Estimation of bubble cavitation rates in a symmetrical Lennard-Jones mixture by NVT seeding simulations

Understanding the formation of bubbles from metastable liquids is a problem of high interest in many fields, such as sonochemistry¹ or geology.² Experiments often provide information about the cavitation rates, i.e., about the number of critical bubbles formed per unit of time and per unit of volume, but not about the microscopic mechanism that leads to their formation. With the aim of understanding this process, the Classical Nucleation Theory^3,4 (CNT) was proposed a century ago, initially to rationalize the condensation of supersaturated vapors, but later it was also applied to other transitions, including the cavitation of bubbles.⁵ 

Molecular simulations have also been extensively used to study the nucleation of bubbles.^6–14 Experiments are often performed at conditions in which the nucleation rates are too low to be computed with the typical simulation times; therefore, special simulation techniques such as umbrella sampling,^15–18 forward-flux sampling,^6,7,19 and transition path sampling^20–22 have been employed. Even though these methods are rigorous, they are often technically challenging and computationally demanding. As an alternative, during the last few years, the seeding method has been gaining popularity, mainly due to the lower computational cost and easier implementation than the preceding methods, with a focus not only on bubble cavitation^9–11,13 but also on other transitions, such as crystallization or condensation.^23–27 This method consists in inserting a seed of the stable phase within the metastable liquid and searching for the conditions for which the inserted seed is critical, i.e., it stands at the top of the free energy barrier that separates the metastable and stable phases. For one-component systems, this method can be implemented both in the NpT ensemble, in which the critical nucleus is in unstable equilibrium, and in the NVT ensemble, in which it is in stable equilibrium, thus enabling the study of the properties of the critical cluster over longer times.^10,11,13 The information about the critical cluster size is then combined with expressions from CNT to calculate the cavitation rates. The application of NVT and NpT seeding to study the cavitation of bubbles in the truncated and force-shifted Lennard-Jones (TFS-LJ) liquid has provided rates consistent with those obtained from other simulation techniques that do not depend on CNT.^9–11 

The main goal of this article is to extend the method of NVT seeding to study the nucleation of bubbles in binary systems. The applicability of NVT seeding has already been demonstrated for other transitions besides cavitation, in particular for crystal nucleation, but always for one-component systems.²⁸ NpT seeding was not used to study bubble nucleation in multi-component systems, but it has been successfully applied to study the nucleation of crystals in several multi-component systems, such as, for example, the crystallization of NaCl and ice from aqueous solutions^29–31 and the nucleation of methane hydrates.³² In this work, the applicability of NVT seeding to multi-component systems will be investigated by studying bubble nucleation in a partially miscible binary symmetrical LJ mixture. This is one of the simplest mixtures one can think of, and this makes it an ideal model system in which to investigate the applicability of NVT seeding before addressing nucleation in more complex binary systems.

Simulation studies of cavitation of bubbles from binary mixtures are rather scarce, and often this phenomenon has been studied theoretically, for example, in the framework of density functional theory or in the framework of CNT.^33–36 However, the extension of CNT to multi-component systems has not been straightforward.^37–39 Several of the proposed formulations lead to unphysical results, especially for non-ideal mixtures, which has been ascribed to the fact that the effect of the curvature on the interfacial free energy cannot be neglected.³⁹ The nucleation of bubbles in binary mixtures is very rich and involves many different scenarios depending on the interaction parameters. Curious phenomena have been found, such as, for example, an anomalous decrease in the nucleation rate when temperature increases at constant pressure in mixtures of species of different volatility.³⁴ However, as our main goal is to illustrate the applicability of NVT seeding, we will focus on the case of a partially miscible symmetrical LJ mixture at a temperature at which the liquid is miscible. In these conditions, we expect that the cavitation behavior follows a similar mechanism as that in the pure system.

The model system is a binary symmetrical mixture of equal-sized spheres. The interaction between any pair of molecules is described by the truncated and force-shifted Lennard-Jones model (TSF-LJ),⁶ 

where r is the distance between two molecules, $ULJr$ is the standard 12-6 LJ potential and $ULJ′r$ is its first derivative, and r_c = 2.5σ is the cut-off distance. We have chosen this short-range model to avoid dealing with long range corrections, which are problematic in heterogeneous systems with large variations in density, as the liquid–bubble systems considered in this work. Besides, there are numerous studies of cavitation in pure TSF-LJ that can be compared with the results for the binary symmetrical mixture.^6,7,9–11 Here we consider a mixture of spheres of equal size, σ_AA = σ_BB, and equal mass, m_A = m_B.

The phase diagram of partially miscible symmetrical TSF-LJ mixtures (σ_AA = σ_BB = σ_AB ≡ σ and ϵ_AA = ϵ_BB ≡ϵ, and ϵ_AB/ϵ less than but close to 1) has not yet been reported to the best of our knowledge. However, for the standard LJ potential, the global phase diagram (which delimits regions of parameter space with the same phase diagram topology) has been calculated⁴⁰ from the equation of state proposed by Johnson et al.⁴¹ These authors found that symmetrical LJ mixtures with ϵ_AB/ϵ > 0.75 exhibit a type II diagram according to Bolz’s classification.⁴² This means that in a pressure vs temperature projection, at low temperature, there is a liquid–liquid–vapor three coexistence line that, as T increases, transforms first into a critical azeotropic line and then into a critical liquid–vapor line. According to Ref. 40, for a mixture with ϵ_AB/ϵ = 0.8, T* = k_BT/ϵ = 0.785 (which is the T* for which cavitation was studied for the pure TSF-LJ liquid in previous studies^10–12) falls within the mixing region. As our intention is to study the cavitation of bubbles from a homogeneous binary liquid mixture, we set ϵ_AB/ϵ = 0.9 so that the liquid does not undergo demixing at the considered T*. In the seeding simulations performed in this work, the concentration of the liquid is set to the molar fraction x_A,l = 0.75, where the subscript l is added to differentiate the composition of the liquid from that of the bubble.

