Study of homogeneous bubble nucleation in liquid carbon dioxide by a hybrid approach combining molecular dynamics simulation and density gradient theory

A new method for predicting homogeneous bubble nucleation rates of pure compounds from vapor-liquid equilibrium (VLE) data is presented. It combines molecular dynamics simulation on the one side with density gradient theory using an equation of state (EOS) on the other. The new method is applied here to predict bubble nucleation rates in metastable liquid carbon dioxide (CO 2 ). The molecular model of CO 2 is taken from previous work of our group. PC-SAFT is used as EOS. The consistency between the molecular model and the EOS is achieved by adjusting the PC-SAFT parameters to VLE data obtained from the molecular model. The influence parameter of density gradient theory is fitted to the surface tension of the molecular model. Massively parallel molecular dynamics simulations are performed close to the spinodal to compute bubble nucleation rates. From these simulations the kinetic prefactor of the hybrid nucleation theory is estimated, whereas the nucleation barrier is calculated from density gradient theory. This enables the extrapolation of molecular simulation data to the whole metastable range including technically relevant densities. The results are tested against available experimental data and found to be in good agreement. The new method does not suffer from typical deficiencies of classical nucleation theory concerning the thermodynamic barrier at the spinodal and the bubble size dependence of surface tension, which is typically neglected in classical nucleation theory. In addition, the density in the center of critical bubbles and their surface tension is determined as a function of their radius. The usual linear Tolman correction to the capillarity approximation is found to be invalid.


I. INTRODUCTION
Bubble nucleation is important in many applications. It occurs for example during ultrasonic cleaning or upon cavitation in technical equipment. In the present work, homogeneous bubble nucleation in metastable liquids is studied. That process is not spontaneous, but inhibited by a free energy barrier. In metastable liquids small bubbles stochastically form and break down. If a bubble exceeds the critical size, where the free energy of formation is maximal, it is likely to grow and become a macroscopic bubble. In this paper, a hybrid approach for describing homogeneous bubble nucleation is developed.
* Author to whom correspondence should be addressed. Electronic mail: kai.langenbach@mv.uni-kl.de As homogeneous bubble nucleation requires activation energy, the rate of forming bubbles J can formally be written in the form of an Arrhenius type expression * 0 exp JJ kT       (1) where 0 J is the so-called kinetic prefactor, prefactor 0 J is often simply adopted from some version of CNT e.g. [20][21][22] . This is dangerous, as CNT makes assumptions on the exponential term, which may be inconsistent with the DGT/DFT results.
A different approach is proposed here: DGT is used for determining the thermodynamic factor * exp( / (kT))  , while the kinetic prefactor 0 J is determined from molecular dynamics simulations instead of the usual CNT expressions, thereby hybridizing both theories. Presently studies of nucleation by molecular simulation can only be carried out for states close to the spinodal, where nucleation events occur frequently. Hence, the hybrid nucleation theory (HyNT) developed here can be seen as an extrapolation of nucleation rates determined by molecular dynamics close to the spinodal to state points in the entire metastable range by DGT.
As an example, the bubble nucleation in metastable liquid CO 2 is studied. The molecular model of CO 2 is taken from Merker et al. 25 . It is known to give good results for vapor pressure, liquid and vapor coexistence densities, second virial coefficients as well as transport properties 25 . The parameters of the PC-SAFT EOS are fitted to vapor-liquid equilibrium (VLE) data of Merker et al. 25 in order to ensure consistency between both models. Molecular simulations are carried out for homogenous metastable states and compared to results from PC-SAFT to verify that the PC-SAFT model is valid in the range of interest for nucleation. The influence parameter of DGT is obtained from a fit to surface tension data obtained with the molecular model of Merker et al. 25 taken from 26 . It should be noted that the experimental surface tension of CO 2 is overpredicted by roughly 20% 26 by the model of Merker et al. 25 . That deviation is accepted here, as no molecular model of CO 2 is presently available that would yield better results for the entire set of thermodynamic properties which are relevant for the present application. Finally, direct molecular dynamics simulations of bubble nucleation are carried out for two temperatures and several bulk phase densities close to the spinodal. The bubble formation statistics are evaluated with the approach introduced by Yasuoka and Matsumoto 27 for droplet nucleation and later on applied by Diemand et al. 28 for bubble nucleation. Using molecular simulation, no assumptions are made regarding the separation of kinetic prefactor and thermodynamic factor. However, in the range, where direct molecular simulation of nucleation is possible, the kinetic prefactor is dominant. This dominance is stronger the closer simulations are performed to the spinodal. Therefore, the prefactor of HyNT is calculated from the simulated nucleation rate closest to the spinodal. The other simulated nucleation rates are used to test the accuracy of the extrapolation using HyNT. The results of HyNT for homogeneous bubble nucleation are further compared to experimental data for the bubble nucleation of CO 2 29,30 .

