Oscillatory thermocapillary instability of a film heated by a thick substrate

TL;DR

This paper investigates oscillatory thermocapillary instabilities in heated liquid films, revealing conditions under which these instabilities occur, especially emphasizing the role of substrate thickness and thermal properties in their emergence.

Contribution

It introduces a nonlinear PDE model coupling film dynamics with substrate heat transfer, analyzing oscillatory instability conditions via a non-self adjoint eigenvalue problem.

Findings

Oscillatory instabilities are more prominent with insulating substrates.

Large substrate thicknesses can lead to complex eigenvalues indicating oscillatory behavior.

The model provides a comprehensive stability map for different parameters.

Abstract

In this work we consider a new class of oscillatory instabilities that pertain to thermocapillary destabilization of a liquid film heated by a solid substrate. We assume the substrate thickness and substrate-film thermal conductivity ratio are large so that the effect of substrate thermal diffusion is retained at leading order in the long-wave approximation. As a result, system dynamics are described by a nonlinear partial differential equation for the film thickness that is nonlocally coupled to the full substrate heat equation. Perturbing about a steady quiescent state, we find that its stability is described by a non-self adjoint eigenvalue problem. We show that, under appropriate model parameters, the linearized eigenvalue problem admits complex eigenvalues that physically correspond to oscillatory (in time) instabilities of the thin film height. As the principal results of our work, we provide a complete picture of the susceptibility to oscillatory instabilities for different model parameters. Using this description, we conclude that oscillatory instabilities are more relevant experimentally for films heated by insulating substrates. Furthermore, we show that oscillatory instability where the fastest-growing (most unstable) wavenumber is complex, arises only for systems with sufficiently large substrate thicknesses.

Tables

Table 1. Table 1: Physical properties for silicone oil, copper and PMMA (resp. Hintz et al. ( 2001 ) , Araki et al. ( 1992 ) , and Assael et al. ( 2005 ) ). The two values of 𝒬 𝒬 \mathcal{Q} are obtained by pairing each substrate with the silicone oil and estimating the heat transfer coefficient.

		silicone oil
$\mu$	(kg/m$\cdot$s)	4.94$\times {10}^{-4}$
${\kappa }_f$	(W/m$\cdot$K)	0.1
${\sigma }_0$	(kg/s2)	1.59$\times {10}^{-2}$
$\gamma$	(kg/s${}_{}{}^2\cdot$K)	6.4$\times {10}^{-5}$
${\theta }_0$	(K)	293

Equations

$$\rho \left({\partial }_t\mathbf{v}+\mathbf{v}\cdot \nabla \mathbf{v}\right)$$

$$\nabla \cdot \mathbf{v}$$

$${\partial }_t\theta +\mathbf{v}\cdot \nabla \theta$$

$${\partial }_{t}h+\mathbf{v}\cdot \nabla h=w\text{at}y=h,$$

$${\theta }_{\mathrm{h}}={\theta }_{\mathrm{h}}(x,t)\equiv \theta (x,h,t),$$

$$p_{\mathrm{g}}-p+\text{mathbfsf}\mathbf{T}\cdot \mathbf{n}\cdot \mathbf{n}=-2\mathcal{H}\sigma ({\theta }_{\mathrm{h}})$$

$$\text{mathbfsf}\mathbf{T}\cdot \mathbf{n}\cdot \mathbf{t}=\nabla \sigma ({\theta }_{\mathrm{h}})\cdot \mathbf{t}$$

$$\mathbf{n}=\frac{\mathbf{k}-{\partial }_{x}h\mathbf{i}}{\sqrt{1+{\left({\partial }_{x}h\right)}^2}}\mathbf{t}=\frac{\mathbf{i}+{\partial }_{x}h\mathbf{k}}{\sqrt{1+{\left({\partial }_{x}h\right)}^2}},$$

$${\kappa }_{\mathrm{f}}\nabla \theta \cdot \mathbf{n}+q\left({\theta }_{\mathrm{h}}-{\theta }_g\right)=0\text{at}y=h,$$

$$\theta =\psi$$

$${\partial }_y\theta =\kappa {\partial }_z\psi$$

$$\kappa ={\kappa }_{\mathrm{s}}/{\kappa }_{\mathrm{f}}$$

$${\partial }_t\psi ={\chi }_{\mathrm{s}}{\nabla }_s^2\psi ,$$

$$\psi ={\psi }_{\mathrm{b}}\text{at}z=-d,$$

$$\mathcal{B}=\frac{q\overline{h}}{{\kappa }_{\mathrm{f}}},\mathcal{S}=\frac{qd}{{\kappa }_{\mathrm{s}}},\mathcal{C}=\frac{\gamma \Delta }{{\sigma }_0}\frac{{\kappa }_{\mathrm{s}}^2}{{\kappa }_{\mathrm{f}}^2},\mathcal{Q}=\frac{q\mu {\chi }_{\mathrm{s}}}{{\sigma }_0}\frac{{\kappa }_{\mathrm{s}}^2}{{\kappa }_{\mathrm{f}}^3},$$

$$\left.\begin{aligned} \begin{array}[]{lllll}x^{\prime}=d,&z^{\prime}=d,&y^{\prime}=\dfrac{\kappa_{\rm f}}{q},&t^{\prime}=\dfrac{d^{2}}{\chi_{\rm s}},&h^{\prime}=\dfrac{\kappa_{\rm f}}{q},\\ u^{\prime}=\dfrac{x^{\prime}}{t^{\prime}}=\dfrac{\chi_{\rm s}}{d},&w^{\prime}=\dfrac{y^{\prime}}{t^{\prime}}=\dfrac{\chi_{\rm s}\kappa_{\rm f}}{d^{2}q},&p^{\prime}=\dfrac{q^{2}\,\mu\,\chi_{\rm s}}{\kappa_{\rm f}^{2}},&\theta^{\prime}=\Delta,&\psi^{\prime}=\Delta.\end{array}\end{aligned}\right\}$$

$$\left.\begin{aligned} \begin{array}[]{ccrlcllcl}u&=&u^{\prime}&(U_{0}&+&\varepsilon\,U_{1}&+\dots+\varepsilon^{n}\,U_{n})&&\\ w&=&w^{\prime}&(W_{0}&+&\varepsilon\,W_{1}&+\dots+\varepsilon^{n}\,W_{n})&&\\ p&=&p^{\prime}&(P_{0}&+&\varepsilon\,P_{1}&+\dots+\varepsilon^{n}\,P_{n})&+&p_{\rm g}\\ h&=&h^{\prime}&(H_{0}&+&\varepsilon\,H_{1}&+\dots+\varepsilon^{n}\,H_{n})&&\\ \theta&=&\theta^{\prime}&(\mathit{\Theta}_{0}&+&\varepsilon\,\mathit{\Theta}_{1}&+\dots+\varepsilon^{n}\,\mathit{\Theta}_{n})&+&\theta_{\rm g}\\ \psi&=&\psi^{\prime}&(\mathit{\Psi}_{0}&+&\varepsilon\,\mathit{\Psi}_{1}&+\dots+\varepsilon^{n}\,\mathit{\Psi}_{n})&+&\theta_{\rm g}\end{array}\end{aligned}\right\}$$

$${\partial }_Y^{2}U-{\partial }_{X}P=0$$

$${\partial }_{Y}P=0$$

$$P=-({\mathcal{S}}^2\mathcal{Q}{)}^{-1}{\partial }_X^{2}H$$

$${\partial }_{Y}U=-\mathcal{C}{\mathcal{Q}}^{-1}{\partial }_X{\Theta }_{\mathrm{H}}$$

$$U=(\mathcal{Q}{\mathcal{S}}^2{)}^{-1}H{\partial }_X^{3}H\left(YH-\frac{1}{2}{Y}^2\right)-\mathcal{C}{\mathcal{Q}}^{-1}Y{\partial }_X{\Theta }_{\mathrm{H}},$$

$$0={\partial }_{T}H+{\partial }_X{\int }_0^{H}UdY.$$

$${\partial }_{T}H+(\mathcal{Q}{\mathcal{S}}^2{)}^{-1}{\partial }_X\left\{\frac{1}{3}{H}^3{\partial }_X^{3}H-\frac{1}{2}(\mathcal{C}{\mathcal{S}}^2)H^2{\partial }_X{\Theta }_{\mathrm{H}}\right\}=0.$$

$${\partial }_Y^2\Theta =0$$

$${\partial }_T\Psi -{\partial }_X^2\Psi -{\partial }_Z^2\Psi =0$$

$${\partial }_Y\Theta +\Theta =0$$

$$\Theta -\Psi =0$$

$$\mathcal{S}{\partial }_Y\Theta -{\partial }_Z\Psi =0$$

$$\Psi =1$$

$$\left.\begin{aligned} \begin{array}[]{@{~}r@{~}c@{~}c@{~}ccl}\leavevmode\nobreak\ \lx@intercol\hfil H(X,T)=\leavevmode\nobreak\ &\mathcal{B}\hfil\leavevmode\nobreak\ &+\hfil\leavevmode\nobreak\ &\delta\;\hat{H}\cos{(\beta{X})}\exp(\Omega{T})&+&\mathcal{O}(\delta^{2})\\ \leavevmode\nobreak\ \lx@intercol\hfil\mathit{\Theta}_{\rm H}(X,T)=\leavevmode\nobreak\ &\bar{\mathit{\Theta}}_{\rm H}\hfil\leavevmode\nobreak\ &+\hfil\leavevmode\nobreak\ &\delta\;\hat{\mathit{\Theta}}_{\rm H}\cos{(\beta{X})}\exp(\Omega{T})&+&\mathcal{O}(\delta^{2})\\ \leavevmode\nobreak\ \lx@intercol\hfil\mathit{\Theta}(X,Y,T)=\leavevmode\nobreak\ &\bar{\mathit{\Theta}}({Y})\hfil\leavevmode\nobreak\ &+\hfil\leavevmode\nobreak\ &\delta\;\hat{\mathit{\Theta}}({Y})\cos{(\beta{X})}\exp(\Omega{T})&+&\mathcal{O}(\delta^{2})\\ \leavevmode\nobreak\ \lx@intercol\hfil\mathit{\Psi}(X,Z,T)=\leavevmode\nobreak\ &\bar{\mathit{\Psi}}({Z})\hfil\leavevmode\nobreak\ &+\hfil\leavevmode\nobreak\ &\delta\;\hat{\mathit{\Psi}}({Z})\cos{(\beta{X})}\exp(\Omega{T})&+&\mathcal{O}(\delta^{2}).\end{array}\end{aligned}\right\}$$

$${\bar{\Theta }}^{\prime \prime }=0$$

$${\bar{\Psi }}^{\prime \prime }=0$$

$${\bar{\Theta }}^{\prime }+\bar{\Theta }=0$$

$$\bar{\Theta }-\bar{\Psi }=0$$

$$\mathcal{S}{\bar{\Theta }}^{\prime }-{\bar{\Psi }}^{\prime }=0$$

$$\bar{\Psi }=1$$

Full text

1 Introduction

The tendency of thin liquid films to destabilize and form wavy patterns is an important area of research for a wide range of applications. For some applications, such as coatings and glass manufacturing, one may wish to operate under conditions that avoid these instabilities. In others, such as multiphase heat/mass transfer technology and nanoscale patterning of liquid metals/polymers, precise control of the emerging wave pattern is of utmost concern. In either case, the parametric conditions of interest can be determined, most simply, by applying the long-wave approximation to the governing nonlinear equations, see Oron et al. (1997) and Craster and Matar (2009). In the long-wave approach, physical effects such as gravity, mean surface tension, thermocapillarity, solutocapillarity, and electromagnetism can be easily be accounted for, and one typically obtains a single nonlinear partial differential equation (PDE) for the spatiotemporal evolution of the local film thickness. This method assumes the film dynamics are non-inertial and governed by a (first order in time) nonlinear PDE.