Molecular dynamics simulations were performed with the GROMACS 4.6.7 package (double precision).⁴³ To implement the chosen model system, we used the argon parameters, i.e., the diameter of the particles is set to σ = 3.405 Å, the energy parameter to ϵ = 0.996 078 kJ mol⁻¹, and the mass is m = 39.946 u. In the remainder of this article, all the quantities will be presented in LJ reduced units: the volume V* = V/σ³, the density ρ* = (N/V)σ³ (N is the number of particles in the system), the temperature T* = k_BT/ϵ, the pressure p* = pσ³/ϵ, and the time t* = t/τ, with $τ=mσ2/ϵ=$ 2.156 ps. In these units, p* = 1 corresponds to p = 419 bar.

The pressure was controlled using an isotropic Parrinello–Rahman barostat,⁴⁴ with a relaxation time of t = 2 ps (t* ≈ 0.92) and a compressibility of 7 × 10⁻⁵ bar⁻¹. For temperature coupling, velocity rescaling was used with a stochastic term⁴⁵ and a relaxation time of t = 1 ps (t* ≈ 0.46). The temperature was fixed at T* = 0.785. The equations of motion were integrated using the velocity-Verlet algorithm with a time step of t = 0.005 ps (t* ≈ 0.002).

According to CNT, the nucleation rate, i.e., the number of nuclei from the stable phase that form in the metastable phase per unit of volume and per unit of time at constant pressure, temperature, and composition, can be calculated from^{3,4,37,39,46,47}

where k_B is the Boltzmann constant, ΔG_c is the Gibbs free energy barrier for the formation of a critical bubble, and K is the kinetic prefactor. ΔG_c can be computed as

where V and A are, respectively, the volume and the area of the surface of the critical bubble, γ is the interfacial free energy of the critical bubble, and Δp = p_bubble − p_l is the difference in pressure between the critical bubble and the surrounding liquid. Assuming a spherical critical bubble of radius R_c, and using the Laplace equation

the free energy barrier can also be expressed as

Note that this expression for ΔG_c is formally identical to that for pure systems.^9,46 The main difference with respect to mono-component systems is the procedure used to calculate Δp. In equilibrium, the chemical potential of the liquid and the bubble is the same. In previous studies of pure systems, the pressure inside the bubble was estimated by integrating the equation of state of a bulk vapor to find the pressure at which the chemical potential of the bubble is the same as that of the surrounding liquid.^9–11 This calculation is a bit more complex in the mixture because we have to find the pressure and composition at which the chemical potentials of the two components are equal in both phases. It is also important to point out that the radius of the critical cluster R_c that enters in Eq. (5) is the one that satisfies the Laplace equation [Eq. (4)], i.e., the radius at the surface of tension where the value of γ is a minimum with respect to formal changes of the dividing surface.⁴⁸ The kinetic factor can be calculated from^34,46

where m_A and m_B are the masses of the particles of types A and B, x_A,v and x_B,v are their vapor molar fractions, and ρ_l is the number density of the metastable liquid. In our system, the mass of both particle types is the same (m_A = m_B = m), and the kinetic prefactor simplifies down to the same expression as that used for one-component systems,⁵ 

In the seeding method,²⁷ a spherical seed from the stable phase is inserted into the equilibrated metastable phase. For one-component systems, simulations of this heterogeneous system can be performed either in the NpT or NVT ensembles.^10,28 In pure systems, the critical cluster corresponds to a maximum in G at constant N, p, and T and to a minimum in F at constant N, T, and V (corresponding at the same p).^10,49 Therefore, it is possible to have the critical bubble in stable equilibrium with the liquid, allowing one to study its properties over longer periods of time. The main goal of this paper is to show the applicability of NVT seeding for binary systems. In our simulations, the number of particles of each type is kept constant, so strictly speaking, we are working in the N_AN_BVT ensemble or the NVTx_A ensemble, with x_A being the overall molar fraction of component A. However, for simplicity, we will use the notation NVT (and NpT), keeping in mind that simulations are performed at constant overall composition. The critical bubble size obtained from seeding will be used together with Eq. (5) to estimate the free energy barrier and the nucleation rate using Eq. (2). As the results of seeding are very sensitive to the size of the critical cluster, it is extremely important to choose a criterion to determine the bubble sizes that gives rates consistent with those derived from methods that do not depend on it (Refs. 9 and 10). In this work, we will use brute force simulations to test two criteria to define the bubble sizes.

The first step to study the cavitation in the symmetrical TSF-LJ mixture is to locate the liquid–vapor coexistence point at the chosen concentration, x_A,l = 0.75. To facilitate the comparison with cavitation in the pure TSF-LJ liquid, simulations are performed at the same reduced temperature as in previous studies (T* = 0.785).^10,11 The coexistence point was evaluated by two independent routes: by direct coexistence simulations^50–52 and by calculation of the chemical potential of the two components in the two phases and imposing the condition of phase equilibrium.

In order to locate the liquid–vapor equilibrium using direct coexistence simulations, we first equilibrate a liquid composed of 2370 particles of type A and 790 particles of type B (i.e., a liquid with a molar fraction x_A = 0.75) at T* = 0.785 in the NpT ensemble. The approximate simulation box dimensions are 4.32 × 4.32 × 8.64 nm³. Once the liquid is equilibrated, we introduce an empty space next to the liquid of length about twice its longer dimension (the size of the simulation box is 4.32 × 4.32 × 30 nm³ when the empty space is inserted). This system is then allowed to evolve in the NVT ensemble.

