Free energy calculations along entropic pathways. II. Droplet nucleation in binary mixtures

The nucleation of liquid droplets is a ubiquitous phenomenon, central to many applications in chemistry, physics, and atmospheric sciences.^1–49 Nucleation from a single component system can be rationalized in terms of a change in a single intrinsic variable, such as the chemical potential or pressure. In this case, the supersaturation of the parent vapor phase is simply characterized by a given value of the pressure^35,48–53 (or equivalently of the chemical potential) that departs from the pressure at coexistence. On the other hand, nucleation from a mixture involves a myriad of pathways arising from the larger dimension of the system with, as additional intrinsic variables, the mole fractions for each of the components.^54–70 This also makes the definition of an appropriate reaction coordinate for the system especially challenging since, for instance, the total number of particles in the cluster is not the only significant variable, as the number of particles of each type also needs to be taken into account to fully characterize the nucleation process. Here, we propose and implement a new approach that uses the entropy S, which captures the interplay between size increase and molecular selectivity during nucleation, as the reaction coordinate for the process.

The aim of this work is to shed light on the entropic pathways followed during the nucleation of two binary mixtures. The first example we consider is a binary mixture of two highly miscible gases, specifically the Ar–Kr mixture. The second example involves a mixture of carbon dioxide with an alkane (here, as an example we consider C₂H₆) that is of technological relevance for the oil industry and for separation applications. We develop here the $μ1μ2VT–S$ method, where $μ1$ and $μ2$ are the chemical potentials for the two components of the mixture, and V and T are the volume and temperature of the system, to simulate the nucleation of a liquid droplet for these two binary mixtures. In Paper I,⁷¹ we discussed how, in the case of a single component system, the $μVT–S$ approach provided a direct connection with the classical nucleation theory,⁷² as any arbitrary value of $μ$ (or, equivalently, of the supersaturation $Δμ$⁠) could be applied, and allowed the system to overcome the free energy barrier of nucleation. In practice, this was achieved by using the entropy of the system, which can be readily calculated in the grand-canonical ensemble, as the reaction coordinate and by driving the system along an entropic pathway using the umbrella sampling simulation technique. Extending this approach to the case of mixtures is especially appealing since the grand-canonical ensemble naturally allows the number of molecules of each component to vary as nucleation proceeds. It is therefore very well suited to shed light on the impact of the choice of a given supersaturation (i.e., $Δμ1$ and $Δμ2$ for a binary mixture) on the free energy barrier of nucleation as well as on the interplay between the size of the droplet and its composition as nucleation takes place.

The paper is organized as follows. In Sec. II, we present the simulation method as well as the molecular models used in this work. We explain how we set up the $μ1μ2VT–S$ simulations, detailing how we proceed with the choice of chemical potentials, supersaturations, the range of entropies to be sampled, and the relation between the conditions of nucleation and the composition of the bulk. We then discuss the results obtained during the simulations of the nucleation process from supersaturated vapor phases of Ar–Kr and CO₂–C₂H₆. For all systems, we determine the free energy profile of nucleation and show how the choice of the conditions of nucleation impacts the height of the free energy barrier. We also focus on the analysis of the nucleation mechanism and on unraveling the interplay between size and composition during the formation of the liquid droplet, before finally drawing the main conclusions from this work in Sec. IV.

We extend to the case of mixtures the simulation method developed in Paper I.⁷¹ The approach proposed is termed as $μVT–S$ and consists in sampling configurations of the system around a value of the entropy S₀ in the grand-canonical (⁠$μVT$⁠) ensemble. As discussed in the first part of this series, the $μVT–S$ simulation method provides a direct connexion with the classical nucleation theory,⁷² since simulations of the nucleation process can be carried out for any value of the supersaturation $Δμ$⁠. Another advantage of this method is that it allows the calculation of the entropy of the system S. In the case of binary mixtures, we therefore define the method as $μ1μ2VT–S$ and evaluate the entropy during the simulations through