The principal phenomenon that a single-equation long-wave model cannot describe is the emergence of instabilities that are oscillatory in time, i.e., overstability (see Nepomnyashchy et al. (2001), chapter 5). Whereas single-equation film models predict monotonic perturbations that grow or decay exponentially in time, oscillatory instabilities can only be observed in systems that describe the interaction between processes that occur on distinct time scales. Thus, oscillatory instabilities are commonly obtained from Orr-Somerfeld type analyses of governing equations of motion that retain inertial effects and diffusive time scales. Wide-ranging examples that highlight the emergence of oscillatory instabilities in fluid layers include work by Sternling and Scriven (1959), Takashima (1981), Anderson and Worster (1996), and Rednikov et al. (1998). A common theme to these works is the level of analytical difficulty; each obtains a linear dispersion relation (describing system stability) that is transcendental and implicit in the perturbation growth rates. Combined with large parametric spaces and the fact that oscillatory perturbations necessarily reside in the complex plane, concise description of the emergence of oscillatory instability can be a challenging task. Alternatively, the long wave approximation offers a convenient means to couple free surface deformation to other time-dependent physical processes of interest.

Several authors have investigated oscillatory instabilities of thin liquid films in the context of the long-wave approximation. In many cases, e.g. Podolny et al. (2005) and Bestehorn and Borcia (2010), such instabilities originate from the coupling between the local thickness and bulk concentration of a film composed of a binary mixture. In addition to the bulk concentration dynamics, Morozov et al. (2014) investigated oscillatory instability with the added effect of absorption/desorption kinetics between interfacial and bulk film surfactant concentration. In other cases, oscillatory instabilities have been uncovered in multiple stacked layers of films, as described theoretically by coupled sets of film thickness evolution equations (Nepomnyashchy and Simanovskii (2007), Beerman and Brush (2007)). Multi-layer film configurations do not, however, guarantee oscillatory modes: for example, such instabilities were not obtained by Pototsky et al. (2005) who investigated the dewetting dynamics of isothermal, ultrathin bilayers. Of particular interest to the present work are oscillatory instabilities reported by Shklyaev et al. (2012) in a model of thin-film thermocapillary destabilization from below. While there are similarities between that work and the present, we point out one important difference: in Shklyaev et al. (2012), the instability is driven by imposing a heat flux at the film-substrate interface; instead, in the present work we consider the full time-dependent heat-transfer in the substrate. We also note that each of these works on oscillatory instabilities of thin liquid films obtains low-order polynomial equations for the perturbation growth rates (in contrast to the transcendental, implicit dispersion we obtain in the present work).

The problem we investigate in this work is the deformational thermocapillary instability. This classic long-wavelength instability was first introduced by Scriven and Sternling (1964) and later verified experimentally by VanHook et al. (1997). In short, thermocapillary stresses that destabilize the free surface are generated by heating a film from below (transverse heating). For sufficiently thin layers, these stresses can surpass capillary stabilization and deform an initially flat film. Recently, Dietzel and Troian (2009) connected this mechanism with the formation of nanopillars ($\sim$10 $\mu$m spacing) on ultrathin ($\sim$100 nm) polymer films. Continuing work on thermocapillary patterning in (ultra)thin polymer films has been reviewed by Singer (2017). A patterning application that directly motivates our study is pulsed-laser dewetting of nanometric liquid metal films. Experiments by Trice et al. (2007) demonstrated dewetting pattern wavelengths that were commensurate with the predictions of long-wavelength thermocapillary modes. Driven by such results, several workers have developed and investigated theoretical models for the pulsed-laser process. Atena and Khenner (2009) augmented a long-wave theory with pulsed laser irradiation to describe the thermocapillary dewetting of liquid cobalt on silicon oxide substrates. Notably, they assumed that the substrate was thin, thereby ensuring model dynamics could be described by a (first order in time) single nonlinear PDE for the film thickness. As a result, oscillatory instabilities do not arise in their model.

Recently, oscillatory modes for pulsed-laser thermocapillary dewetting of liquid metal films were uncovered by Dong and Kondic (2016) and Seric et al. (2018). These authors made observations primarily via nonlinear simulations of a model that couples the film PDE to the full heat equation for the substrate. These works have not precisely characterized the emergence of oscillatory instabilities, in particular, because the task is complicated by the parameter space introduced by laser heating. Thus, in the present work, we investigate the emergence of oscillatory instability for the simpler problem: a film heated by a thick solid substrate. To do so, we initially follow the work by Saeki et al. (2011, 2013) that considered the linear analysis of a coupled film-substrate model, which induces thermocapillary film deformation driven by laser heating. In the present work we follow their asymptotic assumptions so that the full heat equation of the substrate is retained at leading order in the long-wave expansion of the governing equations. Effectively, we assume that the substrate-film thermal conductivity and thickness ratios are large. Although we obtain a dispersion relation that is similar to that of Saeki et al. (2013), it is important to note that they did not observe oscillatory modes. This may be due to a limited examination of model parameter values in their investigation.

The manuscript proceeds as follows: in §2, we present the dimensional equations of motion and boundary conditions for a deformable liquid layer heated by a thick substrate. In §3 we introduce a long-wave asymptotic expansion and derive an evolution equation for the film thickness that is nonlocally coupled to the diffusive (time-dependent) heat conduction problem in the substrate. In this section we also introduce a unique nondimensionalization that casts the nonlocal model in terms of four dimensionless parameters: (1) $\mathcal{B}$, characterizing the mean thermal thickness of the film; (2) $\mathcal{S}$, characterizing the thermal thickness of the substrate; (3) $\mathcal{C}$, characterizing the imposed temperature difference; (4) $\mathcal{Q}$, depending only on material properties. In the following section, §4, we perform a linear analysis of the nonlocal model and demonstrate that its stability is governed by a generalized two-point boundary value problem that is not self-adjoint. Solution of this problem yields the (implicit) dispersion relation that sets the course of investigation for the remainder of the paper. In section §5 we characterize the root structure of the dispersion relation and introduce the numerical contintuation methods we use to track its roots as functions of the perturbation wavenumbers. Then, in section §6 we classify the two characteristic pathways by which oscillatory instability manifests itself. Finally, in section §7, we provide a complete picture of the emergence of oscillatory instabilities within the considered parameter space.

2 Dimensional equations

Here we introduce equations that describe the fluid and temperature dynamics of the laterally infinite, two dimensional film-substrate system depicted schematically in figure 1. The film is composed of a Newtonian, incompressible liquid with average thickness $\overline{h}$, density $\rho$, dynamic viscosity $\mu$, kinematic viscosity $\nu=\mu/\rho$, thermal conductivity $\kappa_{\rm f}$, and thermal diffusivity $\chi_{\rm f}$. Neglecting gravity, we have

[TABLE]

where $\bm{v}\equiv{}\{u(x,y,t),w(x,y,t)\}$, $p(x,y,t)$, and $\mathit{\theta}(x,y,t)$ are the film velocity, pressure and temperature fields, respectively. With $\bm{\nabla}=\{\partial_{x},\partial_{y}\}$, equations (2.1-2.3) govern the evolution of these fields in time $t$ on the horizontal domain $x\in(-\infty,\infty)$ and the vertical domain $y\in{[0,h]}$ where $h=h(x,t)$ is the local, instantaneous film thickness.

The dynamics of this system are decoupled from those of the gas phase by assuming that the ratios between the liquid and gas phase densities, viscosities, and thermal diffusivities are large. Accordingly, at the free surface, we have the kinematic condition

[TABLE]

which states that the speed of the free surface is equal to the velocity of the fluid. Using

[TABLE]

to denote the free surface temperature, the normal and tangential stress balances that hold at the free surface are

[TABLE]

respectively, with gas pressure $p_{\rm g}$, rate of deformation tensor $\bm{\mathbfsf{T}}=\mu\left[\bm{\nabla}\bm{v}+\left(\bm{\nabla}\bm{v}\right)^{\intercal}\right]$ in the liquid phase, surface normal and tangent unit vectors,

[TABLE]

and twice mean curvature $2{\mathcal{H}}=-\bm{\nabla}\cdot\bm{n}$. We consider fluids whose surface tension decreases linearly with temperature according to $\sigma(\mathit{\theta}_{\rm h}{}{})=\sigma_{0}-\gamma\left(\mathit{\theta}_{\rm h}{}{}-\mathit{\theta}_{0}\right)$ where $\gamma=-d\sigma/d\mathit{\theta}_{\rm h}{}{}$ is positive and $\mathit{\theta}_{0}$ is a reference temperature.

Variations in $\theta_{\rm h}$ leading to thermocapillary destabilization are driven by the heat exchanged with the bounding gas phase. This process is modeled using Newton’s Law of Cooling, viz.,

[TABLE]

where $\theta_{\rm g}$ is the uniform gas temperature and $q$ is the empirical rate of heat transfer between the surface and the gas.

The film temperature evolves according to (2.3), and, at $y=0$, the film is in thermal contact with a rigid substrate of temperature $\psi$, thermal conductivity $\kappa_{\rm s}$, and diffusivity $\chi_{\rm s}$. No-slip and no-penetration enforce $\bm{v}=\bm{0}$, and, continuity of temperature and heat flux require

[TABLE]

where

[TABLE]

is the conductivity ratio. Here, the vertical domain of $\psi$ is assigned to a second vertical coordinate $z\in{}[-d,0]$ in anticipation that two vertical length scales will be introduced in the asymptotic analysis of the thick substrate case. The evolution of $\mathit{\psi}$ throughout the substrate is governed by

[TABLE]

where $\bm{\nabla}_{s}=\{\partial_{x},\partial_{z}\}$ is defined with respect to $z$ and $\chi_{\rm s}=\kappa_{\rm s}(c_{s}\rho_{s})^{-1}$ is substrate thermal diffusivity.

Opposite the film, we assume the substrate is in perfect thermal contact with a blackbody of uniform temperature $\mathit{\mathit{\psi}}_{\rm b}$ and impose a Dirichlet condition,

[TABLE]

and define the temperature difference $\Delta\equiv\mathit{\theta_{\rm g}}-\mathit{\psi_{\rm b}}$. At the cost of introducing a second heat transfer coefficient, a mixed boundary condition accounting for interfacial resistances to heat transfer could also be imposed at $z=-d$. By assuming instead that the blackbody transfers heat efficiently to the substrate, the parametric burden of the model is lessened.

3 Dimensionless asymptotic model

In this section we perform a formal asymptotic expansion of the model in Section 2 that describes the evolution of long wavelength disturbances driven by thermally diffusive substrates. The asymptotic model will be written to depend on four dimensionless quantities,

[TABLE]

where the Biot numbers $\mathcal{B}$ and $\mathcal{S}$ measure the thermal thickness of the film and substrate, respectively, $\mathcal{C}$ measures the imposed temperature difference, and $\mathcal{Q}$ measures the combined effects of the film viscosity and substrate thermal diffusivity. These groups arise in the dimensionless asymptotic model if we choose the characteristic scales

[TABLE]

