Heterogeneous bubble nucleation dynamics

Abstract

Heterogeneous nucleation is the most effective mechanism for the inception of phase transformation. Solid walls and impurities act as a catalyst for the formation of a new thermodynamic phase by reducing the activation energy required for a phase change, hence enhancing nucleation. The formation of vapour bubbles close to solid, ideally flat, walls is addressed here by exploiting a mesoscale description that couples diffuse interface modelling of the two-phase vapour–liquid system with fluctuating hydrodynamics, extending previous work by the authors on homogeneous nucleation. The technical focus of this work is to directly account for hydrophobic or hydrophilic walls through appropriate boundary conditions compliant with the fluctuation–dissipation balance, a crucial point in the context of fluctuating hydrodynamics theory. This methodology provides access to the complete dynamics of the nucleation process, from the inception of multiple bubbles up to their long-time macroscopic expansion, on time and spatial scales unaffordable by standard techniques for nucleation, such as molecular dynamics. The analysis mainly focuses on the effect of wall wettability on the nucleation rate, and, albeit qualitatively in agreement with classical nucleation theory predictions, it reveals several discrepancies to be ascribed to layering effects in the liquid close to the boundary and to bubble–bubble interactions. In particular, it is found that, close to moderately hydrophilic surfaces, the most probable nucleation events occur away from the wall through a homogeneous mechanism.

1. Introduction

Equilibrium systems exhibit thermal fluctuations. At the macroscopic level, in the so-called thermodynamic limit, such fluctuations are usually irrelevant. However, for small systems, below the micrometre scale, fluctuations are always crucial. Given the equilibrium configuration maximizing the system entropy under the relevant constraints, the probability density function (p.d.f.) of a fluctuation can be taken to be proportional to the exponential of the entropy decrease brought about by the fluctuation (Einstein 1956). The p.d.f. gives access to the equilibrium correlation function of the fluctuating field and provides information on the ensuing relaxation dynamics, which must obey the relevant fluctuation–dissipation balance (see e.g. De Groot & Mazur 2013). Within this general framework, under special circumstances, the effect of fluctuations emerges to the macroscale and determines the large-scale dynamics. In these most interesting cases, a (mesoscale) model able to bridge the gap between the macroscopic scale of interest and the microscopic level governed by fluctuations is urgently needed. After Landau & Lifshitz's seminal work (see the reprinted version in Landau & Lifshitz (1980)), several coarse-grained models have been proposed to include thermal fluctuations in hydrodynamic descriptions, at both continuum (Fox & Uhlenbeck 1970) and discrete (Español, Serrano & Öttinger 1999; Español, Anero & Zúñiga 2009) levels. These works contributed significantly to the growing field of ‘fluctuating hydrodynamics’, fostering interest in the numerical solution of the related stochastic partial differential equations (SPDEs) (Donev et al. 2010, 2014; Balboa et al. 2012; Delong et al. 2013). Mesoscale approaches, besides playing an important role in the theory of fluids, have significant impact in cases where stochastic fluctuations are crucial, e.g. for microfluidic design, biological systems like lipid membranes (Naji, Atzberger & Brown 2009), Brownian engines and molecular motors (Peskin, Odell & Oster 1993; Detcheverry & Bocquet 2012).

Among the variety of fluctuation-induced effects, nucleation, the incipit of phase transitions in metastable fluids, is ubiquitous. It occurs, for instance, during bubble cavitation (Brennen 2013) and boiling (Carey 2018) as well as in freezing rain (Cao et al. 2009) and in crystal formation more generally (Lutsko 2019). Metastability implies that the transition is an activated process, i.e. an amount of energy is required to overcome the barriers separating different metastable basins. In particular, a liquid sustains a tensile condition reaching pressures well below the equilibrium vapour tension without forming a vapour phase. Equivalently, metastability is observed in superheated fluids, where bubble formation occurs far above the boiling temperature. In these conditions, bubble nucleation must be interpreted in statistical terms, with the probability of phase change related to the level of superheating or stretching of the liquid (Brennen 2013). However, notwithstanding the extremely large tensile strength of ultra-pure water, which is able to sustain up to $1$ kbar tensions in specially designed experimental conditions (El Mekki Azouzi et al. 2013), bubbles are very common in real life. In fact, bubble nucleation is most often favoured by impurities or dissolved gas nuclei, which strongly lower the energy barrier and exponentially enhance the bubble formation rate. A similar effect is induced by solid boundaries. Close to the wall, the reduction of the energy barrier to nucleation depends on the wetting properties of the surface as measured by the contact angle. For instance, recent experimental works show how the wettability of ultra-smooth surfaces influences the onset temperature for pool boiling in superheated liquids (Bourdon et al. 2012, 2015; Malavasi et al. 2015). In all the aforementioned cases, energy barrier crossing is due to stochastic fluctuations, which activate the phase transition (Jones, Evans & Galvin 1999; Kashchiev & Van Rosmalen 2003; Lohse & Prosperetti 2016). The accurate and thorough modelling of thermodynamics (i.e. metastability) and thermal fluctuations is then crucial to correctly reproduce nucleation dynamics.

With few exceptions (see e.g. Allen, Valeriani & ten Wolde 2009; Marchio et al. 2019), our current understanding of nucleation is built on quasi-static descriptions such as classical nucleation theory (CNT) (Blander & Katz 1975) or more recent extensions thereof (Lutsko & Durán-Olivencia 2015). CNT provides the basic interpretative framework in terms of simple energetic arguments and allows the estimation of energy barriers and nucleation rates in homogeneous conditions. Although introducing substantial simplifications, e.g. the shape of the embryo or the existence of a sharp interface between phases that retain their bulk thermodynamic properties despite their nanometric size, CNT is a powerful and predictive tool for the investigation of nucleation phenomena. It is easily extended to heterogeneous conditions, e.g. solid boundaries (Ward et al. 1983). From a thermodynamic standpoint, there are two major contributions to the work of vapour bubble formation: (i) the (positive) energy spent to form the liquid–vapour interface, the product of surface tension and interface area; and (ii) the (negative) work released in transforming liquid into vapour, proportional to bubble volume and vapour–liquid pressure difference. The (positive) surface energy prevails at small radii, whereas at larger radii the negative volume contribution is dominant. The maximum work is achieved at the critical radius, which, in the context of CNT, corresponds to the energy barrier to be overcome for activating the phase change. As intuitively expected on geometrical grounds (see § 2 for a more complete discussion), the surface contribution is lowered for wall-attached bubbles, since the surface of a spherical cap is smaller than the corresponding full sphere, entailing a larger probability of nucleating a bubble at the surface than in the bulk.

For vapour condensation and solidification from dilute solutions, the embryos grow by aggregation of mother-phase molecules into the cluster, a process described in terms of kinetic theory and attachment/detachment rates (Oxtoby 1992). Bubbles are more elusive (Shen & Debenedetti 2003) since the embryo, basically a region depleted of molecules, is more like a void than a molecular cluster. Since CNT describes microscopic objects in macroscopic terms and deals with embryos as if they had the same uniform thermodynamic properties as the stable phase (vapour for bubble nucleation), it cannot predict the energy barriers vanishing at spinodal conditions. Other more fundamental approaches like density functional theory (DFT) (Oxtoby & Evans 1988; Shen & Debenedetti 2001) and CNT extensions (Lutsko & Durán-Olivencia 2015; Lutsko 2018) have been proposed, addressing and correcting some CNT artifacts. Brute-force molecular dynamics (MD) simulations, although in principle representing a powerful tool (Novak, Maginn & McCready 2007; Diemand et al. 2014), can hardly cope with metastability, given the large disparity between atomistic and nucleation time scales. Specialized rare-event techniques (Bolhuis et al. 2002; Allen, Frenkel & ten Wolde 2006; Menzl et al. 2016) overcome the problem but are still limited to equilibrium, quasi-static conditions excluding dynamic effects. Forward flux methods may (Wang, Valeriani & Frenkel 2008) partially remove the limitation of quasi-static approaches but share the problem of being confined to extremely small spatial and time scales. Indeed, in general, atomistic approaches, given their extremely high computational cost, are confined to small systems of a few nanometres in extension, in many cases very far from the target technological application.

Recently, a novel approach in the context of continuum mechanics, based on a diffuse interface description of the two-phase vapour–liquid system endowed with thermal fluctuations in the spirit of fluctuating hydrodynamics (Landau & Lifshitz 1980), has been exploited to address the bubble nucleation process (Gallo, Magaletti & Casciola 2017, 2018; Gallo et al. 2020) in homogeneous conditions and the boiling process close to neutrally wettable surfaces (Magaletti, Georgoulas & Marengo 2020). The model is robust from both the thermodynamic and the statistical points of view. For equilibrium systems, its deterministic ingredients can be traced back to the more fundamental DFT description with a squared-gradient approximation of the excess energy (Lutsko 2011) and, in general, can be framed in the context of non-equilibrium statistical mechanics, e.g. the GENERIC framework (Öttinger & Grmela 1997; Espanol 2001). The approach has been shown to be able to capture the rich hydrodynamics of a multiphase system, such as phase transformation, latent heat release, shock emission and topological changes (Magaletti, Marino & Casciola 2015; Magaletti et al. 2016). The square-gradient approximation was also used to address phase change and spreading of droplets (Teshigawara & Onuki 2010) or the thermodynamics of boiling (Laurila et al. 2012). The stochastic ingredients introduced to model thermal fluctuations obey the fluctuation–dissipation balance (Fox & Uhlenbeck 1970) and reproduce the statistics of the fluctuating fields (Chaudhri et al. 2014).

The aim of this work is to extend our previous studies to heterogeneous nucleation by focusing on spontaneous vapour bubble nucleation on a flat solid surface. A detailed derivation of the model and related boundary conditions is presented, and results of numerical simulations are produced to demonstrate the potential of the proposed approach. In doing so, some of the material discussed by the authors in previous papers – in particular, the diffuse interface model for non-isothermal fluids with a general equation of state (Magaletti et al. 2015, 2016) and fluctuating hydrodynamics for homogeneous nucleation (Gallo et al. 2017, 2018, 2020) – is reviewed to present a comprehensive description of the model. The paper is structured as follows. In § 2 a brief review of CNT is provided, with a particular focus on heterogeneous conditions and nucleation in the microcanonical $NVE$ ensemble (constrained number, volume and energy). In order to correctly take into account the solid–fluid interaction that determines the wetting properties of the surface in the diffuse interface modelling, an analytic form of the solid-wall free energy is derived in § 3. Its expression connects bulk-phase thermodynamics and contact angle with the fluid density at the wall. Next, § 4 is devoted to the Navier–Stokes equations with capillarity, with particular attention paid to the constitutive relations induced by the new boundary terms. The new energetic contribution associated with the solid–fluid interaction is found to modify the fluctuation statistics (§ 5), still preserving the fluctuation–dissipation balance in the general setting (§ 6). In (§ 7) explicit expressions are provided for a particular case. The whole model is assembled together in § 8, combining deterministic and stochastic contributions. In § 9 numerical validation against theoretical statistical properties of thermal fluctuations is provided and a detailed analysis of heterogeneous nucleation for changing surface wetting characteristics is presented. Conclusions and the future perspective are finally drawn in § 10.

2. Classical nucleation theory

Classical nucleation theory (CNT) (Ward et al. 1983; Kashchiev 2000; Brennen 2013) provides the fundamental understanding of bubble nucleation in a metastable liquid, for both homogeneous (bubble forming inside the bulk liquid) and heterogeneous conditions (bubble forming in contact with an extraneous phase, usually a solid with given geometry and chemical properties). Although classical, it is reviewed here to put the new material discussed in the paper in the proper perspective.

For a vapour bubble nucleating on a flat solid surface with prescribed contact angle $\phi$ (and neglecting gravity), CNT models the bubble as a spherical cap of radius $R$ on top of the flat solid surface (figure 1a). The system is described in terms of free energy difference with respect to the pure liquid,

(2.1)\begin{equation} {\rm \Delta} \varOmega (R, \phi ) = -{\rm \Delta} p \,V_V(R, \phi ) + \gamma_{LV} A_{LV}(R, \phi ) + \gamma_{SV} A_{SV}(R, \phi ) + \gamma_{LS} A_{LS}(R, \phi ), \end{equation}

where liquid (outside) and vapour (inside) are assumed to be in their respective bulk states and are separated by sharp interfaces from the neighbouring phases. In the above expression, the subscripts $L$, $V$ and $S$ refer to liquid, vapour and solid phases, respectively, ${\rm \Delta} p = p_V - p_L$ is the vapour–liquid pressure difference (the Laplace pressure), $V_V$ is the bubble volume, $A_{LV}$, $A_{SV}$ and $A_{LS}$ are the liquid–vapour, solid–vapour and liquid–solid interface areas, respectively, with $\gamma _{LV}$, $\gamma _{SV}$ and $\gamma _{LS}$ the corresponding surface energies. The other relevant parameter is the equilibrium (or Young) contact angle $\phi = \cos ^{-1}((\gamma _{SV} - \gamma _{LS})/\gamma _{LV})$ sketched in figure 1(a), which, according to the standard convention, is measured from the liquid–solid interface, i.e. $\phi < {\rm \pi}/2$ means hydrophilic. The contact angle allows the relevant geometric quantities to be re-expressed as $A_{SV} = {\rm \pi}R^2 \sin ^2 \phi$, $A_{LV} = 2 {\rm \pi}R^2 (1+ \cos \phi )$, $A_{LS} = A_w - A_{SV}$ and $V_V (R, \phi ) = V_V (R,{\rm \pi} ) \varPsi (\phi )$, where $A_w$ is the total area of the solid surface and $\varPsi (\phi )=\frac{1}{4} (1+ \cos \phi )^2 (2 - \cos \phi )$, representing a geometric factor that accounts for the contact angle. As $\phi \to 0$ the free energy reduces to the homogeneous case where the work required to form a bubble with radius $R$ is ${\rm \Delta} \varOmega _{hom} = -{\rm \Delta} p \,V_V (R, 0) + \gamma _{LV} A_{LV}(R, 0)$.

Figure 1. (a) Sketch of a vapour bubble on a flat solid surface, illustrating some geometrical parameters. (b) Energy landscapes of a vapour bubble as a function of the radius, at different contact angles $\phi$. The solid line shows the homogeneous case, as a reference. Both energy and radius are rescaled with their critical values, ${\rm \Delta} \varOmega ^*$ and $R^*$, respectively.

It should be noted that the energy required to form a spherical cap at the wall is the fraction $\varPsi (\phi )$ of that required to nucleate a bubble with the same radius in the bulk,

(2.2)\begin{equation} {\rm \Delta} \varOmega (R, \phi ) = {\rm \Delta} \varOmega_{hom} (R ) \varPsi(\phi). \end{equation}

Since the free energy consists of two contributions – a volume term scaling with $R^3$ and a surface term scaling with $R^2$ – a free energy maximum (critical state) is attained at the critical radius

(2.3)\begin{equation} R^*=\frac{2\gamma_{LV}}{{\rm \Delta} p}, \end{equation}

which is intrinsic to the fluid and independent of surface wettability. The corresponding free energy barrier, i.e. the free energy difference between the critical and pure liquid states, is

(2.4)\begin{equation} {\rm \Delta} \varOmega^*= {\rm \Delta} \varOmega (R^*,\phi ) = {\rm \Delta} \varOmega_{hom} ^* \varPsi(\phi ) = \frac{16}{3} {\rm \pi}\frac{\gamma_{LV}^3}{{\rm \Delta} p ^2} \varPsi(\phi). \end{equation}

At variance with the critical radius, which is the same under both homogeneous and heterogeneous conditions, irrespective of the contact angle, the energy barrier ${\rm \Delta} \varOmega ^*$ is lowered by the contact-angle-dependent factor $\varPsi (\phi ) \leq 1$ with respect to nucleation in the bulk. The free energy profiles ${\rm \Delta} \varOmega (R)$ are shown in figure 1(b) for two contact angles, hydrophilic and hydrophobic, respectively, and compared with the homogeneous case. The role of the surface in reducing the free energy barrier below ${\rm \Delta} \varOmega ^*_{hom}$ is apparent, and nucleation over a hydrophobic wall is favoured, as expected.

Together with the energy barrier, the other crucial quantity in nucleation processes is the nucleation rate, i.e. the normalized number of supercritical bubbles formed per unit time. In the heterogeneous context, the nucleation rate is normalized per unit surface (as opposed to unit volume, as used in homogeneous conditions). The classical expression for the nucleation rate is (Blander & Katz 1975; Debenedetti 1996)

(2.5)\begin{equation} J_{BK} = n^{2/3}_L \frac{(1+\cos\phi)}{2}\sqrt{\frac{2\gamma_{LV}}{{\rm \pi} m}} \exp\left(- \frac{{\rm \Delta} \varOmega^*}{k_B \theta}\right) , \end{equation}

