Toward a Molecular Theory of Homogeneous Bubble Nucleation: I

Article pubs.acs.org/JPCB

Toward a Molecular Theory of Homogeneous Bubble Nucleation: I. Equilibrium Embryo Definition Korosh Torabi and David S. Corti*,† †

School of Chemical Engineering, Purdue University, West Lafayette, Indiana 47907-2100, United States ABSTRACT: We propose an equilibrium embryo definition for homogeneous bubble nucleation within the pure component superheated Lennard-Jones liquid. The suggested embryo definition serves as an improvement to a previous (n,v) embryo model (Uline and Corti 2007). In that model, the constrained equilibrium between the bubble and the surrounding superheated liquid was maintained by placing n particles within a spherical volume v, without concern for the redundant counting of configurations and the relevance of the model to the dynamics of a nucleation process. To eliminate this redundancy, while only considering those embryos that should most likely appear at a transitional state, we now define the volume of the embryo via the use of a shell particle and only include n particles inside the volume that are deemed to be “vapor-like”. The underlying free energy surface of formation of the new (n,v) embryo model is generated via Monte Carlo simulation and also approximately by a suggested density functional theory method. The resulting surface implies that small and locally confined regions of near-zero density serve as precursors initiating homogeneous bubble nucleation. Furthermore, the nonredundant counting of the embryo configurations yields a well-defined and sharp conduit in the free energy surface that directs the initially formed embryos toward the critical nucleus. We discuss how the suggested equilibrium embryo model aids in both the identification of the transitional configurations and calculation of the average number density of the critical nuclei within the superheated liquid phase, which is the focus of the companion paper (DOI 10.1021/jp404151h).

I. INTRODUCTION Homogeneous bubble nucleation refers to the process by which a first-order phase transition occurs within a bulk metastable liquid phase in the absence of impurities or external perturbations. Due to the effectively random movements of the particles comprising the liquid, regions with densities both higher and lower than the average bulk density spontaneously form throughout the liquid. Since a metastable phase is not the lowest free energy state for the conditions imposed on the system, the possibility therefore arises that some of the spontaneously formed low density regions, or bubble embryos, initiate a phase transition. By surpassing some critical size, the bubble embryo then rapidly grows, eventually forming a macroscopic region of the new vapor phase.1 Since local spontaneous density fluctuations provide the trigger for the inevitable phase change within the metastable phase, a molecular-level study of homogeneous nucleation requires one to focus on a small but representative portion of the bulk metastable phase. The volume of this system must be large enough to accommodate a typical fluctuation that is relevant to homogeneous nucleation in a bulk fluid yet small enough (in terms of the total number of particles contained in the subsystem) to be computationally tractable with, say, a molecular dynamics (MD) simulation. For a liquid that is not too highly superheated, the dynamical trajectory that describes the time evolution of all the particles comprising such a system spends a large amount of time within a region of phase space that constitutes the basin of attraction of the metastable liquid phase. As long as the system remains within this region, one can © XXXX American Chemical Society

speak of an equilibrium metastable phase with well-defined thermodynamic properties. However, after a much greater time span, a trajectory of the system will eventually escape from the basin of attraction of the metastable liquid. As this trajectory exits from the basin, at least one bubble embryo will have grown beyond some critical size (the so-called critical bubble). This post-critical embryo initiates the irreversible phase transition, while the system can no longer be represented by an equilibrium state. The clear separation of the equilibration time scale of the metastable system (fluctuations that remain within the basin of attraction) and the nucleation time scale (average lifetime of the metastable system) is merely a necessary condition for the existence of the metastable phase. Such separation of time scales is typically met for a pure liquid that is not too superheated and in which homogeneous nucleation is the primary mechanism of the phase transition. However, by increasing the degree of superheating, i.e., moving the system toward the thermodynamic limit of stability or spinodal point, these time scales are not so easily separable. At near-spinodal conditions, the lifetime of the metastable fluid is too short to form an equilibrium phase; the phase transition occurs almost immediately via a process called spinodal nucleation.1 Unless one is dealing with highly superheated liquids, the separation of time scales makes the study of homogeneous Received: April 26, 2013 Revised: July 10, 2013

A

dx.doi.org/10.1021/jp404149n | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

An obvious example of an equilibrium-based approach is the well-known and still widely used phenomenological description of nucleation, called classical nucleation theory (CNT).1,11,12 Within CNT, bubble embryos, regardless of their molecular size, are described by a bulk uniform vapor phase that is separated from the surrounding uniform metastable liquid phase via a sharp interface to which a surface energy is assigned. In CNT, a bubble can be modeled as comprising n particles residing within a spherical volume v, with n and v being the chosen order parameter (in some cases,12 the internal pressure may replace n). The critical bubble is located at the saddle point of the (n,v) free energy surface, and its corresponding reversible work of formation denotes the free energy barrier of the nucleation process.12 The value of the free energy barrier is also required to determine the number density of critical bubbles throughout the metastable liquid, which along with a kinetic prefactor that describes how fast, on average, the critical bubble crosses over the barrier, allows one to calculate the rate of homogeneous bubble nucleation. CNT explains, at least qualitatively, the various competing factors that give rise to a free energy barrier, though its limitations have been known for some time.1 Yet, CNT does provide a convenient framework for understanding bubble nucleation. With this in mind, Uline and Corti developed a molecularlevel description of an (n,v) equilibrium embryo for both bubble and droplet nucleation.13,14 In their model, a rigid, diathermal, and impermeable boundary was employed as a constraint to separate an embryo with a given n and v from the surrounding metastable phase. Unlike CNT, this model accounted for the inhomogeneity that develops within both the embryo and the surrounding metastable fluid. The free energy surfaces and the related density profiles of the (n,v) equilibrium embryos were generated via Monte Carlo (MC) umbrella sampling simulations and separately via a novel constrained-density-functional theory (DFT) model for the pure component Lennard-Jones (LJ) fluid. With this embryo definition, a new and interesting picture of nucleation emerged. Namely, a locus of instabilities appeared along the free energy surfaces, beyond which the surrounding metastable liquid or vapor phase could not maintain its liquid-like or vapor-like densities. Each surface was also found to exhibit a flat ridge. Points on the ridge were found, however, to be representative of the same physical embryo, meaning the same configurations were being redundantly assigned to different (n,v) points. The DFT analysis was subsequently refined by Uline et al.,10 albeit in an ad hoc manner, in an attempt to eliminate this redundant counting of the configurations. While improvements were made, the current authors have since noted that additional changes to the chosen (n,v) equilibrium-based model were still needed. If the equilibrium embryo is going to faithfully describe the actual dynamics of a nucleation event, the embryo definition must ensure that the full configuration space of the metastable liquid phase be mapped uniquely (or nonredundantly) onto (n,v) order parameter space. How that mapping should occur is the focus of the current paper. As will be explained in more detail in the companion paper (DOI 10.1021/jp404151h), if such a mapping is not performed properly, the number density of transitional nuclei within the bulk metastable phase will not be correctly determined, leading to erroneous predictions of the nucleation rate. A further necessary modification to the present embryo definition, which prevents sampling of the configurations that are not within the

