Linear stability theory and molecular simulations of nanofilm dewetting with disjoining pressure, strong liquid–solid slip and thermal fluctuations

Abstract

The dewetting of thin nanofilms is significantly affected by thermal fluctuations, liquid–solid slip and disjoining pressure, which can be described by lubrication equations augmented by appropriately scaled noise terms, known as stochastic lubrication equations. Here molecular dynamics simulations along with a newly proposed slip-generating method are adopted to study the instability of nanofilms with arbitrary slip. These simulations show that strong-slip dewetting is distinct from weak-slip dewetting by faster growth of perturbations and fewer droplets after dewetting, which cannot be predicted by the existing stochastic lubrication equation. A new stochastic lubrication equation considering the strong-slip boundary condition is thus derived using a long-wave approximation to the equations of fluctuating hydrodynamics. The linear stability analysis of this equation, i.e. surface spectrum, agrees well with molecular simulations. Interestingly, strong slip can break down the usual Stokes limits adopted in weak-slip dewetting and bring the inertia into effect. The evolution of the standard deviation of the film height $W^2(t)={\overline {h^2}-{\bar {h}}^2}$ at the initial stage of the strong-slip dewetting is found to be $W\sim t^{1/4}$ in contrast to $W\sim t^{1/8}$ for the weak-slip dewetting.

1. Introduction

Nanometric thin liquid films deposited on substrates exist in a host of applications such as in lubricants (Jhon et al. 1999), coatings (Weinstein & Ruschak 2004) and microfluidics (Stone, Stroock & Ajdari 2004). The reliability of those applications depends heavily on our understanding of their stability mechanism, which is usually investigated in the context of thin-film flows (Oron, Davis & Bankoff 1997; Craster & Matar 2009). Thin-film flows are characterised by the disparity of length scale in different dimensions, i.e. the ratio of film height $h$ to characteristic film length $\lambda$ is very small: $\chi = h/\lambda \ll 1$. This allows the adoption of a long-wave theory to derive lubrication equations from the full governing equations and boundary conditions, reducing the dimensionality and complexity of the problems (Oron et al. 1997; Craster & Matar 2009).

Polymeric or metallic films on substrates with thicknesses below 100 nm have been observed to undergo spontaneous rupture and dewetting (Xie et al. 1998; Seemann, Herminghaus & Jacobs 2001; Becker et al. 2003; González, Diez & Sellier 2016). The dewetting mechanism in these films may be complicated due to the contamination of defects in the liquid. However, the primary dewetting mechanism for homogeneous liquid films is usually called spinodal dewetting. In this process, disjoining pressure, as a result of intermolecular forces between liquid and solid, leads to the instability of films. Basically, from the classical perspective, thermally excited capillary waves can be amplified by the disjoining pressure, but in competition with the restoring force of surface tension, such that disturbances above a critical wavelength can grow and lead to film rupture (Vrij & Overbeek 1968).

For interfacial flows at the nanoscale, thermal fluctuations can play an important role in the instability process (Moseler & Landman 2000; Grün, Mecke & Rauscher 2006; Zhang, Sprittles & Lockerby 2019). Thermal fluctuations can spontaneously generate thermal capillary waves (TCWs) and roughness on the free surface of a liquid film at rest. The magnitude of thermal roughness is usually proportional to thermal length $\sqrt {k_BT/\gamma }$ ($k_B$, $T$ and $\gamma$ are the Boltzmann constant, temperature and surface tension, respectively) (Buff, Lovett & Stillinger 1965; MacDowell 2017). Though the roughness is small and usually on the scale of nanometres, it becomes comparable to the size of films when the film thickness goes down to several nanometres. Note that micrometre roughness can also be obtained and, thus, observed optically in real space using ultra-low surface tension mixtures ($\gamma \sim 10^{-6}\ {\rm N}\ {\rm m}^{-1}$) (Aarts, Schmidt & Lekkerkerker 2004). In the equilibrium state, the amplitude of TCWs, known as the static spectrum, can be described by the renowned capillary wave theory (Buff et al. 1965; Höfling & Dietrich 2015; MacDowell 2017; Höfling & Dietrich 2024). Recently, an extension of the capillary wave theory has been proposed utilising a Langevin equation to describe the transient dynamics of non-equilibrium TCWs and their approach to thermal equilibrium (Zhang, Sprittles & Lockerby 2021). This advancement has led to the identification of a universality class governing the roughening behaviour of film surfaces (Zhang et al. 2021).

The increasing importance of thermal fluctuations as the film height decreases may make the deterministic description of hydrodynamics at the nanoscale break down. For example, the breakup of liquid nanojets in molecular dynamics (MD) simulations (Moseler & Landman 2000; Zhao, Sprittles & Lockerby 2019) and experiments (Hennequin et al. 2006) shows a double-cone rupture profile, in contrast to the long-thread profile predicted by the deterministic lubrication equation (Eggers & Dupont 1994). Moseler & Landman (2000) pioneered in showing that the deficiency of this deterministic lubrication equation for describing nanojet dynamics can be remedied by adding a noise term of appropriate strength to the equation, which leads to a stochastic lubrication equation for nanojets. Eggers (2002) later showed that the evolution of the minimum neck radius is accelerated by thermal fluctuations, leading to $h_{min}\propto (t_0-t)^{0.418}$ (where $t_0$ is the rupture time) in contrast with $h_{min}\propto (t_0-t)$ for the deterministic pinching. This noise-dominated breakup for nanojets has been observed in experiments using ultra-low surface tension mixtures (Hennequin et al. 2006).

For nanofilm rupture, Grün et al. (2006) and Davidovitch, Moro & Stone (2005) independently derived the same stochastic lubrication equation for liquid films on no-slip substrates. The numerical solution to this equation (Grün et al. 2006) can resolve the discrepancy in dewetting time between experimental results (Becker et al. 2003) and the solution to the deterministic counterpart. Subsequently, the rupture of thin films with the effects of thermal fluctuations has been widely investigated by numerical solutions to the stochastic film equation (Grün et al. 2006; Nesic et al. 2015; Diez, González & Fernández 2016; Durán-Olivencia et al. 2019; Shah et al. 2019; Zitz, Scagliarini & Harting 2021). These studies have consistently demonstrated that thermal fluctuations indeed accelerate the rupture process. The application of this stochastic film equation is not restricted to nanofilm dewetting. It has been extended to study, for example, nanodroplet spreading under an elastic sheet (Carlson 2018), curvature-induced film drainage (Shah et al. 2019) and mediated diffusion of particles confined in channels with a fluctuating wall (Marbach, Dean & Bocquet 2018).

In addition to numerical solutions, linear stability analyses of the stochastic film equation have also been studied a lot, which allows us to obtain the evolution of the capillary spectra of surface waves (Mecke & Rauscher 2005; Fetzer et al. 2007; Zhang et al. 2019; Zhao et al. 2019). The analytical spectra show thermal fluctuations can massively amplify the growth of waves, shift the critical wavenumber to a larger value and cause the dominant wavelength to evolve in time (in contrast to a constant value predicted by the deterministic lubrication equation). These interesting findings have been validated both in MD simulations (Zhang et al. 2019) and experiments (Fetzer et al. 2007).

The original stochastic film equation mentioned previously adopts the classical no-slip boundary condition. As the flow scale reaches nanometres, surface effects such as liquid–solid slip can have significant effects on flow behaviours (see the reviews by Lauga, Brenner & Stone 2007; Bocquet & Charlaix 2010). Obviously, nanofilm flows can be significantly affected by slip as well, since the ratio of slip length $b$ to the film thickness $h$ can get close to unity or even much larger than unity (Bäumchen & Jacobs 2009). In fact, in the deterministic setting, the introduction of slip to the deterministic film equation has been studied extensively for various phenomena, such as in coating a plate (Liao, Li & Wei 2013), droplet spreading (Savva, Kalliadasis & Pavliotis 2010), film rupture (Martínez-Calvo, Moreno-Boza & Sevilla 2020) and falling films down a slippery plate (Ding & Wong 2015). However, only recently has the no-slip stochastic film equation been generalised to consider slip (Zhang, Sprittles & Lockerby 2020), which is non-trivial and requires the usage of the Green–Kubo-type expression (Bocquet & Barrat 1994) that relates slip length to the random stress tensor at the wall. The derived slip equation is validated by the well-controlled molecular simulations (Zhang et al. 2020).

However, the derived slip equation (Zhang et al. 2020) is limited to the case of weak slip $b/h \approx 1$. In many cases, the slip length can be as large as micrometres so that $b/h\gg 1$, such as flow over graphene sheets (Falk et al. 2010), flow over engineered hydrophobic materials (Rothstein 2010) and flow over substrates in presence of gas cavities and surface nanobubbles (Lohse & Zhang 2015). For polymer liquids, increasing molecular weights can also increase the slip length up to micrometres (Bäumchen, Fetzer & Jacobs 2009). In fact, the dewetting of polymeric films on dodecyltrichlorosilane substrate (DTS), where slip length can be up to $1~{\rm \mu}$m, has been examined extensively in the deterministic framework (Kargupta, Sharma & Khanna 2004; Fetzer et al. 2005; Münch, Wagner & Witelski 2005; Bäumchen et al. 2014). Notably, different levels of slip give rise to different deterministic lubrication models (Münch et al. 2005). However, as we mentioned earlier, the effects of thermal fluctuations on nanofilms are significant, so a generalisation of our current understanding of the instability of strong-slip films to consider thermal fluctuations is essential.

