2. Results

The two-phase system in the isothermal condition is described here by following the DI approach through the (Landau) free energy functional

(2.1) \begin{equation} \varOmega \left [\rho \right ] = \int _{V} \left [ f_b\left (\rho \right )+ f_c \left ({\boldsymbol{\nabla }}\rho \right ) - \mu _{eq} \rho \right ] \,{\textrm{d}}V + \oint _{\partial V} f_w \left (\rho \right ) \,{\textrm{d}}S \, , \end{equation}

where the van der Waals’ square gradient approximation is used to express the non-local excess free energy contribution as $f_c = (\lambda /2)\vert {\boldsymbol{\nabla }}\rho \vert ^2$ , while $f_b$ is the classical Helmholtz bulk free energy density and $\mu _{eq}$ the equilibrium chemical potential, see the Supplementary materials for additional information is available at https://doi.org/10.1017/jfm.2025.10591. The capillary coefficient $\lambda$ controls the liquid–vapour interface properties, namely the surface tension and the interface thickness. The temperature dependence of all the terms has been neglected for the ease of notation. Finally, the surface contribution $f_w$ arises as a mean field approximation of the fluid–wall interactions and accounts for the wetting properties of the surface. Its explicit expression, derived in Gallo, Magaletti & Casciola Reference Gallo, Magaletti and Casciola2021,

(2.2) \begin{equation} f_w(\rho ) = f_w\left (\rho _V\right ) -\cos \theta \, \int _{\rho _V}^\rho \sqrt {2\lambda \left [w_b\left (\tilde {\rho }\right ) - w_b\left (\rho _V\right )\right ]} \,{\textrm{d}} \tilde {\rho } \, , \end{equation}

has been shown to recover the well-known Young condition for the equilibrium contact angle, $\theta$ . In the previous expression, $f_w(\rho _V) = \gamma _{SV}$ is the surface energy at the solid–vapour interface, $\rho _V$ is the vapour density and $w_b = f_b - \mu _{sat}\rho$ , with $\mu _{sat}$ the chemical potential at the saturation condition ( $\mu = \partial f_b/\partial \rho$ ). The starting point for deriving the form of the wall free energy is based on geometric considerations regarding the orientation of the interface normal, aligned with the density gradient, and its contact with the solid surface, along with the expression for surface tension from DI theory, see Gallo et al. Reference Gallo, Magaletti and Casciola2021 for details. In this work, the modified Benedict–Webb–Rubin equation of state (Johnson, Zollweg & Gubbins Reference Johnson, Zollweg and Gubbins1993), mimicking the behaviour of a Lennard-Jones fluid, has been used as $f_b(\rho )$ . Here $f_b(\rho )$ is rendered dimensionless using the reference energy $f_b{_{\textit{ref}}} = \epsilon /\sigma ^3$ where $\sigma = 3.4\times 10^{-10} \, \text{m}, \, \epsilon = 1.65\times 10^{-21} \, \text{J}$ . The density field is non-dimensionalised with the reference density, $\rho _{\textit{ref}} = m/\sigma ^3$ , with $ m = 6.63\times 10^{-26} \, \text{kg}$ . The capillary coefficient is set to $\lambda = 5.224$ with its reference value $ \lambda _{\textit{ref}} = \sigma ^5 \epsilon /m^2$ ensuring consistency with surface tension values derived from Monte Carlo simulations (Gallo, Magaletti & Casciola Reference Gallo, Magaletti and Casciola2018; Magaletti et al. Reference Magaletti, Gallo, Perez, Carrillo and Kalliadasis2022).

The minimisation of (2.1) leads to the following Euler–Lagrange equation:

(2.3) \begin{equation} \dfrac {\delta \varOmega }{\delta \rho } = \mu -\lambda {\nabla} ^2 \rho - \mu _{eq} = 0 \,, \end{equation}