where $n_L$ is the number density of liquid molecules and $m$ their mass.

2.1. Critical bubble in the grand canonical ensemble

After retracing the main features of the classical theory of the nucleation of vapour bubbles in a liquid, the analysis is now specialized for the grand canonical $\mu V \theta$ ensemble (with parameters the chemical potential $\mu$, system volume $V$ and temperature $\theta$). Here, at fixed volume, the two different phases (liquid and vapour) exhibit the same chemical potential and the same temperature. The analysis is explicitly conducted for the homogeneous case, the extension to the heterogeneous case being simply obtained using the geometrical function $\psi (\phi )$ that rescales the energy barrier and the critical volume. The grand canonical ensemble is particularly suitable for addressing vapour bubble nucleation, since it imposes only mechanical and thermodynamical equilibrium between phases, without enforcing macroscopic constraints. Within this setting, the liquid and vapour densities, $\rho _L$ and $\rho _V$, respectively, follow from the temperature and chemical potential, $\mu _L (\rho _L, \theta ) = \mu$ and $\mu _V (\rho _V, \theta ) = \mu$, and the critical vapour bubble radius is determined by (2.4).

2.2. Critical bubble in the microcanonical ensemble

The microcanonical $NVE$ ensemble requires, by definition, mass, $M$ (or equivalently number of particles $N$), volume $V$ and energy $E$ to be constrained. Clearly, when the system is large in comparison with the typical bubble size, the $\mu V \theta$ and the $N V E$ ensembles are equivalent. However, for confined and crowded systems with many bubbles, the nucleation processes can be different. In fact, this difference will turn out to be important to understand the numerical results in § 9, suggesting the need to review $NVE$ nucleation explicitly.

In the $NVE$ ensemble the equations to be solved to determine the critical bubble are more complicated than in the previous case:

(2.6)\begin{equation} \left.\begin{array}{c@{}} R^* = \dfrac{2\gamma_{LV}}{p_V (\rho_V, \theta) - p_L (\rho_L, \theta) }, \quad \mu_V(\rho_V, \theta ) = \mu_L (\rho_L, \theta) , \\ U_L (\rho_L, \theta) + U_V (\rho_V, \theta ) + E_c = E, \quad V_V \rho_V + (V -V_V)\rho_L = M . \end{array}\right\} \end{equation}

The unknowns here are the liquid and vapour densities, the temperature and the bubble radius. In the equation expressing the energy constraint, $E_c$ is the interfacial (capillary) energy, which in the case of a single nucleated bubble is given by $\gamma _{LV} 4 {\rm \pi}R^{*2}$. Here, it is left as an additional parameter for reasons to be discussed in § 9. The equations of state for the internal energy of the liquid and vapour $U_{L/V}$ and the corresponding chemical potentials $\mu _{L/V}$ and pressures $p_{L/V}$ are assumed to be given.

The ‘free energy landscape’, for different levels of confinement, is reported in figure 2. The thermodynamic potential of the $NVE$ ensemble, i.e. the entropy, is plotted as a function of the bubble radius, taken as reaction coordinate for the process (stability corresponds to entropy maximum). The overall picture corresponds quite well to the phenomenology observed in the $\mu V \theta$ system: at extreme confinement the liquid remains stable, similar to what is found in Marti et al. (2012) and Vincent & Marmottant (2017). When the available volume is large enough, a barrier separates the metastable liquid from the equilibrium state featuring a bubble. The radius of the equilibrium bubble is also reported in figure 2 as a function of confinement. The equilibrium radius increases monotonically, consistently with the familiar result that full vapour is the equilibrium state in the ${\mu V \theta }$ ensemble. The critical radii in the $NVE$ and ${\mu V \theta }$ ensembles are substantially identical as long as the confinement is not extreme, as in the simulations to be discussed in § 9.

Figure 2. Main panel: The $NVE$ CNT for homogeneous fluids, showing entropy profiles ${\rm \Delta} S \,\theta _0/{\rm \Delta} \varOmega ^*$ (normalized with the $\mu V \theta$ free energy (grand potential) barrier ${\rm \Delta} \varOmega ^*$ at $\theta = \theta _0$) versus bubble radius $R$ (normalized with the $\mu V \theta$ critical radius $R^*$) at changing confinement $V/V^*$, where $V^* = \frac 43 {\rm \pi}R^{*3}$ is the critical bubble volume. For volumes below $V/V^* = 481$ no critical bubble exists: confinement is strong enough to stabilize the liquid. Above this limiting confinement, a barrier separates the metastable liquid from the stable configuration featuring a finite size bubble (maximum in the entropy, explicitly appearing in the plot only for relatively small confinement volumes). Please note that with entropy as thermodynamic potential, the sign of stable and unstable extrema are reversed with respect to the more usual cases involving free energies. Inset: Radius of the equilibrium bubble versus confinement (log scale). As the confinement is reduced (available volume increased), the size of the equilibrium bubble grows larger, to eventually grow unbounded for infinite volume to meet the $\mu V \theta$ prediction. Note that the curve in the inset starts at a finite (small) $V/V^*$.

3. Thermodynamics of liquid–vapour systems at solid walls

The purpose of the present section is modelling hydrophilic/hydrophobic walls in the context of a capillary fluid à la van der Waals. In particular, the surface free energy density at the wall is explicitly determined consistently with the thermodynamics of the bulk fluid and the prescribed equilibrium contact angle. To this end, the van der Waals square-gradient approximation of the (Helmholtz) free energy functional (van der Waals 1979) is extended as

(3.1)\begin{equation} F[\rho,\theta] = \int_{V} f(\rho, \boldsymbol{\nabla}\rho, \theta)\,\textrm{d} V + \oint_{\partial V} f_w (\rho, \theta)\,\textrm{d} S, \end{equation}

where $f = f_b + f_c$, with $f_b$ the free energy density (per unit volume) of the bulk fluid at mass density $\rho$ and temperature $\theta$ and $f_c = (\lambda /2)\boldsymbol {\nabla }\rho \boldsymbol {\cdot } \boldsymbol {\nabla }\rho$ the capillary contribution. As detailed below, in equilibrium, this free energy model describes a smooth density profile transitioning between liquid and vapour on a typical scale $\epsilon$, the thickness of the interface. As discussed in Anderson, McFadden & Wheeler (1998) and Lutsko (2011), this corresponds to a gradient approximation of more general, non-local descriptions, like those exploited in DFT (Tarazona & Evans 1984). The free energy contribution $f_w$ arises from the fluid–wall interactions and accounts for the wetting properties of the surface. Its explicit expression, provided by Cahn (1977) in the context of the isothermal phase field model for immiscible liquids, is extended here to the non-isothermal van der Waals square-gradient model relevant for liquid–vapour phase transitions.

It is instrumental to introduce the entropy functional $S$, obtained as the functional derivative of the free energy with respect to the temperature,

(3.2)\begin{align} S[\rho,\theta]&= \int_{V} - \frac{\delta F}{\delta \theta}\,\textrm{d} V = - \int_{V} \frac{\partial f_b}{\partial \theta} \,\textrm{d} V - \oint_{\partial V} \frac{\partial f_w}{\partial \theta}\,\textrm{d} S \nonumber\\ &= \int_{V} s_b(\rho,\theta) \,\textrm{d} V + \oint_{\partial V} s_w(\rho, \theta) \,\textrm{d} S , \end{align}

where the third equality holds for a temperature-independent $\lambda$ (assumed here for the only purpose of simplifying the exposition). The last identity follows from the definition of the bulk entropy density $s_b$ and after the introduction of the surface entropy density $s_w$.

For a closed and isolated thermodynamic system of given energy $E_0$ and mass $M_0$, the constrained entropy functional ($S_c$) reads

(3.3)\begin{equation} S_c = S + l_1 \left(M_0 - \int_{V} \rho\, \textrm{d} V \right) + l_2(E_0 - U), \end{equation}

where $l_1$ and $l_2$ are the Lagrange multipliers enforcing mass and energy constraints. The internal energy $U$ is

(3.4)\begin{align} U& = F + \theta S = \int_{V} ( u_b(\rho,\theta) + u_c (\boldsymbol{\nabla}\rho) ) \, \textrm{d} V + U_w [\rho, \theta ]\nonumber\\ &= \int_{V} \left(u_b(\rho,\theta) +\tfrac12\lambda\boldsymbol{\nabla}\rho\boldsymbol{\cdot}\boldsymbol{\nabla} \rho\right) \textrm{d} V + \oint_{\partial V} (\,f_w + \theta s_w) \, \textrm{d} S , \end{align}

with $u_b = f_b - \theta \partial f_b/\partial \theta$ and $u_c =(\lambda /2) \boldsymbol {\nabla }\rho \boldsymbol {\cdot }\boldsymbol {\nabla }\rho$. By maximizing the constrained entropy,

(3.5)\begin{align} \delta S_c [ \rho,\theta] &= \delta \int_{V}\, (s_b - l_2u (\rho, \boldsymbol{\nabla}\rho,\theta) - l_1 \rho ) \,\textrm{d} V \nonumber\\ &\quad + \delta \oint_{\partial V} [s_w - l_2 (\,f_w + \theta s_w )]\,\textrm{d} S = 0, \end{align}

the Lagrange multipliers are identified as $l_1 = - (\mu ^b_c - \lambda \nabla ^2 \rho )/ \theta$ and $l_2= 1/\theta$, where $\mu ^b_c = \partial f_b/\partial \rho$ is the bulk chemical potential. It follows that, at equilibrium, the temperature and the (generalized) chemical potential $\mu _c$ must be constant:

(3.6)\begin{gather} \theta = \textrm{const.} = \theta_{eq}, \end{gather}

(3.7)\begin{gather}\mu_c = \mu^b_c - \lambda \nabla^2 \rho= \textrm{const.}= \mu^{eq}_c. \end{gather}

The boundary term produces the additional requirement

(3.8)\begin{equation} \left.\left(\lambda\boldsymbol{\nabla}\rho \boldsymbol{\cdot} \hat{\boldsymbol{n}} + \frac{\partial f_w}{\partial \rho}\right) \right\vert_{\partial V} =0 , \end{equation}

where $\hat {\boldsymbol{n}}$ is the outward normal, to be read as a (nonlinear) boundary condition for the density.

The above equilibrium conditions provide the relationship between density distribution and system thermodynamics. It is first illustrated for a flat interface away from solid walls (and constant $\lambda$) to set the stage for the discussion of the boundary condition at the solid surface. For a flat interface the inhomogeneous direction is $\hat {\boldsymbol s}$ (i.e. $\rho = \rho (s)$) and the equilibrium condition (3.7) reads

(3.9)\begin{equation} \mu_c = \mu_c^b(\rho, \theta)-\lambda \frac{\textrm{d}^2\rho}{\textrm{d} s^2} = \mu_{eq}. \end{equation}

After multiplying (3.9) by $\textrm {d}\rho /\textrm {d} s$ and integrating between $\rho _\infty =\rho _V$ (the saturation vapour density) and $\rho$, one has

(3.10)\begin{equation} {w}_b(\rho,\theta) - { w}_b(\rho_V(\theta)) = \frac{\lambda}{2}\left(\frac{\textrm{d}\rho}{\textrm{d} s}\right)^2, \end{equation}

with ${w}_b = f_b - \mu _{eq}\rho$ the bulk Landau free energy density (grand potential). The grand potential of this unbounded inhomogeneous system is defined, as usual, as the Legendre transform of the free energy,

(3.11)\begin{equation} \varOmega = F - \int_{V} \rho \frac{\delta F}{\delta \rho} \, \textrm{d} V = \int_{V} {w} \, \textrm{d} V, \end{equation}

where $w=f_b + (\lambda /2) \boldsymbol {\nabla } \rho \boldsymbol {\cdot } \boldsymbol {\nabla } \rho - \mu _c \rho$ is the grand potential density (it is worth stressing that $w_b$ and $w$ are different, unless the system is homogeneous, the latter including interfacial effects). The surface tension is the excess grand potential density is

(3.12)\begin{align} \gamma_{LV} &= \int_{-\infty}^{s_i} ( {w}[\rho,\theta] - {w}[\rho_V] )\,\textrm{d} s \nonumber\\ & \quad + \int^{\infty}_{s_i} ( {w}[\rho,\theta] - {w}[\rho_L] ) \, \textrm{d} s =\int_{-\infty}^{\infty} ( {w}[\rho,\theta] - {w}[\rho_V] ) \, \textrm{d} s , \end{align}

where $s_i$ denotes the position of the Gibbs dividing surface, and the equilibrium property $w[\rho _V] = w[\rho _L]$ has been used. The definition of ${ w}[\rho ]$, (3.11) and the equilibrium condition (3.9) provide

(3.13)\begin{align} \gamma_{LV} &= \int_{-\infty}^{\infty} \left[ { f}_b + \frac{1}{2} \lambda \left(\frac{\textrm{d} \rho}{\textrm{d} s} \right)^2 - \mu_{eq} \rho - {w}_b(\rho_V) \right] \textrm{d} s \nonumber\\ &= \int_{-\infty}^{\infty} \left[ {w}_b + \frac{1}{2} \lambda \left(\frac{\textrm{d} \rho}{\textrm{d} s} \right)^2 - { w}_b(\rho_V) \right] \textrm{d} s, \end{align}

which, combined with (3.10), yields

(3.14)\begin{equation} \gamma_{LV} = \int_{-\infty}^{+\infty}\lambda \left(\frac{\textrm{d}\rho}{\textrm{d} s}\right)^2\, \textrm{d} s = \int_{\rho_V}^{\rho_L} \sqrt{2\lambda ({w}_b(\rho,\theta) - {w}_b(\rho_V))}\,\textrm{d} \rho.\end{equation}

It should be noted that, besides the capillary coefficient $\lambda$, the surface tension depends only on the bulk grand potential density in the coexistence mass density range $\rho _V(\theta ) \le \rho \le \rho _L(\theta )$.

The model is complete after determining the explicit expression for the fluid–surface interaction term $f_w$. The resulting expression will generalize the one proposed by Cahn (1977) for two immiscible fluids of equal density that has been exploited in several phase-field-based numerical simulations (see e.g. Jacqmin 2000; Yue, Zhou & Feng 2010; Sartori et al. 2015). The starting point to derive the proper expression for $f_w$ is the geometrical relation $\hat {\boldsymbol s}\boldsymbol {\cdot } \hat {\boldsymbol n} = \cos \phi$, involving the (contact) angle between the wall-normal and interface-normal directions (see figure 3), and (3.8), which is rearranged as

(3.15)\begin{equation} \frac{\textrm{d} f_w}{\textrm{d} \rho} + \lambda \frac{\textrm{d} \rho}{\textrm{d} s} \cos \phi = 0 , \end{equation}

using $\boldsymbol {\nabla }\rho = (\textrm {d} \rho /\textrm {d} s ) \hat {\boldsymbol s}$.

Figure 3. Sketch of the capillary forces originating at the triple contact line. The vector normal to the wall is indicated with $\hat {\boldsymbol {n}}$, while $\hat {\boldsymbol {s}}$ indicates the vector normal to the interface, in the direction of the density gradient. The geometrical relation $\hat {\boldsymbol s}\boldsymbol {\cdot }\hat {\boldsymbol n} = \cos \phi$ holds between the vectors and the contact angle $\phi$.

At equilibrium the interface-normal density variation follows from (3.10), allowing one to integrate the above equation as

(3.16)\begin{equation} f_w(\rho,\theta) = -\cos \phi \int_{\rho_V}^\rho \sqrt{2\lambda (w_b (\tilde{\rho},\theta ) - w_b (\rho_V ) )} \,\textrm{d} \tilde{\rho} + f_w (\rho_V ) . \end{equation}

This form of $f_w$ is consistent with the surface free energy of a pure vapour with density $\rho _V$ in contact with the wall, given by the solid–vapour surface tension, $f_w(\rho _V) = \gamma _{SV}$. Similarly, for a pure liquid with density $\rho _L$, (3.16) yields

(3.17)\begin{equation} f_w(\rho_L) = \gamma_{LS} = -\cos \phi \int_{\rho_V}^{\rho_L} \sqrt{2\lambda(w_b(\tilde{\rho},\theta) - w_b(\rho_V))} \,\textrm{d} \tilde{\rho} + \gamma_{SV} = -\gamma_{LV}\cos \phi + \gamma_{SV}, \end{equation}

which is the well-known Young equation for the equilibrium contact angle.