A. Hybrid Nucleation Theory
In hybrid nucleation theory (HyNT), with the nucleation rate formally written as in Eq. (1), the free energy *  of forming a critical bubble nucleus in a surrounding mother phase of density h  and temperature T is calculated from DGT.
DGT 23,31,32 can be used in conjunction with a suitable EOS to investigate systems with a heterogeneous density distribution. In this work the PC-SAFT EOS 24 is used, which has been shown before to yield good results for both bulk fluid phase behavior and interfacial properties when combined with DGT e.g. 33 . For a heterogeneous one-component system, the grand potential is written as a functional of the spatially varying density   x  and a function of the temperature T according is the grand potential per volume. Here, 0 f is the homogeneous free energy density calculated from the PC-SAFT EOS at the local density () x  ,  is the chemical potential calculated for the density of the mother phase h  , p is the pressure, and  is the so-called influence parameter. Since we are interested here in the grand potential of critical bubbles, which corresponds to the free energy of forming them in the mother phase of the molar density h  , radial symmetry is assumed. In spherical coordinates, the grand potential is given by where r is the radial coordinate. Functional minimization leads to the Euler-Lagrange equation 2 (9) The influence parameter  is obtained from a fit to data for the planar surface tension 35 using the expression 31 where '  and ''  are the molar liquid and gas density in equilibrium respectively.
The kinetic prefactor in Eq. (1) is defined here as where k 0 is the kinetic part determined from molecular simulations close to the spinodal as elaborated in paragraph IV and *  is the unit cavity density defined here as (12) where av N is Avogadro's number.

B. Classical nucleation theory
and Döring 38 , and Zeldovich 39 . CNT has been modified empirically and semi-empirically in many ways, cf.  (13) where m is the mass of the CO 2 molecule. Since the surface tension of the planar interface   is used, this is also called the capillarity approximation. We use PC-SAFT to connect the pressure, temperature and density.

A. Molecular simulation
The simulations were carried out using the ls1 mardyn 40 program, which was modified by the authors to evaluate bubble nucleation statistics. The CO 2 model is the rigid three Lennard-Jones sites + point quadrupole model of Merker et al. 25 . It is given in Appendix B for completeness. In all simulations a cutoff at r c = 13.2 Å was used. The long range correction for the Lennard-Jones force field was treated with the method proposed by Lustig 41 , while the quadrupole does not need such correction since it vanishes by angular integration. In the method of Lustig 41 , a center of mass cutoff is employed adding constant mean field correction terms to the potential energy and the virial, where contributions from individual interaction sites are determined by angular integration. This is a simplification for bubbles, however, as long as the bubbles are small it should be a good approximation. The leapfrog integrator 42 was used with a time step of 2.5 fs. Simulations were conducted in the NVT ensemble. The occurrence of gaseous regions was sampled on a regular grid. The number of grid points M was of the order of N/4 in all cases, see Appendix C for details. A point was considered to be gaseous, if N g,max = 5 or less molecular centers of mass were in a sphere of the radius r b = 6.88 Å around the point. This local density criterion corresponds to a threshold density of roughly 6.1 mol l -1 . By variation of both the number of molecular centers N g,max and the radius r b , it was established that the results are not particularly sensitive to the particular choices, as long as the corresponding density is gaseous, i.e. below the critical density. This is the case for the present choice of parameters. A cluster c is an aggregation of unit cavities. The cluster size is determined by connecting gaseous grid points that are immediate neighbors and counting the number M c of connected grid points per cluster. The volume of these clusters corresponding to bubble volume is calculated via where V is the simulation box volume. The number of such bubbles exceeding a certain threshold volume V min was counted and used to determine the nucleation rate as explained below.
A crystal like structure, with superimposed random displacements of the individual molecules, was used as the initial configuration. No pre-equilibration was conducted. In this way two series of simulations close to the liquid spinodal were carried out in which for T = 220 K and T = 280 K several densities were studied (cf. Table 4 and Table 5 in Appendix C).
In each simulation the bubble nucleation rates are determined using the Yasuoka-Matsumoto 27 Tables 4 and Table 5 in Appendix C. Additional information on the simulations is supplied in Table 6 of Appendix C. or smaller, so that the influence of the finite threshold volume, which was at least 1 nm 3 in all cases, can be neglected and FSE 1 is absent. For T = 220 K and ρ h = 23.2 mol l -1 , the nucleation rate is comparably large, so that at large values of V min , a decay of the rate of formation is found due to the influence of the finite simulation volume. Hence, FSE 2 is present. In this case, the macroscopic nucleation rate corresponds to the rate of formation for small bubbles. It is labeled as case A (FSE 1: no; FSE 2: yes). For T = 220 K and ρ h = 23.9 mol l -1 also FSE 2 is absent, since J is small, and no finite-size effects are observed. A flat diagram is obtained, and any of the values can be used directly to estimate the macroscopic nucleation rate. For obtaining a nucleation rate from a given simulation not only a single threshold value V min is used. Rather, results obtained using several values of V min from the region in which the influence of both FSE 1 & FSE 2 is small are used. The highest and lowest of these values for one simulation are considered the confidence interval and their geometric mean is the J reported in Table 1. All values for J are given in Table 4 where those inside the confidence interval are underlined. comparably small systems. Particle numbers were between 800 and 1500. These simulations were conducted to capture the relation between density and pressure for the employed CO 2 model at T = 220 K and T = 280 K. Both stable and metastable states were studied. The same local density criterion as for the large scale simulations of bubble nucleation was applied to a regular grid with size M = 10, where only simulations were taken into account with more 999.95 grid points on average belonging to a single metastable or stable phase. The simulation results are listed in Table 7.