nucleation a computationally challenging problem, even for simple liquids. In order to observe a nucleation event, one needs to generate, for example, an extremely long MD simulation trajectory. To reduce the computational expense of waiting for the phase transferring trajectories to happen spontaneously, one can use the method of transition path sampling (TPS).2−4 TPS is a method by which an ensemble of transitional trajectories of a rare event is created in a systematic fashion, and has been used to generate nucleating trajectories within the Ising lattice model.5,6 Yet, implementing TPS for a more realistic system of Lennard-Jones particles is still a computationally expensive task. Another algorithm for systematically harnessing transposal trajectories of a rare event is the forward flux sampling (FFS).7 Wang et al. have implemented FFS to the bubble nucleation process within a LJ liquid.8 More detailed discussion of FFS and its application to a nucleation problem is presented in paper II (10.1021/jp404151h). What makes the development of a fully dynamical molecularlevel analysis of nucleation quite challenging is the exceedingly large number of degrees of freedom required to specify the local density fluctuations (i.e., coordinates of all the particles comprising the system) along a given trajectory in phase space. To bypass this difficulty, one typically views the local density fluctuations of the system through the lens of a more manageable number of collective variables that are assumed to be relevant to the process of nucleation. By employing these collective variables, also referred to as order parameters, the full degrees of freedom of the system are therefore mapped onto a significantly smaller number of variables. In this way, a point in order parameter space becomes representative of all the configurations of the system that meet the criteria dictated by the chosen order parameters. By working in this reduced space, certain key physical characteristics of the transitional nuclei become more readily apparent. Furthermore, the use of order parameters facilitates the computational procedures involved in a molecular simulation study of nucleation. Mapping the configurations of the system onto a chosen order parameter space can be accomplished by constraining the partition function to comply with the desired values of the order parameters while integrating over the remaining degrees of freedom available to the system. As such, the concept of an ensemble average or equilibrium embryo naturally arises. Moreover, by equilibrating or ensemble averaging the system under a set of internally imposed constraints dictated by the order parameters, a free energy, or minimum work of embryo formation, can then be defined. Within this framework, the nucleation rate therefore becomes the average rate at which trajectories that describe the evolution of average embryos cross over the free energy barrier. Such an equilibrium-based method invokes what has been referred to as “inversion of order of averaging”.9 Although each dynamical evolution of the system in the configuration space can be uniquely mapped (in the case of a properly defined embryo) onto a trajectory in the order parameter space, this coarse-grained trajectory may not necessarily describe how the system actually evolves during a nucleation event.10 Only if the actual system trajectory remains near the same point in the order parameter space for a long enough time, thereby sampling a sufficient number of configurations consistent with the specific set of order parameters, will the time evolution of the average embryo be the same as (or close enough to) the dynamics of the actual nucleation event. B



Article

Figure 1. Illustrations of how a given configuration can be assigned to multiple points in the (n,v) order parameter space. (a) The dashed lines represent different (n,v) embryos that are assigned to the same snapshot. (b) Even after distinguishing the vapor-like particles (open circles) from the liquid-like particles (dark circles), there still exits a degree of freedom in choosing the boundary of the embryo, as shown by the dotted circles. To remove this degree of freedom, the boundary of the embryo is always taken to coincide with the shell particle (red circle), which is defined to be the closest liquid-like particle to the center of the embryo.

II. (n,v) EQUILIBRIUM EMBRYO DEFINITION REVISITED The (n,v) equilibrium embryo definition originally introduced by Uline and Corti13 is the ensemble average of all the configurations of the system that happen to have n particles inside a spherical volume v located at a fixed position within the bulk metastable liquid phase. With this definition, any configuration of the n particles, regardless of the vapor-like character of those particles, is assigned to a certain (n,v) embryo. Such an “all or nothing” definition leads to a redundant counting of embryo states, such that any configuration could possibly be assigned to different points in the (n,v) order parameter space (see Figure 1). As a result, a flat ridge appeared in the free energy surface. These ridge points, being the maxima of v for each n, were found to be described by the same continuous density profile, and so had the same free energy. Thus, the ridge points are not distinct embryos, and are simply generated by altering the location of v along the given density profile and determining the corresponding value of n. This problem does not appear only along the ridge. Other large enough (n,v) embryos contain significant portions of liquid-like densities, and also should not be counted as distinct fluctuations (see ref 10 for further details). To eliminate this problem of redundant counting of configurations, the (n,v) embryo definition requires an additional constraint: the embryo volume must contain only “vapor-like” particles. This constraint also serves an additional purpose. Bubbles with a large number of “liquid-like” particles are not representative of those low density regions assumed to be relevant to the process of bubble nucleation, and should not be included in the free energy surface. The criterion used to distinguish between vapor-like and liquid-like particles is arbitrary, though not completely so. ten Wolde and Frenkel18 showed via molecular simulation for a pure component LJ fluid that the particles within a bulk liquid phase and a bulk vapor phase at coexistence have notably different numbers of particles neighboring them. Particles are defined to be neighbors if they are located within a distance of 1.5σ from one another (where σ is the LJ particle diameter). Inside the liquid phase, particles are hardly ever found to have less than five neighboring particles, with the average number of neighbors being more than ten. Particles inside the vapor phase almost always have less than five neighboring particles. Here, as was also adopted in ref 8, we utilize the same criterion. For a given configuration, we label all

boundaries of the metastable liquid phase, is presented in the companion paper (DOI 10.1021/jp404151h). The new equilibrium embryo definition presented here uniquely defines the volume of the embryo via a “shell particle”9,15−17 and also ensures that the n particles residing in the volume v are all “vapor-like”.8,18 Since this modified embryo definition requires that the location of all particles be known, we turn to MC simulation to generate the (n,v) free energy surface. In the companion paper (DOI 10.1021/jp404151h), we show how this free energy surface can be used to properly calculate the number density of the transitional nuclei within the superheated liquid phase. This in turn, when supplemented with the obtained kinetic data of the time evolution of the transitional nuclei into macroscopic bubbles, leads to a prediction of the rate of homogeneous bubble nucleation. Since molecular simulation is now needed to apply the new embryo definition, the generation of the free energy surface at any state point is time-consuming. Although the new bubble definition cannot be rigorously incorporated into DFT, we nevertheless discuss how the shell particle can be accounted for in DFT calculations, invoking similar approximations already present in the standard DFT analysis. We also show how a free energy penalty that mimics the constraint needed to ensure the vapor-like trait of the particles contained within the bubble can be easily incorporated into DFT. This modified DFT method yields qualitatively correct trends, as compared to the simulation results. As such, DFT still has the potential to reliably, and more rapidly, explore bubble nucleation over a broad range of conditions. The outline of the present paper is as follows. In Section II, the two needed additional constraints to the (n,v) embryo definition of Uline and Corti13 are presented. In Section III, we discuss the molecular simulation procedures needed to generate the free energy surface of bubble formation within the superheated LJ liquid. A resulting free energy surface along with additional comment upon the nucleation mechanism implied by the new embryo definition will also be provided. The approximate incorporation of the new additional constraints into DFT is presented in Section IV. Conclusions are provided in Section V. How the new embryo definition and its corresponding free energy surface can be used to calculate the rate of homogeneous bubble nucleation is the focus of the companion publication. C