The scales for the film vertical coordinate and its thickness, $y^{\prime}=h^{\prime}=\kappa_{\rm f}/q$, are chosen to ensure that $\mathcal{B}$ arises as the mean value of the local dimensionless film thickness. A different vertical scale, $z^{\prime}=d$, is the natural choice for the substrate, so that the parameter $\mathcal{S}$ will enter into the non-dimensional version of the heat flux condition (2.11). Taking $x^{\prime}=d$ and $t^{\prime}=d^{2}/\chi_{\rm s}$ as the lateral length and time scales of the substrate, the lateral film velocity and pressure scales that follow from these choices are as given in (3.2), which leads to the emergence of the quantities $(\mathcal{S}^{2}\mathcal{Q})^{-1}$ and $\mathcal{C}\,\mathcal{Q}^{-1}$ in the dimensionless normal and tangential stress balances, respectively. Lastly, the film transverse velocity scaling $w^{\prime}$ is different than $u^{\prime}$ as a result of the different scaling choices for $x^{\prime}$ and $y^{\prime}$, and the (subsequent) nondimensionalization of the continuity equation (2.2).

To perform the asymptotic expansion, we first define the aspect ratio parameter $\varepsilon\equiv{}y^{\prime}/x^{\prime}=(\mathcal{S}\kappa)^{-1}$ and require that $\varepsilon\ll{1}$. Having set $z^{\prime}=x^{\prime}$, satisfying $\epsilon\ll{1}$ ensures we consider systems with mean film thicknesses that are small compared to both its lateral variations and the substrate thickness. This definition of $\epsilon$ (i.e., with respect to two system dimensions) contrasts conventional long-wavelength analyses that define $\epsilon$ as the ratio of the film thickness to a characteristic horizontal wavelength. As a result, in the current problem, $\epsilon$ arises naturally in the model following nondimensionalization with (3.2). The formal expansions of the dependent variables take the form

[TABLE]

where the variables subscripted with $n=0,1,...$ are dimensionless and assumed to be $O(1)$ in magnitude.

Substituting (3.3) into the governing equations (2.1-2.14), we retain the leading order terms, drop the [math] subscripts on the dependent variables, and obtain dimensionless long-wavelength, thick substrate equations and boundary conditions. Attending first to the film equations of motion, we have, from (2.1),

[TABLE]

a boundary value problem for $U$ and $P$ that is closed by applying the conditions $U=0$ at $Y=0$, and, from (2.6) and (2.7),

[TABLE]

at the free surface.

This boundary value problem (3.4)-(3.7) describes viscous, locally-parallel flows that may be driven by capillary normal stresses or thermocapillary tangential stresses at the free surface. Because the leading order vertical pressure gradient is equal to zero via (3.5), the horizontal pressure gradient appearing in (3.4) is independent of $Y$ and is evaluated using the interfacial value specified by (3.6). Solution of the boundary value problem for the horizontal velocity $U(Y)$ yields

[TABLE]

where $\mathit{\Theta}_{\rm H}(X,T)=\mathit{\Theta}(X,H(X,T),T)$ is the temperature at the free surface. From (2.2), we use the dimensionless equation for continuity to rewrite the kinematic condition as

[TABLE]

Evaluating the integral (3.9) using (3.8), we obtain a nonlinear partial differential equation for the spatiotemporal evolution of $H(X,T)$,

[TABLE]

This equation represents a standard model for the dynamics of a liquid film subject to capillary stabilization and thermocapillary destabilization. Via the free surface temperature $\mathit{\Theta}_{\rm H}(X,T)$, (3.10) is coupled to the long-wave counterparts to equations (2.9)-(2.14), viz.,

[TABLE]

which are subject to

[TABLE]

To summarize, equations (3.10)-(3.16) represent an asymptotic model that couples, via $\mathit{\Theta}_{\rm H}(X,T)$, a nonlinear partial differential equation for the film thickness $H(X,T)$ to a thermal boundary value problem for temperature profiles $\Theta(X,Y,T)$ and $\Psi(X,Z,T)$ in the film and substrate, respectively. The model can be recast without $\Theta(X,Y,T)$, given that (3.11) prescribes profiles $\Theta(X,Y,T)$ that are linear in $Y$. However we find it easier to present the linear analysis that follows by first perturbing the system as written in (3.10)-(3.16). We also note that the model can be recast to include conventional capillary and Marangoni numbers if (3.10)-(3.16) are instead nondimensionalized with respect to the viscous scales of the film. However, we find the parameter set $(\mathcal{C},\mathcal{Q},\mathcal{B},\mathcal{S})$ is most conducive to a complete presentation of oscillatory instabilities.

4 Linear analysis

In this section we present a linear stability analysis of small perturbations to a steady state solution of (3.10)–(3.16). The steady state solution, which we will also refer to as the basic state, consists of a horizontally uniform (i) flat film of constant height, and (ii) temperature profile that depends linearly on the vertical $Y$ and $Z$-coordinates. Notably, we demonstrate that the resulting linear equations can be cast as a generalized eigenvalue problem that is not self adjoint (in the standard $L^{2}$ inner product). The key result of the linear analysis is the determination of the dispersion relation that characterizes the perturbation growth rate $\Omega$ implicitly in terms of the wavenumber $\beta$. Approximate and numerical assessment of system stability as governed by the dispersion relation then sets the course of investigation for the remainder of the paper.

To proceed with the linear analysis we introduce a normal-mode perturbation to a steady state solution of (3.10)–(3.16),

[TABLE]

Here $\delta\ll{1}$ is the real amplitude of a horizontally-periodic perturbation of real wavenumber $\beta$ and complex growth rate $\Omega$; while we choose the functions $(\mathcal{B},\bar{\mathit{\Theta}}_{\rm H},\bar{\mathit{\Theta}}({Y}),\bar{\mathit{\Psi}}({Z}))$ to be a steady solution of (3.10)–(3.16). To determine the steady solutions first note that the functions $H(X,T)=\mathcal{B}$ and $\mathit{\Theta}_{\rm H}(X,T)=\bar{\mathit{\Theta}}_{\rm H}$ are both constants. The remaining equations (3.11)–(3.16) then govern the basic state temperature profiles for $\bar{\mathit{\Theta}}(Y)$ and $\bar{\mathit{\Psi}}(Z)$:

[TABLE]

subject to the interface and boundary conditions

[TABLE]

Solving the linear equations (4.2)–(4.7) yields the complete steady solution:

[TABLE]

Note that the value of $\bar{\mathit{\Theta}}_{\rm H}$ in (4.8) is determined from $\bar{\mathit{\Theta}}(Y)$ via $\bar{\mathit{\Theta}}_{\rm H}=\bar{\mathit{\Theta}}(\mathcal{B})$, since $\bar{\mathit{\Theta}}_{\rm H}$ is defined as the temperature profile $\bar{\mathit{\Theta}}(H)$ at $Y=H=\mathcal{B}$. Together equations (4.8) define the film and substrate temperatures of a horizontally-uniform basic state as linear functions of their respective vertical coordinates.

We now move to compute the $O(\delta)$ perturbation about the basic state. First note that the dependent variable $\mathit{\Theta}_{\rm H}(X,T)$ is just the value of $\mathit{\Theta}(X,Z,T)$ evaluated at the surface $Z=H$, i.e. $\mathit{\Theta}_{\rm H}(X,T)=\mathit{\Theta}(X,H,T)$. Hence, the perturbation variables $\bar{\mathit{\Theta}}_{\rm H}$ and $\hat{\mathit{\Theta}}({Y})$ for $\mathit{\Theta}_{\rm H}(X,T)$ and $\mathit{\Theta}(X,Z,T)$ are coupled. To them we (i) Taylor expand $\mathit{\Theta}(X,H,T)$ about the base state value $H=\mathcal{B}$ in powers of $\delta$, and (ii) equate the $O(\delta)$ terms in $\mathit{\Theta}(X,H,T)$ with those of $\mathit{\Theta}_{\rm H}(X,T)$. We then obtain the relation:

[TABLE]

where we have used the fact (from (4.8)) that $\bar{\mathit{\Theta}}^{\prime}(\mathcal{B})=-\bar{\mathit{\Theta}}_{\rm H}$. Equation (4.9) will be used to eliminate the variable $\hat{\mathit{\Theta}}_{\rm H}$ from the linear stability analysis.

To obtain the linearized equations about the basic state, we substitute the ansatz (4.1) into the long-wavelength model given by (3.10)–(3.16). Collecting the $O(\delta)$ terms, and using (4.9) to eliminate $\hat{\mathit{\Theta}}_{\rm H}$, gives rise to equations for $\hat{H}$, $\hat{\mathit{\Theta}}(Y)$, and $\hat{\mathit{\Psi}}(Z)$:

[TABLE]

In (4.11) we have introduced

[TABLE]

which plays the role of a (complex-valued) wavenumber in the $Z$-direction for perturbations confined to the substrate domain. Note that the sign convention assumed in (4.12) is intentional for the subsequent stability analysis. Equations (4.10)–(4.11) are also subject to the film dispersion relation

[TABLE]

where

[TABLE]

and the boundary conditions are

[TABLE]

To obtain non-zero solutions to (4.10)–(4.17), we first recast the system as an eigenvalue problem for $\hat{\mathit{\Psi}}(Z)$ by eliminating the variables $\hat{H}$ and $\hat{\mathit{\Theta}}(Y)$. To first eliminate $\hat{\mathit{\Theta}}(Y)$, we solve equation (4.10), writing $\hat{\mathit{\Theta}}(Y)=\hat{\mathit{\Theta}}^{\prime}\,Y$+$\hat{\mathit{\Theta}}(0)$ for constants $\hat{\mathit{\Theta}}^{\prime}$ and $\hat{\mathit{\Theta}}(0)$. Inserting the solution for $\hat{\mathit{\Theta}}(Y)$ into the two boundary conditions (4.13)–(4.14) allows one to solve for the constants $\hat{\mathit{\Theta}}^{\prime}$ and $\hat{\mathit{\Theta}}(0)$ in terms of $\hat{H}$ only. Writing $\hat{\mathit{\Theta}}(Y)$ in terms of $\hat{H}$, the interface conditions (4.15)–(4.16) then take the form:

[TABLE]

The variable $\hat{H}$ can be eliminated in the interface equations (4.18)–(4.19), yielding a boundary condition for $\hat{\mathit{\Psi}}(Z)$ at $Z=0$. The resulting boundary condition at $Z=0$, together with the ODE (4.11), and boundary condition (4.17) at $Z=-1$, gives rise to the following problem for eigenvalues $\lambda^{2}$ and eigenfunctions $\hat{\mathit{\Psi}}(Z)$:

[TABLE]

with real constants

[TABLE]

Note that (4.12) has been used to replace $\Omega$ in terms of $\lambda$ in (4.20). The problem (4.20) is irregular in the sense that the eigenvalue $\lambda^{2}$ appears in both the boundary condition as well as the domain equation. To solve for the eigenvalues, we write the general solution for $\hat{\Psi}(Z)$ as $\hat{\mathit{\Psi}}(Z)=c_{1}\lambda^{-1}\sinh{(\lambda{Z})}+c_{2}\cosh{(\lambda{Z})}$, and require that it satisfies the two boundary conditions in (4.20). We include the extra factor of $\lambda^{-1}$ in the ansatz so that $\lim_{\lambda\rightarrow 0^{+}}\hat{\mathit{\Psi}}(Z)=c_{1}Z+c_{2}$ solves the ODE (4.20) when $\lambda=0$ (this will then allow for the simultaneous treatment of $\lambda=0$ and $\lambda\neq 0$ in the subsequent calculations). Substitution then requires that the following determinant vanish, viz.,

[TABLE]

In the equation (4.22), the $\lambda^{-1}\tanh{\lambda}$ term has a removable singularity at $\lambda=0$ (with the limit value of $1$ when $\lambda\rightarrow 0$). Equation (4.22) may then be compactly written as an implicit function relating $\lambda$ and $\beta$:

[TABLE]

where

[TABLE]

Any solution $(\lambda,\beta)$ to equation (4.22), or equivalently (4.23), then determines $\Omega$ via equation (4.12). As a result, equation (4.23) defines an (implicit) dispersion relation since it describes the values of $\lambda$ (and hence $\Omega$), in terms of $\beta$, for which non-zero solutions $\hat{\Psi}(Z)$ exist. We will therefore refer to $f(\lambda,\beta)=0$ as the dispersion relation. For values of $\lambda\neq 0$ and $a_{1}+\lambda^{2}b_{1}\neq 0$, the dispersion relation $f(\lambda,\beta)=0$ can be recast into a form that is more commonly encountered in linear stability analyses of thin film models,

[TABLE]

In Appendix §A, this form is used to readily obtain a dispersion relation for thin substrates. The third term in (4.25) describes thermocapillarity and we note that it vanishes in situations that render the free surface isothermal (this occurs if $\mathcal{C}=0$).

To conclude the solution of the linearized system, we solve for the substrate temperature eigenfunction $\hat{\Psi}({Z})$, and the film temperature eigenfunction $\hat{\mathit{\Theta}}({Z})$ in terms of the film perturbation amplitude $\hat{H}$. For a fixed $\beta$, take $\lambda$ as a root of the dispersion relation (4.24) and fix $\Omega$ via (4.12). Then the ansatz (4.1) solves the linearized equations, with eigenfunction profiles given by:

[TABLE]

5 Root structure of the dispersion relation

In this section, we compute the growth factors $\Omega(\beta)=\lambda(\beta)^{2}-\beta^{2}$ by solving the dispersion relation $f(\lambda,\beta)=0$ for $\lambda$ in terms of $\beta$, and applying (4.12). The growth factors $\Omega(\beta)$ are important as they dictate the stability of the basic steady state solution, and can be used to investigate the physical regimes having (qualitatively) different linear instabilities.

Closed form solutions for the implicit functions $\lambda(\beta)$ (and hence $\Omega(\beta)$) satisfying the dispersion relation $f(\lambda,\beta)=0$ cannot be determined and must instead be investigated numerically. Thus, in the work that follows, we adopt a continuation method (see Boyd (2014)) and use $\beta\in[0,\infty)$ as the continuation parameter. Starting with the value $\beta=0$ and $\lambda(0)$, we will track the implicit solutions $\lambda(\beta)$ to the dispersion relation (note that there are infinitely many) as continuous functions of $\beta$. For notational purposes we will refer to the solutions $\lambda(\beta)$ as roots to the dispersion relation. In addition, we will compute the asymptotic behavior of the functions $\lambda(\beta)$ for both small and large $\beta$. Together, the asymptotic calculations and numerical continuation will provide a comprehensive picture of the values $\Omega(\beta)$, for any given set of physical parameters $(\mathcal{B},\mathcal{S},\mathcal{Q},\mathcal{C})$. This will then enable an investigation into the different physical behaviors captured by the model.

5.1 The continuation method

We first remark on the symmetries of the dispersion relation, which will help to simplify the computation of the implicit functions $\lambda(\beta)$. Note that for any fixed $\beta$, the dispersion relation satisfies:

[TABLE]

With the above symmetries in mind, the continuation method may be restricted to the first quadrant of the $\lambda$-complex plane (equivalently, to the upper half of the $\Omega$-complex plane). Solutions $\lambda(\beta)$ may then be extended to the remaining three quadrants by symmetry.

We initialize the continuation method at $\beta=0$, for which the roots $\lambda(0)$ satisfy:

[TABLE]

and $\mathcal{R}=\mathcal{S}^{-1}(1+\mathcal{B})$ is a positive constant. The initialization value $\beta=0$ is useful as we may enumerate exactly all of the roots to equation (5.2) as follows.

First, observe that all nonzero roots to equation (5.2) are purely imaginary. To show this, it is sufficient to write $\lambda(0)=\xi+i\zeta$ for real values $\xi,\zeta$ and verify that there are no solutions to (5.2) for values $\xi>0$ and $\zeta\geq 0$ (by symmetry we may restrict to the first quadrant). Equating the real and imaginary parts of (5.2), the values $(\xi,\zeta)$ must satisfy the simultaneous equations:

[TABLE]

Note that, since $\sinh\xi+\mathcal{S}\,\xi\cosh\xi>0$ (similarly $\cosh\xi+\mathcal{S}\,\xi\sinh\xi>0$) equations (5.3)–(5.4) imply that if $\sin\zeta=0$ then $\cos\zeta=0$ (or if $\cos\zeta=0$ then $\sin\zeta=0$) — which is not possible. Hence, $\zeta$ cannot satisfy $\sin\zeta=0$, or $\cos\zeta=0$, and we are free to divide (5.3) by (5.4) to obtain the following (necessary) equation for a root:

[TABLE]

Equation (5.5), however, has no solutions for $\xi>0$ since the left-hand side is (strictly) positive and the right-hand side is non-positive. Hence, $\xi=0$, which shows that the roots to (5.2) must have the form $\lambda(0)=i\zeta$.

We can now enumerate the roots of (5.2) as $\lambda_{\pm n}(0)=\pm i\zeta_{n}$, where $0=\zeta_{0}<\zeta_{1}<\ldots$, and the values of $\zeta_{n}$ are the non-negative solutions to the equation

[TABLE]

Note that in equation (5.2), the value $\lambda_{\pm 0}(0)=0$ is a double root, and can be understood by considering $\lambda_{+0}(0)$ and $\lambda_{-0}(0)$ as two distinct roots. With this convention, writing $\lambda_{\pm n}(0)=\pm i\zeta_{n}$ then captures all roots of (5.2), including multiplicity.

The roots $\zeta_{n}$ are presented graphically in figure 2 as the intersections of the functions $g_{1}(\zeta)=\tan{\zeta}$ and $g_{2}(\zeta)=-\mathcal{R}\,\zeta$. In the limit of small $\mathcal{R}\ll 1$ (resp. large $\mathcal{R}\gg 1$), the roots $\zeta_{n}$ asymptotically approach the zeros of $\cos{\zeta_{n}}$ (resp. $\sin{\zeta_{n}}$). In the asymptotic limit $n\rightarrow\infty$, the roots $\zeta_{n}\rightarrow n\pi-\pi/2$.

We now restrict attention to the roots $\lambda_{+n}(0)$, $n\geq 0$, initialized to the upper-half plane (and for brevity drop the $+$ in the subscript of $\lambda_{+n}(\beta)$), since the remaining roots are negatives by symmetry of (5.1). With the roots $\lambda(\beta)$ of (4.24) initialized to $\lambda_{n}(0)=i\zeta_{n}$, we continuously vary $\beta\in[0,\infty)$ and track the roots $\lambda_{n}(\beta)$ as functions of $\beta$. In our subsequent linear stability analysis, the $\lambda_{0}(\beta)$ root (initialized to $\lambda_{0}(0)=0$) will play a particularly important role. As a result, we will refer to $\lambda_{0}(\beta)$ as the film root, and (from now on) write $\lambda_{\rm f}(\beta)$. The phrase “film root” is motivated by the fact that the corresponding complex frequency $\Omega_{\rm f}(\beta)=\lambda_{\rm f}(\beta)^{2}-\beta^{2}$ is analogous to the frequencies $\Omega(\beta)$ given by a free thin film equation (see for instance §A). The remaining roots are initialized to $\lambda_{n}(0)=i\zeta_{n}$ for $n=1,2,\ldots$.

As a technical point, we stress that the roots $\lambda_{\rm f}(\beta)$ or $\lambda_{n}(\beta)$ are only (uniquely) identifiable by their initial values $i\zeta_{n}$ within an interval $0\leq\beta\leq\beta_{\rm coll}$ for which no collision has occurred. Once a collision occurs, i.e. two (or more) roots collide at a value $\beta=\beta_{\rm coll}$, it is generally not possible to identify uniquely two (or more) post-collision roots $\lambda(\beta)$ at values $\beta>\beta_{\rm coll}$ with their initial values $i\zeta_{n}$.

As a final remark on the numerical computations, we follow a standard continuation approach: at each step, using the value $\lambda(\beta)$ as an initial guess, we use Newton’s method to compute $\lambda(\beta+\Delta\beta)$, where $\Delta\beta$ is the increment (chosen adaptively to ensure convergence at each step). As a practical detail, to enable the method to find complex valued solutions, we initialize the Newton algorithm with a value that does not lie strictly on either the real or complex axis by perturbing the initial guess via $\lambda(\beta)+\epsilon(1+i)$, for $\epsilon\ll 1$. This is to avoid having Newton iterates become trapped to the (invariant) real or complex axis.

5.2 Asymptotic behavior of the roots for small $\beta$

The previous section demonstrates that the values $\lambda_{n}(0)=i\zeta_{n}$ are purely imaginary. Hence, at $\beta=0$, the growth rates $\Omega$ lie along the negative real axis: $\Omega_{\rm f}(0)=0$ (corresponding to the film root $\lambda_{\rm f}(0)=0$), and $\Omega_{n}(0)=-\zeta_{n}^{2}<0$ (for the roots $\lambda_{n}(0)=i\zeta_{n}$, $n\neq 0$); see (4.12). The purpose of this section is to examine how the values $\Omega_{n}(\beta)$ and $\Omega_{\rm f}(\beta)$ change for small values of $0\leq\beta\ll 1$.

We first compute the small $\beta$ behavior of the film root $\lambda_{\rm f}(\beta)$ and corresponding growth rate $\Omega_{\rm f}(\beta)$. This can be done by expanding (4.24) in powers of $\lambda$ (about $\lambda_{\rm f}(0)=0$) to obtain:

[TABLE]

where, recall, $a_{1},a_{2},b_{1},b_{2}$ depend on $\beta$ via (4.21). Truncating the expansion (5.7) at $\mathcal{O}(\lambda^{4})$, and setting $f(\lambda,\beta)$ to zero, yields an approximate solution for $\lambda_{\rm f}(\beta)$, valid at small $\beta$:

[TABLE]

Here, the Taylor coefficients $\eta_{1}$ and $\eta_{3}$ are

[TABLE]

respectively, where $\eta_{3}$ has been written compactly using the definition of $\eta_{1}$. The above calculation shows that the double root at $\lambda_{\rm f}(0)=0$ splits immediately into two real nonzero roots given (approximately) by (5.8). As a convention, we use $\lambda_{\rm f}(\beta)$ to denote the positive branch in (5.8). Via (4.12) we then have the small-wavenumber expansion of the corresponding growth rate $\Omega_{\rm f}(\beta)$, viz.,

[TABLE]

Here, the first term inside the square brackets describes capillary stabilization, and the other terms including $\bar{\mathit{\Theta}}_{\rm H}$ pertain to thermocapillary effects. Notably, this expression predicts that thermocapillarity acts both to destabilize small wavenumbers and stabilize large ones, in contrast to the strict thermocapillary destabilization observed for thin substrates (see Appendix A). In particular, (5.10) shows that the relative importance of thermocapillary stabilization will increase for large $\mathcal{S}$ (substrate thickness), large $\mathcal{C}$ (imposed temperature difference), or small $\mathcal{Q}$ (diffusive effects).

