1. Introduction
Droplets suspended and transported in flow media are central in several emerging applications encompassing biomedical technology (Gañán-Calvo et al. Reference Gañán-Calvo, Montanero, Martín-Banderas and Flores-Mosquera2013; Mazutis et al. Reference Mazutis, Gilbert, Ung, Weitz, Griffiths and Heyman2013; Pan et al. Reference Pan, Mukherjee, Rahaman, Samanta and Dasgupta2013; Tasoglu et al. Reference Tasoglu, Gurkan, Wang and Demirci2013; Vladisavljević et al. Reference Vladisavljević, Khalid, Neves, Kuroiwa, Nakajima, Uemura, Ichikawa and Kobayashi2013; Zhao et al. Reference Zhao, Chen, Yue, French, Rufo, Benkovic and Huang2013; Rosenfeld et al. Reference Rosenfeld, Lin, Derda and Tang2014; Schlicht & Zagnoni Reference Schlicht and Zagnoni2015), industrial processes (Karabelas Reference Karabelas1977; Kaushal & Tomita Reference Kaushal and Tomita2002), flow cytometry (Adan et al. Reference Adan, Alizada, Kiraz, Baran and Nalbant2017) and beyond. As droplets migrate in confined fluidic pathways, they are likely to confront obvious variations in their thermal ambience, either purposefully imposed or inherent to the system itself. For instance, when droplets are used as drug carriers, they need to migrate to the specific target areas inside the human body by invading complex non-isothermal environments inside physiological pathways (Lawson & Chughtai Reference Lawson and Chughtai1963; Werner & Buse Reference Werner and Buse1988). However, in engineered flow-focusing experiments (Baroud et al. Reference Baroud, Delville, Gallaire and Wunenburger2007a; Baroud, de Saint Vincent & Delville Reference Baroud, de Saint Vincent and Delville2007b; Nguyen et al. Reference Nguyen, Ting, Yap, Wong, Chai, Ong, Zhou, Tan and Yobas2007; Stan, Tang & Whitesides Reference Stan, Tang and Whitesides2009), micro-fabricated thermal units are known to facilitate controlled droplet motion. As the droplets migrate, deformation of the droplet–carrier fluid interface is likely to trigger non-trivial alterations in the incipient flow field and a consequent alteration in the droplet morpho-dynamics via complex two-way thermo-fluid coupling, an aspect that may not be trivially captured by simplified theoretical frameworks. A related area of research is the thermo-capillary dynamics of sessile droplets (Karapetsas et al. Reference Karapetsas, Sahu, Sefiane and Matar2014; Pradhan & Panigrahi Reference Pradhan and Panigrahi2015) which is governed by three-phase contact line phenomena mediated by fluid motion inside the droplet. It should be noted that the physical problem addressed in this paper involves a droplet which is completely immersed in another fluid and experiences a dynamically evolving thermal field on its entire periphery, as opposed to a sessile droplet which has a footprint on the solid surface and hence, the thermal gradient imposed on the solid surface directly acts on the droplet. This renders the dynamics of droplets in the scenario of a freely moving droplet physically different as compared to that of a sessile droplet, and hence the physics of migration as observed for the latter cannot be trivially presumed to govern the former.
Motion of a droplet in a thermal field is known to be strongly influenced by the interfacial stress stemming from the gradient of surface tension with temperature, culminating in an interfacial tangential stress, alternatively interpreted as the Marangoni (thermo-capillary) effect (Subramanian Reference Subramanian1983). Starting from the early work of Young, Goldstein & Block (Reference Young, Goldstein and Block1959), several researchers have studied the thermally modulated dynamics of droplets and bubbles in a quiescent fluid medium, as governed by Marangoni effects. In a classical text titled ‘Film notes for surface tension in fluid mechanics’, Trefethen (Reference Trefethen1969) describes an experiment where air bubbles propel themselves within a capillary tube, along the direction of the externally imposed temperature gradient. They explained this phenomenon qualitatively by analysing the variation in the thermally induced surface tension of the bubble and the consequent fluid motion due to viscous drag. Nadim, Haj-Hariri & Borhan (Reference Nadim, Haj-Hariri and Borhan1990) reported thermo-capillary migration of a slightly deformed droplet in a quiescent fluid, where the droplet deformation occurs either due to surfactants or due to fluid inertia. They estimated the droplet migration speed by expanding the droplet shape using a regular perturbation approach and concluded that in the presence of droplet deformation, the migration speed can either increase or decrease, depending on the relevant physical property ratios of the carrier and suspending fluid phases. Later, Raja Sekhar and co-workers (Choudhuri & Raja Sekhar Reference Choudhuri and Raja Sekhar2013; Sharanya & Raja Sekhar Reference Sharanya and Raja Sekhar2015) studied the thermo-capillary migration of a spherical droplet in the presence of a background flow and inferred that the effects of non-uniform temperature field and background flow can be linearly superimposed. Recently, Das, Mandal & Chakraborty (Reference Das, Mandal and Chakraborty2018); Das et al. (Reference Das, Mandal, Som and Chakraborty2017) studied the migration of a surfactant-coated droplet in a non-isothermal pressure driven flow. They used a perturbation approach for studying the thermo-capillary migration in the regimes of diffusion-dominated and advection-dominated transport of surfactants along the droplet interface. It was found that the droplet migrates faster both along the axial and cross-stream wise directions when the temperature field is increasing along the direction of flow. Conversely, a linearly decreasing temperature field in the flow direction was found to slow down the droplet migration along both the axial as well as cross-streamline directions. In a related work by the same authors, the effect of asymptotically small shape deformation was also taken into account (Das & Chakraborty Reference Das and Chakraborty2018), and it was observed that when the temperature decreases along the direction of flow, a reversal in the direction of cross-streamline migration of the droplet occurs beyond a certain magnitude of the temperature gradient.
Despite offering significant insights on thermally modulated droplet transport, the above results appear to be physically somewhat restrictive because of the following simplifications made for analytical tractability: (a) shape deformation of the droplet during its dynamic evolution is neglected (Das et al. Reference Das, Mandal and Chakraborty2018) or at the best considered to be asymptotically small (Das & Chakraborty Reference Das and Chakraborty2018); and (b) the externally imposed thermal field is directly super-imposed on the droplet, obviating a possible two-way coupling mediated by the balance of viscous stress and capillary stress at the dynamically evolving interface.
Here, we show that by relaxing the above assumptions via an explicit interface-tracking approach, one may confront a physical scenario that conflicts with results from reported theories where only one-way coupling between the externally imposed thermal field and the interfacial stress field is considered. Towards establishing this proposition, we use a boundary integral based formalism that essentially converts the physical three-dimensional problem into a set of two-dimensional integral equations over the droplet surface. For illustration, we consider a droplet initially suspended in an externally imposed parabolic flow (Griggs, Zinchenko & Davis Reference Griggs, Zinchenko and Davis2007; Janssen & Anderson Reference Janssen and Anderson2008) and a linearly varying temperature gradient along the direction of the flow. In such a scenario, the boundary integral equations can further be simplified into a set of integrals over the droplet interface by including additional corrections to the Green's function (Pozrikidis et al. Reference Pozrikidis1992). Nevertheless, the transport phenomenon addressed herein is elusively more complex as against an intuitive two-dimensional paradigm, as attributed to an inherent three-dimensionality of the dynamically evolving droplet–carrier fluid interface. Consequently, a three-dimensional consideration of the interface is warranted to bring out the correct interfacial forces.
As a decisive advantage, as opposed to the established diffuse interface methods that artificially smear out macroscopically sharp interfaces to gain computational advantages (Anderson, McFadden & Wheeler Reference Anderson, McFadden and Wheeler1998; Liu & Zhang Reference Liu and Zhang2010; Teigen et al. Reference Teigen, Song, Lowengrub and Voigt2011) of avoiding mesh entanglement during the interface deformation at the expense of using soft parameters in the model that cannot be directly mapped with the physical data, our approach enables an explicit representation of the interface without invoking any artificial model parameters, and hence holds the capability of extremely accurate depiction of the interfacial forces. Further, our boundary integral based formulation is established to be capable of representing a three-dimensional interface via a derived system of two-dimensional equations, establishing computational tractability amidst the dynamical complexity. In contrast to the predictions from analytical theory, we show that the axial and cross-stream migration speeds of a droplet increase in the presence of a linearly decreasing temperature field and conversely, an increasing temperature field slows down the migration velocity along both axial and cross-stream directions. This paradigm of two-way coupling between the imposed flow and other field mediated interfacial effects is likely to be decisive in a wide class of problems encompassing droplet and bubble dynamics, offering interfacial descriptions with high level of quantitative accuracy as compared to established theories.
2. Problem formulation
2.1. System description
We consider a Newtonian droplet with density 
$\rho$, viscosity 
$\mu$ and thermal conductivity 
$k_{i}$. Initially, it holds a spherical shape with a radius 
$a$. It is suspended in another Newtonian medium having identical density and viscosity, and a thermal conductivity equal to 
$k_{e}$. The surface tension, 
$\bar {\sigma }$, varies along the droplet–carrier fluid interface, due to local variations in temperature. For simplicity, yet without sacrificing the essential physics, we consider a linear decrease of 
$\bar {\sigma }$ with increase in the interfacial temperature 
$\bar {T}_{s}$, with the reference value of 
$\bar {\sigma }_{0}$ at temperature 
$\bar {T}_{0}$.
We fix up a Cartesian co-ordinate system (
$\bar {x}$,
$\bar {y}$,
$\bar {z}$) at the channel centreline, on the inlet boundary. Initially, the droplet centroid is positioned with a slight offset from the channel centreline at (
$\bar {X}_{c},\bar {Y}_{c},\bar {Z}_{c}$). The imposed free stream velocity, 
$\boldsymbol {\bar {V}}_{\infty }$, is taken to be a parabolic profile with centreline velocity 
$\bar {V}_{c}$ and the externally imposed temperature field is denoted by 
$\bar {T}_{\infty }$. The domain is bounded along the z-direction by two parallel plates separated by a distance 
$\bar {H}$. The velocity of the droplet centroid is denoted by 
$\bar {\boldsymbol {U}}$, having components equal to 
$\bar {U}_{x}$, 
$\bar {U}_{y}$ and 
$\bar {U}_{z}$ along the 
$x$, 
$y$ and 
$z$ directions, respectively. The temporal co-ordinate is denoted by 
$t$.
2.2. Normalization of the parametric space
Normalization of the parametric space is essential towards gaining a generic insight into the physical situation disregarding absolute values of the chosen physical parameters. Towards this, we choose the reference length scale of the system to be the un-deformed radius of the droplet (
$a$), the velocity scale as the centreline velocity 
$\bar {V}_{c}$ of the imposed flow, characteristic time scale to be 
$t_c$ and the characteristic temperature difference (
$\Delta T_{c}$) as the temperature difference across the length scale of the droplet (
$|\bar {G}|a$), where 
$|\bar {G}|$ denotes the axial gradient of the externally applied temperature field (see figure 1). The non-dimensional temperature field in the domain is expressed relative to the reference temperature 
$\bar {T}_{0}$. Finally, we choose the thermal conductivity of the carrier fluid 
$k_{e}$ to be the reference value for thermal conductivity. Accordingly, the following relevant dimensionless parameters emerge.
Figure 1. Schematic of the physical set-up where a deformable droplet is suspended in a pressure-driven flow in a parallel plate channel. For the purpose of illustration, we consider an imposed temperature profile decreasing along the flow direction. The model framework, however, is generic to accommodate any imposed variations in the thermal field.
(i) Capillary number
$Ca = {\mu \bar {V}_{c}}/{\bar {\sigma }_{0}}$, denoting the relative importance of viscous and surface tension forces.(ii) Marangoni number
$Ma_{T} = {\beta \overline {|G|}a|}/{\mu \bar {V}_{c}}$, denoting the magnitude of temperature-gradient driven interfacial flow as compared to the strength of the imposed flow.(iii) Thermal Péclet number
$Pe_{T} = {\bar {V}_{c}a}/{\alpha _{e}}$, denoting the relative magnitude of advective to diffusive transport of thermal energy. Here, $\alpha _{e}$ refers to the thermal diffusivity of the carrier phase.(iv) Grashof number
$Gr = {g\gamma _{e}\rho _{e}^{2}\Delta Ta^{3}}/{\mu _{e}^{2}}$, denoting the relative magnitude of thermal buoyancy to viscous force.(v) Reynolds number
$Re = {\rho _{e}\bar {V}_{c}a}/{\mu _{e}}$, denoting the ratio of the imposed flow velocity to a diffusive velocity scale due to viscous dissipation and the Strouhal number is denoted by $S = {t_c}/{(a/\bar {V}_{c})}$, such that the relative importance of Eulerian acceleration is given by the ratio $Re/S$.(vi) The ratio of the thermal conductivities of the suspending and carrier fluid phases,
$\delta$.(vii) Confinement ratio
$Wc = {2a}/{H}$, denoting the extent of occupancy of the droplet relative to the channel height.
We use the boolean variable 
$\zeta$ to denote the direction of the temperature gradient. For a linearly increasing temperature field, 
$\zeta = 1$ and for a linearly decreasing temperature field, 
$\zeta = -1$. After normalization, all the variables are denoted by 
$v$ instead of 
$\bar {v}$, where 
$v$ refers to any physical variable. Further, we denote the value of a physical variable inside the suspending phase using a subscript ‘i’ (
$v_{i}$) and corresponding value for the carrier phase using a subscript ‘e’ (
$v_{e}$). All vector quantities are denoted by 
$\boldsymbol {v}$ and second-order tensors are denoted using 
$\boldsymbol{\mathsf{v}}$.
2.3. Assumptions
The problem considered here is a complex multi-physics problem involving multi-phase fluid mechanics and thermal transport. The dynamically evolving deformable interface further complicates the solution by imposing a two-way coupling between the fluid flow and heat transfer mediated by a topologically evolving phase boundary. To simplify the solution procedure without sacrificing the essential physics of interest here, the following assumptions are made.
(i) The fluid flow is dominated by viscous effects (
$Re \ll 1$) and the magnitude of Eulerian acceleration is assumed to be negligible in comparison to viscous and pressure-gradient terms ($Re/S \ll 1$). These assumptions are quite standard in the literature of low-Reynolds-number hydrodynamics and are commonly referred to as the ‘creeping-flow’ approximation (Leal Reference Leal2007). A detailed justification of this assumption along with its applicability to the present problem is discussed in Appendix B.1.(ii) Thermal energy transport occurs predominantly via diffusion, so that the thermal Péclet number is negligibly small (
$Pe_{T} \ll 1$). This is consistent with several physical systems of practical interest (Nallani & Subramanian Reference Nallani and Subramanian1993; Chen et al. Reference Chen, Lu, Yang and Maa1997).(iii) Thermal buoyancy effects are not important, considering
$Gr \ll 1$.(iv) The surface tension is assumed to vary linearly with the interfacial temperature (Young et al. Reference Young, Goldstein and Block1959; Gittens Reference Gittens1969), so that
(2.1)where\begin{equation} \bar{\sigma} = \bar{\sigma}_{0} - \beta(\bar{T}_{s} - \bar{T}_{0}), \end{equation}$\beta$ refers to the sensitivity of the interfacial tension with change in temperature.(v) The domain is assumed to be unbounded along the
$x$ and $y$ directions. This physically replicates a set-up where the inlet and outlet of the channel lie far apart from the droplet surface and the walls bounding the fluid domain along the $y$ direction lie far apart from each other in comparison to the walls along the $z$ direction. Thus, the extent of the channel along the $z$ direction determines the effective role of the confinement.
2.4. Governing equations and boundary conditions
2.4.1. Thermal problem
In the absence of any heat source/sink within the domain and negligible thermal convection, the thermal problem is governed by the homogeneous Laplace equation, as given by
\begin{equation} \left. \begin{gathered} \displaystyle{\nabla}^{2}T_{i} = 0,\\ \displaystyle{\nabla}^{2}T_{e} = 0. \end{gathered} \right\} \end{equation}The ambient temperature field, which exists in the absence of the droplet, is given by
\begin{equation} T_{\infty} = \zeta z, \end{equation}
where 
$\zeta$ can have a value of 1 or 
$-1$ depending on whether the temperature increases or decreases along the direction of flow.
The temperature and heat flux are continuous across the interface of the droplet. This can be expressed mathematically as
\begin{equation} \left. \begin{gathered} \displaystyle T_{i} = T_{e}, \\ \displaystyle \delta(\boldsymbol{\nabla} T_{i}\boldsymbol{\cdot}\hat{\boldsymbol{n}}) = \boldsymbol{\nabla} T_{e}\boldsymbol{\cdot}\hat{\boldsymbol{n}} \end{gathered} \right\} \quad {\rm on} \ \varGamma, \end{equation}
where 
$\varGamma$ denotes the interface of the droplet and 
$\hat {\boldsymbol {n}}$ denotes the unit normal vector on 
$\varGamma$ (see figure 1).
The channel walls are kept at the ambient thermal conditions:
\begin{equation} \displaystyle T_{e} = T_{\infty} \quad {\rm on} \ B, \end{equation}
where 
$B$ denotes the rigid walls (see figure 1).
2.4.2. Hydrodynamics
As the Reynolds number of the flow is vanishingly small, the fluid flow is governed by the Stokes flow equations:
\begin{equation} \left.
\begin{gathered}
\displaystyle\mu{\nabla}^{2}\boldsymbol{u}_{i} =
\boldsymbol{\nabla} p_{i}, \quad
\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_{i} =
0,\\ \displaystyle\mu{\nabla}^{2}\boldsymbol{u}_{e} =
\boldsymbol{\nabla} p_{e}, \quad
\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_{e} =
0. \end{gathered} \right\}
\end{equation}
Here, 
$\boldsymbol {u}_{i}$ and 
$p_{i}$ refer to the non-dimensional velocity and pressure fields, respectively, inside the droplet. The corresponding variables in the carrier phase are denoted by 
$\boldsymbol {u}_{e}$ and 
$p_{e}$. The ambient velocity field in the absence of the droplet is given by
\begin{equation} \boldsymbol{V}_{\infty} = \left(1 - 4\frac{z^2}{H^2}\right)\hat{x}. \end{equation}The traction discontinuity at the interface is balanced by the surface tension forces. It can be expressed as a force balance at the interface:
\begin{equation} (\boldsymbol{\mathsf{S}}_{e} - \boldsymbol{\mathsf{S}}_{i})\boldsymbol{\cdot}\hat{\boldsymbol{n}} = \sigma\hat{\boldsymbol{n}}(\boldsymbol{\nabla}\boldsymbol{\cdot}\hat{\boldsymbol{n}}) - \boldsymbol{\nabla}_{s}\sigma , \end{equation}
where 
$\boldsymbol{\mathsf{S}}_{i}$ and 
$\boldsymbol{\mathsf{S}}_{e}$ denote respectively the hydrodynamic stress tensors for the fluid inside and outside the droplet. Here, 
$\boldsymbol {\nabla }_{s}$ denotes the gradient operator projected along the surface, i.e. 
$\boldsymbol {\nabla }_{s} \equiv (\boldsymbol{\mathsf{I}} - \hat {\boldsymbol {n}}\hat {\boldsymbol {n}})\boldsymbol {\cdot }\boldsymbol {\nabla }$.
The tangential component of velocity is continuous across the interface, which implies
\begin{equation} \boldsymbol{u}_{i} - (\boldsymbol{u}_{i}\boldsymbol{\cdot}\hat{\boldsymbol{n}})\hat{\boldsymbol{n}} = \boldsymbol{u}_{e} - (\boldsymbol{u}_{e}\boldsymbol{\cdot}\hat{\boldsymbol{n}})\hat{\boldsymbol{n}} \quad {\rm on}\ \varGamma. \end{equation}The evolution of the interface is governed by the normal component of the interfacial velocity:
\begin{equation} \frac{\mathrm{d}\kern0.06em \boldsymbol{x}_{\varGamma}}{\mathrm{d} t}\boldsymbol{\cdot}\hat{\boldsymbol{n}} = \boldsymbol{u}_{i}\boldsymbol{\cdot}\hat{\boldsymbol{n}} = \boldsymbol{u}_{e}\boldsymbol{\cdot}\hat{\boldsymbol{n}} \quad {\rm on} \ \varGamma, \end{equation}
where 
$\boldsymbol {x}_{\varGamma }$ denotes the position vector of any point on the interface.
Equation (2.10) is physically reminiscent of two important features, namely, the interface evolution and kinematic boundary condition. The first part of the equation denotes an expression for evolution of the droplet interface, which is akin to several previous approaches that have employed a Lagrangian interface tracking for simulation of deformable droplets (Loewenberg & Hinch Reference Loewenberg and Hinch1996). Since the shape of the droplet changes exclusively due to the normal component of the interfacial velocity, we have considered only the normal component of the interfacial velocity in the first part of (2.10). The second part of (2.10) represents the kinematic boundary condition under the assumption of no phase change (Leal Reference Leal2007). Summarily, in the absence of any phase change processes, (2.10) accurately describes both the kinematic boundary condition and evolution of the droplet interface. This could also be interpreted as continuity of velocity at the interface as (2.10) implies that at the interface, both the droplet and carrier phases have identical velocities which is equal to the velocity of the interface.
The velocity is assumed to satisfy the no-slip and no-penetration conditions on the surface of the rigid walls, i.e.,
\begin{equation} \boldsymbol{u}_{e} = \boldsymbol{0} \quad {\rm on} \ B. \end{equation}3. Numerical solution
3.1. Formulation of integral equations
We employ two separate three-dimensional boundary element methods (BEMs) for the solution of the thermal and the flow problems. The two methods are coupled via the interfacial force balance (2.8) considering the temperature dependence of surface tension (2.1). Whereas BEM has been extensively used to solve the isothermal fluid dynamics, here, for the first time, we explore the crucial capabilities of BEM in coupling the conservation equations of momentum and thermal energy. In the foregoing discussion, we first describe the BEM for the thermal problem and subsequently provide an outline of the BEM for capturing the flow dynamics.
A general solution of the energy equation within the domain can be constructed by a linear superposition of a set of point singularities on the domain boundaries. The strength of these point singularities is unknown a priori and is obtained from the solution of a Fredholm integral equation. For the temperature field inside the carrier phase, the general solution can be expressed as (Pozrikidis Reference Pozrikidis2002)
$$\begin{align}
T_{e}(\boldsymbol{x}_{0}) &= T_{\infty}(\boldsymbol{x}_{0})
-
\underbrace{\int_{B\cup\varGamma}G_{l}(\boldsymbol{x},\boldsymbol{x}_{0})\{
\boldsymbol{\nabla} T_{e}(\boldsymbol{x})
\boldsymbol{\cdot} \hat{\boldsymbol{n}}(\boldsymbol{x})
\} \,\mathrm{d}S(\boldsymbol{x})}_{Single \ Layer \
Operator} \nonumber\\ &\quad+
\underbrace{\int_{B\cup\varGamma}T_{e}(\boldsymbol{x})\{
\hat{\boldsymbol{n}}(\boldsymbol{x}) \boldsymbol{\cdot}
\boldsymbol{\nabla}
G_{l}(\boldsymbol{x},\boldsymbol{x}_{0})\}
\,\mathrm{d}S(\boldsymbol{x})}_{Double \ Layer Operator}.
\end{align}$$
Here, 
$\boldsymbol {x}_{0}$ denotes any point lying inside the domain occupied by the carrier fluid and 
$G_{l}(\boldsymbol {x}, \boldsymbol {x}_{0})$ denotes the Green's function for the Laplace equation, which can be expressed as
\begin{equation} G_{l}(\boldsymbol{x},\boldsymbol{x}_{0}) = \frac{1}{4{\rm \pi}(|\boldsymbol{x}-\boldsymbol{x}_{0}|)}, \end{equation}
where 
$|\boldsymbol {x}-\boldsymbol {x}_{0}|$ is the magnitude of the Euclidean distance between the points 
$\boldsymbol {x}$ and 
$\boldsymbol {x}_{0}$. Here, we introduce two integral operators, the single- and double-layer operator. When the integrand constitutes a singular Green's function multiplied by a continuous function, the integral operator is referred to as the single-layer operator. However, a double-layer operator involves a continuous function multiplied by the normal component of the gradient of the Green's function, as can be seen in (3.1). This notation will be used to refer to these integrals in all subsequent discussions.
It can be observed that the surface integrals in (3.1) are carried out over the surface of the droplet and the rigid walls. This can be simplified by carefully modifying the Green's function such that it vanishes on the walls (Pozrikidis et al. Reference Pozrikidis1992):
\begin{equation} G_{l}^{2W} = 0 \quad {\rm and} \quad \boldsymbol{\nabla} G_{l}^{2W} = 0 \quad \forall \boldsymbol{x} \in B, \end{equation}
where 
$G_{l}^{2W}$ represents the modified Green's function. Details on the computation of 
$G_{l}^{2W}$ is given in Appendix A.1. In (3.1), we neglect the contribution of the surface integrals from the surfaces bounding the domain along the 
$x$ and 
$y$ directions, with the assumption that they lie far away from the droplet interface. Using the modified Green's function, (3.1) can be simplified as
$$\begin{align} T_{e}(\boldsymbol{x}_{0}) &= T_{\infty}(\boldsymbol{x}_{0}) - \int_{\varGamma}G_{l}^{2W}(\boldsymbol{x},\boldsymbol{x}_{0})\{ \boldsymbol{\nabla} T_{e}(\boldsymbol{x})\boldsymbol{\cdot}\hat{\boldsymbol{n}}(\boldsymbol{x}) \} \,\mathrm{d}S(\boldsymbol{x}) \nonumber\\ &\quad + \int_{\varGamma}T_{e}(\boldsymbol{x})\{ \hat{\boldsymbol{n}}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{\nabla} G_{l}^{2W}(\boldsymbol{x},\boldsymbol{x}_{0}) \}\,\mathrm{d}S(\boldsymbol{x}). \end{align}$$
We consider that the point 
$\boldsymbol {x}_{0}$ is being brought vanishingly close to the surface of the droplet from the outer side, and consider the principal value of the surface integrals:
$$\begin{align} \frac{T_{s}(\boldsymbol{x}_{0})}{2} &= T_{\infty}(\boldsymbol{x}_{0}) - \int_{\varGamma}G_{l}^{2W}(\boldsymbol{x},\boldsymbol{x}_{0})\{ \hat{\boldsymbol{n}}(\boldsymbol{x})\boldsymbol{\cdot} \boldsymbol{\nabla} T_{e}(\boldsymbol{x}) \} \,\mathrm{d}S(\boldsymbol{x}) \nonumber\\ &\quad + \int_{\varGamma}^{PV} T_{s}(\boldsymbol{x})\{ \hat{\boldsymbol{n}}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{\nabla} G_{l}^{2W}(\boldsymbol{x},\boldsymbol{x}_{0}) \} \,\mathrm{d}S(\boldsymbol{x}), \end{align}$$
where 
$T_{s}$ denotes the temperature field on the interface of the droplet and the superscript 
$PV$ represents the principal value of the double-layer operator.
We next consider the boundary integral representation for the temperature field inside the droplet. The temperature inside the droplet can be expressed as
$$\begin{align} T_{i}(\boldsymbol{x}_{0}) &= \int_{\varGamma}G_{l}^{2W}(\boldsymbol{x},\boldsymbol{x}_{0})\{ \hat{\boldsymbol{n}}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{\nabla} T_{i}(\boldsymbol{x}_{0})\}\,\mathrm{d}S(\boldsymbol{x}) \nonumber\\ &\quad - \int_{\varGamma}T_{i}(\boldsymbol{x})\{ \hat{\boldsymbol{n}}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{\nabla} G_{l}^{2W}(\boldsymbol{x},\boldsymbol{x}_{0}) \} \,\mathrm{d}S(\boldsymbol{x}), \end{align}$$
where 
$\boldsymbol {x}_{0}$ denotes any point inside the droplet. On letting the point 
$\boldsymbol {x}_{0}$ approach the droplet interface, we obtain
$$\begin{align} \frac{T_{s}(\boldsymbol{x}_{0})}{2} &= \int_{\varGamma}G_{l}^{2W}(\boldsymbol{x},\boldsymbol{x}_{0})\{ \hat{\boldsymbol{n}}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{\nabla} T_{i}(\boldsymbol{x}_{0})\}\,\mathrm{d}S(\boldsymbol{x}) \nonumber\\ &\quad - \int_{\varGamma}^{PV}T_{s}(\boldsymbol{x})\{ \hat{\boldsymbol{n}}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{\nabla} G_{l}^{2W}(\boldsymbol{x},\boldsymbol{x}_{0}) \} \,\mathrm{d}S(\boldsymbol{x}). \end{align}$$ We eliminate the normal derivatives of the temperature gradient by multiplying (3.7) with 
$\delta$ and adding to (3.5). The normal derivative terms cancel out due to the continuity of the heat flux across the interface (2.4). This leads to the final form of the boundary integral equation used for computing the interfacial temperature:
\begin{equation} \frac{1+\delta}{2}T_{s}(\boldsymbol{x}_{0}) = T_{\infty}(\boldsymbol{x}_{0}) + (1 - \delta)\int_{\varGamma}^{PV}T_{s}(\boldsymbol{x})\{ \hat{\boldsymbol{n}}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{\nabla} G_{l}^{2W}(\boldsymbol{x},\boldsymbol{x}_{0})\}\,\mathrm{d}S(\boldsymbol{x}), \end{equation}
where 
$\boldsymbol {x}_{0}$ denotes the position vector of any point on the droplet interface.
The interfacial velocity of an equi-viscous droplet can be obtained using the well-established equations from the boundary integral representation of Stokes’ flow in the presence of deformable droplets (Pozrikidis et al. Reference Pozrikidis1992; Janssen & Anderson Reference Janssen and Anderson2008):
\begin{equation} \boldsymbol{u}_{s}(\boldsymbol{x}_{0}) = \boldsymbol{u}_{\infty} - \frac{1}{8{\rm \pi}}\int_{\varGamma}[(\boldsymbol{\mathsf{S}}_{e} - \boldsymbol{\mathsf{S}}_{i})\boldsymbol{\cdot}\hat{\boldsymbol{n}}(\boldsymbol{x})]\boldsymbol{\cdot}\boldsymbol{\mathsf{G}}_{s}^{2W}(\boldsymbol{x},\boldsymbol{x}_{0})\,\mathrm{d}S(\boldsymbol{x}), \end{equation}
where 
$\boldsymbol{\mathsf{G}}_{s}^{2W}$ refers to the Green's function for the Stokes equation modified to satisfy the no-slip condition at the walls (Griggs et al. Reference Griggs, Zinchenko and Davis2007). The computation of the modified Green's function for the Stokes equation is explained in Appendix A.2. The traction jump appearing in (3.9) can be obtained by simplifying the boundary condition for force balance at the interface (2.8):
\begin{equation} (\boldsymbol{\mathsf{S}}_{e}-\boldsymbol{\mathsf{S}}_{i})\boldsymbol{\cdot}\hat{\boldsymbol{n}}(\boldsymbol{x}) = \frac{1}{Ca}[\sigma(\boldsymbol{x})\kappa_{m}(\boldsymbol{x})\hat{\boldsymbol{n}}(\boldsymbol{x}) - \boldsymbol{\nabla}_{s}\sigma(\boldsymbol{x}) ], \end{equation}
where 
$Ca$ denotes the capillary number.
3.2. Discretization in space and time
3.2.1. Mesh generation
The boundary integral equations for the interfacial velocity (3.9) and temperature (3.8) are solved on a dynamically evolving interface which is triangulated into an unstructured grid. This grid is generated by successively dividing an icosahedron into triangles (Yon & Pozrikidis Reference Yon and Pozrikidis1998). At each step, every triangle is subdivided into four new triangles, such that the total number of triangles increases by four times in each iteration. Finally, the points are projected onto a unit sphere. The programs for performing this procedure were adopted from the open-source package BEMLIB (Pozrikidis Reference Pozrikidis2002). Henceforth, the curved triangles will be referred to as elements, which are uniquely identified using a set of six points, three on its vertices and one on the midpoint of each edge. In this study, we found that 1280 elements were sufficient to obtain grid independent results.
3.2.2. Evaluation of boundary integrals
The surface integrals were evaluated by summing the integrals up over the individual elements. To perform the numerical integrations, a curved element is mapped to a planar right-angle triangle in an orthogonal parametric space and all the relevant quantities are mapped using a second-order isoparametric mapping (Yon & Pozrikidis Reference Yon and Pozrikidis1998). When the point of singularity 
$\boldsymbol {x}_{0}$ lies on an element over which the integral is evaluated, that element is referred to as a singular element. For a singular element, the single-layer operator is evaluated using a polar integration rule (Pozrikidis et al. Reference Pozrikidis1998) and the double-layer operator is evaluated by subtracting the singular part and integrating it analytically (Pozrikidis et al. Reference Pozrikidis1992). Over non-singular elements, the single-layer and double-layer operators are evaluated using a 7-point Gauss–Legendre quadrature.
3.2.3. Numerical computation of traction discontinuity
In this section, we discuss the methods for computation of the right-hand side of (3.10). The important quantities of interest are: (i) unit normal vector and curvature; (ii) surface tension; and (iii) surface gradient of interfacial tension. The computation of the mean curvature is the most sensitive part of the simulation, as inaccuracies in the curvature calculation give rise to numerical instabilities, especially when large deformations are encountered. To mitigate this difficulty, the normal vector and curvature are computed using a local parabolic interpolation of the surface at each point (Zinchenko, Rother & Davis Reference Zinchenko, Rother and Davis1997). Once the interfacial temperature is known at every point, the surface tension is computed using (2.1). The surface gradient of the interfacial tension is obtained by an indirect method which involves the solution of a 
$3\times 3$ matrix equation relating the surface gradient in the Cartesian co-ordinates with the surface gradient in the local parametric space of an element (Yon & Pozrikidis Reference Yon and Pozrikidis1998).
3.2.4. Advancing in time
Once the interfacial velocity is computed at every point, the interface is advanced using a second-order accurate Runge–Kutta method (RK2), which can be expressed as follows:
\begin{equation} \boldsymbol{x}_{n+1} = \boldsymbol{x}_{n} + \frac{\Delta t}{2}\{ \boldsymbol{V}_{n} + \boldsymbol{V}_{n+1/2} \}, \end{equation}
where 
$\boldsymbol {x}_{n+1}$ and 
$\boldsymbol {x}_{n}$ represent the position of the interface at the time step 
$n$ and 
$n+1$, respectively, 
$\Delta t$ denotes the size of the time step, 
$\boldsymbol {V}_{n}$ is the velocity of the interfacial points at the 
$n$th time step and 
$\boldsymbol {V}_{n+1/2}$ refers to the velocity of the interfacial points at an intermediate time step where the position of the interface is approximated using the velocity (
$\boldsymbol {V}_{n}$) and position (
$\boldsymbol {x}_{n}$) at the 
$n$th time step. The velocities of the interfacial marker points (
$\boldsymbol {V}_{n},\boldsymbol {V}_{n+1/2}$) are computed by a linear superposition of the normal component of the interfacial velocity obtained from BEM and an artificial tangential velocity referred to as the mesh-relaxation (MR) velocity. The MR velocity is imposed to maintain the stability of the mesh for long simulation times. The calculation of this velocity is based on an iterative technique for minimizing the mesh distortion energy (Zinchenko et al. Reference Zinchenko, Rother and Davis1997). It should be noted that the addition of this artificial tangential velocity to the interfacial marker points does not alter the physics of the problem as it does not lead to any additional deformation of the interface. However, the inclusion of the MR is extremely important because, without the same, the simulation convergence in a background pressure-driven flow is not achieved for all scenarios.
3.3. Summary of the numerical method
The overall numerical approach may be summarized as follows.
(i) Initialization and mesh generation.
(ii) Computation of geometric quantities like normal vector and curvature.
(iii) Solution of interfacial temperature from BEM.
(iv) Computation of surface tension and interfacial gradient of surface tension.
(v) Computation of interfacial velocity from BEM.
(vi) Calculation of MR velocity.
(vii) Advection of interfacial marker points using RK2 method (3.11).
(viii) Go to step (ii) and iterate for subsequent time steps.
Further details regarding the computational implementation of the numerical algorithm are outlined in Appendix B.2.
4. Results and discussions
For illustration, we consider the migration of a neutrally buoyant deformable droplet having the same viscosity as the carrier fluid in a linearly varying temperature field. We consider two separate cases: (i) linearly increasing temperature profile (
$\zeta$ = 1); and (ii) linearly decreasing temperature profile (
$\zeta = -1$) along the flow direction. The magnitude of the temperature gradient is quantified using the non-dimensional Marangoni number (
$Ma_{T}$). It has been noticed that the physically acceptable range for 
$Ma_{T}$, within the limit of droplet breakup, is different for 
$\zeta = 1$ and 
$\zeta = -1$. For 
$\zeta = 1$, increasing the value of 
$Ma_{T}$ beyond 
$0.05$ leads to excessive droplet deformation, and a further increase in 
$Ma_{T}$ leads to breakup. This phenomenon can be understood qualitatively by considering the fact that the surface tension of an interface decreases with an increase in temperature. Hence, as the droplet migrates towards a region of higher temperature, the average temperature on the droplet interface also increases. Due to this, the average surface tension on the droplet interface decreases, leading to larger deformations. Quantitatively, this effect is captured by an instantaneous capillary number (
$Ca_{i} = {\mu V_{c}}/{\sigma _{avg}}$) which varies as the position of the droplet evolves in time. Here, 
$\sigma _{avg}$ refers to the average value of the surface tension at the interface at any given point in time. For 
$\zeta = 1$, 
$Ca_{i}$ increases with time. However, 
$Ca_{i}$ decreases with time for 
$\zeta = -1$. For 
$\zeta = -1$, we have found that the maximum value of 
$Ma_{T} \sim 0.10$. We have arrived at this estimate by taking experimentally obtained values of 
$\beta$ (Yakhshi-Tafti, Kumar & Cho Reference Yakhshi-Tafti, Kumar and Cho2011) and the maximum value of 
$\Delta T_{c}$ that obviates any perceptible role of buoyancy. The distinguishable interfacial phenomenon, despite holding the same magnitude of the imposed temperature gradient, occurs in cases of the positive and the negative temperature gradients because of its interconnection with the shape deformation, which had not been considered in the reported models. However, in the subsequent discussions, we show that this effect can have a strong influence on the physical outcome, so much so that it contradicts the reported theoretical predictions. The difference in magnitude of Marangoni number for the two directions is due to the physical asymmetries associated with positive and negative temperature gradients. In the case of a positive temperature gradient, the average capillary number increases, which leads to greater droplet deformation with time; whereas for a negative temperature gradient, the average surface tension increases and the deformation of the droplet does not increase with time. To avoid situations of perpetual enhancements in the droplet deformation to an extent that the same approach the conditions for droplet breakup, which does not conform to the focal theme of this work, we have limited the Marangoni number to the case of a positive temperature gradient. Since alteration of the directions of the temperature gradient by itself disrupts the physical symmetry, adhering to the same order of magnitude of it for the two scenarios, thus, does not necessarily conform to the same physics. Consistent with the focal theme of the present work that aims to establish significant differences between the droplet dynamics, as predicted by asymptotic theory and the same obtained via numerical computation, we therefore purposefully attempted to remain restrictive to the parametric limits that ensure droplet migration without approaching physical disintegration for which the qualitative physics by itself would change altogether and the essential basis of comparing with the asymptotic theory loses its fundamental proposition. The key to this proposition has been to highlight the role of a dynamic two-way coupling between the thermal and the flow fields at the interface, and exemplify serious deficits in the reported theoretical propositions that superimpose directly the external thermal field instead.
We begin by validating our numerical method with the existing literature in § 4.1, where we validate both the trajectory and shape of a deformable droplet migrating in an isothermal pressure-driven flow. Subsequently, in § 4.2, we discuss the effect of the temperature gradient by comparing the trajectory and axial migration velocity of a droplet for a positive temperature gradient, negative temperature gradient and isothermal flow. These results are discussed in detail using physical arguments supported by additional numerical experiments. In § 4.3, we provide a direct comparison of the results obtained from the present numerical simulations with the results obtained from reported analytical theories, for a valid range of physical parameters obeying the limits of the assumptions made in analytical theories. This shows that the conclusions drawn from the present simulations hold true in a broad sense, even within the parametric limits of the analytical theories. Finally, we investigate the quantitative variations in migration trajectory and the axial component of migration velocity, due to variations in the confinement ratio 
$Wc$ and magnitude of the thermal Marangoni number 
$Ma_{T}$ in §§ 4.4 and 4.5, respectively.
4.1. Validation
Figure 2 depicts a comparison of the droplet migration trajectory obtained from the present method with the results obtained by Griggs et al. (Reference Griggs, Zinchenko and Davis2007). In this problem, the evolution of an initially spherical droplet in an isothermal pressure driven flow between two parallel plates is investigated. Initially, we place the droplet at an off-centreline position (
$Z_{c} = -0.10$) and follow its motion. We observe that the trajectories obtained from the two separate numerical methods show nearly identical characteristics. Only a slight deviation may be attributed to the difference in the methods of numerical integration employed for computing the single-layer and double-layer integrals. The integration quadrature used in the present study has been found to yield higher order accuracy and also requires fewer discrete surface elements for obtaining grid independent results. In the insets, we have shown the shape of the droplet at various axial locations. These shapes also match well with those presented by Griggs et al. (Reference Griggs, Zinchenko and Davis2007).
Figure 2. Variation in the transverse location of the droplet centroid with the axial position of the centroid for an isothermal flow. The surface plots show the deformed shape of the droplet at various locations, as demarcated by the respective arrows. The values of relevant physical parameters are 
$Ca = 0.33$, 
$W_{c} = 0.60$.
The results presented in figure 2 serve as a baseline for studying the thermocapillary dynamics of a droplet, which is discussed in subsequent sections. Therefore, it is imperative to identify certain key features from this figure. We notice that the droplet deforms maximally near the beginning of its trajectory, where it has a highly asymmetrical shape. Subsequently, the shape evolves to attain a greater degree of symmetry with respect to the axis of the imposed flow. When the droplet reaches its equilibrium position at the centreline of the flow, it attains a ‘bullet-like’ shape where the rear section of the droplet has a higher radius of curvature than the front section.
In the investigations presented in the subsequent sections, our approach remains the same, i.e. we initially place a spherical droplet at an off-centreline position and observe its evolution in time. For all the cases reported in this article, we fix the value of the conductivity ratio 
$\delta$ to 0.10 as we observed that the change in 
$\delta$ does not lead to any significant differences in the migration dynamics of the droplet within the physically realistic limits of 
$\delta$.
Before proceeding to the analysis for a non-isothermal flow, we illustrate a comparison of the numerical results with analytical theory for the limiting case of a droplet suspended in an unbounded circular Poiseuille flow, under the assumption of negligible droplet deformation. Mathematically, this translates to a very low value of the capillary number and confinement ratio i.e. 
$Ca \ll 1$ and 
$W_c \ll 1$. To quantify the radial location of the droplet centroid, we introduce the variable 
$e = \sqrt {Z_c^2 + Y_c^2}$, denote the radius of the circular channel by 
$R$ and represent the radial migration velocity of the droplet by 
$U_r$ such that a positive value of 
$U_r$ denotes migration away from the channel centreline. Interestingly, analytical theory for an equiviscous droplet in a unbounded Poiseuille flow predicts a migration velocity of the droplet which is directed away from the channel centreline (Haber & Hetsroni Reference Haber and Hetsroni1971). In figure 3, we show that the results from our numerical method exhibit a good match with the analytical predictions. The deviation between the two approaches becomes significant only for higher values of 
$e/R$. This can be attributed to an increase in the magnitude of fluid shear with an increase in the distance of the droplet centroid from the channel centreline, leading to a higher degree of droplet deformation which violates the assumption of negligible droplet deformation made in analytical theory.
Figure 3. Transverse migration velocity of a droplet as a function of lateral position of the centroid when the droplet is suspended in a circular Poiseuille flow. The values of relevant physical parameters are 
$Ca = 0.01$, 
$W_{c} = 0.10$.
4.2. Effect of temperature gradient
The migration trajectory of a droplet and the axial component of centroid velocity (
$U_{x}$) are plotted in figures 4(a) and 4(b), respectively, for three cases: 
$\zeta = 1$, 
$\zeta = 0$ (isothermal flow) and 
$\zeta = -1$. From figure 4(a), we can observe that in the presence of a linearly decreasing temperature field (i.e. negative temperature gradient along the channel length), the droplet migrates faster towards the centreline as compared to an isothermal flow and in the presence of a linearly increasing temperature field (i.e. positive temperature gradient along the channel length), the cross-streamline migration of the droplet is slowed down. Similarly, the axial component of the velocity also increases in the presence of a linearly decreasing temperature profile and decreases in the presence of a linearly increasing temperature field, as compared to an isothermal flow. These results contradict the prediction of recently reported analytical theories which suggest that a linearly increasing temperature field will speed up both the axial and cross-stream migration of the droplet due to the Marangoni effect (Das et al. Reference Das, Mandal, Som and Chakraborty2017, Reference Das, Mandal and Chakraborty2018). This contradiction is rooted in the fact that the reported analytical theories assume a spherical interface and solve for the temperature profile at a given transverse location. This temperature profile is used to calculate the migration velocity of the droplet at that transverse location using Faxen's laws (Kim & Karrila Reference Kim and Karrila2013). However, in reality, when the droplet migrates in a non-isothermal environment, the average local temperature over the surface of the droplet also changes. This leads to a variation in the average value of the interfacial tension (
$\sigma _{avg}$). When the droplet migrates along a linearly decreasing temperature field (
$\zeta = -1$), 
$\sigma _{avg}$ increases with time and analogously, when 
$\zeta = +1$, 
$\sigma _{avg}$ decreases with time. Here we show that the variation in 
$\sigma _{avg}$ causes a change in the magnitude of the component of the viscous drag force on the droplet (
$\boldsymbol {F}_{D}$) which is exclusively due to the effect of interfacial tension, and this change is an order of magnitude larger than the force due to the Marangoni effect. It should be pointed out that 
$\boldsymbol {F}_{D}$ is not the total drag force on the droplet, but it is that component of the drag force which arises exclusively due to the effect of interfacial tension of the droplet interface. It has been shown that an equiviscous droplet migrates towards the centreline in a confined Poiseuille flow due to a confluence of wall-induced lift forces and deformability of the droplet interface (Griggs et al. Reference Griggs, Zinchenko and Davis2007; Janssen & Anderson Reference Janssen and Anderson2008). Since the viscous drag opposes this motion, it implies that a reduction in magnitude of the drag forces facilitates a faster cross-streamwise migration of the droplet.
Figure 4. (a) Variation of transverse centroid position with axial centroid location, and (b) variation of axial centroid velocity with time, when 
$Ca = 0.10$, 
$\delta = 0.10$ and 
$W_{c} = 0.50$.
To estimate the magnitude of 
$\boldsymbol {F}_{D}$, we refer to the work by Haber & Hetsroni (Reference Haber and Hetsroni1971), where they considered the migration of a deformable droplet in a Poiseuille flow in the limit of vanishingly small values of the capillary number (
$Ca$). In equation (29) of their paper, they present an expression for the transverse component of the drag force (
$F_{Dz}$) acting away from the centreline. On analysing this expression, we find that if fluid viscosities, channel width and centreline velocity of the flow is kept constant, 
$F_{Dz} \sim {a^{4}}/{\sigma _{avg}}$, where 
$a$ denotes the droplet radius. As we are not considering very large deformations of the droplet interface, we can approximate 
$a^2 \sim A$, where 
$A$ denotes the total surface area of the droplet and arrive at the following scaling relation:
\begin{equation} F_{Dz} \sim \frac{A^2}{\sigma_{avg}}. \end{equation} From (4.1), it can be inferred that the magnitude of 
$F_{Dz}$ varies inversely with the average interfacial tension 
$(\sigma _{avg})$, which changes significantly during the course of droplet migration. This change is proportional to the magnitude of the temperature gradient multiplied by the length scale of axial migration of the droplet, i.e. 
$\Delta \sigma _{avg} \sim |G|\Delta X_{c}$. In addition to the viscous drag force, there is another force due to the Marangoni effect which acts along a direction opposite to the viscous drag force. In the present case, Marangoni force results from a gradient in interfacial tension due to temperature difference across the droplet interface. Hence, the magnitude of Marangoni force scales as the temperature gradient multiplied by the length scale of the droplet, i.e. 
$F_{Ma} \sim |G|a$, where 
$F_{Ma}$ is representative of the magnitude of the Marangoni force. It can be seen from figures 2 and 4 that 
$\Delta X_{c} \gg a$, which implies that the magnitude of 
$\Delta \sigma _{avg}$ is greater than 
$F_{Ma}$. It can be seen from (2.1) that for a droplet migrating along a linearly decreasing temperature gradient, 
$\sigma _{avg}$ increases during the course of migration, which contributes towards lowering the value of 
$F_{Dz}$ and analogously, 
$F_{Dz}$ increases when the droplet migrates along an increasing temperature gradient. In addition to 
$\sigma _{avg}$, 
$F_{Dz}$ also depends on the total surface area of the droplet. In figure 5, we depict the evolution of the total surface area of the droplet with time and it can be seen that the total surface area is smaller for a negative temperature gradient as compared with the isothermal case and is larger for the case of a positive temperature gradient. Since 
$F_{Dz} \sim A^2$, a decrease in 
$A$ will also contribute towards a decrease in the magnitude of the drag force.
Figure 5. Variation of total surface area of the droplet with time when 
$Ca = 0.10$, 
$\delta = 0.10$, 
$W_{c} = 0.50$, with 
$Ma_{T} = 0.20$ for a negative temperature gradient and 
$Ma_{T} = 0.0125$ for a positive temperature gradient.
Therefore, the combined effects of change in 
$\sigma _{avg}$ and 
$A$ leads to a decrease in the viscous drag force acting on the droplet when it migrates along a linearly decreasing temperature profile. This has been confirmed by plotting the temporal evolution of 
${A^2}/{\sigma _{avg}}$ in figure 6, where it can be seen that the value of 
${A^2}/{\sigma _{avg}}$ decreases with time for a negative temperature gradient and increases for the case of a positive temperature gradient. Moreover, it has been shown using scaling arguments that the change in magnitude of viscous drag force is much greater than the Marangoni force. Although Haber & Hetsroni (Reference Haber and Hetsroni1971) did not obtain any correction to the axial component of the drag force due to interfacial tension, due to their assumption of vanishingly small droplet deformation, it can be expected that the axial component of the viscous drag also follows a similar scaling law as the underlying mechanism behind both the axial and cross-streamwise components of drag force is similar. This explains the change in magnitude of both the axial and cross-streamwise components of the droplet migration velocity in the presence of a non-isothermal temperature field.
Figure 6. Variation of 
${A^2}/{\sigma _{avg}}$ with time when 
$Ca = 0.10$, 
$\delta = 0.10$, 
$W_{c} = 0.50$, with 
$Ma_{T} = 0.20$ for a negative temperature gradient and 
$Ma_{T} = 0.0125$ for positive temperature gradient.
In figures 7 and 8, we show the interfacial variation in temperature and surface tension, respectively, for the results corresponding to figure 4. For constructing the surface plots, we compute the average value of a variable over the surface and plot the deviation at each point. The deviations increase from a negative to positive value as the colour changes progressively from blue to yellow (visualized by an equivalent change in shade in greyscale). The dark blue regions indicate a lower value of the variable as compared to the average and the bright yellow regions indicate an above average value for the variable. The plots in figure 7(a i–a iv) denote the temporal evolution of the interfacial temperature for a droplet migrating in a negative temperature gradient. Similarly, the plots in figure 7(b i–b iv) denote the temporal evolution of the interfacial temperature for a droplet migrating in a positive temperature gradient. We observe that in figure 7(b i–b iv), the front section of the droplet has a higher temperature than the rear section, as it is migrating along a path of increasing temperature. Similarly, in figure 7(a i–a iv), we see that the front section has a lower value of the temperature as the droplet is moving along a path of decreasing temperature. By comparing figure 7(a i–a iv) with figure 7(b i–b iv), one may notice that the shapes for figure 7(b i–b iv) show a greater degree of deformation as compared to those for figure 7(a i–a iv). This difference is most prominent between figures 7(a iv) and 7(b iv), even though the value of 
$Ma_{T}$ corresponding to figure 7(b iv) is less than one tenth as compared to that for figure 7(a iv). This again highlights the physical asymmetry of the scenarios between a spatially increasing and decreasing temperature gradient of the same magnitude. In figure 8(a i–a iv), we observe that the front section of the droplet has a higher value of surface tension when the same migrates along a linearly decreasing temperature field. This is expected, as the surface tension of the interface increases with a decrease in the interfacial temperature. Following a similar point of view, we can explain the temporal evolution of the surface tension for the linearly increasing temperature field, depicted in figure 8(b i–b iv).
Figure 7. Variation of interfacial temperature over the droplet surface for (a i–a iv) the negative temperature gradient, and (b i–b iv) the positive temperature gradient when 
$Ca = 0.10$, 
$\delta = 0.10$, 
$Wc = 0.50$, with 
$Ma_{T} = 0.20$ for linearly decreasing and 
$Ma_{T} = 0.0125$ for linearly increasing temperature profiles.
Figure 8. Variation of interfacial tension over the droplet surface for (a i–a iv) the negative temperature gradient and (b i–b iv) the positive temperature gradient, when 
$Ca = 0.10$, 
$\delta = 0.10$, 
$Wc = 0.50$, with 
$Ma_{T} = 0.20$ for decreasing and 
$Ma_{T} = 0.0125$ for the linearly increasing temperature profile.
To summarize, we have investigated the axial and cross-streamline migration of a droplet in a non-isothermal flow and have shown that, in addition to the interfacial gradients in surface tension, one also needs to account for a variation in the average value of the surface tension. In fact, this variation in the average value of the surface tension dictates the nature of migration trajectory of the droplet. This has been explained in detail by considering the variation of the total viscous drag on the droplet and it has been shown that our arguments are in line with the classical analytical results for a deformable viscous droplet.
4.3. Direct comparison with analytical theory
In this section, we present a direct one-to-one comparison between the results obtained from the present study and the analytical study by Das & Chakraborty (Reference Das and Chakraborty2018), where they computed the migration velocity of a droplet by considering small deviations of the shape from a sphere and a one-way coupling between the thermal field and fluid flow.
In figure 9, we plot the axial migration velocity for a droplet located initially at an off-centreline position, and compare the results obtained from the present numerical simulations and reported analytical theory. We observe that there is a clear deviation regarding the effect of the temperature gradient. In figure 9(a), we see that the present results predict a decrease in magnitude of axial migration velocity on application of a positive temperature gradient along the direction of the flow and an increase in the migration velocity on application of a negative temperature gradient. However, it can be seen from figure 9(b) that the analytical results predict an increase in the migration velocity on application of a positive temperature gradient and a decrease in the axial migration velocity upon application of a negative temperature gradient.
Figure 9. Comparison of the axial component of migration velocity with transverse centroid location from (a) numerical simulations and (b) analytical calculations when 
$Wc = 0.10$, 
$Ca = 0.10$, 
$\delta = 0.10$, with 
$Ma_{T} = 0.20$ for 
$\zeta = -1$ and 
$Ma_{T} = 0.0125$ for 
$\zeta = 1$.
In figure 10, the cross-stream migration velocity of the droplet is plotted and here also we see a clear contradiction between the results obtained from the existing analytical theory and the present simulations. The analytical theories predict an increase in the transverse migration velocity on application of a positive temperature gradient, whereas the completely opposite outcome is predicted by the numerical simulations. Similarly, on application of a linearly decreasing temperature gradient, the transverse migration velocity is supposed to decrease according to analytical theory but the present simulations show that it actually increases.
Figure 10. Comparison of transverse component of migration velocity with transverse centroid location from (a) numerical simulations and (b) analytical calculations when 
$Wc = 0.10$, 
$Ca = 0.10$, 
$\delta = 0.10$, 
$Ma_{T} = 0.20$ for 
$\zeta = -1$ and 
$Ma_{T} = 0.0125$ for 
$\zeta = 1$.
4.4. Effect of confinement ratio
In figure 11, we have compared the trajectory and axial migration velocity of a droplet for various values of confinement ratio (
$Wc$). In figure 11(a i), we have compared the migration trajectories for a droplet suspended in a linearly decreasing temperature field and the corresponding result for a linearly increasing temperature field is depicted in figure 11(b i). From figures 11(a i) and 11(b i), we observe that the extent of centreline migration increases as the value of the confinement ratio is increased. This effect can be primarily attributed to the fact that a larger droplet experiences a greater magnitude of lift forces from the rigid walls which tend to push it towards the centreline. If we closely compare figures 11(a i) and 11(b i), we can note that the difference between the migration trajectories for 
$Wc = 0.50$ and 
$Wc = 0.60$ is more prominent when the droplet migrates along a linearly increasing temperature field. This is because of the fact that since a larger droplet extends across a longer distance, the Marangoni effects due to the variation in surface tension become more prominent. In the case of a positive temperature gradient, it aids in increasing the magnitude of cross-streamline migration, whereas in the presence of a negative temperature gradient, it slows down the migration speed. Hence, although the Marangoni effect does not dictate the overall dynamics of a droplet, it explains the differences in the migration trajectory with confinement ratio between a linearly increasing and decreasing temperature field.
Figure 11. Variation of migration trajectory (a i,b i) and axial migration velocity (a ii,b ii) with confinement ratio for a droplet suspended in a linearly decreasing (a i,a ii) and linearly increasing (b i,b ii) temperature field, when 
$Ca = 0.10$, 
$\delta = 0.10$ and 
$Ma_{T} = 0.0125$.
In figures 11(a ii) and 11(b ii), we compare the axial migration velocity for different values of confinement ratio for 
$\zeta = -1$ and 
$\zeta = 1$, respectively. In both cases, we observe that the magnitude of the axial migration velocity decreases with an increase in the confinement ratio. This can be explained by noting that a smaller droplet experiences a lower magnitude of viscous drag from the surrounding flow, due to a lower value of its cross-sectional area. In this case also, we observe that the Marangoni effect leads to a higher order correction to the axial migration velocity. This can be seen by closely examining the differences in the axial migration velocity between 
$Wc = 0.50$ and 
$Wc = 0.60$ for figures 11(a ii) and 11(b ii). Here, we can observe that the difference between the values of 
$U_{x}$ for 
$Wc = 0.60$ and 
$Wc = 0.50$ is less in figure 11(a ii) than in figure 11(b ii). This is because for 
$\zeta = -1$, the Marangoni effect tends to slow down the axial migration velocity of the droplet, and this effect gains more prominence when the size of the droplet increases. Since the axial velocity of a droplet decreases with an increase in 
$Wc$ due to viscous drag, the Marangoni effect further slows down the droplet leading to a greater difference in the migration speeds on increasing 
$Wc$. However, when 
$\zeta = 1$, the Marangoni effect increases the axial migration velocity, leading to a lesser magnitude of decrease in the axial migration speed on increasing 
$Wc$.
Therefore, we observe that although the magnitude of change in velocity due to the Marangoni effect is small in comparison to that due to viscous drag, it explains the differences in the variation of droplet migration with confinement ratio for an increasing and decreasing temperature profile.
4.5. Effect of Marangoni number
In this section, we quantify the difference in thermo-capillary dynamics of a deformable droplet when subjected to alterations in the thermal Marangoni number, 
$Ma_{T}$. In figure 12(a), we show the migration trajectories for two different values of 
$Ma_{T}$ in a linearly decreasing temperature profile. Here, we observe that as the value of 
$Ma_{T}$ increases, the droplet migrates faster towards the centreline. This effect can be explained by considering the fact that in the case of 
$\zeta = -1$, the increase in the magnitude of the average interfacial tension is greater for a larger value of 
$Ma_{T}$, leading to a greater reduction in the viscous drag. To confirm this, we have plotted 
${A^2}/{\sigma _{avg}}$ for 
$Ma_{T} = 0.20$ and 
$Ma_{T} = 0.0125$ in figure 13. As we have shown in an earlier section that the drag force scales with 
${A^2}/{\sigma _{avg}}$, it implies from figure 13 that the decrease in magnitude of viscous drag increases as the value of 
$Ma_{T}$ is increased. Due to a reduction in the viscous drag on increasing the value of 
$Ma_{T}$, the migration speed of the droplet increases. This is also reflected in figure 12(b), where we observe an increase in the magnitude of the axial migration velocity due to an increase in 
$Ma_{T}$ when 
$\zeta = -1$. However, we observe that the change in magnitude of the axial migration velocity is not very large as compared to the change observed by varying the confinement ratio (
$Wc$) or by switching the direction of the temperature gradient. This can be explained by considering the fact that the magnitude of decrease in viscous drag, due to the increase in 
$Ma_{T}$, is very small as compared to the variation in drag forces due to changes in 
$Wc$ and 
$\zeta$.
Figure 12. Variation of migration trajectory (a) and axial migration velocity (b) with 
$Ma_{T}$, for a droplet suspended in a linearly decreasing temperature field when 
$Ca = 0.10$, 
$\delta = 0.10$ and 
$Wc = 0.50$.
Figure 13. Variation of 
${A^2}/{\sigma _{avg}}$ with time for two different values of 
$Ma_T$ when 
$\zeta = -1$, 
$Ca = 0.10$, 
$\delta = 0.10$ and 
$Wc = 0.50$.
5. Conclusions
We have presented a robust three-dimensional computational framework for investigating the thermo-capillary dynamics of a deformable droplet in a confined fluidic environment. This method allows the imposition of any external velocity and temperature field, in a framework that is consistent with the boundary conditions relevant to wall-bounded flows. We have used this framework for analysing the motion of a deformable droplet in pressure-driven flow, in the presence of a temperature field varying linearly along the axis of the flow. The results obtained from this study indicate drastic alteration in the migration dynamics of a droplet as compared to the predictions from reported analytical theories which consider asymptotically small shape deformation of the interface. In particular, we observe that when a droplet migrates towards a colder region, its migration speed increases both along the axial and cross-streamwise directions, and alternatively, when a droplet migrates towards a hotter region, its velocity decreases in magnitude along both the directions. This is in contrast to the commonly portrayed consequences of the Marangoni effect which predict a decrease in the migration speed along both the axial and cross-streamwise directions when a droplet migrates in a linearly decreasing temperature field. The apparently non-intuitive physical phenomenon has been explained by considering the change in viscous drag over the surface of the droplet, due to variation in the average interfacial tension of the droplet–carrier fluid interface. A direct one-to-one comparison has been provided between the results obtained from the present numerical method and reported analytical theories, and clear contradictions have been observed between the two. The reasons for these contradictions have been pinpointed and it has been observed that the predictions obtained from the current study also hold true in a broader sense, which encompasses the limiting cases addressed by the analytical theories. Finally, we have deciphered the effects of confinement ratio and thermal Marangoni number on the migration characteristics of the droplet and we conclude that alterations in the confinement ratio can exert significant influence over the quantitative nature of droplet migration. In comparison, the variations due to change in thermal Marangoni number are smaller in magnitude.
While the analysis presented in this work has shed new insights on the exclusive two-way coupling between the thermal and the flow field, as mediated by interfacial stress balance on the droplet surface, severe droplet deformations hallmarking their coalescence and break-up phenomena are not considered in this work. Those considerations, amidst possible entanglements of the Lagrangian grid points located on the interface, need special re-meshing considerations (Cristini, Bławzdziewicz & Loewenberg Reference Cristini, Bławzdziewicz and Loewenberg2001) and physically refer to unstable dynamical regimes, precluding the controllability of precise migration of an identified droplet without further disintegration or re-integration, which are beyond the purview of the central theme of this work.
Acknowledgements
D.P.P. acknowledges the help of Dr S. Ray in proofreading the manuscript. D.P.P. and P.K. are grateful to Mr S. Santra, Mr T.N. Banuprasad and Dr S. Das for providing valuable feedback related to the development of the mathematical model. S.C. acknowledges the Department of Science and Technology, Government of India, for the Sir J.C. Bose National Fellowship and for their Sponsorship under the National Supercomputing Mission (NSM) Project titled: ‘Highly Parallelized Solver for Multiphase Micro-Hydrodynamics’.
Declaration of interests
The authors report no conflict of interest.
Appendix A
In this section, we discuss the computation of the Green's function for both the thermal energy and the momentum equation in a domain bounded by two parallel plates. The premise behind constructing both the Green's functions is to superimpose the free-space Green's function with additional corrections which negate it on the boundary. The computation of the Green's function for the Stokes equation is complicated by the fact that it is a second-order tensor and additional care has to be taken to ensure that all the components vanish at the wall.
A.1. Computation of Green's function for the energy equation
To construct the Green's function for the energy equation, we consider successive reflections of the point of singularity over the two rigid walls. The first sets of two reflections on each wall are then taken as separate points of singularity, and the resulting Green's function is subtracted from the free space Green's function. This leads to a cancellation of the Green's function on the two walls, due to the free space component of Green's function. However, it gives rise to additional contributions at the wall, for example, the image on the upper wall will give rise to a contribution at the lower wall, and analogously, the image at the lower wall will have non-zero contribution on the upper wall. Therefore, additional reflections of these two points have to be considered on the two walls to nullify this contribution, and the new set of images will in turn have additional contributions on the other wall (other than the wall over which they were originally reflected). This gives rise to an infinite sum, whose convergence can be physically understood from the fact that the contribution due to each successive family of reflections is less in magnitude, since the Green's function scales inversely with distance. The modified Green's function can be expressed as
$$\begin{gather} G_{l}^{2W} = \underbrace{\frac{1}{4{\rm \pi}|\boldsymbol{x} - \boldsymbol{x}_{0}|}}_{free \ space\ contribution} + \underbrace{\sum_{k={-}\infty,k\neq0}^{\infty} \frac{1}{4{\rm \pi}\sqrt{(x - x_{0})^2 + (y - y_{0})^2 + (z - (z - (2Hk + z_{0})))^2}}}_{even \ set \ of \ reflections} \nonumber\\ - \underbrace{\sum_{k={-}\infty}^{\infty} \frac{1}{4{\rm \pi}\sqrt{(x - x_{0})^2 + (y - y_{0})^2 + (z - (z - ((2k-1)H - z_{0})))^2}}}_{odd \ set \ of \ reflections}, \end{gather}$$
where 
$k$ belongs to the set of integers. To implement this numerically, we consider the leading order contributions from the odd and even set of reflections and interpolate the higher order terms from a pre-processed array containing the values of the higher order summation terms at several discrete points belonging to the computational domain. This was found to be a fast and accurate method for the computation of Green's function and its derivatives. To calculate the derivatives, we simply differentiate (A1) and follow the same procedure outlined above.
A.2. Computation of Green's function for the Stokes equation
The Green's function was computed using the formulae given in Janssen & Anderson (Reference Janssen and Anderson2007). To illustrate the numerical computation, we take the example of one component of the Green's function tensor, and the procedure for computing all other components is identical:
\begin{equation} G_{s,zz}^{2W} = \underbrace{G_{s,zz}^{\infty}}_{free \ space \ component} + \underbrace{\int_{0}^{\infty}{\rm J}_{0}(qs)t_{1nn}(q,z,z_{0})dq}_{correction \ due \ to \ walls}, \end{equation}
where 
$\textrm {J}_{0}(qs)$ denotes the Bessels's function, 
$t_{1nn}$ is a function of the 
$z$ coordinate of the point of evaluation and point of singularity, and 
$s$ denotes the distance between the point of singularity and point of evaluation in the 
$XY$ plane. To evaluate the improper integral, 
$t_{1nn}$ is decomposed into two components:
\begin{equation} t_{1nn}(q,z,z_{0}) = \hat{t}_{1nn}(z,z_{0}) + \tilde{t}_{1nn}(q,z,z_{0}), \end{equation}
where 
$\hat {t}_{1nn}(z,z_{0}) = \lim _{q\to \infty }t_{1nn}(q,z,z_{0})$. The integral arising from 
$\hat {t}_{1nn}(z,z_{0})$ is evaluated analytically and the other integral is found to be a smooth function of 
$z,z_{0}$ and 
$s^{2}$; therefore, it is pre-computed on a set of discrete points identified by these three variables and interpolated at any given point using a tri-linear interpolation method. We found this method to be simple and accurate for obtaining the Green's function for both the Stokes flow and energy equation.
Appendix B
B.1. Justification of the creeping flow assumption
Under the creeping flow assumption, fluid flow is modelled using the Stokes equation which disregards the inertial terms present in the Navier–Stokes equations. Mathematically, this corresponds to 
$Re \ll 1$ and 
$Re/S \ll 1$, and as a result of this, time does not explicitly appear in the equations of fluid dynamics. This corresponds to a ‘quasi-steady’ state where time acts more like a parameter than an independent variable (Leal Reference Leal2007). Now, to see why this is a good approximation for the present case which involves an unsteady process of droplet deformation, we need to compare two distinct time scales.
(i) Time scale for significant droplet deformation (
$t_1$): this refers to the time required for a change in the droplet deformation which is comparable to the characteristic length scale of the problem, which is the droplet radius ($a$). For the current problem, this can be estimated by analysing the changes in the droplet shape during the course of its migration. For example, it can be seen from the insets in figure 2 that significant droplet deformation occurs across ${\sim }10$ non-dimensional units of $X_{c}$ (see the first and third insets of the droplet shape in figure 2). Since the axial migration velocity of the droplet $U_{x} \sim 1$ (see figures 9 and 11a ii,b ii), the dimensional scaling law for $t_1$ can be expressed as $t_1 \sim {10a}/{\bar {V}_{c}}$. Note that here the lengths are non-dimensionalized by the droplet radius $a$ and velocities are non-dimensionalized using the centreline velocity for the imposed flow $\bar {V}_{c}$, as explained in detail in § 2.2 of this manuscript.(ii) Time scale for evolution of fluid flow from one steady state to another (
$t_2$): in the low-Reynolds-number limit, this is equal to the diffusive time scale which is ${\sim }{a^2}/{\nu }$ (Leal Reference Leal2007), where $\nu$ refers to the kinematic viscosity of the fluids (both the droplet and carrier fluids are considered to be of the same kinematic viscosity in this problem).Therefore, the ratio of the two time scales is
${t_{2}}/{t_{1}} \sim {\bar {V}_{c}a}/{10\nu }$ which is one order of magnitude lower than the Reynolds number. Now, as $Re \ll 1$, it can be concluded that ${t_{2}}/{t_{1}} \ll 1$. Hence, the fluid flow takes negligibly small time to adjust with respect to the time scale for droplet deformation and can be effectively treated to be in a steady state for each deformed configuration of the droplet.
B.2. Some details regarding computational implementation of the numerical algorithm
The computation of interfacial temperature requires the solution of a system of linear equations, and is the most time-consuming part of the simulation. If 
$N$ denotes the number of points and 
$M$ denotes the number of discrete elements over the surface of the droplet, then the total number of steps required for constructing the coefficient matrix for the interfacial temperature is 
$\sim N^{2}M$. The resulting system of linear equations is solved using the DGESV subroutine in LAPACK (Anderson et al. Reference Anderson1999), which is parallelized in the ATLAS (Whaley & Petitet Reference Whaley and Petitet2005) package. The computation of interfacial velocity takes 
$\sim NM$ steps. The code was developed in-house using the FORTRAN programming language and the mesh discretization and integration quadrature was adopted from the BEMLIB (Pozrikidis Reference Pozrikidis2002) package. The code was fully parallelized using OPENMP (Dagum & Menon Reference Dagum and Menon1998) directives and run on a CPU having the i7-8700k processor. Each simulation required approximately 4 GB of RAM and took nearly eight days for completion.
Open access