expressing the chemical equilibrium in terms of a constant (generalised) chemical potential. This equilibrium condition is supplemented with the boundary condition $\partial f_w/\partial \rho + \lambda \partial \rho /\partial n = 0$ at the solid surface. When the prescribed equilibrium chemical potential corresponds to a metastable condition, $\mu _{\textit{spin}} \lt \mu _{eq} \lt \mu _{sat}$ , (2.3) has three solutions: (i) the homogeneous liquid $\rho ({\boldsymbol{x}}) = \rho _L$ ; (ii) the homogeneous vapour $\rho ({\boldsymbol{x}}, t) = \rho _V$ ; (iii) a two-phase solution with a bubble of a given radius surrounded by the metastable liquid, the critical nucleus representing the transition (or critical) state, i.e. the density configuration $\rho _{crit}({\boldsymbol{x}})$ that must be reached to trigger the complete transition from the metastable liquid to the more stable vapour state. Several previous studies have addressed nucleation within frameworks similar to the one employed in the present work (Talanquer & Oxtoby Reference Talanquer and Oxtoby1996; Baidakov Reference Baidakov2016). However, these approaches have been limited to the calculation of critical nucleus profiles and the subsequent estimation of nucleation rates, without resolving the full transition pathway using a rare-event technique. We address here the problem of evaluating the non-trivial solution of case (iii) by exploiting the powerful string method (E et al. Reference E, Ren and Vanden-Eijnden2007; Magaletti et al. Reference Magaletti, Gallo and Casciola2021), with the twofold objective of determining both the density field of the critical nucleus and the complete MEP, which describes the system configurations $\rho (s,{\boldsymbol{x}})$ along a suitable reaction coordinate $s$ advancing through the transition. As detailed in Appendix B, the MEP corresponds to the MLP for a gradient dynamics and thus provides a physically grounded description of the liquid–vapour phase transition. This path can be visualised as a curve (the string), parametrised with $s$ , in the infinite configurational space of the density fields, with the configuration of maximum energy along the MEP representing the transition state. For the homogeneous case, the bubble radius, in each string image along the path, is extracted a posteriori from the density field. In the spherical case, a natural radius definition is

(2.4) \begin{align} R(s) = \int _0^\infty r (\partial \rho (r,s)/\partial r)^2 \,{\textrm{d}}r/\int _0^\infty (\partial \rho (r,s)/\partial r)^2 \,{\textrm{d}}r \, , \end{align}

Figure 1.

(a) Comparison of energy landscapes for homogeneous and heterogeneous bubble nucleation predicted by CNT (dashed lines) and DI model (solid lines). The heterogeneous case corresponds to a neutrally wetting surface ( $\theta = 90^\circ$ ). All curves are computed at $T=1.2$ and $\mu _{lev}=0.2$ . (b) Transition path of homogeneous nucleation projected onto the two-coordinate space $\{\rho _{av}, R\}$ at $T=1.2$ , $\mu _{lev}=0.9$ , Arrows identify the phase change direction, red DI and blue CNT, respectively. Transition states from CNT and DI models are marked with full and empty circles, respectively. Red squares indicate $\rho _{av}$ and $R$ during a FH simulation. Inset: transition paths at varying metastabilities, with circles marking the corresponding transition states. The MEPs of heterogeneous nucleation in $\{\rho _{av}, V\}$ space for different surface wettabilities: (c) $T=1.2$ , $\mu _{lev}=0.2$ ; (d) $T=1.2$ , $\mu _{lev}=0.6$ .

which, as will be discussed below, provides a good indicator for the bubble radius converging to the classical one used in sharp interface models. In classical theories, the bubble radius (or its volume, equivalently) is used as the reaction coordinate to describe the progress of the nucleation process, both in homogeneous and heterogeneous conditions. Here CNT is often used to quantitatively depict the energy landscape of the system in terms of the bubble size, with the energy maximum representing the barrier for nucleation and its corresponding critical radius. In figure 1(a), the landscapes are plotted in terms of the normalised grand potential $\Delta \varOmega /(k_B T) = (\varOmega [\rho ] - \varOmega [\rho _L])/(k_B T)$ as a function of the bubble radius $R$ . Both homogeneous and heterogeneous CNT predictions, at non-dimensional temperature $T=1.2$ ( $T\simeq 1.31$ being the critical temperature) and metastability level $\mu _{lev}= (\mu _{eq}-\mu _{sat})/(\mu _{\textit{spin}}-\mu _{sat}) = 0.2$ , are compared with those obtained by applying the string method to the DI modelling of the two-phase system, as previously described. For $T=1.2$ , we have $\mu _{sat} = -3.9506$ , $\mu _{\textit{spin}} = -4.0491$ , $\rho _{Lsat} = 0.5669$ and $\rho _{Lspin} = 0.4798$ . The heterogeneous landscape in the plot refers to the specific case of a neutrally wetting surface ( $\theta =90^o$ ). Despite the low metastability level – corresponding to a condition close to saturation – the discrepancy between the different energy barrier predictions is apparent in both cases, with a difference of the order of $6\,\%$ , in line with the results in Shen & Debenedetti Reference Shen and Debenedetti2001. The differences in the nucleation barrier are minor, and the critical radii are accurate when approaching saturation conditions. However, discrepancies in the barrier become significantly more pronounced near the spinodal, where CNT predicts a finite barrier, whereas it should vanish. This aspect will be discussed in detail later in the paper. The major difference, however, is the behaviour at low energies where the DI model predicts a non-bijective correspondence between the radius and the energy, suggesting that the radius alone is insufficient to fully describe the nucleation process.