B. PC-SAFT + DGT
CO 2 is modeled here with the original PC-SAFT 24 without any polar terms [43][44][45] . Hence, the PC-SAFT CO 2 model has three parameters: the number of segments m, the segment diameter σ, and the interaction energy ε. They were fitted to the molecular simulation data of Merker 25 for the vapor pressure, saturated liquid density, and saturated vapor volume, who used  Table 2, a comparison between the molecular simulation data for   6 shows results for the thermal behavior of CO 2 obtained from the molecular simulations and from PC-SAFT.
Besides data for the vapor-liquid equilibrium, which agree very well, data for isotherms at 220 and 280 K is shown. That data covers both the stable, metastable and unstable region. Molecular simulation data is available for stable and metastable states.
It agrees very well with the PC SAFT predictions.

A. Evaluation of kinetic prefactor and comparison of nucleation rates from theory and simulation
The prefactor k 0 of Eq. (11) as used in Eq. (1) is fitted such that to the nucleation rate of Table 1 for the lowest density above the spinodal density (obtained from PC-SAFT) for the respective temperature and that from the HyNT nucleation rate agree. The kinetic prefactor is temperature dependent. We choose a simple Arrhenius type dependence with coefficients a = -18.36 and b = 8.611 to describe that dependence. This is somewhat arbitrary, but it has been argued in literature 14,47 that the kinetic prefactor depends on several quantities (e.g. viscosity, vapor pressure) that themselves have an exponential temperature dependence, hence our choice. Another sensible option would have been a power law dependence, however, with only two measured data points for k 0 both choices seemed likewise plausible. In FIG. 8 bubble nucleation rates for T = 220 K and T = 280 K from molecular simulation, CNT according to Eq. (13) and HyNT are shown.
The molecular simulation results show that the nucleation rate close to the spinodal is larger for lower temperatures. Since HyNT is fitted to that nucleation rate it shows the same temperature dependence. CNT predicts the inverse trend. With increasing density, nucleation rates drop super-exponentially. For T = 280 K, HyNT shows a very good agreement with simulation. CNT overestimates the simulation data by roughly two orders of magnitude at that temperature, even though the same EOS and surface tension data is used. For T = 220 K the agreement between HyNT and simulation is still reasonable, however, not as good as for the higher temperature. This may be due to the simplifications in DGT, namely the assumptions that the influence parameter κ is a constant and the truncation of the series expansion (cf. Eq.

B. Comparison to experimental data
For CO 2 only limited experimental data on metastable states and nucleation is available. Straub 29 has measured the limit of supercooling for isochoric cooling of liquid CO 2 below the coexistence temperature. This was also done by Huang et al. 30 .
In both cases measurements were made close to the critical temperature of CO 2 . In order to compare results from such experimental data to any nucleation theory both the sample volume and the relaxation time, i.e. the average time until the first macroscopic bubble is observed, need to be known. Neither Straub 29 nor Huang et al. 30

C. Discussion
The HyNT method presented above relies both on molecular simulation and an EOS + DGT. The employed molecular model and the EOS should be consistent and they should describe the available experimental data on VLE and surface tension of the studied fluid well. In CNT an EOS describing the VLE and surface tension data is needed so that both HyNT and CNT require the same input. As shown above, CNT predicts the wrong temperature dependence of the nucleation rate J, whereas in HyNT that problem does not occur. Nucleation rate is calculated from two different factors, the kinetic prefactor J . The former is calculated using DGT + PC-SAFT and the latter is fitted to simulated nucleation rates for two temperatures, as close to the spinodal as possible, because there the thermodynamic barrier vanishes. HyNT shows a good agreement with nucleation rates determined from simulation, also for simulations that were not included in the fit. CNT in the version of Blander and Katz 4 , using the same EOS, shows the wrong temperature dependence of the nucleation rate. This is explained in terms of the non-zero limiting nucleation barrier height of CNT at the spinodal. Hence, CNT should not be used to interpret data from molecular simulation close to the spinodal.
Further improvements of HyNT might be possible using a density functional instead of the DGT used here, however, at the expense of increased complexity. Results for the isochoric limit of supercooling for liquid CO 2 obtained from HyNT are compared to experimental data 29,30 and show good agreement. However, also CNT agrees well with that data, which is only available in a small temperature range close to the critical point and therefore not suitable for discriminating the results of both theories. The new theory should be tested for temperatures closer to the triple point of CO 2 preferably using direct bubble nucleation rate measurements, if such data becomes available. Additionally, the linear Tolman correction to the surface tension is found to be qualitatively wrong for very small bubbles, i.e. in the regime, where molecular simulation takes place. The qualitative dependency of the density in the center of critical bubbles on the bubble size that was observed earlier for the Lennard-Jones truncated and shifted fluid is confirmed here for CO 2 .

APPENDIX A: Numerical Procedure for Determining Radial Density Profiles
Density profiles are calculated according to the differential Eq. (4). This equation is discretized using finite differences with a step width of h = 0.2 Å using a central difference scheme. Rewriting Eq. (4) for point i in the finite difference scheme results in for i > 0 and in    4) is exactly zero. This means that at vanishing derivative at some large r, the homogenous density h  must be recovered for the non-trivial solution of this boundary value problem. Using the iteration formulas above, we guess a starting density 0  and calculate the profile until the absolute change of density from i to i+1 divided by the total density difference between center of the bubble and homogeneous phase is smaller than 10 -6 . If the density at that point n is larger than the homogeneous phase, 0  is reduced and if larger it is increased and the density profile recalculated. After that initial step, the bisection method is used to find the Using this method, the calculation of profiles is limited by the machine precision. If the lower and upper bound in the bisection method differ by less than double the machine precision, no better solution is possible. This occurs approaching the binodal, where the radius of the bubble grows to infinity. There have been suggestions in literature 50 to work around this, however, we do not apply such methods here, since we are mainly interested in the range closer to the spinodal.