So far, experimental studies on the effects of thermal fluctuations on thin-film flows are limited due to the technical difficulties associated with the spatiotemporal scale. Though the experiments of the dewetting of polymer nanofilms on no-slip SiO$_{2}$-coated silicon wafers have demonstrated the great effects of thermal fluctuations (Fetzer et al. 2007) on the growth of surface perturbations, dewetting of polymer films with strong slip and thermal fluctuations have not been considered in experiments. As such, MD simulations are a natural and convenient tool to investigate the thin-film problem at the nanoscale as thermal fluctuations are inherent in MD simulations.

In this work, MD simulations are employed to simulate the rupture of nanofilms on substrates with a strong slip, in comparison with the small-slip rupture. A new simulation strategy is proposed to generate strong slip in molecular simulations since the classic molecular simulations are limited to weak slip. We obtain the evolution of film surface spectra, rupture time and number of droplets after film rupture from molecular simulations. A new stochastic lubrication equation considering the strong slip is derived from fluctuating hydrodynamics (FH). A linear stability analysis of this stochastic lubrication equation leads to the analytical spectra which are validated against the MD results to establish the applicability of the new theory to predict future experiments.

This paper is organised as follows. In § 2, the stochastic lubrication equation for the strong-slip dewetting is derived from the equations of FH using a long-wave approximation. In § 3, a linear stability analysis of the newly derived stochastic equation is performed to obtain the surface spectrum. In § 4, molecular simulations of the rupture of nanofilms with the method to generate strong slip are presented. Section 5 compares the new model with molecular simulation results, and discusses new findings. In § 6, we summarise the main contributions of this work and outline future directions of research.

2. Stochastic lubrication equation for films with strong slip

2.1. Governing equations and boundary conditions

As shown in figure 1, a molecularly thin liquid film is deposited on a solid surface and destabilised by both disjoining pressure $\phi$ and thermal fluctuations $\psi$. The film is quasi-two-dimensional (quasi-2-D) confined by its size ($L_x$, $L_y$, $h$) where $L_x\gg L_y$ and $L_x\gg h$. Without slip, the stochastic thin film equation is derived in detail by Grün et al. (2006) using a long-wave approximation ($\chi =h_0/\lambda \ll 1$) to FH (Landau & Lifshitz 1959). The no-slip stochastic equation has been extended to consider weak slip $b\sim O(h)$ (Zhang et al. 2020). Here we present the derivation of a new equation considering $b\gg h$.

Figure 1. Sketch of a (quasi-two-dimensional) molecularly thin liquid film on a slippery substrate. Here $h(x,t)$ is the film thickness, $\lambda$ is the characteristic length, $u$ is the horizontal velocity, $\phi$ is the thermal fluctuations, $\psi$ is the disjoining pressure and $b$ is the slip length. The film has a small depth $L_y$ into the page.

The governing equations for this problem are given by equations of FH, where thermal fluctuations are modelled by an additional random stress tensor. The (incompressible) continuum equation and momentum equations are

(2.1)\begin{gather} {\partial_x}u+{\partial_z}w=0, \end{gather}

(2.2)\begin{gather} \rho \left({\partial_t}u+u{\partial_xu}+w{\partial_zu} \right)=-{\partial_x p}+\mu \left( {\partial_{xx}u}+{\partial_{{zz}}u} \right) +{\partial_x}{{\psi }_{xx}}+{\partial_z}{{\psi }_{zx}}, \end{gather}

(2.3)\begin{gather} \rho \left( {\partial_t}w+u {\partial_x w}+w{\partial_zw} \right)=-{\partial_z p}+\mu \left( {\partial_{{xx}}w}+{\partial_{{zz}}}w \right) +{\partial_x}{{\psi }_{xz}}+{\partial_z}{{\psi }_{zz}}. \end{gather}

Here $u$ and $w$ are the $x$ and $z$ components of velocity, and $\psi$ is a 2-D random stress tensor with zero mean and covariance given by

(2.4)\begin{equation} \left\langle {{\psi }_{ij}}(\boldsymbol{x},t) {{\psi }_{lm}}(\boldsymbol{x}',t') \right\rangle= \frac{2\mu {{k}_{B}}\theta}{L_y}\left( {{\delta }_{il}}{{\delta }_{jm}}+{{\delta }_{im}}{{\delta }_{jl}} \right) \delta ({x}-{x}')\delta ({z}-{z}')\delta(t-t'). \end{equation}

Here $\theta$ is the temperature. The factor $1/L_y$ appears because the films are quasi-2-D ($L_y \ll L_x$), allowing all variables of interest to be averaged over the $y$ direction (Zhang et al. 2020).

For boundary conditions, at $z=h$, we have the dynamic condition and kinematic condition:

(2.5a,b)\begin{equation} (\boldsymbol{\sigma}+\boldsymbol{\psi})\boldsymbol{\cdot}\boldsymbol{n}=-[\gamma\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{n}+\phi(h)]\boldsymbol{n}, \quad {\partial_t h}+u{\partial_x h}=w, \quad \text{at } z=h. \end{equation}

Here $\boldsymbol {\sigma }=-p{{\delta }_{ij}}+\mu (\partial _{{x}_{j}}{u}_{i}+\partial _{{x}_{i}}{{u}_{j}} )$ is the hydrodynamic stress tensor (here $i=(x,z)$ and $j=(x,z)$), $\gamma$ is the surface tension, $\phi (h)$ is the disjoining pressure and $\boldsymbol {n}$ is the outer normal vector at the surface $\boldsymbol {n}={( -\partial _x h,1 )}/{\sqrt {1+{{(\partial _x h )}^{2}}}}$.

At $z=0$, the impermeable condition and Navier's slip boundary condition are separately given by

(2.6a,b)\begin{equation} w=0, \quad u=b\frac{\partial u}{\partial z}, \quad \text{at } z=0, \end{equation}

where $b$ is the slip length. The covariance of the random shear stress at the wall is given by (Bocquet & Barrat 1994; Zhang et al. 2020):

(2.7)\begin{equation} \left\langle\psi_{zx}{{|}_{z=0}}(x,t)\psi'_{zx}{{|}_{z=0}}(x',t')\right\rangle=\frac{2\mu {{k}_{B}}\theta}{b L_y}\delta (x-x')\delta (t-t'). \end{equation}

Navier's slip condition is chosen because it has been validated extensively by both experiments and MD simulations (see the reviews by Lauga et al. 2007; Bocquet & Charlaix 2010). However, there are many other forms of slip boundary conditions, as discussed by Sibley et al. (2015).

2.2. A long-wave approximation

To derive a lubrication equation, (2.1)–(2.6a,b) have to be scaled based on the dominant mechanism of momentum balance, which varies with the level of slip length (Münch et al. 2005). Classically, for the weak-slip case, the momentum balance happens in the horizontal direction $\partial _x p \sim \mu \partial _{zz}u$, which means $ph_0/(\mu u_0) \sim 1/\chi$, where $u_0$ is the characteristic velocity. At the free surface, surface tension has to be balanced with pressure $p=-\gamma \partial _{xx} h$, which means $\gamma /(\mu u_0)\sim 1/\chi ^3$. At the solid surface, the order of slip length is $b\sim h_0$. Therefore, the pressure term, surface tension term, and slip length will be scaled to $P=\chi ph_0/(\mu u_0)$, $\varGamma =\chi ^3 \gamma /(\mu u_0)$, and $B=b/h_0$.

However, for the strong-slip (including free slip) flow, the velocity profile essentially becomes uniform (plug flow) instead of being parabolic so that the momentum balance happens in the vertical direction $\partial _z p \sim \mu \partial _{zz}w$ (Münch et al. 2005), which leads to the scaling $ph_0/(\mu u_0) \sim \chi$ and $\gamma /(\mu u_0)\sim 1/\chi$. These scalings need the slip length to be strong and $b\sim h_0\chi ^{-2}$ (Münch et al. 2005).

As for the scaling of the random stress tensor, it may be scaled the same as the scale of their deterministic counterparts (Grün et al. 2006). For example, $\psi _{xx}\sim \mu \partial _x u \sim \mu u_0/\lambda$ and $\psi _{zz}\sim \mu \partial _z w \sim \mu u_0/\lambda$. Note that, in this way, two scaling exists for $\psi _{zx}$: $\psi _{zx}\sim \mu \partial _x w \sim \mu \chi u_0/\lambda$ and $\psi _{zx}\sim \mu \partial _z u \sim \mu u_0/(\chi \lambda )$. The former is used to keep $\psi _{zx }$ at a lower order. In summary, the following scalings are used:

(2.8)\begin{align} \left. \begin{gathered} X=\frac{x}{\lambda },Z=\frac{z}{{{h}_{0}}},H=\frac{h}{{{h}_{0}}},B=\frac{b}{{{h}_{0}}\chi^{-2}},U=\frac{u}{{{u}_{0}}},W=\frac{w}{\chi {{u}_{0}}},T=\frac{{{u}_{0}}t}{\lambda },\varGamma =\frac{\chi \gamma }{\mu {{u}_{0}}},\\ \left(P,\varPhi \right)=\frac{{{h}_{0}}}{\chi {{u}_{0}}\mu }\left(\,p,\phi \right),\left( {{\varPsi }_{xx}},{{\varPsi }_{zz}} \right)=\frac{\lambda }{{{u}_{0}}\mu }\left( {{\psi }_{xx}},{{\psi }_{zz}} \right),\left( {{\varPsi }_{xz}},{{\varPsi }_{zx}} \right)=\frac{{\lambda}}{{{\chi }}{{u}_{0}}\mu }\left( {{\psi }_{xz}},{{\psi }_{zx}} \right). \end{gathered} \right\} \end{align}