New insights on the transition pathway are gained when the homogeneous nucleation MEP is plotted using a two-coordinate system, $\{\rho _{av}, R\}$ , with the average density inside the bubble a posteriori evaluated along the string as $\rho _{av} = 3/(4 \pi R^3) \int _{0}^{R} \rho (r) 4\pi r^2 \,{\textrm{d}}r$ , see figure 1(b). The DI model reveals a non-classical C-shape nucleation pathway. In line with the findings on crystal nucleation (Lutsko Reference Lutsko2019; Lutsko & Lam Reference Lutsko and Lam2020), the process starts with a spatially extended density variation of small intensity, hence characterised by a large radius and $\rho _{av} \simeq \rho _L$ . Successively, the embryo spatially localises, and its density variation increases. After the transition state is reached, the bubble further grows and expands, see the red arrows indicating the direction of the transition. The characteristic C-shaped form of the transition pathway is not sensitive to the specific definition of the radius and can also be observed when using the equimolar radius, as shown in the Supplementary materials. The position of the transition state along the curve strongly depends on the degree of metastability $\mu _{lev}$ , as shown in the inset in figure 1(b): a large and low-density critical embryo characterises conditions close to saturation (small $\mu _{lev}$ ), in agreement with CNT; when approaching spinodal conditions the critical embryo is instead characterised by small size and a high average density. The C-shape pathway is also confirmed by dynamic simulations with the spherically averaged FH model proposed by these authors in Gallo et al. Reference Gallo, Magaletti, Cocco and Casciola2020, see Appendix A for details. Due to thermal fluctuations, the system explores a pseudotube in the $\{\rho _{av}, R\}$ space around the zero-temperature MEP provided by the string method. Fluctuations of the average density intensify as the vapour embryo localises in a small region and, consequently, triggers the formation of a low-density cluster that successively expands beyond the transition state. Please note that fluctuations increase when reducing the average (bubble) volume (Gallo Reference Gallo2022). We observe that the MEP is perfectly followed by brute force FH simulations, as demonstrated in figure 1(b). This shows unequivocally that, among all the possible thermally induced fluctuations (always present also in stable liquids), the relevant ones for triggering the transition from the metastable liquid state are indeed those predicted by the MEP we have calculated. In addition, despite the nucleation path in FH being expected to differ from the MEP (Grafke, Grauer & Schäfer Reference Grafke, Grauer and Schäfer2015; Yao & Ren Reference Yao and Ren2022), the bubble nucleation mechanism appears to be well described by free-energy calculations.

Heterogeneous nucleation follows a similar pathway (with the volume $V$ instead of the radius) with a small but spatially extended density variation as a precursor of the actual bubble formation. Results at different wettabilities and degrees of metastability are plotted in figure 1(c,d). As expected, the volume of the critical nuclei decreases as the contact angle increases, namely, hydrophilic surfaces require a larger critical bubble to initiate the phase transition compared with hydrophobic walls. In addition, the critical volume decreases when the metastability level is increased.

Figure 2.

(a) Normalised free-energy barrier $\Delta \varOmega ^\star /k_B T$ versus $\mu _{lev}$ . Dashed lines refer to CNT, while solid lines refer to the DI. Inset: the left-hand axis reports $R$ versus $\mu _{lev}$ , the right-hand axis depicts $R$ normalised with interface thickness $l_{10-90}$ versus $\mu _{lev}$ . (b) Density profiles along the transitions, precritical and postcritical profiles are depicted in blue and red, respectively. The critical profile is reported in black ( $\mu _{lev} = 0.8$ , $T=1.20$ ). (c) Here $(\partial \rho /\partial r)^2$ normalised with its $L_2$ norm versus radial coordinate. The panel refers to precritical states, while the inset refers to postcritical conditions. (d) Energy landscape as a function of the tuple $(\rho _{av}, R)$ , $\mu _{lev} = 0.2$ , $T=1.20$ . All cases refer to homogeneous nucleation, the reaction coordinate $s$ increasing directions are also indicated.