Once the system reaches equilibrium, the coexistence pressure is estimated from the normal component of the pressure tensor, which for the conditions just described is $plv*=$ 0.0334(3). The equilibrium composition of each phase is calculated by measuring the total and partial density profiles along the direction perpendicular to the interface. We obtained that, after equilibration, the concentration of particles of type A in the liquid phase increased up to about x_A,l = 0.77(1) and that of the gas phase decreased down to about x_A,v = 0.69(1). This enrichment of component A in the liquid phase can be easily understood by energetic arguments. As crossed interactions are lower than those between like particles, and interactions are much more relevant in the liquid phase than in the gas phase (in which particles are on average far away from each other), it is energetically favorable to reduce the number of crossed interactions in the liquid by releasing to the vapor a larger proportion of particles of type B.

As all the cavitation studies presented in this work were performed for a liquid with composition x_A,l = 0.75, we repeated the direct coexistence simulations to try to locate the liquid–vapor transition when the liquid has a composition closer to this value. As the concentration of A in the liquid, x_A,l, tends to increase when brought into contact with vacuum, we now start the direct coexistence simulations from a liquid phase containing 2323 particles of type A and 837 particles of type B, so that the initial composition of the liquid is x_A = 0.735, i.e., somewhat below the target x_A,l = 0.75. This system equilibrates at a slightly lower coexistence pressure than before, $plv*=0.0341(3)$⁠, but with the desired liquid equilibrium concentration of A, namely, x_A,l = 0.75(1) and x_A,v = 0.68(1) in the vapor. Finally, we also checked that the coexistence pressure coincides with that obtained using the equations of state (EOS) of the bulk liquid and vapor phases together with the densities and compositions of the liquid and vapor phases extracted from the bulk region of each phase (i.e., sufficiently away from the interface) in the direct coexistence simulations. Using this procedure, we get $pl*=0.033(1)$ from the EOS of the bulk liquid and $pv*=0.034(1)$ from the EOS of the bulk vapor. Both values are in very good agreement with each other and also with the estimate from the normal component of the pressure tensor.

The interfacial free energy, γ, at coexistence can be calculated from the average pressure tensor,^48,53–58

where L_z is the length of the simulation box in the direction normal to the interface, $p̄zz$ is the average of the normal component of the pressure, and $p̄xx$ and $p̄yy$ are the averages of the tangential components of the pressure. The values of the average components of the pressure tensor do not depend on any particular formulation of the local pressure tensor.^53,57,58 The factor 2 in the denominator of Eq. (8) takes into account that two liquid–vapor interfaces are formed by the periodicity of the simulation box.

The coexistence conditions as well as the interfacial free energy for the TSF-LJ mixture are gathered in Table I. For completeness, we also performed simulations for the pure TSF-LJ system at T* = 0.785, and these results are compared with data from previous studies^6,11 in the same Table. The first observation is that the liquid–vapor coexistence and the interfacial free energy obtained by us for the pure system are fully consistent with those found in the literature.^6,11 Comparing now the data for the pure and the symmetrical TSF-LJ mixtures at the same temperature, it can be seen that the liquid–vapor coexistence pressure of the symmetrical mixture is higher than that of the pure TSF-LJ system, which is in keeping with the trends observed in Ref. 40. Finally, the interfacial free energy is appreciably lower in the symmetrical TSF-LJ mixture than in the pure system. The same qualitative behavior was found in previous work,⁵⁹ but a quantitative comparison cannot be made as this work considered a mixture with a different composition (x_A = 0.5) and used a different truncation of the interatomic potential.

An alternative method to calculate the liquid–vapor coexistence point is to find the thermodynamic conditions at which the chemical potential of the two components is equal in both phases at the same temperature T and at the same pressure p,

where μ_i,α(p, T, x_i,α) and x_i,α are, respectively, the chemical potential and molar fraction of component i (i being an A or B component) in phase α (α being the liquid or vapor phase).

The excess chemical potential of components A and B in a liquid mixture with x_A,l = 0.75 and in the vapor with x_A,v = 0.60, 0.65, and 0.70 as a function of pressure were calculated using the particle test Widom method^60,61 as implemented in GROMACS.⁴³ For that aim, we performed NpT simulations of each phase in a cubic simulation box containing 1000 particles. Typically, the chemical potential was averaged over 2000 configurations, collected over a MD simulation run consisting of 2 × 10⁷ steps, after reaching equilibrium. In each configuration, 10 million particle insertion tests were performed. The total chemical potential was obtained by adding to the excess chemical potential the ideal contribution, which for component i is given by μ_i,id/(k_BT) = ln(x_iρσ³), where x_i is its molar fraction. The thermal de Broglie length is the same in both phases, and its value does not change the chemical potential difference between them. Here it was assigned to the arbitrary value σ.

The chemical potential of the two components in the liquid and vapor phases is shown in Fig. 1(a). As can be seen, the chemical potential of component A in the liquid with x_A,l = 0.75 (black lines and symbols) crosses the chemical potential of component A in the vapor with concentrations x_A,v = 0.60 (red), 0.65 (blue), and 0.70 (green) at decreasing pressures, and the opposite occurs for the minority component B (i.e., the chemical potential of the liquid crosses that of the gas phase at increasing pressures as x_A,v increases). These crossing points fulfill the condition that the chemical potential of the corresponding component is equal in the liquid and vapor phases. The crossing points for each component can be represented as a function of the vapor composition, x_A,v, obtaining two lines, one for each component, with opposite slopes [see Fig. 1(b)]. Only at the crossing point between these two lines does the equality of chemical potentials of each component occur at the same composition of the vapor (x_a,v) and at the same pressure. From this, we obtain that, at T* = 0.785, the symmetrical TSF-LJ liquid with composition x_A,l = 0.75 coexists with the vapor phase with composition x_A,v = 0.67(1) at $plv*=0.0341(3)$⁠. As can be seen in Table I, these results are in full agreement with those obtained by the direct coexistence route.