where $μ1$ and $μ2$ are the chemical potentials for the binary mixture, N₁ and N₂ are the number of atoms/molecules for each of the two components, and U is the internal energy for the entire system given by, in the case of a binary mixture of atoms,

where U_pot is the potential energy for the system. In the case of linear molecules, we add a contribution of k_BT per molecule to account for the rotational degrees of freedom.

The $μ1μ2VT–S$ simulation relies on gradually decreasing the entropy of the system through the application of a bias potential. This bias potential is defined within the framework of the umbrella sampling technique as

in which S₀ is the target value for the entropy, S is the current value of the entropy of the system, and k is a spring constant. This bias potential is then added to the potential energy of the system, and the total potential energy is used in the conventional Metropolis criteria for the acceptance of the different types of Monte Carlo (MC) steps. For the Ar–Kr mixture, MC steps include the insertion (12.5% of the attempted MC steps) or deletion (12.5% of the attempted MC steps) of atoms as well as the translation of a single atom (75% of the attempted MC steps). In the case of the C₂H₆–CO₂ mixture, we have the following rates: rotation (37.5% of the attempted MC steps), translation (37.5%), insertion (12.5%), and deletion (12.5%).

Successive umbrella simulations with decreasing values for the target entropy (S₀) are carried out to achieve the formation of a liquid droplet of a critical size. During each of these simulations, histograms for the number of times a given entropy interval is visited are collected, allowing for the calculation of the free energy profile associated with the nucleation process.^70,73–76

The simulation models used in this work for Ar and Kr are based on the Lennard-Jones potential,⁷⁷ with the following parameters for Ar: $σAr$ = 3.3952 Å  and $ϵAr/kB$ = $116.79K$ and for Kr: $σKr$ = 3.6274 Å and $ϵKr/kB$ = 162.58 K.

For the CO₂–C₂H₆, we use a force field that models the dispersion-repulsion interactions through an exp-6 functional form,^78–81 

In Eq. (4), r_m is the distance for which the potential reaches a minimum, r_max is the smallest positive distance for which du(r)/dr = 0, and $ϵ$ and $α$ are two potential parameters. As discussed by Errington and Panagiotopoulos,⁷⁹ it is convenient to discuss the potential parameters in terms of $σ$ (i.e., the distance for which u(r) = 0, obtained numerically by solving the equation $u(σ)=0$⁠) rather than in terms of r_m. In the case of C₂H₆, we use an united atom-type force field and model the molecule with two exp– 6 sites, each site standing for a CH₃ group. We use the following set of parameters: $ϵCH3/kB=129.64$ K, $σCH3=3.679$ Å, and $αCH3=16$⁠. In the case of CO₂, in addition to the repulsion-dispersion interactions, a Coulombic term is added to account for the quadrupolar nature of CO₂. The CO₂ molecule is modeled with a distribution of three exp – 6 sites and three point charges located on each of the atoms. We use the following parameters for the exp – 6 sites: $ϵC/kB=29.07$ K, $σC=2.753$ Å and $αC=14$⁠, $ϵO/kB=83.20$ K, $σO=3.029$ Å, and $αO=14$⁠. We also have q_C = 0.6466 e and q_O = $−03233$ e. Both molecules are considered to be rigid, with a distance between the two CH₃ exp – 6 sites set to 1.839 Å for ethane and a length of the C–O bond fixed to 1.1433 Å for CO₂. The calculation of the interaction energy is performed for distances up to 13.5 Å for both systems with electrostatic interactions calculated with the Ewald sum method beyond that distance using the parameters given in previous work.⁸² In line with previous simulation work on nucleation,³⁵ we do not include any tail corrections for the dispersion-repulsion interactions beyond the cutoff distance.

We carry out $μ1μ2VT–S$ simulations for the Ar–Kr mixture at T = 148.15 K in cubic cells with an edge of 100 Å  (the usual periodic boundary conditions are applied). We take advantage of the Expanded Wang-Landau (EWL) method we recently developed^82–85 to obtain very accurate estimates for $μ1$ and $μ2$⁠, both at the vapor-liquid coexistence and for supersaturated vapor phases, which will be the parent phase for the nucleation events. The values are presented in Table I.