In the case of homogeneous nucleation, figure 2(a) reports the nucleation barrier normalised by $k_B T$ as a function of the degree of metastability. As shown, CNT captures well the qualitative trend of the barrier with increasing metastability. However, as the system approaches the spinodal limit, $\mu _{{lev}} \to 1$ , CNT predicts a large limit barrier of approximately $ \Delta \varOmega ^\star _{CNT} \simeq 18 k_B T$ , whereas the DI approach yields the expected vanishing barrier $\Delta \varOmega ^\star _{DI} \simeq 0$ . This discrepancy arises from the gradual breakdown of the sharp interface assumption as the liquid–vapour interfacial thickness becomes comparable to the typical bubble radius. This effect is illustrated in the inset of figure 2(a) (red curve right-hand axis), which shows the ratio of the bubble radius to the interface thickness $l_{10-90}$ , defined as the width of the region over which the density transitions from $\rho _{10} = 0.1 \rho _L + 0.9 \rho _V$ to $\rho _{90} = 0.9 \rho _L + 0.1 \rho _V$ (Caupin Reference Caupin2005). On the right-hand axis, CNT (dashed line) and DI (solid line) radii are depicted. The DI predicts a non-monotonic behaviour of the critical radius as the spinodal is approached, in contrast to the monotonic trend obtained from CNT. Nevertheless, the actual values of the critical radii remain in reasonably good agreement over the range of metastabilities considered. Figure 2(b) shows the transition pathway for homogeneous nucleation in terms of the density profiles $\rho (r)$ along the reaction coordinate $s$ of the string. Precritical profiles are shown in blue, postcritical ones in red, and the critical (saddle-point) profile in black. Based on the structure of these profiles, we evaluate the quantity $(\partial \rho / \partial r)^2 = \rho '(r)^2$ , normalised by its $L_2$ norm, $||\rho '(r)||_{L_2}$ . This normalised profile offers insight into the interfacial region and, in the sharp interface limit, serves as an indicator of the interface location. Specifically, as the interface thickness tends to zero, $\rho '^2/||\rho '(r)||_{L_2}^2 \to \delta (r - R)$ , converging to the sharp interface limit, where $\delta$ denotes the Dirac delta function. Within the DI framework, this provides a generalised representation of the interface. This normalised gradient is reported in figure 2(c), where the panel shows precritical profiles and the inset includes postcritical profiles and the saddle point. As indicated by the arrow representing the direction of the transition (increasing $s$ ), the initial profiles are broad, with maxima located at larger distances $r$ and only modest density variations (red curves in figure 2 b), corresponding to a slightly rarefied liquid. As the transition proceeds, the maxima shift towards smaller radii and eventually increase in amplitude near the saddle point and beyond (blue curves in figure 2 b), consistent with the complex C-shaped transition pathway observed in the $(R, \rho _{{av}})$ plane. To better highlight the nucleation mechanism, figure 2(d) also reports the free energy landscape as a function of the two variables $R$ and $\rho _{{av}}$ , with the direction of the reaction coordinate $s$ indicated.

Figure 3.

(a) Density fields along the MEP, showing nucleation progress from (i,ii) to (vii,viii). Panels (a iii) and (a iv) correspond to the transition states: (iii) hydrophilic case at $T=1.2$ , $\mu _{lev}=0.5$ , $\theta =30^\circ$ ; (iv) hydrophobic case at $T=1.2$ , $\mu _{lev}=0.5$ , $\theta =110^\circ$ . (b) Energy barrier ratio between heterogeneous and homogeneous nucleation as a function of contact angle, obtained using the string method with the DI model. Symbols indicate different levels of metastability; the solid black curve shows the CNT prediction, based solely on the geometrical factor $\varPsi$ . (c) Mean first passage time for homogeneous nucleation versus metastability. The solid line is the DI model prediction, while red squares refer to Corrected CNT prediction (Menzl et al. Reference Menzl, Gonzalez, Geiger, Caupin, Abascal, Valeriani and Dellago2016), and the blue triangle corresponds to brute force FH simulations. (d) The DI model prediction of heterogeneous nucleation’s mean first passage time. Each curve corresponds to a different contact angle.