As mentioned earlier, one of the difficulties of studying nucleation by seeding is that this method depends critically on having a good criterion to estimate the size of the critical nucleus. Typically, this problem has been sorted out by comparing the nucleation rates obtained with seeding with those coming from more rigorous simulation techniques, such as umbrella sampling, forward-flux sampling, or brute force simulations.^9,10,27 Given that we are not aware of previous estimates of the bubble cavitation rate for the symmetrical TSF-LJ mixture, we started by studying cavitation by brute force simulations. At high metastabilities, the free energy barrier for nucleation decreases, and spontaneous cavitation becomes accessible within the typical length and time scales of simulations.

For that aim, we performed NpT simulations of a liquid with composition x_A,l = 0.75, containing 18 960 particles of type A and 6320 of type B. At these conditions, we were able to observe spontaneous cavitation only at negative pressures. For each considered pressure, we ran 10 independent simulations. Cavitation events can be easily detected by an abrupt increase in density. The cavitation rate can be estimated from these simulations as

where ⟨t⟩ is the average time for observing spontaneous cavitation, which is calculated simply as the average of the time for cavitation in the 10 independent simulations performed for each pressure, and $V$ is the average volume of the simulation box before the cavitation event. We checked that 10 simulations are enough to obtain an accurate estimate of the nucleation rate with an error significantly lower than one order of magnitude by comparing the rate obtained from 20 independent simulations with those obtained from only 10. The cavitation rates obtained by brute force simulations are given in Table II.

As seen in Sec. II B, the estimation of the cavitation rate requires to calculate the free energy barrier ΔG_c and the kinetic prefactor K. The cavitation free energy barrier can be obtained from the size of the critical cluster, R_c, and the pressure difference between the bubble and the surrounding liquid, Δp [Eq. (5)]. In the next subsections, we will see how to estimate these properties. It is important to recall that the expression of the free energy barrier ΔG_c, Eq. (5), was derived under the assumption that the composition of the metastable phase remains constant during the bubble nucleation. This is true in the thermodynamic limit, where the size of the critical bubble is negligible compared to the size of the system (so that the formation of the bubble causes only an infinitesimal change in the composition of the external liquid phase). In previous studies of the nucleation of droplets and crystals from binary mixtures, simulations were performed in the isothermal-isobaric semigrand ensemble (in which p, T, N, and the chemical potential difference between the two species are kept constant).^18,47,62 As will be shown later, in the liquid–bubble simulations performed in this work, the composition of the liquid remains almost unchanged. The reason is that, even though the composition of the bubble is different from that of the metastable liquid, the number of particles in the bubble is so small compared to that in the surrounding liquid that it barely changes its composition. Therefore, the use of the NVT ensemble, keeping the overall composition constant, is justified in this case.

The protocol used to construct initial configurations of the liquid in coexistence with bubbles is as follows: first, we equilibrate the liquid with a composition x_A,l = 0.75 by performing NpT simulations at T* = 0.785. In order to be able to accommodate bubbles of different sizes, we considered two system sizes: a smaller one with N = 10 000 particles and another one with N = 130 000 that equilibrated during 1 × 10⁶ MD steps. Bubbles of different sizes are then generated by removing all the particles from a spherical region defined in the equilibrated liquid simulation box. The sizes of the bubbles were chosen so that the ratio between the volumes of bubble and simulation box was within 5%–15% (see Table V). Previous work suggests that within these limits seeding simulations are not significantly affected by finite size effects.⁹ 

These systems containing the artificially constructed bubbles are then allowed to evolve in the NVT ensemble for 2 × 10⁶ MD steps. Initially, bubbles grow or shrink until an equilibrium with the surrounding liquid is established. Note that not all the bubbles will equilibrate, depending on the choice of the total number of particles, size of the simulation box, and total density.⁶³ In those systems in which the bubble remains stable after equilibration, i.e., it does not completely dissolve in the timescale of our simulation, we estimate their average radius in simulations consisting of 100–400 × 10⁶ MD steps.

As mentioned earlier, in seeding, nucleation rates can change significantly depending on the criterion used to define the cluster size.^9,29 The radius that enters a CNT is that given by the surface of tension, i.e., the surface on which the surface tension acts.^48,63–65 However, the location of that interface is not known.^66,67 For that reason, in previous implementations of the seeding method, the radius of the critical nucleus is estimated using a criterion that yields nucleation rates consistent with those obtained from rigorous simulations that do not depend on the critical cluster size, such as umbrella sampling or brute force simulations.^9,11,31

Here we employed two criteria used in previous studies to define the bubble size.^9,11 Later, we will compare the cavitation rates derived from these two alternative criteria to see which one gives results in closer agreement with those from brute force simulations. In the first criterion, the radius of the bubble is taken as the equidensity radius,⁹ which is defined as the distance at which the density is equal to the average of the liquid and gas densities measured sufficiently away from the interface (where radial densities remain constant). The equidensity radius is estimated from the radial overall density profile measured from the center of the bubble, as shown in Fig. 2(a). Following previous work,¹⁰ the center of the bubble in each instantaneous configuration is located by searching for the positions of the minima in the density profiles calculated along the x, y, and z directions using a bin width of about 0.7σ. The radial density profile is calculated using a bin of width 0.15σ in the radial direction and is typically averaged over 1000 independent configurations stored every 20 000 steps, so that the typical length of the simulations used for these calculations is 20 × 10⁶ MD steps. The radial overall density profile is then fitted to the following sigmoid function:

that has four fitting parameters: ρ_v,dp and ρ_l,dp provide estimates of the densities of the vapor and liquid phases (the subscript dp stands for density profile), α is a measure of the width of the interfacial region, and R_ed is the distance at which the density adopts the average liquid and vapor densities. The values for ρ_v,dp and ρ_l,dp obtained from this fit are coherent (within their statistical uncertainty) with the values inferred from averaging the density profiles of the vapor and the liquid regions, respectively, in Fig. 2. Notice that the error in the estimation of the vapor density ρ_v,dp is higher than that of the liquid density ρ_v,dp since there are only a few particles in the bubble, resulting in large fluctuations in the density.

Besides the radial density profile, we also calculated the radial concentration profile by dividing the density of particles of type A by the total density for each radial distance. A sample profile is shown in Fig. 2(b). The concentration of particles of type A is higher in the liquid phase than in the vapor phase, similarly to what was observed in the direct coexistence simulations. In this case, as the volume of the gas phase is relatively small, the composition of the liquid phase remains fairly close to the starting composition, namely, x_A,l ≃ 0.75. However, the molar fraction of A in the vapor phase is again slightly lower than that in the surrounding liquid (x_A,v ≈ 0.66). As explained earlier, this is due to the fact that the strength of the crossed interactions is 10% lower than that of like-interactions. In general, in all the NVT seeding simulations performed in this work, the concentrations of the bubbles are within x_a,v ≈ 0.670–0.673, and that of the surrounding liquid is typically within x_A,l ≈ 0.750–0.753.

Alternatively, one can also resort to the equimolar Gibbs dividing surface as a criterion to define the bubble radius.^9,66 The Gibbs dividing surface is obtained by assuming that the densities of the liquid and gas phases are constant up to the dividing surface with radius R_G. In multi-component systems, a Gibbs dividing surface can be defined for each component. Here we define the Gibbs radius using the majority component A,

where V_bubble and V are, respectively, the volume of the bubble and of the simulation box. Substituting the volume of the bubble by

and reorganizing we obtain that the Gibbs radius can be estimated as

where L is the edge of the simulation box. As mentioned earlier, we could have also defined the Gibbs radii using the minoritary component, B. These two Gibbs radii are only identical when the composition in the liquid and vapor phases are the same. As can be seen in Tables I and III, in our system, this is neither true at coexistence nor in states where nucleation is studied. As discussed earlier, as a consequence of the reduced crossed interactions, the minoritary component tends to go to the vapor phase, where interactions are less relevant. From Eq. (12), it is easy to see that this leads to R_G,B > R_G,A. However, as can be seen in Table IV, this effect is mild in our system. R_G,B is only about 2% bigger than R_G,A. Consequently, nucleation rates calculated with R_G,B are 0–6 orders of magnitude smaller than those calculated with R_G,A, but, as will be shown later, this difference is within the uncertainty of our calculations. Since it seems more sensible to use the majority component, we have chosen to use R_G,A in our calculations. In what follows, for simplicity in the notation, this radius will be denoted R_G.

As can be seen in Fig. 3(a) and Tables V and VI, the Gibbs radius is systematically larger than the equidensity radius. Taking into account Eqs. (2) and (5), this means that the Gibbs radius will lead to a higher cavitation free energy barrier and, consequently, to a lower cavitation rate than the equidensity radius. The same trend was also observed in previous seeding studies of cavitation^9–11 [see Fig. 3(a)] and droplet nucleation¹¹ in the pure TSF-LJ system. Gibbs radii are not provided in Refs. 10 and 11, but they can be estimated from the data reported in Ref. 11. The results for the pure TSF-LJ system using the Gibbs radii are also included in Figs. 3 and 5.

Following previous work on one-component systems,^10,11,28 we checked what happens to the stabilized bubbles in the NVT ensemble when simulated in the NpT ensemble. For that aim, we launched 36 simulations in the NpT ensemble starting from bubble–liquid configurations taken along the trajectory of the system with index 9 in Table III at the corresponding $pl*$ $(pl*=0.0031)$⁠. The outcome of those simulations is shown in Fig. 4, where the evolution of the total density of the system in all the NpT simulations is monitored. As can be seen, in about half of the trajectories, the bubble dissolves, and in the other half, the bubble grows irreversibly to the vapor phase, indicating that the cluster is indeed critical.

One might intuitively expect that by virtue of the Gibbs phase rule, for binary systems, it would be possible to have a critical cluster in stable equilibrium at constant p and T. Indeed, this has been exploited to not only have a planar interface of a binary system at equilibrium in the NpT ensemble^68,69 but also to stabilize spherical crystal seeds of NaCl⁶⁸ in coexistence with a salt solution and of methane hydrate in coexistence with water with some methane dissolved.⁷⁰ However, as we have just seen, the critical bubbles in the binary symmetrical TSF-LJ mixture cannot be stabilized in the NpT ensemble. Therefore, why is it possible to stabilize the critical seed in the NpT ensemble for some multi-component systems but not for others? Further work is needed to clarify this issue, but it is important to highlight that in all cases in which stable clusters in the NpT ensemble were achieved, they had a constant composition, given by the stoichiometry of the crystals (the NaCl rock salt and the methane hydrate), whereas in our case, the composition of the bubbles is not fixed. Small fluctuations in the composition of the clusters change the free energy barrier for bubble formation, making the cluster precritical or postcritical at the specified p and T and, therefore, leading to their dissolution or growth. As for the one-component systems, we found that the bubbles only remain stable in simulations at constant volume.