The expression (3.16) for the surface energy differs from the usual recipes found in the literature. Seppecher (1996), for example, proposed a linear dependence on the density, while Sibley et al. (2013) considered a third-order polynomial (see also Shen et al. 2017). The cubic ansatz is particularly tricky, since it turns out to be consistent with the present equation (3.16) for the special case of a quartic, Ginzburg–Landau, double-well potential $\omega _b(\rho ) = \alpha (\rho -\rho _L)^2(\rho -\rho _V)^2$, which is the most common free energy model for the Cahn–Hilliard equation (Anderson et al. 1998; Jacqmin 2000; Carlson, Do-Quang & Amberg 2010; Magaletti et al. 2013), as also noted in Laurila et al. (2012), where the authors comment that for a van der Waals fluid the third-order polynomial is just an approximation. Equation (3.16) suggests that $f_w$ cannot in general be assigned independently of fluid bulk free energy $f_b(\rho )$, a prescription that is often ignored in applications. In figure 4 we show the equilibrium configurations of a two-dimensional vapour bubble laying on a flat wall with two different wettabilities. In both cases the contact angle $\phi$ has been imposed through the boundary condition for the density field, (3.8). The critical density isoline forms the prescribed contact angle with the solid wall, as expected.

Figure 4. (a) Adsorption/depletion of the density field occurring in the wall-normal direction, $z$, as a function of the contact angle. In the main panel the bulk liquid density is $\rho _b=0.48$. All the quantities are expressed in Lennard-Jones units, as will be explained in § 9. Hydrophobic walls, $\phi > 90^\circ$, show a reduction of the density close to the wall; the opposite behaviour is provided by hydrophilic walls. The inset shows that this layering effect is amplified when the degree of metastability is increased. (b) Vapour–liquid contact angle for hydrophobic (top) and hydrophilic (bottom) solid surfaces.

The proposed fluid–solid free energy describes density layering at the solid surface for non-neutrally wettable surfaces ($\cos \phi \ne 0$). Layering is a common feature for liquids in contact with solid walls. At nanoscale, oscillations of the density field close to the wall are known to occur, as observed through X-ray scattering experiments (Huisman et al. 1997; Yu et al. 2000) or MD simulations (Sikkenk et al. 1987). They are well captured by the intrinsically non-local DFT (Tarazona & Evans 1984; Tarazona, Marini Bettolo Marconi & Evans 1987; Evans, Stewart & Wilding 2017). The possible effect on nucleation of these near-wall density oscillations cannot be addressed in the present framework, which, due to the intrinsic limitations of purely local theories, should be understood as a coarse-grained description of the actual phenomenology, resulting in a monotonic density layering at the wall. A related approach is described in (Carey & Wemhoff 2005), where the fluid–solid interaction is accounted for through an interaction potential between solid and fluid particles, leading, through a mean-field theory, to a disjoining pressure, which induces a similar monotonic density stratification at the wall.

In the hydrophilic case, the liquid–solid interaction accumulates fluid at the wall, resulting in a local increase of density. Conversely, the lower affinity of hydrophobic walls produces a depletion region with layers of the order of a few nanometres – see figure 4 reporting the density profiles obtained by applying the boundary condition (3.8) for different wetting properties of the surface. This aspect is important for heterogeneous nucleation, since it facilitates the process at hydrophobic walls while discouraging it on hydrophilic surfaces (see the more detailed discussion in § 9).

4. The capillary Navier–Stokes model

Aim of the present section is deriving constitutive relations for stresses and energy fluxes consistent with the Clausius–Duhem formulation of thermodynamic irreversibility (De Groot & Mazur 2013) for non-homogeneous, wall-bounded capillary systems – see Magaletti et al. (2016) for the case of a homogeneous fluid.

The theory explicitly takes into account capillary effects occurring in the fluid, not only at the vapour–liquid interface but also in the stratified layer close to a substrate (typically, a solid wall). Stratification at a substrate is induced by the combination of capillarity (the $\lambda |\boldsymbol {\nabla } \rho |^2$ contribution to the free energy (3.1)) and surface free energy $f_w(\rho ,\theta )$ included in the boundary contribution. Associated with the surface free energy term there is a surface entropy with density $s_w = - \partial f_w/\partial \theta$. In order to account for possible additional entropy production terms arising from the wall potential, the system can be thought of as comprising fluid, solid and a concentrated, zero-thickness layer separating fluid and solid. In total generality the layer may possess a (concentrated) mass density $\rho _w$ (units of mass per unit area) and a velocity $\boldsymbol {v}_w$, hence a kinetic energy density $\frac 12 \rho _L \vert \boldsymbol {v}_w \vert ^2$, and a total surface energy density, the sum of internal, $u_w = f_w + \theta \partial f_w/\partial \theta$, and kinetic energy density. The surface free energy density is assumed to depend on the mass density of the fluid in contact with the layer $\rho$, on the concentrated mass density associated with the layer $\rho _w$ and on temperature. The equations for the layer establish mass, momentum and energy conservation, and will be used in the limit of vanishing layer mass density, $\rho _w$, such that, in the limit, the surface potential recovers the form $f_w(\rho ,\theta )$. In order not to burden the discussion, the governing equations inside the layer will be described in detail in appendix A. Here, only the results of the procedure are reported, consisting of a relaxation condition for the contact angle in the diffuse capillary context.

For a capillary fluid enclosed in a volume ${\mathcal {D}}$, the conservation laws for mass $M$, momentum $\boldsymbol {P}$ and total energy $E$ are (the inclusion of mass forces is straightforward)

(4.1)\begin{equation} \left.\begin{array}{c@{}} \displaystyle\dfrac{\textrm{d} M}{\textrm{d} t} = \dfrac{\textrm{d} }{\textrm{d} t}\int_{\mathcal{D}} \rho \, \textrm{d} V= 0 , \\ \displaystyle\dfrac{\textrm{d} \boldsymbol{P}}{\textrm{d} t} = \dfrac{\textrm{d} }{\textrm{d} t}\int_{\mathcal{D}} \rho\boldsymbol{v} \, \textrm{d} V = \int_{\partial {\mathcal{D}}} \boldsymbol{\varSigma}\boldsymbol{\cdot} \boldsymbol{n} \, \textrm{d} S , \\ \displaystyle\dfrac{\textrm{d} E}{\textrm{d} t} = \dfrac{\textrm{d} }{\textrm{d} t}\int_{\mathcal{D}} e \, \textrm{d} V = \int_{\partial {\mathcal{D}}} (\boldsymbol{\varSigma} \boldsymbol{\cdot} \boldsymbol{v}-\boldsymbol{q})\boldsymbol{\cdot} \boldsymbol{n} \, \textrm{d} S , \end{array}\right\} \end{equation}

where $\boldsymbol {n}$ is the outward unit vector normal to the boundary, $\boldsymbol {v}$ is the fluid velocity, $e=u+\frac 12 \rho \vert \boldsymbol {v}\vert ^2$ is the total energy density, with $u=u_b+u_c$ the internal energy density defined in (3.4), and $\boldsymbol {\varSigma }$ and $\boldsymbol {q}$ the are stress tensor and heat flux, respectively, to be determined in the following. The local form of the conservation laws is obtained by applying Green's theorem, and take the classical form

(4.2)\begin{equation} \left.\begin{array}{c@{}} \displaystyle \dfrac{\partial \rho}{\partial t}+ \boldsymbol{\nabla}\boldsymbol{\cdot}\left( \rho\boldsymbol{v}\right) = 0 , \\ \displaystyle \dfrac{\partial \rho\boldsymbol{v}}{\partial t}+ \boldsymbol{\nabla}\boldsymbol{\cdot}\left( \rho\boldsymbol{v}\otimes\boldsymbol{v}\right) = \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\varSigma} ,\\ \displaystyle \dfrac{\partial e}{\partial t}+ \boldsymbol{\nabla}\boldsymbol{\cdot}( \boldsymbol{v} e) = \boldsymbol{\nabla}\boldsymbol{\cdot}( \boldsymbol{\varSigma}\boldsymbol{\cdot}\boldsymbol{v} -\boldsymbol{q}) . \end{array}\right\} \end{equation}

Under the assumption of local equilibrium, the above equations can be manipulated, as usual in the context of non-equilibrium thermodynamics (De Groot & Mazur 2013), to obtain an evolution equation for the entropy density. Namely, the internal energy evolution is first derived by subtracting the kinetic energy contribution from the total energy density,

(4.3)\begin{equation} \rho\frac{\textrm{D} \tilde{u}}{D t } - \boldsymbol{\varSigma} \boldsymbol{:} \boldsymbol{\nabla} \boldsymbol{v} + \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{q} = 0, \end{equation}

where $\textrm {D}/\textrm {D}t = \partial /\partial t + \boldsymbol {v}\boldsymbol {\cdot } \boldsymbol {\nabla }$ is the material derivative and $\tilde {u} = u/\rho$ is the specific internal energy. By definition $\tilde {u} = \tilde {f} + \theta \tilde {s}$, with $\tilde {f} = (\,f_b + (\lambda /2)\vert \boldsymbol {\nabla }\rho \vert ^2)/\rho$ – see (3.1) and comments below – and hence its differential reads

(4.4)\begin{equation} {\mathrm d} \tilde{u} = \dfrac{\partial \tilde{f}}{\partial \rho}{\mathrm d}\rho + \dfrac{\partial \tilde{f}}{\partial \boldsymbol{\nabla}\rho}\boldsymbol{\cdot}{\mathrm d}\boldsymbol{\nabla}\rho + \theta\, {\mathrm d}\tilde{s}. \end{equation}

The partial derivatives of the specific free energy are derived from its definition, (3.1), providing

(4.5)\begin{equation} \rho \frac{\textrm{D} \tilde{u}}{\textrm{D} t} = \frac{1}{\rho}\left(p-\frac{\lambda}{2}\vert\boldsymbol{\nabla}\rho\vert^2\right) \frac{\textrm{D} \rho}{\textrm{D} t} + \lambda\boldsymbol{\nabla}\rho\boldsymbol{\cdot} \frac{\textrm{D} \boldsymbol{\nabla}\rho}{\textrm{D}t} + \rho\theta\frac{\textrm{D} \tilde{s}_b}{\textrm{D} t} , \end{equation}

with $p$ the pressure. The material derivative of the density gradient (second term on the right-hand side of (4.5)) is evaluated from the mass conservation equation (Magaletti et al. 2016),

(4.6)\begin{equation} \lambda \boldsymbol{\nabla}\rho\boldsymbol{\cdot} \frac{\textrm{D} \boldsymbol{\nabla}\rho}{\textrm{D}t} = \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\lambda \boldsymbol{\nabla}\rho \frac{\textrm{D} \rho}{\textrm{D}t}\right) - \lambda \nabla^2 \rho \boldsymbol{\nabla}\rho \frac{\textrm{D} \rho}{\textrm{D}t} - \lambda \boldsymbol{\nabla}\rho \otimes \boldsymbol{\nabla}\rho \boldsymbol{:} \boldsymbol{\nabla}\boldsymbol{v}. \end{equation}

After substitution of (4.5) and (4.6) into (4.3), a few more elementary steps isolate the entropy time derivative as

(4.7)\begin{align} \rho \frac{\textrm{D} \tilde{s}}{\textrm{D} t} &= -\boldsymbol{\nabla}\boldsymbol{\cdot}\left[\frac{1}{\theta} \left(\lambda\boldsymbol{\nabla}\rho \frac{\textrm{D} \rho}{\textrm{D}t} + \boldsymbol{q}\right) \right] - \left[ \lambda \boldsymbol{\nabla}\rho \frac{\textrm{D} \rho}{\textrm{D}t}+ \boldsymbol{q}\right] \boldsymbol{\cdot}\frac{\boldsymbol{\nabla}\theta}{\theta^2} \nonumber\\ &\quad + \left[\boldsymbol{\varSigma} + \left(p - \frac{\lambda}{2}\lvert \boldsymbol{\nabla}\rho\rvert^2 - \lambda \rho\nabla^2 \rho \right) {\boldsymbol{\mathsf{I}}} + \lambda \boldsymbol{\nabla} \rho\otimes\boldsymbol{\nabla}\rho\right]\boldsymbol{:} \frac{\boldsymbol{\nabla}\boldsymbol{v}}{\theta} , \end{align}

where the divergence of the entropy flux (first term on the right-hand side) and two entropy production terms appear. By requiring a positive entropy production (Clausius–Duhem inequality) and by restricting the analysis to linear constitutive prescriptions for stress

(4.8)\begin{align} \boldsymbol{\varSigma} &= \boldsymbol{\varSigma}_{rev} + \boldsymbol{\varSigma}_{v} \nonumber\\ &= \left(-p + \frac{\lambda}{2}\vert\boldsymbol{\nabla}\rho\vert^2 + \lambda\rho\nabla^2\rho\right) {\boldsymbol{\mathsf{I}}} - \lambda\boldsymbol{\nabla}\rho\otimes\boldsymbol{\nabla}\rho + \eta\left[(\boldsymbol{\nabla}\boldsymbol{u} + \boldsymbol{\nabla}\boldsymbol{u}^\textrm{T}) - \frac{2}{3}\boldsymbol{\nabla} \boldsymbol{\cdot}{u}{\boldsymbol{\mathsf{I}}}\right] \end{align}

and energy flux

(4.9)\begin{equation} \boldsymbol{q} = \boldsymbol{q}_c + \boldsymbol{q}_h = - \lambda \boldsymbol{\nabla}\rho\frac{\textrm{D} \rho}{\textrm{D}t} - k \boldsymbol{\nabla}\theta, \end{equation}

capillary effects are explicitly included in the model. The irreversible part of the stress is the classical Newtonian viscous tensor, where $\eta >0$ is the dynamic viscosity. Since viscosity is not expected to play a crucial role in nucleation (e.g. in CNT it is completely neglected), the simplest choice $-2 \eta /3$ is made here for the second viscosity coefficient. A general value can be, however, straightforwardly included in the model as needed, e.g. to describe damping in nonlinear bubble oscillations (Scognamiglio et al. 2018). Similarly, the classical Fourier's heat flux, with $k>0$ the thermal conductivity, provides the irreversible part of the energy flux. It is worth noting that the capillary (square density gradient) term in the free energy leads to capillary stresses included in the reversible part of the dynamics (no entropy production), as expected from the general thermodynamic definition of surface tension.

The field equations for a capillary fluid, (4.2), (4.8) and (4.9), require an additional boundary condition for the density field. The equilibrium considerations of § 3 suggest the density normal derivative as the appropriate prescription, (3.8). This result can be extended to non-equilibrium conditions through a procedure strictly similar to that just reviewed for the bulk, applied in this case to the layer adjoining the fluid and solid surface. The layer is described in terms of suitable fields defined on the manifold supporting the layer that are energetically, mechanically and kinematically coupled with the neighbouring fluid. By prescribing that each entropy production term is non-negative (Clausius–Duhem inequality combined with Curie principle), linear constitutive laws for the layer stress tensor ${\boldsymbol{\mathsf{Q}}}_{\boldsymbol\pi}$ and for the layer tangential heat flux $\boldsymbol {q}_w$ are derived (see appendix A for derivation details).

The requirement that the entropy production associated with the layer should, in general, be non-negative leads to relaxing the equilibrium wetting condition. In this more general setting, the prescription of the equilibrium (Young) contact angle, (3.8), is replaced by the boundary equation

(4.10)\begin{equation} \dfrac{\textrm{D} \rho}{\textrm{D} t} = - D_w\left( \frac{\partial f_w}{\partial \rho} + \lambda \dfrac{\partial \rho}{\partial n} \right) , \end{equation}

where $D_w$ is a mobility coefficient determining the relaxation time scale. Notice that the surface free energy $f_w$, as an equilibrium property of the system, keeps the form (3.16) obtained in § 3. This equation describes the relaxation of the contact angle towards equilibrium and can be understood as a generalization to the non-isothermal case and to general equations of state for the bulk fluid of the contact angle dynamics proposed in Jacqmin (2000), Carlson, Do-Quang & Amberg (2009, 2011) and Ren & E (2011). In what follows, the boundary term associated with (4.10) will be assumed to provide a negligible contribution to the entropy, implying the equilibrium Young's law (3.15).

5. Equilibrium thermal fluctuations of a wall-bounded capillary fluid

Fluids at a mesoscopic scale are subject to thermal fluctuations that need to be included in the hydrodynamic equations. The purpose of this section is discussing thermal fluctuations for a capillary fluid in contact with a solid surface. What follows falls within the theoretical framework of fluctuating hydrodynamics, first proposed by Landau & Lifshitz (1980 reprint) and systematically developed since then (Fox & Uhlenbeck 1970; Zubarev & Morozov 1983; Español et al. 2009). The theory was retraced and exploited in the context of liquid–vapour systems far from boundaries to address homogeneous bubble nucleation in a recent work (Gallo et al. 2018). The extension of this theory to heterogeneous systems bounded by solid walls is presented here. The focus here is on a fluid enclosed between two fixed rigid solid walls to provide the static correlation functions of the fluid-phase fluctuations. In comparison with the unbounded cases, the interaction with the solid influences the fluctuating fluid field through the boundary conditions.