We start by determining the conditions for the vapor-liquid coexistence at the coex I point. In practice, this is done using the EWL method by finding numerically the values for $μ1$ and $μ2$ which lead to equal probabilities for the vapor and the liquid phase (see more details in previous work⁸²). The EWL method also allows us to obtain all thermodynamic properties for the mixture including the mole fractions, pressure, and the entropies for the two coexisting phases. These entropies, which are also given in Table I, provide an idea of the range of entropies that need to be sampled for the system to undergo the vapor $→$ liquid transition. From the coexistence point, we can increase the value of the two chemical potentials $μ1$ and $μ2$ or, in other words, create a supersaturated vapor that will serve as a starting point for the nucleation process. We list in Table I the two sets of supersaturations (system 1 and system 2) generated from the coexistence point coex I. As can be seen from Table I, increasing $Δμ$ brings the supersaturated vapor more deeply into the liquid domain of the phase diagram, resulting in a larger value for the pressure and a lower value for the entropy of the liquid. We proceed along the same lines from the second coexistence point (coex II) to define two supersaturated vapors (system 3 and system 4), this time with a mole fraction in the liquid (x_Kr set to 0.250). We finally add that there are other ways of determining the chemical potential at the vapor-liquid coexistence and for supersaturated vapors.^78,86–97

We plot in Fig. 1 the successive umbrella sampling windows carried out during the $μ1μ2VT–S$ simulations on system 1. Fig. 1 shows the histograms corresponding to the probability according to which a given entropy interval is visited during the $μ1μ2VT–S$ simulations. Each of the umbrella sampling windows (labeled with an index i) is obtained by imposing through Eq. (3) a different value of the target entropy S_0,i. During the nucleation of a liquid droplet, the system goes from a supersaturated vapor with a low density, and thus of high entropy, to a system containing a liquid droplet of a critical size, i.e., to a much more dense system of lower entropy. To observe the formation of the liquid droplet, we therefore carry out successive umbrella sampling windows for decreasing values of the target entropy S_0,i and obtain the histograms, shown in Fig. 1, that cover the entire nucleation process. The progress of the system towards the formation of a liquid droplet can be followed by monitoring the relative location of S_max,i, the entropy for which the histogram p_i(S) reaches its maximum, and of S_0,i, the target entropy for the ith umbrella sampling window. We start with the window located to the right of Fig. 1, associated with the largest value of S_0,1 which corresponds to a very dilute vapor (S_0,1 = 0.11 (kJ/mol)/K). As shown in Fig. 1, for the first window (starting from the right), we have S_max,1 = 0.1097 (kJ/kg)/K, which is less than the target value S_0,1 = 0.11 (kJ/kg)/K. This means that the target value S_0,1 is greater than the entropy of the metastable supersaturated vapor for the choice of $(μ1,μ2)$ made for system 1. Gradually decreasing the target entropy for the next windows allows us to find the value S_0,i which coincides with the maximum for p_i(S). This occurs here for S_max,i = S_0,i = 0.1091 (kJ/kg)/K. At this point, we have the metastable supersaturated vapor. Then, during the next few umbrella sampling windows, we observe a change in behavior as the histograms p_i(S) now lag behind the target value for the entropy with $Smax,i>S0,i$⁠. This corresponds to the fact that the system has to overcome the free energy cost in forming the liquid droplet. Later on, for S_0,i = 0.103 (kJ/kg)/K, we find again that the maximum for p_i(S) coincides with S_0,i, indicating that we have reached the top of the free energy barrier of nucleation and that a liquid droplet of a critical size has formed. For target values of the entropy greater than 0.103 (kJ/kg)/K, we observe again a change in behavior as the histograms p_i(S) run ahead of the target value for the entropy with $Smax,i<S0,i$⁠, corresponding to the spontaneous growth of the liquid droplet.