Once we have estimated the bubble size, we proceed to calculate the difference in pressure between the bubble and the surrounding liquid, $Δp*=pbubble*−pl*$⁠, which is the other quantity needed to estimate ΔG_c [Eq. (5)]. In this work, we followed two different routes to achieve this aim. In both approaches, the pressure of the liquid that surrounds the bubble, $pl*$⁠, was taken as the virial pressure in the NVT simulations, as provided by GROMACS.⁴³ Note that as we are dealing with a heterogeneous system in which the bubble is in coexistence with the liquid, the pressure is not uniform in the system. However, the average of the pressure tensor in the whole system coincides with that of the external phase, i.e., with that of the liquid.⁷¹ The pressure of the liquid in all the seeding simulations is reported in Table III. The difference between the two methods relies on how to estimate the pressure inside the bubble, $pbubble*$⁠.

In the first approach, the pressure inside the bubble is taken as the mechanical pressure, which can be estimated from the density and composition inside the bubble together with an EOS for the bulk vapor phase obtained from simulations $[pbubble*=p*(ρv*,xA,v)]$⁠. We checked that the EOS of the bulk vapor shows a mild dependence on composition within the range x_A,v = 0.670–0.673. This is rather convenient, because even though the composition of the bubbles is not exactly the same in all the seeding simulations, it remains within these boundaries, and the use of the same bulk EOS for all the bubbles is justified. Here, the EOS was calculated by performing NpT MD simulations of a vapor with composition x_A,v = 0.672 and sweeping pressures within p* = 0.002–0.033.

In the second approach, we use the thermodynamic route to estimate $pbubble*$⁠. The equilibrium established between the bubble and the surrounding liquid implies that the chemical potential of each component is the same in the two phases. The pressure inside the bubble is the one that fulfills this condition, which can be obtained from the chemical potentials calculated in Sec. III A 2 and shown in Fig. 1(a). The procedure is as follows: first, we estimate the chemical potentials of components A and B in a liquid with x_A,l = 0.75 at pressure $pl*$⁠, namely, $μA,l(pl*,xA,l)$ and $μB,l(pl*,xA,l)$⁠. The gray vertical line in Fig. 1(a) marks the pressure of the surrounding liquid, $pl*$⁠, for the system with index 10 (Table III), and the two horizontal gray lines mark the chemical potentials of components A and B in such a liquid. Second, we look for the crossing points of these horizontal lines with the chemical potential lines of the corresponding component in the vapor phase at several vapor compositions (x_A,v = 0.60, 0.65, and 0.70). At these crossing points, the chemical potential of each component is the same in the liquid and in the vapor at concentration x_A,v, i.e., they fulfill $μA,l(pl*,xA,l)=μA,v(pv*,xA,v)$ (purple symbols) and $μB,l(pl*,xA,l)=μB,v(pv*,xA,v)$ (red symbols). These crossing points are represented and fitted to straight lines in Fig. 1(c). Only at the crossing point between these two lines does the equality of chemical potentials of the two components in both phases occur at the same composition of the vapor phase, x_A,v, and at the same vapor pressure $pv*$⁠, giving the equilibrium pressure inside the bubble $pbubble*$ and its composition. For the system labeled with index 10 and shown in Fig. 1(c), the pressure and composition of the bubble in equilibrium with the liquid are $pbubble*=0.0316$ and x_A,v = 0.671, which compare well with the values $pbubble*=0.0309$ and x_A,v = 0.656 obtained with the first described method, i.e., using estimates of the density and composition inside the bubble together with the simulated EOS of the bulk vapor.

The results of Δp* obtained from the mechanical and thermodynamic routes are plotted jointly in Fig. 3(b) for comparison. As can be seen, both methods give results that are consistent with each other, but the chemical potential route provides less noisy data than the EOS route. Consistent with previous findings in the pure system, this means that the mechanical pressure is equivalent to the thermodynamic pressure.^11,72 Note that this does not always hold, as has been shown in a recent study of hard-sphere crystal nucleation.²⁸ The higher statistical uncertainty in the mechanical pressure is expected, given the difficulty of accurately estimating the density inside the bubble. The number of vapor particles is very small, and fluctuations in the number of particles lead to relatively large changes in the density. Besides, only the results obtained from the thermodynamic route guarantee that the chemical potentials are the same in the bubble and the surrounding liquid. For those reasons, the values of Δp* obtained from the chemical potential route are the ones that will be used in the remainder of the article (these data are tabulated in Table III). Curiously, as can be seen in Fig. 3(b), the difference in pressure between the bubble and surrounding liquid is very similar in the pure and the symmetrical binary TSF-LJ mixture when compared at the same stretching $pl*−plv*$⁠, despite the bubbles being much bigger in the pure system.

The density of the metastable liquid, $ρl*$⁠, needed to calculate the kinetic factor K, Eq. (7), is shown in Fig. 3(b). Again, we tested two possible ways to estimate this density. In the first one, we used the EOS of the bulk liquid with x_A,l = 0.75 to obtain the density at the corresponding liquid pressure $pl*$⁠, obtained from the virial pressure in seeding simulations. In the second approach, the density is directly measured from the seeding simulations using the sigmoid fit. As can be seen in Fig. 3(c), the agreement between both sets of data is fairly good, except for the point at the lower pressure and some of the points at the highest pressures. It is also worth pointing out that the seeding points that exhibit the largest departures from the EOS are precisely those with a higher V_bubble/V_box ratio (see Table V), which might be indicative that these small deviations arise from finite size effects. In any case, differences are always lower than 0.3%, which can be taken as an indication that finite size effects are not very significant in any of our seeding simulations.