In figure 3, the small-$\beta$ approximation (5.10) for $\Omega_{\rm f}(\beta)$ is compared to exact root branches calculated via numerical continuation. In each panel, one comparison is made to demonstrate a parameter set for which the agreement is qualitatively good; both the maximum growth rates and the cutoff wavenumbers are adequately predicted by (5.10). The surprisingly good prediction for the cutoff wavenumber is due to the original expansion in (5.7) being based on small $\lambda$: since $\lambda^{2}=\Omega+\beta^{2}$, we expect the resulting approximation (5.10) to be good near the origin $\beta=\Omega=0$ and at the cutoff where $\beta^{2}\ll 1$ still and $\Omega=0$. In the region between $\beta=0$ and the $O(1)$ cutoff values of the growthrate $\Omega$, the small-$\lambda$ assumption made to obtain (5.7) is certainly violated. This is reflected in the poor agreement observed in this region between the asymptotic approximation (5.10) and the numerically-calculated curves in figure 3.

Each branch in figure 3 is presented from $\beta=0$ to a critical value of $\beta_{\rm coll}$, beyond which the real film root ceases to exist ($\beta_{\rm coll}$ will be defined more precisely below). Before focusing on these critical points in the next section, we note that several key physical behaviors can be inferred from figure 3. First, we see in panel (a) that larger imposed temperature differences (larger $\mathcal{C}$) increase both the unstable band of wavenumbers and their associated growth rates. Panel (b) then shows that diffusive effects, as measured by $\mathcal{Q}$, suppress the growth rates of instability without significantly changing the bandwidth of unstable wavenumbers. The loss of agreement in panel (c) occurs because, thinner films, via small values of $\mathcal{B}$, promote the relative importance of the substrate thermal process.

Panel (d) shows the dependence on $\Omega_{\rm f}(\beta)$ on the parameter $\mathcal{S}$. We note that the increase in the growth rates with $\mathcal{S}$ shown in panel (d) results from having scaled time with respect to $d^{2}$, see (3.2). The (approximately linear) increase of $\Omega_{\rm f}$ with $\mathcal{S}$ shown here in fact corresponds to a linear decrease in the dimensional growth rates–due primarily to the increased thermal resistance of thicker substrates.

5.3 The roots for $n\neq 0$

To examine the roots $\lambda_{n}(\beta)$ and $\Omega_{n}(\beta)$ with $n\neq 0$ in the small-$\beta$ limit, we substitute the power series,

[TABLE]

into the dispersion relation with $\zeta_{n}$ as given by (5.6) (see also figure 2) and unknown coefficients $\sigma_{n}$ (which will turn out to be real). Such an expansion in even powers of $\beta$ is justified by the symmetry relations (5.1) and by the fact that the values $\lambda_{n}(\beta)$ are simple near $\beta=0$. Note that this is in contrast to the roots $\pm\lambda_{\rm f}(\beta)$ that emerge from the double root $\lambda_{\rm f}(0)$, which does not have an even power series at $\beta=0$ (double roots generally split via a square root dependence on the continuation parameter).

Expanding $f(\lambda_{n}(\beta),\beta)=0$ (see (4.24)) for small $\beta$, and setting the $O(\beta^{2})$ term to zero yields

[TABLE]

In this expression, $\tan{\zeta_{n}}$ has been replaced with $-\mathcal{R}\,\zeta_{n}$ per (5.6) (recall that $\mathcal{R}=\mathcal{S}^{-1}(1+\mathcal{B})$). Substituting the expressions (5.11)–(5.12) into (4.12) yields

[TABLE]

At zero wavenumber, all these roots are real and negative: $\Omega_{n}(0)=-\zeta_{n}^{2}$. As $\beta$ increases from zero, if $2\,\zeta_{n}\,\sigma_{n}>1$, then the roots (initially) move along the negative real axis in the complex plane towards the right-half plane (RHP) (we will also use LHP to denote the left-half plane). Alternatively if $2\,\zeta_{n}\,\sigma_{n}<1$, then the roots move to the left along the negative real axis (see figure 4).

We now remark that the product $\zeta_{n}\,\sigma_{n}>0$ is (i) always positive, (ii) monotonically decreases with increasing values of $n$ (due to the fact that the values $\zeta_{n}$ monotonically increase with $n$), and (iii) $\lim_{n\rightarrow\infty}\zeta_{n}\,\sigma_{n}\rightarrow 0$. As a result, only a finite number of the values $\Omega_{n}(\beta)$ (closest to the origin) can have a positive $O(\beta^{2})$ coefficient, and hence initially move towards the unstable RHP. All other roots move (at small $\beta$) farther into the LHP. Figure 4 demonstrates, via the continuation method, the motion of the roots $\{\Omega_{\rm f}(\beta),\Omega_{n}(\beta)\}$, with varying $\beta$. The plot shows the film root ($\Omega_{\rm f}(0)=0$) and $\Omega_{n}(0)=-\zeta_{n}^{2}$ roots at $\beta=0$. The arrows for the markers denote the numerically computed directions in which the roots move as $\beta>0$ increases. The arrows also coincide with the small-$\beta$ approximations (5.10) and (5.13). Note that the figure also demonstrates that the roots become complex only after a collision, and generally move to the left (becoming more stable) with increasing $\beta$.

6 Oscillatory instability classification

In this section, we investigate how the frequencies $\Omega_{\rm f}(\beta)$ and $\Omega_{n}(\beta)$ move in the complex plane (as functions of the continuation parameter $\beta$), and lead to oscillating in time solutions of the linearized equations (4.13)–(4.17) with exponentially growing amplitudes. A value $\Omega$ (that satisfies the dispersion relation) is oscillatory unstable if $\Omega$ lies in the strict RHP, and does not lie along the real axis. That is, $\Omega$ satisfies:

[TABLE]

Combined with the symmetry observations from the previous section (§5), conditions (6.1) place restrictions on how exactly a root $\Omega(\beta)$ can become oscillatory unstable as the wavenumber $\beta$ increases from [math].

First, condition (ii) in (6.1) requires that (as $\beta$ varies) two frequencies $\Omega(\beta)$ must collide at some value of $\beta$ — that is, there is a value of $\beta$ for which two of the frequencies $\{\Omega_{\rm f}(\beta),\Omega_{n}(\beta)\}$ are equal. This is because the complex frequencies $\{\Omega_{\rm f}(0),\Omega_{n}(0)\}$ are simple (at $\beta=0$), move continuously with $\beta$, and cannot leave the real axis as long as they remain simple (due to conjugate symmetry of $\Omega(\beta)$, see (5.1)). Hence, a necessary condition for (ii) is a collision (double frequency) at some $\beta$.

Second, the results from section 5 show that the frequencies $\{\Omega_{\rm f}(0),\Omega_{n}(0)\}$ are on the non-positive real axis, and only a finite number of the largest roots initially move towards the RHP (as $\beta$ increases). Hence, the two largest roots $\Omega_{\rm f}(\beta)$ (film root) and $\Omega_{1}(\beta)$ are the most likely candidates to collide and satisfy (6.1). In other words, oscillatory instabilities most likely arise from the motion and collisions of $\Omega_{\rm f}(\beta)$ and $\Omega_{1}(\beta)$.

Treating the continuation parameter $\beta$ as a bifurcation parameter (holding the parameters $(\mathcal{C},\mathcal{Q},\mathcal{B},\mathcal{S})$ fixed), we now classify how exactly $\Omega_{\rm f}(\beta)$ and $\Omega_{1}(\beta)$ collide (bifurcate) and satisfy (6.1) (lead to oscillatory instabilities). We define the (first) collision of $\{\Omega_{\rm f}(\beta),\Omega_{1}(\beta)\}$ to occur at the value $\beta_{\rm coll}$, and denote

[TABLE]

Note that the values $(\beta_{\rm coll},\Omega_{\rm coll})$ are readily identifiable: in addition to the dispersion relation $\Omega_{\rm coll}=\lambda_{\rm coll}^{2}-\beta_{\rm coll}^{2}$, with $f(\lambda_{\rm coll},\beta_{\rm coll})=0$, the values $(\beta_{\rm coll},\Omega_{\rm coll})$ satisfy the condition of a double root required by the implicit function theorem: $\partial_{\lambda}f(\lambda_{\rm coll},\beta_{\rm coll})=0$. For values of $\beta\in[0,\beta_{\rm coll}]$, the frequencies $\Omega_{\rm f}(\beta),\Omega_{1}(\beta)\in\mathbb{R}$; however, for $\beta>\beta_{\rm coll}$ (after the collision), the roots appear as complex conjugate pairs that we denote as $\Omega_{\rm f,1}(\beta)=\Omega_{\rm r}(\beta)+i\Omega_{\rm i}(\beta)$ and $\Omega_{\rm f,1}^{*}(\beta)$, with $\Omega_{\rm r}(\beta),\Omega_{\rm i}(\beta)\in\mathbb{R}$. We further denote $\beta_{\rm im}>\beta_{\rm coll}$ (if it exists) as the first value at which the frequency $\Omega_{\rm f,1}(\beta)$ crosses the imaginary axis, i.e. $\Omega_{\rm r}(\beta_{\rm im})=0$. Summarizing the notation, we have:

[TABLE]

We now identify two characteristic ways for the roots $\Omega_{\rm f}(\beta),\Omega_{1}(\beta)$ to satisfy conditions (6.1) and give rise to oscillatory instabilities.

Type I: $\Omega_{\rm coll}>0$. The top-left subfigure 5 highlights how the frequencies $\Omega_{\rm f}(\beta),\Omega_{1}(\beta)$ move from locations $a$–$c$ at wavenumbers $\beta<\beta_{\rm coll}$; collide at location $d$ ($\beta=\beta_{\rm coll}$ with $\Omega_{\rm coll}>0$); are unstable and satisfy (6.1) at any point $e$ between $d$ and $f$; cross over into the LHP at $f$; and are stable at points $g$ in the LHP. The behavior of the perturbations (4.1)qualitatively changes at different wavenumbers $\beta$ through a series of bifurcations. To provide a visual characterization of the linearized dynamics of the perturbations (4.1) at different points $a$–$g$, we plot the phase plane trajectories in the eigenmodes with frequencies $\Omega_{\rm f}(\beta),\Omega_{1}(\beta)$ (when $0\leq\beta\leq\beta_{\rm coll}$) or frequencies $\Omega_{\rm f,1}(\beta),\Omega_{\rm f,1}^{*}(\beta)$ (for $\beta>\beta_{\rm coll}$). For wavenumbers $\beta\leq\beta_{\rm coll}$, a perturbation in $H(X,T)$ with amplitude $\delta\hat{H}_{\rm f}$ in frequency $\Omega_{\rm f}(\beta)$, and amplitude $\delta\hat{H}_{1}$ in frequency $\Omega_{1}(\beta)$, evolves as:

[TABLE]

Figure 5$(a-d)$ plots the phase plane dynamics in the $h_{1}(T)$–$h_{\rm f}(T)$ plane. For wavenumbers $\beta>\beta_{\rm coll}$, the frequencies are complex $\Omega_{\rm f,1}(\beta)$,$\Omega_{\rm f,1}(\beta)^{*}$ and we write the perturbation as

[TABLE]

where $\hat{H}_{\rm R}$ and $\hat{H}_{\rm I}$ are the two amplitudes of the perturbation. The panels $e$–$g$ in figure 5) plot the phase plane trajectories in the $h(T)$–$dh(T)$ plane. Here we use the short-form notation $dh(T)$ for the second linearly independent term in (6.5) since it is proportional to $\frac{dh}{dt}$. The subfigures show the qualitatively different phase plane behavior, and emergence of oscillatory instabilities, as the bifurcation parameter $\beta$ varies. The boundary in the parameter space $(\mathcal{C},\mathcal{Q},\mathcal{B},\mathcal{S})$, for Type I behavior to occur, must satisfy (as a necessary condition) $\Omega_{\rm coll}=0$.

