Stokes theory of thin-film rupture

TL;DR

This paper investigates the flow dynamics leading to thin-film rupture due to van der Waals forces, revealing a universal self-similar solution that differs from lubrication theory predictions.

Contribution

It demonstrates that near rupture, the flow transitions from lubrication theory to a universal Stokes flow solution, providing new insights into thin-film rupture dynamics.

Findings

Lubrication theory predicts $h_{min} \,\propto\, \tau^{1/5}$.

Near rupture, flow follows a universal self-similar solution with $h_{min} \,\propto\, \tau^{1/3}$.

The opening angle near rupture is approximately 37 degrees.

Abstract

The structure of the flow induced by the van der Waals destabilization of a non-wetting liquid film placed on a solid substrate is unraveled by means of theory and numerical simulations of the Stokes equations. Our analysis reveals that lubrication theory, which yields $h_{\text{min}} \propto \tau^{1/5}$ where $h_{\text{min}}$ is the minimum film thickness and $\tau$ is the time until breakup, cannot be used to describe the local flow close to rupture. Instead, the slender lubrication solution is shown to experience a crossover to a universal self-similar solution of the Stokes equations that yields $h_{\text{min}} \propto \tau^{1/3}$, with an opening angle of $37^{\circ}$ off the solid.

Equations

$$\nabla \cdot \mathbf{u}=0,0=-\nabla \varphi +\nabla \cdot \mathbf{T},$$

$$\mathbf{T}\cdot \mathbf{n}+\left(\nabla \cdot \mathbf{n}\right)\mathbf{n}$$

$$\mathbf{n}\cdot \left(\partial {\mathbf{x}}_s/\partial t-\mathbf{u}\right)$$

$$h_t+{\left(h^3{h}_{xxx}/3+h^{-1}{h}_x\right)}_x=0,$$

$$\left[u^{\ast },v^{\ast },p^{\ast }\right]$$

$$h^{\ast }$$

$$\left[u^{\ast },v^{\ast }\right]$$

$$p^{\ast }$$

$$h^{\ast }$$

$$\left[u^{\ast },v^{\ast }\right]$$

$$p^{\ast }$$

$$h^{\ast }$$

$${\omega }_{\text{ Stokes}}$$

$${\omega }_{\text{ Lub}}$$

$${\omega }_{\text{ Stokes}}=k^2\left[\frac{1}{(h_o/a)}-\frac{(h_o/a{)}^3}{3}{k}^2\right]-\frac{9(h_o/a)}{5}{k}^4+\left[\frac{23(h_o/a{)}^3}{7}+\frac{3(h_o/a{)}^5}{5}\right]k^6+O\left(k^8\right),$$

$$t_{\mathrm{R}}^{\text{Lub}}$$

$$\frac{t_{\mathrm{R}}^{\text{Stokes}}}{t_{\mathrm{R}}^{\text{Lub}}}$$

$$x={\tau }^{1/3}\xi ,y={\tau }^{1/3}\eta ,h={\tau }^{1/3}f(\xi ),u={\tau }^{-2/3}U(\xi ,\eta ),v={\tau }^{-2/3}V(\xi ,\eta ),p={\tau }^{-1}P(\xi ,\eta ),$$

$$U_{\xi }+V_{\eta }=0,$$

$$U_{\xi \xi }+U_{\eta \eta }=P_{\xi }-3f^{-4}{f}_{\xi },$$

$$V_{\xi \xi }+V_{\eta \eta }=P_{\eta },$$

$$(1+f_{\xi }^2)P-2V_{\eta }+2(V_{\xi }-f_{\xi }{U}_{\xi }+U_{\eta })f_{\xi }=0,$$

$$(1-f_{\xi }^2)(V_{\xi }+U_{\eta })+2(V_{\eta }-U_{\xi })f_{\xi }=0,$$

$$f/3-(\xi /3+U)f_{\xi }+V=0,$$

$$F^{\text{(iv)}}+\left[4+2\lambda (2+\lambda )\right]F^{\prime \prime }+{\lambda }^2(2+\lambda {)}^{2}F=0,$$

$$(\lambda +1{)}^2\cos (2{\theta }_o)+\cos [2{\theta }_o(\lambda +1)]-\lambda (\lambda +2)=0,$$

$$Re=\frac{\rho h_{\text{min}}^{\ast }}{\mu }\frac{\mathrm{d}{h}_{\text{min}}^{\ast }}{\mathrm{d}{\tau }^{\ast }}\Rightarrow Re=\frac{\rho a\sigma }{{\mu }^2}{h}_{\text{min}}\frac{\mathrm{d}{h}_{\text{min}}}{\mathrm{d}\tau },$$

$$Re=0.034\frac{\rho \sqrt{A\sigma }}{{\mu }^2}{\tau }^{-1/3}\to \infty \text{as}\tau \to 0.$$

$${\tau }_I=3.9\times 10^{-5}\frac{{\rho }^3{\left(A\sigma \right)}^{3/2}}{{\mu }^6}.$$

$$\varphi (h)=h^{-3}-\frac{48\pi ϵ}{A{a}^6}{h}^{-9},$$