A further comparison with heterogeneous CNT is shown in figure 3(b) where the ratio between the energy barriers of the heterogeneous and homogeneous cases, $\Delta \varOmega _{het}^*/\Delta \varOmega _{hom}^*$ , is reported as a function of the contact angle. The CNT again provides a theoretical prediction, with this ratio equal to the purely geometrical function $\varPsi (\theta ) = 1/4(2 - \cos (\phi ))(1+ \cos (\phi ))^2$ . The numerical results obtained with the string method applied to the DI model show, instead, a behaviour strongly dependent also on the degree of metastability. Close to saturation, $\mu _{lev}=0.2$ , CNT predictions are well reproduced by the DI model. Mesoscale properties of the critical bubble start to be effective at larger metastability levels. When the interface thickness is comparable to the critical bubble dimension, at $\mu _{lev}=0.5$ and even more apparently at $\mu _{lev}=0.8$ , the deviation from CNT is more pronounced. It is worth noticing that at high metastability, the hydrophobic surface can anticipate the spinodal condition. The energy barrier at $\mu _{lev}=0.8$ and $\theta =110^{\circ}$ is actually zero, suggesting a barrierless mechanism of phase transition typical of the spinodal decomposition. This result is consistent with the findings of Talanquer & Oxtoby Reference Talanquer and Oxtoby1996, where the so-called surface spinodal was identified using a different form of fluid–solid free-energy. Moreover, on the hydrophilic side of the plot, the case at $\mu _{lev}=0.8$ shows a heterogeneous energy barrier independent of the contact angle up to $60^{\circ}$ and with a value slightly smaller than the homogeneous barrier. To understand this peculiar behaviour, it is instrumental to look at the full-density fields during the transition progress, shown in figure 3(a). The highly hydrophilic surface, even at a moderate metastability level, induces the localisation of the density variation in a region close to the wall but not directly in contact with it. As a consequence, the critical nuclei are almost spherical, resembling homogeneous nucleation. This behaviour has also been observed in full three-dimensional FH simulations (Gallo et al. Reference Gallo, Magaletti and Casciola2021), in MD simulations (Nagayama, Tsuruta & Cheng Reference Nagayama, Tsuruta and Cheng2006; Chen et al. Reference Chen, Zou, Wang, Han and Yu2018; Zou et al. Reference Zou, Gupta and Maroo2018) and it could be ascribed to the mesoscale interaction between the vapour–liquid interface and the fluid–solid layering induced by the strong attraction of the hydrophilic surface. The configurations along the MEP show that, after the transition state, the bubble grows at the solid surface, recovering the expected behaviour always observed in experiments. Conversely, in the hydrophobic case, the embryos always sit at the solid wall due to the high affinity with the vapour.

A key observable in stochastic processes such as nucleation is the mean waiting time for forming a critical nucleus. This time scale is usually prohibitively long for direct brute-force simulations (e.g. MD or FH), making theoretical approaches essential. Kramers’ theory provides a fundamental framework for estimating the characteristic nucleation rate, see Hänggi et al. Reference Hänggi, Talkner and Borkovec1990 for a general discussion on barrier-crossing problems and Gallo et al. Reference Gallo, Magaletti, Cocco and Casciola2020 for bubble nucleation application. The theory requires the free energy profile and the diffusion coefficient of the fluctuating vapour nucleus embedded in the liquid. Having identified the MEP, described through parameterisation with the curvilinear abscissa $s$ , we can conjecture the simplest dynamics for the reaction coordinate to determine the nucleation times. Such dynamics sees the thermodynamic system as a Brownian walker in the energy landscape $\varOmega = \varOmega (s)$ ,

(2.5) \begin{equation} \dfrac {{\textrm{d}}{s}}{{\textrm{d}}{t}} = -\alpha \dfrac {{\textrm{d}} \varOmega (s)}{{\textrm{d}} s} + \sqrt {2 D}\eta (t) \end{equation}

with $\alpha$ a friction coefficient, $D = k_B T \alpha$ the diffusion and $\eta (t)$ a white noise with zero mean and $\langle \eta (t)\eta (q)\rangle = \delta (t-q)$ . The friction coefficient is estimated using the following procedure: the critical bubble $ \rho _{{crit}} ({\boldsymbol{x}}, s^\star )$ is identified and perturbed along the directions of the two stable basins, $\rho _V$ and $ \rho _L$ . The perturbed configurations are then allowed to evolve according to the isothermal NSK equations

(2.6) \begin{align} \dfrac {\partial \rho }{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\rho {\boldsymbol{v}}) &= 0, \\ \nonumber \dfrac {\partial \rho {\boldsymbol{v}}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\rho {\boldsymbol{v}} \otimes {\boldsymbol{v}}) &= - \rho \boldsymbol{\nabla }\left (\dfrac {\delta \varOmega }{\delta \rho }\right ) + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\varSigma }, \end{align}

