Toward a Molecular Theory of Homogeneous Bubble Nucleation: II

Article pubs.acs.org/JPCB

Toward a Molecular Theory of Homogeneous Bubble Nucleation: II. Calculation of the Number Density of Critical Nuclei and the Rate of Nucleation Korosh Torabi and David S. Corti* School of Chemical Engineering, Purdue University, West Lafayette, Indiana 47907-2100, United States ABSTRACT: In the present paper, we develop a method to calculate the rate of homogeneous bubble nucleation within a superheated L-J liquid based on the (n,v) equilibrium embryo free energy surface introduced in the first paper (DOI: 10.1021/ jp404149n). We express the nucleation rate as the product of the concentration of critical nuclei within the metastable liquid phase and the relevant forward rate coefficient. We calculate the forward rate coefficient of the critical nuclei from their average lifetime as determined from MD simulations of a large number of embryo trajectories initiated from the transitional region of the metastable liquid configuration space. Therefore, the proposed rate coefficient does not rely on any predefined reaction coordinate. In our model, the critical nuclei belong to the region of the configuration space where the committor probability is about one-half, guaranteeing the dynamical relevance of the proposed embryos. One novel characteristic of our approach is that we define a limit for the configuration space of the equilibrium metastable phase and do not include the configurations that have zero committor probability in the nucleation free energy surface. Furthermore, in order to take into account the transitional degrees of freedom of the critical nuclei, we develop a simulationbased approach for rigorously mapping the free energy of the (n,v) equilibrium embryos to the concentration of the critical nuclei within the bulk metastable liquid phase.

I. INTRODUCTION In paper I,44 we introduced an equilibrium embryo definition for the process of homogeneous bubble nucleation based on the two order parameters v (volume of a sphere located at a fixed point within the metastable phase) and n (number of vapor-like particles that fall inside the sphere). The (n,v) equilibrium embryo is defined to be the ensemble of all the configurations that have strictly n vapor-like particles (and no liquid-like particles) inside the spherical volume v and a shellparticle that always coincides with the boundary of v (the shellparticle is defined to be the closest liquid-like particle to the center of the embryo). This embryo definition allows for any configuration to correspond to one and only one point in the (n,v) order parameter space. The free energy surface of the equilibrium embryos in the (n,v) space was generated via several Monte Carlo (MC) umbrella sampling algorithms discussed in paper I.44 The main purpose of introducing an equilibrium embryo definition is to determine the corresponding free energy surface. In turn, this surface allows one to identify the critical nuclei and to estimate their average number density within the bulk metastable phase. Yet, if one is interested in calculating the rate of nucleation, an additional piece of information, and one that is fully dynamic, must also be obtained. As is the case in the transition state theory (TST), one must estimate the number density of the transition states as well as the corresponding kinetic rate factor, which is a measure of how fast the critical nuclei cross over the free energy barrier or evolve toward the © 2013 American Chemical Society

formation of the new phase. In the case of a diffusive barrier crossing, one also needs to determine the transmission coefficient, which is the fraction of trajectories, out of all those trajectories (i.e., pathways along the free energy surface) initially passing over the free energy barrier, that eventually do not return back over the barrier and toward the bulk metastable liquid phase. Employing the equilibrium free energy surface introduced in paper I44 for subsequent calculation of the nucleation rate is the main focus of the present paper. But before proceeding to discuss our proposed approaches for determining the nucleation rate, we first draw the reader’s attention to several key issues that arise when utilizing the framework of the TST for any equilibrium-based molecular theory of homogeneous nucleation, some of which have not been fully considered in the nucleation literature. First, we address the intrinsic difference between a homogeneous nucleation event and a chemical reaction. When dealing with a chemical reaction, the underlying free energy surface describes the equilibrium concentrations of the reactants, various intermediate and transition states, as well as products (i.e., the reaction is reversible and so at equilibrium the net rate of reaction is zero). In homogeneous nucleation, however, the product side of the free energy surface is Received: April 26, 2013 Revised: September 6, 2013 Published: September 10, 2013 12491

dx.doi.org/10.1021/jp404151h | J. Phys. Chem. B 2013, 117, 12491−12504

The Journal of Physical Chemistry B

Article

nonexistent (i.e., the “reaction” is strictly irreversible). This should not be confused with the possibility of having two phases at equilibrium with each other at coexistence conditions. Homogeneous nucleation refers to the process of formation of embryos of the new phase that grow into a macroscopic new phase, and consequently leave the bulk metastable phase. To have this process run in the reverse direction, meaning that macroscopic regions of the stable phase re-enter and phase transfer back into the bulk metastable phase, is irrelevant. Therefore, one is not able to envision any kind of equilibrium for trajectories that travel between the reactants and products for the process of homogeneous nucleation within a metastable phase. In this case, a meaningful free energy surface must represent just the distribution of the embryos within the equilibrium metastable phase. Trajectories in both the forward and reverse direction are only possible within such a distribution (and not between the metastable and stable phases). Second, as opposed to the case of a chemical reaction, the concentration (number per unit volume) of critical nuclei within the bulk metastable phase cannot always be readily calculated from the free energy surface. This problem is more pronounced for the equilibrium embryo definitions that are based on configurations or density fluctuations about a fixed point in space. For these types of definitions, the equilibrium embryos are not entities that can be counted within a given volume of the bulk metastable phase (as some of the degrees of freedom of the embryo particles may be redundantly counted as the embryo is translated throughout the system). Such definitions instead describe events with a certain probability of occurrence at a given point in space. Thus, the question always arises as to what would be the equilibrium concentration of entities corresponding to that embryo within the bulk metastable phase. Estimating the number density of critical nuclei has been a matter of controversy for this type of embryo definition, which has been employed within density functional theory (DFT) models,1 classical nucleation theory (CNT) models,2 and also our (n,v) embryo model presented in paper I.44 There has been an ongoing debate over the way in which the free energies for such models3−7 should be mapped into number densities within the bulk metastable phase. Several ad hoc and after-the-fact solutions to this problem have been suggested and critiqued3,8−13 in the literature. As opposed to CNT and DFT models that do not properly account for the translational degrees of freedom of the nuclei, there are several molecular-based approaches, mainly devised for the process of homogeneous droplet nucleation, in which using the center of the mass of the embryo as the reference point14,15 or introducing a particular liquid cluster definition16 does result in embryos that are indeed countable entities. Models based on the center of mass or a nearest-neighbor liquid cluster definition cannot, however, be readily applied to bubble nucleation, as they are inconsistent with various configurations that correspond to a low-density vapor bubble. Finally, finding a reliable reaction coordinate for the homogeneous bubble nucleation is not always a straightforward task. Complicating matters further is the acknowledgment that the rate of nucleation is dependent on the choice of order parameters. Therefore, any uncertainty in the chosen order parameters leads to questions about the reliability and accuracy of the calculated nucleation rate. Moreover, since barrier crossing in a nucleation process is highly diffusive, estimating the transmission coefficient includes a relatively large error.17

In this paper, we address the above issues in the context of the (n,v) embryo definition that we proposed in paper I,44 while also discussing how the corresponding free energy surface can be used, along with additional dynamical information, to calculate the rate of the homogeneous bubble nucleation. Here, in section II, we introduce an additional constraint within our MC simulation algorithms in order to prevent the sampling of configurations that are outside of the configuration space of the superheated liquid phase. In this way, we make sure the free energy surface that corresponds to the number density of the nuclei within the metastable phase is generated based on the configurations that strictly belong to the configuration space of the equilibrium superheated liquid phase. This in turn paves the way toward the proper calculation of the critical nuclei concentration within the bulk metastable phase. In section III, we explain why the probability of a given (n,v) embryo (which is defined with reference to a fixed point in space) is equal to the average volume fraction available for insertion of the corresponding (n,v) sphere within the bulk metastable phase. In the same section, we devise a MC simulation algorithm which enables us to calculate the number density of the critical nuclei as countable entities based on the information provided by the (n,v) embryo free energy surface. In section IV, we introduce a molecular dynamics (MD) simulation-based method to calculate the kinetic rate factor as a measure of how fast the critical nuclei move in the forward direction toward the formation of a macroscopic bubble. Instead of using a particular reaction coordinate and calculating the time derivative of the reaction coordinate of the critical nuclei, we use a large number of MD simulation trajectories initiated from uncorrelated critical configurations to calculate the average lifetime of the critical nuclei. The product of the number density of critical nuclei and the kinetic rate factor yields the rate of homogeneous bubble nucleation, which is expressed in terms of the number of bubbles formed per unit time per unit volume of the superheated liquid phase. In section V, we acknowledge that the (n,v) order parameter is not a unique choice for an equilibrium embryo definition (although in light of CNT, it is an intuitive starting point for studying bubble nucleation). Other embryo definitions would, of course, lead to different free energy landscapes. Ultimately, the quality of an equilibrium embryo definition depends on its ability to appropriately describe the inherently dynamical process of homogeneous nucleation. Hence, there are some restrictions on what constitutes a good set of order parameters. The critical nuclei predicted by the chosen set of order parameters in an equilibrium-based model should represent the transitional configurations of the nucleation process (i.e., configurations with a committor probability of approximately one-half when the system is allowed to dynamically evolve.18 In this work, we define the committor probability as the probability that the dynamical trajectories initiated from a given configuration fall back to the basin of the metastable phase). To that end, and with recourse to a previous dynamical study,19 we also discuss why embryos defined with our chosen “local order” parameter do appear to be relevant to the dynamical process of homogeneous bubble nucleation. In contrast, a “global order” parameter20 that results in a system spanning and ramified critical embryo does not seem to be relevant to the transition configurations within the superheated liquid phase. Finally, conclusions are presented in section VI. 12492