$${\int }_{\Omega }\left[\tilde{p}\nabla \cdot \mathbf{u}-(p+\varphi )\nabla \cdot \tilde{\mathbf{u}}+\left(\nabla \mathbf{u}+\nabla {\mathbf{u}}^{\mathrm{T}}\right):\nabla \tilde{\mathbf{u}}\right]\mathrm{d}\Omega +{\int }_{{\Gamma }_s}{\nabla }_s\cdot \tilde{\mathbf{u}}\mathrm{d}S=0,$$

$$h_t+{\left(h^3{h}_{xxx}/3+h^{-1}{h}_x\right)}_x=0,t>0,0<x<\pi /k,$$

$$(h_t,0)+{\partial }_x\mathbf{G}=\mathbf{g},\text{with}\mathbf{G}=\left[h_x/h+h^3{h}_{xxx}/3,h_x\right],\mathbf{g}=(0,h_{xxx}),$$

$${\int }_{\Omega }\left[\tilde{P}\nabla \cdot \mathbf{U}-\left[P+f^{-3}\right]\nabla \cdot \tilde{\mathbf{U}}+\left(\nabla \mathbf{U}+\nabla {\mathbf{U}}^{\mathrm{T}}\right):\nabla \tilde{\mathbf{U}}\right]\mathrm{d}\Omega =0,$$

$$\left(-P\mathbf{I}+\nabla \mathbf{U}+\nabla {\mathbf{U}}^{\mathrm{T}}\right)\cdot \mathbf{n}=0,$$

$${\int }_{{\Gamma }_s}\tilde{\delta }\left[f/3-(\xi /3+U)f_{\xi }+V\right]\mathrm{d}\Gamma =0,$$

Full text

††preprint: Accepted for publication in Phys. Rev. Fluids

I Introduction

Thin liquid films are ubiquitous in nature and everyday life, and they play important roles in many engineering processes, medical and physiological contexts, and in geophysics and biophysics, among other fields. As a consequence, many research efforts have been devoted to study the structure, stability and nonlinear dynamics of liquid films. For instance, falling films are relevant in coating processes Huppert (1982a); Kalliadasis et al. (2011), gravity currents are important in geological phenomena  Pattle (1959); Huppert (1982b, 1986), and pulsed-laser-heated thin films are used in patterning and plasmonic applications Makarov et al. (2016); Hughes et al. (2017); Kondic et al. (2019). Tear-film dynamics Braun (2012); Fuller and Vermant (2012); Hermans et al. (2015); Dey et al. (2019) and surfactant replacement therapy Jensen and Grotberg (1992); Halpern and Grotberg (1992); Jensen and Grotberg (1993); Halpern and Grotberg (1993); Grotberg (1994); Halpern et al. (1998); Cassidy et al. (1999) are examples where the stability of surfactant-driven films is crucial for the healthy functioning of the eye and lungs. The reader is referred to De Gennes (1985); Oron et al. (1997); Bonn et al. (2009); Craster and Matar (2009); Blossey (2012); Kondic et al. (2019) for detailed and excellent reviews of the vast amount of scientific literature devoted to liquid films.

The present work is of fundamental character, and deals with a non-wetting ultrathin liquid film placed on a solid substrate. The liquid film becomes unstable to infinitesimal surface waves when its thickness becomes smaller than about 100 nm, leading to a spinodal dewetting pathway coexistent with hole nucleation Wyart and Daillant (1990); Reiter (1992); Oron et al. (1997); Reiter et al. (1999); Seemann et al. (2001); Becker et al. (2003); Thiele (2007); Craster and Matar (2009); Bonn et al. (2009); Blossey (2012). The spontaneous growth of perturbations takes place when the destabilizing van der Waals (vdW) forces exceed the stabilizing surface tension force, provided that the disturbance wavenumber is below a certain cut-off (Scheludko, 1962; Vrij, 1966; Ruckenstein and Jain, 1974). Previous theoretical efforts to describe the nonlinear dynamics leading to film rupture were based on lubrication theory (Williams and Davis, 1982; Burelbach et al., 1988; Zhang and Lister, 1999; Blossey, 2012), which assumes that the longitudinal length scale is much larger than the film thickness, and provides models with simpler mathematical structure than the Navier-Stokes equations. Indeed, while the latter must be solved as a free boundary problem where the film thickness $h$ is part of the solution, the former leads to a partial differential equation (PDE) for $h$ as a function of the relevant spatial coordinates and time, that is physically transparent and more amenable to analysis.

In this paper, through dimensional arguments, numerical computations and similarity theory, we reveal that the Stokes flow close to the rupture singularity is not slender, and provides results markedly different from previous ones based on lubrication theory Zhang and Lister (1999). Our local theory, motivated by dimensional analysis, is inspired by pioneering studies of singularities in free-surface flows without resorting to one-dimensional approximations of the equations of motion, in contexts like the breakup of liquid jets (Papageorgiou, 1995; Day et al., 1998), or the ejection of jets from Faraday waves (Zeff et al., 2000).

II Formulation

Consider a planar film of Newtonian liquid of viscosity $\mu$ coating a solid surface that spans the $(x,z)$ plane. As sketched in Fig. 1 the film, of initial height $h_{o}$, is surrounded by a passive gaseous atmosphere at constant pressure, such that the gas-liquid interface has a surface tension $\sigma$, and is described by the function $y=h(x,t)$. The liquid film, initially at rest, becomes unstable due to the long-range vdW forces, whose collective effects are modeled through a disjoining pressure, $A/(6\pi h^{3})$, with associated Hamaker constant $A$ Hamaker (1937), being $A>0$ in the non-wetting case considered herein. The latter intermolecular force model, which considers only non-retarded vdW interactions, is the simplest one among a hierarchy of existing models to rationalize the experimental observations Thiele et al. (2001); Seemann et al. (2001); Israelachvili (2011). Note that, as argued in Appendix A, liquid inertia is negligible under realistic experimental conditions. The cartesian components of the liquid velocity field in the $(x,y)$ directions are $(u,v)$, and the pressure field is $p$.