Adopting these scalings (2.8), the rescaled and dimensionless continuum and momentum equations are

(2.9)\begin{gather} {\partial_X U}+{\partial_Z W}=0, \end{gather}

(2.10)\begin{gather} \chi {Re}\left({\partial_T U}+U{\partial_X U}+W{\partial_Z U} \right)=-{{\chi }^{2}}{\partial_X P}+ {{\chi }^{2}}\partial_{{XX}}U+{\partial_{{ZZ}}}U+{{\chi }^{2}}{\partial_X {{\varPsi }_{xx}}}+{{\chi }^{2}}{\partial_Z {{\varPsi }_{zx}}}, \end{gather}

(2.11)\begin{gather} \chi {Re}\left({\partial_T W}+U{\partial_X W}+W{\partial_{Z} W} \right)=-{\partial_Z P}+{{\chi }^{2}}{\partial_{{XX}}W}+{\partial_{{ZZ}}W} + {{\chi }^{2}}{\partial_X {{\varPsi }_{xz}}}+{\partial_Z {{\varPsi }_{zz}}}. \end{gather}

Here the Reynolds number is ${Re}={\rho u_0 h_0 }/{\mu }$. As the characteristic velocity is $u_0\sim \chi {\gamma }/{\mu }$, ${Re}$ can be also written as ${Re}=\chi {\rho \gamma h_0 }/{\mu ^2}=\chi {La}$, where $La={\rho \gamma h_0 }/{\mu ^2}$ is the Laplace number. The ${La}$ is assumed to be of order one.

In terms of the rescaled boundary conditions,

(2.12a) $$\begin{gather} -P{\,+\,}\frac{1}{1{\,+\,}{{\chi }^{2}}{{\left( {{\partial }_{X}}H \right)}^{2}}}\left\{ 2\left[ 1{\,-\,}{{\chi }^{2}}{{\left( {{\partial }_{X}}H \right)}^{2}} \right]{{\partial }_{Z}} W{\,-\,}2{{\partial }_{X}}H\left( {{\partial }_{Z}}U+{{\chi }^{2}}{{\partial }_{X}}W \right){\,+\,}{{\chi }^{2}}{{\left( {{\partial }_{X}}H \right)}^{2}}{{\varPsi }_{xx}} \right. \nonumber\\ -\left.{{\chi }^{2}}{{\partial }_{X}}H\left( {{\varPsi }_{xz}}+{{\varPsi }_{zx}} \right)+{{\varPsi }_{zz}} \right\} =\varGamma \frac{{{\partial }_{XX}}H}{{{\left[ 1+{{\chi }^{2}}{{\left( {{\partial }_{X}}H \right)}^{2}} \right]}^{3/2}}}+\varPhi, \quad \text{at } z=H,\end{gather}$$

(2.12b) $$\begin{gather} 2{{\chi }^{2}}{\partial_X H}\left[ \left( {\partial_Z W}-{\partial_X U} \right)+\frac{1}{2}\left( {{\varPsi }_{zz}}-{{\varPsi }_{xx}} \right) \right]\nonumber\\ +\left[ 1-{{\chi }^{2}}{{\left({\partial_X H} \right)}^{2}} \right]\left( {\partial_Z U}+{{\chi }^{2}}{\partial_X W}+{{\chi }^{2}}{{\varPsi }_{zx}} \right)=0\ \text{at}\ z=H , \end{gather}$$

(2.12c)$$\begin{gather}{\partial_T H}+U{\partial_X H}=W ,\quad \text{at } z=H, \end{gather}$$

(2.12d)$$\begin{gather}W=0,\quad \text{at } z=0, \end{gather}$$

(2.12e)$$\begin{gather}U=\frac{B}{{{\chi }^{-2 }}}{\partial_Z U},\quad \text{at } z=0. \end{gather}$$

The rescaled equations can be approximately solved by the perturbation expansion of $U,W,P,\varPsi,H$:

(2.13)\begin{equation} \left( U,W,P,\varPsi,H \right)=\left( {{U}_{0}},{{W}_{0}},{{P}_{0}},{{\varPsi }_{0}},H_0 \right)+{{\chi }^{2}}\left( {{U}_{1}},{{W}_{1}},{{P}_{1}},{{\varPsi }_{1}},H_1 \right)+\cdots. \end{equation}

Then, at leading orders of governing equations, one can get

(2.14)\begin{gather} {\partial_{{ZZ}}U_0}=0 , \end{gather}

(2.15)\begin{gather} {\partial_{{ZZ}}W_0}-{\partial_Z P_0}+{\partial_Z {{\varPsi }_{zz0}}}=0, \end{gather}

(2.16)\begin{gather} {\partial_X {{U}_{0}}}+{\partial_Z {{W}_{0}}}=0. \end{gather}

For boundary conditions, their leading-order forms are

(2.17a)\begin{gather} -{{P}_{0}}+2\left( {\partial_Z {{W}_{0}}}-{\partial_X H_0}{\partial_Z {{U}_{0}}}\right)+{{\varPsi }_{zz0}}=\varGamma {{\partial_{XX} }}H_0+\varPhi, \quad \text{at } z=H, \end{gather}

(2.17b)\begin{gather} {\partial_Z {{U}_{0}}}=0, \quad \text{at } z=H, \end{gather}

(2.17c)\begin{gather} {{W}_{0}}={\partial_T {{H}_{0}}}+{{U}_{0}}{\partial_X {{H}_{0}}},\quad \text{at } z=H, \end{gather}

(2.17d)\begin{gather} W_0=0, \quad \text{at } z=0, \end{gather}

(2.17e)\begin{gather} \partial_Z U_0=0, \quad \text{at } z=0. \end{gather}

Using leading-order equations, we find

(2.18)$$\begin{gather} {{U}_{0}} \equiv {{U}_{0}}(X,T), \end{gather}$$

(2.19)$$\begin{gather}{{W}_{0}} =-{\partial_X {{U}_{0}}}Z, \end{gather}$$

(2.20)$$\begin{gather}{{P}_{0}}=-2{\partial_X U_0}-\varGamma{\partial_{XX}H_0}-\varPhi +{{\varPsi }_{zz0}} , \end{gather}$$

which leads to the first equation of the desired stochastic lubrication equation:

(2.21)\begin{equation} {\partial_T {{H}_{0}}}+{\partial_X \left( {{H}_{0}}{{U}_{0}} \right)}=0. \end{equation}

In the next order, the governing equation that will be used is,

(2.22)\begin{equation} {La}\left({\partial_T {{U}_{0}}}+{{U}_{0}}{\partial_X {{U}_{0}}}\right)={{{\partial}_{XX}}{{U}_{0}}}-{\partial_X {{P}_{0}}}\ + {{\partial_{ZZ} }}{{U}_{1}}+{\partial_X {{\varPsi }_{xx0}}}+{\partial_Z {{\varPsi }_{zx0}}} \end{equation}

and the boundary condition that will be needed is

(2.23)\begin{gather} -4{\partial_X H_0}{\partial_X U_0}+{\partial_X H_0}\left( {{\varPsi }_{zz0}}-{{\varPsi }_{xx0}} \right)+ {\partial_Z {{U}_{1}}}-H_0{{{\partial_{XX} }}U_0}+{{\varPsi }_{zx0}}=0, \quad \text{at}\ z=H, \end{gather}

(2.24)\begin{gather} {{U}_{0}}=B{\partial_Z {{U}_{1}}},\quad \text{at}\ z=0. \end{gather}

Integrating (2.22) from $0$ to $H_0$ and using boundary conditions (2.23) and (2.24) leads to

(2.25) \begin{align} & H_0\left[ {La}\left({\partial_T {{U}_{0}}}+{{U}_{0}}{\partial_X {{U}_{0}}}\right) \right] =H_0\left[ 3{{{\partial_{XX} }}{{U}_{0}}}+{\partial_X }\left( \varGamma {{{\partial }_{XX}}H_0}+\varPhi \right) \right]\nonumber\\ &\qquad +\left({\partial_Z {{U}_{1}}}+{{\varPsi }_{zx0}} \right){{|}_{Z=H_0}}-\left({\partial_Z {{U}_{1}}}+{{\varPsi }_{zx0}} \right){{|}_{Z=0}}+\int_{0}^{H_0}{{\partial_X }\left( {{\varPsi }_{xx0}}-{{\varPsi }_{zz0}} \right)\,{\rm d}Z} \nonumber\\ &\quad =H_0{\partial_X }\left(\varGamma{\partial_{XX}H_0}+\varPhi \right)+4{\partial_X }\left( H_0{\partial_X {{U}_{0}}} \right)\notag\\ &\qquad -\frac{{{U}_{0}}}{B}+{\partial_X}\int_{0}^{H_0}{\left( {{\varPsi }_{xx0}}-{{\varPsi }_{zz0}} \right)\,{\rm d}Z}-{{\varPsi }_{zx0}}{{|}_{Z=0}} . \end{align}