dx.doi.org/10.1021/jp404151h | J. Phys. Chem. B 2013, 117, 12491−12504

The Journal of Physical Chemistry B

Article

II. LIMITING THE CONFIGURATION SPACE OF THE METASTABLE LIQUID PHASE A. Dynamical Limits of a Metastable Phase Space. As the equilibrium metastable phase does not correspond to the lowest free energy for the given external conditions, there are various microstates (e.g., configurations containing a large region of vapor-like density spontaneously formed within the superheated liquid) from which a dynamical trajectory quickly continues with a phase transition toward the stable vapor phase (i.e., these microstates essentially have no chance of staying, or committing, to the metastable liquid phase). Such microstates are not members of the phase space of the equilibrium metastable liquid. Consequently, only a limited region of phase space, the so-called basin of attraction (Figure 1), defines a

within the embryo. Since the vapor particles confined inside the impermeable boundary (or embryo volume) attract the surrounding liquid particles, for larger values of n, larger embryo volumes can be explored before the system becomes unstable. The dynamical limit of the metastable phase space should, however, be determined directly from dynamical trajectories of the metastable system. The (n,v) limit of stability within the MC simulation is, by default, determined within configuration space only, and is not based on any dynamical information. This clean separation of configuration and momentum space is implicitly invoked when utilizing the MC algorithm, such as was done for the MC sampling method developed in paper I44 for obtaining the free energy surface of our chosen (n,v) embryo definition. Configurations belonging to small nonbarrier (n,v) embryos, when allowed to evolve dynamically (via separate MD simulation and after the removal of the impermeable boundary used to create the equilibrium (n,v) embryo19), almost always return to the liquid phase. But for configurations that are near the barrier, the previous MC method may label a given (n,v) embryo as being unstable (stable) even though various configurations belonging to that embryo, when placed within a dynamical simulation, yield nonzero (zero) committor probability. As a matter of fact, for large values of n, the MC simulation may allow for the sampling of configurations that should be considered outside of the metastable basin (i.e., when the impermeable boundary is removed and initial momenta are assigned to all the particles, essentially all the dynamical trajectories that are generated quickly trigger a phase transition). And on the other hand, for small values of n, as we get closer to the transitional region of the order parameter space, the MC simulation might render a certain (n,v) embryo as unstable, while it contains various configurations that would, in a dynamic simulation, still be inside the metastable basin. In this case, some important transitional states are missing from the equilibrium free energy surface. The above discussion indicates that our embryo definition must be further modified. Specifically, an additional constraint must be imposed to allow us to expand the free energy surface up to the end of the transitional region, while at the same time preventing the sampling of configurations that do not belong to the basin of the metastable phase. This refinement of our embryo definition requires, however, input from dynamical simulations (and so its introduction was delayed until now, and not included in paper I,44 which focused solely on the equilibrium-based embryo definition). Determining the dynamical limit of the metastable phase, thereby identifying those states to include in the metastable basin, is, however, a computationally prohibitive task. For each given configuration (initial particle positions), separate dynamical trajectories should be followed after being initiated with a broad range of initial particle velocities. If most of the trajectories initiated from a given configuration remain within the metastable basin, that configuration is included in the equilibrium phase space of the metastable liquid; otherwise, it is deemed to be outside the metastable phase space limit. Such a fully dynamical sorting of the configurations is computationally intractable. In the following section, we introduce an approximate method to modify the previously obtained (n,v) free energy surface such that only (or approximately so) those configurations within the dynamical limits of the metastable liquid are considered.

Figure 1. Two-dimensional representation of the configuration space of a hypothetical metastable system with degrees of freedom q1 and q2. The shaded area (blue) for which the committor probability is between zero and one are transition states that contain the critical nuclei of the new phase within the metastable system. The phase space of the metastable system is limited to the region in which the committor probability is greater than zero.

thermodynamic equilibrium for a metastable phase. In that sense, a metastable phase is considered to be a uniform phase at thermodynamic equilibrium rather than a heterogeneous phase, including both the liquid and regions of growing vapor phase (bubbles). As soon as a trajectory steps outside the metastable basin, where the committor probability is zero, it enters a different region of the phase space where dynamical trajectories quickly follow pathways that depart ever further from the metastable basin via the continual growth of the new phase. Since the purpose of generating a free energy surface is to estimate the concentration of the critical nuclei within the bulk metastable liquid (transitional states with a committor probability of around 0.5), we should determine the free energy surface based on only those microstates that are members of the phase space of the metastable liquid. Including the configurations with committor probability of zero within the free energy surface inevitably leads to overestimation of the critical nuclei concentration within the metastable liquid. The dynamical limit of the metastable phase space is related to the (n,v) limit of stability that emerges in the free energy surface developed in paper I44 and also within previously suggested models that used similar isothermal−isobaric MC algorithms.21,22 As is apparent in Figure 3 of paper I,44 the free energy surface can only be extended up to certain embryo radii, beyond which the NPT MC simulation becomes unstable (the size of the simulation box diverges from its average value and the metastable liquid surrounding the embryo can no longer maintain its liquid-like density). Such a limit of stability depends on the number of vapor-like particles n confined 12493

dx.doi.org/10.1021/jp404151h | J. Phys. Chem. B 2013, 117, 12491−12504

The Journal of Physical Chemistry B

Article

B. Finding the Proper Constraint That Limits the Free Energy Surface to the Metastable Basin. As explained above, the most precise limit of the metastable state is defined by the committor probability analysis. Examining the committor probability for every single microstate is not possible, however, so we instead look for some additional nondynamical characteristic of the equilibrium embryo that correlates with the corresponding committor probabilities. While there is some loss of rigor by doing so, as we again consider only the configuration space of the system, this method is still computationally within reach. Furthermore, a set of subsequent dynamic simulations suggests that the resulting free energy surface is much closer to containing only those states within the phase space of the metastable liquid. We now consider the size of the largest connected volume within the simulation cell that is empty of liquid-like particles (in short, we call it the void volume, although it may contain vapor-like particles). For states near the free energy barrier, this void volume will invariably contain the spherical volume v, and may extend beyond the embryo’s impermeable boundary. It is has been demonstrated that there exists a one-to-one correspondence between the bubble size and the committor probability.23,24 Therefore, we assume the size of this void space is the appropriate indicator of whether or not a given configuration, when placed within a dynamic simulation, will lead to a phase transition. By limiting its value during the MC simulations, we avoid sampling configurations with a zero committor probability and at the same time we are able to reach those nonzero committor probability states that previously could not be sampled. The definition of the void space volume and a procedure to calculate its value for a given configuration within the simulation cell is as follows. Similar to what has been done in the bubble nucleation study by Wang et al.,23 we divide the whole space of the simulation box into a fine mesh of cubical subcells and identify the subcells that do not contain any liquidlike particles and label them as void subcells (see paper I44) for the definition of a liquid-like particle). We use a relatively small mesh size of 0.5σ (where σ is the Lennard-Jones particle diameter) which clearly captures the outline of the void space. In order to exclude the interstitial void space between the liquid particles, we then label all 26 subcells surrounding a liquid-like subcell as liquid-like. With this recipe, we are able to reasonably estimate the void volume formed inside the MC simulation box due to presence of the impermeable boundaries representing the (n,v) embryos (see Figure 2). However, the question remains as to what should be the limiting value of the void volume in order to keep the MC simulations within the limits of the metastable basin. To help provide an answer, we generated histograms of the void space volumes sampled within MC simulations (the MC algorithm with impermeable spherical boundaries was introduced in paper I44) containing several (n,v) embryos suspected of being located near to or inside the transitional region. As mentioned earlier, the previously obtained free energy surface most likely will not extend into the transition region for embryos with a small number of particles (n < 5). For example, as seen in Figure 3 of paper I,44 the limit of stability for a cavity (n = 0) appears at λ ≈ 2.6 (where λ is the radius of the (n,v) embryo, =4/3πλ3). When sampled configurations at this cavity radius are placed within a series of MD simulations,19 almost all these starting configurations yield committor probabilities of one. But for larger values of n, the NPT MC simulation used to generate

Figure 2. Schematic of the void space volume (encompassed by the dashed line) formed about the impermeable spherical boundary (solid line) that defines the volume of the equilibrium embryo. The filled circles represent the liquid-like particles within the metastable liquid phase, while the empty circles are the vapor-like particles. Note that as we impose no constraint on the particles outside the embryo there is always a chance that a vapor-like particle appears outside of the impermeable boundary.

the free energy surface is stable enough so that larger embryos can indeed sample some configurations that appear to reside outside of the basin of the metastable phase. Figure 3 shows the