To provide additional insight into the evolution of the average properties of the system for different umbrella sampling windows, we show in Fig. 2 the potential energy and the number of atoms for different values of the target entropy S_0,i = 0.107 (kJ/mol)/K and S_0,i = 0.104 (kJ/mol)/K during a production run. These plots show the impact of decreasing the target value for the entropy on the system. For example, comparing the results of two different umbrella sampling windows, 0.107 (kJ/mol)/K and 0.104 (kJ/mol)/K, we observe a decrease in the potential energy by about 60%. This happens simultaneously with an increase in the number of atoms of each component in the mixture by about 15% for Ar and by 20% for Kr. This confirms that the decrease in the target value for the entropy occurs with a greater organization and density of the system, which we will analyze in depth in Sec. III. We finally add that, for each umbrella sampling window, we run an equilibration run of $100×106$ MC steps, followed by a production run of $200$ × $106$ MC steps. Throughout the $μ1μ2VT–S$ simulations, we also check that the acceptance rate for the insertion/deletion steps remains high enough to ensure an accurate sampling of the configurations of the systems. For instance, in the case of system 1 and for the configurations of highest density (umbrella sampling window with a critical liquid droplet for S_0,i = 0.103 (kJ/mol)/K), the acceptance rates for the insertion/deletion steps are of 45.8%.

Simulations of droplet nucleation for the C₂H₆–CO₂ mixture are performed at T = 263.15 K in cubic cells, with an edge of 100 Å, and with the usual periodic boundary conditions. We use EWL simulations to determine the other input parameters for the $μ1μ2VT–S$ simulations and list in Table II the sets of chemical potentials $(μ1,μ2)$ used for these simulations.

We identify a first state point, coex III, leading to vapor-liquid coexistence for the C₂H₆–CO₂ mixture. As previously discussed, this is achieved by finding numerically $μ1$ and $μ2$ such that the liquid and the vapor phases are equally probable. The EWL simulations also yield the mole fractions as well as the entropies of the two coexisting phases, bracketing the range of entropies needed to sample the nucleation of the liquid droplet. From coex III, we increase the chemical potentials of the two mixture components by $Δμ1$ and $Δμ2$ to obtain thermodynamic conditions located in the domain of the liquid in the phase diagram and, for which, we can observe a supersaturated vapor. We repeat this step for two different supersaturations leading to systems 5 and 6 in Table II (as for Ar–Kr, the choices for $Δμ1$ and $Δμ2$ are made such that the liquid mole fraction in the second component, here CO₂, remains constant). We choose another coexistence point, coex IV, and two corresponding supersaturated vapors, system 7 and system 8. For each umbrella sampling window, an equilibration run of $50×106$ MC steps is run, followed by a production run of $100×106$ MC steps. As with the Ar–Kr mixture, we check that the acceptance rates for the insertion/deletion steps remain high enough to ensure an accurate sampling. This is the case here for all umbrella sampling windows. For instance, in the case of system 5 and for the window associated with the highest density (umbrella sampling window with a critical liquid droplet for S_0,i = 0.155 (kJ/mol)/K), the acceptance rates for the insertion/deletion steps are of 36.9% for C₂H₆ and of 46.6% for CO₂.

We start by analyzing the free energy barriers obtained for the Ar–Kr mixture at T = 148.15 K. The left panel of Fig. 3 shows the free energy profiles for two different supersaturations, system 1 and system 2, at x_Kr = 0.44. The supersaturation has a direct effect on the height of the free energy barrier, with a barrier of $25±2kBT$ for system 1 and of $15±2kBT$ for a system 2. $Δμ1$ and $Δμ2$ are greater for system 2, which means that the parent supersaturated vapor of system 2 is located more deeply into the domain of the liquid in the phase diagram, and, as such, that the nucleation of a liquid droplet occurs more easily. This also leads to a lower free energy of nucleation than for system 1. We find that the range of entropies spanned during the nucleation process is also impacted by the amount of supersaturation. For instance, looking at system 2 in Fig. 3, we find that, for the higher supersaturation, the passage from the supersaturated vapor to a system containing a droplet of a critical size occurs over entropies between 0.1085 (kJ/mol)/K,

