1. Introduction
The dynamics of spreading and wetting of droplets on solid substrates have a spectrum of applications in industrial processes, such as painting and surface coatings on solid substrates, spraying, lubrication, ink-jet printing and many more (Daniel & Berg Reference Daniel and Berg2006; Bonn et al. Reference Bonn, Eggers, Indekeu, Meunier and Rolley2009). Furthermore, the recent advancements in the lab-on-a-chip technology have motivated researchers to find optimal strategies for effective manipulation of microdroplets in view of their importance in broad-ranging bio-microfluidic processes, specifically in drug delivery of ‘microbicides’ as anti-HIV agents (Hu & Kieweg Reference Hu and Kieweg2012), treatment of dry-eye syndrome caused by early onset of tear film rupture, analytical detection, digital microfluidics and encapsulation of biological cells (Stone, Stroock & Ajdari Reference Stone, Stroock and Ajdari2004; Braun Reference Braun2012; Anna Reference Anna2016), to name a few.
The spreading behaviour of sessile drops in the above-mentioned practical settings is highly influenced by the surface tension at the liquid–air interface. Surface tension manipulation using a temperature gradient at the solid substrate has been recognized as an effective control mechanism for the movement and wetting of droplets along the surface. This phenomenon is known as the thermocapillary or thermal Marangoni actuation of droplets (Young, Goldstein & Block Reference Young, Goldstein and Block1959; Subramanian & Balasubramaniam Reference Subramanian and Balasubramaniam2001). The complex interfacial dynamics of droplet spreading over non-isothermal substrates have prompted a number of experimental and theoretical investigations over the last few decades (Ehrhard & Davis Reference Ehrhard and Davis1991; Brzoska, Brochard-Wyart & Rondelez Reference Brzoska, Brochard-Wyart and Rondelez1993; Pratap, Moumen & Subramanian Reference Pratap, Moumen and Subramanian2008; Gomba & Homsy Reference Gomba and Homsy2010; Sui Reference Sui2014; Chaudhury & Chakraborty Reference Chaudhury and Chakraborty2015; Sui & Spelt Reference Sui and Spelt2015; Mac Intyre et al. Reference Mac Intyre, Gomba, Perazzo, Correa and Sellier2018; Dominguez Torres et al. Reference Dominguez Torres, Mac Intyre, Gomba, Perazzo, Correa, Lopez-Villa and Medina2020; Xu et al. Reference Xu, Chen, Liu, Sahu and Ding2021). The theoretical study of Ehrhard & Davis (Reference Ehrhard and Davis1991) employed the lubrication approximation of the Navier–Stokes equations by exploiting the disparity between the length scale of the droplet along the spreading direction and the drop height. They concluded that cooling the substrate along the length augments the spreading action. Furthermore, longitudinal heating opposes the capillary driving force responsible for isothermal spreading, emphasizing the importance of the heat transfer mechanism in regulating droplet motion. In a later work, Ehrhard (Reference Ehrhard1993) confirmed his numerical results and prescribed a power law for the spreading rate in both isothermal and non-isothermal conditions. Many years later, Chaudhury & Chakraborty (Reference Chaudhury and Chakraborty2015) discovered that thermocapillary-dominated spreading is fundamentally distinct from capillary-dominated spreading and established two scaling regimes for the transient spreading width of the droplet. In addition, the variation of surfactant concentration on the droplet surface was also reported to significantly alter the droplet migration and deformation characteristics (Leal Reference Leal2007; Poddar, Bandopadhyay & Chakraborty Reference Poddar, Bandopadhyay and Chakraborty2019a; Poddar et al. Reference Poddar, Mandal, Bandopadhyay and Chakraborty2019b)
Any attempt to resolve the contact line dynamics is challenged by the diverging stress condition near the contact line, known as the ‘contact line paradox’ (Snoeijer & Andreotti Reference Snoeijer and Andreotti2013). Several alternate routes had been suggested in the literature to address this moving contact line singularity within the continuum hydrodynamics, with proper consideration of the solid–fluid interaction (Snoeijer & Andreotti Reference Snoeijer and Andreotti2013). Among these, the models based on the Navier slip boundary condition at the solid–liquid interface relieve the singularity but come with the additional complexity of tracking the dynamic contact angle (Ehrhard & Davis Reference Ehrhard and Davis1991; Ren, Trinh & Weinan Reference Ren, Trinh and Weinan2015) in numerical calculations. A model based on molecular interaction potential, on the other hand, alleviates the contact line singularity by incorporating a realistically observed (Starov, Kalinin & Chen Reference Starov, Kalinin and Chen1994; Popescu et al. Reference Popescu, Oshanin, Dietrich and Cazabat2012) thin precursor film surrounding the drop. The latter approach not only eliminates explicit tracking of the dynamic contact angle but also allows for specifying an independent apparent contact angle. The precursor film approach was shown to be advantageous in a range of studies on wetting hydrodynamics under thermocapillary actuation (Gomba & Homsy Reference Gomba and Homsy2009, Reference Gomba and Homsy2010; Mac Intyre et al. Reference Mac Intyre, Gomba, Perazzo, Correa and Sellier2018). The numerical investigation of Gomba & Homsy (Reference Gomba and Homsy2010) used the same approach to classify the diversified spreading behaviours of a partially wetting drop in microgravity situations based on widely contrasting drop morphologies induced by a combination of thermal Marangoni stress and conjoining–disjoining pressure. Subsequently, Mac Intyre et al. (Reference Mac Intyre, Gomba, Perazzo, Correa and Sellier2018) analysed the effects of gravity and different forms of the molecular interaction potential. Their results provided a comprehensive understanding of the earlier experimental observations (Brzoska et al. Reference Brzoska, Brochard-Wyart and Rondelez1993) on droplet elongation and a saturation droplet migration speed under specific conditions.
Many industrial (Rafaï, Bonn & Boudaoud Reference Rafaï, Bonn and Boudaoud2004; Wang et al. Reference Wang, Lee, Peng and Lai2007a,Reference Wang, Zhang, Lee and Pengb) and biological (Hwang, Litt & Forsman Reference Hwang, Litt and Forsman1969; Walburn & Schneck Reference Walburn and Schneck1976; Yu et al. Reference Yu, Malcolm, Woolfson, Jones and Andrews2011) fluids exhibit complex mechanical responses when subjected to external forces, and a constant viscosity can no longer describe their rheology. A significant volume of experimental (Carré & Eustache Reference Carré and Eustache2000; Rafaïet al. Reference Rafaï, Bonn and Boudaoud2004; Wang et al. Reference Wang, Lee, Peng and Lai2007a,Reference Wang, Zhang, Lee and Pengb; Liang et al. Reference Liang, Wang, Lee, Peng and Su2009) as well as theoretical (Starov et al. Reference Starov, Kalinin and Chen1994; Flitton & King Reference Flitton and King2004; Miladinova, Lebon & Toshev Reference Miladinova, Lebon and Toshev2004) research was devoted to the spreading characteristics of isothermal droplets featuring time-independent and inelastic rheological responses. One essential rheological aspect of these fluids is the nonlinear variation of the shear stress with the strain rate (Poddar et al. Reference Poddar, Bandopadhyay and Chakraborty2019a,Reference Poddar, Mandal, Bandopadhyay and Chakrabortyc). Moreover, a branch of these fluids, known as the shear-thinning fluids, were found to possess an exclusive quality of suppressing the diverging stress condition at the contact line (Weidner & Schwartz Reference Weidner and Schwartz1994), a common problem associated with Newtonian fluids. Rafaïet al. (Reference Rafaï, Bonn and Boudaoud2004) were the first to present experimental results on the spreading dynamics of polymeric solutions that are shear thinning in nature. They related the power-law index of such fluids with the polymeric concentration and discovered that the shear-thinning fluids spread at a slower rate than Newtonian fluids. Subsequently, the work of Wang et al. (Reference Wang, Lee, Peng and Lai2007a) provided crucial insights into the spreading dynamics of both shear-thinning and shear-thickening fluids. Combining experiments and a thin film-based theory, they provided important information about film thickness and dynamic contact angle for both the capillary and gravity-dominated regimes. Later, Starov et al. (Reference Starov, Tyatyushkin, Velarde and Zhdanov2003) adopted a self-similar solution approach to establish theoretical spreading rules for both shear-thinning and shear-thickening fluids. In another work, Kheyfets & Kieweg (Reference Kheyfets and Kieweg2013) identified that shear-dependent rheology has a significant impact on the gravity-driven spreading of microbicide gel on the vaginal epithelium.
Despite considerable interest in the literature on the capillary-dominated spreading of non-Newtonian droplets, no attempt has been directed so far to investigate the interplay between the thermocapillary driving force and the non-Newtonian rheology of a sessile droplet. The intricate interfacial dynamics associated with the Marangoni flow are likely to interact with the shear-dependent flow conditions in a unique fashion, calling for further investigations. In the present work we address this issue through numerical simulations and seek to gain a fundamental understanding of thermocapillary spreading of partially wetting non-Newtonian droplets by employing an inelastic, time-independent power-law model for fluid rheology. Our hydrodynamic analysis, premised on a thin film model, unveils the dramatic influence of the non-Newtonian rheology on the spreading behaviour as well as the morphology of non-isothermal drops.
We have analysed the final shape and migration behaviour of the droplets for various practical combinations of the dimensionless thermocapillary strength, the power-law index and the equilibrium contact angle. Accordingly, the spreading behaviours have been categorized into three distinct regimes: Marangoni film, droplet and transition regimes. Even though a similar three regimes were reported previously for Newtonian droplets (Gomba & Homsy Reference Gomba and Homsy2010), the current work stands apart in bringing out the substantial changes in droplet spreading that are exclusive to non-Newtonian fluids. Variation in the power-law index 
$(n)$, quantifying the fluid's shear rate dependency, has been found to have a two-step influence on spreading. Firstly, it creates severe distinctive modulations in droplet form, height attained, time scales for transient spreading and many other spreading attributes. Secondly, the parameter 
$(n)$ serves as a switching mechanism between different spreading regimens. We envisage that the fundamental understanding of thermocapillary-driven non-Newtonian sessile drops obtained from the present study will pave the way for designing efficient microfluidic devices that can handle complex biofluids.
2. Mathematical formulation
The physical system considered here consists of a two-dimensional complex fluid droplet with density 
$\tilde {\rho }$ and viscosity 
$\tilde {\mu }$ that is deposited on a non-isothermal solid surface. We have defined a planar Cartesian coordinate system 
$(\tilde {x},\tilde {y})$ as shown in figure 1. Following continuum hydrodynamics, the governing equations for the flow field 
$(\tilde {p},\boldsymbol {\tilde {v}})$ are the continuity and Cauchy momentum equations
\begin{equation} \tilde{\boldsymbol{\nabla}} \boldsymbol{\cdot} \boldsymbol{\tilde{v}} = 0 \quad \text{and} \quad \tilde{\rho} \left(\frac{\partial \boldsymbol{\tilde{v}}}{\partial \tilde{t}}+\boldsymbol{\tilde{v}} \boldsymbol{\cdot} \tilde{\boldsymbol{\nabla}}\boldsymbol{\tilde{v}}\right) = \tilde{\boldsymbol{\nabla}} \boldsymbol{\cdot} \tilde{\boldsymbol{\mathsf{T}}}. \end{equation}
Here, 
$\tilde {\boldsymbol{\mathsf{T}}} = -p \mathbb {I} + \boldsymbol {\tilde {\tau }}$ denotes the stress tensor, and the deviatoric stress 
$\boldsymbol {\tilde {\tau }}$ is related to the strain rate tensor as 
$\boldsymbol {\tilde {\tau }}=2 \tilde {\mu } \tilde {\boldsymbol{\mathsf{D}}}$. The shear rate dependence of the fluid apparent viscosity has been captured by the power-law constitutive relation, given by (Deen Reference Deen1998; Garg et al. Reference Garg, Kamat, Anthony, Thete and Basaran2017)
\begin{equation} \tilde{\mu}\left(\tilde{\dot{\gamma}}\right)= \mu_0 \left| 2 \tilde{\mathcal{M}} \tilde{\dot{\gamma}}\right|^{n-1}, \end{equation}
where 
$1/\tilde {\mathcal {M}}$ denotes the characteristic deformation rate and 
$\mu _0$ is the constant viscosity of the Newtonian base fluid, which is recovered from the above equation for 
$n=1$. The exponent 
$n$ is known as the power-law index, where 
$n >1$, 
$n =1$ and 
$n<1$ stand for the shear-thickening (or dilatant), Newtonian and shear-thinning (or pseudoplastic) fluids, respectively. In addition, 
$\tilde {\dot {\gamma }}$ denotes the second invariant of the strain rate tensor 
$\tilde {\boldsymbol{\mathsf{D}}}$, i.e.
\begin{equation} \tilde{\dot{\gamma}} = \left[\tfrac{1}{2} (\tilde{\boldsymbol{\mathsf{D}}} \colon \tilde{\boldsymbol{\mathsf{D}}})\right]^{1/2}, \end{equation}where
\begin{equation} \tilde{\boldsymbol{\mathsf{D}}}=\tfrac{1}{2} \left[(\tilde{\boldsymbol{\nabla}} \tilde{\boldsymbol{v}})+(\tilde{\boldsymbol{\nabla}} \tilde{\boldsymbol{v}})^\mathrm{T}\right]. \end{equation}
It is worth noting that the term 
$\tilde {\mathcal {M}}$ appearing in (2.2) is different from the consistency index that is commonly used in the power-law equation (Bird, Armstrong & Hassager Reference Bird, Armstrong and Hassager1987). Unlike the conventional form, the present equation (Garg et al. Reference Garg, Kamat, Anthony, Thete and Basaran2017) eliminates the dependency of the reference scales of time, velocity and pressure on 
$n$ (Gorla Reference Gorla2001), thereby facilitating direct comparison of spreading characteristics of non-Newtonian and Newtonian droplets.
Figure 1. Schematics of a two-dimensional non-Newtonian sessile droplet with height 
$\tilde {h}(\tilde {x},\tilde {t})$ and base width 
$\tilde {w}(\tilde {t})$ on a flat substrate having a linear temperature distribution 
$\tilde {T}_w(\tilde {x})$. The temperature gradient at the drop–air interface creates a surface tension gradient, triggering Marangoni flow, and subsequently, the drop moves along the substrate. A constant film of equilibrium thickness 
$\tilde {h}_f$ surrounds the droplet. The drop has partially wetting characteristics with a contact angle 
$\theta$ at the solid–liquid interface.
The applicability of the power-law constitutive relation in two limiting flow conditions has often been questioned. Firstly, it predicts infinite fluid viscosity for shear-thinning fluids 
$( n<1 )$ at locations in the flow domain where the shear rate (and, consequently, the shear stress) becomes zero (Myers Reference Myers2005). Although criticized for its unphysical viscosity predictions at low shear rates (Acrivos, Shah & Petersen Reference Acrivos, Shah and Petersen1960; Myers Reference Myers2005), the same model was used for capturing the essential physics of the stress-dependent rheology of thin films under a broad spectrum of practical circumstances, such as the spreading of a drop over a solid substrate (Starov et al. Reference Starov, Tyatyushkin, Velarde and Zhdanov2003), a droplet flowing down an inclined plane (Perazzo & Gratton Reference Perazzo and Gratton2003; Miladinova et al. Reference Miladinova, Lebon and Toshev2004; Hu & Kieweg Reference Hu and Kieweg2012; Ruyer-Quil, Chakraborty & Dandapat Reference Ruyer-Quil, Chakraborty and Dandapat2012; Noble & Vila Reference Noble and Vila2013), blade coating (Ross, Wilson & Duffy Reference Ross, Wilson and Duffy1999), rupture of thin films (Garg et al. Reference Garg, Kamat, Anthony, Thete and Basaran2017), vaginal delivery of microbicides (Kheyfets & Kieweg Reference Kheyfets and Kieweg2013) and so on. Along with this, in the present work, the Marangoni stress condition at the free surface alleviates the anomalous infinite viscosity condition due to the zero shear boundary conditions employed in a host of earlier studies (Acrivos et al. Reference Acrivos, Shah and Petersen1960; Carré & Eustache Reference Carré and Eustache2000; Flitton & King Reference Flitton and King2004; Noble & Vila Reference Noble and Vila2013). Secondly, according to the ‘contact line paradox’ for Newtonian fluids, a diverging shear rate is predicted at the moving contact line (Snoeijer & Andreotti Reference Snoeijer and Andreotti2013). Contrastingly, the viscosity mathematically goes to zero for shear-thinning fluids, thus avoiding the thermodynamically inconsistent diverging dissipation condition. However, laboratory experiments suggest that the viscosity asymptotically attains the value of the solvent viscosity instead of becoming zero. Rafaïet al. (Reference Rafaï, Bonn and Boudaoud2004) showed that the length scale for this transition is 
${\sim }100\,\mathrm {nm}$, which is of the same order of magnitude as the zone of intermolecular forces. Now, the molecular film thickness 
$\tilde {h}_*$ chosen in the present study (details in § 3) remains above the viscosity transition length scale, rendering the power-law model both physically and mathematically consistent in the present scenario.
The liquid–air interface at 
$\tilde {y} = \tilde {h}(\tilde {x},\tilde {t})$ is subject to the normal and tangential stress balance conditions
$$\begin{gather} \boldsymbol{\hat{n}}\boldsymbol{\cdot} \boldsymbol{\tilde{\boldsymbol{\mathsf{T}}}} \boldsymbol{\cdot} \boldsymbol{\hat{n}} = \tilde{\kappa} \tilde{\gamma}, \end{gather}$$
$$\begin{gather}\boldsymbol{\hat{n}}\boldsymbol{\cdot} \tilde{\boldsymbol{\mathsf{T}}} \boldsymbol{\cdot} \boldsymbol{\hat{t}} = \boldsymbol{\hat{t}} \boldsymbol{\cdot} \tilde{\boldsymbol{\nabla}} \tilde{\gamma}, \end{gather}$$
where 
$\boldsymbol {\hat {n}}$, 
$\boldsymbol {\hat {t}}$ are the outward unit normal and tangential vectors at the interface, respectively; and 
$\tilde {\kappa } = -\tilde {\boldsymbol {\nabla }} \boldsymbol {\cdot } \boldsymbol {\hat {n}}$ denotes the mean curvature of the interface.
At the interface between the drop and the solid surface 
$(\tilde {y}=0)$ the no-slip and the no-penetration conditions apply, i.e.
\begin{equation} \tilde{u}(\tilde{x},0,\tilde{t})=0 \quad \text{and} \quad \tilde{v}(\tilde{x},0,\tilde{t})=0 , \end{equation}
where 
$\tilde {u}, \tilde {v}$ are the fluid velocity components parallel and normal to the solid surface. Now, the kinematic boundary condition at the free surface is given by
\begin{equation} \frac{{\rm D}\tilde{f}}{{\rm D}\tilde{t}}=\frac{\partial \tilde{f}}{\partial \tilde{t}} + \boldsymbol{\tilde{v}} \boldsymbol{\cdot} \tilde{\boldsymbol{\nabla}} f = 0, \end{equation}
where 
$\tilde {f}(\tilde {x},\tilde {y},\tilde {t})=\tilde {y}-\tilde {h}(\tilde {x},\tilde {t})$.
The molecular interaction between the solid and liquid has been modelled with a disjoining–conjoining pressure term (Derjaguin & Kusakov Reference Derjaguin and Kusakov1936; Mitlin & Petviashvili Reference Mitlin and Petviashvili1994) 
$\tilde {\varPi }(\tilde {h})$ of the form
\begin{equation} \tilde{\varPi}=\tilde{\mathcal{K}} \left[\left(\frac{\tilde{h}_*}{\tilde{h}}\right)^{3}- \left(\frac{\tilde{h}_*}{\tilde{h}}\right)^{2}\right], \end{equation}
where 
$\tilde {h}_*$ is the energetically favoured thickness of the molecular film (Gomba & Homsy Reference Gomba and Homsy2009). Thus, the total pressure within the drop can be written as
\begin{equation} \tilde{p}={-}\tilde{\gamma}\frac{\partial^2 \tilde{h}}{\partial \tilde{x}^2}- \tilde{\varPi}, \end{equation}
where the first term represents the Laplace pressure as a consequence of the interface curvature, and the second term stands for the disjoining pressure 
$\tilde {\varPi }(\tilde {h})$ due to van-der-Waals interaction at the solid–liquid interface. This model alleviates the singularity at the contact line of the drop by incorporating an experimentally observed thin 
$({\sim }10\unicode{x2013}100\,\text {nm})$ precursor film ahead of the drop (Starov et al. Reference Starov, Kalinin and Chen1994; Hoang & Kavehpour Reference Hoang and Kavehpour2011; Popescu et al. Reference Popescu, Oshanin, Dietrich and Cazabat2012). Also, for the partially wetting systems, the equilibrium contact angle is related to the intermolecular forces through the term 
$\tilde {\mathcal {K}}$ as 
$\tilde {\mathcal {K}} = 2 \widetilde {\gamma _0}(1-\cos (\theta _e)) / \tilde {h}_*$.
It is worth mentioning that mathematical calculations predict a diverging shear rate at the contact line, i.e. 
$\tilde {\dot {\gamma }} \to \infty$. Thus, the power-law model for the shear-thinning fluids 
$( n<1 )$ gives a zero viscosity (see (2.2)) in the contact line zone so that the viscous dissipation remains finite (Weidner & Schwartz Reference Weidner and Schwartz1994; Rafaïet al. Reference Rafaï, Bonn and Boudaoud2004), keeping the model thermodynamically consistent. Interestingly, this suppression of contact line singularity occurs without using a disjoining pressure or slip model, and the contact line movement remains possible. However, accommodation of the partially wetting conditions 
$(\theta _e > 0)$ necessitates considering the conjoining–disjoining pressure 
$(\tilde {\varPi })$. Moreover, the present model facilitates developing a general analysis framework for shear-thinning, shear-thickening and Newtonian fluids.
The surface tension at the liquid–air interface has been considered to vary with the temperature as (Leal Reference Leal2007)
\begin{equation} \tilde{\gamma} = \tilde{\gamma}_0-\tilde{\sigma} (\tilde{T}-\tilde{T}_0), \end{equation}
where 
$\tilde {\gamma }_0$ represents the surface tension at 
$\tilde {T}=\tilde {T}_0$ and the positive constant 
$\sigma$ denotes the temperature gradient of surface tension. This thermally induced gradient in surface tension causes Marangoni stress at the interface, given as
\begin{equation} \hat{\tau}=\frac{\partial \tilde{\gamma}}{\partial \tilde{x}} = \frac{\partial \tilde{\gamma}}{\partial \tilde{T}} \frac{\partial \tilde{T}}{\partial \tilde{x}}. \end{equation}In order to capture the spreading characteristics of the drop under the realm of long-wave/lubrication approximation (Leal Reference Leal2007; Eddi, Winkels & Snoeijer Reference Eddi, Winkels and Snoeijer2013; Chaudhury & Chakraborty Reference Chaudhury and Chakraborty2015), the following conditions in terms of the lubrication parameters are to be satisfied (Deen Reference Deen1998), i.e.
\begin{equation} (\tilde{h}_c/\tilde{w})^2 \ll 1 \quad \text{and}\quad Re (\tilde{h}_c/\tilde{w}) \ll 1, \end{equation}
where 
$\tilde {h}_c$ is the maximum droplet height, 
$\tilde {w}$ is the equilibrium foot width of the droplet and 
$Re$ is the Reynolds number, defined as 
$Re=\tilde {\rho } \tilde {U} \tilde {h}_c/\tilde {\mu }_0$. Here, 
$\tilde {U}$ represents the characteristic spreading speed of the drop, 
$\tilde {\mu }_0$ is the viscosity of the Newtonian base fluid and 
$\tilde {\rho }$ is the density of the fluid. Typical values of different dimensional quantities can be estimated from the reported experiments of Pratap et al. (Reference Pratap, Moumen and Subramanian2008): 
$\tilde {\mu }_0 = 10^{-3}$ Pa s, 
$\tilde {\rho } = 1000\,{\rm kg}\,{\rm m}^{-3}$, 
$\tilde {h}_c \le 0.3\,\text {mm}$, 
$1.2\,\text {mm}\le \tilde {w} \le 3.2\,\text {mm}$ and 
$200 \le \tilde {U} \le 800\,\mathrm {\mu }{\rm m}\,{\rm s}^{-1}$. We have considered 
$\tilde {w}=2 \tilde {R}$, where 
$\tilde {R}$ is the droplet footprint radius obtained from literature. Consequently, the lubrication parameters are found within the acceptable limits, i.e. 
$(\tilde {h}_c/\tilde {w})^2 \lesssim 0.063$ and 
$Re \tilde {h}_c/\tilde {w} \lesssim 0.06$, thereby justifying the lubrication approximation in the present work.
Under the lubrication approximation, the inertia term in the Cauchy momentum equation (
$\tilde {\boldsymbol {v}}\boldsymbol {\cdot } \boldsymbol {\nabla } \tilde {\boldsymbol {v}}$ in (2.1)) becomes negligible. Such a flow regime is dominated by the viscous and capillary effects. This approximation can be exploited to derive a single time-dependent differential equation for droplet height 
$\tilde {h}(\tilde {x},\tilde {t})$. Ideally, the lubrication theory was developed for vanishing contact angles at the drop–solid interface. Thus, its application in partial wetting conditions, especially when the contact angle is 
$40^{\circ }$ or higher, is questionable since the ratio 
$\tilde {h}/\tilde {w}$ may not be necessarily small. However, results of lubrication theory were compared against the full Navier–Stokes simulations for a multitude of flow problems involving contact angles as high as 
$40^{\circ }$, and, only small deviations were reported (Goodwin & Homsy Reference Goodwin and Homsy1991; Mitlin & Petviashvili Reference Mitlin and Petviashvili1994; Schwartz & Eley Reference Schwartz and Eley1998; Perazzo & Gratton Reference Perazzo and Gratton2004; Diez & Kondic Reference Diez and Kondic2007; Chaudhury & Chakraborty Reference Chaudhury and Chakraborty2015). The analysis of Perazzo & Gratton (Reference Perazzo and Gratton2004) revealed that, for high contact angles, the lubrication theory is not so accurate in predicting the detailed velocity field. However, the global flow properties like the cross-section of the thin film, volumetric flow rate 
$Q = \int _{0}^{h} u\,{{\rm d}y}$ (2.18) and average velocity are in close agreement with the full Navier–Stokes solution. They attributed this dramatic agreement to the redistribution of momentum without alterations in its total value. As a result, the inaccuracies of velocity distribution cancel out, leading to almost accurate average values.
We consider small drop sizes for which the surface tension dominates the gravity-driven deformation and spreading of the drop (Leal Reference Leal2007). The gravity effect is neglected when the Bond number 
$(Bo)$, defined as the ratio between the horizontal gradient of hydrostatic pressure and the capillary force, i.e. 
$Bo = \tilde {\rho } g \tilde {h}_c \tilde {R}/\tilde {\gamma }_0$, remains much less than unity. Pratap et al. (Reference Pratap, Moumen and Subramanian2008) performed experiments for a wide range of surface tension gradients and drop sizes 
$0.6\,{\rm mm}\le \tilde {R}\le 1.6$ mm, 
$\tilde {h}_c \le 0.3$ mm and found 
$Bo$ in the range of 
$10^{-3}\unicode{x2013}10^{-2}$. Similarly, Schwartz & Eley (Reference Schwartz and Eley1998) performed numerical simulations with and without the inclusion of gravity effects and found that, for a footprint radius 
$\tilde {R} \le 4$ mm, the gravity effects are unimportant for spreading on a horizontal substrate. On the other hand, the assumption of negligible gravity effect has practical relevance in microgravity conditions that are observed in many of the experiments conducted in orbiting space aircrafts (Subramanian & Balasubramaniam Reference Subramanian and Balasubramaniam2001).
Under the lubrication approximation, it is reasonable to consider that the thermal Péclet number is very small and conduction remains the dominant heat transfer mode within the drop (Oron & Rosenau Reference Oron and Rosenau1994; Gomba & Homsy Reference Gomba and Homsy2010). Thus, the thermal field in the drop fluid is governed by the energy equation
\begin{equation} \frac{\partial^2 \tilde{T}}{\partial \tilde{y}^2}=0, \end{equation}which is complemented by the boundary conditions
$$\begin{gather} \text{at }\tilde{y}=0:\quad \tilde{T}=\tilde{T}_w, \end{gather}$$
$$\begin{gather}\text{at }\tilde{y}=\tilde{h}: \quad -\tilde{k}_l \frac{\partial \tilde{T}}{\partial \tilde{y}} = \mathcal{H}_\infty(\tilde{T}-\tilde{T}_\infty). \end{gather}$$
Here, 
$\tilde {k}_l$ is thermal conductivity of the liquid, 
$\tilde {T}_\infty$ is the ambient temperature and 
$\mathcal {H}_\infty$ is the convective heat transfer at the interface. Thus, the temperature field within the drop is obtained as
\begin{equation} \tilde{T}= \tilde{T}_w-\frac{( \tilde{T}_w- \tilde{T}_\infty)\mathcal{H}_\infty \tilde{y} / \tilde{k}_l}{1+\mathcal{H}_\infty \tilde{h} / \tilde{k}_l}. \end{equation}
Again using the lubrication approximation, it has been found that 
$\mathcal {H}_\infty \tilde {y} /\tilde {k}_l,\mathcal {H}_\infty \tilde {h} /\tilde {k}_l \ll 1$ (Ehrhard & Davis Reference Ehrhard and Davis1991; Chaudhury & Chakraborty Reference Chaudhury and Chakraborty2015), leading to 
$\tilde {T} \approx \tilde {T}_w(\tilde {x}).$ Motivated by the high thermal conductivity of the widely used substrate materials, e.g. silicon (Brzoska et al. Reference Brzoska, Brochard-Wyart and Rondelez1993), we assume a linearly varying surface temperature distribution, such that
\begin{equation} \frac{\partial \tilde{T}_w}{\partial \tilde{x}}= \varGamma. \end{equation}
Substituting the above expression for temperature gradient (2.16) and using (2.10) in (2.11), we get the expression for Marangoni stress as 
$\hat {\tau }= -\sigma \varGamma$, which is 
$> 0$. Therefore, the tangential stress boundary condition (2.5b) takes the form
\begin{equation} \tilde{\tau}_{\tilde{x}\tilde{y}}|_{\tilde{h}} = \hat{\tau}. \end{equation} Next, we non-dimensionalize the variables by choosing the following characteristic scales: 
$x_c=y_c=a$, 
$t_c= \mu x_c^4/\gamma _0 h_c^3= \mu _c a/\gamma _0$, 
$p_c=\gamma _0/a$, 
$\mu _c=\mu _0$, 
$u_c= x_c/t_c=\gamma _0/\mu _0$ and 
$\mathcal {M}_c=y_c/u_c$. Here, 
$a = \sqrt {\gamma _0/\rho g}$ is the capillary length scale. Note that the present characteristic time scale is the capillary time scale and it differs from its Newtonian counterpart defined in many earlier studies (Gomba & Homsy Reference Gomba and Homsy2010; Chaudhury & Chakraborty Reference Chaudhury and Chakraborty2015) only by a factor of three. In what follows, we drop the ‘
$\tilde {}$’ from different quantities and work only with their dimensionless forms.
Using the continuity (2.1a) and the no-penetration condition (2.6b) in (2.7), we can derive the following form of the kinematic boundary condition (Leal Reference Leal2007):
\begin{equation} \frac{\partial h}{\partial t} + \frac{\partial Q}{\partial x} = 0. \end{equation}
Here the flow rate is defined as 
$Q= \int _{0}^{h} u\,{{\rm d}y}$. Upon using the constitutive relation introduced in (2.2) and using the no-slip condition (2.6a,b), the 
$x$ component of the flow velocity 
$(u)$ is obtained as (details in § A)
\begin{equation} u = \frac{\mathcal{M}^{1/n-1}}{p_x} \frac{n}{n+1} \left[\left|p_x(y-h)+\beta\right|^{{1}/{n}+1} - \left|\beta-hp_x\right|^{{1}/{n}+1} \right], \end{equation}
where 
$p_x$ denotes the axial pressure gradient 
$\partial p/\partial x$. Also, we have introduced a dimensionless thermocapillary number 
$\beta$, defined as 
$\beta = \hat {\tau }/p_c =\hat {\tau } a/{\gamma _0}$. It was reported that for the thin film dynamics of Newtonian drops under thermocapillary action, the streamwise flow velocity is a linear superposition of two driving forces, namely a Poiseuille flow induced by the pressure gradient 
$p_x$ and a shear flow originating from the Marangoni stress 
$\hat {\tau }$ at the free surface (Brzoska et al. Reference Brzoska, Brochard-Wyart and Rondelez1993). In strong contrast, the velocity profile in (2.19) suggests a nonlinear coupling of these driving mechanisms, the resultant of which is dictated both by the power-law index 
$(n)$ and the thermocapillary number 
$(\beta )$.
Consequently, the flow rate is given by
\begin{align} Q &= \mathcal{M}^{{1}/{n}-1} \left[\frac{n^2}{(n+1)(2n+1)}\frac{1}{p_x^2} \left\{\beta^{{1}/{n}+2} - \left|\beta-h p_x\right|^{{1}/{n}+1} \left(\beta- h p_x\right)\right\} \right. \nonumber\\ &\quad \left.- \frac{n}{(n+1)}\frac{h}{p_x} \left|\beta-h p_x\right|^{{1}/{n}+1} \right]. \end{align}
Using the above form of 
$Q$ in (2.18) yields the following highly nonlinear partial differential equation:
\begin{align} & \frac{\partial h}{\partial t} + \mathcal{M}^{{1}/{n}-1} \frac{\partial }{\partial x} \left[\frac{n^2}{(n+1)(2n+1)}\frac{1}{p_x^2} \left\{\beta^{{1}/{n}+2} - \left|\beta-h p_x\right|^{{1}/{n}+1} \left(\beta- h p_x\right) \right\} \right.\nonumber\\ &\quad \left.- \frac{n}{(n+1)}\frac{h}{p_x} \left|\beta-h p_x\right|^{{1}/{n}+1} \right] = 0. \end{align}
Equation (2.21) serves as the governing equation for the evolution of the free surface height 
$h(x,t)$. Substituting 
$n = 1$ in the above equation recovers the form of the same equation derived earlier for Newtonian fluids (Gomba & Homsy Reference Gomba and Homsy2010). It is noteworthy that the modulus sign in different terms of (2.21) has to be retained throughout the calculations (Ross et al. Reference Ross, Wilson and Duffy1999) to capture both the physical scenarios: 
$\beta >hp_x$ and 
$\beta < h p_x$. Moreover, the nonlinear contribution of the term 
$(\beta -hp_x)$ in the equation embodies a stark contrast in the coupling behaviour of the Marangoni and viscous stresses against the backdrop of the linear superposition of different effects in the Newtonian counterpart of the same equation (Gomba & Homsy Reference Gomba and Homsy2010; Mac Intyre et al. Reference Mac Intyre, Gomba, Perazzo, Correa and Sellier2018).
Here, the dimensionless pressure profile within the thin film is
\begin{equation} p={-}\frac{\partial^2 h}{\partial x^2}- \varPi(h), \end{equation}
where the dimensionless disjoining pressure 
$\varPi (h)$ is given as
\begin{equation} \varPi=\mathcal{K} \left[ \left(\frac{h_*}{h}\right)^3 - \left(\frac{h_*}{h}\right)^2\right], \end{equation}
with 
$\mathcal {K}= ({2(1-\cos (\theta _e))}/{\tilde {h}_*}) ({x_c^2}/{h_c^2})$.
It is to be noted that the shear-dependent rheology influences both the capillary and thermocapillary effects. Although (2.10) suggests that the surface tension is not affected by the power-law viscosity through the parameter 
$n$, the capillary pressure term (2.9) inherently contains the effect of 
$n$ through the droplet curvature 
$\partial ^2 h/\partial x^2$. Along similar lines, the thermocapillary strength, represented by the parameter 
$\beta$, is independent of 
$n$. Nevertheless, the actual influence of 
$\beta$ on the fluid flow and drop profile is not independent from 
$n$, as described in (2.19) and (2.21), respectively.
3. Solution methodology
The governing equation (2.21) coupled with (2.22) has been solved numerically using the finite element package COMSOL Multiphysics
$\circledR$ v.6.0. Periodic boundary conditions have been applied for the variables 
$h$ and 
$p$ at the two ends of the computational domain (Schwartz et al. Reference Schwartz, Roy, Eley and Petrash2001; Mac Intyre et al. Reference Mac Intyre, Gomba, Perazzo, Correa and Sellier2018). The length of the computational domain has been chosen between 
$200$ and 
$400$, as required to capture the essential physics of droplet spreading and deformation. The PARDISO (parallel direct sparse solver for clusters) solver has been deployed for solving the system of linear algebraic equations arising from the discretization. The pressure gradient term 
$p_x$ in the denominator of the second term of (2.21) poses an additional numerical hurdle for the solver. This issue has been tackled by adding a very small number 
$eps \approx 10^{-15}$ to 
$p_x$ whenever it reaches an exact zero. The converged simulation results for the Newtonian fluids 
$(n = 1)$ have been compared with the finite-difference results of Gomba & Homsy (Reference Gomba and Homsy2010) and finite element results of Mac Intyre et al. (Reference Mac Intyre, Gomba, Perazzo, Correa and Sellier2018) in figure S-1(a–c) of the supplementary material available at https://doi.org/10.1017/jfm.2022.900. An excellent agreement between the present numerical simulations and the known results can be found in these figures for three distinct values of the equilibrium contact angle (
$\theta _e = 5^{\circ }, 10^{\circ }$ and 
$30^{\circ }$) and different times 
$(t)$. It has been observed that an increased shear-thinning 
$(n<1)$ behaviour takes lesser computational time than the Newtonian case, and the converse happens for an intensifying shear-thickening 
$(n>1)$ nature of the fluid.
The choice of a realistic value of the molecular film thickness 
$\tilde {h}_*\sim 10\,\mathrm {nm}$ demands a computational mesh size in the same range, i.e. 
${\rm \Delta} x \approx 10^{-6}\unicode{x2013}10^{-5}$ to ensure sufficient spatial resolution in the moving contact line region (Diez & Kondic Reference Diez and Kondic2001; Gaskell et al. Reference Gaskell, Jimack, Sellier and Thompson2004; Gomba & Homsy Reference Gomba and Homsy2010); however, this increases the computational cost to an unacceptable level. To optimize the computational performance, Gomba and co-workers (Gomba & Homsy Reference Gomba and Homsy2010; Mac Intyre et al. Reference Mac Intyre, Gomba, Perazzo, Correa and Sellier2018) showed that a choice of 
$h_*\approx 0.01$ provides optimized computational performance without affecting the overall flow behaviour in the thermocapillary spreading of Newtonian drops. Accordingly, we choose 
$h_*={\rm \Delta} x=0.01$ throughout the simulations.
The breakup of a thin sessile drop is characterized by its rupture into multiple drops connected by ultra-thin films (Gomba & Homsy Reference Gomba and Homsy2009; Mac Intyre et al. Reference Mac Intyre, Gomba, Perazzo, Correa and Sellier2018). Both surface tension and intermolecular forces at the interface play crucial roles in dictating the time and intensity of this rupture. The presently reported numerical simulations account for these physical aspects by using a widely adopted form (Derjaguin & Kusakov Reference Derjaguin and Kusakov1936; Mitlin & Petviashvili Reference Mitlin and Petviashvili1994; Gomba & Homsy Reference Gomba and Homsy2009) of molecular interaction potential (2.8) that models both attractive (destabilizing) and repulsive (stabilizing) forces at the interface. This consideration gives rise to an equilibrium film thickness 
$\tilde {h}_*$, enabling the numerical framework to provide a solution of the lubrication equation (2.21) for all time and, thus, model the formation of multiple connected drops.
The initial condition has been chosen as the steady-state profile of the drop shape without temperature gradient 
$(\beta =0)$ as derived by Gomba & Homsy (Reference Gomba and Homsy2009). This choice of initial condition is motivated by the physically consistent drop shape parameters, such as the initial apparent contact angle 
$(\theta _e)$ between the liquid and the solid surface, the maximum height of the drop surface 
$(h_{max})$ and the cross-sectional area 
$(\mathcal {A})$ of the drop, which can be related to the absorbed film thickness around the drop 
$(h_f)$ in its equilibrium condition. Intending to follow a consistent comparison basis for the results with different power-law indices 
$(n)$, we have chosen identical initial conditions for drops with different 
$n$.
4. Results and discussions
In this section we demonstrate the results obtained from the numerical simulations and focus on the intricate interplay between the thermocapillary actuation mechanism and the shear-dependent fluid rheology. The values of the different parameters used to perform these simulations are chosen in coherence with the relevant experimental literature. Experiments related to thermocapillary actuation of non-Newtonian droplets reported the surface tension at the reference temperature 
$(\tilde {\gamma }_0)$ in the range of 
$2 \times 10^{-2}\unicode{x2013}5 \times 10^{-2}\,{\rm N}\,{\rm m}^{-1}$ (Wang et al. Reference Wang, Lee, Peng and Lai2007a,Reference Wang, Zhang, Lee and Pengb). Similarly, the measurements of Chen et al. (Reference Chen, Troian, Darhuber and Wagner2005), Brzoska et al. (Reference Brzoska, Brochard-Wyart and Rondelez1993) and Pratap et al. (Reference Pratap, Moumen and Subramanian2008) suggest that the temperature gradient of surface tension is in the range 
$\tilde {\sigma } \approx 10^{-5}\unicode{x2013}10^{-4}\,{\rm N}\,{\rm m}^{-1}\,^{\circ }{\rm C}$, and the longitudinal temperature gradient is in the range of 
${\partial \tilde {T}}/{\partial \tilde {x}} \approx 10^{2}\unicode{x2013}10^{3}\,^{\circ }{\rm K}\,{\rm m}^{-1}$. We choose the parameter 
$\tilde {\mathcal {M}} = \tilde {\mu }_0$ with regard to recovering the Newton's law of viscosity from the constitutive relation in (2.2) in the limit of 
$n = 1$. The capillary length scale is calculated to be in the range 
$a \approx 10^{-4}\unicode{x2013}10^{-3}$ m. Based on different practical ranges of the dimensional quantities, we have chosen the different dimensionless parameters in the following ranges: 
$\beta \approx 0.001 \ \mathrm {to} \ 0.03$ and 
$\mathcal {M} \approx 14$. The cross-sectional area of the two-dimensional droplets has been approximated from the experimentally measured radius of the three-dimensional droplets to obtain 
$\mathcal {A} \approx 10$ (Brzoska et al. Reference Brzoska, Brochard-Wyart and Rondelez1993; Gomba & Homsy Reference Gomba and Homsy2010). The values of the power-law index 
$n$ are considered in the range 0.4–1.6, motivated by the rheological responses of various natural as well as synthetic polymeric solutions and particulate suspensions (Rafaïet al. Reference Rafaï, Bonn and Boudaoud2004; Wang et al. Reference Wang, Lee, Peng and Lai2007a,Reference Wang, Zhang, Lee and Pengb; Kheyfets & Kieweg Reference Kheyfets and Kieweg2013).
4.1. Summary of non-Newtonian spreading regimes
In figure 2 we present the different spreading regimes observed due to variations in the power-law index 
$n$, equilibrium contact angle 
$\theta _e$ and the thermocapillary number 
$\beta$. Figure 2 conveys the key message of this work. In the subsequent sections we describe the different spreading regimes and provide physical justifications for the central findings. The Newtonian spreading regimes reported earlier (Gomba & Homsy Reference Gomba and Homsy2010) have been highlighted in the above regime maps. Classification of the dynamical phase regimes has been achieved after reflecting on the results of the numerical simulation for a wide range of parameter sets 
$(\beta, n, \theta _e)$. Similar to Gomba & Homsy (Reference Gomba and Homsy2009), three regimes of droplet spreading were identified, viz., the Marangoni film regime, transition regime and the droplet regime, the distinguishing features of which are mentioned below. For numerical demarcation, if the speeds of leading and trailing edges of the drop have a relative difference of less than 1 % throughout the period of droplet spreading, then such a drop is categorized under the droplet regime. On the contrary, when the relative difference in the speeds is higher than the said tolerance, and the droplet gradually splits into multiple parts, it has been identified as the transition regime. Usually, the droplet develops a linear profile at the beginning, and after further spreading, the film ruptures and splits into multiple droplets from its rear. The complex transient behaviour and film rupture pose a difficulty in characterizing the droplet spreading in the transition regime. Finally, the droplets under the Marangoni film regime are those having a high speed difference at the two edges, but do not split into multiple droplets, and have the characteristic linear profile. In this case, the leading front of the droplet moves at a faster (identified by the said tolerance) velocity than the rear; hence, the width of the droplet increases with time. The trailing edge has a characteristic capillary ridge and exhibits a linear profile, with its slope decreasing over time.
Figure 2. Regime maps depicting the various modes of droplet spreading in the 
$n-\theta _e$ plane (
$\theta _e$ in degrees) for thermocapillary numbers (a) 
$\beta = 0.001$, (b) 
$\beta = 0.007$ and (c) 
$\beta = 0.03$. The points inside the black dashed rectangles denote the Newtonian (
$n = 1$) results reported by Gomba & Homsy (Reference Gomba and Homsy2010). Due to difference in the non-dimensionalization schemes, the parameter 
$\beta$ is related to the parameter 
$B$ defined in Gomba & Homsy (Reference Gomba and Homsy2010) by the relation 
$\beta = 2B/3$.
Furthermore, we have checked the regimes reported in figure 2 with 
$h_* = 0.005$ and found no change in the regime boundaries. Lowering the 
$h_*$ values has an effect only in the transition regime (defined in § 4.4), where the number of split drops and their shapes vary and have no effect on the spreading behaviour of droplets in other regimes, thereby justifying the choice of 
$h_*$ in the current numerical simulations much above the experimentally observed molecular film thickness.
4.2. Non-Newtonian effects on Marangoni films
Figure 3(a) illustrates the axial variation of the free surface profile at time 
$t=3\times 10^4$ and for different values of the power-law index 
$n$. It is observed that a shear-thickening 
$(n= 1.4)$ drop spreads much faster than a Newtonian drop, while a shear-thinning droplet 
$(n = 0.6)$ moves at a slower pace. The observed nature of relative speeds for different non-Newtonian drops is against the intuitive prediction that the speed of the contact line should reduce due to a corresponding decrease in viscous resistance offered by shear-thinning fluids. Nevertheless, the experimental observations of Rafaï et al. (Reference Rafaï, Bonn and Boudaoud2004) under isothermal conditions support our findings. They pointed out that the flow near the contact line region gets accelerated due to the amplified shear-thinning effect. At the same time, this acceleration mechanism is counteracted by the escalating shear rates in this zone, causing faster displacement of the corner points in the contact line region. Thus, the contact angle is lowered, and subsequent suppression of the driving mechanism for spreading results. As a cumulative effect, the contact line moves at a slower pace.
Figure 3. (a) Droplet profiles at time instant 
$t=3\times 10^4$ for different power-law indices 
$n$. (b) Scaled droplet height (
$h \beta ^{1/n}/\mathcal {M}^{1-1/n}$) versus axial distance relative to the rear contact point 
$(x-x_0)$ at large times. Other parameters are 
$\mathcal {A}=10,\beta =0.007, \theta _e=5^{\circ }$. The dotted lines indicate the asymptotic linear profile in each case. (c) Variation of the normalized slope of the droplet shape in the asymptotic regime on the 
$(n, \beta )$ plane.
A distinguishing feature of drop shape found in the above figure is a linear bulk region, bearing similarity with the formation of a capillary ridge in Marangoni films (Kalliadasis, Kiyashko & Demekhin Reference Kalliadasis, Kiyashko and Demekhin2003; Sauleda et al. Reference Sauleda, Hsieh, Xu, Tilton and Garoff2022). The droplet's shape near the advancing front becomes steeper with the augmented shear-thickening nature of the fluid. The shear-thinning drop (
$n=0.6$) has been found to show a significant deviation from the original curvature, while its Newtonian (
$n=1$) and shear-thickening (
$n=1.4$) counterparts have already developed a linear shape in the bulk region. The physical mechanism behind this behaviour is described below.
4.2.1. Asymptotic solution
To gain physical insight into the spreading kinematics, we appeal to (2.20). The relative importance of different terms can be estimated by defining a variable 
$\xi = {hp'}/{\beta }$, which gives the following alternative form of the flow rate:
\begin{align} Q &= \mathcal{M}^{{1}/{n}-1}\beta^{1/n}h^2\left[\frac{n^2}{(n+1)(2n+1)} \frac{1}{\xi^2} \left\{1 - |1-\xi|^{{(n+1)}/{n}} \left(1- \xi\right) \right\} \right. \nonumber\\ &\quad \left.- \frac{n}{(n+1)}\frac{1}{\xi}|1-\xi|^{{(n+1)}/{n}}\right]. \end{align}
Now, far from the leading and trailing edges of the drop the curvature is negligible, giving rise to a negligible contribution of the capillary pressure term in the total pressure. Not only that, the term 
$\mathcal {K}$ in (2.23) turns out to be negligible due to a small value of the equilibrium contact angle 
$\theta _e$. This indicates a diminishing disjoining–conjoining pressure term in (2.22). Consequently, the limit 
$\xi ={hp'}/{\beta } \to 0$ is reached, which is also verified from the numerical simulations. Thus, a binomial expansion of 
$|1-\xi |^{{(n+1)}/{n}}$ allows us to write the asymptotic form of (2.21) in the limit 
$\xi \to 0$ as
\begin{equation} \frac{\partial h}{\partial t} + \frac{\partial}{ \partial x} \left(\mathcal{M}^{1/n-1}\beta^{1/n}\frac{h^2}{2} \right) = 0, \end{equation}signifying a Marangoni stress dominated spreading regime. Following Gomba & Homsy (Reference Gomba and Homsy2010), a self-similar solution of the above equation can be derived as
\begin{equation} h=\frac{x-x_0}{\mathcal{M}^{1/n-1}\beta^{1/n} t}, \end{equation}
where 
$x_0$ is the initial position of the rear contact angle. Thus, in this linear regime the droplet height can be represented using a scaled variable 
$h \beta ^{1/n}/\mathcal {M}^{1-1/n}$ that superimposes the dynamics for different power-law fluids and thermocapillary numbers, as shown in figure 3(b). The figure also shows that, with the increase in time, the drop height reduces, further indicating the dominance of the linear profile for a larger portion of the drop volume.
It is also inferred from (4.3) that a point on the free surface with height 
$h$ has a velocity 
$U= {\partial x}/{\partial t} = \mathcal {M}^{1/n-1}\beta ^{1/n} h$, highlighting the effect of the power-law index on the drop speed. The asymptotic free surface profile 
$h(x,t)$ depicted by (4.3) is shown with dashed lines for all the cases shown in figure 3(a). It is prominent that the accuracy of the asymptotic prediction spans a wider portion of the drop with increasing values of the power-law index 
$n$. This indicates a long-lasting contribution of the capillary pressure for the shear-thinning fluids compared with other fluids and explains the greater resistance to shape change offered by shear-thinning fluids.
The slope of the droplet in the asymptotic regime 
$(\mathcal {S})$ can be approximated from (4.3) as 
$\mathcal {S} = \partial h/\partial x = \mathcal {M}^{1-1/n} \beta ^{-1/n} t^{-1}$. The slope 
$\mathcal {S}$ for the non-Newtonian drops, normalized by the slope of a corresponding Newtonian drop at the same time 
$(\mathcal {S}_N)$, has been shown as a function of both the thermocapillary number 
$\beta$ and power-law index 
$n$ in figure 3(c). It shows that increasing the strength of the Marangoni stress, quantified by increasing 
$\beta$, causes a decrease in slope, representing the higher spreading rates with high thermocapillary strength for the same value of 
$n$. On the other hand, a decreasing slope predicts reducing values of the maximum free surface height 
$h$ with increasing 
$\beta$, a finding in accordance with the qualitative observations related to Newtonian Marangoni films (Kalliadasis et al. Reference Kalliadasis, Kiyashko and Demekhin2003).
A triangular approximation of the drop area 
$\mathcal {A}$ in the bulk region yields
\begin{equation} \mathcal{A}=\frac{w (x_f-x_r)}{2\mathcal{M}^{1/n-1}\beta^{1/n}t}= \frac{w^2}{2\mathcal{M}^{1/n-1}\beta^{1/n}t}, \end{equation}
with 
$x_f$ and 
$x_r$ being the locations of the front (or advancing) and rear (or receding) contact lines, and 
$w$ denoting the base width of the drop, defined as 
$w=x_f-x_r$. From (4.4) the width can be expressed as
\begin{equation} w= C(\beta, \mathcal{A}, \mathcal{M}, n) t^{\alpha}, \end{equation}
where 
$C= (2\mathcal {A}\mathcal {M}^{1/n-1}\beta ^{1/n})^{1/2}$ and 
$\alpha =1/2$ are the prefactor and the spreading exponent for droplet spreading. The asymptotic expression for 
$w(t)$ is observed to match well with the actual numerical data, as shown in figure 4(a). The functional form of 
$w(t)$ suggests that in the thermocapillary regime the exponent for spreading 
$(\alpha )$ does not depend on the power-law index, whereas the prefactor 
$(C)$ is a function of 
$n$. This observation is in contrast with the spreading characteristics of power-law fluids in the capillary regime (King Reference King2001; Starov et al. Reference Starov, Tyatyushkin, Velarde and Zhdanov2003; Rafaïet al. Reference Rafaï, Bonn and Boudaoud2004) where the spreading exponent also depends on 
$n$. This contrast is a consequence of the fact that the thermocapillary-driven spreading is intrinsically different from the capillary regime of spreading, distinguished by the asymmetry in the drop shape about a central plane, even for a two-dimensional droplet.
Figure 4. (a) Time evolution of the droplet base width 
$w$ for different 
$n$. The dotted lines indicate the corresponding asymptotic curve in each case obtained from (4.5). Other parameters are 
$\mathcal {A}=10,\beta =0.007$ and 
$\theta _e=5^{\circ }$. (b) Contour map for the droplet foot width at 
$t=2.1\times 10^5$ normalized with respect to its initial value 
$(w/w_0)$ on the 
$n$–
$\beta$ plane. Time evolution of the velocity of the peak height (c) and its axial position (d), for a different power-law index 
$n$. Insets show time variation of 
$h_{peak}$ (c) and 
$x_{peak}$ (d). A time interval 
${\rm \Delta} t= 1000$ has been used in the calculations of the peak velocity components from the numerical results.
In order to unveil the time scalings of a non-Newtonian, non-isothermal spreading drop, we have plotted the time evolution of the droplet base width 
$w$ for changing values of the parameter 
$n$ in figure 4(a). Unlike the thermocapillary regime, where we have derived 
$\alpha = 1/2$, the capillary regime spreading characteristics were reported as 
$w \sim t^\alpha$ with 
$\alpha = n/(5+2n)$ (Starov et al. Reference Starov, Tyatyushkin, Velarde and Zhdanov2003). This law dictates the spreading in early times where the linear profile assumption is not applicable, and the drop curvature effect has a significant role to play. In this connection, figure 4(a) depicts that the critical time for the onset of the large time behaviour (4.5) is higher for shear-thinning fluids, signifying prolonged capillary-dominated spreading characteristics. Conversely, the pure thermocapillary dominance commences much earlier for shear-thickening drops. The spreading extent of the drop has been quantified by the droplet foot width normalized with respect to its initial value, i.e. 
$w/w_0$ (obtained from (4.5)), is shown in figure 4(b). The increasing nature of 
$w/w_0$ with both 
$\beta$ and 
$n$ reflects that both the increasing thermocapillary intensity and the shear-thickening nature induce greater speeds to the advancing contact line 
$x_f$ than the receding contact line 
$x_r$.
The velocity of the peak point of the free surface profile 
$(x_{peak},h_{peak})$ sheds light onto the time dependence of the spreading. The velocity components of the peak point 
$(U_{peak},V_{peak})$, as delineated in figures 4(c) and 4(d), have been approximated as differences between successive elements of the 
$x_{peak}(t)$ and 
$h_{peak}(t)$ data obtained from the numerical simulations. A time interval 
${\rm \Delta} t= 1000$ has been used in these calculations. The initial rise in peak (see inset of figure 4c) signifies a critical competition between the capillary-dominated spreading and viscous effect within the early time zone. As time increases, the thermocapillary effect becomes significant and the capillary driving force is suppressed. In effect, a corresponding decrease in the peak position is observed. The negative values of vertical peak velocity 
$V_{peak}$ for different fluids physically signify a falling droplet height with time as it spreads along the surface. Finally, the peak saturates to a constant height, represented by its constant velocity 
$(U_{peak},V_{peak})$ after a certain time of spreading. For shear-thickening fluids, the initial viscous prominence is more significant due to increasing velocity gradients, and the competitive behaviours of the viscous and thermocapillary effects are observed up to a much longer time than that of the Newtonian or shear-thinning fluids. This is reflected in the sharp gradients of 
$U_{peak}(t)$ and 
$V_{peak}(t)$ curves till much longer times for shear-thickening fluids. These sharp changes mark the transition from the early time capillary regime to the large time thermocapillary regime (Chaudhury & Chakraborty Reference Chaudhury and Chakraborty2015). In this connection, a qualitatively similar sharp transition was reported earlier by Ren et al. (Reference Ren, Trinh and Weinan2015) for contact line velocity in the capillary spreading of a Newtonian drop.
4.2.2. Onset of new Marangoni film regimes
In addition to causing characteristic modulations to the drop profiles within the Marangoni film regime itself, the shear-dependent rheology is capable of instigating Marangoni films in place of a different spreading regime observed for Newtonian fluids. One such case is demonstrated in figure S-2 of the supplementary material, which depicts the gradual development of the profile of a shear-thickening drop with 
$n =1.4$ into a Marangoni film from a transition regime spreading of drops having 
$n \lesssim 1.1$ (refer to figure 2a).
Again, comparing the regime maps in figure 2(a–c), it is found that the Marangoni film regime emerges at a high value of the power-law index 
$(n_{cr})$, which is dependent on the parameter 
$\beta$. On the other hand, the Marangoni film zone is restricted upto a lower range of 
$\theta _e$ for lower values of 
$\beta$. For example, the Marangoni zone has been observed till 
$\theta _e=5^{\circ }$ for 
$\beta = 0.001$ and 
$0.007$, however, the same zone spreads upto 
$\theta _e = 10^{\circ }$ for 
$\beta = 0.03$. The physics responsible for such behaviour can be explained with due consideration that the effect of disjoining pressure is suppressed to a greater extent for high 
$\beta$. At the same time, the competition between the disjoining pressure and the Marangoni effect is severely influenced by the shear dependence of the viscosity, giving rise to critical power-law indices 
$n_{cr}$.
4.3. Non-Newtonian effects in uniform droplet motion
For a high contact angle value, the curvature effect is significant, and the droplet profile can no longer be approximated as a Marangoni film. This can be perceived as the manifestation of a strengthened disjoining–conjoining pressure at high contact angles that competes with the Marangoni stress in determining the curvature. Figure 5(a) shows the droplet profiles at 
$t=3\times 10^4$ with 
$\theta _e=40^{\circ }$ for different power-law indices 
$(n)$ and thermocapillary numbers 
$(\beta )$. The numerical simulations in this parametric zone predict that the droplet advances with a constant velocity 
$U$ and retains almost the same shape during the movement. Such a condition implies equal velocities of the drop's advancing and receding contact lines. We show that the non-Newtonian rheology strongly influences the uniform droplet migration speed, as evident from the different locations attained for fluids having different power-law indices. Compared with a Newtonian fluid, the uniform droplet speed 
$U$ enhances for shear-thickening fluids, while a speed reduction occurs for shear-thinning fluids. It is to be noted that the analytical approximation of drop speed derived earlier for Newtonian fluids (Brzoska et al. Reference Brzoska, Brochard-Wyart and Rondelez1993; Pratap et al. Reference Pratap, Moumen and Subramanian2008) is not feasible in the present scenario owing to the nonlinearity associated with the power-law model.
Figure 5. (a) Axial position 
$(x)$ versus height of the droplet 
$(h)$ at 
$t=3\times 10^4$ for different power-law indices 
$n$. The dotted and solid lines correspond to the cases with 
$\beta =0.01$ and 
$\beta = 0.007$, respectively. (b) Variation of the constant droplet velocity 
$U$ with 
$n$ for a different thermocapillary number 
$\beta$. Other parameters are 
$\mathcal {A}=10$ and 
$\theta _e=40^{\circ }$.
An intricate interplay between the thermocapillary effect and the shear-dependent viscous effects in deciding the droplet migration speed is portrayed in figure 5(b). It shows that the rate of increase in 
$U$ with 
$n$ is higher for 
$n<1$, while in the region 
$n>1$, the droplet gradually reaches a constant velocity 
$U$. In addition, for high 
$\beta$, the velocity has an overall increase, but it saturates faster than in lower 
$\beta$ cases. Such observations can be substantiated by considering that the fluid rheology affects the mechanisms of disjoining pressure and the Marangoni stress differently for varied strengths of the parameter 
$\beta$, and the resultant effect is reflected accordingly.
4.3.1. Onset of new droplet regimes
For some high values of 
$\beta$ (e.g. 
$\beta =0.03, \theta _e = 30^{\circ }$ in figure 2c), the droplet shows uniform motion for some shear-thinning fluids, while the transition regime sets in at a critical high value of 
$n_{cr}$, which is again dependent on the equilibrium contact angle 
$\theta _e$. Figure S-3 of the supplementary material shows one such instance. In this case, the drop peak initially falls, but retains a steady shape after a certain time. Here, the early time behaviour is greatly influenced by the intermolecular force at the drop–substrate interface and the drop shape tends to flatten. However, with an increase in time, the thermocapillary force is modified by the shear-thinning effect to such an extent that the intermolecular force is overpowered. Consequently, the destabilization of the thin film at the interface is suppressed, and the droplet shape remains intact.
The regime maps in figure 2 illustrate a common trend that for all the 
$\beta$ values considered, decreasing the power-law index has the effect of pushing the drop into the ‘droplet’ regime. It is observed from figure 2(a) that beyond the critical equilibrium contact angle 
$\theta _{cr}=10^{\circ }$ the droplet regime is the only regime that can occur with 
$\beta =0.001$. However, this critical angle becomes 
$\theta _{cr} = 20^{\circ }$ for 
$\beta =0.007$, and 
$\theta _{cr} > 40^{\circ }$ for 
$\beta =0.03$, as depicted by figures 2(b) and 2(c), respectively. This observation suggests that the dominance of the disjoining pressure in preserving the drop shape is diversely intervened by the fluid rheology depending on the strength of the Marangoni stress.
4.4. Non-Newtonian effects in the transition zone
It has been discussed in the preceding sections how the complex rheology can influence the different governing mechanisms responsible for the Marangoni film or the droplet regime of spreading. Now, it is found from figure 2(a–c) that the above two spreading regimes are separated by another spreading zone, defined here as the ‘transition’ regime. This spreading regime is characterized by rupturing of the main droplet body into multiple sister droplets or the advancement of the drop with an irregular shape without rupture, which cannot be classified either as a ‘droplet’ or ‘Marangoni film’. In this intermediate regime both the Marangoni stress and the disjoining–conjoining pressure terms turn out to be of comparable magnitude. As a result, the droplet eventually ruptures and splits into multiple smaller droplets that follow the parent droplet body. This being the case for a Newtonian droplet, the shear-dependent rheology modulates both these mechanisms in a complex and interconnected fashion for non-Newtonian droplets. Figure 6(a) describes the time evolution of a shear-thickening droplet for 
$\theta _e = 20^{\circ }$ and 
$\beta = 0.007$ in the transition regime where the droplet ruptures into multiple smaller droplets.
Figure 6. Time evolution of the droplet in the transition regime (a) 
$\theta _e= 20^{\circ }, \beta =0.007$ and 
$n =1.6$; (b) 
$\theta _e= 5^{\circ }, \beta =0.001$ and 
$n =0.6$; and (c) 
$\theta _e= 5^{\circ }, \beta =0.001$ and 
$n = 1$.
Under certain circumstances, the shear-thinning fluids can also show transition behaviour (see figure 6b), which is characteristically different from the corresponding transition behaviour for a Newtonian drop (refer to figure 6c). These figures also suggest that at low 
$\beta$ and 
$\theta _e$ rupturing of the droplet is delayed for a shear-thinning fluid. Notably, the drop behaves like a single entity with irregular shape for 
$n = 0.6$ even upto a large time 
$t= 6 \times 10^6$, as elucidated in figure 6(b). For Newtonian fluids the transition regime commences at a fixed intermediate contact angle 
$(\theta _{tr})$ when the thermocapillary number 
$(\beta )$ is held fixed. However, the non-Newtonian droplets show a more involved transition characteristic such that the critical transition contact angle 
$\theta _{tr}$ depends on the power-law index 
$n$. Figure 2(a) shows that, for a low thermocapillary number 
$\beta = 0.001$, the transition region appears in the phase map for a narrow range of contact angles in the zone 
$\theta _e \lesssim 10^{\circ }$. In comparison to this, the high 
$\beta$ cases in figures 2(b) and 2(c) show broader ranges of 
$\theta _e$, resulting in transition behaviour, ranging upto 
$\theta _e \approx 40^{\circ }$ in the last case. It is also revealed from the regime map in figure 2(a) that, for 
$\theta _e=5^{\circ }$, the droplet eventually ruptures till 
$n \lesssim 1.1$ instead of developing Marangoni films. Thus, in this zone, the effect of the power-law index on Marangoni stress plays a minor role compared with its effect on the disjoining pressure. On the other hand, for 
$n > 1.1$, the shear-thickening effect makes the Marangoni stress so strong that the droplet mimics the behaviour of a Marangoni film.
4.4.1. Onset of new transition regimes
Figures 7(a) and 7(b) highlight two contrasting spreading characteristics with varying power-law index 
$n$ for moderate 
$(\beta =0.007)$ and high 
$(\beta = 0.03)$ thermocapillary numbers, respectively. While, for 
$\beta =0.007$, the emergence of transition characteristics with shear-thickening fluids is preceded by droplet characteristics with lower 
$n$ values, the 
$\beta =0.03$ case exhibits a transition behaviour for all ranges of 
$n$ considered. Our simulations disclose that in this parametric zone the shear-thickening fluids exhibit a greater tendency of droplet rupture than the other fluid types. This is reflected in figure 7(b), which shows a greater number of ruptured droplets for 
$n = 1.4$ than 
$n =1$ and 0.6 cases.
Figure 7. Axial position 
$(x)$ versus height of the droplet 
$(h)$ for different power-law indices 
$n$, at (a) 
$\beta = 0.007$ and (b) 
$\beta = 0.03$. Other parameters are 
$t=4\times 10^4, \mathcal {A}=10$ and 
$\theta _e=20^{\circ }$.
Another type of new transition characteristics emerge in figure 2(c) for shear-thinning fluids 
$( 0.4 \lesssim n \lesssim 0.7 )$ with 
$\theta _e = 10^{\circ }$ and 
$\beta = 0.03$. It shows that 
$n$ values beyond this range facilitate the development of Marangoni films. Here, we find a tendency of early droplet rupture near the trailing edge as the parameter 
$n$ is lowered in the zone 
$0.4 \lesssim n \lesssim 0.7$. An example of this behaviour is portrayed in figure S-4 of the supplementary material. This observation is in sharp contrast with the pronounced ‘transition’ tendency of shear-thickening fluids as observed in the high contact angle zone 
$\theta _e \gtrsim 25^{\circ }$ in figure 2(c). These opposing effects of power-law index 
$n$ on the onset of transition regime for a distinct parameter set 
$(\theta _e,\beta )$ reiterates the contrasting influence of fluid rheology on the competing mechanisms, namely the conjoining–disjoining pressure and the Marangoni stress.
5. Conclusions and remarks
We have developed a mathematical model to explore the rich interfacial dynamics of a thermal gradient-driven sessile drop having shear-dependent viscosity, relevant to diverse industrial and biological flow scenarios. The momentum equation coupled with the incomprehensibility condition has been solved using the lubrication approximation, along with a conjoining–disjoining pressure model to tackle the contact line singularity. The flow field has been shown to be a nonlinear composition of a power-law Poiseuille flow and a shear flow triggered by the interfacial driving force of Marangoni flow. Consequently, the numerical solution of the governing equations has been found to agree reasonably well with an asymptotic analytical solution derived in limiting scenarios.
A systematic analysis of the droplet profiles for varying droplet wettability, thermocapillary strength and the shear dependence of viscosity enabled us to categorize the droplet spreading features into three distinct regimes, namely (i) the Marangoni film regime, where the droplet attains a linear profile similar to a capillary ridge; (ii) the droplet regime, where the droplet migrates at a constant speed with almost the same shape; and (iii) the transition regime, where the droplet either ruptures or moves with an irregular shape.
The reported results related to the Marangoni film regime reveal that at a low value of the equilibrium contact angle, the Marangoni effect overpowers the weakened disjoining pressure, resulting in negligible curvature of the droplet. In addition, it has been inferred from the regime maps for both low and high 
$\beta$ (figure 2a,c) that a high value of the power-law index is responsible for a greater tendency for the conversion to the Marangoni film regime from the transition regime. It has been found that a shear-thickening droplet has a steeper advancing front, whereas its shear-thinning counterpart retains a significant curvature for a much longer time. Moreover, a shear-thinning (thickening) droplet spreads at a slower (faster) rate than a Newtonian droplet. This behaviour has been linked to opposing causes of flow acceleration at the contact line area and contact angle reduction. We also discovered that the capillary- and thermocapillary-driven spreading mechanics are prominent at various time scales for the non-Newtonian fluids studied. The thermocapillary-dominated spreading behaviour sets in at a much earlier time for shear-thickening fluids (figure 4a). Furthermore, the time variation of the droplet peak (figure 4c,d) unravels that the shear-thickening fluid property makes the initial viscous stresses a dominant participant till a much longer time during its interaction with the thermocapillary stresses, causing sharp gradients in the time evolution of peak point velocity.
The droplet regime occurs at a high equilibrium contact angle 
$\theta _e$. In this regime a strengthened disjoining–conjoining pressure prevents the droplet from attaining a linear shape, and the droplet migrates at a uniform speed 
$U$ with almost the same shape. The regime maps illustrate that decreasing the power-law index has the effect of retaining the drop within the ‘droplet’ regime. In effect, the onset of this regime does not commence at a fixed value of 
$\theta _e$ when the rheology deviates from that of a Newtonian fluid. While some shear-thinning fluids cause uniform droplet motion, the transition regime sets in beyond a critical power-law index. While certain shear-thinning fluids induce uniform droplet motion, the transition regime kicks in after a threshold power-law index. Not only this, the minimum equilibrium contact angle 
$(\theta _{cr})$ required for preserving the droplet shape shows a steady increase with 
$\beta$. For example, 
$\theta _{cr} = 10^{\circ },20^{\circ }$ and 
$40^{\circ }$ with 
$\beta =0.001, 0.007$ and 
$0.03$, respectively. We have further highlighted an enhancement (attenuation) of 
$U$ in the 
$n>1$ (
$n<1$) cases (figure 5b). An intricate interplay between the thermocapillary effect and the shear-dependent viscous effects decides the droplet speed. For high values of 
$\beta$, the droplet velocity has an overall increase, but it gradually reaches a constant value at a faster rate. On the other hand, 
$U$ increases with the power-law index 
$n$. The rate of increase of 
$U$ with 
$n$ is high for shear-thinning fluids, while saturation of 
$U$ has been observed with the shear-thickening rheology.
In the intermediate transition regime, the Marangoni stress and the disjoining pressure terms are comparable. The nonlinear interaction of these effects in the case of non-Newtonian fluids greatly influences the drop spreading characteristics. In conditions with low 
$\beta$ and 
$\theta _e$ values, a delayed rupturing of the droplet takes place for shear-thinning fluids (figure 6b). The regime maps show that the transition zone is extended for a larger span of 
$\theta _e$ with an increase in 
$\beta$. In addition, an overall higher tendency of shear-thickening fluids has been revealed due to high chances of droplet rupture. In contrast, for some low 
$\theta _e$ and high 
$\beta$ combinations, the rising power-law indices do not ensure transition characteristics, rather the Marangoni film is encountered (see, e.g. figure S-4).
To summarize, we have unraveled the unique characteristics of the different spreading regimes that are engendered due to a highly nonlinear impact of the non-Newtonian rheology on the competitive effects of Marangoni stress and conjoining–disjoining pressure. It is worth noting that the results of the power-law model should be used cautiously when the shear rate is low. For a more comprehensive study in the future, focusing on the shear rate dependent fluid rheology, Carreau or other such models (Bird et al. Reference Bird, Armstrong and Hassager1987; Myers Reference Myers2005) may be employed to overcome the shortcomings of the power-law model by incorporating finite zero and infinite shear rate viscosities. Furthermore, an enhanced understanding of other non-Newtonian effects, for example, normal stresses in fluids (Rafaïet al. Reference Rafaï, Bonn and Boudaoud2004) remains an exciting arena to explore and may complement the different aspects of the present work. We conjecture that the incorporation of gravity effects in the present problem set-up may have interesting consequences like enhanced spreading rates (Karapetsas et al. Reference Karapetsas, Sahu, Sefiane and Matar2014), elongation of droplets (Brzoska et al. Reference Brzoska, Brochard-Wyart and Rondelez1993), etc.
The present model stands as a first attempt to identify the complicated interplay between the thermocapillary effects and the non-Newtonian rheology on droplet spreading, considering the two-dimensional droplet as the representation of the symmetry plane passing through the maximum thickness of the original three-dimensional drop, similar to a host of earlier works on Newtonian droplets (Ehrhard & Davis Reference Ehrhard and Davis1991; Gomba & Homsy Reference Gomba and Homsy2010; Sui Reference Sui2014; Chaudhury & Chakraborty Reference Chaudhury and Chakraborty2015; Mac Intyre et al. Reference Mac Intyre, Gomba, Perazzo, Correa and Sellier2018). However, inferences drawn from the current study should be used with caution when comparing to experiments due to the assumed two-dimensional shape in simulations, whereas actual drops are three-dimensional. A two-dimensional droplet model was found to describe the experimental results with great accuracy for the uniform droplet regime (Chen et al. Reference Chen, Troian, Darhuber and Wagner2005), supporting our considerations. Conversely, due to the presence of another radius of curvature in three dimensions, the capillary stress is expected to be different for three-dimensional drops, particularly in the Marangoni film regime. In addition, the droplet rupture in two- and three-dimensional scenarios could differ, indicating the need for detailed three-dimensional simulations to validate the droplet characteristics in the transition regime. We intend to address these complicated issues arising in a three-dimensional droplet in future studies.
Given the flexibility of tuning both the shear-thinning and shear-thickening properties of polymeric fluids by using different combinations of the base fluid and the solute (Wang et al. Reference Wang, Zhang, Lee and Peng2007b), the current theoretical insights may open up new avenues for experimental exploration. Besides, the reported distinct droplet behaviours due to shear-dependent fluid viscosity may be correlated for accurate extraction of the complex rheological properties of the working fluid. Finally, the technological advancements in thermal gradient-based control in droplet microfluidics (Nguyen, Pang & Huang Reference Nguyen, Pang and Huang2006; Won, Lee & Song Reference Won, Lee and Song2017) may be coupled with the fundamental understanding of the complex rheological responses of the spreading droplets developed in the present work toward designing efficient microfluidic systems used for enhanced manipulation of complex biofluids.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2022.900.
Acknowledgements
The authors would like to acknowledge J.R. Mac Intyre (Aalto University, Finland) for useful discussions.
Funding
A.P. acknowledges the support provided by the Indian Institute of Technology (ISM), Dhanbad though the project: FRS(169)/2021-2022/ME.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Detailed derivation of flow velocity and flow rate
The two-dimensional Cauchy momentum equation (2.1a) takes the form
\begin{equation} \frac{\partial \tilde{\tau}_{xy}}{\partial \tilde{y}}=\frac{\partial \tilde{p}}{\partial \tilde{x}}. \end{equation}
Integrating with respect to 
$\tilde {y}$ from 
$\tilde {y}=\tilde {h}$ to 
$\tilde {y}$ gives
\begin{equation} \tilde{\tau}_{xy}\left|_{\tilde{y}}- \tilde{\tau}_{xy}\right|_{\tilde{h}}=\tilde{p}_x \tilde{y} - \tilde{h} \tilde{p}_x. \end{equation}
Since 
$\tilde {\tau }_{xy}|_{\tilde {h}} = \hat {\tau }$ (2.17), we can write
\begin{equation} \tilde{\tau}_{xy}|_{\tilde{y}}=\tilde{p}_x \tilde{y} + c_1,\quad \text{where}\ c_1={-}\tilde{h} \tilde{p}_x+\hat{\tau}. \end{equation}Applying the lubrication approximation, the power-law constitutive equation (2.2) can be used in the above equation to yield
$$\begin{gather} \mu_0 \tilde{\mathcal{M}}^{n-1} \left|\frac{\partial \tilde{u}}{\partial \tilde{y}}\right|^{n-1} \frac{\partial \tilde{u}}{\partial \tilde{y}}=\tilde{p}_x \tilde{y}+c_1 \end{gather}$$
$$\begin{gather}\text{or}\quad \frac{\partial \tilde{u}}{\partial \tilde{y}}=\frac{\tilde{\mathcal{M}}^{1/n-1}}{\mu^{1/n}_0}\left|\tilde{p}_x \tilde{y}+c_1\right|^{{1}/{n}-1}(\tilde{p}_x \tilde{y}+c_1). \end{gather}$$
Again integrating from 
$\tilde {y}=0$ to 
$\tilde {y}$ and using the no-slip condition lead to
\begin{align} \tilde{u} &= \frac{\tilde{\mathcal{M}}^{1/n-1}}{\mu^{1/n}_0} \int_{\tilde{y}=0}^{\tilde{y}}\left|\tilde{p}_x \tilde{y}+c_1\right|^{{1}/{n}-1}(\tilde{p}_xy+c_1)\,{\rm d} \tilde{y} \end{align}
\begin{align} &= \frac{\tilde{\mathcal{M}}^{1/n-1}}{\mu^{1/n}_0} \left[\frac{\left|\tilde{p}_x \tilde{y}+c_1\right|^{{1}/{n}+1} - \left|c_1\right|^{{1}/{n}+1}}{\left(\frac{1}{n}+1\right) \tilde{p}_x} \right] \end{align}
\begin{align} &= \frac{\tilde{\mathcal{M}}^{1/n-1}}{\tilde{p}_x\mu^{1/n}_0} \frac{n}{n+1} \left[\left|\tilde{p}_x\tilde{y}+c_1\right|^{{1}/{n}+1} - \left|c_1\right|^{{1}/{n}+1} \right]. \end{align}Thus, the expression for flow rate becomes
\begin{align} \tilde{Q}= \int_{0}^{\tilde{h}} \tilde{u}\,{\rm d}\tilde{y} &= \frac{\tilde{\mathcal{M}}^{1/n-1}}{\mu^{1/n}_0} \frac{n^2}{(n+1)(2n+1)} \frac{1}{\tilde{p}_x^2} \left\{ \hat{\tau}^{{1}/{n}+2} - \left|\hat{\tau}-\tilde{h} \tilde{p}_x\right|^{{1}/{n}+1} \left(\hat{\tau}- \tilde{h} \tilde{p}_x\right) \right\} \nonumber\\ &\quad -\frac{\tilde{\mathcal{M}}^{1/n-1}}{\mu^{1/n}_0} \frac{n}{(n+1)} \frac{1}{\tilde{p}_x} \left\{\left|\hat{\tau}-\tilde{h} \tilde{p}_x\right|^{{1}/{n}+1} \tilde{h}\right\}. \end{align}