where $ \rho$ is the density, $ {\boldsymbol{v}}$ is the velocity field and $\boldsymbol{\varSigma } = \eta ( \boldsymbol{\nabla }{\boldsymbol{v}} + (\boldsymbol{\nabla }{\boldsymbol{v}})^T ) + ( \zeta - 2/3 \eta ){} (\boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{v}}) \textit{I}$ is the viscous stress tensor, $\eta , \zeta$ the two viscosity coefficients. We adopt the viscosity expression proposed by Rowley & Painter Reference Rowley and Painter1997 to ensure full consistency with the Lennard-Jones properties. By monitoring the energy variation $\delta \varOmega$ occurring as the system transitions between two nearby configurations $\rho ({\boldsymbol{x}}, s^\star )$ and $\rho ({\boldsymbol{x}}, s^\star \pm \delta s)$ , with $\delta s = ||\rho ({\boldsymbol{x}}, s^\star \pm \delta s) - \rho ({\boldsymbol{x}}, s^\star ) ||_{L^2}$ and the time required for this transition $\delta t$ (as measured in NSK equations), the friction $ \alpha$ is determined as $\alpha = - (\delta s/\delta t)/(\delta \varOmega /\delta s)$ . Specifically, as detailed in Appendix B, the string is a collection of fields $\rho ({\boldsymbol{x}}, s)$ . We aim to estimate the friction coefficient at the saddle point. Starting from the saddle point, defined by the field $\rho _{{crit}}({\boldsymbol{x}}) = \rho ({\boldsymbol{x}}, s^\star )$ , we introduce perturbations in the directions of the basins of the liquid and vapour by considering the configurations just before and after the saddle point, $\rho ({\boldsymbol{x}}, s^\star \pm \delta s)$ . These are then evolved under the NSK equations towards the respective basins of liquid and vapour. By measuring the time taken to reach the subsequent configurations $\rho ({\boldsymbol{x}}, s^\star \pm 2 \delta s)$ , we obtain two values of $\alpha$ , namely $\alpha _r = - (\delta s/\delta t_r)/(\delta \varOmega _r/\delta s)$ and $\alpha _l = - (\delta s/\delta t_l)/(\delta \varOmega _l/\delta s)$ , corresponding to the right- and left-hand sides of the saddle point, with $\delta \varOmega _r = \varOmega (s^\star + 2 \delta s) - \varOmega (s^\star + \delta s)$ and $\delta \varOmega _l = \varOmega (s^\star - 2 \delta s) - \varOmega (s^\star - \delta s)$ , and $\delta t_{r/s}$ the time measured from NSK numerical simulations. We then define the effective diffusion coefficient at the saddle as $\alpha ^\star = (\alpha _r + \alpha _l)/2$ . The parameter $\alpha$ has physical dimensions of energy times density squared per unit time, while $s$ has the dimensions of a density. This is only an estimate, as the MEP does not represent the MLP of the NSK dynamics, given that the two systems evolve in different spaces and their trajectories are not necessarily close (Grafke & Vanden-Eijnden Reference Grafke and Vanden-Eijnden2019; Zakine & Vanden-Eijnden Reference Zakine and Vanden-Eijnden2023). However, numerical results suggest that the forecast is remarkably accurate, as discussed subsequently. The procedure was repeated for both the homogeneous and heterogeneous cases. After estimating the friction coefficient, the mean first passage time $\tau$ was evaluated by using Kramers’ theory (Kramers Reference Kramers1940; Schulten, Schulten & Szabo Reference Schulten, Schulten and Szabo1981),

(2.7) \begin{equation} \tau = \int _{B_{L}} \exp \left (-\dfrac {\Delta \varOmega }{k_BT}\right ) {\textrm{d}}s \int _{{B}_{SP}} \dfrac {1}{D(s)} \exp {\left (\dfrac {\Delta \varOmega }{k_BT}\right )} {\textrm{d}}s, \end{equation}

where $B_{L}$ and ${B}_{SP}$ are the metastable and saddle-point basins. The two integrals can be evaluated using the classical saddle-point approximation, after noting that the free energy profiles’ curvatures ( $\varOmega ''(s)$ ) are positive and negative in the two basins, respectively. The second integral can be expanded as

(2.8) \begin{align} &\int _{{B}_{SP}} \dfrac {1}{D(s)} \exp {\left (\dfrac {\Delta \varOmega }{k_BT}\right )} {\textrm{d}}s \sim \exp \left (\dfrac {\Delta \varOmega ^\star }{k_BT}\right ) \int _{-\infty }^{\infty } \exp \left (-\dfrac {(s-s^\star )^2}{2\nu ^2}\right ) \nonumber \\&\quad \times \, \dfrac {1}{D^\star + D^{\prime}_\star (s-s^\star ) + 1/2 D^{\prime\prime}_\star (s-s^\star )^2 + \ldots } d(s-s^\star ) , \end{align}

with $\nu ^2 = k_B T /|\varOmega ''(s^\star )|$ . Our data show that the diffusion coefficient is gently varying in the neighbourhood of the saddle point, leading to the conclusion that, within an order one prefactor,

(2.9) \begin{equation} \tau \sim \dfrac {2\pi }{ \alpha (s^\star ) \sqrt {\varOmega ''(0) |\varOmega ''(s^\star )|}} \exp \left (\dfrac {\Delta \varOmega ^\star }{k_B T}\right )\!. \end{equation}

