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),

(2.1a)$$\begin{gather} \tilde{\boldsymbol{\nabla}}\boldsymbol{\cdot}\tilde{\boldsymbol{u}}^j =0, \end{gather}$$

(2.1b)$$\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}$$

(2.1c)$$\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)

(2.2)\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

(2.3a)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}^j =0, \end{gather}$$

(2.3b)$$\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}$$

(2.3c)$$\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

(2.4)\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

(2.5a)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}^j =0, \end{gather}$$

(2.5b)$$\begin{gather}-\boldsymbol{\nabla} p^j + ({\mu^j}/{\mu^e})\,\nabla^2 \boldsymbol{u}^j=0, \end{gather}$$

(2.5c)$$\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

(2.6)\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.1. Thermal field

Both inside and outside of the droplet satisfy the Laplace equation (2.5c) with a complete general solution in the form

(3.1a,b)\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

(3.2a,b)\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,

(3.3)\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

(3.4a,b)\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

(3.5a,b)\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:

(3.6a,b)\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}

(3.7a,b)\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

(3.8)\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)

(3.9)\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:

(3.10a)$$\begin{gather} O(1){:}\quad \boldsymbol{\nabla}_s^2 \varGamma_0=0, \end{gather}$$

(3.10b)$$\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}$$

(3.10c)$$\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

(3.11)\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

(3.12) \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

(3.13a) \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}

(3.13b) \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

(3.14a) \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}

(3.14b) \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

(3.15)\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)

(3.16)\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

(3.17a)$$\begin{gather} O(1){:}\quad \boldsymbol{\nabla}_s\boldsymbol{\cdot}\boldsymbol{u}_{s,0}=0, \end{gather}$$

(3.17b)$$\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

(3.18a)$$\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}$$

(3.18b)$$\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

(3.19a,b)\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

(3.20a,b)\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).