for which the free energy reaches a minimum for the entropy of the metastable parent phase, and 0.1055 (kJ/mol)/K, for which the free energy reaches a maximum corresponding to the formation of a droplet of a critical size. On the other hand, for system 1 (lower supersaturation), we find that liquid nucleation takes place over a much broader entropy range, with the entropic pathway ranging from 0.1091 (kJ/mol)/K to 0.103 (kJ/mol)/K. This can be attributed to the combination of two effects. For a higher supersaturation, the density for the parent supersaturated vapor is larger and, therefore, its entropy is lower. Furthermore, at high supersaturation, the thermodynamic conditions lie further inside the liquid domain. Thus, the critical size for the liquid droplet becomes smaller, resulting in a system at the top of the free energy barrier that has higher entropy. Both effects account for the narrowing, at high supersaturation, of the entropy range spanned during the droplet nucleation along the entropic pathway.

We now turn to the second set of simulations carried out for system 3 and system 4 at x_Kr = 0.25. The free energy profiles obtained for the two supersaturations are shown in the right panel of Fig. 3. The behavior observed for x_Kr = 0.25 exhibits similar qualitative features as for x_Kr = 0.44. We find that the height of the free energy barrier of nucleation increases as we go from system 4 to system 3 (as the supersaturation decreases), with the free energy of nucleation increasing from $7±1kBT$ to $21±2kBT$⁠. The range of entropies spanned during nucleation is also found to increase as the supersaturation decreases. For system 4, we observe that the entropy of the parent phase is of 0.1047 (kJ/mol)/K while the entropy at the top of the free energy barrier is of 0.103 (kJ/mol)/K. For a lower supersaturation (system 3), the entropy of the supersaturated vapor is 0.1055 (kJ/mol)/K while it is of 0.100 (kJ/mol)/K at the top of the free energy barrier. This means that a significantly broader range of entropies is sampled during the nucleation process. There are, however, notable differences between the two plots, which result from the interplay between the pressure and the chemical compositions of the system. Considering the results at fixed supersaturation, we find that the free energy of nucleation decreases by 16% (going from system 1 to system 3), while it decreases by 53% between system 2 and system 4. Similarly, the dependence of the height of the free energy barrier on supersaturation is also shown to be impacted, with a decrease by 40% from system 1 to system 2 at x_Ar = 0.56 and a much larger decrease by 66% from system 3 to system 4 at x_Ar = 0.75. Both findings can be attributed to the fact that liquid mixtures with a lower fraction of Kr are obtained for higher pressures. This, in turn, implies that the corresponding parent phases will be supersaturated vapors of larger densities and, therefore, result in lower free energy barriers of nucleation.

The plot shown in Fig. 3 for the free energy profile of nucleation is a projection of the multi-dimensional free energy surface, as discussed in prior work on the nucleation in binary systems.⁹⁸ Here, we do not take into account any non-isothermal effect during nucleation.^99,100 To interpret further the results obtained for the free energy of nucleation as a function of S, we now focus on the interdependence between the entropy S, the droplet size, and its composition. As we have seen from Fig. 2, the number of atoms steadily increases as the target entropy S_0,i is decreased. To check that the decrease in entropy undergone by the system leads to an increased organization within the system and to the formation of a liquid droplet, we show in Fig. 4 snapshots of the system, obtained during umbrella sampling simulations for decreasing values for the target entropy. As S_0,i decreases, the droplet size increases, with a larger number of atoms being incorporated to the droplet as the target entropy decreases from 0.1075 (kJ/mol)/K, to 0.105 (kJ/mol)/K, and finally to 0.103 (kJ/mol)/K. To assess further this point, we perform a detailed analysis of the size and composition of the droplet along