Type II: $\Omega_{\rm coll}\leq 0$ and $\Omega_{\rm r}(\beta)>0$ for some $\beta>\beta_{\rm coll}$. Figure 6 highlights (via arrows) the motion of the roots $\Omega_{\rm f}(\beta)$ (red), $\Omega_{1}(\beta)$ (black), and post collision roots $\Omega_{\rm f,1}(\beta)$ (black with white dashed line) with varying $\beta$. The top-left panel in Figure 7 shows again the motion of the frequencies $\Omega_{\rm f}(\beta)$, $\Omega_{1}(\beta)$, and $\Omega_{\rm f,1}(\beta)$. The frequencies $\Omega_{\rm f}(\beta)$, $\Omega_{1}(\beta)$ move from locations $a$–$c$ and collide at $d$ with $\Omega_{\rm coll}\leq 0$. The values $\Omega_{\rm f,1}(\beta)$ are complex at location $e$, and then move into the right-half plane $\Omega_{\rm r}(\beta)\geq 0$ at $f$–$g$; thereby satisfying conditions (6.1). The subfigures also show the phase plane trajectories for perturbations to the film height (in a fashion completely analogous to figure 5) given by equations (6.4)–(6.5). The boundary in the parameter space $(\mathcal{C},\mathcal{Q},\mathcal{B},\mathcal{S})$, for Type II behavior, requires (as a necessary condition) that $\max_{\beta\geq\beta_{\rm coll}}\Omega_{\rm r}(\beta)=0$.

Having criteria for the boundaries in the parameter space $(\mathcal{C},\mathcal{Q},\mathcal{B},\mathcal{S})$ of Types I or II oscillatory instabilities will become useful in the following section. Specifically, we will use these conditions to help plot phase diagrams and identify model parameters and experimental conditions that may yield oscillatory instabilities.

The underlying distinction between Types I and II instabilities occurs from viewing $\beta$ as a bifurcation parameter. Each root $\Omega$ for any value of $\beta$ is associated with a linear dynamical system for the variables in equations (4.1). Type I oscillatory instabilities occur from one bifurcation when the roots collide at $\beta=\beta_{\rm coll}$. Meanwhile, Type II oscillatory instabilities occur through two bifurcations: the first at $\beta=\beta_{\rm coll}$ when the roots collide, and the second at the value $\beta=\beta_{\rm im}$, when the roots enter into the RHP (e.g., see, Strogatz (2015), chapter 8). Categorizing oscillatory instabilities as Type I or Type II through different bifurcations provides criteria that we will use in §7 to systematically determine which parameter values $(\mathcal{C},\mathcal{Q},\mathcal{B},\mathcal{S})$ give rise to oscillatory instabilities. Type I and Type II instabilities are physically distinguishable through their bands of unstable wavenumbers: Type I does not contain a band of wavenumbers with values $\beta<\beta_{\rm coll}$ where $\Omega_{\rm f}(\beta)$ is stable (figure 5 at points $a$, $b$, $c$ and $d$ are unstable); in contrast, Type II does contain an interval of wavenumbers with values $\beta<\beta_{\rm coll}$ where $\Omega_{\rm f}(\beta)$ is stable (figure 7 at points $c$ and $d$ are unstable).

In addition to classifying Types I and II oscillatory instabilities, we further distinguish whether a set of parameters $(\mathcal{C},\mathcal{Q},\mathcal{B},\mathcal{S})$ gives rise to global oscillatory instabilities. A set of parameters $(\mathcal{C},\mathcal{Q},\mathcal{B},\mathcal{S})$ is said to be globally oscillatory unstable if there is a $\beta_{\rm c,max}$ and $\Omega_{\rm c,max}=\Omega(\beta_{\rm c,max})$ satisfying the dispersion relation (4.23) that is oscillatory unstable, and has the largest growth rate $\textrm{Re}(\Omega_{\rm c,max})\geq\textrm{Re}(\Omega(\beta))$ for all $\beta$, $\Omega(\beta)$ (satisfying the dispersion relation (4.23)). Global oscillatory unstable parameter values are physically significant because they represent experimental situations where the most unstable perturbation to the linearized system (4.1) is oscillatory unstable (and hence the most likely to be observed).

The bottom panel of figure 4 provides additional information for the motion (as functions of $\beta$) of the frequencies $\Omega_{\rm f}(\beta),\Omega_{n}(\beta)$. This highlights the fact that the complete motion of the roots is quite complicated. In the figure, bold lines show the trajectories of the roots in the complex plane over an interval $0\leq\beta\leq\beta_{*}$ ($\beta_{*}>0$ is chosen somewhat arbitrarily). The bold curves trace out several collisions of the roots (including collisions made by roots $\Omega_{n}(\beta)$ for $n>1$), and highlight that they always appear in complex conjugate pairs. The arrows and red circles show roots at the value $\beta=\beta_{*}$, with the arrows indicating the direction of motion for the roots at $\beta=\beta_{*}$. Although several of the roots depicted at $\beta=\beta_{*}$ are moving in the positive direction along the real axis, we note that each of these roots ultimately reverses direction and moves into the left-half plane for sufficiently large $\beta$.

In the numerical continuation of the roots $\{\Omega_{\rm f}(\beta),\Omega_{n}(\beta)\}$, we have always observed that oscillatory instabilities arise as Type I or Type II, as described in this section. It may be possible that oscillatory instabilities occur from the collision of other roots (for instance, a collision including $\Omega_{2}(\beta)$); however we did not observe this in any of our investigations. If the only possible mechanism to obtain oscillatory instabilities is through Type I or Type II, and the largest growth rate $\Omega_{\rm c,max}$ occurs at a value $\Omega_{\rm f}(\beta)$ (which we also numerically observe to be the case in our studies), then we may simplify the condition for global oscillatory instabilities to growth rates computed in terms of $\Omega_{\rm f}(\beta)$ and $\Omega_{\rm f,1}(\beta)$ by defining:

[TABLE]

The condition for global oscillatory instabilities is then

[TABLE]

To characterize oscillatory instabilities, we will also make use of the most unstable wavenumber, $\beta_{\rm max}$ defined as:

[TABLE]

or alternatively written as the argument of the maximum $\beta_{\rm max}=\textrm{argmax}_{\beta\geq\beta_{\rm coll}}\Omega_{\rm r}(\beta)$. We will also denote the imaginary frequency of the most unstable wave number as $\Omega_{\rm i,max}=\Omega_{\rm i}(\beta_{\rm max})$.

In practice, the maximization $\max_{\beta\geq\beta_{\rm coll}}\Omega_{\rm r}(\beta)$ in ondition (6.7) implies that we maximize the real value of the root $\Omega_{\rm f,1}(\beta)$ over a suitably large range of $\beta$ values (by taking $\beta$ large enough $\Omega_{\rm r}(\beta)$ will eventually become negative). Replacing the inequality $(\leq)$ in (6.7) with an equality $(=)$ then provides a condition for the boundary of the global oscillatory instability region.

7 Emergence of oscillatory instabilities

The purpose of this section is to explore the model parameter space $(\mathcal{C},\mathcal{Q},\mathcal{B},\mathcal{S})$ and characterize which parameter regions give rise to oscillatory instabilities. This will provide a guide for experimental scenarios in which one may likely see oscillatory instabilities. To compute these regions we use the formulas for the boundaries of the parameter regions, satisfied by Type I and Type II instabilities, developed in §6. As a general guide, we also introduce a heuristic value

[TABLE]

as the rate of change of the real value of the roots $\Omega_{\rm f,1}(\beta)=\Omega_{\rm r}(\beta)+i\Omega_{\rm i}(\beta)$ immediately after the collision $\beta\rightarrow\beta_{\rm coll}^{+}$. A positive value $\partial_{\beta}\Omega_{\rm r}>0$ (resp. $<$) implies the roots $\Omega_{\rm f,1}(\beta)$ move towards the right (resp. left) in the complex plane as $\beta$ increases past $\beta_{\rm coll}$. Knowing whether the roots $\Omega_{\rm f,1}(\beta)$ move towards the left (more stable, viz. figures 4 and 5) or right (more unstable, viz. figures 6 and 7) in the complex plane after the collision is a useful heuristic when identifying regions of global oscillatory instabilities. Specifically, numerical evidence shows that parameter values $(\mathcal{C},\mathcal{Q},\mathcal{B},\mathcal{S})$ that have $\partial_{\beta}\Omega_{\rm r}<0$ (solutions initially become more stable after collisions) do not exhibit global oscillatory instabilities.

This section is organized as follows: in §7.1 we plot and detail the behavior of the phase diagram for oscillatory instabilities, with a focus on the parameters $(\mathcal{C},\mathcal{B},\mathcal{S})$. In §7.2 we examine the effect of the parameter $\mathcal{Q}$ (material property dependent) on the behavior of the phase diagram. The results from §7.2 will help guide realistic choices of material properties and experimental conditions for observing oscillatory instabilities. Guided by the results in §7.1–§7.2, in §7.3 we discuss materials that give rise to reasonable model parameter values for observing oscillatory instabilities.

7.1 Material phase diagrams and oscillatory instabilities

In this section we plot phase diagrams that show for which model parameters oscillatory instabilities occur. Our approach for plotting the diagrams is motivated by experimental considerations. The parameter $\mathcal{Q}$ depends on the material properties, and is the most difficult to change in experiments (it requires changing the substrate or fluid materials in the experiment). The values of $\mathcal{B},\mathcal{S}$ can by varied by modifying the thickness of the film ($\mathcal{B}$) and substrate ($\mathcal{S}$), while $\mathcal{C}$ may be varied easily by modifying the temperature difference across the film and substrate. Since we have four parameters $(\mathcal{C},\mathcal{Q},\mathcal{B},\mathcal{S})$ we adopt the following approach to visualize the phase diagrams: we fix a value of $\mathcal{Q}$, and then plot the phase diagram in the $\mathcal{B}\,$–$\,\mathcal{S}$ plane for different values of $\mathcal{C}$. This is equivalent to plotting cross-sections of the three dimensional phase diagram $(\mathcal{C},\mathcal{B},\mathcal{S})$ (holding $\mathcal{Q}$ constant). For the purpose of developing better intuition, one may think of $\mathcal{B}$ and $\mathcal{S}$ being the thicknesses of the film and substrate (respectively) and $\mathcal{C}$ the imposed temperature difference.

Figure 8 plots the $\mathcal{B}\,$–$\,\mathcal{S}$ phase diagram for values of $0.425\leq\mathcal{C}\leq 0.57$, holding $\mathcal{Q}=5.114\times{10}^{-4}$ fixed. The subfigures in 8 reveal an onset parameter value $\mathcal{C}_{\rm o}=0.5$ such that for $\mathcal{C}<\mathcal{C}_{\rm o}$ there is no region of oscillatory instability in the $\mathcal{B}\,$–$\,\mathcal{S}$ phase diagram; while values $\mathcal{C}>\mathcal{C}_{\rm o}$ give one connected region (shown in color) of oscillatory instability ($\mathcal{Q}$ is chosen to three decimals so that $\mathcal{C}_{\rm o}$ is a single decimal). At $\mathcal{C}=0.5$, the region of oscillatory instability emerges from a single point $(\mathcal{B}_{\rm o},\mathcal{S}_{\rm o})=(1.12,3.21)$ (labeled D in the figure). The different shadings in figure 8 provide details on how the roots $\Omega_{\rm f}(\beta),\Omega_{1}(\beta)$ become unstable (i.e. Type I or II), as well as the sign of the heuristic quantity $\partial_{\beta}\Omega_{\rm coll}$. The region inside the dashed curves indicates where $\partial_{\beta}\Omega_{\rm coll}>0$. When investigating the four dimensional parameter space $(\mathcal{C},\mathcal{Q},\mathcal{B},\mathcal{S})$, the heuristic $\partial_{\beta}\Omega_{\rm coll}>0$ is helpful in identifying regions that have global oscillatory instabilities, as they numerically appear inside the heuristic, see e.g. figure 10 (note that figure 8 shows no regions of global oscillatory instability). The boundary of the heuristic is easy to compute and can then be used to restrict the region where a refined search for global oscillatory instabilities can be done. In addition to $\mathcal{C}_{\rm o}$, we introduce $\mathcal{C}_{\rm g}$ as the critical parameter value for which global oscillatory instabilities occur, i.e. global oscillatory instabilities occur when $\mathcal{C}>\mathcal{C}_{\rm g}$, while for values of $\mathcal{C}<\mathcal{C}_{\rm g}$ all oscillatory instabilities are non-global (such as those in figure 8). A key observation from the figure is that oscillatory instabilities do not occur at low temperatures ($\mathcal{C}<\mathcal{C}_{\rm o}$).