Figure 3. Histogram of the number of void subcells for the MC umbrella sampling simulation of an embryo of n = 10 and λ = 3.2 generated within the Lennard-Jones metastable liquid at Pσ3/ε = −0.33 and kBT/ε = 0.8.

histogram of the void space volume for an equilibrium embryo of n = 10 and λ = 3.2 formed within the Lennard-Jones (LJ) superheated liquid at the same metastable conditions of paper I (ref 44; kBT/ε = 0.8 and Pσ3/ε = −0.33, where kB is the Boltzmann constant and ε and σ are LJ energy well depth and particle diameter, respectively). For this embryo, when randomly selected configurations are placed within the MD simulation, almost all the starting configurations have committor probabilities of less than 10%. In other words, almost all the sampled configurations for this large embryo do not belong to the basin of the metastable phase. As shown in Figure 3, the majority of the configurations sampled for this embryo contain more than 2200 void subcells. Thus, we take the limit of the void space volume to be 2200 subcells. Below, we report on the sensitivity of the (n,λ) free energy surface to this value of the imposed limit. C. Applying the Configuration Space Limit to the (n,v) Free Energy Surface. Assuming the void space volume limit is an effective tool to avoid sampling configurations that are not members of the metastable basin, we now have to implement that limit within the simulation procedures that we developed to generate the (n,v) free energy surface. In paper I (ref 44; Section IV.B), we introduced a MC umbrella sampling 12494

dx.doi.org/10.1021/jp404151h | J. Phys. Chem. B 2013, 117, 12491−12504

The Journal of Physical Chemistry B

Article

Our present (n,v) embryo definition ensures a unique one-toone mapping of the configuration space onto the order parameter space. As such, P(n,v) is properly normalized, that is ∑n∞= 0∑i ∞= 1P(n,λi) = 1, but P(n,v) is the probability of bubble formation about a fixed point in space, designated as the geometric center of the (n,v) embryo. With this definition, embryos are defined as events that may or may not happen at a certain point within the metastable liquid phase. However, while keeping all the particles for a given configuration fixed (both inside and outside the embryo), the geometric center of the embryo may be translated a small distance without generating a new (n,v) embryo. In other words, the same (i.e., redundant) configuration has been assigned more than once to a given (n,v) embryo.6,7,10−12,25 The current bubble definition therefore does not produce countable entities in the sense that it is meaningless to ask how many distinct embryos with a particular n and v are contained within a given snapshot of the bulk metastable phase. The removal of these translational redundancies, or degrees of freedom, is relatively straightforward when dealing with dense liquid clusters, whereby the center of mass of the cluster is forced to coincide with the geometric center of the embryo volume.14,15 Doing likewise for bubbles or low-density embryos is more problematic. So far, the free energies of formation that we have determined are directly related to the probability of realizing a particular (n,v) embryo at a given point in space, as P(n,v) = exp(−W(n,v)/kBT). We therefore need to devise a systematic method of converting this information provided by the free energy surface into the number density or concentration of the critical nuclei (as countable entities) within the bulk metastable phase. A. An Alternative Interpretation of the (n,v) Equilibrium Embryo Free Energy Surface. From another point of view, the probability of finding a given (n,v) embryo at fixed point within the bulk metastable phase, P(n,v), can be interpreted as follows. Consider the (n,v) embryo definition as a lens that can be moved around within a given snapshot of the bulk metastable phase in order to find all the points that are the center of that particular (n,v) embryo. With this thought experiment, one realizes that P(n,v) is also equal to the average value of the available volume for successful insertion of the geometric center of the given (n,v) embryo per unit volume of the bulk metastable liquid. Non-redundant counting of the configurations within the (n,v) embryo definition ensures that there are no overlaps between the volumes available for different (n,v) embryos within a given snapshot of the system, and therefore, the summation of the available volume fractions for all possible (n,v) embryos is merely 1. Hence, by considering the volume available for insertion of critical (n,v) embryo that does not generate a distinct nucleus (defined as a countable entity), we can determine how to modify P(n,v) in order to correct for the translational degrees of freedom not accounted by the embryo definition. As such, we will be able to convert P(n,v) into the number density of the critical nuclei. In the following section, we demonstrate how to estimate the average volume available for insertion of the critical (n,v) embryo (event) within a single critical nucleus (entity). B. An MC Simulation Algorithm to Map (n,v) Probabilities into a Number Density. What enables us to calculate the number density of the critical nuclei as countable entities is the fact that regions of the bulk metastable phase that can accommodate the center of the critical (n,v) embryo are

procedure in which we employed impermeable spherical boundaries to sample different (n,v) equilibrium embryo configurations. For such simulations, as we march toward larger radii, the chance of sampling configurations outside the limits of the metastable phase increases. In order to generate a free energy surface exclusively for those states within the metastable basin, we use the same MC umbrella sampling procedure, but now a particle move or a change in the simulation cell volume is rejected if the resulting total number of the void subcells exceeds the limit of 2200. The updated free energy surface is provided in the next section (see Figure 7). Figure 4 shows the free energy curve for a cavity (n = 0) with the void volume limit of 2200. As can be seen, the free energy

Figure 4. Free energy surface of embryos with n = 0 for three different void volume limits of 2000, 2200, and 2600 subcells, generated within the Lennard-Jones metastable liquid at Pσ3/ε = −0.33 and kBT/ε = 0.8. All three curves are extended up to transitional region. The curve corresponding to the 2600 subcell limit is further extended to show how the embryo probability P(n,λ) rapidly drops to zero for embryos beyond the transitional region.

sharply rises as we go beyond the transitional region. As a direct consequence of the imposed limitation on the configuration space, the probability of forming an (n,v) embryo outside of the limits of the metastable basin rapidly drops to zero (as P(n,v) = exp(−W(n,λ)/kBT), where v = 4πλ3/3). Also as evident in Figure 4, the calculated free energies are not very sensitive to the value of the void volume limit imposed, especially for those embryos of interest found before and inside the transitional region.

III. CALCULATING THE NUMBER DENSITY OF THE CRITICAL NUCLEI In order to express the rate of homogeneous nucleation in terms of the number of bubbles formed per unit time per unit volume, one first needs to estimate the average number density or concentration of the critical nuclei (defined as countable entities) within the bulk metastable phase. As mentioned in the Introduction section, any model that resorts to a fixed point in space within its embryo definition (e.g., CNT and DFT) neglects the transitional degree of freedom of the nuclei. Therefore, within such models, the number density of the critical nuclei cannot be so readily expressed in terms of the free energy surface of the embryos formed about a fixed point in space. 12495

dx.doi.org/10.1021/jp404151h | J. Phys. Chem. B 2013, 117, 12491−12504

The Journal of Physical Chemistry B

Article

the number of vapor-like particles included in the vapor embryo. Based on the above assumptions, we conclude that the transitional configurations within the metastable phase correspond to equilibrium embryos of radius λ‡ = 3.2 and arbitrary values of n. Therefore, by moving around a lens of radius λ‡ (regardless of the value of n) within a given snapshot of the bulk metastable liquid, we are able to identify critical nuclei within the system and calculate the value of the volume quantum corresponding to the critical radius λ‡. Once the volume quantum ⟨v⟩ is known, we can calculate the number density of the critical nuclei ρ‡ as follows

scarce, uncorrelated and spatially separated from each other in space. This is equivalent to assuming that the chance of having two or more correlated critical nuclei forming close to each other in space is negligible. This is a valid assumption for nucleation within a metastable phase not too close to the spinodal region; as found below, the extremely low number density of the critical nuclei that we calculate for the conditions of our study validates this assumption. The direct consequence of the above assumption is that the total available volume for insertion of the geometric center of the critical (n,v) embryo within a given snapshot of the system comes in the form of locally concentrated, and hence countable, “lumps” at distinct and separate regions within the bulk metastable phase (see Figure 5). Each of these volume lumps is associated with a

∞

‡

ρ =

∑n = 0 P(n , λ‡ = 3.2) < v>

(1)