APPENDIX B: Molecular Model for CO 2
The molecular model for CO 2 by Merker et al. 25 consists of three linearly arranged Lennard-Jones sites, one for each oxygen and the carbon and a point quadrupole on the carbon aligned with the molecular axis. The parameters of the model are the Lennard-Jones diameters of the oxygen σ O and carbon σ C as well as their interaction energies ε O and ε C The distance between carbon and oxygen r CO and the quadrupole moment Q, all of which are listed in Table 3. In Table 4 and Table 5 the nucleation rate from molecular simulation, evaluated as described in paragraph III.A, is given for the two investigated temperatures as indicated in the heading for different densities close to the spinodal. The numbers of particles used for the different simulations as well as the numbers of grid points for bubble evaluation in one dimension are given in Table 6. Simulation results for metastable and stable p-ρ-T behavior of the CO 2 model are given in Table 7.  17.575 9.7 • 10 33 5.9 • 10 32 5.0 • 10 31 n/a n/a n/a n/a n/a 17.650 7.8 • 10 33 3.6 • 10 32 n/a n/a n/a n/a n/a n/a  Table 7. Thermal properties of homogeneous states of CO 2 as described by the model of Merker et al. 25 : Pressure p as a function of temperature T and density ρ. The number in parenthesis is the uncertainty of the last digit of the number reported for the pressure. Due to data loss, no uncertainty for the simulations at T = 280 K can be reported here.
ρ / mol l -1 p / MPa    Both data for the vapor-liquid equilibrium and data for isotherms at 220 K and 280 K in the stable metastable and unstable range are shown. Full line: binodal PC-SAFT; dashed line: isotherm PC-SAFT; squares: VLE molecular simulation Merker et al. 25 ; symbols are data for the isotherms obtained from molecular simulation (diamonds: 220 K; circles: 280 K). Panel a) shows the gas side while panel b) shows the entire studied range.