The qualitative differences in how the roots become unstable (as $\beta$ varies) are shown in figure 9. Each part of figure 9 plots $\Omega_{\rm f}(\beta)$, $\Omega_{1}(\beta)$, and the real and imaginary values of $\Omega_{\rm f,1}(\beta)$ (defined in (6.3)) for parameter values that capture the behavior at locations A–K in figure 8. Specifically, in figure 9, panels A–C first show the change in sign of the heuristic quantity $\partial_{\beta}\Omega_{\rm coll}$; panels E, G, and J correspond to Type I instabilities; panels H and K to Type II instabilities. Meanwhile, panels D, F, H and I provide a comprehensive survey of parameter values that lie on the boundaries of Type I or II instabilities. Type I oscillatory instabilities have one continuous band of unstable wavenumbers $0\leq\beta<\beta_{\rm im}$, of which $[0,\beta_{\rm coll}]$ is monotonically unstable and $(\beta_{\rm coll},\beta_{\rm im})$ is oscillatory unstable. Type II oscillatory instabilities have two continuous bands of unstable wavenumbers, separated by a gap of stable wavenumbers that includes the interval $(\beta_{\rm coll},\beta_{\rm im}]$. Together, the panels in figure 9 characterize all the different possibilities for (possibly oscillatory) instability development. Note that figure 9 does not admit global oscillatory instabilities: all the subfigures are plotted for values of $\mathcal{C}<\mathcal{C}_{\rm g}$, where $\mathcal{C}_{\rm g}$ is the critical parameter such that global oscillatory instability can only occur for $\mathcal{C}>\mathcal{C}_{\rm g}$ (oscillatory instabilities for $\mathcal{C}<\mathcal{C}_{\rm g}$ are non-global).

Figure 10 continues the phase diagrams in figure 8 to the values $1.5635\leq\mathcal{C}\leq 5$ (again with $\mathcal{Q}{}=5.114\times{10}^{-4}$). The value of $\mathcal{C}_{\rm g}=1.5635$ in the first panel of figure 10 is significant: it is the onset value for global oscillatory instabilities, which emerge at the point labeled L. For $\mathcal{C}>\mathcal{C}_{\rm g}$, there are regions of parameter values $(\mathcal{B},\mathcal{S})$ (red shading, figure 10) that are globally oscillatory unstable. The value $\mathcal{C}=5$ in figure 10 highlights that global oscillatory unstable modes can be either Type I or Type II. Figure 11 plots the roots $\Omega_{\rm f}(\beta),\Omega_{1}(\beta)$ and $\Omega_{\rm f,1}(\beta)$ for behavior indicative of parameter values L–Q in figure 10. Subfigures 11 L, N and O correspond to parameter values on the boundary of the global oscillatory instability region (L being Type I, O being Type II and N being on the boundary of Type I). The remaining panels M, P and Q show points for parameter values in the interior of the global oscillatory instability region.

Lastly, we characterize the phase diagrams as $\mathcal{C}$ (which one may think of as the temperature difference) becomes large. Figure 12 fixes $\mathcal{Q}=0.006$ and plots contours (corresponding to the curve $\Omega_{\rm coll}=0$) in subpanels (a), (d) that enclose Type I instabilities; contours in subpanels (b), (e) that define the heuristic (corresponding to the curve $\partial_{\beta}\Omega_{\rm coll}=0$); and contours in subpanels (c), (f) that enclose regions of global oscillatory instabilities (corresponding to the curve $\Omega_{\rm f,max}=\Omega_{\rm r,max}$). In the subpanels the contours are enumerated 1–8 for convenience, and correspond to increasing values of $\mathcal{C}$, with the smallest value of $\mathcal{C}$ labelled $1$ and the largest labelled $8$ (the numerical values of $\mathcal{C}$ are stated in the caption). We define $\mathcal{C}_{\rm s}$ to be the (smallest) onset value of $\mathcal{C}$ for which the zero contour of the heuristic $\partial_{\beta}\Omega_{\rm coll}$ emerges in the $\mathcal{B}\,$–$\,\mathcal{S}$ plane; and $(\mathcal{B}_{\rm s},\mathcal{S}_{\rm s})$ as the point at which the heuristic $\partial_{\beta}\Omega_{\rm coll}=0$ first emerges. The bottom panels (d)–(f) of figure 12 show the onset coordinate values $(\mathcal{B}_{\rm o},\mathcal{S}_{\rm o})$ for oscillatory instabilities (subpanel (d)); $(\mathcal{B}_{\rm o},\mathcal{S}_{\rm o})$ for the heuristic (subpanel (e)); and $(\mathcal{B}_{\rm g},\mathcal{S}_{\rm g})$ for global oscillatory instabilities (subpanel (f)).

Subfigures (a), (d) show that the contours 1–8 are nested: given two contours $1\leq j,k\leq 8$ with $\mathcal{C}_{j}>\mathcal{C}_{k}$, then contour $\mathcal{C}_{j}$ encloses contour $\mathcal{C}_{k}$. This fact implies that the regions of oscillatory instability in the phase diagrams $\mathcal{B}\,$–$\,\mathcal{S}$ become large as $\mathcal{C}$ increases. The regions of global oscillatory instabilities in subfigures (c), (f) are not nested, but still grow in size as $\mathcal{C}$ increases. How the oscillatory unstable regions change as $\mathcal{C}$ varies is significant. Generally speaking, contours 5–8 (subpanel (a)) show that these contours expand in all directions with increasing $\mathcal{C}$, and eventually cover the entire $\mathcal{B}\,$–$\,\mathcal{S}$ plane. This suggests that (for this value of $\mathcal{Q}$), any pair $(\mathcal{B},\mathcal{S})$ will lead to oscillatory instabilities for a sufficiently large $\mathcal{C}$ (temperature difference) value. In contrast, the existence of global oscillatory unstable regions (subpanels (c), (f)) in the $\mathcal{B}\,$–$\,\mathcal{S}$ plane depends on $\mathcal{C}$:

• •

There are values of $(\mathcal{B},\mathcal{S})$ (for instance, if the ratio $\mathcal{S}/\mathcal{B}$ is sufficiently small) that will not be globally oscillatory unstable for any value of $\mathcal{C}$.

• •

An arbitrary set of parameters $(\mathcal{B},\mathcal{S})$ (material thicknesses) will only be subject to a global oscillatory instability for a range of $\mathcal{C}$ (temperature differences) values (between a minimum value $\mathcal{C}_{\rm g}$ and some maximum value).

Similar conclusions can be drawn that describe the change in the zero contours of the heuristic $\partial_{\beta}\Omega_{\rm coll}$ with $\mathcal{C}$ (subpanels (b), (e)). These observations provide a valuable guide for choosing materials that yield realistic experimental setups for observing Type II and global oscillatory instabilities.

7.2 The effect of $\mathcal{Q}$ on the behavior of the phase diagrams

As previously stated, the value of $\mathcal{Q}$ depends on the material properties in the experiment and cannot be modified by the experimental setup (i.e. experimental geometry or temperature). Therefore, identifying values of $\mathcal{Q}$ that give rise to oscillatory instabilities (and in particular, global oscillatory instabilities) will guide the choice of experimental materials.

Numerical experiments show that modifying the value of $\mathcal{Q}$ does not change the qualitative behavior of the phase diagrams in §7.1. Changes in the value of $\mathcal{Q}$ can, however, result in a (potentially significant) quantitative change in the onset values $\mathcal{C}_{\rm o}$, $\mathcal{C}_{\rm g}$ (temperature difference), as well as the locations in the $\mathcal{B}\,$–$\,\mathcal{S}$ plane (i.e. thicknesses of the film and substrate) for which the regions of oscillatory and global oscillatory instabilities emerge: $(\mathcal{B}_{\rm o},\mathcal{S}_{\rm o})$ and $(\mathcal{B}_{\rm g},\mathcal{S}_{\rm g})$.

Figure 13 plots the onset values $\mathcal{C}_{\rm o},\mathcal{B}_{\rm o},\mathcal{S}_{\rm o}$ (solid line) for oscillatory instabilities, as well as the onset values $\mathcal{C}_{\rm g},\mathcal{B}_{\rm g},\mathcal{S}_{\rm g}$ (dashed line) for global oscillatory instabilities, as functions of $\mathcal{Q}$. Onset values $\mathcal{C}_{\rm s},\mathcal{B}_{\rm s},\mathcal{S}_{\rm s}$ of the emergence of the zero of the heuristic (dotted lines) are also given. Figure 13 restricts the range of $0\leq\mathcal{Q}\leq 0.1$ to an experimentally feasible range. The plots are significant since $\mathcal{C}_{\rm o}$ and $\mathcal{C}_{\rm g}$ are the minimum values for which oscillatory and global oscillatory instabilities occur. In addition, the values of $\mathcal{B}_{\rm g}$ and $\mathcal{S}_{\rm g}$ provide information on choosing $\mathcal{B}$ and $\mathcal{S}$. Guided by the qualitative behavior in figure 8, the diagrams show that choosing $(\mathcal{B},\mathcal{S})$ close to $(\mathcal{B}_{\rm g},\mathcal{S}_{\rm g})$ will likely yield global oscillatory instabilities for some range of $\mathcal{C}>\mathcal{C}_{\rm g}$.

As a computational remark, the values $\mathcal{C}_{\rm g}$ and $\mathcal{C}_{\rm o}$ are calculated by minimizing the value of $\mathcal{C}$ in the region of $(\mathcal{C},\mathcal{B},\mathcal{S})$ parameter space that satisfies the Type I instability criterion (i.e. satisfy the condition $\Omega_{\rm coll}=0$) or the global oscillatory instability criterion (condition (6.7)).

7.3 Experimental considerations for oscillatory instabilities

In this section, we discuss oscillatory instabilities in the context of an experimental setting. This will shed light on the physical mechanism for oscillatory instabilities. In particular, we will contrast two cases: a substrate that is a conductor (having a large thermal conductivity); and a substrate that is an insulator (having a small thermal conductivity). We conclude that oscillatory instabilities are far more likely to be observed for films heated by substrates that are insulators.

To draw this conclusion, we estimate the rate of heat transfer between the surface and the gas $q=5$ W/m$\cdot{}$K (see table 1) and consider low-viscosity silicone oil films for all the cases presented in this section. Physical properties of the silicone oil, copper (a conductor), and PMMA (poly(methyl methacrylate), an insulator) are given in Table 1. We first consider films of silicone oil heated by copper substrates. The results from the previous section will demonstrate that oscillatory instabilities for silicone oil-copper systems are not likely in experimentally feasible conditions.

