Theoretical prediction of the homogeneous ice nucleation rate: disentangling thermodynamics and kinetics

TL;DR

Addressing the long-standing mismatch in homogeneous ice nucleation rates, the paper develops a framework to separate thermodynamic and kinetic factors in nucleation from undercooled water. It computes the nucleation free-energy profile for a perfect ice Ih nucleus using umbrella sampling with a global CV and a Gibbs dividing surface, then adds a stacking-disorder correction based on a stacking-fault free energy $\gamma_{\rm sf}(T)$ to capture Ic/Ih competition. The kinetic prefactor is obtained from a stochastic Langevin model fitted to umbrella-sampling data, enabling estimation of the rate $J$ via $J = (1/v_{\text{l}}) Z f^+ \exp(-G^{\star}/(k_B T))$. Applying the method to the mW water model at $T=240$ K and $P=1$ bar yields $J \approx 0.3~\mathrm{m^{-3}s^{-1}}$, nine orders of magnitude larger than prior seeding-based results, highlighting the importance of curvature corrections and stacking disorder. The framework provides a path toward ab initio, cross-validated predictions of homogeneous nucleation rates and can be extended to other molecular crystals and complex nucleation pathways.

Abstract

Estimating the homogeneous ice nucleation rate from undercooled liquid water is at the same time crucial for understanding many important physical phenomena and technological applications, and challenging for both experiments and theory. From a theoretical point of view, difficulties arise due to the long time scales required, as well as the numerous nucleation pathways involved to form ice nuclei with different stacking disorders. We computed the homogeneous ice nucleation rate at a physically relevant undercooling for a single-site water model, taking into account the diffuse nature of ice-water interfaces, stacking disorders in ice nuclei, and the addition rate of particles to the critical nucleus.We disentangled and investigated the relative importance of all the terms, including interfacial free energy, entropic contributions and the kinetic prefactor, that contribute to the overall nucleation rate.There has been a long-standing discrepancy for the predicted homogeneous ice nucleation rates, and our estimate is faster by 9 orders of magnitude compared with previous literature values. Breaking down the problem into segments and considering each term carefully can help us understand where the discrepancy may come from and how to systematically improve the existing computational methods.

Figures (4)

• Figure 1: The light blue curves are the free energy profiles as a function of the collective variable $\Phi$ for 8 sets of umbrella sampling simulations, and the dark blue curve is the averaged $\tilde{G}(\Phi)$ from these runs. Each dotted light red curve is the free energy profile of a perfect Ih nucleus extracted from each set of simulations, the dotted dark red curve is the averaged result, and the orange curve indicates a CNT fit using Eqn. \ref{['eq:cnttol']}. The inset shows a snapshot of an ice Ih nucleus embedded in liquid water.
• Figure 2: Upper panel: The stacking fault free energy per area as a function of temperature. The black line is the estimate from the potential energy difference at 0 K, the blue curve represents the harmonic approximation (HAR), and the red dots show the results from thermodynamic integration that considers anharmonicity (ANH). Statistical uncertainties are indicated by the error bars. Middle panel: the free energy difference $\Delta G_{\rm sf}$ between a pure ice Ih nucleus and a one that has the same size and a mixed stacking disorder at temperatures 230 K, 240 K and 250 K as predicted by the analytic model illustrated in the inset. Lower panel: The free energy profile as a function of cubicity for nuclei of three different sizes at 240 K.
• Figure 3: The red curve is the free energy profile $G_{\rm{Ih}}$ of a pure Ih nucleus, the green curve is $G_{\rm{Ic}}$ of a pure Ic nucleus, and the black curve indicates the free energy profile of an ice nucleus that can have a mixed stacking order. The width of the curves indicates the statistical error in the free energy estimation, computed from the error of the mean from independent simulation runs.
• Figure 4: The orange, green and blue curves are the spectra $\omega S(\omega)$ for bulk solid, bulk liquid, and a solid-liquid system that contains a solid critical nucleus ($n^{\star}=550$) at 240 K and 1 bar, respectively. The red curve is the fitting curve using Eqn. \ref{['eq:somega']} with parameters $m_f=3\times10^{-5}~{\rm ps^2 \,kJ/mol}$, $\kappa_f=0.03~{\rm kJ/mol}$, $\gamma_f=0.0025~{\rm ps \,kJ/mol}$, and $\gamma=0.06~{\rm ps \, kJ/mol}$.