Equilibrium thermal fluctuations were originally described by Einstein (1956), who provided the static correlation functions as a result of the entropy deviations ${\rm \Delta} S$ from equilibrium, where the entropy $S_{eq}$ is a maximum. For the system described in § 4, the entropy consists of two contributions: the fluid, $S_f$, and the layer entropy, $S_w$, respectively. In § 3 the entropy deviation was expressed as a functional of mass density, $\delta \rho (\boldsymbol x,t)$, velocity, $\delta {\boldsymbol v}(\boldsymbol x,t)$, and temperature, $\delta \theta (\boldsymbol x,t)$, fluctuations (see (3.3)).

Assuming small fluctuations with respect to the mean field, the entropy – a functional of the mass density, velocity and temperature fields – can be Taylor-expanded to second order around equilibrium (second order at least is required, since at equilibrium the first variation vanishes altogether). Explicitly, using the variables $\boldsymbol {U} = (\rho , \boldsymbol {\nabla }\rho , \theta , \boldsymbol {v})^\textrm {T}$, where field derivatives higher than the density gradient will not be needed since the free energy (3.1) is quadratic in the density gradients, the entropy functional is

(5.1)\begin{align} {\rm \Delta} S_c &= S_c -S_{eq} = \int_{V}{\rm \Delta} s_{c f} (\rho, \boldsymbol{\nabla}\rho, \theta, \boldsymbol v)\, \textrm{d} V + \oint_{\partial V}{\rm \Delta} s_{c w} (\rho, \theta, \boldsymbol v)\, \textrm{d} S \nonumber\\ &= \left.\int_{V} \frac{1}{2}\,\sum_{i,j} \frac{\partial^2 {\rm \Delta} s_c }{\partial U_i \partial U_j}\right\vert_{eq} \delta U_i \delta U_j \, \textrm{d} V + \left.\oint_{\partial V} \frac{1}{2}\,\sum_{i,j} \frac{\partial^2 {\rm \Delta} s_{c w} }{\partial U_i \partial U_j}\right\vert_{eq} \delta U_i \delta U_j \, \textrm{d} S \nonumber\\ &= {\rm \Delta} S_{c f} + {\rm \Delta} S_{c w} . \end{align}

Exploiting the well-known relations (the subscript b stands for homogeneous bulk conditions where capillary contributions do not appear)

(5.2)\begin{equation} \left.\begin{array}{c@{}} \displaystyle \mathrm{d} s_b = \dfrac{1}{\theta}\,\mathrm{d} u_b - \dfrac{\mu_b}{\theta}\,\mathrm{d}\rho ,\\ \displaystyle \mathrm{d} u_b = \rho c_{v}\,\mathrm{d}\theta+\left(\mu_b+\theta\left. \dfrac{\partial s_b}{\partial\rho}\right\vert_\theta\right)\, \mathrm{d}\rho, \\ \displaystyle \mathrm{d}\mu_b = \dfrac{c_{T}^{2}}{\rho}\,\mathrm{d}\rho+\left(\dfrac{1}{\rho} \left.\dfrac{\partial p}{\partial\theta}\right\vert_\rho-\dfrac{s_b}{\rho}\right)\, \mathrm{d}\theta , \end{array}\right\} \end{equation}

where $c_{v}$ is the constant-volume specific heat and $c_{T}$ is the isothermal speed of sound, (5.1) can be rearranged as

(5.3)\begin{align} {\rm \Delta} S_c &= -\frac{1}{2} \int_V \left[ \frac{c_{T, eq}^{2}}{\theta_{eq}\rho_{eq}}\delta\rho^{2}-\frac{\lambda}{\theta_{eq}} \delta\rho (\nabla^2\delta\rho ) \right] \textrm{d} V \nonumber\\ &\quad - \frac{1}{2} \int_V \left[ \frac{\rho_{eq}}{\theta_{eq}}\delta\boldsymbol{v}\boldsymbol{\cdot} \delta\boldsymbol{v}+\frac{\rho_{eq}c_{v,eq}}{\theta_{eq}^{2}} \delta\theta^{2}\right] \textrm{d} V \nonumber\\ &\quad - \frac{1}{2} \oint_{\partial V} \frac{1}{\theta_{eq}}\left[ \lambda \delta\left( \frac{\partial \rho}{\partial n}\right)_{eq} \delta\rho + \left.\frac{\partial^2 f_w}{\partial \rho^2}\right\vert_{eq} \delta\rho^2\right] \textrm{d} S - \frac{1}{2} \oint_{\partial V} \frac{1}{\theta_{eq}} \left.\frac{\partial s_w}{\partial \theta}\right\vert_{eq} \delta\theta^2 \, \textrm{d} S . \end{align}

The first two rows are the volume contributions already discussed in

Gallo et al. (2018). The new terms on the third line arise from the solid–liquid interaction due to the layer entropy. The first integral in the third row of (5.3) has two different contributions: the first one arises from the integration by parts of the square gradient, and the second comes from the solid–liquid energy $f_w$. As a consequence of the static boundary condition, (3.8), the whole first integral vanishes (see the comments at the end of § 4 and appendix A):

(5.4)\begin{equation} \oint_{\partial V} \lambda \delta\left( \frac{\partial \rho}{\partial n}\right) \delta\rho + \frac{\partial^2 f_w}{\partial \rho^2} \delta\rho^2 \, \textrm{d} S = \oint_{\partial V} \frac{\delta}{\delta \rho} \left[\left(\lambda\boldsymbol{\nabla}\rho \boldsymbol{\cdot} \hat{\boldsymbol{n}} + \frac{\partial f_w}{\partial \rho}\right)\right] \delta\rho \, \textrm{d} S =0 . \end{equation}

Fluctuating field correlations are most easily obtained after recasting the entropy in a form more appropriate to operational manipulations:

(5.5)\begin{align} {\rm \Delta} S_c & = -\frac{1}{2} \int_{V} \int_{V} \textrm{d} V_{\boldsymbol x} \, \textrm{d} V_{\tilde {\boldsymbol x}} \left\{\delta {\boldsymbol v}({\boldsymbol x}) \frac{\rho_{eq}}{\theta_{eq}} \delta({\boldsymbol x} - {\tilde {\boldsymbol x}} ) \boldsymbol{\cdot} \delta {\boldsymbol v}({\tilde {\boldsymbol x}}) \right. \nonumber\\ &\quad +\delta\rho({\boldsymbol x}) \left[\frac{c_{T,eq}^{2}}{\theta_{eq}\rho_{eq}} \delta({\boldsymbol x} - {\tilde {\boldsymbol x}} ) -\frac{\lambda}{\theta_{eq}} \nabla_{\boldsymbol x}^2 \delta({\boldsymbol x} - {\tilde {\boldsymbol x}} ) \right] \delta\rho({\tilde {\boldsymbol x}}) \nonumber\\ &\quad \left.+ \, \delta \theta({\boldsymbol x}) \frac{\rho_{eq} c_{v,eq}}{\theta_0^2} \delta({\boldsymbol x} - {\tilde {\boldsymbol x}} ) \delta \theta({\tilde {\boldsymbol x}}) \right\} \nonumber\\ &\quad -\frac{1}{2} \oint_{\partial V}\oint_{\partial V} \textrm{d} S_{ {\boldsymbol{ x}}_w} \, \textrm{d} S_{\tilde {\boldsymbol x}_w}\left[ \delta \theta({\boldsymbol x}_w) \left.\frac{\partial s_w}{\partial \theta} \right\vert_{eq}\delta_w({\boldsymbol x}_w - {\tilde {\boldsymbol x}_w} ) \delta \theta({\tilde {\boldsymbol x}_w})\right] . \end{align}

In the above equation, $\delta ({\boldsymbol x} - {\tilde {\boldsymbol x}})$ is the ordinary three-dimensional Dirac delta function and $\delta _w({\boldsymbol x}_w - {\tilde {\boldsymbol x}_w} )$ denotes the Dirac delta function on the manifold,

(5.6)\begin{equation} \oint_{\partial V} f({\tilde {\boldsymbol x}_w}) \delta_w({\boldsymbol x}_w - {\tilde {\boldsymbol x}_w} ) \, \textrm{d} S_{\tilde {\boldsymbol x}_w} := \int_{{\mathcal{S}}} \, \textrm{d} S_{\tilde {\boldsymbol \xi}} \,\, f({\tilde {\boldsymbol \xi}}) \delta_{2D}({\boldsymbol \xi} - {\tilde {\boldsymbol \xi}} ), \end{equation}

where ${\boldsymbol x}_w = {\boldsymbol x}_w ({\boldsymbol \xi }) \in \mathbb {R}^3$, ${\boldsymbol \xi } = (\xi ^1,\xi ^2) \in {\mathcal {S}} \subset \mathbb {R}^2$ is the parametric equation (locally) describing the boundary and $\delta _{2D}$ is the ordinary two-dimensional Dirac delta function. Integration by parts is used twice to move the Laplacian $\nabla ^2$ from the density fluctuation $\delta \rho$ to the Dirac delta function. In a compact operator form, (5.5) reads

(5.7)\begin{equation} {\rm \Delta} S_c = {\rm \Delta} S_{c f} + {\rm \Delta} S_{c w} = -\tfrac{1}{2} \int_{V} {\boldsymbol \varDelta}^{\dagger}\boldsymbol{\mathcal{H}}_{\boldsymbol{f}} \,{\boldsymbol \varDelta} \,\textrm{d} V -\tfrac{1}{2} \oint_{\partial V} {\delta \theta}^{\dagger}{\mathcal{H}}_w \,{\delta\theta} \, \textrm{d} S , \end{equation}

with $\boldsymbol {\varDelta } = (\delta \rho ,\delta \boldsymbol {v},\delta \theta )$ denoting the fluctuating fields, and $\boldsymbol{\mathcal{H}}_{\boldsymbol{f}}$ and ${\mathcal {H}}_w$ are diagonal, positive definite operators, defined as

(5.8)\begin{gather} (\boldsymbol{\mathcal{H}}_{\boldsymbol{f}} {\boldsymbol \varDelta})({\boldsymbol x}) = \int_V {\boldsymbol{\mathsf{H}}}_f({\boldsymbol x}, {\tilde {\boldsymbol x}}) {\boldsymbol\varDelta}({ \tilde{\boldsymbol x}} )\,\textrm{d} V_{ \tilde{\boldsymbol x}} = \int_V \hat{\boldsymbol{\mathsf{H}}}_f({\boldsymbol x}) \delta({\boldsymbol x}-{\tilde {\boldsymbol x}}) {\boldsymbol\varDelta}({ \tilde{\boldsymbol x}} )\,\textrm{d} V_{ \tilde{\boldsymbol x}}, \end{gather}

(5.9)\begin{gather}({\mathcal{H}}_w {\delta \theta})({\boldsymbol x}_w) = \oint_{\partial V} {\boldsymbol{\mathsf{H}}}_w({\boldsymbol x}_w, {\tilde {\boldsymbol x}_w}) {\delta\theta}({\tilde {\boldsymbol x}_w}) \, \textrm{d} S_{\tilde {\boldsymbol x}_w} = \oint_{\partial V} \hat{\boldsymbol{\mathsf{H}}}_w({\boldsymbol x}_w) \delta_w({\boldsymbol x}_w-{\tilde {\boldsymbol x}_w}) {\delta \theta }({\tilde {\boldsymbol x}_w}) \, \textrm{d} S_{\tilde {\boldsymbol x}_w}, \end{gather}

where

(5.10)\begin{equation} \boldsymbol {\hat H}_f = \left( \begin{array}{@{}ccc@{}} \dfrac{c_{T,eq}^{2}}{\theta_{eq}\rho_{eq}}-\dfrac{\lambda}{\theta_{eq}}\nabla_{\boldsymbol x}^2 & 0 & 0 \\ 0 & \dfrac{\rho_{eq}}{\theta_{eq}} {\boldsymbol{\mathsf{I}}} & 0 \\ 0 & 0 & \dfrac{\rho_{eq}c_{v,eq}}{\theta_{eq}^{2}} \end{array} \right) \end{equation}

involves differential operators with ${\boldsymbol{\mathsf{I}}}$ the $3 \times 3$ identity matrix and

(5.11)\begin{equation} \hat{\boldsymbol{\mathsf{H}}}_w = \left.\frac{\partial s_w}{\partial \theta} \right\vert_{eq} . \end{equation}

The probability distribution functional for the fluctuating fields $\boldsymbol {\varDelta }$ (Einstein 1956),

(5.12)\begin{equation} P_{eq}[\boldsymbol{\varDelta}]=\frac{1}{Z}\exp\left({\frac{{\rm \Delta} S_c}{k_{B}}}\right) = \frac{1}{Z_f}\exp\left({\frac{{\rm \Delta} S_{c f}}{k_{B}}}\right) \frac{1}{Z_w}\exp\left({\frac{{\rm \Delta} S_{c w}}{k_{B}}}\right), \end{equation}

can be approximated by exploiting the second-order approximation, (5.7),

(5.13)\begin{equation} P_{eq}[{\boldsymbol \varDelta}]=\frac{1}{Z_f}\exp\left({-\frac{1}{2k_{B}} \int_{V}{\boldsymbol \varDelta}^{\dagger}\boldsymbol{\mathcal{H}}_{\boldsymbol{f}}\,{\boldsymbol \varDelta}\,\textrm{d} V} \right) \frac{1}{Z_w}\exp\left({-\frac{1}{2k_{B}}\oint_{\partial V} \delta\theta^{\dagger}{\mathcal{H}}_w\,\delta\theta\,\textrm{d} S} \right). \end{equation}

Owing to the exponential form of the probability distribution functional, the diagonal nature of the operator $\hat {\boldsymbol{\mathsf{H}}}_f$ together with the local dependence on temperature of the entropy, $P_{eq}$ can be factorized as

(5.14)\begin{equation} P_{eq}[{\boldsymbol \varDelta}] = P_f[\boldsymbol{\varDelta}] P_w[\delta\theta] = P_f[\delta\rho] P_f[\delta\boldsymbol{v}] P_f[\delta\theta] P_w[\delta\theta] , \end{equation}

which implies that the fluctuations of temperature in the fluid and in the layer are statistically independent. Focusing on the fluid domain, the relevant property of the probability distribution functional is the correlation function

(5.15)\begin{equation} {\boldsymbol{\mathsf{C}}}_f(\boldsymbol{x}, \tilde{\boldsymbol{x}} ) = \langle\boldsymbol{\varDelta}(\boldsymbol{x})\otimes \boldsymbol{\varDelta}^{\dagger}( \tilde{\boldsymbol{x}} )\rangle = \frac{1}{Z}\int {\mathcal{D}}\boldsymbol{\varDelta} \boldsymbol{\varDelta}(\boldsymbol{x})\otimes \boldsymbol{\varDelta}^{\dagger}( \tilde{\boldsymbol{x}} ) P_{eq} [\boldsymbol{\varDelta} ], \end{equation}

where the integral over all the possible field fluctuations is understood on the right-hand side (path integral). It can be evaluated in closed form given the quadratic approximation for the entropy, which leads to Gaussian path integrals. The final expression is (see Gallo et al. (2018) for details)

(5.16)\begin{equation} {\boldsymbol{\mathsf{C}}}_f (\boldsymbol{x}, \tilde{\boldsymbol{x}} ) = k_B {\boldsymbol{\mathsf{H}}}_f^{-1} (\boldsymbol{x}, \tilde{\boldsymbol{x}}), \end{equation}

implying the following identity:

(5.17)\begin{equation} \int_V {\boldsymbol{\mathsf{H}}}_f (\boldsymbol{x}, \hat{\boldsymbol x}) {\boldsymbol{\mathsf{C}}}_f (\hat{\boldsymbol x}, \tilde{\boldsymbol{x}} ) \, \textrm{d} V_{\hat{\boldsymbol x}} = {\boldsymbol{\mathsf{I}}} k_B \delta (\boldsymbol{x} - \tilde{\boldsymbol{x}}). \end{equation}

From a practical standpoint, the correlation functions are obtained by solving the equivalent equation

(5.18)\begin{equation} \hat{\boldsymbol{\mathsf{H}}}_f(\boldsymbol{x}) {\boldsymbol{\mathsf{C}}}_f(\boldsymbol{x}, \tilde{\boldsymbol{x}}) = {\boldsymbol{\mathsf{I}}} k_B \delta (\boldsymbol{x} - \tilde{\boldsymbol{x}}) \end{equation}

(recall that ${\boldsymbol{\mathsf{H}}}_f (\boldsymbol { x}, \hat {\boldsymbol x})$ and $\hat {\boldsymbol{\mathsf{H}}}_f(\boldsymbol { x})$ are slightly different objects, (5.8)). The velocity and temperature correlations at equilibrium are straightforwardly obtained as