The Leibniz integral rule is used here. The covariance of $\varPsi _{zx0}{|}_{Z=0}$ is given by the Green–Kubo expression for slip length $\langle \varPsi _{zx}{{|}_{Z=0}}(X,T)\varPsi '_{zx}{{|}_{Z=0}}(X',T')\rangle =({2 {{k}_{B}}\theta }/{\chi ^2 \mu u_0 b L_y})\delta (X-X')\delta (T-T')$ (Zhang et al. 2020). In terms of the integral of the white noise in (2.25), as ${\varPsi }_{xx0}$ is uncorrelated with ${\varPsi }_{zz0}$, they can be combined as $\sqrt 2{\varPsi }_{xx0}$. Now using the theorem provided in the appendix of Zhang et al. (2020), the integral of the white noise is

(2.26)$$\begin{gather} {\partial_X }\int_{0}^{H_0}{\left( {{\varPsi }_{xx0}}-{{\varPsi }_{zz0}} \right)\,{\rm d}z} =\sqrt{2}{\partial_X }\int_{0}^{H_0}{ {{\varPsi }_{xx0}} \, {\rm d}z} \nonumber\\ =\sqrt{2}{\partial_X }\left( \sqrt{H_0}\mathrm{Re} \right) , \end{gather}$$

where $\langle {{\mathrm {Re}(X,T)}}{{\mathrm {Re}(X',T') }}' \rangle =({4 {{k}_{B}}\theta }/{\mu u_0 h_0 L_y})\delta ( X-X' )\delta ( T-T' )$.

Putting (2.21) and (2.25) together and returning to their dimensional forms, we arrive at the stochastic lubrication equation for films with strong slip (referred to as S-S model hereafter)

(2.27a)$$\begin{gather} {\partial_t h}+{\partial_x \left( hu \right)}=0, \end{gather}$$

(2.27b)$$\begin{gather}\rho \left({\partial_t u}+u{\partial_x u}\right)={\partial_x}\left( -\gamma {{{\partial_{xx} }}h}+\phi \right)+\frac{4\mu }{h}{\partial_x }\left(h{\partial_x u}\right)-\frac{\mu }{h}\frac{u}{b}+\frac{1}{h}\left[ 2 {\partial_x }\left( \sqrt{h}{{\xi }_{1}} \right)-\frac{\sqrt b}{b}{{\xi }_{2}}, \right] \end{gather}$$

where the covariance of the noise term is $\langle {{\xi _i(x,t) }}{{\xi _j (x',t')}} \rangle =({2\mu {{k}_{B}}\theta }/{L_y})\delta _{ij}\delta ( x-x' )\delta ( t-t' )$.

The above-derived stochastic lubrication equation validates for the strong slip length on the order of $b\sim h_0\chi ^{-2}$. In terms of the weak slip $b\sim h_0$, a similar long-wave approximation to FH equations (see Zhang et al. (2020) for details) can result in the weak-slip stochastic lubrication equation (referred to as W-S model hereafter)

(2.28)\begin{equation} {\partial_t h}=\frac{1}{\mu }{\partial_x }\left[\left(\frac{1}{3}{{h}^{3}}+b{{h}^{2}}\right){\partial_x}\left(-\gamma{\partial_{xx} h}+\phi\right)+\sqrt{\frac{1}{3}{{h}^{3}}+b{{h}^{2}}}\xi_3 \right], \end{equation}

where the noise $\xi _3$ has zero mean and covariance $\langle \xi _3(x,t)\xi _3(x',t')\rangle =({2\mu k_B \theta }/{L_y})\delta (x-x')\delta (t-t')$. As the W-S model works for the no-slip case ($b=0$), the S-S model validates for the free-slip case ($b=\infty$) where the terms containing $b$ in (2.27) vanish. Interestingly, this free-slip model can be applied to study free films such as foam films (Erneux & Davis 1993; Vaynblat, Lister & Witelski 2001). Note that the free-slip model has not been reported previously.

3. Surface spectrum for films with strong slip

With the newly developed S-S model (2.27), the spectrum of surface waves in linear stages is derived here. We first rewrite the noise $\xi$ in terms of the white noise $N$ with unit variance $\langle {{N_i}}{{N_j }}' \rangle =\delta _{ij}\delta ( x-x' )\delta ( t-t' )$, and linearise (2.27) with $h=h_0+\tilde {h}$, $u=0+\tilde {u}$ and $N=0+\tilde {N}$:

(3.1)\begin{equation} \frac{{{\partial }^{2}}\tilde{h}}{\partial {{t}^{2}}}=\frac{\mu }{\rho }\frac{\partial }{\partial t}\left[ 4\frac{{{\partial }^{2}}\tilde{h}}{\partial {{x}^{2}}}-\frac{{\tilde{h}}}{{{h}_{0}}b} \right]-\frac{{{h}_{0}}}{\rho }\left[ \frac{\partial \phi }{\partial h}\frac{{{\partial }^{2}}\tilde{h}}{\partial {{x}^{2}}}+\gamma \frac{{{\partial }^{4}}\tilde{h}}{\partial {{x}^{4}}} \right]-f_1\frac{{{\partial }^{2}}{\tilde{N}_{1}}}{\partial {{x}^{2}}}+f_2\frac{\partial {\tilde{N}_{2}}}{\partial x}, \end{equation}

where the factors ${{f}_{1}}=\sqrt {8\mu {{k}_{B}}\theta {{h}_{0}L_y}}/{(\rho L_y) }$ and ${{f}_{2}}=\sqrt {{2\mu {{k}_{B}}\theta }{b}L_y}/{(\rho b L_y) }$. Note that the second derivative of $\tilde {h}$ with respect to $t$ comes from putting the linearised (2.27a) into the linearised (2.27b) to eliminate the variable $\tilde {u}$. Then a Fourier transform of (3.1) is performed using $\hat {h}(q,t)=\int _{-\infty }^{\infty }\tilde {h}(x,t){{\textrm {e}}^-}^{\textrm {i}qx}\,\textrm {d}\kern0.06em x$ and $\hat {N}(q,t)=\int _{-\infty }^{\infty }\tilde {N}(x,t){{\textrm {e}}^-}^{\textrm {i}qx}\,\textrm {d}\kern0.06em x$ to get

(3.2)\begin{equation} \frac{{{\partial }^{2}}\hat{h}}{\partial {{t}^{2}}}=-C\frac{\partial \hat{h}}{\partial t}-D\hat{h}+f_1{{q}^{2}}{{\hat{N}}_{1}}+f_2 q{\rm i}{{\hat{N}}_{2}}. \end{equation}

Here

(3.3a,b)\begin{equation} C=\frac{\mu }{\rho }\left( 4{{q}^{2}}+\frac{1}{{{h}_{0}}b} \right),\quad D=\frac{{{h}_{0}}}{\rho }\left[ \frac{\partial \phi }{\partial h}{{q}^{2}}+\gamma {{q}^{4}} \right]. \end{equation}

The solution of (3.2) can be represented as the linear superposition of two contributions (Zhang et al. 2019; Zhao et al. 2019)

(3.4)\begin{equation} \hat h={\hat h}_{{det}}+{\hat h}_{{sto}}, \end{equation}

where ${{\hat {h}}_{{det}}}$ is the solution to the deterministic part of (3.2) and ${{\hat {h}}_{{sto}}}$ is the contribution purely caused by thermal fluctuations. To find ${{\hat {h}}_{{det}}}$, the deterministic equation

(3.5)\begin{equation} \frac{{{\partial }^{2}}\hat{h}}{\partial {{t}^{2}}}+C\frac{\partial \hat{h}}{\partial t}+D\hat{h}=0, \end{equation}

is solved by Laplace transform. Using $g(q,s)=\int _{0}^{\infty }{\hat {h}(q,t){{\textrm {e}}^-}^{\textrm {i}ts}\,\textrm {d}}t$ and assuming ${\partial \hat {h}}/{\partial t}{{|}_{t=0}}$, one can get

(3.6)\begin{equation} g=\frac{s+C}{{{s}^{2}}+Cs+D}\hat{h}(q,0), \end{equation}

whose inverse Laplace transform is

(3.7)\begin{equation} {{\hat{h}}_{\det }(q,t)}=\hat{h}(q,0)\left[ \frac{{{{\rm e}}^{{{\omega }_{1}}t}}+{{{\rm e}}^{{{\omega }_{2}}t}}}{2}+\frac{C}{2\sqrt{{{C}^{2}}-4D}}\left( {{{\rm e}}^{{{\omega }_{1}}t}}-{{{\rm e}}^{{{\omega }_{2}}t}} \right) \right]. \end{equation}

Here $\omega _{i=1,2}=({-C\pm \sqrt {C^2-4D}})/{2}$ is the solution to $s^2+Cs+D=0$. Notably, $\omega _1$ is the dispersion relation (growth rate) of the deterministic lubrication equation, namely, (2.27) without the noise terms.

To obtain $\hat h_{sto}$, one has to determine the impulse response of the linear system. Using the Laplace transform of ${{{\partial }^{2}}\hat {h}}/{\partial {{t}^{2}}}+C({\partial \hat {h}}/{\partial t})+D\hat {h}=\delta$, and assuming $\hat {h}(q,0)=0$, it is found

(3.8)\begin{equation} g=\frac{1}{{{s}^{2}}+Cs+D}. \end{equation}

The impulse response is, thus, the inverse Laplace transform of (3.8)

(3.9)\begin{equation} F(q,t)={{\hat{h}}}=\frac{{{{\rm e}}^{{{\omega }_{1}}t}}-{{{\rm e}}^{{{\omega }_{2}}t}}}{\sqrt{{{C}^{2}}-4D}}. \end{equation}

Now with thermal fluctuations $f_1{{q}^{2}}\hat {N}_1$ and $f_2 qi\hat {N}_2$ as input, we find

(3.10)\begin{equation} {{\hat{h}}_{sto}}=f_1 q^2{{{}_{_{}}}\int_{0}^{t}{\hat{N}_1\left( q,t-\tau \right)}F(q,\tau )\,{\rm d}\tau } +f_2 q {\rm i}{{{}_{_{}}}\int_{0}^{t}{\hat{N}_2\left( q,t-\tau \right)}F(q,\tau )\,{\rm d}\tau }. \end{equation}

As $\hat {h}$ is both a random and complex variable, the root mean square (r.m.s.) of its norm is sought, namely, surface spectrum,

(3.11)\begin{equation} {{\left| \hat{h}(q,t) \right|}_{{rms}}}=\sqrt{\overline{{{\left| {{\hat{h}}_{{det}}}+{{\hat{h}}_{{sto}}} \right|}^{2}}}} =\sqrt{{{\left| {{\hat{h}}_{{det}}} \right|}^{2}}+\overline{{{\left| {{\hat{h}}_{{sto}}} \right|}^{2}}}}, \end{equation}

where from (3.7)

(3.12)\begin{equation} {\left| {{\hat{h}}_{{det}}} \right|}^{2}=\overline{{{\left| \hat{h}(q,0) \right|}^{2}}}{{\left[ \frac{{{{\rm e}}^{{{\omega }_{1}}t}}+{{{\rm e}}^{{{\omega }_{2}}t}}}{2}+\frac{C}{2\sqrt{{{C}^{2}}-4D}}\left( {{{\rm e}}^{{{\omega }_{1}}t}}-{{{\rm e}}^{{{\omega }_{2}}t}} \right) \right]}^{2}} , \end{equation}

and from (3.10)

(3.13)\begin{align} \overline{{{\left| {{\hat{h}}_{sto}} \right|}^{2}}}&={{{\left| \,{f_1 q^2}\int_{0}^{t}{\hat{N}_1\left( q,t-\tau \right)}F(q,\tau )\,{\rm d}\tau \right|}^{2}}} +{{{\left| \,{f_2 q {\rm i}}\int_{0}^{t}{\hat{N}_1\left( q,t-\tau \right)}F(q,\tau )\,{\rm d}\tau \right|}^{2}}} \nonumber\\ & ={{{\left| \,{f_1 q^2} \right|}^{2}}\int_{0}^{t}{{{\left| \hat{N}_1\left( q,t-\tau \right) \right|}^{2}}}F{{(q,\tau )}^{2}}\,{\rm d}\tau }+{{{\left| \,{f_2 q {\rm i}} \right|}^{2}}\int_{0}^{t}{{{\left| \hat{N}_1\left( q,t-\tau \right) \right|}^{2}}}F{{(q,\tau )}^{2}}\,{\rm d}\tau } \nonumber\\ & =\left( {{\left| \,{f_1 q^2} \right|}^{2}}+{{\left| \,{f_q q {\rm i}} \right|}^{2}} \right){{L}_{x}}\int_{0}^{t}{F{{(q,\tau )}^{2}}}\,{\rm d}\tau \nonumber\\ & =\frac{2\mu {{k}_{B}}T{{L}_{x}}{{q}^{2}}}{{{\rho }^{2}}{{L}_{y}}\left( {{C}^{2}}-4D \right)}\left( 4{{h}_{0}}{{q}^{2}}+\frac{1}{b} \right)\left[ \frac{{{{\rm e}}^{2{{\omega }_{1}}t}}-1}{2{{\omega }_{1}}}+\frac{{{{\rm e}}^{2{{\omega }_{2}}t}}-1}{2{{\omega }_{2}}}+\frac{2\left( {{{\rm e}}^{-Ct}}-1 \right)}{C} \right]. \end{align}

Here we have used $\overline {{{| \hat {N}( q,t ) |}^{2}}}={{L}_{x}}$, due to the finite length of the discrete Fourier transform used in MD simulations. Thus, we derive the spectrum of surface waves of a bounded film with strong slip as

(3.14)\begin{align} S\left(q,t\right) & =\left\{\overline{{{\left| \hat{h}(q,0) \right|}^{2}}}{{\left[ \frac{{{{\rm e}}^{{{\omega }_{1}}t}}+{{{\rm e}}^{{{\omega }_{2}}t}}}{2}+\frac{C}{2\sqrt{{{C}^{2}}-4D}}\left( {{{\rm e}}^{{{\omega }_{1}}t}}-{{{\rm e}}^{{{\omega }_{2}}t}} \right) \right]}^{2}} \right.\nonumber\\ & \quad \left.+\,\frac{2\mu {{k}_{B}}T{{L}_{x}}{{q}^{2}}}{{{\rho }^{2}}{{L}_{y}}\left( {{C}^{2}}-4D \right)}\left( 4{{h}_{0}}{{q}^{2}}+\frac{1}{b} \right)\left[ \frac{{{{\rm e}}^{2{{\omega }_{1}}t}}-1}{2{{\omega }_{1}}}+\frac{{{{\rm e}}^{2{{\omega }_{2}}t}}-1}{2{{\omega }_{2}}}+\frac{2\left( {{{\rm e}}^{-Ct}}-1 \right)}{C} \right]\right\}^{1/2}. \end{align}

Equations (3.14) and (2.27) are the main contributions of this work. For the W-S model, the analytical spectrum is (Zhang et al. 2020)

(3.15)\begin{equation} S\left(q,t\right)=\sqrt{\overline{{{\left| \hat{h}(q,0) \right|}^{2}}}{{{\rm e}}^{2{{\omega }_{3}}(q)t}}+\frac{{{L}_{x}}}{{{L}_{y}}}\frac{{{k}_{B}}T}{\gamma {{q}^{2}}+\left.\left( {\rm d}\phi /{\rm d}h \right)\right.|_{h_0}}\left[1-{{{\rm e}}^{2{{\omega }_{3}}(q)t}}\right]}, \end{equation}

where the dispersion relation is

(3.16)\begin{equation} {{\omega }_{3}}=-\frac{h_0^3+bh_0^2}{3\mu }\left( \gamma {{q}^{4}}+\left.\frac{{\rm d} \phi }{{\rm d} h}\right.| _{{{h}_{0}}}{{q}^{2}}\right). \end{equation}

4. MD simulations and a new slip-generating method

MD simulations are performed to simulate the instability of nanofilms on the strong-slip solid and the weak-slip solid. The open-source MD code LAMMPS (Plimpton 1995) is adopted. As shown in figure 2(a), the liquid film is composed of liquid argon (represented in orange) and it is simulated with the standard Lennard-Jones (LJ) 12-6 potential:

(4.1)\begin{equation} U({{r}_{ij}}) = \begin{cases} 4\varepsilon \left[ {{\left( \dfrac{\sigma }{{{r}_{ij}}} \right)}^{12}}-{{\left( \dfrac{\sigma }{{{r}_{ij}}} \right)}^{6}} \right], & \text{if}\ {{r}_{ij}}\le {{r}_{c}},\\ 0, & \text{if}\ {{r}_{ij}}>{{r}_{c}}. \end{cases} \end{equation}

Here $r_{ij}$ is the distance between two atoms and $ij$ represent pairwise particles. The energy parameter $\varepsilon$ is $1.67\times {10^{-21}}$ J and the length parameter $\sigma$ is $0.34\ \textrm {nm}$. We use $r_c=5.5 \sigma$ as the cutoff distance, beyond which the interaction vanishes.

Figure 2. Rupture of thin liquid films with strong slip in molecular simulations. (a) Initial setting of a liquid film with a flat free surface. The film has a small depth $L_y$ into the page. (b) Growth of perturbations and film rupture. (c) Droplet formation after the film ruptures. Note that the wall coloured in blue only denotes the position of the solid boundary as the method of a virtual wall is used in simulations.

The temperature of this system is kept at $T=85$ K using the Nosé–Hoover thermostat. At this temperature, the mass density is $\rho =1.4\times 10^3\ \textrm {kg}\ \textrm {m}^{-3}$ and the number density is $n=0.83/\sigma ^3$. The density of the vapour phase is about $(1/400) \rho$ so that the effects of vapour on the film are neglected. The surface tension of liquid is $\gamma = 1.52\times {10^{-2}}\ \textrm {N}\ \textrm {m}^{-1}$ and the dynamic viscosity is $\mu = 2.87\times {10^{-4}}\ \textrm {kg}\ (\textrm {ms})^{-1}$ (Zhang et al. 2020). The time step for all simulations is $0.004\sqrt {\varepsilon /(m\sigma ^2)}$ where $m$ is the atomic mass of argon.

The initial dimensions of the liquid film (see figure 2a) are chosen as $L_x=313.90\ \textrm {nm}$, $L_y=3.14\ \textrm {nm}$. The height of the film varies for three different cases ($h_0=1.2\ \textrm {nm}$, 1.6 nm and 3.14 nm). Therefore, the film is thin ($L_x\gg h_0$), and the film is quasi-2-D ($L_x\gg L_y$) to save computational costs. To prepare the initial configuration of the liquid film, the liquid film slab with the desired size is cut from a periodic box of liquid atoms, which makes the film surface flat initially.

Conventionally, the solid wall, serving as the boundary condition for the film, is simulated using real solid atoms (real wall) (Willis & Freund 2009; Zhang et al. 2019, 2020, 2021). In our previous work (Zhang et al. 2019, 2020, 2021), for a real wall, the solid is platinum with its isotropic $\langle 100\rangle$ surface in contact with the liquid. The liquid–solid interactions are modelled by the same 12-6 LJ potential with $\sigma _{ls}=0.8\sigma$ for the length parameter and $\varepsilon _{ls} =k\varepsilon$. Conventionally, one may vary the energy parameter $\varepsilon _{ls}$ to obtain different levels of slip length. For example, $b=0.68\ \textrm {nm}$ using $k=0.65$ whereas $b=8.8\ \textrm {nm}$ using $k=0.2$ (Zhang et al. 2021). The disjoining pressure and contact angles of liquid argon on the solid are also changed in this way. However, the slip length obtained using the real solid is usually in the range of a few tens of nanometres (Bocquet & Charlaix 2010) in MD simulations, which cannot match the micrometre slip length present in experiments. The origin of strong slip in experiments, as discussed earlier in the introduction, is usually complicated and cannot be straightforwardly simulated in previous MD simulations using LJ potentials. The use of a real solid wall with cutoff-distance-limited intermolecular interactions between liquid and solid also means that the disjoining pressure due to the solid vanishes when the film height is larger than the cutoff distance (Willis & Freund 2009; Zhang & Ding 2023), which is undesired.

Here we propose a new approach to allow us to generate any values of slip length simply in MD simulations and allow disjoining pressure to be effective at infinite distances physically. Instead of simulating a real wall, we apply a force to fluid atoms to mimic the fluid–substrate interactions and the force acts as the virtual wall, based on the work of (Steele 1973; Barrat & Bocquet 1999; Hadjiconstantinou 2021). The (total) force $\boldsymbol {f}$ for a fluid atom interacting with a face-centred-cubic solid substrate by the LJ potential has been calculated analytically by Steele (1973) and it is adopted here with some modifications (see Appendix A for detailed discussions):

(4.2)\begin{equation} \left. \begin{aligned} \boldsymbol{f} & ={{f}_{x}}\boldsymbol{e_x}+{{f}_{z}}\boldsymbol{e_z},\quad \text{with} \\ {{f}_{x}} & =\frac{{{(2{\rm \pi} )}^{2}}\varepsilon}{100{\ell}_{s}^{3}}\left[ \frac{ {{\sigma }^{12}}}{30}{{\left( \frac{{\rm \pi} }{{{\ell}_{s}}z} \right)}^{5}}{{{\rm K}}_{5}}\left( \frac{2{\rm \pi} }{{{\ell}_{s}}}z \right)-2{{\sigma }^{6}} {{\left( \frac{{\rm \pi} }{{{\ell}_{s}}z} \right)}^{2}}{{{\rm K}}_{2}}\left( \frac{2{\rm \pi} }{{{\ell}_{s}}}z \right) \right]\sin \left( \frac{2{\rm \pi} }{{{\ell}_{s}}}x \right),\\ {{f}_{z}} & =-\frac{{\rm d}U_0 \left(z\right)}{{\rm d}z}. \end{aligned} \right\} \end{equation}