To calculate the volume quantum ⟨v⟩, we developed a new NPT MC simulation procedure that employs a spherical boundary of radius 3.2, fixed at the center of the simulation box, which is permeable only to vapor-like particles. In this way, we allow the system to freely sample all of the configurations available to embryos of radii 3.2 and any number of vapor-like particles according to their statistical weight. We also use the void subcell limit of 2200 (discussed in section II) to make sure the sampled configurations belong to the metastable phase configuration space and to keep the NPT MC simulation stable. For this simulation, we used 5000 particles, and after 30 000 equilibration cycles, we started saving uncorrelated configurations every 300 MC cycles (each cycle consists of 5000 particle move attempts and one box volume move attempt). Within a given configuration, we randomly pick a point near the center of the simulation box and evaluate n and λ for that point (in this way we are moving around and labeling different points by our (n,λ) lens). This is done by first finding the shell particle (closet liquid-like particle to the center) and then counting the number of particles that fall inside a sphere centered at the randomly chosen point with a radius that is defined by the shell particle (see Figure 6). For each saved configuration, we repeat this procedure for 50 random points and are all located inside a smaller spherical test-volume of radius 2.0 placed at the center of the simulation box, which coincides with the center of the spherical volume of radius λ‡ = 3.2. (In this way, we are allowing for the translation of the center of the embryo to nearby points.) We repeat this procedure for 10 000 uncorrelated configurations and count the number of times for which the random test points happen to be the center of an embryo of radius 3.2 (or to be more exact a radius between 3.15 and 3.25, defined by the shell particle). For a simulation run, we counted 50 418 embryos out of the 500 000 random points. Therefore, the volume quantum available for critical nucleation sites (n = 0,1, ..., and λ = 3.2) is equal to ⟨v⟩ = (50 418/500 000)(4π/3)(2.0)3 = 3.379, where (4π/3) (2.0)3 is the volume of the spherical test region at the center of the embryo. Here, we would like to mention that the test-volume should be large enough to contain the whole volume available for the center of the (n,λ = 3.2) configurations. In other words, the chance of having an available point that falls outside the test volume should be (essentially) zero. To make sure that our test volume of radius 2.0 meets this criterion, we generated a histogram of the points that accept a sphere of radius 3.2 or greater based on their distance from the center of the embryo. The generated histogram fell to zero well before a distance of 2.0 was reached. This result indicates that the test volume was

Figure 5. Schematic of the volume available (shown in red) for the insertion of the geometric center of the critical (n,v) embryo (depicted with the dashed line in the upper right magnified picture). The available volume for the critical (n,v) embryo comes as locally concentrated, and so countable, lumps at distinct regions within the bulk metastable phase. As shown in the lower right picture, this assumption does not hold for embryos of much smaller radii, in which the volume available spans various interconnected regions in space.

distinct critical nucleus (as a countable entity). We define the average volume available, within a single critical nucleus, for the center of a critical embryo (an event) as the volume quantum ⟨v⟩. To identify the critical nuclei within the bulk metastable phase, first we have to find the values of n and λ that correspond to the transitional configurations within the system. Previously, via a MD simulation committor probability analysis, we showed that at the conditions of our study (kBT = 0.8 and ε Pσ3/ε = −0.33) a cavity of λ = 3.2 has a committor probability of about one-half.19 This means that the configurations sampled with a void impermeable boundary of λ = 3.2 correspond to transitional configurations of the metastable liquid phase. At this point, we assume that the radius of the embryo is what determines whether the configuration is transitional or not. In this way, we are assuming that the number of vapor-like particles inside the embryo (which are relatively separated from one another and the surrounding liquid-like particles) does not drastically change the committor probability of the embryo. Later, we show this assumption is consistent with the results of our work. Additionally, the results reported by Meadley et al.24 shows that there is a one-to-one correspondence between the vapor-like volume and the committor probability, regardless of 12496

dx.doi.org/10.1021/jp404151h | J. Phys. Chem. B 2013, 117, 12491−12504

The Journal of Physical Chemistry B

Article

fit, which is used for interpolation between the data generated for different values of n. The data confirms that ∑n∞= 0W(n,λ)/ kBT can be safely truncated at n = 15 and shows that the probability has a maximum around n = 7. These data are also consistent with what was found during the simulation for obtaining the volume quantum. In this simulation, in which the number of particles inside the embryo was allowed to vary, a resulting histogram of n indicates that n = 7 was the most likely embryo (see Figure 8).

Figure 6. Schematic representation of how the average volume quantum ⟨v⟩ within a critical nucleus is determined. For a given a snapshot of the system, randomly chosen points (depicted as a cross mark) are selected inside a spherical test volume of radius 2.0 (shown as the dashed area). The test volume is located at the center of the spherical boundary of radius 3.2 (depicted as the double-lined circle), which is held fixed at the center of the simulation box. This boundary allows only for the vapor-like particles to enter. For each test point, we determine if there is at least one liquid-like shell-particle within a spherical shell of thickness 0.1 and radius 3.2 centered at the test point, and if any liquid-like particle is found within a distance of 3.2 from the test point. As an example, the off-centered dashed circle represents a successful embryo of n = 2 vapor-like particles and radius λ = 3.2 defined by the shell-particle.

Figure 8. Histogram of the number of vapor-like particles inside an embryo with the critical radius of λ = 3.2. This plot shows the number of times that each value n was visited out of 50 418 successful embryo insertions of the algorithm explained in section III.B.

From the data presented in Figure 7, we find that n=0

large enough, and so the volume quantum was calculated accurately. Now that we have estimated the volume quantum ⟨v⟩, the only term that needs to be evaluated is ∑n ∞= 0P(n,λ = 3.2) (which is equal to ∑n∞= 0 exp[−W(n,λ)/kBT]). This summation can be evaluated using the modified free energy surface generated after imposing the proper metastable configuration space limit introduced in section II. Figure 7 shows the modified free energy surface extended up to a radius of λ = 3.2 for several values of n between 0 and 15 (values of P(n,λ = 3.2) for n larger than 15 are negligible). The inset in Figure 7 shows a plot of W(n,λ = 3.2)/kBT along with a third order polynomial

∑ P(n , λ = 3.2) = 1.9415 × 10−24 (2)

∞

With the value of ⟨v⟩ = 3.379 calculated above, the dimensionless number density of the critical nuclei is therefore equal to ∞

∑n = 0 P(n , λ = 3.2)

‡

ρ =

⟨v⟩

= 5.7458 × 10−25

(3)

All of the above calculations are based on dimensionless quantities for the metastable LJ liquid at T(kB/ε) = 0.80 and P(σ3/ε) = −0.33. Using the LJ parameters of σ = 3.405 Å and ε/kB = 125.2 K for argon,26 we have ρ‡ = 1.46 × 104 m−3 or about 14 critical nuclei per one liter of the superheated liquid argon at T = 100 K and P = −145 bar. Such a small number density of the critical nuclei justifies our earlier assumption that the critical nuclei are uncorrelated, that are located far away from each other within the metastable liquid.

IV. DETERMINING THE RATE OF HOMOGENEOUS BUBBLE NUCLEATION The rate of homogeneous nucleation is defined as the number of nuclei of the new phase that successfully move beyond the nucleation free energy barrier, per unit time, per unit volume of the bulk metastable phase. Based on transition state theory (TST), the rate of a reaction could be written as the net flux of reactants in the forward direction along the reaction coordinate.27 Accordingly, the nucleation rate J could be written as J = k × ρ‡

Figure 7. Free energy surface of formation of (n,v) embryos within the LJ superheated liquid at a reduced temperature of kBT/ε = 0.8 and a reduced pressure of Pσ3/ε = −0.33 extended to the critical radius of 3.2 consistent with the metastable phase limit introduced in section II. The inset is a plot of W(n,λ = 3.2) kBT for different values of n. The solid line represents a 3rd order polynomial fit to the data.

(4)

‡

where ρ is the density of nuclei at the critical value of the reaction coordinate and the forward rate coefficient k is the average velocity of the trajectories in the forward direction at the critical value of the reaction coordinate multiplied by the 12497

dx.doi.org/10.1021/jp404151h | J. Phys. Chem. B 2013, 117, 12491−12504

The Journal of Physical Chemistry B

Article

absence of the spherical boundary and initial velocities were randomly assigned to all particles based on a Maxwell− Boltzmann distribution corresponding to the temperature of the system. The MD time step for these simulations was chosen to be 0.002 σ(m/ε)1/2, where m is the mass of a LJ particle. Trajectories started from these initial configurations were found to either fall back to the basin of the metastable liquid or to instigate a phase transition. We monitored the overall density of the simulation box to identify in which direction a trajectory was moving. Each MD simulation was terminated as soon as the overall density reached the value of 0.74 (density of the bulk metastable phase) or became smaller than 0.68 (the value at which the system was always found to definitely move toward phase transition, as was also found in ref 19). Figure 9 shows a plot of the overall density versus MD time step for 10 different MD simulation runs, out of which 4 were