(5.19)\begin{gather} {\boldsymbol{\mathsf{C}}}_f^{{\delta \boldsymbol{v} \delta \boldsymbol{v}}} (\boldsymbol{x}, \tilde{\boldsymbol{x}}) = \frac{k_B \theta_{eq}} {\rho_{eq}} \delta(\boldsymbol{x} - \tilde{\boldsymbol{x}}) {\boldsymbol{\mathsf{I}}}, \end{gather}

(5.20)\begin{gather}{\boldsymbol{\mathsf{C}}}_f^{\delta \theta \delta \theta}(\boldsymbol{x}, \tilde{\boldsymbol{x}}) = \frac{k_B \theta_{eq}^2} {c_{v,eq}} \delta (\boldsymbol{x} - \tilde{\boldsymbol{x}}), \end{gather}

respectively. The density correlation involves instead a differential operator, calling for the solution of the differential equation

(5.21)\begin{equation} \left(\frac{c_{T,eq}^{2}}{\theta_{eq}\rho_{eq}} -\frac{\lambda}{\theta_{eq}}\nabla_{\boldsymbol x}^2\right) { {\boldsymbol{\mathsf{C}}}_f^{{\delta \rho \delta \rho}}} (\boldsymbol{x}, \tilde{\boldsymbol{x}}) = k_B \delta (\boldsymbol{x} - \tilde{\boldsymbol{x}}), \end{equation}

complemented with the boundary conditions specifying the wall wettability. Since the linearized form of (3.8) reads

(5.22)\begin{equation} \left.\lambda\frac{\partial \delta\rho}{\partial n} + \frac{\partial^2 f_w}{\partial \rho^2}\right\vert_{eq} \delta\rho = 0, \end{equation}

the boundary conditions for the density–density correlation ${\boldsymbol{\mathsf{C}}}_f^{\delta \theta \delta \theta } = \langle \delta \rho (\boldsymbol {x}) \delta \rho (\tilde {\boldsymbol {x}}) \rangle$, i.e. the Green's function in (5.21), is

(5.23)\begin{equation} \left.\lambda\frac{\partial {{\boldsymbol{\mathsf{C}}}_f^{{\delta \rho \delta \rho}}} (\boldsymbol{x}, \tilde{\boldsymbol{x}})}{\partial n_{\boldsymbol x}} + \frac{\partial^2 f_w}{\partial \rho^2}\right\vert_{eq} { {\boldsymbol{\mathsf{C}}}_f^{{\delta \rho \delta \rho}}}(\boldsymbol{x}, \tilde{\boldsymbol{x}}) = 0. \end{equation}

At variance with the unbounded system (Gallo et al. 2018), a closed-form solution of (5.21) is not readily available in the presence of the solid wall, because of the inhomogeneity of the coefficients ($C_{T,eq} = C_{T,eq}(\boldsymbol {x})$) in the linear equation for the Green's function, due to the stratification of the equilibrium density at the wall, $\rho _{eq} = \rho _{eq}(\boldsymbol {x})$, which exists whenever the equilibrium contact angle differs from $90^\circ$ (neutral wettability).

6. Fluctuation–dissipation balance for wall-bounded capillary fluids: general theory

In a nutshell, fluctuating hydrodynamics amounts to including suitable stochastic forcing terms into the linearized version of the field equations. The noise is assumed to be generated by a linear operator (noise operator) acting on a zero-mean, delta-correlated-in-time Gaussian process (white noise). The equilibrium solution of the resulting system of SPDEs is a Gaussian field, completely determined by the correlations. After the noise operator is determined in such a way to generate the equilibrium correlations independently defined in terms of the entropy functional, the resulting SPDEs will sample the Gaussian equilibrium probability functional $P_f[\boldsymbol {\varDelta }]$ discussed in the previous section. In other words, once the equilibrium correlation functions of the fluctuating fields are known, the fluctuation–dissipation balance (FDB) provides the explicit form of the stochastic forcing. Eventually, the capillary Navier–Stokes equations (4.2) and (4.8) completed with the additional stochastic contributions (4.9) result in a model that could be dubbed the capillary Landau–Lifshitz–Navier–Stokes (CLLNS) system. The purpose of the present section is to extend this model to the case of a wall-bounded capillary fluid. As will be shown in detail, the FDB can be satisfied keeping the noise in the form of a standard white noise. The noise terms act on momentum and energy and, also in the presence of walls, keeps the divergence form, consistently with the expected conservation properties. Special care is devoted to include the appropriate boundary conditions in the relevant operators, a crucial point in order to avoid the need to modify the structure of the noise.

The three-dimensional CLLNS system augmented with noise is

(6.1)\begin{equation} \left.\begin{array}{c@{}} \dfrac{\partial\rho}{\partial t}+ \boldsymbol{\nabla}\boldsymbol{\cdot}( \rho\boldsymbol{v}) = 0 , \\ \dfrac{\partial\rho\boldsymbol{v}}{\partial t}+ \boldsymbol{\nabla}\boldsymbol{\cdot}( \rho\boldsymbol{v}\otimes\boldsymbol{v}) = \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\varSigma} + \boldsymbol{f}_\varSigma ,\\ \dfrac{\partial E}{\partial t}+ \boldsymbol{\nabla}\boldsymbol{\cdot}( \boldsymbol{v} E) = \boldsymbol{\nabla}\boldsymbol{\cdot}( \boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\varSigma} -\boldsymbol{q} ) + \boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{f}_\varSigma + f_q , \end{array}\right\} \end{equation}

where $E$ is the total energy density, $E = {\mathcal {U}} + \frac 12 \rho \vert \boldsymbol {v}\vert ^2 + \frac 12 \lambda \vert \boldsymbol {\nabla }\rho \vert ^2$, with ${\mathcal {U}}$ the internal energy density. In the momentum and energy equations, respectively, $\boldsymbol \varSigma$ and $\boldsymbol q$ are the classical deterministic stress tensor and energy flux (derived in § 4), respectively. The boundary conditions consist of the no-slip/impermeability condition for velocity, zero heat flux and prescribed wall wettability at the solid surface:

(6.2a–c)\begin{equation} \rho \boldsymbol{v} = 0, \quad \boldsymbol{q}_h \boldsymbol{\cdot} \hat{\boldsymbol{n}} = 0, \quad \lambda\boldsymbol{\nabla}\rho\boldsymbol{\cdot}\hat{\boldsymbol{n}} + \frac{\partial f_w}{\partial \rho} = 0. \end{equation}

The stochastic forces, whose statistical properties are to be inferred from the fluctuation–dissipation theorem, are $\boldsymbol {f}_\varSigma$ and $f_q$. After introducing the equilibrium state $\boldsymbol {U}_{eq} =(\rho _{eq}(\boldsymbol {x}), \boldsymbol {0}, \theta _{eq})$ and denoting the fluctuating fields by $\boldsymbol {\varDelta } = (\delta \rho ,\delta \boldsymbol {v},\delta \theta )$, as previously, the state is described by $\boldsymbol {U} = \boldsymbol {\varDelta } + \boldsymbol {U}_{eq}$. With this notation, equations (6.1) written in a compact form are

(6.3)\begin{equation} \left.\begin{array}{c@{}} \dfrac{\partial \boldsymbol{U}}{\partial t} = {\boldsymbol N}[ \boldsymbol{U} ] + \boldsymbol{f}, \quad \boldsymbol{x}\in {\mathcal{D}}, \\ \boldsymbol B [\boldsymbol U] = 0 \quad \boldsymbol{x} \in \partial {\mathcal{D}} , \end{array}\right\} \end{equation}

with $\boldsymbol B$ the operator that specifies the boundary conditions. Linearization around the equilibrium state $\boldsymbol {U}_{eq}$ leads to

(6.4)\begin{align} \left.\begin{array}{c@{}} \dfrac{\partial\delta\rho}{\partial t} = - \boldsymbol{\nabla}\boldsymbol{\cdot}( \rho_{eq} \delta\boldsymbol{v}) , \\ \dfrac{\partial\rho_{eq}\delta\boldsymbol{v}}{\partial t} = - \rho_{eq} \boldsymbol{\nabla}\left( \dfrac{c_{eq}^2}{\rho_{eq}} - \lambda\nabla^2 \right)\delta\rho + \eta\nabla^2\delta{\boldsymbol v} + \dfrac{1}{3}\eta\boldsymbol{\nabla}\boldsymbol{\nabla}\boldsymbol{\cdot} \delta\boldsymbol v - \boldsymbol{\nabla}(\partial_\theta p\vert_{eq} \delta\theta) + \boldsymbol{f}_\varSigma ,\\ \dfrac{\partial c_{v,eq}\rho_{eq} \delta\theta}{\partial t} = - \theta_{eq} \partial_{\theta}p\vert_{eq}\boldsymbol \nabla \boldsymbol{\cdot} \delta\boldsymbol{v} + k \nabla^2\delta\theta + f_q . \end{array}\right\} \end{align}

The linearized system can be rewritten in conservative (divergence) form with the change of variables ${\boldsymbol \varDelta }' = (\delta \rho , \rho _{eq} \delta \boldsymbol {v}, c_{v,eq}\rho _{eq} \delta \theta ) = (\delta \rho , \delta {\boldsymbol \pi }, \delta \iota )$,