Here $\boldsymbol {e}_x$ and $\boldsymbol {e}_z$ are unit vectors in the $x$ and $z$ direction, $\ell _s$ is the lattice spacing and $\textrm {K}_5$ and $\textrm {K}_2$ are the modified Bessel functions of the second kind. The applied force in the $x$ direction $f_x$ decays rapidly with $z$: $f_x$ is closely related to slip length and its sinusoidal-form force represents the energy corrugation of the solid surface. A smaller lattice spacing $\ell _s$ results in a smooth surface energy distribution and then a larger slip length (Thompson & Robbins 1990; Barrat & Bocquet 1999; Hadjiconstantinou 2021). In our MD simulations, the functions $\textrm {K}_5$ and $\textrm {K}_2$ are approximated by $(({{\rm \pi} }/{2x})^{1/2})\textrm {e}^{-x}$ for simplicity. Here $U_0 (z)$ is the total interaction energy exerted on a liquid molecule by the (continuous) substrate. As shown in Appendix A, $U_0(z)$ is related to the function of disjoining pressure as $U_0(z)=-\phi (z)/n$. Note that disjoining pressure is a function of film height $h$. One has to replace $h$ with $z$ while the form of the function itself is the same.

Here $\phi$ takes the usual form (Israelachvili 2011):