Article

Figure 2. The reversible work of forming a bubble with a given n and radius λ within the superheated LJ liquid (see text for bulk conditions) generated from the probabilities calculated within the uniform phase simulation. The dashed line is the product of the bulk density and the volume of the shell volume, or 4πλ2Δλ ρbulk.

the particles that have five or fewer neighbors as vapor-like. All other particles are labeled as liquid-like. Hence, when determining the equilibrium properties of the bubble, we now average over only those configurations in which n vapor-like particles are located inside the embryo volume. While the above constraint removes those (n,v) bubbles that contain portions of the liquid, another potential redundancy still remains when mapping a given configuration to a particular volume in the order parameter space. As illustrated in Figure 1b, the volume v can be assigned to any value between the closest liquid-like particle and the farthest vapor-like particle from the origin of the embryo. In order to remove this volume redundancy, the volume of the sphere is chosen to always coincide with the liquid-like particle closest to the center of the embryo. The volume is therefore uniquely defined via a liquidlike “shell-particle”, which can reside anywhere within the spherical shell of infinitesimal thickness around the defining volume of the embryo. A shell particle was also utilized to eliminate the redundant counting of volume states in previous studies of homogeneous droplet nucleation within a supercooled vapor.9

bubbles forming within the pure component superheated LJ liquid. The particles interact through the following truncated and shifted LJ potential ⎧ ⎡⎛ ⎞−12 ⎛ ⎞−6 ⎤ ⎪ 4ε⎢⎜ r ⎟ − ⎜ r ⎟ ⎥ − us r /σ ≤ 4 ⎝σ ⎠ ⎦ uLJ (r ) = ⎨ ⎣⎝ σ ⎠ ⎪ ⎪ ⎩ 0 r /σ > 4

us = 4ε(4−12 − 4−6)

where ε is the LJ well-depth. For all the following MC simulations, the system is maintained at a reduced temperature of kBT/ε = 0.8 and the same reduced pressure of Pσ3/ε = −0.33 (the reduced spinodal and binodal pressures at this temperature are −0.60 and 0.006, respectively19). The reduced bulk density of the liquid at these conditions is ρbulkσ3 = 0.745. In the rest of the paper, all parameters are reported in dimensionless form, unless specifically mentioned, using the scales of length, energy, temperature, and pressure of σ, ε, ε/kB, and ε/σ3, respectively. III.A. Uniform Liquid Phase Local Density Fluctuations. For an (n,v) embryo that does not have a large free energy of formation, its frequency of occurrence within the metastable liquid can be calculated with high enough accuracy simply by simulating a uniform phase and periodically examining random snapshots of the equilibrated system. From eq 1, values of W around 13−14 kBT or less correspond to probabilities greater than about 10−6, which can be evaluated within a standard bulk simulation of the metastable liquid. To generate the required snapshots of the equilibrated system, we first ran an isothermal−isobaric MC simulation for 40 000 equilibration cycles with 1000 LJ particles confined within a cubical volume with periodic boundary conditions. Each MC cycle consisted of 1000 attempted particle moves, with particles chosen at random, and one attempted volume move. The sizes of the maximum allowed particle and volume moves were adjusted for a 50% acceptance rate only during the equilibration portion of the simulation. After the equilibration cycles, we ran 60 000 sampling cycles. At the end of each cycle, the current snapshot of the system was used to identify the (n,v) embryos. Since we are simulating a uniform system, the center of the embryo may be located at any arbitrary point within the simulation box. To improve the statistical significance of the results, we randomly chose 50 different locations within each selected snapshot to serve as the

III. (n,v) FREE ENERGY SURFACE VIA MOLECULAR SIMULATION With our modified (n,v) embryo definition now requiring explicit knowledge of all particle positions in order to identify both vapor-like and liquid-like particles and the location of the shell particle, we turn to molecular simulation to determine the free energy surface of bubble formation. The free energy, or the reversible work of formation W, is defined as W (n , λ) = −kBT ln[P(n , λ)]

(2)

(1)

where kB is Boltzmann’s constant, T is the absolute temperature, and P(n, λ) is the probability of finding, at a given point within the metastable phase, a configuration with n particles inside a volume of radius λ that satisfies the following two conditions: (1) the closest liquid-like particle is located at a distance λ from the center of the embryo and (2) the spherical boundary of radius λ contains only n vapor-like particles. P(n, λ) is the same for any arbitrary point inside the bulk metastable liquid phase, chosen as the center of the embryo. While the suggested embryo definition and the molecular simulation procedures to generate the related free energy surface are applicable to several model liquids, we consider here D



Article

the free energy surface of different homogeneous nucleation models.13,14,18,21,22 In the present study, we employ an umbrella sampling procedure that only implicitly utilizes a biasing potential, merely limiting the accessible configurations of the embryo within some prescribed window of allowed volume states (similar to what was done for the Ising model20). To illustrate the method, we first discuss how to determine the free energy difference between two cavities, or embryos of (n = 0, λ1) and (n = 0, λ2). A void impermeable spherical boundary of fixed radius λ1 is placed at the center of the simulation box. The superheated liquid surrounding the λ1 cavity is then allowed to come to equilibrium within an isothermal−isobaric MC simulation. All particle moves that lead to a particle crossing the boundary of radius λ1 are rejected. After 20 000 equilibration cycles are completed, and within random snapshots of the system chosen after every 10 MC cycles, the closest liquid-like particle to the center of the simulation box is located. A histogram of different cavity radii are generated, allowing us to calculate the relative probability of the two embryos of (n = 0, λ1) and (n = 0, λ2) for λ2 > λ1. Our method makes use of a biasing potential that is infinite for any configuration that has an unwanted particle inside λ1 and zero for all other configurations that are being sample during the simulation. Hence, the calculated relative probabilities do not need to be corrected for the bias imposed, though we need to know the absolute free energy for a cavity of for at least on value of λ. We also do not need to strictly impose an upper bound for each sampling window, since larger embryos are naturally less probable. To get the best possible statistics, the separation between λ1 and λ2 should not be too large. If this separation is too small, however, reasonable statistics with fewer samplings are obtained, but more sampling windows are needed to climb up the free energy surface. For the system at the particular state of our study, we found that a radial interval of 0.1 yielded reasonably accurate frequency ratios with not more than 100 000 MC sampling cycles. When dealing with embryos of n > 0, the umbrella sampling procedure explained above, needs to be slightly modified. First, the MC simulation should only sample configurations with vapor-like particles inside the boundary of the embryo. To achieve this, we begin the MC simulation from a random configuration in which, there are n vapor-like particles inside the impermeable spherical boundary of desired radius. Then the MC simulation proceeds as before, but now all particle moves that cause any of the n particles to have more than 5 neighbors are rejected and an additional trial boundary move is performed. Maintaining a fixed impermeable sphere of radius λ1, as was done for the n = 0 embryos, improperly prevents the n vapor-like particles inside the embryo from occupying the volume between λ1 and λ2. In order to allow these particles to sample the entire volume allocated to them, we therefore allow the radius of the impermeable boundary to change during the simulation. At the end of each MC cycle, a switching attempt in the radius of the impermeable boundary, from λ1 to λ2 or vice versa depending on the initial state of the boundary, is performed and the trial move is accepted only if no particles crosses the surface of the sphere (thereby keeping n fixed). In this way, the relative probability of two embryos with different radii but the same n can be determined as the ratio of the number of configurations with the boundary at λ1 and a shell particle within λ1 + Δλ to those with the boundary at λ2 and a shell particle within λ2 + Δλ.