(6.5)\begin{equation} \left.\begin{array}{c@{}} \dfrac{\partial \boldsymbol{\varDelta}'}{\partial t} = \boldsymbol{\mathcal{L}} \boldsymbol{\varDelta}' + \boldsymbol{f}, \quad \boldsymbol{x} \in {\mathcal{D}}, \\ \boldsymbol{\mathcal{B}} \boldsymbol{\varDelta}' = 0, \quad \boldsymbol{x} \in \partial {\mathcal{D}} , \end{array}\right\} \end{equation}

where the (linearized) hydrodynamic operator $\boldsymbol {\mathcal {L}}$ is

(6.6)\begin{equation} \boldsymbol{\mathcal{L}}=\left(\begin{array}{@{}ccc@{}} 0 & - \boldsymbol \nabla \boldsymbol{\cdot} & 0\\ -\rho_{eq}\boldsymbol \nabla \left( \dfrac{c_{T,eq}^2}{\rho_{eq}} - \lambda \nabla^2 \right) & {\eta}\left( \nabla^2 + \dfrac{1}{3} \boldsymbol \nabla \boldsymbol \nabla \boldsymbol{\cdot} \right)\dfrac{1}{\rho_{eq}} & -\boldsymbol \nabla \left(\partial_{\theta}p\vert_{eq} \dfrac{1}{\rho_{eq} c_{v,eq}} \right) \\ 0 & -\theta_{eq} \partial_{\theta}p\vert_{eq} \boldsymbol \nabla \boldsymbol{\cdot} \dfrac{1}{\rho_{eq}} & k \nabla^2 \dfrac{1}{\rho_{eq}c_{v,eq}} \end{array}\right) , \end{equation}

and $\boldsymbol {\mathcal {B}}$ is the (linearized) boundary conditions operator

(6.7a,b)\begin{equation} \left.\lambda\frac{\partial \delta\rho}{\partial n} + \frac{\partial^2 f_w}{\partial \rho^2}\right\vert_{eq} \delta \rho = 0,\quad \delta{\boldsymbol \pi} = \boldsymbol {0} , \quad \lambda \frac{\partial \rho_{eq}}{\partial n} \boldsymbol{\nabla}\boldsymbol{\cdot}\delta {\boldsymbol \pi} + k \frac {\partial}{\partial n} \left( \frac{\delta {\iota}}{c_{v eq}\rho_{eq}} \right)= 0 . \end{equation}

The stochastic contribution can be written through the linear operator $\boldsymbol {\mathcal {K}}$ acting on independent white noise processes ${\boldsymbol \varXi }$, as $\boldsymbol f = \boldsymbol {\mathcal {K}} {\boldsymbol \varXi }$, where

(6.8)\begin{equation} \boldsymbol{\mathcal{K}} = \left(\begin{array}{ccc} 0 & 0 & 0\\ 0 & \boldsymbol{\sigma}_{\boldsymbol\pi} & 0\\ 0 & 0 & \sigma_{\iota} \end{array}\right) . \end{equation}

In three dimensions ${\sigma _{\pi }}$ is a $3 \times 3$ matrix whose components are scalar linear operators to be determined. The vector ${\boldsymbol \varXi } = (\varXi^\rho , {\boldsymbol \varXi }^\pi , \varXi^\iota )^\textrm {T}$ collects the white noise processes for the fluctuating fields, with ${\boldsymbol \varXi }^\pi = ( \varXi^{\pi _x}, \varXi^{\pi _y}, \varXi^{\pi _z} )^\textrm {T}$. The Gaussian process ${\boldsymbol \varXi }$ has zero mean and correlations given by

(6.9)\begin{equation} \langle{\boldsymbol \varXi}( \tilde{\boldsymbol{y}},s) \otimes {\boldsymbol \varXi}^\textrm{T}(\hat{\boldsymbol y},q)\rangle = {\boldsymbol{\mathsf{I}}} \delta(\tilde{\boldsymbol{y}} - \hat{\boldsymbol y}) \delta(s-q) , \end{equation}

with ${{\boldsymbol{\mathsf{I}}}}$ now indicating a $( 5 \times 5 )$ identity matrix in ${\boldsymbol \varXi }$ space, $\tilde {\boldsymbol {y}}$ and $\hat {\boldsymbol y}$ denoting space points, and $s$ and $q$ time instants. It may be worth recalling that spatial and temporal Dirac delta functions carry dimensions of reciprocal space and time, respectively. Hence the noise ${\boldsymbol \varXi }$ has dimensions of reciprocal square root of time per volume.

Introducing the operator ${\boldsymbol{\mathsf{L}}}$, which accounts for the differential operator (6.6) completed with the boundary conditions (6.7a,b), problem (6.5) is compactly rewritten as $\textrm {d} \boldsymbol {\varDelta }'/\textrm {d} t = {\boldsymbol{\mathsf{L}}} \boldsymbol {\varDelta }' + \boldsymbol {f}$ (in other words, the action of the operator $\boldsymbol {\mathcal {L}}'$ is restricted to the subspace of fluctuations that satisfy the boundary conditions). The solution can then be expressed in terms of the propagator

(6.10)\begin{equation} {\mathcal{P}}_t = \exp( {\boldsymbol{\mathsf{L}}} t). \end{equation}

Considering that ${\mathcal {P}}_{t} {\mathcal {P}}_{-t} = {\mathcal {U}}$, with ${\mathcal {U}}$ the relevant identity operator, and using the propagator as functional integrating factor, system (6.5) becomes

(6.11)\begin{equation} \frac{\partial}{\partial t} ({\mathcal{P}}_{-t} \boldsymbol{\varDelta}' ) = {\mathcal{P}}_{-t} \,\, \boldsymbol{f}(t) , \end{equation}

where $\textrm {d} ({\mathcal {P}}_{-t}) / \textrm {d} t = - {\boldsymbol{\mathsf{L}}} {\mathcal {P}}_{-t}$ and the obvious commutation ${\boldsymbol{\mathsf{L}}} {\mathcal {P}}_{-t} = {\mathcal {P}}_{-t} {\boldsymbol{\mathsf{L}}}$ has been used. This equation can be now integrated in time, leading to the so-called mild solution (Da Prato 2012)

(6.12)\begin{equation} \boldsymbol{\varDelta}'(\boldsymbol{x},t) =\int_{0} ^{t}{\mathcal{P}}_{t-s}\,\,\boldsymbol{f}(s)\,\textrm{d} s + {\mathcal{P}}_{t} \boldsymbol{\varDelta}'_{0} , \end{equation}

where the last term retains memory of the initial conditions and vanishes for large times since ${\boldsymbol{\mathsf{L}}}$ is a dissipative operator. Consistently, the equilibrium correlations take the explicit form

(6.13)\begin{equation} \langle \boldsymbol{\varDelta}'(\tilde{\boldsymbol{x}},t)\otimes \boldsymbol{\varDelta}^{\prime\dagger}(\hat{\boldsymbol x},t)\rangle =\int_{0}^{t}{\mathcal{P}}_{t-s} \, {\boldsymbol{\mathsf{Q}}}\,{\mathcal{P}}_{t-s}^\dagger \, \textrm{d} s , \end{equation}

where the stochastic force correlation

(6.14)\begin{equation} \boldsymbol{\mathcal{Q}} = \langle \, \boldsymbol{f}(\boldsymbol{x}, s)\otimes \,\boldsymbol{f}^\dagger( \tilde{\boldsymbol{x}}, q)\rangle = {\boldsymbol{\mathsf{Q}}}(\boldsymbol{x}, \tilde{\boldsymbol{x}}) \delta(s-q) \end{equation}

is straightforwardly expressed in terms of the noise operator $\boldsymbol{\mathcal{K}}$ and the stochastic process ${\boldsymbol \varXi }$. Note that it should be stressed that the delta correlation in time for the stochastic forces is always adopted in fluctuating hydrodynamics based on the assumption of a clear scale separation between molecular and hydrodynamic time scales (Markovian dynamics) – see e.g. Duque-Zumajo et al. (2019) for additional comments.

In view of producing the explicit correlation (6.13), the Hermitian non-singular operator ${\boldsymbol{\mathsf{E}}}^{-1}$ is assumed to exist and is able to factorize ${\boldsymbol{\mathsf{Q}}}$ as

(6.15)\begin{equation} {\boldsymbol{\mathsf{Q}}}=-( {\boldsymbol{\mathsf{L}}} {\boldsymbol{\mathsf{E}}}^{-1} + {\boldsymbol{\mathsf{E}}}^{-1} {\boldsymbol{\mathsf{L}}}^\dagger ). \end{equation}

In this case the integrand in (6.13) would be the $s$ derivative of ${\mathcal {P}}_{t-s}\,{\boldsymbol{\mathsf{E}}}^{-1} \,{\mathcal {P}}_{t-s}^\dagger$, implying that the exact integral

(6.16)\begin{equation} \lim_{t\rightarrow\infty}\langle \boldsymbol{\varDelta}'\otimes \boldsymbol{\varDelta}^{\prime\dagger} \rangle = {\boldsymbol{\mathsf{E}}}^{-1} = {\boldsymbol{\mathsf{C}}}'_f . \end{equation}

Hence the operator ${\boldsymbol{\mathsf{E}}}^{-1}$ does indeed exist and is the correlation matrix ${\boldsymbol{\mathsf{C}}}_f$ (5.16), resulting in

(6.17)\begin{equation} {\boldsymbol{\mathsf{C}}}'_f = \left( \begin{array}{ccc} {\boldsymbol{\mathsf{C}}}_f^{{\delta \rho \delta \rho}}(\boldsymbol{x}, \tilde{\boldsymbol{x}}) & 0 & 0 \\ 0 & {k}_B \theta_{eq} \rho_{eq} \delta (\boldsymbol{x} - \tilde{\boldsymbol{x}}) {\boldsymbol{\mathsf{I}}} & 0 \\ 0 & 0 & {k}_B \theta_{eq}^{2} \rho_{eq} c_{v,eq} \delta (\boldsymbol{x} - \tilde{\boldsymbol{x}}) \end{array} \right), \end{equation}

By combining (6.15) and (6.16) the stochastic force correlation follows:

(6.18)\begin{equation} {\boldsymbol{\mathsf{Q}}} =-( {\boldsymbol{\mathsf{L}}}{\boldsymbol{\mathsf{C}}}'_f + {\boldsymbol{\mathsf{C}}}'_f{\boldsymbol{\mathsf{L}}}^{\dagger}) = ( \boldsymbol{\mathcal{M}} + \boldsymbol{\mathcal{M}}^{\dagger}) = 2 {k}_B \boldsymbol{O}, \end{equation}

where $\boldsymbol {\mathcal {M}} = - {\boldsymbol{\mathsf{L}}} {\boldsymbol{\mathsf{C}}}'_f$ and $\boldsymbol {O}$ is the Onsager matrix. Relationship (6.18) is the form the celebrated FDB takes for the present system highlighting the connection between fluctuation intensity and dissipation mechanisms. The physical interpretation is that, in thermodynamic equilibrium, the response of a system to a perturbation is equivalent to the evolution of a spontaneous fluctuation. One is then enabled to infer non-equilibrium properties of a physical system from equilibrium properties.

7. Fluctuation–dissipation balance: explicit expression for neutrally wettable surfaces

The present section is devoted to obtain the explicit form of the stochastic forces following the general procedure outlined in the previous section. For the sake of simplicity, the contact angle $\phi = {\rm \pi}/2$ is assumed, in order to deal with homogeneous equilibrium fields. In this case the boundary conditions operator $\boldsymbol {\mathcal {B}}$, besides enforcing vanishing velocity ($\delta \boldsymbol{\pi} = \boldsymbol 0$) and zero heat flux at the walls ($\partial \delta \iota /\partial n = 0$), also requires vanishing density normal derivative ($\partial \delta \rho /\partial n = 0$) (see (6.7a,b)). These boundary conditions describe the evolution of a closed and isolated system, in the spirit of the $NVE$ ensemble of statistical mechanics.

Compared with the general case, the above assumptions lead to homogeneous equilibrium fields (no stratification at the wall) with the advantage of facilitating the mathematical treatment. The extension to the case of generically wettable surfaces – albeit of greater complexity from the mathematical point of view – is not expected to alter the structure of stochastic forces. In fact, the wall wettability does not introduce an additional dissipation mechanism, and, as such, does not produce entropy. Hence a general solid–fluid energy is expected not to alter the fluctuation–dissipation balance.

Following the general procedure illustrated in § 6, the operators $\sigma _{\boldsymbol \pi }$ and $\sigma _{\iota }$ can be identified from the FBD (6.18),

(7.1)\begin{equation} \boldsymbol{\mathcal{Q}} = \boldsymbol{\mathcal{K}} \langle {\boldsymbol \varXi} {\boldsymbol \varXi}^{\dagger} \rangle \boldsymbol{\mathcal{K}}^{\dagger} = 2 {k}_B \boldsymbol{O} \delta(s-q) , \end{equation}

leading to

(7.2)\begin{equation} \boldsymbol{\mathcal{K}} \boldsymbol{\mathcal{K}}^{\dagger} = 2 {k}_B \boldsymbol{O} = \boldsymbol{\mathcal{M}} + \boldsymbol{\mathcal{M}}^{\dagger} . \end{equation}

It should be stressed that the operator entering the definition of $\boldsymbol {\mathcal {M}}$ is ${\boldsymbol{\mathsf{L}}}$. In practice, to solve the problem, one may start from the differential operator $\boldsymbol {\mathcal {L}}$, (6.6), and the correlation matrix ${\boldsymbol{\mathsf{C}}}_f$, (5.16). Since the correlation matrix satisfies the boundary conditions, this leads to the same result that the actual operator ${\boldsymbol{\mathsf{L}}}$, i.e. $\boldsymbol {\mathcal {L}}$ completed with boundary conditions, would imply. One finds

(7.3)\begin{equation} \boldsymbol{\mathcal{M}} = \left(\begin{array}{ccc} 0 & m_{12} & 0 \\ m_{21} & \boldsymbol m_{22} & m_{23} \\ 0 & m_{32} & m_{33} \end{array} \right) . \end{equation}

The entries of the matrix $\boldsymbol {\mathcal {M}}$ are

(7.4)\begin{gather} m_{12} = \boldsymbol \nabla \boldsymbol{\cdot} [ \rho_{eq} {k}_B \theta_{eq} \delta ( \boldsymbol x - \hat{\boldsymbol x} ) {\boldsymbol{\mathsf{I}}} ], \end{gather}

(7.5)\begin{gather}m_{21} = \rho_{eq} \boldsymbol \nabla [ {k}_B \theta_{eq} \delta ( \boldsymbol x - \hat{\boldsymbol x} ) ], \end{gather}

(7.6)\begin{gather}m_{23} = \boldsymbol \nabla [ {k}_B \theta_{eq}^ 2 \partial_\theta p\vert_{eq} \delta ( \boldsymbol x - \hat{\boldsymbol x} ) ], \end{gather}

(7.7)\begin{gather}m_{32} = \partial_\theta p\vert_{eq} \boldsymbol \nabla \boldsymbol{\cdot} [ {k}_B \theta_{eq}^ 2 \delta ( \boldsymbol x - \hat{\boldsymbol x} ) {\boldsymbol{\mathsf{I}}} ], \end{gather}

(7.8)\begin{gather}\boldsymbol m_{22} = - \eta ( {\boldsymbol{\mathsf{I}}} \nabla^2 + \tfrac{1}{3} \boldsymbol \nabla \otimes \boldsymbol \nabla ) {k}_B \theta_{eq} \delta ( \boldsymbol x - \hat{\boldsymbol x} ), \end{gather}

(7.9)\begin{gather}m_{33} = - k_0 \nabla^2 {k}_B \theta_{eq}^2 \delta ( \boldsymbol x - \hat{\boldsymbol x} ). \end{gather}

Thus, the sum of $\boldsymbol {\mathcal {M}}$ and its hermitian conjugate $\boldsymbol {\mathcal {M}}^{\dagger}$ provides the explicit expression for the square of the unknown matrix operator $\boldsymbol {\mathcal {K}}$, (6.8), i.e. the explicit form of the FDB,

(7.10)\begin{equation} \boldsymbol{\mathcal{K}} \boldsymbol{\mathcal{K}}^{\dagger} = \boldsymbol{\mathcal{M}} + \boldsymbol{\mathcal{M}}^{\dagger} = \left(\begin{array}{ccc} 0 & 0 & 0 \\ 0 & 2 \boldsymbol m_{22} & 0\\ 0 & 0 & 2 m_{33} \end{array} \right). \end{equation}

In deriving (7.10), one needs to take into account that the differential operators $\boldsymbol {\nabla }$ and $\nabla ^2$ are effectively (and not purely formally) anti- and self-adjoint, respectively, thanks to the boundary conditions ($\boldsymbol {\nabla }^\dagger = - \boldsymbol {\nabla }\boldsymbol {\cdot }$ and $(\nabla ^2)^\dagger = \nabla ^2$; appendix B). After recalling (6.8), determining $\boldsymbol {\mathcal {K}}$ amounts to solving (7.10) component-wise,

(7.11)\begin{gather} \sigma_{\iota} \sigma^{\dagger} _{\iota} = -2 {k}_B \theta_{eq}^2 k \nabla^2 \delta ( \hat{\boldsymbol x} - \tilde{\boldsymbol x} ), \end{gather}

(7.12)\begin{gather}\boldsymbol{\sigma_\pi} \otimes \boldsymbol{\sigma}_{\pi}^{\dagger} = - 2 \eta {k}_B \theta_{eq} ( {\boldsymbol{\mathsf{I}}}\, \nabla^2 + \tfrac{1}{3} \boldsymbol \nabla \otimes \boldsymbol \nabla ) \delta( \hat{\boldsymbol x} - \tilde{\boldsymbol x} ). \end{gather}

The explicit solution for the stochastic forces is

(7.13a,b)\begin{equation} \boldsymbol{f}_\varSigma = \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\delta \varSigma}, \quad \boldsymbol{f_q} = \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol {\delta q} , \end{equation}

where $\boldsymbol {\delta \varSigma }$ and $\boldsymbol {\delta q}$ are

(7.14)\begin{gather} \boldsymbol{\delta \varSigma} = \sqrt{2\eta {k}_B\theta_{eq}} \,\tilde{{\boldsymbol \varXi}}^{\boldsymbol\pi} - \tfrac{1}{3}\sqrt{2\eta {k}_B\theta_{eq}} \,\textrm{Tr} (\tilde{{\boldsymbol \varXi}}^{\boldsymbol\pi}){\boldsymbol{\mathsf{I}}}, \end{gather}

(7.15)\begin{gather}\boldsymbol{\delta q} = \sqrt{2k{k}_B\theta_{eq}^2}\,\boldsymbol{{\boldsymbol \varXi}}^{\iota}, \end{gather}

with ${\tilde {{\boldsymbol \varXi }}^{\boldsymbol\pi} } =( {{\boldsymbol \varXi }}^{\boldsymbol\pi} + {{\boldsymbol \varXi }}^{{\boldsymbol\pi} \,\textrm {T}} )/\sqrt {2}$ a stochastic symmetric tensor field, and ${{\boldsymbol \varXi }}^{\iota }$ a stochastic vector, with correlations

(7.16)\begin{gather} \langle{\varXi^\pi}_{\alpha\beta}(\hat{x},\hat{t})\, {\varXi^\pi}_{\gamma\delta}(\tilde{x}, \tilde{t})\rangle = \delta_{\alpha\gamma}\delta_{\beta\delta}\delta(\hat{x}-\tilde{x}) \delta(\hat{t}-\tilde{t}), \end{gather}

(7.17)\begin{gather}\langle{\varXi^{\iota}}_{\alpha}(\hat{x},\hat{t})\, {\varXi^{\iota}}_{\beta}(\tilde{x},\tilde{t})\rangle = \delta_{\alpha\beta}\delta(\hat{x}-\tilde{x}) \delta(\hat{t}-\tilde{t}). \end{gather}

It is immediately shown that expressions (7.14) and (7.15) satisfy (7.11) and (7.12), namely

(7.18)\begin{equation} \langle \sigma_{\iota} \boldsymbol{\varXi}^{\iota} \boldsymbol{\varXi}^{\iota {\dagger}}\sigma^{\dagger} _{\iota} \rangle = \langle {\boldsymbol{\nabla}}_{\hat{\boldsymbol{x}}} \boldsymbol{\cdot} \boldsymbol{\delta q} ( \hat {\boldsymbol x}, t ) \, {\boldsymbol{\nabla}}_{\tilde{\boldsymbol{x}}} \boldsymbol{\cdot} \boldsymbol{\delta q} ( \tilde {\boldsymbol x}, t ) \rangle = -2 {k}_B \theta_{eq}^2 k \nabla^2 \delta ( \hat{\boldsymbol x} - \tilde{\boldsymbol x} ) \end{equation}

and

(7.19)\begin{align} \langle \boldsymbol{\sigma_{\boldsymbol\pi}} {\boldsymbol \varXi}^{\boldsymbol\pi} \otimes {\boldsymbol \varXi}^{{\boldsymbol\pi} {\dagger}} \boldsymbol{\sigma}_{\boldsymbol{\pi}}^{\dagger} \rangle & = \langle {\boldsymbol{\nabla}}_{\hat{\boldsymbol{x}}} \boldsymbol{\cdot} \boldsymbol{\delta \varSigma}( \hat {\boldsymbol x}, t ) \otimes {\boldsymbol{\nabla}}_{\tilde{\boldsymbol{x}}} \boldsymbol{\cdot} \boldsymbol{\delta \varSigma}( \tilde {\boldsymbol x}, t ) \rangle \nonumber\\ &= - 2 \eta {k}_B \theta_{eq} ( {\boldsymbol{\mathsf{I}}} \, \nabla^2 + \tfrac{1}{3} \boldsymbol \nabla \otimes \boldsymbol \nabla ) \delta( \hat{\boldsymbol x} - \tilde{\boldsymbol x} ). \end{align}

The covariance of the stochastic process corresponding to the fluctuating stress is

(7.20)\begin{equation} \langle \boldsymbol{\delta \varSigma}(\hat{x},\hat{t})\otimes \boldsymbol{\delta \varSigma}^\dagger(\tilde{x},\tilde{t})\rangle= {\boldsymbol{\mathsf{Q}}}^\varSigma\delta(\hat{x}-\tilde{x}) \delta(\hat{t}-\tilde{t}), \end{equation}

with

(7.21)\begin{equation} {\boldsymbol{\mathsf{Q}}}^\varSigma_{\alpha\beta\nu\eta}=2{k}_B\theta\eta ( \delta_{\alpha\nu}\delta_{\beta\eta}+\delta_{\alpha\eta}\delta_{\beta\nu}- \tfrac{2}{3}\delta_{\alpha\beta}\delta_{\nu\eta}). \end{equation}

Analogously, the covariance of the fluctuating heat flux is

(7.22)\begin{equation} \langle \boldsymbol{\delta q}(\hat{x},\hat{t}) \otimes\boldsymbol{\delta q}^\dagger(\tilde{x},\tilde{t})\rangle = {\boldsymbol{\mathsf{Q}}}^q\delta(\hat{x}-\tilde{x}) \delta(\hat{t}-\tilde{t}), \end{equation}

with

(7.23)\begin{equation} {\boldsymbol{\mathsf{Q}}}^q_{\alpha\beta}=2{k}_B\theta_{eq}^2 k\delta_{\alpha\beta}. \end{equation}

Several comments are in order: (i) The stochastic forcing terms are in divergence form for momentum and energy (temperature) equations, and absent for mass, consistently with integral conservation principles. (ii) The correlation between the thermodynamic force of different tensor ranks is zero, consistently with the Curie–Prigogine principle, i.e. $(\langle \boldsymbol{\delta} \boldsymbol {q}^\dagger (\tilde {x},\tilde {t})\otimes \boldsymbol{\delta} \boldsymbol {\varSigma }(\hat {x},\hat {t})\rangle =0)$. (iii) The formal expressions of the FDB, (6.18), and of the stochastic forces, (7.13a,b), are exactly the same as obtained in free space, Gallo et al. (2018), i.e. they are independent of the specific boundary conditions, as long as no additional entropy production mechanisms is implied. (iv) The presence of the wall is intrinsic to the functional representation of the operators, hence the statistical properties of the fluctuating fields will be affected by the boundary conditions, as already pointed out when commenting on the density correlation equation, (5.21).

As a final remark before closing the section, the reader should be aware that the FDB, as derived here, is strictly valid in the limit of small fluctuations about the equilibrium state, as follows from the original Landau–Lifshitz approach, which starts from the linearized Navier–Stokes equations. As common, the extension to the nonlinear case is based on the local equilibrium assumption (Bogoliubov 1946; De Zarate & Sengers 2006). The assumption consists in assuming that out-of-equilibrium fluctuation statistics can be obtained from an equilibrium state corresponding to the local hydrodynamic field. A more complete derivation of the nonlinear fluctuating hydrodynamics model for a simple fluid has been provided by Español (1998), obtaining the same noise terms as discussed here.

8. The capillary Landau–Lifshitz–Navier–Stokes model

One of the central results of the previous section is that the stochastic terms keep the same form they have for a homogeneous capillary system (the interested reader is referred to Gallo et al. (2017) for further details). The solid wall imposing no-slip conditions and zero heat transfer and the solid–liquid free energy contribution $f_w$ do not affect the FDB (although they do affect the fluctuations). This result is deeply rooted in the van der Waals thermodynamics of non-homogeneous fluid systems, where capillary effects (both the square-gradient volume term and the wall contribution, (3.1)) are included as free energy contributions. As a result, the capillary stresses are generated by the reversible part of the dynamics, without any entropy production. On the other hand, the fluctuation–dissipation balance asserts that only the dissipative part of the dynamics (irreversible) contributes to stochastic fluxes, and thus there are no capillary contributions to the stochastic part of the equations. As a consequence, the form of the stochastic fluxes in both the bounded and unbounded CLLNS model is exactly the same of the simpler classical fluctuating Landau–Lifshitz–Navier–Stokes (LLNS) model (Fox & Uhlenbeck 1970).

The full CLLNS model reads

(8.1)\begin{equation} \left.\begin{array}{c@{}} \dfrac{\partial\rho}{\partial t}+ \boldsymbol{\nabla}\boldsymbol{\cdot} (\rho\boldsymbol{v}) = 0 , \\ \dfrac{\partial\rho\boldsymbol{v}}{\partial t}+ \boldsymbol{\nabla}\boldsymbol{\cdot} (\rho\boldsymbol{v}\otimes\boldsymbol{v}) = \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\varSigma} + \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\delta \varSigma} ,\\ \dfrac{\partial E}{\partial t}+ \boldsymbol{\nabla}\boldsymbol{\cdot} (\boldsymbol{v} E) = \boldsymbol{\nabla}\boldsymbol{\cdot} [ (\boldsymbol{\varSigma} + \boldsymbol{ \delta \varSigma})\boldsymbol{\cdot} \boldsymbol{v} -\boldsymbol{q} + \boldsymbol{\delta q} ], \end{array}\right\} \end{equation}