transmission coefficient which is the fraction of the forward trajectories crossing the critical value that do not cross back into the metastable basin. Based on the above idea, ten Wolde and Frenkel developed a computer simulation algorithm to estimate the rate of droplet nucleation in a subcooled LJ vapor phase.17 They took the number of particles in the largest liquid cluster of the system as their reaction coordinate. Using a modified form of the Bennet-Chandler scheme28 and by following 600 trajectories that were started from 300 independent configurations near the top of the free energy barrier, they calculated a transmission coefficient to be on the order of 10−2. Such a small value for the transmission coefficient is a signature of a highly diffusive barrier crossing. Other studies also support the idea of a highly diffusive barrier crossing mechanism for homogeneous nucleation.17,29,30 In this regime, one can use Kramers’ theory to express the rate coefficient in terms of the diffusivity of a reactive trajectory at the top of the free energy barrier.31 (Diffusivity is related to the mean squared displacement of a trajectory starting at the top of the free energy barrier.) This can be done using either a simple Einstein model or a more complicated Smoluchowski equation, which takes into account the shape of the free energy surface.32 The TST-based models mentioned above necessarily require a predefined reaction coordinate. The shape of the free energy surface is directly dictated by the choice of the embryo definition and a different set of reaction coordinates will result in different embryo distributions at the top of the barrier.33 Systematic methods have been devised to identify the (apparently) relevant reaction coordinates for complex systems.18,34 However, these methods are in general computationally expensive when applied to a system that has a very large number of degrees of freedom and for which generating the reactive trajectories are costly. Therefore, most of the nucleation studies for such systems are based on assuming a priori the reaction coordinates. Although the approach we have been following so far utilizes a free energy surface based on the chosen set of order parameters n and v, we may still determine the kinetics of the bubble nucleation without calculating the velocity of phase transferring trajectories along a predefined reaction coordinate. To do so, we assume that the transitional nuclei instigate a phase transition on a completely random and uncorrelated basis. In other words, as the barrier crossing is highly diffusive, we may assume, for all the critical nuclei, there is an average probability of phase transition per unit time which is simply the forward rate coefficient k as defined in eq 4. In this section, we present a simulation-based procedure to calculate k. If the assumption that the phase transition is randomly triggered by each critical nucleus is valid, the survival probability of a phase transferring critical nucleus, Ps(t), or the probability that a critical nucleus has not yet instigated a phase transition by time t should be of the form e−kt. To test this relation, we followed the time evolution of 200 uncorrelated snapshots representative of the critical nuclei. In the previous section, we presented a MC simulation procedure that employed a semi-impermeable sphere (only vapor-like particles were allowed to enter) of radius λ = 3.2 to generate critical nuclei configurations. Here, we collect 200 uncorrelated configurations generated during the same MC simulation, and use them as the initial configurations within a subsequent NPT MD simulations performed at the conditions of our system (Pσ3/ε = −0.33 and kBT/ε = 0.80). Similar to what has been done in ref 19, NPT MD simulations were performed in the

Figure 9. Plot of the overall density of the simulation box versus MD time step for 10 different MD simulation runs initiated from 10 uncorrelated transitional configurations (simulation outputs are included in consecutive order, and do not represent one single simulation run). The lifetimes of the phase transferring trajectories are each marked with its corresponding Δt. The horizontal dashed line marks the bulk metastable liquid density.

found to initiate a phase transition. Among 200 separate MD simulation runs, we observed 56 phase transferring trajectories. (The remaining 144 trajectories fell back into the basin of the metastable liquid, i.e., no phase transition was observed. This implies that the committor probability of our transitional configurations generated with the spherical boundary of λ = 3.2 is about 72%.) In order to calculate the survival probability as a function of time, we determine the lifetime of each of the 56 phase transferring trajectories. We define the lifetime to be the time it takes for a trajectory to reach the overall density of 0.68. Making use of the lifetimes of the 56 phase trajectories that led to phase separation, we calculate the survival probability Ps(t) as the fraction of those trajectories that have not yet reached the overall density of 0.68 by the given time t. Based on the assumption of a first-order reaction or randomly instigated phase transitions for the critical nuclei, the survival probability should display an exponential decay. As evident in Figure 10, the data of survival times appears to deviate from exponential behavior only for the initial 10 000 MD time steps. This deviation can be attributed to the way in which the lifetimes are calculated. As a reminder, we have defined the lifetime to be the time it takes for a trajectory to reach the overall density of 0.68. This density seems to be beyond the limits of the configuration space of the metastable phase where the chance of the system going back to the basin of the liquid phase is (essentially) zero, but the embryo begins within the transitional region, where there is roughly a 50% chance that a trajectory will initiate a 12498

dx.doi.org/10.1021/jp404151h | J. Phys. Chem. B 2013, 117, 12491−12504

The Journal of Physical Chemistry B

Article

The factor of 0.28 included in the rate coefficient is one minus the committor probability of the critical nuclei at λ = 3.2. Although the calculated nucleation rate should in principle depend on the chosen value of the critical radius, the above formulation for the rate of nucleation has the advantage of allowing for some freedom in choosing the radius of the critical embryo used in the rate calculations. Note that finding the radius at which the committor probability is exactly 50% requires the generation of a large number of MD trajectories at different radii, and so is a computationally expensive task. However, as long as the committor probability of the transitional ensemble is around 50% (and not too close to neither zero nor one, which would yield a large numerical error), choosing a slightly smaller (larger) radius results in a larger (smaller) number density of the critical nuclei while at the same time generating a smaller (larger) probability of a phase transition (i.e., one minus the committor probability). The product of these two quantities is what ultimately appears in the rate equation, so these two effects tend to cancel each other out. For more explanation on this, we refer the reader to the work of ten Wolde et al.17 where a similar relation between the size of the critical nucleus and the transmission coefficient holds within their TST-based model of homogeneous droplet nucleation. Based on the forward rate coefficient estimated above and the number density of the critical nuclei obtained in the previous section, we calculate the dimensionless nucleation rate to be J = k × ρ‡ = 0.028 × 5.7458 × 10−25 = 1.6 × 10−26. As a reminder, this rate is calculated for the LJ liquid at Pσ3/ε = −0.33 and TkB/ε = 0.80. At this temperature, the binodal and spinodal pressures are 0.006 and −0.6, respectively.37 This corresponds to the system being about halfway within the metastable region. For comparison, the nucleation rate calculated based on the CNT formulation of Blander and Katz38 for the LJ liquid at the same conditions is on the order of 10−29, which is 3 orders of magnitude smaller than our results. There is a great uncertainty associated with the CNT prediction due to the high sensitivity of the nucleation rate to the value of the surface tension and the unavailability of accurate values of the surface tension for such high curvatures (that is, for molecular-sized bubbles). The value of the dimensionless surface tension used for the above CNT prediction is 0.68, which was adapted from the DFT results reported by Corti et al.39 for the critical nucleus at a dimensionless temperature of 0.8 and metastability parameter of 50% (defined based on chemical potential differences at the binodal and spinodal). Using the LJ parameters for argon26 of σ = 3.405 Å, ε/kB = 125.2 K, and m = 6.636 × 10−26 kg, the characteristic time scale corresponds to σ(m/ε)1/2 = 2.11 ps. Hence, the forward rate coefficient of the critical nuclei is 0.028 × (1/2.11 ps) = 9.6 × 1011 s−1 . Substituting this number into eq 4, along with the previously obtained number density of the nucleation sites of ρ‡ = 1.46 × 104 m−3, we estimate the rate of homogeneous bubble nucleation for superheated liquid argon at T = 100 K, P = −145 bar to be J = 4.84 × 107 cm−3 s−1. As a comparison, we reference the work of Wang et al.,23 which reports, to the best of our knowledge, some of the very few available data on the rate of homogeneous bubble nucleation in the LJ liquid. For a slightly superheated LJ liquid at Pσ3/ε = 0.026 and TkB/ε = 0.855, they calculated a dimensionless nucleation rate on the order of 10−14 (or 1019 cm−3 s−1 using the LJ parameters for argon). The saturation

Figure 10. Logarithmic-linear plot of the survival probability (cumulative distribution of the lifetimes) of critical nuclei, where Ps(t) is the probability that a critical nucleus has not yet instigated a phase transition by time t. After the initial relaxation period, the plot shows a linear tail in support of the first-order rate assumption.

phase transition. As a result, there is an initial relaxation time within which the nuclei “diffuse back and forth” around the barrier before they finally “decide” to move forward and toward the defined phase transition boundary (i.e., system density of 0.68). Therefore, Ps(t) includes an initial relaxation period in which the trajectories are diffusing throughout the transitional region of the configuration space until they reach a steady-state distribution. After this initial period, a steady-state regime is reached in which the number of trajectories reaching the detection point (i.e., system density of 0.68) equals the number of trajectories leaving the transitional ensemble in the forward direction. The linear tail of the logarithmic plot of Ps(t), appearing after some initial time period, supports this description (see Figure 10). In a MD simulation study of homogeneous bubble nucleation, Watanabe et al.35 reported a similar behavior for the lifetime of a metastable LJ liquid. They showed that the survival probability of a metastable system follows a Poisson distribution with an additional relaxation time t0 . They calculated the average lifetime of the metastable liquid τ by assuming the following form for the cumulative distribution function of the lifetimes ⎧ 0 ,t < t0 ⎪ F (t ) = ⎨ ⎛ t − t0 ⎞ ⎟ ,t ≥ t ⎪1 − exp⎜ 0 ⎝ τ ⎠ ⎩

(5)

The above function is the same as what fits our data in Figure 10. Furthermore, an analogous effect was also observed in an MD simulation study of droplet nucleation within a subcooled LJ vapor by Chkonia et al.36 At conditions not too close to the spinodal, an initial relaxation period was also found to precede the exponential decay regime of the survival probability. Going back to our calculation of the rate coefficient k, the slope of the straight line that fits survival probability data in Figure 10 (after the initial relaxation period) is 2 × 10−4 (MD step)−1. Therefore, the dimensionless forward rate coefficient of the critical nuclei, which is the probability that the transitional nuclei start a phase transition per unit time, is k = ((2 × 10−4)/ (MD step)) × ((MD step)/(0.002)) × 0.28 = 0.028. 12499

dx.doi.org/10.1021/jp404151h | J. Phys. Chem. B 2013, 117, 12491−12504