geometric center, or origin, of the embryo. About each origin, we first determined the radial distance to the closest liquid-like particle. Next, we counted the number of particles that fell inside the hypothetical sphere constructed around the center with a radius taken to be the distance to the first liquid-like particle. (The closest liquid-like particle is the shell particle that defines the system volume and is not included in the count of the n vapor-like particles that reside within the sphere.) In this way, all the particles that fall inside the sphere are, by construction, vapor-like. Despite being a continuous variable, the radius was nevertheless discretized into bins of finite thicknesses. The first bin was taken to be a complete sphere of radius 0.04, while the rest were spherical shells of thickness Δλ = 0.04 (built around the first sphere). In this way, we determined the probability of having n vapor-like particles within a spherical volume v = 4πλ3/3 centered about a given point with a liquid-like shell particle that falls inside a spherical shell of thickness Δλ about the sphere. The results of this direct sampling method are presented in Figure 2, which plots W(n, λ) in units of kBT versus radius λ for embryos with n = 0, 1, and 2. Unlike what was obtained from previous (n,v) embryo models that did not include the shell particle and vapor-like constraints,13 the minima no longer appear at W = 0. This is, however, not unexpected, as W = 0 implies that the probability of occurrence is unity, which should not be the case with this new embryo definition. For example, consider n = 0, or a cavity. The likelihood of successfully inserting the cavity such that no particles reside inside and a liquid-like particle is located within the finite spherical shell will increase as the radius decreases, but will not reach unity. For small enough volumes, the probability of inserting a cavity and finding a particle inside the spherical shell is approximately equal to 4πλ2 Δλρbulk (which is the dashed line in Figure 2; for a continuous radius, W diverges at λ = 0 but does not in Figure 2 since we determine probabilities for nonzero volumes about a zero radius). In addition, the value of W at a minimum increases with an increase in n. At least for small embryos, and for volumes around these minima, we expect that another vapor-like particle will not be so easily accommodated, leading to an increase in the free energy as n increases. Interestingly, the increase in the free energy between n = 1 and n = 2 is not as large as the increase between n = 0 and n = 1. Due to their large free energies, the above uniform phase random sampling procedure becomes statistically inaccurate, or not possible, for embryos with λ > 1.5 or n > 2. Much longer runs are needed to determine the frequencies of occurrence for these embryos, so much so that the method becomes computationally inefficient. To determine W(n,λ) for larger embryos, an umbrella sampling method particularly suited to studying low probability fluctuations is needed. III.B. Umbrella Sampling MC Simulations. To determine values of W beyond ∼14 kBT, an umbrella sampling or biasing simulation method20 must be employed. This is usually done by adding a biasing potential to the potential energy in order to force the system to sample more frequently a particular region of phase space. To sample relatively unlikely states, relative frequencies of occurrence are typically determined within a restricted range, or window, of the chosen order parameter. By correcting for the effect of the biasing potential in each window, along with knowledge of the probabilities for at least one unbiased window, the unbiased free energy over a large range of order parameter values can be constructed. Different variations of the umbrella sampling technique have been used to generate E



Article

In the above, N − n is the number of particles in the surrounding liquid phase that are contained in the volume (V − v) external to the embryo (though we divide by (N − n − 1)! because the shell particle defining the volume v is distinguishable from the other particles; by definition, the shell particle is located at a point within the spherical shell of volume 4πλ2Δλ), V is the total volume of the entire system, r(N) denotes the position vectors of all the N particles of the system, U(r(N)) is the potential energy of interaction among all the particles, d3nr is, for example, the volume elements of the positions of the n particles, P is the imposed external pressure, and ξ(r(N)) is a constraint function that has a value of unity for configurations in which all the n particles inside the spherical volume v are vapor-like, and zero otherwise. Δn+1 is defined similarly to Δn in eq 5, though with an additional vapor-like particles inside v while the number of particles within the surrounding volume V − v is still equal to N − n. As such, the ratio of these two partition functions can be expressed as

One important feature of the above simulation procedure is that for any given value of n there is a limiting value of λ beyond which the isothermal−isobaric MC simulation becomes unstable. This instability is of the same nature as what has been previously reported.10,13,23 By increasing the size of the impermeable spherical boundary past a critical size, the metastable phase surrounding the sphere no longer maintains a liquid-like density profile. Consequently, within an isothermal−isobaric simulation, the size of the simulation box diverges and the overall density of the system rapidly drops. This instability is an indication that the simulation is sampling configurations that are not members of the equilibrium metastable phase. A more in depth discussion of this issue is provided in the companion paper, where we demonstrate how the application of an additional constraint properly limits a simulation to the configuration space of the metastable phase. III.C. Particle Insertion Simulation. With the umbrella sampling procedure explained in the previous section, we are able to calculate the Gibbs free energy difference between embryos with different values of v but the same value of n. However, we also need to determine the free energy difference between embryos with different n for at least one fixed value of v. Since our free energy surface is generated for a closed system, W(n + 1, v) − W(n,v) is equal to the sum of the reversible work required to insert an additional vapor-like particle inside an (n,v) embryo (keeping the number of particles outside of the embryo fixed) plus the reversible work required to remove a particle from the surrounding liquid phase. The latter term is simply the negative of the chemical potential of the surrounding fluid, which is also the chemical potential of the uniform metastable liquid, μbulk, at the imposed temperature and pressure (as chemical equilibrium is maintained throughout the surrounding fluid). Using the well-known Widom insertion method,24 μbulk can be calculated within a separate simulation of the bulk fluid phase via the following equation ⎡ 1 ⎤ ⎥ exp( U ) μ bulk = −kBT ln⎢ ⟨ − β ⟩ ins bulk ⎢⎣ ρbulk Λ3 ⎥⎦