Together, the silicone oil–copper system yields a parameter value $\mathcal{Q}=2.88\times{10}^{3}$, which is four orders of magnitude larger than the range plotted in figure 13. The magnitude of $\mathcal{Q}$ is large primarily because the thermal conductivity ratio $\kappa_{\rm s}/\kappa_{\rm f}$ is large. This value of $\mathcal{Q}$ yields onset parameters: $(\mathcal{C}_{\rm o},\mathcal{B}_{\rm o},\mathcal{S}_{\rm o})=(1.187\times{}10^{5},1.03\times{}10^{4},8.65\times{}10^{3})$. Using these onset values as a rough guide to estimate experimental conditions yields the parameter values: $(\Delta,\bar{h},d)=(1.85\textrm{ K},2.05\times{}10^{2}\textrm{ m},6.90\times{}10^{5}\textrm{ m})$. Although the temperature difference ($\Delta$) is feasible, the thicknesses are clearly not.

We now shift our focus to a substrate material that does lead to oscillatory instabilities under experimentally feasible experimental conditions. PMMA is a readily-available insulating material with low thermal conductivity and diffusivity. Again using the properties from Table 1, the silicone oil-PMMA system has a value of $\mathcal{Q}=6.45\times{10}^{-7}$. Substituting this value of $\mathcal{Q}$ into figure 13 yields the onset parameter values, which can then be used to estimate experimental conditions:

[TABLE]

The above dimensional variables (temperature and thicknesses) provide a guide for predicting the range of experimental values for which oscillatory instabilities occur. Figure 14 presents results for the silicone oil-PMMA system with different experimental parameters $\Delta,\bar{h},d$. The figure varies the parameter $0\leq\mathcal{C}\leq 0.3$ which corresponds to a dimensional temperature (difference) range $0\leq\Delta\leq 26.5$ K. In particular, the top row of figure 14 plots the real values of several important unstable frequencies versus $\mathcal{C}$ at different $\mathcal{B},\mathcal{S}$ values: $\Omega_{\rm f,max}$ (maximum real growth rate), $\Omega_{\rm c,max}$ (maximum real growth rate of oscillatory instabilities) and $\Omega_{\rm coll}$. Note that these values are defined in §6. The middle row in figure 14 plots the range of imaginary values $\Omega_{\rm i}$ that are oscillatory unstable. Having information on the possible imaginary values of the complex frequencies that are oscillatory unstable is useful, since these frequencies may be excited via parametric resonance by external forcing. Lastly, the bottom row plots the band of wavenumbers that are unstable, and oscillatory unstable. Plotting the unstable wavenumbers provides information on the length scales (thereby influencing which experimental domain sizes one can utilize) that lead to instability.

Figure 14 demonstrates the effect of temperature difference ($\mathcal{C}$) on oscillatory instabilities. Specifically, the figure plots frequency values $\Omega$ (real parts are in subfigures (a) and (d)) and imaginary parts in subfigures (b) and (e)) and wavenumbers $\beta$ versus $\mathcal{C}$ (in subfigures (c) and (f)). The grey regions correspond to values for which only monotonic (non-oscillatory) instabilities occur, while the pink regions correspond to values at which oscillatory instabilities may occur. Note that since there are an infinite number of roots $\Omega_{n}(\beta)$, the pink regions – which always have oscillatory instabilities – may also contain monotonic instabilities as well as oscillatory instabilities. Figure 14 contrasts the stability behavior of $\Omega$ and $\beta$ versus $\mathcal{C}$ for two different sets of $\mathcal{B}$ and $\mathcal{S}$ values. Here the values of $\mathcal{B}=\mathcal{B}_{\rm o}$ and $\mathcal{S}=\mathcal{S}_{\rm o}$ were chosen to ensure oscillatory instabilities in the subfigures (a–c), while $\mathcal{B}=\mathcal{B}_{\rm g}$ and $\mathcal{S}=\mathcal{S}_{\rm g}$ were chosen to ensure that global oscillatory instabilities occur for a range of $\mathcal{C}$ values in subfigures (d–f). The dashed lines in the subfigures correspond to $\beta_{\rm coll}$ (in (a), (d)) and $\Omega_{\rm coll}$ (in (c) and (f)) and are defined in equation (6.3). The solid lines plot the most unstable wavenumber and frequencies for which global oscillatory instabilities occur; and correspond to the variables $\beta_{\rm max}$, $\Omega_{\rm r,max}$ and $\Omega_{\rm i,max}$ as defined in equation (6.6). Figure 14 shows that the range of unstable wavenumbers and frequencies increase with $\mathcal{C}$, and that in general, large $\mathcal{C}$ values tend to drive oscillatory instabilities.

7.4 Summary

We now recapitulate the most important results of this section. First, this section classifies oscillatory instabilities, for any $(\mathcal{Q},\mathcal{C},\mathcal{B},\mathcal{S})$, as Type I, Type II, and subsequently determines whether they are globally oscillatory unstable. Several important conclusions can be reached:

• •

The choice of materials (i.e. the fluid and the substrate) dictate the parameter $\mathcal{Q}$. The value of $\mathcal{Q}$ then guides which experimental conditions, such as the film and substrate thicknesses, as well as temperature difference, lead to oscillatory instability. Generally speaking, the onset values $\mathcal{S}_{\rm g}$ and $\mathcal{B}_{\rm g}$ provide guides (i.e. order-of-magnitude estimates) that can be used as minimum material thicknesses.

• •

Temperature drives instability. It is well-known that temperature gradients can drive instabilities in fluids (i.e. Rayleigh-Benard convection). This result is also true in the current setting: oscillatory instabilities are more likely when there is a larger temperature difference across the film and substrate. Crucially, we find that oscillatory instabilities, or global oscillatory instabilities arise for $\mathcal{C}$ values (i.e. temperature values) that exceed predefined thresholds $\mathcal{C}>\mathcal{C}_{\rm o}$, or $\mathcal{C}>\mathcal{C}_{\rm g}$, respectively.

• •

Insulating substrates are more likely to give rise to oscillatory instabilities than conducting substrates. The physical reason is that substrates that are thermally conducting transfer heat, and consequently equilibrate their temperatures, over time scales much faster than the characteristic time scales in the thin film. Oscillatory instabilities require thermal coupling between the substrate and the film, and can occur when the natural time scales of the film are on the same order as the time scale governing thermal diffusion in the substrate.

8 Discussion and conclusion

In this work we derived a nonlinear model that couples the thermocapillary dynamics of a liquid film heated by a thermally conductive and diffusive substrate. This was done by assuming a large substrate-to-film thermal conductivity ratio and a substrate thickness that is asymptotically larger than both the mean film thickness and the characteristic lateral disturbance. In order to highlight parameter regimes that are subject to oscillatory instabilities, a scaling was incorporated that grouped the effects of the substrate thermal diffusivity, the imposed temperature difference, the film thickness, and the substrate thickness via four separate dimensional parameters: $(\mathcal{Q},\mathcal{C},\mathcal{B},\mathcal{S})$.

For any set of model parameters, linear stability of the model can be described by the wavenumber-dependent interaction between a perturbation associated with the governing film evolution equation and an infinite number of perturbations associated with the substrate heat equation. The film root coalesces with the root $\Omega_{1}(\beta)$ at certain wavenumbers; at these points the participating roots bifurcate into the complex plane and become oscillatory unstable. To investigate the emergence of these instabilities, complex numerical continuation and optimization algorithms were used in section 7 to describe the emergence of oscillatory instabilities. Notably, we showed that parameter sets subject to a global oscillatory instability occur only for substrate-to-film thickness ratios that are sufficiently large.

We have not provided quantitative predictions of the exact experimental conditions at which one will be able to observe oscillatory thermocapillary instability. This is due to the difficulties associated with prescribing the heat transfer coefficient $q$ that describes the rate of heat transfer between the film and the bounding cold gas layer. In many cases it is probably not possible to determine this parameter prior to running an experiment. One way to handle this issue could be to first conduct experiments on thin substrates, and then use the observed instability wavelength to determine $q$. Still, several workers have highlighted the weakness of Newton’s Law of Cooling, in particular, due to transport effects in the gas layer, see VanHook et al. (1997), or in applications where heat transfer rates are large, see Besson (2012). Therefore, we suggest that the results of our work serve as a foundation upon which more elaborate models can be developed.

To conclude we note that, although most of the oscillatory instabilities we discussed in this work were not global, the coupled model may also serve as a foundation upon which time-periodic excitation could be investigated as a means of driving instability.

Appendix A Thin substrate limit

For sufficiently thin substrates, lateral heat conduction and the thermal diffusivity can be neglected and a single nonlinear PDE can be derived for the evolution of the local film thickness. Instead of re-deriving the long wave-model starting with these these assumptions, we equivalently obtain its dispersion relation by taking the limit of (4.25) as $\sqrt{\Omega+\beta^{2}}\rightarrow{0}$, viz.,

[TABLE]

Effectively we have restricted consideration to film and substrate temperature profiles depending only on the vertical coordinate, as described by the basic state solutions (4.8). In doing so, the full dispersion relation reduces to an explicit expression for strictly real values of ${\Omega}$ in terms of the model parameters.

The dimensional equivalent to (A.1) is obtained by making substitutions (3.1) and (3.2) and solving for the dimensional growth rate $\omega$ as a function of the wavenumber $k$, viz.,

[TABLE]

Aside from $\chi_{\rm s}$ (negligible for thin substrates), this expression for $\omega(k)$ describes the influence of material properties and dimensions on film stability. It is clear that viscosity modifies only the growth rate. Solving for the cutoff wavenumber $k_{\rm c}$ at which $\omega=0$ we obtain

[TABLE]

This wavenumber divides the continuous bands of unstable ($0<k<k_{\rm c}$) and stable ($k>k_{\rm c}$) wavenumbers for a given set of system parameters. Instability described by (A.3) is clearly driven by increasing values of the coefficient $\gamma\,\Delta/\sigma_{0}$ and decreasing film thickness $\overline{h}$.

It is also evident that $\lim_{q\to{0}}k_{\rm c}=\lim_{q\to\infty}k_{\rm c}=0.$ For $q\rightarrow{0}$, the resistance to heat transfer at the film-gas interface becomes infinite and, as a result, the perturbed free surface is uniformly equal to the blackbody temperature $\psi_{\rm b}$. In the absence of variations in the free surface temperature, no thermocapillary stresses arise and perturbations of all wavelengths are stable. Likewise, all values of $k$ are stabilized in the limit $q\rightarrow{\infty}$, which uniformly sets the free surface temperature to the gas temperature $\theta_{\rm g}$. For finite values of $q$, interfacial resistance to heat transfer introduces variations in the free surface temperature that depend locally on the perturbed film thickness. Specifically, with the rate at which heat is removed from the film fixed by $q$, local hot and cold spots form at troughs and crests, respectively, due to their relative proximities to the heating source.

Finally, inspecting the limits

[TABLE]

we see that $k_{\rm c}$ is maximized for situations that effectively transfer the isothermal blackbody temperature directly to the film-substrate interface ($\mathcal{S}\rightarrow{0}$). We conclude by stating that placing a substrate between a film and the blackbody necessarily stabilizes films for finite values of $d$ and $\kappa_{\rm s}$ relative to the case of heating a film directly without a substrate.

Acknowledgments

The research was supported by a fellowship from the New Jersey Institute of Technology Department of Mathematical Sciences (Batson); by NSF CBET–1604351 (Batson, Kondic); by NSF DMS–1815613 (Cummings, Kondic); and by NSF DMS–1719693 (Shirokoff). D. Shirokoff was supported by a grant from the Simons Foundation ($\#359610$).