(4.3)\begin{equation} \phi(h)=\frac{A}{6{\rm \pi} h^3}-\frac{M}{h^9}, \end{equation}

where typical values $A=1.7\times 10^{-20}$ J and $M=0.018\varepsilon \sigma ^6$ are used. This form of disjoining pressure $h^{-3}-h^{-9}$ is obtained by integrating the 12-6 LJ potential over the entire substrate assuming the substrate is continuous, as shown in Appendix A and (Dietrich 1988; Schick 1990; Carey & Wemhoff 2005; Israelachvili 2011; MacDowell et al. 2014). There are many other forms of disjoining potential resulting from different liquid and solid properties (Becker et al. 2003).

We vary the lattice spacing $\ell _s$ and see what the slip length is from independent simulations where a pressure-driven flow goes past a substrate as shown by the MD snapshot in figure 3. The pressure gradient is created by applying a body force $g$ to the fluid. For $\ell _s=0.37\ \textrm {nm}$, one can see that the velocity profile in MD (blue squares in figure 3) is parabolic. However, for $\ell _s=0.10\ \textrm {nm}$, the velocity profile in MD (black triangles in figure 3) is nearly constant (plug flow).

Figure 3. Slip length measured using pressure-driven flows past the substrate in molecular simulations. MD results of velocity (symbols) are fitted with analytical solutions (solid lines) to obtain slip length.

The generated velocity distribution is

(4.4)\begin{equation} u(z)=\frac{\rho g}{2\mu }\left[(z-{{z}_{1}})(2{{z}_{2}}-z_1-z)+2b\left(z_2-z_1\right)\right]. \end{equation}

Here $z_1=0$ and $z_2=9.2\sigma$ are the positions of the wall and the free surface, respectively. This prediction (solid lines in figure 3) is used to fit the MD results (symbols in figure 3) to obtain the slip length. For $\ell _s=0.37\ \textrm {nm}$, a nearly no-slip surface $b=0.2\ \textrm {nm}$ is achieved. For $\ell _s=0.10\ \textrm {nm}$, a strong slip length of $b=400\ \textrm {nm}$ is obtained. Note that the choice of thermostats may have a slight influence on the slip length. For example, for water flows inside carbon nanotubes, the Berendsen and Nosé–Hoover thermostats result in very similar slip length, while a smaller slip length is found under the influence of the Langevin thermostat (Sam et al. 2017). Thus, the same thermostats should be chosen when measuring the slip length from pressure-driven flows and simulating nanofilm dewetting.

To mimic experimental conditions of polymer films on octadecyltrichlorosilane (OTS) and DTS substrates (Lessel et al. 2017), where contact angles (disjoining pressure) are about the same for both cases but slip length is very different (no slip on OTS in contrast to 500 nm slip length on DTS), we thus keep the disjoining pressure (4.3) the same for the simulations of weak-slip nanofilm dewetting and strong-slip nanofilm dewetting. In our simulations, the contact angles of drops after film dewetting for both cases are measured to be about $40^{\circ }$. This can also facilitate the comparison between weak-slip dewetting and strong-slip dewetting in simulations as the slip length is the only variable. Classically, disjoining pressure decreases with increasing slip length and contact angles. This is, however, far from being universal, as discussed previously for the case of polymer films on OTS and DTS substrates. Liquids on hydrophilic substrates with small contact angles can also have large slip length (Rothstein 2010; Ho et al. 2011). As shown in Appendix A, the proposed simulation technique is general and it can include the classic results where slip length and disjoining pressure are coupled and disjoining pressure decreases with increasing slip length and contact angles. As the virtual force is obtained by integrating the LJ potential (especially for $f_x$), this simulation technique is limited for the system where liquid and solid interact with LJ potential.

With the wall modelled by the virtual force, the simulations of nanofilm dewetting are run for 2 ns for the case of $h_0=1.2\ \textrm {nm}$, 10 ns for $h_0=1.6\ \textrm {nm}$ and 100 ns for $h_0=3.14\ \textrm {nm}$.

5. Results and discussion

5.1. Evolving spectra of an unstable film with strong slip

As shown in figure 2, a flat film deposited on the substrate experiences the spontaneous growth of perturbations on its surface, leading to film rupture (see figure 2b) and droplet formation (see figure 2c). To reveal the instability mechanism in the case of the strong-slip dewetting, the evolution of its surface spectra is obtained from MD simulations.

The instantaneous liquid–vapour interface $h(x,t)$, defined by the usual equimolar surface, is extracted from MD simulations (see (Zhang et al. 2019) for methods). A discrete Fourier transform of $h(x,t)$ is performed to obtain the amplitude of surface modes $\hat {h}(q,t)$. The surface spectra are thus calculated from the average (r.m.s.) of 40 independent simulations.

The symbols in figure 4(a) represent the MD spectra for the film with thickness $h=1.2\ \textrm {nm}$ and slip length $b=400\ \textrm {nm}$ at three different times. Compared with the weak-slip case ($b=0.2\ \textrm {nm}$) presented in figure 4(b), one can see that the transient characteristics of the spectra are strongly influenced by the slip length. The spectra for the strong-slip case grow faster than the spectra for the weak-slip case by comparing the blue triangles in figure 4(a) with the wine triangles in figure 4(b). The dominant wavenumber $q_d$ (the one with the maximum amplitude) at different times in the strong-slip dewetting are smaller than those in the weak-slip dewetting (e.g. at $t=0.43$ ns, $q_d\approx 0.25\ \textrm {nm}^{-1}$ for the strong-slip case in figure 4(a) in contrast to $q_d\approx 0.50\ \textrm {nm}^{-1}$ for the weak-slip case in figure 4b).