Δn + 1 1 1 [ = (n + 1) Λ3 Δn

∫ dV ∫v d3nr ∫v d3rins

∫V −v d3(N−n−1)r exp[−β(U(r(N); rins) + PV )] ∫ dV ∫v d3nr ∫V −v d3(N−n−1)r

ξ(r(N ); rins)]/[

exp[−β(U (r(N )) + PV )]ξ(r(N ))]

where rins is the coordinate of the additional particle. Due to the pairwise additivity of the LJ potential, U(r(N);rins) can be written as the sum of U(r(N)) and the potential energy of the interaction of the inserted particle with the rest of the particles in the system, Uins(r(N);rins). Meanwhile, ξ(r(N);rins), which is the constraint function of the system with the additional particle, can be expressed as the product of ξ(r(N)) and an additional constraint function ξins(r(N);rins) that takes a value of unity if after the insertion of the additional particle inside the volume v all the n + 1 particles are vapor-like, and zero otherwise. Therefore, eq 6 can be simplified to

(3)

Δn + 1 1 1 = ⟨ d3rins exp[−β(Uins(r(N ); rins))] (n + 1) Λ3 v Δn ξins(r(N ); rins)⟩(n , v)

∫

where Λ is the de Broglie wavelength, β = 1/kBT, Uins is the change in the potential energy of the system due to the insertion of the test particle and ⟨...⟩bulk indicates the ensemble average over random insertions within the bulk phase. The reversible work required to insert an additional vaporlike particle inside an (n,v) embryo spherical boundary can also be obtained via a molecular simulation procedure closely similar to the standard Widom insertion method. To see how, note that the reversible work of inserting another particle into a given embryo is given by ⎛Δ ⎞ ΔGins = −kBT log⎜ n + 1 ⎟ ⎝ Δn ⎠

=

1 1 n! (N − n − 1)! Λ3N

1 v ⟨exp( −βUins)ξins⟩(n , v) (n + 1) Λ3 (7)

in which ⟨ ⟩(n,v) represents the following volume-average of the inserted particle over the volume v of the embryo ⟨exp(− βUins)ξins⟩(n , v)

∫v d3rins exp[−β(Uins(r(N); rins))]ξins(r(N); rins)⟩(n,v)]/v