Due to the reduced cross interactions in the mixture, the liquid density is lower than in the pure system. Therefore, the kinetic factor K is about 15% lower in the mixture¹¹ both because of the smaller bubble size and the reduced density of the liquid as compared to those of the pure system at the same stretching.

Now that we have estimated the radius of the critical cluster (R_c), the difference of pressure inside and outside the bubble (Δp*), and the density of the metastable liquid, we can proceed to calculate the free energy barrier for cavitation using Eq. (5). The results as a function of stretching relative to the coexistence $(pl*−plv*)$ are shown in Fig. 5(a). As the Gibbs radii are larger than the equidensity radii, the free energy barrier calculated with the Gibbs dividing surface is larger, with the difference increasing with stretching as the critical bubble gets smaller. However, both sets of data support the conclusion that the free energy barrier is significantly lower in the mixture than in the pure TSF-LJ system when compared at the same stretching relative to coexistence. Interestingly, a similar trend was observed in a simulation study on the effect of composition on the nucleation of droplets from a LJ binary mixture.¹⁸ 

Next, we present the cavitation rates obtained from Eq. (2) using the data summarized in Tables V and VI. As can be seen in Fig. 5(b), the small differences in the free energy barrier obtained with the two defined radii are also visible in the nucleation rates. Rates are lower using the Gibbs criteria. Both sets of data give consistent results with brute force simulations, although those derived from the Gibbs radii show better agreement. However, this should be taken with care, because for the pure TSF-LJ system, R_G also coincides better with brute-force simulations,⁹ but at lower metastabilities R_ed provides rates consistent with umbrella sampling.¹¹ 

The interfacial free energy can be obtained from the Laplace equation [Eq. (4)]. The values calculated from the equidensity and Gibbs criteria are shown as a function of the inverse radius and as a function of $pl*−plv*$ in Figs. 5(c) and 5(d), respectively. As can be seen, the extrapolation of both sets of data to an infinite radius gives coherent results with our calculation of the interfacial free energy at coexistence, but the values of γ obtained from R_G exhibit a considerably smaller slope than those derived from R_ed, in this latter case predicting that the cost of creating a bubble–liquid interface is rather insensitive to the curvature of the bubble, adopting values close to those of the planar interface from very small sizes.

The values of γ are compared with those of the pure TSF-LJ system¹¹ in Fig. 5(d). For the pure TSF-LJ system, γ was estimated using R_ed,¹¹ because it gave rates coherent with those from umbrella sampling.⁹ Consistent with the results in the thermodynamic limit, the interfacial free energy is systematically larger for the pure system than for the symmetrical mixture, with a slope intermediate between those derived from R_ed and R_G in the mixture. The lower interfacial free energy of the mixture as compared to that of the pure system explains why cavitation rates are significantly larger in the mixture.

The cavitation rate was fitted using CNT expressions as described in Ref. 10. ΔG_c is a function of R_G and Δp [Eq. (5)], and the kinetic factor K depends on R_G, Δp, and ρ_l. As can be seen in Fig. 3, the dependence of Δp and ρ_l with p_l − p_lv follows a smooth behavior, but that of R_G is more complex. Therefore, it is convenient to rewrite ΔG and A in terms of γ, which shows a nearly linear behavior with p_l − p_lv [Fig. 5(d)]. Using the fits of Δp, γ, and ρ_l as a function of p_l − p_lv provided in the captions to Figs. 3 and 5, we obtain J as a function of p_l − p_lv, which is plotted as a dashed red line in Fig. 5(b). The error in the rate comes from the uncertainty in the estimation of R_G, Δp and ρ_l, with the larger error coming from the critical radius. Assuming that the error in the critical radius is of about 0.3σ and neglecting the uncertainty in the other quantities, we can estimate the error in the rate, which is delimited by the red dotted lines in Fig. 5(b). As can be seen, this fit and its error capture well the variation of the rates with the stretching. If instead of using the fit of γ we use the capillarity approximation (i.e., we assume that adopts the same value as the planar interface), the rate is very similar to that obtained using the rates calculated from R_G, which is expected as we have seen that γ was rather insensitive to the bubble curvature. Differences are much larger for the data derived from R_ed because γ exhibits a larger dependence on curvature.

A theoretical work⁷³ has shown that in the pure LJ system, the nucleation free energy barrier for the formation of bubbles is the same for all temperatures if compared at the same scaled metastability, defined as

where μ_nucl, μ_coex, and μ_spin are the chemical potentials of the liquid at nucleation, coexistence, and spinodal. With this definition, X_m varies between 0 at coexistence and 1 at the spinodal. Interestingly, in subsequent work,¹¹ it was shown that the free energy barrier for droplet and bubble formation as a function of X_m also collapses into a single curve. In that regard, it seems that the free energy barrier for nucleation is the same for bubbles and droplets if compared at the same degree of scaled metastability.

In view of this, it is interesting to compare the free energy barrier for the LJ mixture as a function of X_m with that of the pure system. For a two component system, the driving force for nucleation depends on the chemical potential of the two components. Therefore, we define