The molecular length scale De Gennes (1985), $a=[A/(6\pi\sigma)]^{1/2}$, is taken as the relevant characteristic length scale for the ultrathin liquid films considered in the present work, together with $\mu a/\sigma$, $\sigma/\mu$ and $A/(6\pi a^{3})$ as time, velocity and pressure scales, respectively. The non-dimensional Stokes equations read

[TABLE]

where $\bm{u}=(u,v)$, $\bm{T}=-p\bm{I}+\bm{\nabla}\bm{u}+(\bm{\nabla}\bm{u})^{\rm{T}}$ is the liquid stress tensor and $\phi=h^{-3}$ is the dimensionless vdW potential. The accompanying boundary conditions include the non-slip condition $\bm{u}=0$ at the solid wall $y=0$, and

[TABLE]

at the free surface $y=h(x,t)$, with corresponding parametrization $\bm{x}_{s}$ and unit normal vector $\bm{n}$, accounting for the stress balance and the kinematics of the interface, respectively.

In contrast with Stokes flow, the lubrication approximation provides the much simpler leading-order description

[TABLE]

governing the evolution of the free surface under the small-slope assumption Williams and Davis (1982); Burelbach et al. (1988); Zhang and Lister (1999). Hereafter, subscripts will denote partial derivatives. Finally, note that the dimensionless initial film thickness, $h_{o}/a$, is the only governing parameter.

III Results

The presentation of results is divided into three parts. First, we will show that dimensional analysis demonstrates the existence of a non-slender self-similar solution of the first kind for the Stokes equations, governed by a balance between viscous and vdW forces, for which capillary forces are negligible. Second, we will present several comparisons of the film evolution given by the numerical integration of the Stokes equations, with that predicted by the lubrication model. Finally, we formulate and solve the self-similar Stokes problem, comparing the resulting solution with the rescaled numerical results for the full temporal evolution for times close to rupture.

III.1 Dimensional analysis

Dimensional arguments suggest the existence of a similarity solution of the Stokes equations near film rupture that differs from the lubrication result Zhang and Lister (1999). Indeed, the parametric dependences of the longitudinal velocity, transverse velocity, pressure and film thickness are

[TABLE]

where $\tau_{*}=t^{*}_{\mathrm{R}}-t^{*}$ is the time remaining to rupture, and asterisks denote the dimensional versions of the flow variables. Taking $(\tau_{*},\mu,A)$ as dimensional basis, the Buckingham $\Pi$ theorem provides the reduced functional dependences,

[TABLE]

where $[\xi,\eta,\Pi_{h_{o}}]=(\mu/A)^{1/3}\tau_{*}^{-1/3}[x^{*},y^{*},h_{o}]$ and $\Pi_{\sigma}=\sigma/(A^{1/3}\mu^{2/3})\tau_{*}^{2/3}$. When $\tau_{*}\to 0$, $\Pi_{\sigma}\to 0$ and $\Pi_{h_{o}}\to\infty$, suggesting that as rupture is approached surface tension forces become negligible, and that the local flow becomes independent of $h_{o}$. Note also that the velocity components blow up near rupture, as dictated by (7), a common feature exhibited in most finite-time singularities in free surface flows. Thus, we expect a local self-similar Stokes flow of the form

[TABLE]

Note also that, since the self-similar scalings for $x$ and $y$ are the same, it can be anticipated that the approach of the flow to the singularity cannot be described using lubrication theory, which assumes that the characteristic length in the $x$ direction is much larger than the film thickness.

III.2 Flow evolution

The equations (1)–(3) were numerically integrated with the symmetry conditions $u=v_{x}=0$ imposed at the planes $x=0$ and $x=\pi/k$, and the initial conditions $\bm{u}=0$ and $h(x,0)=h_{o}/a(1-10^{-3}\cos{kx})$ for $0<x<\pi/k$. Here, the wavenumber $k$ was chosen as the wavenumber of maximum amplification deduced from a linear stability analysis of Eqs. (1)–(3) for the Stokes description Jain and Ruckenstein (1976); González et al. (2016), and of Eq. (4) for the lubrication approximation Vrij (1966), which yield, respectively, the dispersion relations

[TABLE]

The Stokes dispersion relation (13) has the small-$k$ expansion

[TABLE]

which, at $O(k^{4})$, differs from the lubrication dispersion relation by the factor $9(h_{o}/a)k^{4}/5$. Note that the latter correction of the lubrication result is important for $h_{o}/a\sim 1$, but becomes negligible for $h_{o}/a\gg 1$. Using the lubrication equation (14), and the expansion (15) of the Stokes dispersion relation truncated at $O(k^{6})$, the film rupture time can be estimated, using the corresponding maximum growth rates, as

[TABLE]