≡ [⟨

(8)

(4)

To calculate the ensemble average in eq 8, we first generated an equilibrated (n,v) embryo via an isothermal−isobaric MC simulation that constrained the system to sample configurations with n vapor-like particles inside the spherical volume v centered at the origin of the simulation cell and a liquid-like shell particle at the surface of the sphere. The simulation was continued during which a large number of uncorrelated configurations of the system were generated. For each configuration, the test particle was inserted at several random

where Δn is the isothermal−isobaric partition function of the system containing an (n,v) embryo that can be written as Δn = 4πλ 2Δλ

(6)

∫ dV ∫v d3nr

∫V −v d3(N−n−1)r exp[−β(U(r(N)) + PV )]ξ(r(N)) (5) F



Article

Table 1. Free Energy Differences between (n, λ) and (n + 1, λ) Embryos Calculated from the Particle Insertion Method Described in section III.C

locations inside the spherical volume v. To expedite the computations, Uins(r;rins) was only calculated when the neighbor constraint was not violated (in the case of constraint violation, the term inside the bracket on the left-hand side of eq 8 is zero). Once this ensemble average was determined, the free energy difference between an (n + 1, v) embryo and an (n,v) embryo was evaluated using

λ = 2.0 n to n+ 1 0 1 2 3 4 5 6 7 8 9

⎛Δ ⎞ G(n + 1, v) − G(n , v) = −kBT ln⎜ n + 1 ⎟ − μ bulk ⎝ Δn ⎠ ⎡ vρ ⟨exp( −βUins)ξins⟩(n , v) ⎤ bulk ⎥ = −kBT ln⎢ × ⟨exp( −βUins)⟩bulk ⎦ ⎣ (n + 1)

(9)

III.D. Simulation Results. For all simulations, the size of the system (total number of particles within the simulation box) was chosen large enough to ensure that a flat tail for the density profile (equal to the metastable liquid bulk density ρbulk = 0.745) was reached within at least a distance of 2.0 from the walls of the simulation box. Since the LJ potential is truncated at a radius of 4.0 (see eq 2), this ensures that the embryo sees around itself (because of periodic boundary conditions) a uniform liquid phase and not periodic image of itself. Of course, a very large system size would automatically meet this condition, though at the cost of computational time. Through a few trial and error, we found that the following system sizes do not introduce any undesired finite size effects: 2000, 3000, 4000, and 5000 particles for radii domains less than 2.0, 2.0 to 2.5, 2.5 to 3.0, and greater than 3.0, respectively. As explained in section III.B, for a given value of n, the free energy curve could be extended all the way to the radius where the simulation becomes unstable. For the conditions of the present study, the limits of instability for n = 0, 1, 2, 5, 8, 10, 12, and 15 were found to be equal to λ = 2.6, 2.6, 2.7, 2.9, 3.1, 3.2, 3.2, and 3.3, respectively. After generating the free energy curves for the values of n mentioned above, we performed test particle insertion simulations (discussed in section III.C) in order to place each free energy curve on its proper relative position. To ensure better statistics, we carried out each insertion simulation at a radius in which the particle insertions were not associated with large potential energy differences (see eq 8). For the range of n = 1 to 9, we used λ = 2.0, and for n = 10 to 14, we used λ = 2.5. For the length of the simulations, 20 000 equilibration cycles and 200 000 sampling cycles with 100 test particle insertions for each snapshot saved at the end of each cycle seem to be long enough to ensure an accurate averaging. This is evaluated by calculating a running average and making sure the value levels off before the sampling is terminated. Free energy differences between (n,v) and (n + 1, v) embryos are listed in Table 1. The resulting free energy surface of the ensemble average (n,v) embryos for a broad range of order parameter space is shown in Figure 3. In Figure 3, we have also included the uniform phase simulation free energies (previously shown in Figure 2). Agreement between the two sets of data for n ≤ 2 serves as validation of the accuracy of the developed simulation procedures. The free energy surface presented in Figure 3, while retaining similar features found using the earlier embryo model of Uline and Corti,13 yields a picture of homogeneous bubble nucleation that is nevertheless different from previous work in certain key aspects. As explained in section III, the problem of the redundant counting of configurations has been completely

to to to to to to to to to to

1 2 3 4 5 6 7 8 9 10

λ = 2.5 ΔG/kBT −1.70 −1.09 −0.71 −0.39 −0.08 0.23 0.55 0.91 1.31 1.75

n to n + 1 10 11 12 13 14

to to to to to

11 12 13 14 15

ΔG/kBT 1.24 2.56 3.94 5.40 6.94

eliminated from the model. Each (n,v) point along the free energy surface therefore represents an ensemble average of those configurations that exclusively belong to that (n,v) point in the order parameter space. As a result, the new free energy surface no longer has a flat ridge, which follows from different (n,v) points being represented by the same configurations. In addition, Figure 3 shows that for a constant value of the radius the free energy exhibits a relatively sharp minimum in the ndirection. This minimum in n for any fixed value of λ corresponds to the most probable number of vapor-like particles that can be found inside the spherical volume of radius λ (that contains no liquid-like particles). Another important feature of the new free energy surface is that embryos are no longer found to contain regions of liquidlike density. This is, of course, an immediate consequence of the vapor-like constraint imposed on the embryo particles. With such constraint, for a given value of n, the volume of the embryo needs to be large enough to allow the particles to maintain their vapor-like characters. By shrinking the radius below some certain value, the configurations that comply with the vapor-like constraint become so scarce that the free energy surface rapidly increases to larger values (for clarity, this diverging behavior is not explicitly shown in Figure 3). In contrast, the (n,v) free energy surface reported in ref 13 did not include the vapor-like constraint; therefore, for any fixed value of n an embryo could be shrunk to volumes that corresponded to bulk liquid and even higher densities. Within such a surface, a valley of zero free energies appeared along which were found all the (n,v) points with a flat density profile equal to the bulk liquid density. All of these embryos were, however, representative of the same physical system, as the interior and exterior density profiles were nothing more than the same uniform and continuous density profile consistent with the bulk metastable liquid phase. Looking at such surface, one infers that nucleation is most likely initiated from an arbitrary point at the bottom of the flat valley of the uniform liquid phase, with the embryo proceeding to climb up the free energy surface by expanding its volume in order to reach the flat ridge. Since all points along the flat ridge also correspond to the same density profile, the net effect of the process is simply the appearance of a critical vapor-like density fluctuation that eventually surmounts a free energy barrier. Although this mechanism does conform to our expectations, we note that nearly all trajectories moving along this free surface start and end at the same fluctuation while sampling different embryos in between. G



Article

Figure 3. Free energy surface of formation of (n,v) embryos within the pure component LJ superheated liquid at reduced temperature of kBT/ε = 0.8 and reduced pressure of Pσ3/ε = −0.33. The closed symbols are the umbrella sampling results and the open symbols are the corresponding values calculated via the uniform phase local density fluctuation simulation method.

basin of attraction. Once that fluctuation appears, and the phase transition begins to occur (moving the liquid further away from the basin of attraction), one cannot speak of an equilibrium metastable liquid phase. Thus, only precritical embryos should appear within a truly equilibrium metastable liquid. Once a postcritical embryo has developed, which is an indication that the phase transition has most likely begun, one should not label this inherently nonequilibrium system as an equilibrium phase. In this regard, the appearance of stability limits prior to the development of any maxima is consistent with the elimination of those fluctuations from an equilibrium free energy surface that would cause the system, more often than not, to move away from its basin of attraction. If we are dealing with a metastable liquid at equilibrium (i.e., no indications of a phase transition taking place), the probability density distribution of bubble sizes should therefore decrease with an increase in volume, and not exhibit a minimum (as the appearance of a minimum indicates that postcritical embryos have developed and are beginning to grow into the new phase). Such a conclusion is consistent with what has been reported by Toxvaerd,25 albeit for droplet nucleation within a metastable vapor phase. As long as the configuration space was limited to the equilibrium metastable phase (i.e., the uniform metastable phase prior to the phase transition occurring), the liquid cluster probability distribution obtained during the molecular simulation was found to be a monotonously decreasing function of the particle number. No clusters beyond a certain maximum particle number were found to appear, and the cluster distribution did not exhibit a minimum.

However, with our new embryo definition, by removing the redundant counting of configurations, the free energy surface now includes only those (n,v) points that are compact regions of vapor-like density. Embryo trajectories no longer begin and end at exactly the same density fluctuation. A well-defined channel, starting from small values of n and small volumes, serves to funnel a range of similar (though not identical) trajectories up toward the top of the free energy barrier. Specifically, the free energy surface of Figure 3 now suggests that bubble nucleation begins with the formation of tiny cavities, regions devoid of particles, as these are the only embryos (λ < 1.0 for n = 0) with free energies of formation very near to zero. Starting from cavities is perhaps not surprising, as small cavities frequently appear within the metastable liquid phase, but we now see that a cavity must expand beyond λ = 1.0 before it becomes favorable for an embryo to include at least one vapor-like particle. Only beyond these sizes does a bubble, i.e., an embryo containing a few vapor-like particles, have a lower free energy of formation than a cavity. Embryos should therefore begin as small cavities and transition into bubbles once they reach a certain size. As such, and as expected, those trajectories that are most likely to reach the free energy barrier will be representative of growing bubbles. In this way, the new (n,v) free energy surface of Figure 3 can be interpreted as a representation of the conduit that channels small cavities all the way up to the larger transitional bubbles. Finally, we also see that the new free energy surface no longer exhibits maxima in volume for a given particle number. Stability limits, which are signals that the (n,v) embryos have become so large that the system has begun to sample those configurations that are not members of the equilibrium metastable liquid phase, now appear before any maxima are seen (and which mostly appeared beyond the maxima in the previous surface). As will be elaborated in more detail in the companion paper, the configuration space of any metastable phase has an inherent limit within which an equilibrium system can be defined. However, here we note that an embryo that surmounts the free energy barrier must be considered as a fluctuation that will most likely instigate the phase transition, i.e., the system has been brought outside of the metastable

IV. CONSTRAINED-DFT FREE ENERGY SURFACE OF (n,v) EMBRYOS Calculating the free energy surface via molecular simulation is a computationally intensive task. Free energy surfaces can, however, be generated more rapidly with density functional theory (DFT), allowing for the possibility of exploring homogeneous nucleation over a much broader range of metastable fluid conditions. Unfortunately, the modified (n,v) embryo definition cannot be rigorously incorporated into DFT. H



Article

on the configuration integral of an (n,v) embryo. As the average density inside the bubble increases, the vapor-like constraint further restricts the configurations that may be sampled by the embryo. As the number of accessible configurations decreases, the probability of observing the corresponding average embryo decreases as well, or the free energy of its formation increases. In our molecular-based embryo definition, only vapor-like particles are allowed to reside within the embryo, i.e., those particles that have no more than five neighbors within a distance of 1.5. Specific configurations of particles are not directly considered with DFT; however, as DFT deals with density profiles and not individual particle locations. The density profile does provide information, in an average sense, on the probability of finding a particle at a given position. Within the random phase approximation29 employed by eq 10, some information on the local density of neighboring particles about a given particle at the origin is also provided. To determine whether the local density profile most likely corresponds to a vapor-like particle at a particular location, we first integrate the density profile within a spherical volume of radius 1.5 centered about the given point of interest. If the result is a number less than six, which approximately corresponds to a particle sitting at the point of interest having, on average, less than five surrounding nearest neighbors, that particle is deemed to be vapor-like and no free energy penalty is imposed. If the result is greater than six, a penalty is added to the free energy, as the particle is deemed to be liquid-like. After averaging the result over all positions just inside the embryo, we suggest a free energy penalty term of the following form

Nevertheless, key aspects of the embryo definition can be approximately accounted for within DFT. In this section we present how such modifications to the constrained-DFT formulation of Uline and Corti13 can be applied. Constrained-DFT begins with the following form of the grand potential functional, Ω,10 Ω[ρ] =

∫ d3rfhs [ρ] + 12 ∫ d3r ∫ ρ(r′) − μ0

∫ d3rρ(r)

d3r′φ2(|r − r′|)ρ(r) (10)

subject to the following constraint

∫ d3rρ(r)θ(λ − |r|) = n

(11)

where f hs is the free energy density of the reference hard sphere fluid calculated using a later modification of the fundamental measure theory of Rosenfeld,26,27 ρ(r) is the density profile that spans both the embryo and the surrounding fluid and is, in general, discontinuous at the radius λ that determines the spherical boundary used in the definition of the (n,v) embryo, φ2 is the attractive part of the LJ potential as designated by Weeks et al.,28 and θ is the Heaviside step function, which ensures that only n particles are located inside the embryo volume. The above DFT formulation of eq 10 and eq 11 does not account for the additional constraints of our new (n,v) embryo definition. In DFT, molecular-level information about different configurations of the system is inherently washed out by being lumped into one single ensemble average density profile. Hence, the shell-particle and vapor-like constraints cannot be rigorously incorporated into DFT. In a previous publication,10 the vapor-like constraint was accounted for in an approximate, and ultimately ad hoc, manner by simply terminating the free energy surface at those (n,v) points where the average density inside the defining volume of the embryo became larger than ρσ3 = 0.296. This value corresponds to the average distance between particles being less than 1.5σ (and so the particles were deemed to be liquid-like), which is equivalent to (v/n)1/3 < 1.5σ. Although this procedure eliminated some embryos that are probably not relevant to the process of bubble nucleation, the free energy surface was, however, abruptly truncated in some regions before local minima were expected to appear. As such, no divergences in the free energy were (or even should have been) found at those limiting volumes for which the constraint was first violated. In order to generate a smoothly varying free energy surface away from the locus of stability and near those volumes where a given constraint should be violated, while still accounting for the vapor-like constraint to some degree in the DFT calculations, we introduce a free energy penalty to those embryos that violate or come close to violating the vapor-like criterion. In this way, we eliminate somewhat more naturally those nonrelevant density fluctuations by making them less probable, or by increasing their free energies of formation. This penalty is again applied after the density profiles are obtained from eqs 10 and 11, in which the vapor-like constraint is not considered, but as long as the resulting free energy surface is similar to the “exact” surface, the predicted probabilities of finding various (n,v) embryos, or what are the most likely trajectories along the free energy surface, will be sufficiently close to what is found in simulation. The free energy penalty term should be based on the effect of the vapor-like constraint

free energy penalty ⎡ ∫ d3rθ(6 − ∫ d3r′ρ(|r − r′|)) ⎤ |r − r |< 1.5σ ⎥ ′ = −αkBT ln⎢ v ⎢ ⎥ v ⎣ ⎦ (12)

where θ is a step function and v is the volume of the embryo. Since the amount by which the above penalty changes the free energy of the system is unknown, eq 12 includes the arbitrary constant α. Although the form of α is not known, the following argument indicates that α should be proportional to n. If N ideal gas particles initially confined to a volume V were then restricted to the smaller volume εV, the configuration integral would change from VN to (εV)N, thereby increasing the free energy by = −kBTN ln ε. While eq 12 does not correspond to a reduction in the volume of the system, it does account for a reduction in the configuration integral that should be shared equally by each particle (as each particle has access to the entire volume of the embryo and is subject to the same constraint). In what follows, we set α = n, which as will be seen below appears to be an appropriate choice. Once the free energies and corresponding density profiles for different (n,v) embryos are generated by the constraint-DFT method of eqs 10 and 11, the free energy surface is modified by adding the free energy penalty of eq 12. Besides that, we also need to account for the shell-particle constraint. While the shell particle can be added directly into eq 10 (corresponding to an additional external field in the system), solving for the corresponding density profile becomes somewhat prohibitive as the resulting density profile is no longer spherically symmetric. The net effect of the shell particle, or how it alters the free energy of the system, can, however, be determined by comparing the configuration integral of the embryo with a shell I



Article

Figure 4. Free energy surface of forming (n,v) embryos within the superheated LJ liquid using the modified form of the constraint-DFT method. The free energy change due to the addition of the shell particle is determined exactly while that following from the application of the vapor-like/ liquid-like constraint is determined approximately using eq 12 with α = n.

particle to the configuration integral without the shell particle. Therefore, we can again first determine the free energy of the system using eq 10, and then include the required increase in the free energy due to the presence of the shell particle. To determine this free energy correction, we begin with the configuration integral for the system without a shell particle Z=

1 n! (N − n)!

ρ+ = [(N − n)

∫V −v d3(N−n−1)r exp[−βU(r(N))]] (15)

A comparison of eqs 14 and 15 indicates that the free energy penalty for having a shell-particle is exactly equal to −kBT ln(4πλ2 Δλ ρ+). Adding this term is equivalent to multiplying the configuration integral of the system without a shell particle by the probability of finding one of the outside particles within the shell volume surrounding the embryo of radius λ. The value of the shell volume thickness was again taken to be the same as what was used in the simulations, or Δλ = 0.1. The free energy surface of forming the (n,v) embryos within the superheated LJ liquid using the modified constraint-DFT method is shown in Figure 4 for P = −0.45 and T = 0.80. The DFT results are reported for the same temperature as the simulations but at a slightly lower pressure, which corresponds to nearly the same bulk metastable liquid density as found in the simulations. As apparent in Figure 4, the DFT generated free energy surface is able to capture the same behavior as was observed in the simulation free energy surface, which is the emergence of a conduit that funnels small regions of vapor-like density up toward the critical embryo or free energy barrier. The size of the critical embryo is in good agreement with the simulation results but the DFT formulation predicts slightly larger radii of instability for larger values of n. The DFT formulation also slightly overestimates the free energy barrier in comparison with the simulation. Nevertheless, the relative shift of the free energies caused by the added constraint penalties seems to be in good agreement with the simulation results. Therefore, the scaling parameter of α = n appears to be a reasonable choice.

∫v d3nr ∫V −v d3(N−n)r exp[−βU(r(N))]

The ratio of the configuration integral for a system with a shell particle to one without a shell particle is therefore given by Zshell = [4πλ 2Δλ(N − n) Z

∫v d3nr ∫V −v d3(N−n−1)r exp[−βU(r(N))]] ∫v d3nr ∫V −v d3(N−n)r exp[−βU(r(N))]]