The Journal of Physical Chemistry B

Article

liquid at ρσ3 = 0.659 and TkB/ε = 0.868 that was 3 orders of magnitude larger than the CNT prediction at the same conditions. While not at the same temperature and pressure of our study, our nucleation rate is arguably consistent with their results, at least with respect to the deviations from the CNT prediction. In addition, we report the recently performed mean first passage time study for the metastable TIP4P water model,41 which predicted the rate of homogeneous bubble nucleation to be roughly 8 orders of magnitude larger than the CNT predictions.

pressure at this temperature is reported to be 0.046. Depending on the value of the surface tension used, they report the calculated nucleation rate to be between 6 and 21 orders of magnitude larger than what is predicted by CNT. The large rate of nucleation predicted in the work of Wang et al. (for a liquid which is only slightly superheated) might be due to one of the key assumptions implicit in the forward flux sampling (FFS) method that they employed. Within FFS, one assumes that small embryos (those that end up crossing the free energy barrier) detected near the basin of the metastable liquid would individually and independently grow into postbarrier bubbles. This is while these embryos have small free energies of formation. Consequently, they have large concentrations and might be located close to each other within the liquid phase. In fact, with the nucleation rate predicted by Wang et al.,23 and assuming that the kinetic rate constant is, e.g., 10 times larger than what we calculated and the probability of forming a critical embryo is consistent with the lowest value of free energy barrier predicted by CNT (36 kBT as reported in ref 22), the number density of embryos close to the basin of the metastable liquid would be approximately 103 σ−3, which is too large to be physically meaningful (this large number density remains quite high even if we chose a rate constant of about 4 orders of magnitude larger than our estimate and therefore seems to be inconsistent with the assumption of independent and uncorrelated formation of nuclei). Such a large number density implies that the phase transferring embryos inevitably join each other as they proceed forward during the early stage of nucleation. And yet, FFS appears to treat these embryos as still being separate and distinct after they cross the first FFS interface. Consequently, many more bubbles are apparently found to be crossing the free energy barrier than what would be found if one began tracking barrier crossings starting just from a critical nucleus. This is because within FFS the rate is calculated as the product of the number of embryos that cross the first FFS interface per unit time per unit volume and the probability that they continue all the way to the other side of the barrier. Therefore, if on average a given number of small embryos crossing the first FFS interface come together as they climb the free energy barrier and eventually form only one post barrier bubble, the rate predicted by FFS needs to be divided by that number. Furthermore, in this case, the various FFS interfaces were defined via the largest bubble that appears in the simulation cell, the identity of which may change as a trajectory approaches the barrier (and so the originally counted embryo that crosses the first interface may not even be the embryo that crosses the barrier). At the same time, we would like to mention that the above deficiencies are not inherent to the FFS algorithm in general (as in many rare event processes such as a chemical reaction, independency of reactive trajectories is not an issue) but they are rather related to its application to the bubble nucleation process. A recent experimental study40 estimates the time scale of nucleating one bubble within superheated water traped inside a microcavity of volume in order of 104 μm3 to be on the order of 1 μs. With a nucleation rate of 1019 cm−3 s−1 for argon estimated by FFS in ref 22, one predicts this time scale to be around 10 fs. With a 10 fs lifetime, it is essentially meaningless to speak of an equilibrium metastable liquid phase. This is while the temperature is only 10% higher than the saturation value. Finally, we also would like to make note of the brute force MD simulation study of Watanabe et al.,35 who calculated a homogeneous bubble nucleation rate for a superheated LJ

V. CHOICE OF EMBRYO DEFINITION: LOCAL VS GLOBAL ORDER PARAMETERS As with any equilibrium-based nucleation model, we began by suggesting an embryo definition relying upon a set of predefined collective variables or order parameters (selected based on some initial physical insight to the process or similar processes) that seem to be relevant to the process of bubble nucleation. This set of order parameters was modified after some initial examination in order to arrive at an improved representation of those embryos that participate in a nucleation event. While arbitrary, the chosen order parameters and the corresponding equilibrium embryo definition should nevertheless be chosen so as to provide a relevant description of the inherently dynamical process of nucleation. Systematic methods have been devised to search for what appears to be a “best set” of order parameters among a set of candidate variables.18,42 These methods are in general computationally expensive and, to the best of our knowledge, have not been implemented for homogeneous bubble nucleation within a three-dimensional fluid, such as one comprising LJ particles. We do note, however, that such methods require some input from the user, as in what candidate variables are to be considered, and so are not completely free from the implicit assumptions invoked by other approaches. Although our chosen (n,v) order parameter has not been explicitly shown to be the best set of order parameters for bubble nucleation, a relatively simpler committor probability analysis can nevertheless be employed to determine the ultimate dynamical relevance of the chosen set of order parameters. In general, it is desired that all of the embryo configurations which are mapped into the same value of the order parameters have approximately the same committor probability. A sharp committor probability histogram therefore indicates the dynamical relevance of the selected set of order parameters.18 If this is the case, the proposed set of order parameters would be able to identify the local density fluctuations that have the potential of initiating a phase transition within the metastable liquid phase. The transition configurations correspond to the points in the order parameter space that have a committor probability value of about one-half. Our chosen equilibrium embryo definition is based on the order parameters n and v, and so is only concerned with the local configuration of particles at a certain location within the bulk metastable phase. In all the simulations developed in order to generate the (n,v) free energy surface, we made sure that the simulation box was large enough to avoid any possible finite size effects. In this way, the nonuniform density profile of any given embryo strictly remains within the confines of the simulation box; all that the embryo “sees” in its surrounding is the bulk liquid phase. Therefore, the suggested (n,v) model serves as way in which those compact low-density regions (not being any larger than a few particles in diameter) that initiate 12500

dx.doi.org/10.1021/jp404151h | J. Phys. Chem. B 2013, 117, 12491−12504

The Journal of Physical Chemistry B

Article

Figure 11. Snapshots of three different critical nuclei are shown, in which the connected void-space volume and the spherical boundary of critical radius (λ = 3.2 used to generate the related critical configurations are provided. For clarity, the vapor-like particles inside the embryo and the surrounding liquid-like particles are not shown. Units of length for the various coordinate axes are the LJ particle diameter σ.

is also the largest vapor nucleus present within the simulation cell) is free to choose any nonspherical shape. Although several other bubble nucleation studies have likewise predicted the critical nuclei to be locally confined regions of vapor-like density,23,24,35 embryo definitions based on a set of local order parameters are not the only choice for studying homogeneous bubble nucleation. In an earlier molecular simulation study of a metastable LJ liquid that employed a nonlocal order parameter, Shen and Debenedetti20 presented a different picture for bubble formation. In an attempt to avoid making a priori assumptions about the internal structure of an embryo, they chose the overall density of the simulation box (or simply the size of the simulation box) as the appropriate order parameter. They calculated the free energy of