and this is the quantity that will be evaluated at coexistence, at the nucleation conditions, and at the mechanical spinodal to estimate X_m. For a pure component system, the spinodal is unique and is defined by the condition (∂p/∂ρ)_T = 0. For a binary mixture, one can define two spinodals, the mechanical one given by (∂p/∂ρ)_T,x = 0 and the diffusional one defined as (∂μ/∂x)_p,T = 0.^74,75 The location of these two spinodals can be quite different.⁷⁴ Since for pure components we found a correlation between the free energy barrier and X_m via the mechanical definition of the spinodal, we shall explore the possibility of a similar correlation for mixtures when also using the mechanical definition of the spinodal. We estimated its location using an approximate method, as described in Ref.11. In particular, we performed brute force simulations as described in Sec. III B at decreasing pressures. The spinodal is located at the highest pressure at which cavitation occurs immediately at the beginning of the simulations, without an induction period. We estimate that the spinodal is located at p* ≈ −0.05. The free energy barriers for bubble formation obtained either with R_ed or R_G are plotted as a function of X_m in Fig. 6, together with the results for the pure TSF-LJ system taken from previous work.¹¹ As can be seen, the free energy barrier for the pure and the binary systems adopt quite similar values when compared at the same scaled metastability. This means that the cost of creating a critical bubble does not depend much on the interactions (at least for symmetrical mixtures with low deviations of the reduced crossed interactions from unity as the ones considered), and it is mostly determined by the scaled metastability.

In this work, we have applied NVT seeding to calculate bubble nucleation rates in a partially miscible binary symmetrical TSF-LJ mixture. The study was conducted at a temperature that falls within the mixing region of the system, for which there are previously available data for the pure TSF-LJ model. One aspect that requires special care when implementing seeding is the choice of an appropriate criterion to determine the bubble size. Here, we used both the equidensity radius and the equimolar (or Gibbs dividing surface) radius. Consistently with previous work,^9,11 the equidensity radii are systematically smaller than the Gibbs radii, leading to a lower free energy barrier and a higher nucleation rate. However, what is the correct value of the critical bubble size? As mentioned before, the correct radius is the one corresponding to the surface of tension, and it is the one that enters into CNT expressions. Therefore, we compared the nucleation rates derived from the equidensity and the Gibbs surface radii with those obtained from brute force simulations. Gibbs radii give rates closer to those obtained by brute force simulations. The same trend was found for bubble nucleation in pure TSF-LJ, but at lower metastabilities, rates calculated with the equidensity radii are closer to those obtained by forward-flux sampling.⁹ For droplet nucleation in the same system, equidensity nucleation rates are also in closer agreement with umbrella sampling results.¹¹ However, a recent study using enhancing sampling methods and transition state theory found results that are in agreement with seeding but obtained from Gibbs radii and concluded that forward-flux sampling gives systematically higher rates.^6,7,12 Therefore, we can conclude that for the LJ model at high metastabilities, close to the region where brute force simulations are feasible, the Gibbs dividing surface radii give the correct nucleation rates. However, at low metastabilities, where special simulation techniques are needed, the situation is not yet clear, as there are conflicting results using rigorous techniques.¹² Further work is needed to clarify why different methods provide different rates at low metastabilities.

In any case, the differences in the cavitation behavior between the pure and symmetrical mixtures are larger than the differences between seeding rates using different criteria for determining the bubble size. Therefore, conclusions about the difference in the cavitation behavior in the two systems can be drawn. At the same stretching, measured as the pressure difference for coexistence, the critical bubble is smaller in the binary system than in the pure system. The reason is that the interfacial free energy is about 25% lower in the mixture than in the pure system, i.e., the free energy cost of creating a bubble is lower in the mixture. This is, at least in part, related to the fact that the critical temperature is lower in the mixture than in the pure system.^40,59 Curiously, when comparing both systems at the same degree of scaled metastability, X = (μ_liq − μ_coex)/(μ_spin − μ_coex), the free energy barrier is similar in both systems. This is an interesting result, as it was found in previous work that free energy barriers at different temperatures for bubble and droplet nucleation merged into a single curve when represented as a function of the scaled metastability.^11,73 Our results suggest that this conclusion can also be extended to mixtures, at least for small deviations of the crossed interactions from those between like interactions, as the one considered here.

In the future work, it would be interesting to extend the NVT method to other systems and transitions, for example, to study droplet formation in strongly non ideal mixtures or crystallization from a binary mixture. In principle, NVT should also be applicable in this case; however, some difficulties might arise. For example, if important fractionation effects occur during nucleation (i.e., the composition of the nucleated and the parent phase are too different), if one wants to keep the composition constant, large simulation boxes are needed or that other ensembles, such as the semigrand canonical ensemble, need to be employed. In NVT seeding, the critical cluster is in equilibrium with the metastable phase, so that the system evolves until the chemical potentials of the two components are the same in the two phases. But if one of the phases is solid, diffusion might be too slow to reach the equilibrium.

It would also be interesting to explore alternative routes to calculate the free energy cost of creating a bubble in the metastable liquid. A possibly viable method is to introduce an external potential that creates a bubble in the liquid and design a reversible path that allows the calculation of the free energy by thermodynamic integration.

This project has been funded by Grant Nos. PID2019-105898GB-C21 and PID2020-115722GB-C21 from the Agencia Estatal de Investigación. C. P. L. thanks the Ministerio de Educación y Formación Profesional (MEFP) for a predoctoral Formación Profesorado Universitario under Grant No. FPU18/03326 and also the Ayuntamiento de Madrid for a Residencia de Estudiantes grant. The authors acknowledge the computer resources and technical assistance provided by the Red Española de Supercomputación (RES).

The authors have no conflicts to disclose.

Cintia P. Lamas: Data curation (equal); Formal analysis (equal); Methodology (equal); Validation (equal); Visualization (equal); Writing – original draft (supporting); Writing – review & editing (equal). Eduardo Sanz: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (supporting); Writing – review & editing (equal). Carlos Vega: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (supporting); Writing – review & editing (equal). Eva G. Noya: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (lead); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.