where the deterministic stress tensor was defined in (4.8) and the deterministic energy flux in (4.9).

The presence of the wall is reflected in the boundary conditions, here the no-slip condition for the fluid velocity at the solid surface and vanishing heat flux, i.e. $\rho \boldsymbol {v} = 0$ and $\boldsymbol {q}_h\boldsymbol {\cdot }\hat {\boldsymbol n} = 0$. The wall wettability is controlled by the condition $\partial \rho / \partial n = g(\theta ,\phi )$, where $g=\cos \phi \sqrt {(2/\lambda) (w_b({\rho },\theta ) - w_b(\rho _V))}$ is positive for hydrophilic walls and negative for hydrophobic ones, (3.16).

8.1. Equation of state

Two additional relations are still needed to close the system of equations (8.1). They are the two equations of state relating thermodynamic pressure and internal energy to density and temperature, $p = p(\rho ,\theta )$ and $u_b = u_b(\rho ,\theta )$. Both follow from the (bulk) free energy density $f_b = f_b(\rho ,\theta )$ specified in (3.1). Here the free energy of a Lennard-Jones fluid (Johnson et al. 1993) will be assumed in order to be able to compare the present results with classical techniques, such as molecular dynamics simulations. Once the appropriate thermodynamic potential is selected, the other thermodynamic properties follow straightforwardly, e.g. the pressure $p = -\rho ^2 \partial (\,f_b/\rho )/\partial \rho$.

For a Lennard-Jones fluid, the relevant expressions are too cumbersome to be repeated here, but it may still be interesting to recall the general aspects of a two-phase, vapour–liquid, system. The phase diagram is reported in figure 5 with binodal and spinodal lines shown in the inset. The binodal, as locus of saturation densities at changing temperature, is obtained from the classical equilibrium conditions stating that identical pressures $p$ and chemical potential $\mu _c$ characterize bulk vapour and liquid. The saturation densities $\rho _{V\, sat}$ and $\rho _{L\, sat}$ are evaluated by (numerically) solving, at given temperature, the nonlinear system of equations

(8.2)\begin{equation} \left.\begin{array}{c@{}} p(\rho_{V\,sat}, \theta) = p(\rho_{L\,sat}, \theta), \\ \mu_c(\rho_{V\,sat}, \theta) = \mu_c(\rho_{L\,sat}, \theta) . \end{array}\right\} \end{equation}

The two spinodals, as locus of the densities ($\rho _{V\, spin}$ and $\rho _{L\, spin}$, respectively) where $\partial p/\partial \rho \vert _\theta = 0$, represent the thermodynamic limit of metastability. All the thermodynamic states placed between the binodal and the spinodal lines are metastable. In these states, nucleation is a thermally activated process requiring an activation energy, as discussed in § 2. The closer the metastable state is to the spinodal, the smaller the activation barrier and the higher the nucleation probability.

Figure 5. Phase diagram for the Lennard-Jones equation of state (Johnson, Zollweg & Gubbins 1993). In the main panel, the isotherm $\theta =1.25$ and the iso-chemical potential $\mu =\mu _{sat}$ with the saturation value are reported with dashed and dash-dotted lines, respectively. The saturation densities are identified as the two points with equal temperature, chemical potential and pressure; the red circle represents the vapour saturation point and the orange circle the liquid one. The other two circles, dark blue and light blue, represent the spinodal points, vapour and liquid, respectively, identified on the isotherm where $\partial p/\partial \rho = 0$. In the inset, the loci of all the saturation and spinodal points at different temperatures are reported in the $\rho$–$\theta$ plane.

9. Numerical results

Heterogeneous vapour bubble nucleation is here simulated by numerically solving the CLLNS system of equations presented in § 8. The configuration consists of a metastable liquid described by the equation of state of a Lennard-Jones fluid (Johnson et al. 1993) which is enclosed in a box with two flat solid surfaces perpendicular to the $z$ direction where no slip, zero heat flux and given wettability are assigned, using periodicity in the $x$–$y$ plane. These boundary conditions globally enforce volume, mass and total energy conservation.

In figure 6 a sketch of the simulation set-up is reported for the reader's convenience. The equations are made dimensionless using the following reference quantities: $\sigma = 3.4 \times 10^{-10}$ m as length, $\epsilon = 1.65 \times 10^{-21}$ J as energy, $m = 6.63 \times 10^{-26}$ kg as mass, $V_r = (\epsilon /m)^{1/2}$ as velocity, $T_r = \sigma /V_r$ as time, $\theta _r = \epsilon /k_B$ as temperature, $\eta _r = \sqrt {m\epsilon }/\sigma ^2$ as shear viscosity, $c_{vr} = mk_B$ as specific heat at constant volume and $k_r = \eta _r c_{vr}$ as thermal conductivity. The dimensionless fields are defined as $\rho ^* = \rho /\rho _r$, $\theta ^* = \theta /\theta _r$ and $\boldsymbol {v}^* = \boldsymbol {v}/V_r$. The dimensionless fluxes, (4.8), (4.9), (7.14) and (7.15), read

(9.1)\begin{align} \boldsymbol{{\varSigma}}^* &= \left( \frac{C}{2}\vert\boldsymbol{\nabla}^*\rho^*\vert^2 + \rho^*\boldsymbol{\nabla}^*\boldsymbol{\cdot} (C\boldsymbol{\nabla}^*\rho^*)\right) {\boldsymbol{\mathsf{I}}} - C\boldsymbol{\nabla}^*\rho^*\otimes\boldsymbol{\nabla}^*\rho^*\nonumber\\ & \quad + \eta^*\left[(\boldsymbol{\nabla}^*\boldsymbol{u}^* + \boldsymbol{\nabla}\boldsymbol{u}^{\textrm{T}*}) - \frac{2}{3}\boldsymbol{\nabla}^*\boldsymbol{\cdot} \boldsymbol{u} ^*{\boldsymbol{\mathsf{I}}}\right], \end{align}

(9.2)\begin{gather} \boldsymbol{q}^* = C \rho^* \boldsymbol{\nabla}^*\rho^*\boldsymbol{ \nabla}^*\boldsymbol{\cdot} \boldsymbol{u}^* - k^* \boldsymbol{\nabla}^*\theta^* , \end{gather}

(9.3)\begin{gather}\boldsymbol{\delta \varSigma}^* = \sqrt{2\eta^*\theta^*} \, \tilde{\boldsymbol{\varXi}}^{\pi*} - \tfrac{1}{3}\sqrt{2\eta^*\theta^*} \,\textrm{Tr}(\tilde{\boldsymbol{\varXi}}^{\pi*}) {\boldsymbol{\mathsf{I}}}, \end{gather}

(9.4)\begin{gather}\boldsymbol{\delta q}^* = \sqrt{2k^*\theta^{*2}}\,\boldsymbol{\varXi}^{\iota*} , \end{gather}

where $C=\lambda \rho _r/\sigma ^2 V_r^2$ is a capillary number. All throughout, the capillary number is fixed at $C=5.244$, in order to reproduce the surface tension expected of a Lennard-Jones fluid (Gallo et al. 2018). Hereafter the symbol $^*$ indicating dimensionless quantities is omitted, for the ease of notation. The system volume $V=750 \times 750 \times 500$ has been discretized on a uniform grid, containing $50$ cells in the $z$ direction and $75 \times 75$ in the $x$–$y$ plane. The dimensionless viscosity and thermal conductivity of a Lennard-Jones fluid (Rowley & Painter 1997) have been used to account for the variability of transport coefficients due to the vapour–liquid density contrast.

Figure 6. Sketch of the simulation set-up.

The system of equations (8.1) has been discretized in the spirit of the method of lines with the numerical scheme proposed by Balboa et al. (2012) and recently adopted and validated by these authors in Gallo et al. (2018). When dealing with stochastic equations, particular attention must be paid to the discrete version of the fluctuation–dissipation balance. In particular, the numerical scheme must reproduce the field statistics as predicted by Einstein–Boltzmann theory. In order to cope with this particular requirement, Balboa et al. (2012) proposed to exploit a staggered grid to spatially discretize the equations. Special attention is here dedicated to the numerical treatment of the boundary conditions at the wall – see details illustrated in appendix C. The time evolution has been performed by means of a second-order Runge–Kutta scheme (Delong et al. 2013), with a constant time step ${\rm \Delta} t = 0.1$, which is small enough to accurately solve acoustic and viscous time scales.

Several metastable conditions have been investigated: this section reports on four different sets of simulations at initial temperature $\theta =1.25$ and five different initial bulk densities, $\rho _L = 0.47, 0.475, 0.48, 0.485, 0.4875$, corresponding to a decreasing degree of metastability, $(\mu _c(\rho _L)-\mu _c^{sat})/(\mu _c^{spin}-\mu _c^{sat})$, with $\mu _c^{sat}$ and $\mu _c^{spin}$ the saturation and spinodal chemical potential, respectively (Shen & Debenedetti 2001). The simulations explore the metastable range of density at the selected temperature, $\rho _L\in [\rho _{L\,spin}, \rho _{L\,sat}]=[0.44, 0.51]$. It may be worth stressing that on decreasing the metastability parameter the barrier height for nucleating vapour bubbles increases (e.g. the saturation line is approached). As a main goal, the effect of changing the wall wettability is discussed and the potential effect of liquid stratification (depletion/accumulation) at the wall is investigated.

In order to correctly identify bubbles and evaluate the nucleation rate in the different conditions, it is crucial to properly distinguish vapour bubbles from subcritical embryos. Here the same technique (the string method; E, Ren & Vanden-Eijnden 2002; Bonella, Meloni & Ciccotti 2012) adopted by Gallo et al. (2018) is exploited, with a slight modification to take the wall into account. In a first step, the critical bubble is evaluated at given metastable conditions ($\rho _L$ and $\theta$). A spherical vapour bubble, with the corresponding density profile $\rho _{crit}(r)$ in the radial direction, is found. Then, the critical volume $V_c$ is rescaled to account for the contact angle using the geometrical function $\psi (\phi )$ described in § 2. This procedure provides the estimate of the critical volume to be used as reference size discriminating bubbles from embryos. A few dedicated computations to determine the critical bubble in the presence of the wall, adapting the string method approach to the case of heterogeneous nucleation, showed that the use of the geometric factor is in excellent agreement with the actual critical volume.

In the second step, a clustering algorithm detects the number of critical bubbles in the domain, as follows: (1) A cell flagging procedure selects the vapour cells based on density ($\rho <\rho _{cut}$ identifies vapour). The threshold density, $\rho _{cut}$, is chosen as the density where the critical profile $\rho _{crit}(r)$ exhibits the maximum slope, corresponding to the interface centre. It turned out that, in the range of metastable states analysed here, $\rho _{cut}$ is in the range $0.35\text {--}0.36$. A sensitivity analysis showed small effects of this threshold value on the results. (2) A region growing procedure clusters vapour cells by examining the neighbourhood of flagged cells. The procedure is iterated until all flagged cells are processed. (3) The volume $V_k$ of the $k$th cell aggregate (cluster), the sum of cell volumes belonging to the cluster, is used as order parameter identifying supercritical bubbles when $V_k \ge V_c$.

Before discussing the nucleation rate, snapshots at two different conditions are shown in figures 7 and 8, at $\rho _L=0.475$ and $\rho _L=0.485$, respectively, to illustrate the evolution at different degree of metastability. Starting from a homogeneous liquid mother phase, thermal fluctuations produce regions of low density with a complex shape. Although energy considerations may suggest that vapour bubbles should preferentially form close to the walls, in both cases vapour embryos are observed all over the domain. After reaching the critical size, a complex dynamics is initiated during which collapses and/or coalescence events take place. In the long run, the bubbles surviving in the condition with smaller metastability, figure 8, are all found attached to the wall. At large degree of metastability (small nucleation barrier) several bubbles persist instead in the bulk liquid.

Figure 7. Snapshots during the nucleation process at thermodynamic condition $\rho _L = 0.475$, taken at times $t = 60\,000$, 80  000, 240  000 and 264  000 (a–d, respectively). The vapour bubbles are represented by density isosurfaces. Close to the bottom wall we report a slice representing the density field with the colour scale and where the regions with high density gradients are indicated by the black isolines.

Figure 8. Snapshots during the nucleation process at thermodynamic condition $\rho _L = 0.485$, taken at times $t=180\,000$, 190  000, 280  000 and 360  000 (a–d, respectively).

In general, bubble number and dimensions depend on the initial metastability, as shown in figure 9, reporting the number of bubbles formed at the wall for the two previously mentioned thermodynamic conditions. The time instants singled out in the panels of figures 7 and 8 are indicated with circles and identified by the letters ${a,b,c,d}$ in figure 9. The evolution shows three main stages.

Figure 9. Time evolution of the number of wall-attached supercritical bubbles at two different thermodynamic conditions $\rho _L=0.475$ (a) and $\rho _L=0.485$ (b). In both panels, the insets report a zoom in the time interval where a linear fitting was applied for calculating the nucleation rate.

Initially the bubble number increases almost linearly with time (constant production rate). At the lower degree of metastability ($\rho _L=0.485$) the constant rate regime follows after a relatively long incubation period with almost no bubble. This waiting time is related to the mean first passage time (Kramers 1940) of classical theory, where nucleation is described as Brownian wandering on a free energy landscape, until the walker escapes the barrier to nucleate the bubble. The constant rate stage persists until a maximum number of bubbles accumulates in the domain (configurations denoted by the letter ${b}$). The higher the degree of metastability, the lower is the energy barrier and, as a consequence, the larger is the maximum number of bubbles in the domain.

The second stage consists of the expansion–coalescence dynamics where bubbles increase their size evolving towards equilibrium while partially coalescing with neighbouring ones. During this phase, the newly formed embryos typically collapse, since the liquid gets compressed by the enlarging vapour phase, given the constraints of constant mass and volume. This kind of event occurs more frequently for the more populated condition at larger degree of metastability.