respectively, where $\varepsilon$ is the initial disturbance amplitude. The approximation (16) for $t_{\mathrm{R}}$ is plotted as a thin line in Fig. 2($d$) for $\varepsilon=10^{-3}$, showing good agreement with the numerical rupture time. From Eq. (17) it is also deduced that the Stokes equations predict a slightly longer lifetime for the liquid film than the lubrication approximation. A detailed description of the numerical techniques can be found in Appendix C.

A representative numerical integration is presented in Fig. 1 for $h_{o}/a=4.34$. The slightly disturbed flat film profile departs from the initial condition by virtue of the destabilizing vdW forces in a self-accelerated process, leading to a rupture singularity in a finite time $t_{\mathrm{R}}$, whose precise computation involved an algebraic fitting procedure that took advantage of the anticipated power-law behavior $h_{\textrm{min}}\propto\tau^{1/3}$ for $\tau\to 0$. Sample computations are depicted in Fig. 2 for several values of $h_{o}/a$. The accompanying instantaneous exponent $n=\mathrm{d}\log_{10}{h_{\textrm{min}}}/\mathrm{d}\log_{10}{\tau}$ reveals the persistent self-similar behavior $h_{\textrm{min}}\to K_{S}\,\tau^{1/3}$ for $\tau\lesssim 0.1$ and all values of $h_{o}/a$, in agreement with dimensional analysis, where $K_{S}=0.665$. The solution of the lubrication equation (4) for $h_{o}/a=43.42$ (dashed line) also exhibits a self-similar behavior Zhang and Lister (1999), but with a different asymptotic law $h_{\textrm{min}}\to K_{L}\,\tau^{1/5}$, where $K_{L}=0.913$. The crossover time between the lubrication and Stokes self-similar solutions can be estimated by equating $K_{L}\tau_{c}^{1/5}=K_{S}\tau_{c}^{1/3}\Rightarrow\tau_{c}=(K_{L}/K_{S})^{15/2}=10.77$, with an associated minimum thickness $h_{\text{min},c}=K_{S}(K_{L}/K_{S})^{5/2}=1.47$. Indeed, the results of Fig. 2 for $h_{o}/a=43.42$ reveal that the evolution of $h_{\textrm{min}}(\tau)$ obtained with the lubrication equation closely follows the Stokes result for $\tau\geq\tau_{1/5}\approx 10^{3}$ with a scaling exponent of $1/5$, followed by a long crossover for $\tau_{1/3}\leq\tau\leq\tau_{1/5}$, and finally reaching the $1/3$ power law for $\tau\leq\tau_{1/3}\approx 0.08$. In terms of the minimum film thickness, the $1/5$-scaling takes place for $h_{\textrm{min}}\geq 3.63$, the corrossover for $0.29\leq h_{\textrm{min}}\leq 3.63$, and the $1/3$-scaling for $h_{\textrm{min}}\leq 0.29$.

The failure of lubrication theory to predict the last stages of the rupture behavior observed in Fig. 2 demands unraveling the local self-similar Stokes flow, presented in the next section.

III.3 Self-similar solution

Dimensional analysis suggests substituting the similarity ansatz

[TABLE]

into (1)–(3) to elucidate the structure of the leading-order flow for $\tau\to 0$. The self-similar Stokes equations read

[TABLE]

which must be integrated in $0<\xi<\infty$, $0<\eta<f(\xi)$, with the boundary conditions $U=V=0$ at the wall $\eta=0$, $U=V_{\xi}=0$ at the symmetry plane $\xi=0$, $0<\eta<f(0)$, and

[TABLE]

at the unknown free surface $\eta=f(\xi)$. Notice that the leading-order contribution of the normal component of the stress balance (2) is $O(\tau^{-1})$, while the capillary pressure is $O(\tau^{-1/3})$. Thus, as anticipated by dimensional arguments, surface tension does not contribute to the normal-stress equilibrium at leading order as $\tau\to 0$. The system of nonlinear elliptic PDEs (19)–(24) describing the local Stokes flow close to rupture is parameter-free. It is interesting to note that problems with similar mathematical structure appear in other free-surface flows like inertial focusing and jet breakup Papageorgiou (1995); Zeff et al. (2000).

The asymptotic description of film rupture is completed by specifying the far-field boundary conditions at $\xi^{2}+\eta^{2}\gg 1$, $0<\eta<f(\xi)$. Inspection of the kinematic boundary condition (24) reveals that $f\propto\xi$ for $\xi\gg 1$ if $U$ and $V$ are subdominant. This suggests that the shape of the free surface sufficiently far from the origin, $r\gg 1$, is a wedge $\theta=\theta_{o}$, where $(r,\theta)$ are polar coordinates such that $\xi-\xi_{o}=r\cos\theta$ and $\eta=r\sin\theta$, and $(V_{r},V_{\theta})$ are the associated radial and polar components of the velocity field. Insight of the far-field behavior was first obtained by numerically integrating (19)–(24) imposing a stress-free boundary condition sufficiently far from the origin, as explained in Appendix C, where a detailed description of the numerical technique employed to solve the similarity problem is provided.

The examination of the radial and polar velocity profiles along rays $\theta=\textrm{constant}$ revealed that $V_{r}\sim V_{\theta}\sim r^{-(1+\lambda)}$ for $r\gg 1$, where $\lambda$ is a positive constant smaller than unity. This suggests a far-field stream function of the form Barenblatt and Zel’Dovich (1972) $\psi=F(\theta)/r^{\lambda}$, such that $V_{r}=r^{-1}\partial\psi/\partial\theta=F^{\prime}r^{-(1+\lambda)}$ and $V_{\theta}=-\partial\psi/\partial r=\lambda Fr^{-(1+\lambda)}$. Since $\psi$ is biharmonic, $F$ is seen to be the solution to the fourth-order linear homogeneous equation