the phase transition within the metastable liquid phase can be systematically identified. Snapshots of some typical critical nuclei are provided In Figure 11, in which both the spherical impermeable boundary (λ = 3.2) and the corresponding connected void-space volume (see section II B) are explicitly shown. This figure reveals that although the void-space volume representing a critical nucleus is not perfectly spherical, it nonetheless maintains a locally compact shape. As is apparent in Figure 11, using a spherical boundary as implemented in the (n,v) embryo definition does not mean that we are restricting the embryos to those of perfectly spherical shapes. For any configuration mapped to a given (n,v) embryo, all particles inside the spherical boundary of volume v should strictly be vapor-like. This is while the corresponding void volume (which 12501

dx.doi.org/10.1021/jp404151h | J. Phys. Chem. B 2013, 117, 12491−12504

The Journal of Physical Chemistry B

Article

nucleation revealed that only compact and locally confined regions of low density could be considered as those fluctuations that ultimately triggered the phase transition within the superheated liquid phase. Second, the global order parameter may fail to capture those density fluctuations that are dynamically relevant to the nucleation process. In our previous work,19 we presented a committor probability analysis of the compact and locally confined cavities generated within a superheated LJ liquid at kBT/ε = 0.8 and Pσ3/ε = −0.33. There, the probability of a cavity instigating a phase transition sharply changed from 0 to 1 upon changing the radius of the cavity from 3.0 to 3.5. This shows that with an embryo definition that considers compact and locally confined regions of low density, we were able to properly identify the transition configurations for bubble nucleation. In the same work, we also performed a committor probability test on a system constrained by the size of the simulation box (without imposing any constraint on the configuration of the particles inside the system), i.e., a global order parameter. We showed that while all the trajectories initiated for a cavity of radius 3.5 instigated a phase transition which corresponds to an initial overall density of 0.68, none of the trajectories starting from those configurations generated by uniformly stretching the system to a density of 0.68 (and without explicitly containing any large cavities) started a phase transition. This indicates that by using only one global order parameter we were unable to capture the dynamical bottleneck of the nucleation process. Kusaka and Oxtoby 43 also demonstrated that in order to properly generate the critical embryo for homogeneous bubble nucleation inside a LJ fluid, one needs to supplement the global order parameter (size of the simulation box) with a second order parameter that drives the system toward the formation of compact and locally confined regions of low density.

the system as a function of the overall density, which was appropriately called a global order parameter, via an NPT MC umbrella sampling technique. A free energy barrier was located, and the critical embryo was assumed to be present within the system at this location. The internal structure of the system was then examined by looking for cavities or regions completely devoid of particles. Upon doing so, they found that as the system climbs the free energy barrier, a large number of small cavities come together where, at the top of the barrier, they form a single, highly ramified cavity that spans the entire simulation box. The critical nucleus that appears in this study with the use of a global order parameter is drastically different from the locally compact embryos that are seen with our equilibrium definition. While the choice of an order parameter for nucleation is ultimately arbitrary, we nevertheless would like to raise some concerns about the use of a global order parameter to describe bubble nucleation in superheated liquids. First, a global order parameter is, by definition, nonlocal. As the overall density of the system is reduced in a uniform manner, such a density change is not confined to a small region of the fluid but rather is propagated throughout the entire bulk liquid via the application of the standard periodic boundary conditions. Consequently, what appears to trigger the phase transition is no longer confined to a small microscopic region. As mentioned before, a similar effect seems to be causing problems within the “boxed MD” simulation method of Meadley et al.24 that employs the same global order parameter of ref 19. They report observation of “bubbles forming connected structures across the periodic boundary conditions of the simulation box” before the critical bubble region is reached. This problem persists for even relatively large systems of 8000 particles. Even for the large system of 12 000 particles where the connected structures are not observed, there still exists the possibility of forming the bubbles near the boundaries of the simulation box where they see their own periodic boundary image. They report a free energy barrier difference of about 7 kBT (4 orders of magnitude difference for the critical bubble concentration) and a 2-fold increase in the critical bubble size by increasing the simulation size from 3000 to 12 000 particles. They do not discuss if the size effects are eliminated for systems larger than 12 000 particles. It is also not clear how the redundancy problem caused by the translational degree of freedom of the embryos3,4,13 is treated within the “boxed MD” simulation algorithm. As explained in section III, in order to calculate the nucleation rate, any equilibrium based model should contain a procedure for converting the free energy values (estimated based on the probability of equilibrium embryos formed within the finite volume of the simulation cell) into the average number of the critical bubbles per unit volume of the bulk liquid phase. For conditions not too close to the critical point nor to the superheated liquid spinodal, where the free energy barrier is at least several multiples of kBT, we expect that the embryos should be uncorrelated and their concentrations (or probabilities of occurrence) should be low. In other words, the density profile generated via the use of the chosen order parameters must approach the bulk liquid value before the boundary of the simulation cell is reached, i.e., the relevant embryos should be confined to local, and not global, regions of the metastable phase. Such embryos were seen in the work of Watanabe et al.35 There, a brute force MD simulation (without resorting to any order parameter) of homogeneous bubble

VI. SUMMARY AND CONCLUSIONS In this paper we demonstrated how the information provided with the (n,v) equilibrium embryo model can be used, with additional modifications based on relevant dynamical information, to calculate the rate of homogeneous bubble nucleation in a metastable liquid. This was done by expressing the rate of nucleation in terms of the product of the number density of the critical nuclei and the corresponding forward rate coefficient. By considering the connected void-space volume (region devoid of liquid-like particles), we showed how the (n,v) equilibrium embryo of critical size serves as a lens for detecting the critical (transitional) nuclei within the superheated liquid. In turn, this lens was used to determine how to convert the probabilities of embryo formation at a fixed point into the number densities of critical nuclei (i.e., the phase transferring nuclei as countable entities). This number density was found to have a small value at the chosen state point of our study (corresponding to about halfway within the metastable region), justifying the needed assumption that the critical nuclei are independent and uncorrelated. We also introduced a MD method for calculating the forward rate coefficient of the transitional nucleation sites. Once a transitional state has been identified, the method assumes that subsequent phase transition events are randomly distributed in time and does not depend on an explicitly predefined reaction coordinate. Within an MD simulation, we calculated the survival times of a large number of critical nuclei (the time required for the critical nuclei to pass over the nucleation barrier and instigate a phase 12502

dx.doi.org/10.1021/jp404151h | J. Phys. Chem. B 2013, 117, 12491−12504

The Journal of Physical Chemistry B

Article

other recently obtained values that were generated by different methods and for conditions close to, relatively speaking, the same location in the metastable region. In addition, we attempted to address several important points that are sometimes ignored within other equilibrium-based theories of nucleation. First, we introduced an equilibrium embryo definition that is capable of correctly identifying and quantifying the transitional bottlenecks of the inherently dynamical process of nucleation. Second, we developed a method of properly relating the equilibrium free energy surface to the number density of the critical nucleation sites within the bulk metastable phase. Last, we accounted for the diffusivenature of the barrier crossing when estimating the rate of nucleation in terms of the number of embryos of the new phase successfully formed per unit time per unit volume of the bulk metastable phase. Finally, one of the novel characteristics of our nucleation model is that the underlying free energy surface is defined based on the configuration space of the metastable liquid phase only. Configurations that have zero probability of committing to the basin of the metastable phase are considered outside the limits of the metastable phase configuration space and therefore are excluded from the free energy surface. This limit is enforced in the MC umbrella sampling simulations via the additional constraint on the void-space volume within the simulation box. As shown in Figures 4 and 7, our suggested free energy surface does not show a maximum in the value of the work, and this is due to the very fact of not including the post barrier density fluctuations in the free energy surface. The nonappearance of a maximum is consistent with what has been demonstrated by Toxvaerd29 via a brute force MD simulation of droplet nucleation within a subcooled LJ vapor. Toxvaerd showed that the cluster distribution of the droplets within the configuration space of the metastable vapor monotonously goes to zero and does not show a minimum as long as the phase transition is not occurring, whereby one can speak of a truly equilibrium metastable vapor phase. Likewise, the embryos included in our suggested free energy surface span the relevant configuration space all the way to the point where the committor probability becomes (essentially) zero and not any further. The critical nucleation sites are defined to be those for which the chance of committing to either side of the barrier is of the same order. Within more conventional equilibrium-based nucleation theories, the critical nucleus is, however, taken to be the embryo that corresponds to the saddle point of (or maximum, in the one-dimensional case) the free energy surface in the order parameter space. In paper I,44 we showed that for complex systems, due to an improper choice of order parameters, there is always a chance that the saddle point of the free energy surface does not indicate the true dynamical bottleneck of the phase transition, as appears to be the case with the global order parameters used for homogeneous bubble nucleation.20

transition). We showed that regardless of a short initial period the cumulative distribution function of the lifetime of the critical nuclei has an exponential shape that can be used to calculate the relevant forward rate coefficient. When combined with the data on the number density and committor probability of the critical nuclei, we estimated the rate of homogeneous bubble nucleation within the LJ liquid at a particular metastable state point. As a summary, the different steps required to calculate the rate of homogeneous bubble nucleation within a superheated liquid at a given temperature and pressure are the following: 1. Generating the (n,v) free energy surface via MC umbrella sampling algorithms developed in paper I44 and subsequent implementation of the dynamical limit of the metastable phase based on the committor probability analysis developed in section II. The committor probability also identifies the transitional region of the free energy surface. 2. Mapping the free energy of formation at a fixed location into the number density of the critical nuclei via the MC algorithm developed in section III. 3. Estimating the rate constant based on the survival times of an ensemble of critical nuclei from the corresponding MD simulation trajectories as explained in section IV. When correcting the free energy surface for the metastable phase limit, we assumed that there was a one-to-one correspondence between the void volume of the embryo (as defined in section II) and the committor probability of the corresponding configurations, irrespective of the number of vapor-like particles within the embryo. With this assumption, we avoid the computationally expensive committor probability calculations for each and every configuration included in the free energy surface. The committor probability values for different void volumes and the insensitivity of the free energy surface to the exact value of the imposed limit (as shown in Figure 4 of section II) strongly suggests that the final results are not greatly influenced by the above assumption. Also, implicit to the number density calculations is the assumption that a single ensemble-averaged volume quantum, as defined and calculated in section III, is a good representation of the average volume available for insertion of a sphere with the critical radius. To ensure the validity of this assumption a large ensemble of uncorrelated critical nuclei must be used in the volume quantum calculations (here we used 50 test insertions for 104 uncorrelated critical nuclei). Finally, the procedure suggested for calculating the rate constant, as introduced in section IV, is valid if a first order kinetics can explain the phase transition of the critical nuclei. This requires the barrier crossing to be highly diffusive, in a way that the lifetime of a critical nucleus in the transitional region is completely random and independent of the trajectory that created it. We note that any of the above assumptions might not be valid for all liquids and in all metastable conditions, though they do appear to be reasonable for the temperature and pressure chosen here. Their validity, however, needs to be carefully checked for other systems of interest. The theoretical and experimental data for homogeneous bubble nucleation rates presently available in the literature is rather scarce and does not correspond to the conditions chosen here, so a meaningful quantitative comparison of our results with existing rate calculations cannot be truly made. Nevertheless, our estimated nucleation rate is comparable to some

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Present Address

Korosh Torabi, Department of Chemistry, Northwestern University, Evanston, Illinois 60208−3113. Notes

The authors declare no competing financial interest. 12503

dx.doi.org/10.1021/jp404151h | J. Phys. Chem. B 2013, 117, 12491−12504

The Journal of Physical Chemistry B

■

Article

(22) Uline, M. J.; Torabi, K.; Corti, D. S. Homogeneous Nucleation and Growth in Simple Fluids. Ii. Scaling Behavior, Instabilities, and the (N,V) Order Parameter. J. Chem. Phys. 2010, 133, 174512. (23) Wang, Z. J.; Valeriani, C.; Frenkel, D. Homogeneous Bubble Nucleation Driven by Local Hot Spots: A Molecular Dynamics Study. J. Phys. Chem. B 2009, 113, 3776. (24) Meadley, S. L.; Escobedo, F. A. Thermodynamics and Kinetics of Bubble Nucleation: Simulation Methodology. J. Chem. Phys. 2012, 137, 074109. (25) Reiss, H.; Bowles, R. K. Some Fundamental Statistical Mechanical Relations Concerning Physical Clusters of Interest to Nucleation Theory. J. Chem. Phys. 1999, 111, 7501. (26) White, J. A. Lennard-Jones as a Model for Argon and Test of Extended Renormalization Group Calculations. J. Chem. Phys. 1999, 111, 9352. (27) Chandler, D. Introduction to Modern Statistical Mechanics; Oxford University Press: New York, 1987. (28) RuizMontero, M. J.; Frenkel, D.; Brey, J. J. Efficient Schemes to Compute Diffusive Barrier Crossing Rates. Mol. Phys. 1997, 90, 925. (29) Toxvaerd, S. Molecular-Dynamics Simulation of Homogeneous Nucleation in the Vapor Phase. J. Chem. Phys. 2001, 115, 8913. (30) ten Wolde, P.-R.; Ruiz-Montero, M. J.; Frenkel, D. Simulation of Homogeneous Crystal Nucleation Close to Coexistence. Faraday Discuss. 1996, 104, 93. (31) Hanggi, P.; Talkner, P.; Borkovec, M. Reaction-Rate Theory 50 Years after Kramers. Rev. Mod. Phys. 1990, 62, 251. (32) Knott, B.; Duff, N.; Doherty, M.; Peters, B. Estimating Diffusivity Along a Reaction Coordinate in the High Friction Limit: Insights on Pulse Times in Laser-Induced Nucleation. J. Chem. Phys. 2009, 131. (33) Auer, S.; Frenkel, D. Numerical Prediction of Absolute Crystallization Rates in Hard-Sphere Colloids. J. Chem. Phys. 2004, 120, 3015. (34) Ma, A.; Dinner, A. A Systematic Method for Identifying Reaction Coordinates in Complex Systems. Protein Sci. 2004, 13, 219. (35) Watanabe, H.; Suzuki, M.; Ito, N. Cumulative Distribution Functions Associated with Bubble-Nucleation Processes in Cavitation. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2010, 82, 051604. (36) Chkonia, G.; Wolk, J.; Strey, R.; Wedekind, J.; Reguera, D. Evaluating Nucleation Rates in Direct Simulations. J. Chem. Phys. 2009, 130, 064505. (37) Johnson, J. K.; Zollweg, J. A.; Gubbins, K. E. The Lennard-Jones Equation of State Revisited. Mol. Phys. 1993, 78, 591. (38) Blander, M.; Katz, J. L. Bubble Nucleation in Liquids. AlChE J. 1975, 21, 833. (39) Corti, D. S.; Kerr, K. J.; Torabi, K. On the Interfacial Thermodynamics of Nanoscale Droplets and Bubbles. J. Chem. Phys. 2011, 135, 024701. (40) Vincent, O.; Marmottant, P.; Quinto-Su, P. A.; Ohl, C.-D. Birth and Growth of Cavitation Bubbles within Water under Tension Confined in a Simple Synthetic Tree. Phys. Rev. Lett. 2012, 108, 184502. (41) Abascal, J. L. F.; Gonzalez, M. A.; Aragones, J. L.; Valeriani, C. Homogeneous Bubble Nucleation in Water at Negative Pressure: A Voronoi Polyhedra Analysis. J. Chem. Phys. 2013, 138, 084508. (42) Ma, A.; Dinner, A. Automatic Method for Identifying Reaction Coordinates in Complex Systems. J. Phys. Chem. B 2005, 109, 6769. (43) Kusaka, I.; Oxtoby, D. Identifying Physical Clusters in Bubble Nucleation. J. Chem. Phys. 1999, 111, 1104. (44) Torabi, K.; Corti, D. S. J. Phys. Chem. B 2013, DOI: 10.1021/ jp404149n.

ACKNOWLEDGMENTS This paper is based upon work supported by the National Science Foundation under Grant No. CHE-0718145.

■

REFERENCES

(1) Oxtoby, D. W.; Evans, R. Nonclassical Nucleation Theory for the Gas-Liquid Transition. J. Chem. Phys. 1988, 89, 7521. (2) Debenedetti, P. G.; Reiss, H. Reversible Work of Formation of an Embryo of a New Phase within a Uniform Macroscopic Mother Phase. J. Chem. Phys. 1998, 108, 5498. (3) Reiss, H.; Kegel, W. K.; Katz, J. L. Resolution of the Problems of Replacement Free Energy, 1/S, and Internal Consistency in Nucleation Theory by Consideration of the Length Scale for Mixing Entropy. Phys. Rev. Lett. 1997, 78, 4506. (4) Reiss, H.; Kegel, W. K.; Katz, J. L. Role of the Model Dependent Translational Volume Scale in the Classical Theory of Nucleation. J. Phys. Chem. A 1998, 102, 8548. (5) Reguera, D.; Reiss, H. The Role of Fluctuations in Both Density Functional and Field Theory of Nanosystems. J. Chem. Phys. 2004, 120, 2558. (6) Kusaka, I.; Oxtoby, D.; Wang, Z. On the Direct Evaluation of the Equilibrium Distribution of Clusters by Simulation. J. Chem. Phys. 1999, 111, 9958. (7) Kusaka, I.; Oxtoby, D.; Wang, Z. On the Direct Evaluation of the Equilibrium Distribution of Clusters by Simulation. II. J. Chem. Phys. 2001, 115, 6898. (8) Weakliem, C. L.; Reiss, H. The Factor 1/S in the ClassicalTheory of Nucleation. J. Phys. Chem. 1994, 98, 6408. (9) Reiss, H.; Kegel, W. K. Replacement Free Energy and the 1/S Factor in Nucleation Theory as a Consequence of Mixing Entropy. J. Phys. Chem. 1996, 100, 10428. (10) Reiss, H. Critique of the Determination of Equilibrium Cluster Distributions by Means of ″Direct Simulation″. J. Mol. Struct. 1999, 485, 465. (11) Reiss, H.; Bowles, R. K. Comparison between Two Methods for Mapping Fluctuations in a Simulation Cell onto a Macrovolume. J. Chem. Phys. 1999, 111, 9965. (12) Reiss, H.; Bowles, R. K. Mapping Volume Scale for Overlapping Clusters. J. Chem. Phys. 2000, 112, 1390. (13) Reiss, H.; Djikaev, Y.; Bowles, R. K. On a Debate over the Simulation and Mapping of Physical Clusters in Small Cells. J. Chem. Phys. 2002, 117, 557. (14) Ellerby, H. M.; Weakliem, C. L.; Reiss, H. Toward a Molecular Theory of Vapor-Phase Nucleation 0.1. Identification of the Average Embryo. J. Chem. Phys. 1991, 95, 9209. (15) Schaaf, P.; Senger, B.; Reiss, H. Defining Physical Clusters in Nucleation Theory from the N-Particle Distribution Function. J. Phys. Chem. B 1997, 101, 8740. (16) ten Wolde, P.; Frenkel, D. Computer Simulation Study of GasLiquid Nucleation in a Lennard-Jones System. J. Chem. Phys. 1998, 9901. (17) ten Wolde, P.; Ruiz-Montero, M.; Frenkel, D. Numerical Calculation of the Rate of Homogeneous Gas-Liquid Nucleation in a Lennard-Jones System. J. Chem. Phys. 1999, 1591. (18) Peters, B.; Trout, B. L. Obtaining Reaction Coordinates by Likelihood Maximization. J. Chem. Phys. 2006, 125, 054108. (19) Torabi, K.; Corti, D. S. Molecular Simulation Study of CavityGenerated Instabilities in the Superheated Lennard-Jones Liquid. J. Chem. Phys. 2010, 133, 134505. (20) Shen, V. K.; Debenedetti, P. G. A Computational Study of Homogeneous Liquid-Vapor Nucleation in the Lennard-Jones Fluid. J. Chem. Phys. 1999, 111, 3581. (21) Uline, M. J.; Torabi, K.; Corti, D. S. Homogeneous Nucleation and Growth in Simple Fluids. I. Fundamental Issues and Free Energy Surfaces of Bubble and Droplet Formation. J. Chem. Phys. 2010, 133, 174511. 12504

dx.doi.org/10.1021/jp404151h | J. Phys. Chem. B 2013, 117, 12491−12504