1. Introduction
The dynamics of suspending droplets and bubbles gained significant attention due to its vast range of applications in biomedical technology, industrial operations, bio-engineering products and microfluidic devices (Karabelas Reference Karabelas1977; Kaushal & Tomita Reference Kaushal and Tomita2002; Stone, Stroock & Ajdari Reference Stone, Stroock and Ajdari2004; Huebner et al. Reference Huebner, Sharma, Srisa-Art, Hollfelder, Edel and Demello2008; Teh et al. Reference Teh, Lin, Hung and Lee2008; Baroud, Gallaire & Dangla Reference Baroud, Gallaire and Dangla2010; Seemann et al. Reference Seemann, Brinkmann, Pfohl and Herminghaus2011; Stan et al. Reference Stan, Guglielmini, Ellerbee, Caviezel, Stone and Whitesides2011; Zhu & Fang Reference Zhu and Fang2013). In the case of cell analysis and droplet PCR, uniformly sized droplets (millimetre or sub-millimetre) are separated from the base sample, and that process depends on the flow geometry, fluid properties and external fields (such as magnetic, electric, surfactant and temperature). Droplets generally exhibit rich science in shape deformation, cross-migration, encapsulation, break-up and merging (Ahn et al. Reference Ahn, Kerbage, Hunt, Westervelt, Link and Weitz2006; Link et al. Reference Link, Grasland-Mongrain, Duri, Sarrazin, Cheng, Cristobal, Marquez and Weitz2006; Anna Reference Anna2016; Bandopadhyay et al. Reference Bandopadhyay, Mandal, Kishore and Chakraborty2016). During the droplet sorting technique, external fields are applied to create unbalanced surface forces to achieve lateral motion and separation.
The migration of droplets can also be managed without using external fields. There are droplets that migrate due to their internal activity, called ‘active droplets’ (Ramos, Cordero & Soto Reference Ramos, Cordero and Soto2020; Grauer et al. Reference Grauer, Schmidt, Pineda, Midtvedt, Löwen, Volpe and Liebchen2021; Rajabi et al. Reference Rajabi, Baza, Turiv and Lavrentovich2021; Shankar, Raju & Mahadevan Reference Shankar, Raju and Mahadevan2022). Studies on Escherichia coli encapsulated in an emulsion droplet showed motility due to internal heat generation by the micro-organism (Rothbaum & Stone Reference Rothbaum and Stone1961), which can be modelled as a droplet subjected to an internal thermal singularity. An important challenge is to find a way to control the droplet propulsion and the dependence based on the nature of singularities present (Marchetti Reference Marchetti2012; Shankar et al. Reference Shankar, Raju and Mahadevan2022).
Are the droplet dynamics not understood in the presence of thermal fields earlier? Partly yes, in the form of uniform or non-uniform thermal gradients in the ambient field, but without any internal thermal heating to mimic living cells cultured situation. First, we overview earlier mathematical findings on thermocapillary-based droplet dynamics, and then we proceed to explain our work on the thermal singularities. Young, Goldstein & Block (Reference Young, Goldstein and Block1959) have found that droplets with small Reynolds number (
$Re\approx 0$) migrate from the cold to the hot side due to thermocapillary forces. Bratukhin (Reference Bratukhin1975) has obtained a flow field solution up to 
$O(Re)$, when a constant temperature gradient is present at the far field. Thompson, DeWitt & Labus (Reference Thompson, DeWitt and Labus1980) have extended the solutions up to 
$O(Re^2)$ and found far-field boundary conditions are difficult to satisfy similar to the Whitehead paradox. Subramanian (Reference Subramanian1983) has used the method of matched asymptotic expansions in terms of the Marangoni number (
$Ma$) for the Stokes flow, and found the droplet velocity up to 
$O(Ma^2)$. Under zero gravity conditions, Balasubramaniam & Chai (Reference Balasubramaniam and Chai1987) have shown that the solution obtained by Young et al. (Reference Young, Goldstein and Block1959) is an exact solution of the Navier–Stokes equations for any value of 
$Re$ as long as 
$Ma$ is small. Subramanian & Balasubramaniam (Reference Subramanian and Balasubramaniam2001) indicated that the droplet migration velocity is unaffected by the surface shear viscosity but depends on the surface dilatational viscosity. Similar to thermocapillary effects, the surfactant concentration variation on the droplet interface can also alter the droplet velocity (Leal Reference Leal2007). For example, a tiny surfactant presence on the droplet interface reduces the velocity (Sadhal & Johnson Reference Sadhal and Johnson1983, Reference Sadhal and Johnson1986). Hanna & Vlahovska (Reference Hanna and Vlahovska2010) have shown that droplets released at an arbitrary location in a Poiseuille flow (one special case of an arbitrary Stokes flow) reach the centreline for large surface Péclet numbers (
$Pe_s$). Schwalbe et al. (Reference Schwalbe, Phelan, Frederick, Vlahovska and Hudson2011) investigated the effects of surface viscosities and found that the presence of interfacial stresses significantly alters the slip velocity and the flow circulation. Dandekar & Ardekani (Reference Dandekar and Ardekani2020) examined the same for small 
$Pe_s$, and found that the droplet cross-stream migration is due to 
$O(Pe_s^2)$ terms. Das, Mandal & Chakraborty (Reference Das, Mandal and Chakraborty2018) have shown that when the interfacial tension is modelled as a linear combination of surfactant and temperature, the droplet exhibits cross-stream migration at 
$O(Pe_s^2)$. Moreover, the direction of this motion changes depending on the thermal Marangoni number 
$Ma_T$. Sharanya, Raja Sekhar & Rohde (Reference Sharanya, Raja Sekhar and Rohde2019) have found that the cross-stream migration increases due to 
$O(Pe_s)$ if the interfacial tension is modelled as a nonlinear combination of surfactant and temperature field. Panigrahi et al. (Reference Panigrahi, Santra, Banuprasad, Das and Chakraborty2021) have shown that the droplet cross-stream migration gets suppressed with the inclusion of shear and dilatational surface viscosities.
To the best of our knowledge, the droplet dynamics with internal thermal singularities (due to the importance of cell separation techniques) has not been modelled yet. In this work, we have investigated a surfactant-laden droplet suspended in an arbitrary non-isothermal Stokes flow subjected to internal thermal singularities in the form of a monopole and a dipole (see figure 1). In our model, we have allowed the droplet thermophysical properties to be different from the outside ambient so that the stress and thermal balance conditions at the interface determine the droplet dynamics. We have varied the monopole location and the dipole orientation inside the droplet, and investigated the overall drag and the corresponding migration velocity. The interfacial tension in our model is a linear function of the temperature and the surfactant concentration. Our investigation highlights that in a specific arbitrary Stokes flow (i.e. a Poiseuille flow), the droplet with internal thermal singularities can show high velocity compared to the results mentioned in the previous works. Further, we have shown that the droplet can show cross-stream migration even if it is released on the Poiseuille flow centreline.
Figure 1. Schematic of a surfactant-laden droplet (of radius 
$a$) moving with velocity 
${\tilde {\boldsymbol {U}}}$ in an ambient Stokes flow. The ambient flow has velocity 
$\tilde {\boldsymbol {v}}_\infty$ and is subjected to a thermal gradient 
$\tilde {\boldsymbol {\nabla }} \tilde {T}_\infty$ (represented as a blue to red shade), coinciding with gravity 
$\boldsymbol {g}$. Based on the droplet centre, a spherical coordinate system 
$(\tilde {r}, {\theta }, {\phi })$ is chosen that connects with a Cartesian coordinate system 
$(\tilde {x}, \tilde {y}, \tilde {z})$ by general transformation rule. The streamwise direction is 
$\tilde {z}$, and the cross-stream directions are 
$\tilde {x}$ (perpendicular to gravity) and 
$\tilde {y}$ (into the plane). In our work, droplets consist of different thermal singularities such as (i) a centred monopole, (ii) an off-centred monopole with a source in the streamwise direction, and (iii) a dipole with an arbitrary orientation. The red circle in the droplet is a hot source, and yellow for a cold sink. Surfactant is represented as a coiled ball on the droplet surface.
In § 2, we explain the mathematical model, corresponding boundary conditions and the parameters involved. In § 3, we outline the solution procedure using the solenoidal decomposition method for small and large 
$Pe_s$. In § 4, we investigate detailed results for droplets suspended in a non-isothermal Poiseuille flow. A summary with future scope is given in § 5.
2. Mathematical formulation
2.1. Description of the mathematical model
In figure 1, we have shown the schematic of a surfactant-laden droplet of radius 
$a$, suspended in an ambient flow that is subjected to an external temperature gradient 
$\tilde {\boldsymbol {\nabla }}\tilde {T}_\infty$. The density, dynamic viscosity, thermal conductivity and thermal diffusivity internal to the droplet are 
$\rho ^i$, 
$\mu ^i$, 
$\kappa ^i$ and 
$\alpha ^i$, respectively. The same quantities external to the droplet are 
$\rho ^e$, 
$\mu ^e$, 
$\kappa ^e$ and 
$\alpha ^e$, respectively. We define two coordinate systems: (1) spherical 
$(\tilde {r},\theta,\phi )$ to find the detailed solutions, and (2) Cartesian 
$(\tilde {x},\tilde {y},\tilde {z})$ to show results on droplet velocity. Here, 
$\theta$ is the angle made by the vector 
$\tilde {\boldsymbol {r}}$ with the 
$\tilde {z}$ axis, and 
$\phi$ is the angle made by the projection of 
$\tilde {\boldsymbol {r}}$ on the 
$\tilde {x}\unicode{x2013}\tilde {y}$ plane with the 
$\tilde {x}$ axis. Consequently, we can express the relationships as 
$\tilde {x}=\tilde {r}\sin {\theta }\cos {\phi }$, 
$\tilde {y}=\tilde {r}\sin {\theta }\sin {\phi }$ and 
$\tilde {z}=\tilde {r}\cos {\theta }$. The gravitational force 
$\boldsymbol {g}$ acts along the 
$\tilde {z}$ direction. The unit vectors along the 
$\tilde {x}$, 
$\tilde {y}$ and 
$\tilde {z}$ directions are 
$\hat {i}$, 
$\hat {j}$ and 
$\hat {k}$, respectively. On a co-moving frame attached to the droplet centre, the mass, momentum and energy equations that describe internal and external flow fields satisfy (Leal Reference Leal2007; Subramanian & Balasubramaniam Reference Subramanian and Balasubramaniam2001),
$$\begin{gather} \tilde{\boldsymbol{\nabla}}\boldsymbol{\cdot}\tilde{\boldsymbol{u}}^j =0, \end{gather}$$
$$\begin{gather}\rho^j\left(\frac{\partial \tilde{\boldsymbol{u}}^j}{\partial \tilde{t}}+\tilde{\boldsymbol{u}}^j\boldsymbol{\cdot}\tilde{\boldsymbol{\nabla}}\tilde{\boldsymbol{u}}^j\right) =\tilde{\boldsymbol{\nabla}}\boldsymbol{\cdot}\tilde{\boldsymbol{{\tau}}}^j, \end{gather}$$
$$\begin{gather}\frac{\partial \tilde{T}^j}{\partial \tilde{t}}+ \tilde{\boldsymbol{u}}^j\boldsymbol{\cdot}\tilde{\boldsymbol{\nabla}} \tilde{T}^j=\alpha^j\,\tilde{\nabla}^2 \tilde{T}^j, \end{gather}$$
where 
$j=i$ (for internal) and 
$e$ (for external). The other variables are as follows: 
$\boldsymbol {\tilde {u}}^j$ is the relative velocity of flow with respect to the droplet centre velocity 
$\tilde {\boldsymbol {U}}$ (i.e. 
$\boldsymbol {\tilde {u}}^j=\boldsymbol {\tilde {v}}^j-\tilde {\boldsymbol {U}}$, where 
$\boldsymbol {\tilde {v}}^j$ is the flow velocity in the laboratory reference frame), 
$\tilde {\boldsymbol {\tau }}^j$ is the stress tensor, 
$\tilde {T}^j$ is the temperature, and 
${\partial }/{\partial \tilde {t}}$ is the partial derivative with respect to time. For Newtonian fluids (as we use), the stress tensor is modelled as 
${{\tilde {\boldsymbol {\tau }}}}^j=-\tilde {p}^j\boldsymbol {I}+\mu ^j[\tilde {\boldsymbol {\nabla }}\tilde {\boldsymbol {u}}^j+({\tilde {\boldsymbol {\nabla }}}\tilde {\boldsymbol {u}}^j)']$, where 
$\tilde {p}^j$ is the pressure, 
$\boldsymbol {I}$ is the identity tensor, 
${\tilde {\boldsymbol {\nabla }}}$ is the gradient operator, and 
$'$ denotes the transpose. The surfactant concentration 
$\tilde {\varGamma }$ at the droplet interface 
$\tilde {r}=a$ satisfies (Stone & Leal Reference Stone and Leal1990; Stone Reference Stone1990)
\begin{equation} \frac{\partial\tilde{\varGamma}}{\partial \tilde{t}}+\tilde{\boldsymbol{\nabla}}_s\boldsymbol{\cdot}(\tilde{\boldsymbol{u}}_s\tilde{\varGamma})=D_s\, \tilde{\boldsymbol{\nabla}}_s^2\tilde{\varGamma}, \end{equation}
where 
$\tilde {\boldsymbol {u}}_s$ is the relative velocity at the surface, and 
$D_s$ is the surfactant diffusivity. The surface gradient operator is 
$\tilde {\boldsymbol {\nabla }}_s \equiv \tilde {\boldsymbol {\nabla }}-\tilde {\boldsymbol {n}}(\tilde {\boldsymbol {n}} \boldsymbol {\cdot } \tilde {\boldsymbol {\nabla }})$, for a surface described by a normal vector 
$\tilde {\boldsymbol {n}}$, and a tangent vector 
$\tilde {\boldsymbol {t}}$. The interface conditions correspond to: (1) vanishing normal velocity, 
${\tilde {\boldsymbol {u}}}^e \boldsymbol {\cdot } \tilde {\boldsymbol {n}}={\tilde {\boldsymbol {u}}}^i\boldsymbol {\cdot } \tilde {\boldsymbol {n}}=0$; (2) continuity in tangential velocity, 
${\tilde {\boldsymbol {u}}}^e \boldsymbol {\cdot } \tilde {\boldsymbol {t}}-\tilde {\boldsymbol {u}}^i\boldsymbol {\cdot }\tilde {\boldsymbol {t}}=0$; and (3) jump in the tangential stress, 
${{\tilde {\boldsymbol {\tau }}}}^e\boldsymbol {\cdot }{\tilde {\boldsymbol {n}} \boldsymbol {\cdot }\tilde {\boldsymbol {t}}}-{\tilde {\boldsymbol {\tau }}}^i\boldsymbol {\cdot }{\tilde {\boldsymbol {n}} \boldsymbol {\cdot }\tilde {\boldsymbol {t}}}=-\tilde {\boldsymbol {\nabla }}_s \tilde {\sigma } \boldsymbol {\cdot } \tilde {\boldsymbol {t}}$. The interfacial tension (
$\tilde {\sigma }$) in our work is modelled as 
$\tilde {\sigma }(\tilde {T})=\tilde {\sigma }_0-\gamma (\tilde {T}-\tilde {T_0})-R_g\tilde {T}_0\tilde {\varGamma }$ for the dilute concentrations and moderate temperature differences (Homsy & Meiburg Reference Homsy and Meiburg1984; Kim & Subramanian Reference Kim and Subramanian1989; Leal Reference Leal2007). The equilibrium interfacial tension (
$\tilde {\sigma }_0$) is considered at reference temperature 
$\tilde {T}_0$, 
$\gamma$ is a positive constant, and 
$R_g$ is the ideal gas constant. The temperature field at the interface satisfies 
$\tilde {T}^e =\tilde {T}^i$ and 
$\kappa ^e ({\partial \tilde {T}^e}/{\partial \tilde {r}}) = \kappa ^i ({\partial \tilde {T}^i}/{\partial \tilde {r}})$. The far-field conditions (at 
$\tilde {r}\to \infty$) for the velocity, pressure and temperature are 
$\tilde {\boldsymbol {u}^e}\to \tilde {\boldsymbol {u}}_\infty =\tilde {\boldsymbol {v}}_\infty -\tilde {\boldsymbol {U}}$, 
$\tilde {p}^e \to \tilde {p}_\infty$ and 
$\tilde {T}^e \to \tilde {T}_\infty$, respectively. Internal to the droplet (at 
$\tilde {r}\to 0$), both the velocity 
$\tilde {\boldsymbol {u}}^i$ and the pressure 
$\tilde {p}^i$ are bounded, and the temperature satisfies, 
$\tilde {T}^i \to \tilde {T}_s$, for a prescribed thermal source that generates temperature 
$\tilde {T}_s$.
2.2. Control parameters and dimensionless governing equations
To non-dimensionalize the variables, we need to define relevant characteristic quantities. For example, in a cylindrical Poiseuille flow of radius 
$\tilde {R}_0$ with the centreline flow velocity 
${U}_c$, we define the dimensionless variables as 
$\boldsymbol {{u}}^j={\tilde {\boldsymbol {u}}^j/ {{U}_c}}$, 
${\boldsymbol {\nabla }}=a\,\tilde {\boldsymbol {\nabla }}$, 
$\boldsymbol {{r}}={\tilde {\boldsymbol {r}}/ a}$, 
${p}^j={\tilde {p}^j/ (\mu ^e {U}_c}/a)$, 
$\sigma =\tilde {\sigma }/\tilde {\sigma }_0$, 
$T^j=(\tilde {T}^j-\tilde {T}_0)/a\,|\tilde {\boldsymbol {\nabla } }\tilde {T}^\infty |$, 
$\varGamma =\tilde {\varGamma }/\tilde {\varGamma }_0$ and 
$t=\tilde {t}/(\tilde {R}_0/{U}_c)$. Here, 
$\tilde {R}_0/{U}_c$ denotes the advection time scale. With these non-dimensional variables, the governing equations (2.1) reduce to
$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}^j =0, \end{gather}$$
$$\begin{gather}Re\,(\rho^j/\rho^e)\left[\frac{a}{\tilde{R}_0}\,\frac{\partial \boldsymbol{u}^j}{\partial t}+(\boldsymbol{u}^j\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}^j)\right]={-}\boldsymbol{\nabla} p^j + (\mu^j/\mu^e)\,\nabla^2 \boldsymbol{u}^j, \end{gather}$$
$$\begin{gather}Pe_T\left[\frac{a}{\tilde{R}_0}\,\frac{\partial T^j}{\partial t}+ (\boldsymbol{u}^j\boldsymbol{\cdot}\boldsymbol{\nabla} T^j)\right]=(\alpha^j/\alpha^e)\,\nabla^2 T^j, \end{gather}$$
where 
$Re\ (\equiv \rho ^e {{U}_c} a / \mu ^e)$ is the Reynolds number, and 
$Pe_T\ (\equiv {{U}_c}a/\alpha ^e)$ is the thermal Péclet number. Note that index 
$j=i$ for the internal flow and 
$e$ for the external flow.
The surfactant transport equation (2.2) in dimensionless form is
\begin{equation} Pe_s \left[\frac{a}{\tilde{R}_0}\,\frac{\partial \varGamma}{\partial t} + \boldsymbol{\nabla}_s \boldsymbol{\cdot}(\varGamma \boldsymbol{u}_s)\right]=\nabla^2_s\varGamma, \end{equation}
where 
$Pe_s\ (\equiv {{U}_c}a/D_s)$ is the surface Péclet number. The interface conditions (at 
$r=1)$ in dimensionless form are: (1) 
${\boldsymbol {{u}}}^e \boldsymbol {\cdot } {\tilde {\boldsymbol {n}}}={\boldsymbol {{u}}}^i \boldsymbol {\cdot } {\tilde {\boldsymbol {n}}}=0$, (2) 
${\boldsymbol {{u}}}^e \boldsymbol {\cdot } {\tilde {\boldsymbol {t}}} -\boldsymbol {{u}}^i\boldsymbol {\cdot }{\tilde {\boldsymbol {t}}}=0$, and (3) 
${\boldsymbol {{\tau }}}^e\boldsymbol {\cdot }{{\tilde {\boldsymbol {n}}}\boldsymbol {\cdot }{\tilde {\boldsymbol {t}}}}-$ 
$\mu \, {\boldsymbol {{\tau }}}^i\boldsymbol {\cdot }{{\tilde {\boldsymbol {n}}}\boldsymbol {\cdot }{\tilde {\boldsymbol {t}}}} = \boldsymbol {\nabla }_s( Ma_T\,T^e + Ma_\varGamma \,\varGamma ) \boldsymbol {\cdot } \tilde {\boldsymbol {t}}$, where 
$\mu =\mu ^i/\mu ^e$, 
$Ma_T={\gamma \,|\tilde {\boldsymbol {\nabla } }\tilde {T}_\infty |\,a}/{\mu ^e {{U}_c}}$ is the thermal Marangoni number, and 
$Ma_\varGamma ={R_g \tilde {T}_0\tilde {\varGamma }_0}/{\mu ^e {{U}_c}}$ is the surfactant Marangoni number. Note that the tangential stress balance condition is simplified using 
$\sigma =1-Ca\,(Ma_T\,T+Ma_\varGamma \, \varGamma )$, where 
$Ca=\mu ^e U_c/\tilde {\sigma }_0$ is the capillary number. The thermal field at the interface satisfies 
${T}^e ={T}^i$ and 
${\partial {T}^e}/{\partial {r}} =\kappa ({\partial {T}^i}/{\partial r})$, where 
$\kappa =\kappa ^i/\kappa ^e$. The far-field conditions (at 
${r} \to \infty )$ are 
$\boldsymbol {{u}}^e \to \boldsymbol {{u}}_\infty =\boldsymbol {{v}}_\infty -\boldsymbol {{U}}$, 
${p}^e \to {p}_\infty$, and 
${T}^e \to {T}_\infty$. Internal to the droplet (as 
$r\to 0$), we have 
${\boldsymbol {u}}^i$ and 
${p}^i$ bounded with 
${T}^i \to {T}_s$, the thermal singularity.
To give an idea about the mentioned dimensionless numbers, we discuss two examples: (1) a thermocapillary droplet (Nallani & Subramanian Reference Nallani and Subramanian1993) of radius 
$50$–
$146\,\mathrm {\mu }{\rm m}$ suspended in an ambient fluid with 
$\rho ^e=955\,\text {kg}\,\text {m}^{-3}$, 
$\mu ^e= 0.0478\,\text {N s}\,\text {m}^{-2}$, 
$\kappa ^e=0.1\,\text {W}\,\text {m}^{-1}\,\text {K}^{-1}$; and (2) a deionized water droplet (Panigrahi et al. Reference Panigrahi, Santra, Banuprasad, Das and Chakraborty2021) suspended in another medium with the properties 
$\rho ^e= 971\,\text {kg}\,\text {m}^{-3}$, 
$\mu ^e= $ 
$0.04855\,\text {N s}\,\text {m}^{-2}$, 
$\alpha ^e=7\times 10^{-8}\,\text {m}^{2}\,\text {s}^{-1}$, 
$D_s=10^{-11}$–
$10^{-8}\,\text {m}^2\,\text {s}^{-1}$ and 
$\tilde {\varGamma }_0=10^{-10}$– 
$10^{-6}\,\text {mol}\,\text {m}^{-2}$. Based on these examples, the corresponding parameter ranges are 
$\rho \ (=\rho ^i/\rho ^e)\sim O(10^{-2})$–
$O(1)$, 
$Re\sim O(10^{-4})$, 
$a/\tilde {R}_0 \sim O(10^{-2})$, 
$\mu \sim O(10^{-2})$–
$O(10^2)$, 
$Pe_T \sim O(10^{-2})$, 
$\alpha \ (=\alpha ^i/\alpha ^e) \sim O(10^{-1})$–
$O(1)$, 
$Pe_s\sim O(10^{-1})$–
$O(10^{2})$, 
$Ma_T\sim O(10^{-2})$ –
$O(10^{1})$, 
$Ma_\varGamma \sim O(10^{-2})$–
$O(10^{2})$ and 
$Ca\sim O(10^{-4})$. Further, this study focuses on droplet dynamics in situations where droplets may have a density different from the surrounding medium. The presence of gravity leads to a buoyancy effect, causing droplets to undergo motion. When gravity acts perpendicular to the flow, the cross-stream velocity of a buoyant droplet is governed primarily by the buoyancy force, and the effects of flow and temperature field are negligible. Consequently, investigating buoyant droplets in such systems is considered less important. However, here we explore specifically droplet migration when gravity acts in the same direction as the imposed flow (as depicted in figure 1). This scenario presents an intriguing area of research to understand the dynamics and behaviour of buoyant droplets under these conditions. Thus we have the Bond number 
$Bo=|\rho ^i-\rho ^e|\,g a^2/\tilde {\sigma }_0$, which is of 
$O(10^{-4})$. Note that the normal stress balance determines the droplet shape (Hetsroni & Haber Reference Hetsroni and Haber1970), and due to small 
$Ca \sim O(10^{-4})$, we have assumed that the droplet is spherical.
Considering 
$Re \sim O(10^{-4})$, we have neglected terms with 
$(a/\tilde {R}_0)\,Re$ since they are significantly smaller than 
$Re$ itself. Similarly, in the energy equation, considering the value of 
$Pe_T \sim O(10^{-2})$, we have disregarded terms with 
$(a/\tilde {R}_0)\,Pe_T$ as they are much smaller compared to 
$Pe_T$. Thus the governing equations become
$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}^j =0, \end{gather}$$
$$\begin{gather}-\boldsymbol{\nabla} p^j + ({\mu^j}/{\mu^e})\,\nabla^2 \boldsymbol{u}^j=0, \end{gather}$$
$$\begin{gather}\nabla^2 T^j =0. \end{gather}$$ With reference to the surfactant transport equation given in (2.2), the order of magnitude analysis of the individual terms indicates that when the advection time scale is restricted as 
$\tilde {R}_0/{U}_c \gg l_c/{U}_c$, the governing equation reaches a quasi-steady state. Accordingly, for the choice 
$l_c=a \ll \tilde {R}_0$ (small droplets), we ensure that 
$a/U_c \ll \tilde {R}_0/{U}_c$. Consequently, we can neglect the transient terms, resulting in the simplified surfactant transport equation
\begin{equation} Pe_s\,\boldsymbol{\nabla}_s \boldsymbol{\cdot}(\varGamma \boldsymbol{u}_s)=\nabla^2_s\varGamma. \end{equation}
We have considered that the surfactant concentration satisfies the condition 
$\int _{\phi =0}^{2{\rm \pi} }\int _{\theta =0}^{{\rm \pi} } \varGamma (\theta,\phi )\sin {\theta }\, {\rm d}\theta \, {\rm d}\phi =4{\rm \pi}$. We can solve (2.6) analytically for a few special cases like small and large 
$Pe_s$. The surface diffusion dominates the transport for 
$Pe_s \ll 1$, and the surface convection dominates the transport for 
$Pe_s \gg 1$.
3. Solution procedure
Irrespective of the values of 
$Pe_s$, the temperature field is independent of the flow field and the surfactant concentration; thus we first solve the temperature equations with the prescribed boundary and interface conditions. Using asymptotic perturbation for small and large 
$Pe_s$, we solve the surfactant concentration and flow fields simultaneously (see figure 2). At each order of perturbation, the Stokes flow equations are solved using the solenoidal decomposition method known as the double-curl method, which is described briefly in § 3.2. In the subsequent subsections, we provide the necessary details about the solution procedure to determine the drag on the surface of the droplet and the corresponding migration velocity.
Figure 2. Schematic of a general procedure to solve for velocity 
$\boldsymbol {u}$ and surfactant concentration 
$\varGamma$ from the temperature field 
$T$, for (a) small-
$Pe_s$ limit, and (b) large-
$Pe_s$ limit. Step-by-step progression in arrows, where each quantity is solved before moving on to the next one. Temperature and previous quantities are utilized at each step, and 
$Ma_{\varGamma }$ is the surfactant-based Marangoni number.
3.1. Thermal field
Both inside and outside of the droplet satisfy the Laplace equation (2.5c) with a complete general solution in the form
\begin{align} T^i =\sum_{n=0}^{\infty} \left(\mathscr{a}_n r^n + \frac{\mathscr{b}_n }{r^{n+1}}\right) \mathcal{S}_n(\theta,\phi)\quad \text{and}\quad T^e =\sum_{n=0}^{\infty} \left( \mathscr{c}_n r^n + \frac{\mathscr{d}_n }{r^{n+1}}\right) \mathcal{S}_n(\theta,\phi), \end{align}
where 
$\mathcal {S}_n(\theta,\phi )=\sum _{m=0}^{n} P_n^m(\cos {\theta })({A}_{nm}\cos {m\phi }+{B}_{nm} \sin {m\phi })$ is the 
$n$th spherical harmonic, 
$P_n^m(\cos {\theta })$ is the associated Legendre polynomial, and 
$A_{nm}$, 
$B_{nm}$ are the known scalars (Sneddon Reference Sneddon2006). In (3.1a,b), 
$\mathscr {a}_n$ and 
$\mathscr {d}_n$ are to be determined using the thermal interface conditions at 
$r=1$, whereas 
$\mathscr {b}_n$ comes from the internal thermal singularity and 
$\mathscr {c}_n$ from the ambient far field, respectively. We consider the chosen thermal singularity and the far-field condition to satisfy the Laplace equation such that
\begin{equation} T_s =\sum_{n=0}^{\infty} \frac{\mathscr{b}_n }{r^{n+1}}\,\mathcal{S}_n(\theta,\phi)\ \text{at } r=0,\quad \text{and} \quad T_{\infty} =\sum_{n=0}^{\infty} \mathscr{c}_n r^n\, \mathcal{S}_n(\theta,\phi)\ \text{at } r\to\infty. \end{equation}
Further, 
$T_s$ varies depending on the type of singularity chosen. For example, a monopole of strength 
$S$ located at the droplet centre (
$r=0$) corresponds to 
$T_s=S/r$, which is the free-space Green's function in 
$\mathbb {R}^3$ of the Laplace equation. If the monopole is off-centred by a distance 
$d$ along the flow direction with 
$0< d<1$, then we have 
$T_s=\sum _{n=0}^{\infty }S({d^n }/{r^{n+1}})\,P_n(\cos {\theta })$. A dipole (which is a combination of a sink and a source of equal strength 
$S$) located at the droplet centre and oriented at an arbitrary angle in spherical coordinates creates a thermal field,
\begin{equation} T_s=\frac{1}{r^2}[P_1^0(\cos{\theta})\,(S \cos\theta)+P_1^1(\cos{\theta})\,(S\sin\theta\cos{\phi}+S \sin \theta\sin{\phi})]. \end{equation}
We have estimated the required non-zero coefficients as follows: for a centred monopole, 
$\mathscr {b}_0 A_{00} =S$; for an off-centred monopole, 
$\mathscr {b}_n A_{n0}={S}d^n$ for all 
$n\geq 0$; and for a dipole, we have 
$\mathscr {b}_1A_{11}=S\sin \theta \cos {\phi }$, 
$\mathscr {b}_1B_{11}=S \sin \theta \sin {\phi }$ and 
$\mathscr {b}_1A_{10}=S \cos \theta$ for fixed 
$\theta$ and 
$\phi$. All other coefficients not explicitly mentioned here are zero.
The far-field temperature gradient is essential in determining the particular flow behaviour and heat transfer. A constant temperature gradient simplifies the thermal boundary conditions at infinity. Further, it gives a more accurate solution by ensuring that the flow field is symmetric and uniform. Following Choudhuri & Raja Sekhar (Reference Choudhuri and Raja Sekhar2013) and Sharanya & Raja Sekhar (Reference Sharanya and Raja Sekhar2015), we can write the temperature field corresponding to the constant far-field temperature gradient as 
$T_\infty =Gr\cos {\theta }$, where 
$G<0$ for a decrease in temperature along the flow direction, and 
$G>0$ otherwise. Consequently, we find that 
$\mathscr {c}_1A_{10}=G$ and the remaining 
$\mathscr {c}_n$ are zero, for all 
$n \neq 1$. By applying the interface conditions, we get the coefficients
\begin{align} \mathscr{a}_n = \frac{(2n+1)\mathscr{c}_n+(n+1)(\kappa-1)\mathscr{b}_n}{(n\kappa+n+1)}\quad \text{and} \quad \mathscr{d}_n = \frac{n(1-\kappa)\mathscr{c}_n+\kappa(2n+1)\mathscr{b}_n}{(n\kappa+n+1)}. \end{align}This determines the thermal field interior and exterior of the droplet due to thermal singularities.
3.2. Double-curl strategy to obtain the velocity field
It is important to note that Lamb's general solution (Lamb Reference Lamb1924) helps to solve any arbitrary steady Stokes flow, in which velocity and pressure fields are represented as a series of an infinite number of harmonic scalar functions. As an alternative, Palaniappan et al. (Reference Palaniappan, Nigam, Amaranath and Usha1992), Raja Sekhar et al. (Reference Raja Sekhar, Rao, Padmavathi and Amaranath1995) and Padmavathi, Raja Sekhar & Amaranath (Reference Padmavathi, Raja Sekhar and Amaranath1998) developed the double-curl method, which is a complete general solution framework to solve the Stokes equations, in which the velocity and the pressure fields are represented in terms of a scalar biharmonic function and a scalar harmonic function. Correspondingly, the velocity and the pressure fields of the Stokes equation are represented as 
$\boldsymbol {u} = \boldsymbol {\nabla } \times \boldsymbol {\nabla } \times (\boldsymbol {r} A ) + \boldsymbol {\nabla }\times ( \boldsymbol {r} B )$ and 
$p = p_\infty + ({\partial }/{\partial r}) (r\,\nabla ^2 A)$. The scalars 
$A$ and 
$B$ satisfy 
$\nabla ^4 A=0$ and 
$\nabla ^2 B=0$. The general solutions for 
$A$ and 
$B$ are expressed in terms of the basis vectors, 
$\{ r^n, { r^{n+2}}, { r^{-(n+1)}}, {r^{-(n-1)}}\}\,\mathcal {R}_n(\theta, \phi )$ and 
$\{ r^n, {r^{-(n+1)}}\}\,\mathcal {T}_n(\theta, \phi )$ for all 
$n\geq 0$, respectively. Here, 
$\mathcal {R}_n(\theta,\phi )=\sum _{m=0}^{n} P_n^m(\cos {\theta })\,(C_{nm}\cos {m\phi }+D_{nm} \sin {m\phi })$ and 
$\mathcal {T}_n(\theta,\phi )=\sum _{m=0}^{n} P_n^m(\cos {\theta })\,(E_{nm}\cos {m\phi }+F_{nm} \sin {m\phi })$ are the spherical harmonics, where 
$C_{nm}$, 
$D_{nm}$, 
$E_{nm}$ and 
$F_{nm}$ are constants (which are determined from the spherical harmonics at the far field). Using the superscripts 
$e$ and 
$i$ to denote external and internal quantities, respectively, the interface conditions (at 
$r=1$) in terms of the scalars are given by 
$A^e=0$, 
$A^i=0$, 
${\partial A^e }/{\partial r} ={\partial A^i }/{\partial r}$, 
$B^e=B^i$, 
${\partial ^2 A^e }/{\partial r^2}-\mu ({\partial ^2 A^i }/{\partial r^2}) =Ma_T\, T^e + Ma_\varGamma \, \varGamma$ and 
$({\partial }/{\partial r})({B^e }/{r}) =\mu ({\partial }/{\partial r})({B^i }/{r})$. We have 
$A^i < \infty$, 
$B^i < \infty$ for 
$r \to 0$, and from the far-field conditions, we infer 
$A^e \to A_\infty$ and 
$B^e \to B_\infty$. In series form, we have
\begin{equation} A_\infty = \sum_{n=0}^{\infty} (\alpha^\infty_n r^n + {\hat\alpha^\infty_n r^{n+2}})\,\mathcal{R}_n(\theta, \phi)\quad \text{and} \quad B_\infty =\sum_{n=0}^{\infty} \chi^\infty_n r^n\,\mathcal{T}_n(\theta,\phi), \end{equation}
where 
$\alpha ^\infty _n$, 
$\hat \alpha ^\infty _n$ and 
$\chi ^\infty _n$ are known coefficients determined by the choice of a particular Stokes flow at the far field. The resultant flow interior (
$r<1$) and exterior (
$r>1$) of the drop follow:
\begin{gather} A^i =\sum_{n=0}^{\infty} (\delta_n r^n + {\hat\delta_n r^{n+2}})\,\mathcal{R}_n(\theta, \phi)\quad \text{and} \quad B^i =\sum_{n=0}^{\infty} \gamma_n r^n\,\mathcal{T}_n(\theta, \phi), \end{gather}
\begin{gather} A^e =A_\infty+\sum_{n=0}^{\infty} \Bigg(\frac{\beta_n }{r^{n+1}}+\frac{\hat\beta_n }{r^{n-1}}\Bigg)\mathcal{R}_n(\theta, \phi)\quad \text{and} \quad B^e =B_\infty+\sum_{n=0}^{\infty} \frac{\zeta_n }{r^{n+1}}\,\mathcal{T}_n(\theta, \phi), \end{gather}
where 
$\delta _n$, 
$\hat \delta _n$, 
$\gamma _n$, 
$\beta _n$, 
$\hat \beta _n$ and 
$\zeta _n$ are the unknown coefficients that are to be determined using the interface conditions. Consequently, we can compute the velocity and pressure fields. The next task is to compute the migration velocity of the droplet. Towards this, we use the force-free condition (Happel & Brenner Reference Happel and Brenner1981) on the droplet, which corresponds to the balance between the drag (
$\boldsymbol {D}^H$) and the buoyancy force as
\begin{equation} \underbrace{\int_{\varOmega} {\boldsymbol{\tau}^e} \boldsymbol{\cdot} \tilde{\boldsymbol{n}}\, {\rm d}{\varOmega}}_{{\substack{\textit{Hydrodynamic drag}}}}+ \underbrace{{\frac{4{\rm \pi}}{3}\,\frac{Bo}{Ca}\,\boldsymbol{\omega}}}_{\textit{Buoyancy}}=0, \end{equation}
where 
${\rm d}\varOmega$ is the droplet surface element, 
$\boldsymbol {\omega }$ denotes the direction of the buoyancy force, and the denser fluid among the droplet and the surrounding phase is determined by the magnitude 
${\omega }$. If 
$\omega =1$, then the drop is denser than the suspending phase (
$\rho ^i>\rho ^e$), and if 
$\omega =-1$, then the surrounding fluid is dense (
$\rho ^i<\rho ^e$).
3.3. A regular perturbation of the surfactant transport and droplet velocity
To compute drag and droplet velocity, we need surfactant distribution and flow field. The surfactant transport equation (2.6) and the scalars 
$A^j$ and 
$B^j$ are coupled via the interface conditions. A complete analytical solution for the entire range of 
$Pe_s$ is not possible, hence we construct an asymptotic solution for 
$Pe_s\ll 1$ and 
$Pe_s\gg 1$ as a regular perturbation problem.
3.3.1. For small surface Péclet number (
$Pe_s \ll 1$)
In the limit 
$Pe_s \ll 1$, the surfactant concentration and velocity can be expanded in terms of a regular perturbation in 
$Pe_s$ as (Pak, Feng & Stone Reference Pak, Feng and Stone2014; Sharanya et al. Reference Sharanya, Raja Sekhar and Rohde2019)
\begin{equation} \left[\varGamma, \boldsymbol{u}^j\right] =\underbrace{\left[ \varGamma_0, \boldsymbol{u}^j_0\right]}_{\textit{Leading order}}+Pe_s \underbrace{\left[\varGamma_1, \boldsymbol{u}^j_1 \right]}_{\textit{First order}}+Pe_s^2 \underbrace{\left[\varGamma_2, \boldsymbol{u}^j_2 \right]}_{\textit{Second order}}+O(Pe_s^3). \end{equation}From (2.6) and (3.9), the surfactant transport equation at various orders follows:
$$\begin{gather} O(1){:}\quad \boldsymbol{\nabla}_s^2 \varGamma_0=0, \end{gather}$$
$$\begin{gather}O(Pe_s){:}\quad \boldsymbol{\nabla}_s^2 \varGamma_1= \boldsymbol{\nabla}_s \boldsymbol{\cdot} (\boldsymbol{u}_{s,0}\varGamma_0), \end{gather}$$
$$\begin{gather}O(Pe_s^2){:}\quad \boldsymbol{\nabla}_s^2 \varGamma_2= \boldsymbol{\nabla}_s \boldsymbol{\cdot} (\boldsymbol{u}_{s,0}\varGamma_1+\boldsymbol{u}_{s,1}\varGamma_0), \end{gather}$$
where 
$\boldsymbol {u}_{s,0}$ and 
$\boldsymbol {u}_{s,1}$ are the surface velocities at the leading order and first order in 
$Pe_s$, respectively. From the first order onwards, the surfactant concentration depends on the previous flow field (or the surface velocity) and the surfactant concentration of the previous order, and appears as a non-homogeneous part to the Laplacian as in (3.10b) and (3.10c). The leading-order surfactant transport equation (3.10a), which further satisfies the mass balance condition, follows 
$\varGamma _0=1$. The leading-order velocity and pressure fields are computed using the double-curl strategy by applying the interface and far-field conditions. The first-order surfactant transport equation is solved by using already computed leading-order results. Further details in terms of spherical harmonics are shown in Appendix A.
The leading-order drag 
$\boldsymbol {D}_0^H$ obtained by using (3.8), in the form of Faxen's laws (Subramanian & Balasubramaniam Reference Subramanian and Balasubramaniam2001) is
\begin{align} { \boldsymbol{{D}}^H_0 }&= 6{\rm \pi}\left({\frac{3\mu +2 }{3(\mu+1)}\,[\boldsymbol{u}_{\infty,0}]_0+\frac{\mu }{6(\mu+1)}\,[\nabla^2\boldsymbol{u}_{\infty,0}]_0}+{\frac{2 Ma_T }{3(\mu+1)(\kappa+2)}\,[\boldsymbol{\nabla} T_\infty]_0}\right.\nonumber\\ &\quad \left.{}+{\frac{2\kappa \mathscr{b}_1\,Ma_T }{3(\mu+1)(\kappa+2)}\,\boldsymbol{Y}_{0}}\right), \end{align}
with zero net torque. The detailed expressions of 
$[\boldsymbol {u}_{\infty,0}]_0$, 
$[\nabla ^2\boldsymbol {u}_{\infty,0}]_0$, 
$[\boldsymbol {\nabla } T_\infty ]_0$ and 
$\boldsymbol {Y}_{0}$ are given in (A2) and (A3). Note that the fourth term in the above equation signifies the role of internal thermal singularity.
The leading-order migration velocity of the drop is found by equating the drag force to zero, which gives
\begin{align} \boldsymbol{U}_0 &=
\underbrace{[\boldsymbol{v}_{\infty,0}]_0+\left(\frac{\mu
}{6\mu+4}\right)[\nabla^2\boldsymbol{v}_{\infty,0}]_0}_{\textit{Ambient flow
contribution}} +\underbrace{\frac{2 Ma_T
}{(3\mu+2)(\kappa+2)}\,[\boldsymbol{\nabla}
T_\infty]_0}_{\textit{Thermal gradient
contribution}} \nonumber\\ &\quad
+\underbrace{\frac{2\kappa \mathscr{b}_1\,Ma_T
}{(3\mu+2)(\kappa+2)}\,\boldsymbol{Y}_{0}}_{{\substack{\textit{Internal thermal}\\ \textit{singularity
contribution}}}}+\underbrace{\boldsymbol{U}_{B,0}}_{{\substack{\textit{Buoyancy} \\ \textit{contribution}}}}. \end{align}
Here, 
$\boldsymbol {U}_{B,0}$ is the velocity contribution due to the buoyancy, which depends on the particular ambient flow (a detailed analysis for Poiseuille flow is provided in Appendix B). The structure in (3.12) allows us to compute the droplet velocity for any arbitrary ambient flow and thermal fields in the presence of thermal singularities interior to the droplet. Hence the expression (3.12) for the migration velocity facilitates a more generalized control of the drop migration based on the choice of various parameters involved. At the leading order, we have 
$\varGamma _0=1$, and 
$\boldsymbol {\nabla }_s\varGamma _0=\boldsymbol {0}$, so the surfactant coating is uniform on the drop surface. Hence there is no influence of surfactant-induced Marangoni stresses on the droplet migration velocity.
Following the same approach, we find the first-order drag 
$\boldsymbol {D}^H_1$ and droplet migration velocity 
$\boldsymbol {U}_1$ as
\begin{gather}
\boldsymbol{D}^H_1 = 2{\rm \pi}\left[-\left(\frac{2+3\mu}{\mu+1}\right)\boldsymbol{U}_1
+\frac{2}{3(\mu+1)}\,Ma_\varGamma\,\boldsymbol{Y}_{1}\right],
\end{gather}
\begin{gather}
\boldsymbol{U}_1 =\frac{2}{3(3\mu+2)}\,Ma_\varGamma\,\boldsymbol{Y}_{1}+\boldsymbol{U}_{B,1},
\end{gather}
where the complete expression of the vector 
$\boldsymbol {Y}_{1}$ is given in (A7).
The second-order drag 
$\boldsymbol {D}^H_2$ and the migration velocity 
$\boldsymbol {U}_2$ are
\begin{gather} \boldsymbol{D}_2^H =
2{\rm \pi}\left[-\left(\frac{2+3\mu}{\mu+1}\right)\boldsymbol{U}_2+\frac{2}{3(\mu+1)}\,Ma_\varGamma\,
\boldsymbol{Y}_{2}\right], \end{gather}
\begin{gather}
\boldsymbol{U}_2=\frac{2}{3(3\mu+2)}\,Ma_\varGamma\, \boldsymbol{Y}_{2}+
\boldsymbol{U}_{B,2}, \end{gather}
where the vector 
$\boldsymbol {Y}_{2}$ is given in (A10). Further, 
$\boldsymbol {U}_{B,1}$ and 
$\boldsymbol {U}_{B,2}$ are due to buoyancy and depend on the background flow. The impact of the non-uniform surfactant distribution is reflected in the migration velocity of the droplet at first and second order. A closer look at the structures of 
$\boldsymbol {Y}_1$ and 
$\boldsymbol {Y}_2$ (refer to Appendix A) reveals that the presence of thermal singularities affects the migration velocity. The results given in (3.11), (3.12), (3.13a), (3.13b), (3.14a) and (3.14b) are more generic and can be used for any arbitrary Stokes flow.
3.3.2. For large surface Péclet number ($Pe_s\gg 1$)
For large 
$Pe_s$ (
$\gg 1$), the surfactant transport equation (2.6) simplifies to
\begin{equation} \boldsymbol{\nabla}_s\boldsymbol{\cdot} (\boldsymbol{u}_s\varGamma)=0, \end{equation}
with a thin mass transfer boundary layer structure very close to the drop. For large surfactant Marangoni number (
$Ma_\varGamma \gg 1$), we can expand the surfactant concentration and velocity fields as a regular perturbation with the form (Hanna & Vlahovska Reference Hanna and Vlahovska2010; Schwalbe et al. Reference Schwalbe, Phelan, Frederick, Vlahovska and Hudson2011)
\begin{equation} [\varGamma,\boldsymbol{u}^j] =[\varGamma_0,\boldsymbol{u}^j_0]+Ma_\varGamma^{{-}1}\,[\varGamma_1,\boldsymbol{u}^j_1]+Ma_\varGamma^{{-}2}\,[\varGamma_2,\boldsymbol{u}^j_2] +O(Ma_\varGamma^{{-}3}). \end{equation}
The leading-order surfactant concentration follows 
$\varGamma _0=1$ (Hanna & Vlahovska Reference Hanna and Vlahovska2010; Schwalbe et al. Reference Schwalbe, Phelan, Frederick, Vlahovska and Hudson2011). Substituting (3.16) in (3.15), the governing equations for surfactant concentration at different orders are
$$\begin{gather} O(1){:}\quad \boldsymbol{\nabla}_s\boldsymbol{\cdot}\boldsymbol{u}_{s,0}=0, \end{gather}$$
$$\begin{gather}O(Ma_\varGamma^{{-}1}){:}\quad \boldsymbol{\nabla}_s\boldsymbol{\cdot} \left( \boldsymbol{u}_{s,0} \varGamma_{1}+\boldsymbol{u}_{s,1} \right)=0. \end{gather}$$The flow field is unaffected by the above expansion except for the tangential stress condition at the interface. In leading and first order, the relevant tangential stress conditions are
$$\begin{gather} O(1){:}\quad (\boldsymbol{\tau}^e\boldsymbol{\cdot}{\tilde{\boldsymbol{n}}}\boldsymbol{\cdot}{\tilde{\boldsymbol{t}}})_0-\mu(\boldsymbol{\tau}^i\boldsymbol{\cdot}{\tilde{\boldsymbol{n}}}\boldsymbol{\cdot}{\tilde{\boldsymbol{t}}})_0 =\boldsymbol{\nabla}_s (Ma_T\, T^e+\varGamma_{1})\boldsymbol{\cdot}{\tilde{\boldsymbol{t}}}, \end{gather}$$
$$\begin{gather}O(Ma_\varGamma^{{-}1}){:}\quad (\boldsymbol{\tau}^e\boldsymbol{\cdot}{\tilde{\boldsymbol{n}}}\boldsymbol{\cdot}{\tilde{\boldsymbol{t}}})_{1}-\mu(\boldsymbol{\tau}^i\boldsymbol{\cdot}{\tilde{\boldsymbol{n}}}\boldsymbol{\cdot}{\tilde{\boldsymbol{t}}})_{1} =\boldsymbol{\nabla}_s\varGamma_{2} \boldsymbol{\cdot} {\tilde{\boldsymbol{t}}}. \end{gather}$$Further, the perturbation expansion does not have any impact on the thermal field. Hence the solution of the thermal field is the same as that presented in § 3.1. The flow field follows a procedure similar to that described for the small surface Péclet number case, except for the tangential stress jump condition. The detailed expressions for the surfactant distribution and the velocity field are given in § A.2.
The leading-order drag and the droplet migration velocity corresponding to (3.16) are
\begin{align} \boldsymbol{D}_0^H =6{\rm \pi}([\boldsymbol{u}_{\infty,0}]_0+\tfrac{1}{6}[\nabla^2 \boldsymbol{u}_{\infty,0}]_0)\quad \text{and} \quad \boldsymbol{U}_0 =[\boldsymbol{v}_{\infty,0}]_0+\tfrac{1}{6}[\nabla^2 \boldsymbol{v}_{\infty,0}]_0+\boldsymbol{U}_{B,0}. \end{align}
We can observe that in the limit as 
$\mu \to \infty$, (3.11) and (3.12) reduce to (3.19a,b). In this context, we refer to Manikantan & Squires (Reference Manikantan and Squires2020), where the large 
$Ma_\varGamma$ limit is discussed in detail. Accordingly, in this limit, due to the resistive behaviour of the surfactant molecules, the surface of the drop becomes completely immobilized. This is referred as a two-dimensional incompressible surface. Stone & Masoud (Reference Stone and Masoud2015) have shown that this incompressible surface induces the subphase fluid to flow in planes parallel to the interface. As a result, this creates a large velocity gradient and enhances the drag compared to a clean interface (no surfactant). By looking into the first two terms in (3.11), we can infer that 
$({3\mu +2})/{3(\mu +1)}$ and 
${\mu }/{(\mu +1)}$ are less than unity. Hence the drag experienced by the droplet in the large 
$Ma_\varGamma$ limit is more compared to that of a droplet with a clean interface.
The first-order drag 
$\boldsymbol {D}^H_1$ and droplet velocity 
$\boldsymbol {U}_1$ are given by
\begin{equation} \boldsymbol{D}_1^H={-}6{\rm \pi}\boldsymbol{U}_1+2{\rm \pi}\boldsymbol{Z}_1,\quad \text{and} \quad \boldsymbol{U}_1=\tfrac{1}{3}\boldsymbol{Z}_{1}+\boldsymbol{U}_{B,1}, \end{equation}
where the complete expression for 
$\boldsymbol {Z}_1$ is given in (A13).
4. Results and discussion
In the preceding section, we covered the general approach for solving an arbitrary Stokes flow. In this section, we investigate a surfactant-laden droplet in a Poiseuille flow as a special case (see figure 3). The Poiseuille flow 
$\boldsymbol {v}_\infty =(0,0,V_z)$ experienced by the droplet is described in the form (Pak et al. Reference Pak, Feng and Stone2014)
\begin{equation} {{V_z}}=1-\left(\frac{r }{R_0}\right)^2\sin^2{\theta}-\left(\frac{b }{R_0}\right)^2-\frac{2rb }{R_0^2}\sin{\theta}\cos{\phi}. \end{equation}
Here, the parameter 
$b$ represents the distance between the centre of the droplet and the centreline of the flow, and 
$R_0$ is the radius of the pipe. For the chosen velocity field (4.1), the associated scalars that describe the far-field conditions are
$$\begin{gather} A_\infty =\Bigg[ \frac{r }{2}\Bigg(1-\frac{b^2}{R_0^2}\Bigg)-\frac{r^3 }{5R_0^2}\Bigg]P^0_1(\cos{\theta}) +\Bigg(\frac{r^3}{30R_0^2}\Bigg)P_3^0(\cos{\theta})+\Bigg(\frac{b r^2 }{9R_0^2}\Bigg)P_2^1(\cos{\theta})\cos{\phi}, \end{gather}$$
$$\begin{gather}B_\infty =\Bigg(\frac{b r}{R_0^2}\Bigg)P_1^1(\cos{\theta})\sin{\phi}. \end{gather}$$
Figure 3. Schematic of a Poiseuille flow (as a special case of the Stokes flow) in a cylinder of radius 
$R_0$ in which a surfactant-laden droplet is released at a distance 
$b$ from the axis. Thermal gradient and gravity in the streamwise direction. Other notation is similar to figure 1.
4.1. Temperature distribution
In figure 4, we have shown temperature isosurfaces on the droplet at 
$r=1$. The meridional coordinate 
$\theta$ is measured along the droplet surface starting from the downstream (
$\theta =0$) to the upstream (
$\theta ={\rm \pi}$) location. In the absence of a far-field temperature gradient (
$G=0$), we observe that a centred monopole generates a uniform temperature on the droplet surface (see figure 4a). For the singularity strength 
$S=2$, we observe the uniform temperature equal to a value 
$2$. With the presence of ambient gradient (
$G=1$), the surface temperature varies, and it is different from a uniform value (see figure 4b). These variations occur primarily due to the surrounding thermal gradient and not due to the internal thermal singularity.
Figure 4. Temperature (
$T$) contours on the droplet surface for two ambient thermal gradients 
$G=0$ and 1. Droplets with (a–d) monopole-based, and ( f–i) dipole-based thermal singularities. On the droplet surface, the 
$\theta$ direction begins at east (E) and ends at north (N) in the clockwise direction covering south (S) and west (W) points. The droplet is subjected to ambient flow experienced as downstream at 
$\theta =0$ and upstream at 
$\theta ={\rm \pi}$. (e,j) Temperature profiles in the 
$\theta$ direction for fixed 
$\phi =0$. Continuous lines for 
$G=0$, and triangles for 
$G=1$. Red for centred (C) monopole, grey for off-centred (OC) monopole, orange for dipole oriented in the streamwise direction, and blue for dipole oriented in the cross-stream direction. Parameters are 
$S=2$, 
$\kappa =1$. Note that the temperature profiles are independent of other parameters such as droplet offset location from the axis 
$b$, 
$Pe_s$, 
$\mu$, 
$Ma_T$, 
$Ma_{\varGamma }$ and others.
In the case of an off-centred monopole (with 
$r=0.5$, 
$\theta =0$, 
$\phi =0$), the temperature is maximum at 
$\theta =0$ and minimum at 
$\theta ={\rm \pi}$ (see figures 4c,d). A similar behaviour is also seen for a dipole oriented in the streamwise direction (see figures 4f,g). However, in the case of a dipole oriented perpendicular to the flow, the maximum surface temperature is observed at 
$\theta ={\rm \pi} /2$ for 
$G=0$ (see figure 4h), and in a direction different from the streamwise and transverse directions for 
$G\neq 0$ (see figure 4i). We have quantified 
$T(r=1,\theta,\phi =0)$ for the presence of a monopole and a dipole in figures 4(e) and 4(j), respectively. Note that the temperature variation on the droplet surface results in additional thermal Marangoni stress due to singularities, and thus affects the droplet migration.
4.2. In the small-$Pe_s$ limit $(Pe_s\ll 1)$
For 
$Pe_s\ll 1$, we have computed surfactant concentration (using (A4) and (A8) in Appendix A), the drag and the droplet migration velocity (see Appendix B). Our generalized calculations match with Pak et al. (Reference Pak, Feng and Stone2014) in the absence of ambient thermal field (
$G=0$) and internal thermal singularities. Even in the presence of thermocapillary effects (without any singularities), our results agree with Das et al. (Reference Das, Mandal and Chakraborty2018) and Sharanya et al. (Reference Sharanya, Raja Sekhar and Rohde2019). In the following two subsubsections, we discuss both surfactant distribution and droplet velocity in the presence of thermal singularities.
4.2.1. Surfactant distribution
The surfactant distribution 
$\varGamma (\theta,\phi )$ for a droplet released on a Poiseuille flow axis (
$b=0$) is shown in figure 5 for different thermal singularities. Here, the ambient thermal field is acting along the imposed Poiseuille flow. The surfactant distribution is relatively high at 
$\theta ={\rm \pi}$ (west) compared to 
$\theta =0$ (east), and the observed is due to the Marangoni stress-induced flow on the droplet surface. Interestingly, the east–west located surfactant distribution moves towards the north–south directions when a dipole is oriented in the cross-stream directions, so we observe surfactant accumulation at the south-west region of the droplet. A similar result, but an increase in non-uniform surfactant distribution, is found when the droplet is released at an offset from the flow centreline (
$b\neq 0$).
Figure 5. Surfactant concentration 
$\varGamma$ on the droplet surface in the presence of (a,e) a centred monopole, (b, f) an off-centred monopole, (c,g) a dipole oriented in the streamwise direction, and (d,h) a dipole oriented in the cross-stream direction. Droplet released on flow centreline/axis (
$b=0$) from (a) to (d), and released at an offset from the flow axis (
$b=1$) from (e) to (h). Colour bar for 
$\varGamma$ in a linear range from 0.99 (blue) to 1.01 (red). Parameters are 
$S=2$, 
$Pe_s=0.1$, 
$G=1$, 
$Ma_{T}=1$, 
$Ma_{\varGamma }=10$, 
$\kappa =1$, 
$\mu =1$, 
$R_0=5$ and 
$Bo=0$.
4.2.2. Droplet migration velocity
In figure 6(a), we have shown the migration speed 
$|U_z|$ in the streamwise direction as a function of singularity strength 
$S$. In the absence of thermal singularity (
$S=0$), we found 
$|{U}_z| \approx 1.05$, which agrees with the results of Das et al. (Reference Das, Mandal and Chakraborty2018) for the chosen parameters. A centred monopole does not affect the migration speed (refer to (B2)), which matches the results of Rednikov & Ryazantsev (Reference Rednikov and Ryazantsev1989) for homogeneous internal heating of a droplet. For 
$S>0$, the speed increases above 
$1.05$ for both the off-centred monopole and the dipole. In the inset of figure 6(a), we have shown 
${U}_z/|{U}_z (S=0)|$. The droplet with the dipole orientation 
$\theta =0$ gives maximum speed compared to any other orientation. In figure 6(b), droplet cross-stream speed 
$|U_x|$ as a function of 
$S$ is shown. Migration velocity increases with increasing 
$S$. Interestingly, we observe that a dipole can reverse the direction of the droplet migration when 
$S$ is greater than a critical value 
$S_c$. In the inset of figure 6(b), we have shown the droplet's relative cross-stream velocity 
$U_x/|U_x(S=0)|$. For small values of 
$S$, the relative speed is negative, which indicates that the drop migrates towards the flow centreline. With an increase in 
$S$, the droplet can migrate away from the centreline. At critical strength 
$S_c$, we find 
$U_x=0$. It is worth noting that 
$S_c$ is a function of the parameters 
$\mu$, 
$\kappa$, 
$Ma_T$, 
$Ma_\varGamma$ and dipole orientation angle 
$\theta$. Figure 7 shows 
$S_c$ as a function of the dipole orientation 
$\theta$. We observe that 
$S_c$ decreases by nearly two orders of magnitude when orientation changes from 
$\theta =0$ to 
${\rm \pi} /2$ with a small variation in 
$Ma_T$. This suggests that when the dipole is oriented in the cross-stream direction (
$\theta ={\rm \pi} /2$), a small 
$S$ is sufficient to perturb the droplet motion from the centreline to away from it. In the insets of figure 7, we have shown 
$S_c$ dependence on 
$\mu$ (figure 7a) and 
$\kappa$ (figure 7b). We observe that 
$S_c$ decreases with both 
$\mu$ and 
$\kappa$; however, we find that 
$\kappa$ dependence corresponds to a significant change in 
$S_c$.
Figure 6. Droplet velocity (in magnitude) as a function of thermal singularity strength 
$S$: (a) streamwise direction, and (b) cross-stream direction. Thick dashed line for the absence of thermal singularity, centred monopole (red dashes) and off-centred monopole (grey line). Dipole with different orientations in continuous lines: 
$\theta =0$ (orange), 
${\rm \pi} /6$ (violet), 
${\rm \pi} /4$ (green), 
${\rm \pi} /3$ (cyan) and 
${\rm \pi} /2$ (blue). Insets for droplet velocity re-scaled (or compared) by the case with no thermal singularity. Note that the droplet moves away from the centreline (or axis) beyond critical singularity strength 
$S_c$ (defined for a dipole). Parameters are 
$Pe_s=0.1$, 
$G=1$, 
$b=1$, 
$Ma_{T}=1$, 
$Ma_{\varGamma }=10$, 
$\kappa =1$, 
$\mu =1$, 
$R_0=5$ and 
$Bo=0$. Here, 
$S=0$ corresponds to no thermal singularity (Das et al. Reference Das, Mandal and Chakraborty2018).
Figure 7. Critical dipole strength 
$S_c$ as a function of dipole orientation angle 
$\theta$ for 
$Ma_T=0.1$ (red), 
$Ma_T=0.3$ (green), and 
$Ma_T=1$ (blue). Parameters are 
$Pe_s=0.1$, 
$G=1$, 
$Ma_{\varGamma }=10$, 
$\kappa =1$, 
$\mu =1$, 
$R_0=5$ and 
$Bo=0$. Insets for different (a) viscosity and (b) thermal conductivity.
4.3. Small-$Pe_s$ versus large-$Pe_s$ limits
The migration velocity of a surfactant-laden droplet (subjected to a streamwise oriented dipole) in the large-
$Pe_s$ limit is given by
\begin{equation} \boldsymbol{U}=\left[1-\frac{2}{3R_0^2}-\frac{b^2}{R_0^2}+\frac{2\omega}{9}\,\frac{Bo}{Ca}\right]\hat{k}- Ma_\varGamma^{{-}1}\left[\frac{(G+\kappa S)b\,Ma_T}{(\kappa+2)R_0^2}+\frac{2b}{3R_0^4}+\frac{\omega b}{9R_0^2}\,\frac{Bo}{Ca}\right]\hat{i}+O(Ma_\varGamma^{{-}1}). \end{equation} In figure 8, we have presented the results for both small- and large-
$Pe_s$ limits. Notably, we can observe that the magnitude of the cross-stream migration velocity is approximately one order higher in the large-
$Pe_s$ limit compared to the small-
$Pe_s$ limit. The behaviour remains consistent for other thermal singularities as well.
Figure 8. Cross-stream velocity (in magnitude) as a function of singularity strength 
$S$ under small-
$Pe_s$ (red) and large-
$Pe_s$ (blue) limits. Parameters are 
$Pe_s=0.1$, 
$G=1$, 
$Ma_{T}=1$, 
$Ma_{\varGamma }=10$, 
$\kappa =1$, 
$\mu =1$, 
$R_0=5$ and 
$Bo=0$.
4.4. A discussion on critical Marangoni number set by thermal singularities for small $Pe_s$
In the earlier sections, we have constructed a generalized result on surfactant-laden droplet motion and its control using thermal singularities. Here, we present an overview of the role played by thermocapillary action and surfactants on droplet motion.
In the case of droplet-based microfluidics, surfactants are added to the continuous or dispersed phase to maintain the stability of the droplet, prevent coalescence and avoid any cross-contamination. However, such addition of surfactant leads to a reduction in droplet speed at first order in the streamwise direction without any cross-stream migration, and a cross-stream migration towards the Poiseuille flow centreline is seen as a second-order effect in 
$Pe_s$ (Pak et al. Reference Pak, Feng and Stone2014). The stabilized droplet that reaches the centreline may control the overall production rate in micro-devices. Note that this technique, with precise control, can be highly advantageous in cell processing applications, optimizing bioreactor efficiency, etc. Now, if we can disperse the droplets from the centreline and at the same time maintain stability, then we can account for more droplets in the device, thus improving the overall production rate. Similarly, an imposed thermal gradient (without any surfactants) results in a leading-order effect in droplet streamwise motion either by speeding up (if the thermal gradient aligns with the flow) or slowing down (if the thermal gradient opposes the flow) without any cross-stream movement. Although manipulating droplet speed is possible, this approach does not achieve a high droplet production rate within a limited device width. An alternative control mechanism that manipulates droplet migration efficiency would be advantageous. In this regard, a thermal singularity helps to control the droplet motion with leading-, first- and second-order effects in the streamwise and cross-streamwise directions. For example, a dipole can help the droplet to speed up or slow down, depending on its orientation.
In this context, we seek whether a critical Marangoni number (
$Ma^{*}_{T,z}$ or 
$Ma^{*}_{\varGamma,z}$) exists at which a droplet halts (
${U_z}={0}$) or reverses its motion. Our focus in this subsection is solely on the streamwise migration velocity of the droplet. The reason behind this emphasis is that it is more feasible to maintain a thermal gradient in the streamwise direction due to the device configurations, instead of in the cross-stream direction. Furthermore, we are interested in exploring the interfacial tension effects, thus neglecting gravity (
$Bo=0$) for the discussion. In the following, we show four cases: (1) when thermal gradient alone is present, without surfactant or the singularities (
$Ma_T\neq 0$, 
$Ma_{\varGamma }=0$ and 
$S=0$); (2) when surfactant alone is present, without thermal gradient or the singularities (
$Ma_T= 0$, 
$Ma_{\varGamma }\neq 0$ and 
$S=0$); (3) when both thermal gradients and surfactant are present, but thermal singularity is absent (
$Ma_T\neq 0$, 
$Ma_{\varGamma }\neq 0$ and 
$S=0$); and (4) when thermal gradients, surfactant and the singularities are present (
$Ma_T\neq 0$, 
$Ma_{\varGamma }\neq 0$ and 
$S\neq 0$). Note that in all the cases, the applied thermal gradient is opposite to the flow direction (i.e. 
$G=-1$) with a specific aim of halting a droplet's motion.
4.4.1. Case 1: $Ma_T\neq 0$, $Ma_{\varGamma }=0$ and $S=0$
We have equated 
${U}_z=0$ (refer to (B2)) and found the critical thermal Marangoni number as
\begin{equation} Ma^*_{T,z}={-}\frac{(3\mu+2)(\kappa+2)}{2G}\Bigg[1-\left( \frac{b}{R_0}\right)^2 \Bigg]+\frac{\mu(\kappa+2)}{2 G R_0^2}. \end{equation} In figure 9(a), we have shown 
$Ma^*_{T,z}$ as a function of droplet offset distance 
$b$ for different values of 
$\mu$. A droplet away from the centreline requires small 
$Ma^*_{T,z}$ to cease its motion. For highly viscous droplets, a higher thermal gradient is needed to cease the droplet. The high thermal conductivity of the droplet corresponds to low thermocapillary stresses on the interfaces, which results in a large 
$Ma^*_{T,z}$.
Figure 9. Critical Marangoni number required to halt droplet released at an offset distance 
$b$ in a Poiseuille flow: (a) for case1, and (b) for case 3. Continuous lines for 
$\kappa =1$ and dashed lines for 
$\kappa =2$. Other parameters are 
$G=-1$ and 
$R_0=5$.
4.4.2. Case 2: $Ma_T = 0$, $Ma_{\varGamma } \neq 0$ and $S=0$
The critical surfactant Marangoni number 
$Ma^*_{\varGamma,z}$ can be found from the equation
\begin{align} 1-\left(\frac{b}{R_0}\right)^2-\frac{2\mu}{(3\mu+2)R_0^2}-\frac{4\,Pe_s}{3(3\mu+2)^2 R_0^2}\,Ma^*_{\varGamma,z}+\frac{4\,Pe^2_s}{3(3\mu+2)^3 R_0^2}(Ma^*_{\varGamma,z})^2=0. \end{align}
Note that this is a quadratic equation in terms of 
$Ma^*_{\varGamma,z}$. In order to obtain a real non-zero solution, the discriminant should be positive. Consequently, simplifying the discriminant of (4.5) yields the condition 
$0 < R_0^2 - b^2 < ({6\mu +1})/({3\mu +2})$. Thus we can infer that unless the droplet is released very close to the wall (i.e. 
$b \approx R_0$), surfactants can indeed reduce the streamwise migration velocity of the droplet, but they cannot bring it to halt completely.
4.4.3. Case 3: $Ma_T \neq 0$, $Ma_{\varGamma } \neq 0$ and $S=0$
We found the critical thermal Marangoni number as (using (B2))
\begin{align}
Ma^*_{T,z}=\frac{({\kappa}+2) [3 (b^2-R^2_0) (3
{\mu}+2)^3+6 {\mu} (3 {\mu}+2)^2+4\, {Ma_\varGamma}\,(3
{\mu}+2)\,{Pe_s}-4\,Ma_\varGamma^2\,Pe_s^2]}{6 G R_0^2
[(3 {\mu}+2)^2-{Ma_\varGamma}\,(3
{\mu}+2)\,{Pe_s}+Ma_\varGamma^2\,Pe_s^2]}.
\end{align} In figure 9(b), we have shown 
$Ma^*_{T,z}$ as a function of 
$b$. As mentioned earlier, the presence of surfactant reduces the droplet speed, and a tiny opposing thermal gradient is enough to halt the droplet completely (compare the 
$\mu =1$ line in figures 9a,b), thereby reducing the critical 
$Ma_T$ required.
In figure 10, we have shown the variation of 
$Ma^*_{T,z}$ as a function of 
$Ma_\varGamma$. For a droplet released either on the centreline or near the walls in a Poiseuille flow, surprisingly, our result reveals that we need larger thermal gradients (or large critical 
$Ma^*_{T,z}$) for 
$Ma_{\varGamma }<20$. Note that (4.6) has two asymptotic limits in terms of small and large values of 
$Ma_{\varGamma }Pe_s$. For smaller values, we have
\begin{align} \left. Ma^*_{T,z}\right|_{Ma_\varGamma Pe_s \to 0}\approx \Bigg({\frac{\kappa+2}{6GR^2_0}}\Bigg)[6\mu+3(3\mu+2)(b^2-R^2_0)+\{2+3(b^2-R^2_0)\}\,Ma_{\varGamma}\,Pe_s], \end{align}
which is a linearly increasing function of 
$Ma_{\varGamma }\,Pe_s$. In the case of large 
$Ma_{\varGamma }Pe_s$, we find
\begin{align} \left.
Ma^*_{T,z}\right|_{{Ma_\varGamma Pe_s} \to \infty}
\approx{-}
\frac{4(\kappa+2)}{6GR_0^2(3\mu+2)^2}\left[1-\frac{1}{(Ma_\varGamma\,Pe_s)^{2}}\left(1+(3\mu+2)^3
\left(\frac{1}{2}+\frac{3(b^2-R_0^2)}{4}\right)\right)
\right].
\end{align}
Both the asymptotic limits are shown in figure 10. Interestingly, at leading order, the small 
$Ma_{\varGamma }Pe_s$ limit depends on 
$b$, whereas at the large limit, it does not. Note that the large limit is non-zero and has a finite positive value close to zero. The peak values observed on critical 
$Ma^*_{T,z}$ curves near 
$Ma_{\varGamma }\approx 20$ decrease with increasing 
$b$.
Figure 10. Critical thermal Marangoni number 
$Ma^*_{T,z}$ as a function of 
$Ma_\varGamma$ for case 3. Low 
$Ma_{\varGamma }$ follows (4.7) (thin dashed lines at different 
$b$), and large 
$Ma_{\varGamma }$ follows (4.8) (thick dashed line). Parameters are 
$\mu =1$, 
$\kappa =1$, 
$Pe_s=0.1$, 
$G=-1$, 
$R_0=5$ and 
$Bo=0$.
4.4.4. Case 4: $Ma_T \neq 0$, $Ma_{\varGamma } \neq 0$ and $S \neq 0$
In the case of thermal singularities, we have generalized 
$Ma^*_{T,z}$ as a function of 
$Ma_{\varGamma }$ and 
$S$ using (B6) for a dipole oriented along (
$\theta =0$) or opposite to (
$\theta ={\rm \pi}$) the streamwise direction. The corresponding equations for 
$Ma^*_{T,z}$ are given by
\begin{align} Ma^*_{T,z}=\frac{({\kappa}+2) [3 (b^2-R_0^2) (3 {\mu}+2)^3+6 {\mu} (3 {\mu}+2)^2+4\, {Ma_\varGamma}\,(3 {\mu}+2)\,{Pe_s}-4\,Ma_\varGamma^2\,Pe_s^2]}{6 (G+\kappa S) R_0^2 [(3 {\mu}+2)^2-{Ma_\varGamma}\,(3 {\mu}+2)\,{Pe_s}+Ma_\varGamma^2\,Pe_s^2]} \end{align}and
\begin{align} Ma^*_{T,z}=\frac{({\kappa}+2) [3 (b^2-R_0^2) (3 {\mu}+2)^3+6 {\mu} (3 {\mu}+2)^2+4\, {Ma_\varGamma}\,(3 {\mu}+2)\,{Pe_s}-4\,Ma_\varGamma^2\,Pe_s^2]}{6 (G-\kappa S) R_0^2 [(3 {\mu}+2)^2-{Ma_\varGamma}\,(3 {\mu}+2)\,{Pe_s}+Ma_\varGamma^2\,Pe_s^2]}, \end{align}
respectively. Note that the mentioned equations are similar to (4.6), except that the term 
$G$ in the denominator is replaced by 
$G\pm \kappa S$. For a dipole oriented in the flow direction, the required critical value can reach infinity if 
$G+\kappa S\to 0$ (since 
$G=-1$). In contrast, for a dipole oriented opposite to the flow, the denominator (in magnitude) increases with an increase in 
$S$, thus the required critical Marangoni number decreases (see figure 11). Note that the functional form and the asymptotic behaviour in case 4 remain the same, as mentioned in case 3.
Figure 11. Critical thermal Marangoni number 
$Ma^*_{T,z}$ as a function of 
$Ma_\varGamma$ for a dipole with different strengths 
$S$ oriented opposite to the streamwise direction (case 4). Other parameters are 
$\mu =1$, 
$\kappa =1$, 
$Pe_s=0.1$, 
$G=-1$, 
$R_0=5$ and 
$Bo=0$.
4.5. Dependency on viscosity ratio in the small-$Pe_s$ limit
Analysing droplet migration in terms of viscosity ratio in microfluidics is crucial for understanding and optimizing microfluidic systems (Baumgartner, Brenn & Planchette Reference Baumgartner, Brenn and Planchette2020). In this context, we have shown 
$U_z$ and 
$U_x$ as functions of 
$\mu$ (refer to (B6) in Appendix B) in figure 12. For a fixed 
$\mu$, 
$U_z$ for a centred monopole and dipole (oriented perpendicular to the flow) coincides with the case without any thermal singularity, and matches the results of Das et al. (Reference Das, Mandal and Chakraborty2018) for the chosen parameter values. Both the off-centred monopole and streamwise oriented dipole increase the migration velocity. For a small viscosity ratio (
$\mu =0.01$), 
$U_z$ is increased by approximately 
$10\,\%$ for an off-centred monopole, and by 
$25\,\%$ for a streamwise oriented dipole when compared to no thermal singularity (Das et al. Reference Das, Mandal and Chakraborty2018). However, a cross-stream oriented dipole does not affect 
$U_z$. For large 
$\mu$, the droplet behaves as a rigid sphere, and the corresponding influence of the Marangoni stress becomes insignificant, resulting in a decrease in 
$U_z$. In figure 12(b), we have shown 
$U_x$ as a function of 
$\mu$. We can observe that the cross-stream migration velocity is influenced significantly by a transversely oriented (
$\theta ={\rm \pi} /2$) dipole in comparison to other singularities.
Figure 12. Droplet migration velocity (a) streamwise (
${U}_z$) and (b) cross-stream (
${U}_x$) versus viscosity ratio 
$\mu$ for 
$S=0$ (dashed line), centred monopole (dotted line), off-centred monopole (grey), dipole oriented streamwise (orange) and dipole oriented transversely (blue). The ambient thermal field acts along the imposed Poiseuille flow, i.e. 
$G=1$. Other parameters are 
$S=2$, 
$Ma_T=1$, 
$Ma_\varGamma =10$, 
$Pe_s=0.1$, 
$\kappa =1$, 
$b=1$, 
$R_0=5$ and 
$Bo=0$.
It is worth mentioning that when considering a centred/off-centred monopole and a streamwise oriented dipole, the cross-stream migration velocity appears at 
$O(Pe_s^2)$. However, in the case of a transversely oriented dipole, the cross-stream migration velocity can be observed from 
$O(1)$ onwards. Further, in this particular configuration (droplet with transversely oriented dipole), the leading-order cross-stream migration is due mainly to the thermocapillary effect (refer to (B6)). From 
$O(Pe_s)$ onwards, the effect of surfactant can be seen (refer to figure 13a), which is different from the earlier 
$O(Pe_s^2)$ corrections (Das et al. Reference Das, Mandal and Chakraborty2018). This is due to the coupling between the thermal and the surfactant-induced Marangoni stresses. Even for a very small singularity strength (
$S=2 \times 10^{-4}$), the leading-order effect can be observed very clearly. In figure 13(b), we classify the regimes with reference to the migration of the droplet. Accordingly, for a specific 
$(Ma_\varGamma, Ma_T)$ pair, one can identify whether the droplet migrates away from or towards the centreline when the other parameters are fixed.
Figure 13. (a) Droplet's cross-stream migration velocity 
${U}_x$ versus viscosity ratio 
$\mu$ for a transversely oriented dipole. Migration velocity for different orders: leading order (triangles), first order (diamonds), second order (circles) and combined (continuous line). Here, 
$Ma_T=1$ and 
$Ma_\varGamma =10$. (b) Regime diagram showing two different patterns of droplet cross-stream motion for 
$\mu =1$. Other parameters are 
$S=2 \times 10^{-4}$, 
$G=1$, 
$Pe_s=0.1$, 
$\kappa =1$, 
$b=1$, 
$R_0=5$ and 
$Bo=0$.
5. Summary and conclusions
In this paper, we have analysed the migration of a droplet in an arbitrary Stokes flow under the influence of thermocapillary and surfactant-induced Marangoni stress for different types of internal thermal singularities. Neglecting the effects of thermal convection, fluid inertia and shape deformation, we have obtained analytical expressions for the temperature field, surfactant concentration and migration velocity using perturbation analysis for the two limiting cases: (i) 
$Pe_s \ll 1$, in which surface diffusion dominates the surfactant transport; and (ii) 
$Pe_s \gg 1$, in which surface convection dominates the surfactant transport. The migration velocity of a surfactant-laden droplet in an ambient Poiseuille flow (as a particular case of an arbitrary Stokes flow) depends strongly upon the type of internal thermal singularity considered. The migration velocity is unaffected by a localized centred monopole heat source, similar to the findings for a droplet subjected to homogeneous internal heating by Rednikov & Ryazantsev (Reference Rednikov and Ryazantsev1989). However, an off-centred monopole or a dipole heat source influences the overall migration of the droplet.
Our analytical results show that in the small-
$Pe_s$ limit, a drop with an off-centred monopole or a streamwise oriented (
$\theta = 0$) dipole experiences cross-stream migration at 
$O(Pe_s^2)$. The effect becomes prominent in the case of an oriented dipole (
$\theta \neq 0$). In this scenario, we observe the cross-stream migration at the leading order, which is due purely to the thermocapillary effects generated by the asymmetric temperature distribution on the droplet surface. Due to the coupling between thermocapillary-induced Marangoni stress and surfactant-induced Marangoni stress, the combined effect is observed at 
$O(Pe_s)$ onwards, which is different from the classical 
$O(Pe_s^2)$ findings by Das et al. (Reference Das, Mandal and Chakraborty2018). Furthermore, when the droplet is positioned at the flow centreline, it can also move across the stream with this configuration. As we increase the strength of thermal singularity, the migration velocity of the droplet gets enhanced. Remarkably, in the presence of an oriented dipole, a critical strength 
$S_c$ exists that determines whether the droplet migrates towards or away from the flow centreline. Further, in the absence of a thermal gradient (due to both ambient temperature as well as thermal singularity), surfactants can reduce the droplet migration velocity, but they are unable to halt it completely unless the droplet is close enough to the wall. However, when both the thermal gradient and surfactants are present, we can easily bring the droplet to a complete halt by applying a dipole oriented opposite to the flow. In comparison to the small-
$Pe_s$ limit, the magnitude of the cross-stream migration velocity is approximately one order higher for large 
$Pe_s$.
In summary, we have developed a general methodology for addressing thermal singularities in droplet dynamics, and provided Faxen's laws for drag and migration velocity. Our research findings indicate that thermal singularities have the potential to control both the streamwise and cross-stream migration velocities of surfactant-laden droplets. The analytical results obtained in this study serve as valuable benchmarks for future computational methods pertaining to droplet dynamics. For future scope, our findings provide essential insights for estimating drag and energy requirements in microbe-encapsulated droplet technologies, thus opening up promising avenues for further exploration in this field.
Acknowledgements
The authors would like to acknowledge Professor S. Chakraborty, Department of Mechanical Engineering, Indian Institute of Technology (IIT) Kharagpur, for his insightful suggestions.
Funding
A.B. gratefully acknowledges the research support from the Council of Scientific & Industrial Research (CSIR), Government of India. R.L. acknowledges research grants from the Board of Research and Nuclear Sciences (BRNS) and the Mathematical Research Impact Centric Support (MATRICS) scheme sponsored by the Government of India to study multiphase flows. We sincerely thank P. Shakti for computing resources developed under National Super Computing Mission at IIT Kharagpur, India.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Drag and migration velocity in the small- and large-$Pe_s$ limits
We have considered a Cartesian coordinate system 
$(x,y,z)$ having its origin at the droplet centre, as shown in figure 1. We are interested in calculating the axial and cross-stream migration velocity in this frame of reference (Subramanian & Balasubramaniam Reference Subramanian and Balasubramaniam2001).
A.1. For $Pe_s \ll 1$
The leading-order drag, as given in (3.11), is
\begin{align} { \boldsymbol{{D}}^H_0 }&= 6{\rm \pi}\left({\frac{3\mu +2 }{3(\mu+1)}\,[\boldsymbol{u}_{\infty,0}]_0+\frac{\mu }{6(\mu+1)}\,[\nabla^2\boldsymbol{u}_{\infty,0}]_0}+{\frac{2\,Ma_T }{{3(\mu+1)(\kappa+2)}}\,[\boldsymbol{\nabla} T_\infty]_0}\right.\nonumber\\ &\quad \left.{}+{\frac{{2\kappa \mathscr{b}_1\,Ma_T }}{{3(\mu+1)(\kappa+2)}}\,\boldsymbol{Y}_{0}}\right). \end{align}The aforementioned compact forms are due to the relationships
\begin{equation} \left. \begin{gathered} {[\boldsymbol{u}_{\infty,\boldsymbol{0}}]}_0 = 2\alpha_1^{\infty (0)} (C^{(0)}_{11}\hat{i}+D^{(0)}_{11}\hat{j}+C^{(0)}_{10}\hat{k}),\\ [\nabla^2 \boldsymbol{u}_{\infty,\boldsymbol{0}}]_0 = 20\hat\alpha^{\infty (0)}_1({C}^{(0)}_{11}\hat{i}+{D}^{(0)}_{11}\hat{j}+C^{(0)}_{10}\hat{k}),\\ {[\boldsymbol{\nabla} T_\infty]}_0 = \mathscr{c}_1({A}_{11}\hat{i}+{B}_{11}\hat{j}+{A}_{10}\hat{k}). \end{gathered} \right\} \end{equation}Note that the vector
\begin{equation} \boldsymbol{Y}_{0}={A}_{11}\hat{i}+{B}_{11}\hat{j}+{A}_{10}\hat{k} \end{equation}
is related to the spherical harmonic 
$\mathcal {S}_n(\theta,\phi )=\sum _{m=0}^{n} P_n^m(\cos {\theta })\,({A}_{nm}\cos {m\phi }+{B}_{nm}$ 
$\sin {m\phi })$.
Now, using leading-order solutions, we solve the first-order surfactant transport equation given by (3.10b). We express the surfactant concentration 
$\varGamma _1$ in terms of surface spherical harmonics as
\begin{equation} \varGamma_1=\sum_{n=1}^\infty \sum_{m=0}^n \left(I_{nm}^{(1)} \cos m\phi+J_{nm}^{(1)} \sin m \phi\right) P_n^m (\cos \theta), \end{equation}
where 
$I_{nm}^{(1)}$ and 
$J_{nm}^{(1)}$ are
$$\begin{gather} I_{nm}^{(1)} ={-}\frac{(2n+1)(n-m)!}{2(n+m)!}\,\frac{1}{{\rm \pi} n(n+1)}\int_{\phi=0}^{2{\rm \pi}}\int_{\theta=0}^{\rm \pi}(\boldsymbol{\nabla}_s\boldsymbol{\cdot}\boldsymbol{u}_{s,0})\,P_n^m (\cos{\theta}) \cos{m\phi}\sin{\theta} \,{\rm d}\theta \,{\rm d}\phi, \end{gather}$$
$$\begin{gather}J_{nm}^{(1)} ={-}\frac{(2n+1)(n-m)!}{2(n+m)!}\,\frac{1}{{\rm \pi} n(n+1)}\int_{\phi=0}^{2{\rm \pi}}\int_{\theta=0}^{\rm \pi}(\boldsymbol{\nabla}_s\boldsymbol{\cdot}\boldsymbol{u}_{s,0})\,P_n^m (\cos{\theta}) \sin{m\phi}\sin{\theta} \,{\rm d}\theta \,{\rm d}\phi, \end{gather}$$which implies
$$\begin{gather} I_{nm}^{(1)} =\left(\frac{2n-1}{\mu+1}\,\alpha_n^{\infty(0)}+\frac{2n+3}{\mu+1}\,\hat\alpha_n^{\infty(0)}\right)C^{(0)}_{nm}-\frac{(\mathscr{c}_n+\kappa \mathscr{b}_n)\,Ma_T}{(\mu+1)(n\kappa+n+1)}\,A_{nm}, \end{gather}$$
$$\begin{gather}J_{nm}^{(1)} =\left(\frac{2n-1}{\mu+1}\,\alpha_n^{\infty(0)}+\frac{2n+3}{\mu+1}\,\hat\alpha_n^{\infty(0)}\right)D^{(0)}_{nm}-\frac{(\mathscr{c}_n+\kappa \mathscr{b}_n)\,Ma_T}{(\mu+1)(n\kappa+n+1)}\,B_{nm}. \end{gather}$$Similarly, the first-order hydrodynamic drag and corresponding migration velocity are given by (3.13a) and (3.13b), where
\begin{equation} \boldsymbol{Y}_1=I_{11}^{(1)} \hat{i}+J_{11}^{(1)}\hat{j}+I_{10}^{(1)}\hat{k}. \end{equation}
Now, let us consider 
$\varGamma _2$ as
\begin{equation} \varGamma_2=\sum_{n=1}^\infty \sum_{m=0}^n (I_{nm}^{(2)} \cos m\phi+J_{nm}^{(2)} \sin m \phi) P_n^m (\cos \theta). \end{equation}
Correspondingly, 
$I_{nm}^{(2)}$ and 
$J_{nm}^{(2)}$ are given by
$$\begin{align} I_{nm}^{(2)}
&={-}\frac{(2n+1)(n-m)!}{2(n+m)!}\,\frac{1}{{\rm \pi}
n(n+1)}\nonumber\\
&\quad\times \int_{\phi=0}^{2{\rm \pi}}\int_{\theta=0}^{\rm \pi}(\boldsymbol{\nabla}_s\boldsymbol{\cdot}(\boldsymbol{u}_{s,0}
\varGamma_1+\boldsymbol{u}_{s,1} \varGamma_0))P_n^m(\cos{\theta}) \cos{m\phi}\sin{\theta} \,{\rm
d}\theta \,{\rm d}\phi,
\end{align}$$
$$\begin{align} J_{nm}^{(2)} &={-}\frac{(2n+1)(n-m)!}{2(n+m)!}\,\frac{1}{{\rm \pi}
n(n+1)}\nonumber\\
&\quad \times \int_{\phi=0}^{2{\rm \pi}}\int_{\theta=0}^{\rm \pi}(\boldsymbol{\nabla}_s\boldsymbol{\cdot}(\boldsymbol{u}_{s,0}
\varGamma_1+\boldsymbol{u}_{s,1} \varGamma_0))P_n^m (\cos{\theta}) \sin{m\phi}\sin{\theta}\,{\rm
d}\theta \,{\rm d}\phi.
\end{align}$$
The second-order hydrodynamic drag and migration velocity are given by (3.14a) and (3.14b). The corresponding compact form is
\begin{equation} \boldsymbol{Y}_2=I_{11}^{(2)} \hat{i}+J_{11}^{(2)}\hat{j}+I_{10}^{(2)}\hat{k}. \end{equation}A.2. For $Pe_s\gg 1$
For large 
$Pe_s$, the first-order surfactant concentration is given by
\begin{equation} \varGamma_{1}=\sum_{n=1}^\infty\sum_{m=0}^n \varLambda_{n}( C_{nm}^{(0)}\cos{m\phi}+D_{nm}^{(0)}\sin{m\phi})P_n^m(\cos{\theta}), \end{equation}where
\begin{align} \varLambda_n=(2n-1)(2n+1)\alpha_n^{\infty{(0)}}+(2n+3)(2n+1)\hat\alpha_n^{\infty (0)}-\frac{(2n+1)(\mathscr{c}_n+\kappa \mathscr{b}_n)\,Ma_T}{(n\kappa+n+1)}. \end{align}Thus the first-order hydrodynamic drag and droplet migration velocity are given by (3.20a,b), where
\begin{equation} \boldsymbol{Z}_1=L_{111}\hat{i}+L_{112}\hat{j}+L_{101}\hat{k}, \end{equation}
and the terms 
$L_{nm1}$, 
$L_{nm2}$ are
$$\begin{gather} L_{nm1}={-}\frac{(2n+1)(n-m)!}{2{\rm \pi}(n+m)!}\int_{\phi=0}^{2{\rm \pi}}\int_{\theta=0}^{\rm \pi}(\boldsymbol{\nabla}_s\boldsymbol{\cdot}\boldsymbol{u}_{s,0})P_n^m (\cos{\theta}) \cos{m\phi}\sin{\theta} \,{\rm d}\theta \,{\rm d}\phi, \end{gather}$$
$$\begin{gather}L_{nm2} ={-}\frac{(2n+1)(n-m)!}{2{\rm \pi}(n+m)!}\int_{\phi=0}^{2{\rm \pi}}\int_{\theta=0}^{\rm \pi}(\boldsymbol{\nabla}_s\boldsymbol{\cdot}\boldsymbol{u}_{s,0})P_n^m (\cos{\theta}) \sin{m\phi}\sin{\theta} \,{\rm d}\theta \,{\rm d}\phi. \end{gather}$$Appendix B. Drag and migration velocity in Poiseuille flow for $Pe_s \ll 1$
Using the mentioned Poiseuille flow field (refer to (4.1)) and thermal conditions, we have computed the surfactant concentration 
$\varGamma$ for different thermal singularities using (A4) and (A8). For a centred monopole, the surfactant concentration is given by
\begin{align} \varGamma &=1+Pe_s
\left[I_{10}^{(1)}\cos{\theta}+I_{21}^{(1)}\cos{\phi}\,
P_2^1 (\cos{\theta})+ I_{30}^{(1)}\,P_3^0
(\cos{\theta})\right]\nonumber\\ &\quad + Pe_s^2
\left[I_{10}^{(2)}\cos{\theta}+I_{21}^{(2)}\cos{\phi}\,
P_2^1 (\cos{\theta}) +I_{30}^{(2)}\,P_3^0 (\cos{\theta})
+I_{11}^{(2)}\cos{\phi}\, P_1^1
(\cos{\theta})\right.\nonumber\\ &\left.\quad +\,I_{20}^{(2)}\,P_2^0 (\cos{\theta})+I_{22}^{(2)}
\cos(2\phi)\, P_2^2 (\cos{\theta})
+I_{31}^{(2)}\cos{\phi}\, P_3^1 (\cos{\theta})\right]
+O(Pe_s^3), \end{align}
where the quantities 
$I_{10}^{(1)}$, 
$I_{21}^{(1)}$, etc. are provided below. After evaluating (3.12), (3.13a), (3.13b), (3.14a) and (3.14b), we obtain the droplet migration velocity as
\begin{align} \boldsymbol{U} &=
\Bigg[1-\left(\frac{b}{R_0}\right)^2-\frac{2\mu}{(3\mu+2)R_0^2}+
\frac{2G\,Ma_T}{(\kappa+2)(3\mu+2)}+\frac{2\omega}{3}\left(\frac{1+\mu}{3\mu+2}\right)\frac{Bo}{Ca}\Bigg]
\hat{k} +Pe_s\Bigg[\frac{2}{3(3\mu+2)}\nonumber\\
&\quad \times\left(Ma_\varGamma\,
I_{10}^{(1)}-\frac{\omega}{3}\,
\frac{Ma_\varGamma}{3\mu+2}\,\frac{Bo}{Ca}
\right)\hat{k}\Bigg] +Pe_s^2 \Bigg[\frac{2}{3(3\mu+2)}
\Bigg(Ma_\varGamma\,I_{10}^{(2)}+\frac{2 \omega\,
Ma_\varGamma^2}{9(3\mu+2)^3}\,\frac{Bo}{Ca}
\Bigg)\hat{k}\nonumber\\
&\quad + \,\frac{2}{3(3\mu+2)}\Bigg(Ma_\varGamma\,I_{11}^{(2)}-\frac{\omega
b
(5\mu+3)\,Ma_\varGamma}{45(3\mu+2)^2(\mu+1)R_0^2}\,\frac{Bo}{Ca}\Bigg)\hat{i}\Bigg]
+O(Pe_s^3). \end{align}
Since a centred monopole does not affect the migration velocity, (B2) is the same as the migration velocity obtained by Das et al. (Reference Das, Mandal and Chakraborty2018).
For an off-centred monopole, the surfactant concentration is
\begin{align} \varGamma &=1+Pe_s
\left[\left(I_{10}^{(1)}-\frac{3\kappa S d\,
Ma_T}{(\kappa+2)(3\mu+2)}\right)\cos{\theta}+I_{21}^{(1)}\cos{\phi}\,
P_2^1 (\cos{\theta})+ I_{30}^{(1)}\,P_3^0
(\cos{\theta})\right]\nonumber\\
&\quad + Pe_s^2
\Bigg[\left(I_{10}^{(2)}+\frac{3\kappa
Sd\,Ma_T\,
Ma_\varGamma}{(\kappa+2)(3\mu+2)^2}\right)\cos{\theta}+I_{21}^{(2)}\cos{\phi}\,
P_2^1 (\cos{\theta}) +I_{30}^{(2)}\,P_3^0
(\cos{\theta})\nonumber\\
&\quad
+\Bigg(I_{11}^{(2)}-\frac{21\kappa Sd b
(5\mu+3)\,Ma_T}{70(\kappa+2)(1+\mu)(3\mu+2)R_0^2}
\Bigg)\cos{\phi}\, P_1^1 (\cos{\theta})\nonumber\\
&\quad +\left(I_{20}^{(2)}-\frac{{\rm \pi}\kappa
Sd\,Ma_T}{7(\kappa+2)(3\mu+2)}\right) P_2^0 (\cos{\theta})
+I_{22}^{(2)} \cos(2\phi)\, P_2^2 (\cos{\theta})\nonumber\\
&\quad +\Bigg(I_{31}^{(2)}+\frac{432\kappa Sd
b\,Ma_T }{1080 R_0^2
(\kappa+2)(3\mu+2)(\mu+1)}\Bigg)\cos{\phi}\, P_3^1
(\cos{\theta})\Bigg]+O(Pe_s^3),
\end{align}and the migration velocity follows:
\begin{align} \boldsymbol{U} &=
\Bigg[1-\left(\frac{b}{R_0}\right)^2-\frac{2\mu}{(3\mu+2)R_0^2}+\frac{2(G+\kappa
Sd)\,Ma_T}{(\kappa+2)(3\mu+2)}+\frac{2\omega}{3}\left(\frac{1+\mu}{3\mu+2}\right)\frac{Bo}{Ca}\Bigg]\hat{k}\nonumber\\
&\quad + Pe_s\left[\frac{2}{3(3\mu+2)}\,Ma_\varGamma
\left(I_{10}^{(1)} -\frac{3\kappa S d\,
Ma_T}{(\kappa+2)(3\mu+2)}-\frac{\omega}{3}\,\frac{Ma_\varGamma}{3\mu+2}\,\frac{Bo}{Ca}
\right)\right]\hat{k}\nonumber\\ &\quad +Pe_s^2
\Bigg[\frac{2}{3(3\mu+2)}\,Ma_\varGamma
\Bigg(I_{20}^{(2)}-\frac{{\rm \pi}\kappa Sd\,
Ma_T}{7(\kappa+2)(3\mu+2)} +\frac{2 \omega\,
Ma_\varGamma^2}{9(3\mu+2)^3}\,\frac{Bo}{Ca}\Bigg)\hat{k}\nonumber\\
&\quad +\, \frac{2}{3(3\mu+2)}\,Ma_\varGamma
\Bigg(I_{11}^{(2)}-\frac{21\kappa Sd b (5\mu+3)\,
Ma_T}{70(\kappa+2)(1+\mu)(3\mu+2)R_0^2}
-\frac{\omega
b (5\mu+3)\,
Ma_\varGamma}{45(3\mu+2)^2(\mu+1)R_0^2}\,\frac{Bo}{Ca}\Bigg)\hat{i}\Bigg]\nonumber\\
&\quad +O(Pe_s^3). \end{align}
For a dipole, the surfactant concentration is
\begin{align} \varGamma &=1+Pe_s
\left[\left(I_{10}^{(1)}-\frac{3\kappa S
\cos\theta\,
Ma_T}{(\kappa+2)(3\mu+2)}\right)\cos{\theta}+\left(-\frac{3\kappa S \sin
\theta\, Ma_T}{(3\mu+2)(\kappa+2)}\right)\cos{\phi}\,P_1^1
(\cos{\theta})\right.\nonumber\\ &\left.\quad +\,I_{21}^{(1)}\cos{\phi}\,P_2^1 (\cos{\theta})+
I_{30}^{(1)}\,P_3^0 (\cos{\theta})\vphantom{\frac{3\kappa S
\cos\theta\, Ma_T}{(\kappa+2)(3\mu+2)}}\right] + Pe_s^2
\left[\left(I_{10}^{(2)}+\frac{3\kappa
S\cos\theta\,Ma_T\,
Ma_\varGamma}{(\kappa+2)(3\mu+2)^2}\right)\cos{\theta}\right.\nonumber\\
&\quad+I_{21}^{(2)}\cos{\phi}\, P_2^1 (\cos{\theta})
+I_{30}^{(2)}\,P_3^0 (\cos{\theta})
+\left(I_{11}^{(2)}-\frac{21\kappa S\cos\theta\,b
(5\mu+3)\,Ma_T}{70(\kappa+2)(1+\mu)(3\mu+2)R_0^2}
\right)\cos{\phi}\nonumber\\ &\quad \times P_1^1
(\cos{\theta})+\left(I_{20}^{(2)}-\frac{{\rm \pi}\kappa
S\cos\theta\, Ma_T}{7(\kappa+2)(3\mu+2)}\right) P_2^0
(\cos{\theta}) +I_{22}^{(2)} \cos(2\phi)\,P_2^2
(\cos{\theta}) \nonumber\\ &\left.\quad +\,\left(I_{31}^{(2)}+\frac{432\kappa
S\cos{\theta}\,b\,Ma_T }{1080 R_0^2
(\kappa+2)(3\mu+2)(\mu+1)}\right)\cos{\phi}\,P_3^1
(\cos{\theta})\right] +O(Pe_s^3),
\end{align}
and the migration velocity is
\begin{align}
\boldsymbol{U} &= \Bigg[\Bigg(1-\Bigg(\frac{b}{R_0}\Bigg)^2-\frac{2\mu}{(3\mu+2)R_0^2}+
\frac{2(G+\kappa
S\cos\theta)\,Ma_T}{(\kappa+2)(3\mu+2)}+\frac{2\omega}{3}\Bigg(\frac{1+\mu}{3\mu+2}\Bigg)\frac{Bo}{Ca}\Bigg)\hat{k}\nonumber\\
&\quad +\,\frac{2\kappa S
\sin{\theta}\,Ma_T}{(\kappa+2)(3\mu+2)}\,\hat{i}\Bigg] +
Pe_s\left[\frac{2}{3(3\mu+2)}\,Ma_\varGamma
\left(I_{10}^{(1)} -\frac{3\kappa S
\cos\theta\,
Ma_T}{(\kappa+2)(3\mu+2)}-\frac{\omega}{3}\,\frac{Ma_\varGamma}{3\mu+2}\right.\right.\nonumber\\
&\quad \times \left.\left.\frac{Bo}{Ca}\right) \hat{k}-
\frac{6\kappa
S \sin \theta\,Ma_T\,
Ma_\varGamma}{3(3\mu+2)^2(\kappa+2)}\,\hat{i}\right]+Pe_s^2
\left[\frac{2}{3(3\mu+2)}\,Ma_\varGamma
\left(I_{20}^{(2)}-\frac{{\rm \pi}\kappa
S\cos\theta\,Ma_T}{7(\kappa+2)(3\mu+2)}\right.\right.\nonumber\\
&\quad +\left.\,\frac{2
\omega\,Ma_\varGamma^2}{9(3\mu+2)^3}\,\frac{Bo}{Ca}\right)
\hat{k}+ \frac{2}{3(3\mu+2)}\,Ma_\varGamma
\left(I_{11}^{(2)}-\frac{21\kappa
S\cos\theta\,b
(5\mu+3)\,Ma_T}{70(\kappa+2)(1+\mu)(3\mu+2)R_0^2}\right.\nonumber\\
&\quad -\left.\left.\,\frac{\omega b
(5\mu+3)\,Ma_\varGamma
}{45(3\mu+2)^2(\mu+1)R_0^2}\,\frac{Bo}{Ca}
\right)\hat{i}\right]+O(Pe_s^3).
\end{align}Here,
$$\begin{gather} I_{10}^{(1)} = \frac{-3 R_0^2G\,Ma_T-2(\kappa+2)}{(\kappa+2)(3\mu+2)R_0^2}-\frac{\omega}{3(3\mu+2)}\,\frac{Bo}{Ca}, \end{gather}$$
$$\begin{gather}I_{21}^{(1)} ={-}\frac{b}{3R_0^2 (1+\mu)}, \end{gather}$$
$$\begin{gather}I_{30}^{(1)} =\frac{1}{6R_0^2(1+\mu)}, \end{gather}$$
$$\begin{gather}I_{10}^{(2)}=\frac{\omega}{3}\,\frac{Ma_\varGamma}{(3\mu+2)^2}\,\frac{Bo}{Ca}+\frac{[3R_0^2G\,Ma_T +2(\kappa+2)]Ma_\varGamma}{(3\mu+2)^2(\kappa+2)R_0^2}, \end{gather}$$
$$\begin{gather}I_{21}^{(2)}={-}\frac{ b\,Ma_\varGamma}{15 R_0^2 (1+\mu)^2}, \end{gather}$$
$$\begin{gather}I_{30}^{(2)}={-}\frac{Ma_\varGamma}{42R_0^2 (1+\mu)^2}, \end{gather}$$
$$\begin{gather}I_{11}^{(2)}={-}V -\frac{b\,[21 R_0^2 (1+\mu)(5\mu+3)G\,Ma_T +(\kappa+2)(70\mu^2+109\mu+40)]}{70(\kappa+2)(1+\mu)^2(3\mu+2)R_0^4}, \end{gather}$$
$$\begin{gather}
I_{20}^{(2)}=\frac{\omega^2}{27(3\mu+2)^2}\left(\frac{Bo}{Ca}\right)^2+\frac{\rm \pi}{21}
[14XI_{10}^{(1)}-6XI_{30}^{(1)}+9 (I_{21}^{(1)})^2+24 I_{30}^{(1)} I_{10}^{(1)}\nonumber\\ +\,4 (I_{30}^{(1)})^2-21WI_{21}^{(1)}], \end{gather}$$
$$\begin{gather}I_{22}^{(2)} ={-}\frac{1}{252}\,\frac{(7\mu+4)b^2}{(1+\mu)^2 R_0^4}, \end{gather}$$
$$\begin{gather}
I_{31}^{(2)} =\frac{2\omega}{45}\,\frac{b}{R_0^2(3\mu^2+5\mu+2)}\,\frac{Bo}{Ca}\nonumber\\+
\frac{b}{R_0^4} \left[
\frac{432(1+\mu)R_0^2 G\,Ma_T +(45\mu^2 +351\mu
+310)(\kappa+2)}{1080(3\mu+2)(\mu+1)^2(\kappa+2)}\right],
\end{gather}$$
and
\begin{equation} X ={-}\frac{3 R_0^2 G\,Ma_T +2(\kappa+2)}{(\kappa+2)(3\mu+2)R_0^2},\quad W=\frac{b}{R_0^2},\quad V =\frac{\omega b}{30R_0^2}\,\frac{5\mu+3}{3\mu^2+5\mu+2}\,\frac{Bo}{Ca}. \end{equation}Appendix C. Comparative analysis of $U_x$ with $\kappa$: small-$Pe_s$ versus large-$Pe_s$ limits
In figure 14, we have shown the variation of cross-stream migration velocity 
$U_x$ with thermal conductivity ratio 
$\kappa$ for both small- and large-
$Pe_s$ limits when the droplet is subjected to a streamwise oriented dipole. In the presence of an ambient thermal field along the Poiseuille flow (
$G=1$), 
$U_x$ is negligibly small for small 
$Pe_s$, so the droplet is dominated by the streamwise migration. However, for large 
$Pe_s$, 
$U_x$ varies significantly, so 
$|U_x|$ increases rapidly with increasing thermal conductivity ratio 
$\kappa$. In the latter case, due to the mass transfer boundary layer around the droplet, the impact of 
$\kappa$ is significant, enhancing the temperature gradient between the droplet and the surrounding medium. Consequently, the stronger thermal gradient increases the Marangoni stress, which enhances the cross-stream migration velocity. Similar behaviour is observed for other thermal singularities also.
Figure 14. For both small- and large-
$Pe_s$ droplets, cross-stream migration velocity versus the thermal conductivity ratio 
$\kappa$ is shown for the case of a streamwise oriented dipole. Other parameters are 
$S=2$, 
$G=1$, 
$Ma_\varGamma =10$, 
$Ma_T=1$, 
$b=1$, 
$R_0=5$ and 
$Bo=0$.
Appendix D. Small $Pe_s$ analysis for a buoyant droplet $(Bo/Ca \neq 0)$
From (B2), (B4) and (B6), we note that though the buoyancy acts in the direction of Poiseuille flow (i.e. along 
$z$ axis), it also affects the cross-stream migration velocity of the droplet, which is further enhanced due to the thermal singularities. We analyse the effect of buoyancy for the following cases: (i) 
$\omega =1$, 
$G=1$, which signifies that the buoyancy and the ambient thermal field act in the flow direction; (ii) 
$\omega =1$, 
$G=-1$, which signifies that the buoyancy acts along the flow, but the temperature field direction is reversed; (iii) 
$\omega =-1$, 
$G=1$, which means that the buoyancy force acts opposite to the flow, and the ambient temperature field is along the flow; and (iv) 
$\omega =-1$, 
$G=-1$, which indicates that the buoyancy and thermal field act opposite to the flow.
Figure 15(a) shows the variation of the cross-stream migration velocity 
$U_x$ with the viscosity ratio 
$\mu$ when 
$Bo/Ca \ne 0$. The magnitude of the cross-stream migration velocity 
$(|{U}_x|)$ is observed to increase significantly compared to the case when 
$Bo/Ca=0$ (as shown in figure 12b). This enhancement is due to the buoyancy force acting in the same direction as the applied temperature gradient and ambient flow. As a result, the Marangoni stress is amplified, causing the droplet to migrate faster towards the flow centreline. Now, when 
$G=-1$ and 
$\omega =1$, from figure 15(b), we can observe that in the absence of a thermal singularity or in the case of a centred monopole (which has no impact on the migration velocity), the droplet moves away from the flow centreline. However, in the presence of an off-centred monopole or a dipole, the combined effect exerted on the droplet becomes stronger, resulting in its movement towards the centreline. Similar effects can also be observed in cases (iii) and (iv).
Figure 15. Droplet cross-stream migration velocity 
${U}_x$ versus viscosity ratio 
$\mu$ for 
$S=0$ (dashed line), centred monopole (dotted line), off-centred monopole (grey), dipole oriented streamwise (orange) and dipole oriented transversely (blue), for (a) 
$G=1$, 
$\omega =1$, and (b) 
$G=-1$, 
$\omega =1$. Other parameters are 
$S=2$, 
$Ma_T=1$, 
$Ma_\varGamma =10$, 
$Pe_s=0.1$, 
$\kappa =1$, 
$b=1$, 
$R_0=5$ and 
$Bo/Ca=1$.
