1. Introduction
The dynamics of thin liquid films is a topic of extensive interest with a number of applications ranging from biomedical (Chen Reference Chen2016; Li & Chu Reference Li and Chu2016) to electronic coatings and nanotechnology (Zhang Reference Zhang2010). The inclusion of thermal effects in thin film dynamics, relevant for many applications, is a mathematically challenging problem. To develop a realistic model one must consider multiple factors, such as the heat supply mechanism(s), possible dependence of material parameters on temperature, heat loss mechanisms and phase changes. When the liquid of interest is placed upon a thermally conductive substrate one must also account for the heat flow within the substrate as well as the interaction between the liquid and the substrate. Numerous models have been developed to address these complications using continuum theory, which in general describes both the thermodynamics and fluid dynamics in terms of partial differential equations (PDEs), derived from first principles. In situations where there is a small aspect ratio (ratio of typical film thickness to typical lateral length scale of interest) long-wave theory (LWT) may be used, which effectively enables the fluid dynamics problem to be reduced to a fourth-order PDE for film thickness. LWT has already proved very valuable in a variety of settings such as liquid crystals, paint coatings, tear films, nanotechnology and many others (see Craster & Matar (Reference Craster and Matar2009) for a comprehensive review). Due to the variety of length and time scales present, the applicability of LWT to the problem of heat conduction in a thin liquid film is not always clear cut. Of the issues outlined above we highlight the following in this work: (i) the influence of temperature on film evolution; (ii) heating/cooling mechanisms; and (iii) the application of LWT to heat conduction.
Various thermal effects that may influence the evolution of the film thickness have been considered in prior work. For an isothermal nanoscale film the primary dewetting mechanism is liquid–solid interaction, often modelled by a disjoining pressure (see Israelachvili (Reference Israelachvili1992) for an extended review). For non-isothermal films, gradients in temperature may give rise to surface tension gradients (thermocapillary or Marangoni effects), which develop when heating from below (Scriven & Sternling Reference Scriven and Sternling1964) and can destabilize the film. The work of Shklyaev, Alabuzhev & Khenner (Reference Shklyaev, Alabuzhev and Khenner2012) finds novel stability thresholds between monotonic and oscillatory instabilities (in both cases, in the linear regime the instability grows as 
$\textrm {e}^{\omega t}$ with growth rate 
$\omega$, and positive real part, 
${\rm Re} (\omega)>0$; but the imaginary part 
${\rm Im} (\omega)$ is zero in the former case and non-zero in the latter) that also account for heat losses from the free surface of the film (referred to here as radiative heat losses). In that work, the film is heated from below via a constant heat flux from a substrate of much lower thermal conductivity. Batson et al. (Reference Batson, Cummings, Shirokoff and Kondic2019) perform a stability analysis similar to that of Shklyaev et al. (Reference Shklyaev, Alabuzhev and Khenner2012), but model the substrate explicitly rather than as a simple boundary condition. They solve a full heat equation for the substrate temperature, and find that oscillatory instabilities arise primarily due to thermal coupling between the film and the substrate. A number of other works have considered the coupling between the evolution of film and substrate temperatures. Saeki, Fukui & Matsuoka (Reference Saeki, Fukui and Matsuoka2011), for example, consider a film/substrate system heated by a laser and find that the rate of change of film reflectivity 
$R$ with thickness 
$h$, 
$\mathrm {d}R/\mathrm {d}h$, may promote either stability or instability of the film depending on the sign of 
$\mathrm {d}R/\mathrm {d}h$. The magnitude of the incident laser energy was earlier shown to influence film thickness evolution by Oron (Reference Oron2000), who showed in particular that increasing the laser energy can partially inhibit film instability.
Another important effect that may influence film evolution is the dependence of material parameters, such as density, thermal conductivity, surface tension, heat capacity and viscosity, on temperature. These relationships are often assumed to be linear, although a strongly nonlinear Arrhenius-type dependence of viscosity on temperature may exist. Oron, Davis & Bankoff (Reference Oron, Davis and Bankoff1997) formulated a thin film model in which viscosity variation is included, and in later work Seric, Afkhami & Kondic (Reference Seric, Afkhami and Kondic2018) found that film evolution is strongly affected by the inclusion of temperature-dependent viscosity. If temperature variations are sufficiently large, the film may undergo a phase change (liquefaction or solidification). This has been considered using a variety of approaches, for example Trice et al. (Reference Trice, Thomas, Favazza, Sureshkumar and Kalyanaraman2007) use a latent heat model to describe such phase change whereas others, such as Seric et al. (Reference Seric, Afkhami and Kondic2018), assume phase change to be instantaneous.
Modelling of heat losses in a liquid film often focuses on the boundary effects, since viscous dissipation can usually be ignored. Radiative heat losses from the liquid to the surrounding medium are typically modelled by a Robin type boundary condition (Oron Reference Oron2000; Atena & Khenner Reference Atena and Khenner2009; Saeki et al. Reference Saeki, Fukui and Matsuoka2011; Shklyaev et al. Reference Shklyaev, Alabuzhev and Khenner2012; Saeki, Fukui & Matsuoka Reference Saeki, Fukui and Matsuoka2013) whereas the heat loss/gain from the film to the substrate has been modelled variously by (i) a constant temperature (Oron & Peles Reference Oron and Peles1998; Oron Reference Oron2000; Saeki et al. Reference Saeki, Fukui and Matsuoka2011), (ii) constant flux (Atena & Khenner Reference Atena and Khenner2009; Shklyaev et al. Reference Shklyaev, Alabuzhev and Khenner2012) or (iii) continuity of temperatures and fluxes, known as perfect thermal contact (Trice et al. Reference Trice, Thomas, Favazza, Sureshkumar and Kalyanaraman2007; Saeki et al. Reference Saeki, Fukui and Matsuoka2011; Dong & Kondic Reference Dong and Kondic2016; Seric et al. Reference Seric, Afkhami and Kondic2018). The choice of boundary conditions plays an important role when formulating and solving equations to describe the heat flow.
In many cases an asymptotic approach may be adopted, giving rise to simplified leading-order temperature equation(s). The work of Saeki et al. (Reference Saeki, Fukui and Matsuoka2011), for example, includes both radiative heat losses and heat transfer at the liquid–solid interface, and gives rise to a depth-averaged (
$z$-direction) equation for film temperature, which retains parametric 
$z$ dependence even when radiative heat losses are ignored. In later work, Saeki et al. (Reference Saeki, Fukui and Matsuoka2013) developed similar leading-order equations for film temperature when the film is optically transparent. In this case the film temperature dependence on 
$z$ is slaved to the inclusion of radiative heat losses. Trice et al. (Reference Trice, Thomas, Favazza, Sureshkumar and Kalyanaraman2007), on the other hand, conclude that using a 
$z$-independent film temperature model is sufficient when radiative heat losses can be neglected and film-to-substrate heat losses are dominant (e.g. when there is a high thermal conductivity ratio between the film and substrate). These previous works demonstrate that boundary conditions play an integral role in the asymptotic formulation of a model and may facilitate simple models that eliminate 
$z$-dependence (e.g. Shklyaev et al. Reference Shklyaev, Alabuzhev and Khenner2012).
Due to the small aspect ratio of the film, a commonly used ‘reduced’ model for heat conduction is one that neglects in-plane diffusion altogether (Trice et al. Reference Trice, Thomas, Favazza, Sureshkumar and Kalyanaraman2007; Dong & Kondic Reference Dong and Kondic2016; Seric et al. Reference Seric, Afkhami and Kondic2018). This model, which we refer to here as (1D), is much simpler than a model that includes full heat diffusion and is typically justified by arguing that in-plane diffusion occurs on a much longer time scale than that of out-of-plane diffusion. Alternative simplified models have also been proposed. The work of Shklyaev et al. (Reference Shklyaev, Alabuzhev and Khenner2012), for example, uses LWT to derive evolution equations for heat conduction that differ significantly from the one-dimensional (1-D) model of Dong & Kondic (Reference Dong and Kondic2016) (and from the asymptotic model considered in this paper). Atena & Khenner (Reference Atena and Khenner2009) also derive leading-order temperature equations that do not rely on the 1-D approximation. More recent work by Seric et al. (Reference Seric, Afkhami and Kondic2018) briefly compares predictions from (1D) with those using a full thermal diffusion model, and suggests that (1D) performs poorly by comparison, though the analysis is far from complete. Despite the extensive literature, the scenarios for which thermal diffusion model (1D) is valid remain unclear. A key objective of the present paper is to present a thermal model for thin film flow that includes in-plane heat conduction at reasonable computational complexity and to compare with both (1D) and with the full heat diffusion model (which serves as a benchmark).
In this paper, we consider films placed upon a thermally conductive substrate and heated by a laser. Heat generation by a laser source is complicated to model and requires in general that one accounts for the optical properties of the film, such as reflectivity, transmittance, and absorption. These properties may depend on refractive indices of the air, film and substrate, as well as the respective extinction coefficients. Again, various modelling approaches have been taken in the literature: we note for example that Saeki et al. (Reference Saeki, Fukui and Matsuoka2011, Reference Saeki, Fukui and Matsuoka2013) present a detailed model for laser energy that involves complicated expressions for the optical properties; whereas Trice et al. (Reference Trice, Thomas, Favazza, Sureshkumar and Kalyanaraman2007) propose a simpler approach (to be discussed later) in which these properties are approximated. An important application of laser heating is pulsed laser-induced dewetting (PLiD) of metal films. The mechanism by which liquid metals evolve into assemblies of droplets has been explored via experiments (Henley, Carey & Silva Reference Henley, Carey and Silva2005), simulations (Dong & Kondic Reference Dong and Kondic2016; Seric et al. Reference Seric, Afkhami and Kondic2018) and theory (Trice et al. Reference Trice, Thomas, Favazza, Sureshkumar and Kalyanaraman2007) with applications ranging from nanowire growth (Kim et al. Reference Kim2009; Ross Reference Ross2010; Shirato et al. Reference Shirato, Strader, Kumar, Vincent, Zhang, Karakoti, Nacchimuthu, Cho, Seal and Kalyanaraman2011), to plasmonics (Halas et al. Reference Halas, Lal, Chang, Link and Nordlander2011) and photovoltaics (Atwater & Polman Reference Atwater and Polman2010); see also Hughes, Menumerov & Neretina (Reference Hughes, Menumerov and Neretina2017) and Makarov et al. (Reference Makarov, Milichko, Mukhin, Shishkin, Zuev, Mozharov, Krasnok and Belov2016) for recent application-centred reviews, and Ruffino & Grimaldi (Reference Ruffino and Grimaldi2019) and Kondic et al. (Reference Kondic, Gonzalez, Diez, Fowlkes and Rack2020) for reviews focusing on molten metal film instabilities. Of late, PLiD has been used to organize nanoparticles into patterns of droplets via Rayleigh–Plateau-type instabilities (Favazza, Kalyanaraman & Sureshkumar Reference Favazza, Kalyanaraman and Sureshkumar2006a; McKeown et al. Reference McKeown, Roberts, Fowlkes, Wu, LaGrange, Reed, Campbell and Rack2012; Ruffino et al. Reference Ruffino, Pugliara, Carria, Bongiorno, Spinella and Grimaldi2012), induced by exposing metal films/filaments upon (typically) Si/SiO
$_2$ substrates to laser irradiation, effectively liquefying the film for tens of nanoseconds. The liquefied film breaks up into droplet patterns, which then resolidify, freezing the patterns in place. Thermal effects are found to be highly relevant, influencing the stability, evolution and final (solidified) configurations of molten metal films (see for example, Trice et al. Reference Trice, Thomas, Favazza, Sureshkumar and Kalyanaraman2007). A number of experimental studies have considered metallic systems such as Co (Favazza et al. Reference Favazza, Kalyanaraman and Sureshkumar2006a,Reference Favazza, Trice, Gangopadhyay, Garcia, Sureshkumar and Kalyanaramanb,Reference Favazza, Trice, Krishna and Kalyanaramanc; Trice et al. Reference Trice, Thomas, Favazza, Sureshkumar and Kalyanaraman2007), Ag (Krishna et al. Reference Krishna, Sachan, Strader, Favazza, Khenner and Kalyanaraman2010), Au (Yadavali, Khenner & Kalyanaraman Reference Yadavali, Khenner and Kalyanaraman2013), Ni (Fowlkes et al. Reference Fowlkes, Kondic, Diez, Gonzalez, Wu, Roberts, McCold and Rack2012) as well as multi-metal systems (Fowlkes, Wu & Rack Reference Fowlkes, Wu and Rack2010). The large variety of experimental work that has been done on nanoscale metal films calls for a firm theoretical foundation, which can both explain existing results and suggest new approaches.
The focus of the present paper is development of a consistent, asymptotically valid, mathematical model that accounts for (i) heat absorption that is influenced by the local value of (time-dependent) film thickness; (ii) in-plane and out-of-plane heat diffusion in a tractable manner; (iii) self-consistent coupling of the heat flow and film evolution; and (iv) thermal variation of material properties, in particular of surface tension and viscosity. LWT is used to reduce modelling of the film evolution to a fourth-order PDE for the film thickness and to develop an asymptotic model for heat conduction. We consider a set-up where the primary heat loss mechanism is through the substrate rather than the liquid–air interface, and the thermal conductivity of the film is much higher than of the substrate, as appropriate for metal films on SiO
$_2$ substrates. We will show that the proposed model (called asymptotic model (A) in what follows) produces accurate results, while avoiding the shortcomings of models that ignore coupling of fluid dynamics and thermal transport and producing results with a reasonable computational effort. It should be emphasized that the use of a more complex model (called full (F) model below) is orders of magnitude more computationally expensive (even for small computational domains, the computing time is measured in days on a modern computational workstation (in a serial mode)). Our asymptotic model provides essentially indistinguishable results at a fraction of the computational cost.
The rest of the paper is organized as follows. In § 2, we formulate a general mathematical model by introducing appropriate scales, the corresponding dimensionless system, and relevant dimensionless parameter groups. We present three different models of heat conduction: a full diffusion model (F), a 1-D diffusion model (1D) and an asymptotic model (A); and we summarize the derivation of the thin film evolution equation (the fluid mechanical model always used), accomplished using LWT and accounting for the possibility of temperature dependence of material parameters. Section 3 contains our main results. In § 3.1, we perform linear stability analysis (LSA) on the film evolution equation to understand the circumstances under which disturbances to the liquid film lead to instabilities, and to predict the manner of film breakup. In § 3.2, we summarize the conditions under which our simulations are carried out, and in § 3.3, we display results comparing the three models for heat conduction. In § 3.4, we restrict attention to the asymptotic model for heat conduction and study how temperature dependence of both viscosity and surface tension influence the results. We find that temperature dependence of the viscosity has the most significant effect on the instability development, while temperature-induced variation of surface tension plays only a minor role. Furthermore, in the physically relevant regime, allowing viscosity to vary with temperature produces films that dewet fully in the liquid phase, while if viscosity is fixed at its melting temperature value the dewetting occurs much closer to the solidification time, which may result in only partial drop formation. We conclude in § 4 with a brief summary and discussion.
2. Model formulation
Consider a molten metal film (assumed initially solid) of characteristic lateral length scale 
$L$, and (nanoscale) thickness 
$H$, heated by a laser, and in contact with a thermally conductive solid substrate of finite thickness, which itself rests upon a much thicker Si slab. The basic set-up is sketched in figure 1. Here, we consider the substrate to be thin, and comparable in size to the film thickness, 
$H$. We define the aspect ratio of the film to be 
$\epsilon =H/L \ll 1$.
Figure 1. Three-dimensional schematic of the film, substrate and laser system. In dimensionless units the mean film thickness is equal to 
$1$, the substrate thickness is given by 
$H_{s}$ and the domain width is 
$2N{\rm \pi}$ (both 
$N=1$ and 
$N=20$ will be used in simulations). The model is presented in three dimensions but for simplicity simulations will be performed only in two dimensions.
In the following, we refer to the in-plane coordinates as 
$x,y$ and the out-of-plane coordinate as 
$z$. For completeness, we present the governing equations for a 3-D system, though the results presented in this paper will be for the 2-D case in which all quantities are independent of 
$y$. We define 
$L$, 
$H$, 
$U$, 
$\epsilon U$, 
$t_{s}$, 
$T_{melt}$, 
$\mu _{f} U/(\epsilon ^2 L)$ and 
$\gamma _{f}$ (where 
$\mu _{f}$ and 
$\gamma _{f}$ are the viscosity and surface tension of the film at melting temperature, 
$T_{melt}$) to be the in-plane length scale, out-of-plane length scale, in-plane velocity scale, out-of-plane velocity scale, time scale, temperature scale, pressure scale and surface tension scale, respectively (the values of the material parameters used are given in table 1). Similar to Gonzalez et al. (Reference Gonzalez, Diez, Wu, Fowlkes, Rack and Kondic2013), we set 
$t_{s}=3 \mu _{f} L/ (\epsilon ^3 \gamma _{f})$, which can be interpreted as a viscous time scale. The in-plane velocity scale is fixed as 
$U=\epsilon ^3 \gamma _{f}/(3 \mu _{f})$ so that 
$t_s=L/U$. The length scale is fixed as 
$L=\lambda _{m}/(2{\rm \pi} )$, where 
$\lambda _{m}$ is the most unstable wavelength obtained from LSA with surface tension and viscosity fixed as 
$\gamma _{f}$ and 
$\mu _{f}$, respectively (see § A.2 for details). We treat the film as an incompressible Newtonian fluid and assume that viscosity is independent of 
$z$. The resultant dimensionless system then comprises the following fluid equations, which hold on 
$0 < z < h$,
\begin{gather} \epsilon^2 {{Re}} \left( \partial_t u + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} u \right) ={-}\partial_x p + \epsilon^2 \boldsymbol{\nabla}_2 \cdot \left( \mathcal{M} \boldsymbol{\nabla}_2 u \right) + \mathcal{M} \partial_z^2 u + \epsilon^2 \partial_y \mathcal{M} \partial_x v, \end{gather}
\begin{gather}\epsilon^2 {{Re}} \left( \partial_t v + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} v \right) ={-}\partial_y p + \epsilon^2 \boldsymbol{\nabla}_2 \cdot \left( \mathcal{M} \boldsymbol{\nabla}_2 v \right) + \mathcal{M} \partial_z^2 v + \epsilon^2 \partial_x \mathcal{M} \partial_y u, \end{gather}
\begin{gather}\epsilon^4 {{Re}} \left( \partial_t w + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} w \right) ={-}\partial_z p + \epsilon^4 \boldsymbol{\nabla}_2 \cdot \left( \mathcal{M} \boldsymbol{\nabla}_2 w \right) + \epsilon^2 \mathcal{M} \partial_z^2 w + \epsilon^2 \boldsymbol{\nabla}_2 \mathcal{M}\cdot \partial_z(u,v), \end{gather}
\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}the following equations of heat conduction,
\begin{gather} \epsilon^2 {Pe}_{f} \partial_t T_{f} = \epsilon^2 \nabla_2^2 T_{f} + \partial_z^2 T_{f} + \epsilon^2 Q,\quad \mbox{for}\ z \in \left( 0,h \right), \end{gather}
\begin{gather}{Pe}_{s}\partial_t T_{s} = \epsilon^2 \nabla_2^2 T_{s} + \partial_z^2 T_{s},\quad \mbox{for}\ z \in \left({-}H_{s},0 \right), \end{gather}and boundary conditions,
\begin{gather} w=\partial_t h + u\partial_x h + v \partial_y h,\quad \mbox{on}\ z=h, \end{gather}
\begin{gather} \left[ \boldsymbol{n} \cdot {\boldsymbol{\mathsf{T}}} \cdot \boldsymbol{n} \right] = 3\gamma \left(\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{n} \right) - 3\varPi(h),\quad \mbox{on}\ z=h, \end{gather}
\begin{gather}\epsilon \left[ \boldsymbol{n} \cdot {\boldsymbol{\mathsf{T}}} \cdot \boldsymbol{t}_1 \right] ={-}3 \partial_x {\gamma},\quad \mbox{on}\ z=h, \end{gather}
\begin{gather}\epsilon \left[ \boldsymbol{n} \boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}} \boldsymbol{\cdot} \boldsymbol{t}_2 \right] ={-} 3\partial_y {\gamma},\quad \mbox{on}\ z=h, \end{gather}
\begin{gather}\boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\nabla} T_{f} = 0,\quad \mbox{on}\ z=h, \end{gather}
\begin{gather}\partial_z T_{f} = \mathcal{K} \epsilon^2 \partial_z T_{s},\quad T_{f} = T_{s},\enspace \mbox{on}\ z=0, \end{gather}
\begin{gather}\boldsymbol{u} = \boldsymbol{0},\quad \mbox{on}\ z=0, \end{gather}
\begin{gather}\partial_z T_{s} = {Bi} \left( T_{s} - T_{a}\right),\quad \mbox{on}\ z={-}H_{s}, \end{gather}
\begin{gather}\partial_x T_{f} = \partial_x T_{s}=0,\quad \mbox{on}\ x={\pm} N {\rm \pi}, \end{gather}
\begin{gather}\partial_y T_{f} = \partial_y T_{s}=0,\quad \mbox{on}\ y={\pm} N {\rm \pi}. \end{gather}
Here, the fluid velocity is given by 
$\boldsymbol {u}=(u,v,w)$, pressure by 
$p$ and film and substrate temperatures by 
$T_{f}$ and 
$T_{s}$, respectively. Subscripts 
${f}$ and 
${s}$ stand for film and substrate, respectively, unless otherwise stated and 
$\boldsymbol {0}=(0,0,0)$. We refer to the gradient operator as 
$\boldsymbol {\nabla }=(\partial _x, \partial _y, \partial _z)$, its in-plane counterpart as 
$\boldsymbol {\nabla }_2=(\partial _x,\partial _y,0)$ and the in-plane Laplacian operator as 
$\nabla _2^2$ (defined by 
$\nabla _2^2 u = \partial _x^2 u + \partial _y^2 u$ for a given scalar function 
$u$). Equations (2.1)–(2.6) are the Navier–Stokes (NS) equations representing conservation of mass and momentum for the film, together with thermal energy conservation in both film (
${0 < z < h}$) and substrate (
$-H_{s} < z < 0$) domains, both of lateral extent 
$2N{\rm \pi}$, 
$-N{\rm \pi} < x,y < N{\rm \pi}$ (for simulations we use either 
$N=1$ or 
$N=20$, but 
$N$ can be any positive integer). The unit vector 
$\boldsymbol {n}$ denotes the outward normal to the film's free surface, 
$z=h$. The equations above introduce the following dimensionless parameters:
\begin{equation} \left.\begin{gathered} {Re} = \frac{\rho_{f}U L}{\mu_{f}},\quad \mathcal{M} = \frac{\mu}{\mu_{f}},\quad \mathcal{K}=\frac{k_{s}}{k_{f}}\epsilon^{{-}2}, \\ {Pe}_{f} = \frac{\left( \rho c \right)_{f} U L }{k_{f} },\quad {Pe}_{s} = \frac{\left( \rho c \right)_{s} U \epsilon H }{k_{s}},\quad {Bi} = \frac{\alpha H}{k_{s}}, \end{gathered}\right\} \end{equation}
which are the Reynolds number, dimensionless viscosity, (scaled) thermal conductivity ratio, Péclet numbers and Biot number, respectively. We assume 
$\epsilon ^2 {Re} \ll 1$; the remaining quantities in (2.17) are assumed 
$O(1)$. For further discussion of the choice of scales and parameter values, see table 1, § 2.1.3, and § A.1 later. The definitions of the laser source term 
$Q$ and the disjoining pressure 
$\varPi (h)$ are given in the discussion below. The material parameters 
$(\rho ,c,k)_{{f},{s}}$ represent the density, specific heat capacity and thermal conductivity of the film and substrate, respectively. We assume that the substrate is optically transparent and does not absorb laser energy.
Table 1. Parameters used for liquid Cu film and SiO
$_2$ substrate. References: (1) Dong & Kondic (Reference Dong and Kondic2016), (2) Gonzalez et al. (Reference Gonzalez, Diez, Wu, Fowlkes, Rack and Kondic2013), (3) McKeown et al. (Reference McKeown, Roberts, Fowlkes, Wu, LaGrange, Reed, Campbell and Rack2012).
In experiments a Si substrate is often used (Henley et al. Reference Henley, Carey and Silva2005; Wu et al. Reference Wu, Fowlkes, Roberts, Diez, Kondic, González and Rack2011; Yadavali et al. Reference Yadavali, Khenner and Kalyanaraman2013) on top of which a native layer of oxide, usually SiO
$_2$, 
$3$–
$4$ nm in thickness, typically exists (though an additional oxide layer, typically about 
$100$ nm thick, may also be deposited). Below the oxide in either case is Si, which has a much higher thermal conductivity, and can therefore be assumed isothermal relative to the SiO
$_2$. Consistently, we consider the (SiO
$_2$) substrate to be positioned on top of a thick layer of much higher conductivity, assumed to be at constant ambient temperature, 
$T_{a}$. We model the heat loss from the top (SiO
$_2$) substrate to the thick Si layer below via a Newton law of cooling at 
$z=-H_{s}$ (2.14) with Biot number, 
${Bi}$ (related to the dimensional heat transfer coefficient, 
$\alpha$). The value of 
${Bi}$ was chosen so that the film melts and solidifies on a time scale comparable to the film evolution (although the value 
${Bi}=2\times 10^{-3}$ is presented, the range 
$7 \times 10^{-4}$–
$7\times 10^{-3}$ was considered). We assume the following form of the heat source, 
$Q$ in (2.5), representing the external volumetric heating due to the laser at normal incidence (see Trice et al. Reference Trice, Thomas, Favazza, Sureshkumar and Kalyanaraman2007; Seric et al. Reference Seric, Afkhami and Kondic2018),
\begin{equation}
\left.\begin{gathered} Q = F(t) [ 1-R(h) ]
\exp{[-\alpha_{f}(h-z )]}, \\ F(t) =
C \exp [ -(t-t_{p} )^2/(2\sigma^2) ],
\\ C = \frac{E_0 \alpha_{f} L^2}{ \sqrt{2 {\rm \pi}}\sigma t_{s}
H k_{f} T_{melt}}, \end{gathered}\right\}
\end{equation}
where 
$C$ is a constant (assumed 
$O(1)$) proportional to the amount of incident energy, 
$E_0$, applied from the laser onto the film, 
$\alpha _{f}^{-1}$ is the (scaled) absorption length for laser radiation in the film and 
$F(t)$ captures the temporal power variation of the laser, taken to be a Gaussian pulse centred at 
$t_{p}$ and of width defined by 
$\sigma =t_{p}/ (2 \sqrt {2 \ln {2}})$. Similar to prior work by a number of authors (Oron Reference Oron2000; Trice et al. Reference Trice, Thomas, Favazza, Sureshkumar and Kalyanaraman2007; Saeki et al. Reference Saeki, Fukui and Matsuoka2011, Reference Saeki, Fukui and Matsuoka2013; Dong & Kondic Reference Dong and Kondic2016; Seric et al. Reference Seric, Afkhami and Kondic2018) the transmittance of laser source heating is modelled via the Bouguer–Beer–Lambert law (see, e.g. Howell, Siegel & Menguc Reference Howell, Siegel and Menguc2010), which in (2.18) is presented as a spatially dependent source term, 
$\exp (-\alpha _{f}(h-z))$. In general the reflectivity of the film, 
$R(h)$, on a transparent substrate, can be determined by solving Maxwell's equations with the appropriate boundary conditions (Heavens Reference Heavens1955). The resultant form is quite cumbersome to work with, however, and instead we approximate 
$R(h)$ by the simple functional form (Seric et al. Reference Seric, Afkhami and Kondic2018)
\begin{equation} R(h)=r_0 \left(1-\exp \left(-\alpha_{r} {h} \right) \right), \end{equation}
where 
$r_0$ and 
$\alpha _{r}$ are dimensionless fitting parameters, determined by a least-squares fit of the approximate 
$R(h)$ to the full expression for reflectivity.
Equations (2.7)–(2.11) are boundary conditions on the free surface, 
$z=h(x,y,t)$, with unit normal 
$\boldsymbol {n}=\boldsymbol {\nabla } (z-\epsilon h )/|\boldsymbol {\nabla } (z- \epsilon h)|$ and tangent vectors 
$\boldsymbol {t}_1$ and 
$\boldsymbol {t}_2$ given by 
$\boldsymbol {t}_1 = (1,0,\epsilon \partial _x h)/\sqrt {1+\epsilon ^2 (\partial _x h)^2}$ and 
$\boldsymbol {t}_2 =(0,1,\epsilon \partial _y h)/\sqrt {1+\epsilon ^2 (\partial _y h)^2}$. The kinematic boundary condition is given by (2.7); (2.8), (2.9) and (2.10) are the dynamic boundary conditions, representing a balance of stress between the liquid and air phases, where 
${\boldsymbol{\mathsf{T}}}$ is the Newtonian stress tensor, 
$\gamma$ is the surface tension, and 
$\varPi (h)$ is the disjoining pressure representing liquid–solid interaction. Many forms of 
${\varPi }({h})$ are used in the literature (for more information regarding the microscopic nature of the disjoining pressure, we refer the reader to Israelachvili (Reference Israelachvili1992).); we take
\begin{equation} {\varPi}({h})= {K} \left[ \left( \frac{h_*}{h} \right)^n - \left( \frac{h_*}{h} \right)^m \right],\quad {K} = \frac{A_{H} L}{6{\rm \pi} \epsilon \gamma_{f} h_*^3 H^3}, \end{equation}
with equilibrium thickness 
$h_*$, exponents 
$n>m>1$ (we use 
$(n,m) = (3,2)$ since these values were shown by Gonzalez et al. (Reference Gonzalez, Diez, Wu, Fowlkes, Rack and Kondic2013) to be appropriate for liquid metals), and Hamaker constant 
$A_{H}$. We assume that the radiative heat loss from the film to the air is small compared to the heat conduction from the film to the substrate. As a consequence, we neglect heat loss through the liquid–air interface and apply (2.11), an insulating boundary condition. Furthermore, we model the (assumed) primary heat loss mechanism through the interface between the film and the substrate at 
$z=0$ by perfect thermal contact via (2.12). At 
$z=0$, we also assume no slip and no penetration of fluid via (2.13). Finally, we assume that the film and substrate are thermally insulated at the lateral ends, 
$x=\pm N{\rm \pi}$ and 
$y=\pm N{\rm \pi}$.
We will now proceed to simplify the full model as outlined above. We begin in § 2.1 with a discussion of the various models for heat conduction, and derive a leading-order asymptotic model that (we will show) compares well with the full heat conduction model. In § 2.2, we discuss the long-wave approximation for thin films and the inclusion of thermal effects in the resultant thin film equation.
2.1. Thermal modelling
In what follows, we present three different models for the inclusion of thermal effects (the fluid dynamics in all cases will be described by the long-wave model, see § 2.2). In § 2.1.1, we give a ‘full’ model for heat conduction, denoted (F), which includes both in-plane and out-of-plane heat diffusion, but omits both viscous dissipation and thermal advection. As discussed in § 1, a number of previous works have utilized a much simpler model that neglects lateral heat diffusion (e.g. Trice et al. Reference Trice, Thomas, Favazza, Sureshkumar and Kalyanaraman2007; Dong & Kondic Reference Dong and Kondic2016; Seric et al. Reference Seric, Afkhami and Kondic2018). Although the exact relevance of such lateral (in-plane) heat transfer has not yet been carefully analysed, prior work by Seric et al. (Reference Seric, Afkhami and Kondic2018) suggests that it may be important. To study and quantify the possible significance, in § 2.1.2 we describe such a ‘one-dimensional’ model for heat conduction, denoted (1D). Finally, in § 2.1.3 we apply LWT to (F) to develop an ‘asymptotic’ model for heat conduction, (A). This model utilizes key assumptions on the non-dimensional parameters introduced in (2.17) to arrive at a system that is simpler than (F), but unlike (1D) retains lateral heat diffusion. These three models will be compared in § 3.3.
2.1.1. Full model
In order to compare models of heat conduction, we must first declare a model that serves as a benchmark. We refer to (2.5), (2.6), (2.11), (2.12) and (2.14)–(2.16) (despite the presence of terms that may appear asymptotically small with respect to 
$\epsilon$ in comparison to other terms), as the full model (F) for heat conduction.
2.1.2. One-dimensional model
Here, we display the model obtained by neglecting in-plane heat conduction in (F), assuming that the term 
$\epsilon ^2 \nabla _2^2 T_{f}$ may be neglected compared with 
$\partial _z^2 T_{f}$ in (2.5) but retaining all other terms. Equation (2.11) is replaced by 
$\partial _z T_{f}=0$ since 
$\boldsymbol {n}=(0,0,1)+O(\epsilon ^2)$. This yields the following 1-D model (1D) for heat conduction:
\begin{gather} \epsilon^2 {Pe}_{f} \partial_t T_{f} = \partial_z^2 T_{f} + \epsilon^2 Q,\quad \mbox{for}\ z \in \left(0,h \right), \end{gather}
\begin{gather}{Pe}_{s} \partial_t T_{s} = \partial_z^2 T_{s},\quad \mbox{for}\ z \in \left({-}H_{s},0 \right), \end{gather}
\begin{gather}\partial_z T_{f} = 0,\quad \mbox{on}\ z=h, \end{gather}
\begin{gather}\partial_z T_{f} = \mathcal{K}\epsilon^2 \partial_z T_{s},\quad \mbox{on}\ z=0, \end{gather}
\begin{gather}T_{f} = T_{s},\quad \mbox{on}\ z=0, \end{gather}
\begin{gather}\partial_z T_{s} = {Bi} \left( T_{s} - T_{a} \right),\quad \mbox{on}\ z={-}H_{s}, \end{gather}
\begin{gather}\partial_x T_{f} =0,\quad \mbox{on}\ x={\pm} N{\rm \pi}, \end{gather}
\begin{gather}\partial_y T_{f} =0,\quad \mbox{on}\ y={\pm} N{\rm \pi}, \end{gather}
where 
$Q$ is given by (2.18). We note that, although the substrate temperature 
$T_{s}$ only diffuses in the out-of-plane direction, 
$z$, it is still functionally dependent on the in-plane coordinates 
$x,y$, due to (2.25) and the dependence of film temperature 
$T_{f}$ on 
$x,y$ (via dependence on film height 
$h$). It follows from (2.25), (2.27) and (2.28) that 
$\partial _x T_{s}=0$ at 
$x=\pm N{\rm \pi}$, and 
$\partial _y T_{s}=0$ at 
$y=\pm N{\rm \pi}$ automatically. For the rest of the paper, we refer to (2.21)–(2.28) as the (1D) model.
2.1.3. Asymptotic model
Next, we formulate a model of intermediate complexity by carrying out further asymptotic analysis. To do so, we first make a number of assumptions about the non-dimensional parameters defined in (2.17) and provide estimates of time scales based on the parameters given in table 1:
(i)
${Pe}_{f}=O(1)$. The term $\epsilon ^2 {Pe}_{f} = [(\rho c)_{f} H^2/k_{f} ]/t_{s}= t_{{D}_{f}}/t_{s}$ appearing in (2.5) is a ratio of two time scales: $t_{{D}_{f}}$, the time scale of diffusion of heat in the film, and $t_{s}$, the time scale of film evolution. Thus, we assume $t_{{D}_{f}}\ll t_{s}$; heat diffuses rapidly through the film, before any significant film evolution can occur. In our set-up, $t_{{D}_{f}}\approx 1.17\ \textrm {ps}$, whereas $t_{s} \approx 26.86\ \textrm {ns}$;(ii)
${Pe}_{s}=O(1)$. Similar to (i) the Péclet number for the solid layer can be written as a ratio of time scales, ${Pe}_{s}=[(\rho c)_{s} H^2/k_{s}]/t_{s}= t_{{D}_{s}}/t_{s}$, where $t_{D_{s}}$ is the time scale of out-of-plane thermal diffusion in the substrate. We assume that this diffusion occurs on a time scale comparable to that of film evolution. Here, $t_{{D}_{s}}\approx 0.147$ ns. Although this is small relative to $t_{s}$, this assumption ensures that the time derivative is retained in (2.22), which is numerically convenient and has a negligible effect on results;(iii)
${Bi}=O(1)$. The Biot number ${Bi}=(H/k_{s})/(1/\alpha )$ can be interpreted as the ratio of internal thermal resistance due to diffusion, $H/k_{s}$, and external thermal resistance, $1/\alpha$, due to convection away from the boundary $z=-H_{s}$. We assume these internal and external thermal resistances are comparable;(iv)
$\mathcal {K} =k_{s}/(\epsilon ^2 k_{f}) = O(1)$; the film has much higher thermal conductivity than the substrate;(v)
$H_{s} = O(1)$, indicating that the substrate thickness is comparable in size to the film thickness. Hence the substrate is also thin.
The difference in length scales in the problem motivates the idea that in-plane and out-of-plane diffusion can occur on different time scales. As a consequence of the thin substrate assumption, (v), the in-plane diffusion is much slower than that of out-of-plane diffusion. The ratio of the film evolution time scale to that of diffusion is therefore much smaller for in-plane diffusion than out of plane. Consequently, in-plane diffusion can be neglected in the substrate (cf. Seric et al. Reference Seric, Afkhami and Kondic2018).
To obtain an asymptotically valid model we assume the following expansions:
\begin{equation} T_{f} = T_{f}^{(0)} + \epsilon^2 T_{f}^{(1)} + \cdots,\quad T_{s} = T_{s}^{(0)} + \epsilon^2 T_{s}^{(1)} + \cdots, \end{equation}so that, on substituting in (2.5), (2.6), (2.11), (2.12), (2.14), (2.15) and using assumptions (i)–(v) listed above, the leading-order model is given by
\begin{gather} \partial_z^2 T_{f}^{(0)} =0,\quad \mbox{for}\ z\in(0,h), \end{gather}
\begin{gather}{Pe}_{s} \partial_t T_s^{(0)} = \partial_z^2 T_{s}^{(0)},\quad \mbox{for}\ z\in({-}H_{s},0), \end{gather}
\begin{gather}\partial_z T_{f}^{(0)} =0,\quad \mbox{on}\ z=h, \end{gather}
\begin{gather}\partial_z T_{f}^{(0)} =0,\quad \mbox{on}\ z=0, \end{gather}
\begin{gather}T_{f}^{(0)} = T_{s}^{(0)},\quad \mbox{on}\ z=0, \end{gather}
\begin{gather}\partial_z T_{s}^{(0)} = {Bi} ( T_{s}^{(0)} - T_{a} ),\quad \mbox{on}\ z={-}H_{s}, \end{gather}
\begin{gather}\partial_x T_{f}^{(0)} = 0,\quad \mbox{on}\ x={\pm} N{\rm \pi}, \end{gather}
\begin{gather}\partial_y T_{f}^{(0)} = 0,\quad \mbox{on}\ y={\pm} N{\rm \pi}. \end{gather}
Equations (2.30)–(2.33) result in a leading-order film temperature that is independent of 
$z$ but still unknown, 
$T_{f}^{(0)}=T_{f}^{(0)}(x,y,t)$. We must therefore proceed to next order in the asymptotic expansion to obtain a closed model for the leading-order film temperature. Collecting terms at next order in (2.5) yields
\begin{equation}
{Pe}_{f} \partial_t T_{f}^{(0)} = \nabla_2^2 T_{f}^{(0)} +
\partial_z^2 T_{f}^{(1)} + F(t) [1-R(h) ] \exp
[- \alpha_{f} ( h-z ) ],
\end{equation}while the boundary conditions (2.11) and (2.12) at the same order are
\begin{gather} \partial_z T_{f}^{(1)} = \boldsymbol{\nabla}_2 h \cdot \boldsymbol{\nabla}_2 T_{f}^{(0)},\quad \mbox{on}\ z=h, \end{gather}
\begin{gather}\partial_z T_{f}^{(1)} = \mathcal{K} \partial_z T_{s}^{(0)},\quad \mbox{on}\ z=0 . \end{gather}
Since 
$T_{f}^{(0)}$ is independent of 
$z$, we can integrate (2.38) from 
$z=0$ to 
$z=h$. Doing so, and applying the boundary conditions (2.39) and (2.40), gives the following evolution equation for leading-order film temperature:
\begin{equation} h \,{Pe}_{f} \partial_t T_{f} = \boldsymbol{\nabla}_2 \boldsymbol{\cdot} \left( h \boldsymbol{\nabla}_2 T_{f} \right) - \mathcal{K} \left(\partial_z T_{s}\right)\vert_{z=0} + h \bar{Q}, \end{equation}
for 
$x,y \in (-N{\rm \pi} ,N{\rm \pi} )$, where 
$\bar {Q}=h^{-1} \int _{0}^{h} F(t) [1-R(h) ] \exp [ -\alpha _{f} (h-z ) ]\,\textrm {d} z$ is the averaged heat source and the superscripts on 
$T_{f}$, 
$T_{s}$ are dropped for convenience, since now only leading-order quantities are considered. Here, 
$\boldsymbol {\nabla }_2 \boldsymbol {\cdot } ( h \boldsymbol {\nabla }_2 T_{f} )$ in (2.41) describes the lateral heat diffusion, while the terms 
$\mathcal {K} \partial _z T_{s}$ and 
$h \bar {Q}$ represent the heat lost from the film due to contact with the substrate and the generation of heat in the film due to the laser source, respectively. The final asymptotic model for heat conduction is (2.41) in the film, together with
\begin{gather} {Pe}_{s}\partial_t T_{s} = \partial_z^2 T_{s},\quad \mbox{for}\ z \in \left({-}H_{s},0 \right), \end{gather}
\begin{gather}T_{f} = T_{s},\quad \mbox{on}\ z=0, \end{gather}
\begin{gather}\partial_z T_{s} = {Bi} \left( T_{s} - T_{a} \right),\quad \mbox{on}\ z={-}H_s, \end{gather}
\begin{gather}\partial_x T_{f} =0,\quad \mbox{on}\ x={\pm} N{\rm \pi}, \end{gather}
\begin{gather}\partial_y T_{f} =0,\quad \mbox{on}\ y={\pm} N{\rm \pi}. \end{gather}
We note that, even though lateral diffusion is neglected in (2.42), the substrate temperature 
$T_{s}$ remains a function of 
$x$, 
$y$ and 
$z$, the in-plane variation entering through the boundary condition (2.43). By the same reasoning as in the previous section, the lateral end insulating conditions on 
$T_{s}$, (2.15) and (2.16), are satisfied vacuously.
In summary, we have formulated an asymptotic model for heat conduction that exploits the natural geometry of the problem as well as the relative sizes of material parameters (assumptions (i)–(v)). This model, denoted (A), has advantages over both (F) and (1D). By integrating over the 
$z$-direction, a closed model is obtained for a leading-order temperature profile that is independent of 
$z$, simplifying the problem significantly. As a consequence, (A) is considerably less computationally demanding than (F). Solving (F) for the temperature profile throughout the evolving film is cumbersome since the domain is deformable (see the appendix for details): model (A) eliminates this complication since film temperature depends only on the in-plane direction(s) and time. Model (A) is also (as we will see) substantially more accurate and faster to compute than (1D).
A number of other authors have developed reduced models for heat transfer within films, which we now briefly highlight and contrast with our model (A). The models presented by Dong & Kondic (Reference Dong and Kondic2016), Seric et al. (Reference Seric, Afkhami and Kondic2018) and Trice et al. (Reference Trice, Thomas, Favazza, Sureshkumar and Kalyanaraman2007) ignore in-plane diffusion in the substrate, similar to (2.42). Furthermore, all use a Dirichlet boundary condition at the bottom of the substrate rather than the Newton law of cooling used here (2.44). Shklyaev et al. (Reference Shklyaev, Alabuzhev and Khenner2012) arrive at a leading-order temperature equation through arguments similar to ours above. Their model also retains the in-plane diffusion term, 
$\boldsymbol {\nabla }_2 \boldsymbol {\cdot } (h\boldsymbol {\nabla }_2 T_{f})$ in (2.41), but considers radiative heat losses through the liquid–air interface to be dominant rather than the heat loss through the substrate. One important difference between our model (A) and that of Shklyaev et al. (Reference Shklyaev, Alabuzhev and Khenner2012) is that in (A) volumetric heating is considered, which depends on the local value of the film thickness. This fully couples the fluid and thermal problems, whereas the heating mode considered by Shklyaev et al. (Reference Shklyaev, Alabuzhev and Khenner2012) (heating from the substrate below) does not depend directly on the film thickness. Atena & Khenner (Reference Atena and Khenner2009) also assume such volumetric heating but consider the case where the internal heat generation is promoted to leading order so that 
$z$-dependence is retained in the film temperature, leading to a more computationally demanding formulation.
2.2. Free-surface evolution
Each of our heat conduction models couples to the film evolution problem, which must be solved simultaneously. Here we briefly summarize the long-wave approximation that we utilize in all our simulations, which effectively reduces the NS equations to a fourth-order PDE for film thickness, 
$h$. To retain maximum generality and reasonable tractability, we allow both viscosity and surface tension (which appears in boundary conditions (2.8)–(2.10)) to vary with temperature but treat material density, specific heat and thermal conductivity as fixed at their respective values at melting temperature. We present forms for both surface tension and viscosity that utilize the average free-surface temperature, defined for our purposes by
\begin{equation} \bar{T}=\frac{1}{\left( 2N{\rm \pi} \right)^2}\int_{{-}N{\rm \pi}}^{N{\rm \pi}} \int_{{-}N{\rm \pi}}^{N{\rm \pi}}(T_{f}\vert_{z=h})\,\mathrm{d}x\,\mathrm{d}y. \end{equation}We assume that surface tension depends linearly on temperature in the following sense:
\begin{equation} \gamma = 1 +\frac{2{Ma}}{3} (\bar{T} - 1) + \epsilon^2 \frac{2{Ma}}{3} {\rm \Delta} T + O(\epsilon^4) = \varGamma + \epsilon^2 \frac{2{Ma}}{3} {\rm \Delta} T + O(\epsilon^4), \end{equation}
where 
${Ma}$ is the Marangoni number, given by 
${Ma} = (3 \gamma _{T} T_{melt})/(2\gamma _{f})$, where 
$\gamma _T=(\gamma _{f}/T_{melt})\,\mathrm {d}\gamma /\mathrm {d}\bar {T} \vert _{\bar {T}=1}$ is the change in surface tension with temperature when the film (on average) is at melting temperature, 
$\bar {T}=1$ (the factors of 
$2/3$ are used for later convenience); and 
${\rm \Delta} T$ is given by
\begin{equation} {\rm \Delta} T=T_{f}\vert_{z=h}-\bar{T}. \end{equation}
Since 
$\bar {T}$ depends only on time, (2.48) can be interpreted as defining a surface tension that varies in time (at leading order) due to variations in the average temperature, and in space (at higher order), due to spatial variations in temperature. This asymptotic form of 
$\gamma$ proposed in (2.48) provides a consistent balance in the normal and tangential stress balances presented below. The temperature dependence of the dimensionless viscosity, 
$\mathcal {M}={\mu }/\mu _{f}$, is modelled by an Arrhenius-type relationship, which we take as
\begin{equation} \mathcal{M} (t) = \exp\left( \frac{E}{RT_{melt}} \left(\frac{1}{\bar{T}}-1 \right) \right), \end{equation}
where 
$R=8.314\ \textrm {J}\ \textrm {K}^{-1}\ \textrm {mol}^{-1}$ is the universal gas constant and 
$E$ is the activation energy (Gale & Totemeier Reference Gale and Totemeier2004).
To leading order in 
$\epsilon ^2$ the normal and tangential stress balances (2.8)–(2.10) are
\begin{gather} p ={-} 3\varGamma \nabla_2^2 h - 3\varPi(h),\quad \mbox{on}\ z=h, \end{gather}
\begin{gather}\mathcal{M} \partial_z (u,v) = 2{Ma} \boldsymbol{\nabla}_2 \left({\rm \Delta} T \right),\quad \mbox{on}\ z=h. \end{gather}To obtain an evolution equation for film thickness, we express conservation of mass in the form
\begin{equation} \partial_t h + \boldsymbol{\nabla}_2 \boldsymbol{\cdot} \left( h \bar{\boldsymbol{u}} \right)=0, \end{equation}
where 
$\bar {\boldsymbol {u}}$ is the film-averaged (in-plane) velocity, 
$\bar {\boldsymbol {u}}=h^{-1} \int _{0}^{h} (u,v)\,\mathrm {d} z$. To determine 
$u$ and 
$v$, we expand pressure and velocity fields in (2.1)–(2.3) to leading order in 
$\epsilon$, assume 
$\varGamma$ is 
$O(1)$ and apply the boundary conditions (2.51) and (2.52), together with the kinematic condition (2.7), to obtain the leading-order velocity profile,
\begin{equation} \left(u,v \right) = \frac{1}{\mathcal{M}}\left[ \left( \frac{z^2}{2}-zh \right) \boldsymbol{\nabla}_2 p + 2z {Ma}\boldsymbol{\nabla}_2 \left( {\rm \Delta} T \right) \right], \end{equation}
and 
$z$-independent pressure, 
$p$. Equation (2.51), therefore, gives the pressure throughout the layer and 
$\boldsymbol {\nabla }_2 p$ is found by taking the gradient of (2.51). After plugging (2.54) into (2.53) we then arrive at the thin film equation,
\begin{equation}
\partial_t h + \boldsymbol{\nabla}_2 \cdot
\left[\frac{1}{\mathcal{M}} ( h^3
\boldsymbol{\nabla}_2 ( \varGamma \nabla_2^2 h +
\varPi(h) ) + h^2 {Ma} \boldsymbol{\nabla}_2 (
{\rm \Delta} T ) ) \right] = 0.
\end{equation}
Following the time derivative term in (2.55), the terms (from left to right) represent the capillary, disjoining pressure and Marangoni terms, respectively. In general, (2.55) describes the evolution of a nanoscale thin film and is applicable for any of our three thermal models (A), (F) or (1D) by using 
${\rm \Delta} T$ (2.49) and 
$\bar {T}$ (2.47) from the appropriate model.
Equation (2.55) is already sufficiently general to incorporate in-plane variation of viscosity. For model (A), this may be accomplished by using 
$T_{f}^{(0)}$ in place of 
$\bar {T}$ in (2.50). This is an additional advantage of (A) that is not immediately shared by (F) or (1D) (including spatial dependence of viscosity is more complex with these models due to the dependence of temperature on 
$z$, as discussed further in the next section).
Although the choice of scales made at the start of § 2 is standard in the long-wave approximation (e.g. Oron et al. Reference Oron, Davis and Bankoff1997), the introduction of heat conduction adds significant complications, and it is important to check for consistency. For example, to retain surface tension to leading order in (2.55), the velocity scale must be such that 
$\varGamma =O(1)$. This leads to the specific choice of time scale 
$t_{s}$, which may be slower than the (nanoseconds) duration of the Gaussian pulse. Further discussion of the choice of scales is provided in § A.1. We note that as a consequence of our chosen scalings, it would be asymptotically consistent to replace (2.41) and (2.42) by their quasi-steady analogues; however, solving the resulting boundary value problems is numerically more complicated (and does not affect results), hence we retain these time derivatives in the formulation.
3. Results
For simplicity, we limit our considerations to two spatial dimensions, eliminating 
$y$-dependence, so that the film's free surface is at 
$z=h(x,t)$. In § 3.1, we perform LSA, which provides a framework for describing instability growth and motivates our choice of initial film profile(s). Section 3.2 outlines the set-up of the simulations, including the initial conditions and numerical procedures. Sections 3.3 and 3.4 show simulation results for both film and thermal evolution. In § 3.3, we compare the thermal models. In § 3.4, we (almost) exclusively use (A) to solve for heat conduction and allow the surface tension and viscosity to vary with temperature. For what follows, we define the spatially averaged film temperature,
\begin{equation} T_{avg}=\frac{1}{2N{\rm \pi}}\int_{{-}N{\rm \pi}}^{N{\rm \pi}} \frac{1}{h}\int_{0}^{h}T_{f}\,\mathrm{d} z\,\mathrm{d}x, \end{equation}
where the film temperature 
$T_{f}$ is found using model (1D), (F) or (A) (leading-order temperature for (A)) (for model (A) this is exactly the 2-D free-surface average given by (2.47). For models (F) and (1D) it is the average temperature of the entire film). The parameters used are as given in table 1, except where specified otherwise.
3.1. Linear stability analysis
To provide insight into the mechanism by which films dewet, we carry out LSA. Consider a uniform film of height 
$h_0$, perturbed as follows:
\begin{equation} h(x,t)=h_0 \left( 1 + \delta \exp({\textrm{i}kx + \beta t})\right), \end{equation}
where 
$k$ is the wavenumber, 
$\beta$ is the growth rate and 
$\delta \ll 1$ is the amplitude. A more complete analysis could also incorporate independent perturbations to temperature profiles, as was done by Shklyaev et al. (Reference Shklyaev, Alabuzhev and Khenner2012); for simplicity we do not take this approach. We also neglect, for now, the influence of thermal gradients on film instability by setting 
${Ma}=0$ in (2.55). LSA then provides the following dispersion relation:
\begin{equation} \beta
(k) = \frac{1}{\mathcal{M}}h_0^3 k^2 ( 2 - \varGamma
k^2 ).
\end{equation}
From (3.3), it is immediately apparent that viscosity sets the time scale of the perturbation growth/decay. The stability of the film, on the other hand, is controlled by the surface tension. For our purposes we only consider perturbations that grow (
$\beta >0$ when 
$k^2<2/\varGamma$). The wavenumber 
$k_{m}$ corresponding to maximum growth is found from (3.3) by setting 
$\partial \beta /\partial k=0$. The wavelength of maximum growth 
$\varLambda _{m}$ and the maximum growth rate 
$\beta _{m}=\beta (k_{m})$ can then be written in the simple form
\begin{equation} \varLambda_{m} = 2{\rm \pi} \sqrt{\varGamma},\quad \beta_{m}\equiv \frac{h_0^3}{\varGamma \mathcal{M}}. \end{equation}
Since increasing temperature decreases 
$\varGamma$ and 
$\mathcal {M}$ (see (2.48) and (2.50)), higher temperatures will lead to smaller 
$\varLambda _{m}$ and larger 
$\beta _{m}$. In what follows next, 
$\varLambda _{m}$ will be used to define simulation geometries. We note that 
$\varLambda _{m}$ and 
$L$ are related when 
$\varGamma =1$ via 
$\varLambda _{m}=\lambda _{m}/L = 2{\rm \pi}$, where the expression for 
$\lambda _{m}$ is given in § A.2.
3.2. Simulation set-up
Here we describe the details of the simulations. The numerical solution of (2.55) is obtained using an approach adapted from Diez & Kondic (Reference Diez and Kondic2002) with uniform grid size, 
${\rm \Delta} x=h_*$ (
$h_*$ is defined in (2.20a,b) and is fixed for all simulations as 
$h_*=0.1$), which is sufficient to ensure accuracy. Models (F), (1D) and (A) are all solved using central difference spatial discretization. Model (F) utilizes a mapping onto a rectangle to account for the moving boundary (this is not needed for models (1D) and (A) due to lack of in-plane diffusion and lack of 
$z$-dependence, respectively); see § C.1 in appendix for details. For (1D), temporal discretization is performed using the Crank–Nicolson scheme, while for (A) an implicit–explicit scheme is used (see § C.2). Model (F) is solved using an alternating direction implicit (ADI) method, treating mixed derivative terms explicitly. Adaptive time stepping is used to ensure a tolerance of 
$10^{-3}$ maximum allowable relative error in temperature and film thickness. Note that the time-stepping tolerance criteria must be satisfied for both film and heat evolution equations in order to proceed with a successful iteration (a point to which we return later). No-flux boundary conditions 
$\partial _x h=\partial _x^3 h=0$ are imposed at 
$x=\pm N{\rm \pi}$ (
$h\bar {\boldsymbol {u}}=0$ from (2.53)). The domain length, 
$2N{\rm \pi} = N\varLambda _{m}(\varGamma =1)$, is now set by fixing 
$N=1$ or 
$N=20$. For 
$N=1$, the initial film profile is set to represent a small perturbation to a uniform film thickness,
\begin{equation} h(x,0) = h_0 \left( 1 - \delta \cos \left( \frac{2 {\rm \pi}x}{\varLambda_{m}} \right) \right), \end{equation}
where 
$\delta =0.01$. We refer to the corresponding simulations as those with domain length 
$\varLambda _{m}$. For 
$N=20$ the following initial film profile is imposed:
\begin{equation} h(x,0) = h_0 \left[ 1 + \delta \sum_{i=1}^{40} A_i \left( \cos \left( 2{\rm \pi} x/\lambda_i \right) + \sin \left( 2{\rm \pi} x/\lambda_i \right) \right) \right], \end{equation}
where the amplitudes 
$A_i$ are randomly chosen in 
$[-1,1]$ and 
$\lambda _i = 2\varLambda _{m}/i$. Similarly, we refer to simulations that use (3.6) as those with domain length 
$20\varLambda _{m}$. For both values of 
$N$, the film and substrate are each initially set to the ambient temperature,
\begin{equation} T_{f}^{(0)}(x,0)=T_{f}(x,z,0)= T_{s}(x,z,0) = T_{a}. \end{equation} The numerical solutions for 
$T_{f}$ and 
$T_{s}$ are found first, with the film static, since the film is initially solid (
$T_{a}<1$). Once the film is melted (we define this shortly), the solutions for 
$h$, 
$T_{f}$, 
$T_{s}$ are then iterated successively. The flow of the numerical algorithm is as follows:
(i) Film solid and static.
(a) Update
$T_{f}$.(b) Update
$T_{s}$.(c) Repeat the previous 2 steps until melted.
(ii) Film melted.
(a) Set
$\varGamma , \mathcal {M}, {\rm \Delta} T$. Update $h(x,t)$.(b) Update
$T_{f}$.(c) Update
$T_{s}$.(d) Repeat the previous 3 steps until re-solidification.
(iii) End.
Once the film is melted, both film evolution and heat conduction are solved at every time step (although the numerical algorithm allows for a less frequent numerical solution of the temperature equation relative to that of 
$h$). The film evolution is coupled to the temperature profile through the material parameters 
$\varGamma$ (surface tension) and 
$\mathcal {M}$ (viscosity), and the Marangoni term 
${\rm \Delta} T$ in (2.55). The film is allowed to evolve (flow) only when the temperature is everywhere greater than the melting temperature: it then evolves according to (2.55). Subsequently, as the laser heat source decays, the film temperature eventually drops below the melting point and the film re-solidifies. All simulations shown in this paper are ended when the average temperature decreases to solidification temperature, 
$T_{avg}=1$. In what follows, we will be using the liquid lifetime (LL), defined as the time interval during which the average film temperature is above melting (
$T_{avg}>1$).
3.3. Model comparison with fixed parameters
We now compare models (F), (1D) and (A) holding the material parameters fixed. As a basic check, we first consider a stationary flat film (
$h=h_0$), with material parameters fixed at the values corresponding to the melting temperature (
$\varGamma =1 ,\mathcal {M}=1$). For such a film, there is no in-plane heat conduction: the temperature 
$T_{f}$ is a function of time 
$t$ only; and models (F), (1D) and (A) all agree.
Figure 2(a) plots average film temperature against 
$h_0$ and time, showing that temperature depends on film thickness in a non-monotonic manner. Figure 2(b) plots the change of temperature with film thickness, 
$\partial T_{f}/\partial h$, evaluated at 
$h=1$ (this value of 
$h$ will be used in later simulations), as a function of time. For early times (
$t<2.95$), 
$\partial T_{f}/\partial h<0$, so that a decrease in film thickness corresponds to an increase in temperature (thinner film is hotter). For later times 
$\partial T_{f}/\partial h>0$ so that a decrease in film thickness leads to a decrease in temperature (thinner film is colder). This non-monotonic behaviour of 
$\partial T_{f}/\partial h$ in time is due to the changing balance between the heating from the source and the heat loss through the substrate (see appendix D for more details).
Figure 2. (a) Contour plots of average film temperature for a static flat film. (b) Rate of change of temperature with film thickness, 
$\partial T_{f}/\partial h$ as a function of time for 
$h=1$. At early times (
$t<2.95$) 
$\partial T_{f}/\partial h<0$ and later (
$t>2.95$) 
$\partial T_{f}/\partial h>0$.
We next consider evolving films, with the initial film profile given by (3.5). Figure 3(a) shows the evolution using model (F), though the behaviour is also representative of (A) and (1D). Melting temperature is reached at 
$t\approx 0.54$; by time 
$t=3.32$ the liquid film begins to evolve appreciably; at 
$t=4.52$ significant film evolution has occurred and the film first reaches the equilibrium film thickness; and at 
$t=6.50$ the film has fully dewetted. Figure 3(b) shows average film temperatures using models (F), (1D) and (A), as well as the film height at the midpoint 
$x=0$. The average film temperatures are in good agreement before dewetting, but afterwards model (1D) begins to deviate significantly from models (F) and (A), which show excellent agreement for the entirety of the simulation. The cooling rate, 
$\mathrm {d}T_{avg}/\mathrm {d}t$, is faster for (1D) than for (F) and (A) since heat cannot diffuse laterally through the film in (1D). This, in turn, produces a film that solidifies sooner (this will be discussed further in § 3.4). Despite this difference, the midpoint film height 
$h(x=0)$, shown here for (F), is similar for all models (see also figure 3a).
Figure 3. (a) Evolution of film thickness when material parameters are fixed and 
$h(x,0)$ is given by (3.5), at a few representative time points. (b) Average film temperature (see (3.1)) and midpoint film thickness 
$h(0,t)$ for the film profiles given in (a). Deviation between the models appears after the film dewets. The material parameters are set to their melting temperature values, 
$\varGamma =1$ and 
$\mathcal {M}=1$.
We next discuss the spatial variation of temperature. Figure 4 shows free-surface temperatures (a–d) and film midpoint temperatures (at 
$x=0$; e–h) at the times displayed in figure 3(a). Since the film thins at its midpoint as 
$t$ increases (see figure 3b), the film domain (
$z>0$) shrinks from (e–h). The lateral spatial variation of temperature is seen to be much weaker for models (F) and (A) than for model (1D) (see (a–d)), due to the inclusion of lateral heat diffusion in (F) and (A); and the substrate/film midpoint temperatures are lower for (1D) than for (F)/(A) in panels (f–h). In particular, we find that the temperature predictions of models (F) and (A) differ by at most 
$0.01$%, whereas (F) and (1D) differ by as much as 30 %. For model (1D), the temperature is initially higher at the film midpoint (
$x=0$) than at the edges, a situation that is reversed at later times (e.g. figure 4b). Using model (1D), therefore, may lead to overestimated temperature gradients, such as in Trice et al. (Reference Trice, Thomas, Favazza, Sureshkumar and Kalyanaraman2007), which may significantly alter the evolution of the film.
Figure 4. (a–d): Free surface (
$z=h$) temperature profiles, 
$T_{f}\vert _{z=h}$ and (e–h) midpoint (
$x=0$) temperature corresponding to times (a,e) 
$t=0.54$, (b,f) 
$t=3.32$, (a,g) 
$t=4.52$ and (d,h) 
$t=6.50$. Here, 
$\varGamma =1$, 
$\mathcal {M}=1$ and the overlapping curves correspond to models (F) and (A). Note the difference in vertical axis scales between panels (a–h). Colour code: (F) (red), (A) (green dashed), (1D) (blue).
To emphasize the lateral temperature variation for all models, figure 5 plots the difference between the free-surface temperatures at the thinnest (
$x=0$) and thickest (
$x=\pm {\rm \pi}$) parts of the film. Consistent with figure 4(a), this temperature difference for model (1D) is much larger than for models (F) and (A). All models show the same trend: initially the thinnest part of the film is hottest, but ultimately the thickest parts are hottest. We attribute this change of behaviour, which occurs relatively early in the film evolution, to the combination of lateral diffusion and the heat loss through the substrate (see (2.33)).
Figure 5. (a) Difference between free-surface temperatures at the thickest part of the film 
$(x=\pm {\rm \pi}$) and the thinnest part 
$(x=0)$ for models (F) (red), (A) (green) and (1D) (blue). (b) Zoom-in of results for models (F) and (A) from (a) to illustrate behaviour more clearly. The black dashed line represents the horizontal axis where 
$T(x={\rm \pi} )=T(x=0)$.
To conclude this section, we summarize our main findings. We have compared models (F), (1D) and (A) and found that (A) provides a much better approximation to (F) than does (1D). After dewetting, the film cools more rapidly with model (1D) than with (A) and (F) due to the neglect of in-plane thermal diffusion. Consequently, the average film temperatures in model (1D) vary significantly from those predicted by models (F) and (A). From a computational point of view, since the film temperature is independent of 
$z$ in model (A), it is significantly more computationally efficient than (F) and even more efficient than (1D). For illustration, we note that for a 
$255 \times 200$ computational grid in 
$x$ and 
$z$, the simulation times for a typical run reported in this section are 
$73.4$, 
$1.3$, 
$4.8$ hours for models (F), (A) and (1D) respectively, on a reasonably fast workstation.
3.4. Variation of material parameters
In this section, we consider the effect of varying surface tension and viscosity with temperature (see (2.48) and (2.50)) and consider the influence of the Marangoni effect, by comparing film evolution with 
${Ma}=0$ and 
${Ma}\neq 0$. In the previous section (as well as in additional tests, not reported here for brevity), we have demonstrated that model (A) provides a good approximation to the full model (F) with considerably less computational effort, and so henceforth, we will focus on exploring the differences between models (A) and (1D). In this section, domain lengths of 
$\varLambda _{m}$ and 
$20\varLambda _{m}$ are simulated, beginning with the former (domain length of 
$\varLambda _{m}$ may be assumed until otherwise stated).
We focus first on the case where both surface tension and viscosity depend on average temperature and are therefore time dependent, 
$\varGamma (t),\mathcal {M}(t)$, but we ignore the Marangoni effect (set 
${Ma} = 0$ in (2.55)). In the subsequent text, any reference to time-dependent surface tension or viscosity, 
$\varGamma (t),\mathcal {M}(t)$, refers solely to time dependence through the average temperature. Figure 6(a) shows the film thickness profiles and figure 6(b) shows the free-surface temperature profiles, 
$T_{f}\vert _{z=h}$, at the thin (
${\bigstar}$) and thicker (
$\blacksquare$) parts of the film. We observe that model predictions differ only after dewetting (B) and therefore the film thickness profiles in (1D) and (A) remain nearly identical. After dewetting, the two models exhibit marked differences in the temperature profiles. The large difference in cooling rates between (A) and (1D) is consistent with figure 3(b). For (1D) this leads to the thin part of the film 
${\bigstar}$ being significantly colder than the thicker parts 
$\blacksquare$ and the difference is exacerbated by the retention of heat in the thick part of the film.
Figure 6. (a) Film thickness profiles (simulated with (F), but representative of (A) and (1D) also) for times A (melting), B (dewetting), C (1D) re-solidification and D (A) re-solidification. The markers in (a) represent 
$x=0$ (
${\bigstar}$) and 
$x={\rm \pi}$ (
$\blacksquare$). (b) Free-surface temperature at 
${\bigstar}$ and 
$\blacksquare$ for (A) and (1D). The temperature profiles agree until the film dewets (B). Then, (1D) temperatures vary significantly at 
${\bigstar}$ and 
$\blacksquare$, whereas (A) produces similar temperatures at both locations. Surface tension and viscosity vary in time, but the Marangoni effect is not included (
${Ma}=0$).
Figure 7 shows the average film temperatures for models (A) and (1D) with 
$\varGamma (t),\mathcal {M}(t)$ as considered in figure 6, for both 
${Ma}=0$ and 
${Ma}\neq 0$. The Marangoni effect seems to have only a very minor influence on 
$T_{avg}$ for both models. Furthermore, the rapid cooling seen in model (1D) leads to a significantly shorter LL (
${LL}=4.73$) than that predicted by model (A) (
${LL}=5.94$).
Figure 7. Average film temperatures using models (A) and (1D) for 
${Ma}=0$ and 
${Ma} \neq 0$, and using 
$\varGamma (t)$ and 
$\mathcal {M}(t)$. Points A–D correspond to those of figure 6(a). Near 
$t=0.54$ (A), the film melts in both models. The models begin to deviate around 
$t=1.72$ (B) when dewetting occurs; from this point until solidification (which for 
${Ma}=0$ occurs at C for (1D), D for (A); the times for the 
${Ma}\neq 0$ cases are similar) the temperature in model (1D) is lower than that in (A). The LLs are 
${LL}\approx 4.73$ (1D) and 
${LL}\approx 5.94$ (A).
Figure 8 compares the results obtained with(out) temperature variation of material properties. From figure 8(a), we immediately conclude that the Marangoni effect is very weak, and will be ignored henceforth (
${Ma=0}$ for the remaining results). The weak Marangoni effect is, at least in part, due to the high thermal conductivity of the film, which sets the thermal time scale and gives rise to the very weak spatial variations in interfacial temperatures seen in figure 4. The second observation is that allowing surface tension to depend on time, 
$\varGamma =\varGamma (t)$, has a small but measurable effect on the results. With this dependence included, the film instability appears to develop faster than in the constant-
$\varGamma$ case (compare the green dashed curve with the red solid curve in figure 8a). The third observation is that the time dependence of viscosity has by far the largest effect on the film instability development, leading to much faster dewetting.
Figure 8. (a) Midpoint film thickness, 
$h(0,t)$, for the cases: surface tension and viscosity fixed, 
$\varGamma =1,\mathcal {M}=1$, with no Marangoni effect, 
${Ma}=0$ (red, solid line); surface tension and viscosity fixed with Marangoni effect, 
${Ma}\neq 0$ (black, dash-dotted line); surface tension varies in time, 
$\varGamma (t)$, with viscosity fixed and no Marangoni (green, dashed line); viscosity varies in time, 
$\mathcal {M}(t)$, surface tension fixed and no Marangoni (blue, dot-dashed line). In all cases, model (A) was used to calculate temperature. (b) Growth rate, 
$\beta$ as a function of wavenumber, 
$k$, using (3.3) with 
$\varGamma =1$ and 
$\mathcal {M}=1$ and no Marangoni effect, 
${Ma}=0$. The blue, green and red dots represent maximum growth rate extracted from the corresponding simulations in (a).
We now revisit the predictions of LSA and compare them to simulation results. Figure 8(b) plots the dispersion curve according to (3.3), for the constant parameter case 
$\varGamma =1$ and 
$\mathcal {M}=1$ (with 
${Ma}=0$). To estimate the maximum growth rate, 
$\beta _{m}$, in our numerical simulations, we assume the film perturbation grows exponentially at early times, consistent with LSA: 
$h=h_0(1+\delta \exp ({\textrm {i}k_{m}x+\beta _{m} t}))$, where 
$k_{m}$ is the corresponding wavenumber (although in practice perturbations of many different wavenumbers exist, the one usually most apparent in the unstable regime is 
$k_{m}$). We then perform a best linear fit of 
$\ln ((h(0,t)-h_0)/(\delta h_0))$ vs 
$t$ for early times. The red, green and blue dots correspond to growth rates 
$\beta _{m}$ extracted from the corresponding colour-coded simulations in figure 8(a). The film with parameters fixed (red) grows at the rate predicted by LSA, whereas the film with time-dependent surface tension (green) grows at a slightly faster rate. The growth rate in the time-dependent viscosity case (blue) is similar to the time-dependent surface tension case, despite the much faster instability development in figure 8(a). This indicates the relevance of the nonlinear part of instability growth.
To highlight further the significance of time-dependent viscosity, figure 9 compares film evolution for the constant (a) and variable (b) viscosity cases. The difference in film evolution is significant, with much larger instability growth rate for time-dependent viscosity.
Figure 9. Comparison of free-surface evolution 
$h(x,t)$ when (a) viscosity is fixed at the melting value, 
$\mathcal {M}=1$, and (b) viscosity varies in time according to average temperature, 
$\mathcal {M}(t)$ (see (2.50)). In both cases, surface tension is fixed at the melting value 
$\varGamma =1$, and all times plotted are prior to film re-solidification. In this and in the figures that follow, 
${Ma}=0$. In (a) we find that re-solidification happens before dewetting (note the vertical axis scale), while the converse is true for (b), indicating the importance of variable viscosity. The times are as follows: 
$t=0$ (red), 
$t=1.47$ (blue), 
$t=1.75$ (orange), 
$t=1.88$ (black dashed), 
$t=1.92$ (light blue), 
$t=2.21$ (green).
To summarize, we have simulated and compared the results with and without the Marangoni effect, with and without time-varying surface tension, and with and without time-varying viscosity. We have found that allowing either surface tension or viscosity to depend on time through the average temperature speeds up the dewetting mechanism. In particular, varying viscosity has the strongest impact on the film dynamics.
For completeness we also consider the possible variation of viscosity with 
$x$. For model (A), since 
$T_{f}^{(0)}(x,t)$ is a function only of the in-plane spatial variable 
$x$ (and time 
$t$), we may replace 
$\bar {T}$ in the definition of 
$\mathcal {M}$ (2.50) by 
$T_{f}^{(0)}$, so that viscosity, which we denote in this case by 
$\mathcal {M}(x,t)$, depends on both space and time, namely
\begin{equation} \mathcal{M} (x,t) = \exp\left( \frac{E}{RT_{melt}} \left( \frac{1}{T_{f}^{(0)}(x,t)}-1 \right) \right). \end{equation}
For model (1D), on the other hand, the calculated film temperature 
$T_{f} (x,z,t)$ depends on both spatial coordinates 
$x$ and 
$z$, so obtaining an analogous model for 
$\mathcal {M}(x,t)$ requires some additional assumptions. We use the free-surface temperature 
$T_{f}\vert _{z=h}$, so that viscosity takes the form
\begin{equation} \mathcal{M} (x,t) = \exp\left( \frac{E}{RT_{melt}} \left( \frac{1}{\left(T_{f}\vert_{z=h}\right)}-1 \right) \right). \end{equation} The film thickness profiles produced in these 
$\mathcal {M}(x,t)$ cases (both (1D) and (A)) are similar to the 
$\mathcal {M}(t)$ cases shown above and are thus omitted from the main text for brevity (an example is shown in appendix B). Figure 10 plots the average film temperature profiles for each of the models; we observe that the temperatures agree up until point A, at which the film dewets for variable viscosity. The inset shows that, with model (A), the variable viscosity cases (
$\mathcal {M}(t), \mathcal {M}(x,t)$) lead to a slightly longer LL than if viscosity is constant, 
$\mathcal {M}=1$. This may be important considering that the films dewet much nearer the resolidification time when 
$\mathcal {M}=1$. This difference, however, is always small relative to that between models (A) and (1D). When viscosity varies with time and/or space, model (1D) predicts LLs that are far shorter. Furthermore, for both models the results for 
$\mathcal {M}(x,t)$ and 
$\mathcal {M}(t)$ are identical. This is, we believe, due to the film temperatures not deviating far from the respective average temperatures (even for model (1D) where 
$x$-variation of temperature is much larger than for (A), the temperature deviates from the average by at most 
$30\,\%$, which appears to be insufficient to cause significant differences between results for 
$\mathcal {M}(t)$ and 
$\mathcal {M}(x,t)$). We conclude, therefore, that to simulate cases where film dewetting occurs much faster than resolidification, the dependence of viscosity on average temperature is the most important effect to include.
Figure 10. Average film temperatures with model (A): 
$\varGamma =1,\mathcal {M}=1$ (red solid line); 
$\varGamma (t), \mathcal {M}=1$ (green dashed line); 
$\varGamma =1, \mathcal {M}(t)$ (cyan dash-dotted line); 
$\varGamma =1$, spatio-temporally varying viscosity 
$\mathcal {M}(x,t)$ (blue dashed double-dotted line). Average film temperatures with model (1D): 
$\varGamma =1, \mathcal {M}(t)$, (magenta dash-dotted line) and 
$\varGamma =1, \mathcal {M}(x,t)$, (black dash-dotted line). In all cases, the domain length is 
$\varLambda _{m}$ from LSA. Point 
${A}$ marks the film dewetting time for the time-varying viscosity 
$\mathcal {M}(t)$ case. Inset: zoomed-in image of solidification point for the model (A) cases. Here, including 
$\mathcal {M}(t)$ or 
$\mathcal {M}(x,t)$ leads to a slightly longer LL than the constant viscosity (
$\mathcal {M}=1$) cases. Contrast with model (1D) where variable viscosity produces a LL that differs significantly from the rest. Simulations are marked 
$1$–
$6$ and placed near the curves as a guide. Simulations 
$1$–
$4$ were simulated with (A) and 
$5$–
$6$ with (1D).
To emphasize the importance of accurately modelling the film viscosity, we present simulation results for domains of length 
$20\varLambda _{m}$, with films subjected to initial random perturbations according to (3.6). Figure 11 plots film thickness for (a) 
$\mathcal {M}=1$ and (b) 
$\mathcal {M}(t)$, both using model (A) with 
$\varGamma =1$ and 
${Ma}=0$. We again see that films with time-dependent viscosity (figure 11b) dewet much faster than those with constant viscosity (figure 11a). The black curve represents the solidification time of (a) (nearly the same as that for (b)) and hence represents the final configurations of both simulations. The main finding is that the 
$\mathcal {M}=1$ case does not fully dewet, whereas the 
$\mathcal {M}(t)$ case does. This difference in dewetting time scales is consistent with the earlier results for domains of length 
$\varLambda _{m}$, except that there the film was completely dewetted at solidification for both cases. We note slight coarsening in figure 11(a,b) with some droplets merging, leading to fewer than 20 resultant (primary) droplets. The film height evolution for model (1D), cases 
$\mathcal {M}(t), \mathcal {M}(x,t)$, and model (A) case 
$\mathcal {M}(x,t)$ is identical to figure 11(b). In appendix B we show a case where models (1D) and (A) lead to diverging results (for a different choice of parameters). We omit the temperature profiles here as the short LL for (1D) and identical 
$\mathcal {M}(t)$ and 
$\mathcal {M}(x,t)$ results are consistent with the results on domain length 
$\varLambda _{m}$.
Figure 11. Evolution of films subject to initial random perturbations for (a) constant viscosity, 
$\mathcal {M}=1$, and (b) time-dependent viscosity, 
$\mathcal {M}(t)$. Colour code: red 
$t=0$; blue 
$t=1.66$; green 
$t=2.30$; orange 
$t=3.31$; black 
$t=6.47$ (re-solidification time for (a)). The domain length was taken to be 
$20\varLambda _{m}$ and all simulations were done using model (A) with 
$\varGamma =1$. The 
$h$-axis is plotted on a log scale to emphasize satellite droplet formation.
To summarize the results on long domains of length 
$20\varLambda _{m}$, we have found that including the time dependence of viscosity permits the films to fully dewet in the span of the LL whereas keeping viscosity fixed at the melting temperature value (
$\mathcal {M}=1)$ produces films that dewet only partially during the liquid phase. Any effects due to spatial variation of viscosity appear to be irrelevant here. Finally, the LL is much shorter for model (1D) than model (A), as expected.
4. Conclusions
To conclude, we have formulated three models for heat conduction in nanoscale thin films on thermally conductive substrates: a full model (F) that accounts for heat conduction in all directions in both film and substrate; an asymptotically reduced model (A) that exploits a disparity in length scales in both film and substrate to derive an equation governing in-plane diffusion of heat within the film coupled to out-of-plane heat diffusion in the substrate; and a 1-D model (1D) that simply neglects any in-plane diffusion in both film and substrate. In all cases, a thin film model is used to describe the associated fluid dynamics. The main finding is that including in-plane diffusion in the thermal modelling influences strongly the film evolution. In particular, neglecting in-plane diffusion is found to amplify (artificially) in-plane thermal gradients and expedite film cooling. We have found that model (A) is significantly more accurate than (1D) while being considerably more computationally efficient than (F). We have also found that when material parameters are allowed to vary in time through the average film temperature, model (A) produces LLs significantly longer than those of model (1D), due to the absence of lateral heat conduction in (1D). Therefore model (A) combines both accuracy and efficiency.
With regard to the individual (dimensionless) material parameters that arise in our models, 
${Ma}$ (Marangoni number), 
$\varGamma$ (surface tension parameter) and 
$\mathcal {M}$ (film viscosity), we find that the variation of viscosity with time has the greatest effect on model outcomes. By including time-dependent viscosity, films exposed to laser heating (on both small and large domains) fully dewet while in the molten state. In contrast, when viscosity is held constant, dewetting occurs much later in the cooling process, which may result in partial droplet formation only. This suggests strongly that time-dependent viscosity is needed to represent accurately experiment-like behaviour. Using a spatio-temporally varying viscosity, 
${\mathcal {M}}(x,t)$, produces essentially identical results to the case where viscosity depends only on time. Introducing time dependence of the surface tension (
$\varGamma (t)$) has a larger effect on the film instability growth rate (increasing it) than does the Marangoni effect (
${Ma}\neq 0$), though the effect is always small, and insignificant when compared to the variation of viscosity with time. The Marangoni effect was found to be negligible in all cases considered.
Although model (A) is found to be useful in the current setting, its validity relies on a number of underlying assumptions. Therefore, its applicability to other problems must be carefully verified prior to use. In this study we have considered the time and space variation of only selected physical parameters (surface tension and viscosity) through temperature, and have not considered how temperature dependence of other material parameters, such as thermal conductivity and density, may influence the results. Furthermore, all simulations presented here are restricted to the two-dimensional geometry: much more significant computational benefit of model (A) is expected in three spatial dimensions. Extending our model and simulations to three dimensions is one of the directions of our future work.
Acknowledgements
The authors are grateful to I. Seric, W. Batson, M. Lam, P. Sanaei and S. Afkhami for fruitful discussions.
Funding
This research was supported by NSF DMS-1815613 (L.J.C., L.K.); by NSF DMS-1615719 (L.J.C., L.K.); by NSF CBET-1604351 (R.H.A., L.K.); and by CNMS2020-A-00110 (L.K.).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Scalings and parameter values
A.1. Discussion of the choice of scales
The choice of scales has important implications for the derivation of both the thin film equation (2.55) and the thermal model (A) (2.41)–(2.46). The choice of time scale is typically based on the fluid flow, 
$L/U$. From the perspective of the thermal model, however, the pulsed laser heating duration is of the order of nanoseconds. With the time scale choice 
$L/U$, 
$L$ and 
$U$ should be chosen consistent with such a thermal time scale, while still retaining surface tension effects in the fluid flow model, known to be important. Therefore, we choose 
$U=\epsilon ^3 \gamma _{m}/(3\mu _{m})$ and scale 
$L$ on the (inverse) wavenumber of maximum growth 
$k_{m}= 2{\rm \pi} /\lambda _{m}$ as 
$L=k_{m}^{-1}$ (from § 3.1) so that surface tension effects and disjoining pressure are retained to leading order. The tradeoff is that the aspect ratio 
$\varepsilon$, consistent with the data (table 1), is rather large, 0.246. However, the time scale is then of the order of nanoseconds, as desired, and we consider this acceptable in order to develop a consistent model.
For the materials in question, the values of 
${Pe}_{f}, {Pe}_{s}, \mathcal {K}$ and 
${Bi}$ are small despite the 
$O(1)$ assumption (table 1). This is primarily a consequence of the size of 
$\epsilon$, which is directly related to the dependence of 
$L$ on 
$\lambda _{m}$. Given this observation, we briefly consider the limit of small 
${Pe}_{f}, {Pe}_{s}, \mathcal {K}$ and 
${Bi}$. Firstly, in the limit 
${Pe}_{f}\to 0$, (2.41) would reduce to a quasi-steady (no time derivative) equation governing temperature 
$T_{f}(x,y,t)$ in the film. The resultant equation is computationally more difficult to solve (and leads to only negligible differences), so we do not adopt this approach. Secondly, in the limit 
${Pe}_{s} \to 0$, the solution to (2.42) would be linear in 
$z$. The numerical solutions given in figure 4 display substrate temperatures that deviate from linear behaviour in 
$z$. We consider this the motivation for retaining 
${Pe}_{s}$ in (2.42) (which leads to better agreement between models (A) and (F)). In the limit 
$\mathcal {K}\to 0$, film temperature no longer depends directly on substrate temperature (although the substrate temperature still depends on the film temperature). Since the primary heat loss mechanism is considered to be through the film–substrate interface, we retain 
$\mathcal {K}$ in (2.41). If the primary heat loss mechanism is elsewhere, it may be possible to drop the term containing 
$\mathcal {K}$ in (2.41). Finally, in the limit 
${Bi} \to 0$, (2.44) becomes an insulating boundary condition, which in turn leads to much higher substrate/film temperatures. In the case where SiO
$_2$ sits on a native layer of Si, it is expected that there is some heat transfer. Therefore, we retain 
${Bi}$ in (2.44).
A.2. Wavelength of maximum growth
When the material parameters are fixed at melting temperature, 
$\varGamma ,\mathcal {M}=1$, the dimensional dispersion relation can be written as
\begin{gather} \beta(k) = \frac{h_0^2 H^2 k^2}{3\mu_{f}} \left( P_0 - \gamma_{f} H h_0 k^2 \right), \end{gather}
\begin{gather}P_0 = \frac{A_{H}}{6{\rm \pi} h_*^3 H^3} \left[ m\left(\frac{h_*}{h_0} \right)^m - n \left( \frac{h_*}{h_0}\right)^n \right], \end{gather}
where (in this section only) 
$\beta (k)$ is the dimensional growth rate and 
$k$ the dimensional wavenumber. The remaining parameters are given in table 1. The wavelength of maximum growth, 
$\lambda _{m}$, can then be found by setting 
$\partial \beta /\partial k=0$ and written as
\begin{equation} \lambda_{m} = 2{\rm \pi}/\sqrt{\frac{A_{H}}{12{\rm \pi} \gamma_{f} h_*^3 H^3 h_0} \left[ m \left( \frac{h_*}{h_0} \right)^m - n \left( \frac{h_*}{h_0}\right)^n \right]}. \end{equation}
For the parameters used in this paper 
$\lambda _{m}=255\ \textrm {nm}$. For the entirety of the paper, this value is used to define a base length scale 
$L=\lambda _{m}/(2{\rm \pi} )$.
Appendix B. Effect of spatially varying viscosity with a larger Biot number
Here we increase the Biot number to 
${Bi}=5.71\times 10^{-3}$ (more than twice the value used in the main text; see table 1). This leads to much earlier resolidification of the film (for (A) resolidification occurs at 
$t\approx 2.8$ whereas for (1D) 
$t\approx 2.7$). We specifically focus on the influence of model choice when the spatio-temporally varying viscosity is used, 
$\mathcal {M}(x,t)$. Figure 12 shows the film evolution for both models (A) (a) and (1D) (b). The black solid line represents the film profiles at the (1D) solidification time (model (A) predicts a larger solidification time but the final solidified film configuration is nearly identical to the black solid line in figure 12a). The times are given in the caption. The main finding is that for sufficiently large 
$\alpha$ the decay of the satellite droplet that forms for both (A) and (1D) is slower in (1D). Essentially, for (1D) the satellite droplet is cold relative to the thicker parts of the film (recall the much more significant 
$x$-variation of temperature observed for (1D)). Since the viscosity varies with space, the centre (
$x=0$) of the satellite droplet approaches the precursor thickness more slowly than the surrounding area. In (A), the temperature variation with 
$x$ is less pronounced and so the satellite droplet has drained more than its (1D) counterpart before reaching the final solidified configuration. In summary, these results demonstrate that the choice of thermal model can influence the final resulting film profiles.
Figure 12. Film thickness evolution for both the asymptotic model (a) and the 1-D model (b) for the spatio-temporally varying viscosity case 
$\mathcal {M}(x,t)$. The 
$y$-axis is plotted on a log scale to emphasize satellite droplet formation which is more prominent with (1D) than with (A). Here the Biot number is 
${Bi}=5.71\times 10^{-3}$. The times are: 
$t=0$ (red), 
$t=1.47$ (green), 
$t=1.84$ (blue), 
$t=2.21$ (magenta), 
$t=2.40$ (orange dashed), resolidification of (1D), 
$t=2.69$ (black).
Appendix C. Numerical schemes
C.1. Numerical solution of model (F)
To solve (2.5) numerically, along with corresponding boundary conditions in model (F), we define the new variables 
$(u,v,\tau )$ as
\begin{equation} u=x,\quad v= \frac{zh_0}{h},\quad \tau=t, \end{equation}which is a time-dependent mapping transforming the deformable domain, describing the film, into a fixed rectangle (see figure 13).
Figure 13. Schematic of mapping used to transform the computational domain for (F). Here, the domain width shown is 
$X=2{\rm \pi}$. Blue points represent the film domain and black points represent the substrate domain.
This approach trades a simplified domain for an increase in complexity of the thermal equation in the film. Equation (2.5) is transformed into the following differential equation:
\begin{equation} \partial_{\tau} T_{f} = \frac{1}{{Pe}_{f}} \partial_u^2 T_{f} + B(u,v) \partial_{uv} T_{f} + A(u,v) \partial_v^2 T_{f} + F(u,v) \partial_v T_{f} + Q(h(u),v,\tau), \end{equation}where
\begin{equation} \left.\begin{gathered} B(u,v)={-}\frac{1}{{Pe}_{f}}\left( \frac{2v \partial_u h}{h} \right), \\ A(u,v)= \frac{1}{{Pe}_{f}} \left[ \left( \frac{v\partial_u h}{h} \right)^2 + \frac{1}{\epsilon^2} \left( \frac{h_0}{h} \right)^2 \right], \\ F(u,v)= \frac{1}{{Pe}_{f}}\left( \frac{vh\partial_u^2 h}{h^2} - \frac{v (\partial_u h)^2}{h^2} \right) + \frac{v \partial_{\tau} h}{h}, \\ Q(h(u),v,\tau) = \frac{F(t)}{{Pe}_{f}} \left[ 1-R(h) \right] \exp{\left[-\alpha_{f}\left(h-hv/h_0 \right)\right]}, \end{gathered}\right\} \end{equation}
are the coefficients, and subscripts 
$u,v,\tau$ denote partial derivatives. To solve (C2) numerically, the ADI method is used, with the term containing 
$B(u,v)$ treated explicitly. A Crank–Nicolson scheme is used to solve (2.6), the heat equation in the substrate. We use a cell-centred grid system
\begin{equation} \left.\begin{gathered} u_i = x_i = x_0 + {\rm \Delta} x \left( i - 1/2 \right),\quad i=1,\ldots, n,\qquad {\rm \Delta} u = {\rm \Delta} x = \frac{\left(x_{max} - x_0 \right)}{n}, \\ z_j ={-}\left(j-1/2 \right) {\rm \Delta} z_s,\quad j=1,\ldots, n_{s},\qquad {\rm \Delta} z_s = \frac{H_{s}}{n_{s} }, \\ v_j = h_0 - \left( j-1 \right) {\rm \Delta} v,\quad j=1,\ldots, n_{f},\qquad {\rm \Delta} v = \frac{h_0}{n_{f}}, \\ T_{i,j}^k \approx \begin{cases} T_{f}(u_i, v_j, t_k), & 1 < j < n_{f}, \\ T_{s}(x_i, z_j, t_k), & n_{f}+1 < j < n_{f} + n_{s}, \end{cases} \\ h_i^k \approx h(u_i, t_k),\quad Q_{i,j}^k \approx Q(u_i, v_j, t_k). \end{gathered}\right\} \end{equation}The numerical system can be then written as
\begin{equation}
\left.\begin{gathered} \frac{T_{i,j}^* - T_{i,j}^k}{{\rm \Delta}
t} = D_{i,j}^* + G_{i,j}^k + M_{i,j}^k + \frac{1}{2}
Q_{i,j}^{k+1/2}, \\ \frac{T_{i,j}^{k+1} -
T_{i,j}^*}{{\rm \Delta} t} = D_{i,j}^* + G_{i,j}^{k+1} +
\frac{1}{2} Q_{i,j}^{k+1/2}, \\ D_{i,j}^* = \begin{cases}
\tfrac{1}{2}{Pe}_{f}^{{-}1} \delta_u^2 T_{i,j}^*, & j\leq
n_{f}, \\ \tfrac{1}{2} \epsilon^2 {Pe}_{s}^{{-}1}
\delta_x^2 T_{i,j}^*, & j\geq n_{f}+1, \end{cases}
\\ G_{i,j}^k = \begin{cases} \tfrac{1}{2} A_{i,j}^k
\delta_v^2 T_{i,j}^k + \tfrac{1}{2} F_{i,j}^k \delta_v
T_{i,j}^k, & j\leq n_{f}, \\ \tfrac{1}{2}
{Pe}_{s}^{{-}1}\delta_z^2 T_{i,j}^k, & j\geq n_{f}+1,
\end{cases} \\ M_{i,j}^k = \begin{cases} B_{i,j}^k
\delta_{uv} T_{i,j}^k, & j\leq n_{f}, \\ 0, & j>n_{f},
\end{cases} \end{gathered}\right\}
\end{equation}
where 
$T_{f}(x_i,z_j,t_k)\approx T_{i,j}^{k}$ is a discretization of the film temperature, and 
$T_{i,j}^{*}$ represents the solution at an intermediate time step. In the interior grid, the spatial derivatives are given by
\begin{equation} \left.\begin{gathered} \delta_u T_{i,j}=\frac{T_{i+1,j}-T_{i-1,j}}{2{\rm \Delta} u},\quad \delta_v T_{i,j}=\frac{T_{i,j-1}-T_{i,j+1}}{2{\rm \Delta} v}, \\ \delta_x^2 T_{i,j} = \delta_u^2 T_{i,j}= \frac{T_{i+1,j}-2T_{i,j}+T_{i-1,j}}{{\rm \Delta} u^2}, \\ \delta_{uv} T_{i,j} = \frac{T_{i+1,j-1}-T_{i-1,j-1}-T_{i+1,j+1}+T_{i-1,j+1}}{4{\rm \Delta} u {\rm \Delta} v}, \\ \delta_z^2 T_{i,j} = \frac{T_{i,j-1}-2T_{i,j} + T_{i,j+1}}{{\rm \Delta} z_s^2},\quad \delta_v^2 T_{i,j}=\frac{T_{i,j+1}-2T_{i,j}+T_{i,j-1}}{{\rm \Delta} v^2}, \end{gathered}\right\} \end{equation}
which are second-order central difference approximations and the source, 
$Q$, is approximated by an average at times 
$t_k$, and 
$t_{k+1}$,
\begin{equation}
Q_{i,j}^{k+1/2} = \tfrac{1}{2} ( Q_{i,j}^{k+1} +
Q_{i,j}^k ).
\end{equation}C.2. Numerical solution of model (A)
We use the same cell-centred grid as § C.1 except for 
$z_j$
\begin{equation} z_j ={-}\left(j-1 \right) {\rm \Delta} z_{s},\quad j=1,\ldots, n_{s},\qquad {\rm \Delta} z_{s} = \frac{H_{s}}{n_{s} - 1/2}, \end{equation}
so that a grid point exists on the liquid–solid interface (a Dirichlet boundary condition is prescribed there). For simplicity, we let 
$T_i^{k} \approx T_{f}(x_i, t_k)$, 
$S_{i,j}^{k} \approx T_{s}(x_i,z_j,t_k)$ and
\begin{equation} \left.\begin{gathered} F_i^{k} = \frac{1}{{Pe}_{f}} \left[ \delta_x^2 T_i^{k} + \left(\frac{\partial_x h}{h}\right)_{i}^{k} \delta_x T_{i}^{k} \right], \\ G_i^{k} ={-}\frac{\mathcal{K}}{{Pe}_{f} h_i^k} \partial_z \left( S \right)_I^{k}, \end{gathered}\right\} \end{equation}
where 
$\partial _z (S)_{I}^k=\partial T_{s}/\partial z \vert _{z=0}(t=t_k)$. To solve (2.41) and (2.42) we use a predictor–corrector Runge–Kutta/Crank–Nicolson scheme. In the prediction phase, we use a forward-Euler scheme to deal with 
$G_i^k$
\begin{equation}
\left.\begin{gathered} \frac{T_{i}^{k+1}-T_{i}^k}{{\rm \Delta}
t}=\frac{1}{2} [ F_i^{k+1} + F_i^{k} ] +
G_i^{k} +
\bar{Q}_{i}^{k+1/2},\quad i= 1, \ldots, n, \\
\frac{S_{i,j}^{k+1} - S_{i,j}^{k}}{{\rm \Delta} t} = \frac{1}{2}
{Pe}_{s}^{{-}1} [ \delta_z^2 S_{i,j}^{k+1} +
\delta_z^2 S_{i,j}^{k} ],\quad j= 2, \ldots,
n_{s},\enspace i = 1, \ldots, n,
\end{gathered}\right\}
\end{equation}
where 
$\bar {Q}_i^{k+1/2}=(\bar {Q}_i^k+\bar {Q}_i^{k+1})/2$. We then correct this prediction by using a Runge–Kutta, order-2, method on 
$G_i^k$ using the prediction 
$\hat {G}_i^k$
\begin{equation}
\left.\begin{gathered} \frac{T_{i}^{k+1}-T_{i}^k}{{\rm \Delta}
t}=\frac{1}{2} [ F_i^{k+1} + F_i^{k} ] +
\frac{1}{2} ( G_i^{k} + \hat{G}_i^{k} ) +
\bar{Q}_{i}^{k+1/2},\quad i = 1, \ldots, n, \\
\frac{S_{i,j}^{k+1} - S_{i,j}^{k}}{{\rm \Delta} t} = \frac{1
}{2} {Pe}_{s}^{{-}1} [ \delta_z^2 S_{i,j}^{k+1} +
\delta_z^2 S_{i,j}^{k} ],\quad j = 2, \ldots,
n_{s},\enspace i = 1, \ldots, n.
\end{gathered}\right\}
\end{equation}C.3. Numerical solution of model (1D)
The numerical scheme used for model (1D) is a simple Crank–Nicolson scheme
\begin{equation}
\left.\begin{gathered} \frac{T_{i,j}^{k+1} -
T_{i,j}^k}{{\rm \Delta} t} =
\frac{1}{2}\epsilon^{{-}2}{Pe}_{f}^{{-}1} ( \delta_z^2
T_{i,j}^k + \delta_z^2 T_{i,j}^{k+1} ) +
Q_{i,j}^{k+1/2},\quad j=1,\ldots,n_{f}, \\
\frac{T_{i,j}^{k+1} - T_{i,j}^k}{{\rm \Delta} t} =
\frac{1}{2}{Pe}_{s}^{{-}1} ( \delta_z^2 T_{i,j}^k +
\delta_z^2 T_{i,j}^{k+1} ),\quad j=n_{f}+1,\ldots,
n_{f}+n_{s}, \end{gathered}\right\}
\end{equation}
where 
$i=1,\ldots , n$ and 
$n, n_{f}$ and 
$n_{s}$ are the same as in § C.1.
C.4. Convergence results
For what follows, we define the discrete 
$L2$ error, 
$E(t)$, where
\begin{equation}
E^2(t) = \frac{\displaystyle\sum_{i=1}^{n}
\left(\sum_{j=1}^{n_{f}}{\rm \Delta} z_i \left|
T^{comp}_{i,j}-T^{bench}_{i,j} \right|^2 + {\rm \Delta} z_{s}
\sum_{j=1}^{n_{s}} \left| T^{comp}_{i,j}-T^{bench}_{i,j}
\right|^2 \right) }{\displaystyle\sum_{i=1}^{n}
\left(\sum_{j=1}^{n_{f}}{\rm \Delta} z_i \left|T^{bench}_{i,j}
\right|^2 + {\rm \Delta} z_{s} \sum_{j=1}^{n_{s}} \left|
T^{bench}_{i,j} \right|^2 \right) },
\end{equation}
and 
$T_{i,j}^{comp}$ is the computed temperature and 
$T_{i,j}^{bench}$ is a benchmark solution which, for the results presented next, we take to be the numerical solution on the finest grid since no exact solution is known.
Figure 14 shows the quadratic convergence in 
$\Delta t$ (a) and 
$\Delta z_{s}$ (b) for each thermal model. Similarly, figure 15(a) shows quadratic 
$\Delta x$ convergence for (F) and (A), whereas figure 15(b) shows quadratic 
$\Delta v$ convergence for models (F) and (1D).
Figure 14. (a) The 
${\rm \Delta} t$ convergence for (F), (A) and (1D) in the stationary curved film case. All models use 
$O({\rm \Delta} t^2)$ schemes. (b) The 
${\rm \Delta} z_{s}$ convergence for (F), (A) and (1D) in the stationary curved film case where 
$h$ is given by (3.5), but time independent. All models use 
$O({\rm \Delta} z_{s}^2)$ schemes.
Figure 15. (a) The 
${\rm \Delta} x$ convergence for (F) and (A) in the stationary curved film case (note that (1D) has no derivatives with respect to 
$x$). Each model uses a 
$O({\rm \Delta} x^2)$ scheme. (b) The 
${\rm \Delta} v$ convergence for (F) and (1D) in the stationary curved film case (note that 
$T_{f}$ is independent of 
$z$ for (A)). Each model uses a 
$O({\rm \Delta} v^2)$ scheme.
Appendix D. Variation of temperature with film thickness
In the case when the film is flat, its temperature is independent of the in-plane variables and conservation of energy may be reduced to an expression for the average film temperature, written as a simple balance of source heating and substrate cooling
\begin{equation} \partial_t T_{avg} = {Pe}_{f}^{{-}1} \left[\frac{-\mathcal{K}}{h} \left(\partial_z T_{s}\right)\vert_{z=0} + \bar{Q}\right]. \end{equation} Figure 16 shows the variation of both (a) source heating (measured by 
${Pe}_{f}^{-1}\bar {Q}$) and (b) substrate cooling (measured by 
$(h {Pe}_{f} )^{-1} ( \partial _z T_{f})\vert _{z=0}$) for three different film thicknesses, 
$h=8,10,12$ nm. Figure 16(c) gives the average film temperatures for each case. As seen in figure 16(a), thinner films retain more energy from the source and, in the absence of cooling, should be hotter. Physically, thicker films reflect more energy and absorb less. Although (to keep the presentation simple) we use 
$\bar {Q}$ here rather than 
$Q$ from (2.18), it can be shown also that 
$\textrm {d} Q/\textrm {d} h<0$ for all times. Figure 16(b) demonstrates that thinner films also cool faster (through the substrate) than thicker ones. When these two effects are combined, we arrive at average film temperatures that are non-monotonic in film thickness 
$h$. In figure 16(c), thinner films are observed to be initially hotter over the early stages of evolution, primarily due the magnitude of the laser heating, which is initially larger than that of cooling. Over the later stages, the source decreases in strength sufficiently so that the trend of cooling with film thickness is dominant and thinner films are colder. Note that this non-monotonic behaviour depends on the relative strengths of the heat source and the cooling term and may change if different forms are used to describe these effects. The explanation given above is the primary reason for the behaviour seen in figure 2(b).
Figure 16. (a) Average source, 
$\bar {Q}$, (b) magnitude of film heat loss through the substrate, 
$\mathcal {K}(h {Pe}_{f})^{-1}\partial _z T_{s}\vert _{z=0}$ and (c) average film temperature, all for flat films of 
$h=8$ nm (blue, solid), 
$h=10$ nm (red, dot-dashed) and 
$h=12$ nm (black, dashed).