[TABLE]

with the no slip condition $F=F^{\prime}=0$ at $\theta=0$, and the vanishing normal, $F^{\prime\prime\prime}+[4+3\lambda(2+\lambda)]F^{\prime}=0$ and tangential, $F^{\prime\prime}-\lambda(2+\lambda)F=0$ stress boundary conditions at $\theta=\theta_{o}$. If $\lambda$ were known, Eq. (25) together with the boundary conditions discussed above would constitute a closed eigenvalue problem for the universal angle $\theta_{o}$. However, $\lambda$ is expected to be determined from the asymptotic matching with a near-field description for $\xi\ll 1$, which is beyond the scope of this work.

We propose the following alternative: nontrivial solutions to (25) satisfying the boundary conditions exist if

[TABLE]

which determines the pairs $(\theta_{o},\lambda)$ classifying the family of allowed far-field solutions of the Stokes equations with a free-surface angle $\theta_{o}$. Thus one may i) extract $\theta_{o}$ from a numerical integration of (19)–(24) with a stress-free far-field boundary condition and then obtain $\lambda$ from (26), and ii) repeat the integration now imposing the far field variables entailed in the description of $F$. This iterative process converges very fast, and is stopped when the successive values of $\theta_{o}$ differ less than a prescribed tolerance, providing the pair of values $(\theta_{o},\lambda)=(37^{\circ},0.19)$.

The universal function $f(\xi)$ is represented in Fig. 3 together with rescaled film profiles $h/\tau^{1/3}$ for $h_{o}/a=4.34$, along with the far-field behavior $f=\xi\tan\theta_{o}$. The local expansion of the function $f$ for $\xi\ll 1$ has the form $f=f(0)+f^{\prime\prime}(0)\xi^{2}+\mathcal{O}(\xi^{4})$ due to the symmetry of the interface, with coefficients $f(0)=0.66$ and $f^{\prime\prime}(0)=0.59$. As a final self-similarity test, Fig. 4 displays isocontours of $U(\xi,\eta)$ obtained from the solution of (19)-(24) (right), and from the temporal integration of (1)–(3) for $\log_{10}{\tau}=-4.62$ (left), while the inset shows the corresponding profiles of $x$– and $y$–velocity along the interface.

IV Conclusions and future prospects

The self-similar solution obtained herein under Stokes flow arises from a balance between viscous and van der Waals forces, while surface tension forces are asymptotically negligible as rupture is approached. Although this result is in agreement with the conclusions reached by Burelbach et al. (1988) for isothermal liquid films using lubrication theory, a more precise analysis by Zhang and Lister (1999) demonstrated that, in fact, the self-similar solution of the lubrication equation (4) for $\tau\to 0$ arises from a balance between viscous, van der Waals and surface tension forces. In particular, it must be emphasized that the self-similar solution found by Zhang and Lister (1999) is inconsistent with the small-slope assumption, since the longitudinal length scales with time-to-rupture as $x\propto\tau^{2/5}$, while the minimum film thickness scales as $h_{\textrm{min}}\propto\tau^{1/5}$. Thus, the slenderness of the film scales as $h_{\textrm{min}}/x\propto\tau^{-1/5}\to\infty$ as $\tau\to 0$, invalidating the long-wavelength assumption near rupture. The latter inconsistency of the lubrication approach was already pointed out by Zhang and Lister (1999), who estimated that the slenderness assumption breaks down when $h^{*}_{\textrm{min}}\sim a$. Our analysis of the Stokes equations has revealed that, in fact, lubrication theory fails to describe the evolution of the film for $h^{*}_{\textrm{min}}\lesssim 4a$. In particular, the shape of the free surface near rupture approaches a wedge $f\to\xi$ with a slope ${\rm d}f/d\xi=\tan{\theta_{o}}\approx 0.75$ for $\xi\gg 1$ according to our Stokes description, while $f\to\xi^{1/2}$ according to Eq. (4) (Zhang and Lister, 1999). Although the leading-order lubrication theory does not describe the last stages of the vdW-induced rupture of thin films correctly, a higher-order theory might yield more accurate results Ruyer-Quil and Manneville (2000). Although the lack of slenderness of the local Stokes solution may invalidate the use of higher-order lubrication theory, this line of research clearly deserves further inquiry.

It is important to note that, regardless of the self-similar nature of the rupture dynamics, lubrication theory predicts an evolution for the liquid film that is markedly different from the Stokes description, as clearly evidenced by Fig. 2. Indeed, this difference appears during the early stages after the onset of the vdW instability, and increases over time. In particular, the 1/5 power law predicted by lubrication theory is only accomplished transiently during a very short intermediate time interval prior to the crossover to the 1/3 power law described here for the first time.