the entropic pathway. The atoms belonging to the droplet are identified through a commonly used geometric criterion.³⁵ For this purpose, we determine the distribution for the number of neighbors within a distance of 5.4 $Å$ of a central atom in the vapor and the liquid mixture. We show in Fig. 5 these distributions. Both distributions are sharply peaked around N_nab=1 for the vapor and around N_nab=10 for the liquid. As shown in Fig. 5, atoms belonging to the vapor always have less than 6 neighbors within a distance of 5.4 Å. This allows us to introduce the following condition to identify a liquid-like atom, or equivalently an atom belonging to the developing droplet, as an atom with at least 6 nearest neighbors within 5.4 Å.

Applying this analysis to the configurations generated during the umbrella sampling windows leads to the determination of the evolution of the size of the droplet during the nucleation process. Furthermore, by keeping track of the identity (either Ar or Kr) of the atoms belonging to the droplet, we can also shed light on the chemical selectivity during the formation of the cluster. This means that the interplay between droplet size and composition during the nucleation process can be directly accessed during the simulations, and, in turn, shed light on the departure in composition of the critical droplet from the bulk composition^55,64 and the possible onset of phase separation in partially miscible mixtures.^57,59,60

We show on the left of Fig. 6 the results obtained for the two supersaturations of systems 1 and 2. In both cases, throughout the nucleation process, the overall size of the droplet is shown to gradually increase as the entropy of the system decreases. The critical size of the droplet is found to be larger at low supersaturation (⁠$Nc=713±60$ for system 1 and $Nc=323$ ± 30 for system 2). It is also reached for a lower value of the entropy at low supersaturation (S_c = 0.103 (kJ/mol)/K for system 1) than at high supersaturation (S_c = 0.1055 (kJ/mol)/K for system 2), in line with the results obtained for the free energy barrier in Fig. 3. A closer inspection of the variation for the number of each type of atoms shows that the composition of the droplet does not remain the same during the entire nucleation process. For small sizes (high entropy), the droplet is richer in Ar with $55%±5%$ of Ar atoms in the droplet at S = 0.1073 (kJ/mol)/K for system 1. There is then a crossover at S = 0.1067 (kJ/mol)/K for which the two types of atoms are equally present. As the entropy further decreases, the droplet becomes richer in Kr atoms (with a fraction of $48%±4%$

in Ar for a droplet of a critical size). Looking at the composition of the droplet past the critical size, we find that the fraction of Ar atoms in the droplet increases again and reaches $49%$ ± $4%$ for S = 0.102 (kJ/mol)/K. The same mechanism is observed for system 2, with a crossover point located at S = 0.1068 (kJ/mol)/K and a fraction of Ar atoms of $47%$ ± $4%$ in the critical droplet. These fluctuations in the nucleus composition through the nucleation process, as well as the departure of the composition of the droplet from the composition of the liquid phase, are consistent with prior work on nucleation, based on a revision of the classical nucleation theory,¹⁰¹ on classical density functional theory calculations,⁵⁵ on molecular simulations,^60,62,64 or on approaches based on macroscopic kinetics.⁶⁸ Our results indicate that the mole fraction in Ar in the nucleus must be less than for the bulk by 0.08 for system 1 and of 0.03 for system 2 with respect to the bulk composition. These deviations, which are moderate due the almost ideal nature of the Ar–Kr mixture as noted by Zeng and Oxtoby,⁵⁵ did not give rise to a significant phase separation effect as reported in simulations of more strongly asymmetric mixtures.^60,102

Looking now at the evolution of the size of the droplet as a function of entropy for the other set of conditions (x_Ar = 0.75), we see that the size of the droplet steadily increases as the entropy of the system decreases. The droplet reaches a size of $629±90$ atoms for the critical droplet in the case of system 3, and a size of $273±50$ for system 4. As for the previous system, the critical droplet has formed once the entropy has reached a critical value S_c = 0.100 (kJ/mol)/K for system 3, a value that is notably lower than its counterpart of 0.103 for system 4. This reflects the fact that the range for the entropies spanned during the nucleation event becomes narrower and narrower as the supersaturation is increased. Unlike for systems 1 and 2, we do not observe any crossover between the mole fractions in Ar and Kr in the droplet, given the very large fraction of Ar in systems 3 and 4. As for systems 1 and 2, we observe, however, fluctuations in the composition of the droplet with the fraction of Ar atoms in the droplet varying between 60% and 80% during nucleation.