d3nr

∫v d3nr ∫V −v d3(N−n)r exp[−βU(r(N))]]

/[

(13)

/[

∫v

(14)

where the shell particle is held fixed at an arbitrary point within the spherical shell of volume 4πλ2Δλ. The various terms in the above equation are similarly defined as in section III.C, though here we do not include the constraint function since the constraints are applied after the fact in DFT. For the system without a shell particle, we note that the density profile of the particles surrounding the embryo will reach a limiting value at the surface of the embryo. This limiting value is the contact density ρ+, at the outside surface of the embryo, and is equal to the probability density that any one of the surrounding (N − n) particles is located at a fixed location within the spherical shell surrounding the embryo irrespective of the positions of the remaining particles (both those inside and outside of the embryo; due to spherical symmetry, this fixed location within the shell is arbitrary). In other words, ρ+ is given by30

V. CONCLUSIONS In the present work, we suggested an equilibrium embryo definition for the description of homogeneous bubble J



Article

(2) Dellago, C.; Bolhuis, P. G.; Csajka, F. S.; Chandler, D. Transition Path Sampling and the Calculation of Rate Constants. J. Chem. Phys. 1998, 108, 1964. (3) Bolhuis, P.; Chandler, D.; Dellago, C.; Geissler, P. Transition Path Sampling: Throwing Ropes over Rough Mountain Passes, in the Dark. Annu. Rev. Phys. Chem. 2002, 291. (4) Bolhuis, P. G.; Dellago, C.; Chandler, D. Sampling Ensembles of Deterministic Transition Pathways. Faraday Discuss. 1998, 110, 421. (5) Pan, A. C.; Chandler, D. Dynamics of Nucleation in the Ising Model. J. Phys. Chem. B 2004, 108, 19681. (6) Pan, A. C.; Rappl, T. J.; Chandler, D.; Balsara, N. P. Neutron Scattering and Monte Carlo Determination of the Variation of the Critical Nucleus Size with Quench Depth. J. Phys. Chem. B 2006, 110, 3692. (7) Allen, R. J.; Valeriani, C.; Rein Ten Wolde, P. Forward Flux Sampling for Rare Event Simulations. J. Phys.: Condens. Matter 2009, 21, 463102. (8) 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. (9) 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. (10) 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. (11) 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. (12) Gunther, L. A Comprehensive Treatment of Classical Nucleation in a Supercooled or Superheated Fluid. Am. J. Phys. 2003, 71, 351. (13) Uline, M. J.; Corti, D. S. Activated Instability of Homogeneous Bubble Nucleation and Growth. Phys. Rev. Lett. 2007, 99, 076102. (14) Uline, M. J.; Corti, D. S. Activated Instability of Homogeneous Droplet Nucleation and Growth. J. Chem. Phys. 2008, 129, 234507. (15) 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. (16) Ellerby, H. M.; Reiss, H. Toward a Molecular Theory of VaporPhase Nucleation 0.2. Fundamental Treatment of the Cluster Distribution. J. Chem. Phys. 1992, 97, 5766. (17) Senger, B.; Schaaf, P.; Corti, D. S.; Bowles, R.; Voegel, J. C.; Reiss, H. A Molecular Theory of the Homogeneous Nucleation Rate. I. Formulation and Fundamental Issues. J. Chem. Phys. 1999, 110, 6421. (18) ten Wolde, P.; Frenkel, D. Computer Simulation Study of GasLiquid Nucleation in a Lennard-Jones System. J. Chem. Phys. 1998, 9901. (19) Johnson, J. K.; Zollweg, J. A.; Gubbins, K. E. The Lennard-Jones Equation of State Revisited. Mol. Phys. 1993, 78, 591. (20) Chandler, D. Introduction to Modern Statistical Mechanics; Oxford University Press: New York, 1987. (21) 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. (22) Punnathanam, S.; Corti, D. S. Homogeneous Bubble Nucleation in Stretched Fluids: Cavity Formation in the Superheated LennardJones Liquid. Ind. Eng. Chem. Res. 2002, 41, 1113. (23) 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. (24) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications, 2nd ed.; Academic: San Diego, 2002. (25) Toxvaerd, S. Molecular-Dynamics Simulation of Homogeneous Nucleation in the Vapor Phase. J. Chem. Phys. 2001, 115, 8913. (26) Rosenfeld, Y. Free-Energy Model for the Inhomogeneous HardSphere Fluid Mixture and Density-Functional Theory of Freezing. Phys. Rev. Lett. 1989, 63, 980.

nucleation within a superheated liquid. The embryo definition was based on the following set of two order parameters: the number of particles n that happen to be inside a spherical volume v fixed at a given location within the liquid phase. The embryos were also identified with the additional criteria that all the n particles were constrained to be vapor-like and only a liquid-like particle was selected to mark the boundary of the spherical volume v. A set of MC simulation algorithms were developed to generate the free energy surface within the (n,v) order parameter space. The resulting free energy surface suggests that the process of homogeneous bubble nucleation begins with the appearance of small void regions that eventually grow into larger but compact and locally confined regions of vapor-like density. These vapor-like regions may grow large enough to reach the free energy barrier, and so become the critical nuclei. The limit of stability of the liquid phase traversed, shortly after the free energy barrier (which instigates the phase transition). In the companion paper (DOI 10.1021/jp404151h) the equilibrium embryo will be employed to determine the rate of bubble nucleation. The nucleation rate requires as input the concentration or number density of critical nuclei within the superheated liquid and how rapidly those nuclei transfer over the free energy barrier. As is discussed in detail in that paper, we show that some additional modifications of the (n,v) embryo are needed, based on dynamic information, before it can be properly used to calculate the nucleation rate. Since the (n,v) embryo definition is based on nonredundant counting of configurations about a fixed point in space, by dividing the total volume of a given snapshot of the system into nonoverlapping subvolumes each labeled with a unique value of the order parameters, we determine how the (n,v) free energy surface can be properly mapped into the average number density of the critical nuclei within the bulk metastable liquid. Furthermore, we discuss why the phase space of the metastable liquid must be further limited, and how the transitional region of the phase space (i.e., the critical nuclei within the metastable liquid phase) can be identified with a molecular dynamics committor probability analysis. This is in contrast to the more conventional equilibrium models of homogeneous nucleation, which associate the fully equilibrium-based critical nuclei to the saddle point of the free energy surface without considering the possible dynamics of the embryos.

■

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.

■

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

■

REFERENCES

(1) Debenedetti, P. G. Metastable Liquids: Concepts and Principles; Princeton University Press: Princeton, NJ, 1996. K



Article

(27) Yu, Y. X.; Wu, J. Z. Structures of Hard-Sphere Fluids from a Modified Fundamental-Measure Theory. J. Chem. Phys. 2002, 117, 10156. (28) Weeks, J. D.; Chandler, D.; Andersen, H. C. Role of Repulsive Forces in Determining Equilibrium Structure of Simple Liquids. J. Chem. Phys. 1971, 54, 5237. (29) Oxtoby, D. W.; Evans, R. Nonclassical Nucleation Theory for the Gas-Liquid Transition. J. Chem. Phys. 1988, 89, 7521. (30) McQuarrie, D. A. Statistical Mechanics; University Science Books: Sausalito, CA, 2000.

L


Toward a Molecular Theory of Homogeneous Bubble Nucleation: I

Recommend Documents