Admittedly, the new self-similar solution presents a number of limitations, which we aim to identify below in terms of direct comparisons with experimental data available in the literature of thin films. Indeed, the self-similar lubrication and Stokes regimes prevail for $h^{*}_{\textrm{min}}>3.63\,a$ and $h^{*}_{\textrm{min}}<0.29\,a$, respectively. Taking, for instance, $A=2.2\times 10^{-20}$ J, $\sigma=3.8\times 10^{-2}$ Jm*-2* as representative values for ultrathin polymer films Becker et al. (2003) provides $a=1.75$ Å. For $h_{o}=4$ nm Becker et al. (2003), the value of $h_{o}/a=22.82$ which, according to Fig. 2$(b)$, corresponds to a case where the self-similar lubrication regime is not established. The self-similar Stokes regime would be reached for $h^{*}_{\textrm{min}}<0.5$ Å, at which the continuum description is not valid. Additionally, taking the values $A=3.90\times 10^{-20}$ J, $\sigma=2.07\times 10^{-2}$ Jm*-2*, given by Reiter et al. (1999) yield $a=3.16$ Å, which for a range of initial thicknesses $30$ nm $<h_{o}<110$ nm, provide values $95<h_{o}/a<348$. Although in the latter experiments the self-similar lubrication regime should be established, the computational cost quickly increases with increasing values of $h_{o}/a$. Although the numerical exploration of higher values of $h_{o}/a$ might prove useful to quantify the differences between the slender theory and the fully two-dimensional theory, it is expected that the evolution would follow the trends shown in Figs. 2(b)–(c) for the highest value of $h_{o}/a$. In particular, for $\tau>\tau_{1/5}$, both the lubrication and the Stokes evolutions would match before departing towards the 1/3-regime for intermediate crossover times $\tau_{1/3}<\tau<\tau_{1/5}$. More importantly, short-range intermolecular forces, not taken into account in the present analysis, would become important at larger values of $h^{*}_{\textrm{min}}$. It is therefore concluded that, in the case of ultrathin polymer films, only the lubrication self-similar regime is experimentally realized, but only for values of $h_{o}/a\gtrsim 100$.

Liquefied metal films have been comparably less studied than polymer films, but there are a few experimental studies trying to characterize their instability and dynamics –see the review by Kondic et al. (2019) and references therein. For instance, uniform thin films of Cu, Ag, and Ni of thicknesses 4 nm $<h_{o}<$ 10 nm, liquefied by laser irradiation over a solid substrate, are studied in González et al. (2013); Krishna et al. (2009); McKeown et al. (2012), respectively. Although there are still some uncertainties in the measurement of the Hamaker constant González et al. (2013, 2016); Kondic et al. (2019), the latter experimental studies provide the values $A^{Cu}=2.58\times 10^{-18}$ J, $A^{Ag}=4.78\times 10^{-18}$ J, and $A^{Ni}=3.8\times 10^{-18}$ J, which, considering $\sigma^{Cu}=1.78$ Jm*-2*, $\sigma^{Ag}=0.93$ Jm*-2* and $\sigma^{Ni}=1.71$ Jm*-2*, yield $h_{o}/a^{Cu}=14.43$, $h_{o}/a^{Ag}=7.66$, and $h_{o}/a^{Ni}=11.65$, for a film thickness of $h_{o}=4$ nm. According to the findings presented herein, lubrication theory would provide an incorrect description for these three particular cases, as evidenced by figure 2($b$). The self-similar Stokes regime would be reached for $h_{\min}^{*,Co}<0.8$ Å, $h_{\min}^{*,Ag}<1.51$ Å, and $h_{\min}^{*,Ni}<1.0$ Å, respectively, thereby compromising the continuum approximation. More importantly, as discussed in detail by González et al. (2016); Kondic et al. (2019), in the case of liquefied metal films the liquid inertia can become relevant during film drainage and rupture. The latter fact will surely alter the self-similar structure obtained herein using the Stokes equations, as already shown by Garg et al. (2017) for liquid films with power-law rheology.

An alternative and promising theoretical approach to characterize the flow of thin films at the nanoscale is performing molecular-dynamics (MD) simulations Fowlkes et al. (2012); Nguyen et al. (2012, 2014). In particular, it was shown in Nguyen et al. (2012) that MD simulations and the weak-slip lubrication model Münch et al. (2005); Fetzer et al. (2005) yield almost identical results in different geometries and flow configurations of metal liquid films when the thickness is about $10$ nm. In particular, these authors stressed that the slip length is essential to obtain consistent predictions between both frameworks. Moreover, Moseler and Landman (2000); Eggers (2002) showed, in the case of liquid nanojets, that the addition of thermal fluctuations to the deterministic lubrication equation is crucial to obtain results in agreement with MD simulations. As for thin liquid films, there are a few theoretical and computational studies comparing the linear and nonlinear dynamics of the stochastic lubrication approximation with its deterministic counterpart Davidovitch et al. (2005); Grün et al. (2006); Fetzer et al. (2007); Nesic et al. (2015). However, the implication of thermally triggered fluctuations on the linear and nonlinear regimes within the Stokes framework remains as an open problem, to be compared with results obtained from the stochastic lubrication theory and MD simulations.

Another poorly understood aspect that deserves further study is the characteristic lifetime of the ultrathin film. According to Eq. (17), the Stokes equations predict slightly longer rupture times than lubrication theory. However, as far as we know, a systematic and rigorous comparison of the theoretically predicted lifetimes with the experimental observations is still missing. Indeed, a mismatch between the rupture times between the experiments and the numerical simulations of the lubrication equation in polymer films was reported in Becker et al. (2003). However, it remains unclear if this disagreement may be attributed to thermal noise Mecke and Rauscher (2005); Grün et al. (2006); Fetzer et al. (2007), to slippage effects Kargupta et al. (2004); Nguyen et al. (2012), or to a poor characterization of the intermolecular interactions Pahlavan et al. (2018).