Here $\tau$ is reported for homogeneous and heterogeneous nucleation in figures 3(c) and 3(d), respectively. Figure 3( c) shows the first passage time as a function of the metastability parameter $ \mu _{{lev}}$ . The time is non-dimensionalised using $ t_r = 2.15 \times 10^{-12}$ s as a reference. The solid black line depicts our prediction, while the red squares represent the mean first passage time evaluated with an adaptation of the nucleation theory proposed by Menzl et al., (Menzl et al. Reference Menzl, Gonzalez, Geiger, Caupin, Abascal, Valeriani and Dellago2016) (see below), and the blue triangle depicts the FH brute force simulations. Kramers’ theory becomes increasingly accurate in the limit of large energy barriers (Hänggi et al. Reference Hänggi, Talkner and Borkovec1990), see Gallo et al. Reference Gallo, Magaletti, Cocco and Casciola2020, where the agreement with brute-force simulations emerges for energy barriers $\sim 6-7 k_B T$ , at a metastability level of roughly $\mu _{lev} = 0.82$ , corresponding to a reduced density of $\rho _L = 0.475$ with $T=1.25$ . In the present case, for the highest metastability level considered ( $\mu _{lev} = 0.8$ ), the energy barrier ( $\Delta \varOmega ^\star \sim 10 k_B T$ ) is well in the range where Kramers’ theory is expected to be robust. This conclusion is supported by 200 FH simulations we performed using the equations reported in Appendix B. The resulting mean first passage for all cases is consistently in the range of $\tau \sim 10^6$ ( $\tau _{{FH}} \sim 0.25 \times 10^6$ , $\tau _{DI} \sim 1.00 \times 10^6$ , $\tau _{cCNT} \sim 3.60 \times 10^6$ ) confirming the robustness of Kramers’ theory even at moderate metastability. Obviously, also in the heterogeneous case, the accuracy of Kramers’ theory may deteriorate for low energy barriers. Such regimes represent only a minor portion of the results presented here, which can, however, be directly investigated via brute-force simulations.

The Menzl et al., theory – applied by the authors to the TIP4P/2005 water model – combines MD microscopic information (vapour embryo surface energy) with the stochastic overdamped Rayleigh–Plesset (RP) equation to estimate the diffusion coefficient $D$ . It is worth noting that the procedure proposed by Menzl et al., is completely general. It consists of correcting the CNT barrier height by either taking the surface energy from MD simulations or by applying a correction using the Tolman length, which is estimated by adjusting the CNT barrier relative to MD results. The diffusion coefficient of the nucleating bubble is estimated via the fluctuation–dissipation theorem applied to an overdamped RP dynamics, requiring that it samples the equilibrium probability distribution function of the radius $R$ , $P(R) \sim \exp (-\Delta \varOmega (R) / k_B T)$ . The diffusion coefficient is then computed in closed form as $D_R = k_B T /(16\pi \eta R)$ . In the present case, which deals with a Lennard-Jones fluid, the DFT barrier automatically incorporates the curvature correction (Blokhuis & Kuipers Reference Blokhuis and Kuipers2006; Wilhelmsen, Bedeaux & Reguera Reference Wilhelmsen, Bedeaux and Reguera2015; Rehner et al. Reference Rehner, Aasen and Wilhelmsen2019; Magaletti et al. Reference Magaletti, Gallo and Casciola2021) while the diffusion coefficient for reaction variable $s$ is estimated using the hydrodynamic equations with the procedure mentioned above. With this approach, we can directly compare the results of the two theories. We found perfect agreement between the two approaches. The theoretical framework developed by Menzl et al., is not directly applicable to heterogeneous nucleation, as the RP equations describe bubble dynamics in an unbounded medium. In contrast, the approach proposed here is more general, as it allows for the determination of the MEP and the estimation of diffusion coefficients under arbitrary conditions. In figure 3(d), the heterogeneous case is presented, highlighting how the wall chemistry plays a crucial role in promoting nucleation by significantly reducing the characteristic transition times. It can be observed that when the contact angle is below approximately $50^\circ$ , the time scales become independent of wettability, approaching the homogeneous limit. This result is consistent with our previous findings of FH for boiling (Gallo et al. Reference Gallo, Magaletti, Georgoulas, Marengo, De Coninck and Casciola2023) and with MD simulations (Zou et al. Reference Zou, Gupta and Maroo2018; Sullivan et al. Reference Sullivan, Dockar and Pillai2025) identifying the homogeneous nucleation as the main transition mechanism when hydrophilic chemistries are considered.

As a final analysis, we examine the velocity fields and the dissipation function of the vapour bubble as it evolves from the saddle point towards the two basins, corresponding to the metastable liquid and the stable vapour. To this end, we introduce the total free energy of the system, defined as the sum of the grand potential and the kinetic energy of the fluid