Turning to the results obtained for the CO₂–C₂H₆ mixture, we plot in Fig. 7 the free energy barrier of nucleation obtained for systems 5 and 6. For a liquid fraction of $xCO2=0.053$ (left panel of Fig. 7), we observe that increasing the supersaturation leads to a decrease in the height of the free energy barrier ofnucleation from $21±2kBT$ (system 5) to $12±1kBT$ (system 6). Similarly, for a mole fraction of $xCO2=0.097$⁠, a lower supersaturation results in a barrier of $12±1kBT$ (system 7), while a higher supersaturation yields a free energy of nucleation of $7±1kBT$ (system 8). The range of entropies spanned along the nucleation pathway is also found to depend strongly on supersaturation, and becomes broader at lower supersaturations. For instance, for system 5, we obtain an entropy for the parent (supersaturated vapor) phase of 0.161 (kJ/mol)/K and an entropy of the system of S_c = 0.151 (kJ/mol)/K, when a droplet of a critical size has formed. Increasing the supersaturation (system 6) leads to a narrower entropy range, most notably as a result of a decrease in the entropy S_c = 0.154 (kJ/mol)/K for which the critical size of the droplet is reached. Similar conclusions apply for systems 7 and 8. We find a narrowing of the entropy range spanned during nucleation at high supersaturation, with the top of the free energy barrier droplet reached for a higher entropy (S_c = 0.157 (kJ/mol)/K) than at low supersaturation (S_c = 0.155 (kJ/mol)/K). The higher entropy required to form a droplet of a critical size at high supersaturation, together with the lower free energy barrier of nucleation obtained at high supersaturation, is most likely due to the smaller size of the critical droplet at high supersaturation,a point that we discuss in Sec. IV in our analysis of the relation between size, selectivity, and entropy during the nucleation process.

We show in Fig. 8 the snapshots of configurations of system 5, obtained throughout the nucleation process. As for the Ar–Kr system, we start from a metastable supersaturated vapor and carry out umbrella sampling windows with decreasing values for the target entropy. During the first few umbrella sampling windows, the target entropy is high enough, so that the droplet is fairly small (see the left of Fig. 8 for S_0,i = 0.157 (kJ/kg)/K). Then, as S_0,i is decreased further,the droplet starts to become larger and larger, and eventually reaches its critical size for S_0,i = 0.151 (kJ/kg)/K. This shows that decreasing the value for the target entropy not only increases the density of the system but also results in an increased level of organization with the formation of the liquid droplet.

As for the Ar–Kr mixture, Fig. 7 shows a projection of the multi-dimensional free energy surface in the (entropy, free energy) plane. We pursue our analysis by characterizing the interdependence between the entropy S, the droplet size, and its composition. We start by determining which molecules in the system have a liquid-like environment and, as such, belong to the droplet. We determine, for each molecule, the distributions, shown in Fig. 9, for the number of neighbors for the vapor and the liquid. Neighboring molecules are defined as being separated by a distance (between the centers of mass of the 2 neighboring molecules) less than 6.4 Å. The distributions obtained for the vapor and liquid reach their maxima for very different numbers of neighbors (for, on average, a single neighbor in the case of the vapor and for 9 neighbors in the case of the liquid). This allows us to define a molecule as having a liquid-like environment and thus belonging to the droplet, if it has 6 or more neighbors within a spherical shell of 6.4 Å.