Natural extensions of our work include for instance the effect of wall slip (Peschka et al., 2019), the study of axisymmetric rupture (Zhang and Lister, 1999; Witelski and Bernoff, 2000), the breakup of free films Vrij and Overbeek (1968); Vaynblat et al. (2001); Champougny et al. (2017), the influence of surfactants Edwards and Oron (1995); Warner et al. (2002), liquid-liquid dewetting Wyart et al. (1993); Pototsky et al. (2004); Ward (2011), the influence of polymer rheology Blossey (2012), and the effect of thermal noise Davidovitch et al. (2005); Grün et al. (2006); Fetzer et al. (2007); Nesic et al. (2015). Inertial effects should also be considered, as required to account for the dynamics of liquid metal films (González et al., 2013, 2016; Kondic et al., 2019). Finally, of particular importance is the inclusion of more detailed models of intermolecular interactions Seemann et al. (2001); Thiele et al. (2001); Israelachvili (2011); Pahlavan et al. (2018), as required to account for the resulting dewetting patterns and their long-term coarsening Becker et al. (2003); Glasner and Witelski (2003).

Acknowledgements.

The authors thank J. Rivero-Rodríguez and B. Scheid for key numerical advice, A.L. Sánchez, C. Martínez-Bazán and J.M. Gordillo for their enduring teaching and encouragement, W. Coenen for carefully reading the manuscript, and H.A. Stone for sharing his insights, and for his kind hospitality at the Complex Fluids Group, Princeton. This research was funded by the Spanish MINECO, Subdirección General de Gestión de Ayudas a la Investigación, through project DPI2015-71901-REDT, and by the Spanish MCIU-Agencia Estatal de Investigación through project DPI2017-88201-C3-3-R, partly financed through FEDER European funds. A.M.-C. also acknowledges support from the Spanish MECD through the grant FPU16/02562 and to its associated program Ayudas a la Movilidad 2018 during his stay at Princeton University.

Appendix A Inertial effects

In this appendix we briefly assess the accuracy of the assumption of negligible inertial effects. The relative importance of liquid inertia compared with the viscous forces is given by the local Reynolds number

[TABLE]

where $\rho$ is the liquid density and $\rho a\sigma/\mu^{2}$ is the Laplace number based on the liquid properties and the molecular length scale $a=[A/(6\pi\sigma)]^{1/2}$. Since $h_{\textrm{min}}=0.665\,\tau^{1/3}$ for $\tau\to 0$, ${\rm d}h_{\textrm{min}}/{\rm d}\tau=0.222\,\tau^{-2/3}$, and the local Reynolds number,

[TABLE]

However, the inertial crossover time, $\tau_{I}$, defined by the condition $Re(\tau_{I})=1$, accomplishes

[TABLE]

Taking $\rho=10^{3}$ kg m*-3*, $A=10^{-19}$ J, $\sigma=0.03$ Jm*-2*, $\mu=1$ kg m${}^{-1}\,$s*-1* provides an inertial crossover time of $\tau_{I}\approx 1.6\times 10^{-22}$, at a minimum radius $0.665\,\tau_{I}^{1/3}\approx 3.6\times 10^{-8}$, which is eight orders of magnitude smaller than the molecular length scale. Similar conclusions where obtained by Zhang and Lister (1999) in their lubrication analysis.

It is thereby deduced that inertia can be safely neglected in the analysis of the van der Waals instability of thin polymer liquid films. By way of contrast, it should be noted that in situations concerning the thinning of liquid metal films, care must be taken in retaining the inertial effects, as discussed in González et al. (2013, 2016); Kondic et al. (2019). For instance, the global Reynolds numbers considered in the second reference scale up to $\sim 10^{4}$, to be compared with that from Becker et al. (2003), on the order of $10^{-15}$, for which the estimations given in the previous paragraph are representative.

Appendix B Near-rupture velocities

Let us now provide an estimation of the typical velocities occurring near rupture. To that end, we take the fluid properties given in A, and a value of $\tau$ at which the short-range repulsive forces balance the long-range attractive vdW forces. The latter balance requires introducing a more realistic version of the potential, not considered in the main text, which according to (Becker et al., 2003) may be written in dimensionless form as

[TABLE]

where the second term includes the strength of the short-range part of the potential, namely, $\epsilon=6.3\times 10^{-76}\,\textrm{Jm}^{6}$. An estimation of the dimensionless precursor film thickness would be given by balancing both terms in (30) resulting in $h_{\textrm{prec}}=\frac{1}{a}(48\pi\epsilon/A)^{1/6}$, which yields a value $h_{\textrm{prec}}\simeq 2.36$. According to 2(a), this value of $h$ would lie well within the crossover region between the 1$/$5 and the 1$/$3 regimes, although the associated value of $\tau$ would hardly be modified. For instance, assuming the validity of the 1/3 power-law, the corresponding value of $\tau\simeq 44.6$ which yields dimensionless velocities on the order of $0.08$ according to (7) assuming $O(1)$ values of the self-similar variables $U$ and $V$. Recasting these into dimensional form provides velocities on the order of $0.08\sigma/\mu=2.4\times 10^{-3}\,\textrm{m}/\textrm{s}$.