(2.10) \begin{equation} H [\rho , {\boldsymbol{v}}] = \varOmega [\rho ] \!+\! K[\rho , {\boldsymbol{v}}] = \int _{V} f_b\left (\rho \right )+ \dfrac {\lambda }{2} | {\boldsymbol{\nabla }}\rho |^2 \!+\! \dfrac {1}{2} \rho |{\boldsymbol{v}}|^2 \!-\! \mu _{eq} \rho \,{\textrm{d}}V \!+\! \oint _{\partial V} f_w \left (\rho \right ) \,{\textrm{d}}S \, . \end{equation}

The dissipation reads

(2.11) \begin{align} \dfrac {{\textrm{d}}H}{{\textrm{d}} t} &= \int _{V} \dfrac {\delta H}{\delta {\boldsymbol{v}}} \boldsymbol{\cdot }\dfrac {\partial {{\boldsymbol{v}}}}{\partial t} + \dfrac {\delta H}{\delta \rho } \dfrac {\partial \rho }{\partial t} \,{\textrm{d}}V = - \int _{V} {{\boldsymbol{v}}} \boldsymbol{\cdot }\left ( \rho {{\boldsymbol{v}}}\boldsymbol{\cdot }\boldsymbol{\nabla }{{\boldsymbol{v}}} + \rho \boldsymbol{\nabla }\left (\dfrac {\delta \varOmega }{\delta \rho }\right ) + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\varSigma } \right ) \nonumber \\&\quad +\, \left ( \left (\dfrac {\delta \varOmega }{\delta \rho }\right ) + \dfrac {1}{2} |{{\boldsymbol{v}}}|^2\right ) \boldsymbol{\nabla }\boldsymbol{\cdot }\rho {{\boldsymbol{v}}} \nonumber \\ &= -\int _V 2\eta \, \boldsymbol{\tilde E}: \boldsymbol{\tilde E} + \dfrac {1}{3} \zeta (\boldsymbol{\nabla }\boldsymbol{\cdot }{{\boldsymbol{v}}})^2 \,{\textrm{d}}V - \oint _{\partial V} \left (\left (\dfrac {\delta \varOmega }{\delta \rho }\right ) + \dfrac {1}{2} |{{\boldsymbol{v}}}|^2\right ) \rho {{\boldsymbol{v}}}\boldsymbol{\cdot }\hat {\boldsymbol{\nu }} {\textrm{d}}S \, , \end{align}

where (2.6) have been invoked and $ \boldsymbol{\tilde E}$ is the deviatoric part of the symmetric part of the velocity gradient ${\boldsymbol{E}} =\text{sym}({\boldsymbol{\nabla }}\otimes {{\boldsymbol{v}}})$ . In the numerical simulations, the system is large enough that the velocity at the boundary is always zero, so the boundary term in the above equation vanishes and ${\textrm{d}}H/{\textrm{d}}t \lt 0$ . It is worth noting that, if the total chemical potential $\delta \varOmega /\delta \rho + 1/2 |{{\boldsymbol{v}}}|^2$ is set to zero, then $H$ remains monotonically decreasing in time as well. The simulations were performed in a domain of length $L = 3000$ (dimensionless units), which corresponds to approximately one micron in physical length. Since the present analysis is carried out for homogeneous nucleation, we exploit spherical symmetry to solve (2.6) (Abbondanza et al. Reference Abbondanza, Gallo and Casciola2023b ). We adopt the same grid space $\Delta r$ , temporal integrator and integration time step $\Delta t$ of stochastic equations (see Appendix A for details). The thermodynamic conditions are $T = 1.20$ and $\mu _{lev} = 0.5$ . In figure 4, the time evolutions of $H(t)$ (red curves) and ${\textrm{d}}H/{\textrm{d}}t$ (blue curves) are shown during the hydrodynamic evolution of the bubble in the precritical and postcritical states, in figures 4(a) and 4(c), respectively. Analogously, figures 4(b) and 4(d) display the corresponding velocity fields. In both configurations, $H(t)$ decreases monotonically over time, with greater dissipation observed in the precritical phase, characterised by a more abrupt change in the velocity fields. It is particularly interesting to note (figure 4 b) that the bubble initially shrinks (black curve $t=10$ and blue curve $t=20$ ), then slightly compresses the liquid and emits a wave that propagates through the liquid with the isotermal speed of sound velocity $c_T = \sqrt {\partial p /\partial \rho \vert _T}$ (purple and green curves at $t=75$ and $t=100$ ). In the postcritical phase, the bubble expands with lower velocities compared with its precritical counterpart.

Figure 4.

(a,c) Time evolution of the dissipation rate ${\textrm{d}}H/{\textrm{d}}t$ (blue) and total free energy $H(t)$ (red) computed from the NSK dynamics during relaxation towards (a) the metastable liquid basin and (c) the stable vapour basin. (b,d) Velocity profiles at selected time instants during the relaxation towards (b) the metastable liquid and (d) the stable vapour state.