We now move on to the analysis of the size of the droplet as a function of the entropy of the system throughout the nucleation process. Fig. 10 shows that, for all systems, the total number of molecules within the cluster increases smoothly as the entropy of the system decreases. Furthermore, we find that the size of the critical droplet decreases as supersaturation is increased. For a liquid mole fraction of $xCO2=0.053$⁠, the critical size for the droplet is of $Nc=715±35$ molecules for system 5 and of $Nc=422±30$ molecules for system 6. The

smaller size of the critical droplet at high supersaturation (system 6) accounts for the higher value of the entropy for which the system reaches the top of the free energy barrier. Similarly, when the liquid mole fraction $xCO2=0.097$⁠, the critical size is of $Nc=494±28$ molecules for a low supersaturation (system 7) and of $Nc=303±21$ molecules for a high supersaturation (system 8). The smaller critical size and higher critical entropy obtained at the higher supersaturation are, once again, found to be consistent with the free energy plot of Fig. 7. Turning to the composition of the droplet, we find that C₂H₆ remains predominant throughout the nucleation process for all systems (see Fig. 10). We also find, however, that the composition of the nucleus depends on its size, as nucleation starts with the formation of a droplet that has a higher CO₂ mole fraction than the bulk. It is around 8% for systems 5 and 6 for droplets containing a total of 50–100 molecules. Similarly, considering the same droplet sizes, it is of about 17% for systems 7 and 8. The fraction of CO₂ then decreases as the size of the droplet increases. For droplets of a critical size, the fraction of CO₂ is of 6% for systems 5 and 6, while it is of 13% for systems 7 and 8. Despite the small sizes of the critical droplets, which contain only a few hundred of molecules, the fractions in the critical droplets are reasonably close to the CO₂ mole fraction of the liquid, and the departure from the bulk compositions (5.3% of CO₂ for systems 5 and 6, and 9.7% of CO₂ for systems 7 and 8) is small. As for binary mixtures of atoms, the departures in the droplet composition with respect to that of the bulk are consistent with prior simulations of droplet nucleation in binary molecular systems (see e.g. recent simulations of the methane-nonane system⁶⁴).

In this work, we propose a new simulation method to study the nucleation process in binary mixtures of atomic fluids (Ar–Kr) and of molecular fluids (C₂H₆–CO₂). The method is based on driving the formation of a liquid droplet through a series of umbrella sampling simulations where the bias potential is a function of the entropy S of the system. The resulting approach is implemented within the grand-canonical ensemble and, since the entropy serves as the reaction coordinate for the nucleation process, it is called $μ1μ2VT–S$⁠. The application of the method to the formation of liquid droplets in binary mixtures sheds light on the interplay between the size of the droplet, its composition, and the supersaturation at which the nucleation process occurs. Our findings show that, at low supersaturation, the range of entropies spanned by the nucleation process becomes broader, as a result of the combined effect of the larger entropy of the metastable supersaturated vapor (parent phase) and of the lower entropy associated with the configurations of the system that contain a liquid droplet of a critical size. These simulations allow us to characterize the critical droplet in terms of a critical value reached by the entropy at the top of the free energy barrier of nucleation. We are also able to obtain the free energy profile along the entropic pathway underlying the formation of the liquid droplet and to calculate the free energy of nucleation as a function of the supersaturation and chemical composition of the system. The analysis of the composition of the droplet shows that the mole fractions fluctuate throughout the nucleation process and depart from the composition of the bulk. This departure is however found to become less and less significant as the size of the droplet increases and its composition starts to conform more and more to that predicted by thermodynamics. Finally, while the $μ1μ2VT–S$ does not yield directly the nucleation rate, the method allows to generate and stabilize configurations of the system close to the top of the free energy barrier. As discussed in previous work,¹⁹ the nucleation rate can be obtained by carrying out additional molecular dynamics simulations, using configurations close to the top of the free energy barrier as a starting point, and following the Bennett-Chandler scheme^103–105 to determine the kinetics of the process. Alternatively, the threshold method of Yasuoka and Matsumoto can also be used to determine the nucleation rate.¹ 

Partial funding for this research was provided by NSF through CAREER Award No. DMR-1052808.