Appendix C Numerical techniques

In this appendix we describe the numerical techniques employed to solve the Stokes equations (1), the lubrication equation (4), and the self-similar Stokes problem (19)–(24). All the equations were written in weak form in terms of the corresponding inner products, upon convenient use of Green’s identities, rendering them amenable for the use of finite elements for the spatial discretization. Comsol COM (1998) was the software of choice for the implementation of such numerical techniques.

C.1 The Stokes equations

The dimensionless Stokes equations were written in weak form using suitable test functions for the velocity and pressure, namely $\tilde{\bm{u}}$ and $\tilde{p}$, as follows:

[TABLE]

where $\Omega$ is the deformable computational domain bounded by $\Gamma_{s}\cup\Gamma_{w}\cup\Gamma_{a}\cup\Gamma_{\infty}$, the surface gradient operator is defined as $\nabla_{s}=(\bm{I}-\bm{n}\bm{n})\cdot\nabla$, $\mathrm{d}\Omega$ is the volume element and $\mathrm{d}S$ is the surface element. The remaining boundary conditions are specified in the sketch of Fig. 5.

Equation (31) was discretized using Taylor-Hood triangular elements for pressure and velocity (and the corresponding test functions) to ensure numerical stability. The use of an arbitrary Lagrangian-Eulerian (ALE) technique for the tracking of the interface allowed us to impose the kinematic boundary condition along $\Gamma_{s}$ by prescribing the normal velocity of the mesh. The displacement of the mesh elements was computed by solving a Laplace equation for the displacement field, namely $\nabla^{2}\bm{q}=0$, $\bm{q}=(q_{x},q_{y})$, with suitable boundary conditions.

The temporal discretization of the system of nonlinear equations was carried out using a variable-step BDF method with $2/4$ variable order. The relative tolerance of the nonlinear solver was set to $10^{-7}$. The initial conditions corresponded to a quiescent state $u=v=0$ and a perturbed interface $h(x,0)=h_{o}/a\left(1-10^{-3}\cos{kx}\right)$. The time-dependent solver was complemented with an automatic remeshing algorithm which redistributed the mesh elements when the deformation of the domain become sufficiently large. The numerical limitations of the computational techniques employed in this study precluded the minimum film thickness $h_{\textrm{min}}=h(0,t)$ to decrease below $\approx 10^{-4}$.

C.2 The lubrication equation

The lubrication equation for the temporal evolution of the film thickness (4), which is herein written again for convenience,

[TABLE]

needs to be integrated with boundary conditions $h_{x}(0,t)=h_{x}(\pi/k,t)=0$, and initial condition $h(x,0)=\frac{10}{9}\frac{h_{o}}{a}\left(1-\frac{1}{10}\cos{kx}\right)$, as used in Zhang and Lister (1999), which was integrated by recasting it into the conservative second-order form

[TABLE]

and discretized using quadratic Lagrange elements over a non-uniform partition of the domain with maximum element size of $10^{-7}$ close to the origin $x=0$. The time-stepper algorithm was the same as that used for the Stokes equations, namely a $2/4$-BDF.

The self-similar scaling of the film shape $h/(t_{\mathrm{R}}-t)^{1/5}=H(\eta)$ with $\eta=x/(t_{\mathrm{R}}-t)^{1/5}$ (note that this $\eta$ is different from that appearing in the main text) recovered the values obtained by Zhang and Lister (1999), namely $H(0)=0.7327$ and $H^{\prime\prime}(0)=0.3010$ as seen in Fig. 6, which served as a benchmark for our numerical scheme.

C.3 The self-similar Stokes problem

The elliptic system of PDEs for the self-similar variables $U$, $V$, $P$ and the a priori-unknown interface shape $f(\xi)$, was discretized using similar techniques as those employed for the complete Stokes equations. The associated weak form is

[TABLE]

where $\bm{U}=U\bm{e}_{\xi}+V\bm{e}_{\eta}$ and the del operator acts on the similarity coordinates $\xi$ and $\eta$. The boundary condition along the unknown free surface is that of vanishing stresses

[TABLE]

which is naturally accommodated in weak form. The ALE method was used again in this time-independent computation but this time with an extra degree of freedom, namely, the vertical displacement $\delta$ of the free surface $\Gamma_{s}$. This displacement was leveraged as a Lagrange multiplier for the imposition of the kinematic boundary condition by solving the weak boundary PDE

[TABLE]

where $\tilde{\delta}$ is the test function for the vertical displacement $\delta$ discretized using first-order Lagrange elements and $\mathrm{d}\Gamma$ is the line element along $\Gamma_{s}$. This additional equation enabled us to use a standard Newton–Raphson root-finding algorithm to iteratively solve for all the unknowns upon a tolerance, fixed to $10^{-5}$, provided a suitable initial guess. Such an initial guess is depicted in Fig. 7, comprising a quadrilateral computational domain of vertical dimensions $\eta_{0}$ and $\eta_{\infty}$, horizontal span $\xi_{\max}$ and initial angle $\theta_{o}^{(0)}$ and vanishing velocities $U=V=0$. A convenient choice for the tentative angle and $\eta_{0}$ was provided by the scaled film shapes obtained from the full temporal evolution of the Stokes equations near the breakup singularity. With this choice of the numerical parameters, the iterative algorithm described in the main text converged in a few iterations.