Finally, during the third stage, the evolution slows down while reaching a more stable condition, where a small number of large vapour bubbles are in equilibrium with the surrounding compressed liquid (see configurations $d$ in the figures.

In the present numerical experiments, the nucleation rate $J$ is accessed during the first stage. The nucleation rate expresses the frequency of bubble formation normalized by the system volume (in the context of homogeneous nucleation) or by the solid surface area (in heterogeneous conditions). From an operative standpoint, there are two main procedures to evaluate the bubble formation frequency: (i) When the simulated system is small and only a single nucleation event occurs, as often happens in MD simulations, the waiting time for the appearance of the bubble is measured. A large number of independent simulations are used to properly sample the statistics of the phenomenon, and the inverse of the mean waiting time is used as bubble nucleation frequency (Menzl et al. 2016). (ii) When the simulated system is large and the nucleation of a multitude of bubbles can be sampled within a single simulation, the frequency is evaluated as the slope of the curve expressing the time evolution of the number of bubbles in the system (Yasuoka & Matsumoto 1998; Diemand et al. 2014). This linear fitting, expressing the number of bubbles formed per unit time, is performed at the beginning of the nucleation process, and the so-called steady-state rate is measured.

Here the second strategy has been exploited, and the slopes have been evaluated as in the inset of figure 9. The calculated nucleation rates are shown in figure 10(a) in comparison with CNT predictions. It might be worth noting that CNT in the $\mu V \theta$ ensemble (§ 2.1) is here used as a reference theory in the entire range of metastability conditions. The difference between the unknown actual rate and the CNT prediction is expected to vanish asymptotically at larger density, close to saturation. The present results are in line with this expectation, with the convergent trend apparent in the figure. Far from saturation, only homogeneous nucleation has been plentifully analysed, finding that the results spread over approximately eight orders of magnitude, at high temperature, and differ by up to 20 orders from CNT (see Diemand et al. (2014) for a summary of the literature data). To the best of our knowledge, quantitative simulations of heterogeneous bubble nucleation rates are not available for a direct comparison. Only few qualitative attempts can be found in Nagayama, Tsuruta & Cheng (2006). The available MD simulations by Novak et al. (2007) have instead been performed in an isothermal–isostress ensemble (${N P_{zz} \theta }$).

Figure 10. Comparison of the nucleation rates between the fluctuating hydrodynamics (FH) numerical results (square symbols) and the classical nucleation theory (CNT) prediction by Blander & Katz (1975) (triangles) in the ${\mu V T}$ setting. (a) Nucleation rates at different degrees of metastability, and constant contact angle $\phi ={\rm \pi} /2$. (b) Different wall wettabilities at thermodynamic condition $\rho _L= 0.48$.

The standard application of CNT refers to the grand canonical (${\mu V \theta }$) context, where the liquid bulk density remains unaltered along the nucleation process in macroscopic systems. The total mass and energy constraints ($NVE$ ensemble) used in the present simulations might lead to crucial differences, particularly when the system is small or extremely crowded with several bubbles. As shown below, this explains the discrepancy between simulation results and predictions found in figure 10.

By elaborating on the $NVE$ model recalled in § 2.2, a system partially filled with a given quantity of vapour to account for nucleation in the presence of multiple bubbles can be considered. In this case, the critical bubble radius and corresponding energy barrier are found by solving equations (2.6) specialized as

(9.5)\begin{equation} \left.\begin{array}{c@{}} R^* = \dfrac{2\gamma_{LV}}{p_V (\rho_V, \theta) - p_L (\rho_L, \theta) } , \\ \mu_V (\rho_V, \theta ) = \mu_L (\rho_L, \theta) , \\ E = U_L (\rho_L, \theta) + \tilde{U}_V + U_{bubble} (\rho_V, \theta) + \tilde{E}_c + \gamma_{LV} A^* , \\ M = \tilde{M}_V + V^* \rho_V + (V -V^* - \tilde{V}_V)\rho_L . \end{array}\right\} \end{equation}

In the above equations, $R^*$, $A^*$ and $V^*$ are the bubble critical radius, area and volume, respectively. In the case of a neutrally wettable solid wall ($90^{\circ }$ contact angle), for example, $A^* = 2{\rm \pi} R^{* 2}$ and $V^* = \frac 23 {\rm \pi}R^{*3}$. Moreover, $U_{bubble}$ is the internal energy of the critical vapour bubble, evaluated through the equation of state, as $u_b(\rho _V,\theta ) V^*$. All terms highlighted with a tilde are free parameters that allow one to take into account the vapour already formed in the system. Specifically, $\tilde {U}_V$ is the total internal energy of the vapour phase, $\tilde {E}_c$ is the (total) surface energy, $\tilde {M}_V$ is the total mass of vapour and $\tilde {V}_V$ its total volume. In comparing numerical results with the predictions of this extended $NVE$ model, the free parameters are measured from the simulation along the evolution. The time-dependent nucleation rate $J_{NVE}(t)$ is evaluated by solving system (9.5) at each instant. Finally, the number of bubbles is obtained by integrating the nucleation rate,

(9.6)\begin{equation} N_B(t) = 2S \int_{t_0}^t J_{NVE}(\tau)\,\textrm{d}\tau , \end{equation}

where $2S$ is the area of the solid walls and $t_0$ the nucleation waiting time, also extracted from the simulation.

The remarkable agreement shown in figure 11 confirms that the constraints on mass and energy strongly affect the nucleation rate, and fully justifies the order-of-magnitude difference in the rate when fluctuating hydrodynamics results are compared with standard CNT in the ${\mu V\theta }$ setting. The modified $NVE$ description envisaged here hinges on input data from the simulation. As such, it is not a closed model of nucleation. Nevertheless, it enlightens the physical mechanisms operating in the system, also capturing the progressive reduction of nucleation rate, due to the progressive filling of the system with nucleated vapour. On the contrary, the constant rate predicted by the CNT ${\mu V\theta }$ leads to a linearly increasing number of bubbles (blue dashed line in figure 11), out of scale with respect to the simulation results.

Figure 11. Comparison of the number of bubbles predicted by the CNT $NVE$ model proposed in (9.5) (solid black curve) with the FH simulation data (solid red curve). The confidence interval in the simulation results (light red band) corresponds to a change of the threshold density, $\rho _{cut}$, in the bubble tracking procedure in the range 0.32–0.38. The blue dashed line corresponds to the linearly growing number of bubbles predicted by the more common CNT model in the ${\mu V \theta }$ setting.

Wall wettability in the specific thermodynamic condition $\rho _L=0.48$ is investigated by changing the contact angle $\phi$, (3.8). The numerical simulations span a range of hydrophobic ($\phi >90^\circ$) and hydrophilic ($\phi <90^\circ$) conditions. The results shown in figure 10(b) are consistent with the expected behaviour, based on energy barrier considerations: the nucleation rate increases when the contact angle is increased, i.e. the more hydrophobic the surface, the more likely is vapour formation close to the wall. Quantitatively, the CNT prediction shows a more pronounced effect with respect to those observed in the numerical simulations.

For a qualitative illustration, snapshots during nucleation at weakly hydrophilic and weakly hydrophobic walls are reported in figure 12. It is clear that the number of embryos formed at the wall is larger for hydrophobic surfaces. Eventually, figure 13, at moderate and strong hydrophilicity, nucleation is mostly observed in the bulk in the same manner as described for homogeneous nucleation.

Figure 12. Snapshots during the nucleation process at the thermodynamic condition $\rho _L=0.48$, and different wall wettabilities (as marked): (a,b) reports the case of a moderately hydrophilic wall; (c,d) a moderately hydrophobic one. (a,b) $\phi =85^{\circ}$ and (c,d) $\phi =96^{\circ}$.

Figure 13. Snapshots of the bubbles formed at the metastable condition $\rho _L=0.48$ as a function of the contact angle in a range of moderately hydrophilic conditions (as indicated). (a) $\phi =80^{\circ}$, (b) $\phi =70^{\circ}$, (c) $\phi =60^{\circ}$ and (d) $\phi =45^{\circ}$.

In quantitative terms, figure 14 reports the number of vapour bubbles detected in the domain for several hydrophilic wettability levels. Figure 14 (a) and (b) compare the number of wall-attached and the total number of bubbles, respectively. When the contact angle decreases, heterogeneous nucleation events are drastically reduced. The sensitivity of bulk nucleation to contact angle observed in the data can be explained by the enforced mass and volume constraints, which makes pressure and density responsive to the overall volume of vapour phase, which is influenced by the number and size of wall bubbles. Separating surface from bulk nucleation events allows the surface and bulk nucleation rates to be compared. By comparing the slope in the constant-rate regime, the homogeneous nucleation rate is found to be constant within the present statistical accuracy, independently of wall wettability. The heterogeneous (surface) nucleation rate reduces at increasing hydrophilicity and, when $\phi < \phi _{crit}$ with $\phi _{crit} \approx 60^{\circ }$, almost all nucleation events take place in the bulk.

Figure 14. Time evolution of the number of supercritical bubbles at different contact angles $\phi$. (a) The number $N_B$ of bubbles nucleated at the solid wall. (b) The total number $N_B^{TOT}$ of bubbles detected in the entire fluid domain. In the inset, the number of bubbles far from the wall, $N_B^{BULK} = N_B^{TOT} - N_B$, is reported for the reader's convenience.

This behaviour can be explained by the following considerations. (1) From a thermodynamic point of view, the hydrophilic nature of the surfaces involves accumulation of liquid to the wall, strongly preventing the formation of vapour nuclei in these zones. (2) The homogeneous nucleation barriers are comparable to the heterogeneous one close to moderately and strong hydrophilic surfaces, but the number of nucleation sites in the bulk volume is clearly larger than those on the walls. As a consequence, the homogeneous nucleation probability is higher. A similar behaviour was also observed in several MD simulations of heterogeneous nucleation, both on flat solid walls (Nagayama et al. 2006; Chen et al. 2018; Marchio et al. 2018) and close to structured surfaces (She et al. 2016). It may be worth stressing that these results contradict the widespread opinion based on CNT (see § 2) that solid surfaces should always act as bubble catalysts.

One of the major advantages of the mesoscale approach considered here consists in the possibility to follow the complete bubble dynamics, from the inception of phase transition to the late stage of macroscopic bubble evolution. The latter is analysed in the following. The cluster analysis described at the beginning of this section enables tracking of individual bubbles. Figure 15 concerns the volume histories of bubbles that survive for the entire simulations. There, several curves (red and blue lines) feature sudden volume jumps that correspond to coalescence events. One of those, see the event labelled $b$ captured in the snapshots of figure 16, involves the coalescence of two wall-attached bubbles. The growth rate of the resulting bubble increases after coalescence, probably due to the fusion process, which, by decreasing the interfacial curvature, tends to favour the expansion.

Figure 15. Time evolution of the volume of individual vapour bubbles at the specific thermodynamic condition $\rho _L=0.475$ and contact angle $\phi ={\rm \pi} /2$. The sudden volume jumps represent bubble–bubble coalescence events. The letters ${a,b,c,d}$ refer to the time instants of the snapshots reported in figure 16.

Figure 16. Snapshots during the nucleation process at the thermodynamic condition $\rho _L=0.48$ and contact angle $\phi ={\rm \pi} /2$ taken at times $t=90\,000$, 165  000, 235  000 and 258  000 (a–d, respectively). The slices help to recognize the coalescence events.

The analysis is completed by the time evolution of mean temperature and pressure inside the bubbles (figure 17). The saturation pressure and initial temperature are indicated as references. Initially, both temperature and pressures are quite noisy, with a large variance due to the small bubble dimensions. The temperature, which initially fluctuates around the liquid temperature, starts decreasing during the expansion stage. Similarly the pressure, which initially fluctuates around the saturation value, decreases together with the temperature. All the bubbles experience a similar trend, with pressure that differs from the saturation pressure up to $10\,\%$.

Figure 17. (a) Time evolution of the mean temperature inside the individual bubbles of figure 16 (letters ${a,b,c,d}$ respectively). (b) Time evolution of the mean pressure inside the same bubbles. Both the initial temperature and the saturation pressure are reported as a reference.

10. Conclusions

This paper deals with the long-standing issue of heterogeneous nucleation in metastable liquids. Common approaches deal with the problem using heuristic models in the context of classical fluid mechanics or, more rigorously, simplified approaches like the classical nucleation theory (CNT) or different kinds of more or less sophisticated simulations based on atomistic descriptions of the fluid. Here a different path is followed, addressing the nucleation process in the context of a formulation where thermal fluctuations are included in the diffuse interface model in a physically consistent way. This leads to a system of SPDEs that extends the capillary Navier–Stokes equations and includes a suitable random forcing to model thermal noise. The diffuse interface model may be seen as originating from the gradient expansion of more general non-local functionals in the context of density functional theory (DFT). The result is an extension to heterogeneous nucleation of a previous proposal based on fluctuating hydrodynamics that was originally used to deal with nucleation in homogeneous conditions.

Technically, the issue was developing a theory of fluctuations in the presence of walls. In this approach, the appropriate boundary condition describing the hydrophobic or hydrophilic walls enters the Onsager operator that determines the noise intensity. Although the theory is in several respects basically formal, specific examples allowed the noise intensity operator to be explicitly determined. As a consequence, it is shown that there is no need to modify the structure of the white noise when boundary conditions are suitably included in the relevant operators. The issue is discussed also for the discrete system, showing how the discrete model can be consistently derived from the original system of SPDEs. The equivalence of the discrete model so derived with previous works that addressed the problem in purely numerical terms by modifying the noise structure at the wall is illustrated.

The proposed model is exploited to produce numerical simulations of bubble nucleation at a solid wall with various degrees of wettability. By inducing layering effects in the liquid close to the boundary, the wettability affects the heterogeneous nucleation rate. At a qualitative level, the numerical data reproduce the predictions of CNT. There are, however, significant differences with respect to the predictions of classical CNT extended to heterogeneous conditions. The usual theory deals with the ${\mu V \theta }$ ensemble and provides an estimate for the mean nucleation rate by considering a single bubble developing in the metastable liquid. Direct comparison with this approach produces a small but significant discrepancy in the rate, which decreases as the metastability level is reduced. In order to better validate the present model, it was found crucial to take into account two basic differences with the assumptions of the classical theory: (i) the simulations were performed keeping the total system energy constant, making the ${NVE}$ ensemble more appropriate; and (ii) the system evolves in time, due to the accumulation of vapour in the closed system. For this reason we adapted the ${NVE}$ CNT to make it able to predict the nucleation rate in a time-dependent state. The comparison of this extended model turns out to be excellent. Although care must be exerted when comparing heterogeneous and homogeneous nucleation rates, since the former is a surface process and the latter a volume process, the present results suggest that, already for moderately hydrophilic walls, homogeneous nucleation may still result in the most effective mechanism, due to the liquid layering at the wall.

The favourable computational cost of this methodology encourages its exploitation for more complex conditions of nucleation in the presence of complex geometries, mean flow and dissolved gas, all aspects left for future works.

One of the potentially interesting novelties obtained as a byproduct of the present analysis concerns a systematic description of the interaction between a capillary fluid and a solid wall in the context of a diffuse interface approach. The equation introduced to model the wall wettability is related to the boundary condition introduced by Cahn in the context of the Cahn–Hilliard model. Here such a relationship is extended to take into account the thermodynamics of a non-isothermal liquid–vapour system showing how to determine the wall contribution to the free energy consistently with the specific equation of state used to describe the bulk fluid. Based on the resulting expression of the wall free energy, it is also shown how to include non-equilibrium effects in the contact angle relaxation dynamics.

Although the present contribution may constitute a basic step to develop mesoscopic numerical models of nucleation for complex systems able to bridge the gap between nucleation and macroscopic bubble dynamics, several important ingredients still need to be considered. Among the most important in view of practical applications is the extension of the numerics to complex geometries in order to take into account wall roughness effects, which are known to be crucial in most cases. The presence of dissolved gas is almost unavoidable in experiments and very important for applications. Both effects are in principle easily included in the model, although at the price of some additional work.

At the microscopic level, slippage is known to occur at liquid–solid interfaces. At mesoscopic scale, this effect is described by introducing a phenomenological slip length in the Navier equation that replaces the no-slip condition at the wall. It turns out that, apart from special systems, such as flow in carbon nanotubes, the slip length is of the order of less than one nanometre. Its effect on nucleation may, however, need to be included in the model. In principle, this can be easily considered, including also non-equilibrium effects on the contact angle dynamics. This will, however, require dedicated consideration, since any new dissipative effect is expected to alter the explicit form of the fluctuation–dissipation balance that determines the noise intensity.

At a more fundamental level, the limit of the model is the very assumption of the square-gradient form for the free energy functional. The consequence is a monotonic density profile, which, although describing the correct trend in a coarse-grained sense, does not capture the details of the density oscillations that are known to occur in liquids both experimentally, via X-ray scattering, and numerically, via MD simulations. More complex functionals derived in the context of DFT can capture these effects, which could, in principle, affect the nucleation rate. A common characteristic of DFT is, however, the non-local nature of the functionals. If included in the model in the spirit of a generalized dynamic DFT, non-locality would make the computations much more demanding, spoiling the the ability of the model to reach macroscopic scales. In this respect, the opinion of the authors is that the square-gradient approximation could be considered as a very good compromise between the ability to account for even the finest details and computational feasibility.