Figure 4. Evolution of the capillary spectra for the film with the thickness $h=1.2\ \textrm {nm}$ at three different times. (a) The film has a strong slip length $b=400\ \textrm {nm}$. Solid lines represent the prediction of the full S-S model whereas dashed lines are the S-S model without the inertial terms. (b) The film has a weak slip length $b=0.2\ \textrm {nm}$.

The proposed stochastic models (S-S model and W-S model) and their analytical spectra are used to predict the MD spectra. In our MD simulations, the initial setting of the film surface is flat $|\hat {h}(q,0)|\approx 0$, so that the deterministic contribution to the surface spectra is $|\hat {h}_{det}|\approx 0$. As elaborately investigated in our previous work (Zhang et al. 2019, 2020, 2021), the capillary spectra for weak-slip film can be predicted well by the W-S model, see e.g. the solid lines in figure 4(b). However, the W-S model does not apply to the case where the slip is strong, since the W-S model with $b=400\ \textrm {nm}$ leads to enormous overpredictions compared with MD results (not shown). This is because the derivation of the W-S model requires the slip length on the order of the film thickness. To predict the spectra of the unstable film with strong slip $b=400\ \textrm {nm}$, the analytical spectra (3.14) of the newly derived S-S model is adopted and it agrees excellently with MD results (see the solid lines in figure 4a). The symbols in figure 5(a) show the MD spectra for the film with a larger thickness $h=1.6\ \textrm {nm}$. Again, our analytical spectra can predict the MD results very well. Note that for large wavenumbers, the spectra at different times simply collapse into the static spectrum $S_s=\sqrt {L_x k_BT/(L_y\gamma q^2)}$.

Figure 5. Evolution of the capillary spectra for the film with the thickness $h=1.6\ \textrm {nm}$ at three different times. (a) The film has a strong slip length $b=400\ \textrm {nm}$. (b) The film has a weak slip length $b=0.2\ \textrm {nm}$.

Another interesting finding is that the assumption of Stokes flow, i.e. negligible inertia, breaks down when the value of slip length is increased from a weak slip to a strong slip. Ignoring the inertial terms in the derived S-S model (the left-hand-side terms of (2.27)), the surface spectra are simply

(5.1)\begin{equation} S\left( q,t \right)=\sqrt{\frac{{{L}_{x}}}{{{L}_{y}}}\frac{{{k}_{B}}T}{\gamma {{q}^{2}}+{\rm d}\phi /{\rm d}h}\left[ 1-{{{\rm e}}^{-2Dt/C}} \right]}, \end{equation}

which however overpredicts the MD results (see the dashed line in figures 4a and 5a). Note that the W-S model where inertia is also absent works well for predicting the MD results of weak-slip dewetting. Therefore, it is the strong slip that brings the inertia into effect. In the weak-slip regime, the Reynolds number is

(5.2)\begin{equation} {Re}_{W-S}=\frac{\rho u_0 h_0}{\mu}=\chi^3 \frac{\rho \gamma h_0}{\mu^2}=\chi^3 La. \end{equation}

Here we have used the momentum balance in the horizontal direction to estimate the characteristic velocity $u_0\sim \chi ^3 \gamma /\mu$ as discussed previously. However, in the strong-slip regime, the ${Re}_{S\textrm -S}$ is

(5.3)\begin{equation} {Re}_{S\text-S}=\frac{\rho u_0 h_0}{\mu}=\chi \frac{\rho \gamma h_0}{\mu^2}=\chi La, \end{equation}

based on the momentum balance in the vertical direction. Therefore, ${Re}_{S\textrm -S}$ of the strong-slip case is $1/\chi ^2$ larger than the weak-slip case. On the other hand, the Laplace number ${La}$ in our simulations is calculated to be $0.3$ which is consistent with the assumption of deriving the strong-slip model where the ${La}\sim 1$. Both conditions (strong slip and ${La}=0.3$) make inertia non-negligible in our simulations.

5.2. Number of droplets

As shown in figure 6, one distinct feature of the strong-slip dewetting (see figure 6b,d) is having fewer droplets after the film ruptures compared with the weak-slip dewetting (see figure 6a,c). For each case with a specific film thickness and slip length (for example, in figure 6(a), $h=1.2\ \textrm {nm}$ and $b=0.2\ \textrm {nm}$), about 40 independent simulations have been run and the probability distribution of the number of formed droplets is calculated and shown in figure 7(a) for $h=1.2\ \textrm {nm}$ and figure 7(b) for $h=1.6\ \textrm {nm}$.

Figure 6. Droplets form after the rupture of the film: (a) for the film with $h=1.2\ \textrm {nm}$ and $b=0.2\ \textrm {nm}$; (b) for the film with $h=1.2\ \textrm {nm}$ and $b=400\ \textrm {nm}$; (c) for the film with $h=1.6\ \textrm {nm}$ and $b=0.2\ \textrm {nm}$; (d) for the film with $h=1.6\ \textrm {nm}$ and $b=400\ \textrm {nm}$.

Figure 7. (a) Distribution of the number of droplets for the film with $h=1.2\ \textrm {nm}$. (b) Distribution of the number of droplets for the film with $h=1.6\ \textrm {nm}$. (c) The decrease of the dominant wavenumber with time predicted by (3.14) (solid lines) to its asymptotic value by the deterministic lubrication equation (dashed lines in the inset). The squares represent the dominant wavenumber at the time of film rupture.

In terms of the thinner film $h=1.2\ \textrm {nm}$, the average number of droplets for weak-slip dewetting is $N_a=14$ (see the blue bars in figure 7(a) and a typical MD snapshot in figure 6a), whereas the number is $N_b=6.5$ for the strong-slip dewetting (see the red bars in figure 7(a) and a typical MD snapshot in figure 6b). For the thicker film with $h=1.6\ \textrm {nm}$, the mean number of droplets for weak-slip dewetting is $N_c=7$ (see the blue bars in figure 7(b) and the MD snapshot in figure 6c) whereas the number is $N_d=3$ for strong-slip dewetting (see the red bars in figure 7(b) and the MD snapshot in figure 6d).

The number of droplets is connected with the dominant wavenumber $q_d$ (or the dominant wavelength $\lambda _d=2{\rm \pi} /q_d$). In the deterministic situation, the dominant wavenumber is constant over time. For the deterministic W-S equation, it is $q_{d}=\sqrt {-({\textrm {d}\phi }/{\textrm {d}h})/({2\gamma })}\approx \sqrt {{A}/({4{\rm \pi} \gamma h_0^4})}$, while for the deterministic S-S equation, it is the solution when $\textrm {d}\omega _1(q)/\textrm {d}t=0$.

However, as found from MD simulations shown in figures 4 and 5, the dominant wavenumber decreases with time, which is surely beyond the prediction from the deterministic lubrication equation. Based on the analytical spectra of the W-S model (3.15) and S-S model (3.14), figure 7(c) shows the theoretical prediction of the evolving dominant wavenumber. One can find that thermal fluctuations induce a dominant wavenumber (see the solid lines) much higher than its corresponding deterministic prediction (see the dashed lines in the inset) and it only gradually decreases to the deterministic predictions over time. In our simulations, the film can rupture far before the dominant wavenumber reaches its asymptotic value.

Therefore, we measure the average time of rupture $t_r$ from MD simulations (shown by the symbols in figure 7c) and use the analytical theory to predict the dominant wavenumber at the time of rupture. For $h=1.2\ \textrm {nm}$ and $b=0.2\ \textrm {nm}$, the rupture time is about $t_r=2$ ns and the dominant wavenumber at this time is $q_d=0.35\ \textrm {nm}^{-1}$. The dominant wavenumber corresponds to $q_d L_x/(2{\rm \pi} )\approx 17$ droplets, in agreement with the $14$ measured directly from MD simulations. For $h=1.2\ \textrm {nm}$ and $b=400\ \textrm {nm}$, the rupture is much faster $t_r=0.75$ ns and the dominant wavenumber at this time is $q_d=0.156\ \textrm {nm}^{-1}$. This dominant wavenumber translates into $7.8$ droplets, in agreement with the $6.5$ measured directly from MD simulations. Therefore increasing slip length from weak slip to strong slip can indeed reduce the number of droplets by decreasing the dominant wavenumber.

Since the strong-slip dewetting and the weak-slip dewetting starts with the same surface condition (flat surface), the deviations for both cases in the dominant wavenumber begin at a very early time (see figure 7c). The smaller dominant wavenumber in the strong-slip dewetting, compared with the weak-slip dewetting, becomes the obvious feature of the strong-slip spectra shown in figure 4(a).

By increasing the film thickness from $h=1.2\ \textrm {nm}$ to $h=1.6\ \textrm {nm}$, the number of droplets is also decreased. For example, for the strong-slip dewetting, the average number of droplets changes from 6.5 to 3 by with increased film thickness as shown in figure 7(b). This can be also predicted by the dominant wavenumbers shown in figure 7(c) where the dominant wavenumber is smaller for a larger film thickness.

5.3. Surface roughening with strong slip

The growth of spectra is equivalent to the surface roughening of the free surface based on Parseval's theorem:

(5.4)\begin{equation} W(t)=\sqrt{\frac{1}{{{L}_{x}}}\left\langle\int_{0}^{{{L}_{x}}}{{{\left( \delta h \right)}^{2}}\,{{\rm d}\kern0.06em x}}\right\rangle}=\sqrt{\frac{1}{2{\rm \pi} {{L}_{x}}}\int_{{2{\rm \pi}}/{{{L}_{x}}}}^{{2{\rm \pi} K}/{{{L}_{x}}}}{S^2\,{\rm d}q}} , \end{equation}

where $W(t)=\sqrt {\overline {h^2}-{\bar {h}}^2}$ is the r.m.s. roughness of the free surface. Here $K$ is the number of bins used to extract the surface profile from MD simulations, which provides an upper bound on the wavenumbers that can be extracted (Zhang et al. 2021).

In our study, the roughening of a flat liquid surface is due to thermal fluctuations (and also disjoining pressure). Surface roughening or growth resulting from other forms of randomness is a common occurrence in nature. Examples include the propagation of wetting fronts in porous media, the growth of bacterial colonies, and the atomic deposition process in the manufacture of computer chips (Barabási & Stanley 1995). These evolutions of the profile of a growing interface can be described by stochastic partial differential equations (SPDE). One of the most famous examples is the Kardar–Parisi–Zhang (KPZ) equation (Kardar, Parisi & Zhang 1986). Scaling relations for surface roughening are then used to analyse the SPDE, which can be summarised as (Barabási & Stanley 1995)

(5.5)\begin{equation} W\left(t\right)\sim L^{\alpha} f(t/L^m), \end{equation}

where $L$ is the system size, $f(v)=v^{\kappa }$ for $v\ll 1$ (during roughness growth), and $f(v)=1$ for $v\gg 1$ (at roughness saturation). The time to transition, between roughness growth and saturation, scales with $t_s \sim L^m$. The three exponents ($\alpha$, $m$ and $\kappa$) define a universality class, and are here related by $\kappa = \alpha /m$. For example, for the one-dimensional KPZ equation, $(\alpha =1/2,m=1/3,\kappa =3/2)$. Basically, the roughness will grow as a power law of time $W\sim t^{\kappa }$ before getting saturated, which is a result of the balance of the deterministic forces (such as surface tension in our study) and stochastic forces (such as thermal motions of molecules in our study).

In our previous work (Zhang et al. 2021), we have shown that for a weak-slip film roughened by thermal fluctuations (with negligible effects of disjoining pressure), a universality class $(\alpha =1/2,m=4,\kappa =1/8)$ exists and $W\sim t^{1/8}$ before the roughness saturation. It is interesting to see how strong slip changes the scaling of surface roughening. We note that the existence of the disjoining pressure breaks the potential balance of surface tension and thermal fluctuations, resulting in the unbounded growth of the surface roughness. However, the initial growth of surface roughness, when the disjoining pressure is weak, may still be described by scaling laws.

As shown in figure 8, we calculate the surface roughening of a film with the thickness $h=3.14\ \textrm {nm}$ for different levels of slip length. In terms of the weak-slip film $b=0.2\ \textrm {nm}$ and without disjoining pressure, the roughness grows with $W\sim t^{1/8}$ and then becomes saturated eventually (see the blue dash line in figure 8). When disjoining pressure is considered, the growth of the roughness becomes unbounded (see the red solid line). However, the initial stage of the growth of roughness ($t<100$ ns) is unaffected by disjoining pressure so that the $W\sim t^{1/8}$ is still valid.

Figure 8. The growth of roughness on films with weak slip and strong slip.

For the strong-slip film $b=400\ \textrm {nm}$ and without disjoining pressure, we find $W\sim t^{1/4}$ as shown by the blue dashed line in figure 8, which can be derived as follows. Here $\alpha$ can be obtained by considering the surface at saturation, i.e. from the static spectrum given by $S_s=\sqrt {L_x k_BT/(L_y\gamma q^2)}$. Substituting the static spectrum into (5.4) leads to

(5.6)\begin{equation} W_s=\sqrt{\frac{1}{2{\rm \pi} L_y}\frac{k_BT}{\gamma}\left(\frac{L_x}{2{\rm \pi}}-\frac{L_x}{2{\rm \pi} K}\right)}. \end{equation}

For large $K$, which is the case in our MD simulations, $W_s$ becomes independent of $K$ and one can find that $W_s=\sqrt {({L_x}/{4{\rm \pi} ^2 L_y})({k_BT}/{\gamma })}\sim {L_x}^{{1}/{2}}$.

An upper estimate on the transition time $t_s\sim L^m$, between growth and saturation, can be estimated from the inverse of the dispersion relation at the largest permissible wavelength ($q={2{\rm \pi} }/{L_x}$). Examining the dispersion relation indicates that there are three time scales for the spectra to reach thermal equilibrium (without disjoining pressure) whose maximum is the right time scale, namely,

(5.7)\begin{align} t_s&=\frac{1}{\min \left(|2\omega_1\left(q=2{\rm \pi}/L_x\right)|, |2\omega_2\left(q=2{\rm \pi}/L_x\right)|, C \left(q=2{\rm \pi}/L_x\right)\right)}\nonumber\\ &=\frac{1}{|2\omega_1\left(q=2{\rm \pi}/L_x\right)|}\nonumber\\ &\sim{L_x^2}, \end{align}

where we have used the condition $C=({\mu }/{\rho })( 4{{q}^{2}}+{1}/{{{h}_{0}}b})\approx {4{{q}^{2}}\mu }/{\rho }$ due to the fact that $4q^2h_0b=4/\chi \gg 1$ (using the long-wave condition $qh_0 \sim \chi$ and the strong-slip condition $bh_0 \sim \chi ^{-2}$).

Therefore, we find the power $m=2$ and $\kappa =\alpha /m=1/4$ for the strong-slip case. In summary, the universality class $(\alpha =1/2,m=2,\kappa =1/4)$ is found for the surface growth of a thin film with strong slip and negligible disjoining pressure. When disjoining pressure is considered, the initial growth of the surface may still be described by $W\sim t^{{1}/{4}}$ (as shown by the blue solid line in figure 8), which is much faster than the $W\sim t^{{1}/{8}}$ of the weak-slip case. These scalings are confirmed in MD simulations (see the symbols in figure 8).

6. Conclusions

In this work, the instability of molecularly thin liquid films on substrates with strong slip has been investigated using both MD simulations and a newly developed stochastic lubrication equation. A new method has been proposed to generate strong slip length in molecular simulations. Our simulations reveal that the strong-slip dewetting is much faster than the weak-slip dewetting and has fewer droplets formed compared with the weak-slip dewetting. Using a long-wave approximation to the equations of FH, a new stochastic lubrication equation considering strong slip and inertia effects has been derived. By the linear stability analysis, the analytical surface spectrum has been obtained and can explain the differences between strong-slip dewetting and weak-slip dewetting observed in simulations. We have further found that inertial effects, which are usually ignored in weak-slip dewetting, come into play in the strong-slip dewetting due to the enhanced velocity by the strong slip.

When the effect of disjoining pressure is negligible, the derived strong-slip stochastic lubrication equation possesses a universality class $(1/2,2,1/4)$ in contrast to the $(1/2,4,1/8)$ of the weak-slip stochastic lubrication equation. In our simulations, the effect of disjoining pressure is important. However, the surface roughening of strong-slip dewetting at the initial stage may be still described by $W\sim t^{1/4}$ in contrast to the $W\sim t^{1/8}$ for weak-slip dewetting.

Experiments of the rupture of polymer nanofilms on no-slip SiO$_{2}$-coated silicon wafers have demonstrated the great effects of thermal fluctuations (Fetzer et al. 2007). The experimental surface spectra can be predicted with the no-slip film stochastic lubrication equation (Grün et al. 2006; Fetzer et al. 2007). However, polymer films on DTS solid can have slip length up to several micrometres (Bäumchen & Jacobs 2009). Therefore, experiments on the rupture of polymer nanofilms on DTS solid can be used to validate the newly developed theory here. However, care should be taken to make sure the dewetting of polymer films is initiated by disjoining pressure since polymer films can be easily contaminated (leading to nucleation dewetting) (Jacobs, Herminghaus & Mecke 1998). As polymers have a very large viscosity $\sim 10^4\ \textrm {kg}\ (\textrm {ms})^{-1}$, the Laplace number is very small so that the effects of inertia can be safely neglected. For metallic films, however, its viscosity is about $\sim 10^{-3}\ \textrm {kg}\ (\textrm {ms})^{-1}$ so that inertial effects can be important (Diez et al. 2016; González et al. 2016).

In the limit of infinite slip, the derived strong-slip stochastic lubrication equation leads to a new model for free films, which can be readily used to study, for example, the symmetrical breakup of a foam film under the effects of thermal fluctuations (Erneux & Davis 1993; Vaynblat et al. 2001). In this work, we focus on the linear stability theory and molecular simulations of nanofilm dewetting, numerical solutions to these stochastic lubrication equations remain to be investigated in the future (see, e.g. Grün et al. 2006; Nesic et al. 2015; Diez et al. 2016; Durán-Olivencia et al. 2019; Shah et al. 2019; Zhao, Lockerby & Sprittles 2020). Future directions may also include the effects of contamination (such as surfactants) (Zhang & Ding 2023) and evaporation (Oron et al. 1997; Craster & Matar 2009) on film stability which often occur at small scales.