1. Introduction
Migration of droplets plays an important role in microfluidic devices (Pamme Reference Pamme2007). Some important applications involving droplet manipulation include cell sorting, separation of particles and drug delivery (Gascoyne, Satayavivad & Ruchirawat Reference Gascoyne, Satayavivad and Ruchirawat2004; Stone, Stroock & Ajdari Reference Stone, Stroock and Ajdari2004; Di Carlo et al. Reference Di Carlo, Irimia, Tompkins and Toner2007; Karimi, Yazdi & Ardekani Reference Karimi, Yazdi and Ardekani2013). In such applications, accurate control of the droplet motion is important for realizing optimum process capabilities. The Reynolds number in such applications is typically much smaller than unity (
$Re \ll 1$) due to the small channel dimensions as well as the droplet size. It is known that a non-deformable droplet placed in a purely viscous Newtonian fluid only translates in the flow direction with a lower velocity than the background flow at the drop centre (Hetsroni & Haber Reference Hetsroni and Haber1970). It does not experience any lateral motion across streamlines due to the reversibility of the Stokes flow. However, the presence of non-linearities can break the fore-aft symmetry and, consequently, the reversible nature of the flow. In such a case, the droplet can exhibit a cross-stream migration as well. Some significant mechanisms through which non-linearities can be introduced in the system are in the form of fluid inertia, particle surface deformation and non-Newtonian rheology of the suspending fluid (Leal Reference Leal2007). Initial works started with the analysis of migration of solid particles in the presence of fluid inertia (Rubinow & Keller Reference Rubinow and Keller1961; Segre & Silberberg Reference Segre and Silberberg1961; Saffman Reference Saffman1965; Cox & Brenner Reference Cox and Brenner1968; Ho & Leal Reference Ho and Leal1974; Vasseur & Cox Reference Vasseur and Cox1976). Subsequently, several works have focused on the cross-stream migration of droplets due to surface deformation (Haber & Hetsroni Reference Haber and Hetsroni1971; Wohl & Rubinow Reference Wohl and Rubinow1974; Chan & Leal Reference Chan and Leal1979; Stan et al. Reference Stan, Guglielmini, Ellerbee, Caviezel, Stone and Whitesides2011; Mandal, Bandopadhyay & Chakraborty Reference Mandal, Bandopadhyay and Chakraborty2015, Reference Mandal, Bandopadhyay and Chakraborty2016) and the non-Newtonian rheology of the suspending fluid (Chan & Leal Reference Chan and Leal1979; Mukherjee & Sarkar Reference Mukherjee and Sarkar2013, Reference Mukherjee and Sarkar2014; Raffiee, Dabiri & Ardekani Reference Raffiee, Dabiri and Ardekani2017; Raffiee, Ardekani & Dabiri Reference Raffiee, Ardekani and Dabiri2019).
Most of the microfluidic systems involving manipulation of droplets are expected to contain surfactants which might be present as contaminants, additives or impurities (Stan et al. Reference Stan, Ellerbee, Guglielmini, Stone and Whitesides2013). The surfactant molecules have two distinct chemical moieties which allows them to adsorb at interfaces between immiscible liquids. In the presence of an external flow, surface advection causes the surfactant molecules to redistribute along the interface. An asymmetric redistribution leads to interfacial tension gradients along the interface which drives a flow in the surrounding fluid known as the Marangoni flow. This phenomenon breaks the fore-aft symmetry of the problem and thereby the reversibility which causes a surfactant-laden droplet to migrate across the streamlines. This is a relatively unexplored phenomenon in literature which compels us to focus on this problem.
Theoretical work in this direction started within the last decade, when Hanna & Vlahovska (Reference Hanna and Vlahovska2010) predicted the cross-stream migration velocity of a surfactant-laden droplet in a Poiseuille flow with negligible surface diffusion (
$Pe_{s} \to \infty$), where 
$Pe_{s}$ represents the surface Péclet number. They performed their analysis for the cases of high Marangoni number and high drop viscosity, by employing a regular perturbation expansion approach to solve the problem. Through their analysis, they found that, for both the cases, the drop always migrates towards the centre of the channel. Later, Pak et al. (Reference Pak, Feng and Stone2014) extended their analysis as they explored the opposite limit of small surface Péclet numbers (
$Pe_{s} \ll 1$). They employed a regular perturbation expansion method with the small parameter as 
$Pe_{s}$ combined with an efficient use of the Lorentz reciprocal theorem to bypass the detailed calculations of the full-flow problem. Interestingly, they found that, the cross-stream velocity is zero up to first order in the surface Péclet number and the first non-zero term only appears at the second order. Again, they observed that the drop migrates towards the centre of the channel with a velocity directly proportional to the distance of the drop from the centreline of the Poiseuille flow. Das, Mandal & Chakraborty (Reference Das, Mandal and Chakraborty2017) considered the droplet to be deformable and predicted the migration velocity of the droplet in terms of an asymptotic expansion with the small parameter as the capillary number (
$Ca$). They found that, for small surface Péclet numbers (
$Pe_{s} \ll 1$), the drop has a tendency to migrate towards or away from the centreline, depending on the viscosity ratio of the fluid inside and outside the drop. More recently, Santra et al. (Reference Santra, Das, Das and Chakraborty2018) experimentally measured the cross-stream migration velocity of a surfactant-laden droplet subjected to an imposed flow in a confined fluidic environment. Panigrahi et al. (Reference Panigrahi, Santra, Banuprasad, Das and Chakraborty2018) explored the effect of interfacial viscosities on the droplet motion in a non-isothermal Poiseuille flow. Gounley et al. (Reference Gounley, Boedec, Jaeger and Leonetti2016) provided a comprehensive numerical analysis regarding the effect of surface viscosities on droplets immersed in a shear flow. However, in their analysis, they neglect the surfactant transport on the surface of the drop and the Marangoni stresses which arise from gradients in the surface tension.
The variation of the surface tension with the surfactant concentration gives rise to Gibbs elasticity which provides a measure of the sensitivity of the surface tension with the surfactant concentration, as well as Marangoni stresses at the interface (Leal Reference Leal2007). In addition, surfactants also give rise to a surface excess rheology in comparison with a clean interface. In particular, the interface between two fluids can exhibit a complex rheological behaviour, which can comprise a surface viscosity as well as surface elasticity. In the present analysis, we neglect the elastic effects on the surface rheology and only focus on the effects of the surface viscosity. The surface viscosity arises from the rate dependent resistance of the interface to velocity gradients in the plane of the interface. We thus expect that interfacial viscosity combined with the interfacial tension will provide a more realistic depiction of the fluid–fluid interface. The simplest model which describes the mechanics of an interface is the Bousinessq–Scriven law (Boussinesq Reference Boussinesq1913; Scriven Reference Scriven1960) which treats the interface as a homogeneous two-dimensional fluid with surface shear and dilatation viscosities. Using the Boussinesq–Scriven model to represent the interfacial viscosities, Schwalbe et al. (Reference Schwalbe, Phelan, Vlahovska and Hudson2011) predicted the cross-stream migration velocity of an eccentrically placed surfactant-covered droplet in a Poiseuille flow with negligible surface diffusion (
$Pe_{s} \to \infty$). In addition, they found that the presence of interfacial stresses significantly altered the slip velocity and the droplet-circulation velocities.
This study focuses on the effect of the interfacial stresses on the droplet dynamics in the limit of small surface Péclet number (
$Pe_{s} \ll 1$). Analysis of this limit will allow us to close the current literature gap and provide a more complete picture regarding the cross-stream migration of a surfactant-covered droplet in an unbounded Poiseuille flow (refer to table 1). In this work, we explore the migration of a surfactant-covered non-deformable droplet in an unbounded Poiseuille flow including the effects of interfacial stresses, at low surface Péclet numbers. We assume that the surfactant is insoluble in the bulk fluid and that a small amount of surfactant is adsorbed on the drop surface. We note that the assumption of a dilute surfactant concentration considerably limits the extension of our analysis to biological flows, where the surfactants are expected to be more tightly packed and the intermolecular interactions between the adsorbed surfactant molecules should be considered (Manikantan & Squires Reference Manikantan and Squires2020). We use the Bousinessq–Scriven law to represent the interface as a two-dimensional fluid with surface shear and dilatation viscosities. This allows us to examine the effect of interfacial rheology on the droplet dynamics. We neglect the elasticity of the interface and assume the interfacial viscosities to be independent of the surfactant concentration. Section 2 summarizes the problem formulation and the governing equations for the problem. Section 3 focuses on the solution methodology and important results from our study.
List of previous works studying the cross-stream migration of an eccentrically placed surfactant-laden droplet in an unbounded Poiseuille flow. Here 
$Pe_{s}$ represents the surface Péclet number and 
$\eta _{\mu }$ and 
$\eta _{\kappa }$ represent the surface shear and dilatation viscosity, respectively.
2. Problem formulation
Figure 1 shows a surfactant-covered neutrally buoyant droplet of radius 
$a$ and viscosity 
$\lambda \eta$ immersed in a fluid of viscosity 
$\eta$ and density 
$\rho$. There are two coordinate systems defined in the manuscript 
$(x,y,z)$ and 
$(r,\theta ,\phi )$, where 
$\theta$ is the angle made by the vector shown as 
$\boldsymbol {r}$ in the figure with the 
$z$-axis and 
$\phi$ is the angle made by the projection of the vector 
$\boldsymbol {r}$ on the 
$x$–
$y$ plane with the 
$x$-axis. The dashed blue line represents the projection operation which transforms the vector 
$\boldsymbol {r}$ to a vector in the 
$x$–
$y$ plane. The migration velocity of the droplet is represented as 
$\boldsymbol {U}$. The droplet is subjected to an external Poiseuille flow which is given by,
\begin{equation} \boldsymbol{U}_{\infty}=U_{c}\left(1-\frac{r_{0}^{2}}{R^{2}}\right)\boldsymbol{e}_{z} \end{equation}
The flow is along the 
$\boldsymbol{e}_{z}$ direction. Here, 
$r_{0}$ is the radial distance measured from the flow axis, 
$U_{c}$ is the fluid velocity at 
$r_{0}=0$ and R is the radial distance from the centreline at which the external flow is stationary. The droplet is eccentrically located at a radial distance of 
$b$ from the flow axis. In the absence of any flow, the surfactant distributes itself uniformly on the droplet surface at equilibrium. The equilibrium surfactant concentration is given by 
$\Gamma _{eq}$. For a clean drop, the interfacial tension is represented by 
$\sigma _{c}$. However, due to the presence of the surfactant, the interfacial tension is lowered. If the surfactant concentration is small, the surface tension can be described by the following linear relationship (Leal Reference Leal2007),
\begin{equation} \sigma=\sigma_{c}- R_{g} T\Gamma \end{equation}
Here, 
$\Gamma$ denotes the surfactant concentration, 
$R_{g}$ is the gas constant and 
$T$ is the ambient temperature. We assume that the surfactant is insoluble in both the fluids inside and outside the drop. The interfacial viscosities are represented by 
$\eta _{\mu }$ and 
$\eta _{\kappa }$, which signify the surface shear and dilatational viscosities, respectively.
Surfactant-laden droplet of radius 
$a$ immersed in a cylindrical Poiseuille flow at a radial distance of 
$b$ from the flow axis. 
$(x,y,z)$ and 
$(r,\theta ,\phi )$ represent the Cartesian and spherical coordinate systems attached to the centre of the droplet, respectively. Here, 
$\theta$ is the angle made by the vector shown as 
$\boldsymbol {r}$ in the figure with the 
$z$-axis, 
$\phi$ is the angle made by the projection of the vector 
$\boldsymbol {r}$ on the 
$x$–
$y$ plane (shown as a blue vector) with the 
$x$-axis. The dashed blue line represents the 
$z$-coordinate of the vector 
$\boldsymbol {r}$.
We solve the problem in the creeping flow limit (
$Re=0$). In this limit, we can neglect the inertial forces exerted on the drop and the velocity field can be assumed to be set up instantaneously. We also assume that the drop is non-deformable. This assumption is valid if the surface tension forces dominate the viscous forces and the capillary number 
$\eta U_{c}/\sigma _{c}$ is small. We use a quasi-steady approximation and assume that the time taken for the surfactant redistribution is very small compared to the time scale associated with the drop migration. Owing to the use of a regular perturbation methodology with 
$Pe_{s}$ as the small parameter, we expect the characteristic velocity scale in the transverse direction to be expressed as 
$Pe_{s}^{n}U_{c}$. Here, 
$n \geq ~1$ and represents the order of approximation at which the transverse migration velocity first becomes non-zero. When interfacial stresses are not considered, Pak et al. (Reference Pak, Feng and Stone2014) showed that the migration velocity is zero up to the first order and thus 
$n=2$. Although, in our problem, we cannot predict the value of 
$n \ {\rm a} \ {\rm priori}$, we expect the magnitude of the velocity scale in the transverse direction to be at least an order of magnitude less than the longitudinal velocity scale 
$U_{c}$ as 
$Pe_{s}\ll 1$. Hence, in order to check the validity of the quasi-steady approximation, the surfactant redistribution time scale (
$a^{2}/D_{s}$) should be compared with the time scale associated with the motion of the droplet in the longitudinal direction (
$a/U_{c}$), where D s is the diffusivity of the surfactant. This implies that 
$Pe_{s} \ll 1$ which is the regime of analysis in our work.
2.1. Governing equations
We express the equations of motion in a frame of reference translating with the drop (refer to the (
$r,\theta ,\phi$) coordinate system in figure 1). The velocity and pressure fields inside and outside the drop are represented by 
$(\boldsymbol {u}_{i},p_{i})$ and 
$(\boldsymbol {u}_{e},p_{e})$, respectively. We make the governing equations dimensionless by using the following scales: length scale 
$(l_{c})=a$, velocity scale 
$(u_{c})=U_{c}$, pressure scale 
$(p_{c})=\eta U_{c}/a$. As inertia is considered to be negligible, the creeping flow dimensionless equations inside and outside the drop are given by,
\begin{gather} -\boldsymbol{\nabla} p_{i}+\lambda \nabla^{2} \boldsymbol{u}_{i}=\boldsymbol{0},\quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}_{i}=0, \end{gather}
\begin{gather}-\boldsymbol{\nabla} p_{e}+\nabla^{2} \boldsymbol{u}_{e}=\boldsymbol{0}, \quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}_{e}=0. \end{gather}
Here 
$\lambda$ represents the viscosity ratio between the droplet and the suspending fluid. Additionally, we have the surfactant transport equation for an insoluble surfactant which describes the steady-state convection and diffusion processes on the interface as follows (Leal Reference Leal2007):
\begin{equation} Pe_{s}\left(\boldsymbol{\nabla}_{s}\boldsymbol{\cdot}(\Gamma \boldsymbol{u_{s}})+\Gamma(\boldsymbol{\nabla}_{s}\boldsymbol{\cdot} \boldsymbol{n})(\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n})\right)=\nabla_{s}^{2}\Gamma \end{equation}
We have made the equation dimensionless by using 
$\Gamma _{eq}$ as the characteristic concentration scale. Here, 
$Pe_{s}=U_{c}a/D_{s}$ is the surface Péclet number, where 
$D_{s}$ represents the diffusivity of the surfactant, 
$\Gamma$ is the surfactant concentration at the interface, 
$\boldsymbol {u_{s}}$ is the interface velocity, 
$\boldsymbol {n}$ is the normal vector at the interface and 
$\boldsymbol {\nabla }_{s}=\boldsymbol {\nabla }-\boldsymbol {n}(\boldsymbol {n}\boldsymbol {\cdot } \boldsymbol {\nabla })$ represents the gradient operator in the plane of the interface. Here, it should be noted that the unsteady term has been neglected in the above equation due to the quasi-steady approximation discussed in the previous section.
2.2. Boundary conditions
Far away from the drop, the velocity field is given by the difference between the background flow and the drop migration velocity as follows:
\begin{equation} \boldsymbol{u}_{e}(r \to \infty)= \boldsymbol{U}_{\infty}-\boldsymbol{U}. \end{equation}At the surface of the droplet, we apply the condition of velocity continuity and the kinematic boundary condition for a non-deforming interface to yield the following:
\begin{gather} \boldsymbol{u}_{i}=\boldsymbol{u}_{e} \quad \mbox{at}\ r=1, \end{gather}
\begin{gather}\boldsymbol{u}_{i} \boldsymbol{\cdot} \boldsymbol{n}=0 \quad \mbox{at}\ r=1. \end{gather}Finally, we write the tangential stress balance at the surface of the drop as follows:
\begin{equation} \boldsymbol{n}\boldsymbol{\cdot} \left(\boldsymbol{\sigma}_{e}-\lambda \boldsymbol{\sigma}_{i}\right) \boldsymbol{\cdot} \left( \boldsymbol{\mathsf{I}}-\boldsymbol{n}\boldsymbol{n}\right)=Ma\boldsymbol{\nabla}_{s}\Gamma-\boldsymbol{\nabla}_{s}\boldsymbol{\cdot} \boldsymbol{\tau}_{s}. \end{equation}
Here, 
$\boldsymbol{\sigma}_{e}, \boldsymbol{\sigma}_{i}$ represent the hydrodynamic stress tensors in the suspending fluid and the droplet, respectively, and 
$Ma=RT\Gamma _{eq}/U_{c}\eta$ is the Marangoni number which signifies the relative strength of the characteristic Marangoni stresses and viscous stresses. We have used (2.2) to express the surface tension in terms of the surfactant concentration. The interfacial stress which arises due to the presence of surface viscosities is represented by 
$\boldsymbol {\tau }_{s}$. We use the Bousinessq–Scriven constitutive law to represent the interfacial stress as follows:
\begin{equation} \boldsymbol{\tau}_{s}= 2Bq_{\mu}\boldsymbol{\mathsf{D}}_{s}+(Bq_{\kappa}-Bq_{\mu})(\boldsymbol{\nabla}_{s} \boldsymbol{\cdot} \boldsymbol{u}_{s})\boldsymbol{\mathsf{I}}_{s} \end{equation}
Here, 
$\boldsymbol{\mathsf{D}}_{s}=\frac {1}{2}\left ((\boldsymbol {\nabla }_{s}\boldsymbol {u}_{s})\boldsymbol {\cdot } \boldsymbol{\mathsf{I}}_{s}+\boldsymbol{\mathsf{I}}_{s}\boldsymbol {\cdot } (\boldsymbol {\nabla }_{s}\boldsymbol {u}_{s})^{T}\right )$, 
$\boldsymbol{\mathsf{I}}_{s}=\boldsymbol{\mathsf{I}}-\boldsymbol {n}\boldsymbol {n}$, where 
$\boldsymbol{\mathsf{I}}$ is the identity tensor, 
$Bq_{\mu }=\eta _{\mu }/\eta a$ is the Bousinessq number for surface shear viscosity and 
$Bq_{\kappa }=\eta _{\kappa }/\eta a$ is the Bousinessq number for surface dilatation viscosity. In the normal direction, the viscous pressure and the stress contributions across the interface are balanced by the capillary pressure and a contribution from the surface tension variation with the surfactant concentration. In particular, the normal stress balance across the interface can be written as follows:
\begin{equation} \boldsymbol{n}\boldsymbol{\cdot} \left(\boldsymbol{\sigma}_{e}-\lambda \boldsymbol{\sigma}_{i}\right) \boldsymbol{\cdot} \boldsymbol{n}=\frac{1}{Ca}\boldsymbol{\nabla}_{s}\boldsymbol{\cdot}\boldsymbol{n}- \frac{\beta \Gamma}{Ca}\boldsymbol{\nabla}_{s}\boldsymbol{\cdot}\boldsymbol{n}. \end{equation}
Here, 
$\beta =RT\Gamma _{eq}/\sigma _{c}$ is the Gibbs elasticity parameter which measures the sensitivity of the variation of the surface tension with the surfactant concentration. Typical values of 
$\beta$ lie between 
$0$ and 
$1$ (Das et al. Reference Das, Mandal and Chakraborty2017). This implies that, if 
$Ca \ll 1$, the difference in the normal stresses can be balanced by a very small contribution from the term 
$\boldsymbol {\nabla }_{s} \boldsymbol {\cdot } \boldsymbol {n}$ and the droplet deformation can be neglected (Leal Reference Leal2007)
Additionally, we note that, due to the kinematic boundary condition (refer to (2.6b)), the second term in the surfactant transport equation (refer to (2.4)) is identically zero.
3. Solution
In order to solve (2.3) and (2.4), we use a regular perturbation expansion approach with 
$Pe_{s}$ as the small parameter. We express the unknown variables in the following form:
\begin{equation} h=h^{(0)}+Pe_{s}h^{(1)}+Pe_{s}^{2}h^{(2)}+O(Pe_{s}^{3}), \quad h=\{\boldsymbol{u}_{i},\boldsymbol{u}_{e},p_{i},p_{e},\Gamma,\boldsymbol{U}\}. \end{equation}
As the velocity fields at every order in 
$Pe_{s}$ satisfy the Stokes equations (refer to (2.3)), we use Lamb's general solution to express the velocity and pressure fields in terms of spherical harmonics. The velocity and pressure fields inside the droplet can be represented in terms of growing spherical harmonics as follows:
\begin{gather} \boldsymbol{u}_{i}^{(q)}=\sum_{n=1}^{\infty}\boldsymbol{\nabla} \times (\boldsymbol{r}\chi_{n})+\boldsymbol{\nabla} \phi_{n}+\frac{n+3}{2(n+1)(2n+3)\lambda}r^{2} \boldsymbol{\nabla} p_{n}-\frac{n}{(n+1)(2n+3)\lambda}\boldsymbol{r}p_{n}, \end{gather}
\begin{gather}p_{i}^{(q)}=\sum_{n=0}^{\infty}p_{n} \end{gather}
Expressions for 
$p_{n}$, 
$\phi _{n}$ and 
$\chi _{n}$ are given in appendix A. Similarly, we can express the velocity field outside the droplet as follows:
\begin{align} \boldsymbol{u}_{e}^{(q)}&=(\boldsymbol{U}_{\infty}-\boldsymbol{U})+ \sum_{n=1}^{\infty}\boldsymbol{\nabla} \times (\boldsymbol{r}\chi_{-n-1})+ \boldsymbol{\nabla} \phi_{-n-1}-\frac{n-2}{2n(2n-1)}r^{2}\boldsymbol{\nabla} p_{-n-1} \nonumber\\ &\quad +\frac{n+1}{n(2n-1)}\boldsymbol{r}p_{-n-1}, \end{align}
\begin{gather} p_{e}^{(q)}=\sum_{n=0}^{\infty}p_{-n-1} \end{gather}
Expressions for 
$p_{-n-1}$, 
$\phi _{-n-1}$ and 
$\chi _{-n-1}$ are given in appendix A. Finally, we express the surfactant concentration at the drop surface in terms of solid spherical harmonics as follows:
\begin{equation} \Gamma^{(q)}=\sum_{n=0}^{\infty}\sum_{m=0}^{n}\left(\Gamma_{n,m}^{(q)}\cos (m\phi)+ \hat{\Gamma}_{n,m}^{(q)}\sin (m\phi)\right)P_{n,m}(\cos \theta). \end{equation}
The unknown coefficients are determined by substituting the velocity fields as described above in the surfactant transport equation (refer to (2.4)) and the boundary conditions at the drop surface (refer to (2.6) and (2.7)). Here, we note that the drop migration velocity 
$\boldsymbol {U}$ is still an unknown and yet to be determined. As the inertia of the drop is neglected in our analysis, the net force exerted on the drop is zero. For a neutrally buoyant drop, this implies that the hydrodynamic force acting on the drop is zero. Consequently, we have an additional equation as follows,
\begin{equation} \boldsymbol{F}_{h}^{(q)}=4{\rm \pi} \boldsymbol{\nabla}\left(r^{3}p_{-2}\right)=\boldsymbol{0}. \end{equation}
Next, we substitute the perturbation expansions of the velocity, pressure and the surfactant fields in the equations of motion and obtain the governing equations and boundary conditions at various orders in 
$Pe_{s}$.
3.1. Solution at leading order
At leading order, the creeping flow and the surfactant transport equations are given by
\begin{gather} -\boldsymbol{\nabla} p_{i}^{(0)}+\lambda \nabla^{2}\boldsymbol{u}_{i}^{(0)}=\boldsymbol{0}, \end{gather}
\begin{gather}-\boldsymbol{\nabla} p_{e}^{(0)}+\nabla^{2}\boldsymbol{u}_{e}^{(0)}=\boldsymbol{0}, \end{gather}
\begin{gather}\nabla^{2}\Gamma^{(0)}=0. \end{gather}We observe that, at leading order, the fluid flow has no effect on the surfactant concentration. Consequently, the surfactant is distributed uniformly along the drop surface. The Stokes equations are subjected to the following boundary conditions:
\begin{gather} \boldsymbol{u}_{i}^{(0)}=\boldsymbol{u}_{e}^{(0)} \quad \mbox{at}\ r=1, \end{gather}
\begin{gather}\boldsymbol{u}_{i}^{(0)} \boldsymbol{\cdot} \boldsymbol{n}=0 \quad \mbox{at}\ r=1. \end{gather}
\begin{gather}\boldsymbol{n}\boldsymbol{\cdot} \left(\boldsymbol{\sigma}_{e}^{(0)}-\lambda \boldsymbol{\sigma}_{i}^{(0)}\right) \boldsymbol{\cdot} \left(\boldsymbol{\mathsf{I}}-\boldsymbol{n}\boldsymbol{n}\right)=-\boldsymbol{\nabla}_{s} \boldsymbol{\cdot}\boldsymbol{\tau}_{s}^{(0)}. \end{gather}We use (3.2) and (3.3) to express the velocity fields inside and outside the drop, respectively. We then obtain the unknown coefficients in the velocity field expansions by satisfying the above boundary conditions and the force-free condition expressed in (3.5). We find that, at leading order, the drop only translates in the direction of the flow without any motion across the streamlines. The translational velocity of the drop at leading order is given as follows:
\begin{equation} \boldsymbol{U}^{(0)}= \left(1-\frac{b^{2}}{R^{2}}-\frac{2a^{2}(2Bq_{\kappa}+3\lambda)}{3R^{2}(2+2Bq_{\kappa}+3\lambda)}\right) \boldsymbol{e}_{z}. \end{equation}
As the surfactant effects influence the droplet motion only at higher orders, Schwalbe et al. (Reference Schwalbe, Phelan, Vlahovska and Hudson2011) obtained the same expression for the drop migration velocity at the leading order for the limit of 
$Pe_{s} \to \infty$, which is the opposite of the regime considered in our analysis. We note that the drop lags the flow for all values of 
$Bq_{\kappa }$ and the migration velocity is independent of the surface shear viscosity. This observation can be supported by studying the leading-order velocity vectors as shown in figures 2(a,c,e). From the flow patterns, we find that shear viscosity (
$Bq_{\mu }$) has a negligible effect on the leading-order velocity field (refer to figures 2a,e), while the dilatation viscosity (
$Bq_{\kappa }$) alters the circulation pattern inside the drop significantly (refer to figure 2c). The detailed expression for the leading-order velocity field is provided in appendix B. Note that the surfactant is bulk insoluble and the surfactant concentration in the bulk is zero which is indicated by the purple colour in figure 2.
Panels (a,c,e) show the leading-order velocity vectors both inside (shown in blue) and outside the drop (shown in black) superimposed with the leading-order surfactant concentration field. Panels (b,d,f) show the first-order velocity vectors both inside (shown in blue) and outside the drop (shown in black) superimposed with the first-order surfactant concentration field for (a,b) 
$Bq_{\kappa }=0,Bq_{\mu }=0$, (c,d) 
$Bq_{\kappa }=100,Bq_{\mu }=0$ and (e,f) 
$Bq_{\kappa }=0,Bq_{\mu }=100$. For all the cases shown, 
$a/R=1,\ b/R=1,\ \lambda =1,\ Ma=1$. The flow direction is from left to right.
3.2. Solution at first order
At first order, the creeping flow and the surfactant transport equations are given by,
\begin{gather} -\boldsymbol{\nabla} p_{i}^{(1)}+\lambda \nabla^{2}\boldsymbol{u}_{i}^{(1)}=\boldsymbol{0}, \end{gather}
\begin{gather}-\boldsymbol{\nabla} p_{e}^{(1)}+\nabla^{2}\boldsymbol{u}_{e}^{(1)}=\boldsymbol{0}, \end{gather}
\begin{gather}\nabla^{2}\Gamma^{(1)}=\boldsymbol{\nabla}_{s} \boldsymbol{\cdot} \boldsymbol{u}_{s}^{(0)}. \end{gather}The above equations are subjected to the following boundary conditions,
\begin{gather} \boldsymbol{u}_{i}^{(1)}=\boldsymbol{u}_{e}^{(1)} \quad \mbox{at}\ r=1, \end{gather}
\begin{gather}\boldsymbol{u}_{i}^{(1)} \boldsymbol{\cdot} \boldsymbol{n}=0 \quad \mbox{at}\ r=1. \end{gather}
\begin{gather}\boldsymbol{n}\boldsymbol{\cdot} \left(\boldsymbol{\sigma}_{e}^{(1)}-\lambda \boldsymbol{\sigma}_{i}^{(1)}\right) \boldsymbol{\cdot} \left(\boldsymbol{\mathsf{I}}-\boldsymbol{n}\boldsymbol{n}\right)=Ma \boldsymbol{\nabla}_{s}\Gamma^{(1)}-\boldsymbol{\nabla}_{s}\boldsymbol{\cdot}\boldsymbol{\tau}_{s}^{(1)}. \end{gather}From the surfactant transport equation, we observe that, at first order, the surfactant will redistribute itself along the interface due to convection of the leading-order velocity along the drop surface. Consequently, the Marangoni stress term appears in the tangential stress continuity boundary condition. After solving the equations of motion, we obtain the equilibrium surfactant concentration as follows:
\begin{align} \Gamma^{(1)} &= \frac{a^{2}}{R^{2}}\left(-\frac{2}{2+2Bq_{\kappa}+3\lambda}P_{1}^{0}(\cos \theta)+ \frac{7}{7+12Bq_{\kappa}+10Bq_{\mu}+7\lambda}P_{3}^{0}(\cos \theta) \right.\nonumber\\ &\quad \left.+\frac{5b}{a(5+6Bq_{\kappa}+4Bq_{\mu}+5\lambda)}P_{2}^{1}(\cos \theta)\cos \phi\right). \end{align}
However, we find that, despite the surfactant redistribution, the drop migrates only in the direction of the flow and the cross-stream migration velocity is zero up to the first order in 
$Pe_{s}$. Here, we note that Pak et al. (Reference Pak, Feng and Stone2014) made a similar observation for a non-viscous interface. The first-order migration velocity is given as follows:
\begin{equation} \boldsymbol{U}^{(1)}=-\frac{a^{2}}{R^{2}}\left(\frac{4Ma}{3(2+2Bq_{\kappa}+3\lambda)^{2}}\right) \boldsymbol{\hat{e}}_{z}. \end{equation} The negative sign indicates that the first-order velocity causes the drop to further lag the flow. This lag is decreased due to the surface dilatation viscosity as 
$Bq_{\kappa }$ appears in the denominator. However, the surface shear viscosity has no effect on the droplet longitudinal speed. The first-order velocity field vectors are displayed in figures 2(b,d,f) for different values of surface viscosities. In all the cases, the velocity vectors outside the drop are oriented in a direction opposite to the external flow, which reduces the longitudinal speed of the droplet as seen from (3.16). We further observe that an increase in the surface shear viscosity increases the resistance of the drop to surface shear forces, which consequently reduces the compression of the streamlines in the vicinity of the drop in comparison with the droplet with zero interfacial viscosity (refer to figures 2b,f). However, surface dilatation viscosity does not cause any significant changes in the streamlines pattern. The detailed expression for the first-order velocity field is provided in appendix B.
From the surfactant concentration fields in figure 2, we observe that, in all the cases, the first-order velocity creates a surplus of surfactant (shown in red) at the rear end of the drop and a deficit at the front end (shown in blue). This generates a Marangoni flow from the rear end to the front end of the droplet as shown in the figure. For a drop without any interfacial viscosities, as the drop is located away from the flow axis, we observe that the peaks in the surfactant concentration are not on the 
$z$-axis passing through the sphere centre (refer to figure 2b), which is the direction of the external flow. Introduction of the dilatation viscosity generates a similar surfactant distribution pattern (refer to figure 2d); however, the shear viscosity causes the surfactant concentration peaks to align with the sphere axis (refer to figure 2f).
3.3. Solution at second order
At second order, the creeping flow and the surfactant transport equations are given by
\begin{gather} -\boldsymbol{\nabla} p_{i}^{(2)}+\lambda \nabla^{2}\boldsymbol{u}_{i}^{(2)}=\boldsymbol{0}, \end{gather}
\begin{gather}-\boldsymbol{\nabla} p_{e}^{(2)}+\nabla^{2}\boldsymbol{u}_{e}^{(2)}=\boldsymbol{0}, \end{gather}
\begin{gather}\nabla^{2}\Gamma^{(2)}=\boldsymbol{\nabla}_{s} \boldsymbol{\cdot} (\Gamma^{(0)}\boldsymbol{u}_{s}^{(1)}+\Gamma^{(1)}\boldsymbol{u}_{s}^{(0)}). \end{gather}The above equations are subjected to the following boundary conditions:
\begin{gather} \boldsymbol{u}_{i}^{(2)}=\boldsymbol{u}_{e}^{(2)} \quad \mbox{at}\ r=1, \end{gather}
\begin{gather}\boldsymbol{u}_{i}^{(2)} \boldsymbol{\cdot} \boldsymbol{n}=0 \quad \mbox{at}\ r=1. \end{gather}
\begin{gather}\boldsymbol{n}\boldsymbol{\cdot} \left(\boldsymbol{\sigma}_{e}^{(2)}-\lambda \boldsymbol{\sigma}_{i}^{(2)}\right) \boldsymbol{\cdot} \left(\boldsymbol{\mathsf{I}}-\boldsymbol{n}\boldsymbol{n}\right)=Ma \boldsymbol{\nabla}_{s}\Gamma^{(2)}-\boldsymbol{\nabla}_{s}\boldsymbol{\cdot}\boldsymbol{\tau}_{s}^{(2)}. \end{gather}
We find that the drop exhibits a cross-stream migration velocity in the lateral direction at this order in 
$Pe_{s}$. The net migration velocity of the drop at this order is given by
\begin{equation} \boldsymbol{U}^{(2)}= -\frac{a^{2}}{R^{2}}\left(\frac{4Ma}{3(2+2Bq_{\kappa}+3\lambda)^{3}}\right) \boldsymbol{e}_{z}-\frac{a^{3}b}{R^{4}}f(Ma,\lambda,Bq_{\kappa},Bq_{\mu})\boldsymbol{e}_{x}. \end{equation}
Here, 
$f(Ma,\lambda ,Bq_{\kappa },Bq_{\mu })$ is given by
\begin{align} \left(\frac{\begin{array}{l} Ma(40+109\lambda+70\lambda^{2}+144Bq_{\kappa}^{2}+80Bq_{\mu}^{2}+4Bq_{\mu}(29+39\lambda)\\ \quad +2Bq_{\kappa}(77+108Bq_{\mu}+102\lambda)) \end{array}}{3(2+3\lambda)(2+2Bq_{\kappa}+3\lambda)(5+6Bq_{\kappa}+ 4Bq_{\mu}+5\lambda)(7+12Bq_{\kappa}+10Bq_{\mu}+7\lambda)}\right) \end{align}
From the above expression, we observe that inclusion of the interfacial stress does not change the direction of the droplet migration and it always moves towards the centre of the channel irrespective of the values of the Bousinessq numbers. Additionally, we find that the cross-stream migration velocity depends on both the interfacial viscosities. If we substitute 
$Bq_{\mu }=Bq_{\kappa }=0$, we obtain the expression for the cross-stream migration velocity obtained by Pak et al. (Reference Pak, Feng and Stone2014) for a non-viscous interface as follows:
\begin{equation} U^{(2)}_{x}=-\frac{a^{3}b}{R^{4}}\left(\frac{Ma(40+109\lambda+70\lambda^{2})}{105(1+\lambda)^{2}(2+3\lambda)^{2}}\right) \end{equation}The detailed expression for the surfactant concentration at this order is provided in appendix B.
Figure 3 shows the effect of the interfacial viscosities on the cross-stream velocity of the droplet (shown by symbol 
$V$). Here, 
$V_{h}$ denotes the cross-stream migration velocity for a droplet with a non-viscous interface (
$Bq_{\kappa }=Bq_{\mu }=0$). We observe that both the viscosities influence the droplet motion in different ways. Surface dilatation reduces the migration velocity of the droplet, and, for high Bousinessq numbers (
$Bq_{\kappa } \gg 1$), the cross-stream migration of the droplet is greatly suppressed. On the other hand, surface shear increases the migration velocity of the drop, which tends to reduce the time taken by the drop to reach the channel centre. Additionally, we find that increasing the shear Boussinessq number beyond a value of 
$Bq_{\mu } \approx 100$ has no effect on the droplet motion and the curve reaches a plateau of 
$V/V_{h} \approx 1.275$. It is experimentally observed that for higher surface pressures, the surface Boussinessq numbers can attain high values (
$Bq_{\kappa },Bq_{\mu } \gg Ma$) (Bhamla et al. Reference Bhamla, Chai, Alvarez-Valenzuela, Tajuelo and Fuller2017), which can affect the droplet dynamics considerably as seen from figure 3. Experimentally measured values of surface shear viscosities are of the order of magnitude of 
${\sim }10^{-3}\ \text {g}\ \text {s}^{-1}$ (Gupta & Wasan Reference Gupta and Wasan1974), while the surface dilatational viscosities are usually higher. In particular, for a number of surfactant solutions (e.g. dodecanediol, lauryl alcohol), it was found that the surface dilatational viscosities are of the order of magnitude of 
${\sim } 10^{3}\ \text {g}\ \text {s}^{-1}$ (van Voorst Vader, Erkens & Van den Tempel Reference van Voorst Vader, Erkens and Van den Tempel1964). Assuming a surfactant solution in water with a micrometre-size droplet, we obtain the estimates for the surface shear and dilatational Bousinessq numbers to be 
$Bq_{\mu }\sim 10^{3}$ and 
$Bq_{\kappa } \sim 10^{9}$, respectively. The Marangoni number (
$Ma$) can be expressed as the ratio of the Gibbs elasticity parameter (
$0<\beta <1$) and the capillary number (
$Ca \ll 1$). We thus expect typical values of Marangoni numbers to be high (
$Ma \sim 10\text {--}1000$). For example, Wang et al. (Reference Wang, Guo, Dabiri, Vlachos and Ardekani2017) performed experiments with surfactant-laden bubbles with a Marangoni number of 
${\sim }30$. From the above estimates and figure 3, we find that, for commonly used surfactant solutions, the effect of the dilatation viscosity is to completely inhibit the cross-stream motion of the droplet, while the surface shear viscosity enhances the migration velocity by a factor of 
${\sim } 1.275$.
Variation of the drop cross-stream migration velocity with (a) Boussinesq number for surface dilatation for 
$Bq_{\mu }=0,\lambda =1$ and (b) Bousinessq number for surface shear for 
$Bq_{\kappa }=0,\lambda =1$.
For experiments involving droplet migration in microchannels, the lift forces acting on the drops due to the fluid inertia and the droplet deformation are considered to be important. We seek to analyse the relative magnitudes of inertial, deformation and surfactant induced lift forces acting on a droplet in a microchannel experiment. Asmolov (Reference Asmolov1999) derived an expression for the inertial lift force acting on particles suspended in a Poiseuille flow. In particular, their analysis is valid for small particle Reynolds numbers denoted as 
$R_{p}$. The theoretical expression for the deformation induced lift forces was given by Chan & Leal (Reference Chan and Leal1979). For small droplets located away from the channel walls (
$a/R=b/R \sim 0.1$), surface Péclet number 
$Pe_{s} \sim 0.1$, viscosity ratio 
$\lambda =1$, particle Reynolds number 
$R_{p} \sim 10^{-3}-1$ and a capillary number 
$Ca \sim 10^{-3}-0.1$, we find the relative magnitudes of the surfactant induced lift force (
$F^{s}$), inertial induced lift force (
$F^{i}$) and deformation induced lift force (
$F^{d}$) to be as follows:
\begin{align}\frac{F^{s}}{F^{i}} \sim 2.6 \times 10^{-6} Ma-2.6 \times 10^{-3}Ma,\quad \frac{F^{s}}{F^{d}} \sim 5 \times 10^{-4} Ma-5 \times 10^{-2}Ma. \end{align}
We thus find that, for large Marangoni numbers, 
$Ma \sim 10^{2}-10^{3}$, the surfactant induced lift forces would be comparable to the inertial and the deformation induced lift forces for small particle Reynolds numbers and capillary numbers.
We can qualitatively understand the influence of the surface viscosities on the cross-stream migration velocity by analysing the surfactant distribution on the droplet surface at this order (refer to figure 4). For the case of a large surface dilatation viscosity (
$Bq_{\kappa }=100$), we observe that the surfactant redistribution is greatly suppressed as the surface concentration appears to be uniform. We expect that a suppression in the surfactant redistribution decreases the strength of the Marangoni flow and, consequently, the cross-stream migration velocity is reduced significantly as seen from figure 3(a). On the other hand, a large surface shear viscosity causes the surfactant to be redistributed, but in a different manner as compared to a non-viscous interface. In particular, the minimum surfactant concentration (shown in yellow-orange) is observed near a region located further away from the local 
$z$-axis passing through the sphere centre which is parallel to the flow direction. We expect this redistribution to generate a stronger Marangoni flow in the cross-stream direction (along the 
$x$-axis), which consequently increases the cross-stream migration velocity of the droplet as seen from figure 3(b).
Second-order surfactant distribution 
$\Gamma ^{(2)}$ for (a) 
$Bq_{\kappa }=0,Bq_{\mu }=0$, (b) 
$Bq_{\kappa }=100,Bq_{\mu }=0$ and (c) 
$Bq_{\kappa }=0,Bq_{\mu }=100$. In all the cases, 
$a/R=1, b/R=1, \lambda =1,Ma=1$ and the external flow is along the 
$z$-axis.
Our work is applicable in the limit of small Reynolds number (
$Re \ll 1$) based on the particle size, small surface Péclet number (
$Pe_{s}\ll 1$) and a small capillary number (
$Ca \ll 1$). In addition, we assume the surfactant concentration to be negligible (
$\Gamma _{eq}\ll 1$), i.e. we neglect any intermolecular interactions between the surfactant molecules which are adsorbed on the droplet surface. As the cross-stream migration occurs as a second-order effect in the surface Péclet number, the effect of the droplet capillary number (i.e. drop deformation) might become important even if the value of the capillary number is small. Hence, for our results to be applicable in practical scenarios, the neglect of the drop deformation is justified if 
$Ca \ll Pe_{s}^{2}$. For surfactant-laden droplets of size 
${\sim } O(\mathrm {\mu } \text {m})$ immersed in water with an average flow speed 
${\sim } 10^{-3}\ \text {m}\ \text {s}^{-1}$, the Reynolds number is 
${\sim }10^{-3}$. Assuming a surface tension of 
${\sim }10^{-3}\ \text {kg}\ \text {s}^{-2}$, the capillary number can be estimated to be 
${\sim } O(10^{-3})$. These requirements have been satisfied in a number of experiments in the literature, including the comprehensive analysis by Stan et al. (Reference Stan, Ellerbee, Guglielmini, Stone and Whitesides2013) which estimated the lift forces acting on a droplet in a microchannel for different pairs of continuous and dispersed phase fluids. The condition of small Péclet number 
$(Pe_{s} \sim 0.1)$ would be satisfied for surfactant diffusivity values which are 
${\sim }10^{-8}\ \text {m}^{2}\ \text {s}^{-1}$. This requirement seems reasonable as typical values of surfactant diffusivity are expected to lie between 
${\sim } 10^{-8}\ \text {m}^{2}\ \text {s}^{-1}$ and 
${\sim }10^{-10}\ \text {m}^{2}\ \text {s}^{-1}$ (van Voorst Vader et al. Reference van Voorst Vader, Erkens and Van den Tempel1964; Das et al. Reference Das, Mandal and Chakraborty2017).
4. Conclusion
We predicted the migration velocity of a surfactant covered droplet in a Poiseuille flow including the interfacial viscosities of the drop. We performed a regular perturbation expansion for low surface Péclet numbers (
$Pe_{s}$) and solved the problem up to the second order in 
$Pe_{s}$. We represented the interface as a two-dimensional homogeneous fluid using the Bousinessq–Scriven constitutive equation and employed Lamb's general solution for representing the velocity fields both inside and outside the droplet. We obtained an expression of the cross-stream migration velocity of the droplet, where the interfacial effects were captured by Boussinessq numbers for surface shear and surface dilatation (
$Bq_{\mu },Bq_{\kappa }$).
We found that, even when interfacial effects are included, the cross-stream migration velocity of the droplet is zero up to the first order in the surface Péclet number and the first non-zero term only appears at the second order. Surface dilatation viscosity restricts the surfactant redistribution on the droplet surface which suppresses the Marangoni flow around the droplet. This causes the cross-stream velocity to decrease with the dilatation viscosity. For a large dilatation viscosity (
$Bq_{\kappa } > 100$), the cross-stream migration almost vanishes completely. On the other hand, surface shear viscosity redistributes more surfactant along the cross-stream axis of the sphere, which drives a stronger Marangoni flow in the cross-stream direction. Consequently, the migration velocity of the drop increases with the surface dilatation viscosity. However, increasing the shear Bousinessq number beyond a certain value (
$Bq_{\mu } \approx 100$) has no effect on the migration velocity and the maximum velocity enhancement is 
$\approx 1.275$ times more than the cross-stream velocity of a droplet with zero interfacial viscosity. We found that only the surface dilatation viscosity affects the longitudinal speed of the droplet. In particular, the dilatation viscosity reduces the lag experienced by the drop due to the Marangoni flows from the droplet rear to the front. Additionally, a high surface shear viscosity suppresses the compression of the streamlines in the vicinity of the droplet. From our analysis, we emphasize that inclusion of interfacial viscosities is important to accurately model the migration characteristics of droplets for experiments involving droplet manipulation.
Acknowledgements
This publication was made possible due to support from NSF grant CBET-1705371. R.D. thanks F. N. Andrews Fellowship for financial support.
Declaration of interests
The authors report no conflict of interest.
Appendix A
Spherical harmonics in Lamb's general solution for flow inside the droplet:
\begin{gather} p_{n}=\lambda r^{n}\sum_{m=0}^{n}\left(A_{n,m}^{(q)}\cos (m\phi)+ \hat{A}_{n,m}^{(q)}\sin (m\phi)\right)P_{n,m}(\cos \theta), \end{gather}
\begin{gather}\phi_{n}=r^{n}\sum_{m=0}^{n}\left(B_{n,m}^{(q)}\cos (m\phi)+\hat{B}_{n,m}^{(q)} \sin (m\phi)\right)P_{n,m}(\cos \theta), \end{gather}
\begin{gather}\chi_{n}=r^{n}\sum_{m=0}^{n}\left(C_{n,m}^{(q)}\cos (m\phi)+ \hat{C}_{n,m}^{(q)}\sin (m\phi)\right)P_{n,m}(\cos \theta). \end{gather}
where, 
$P_{n,m}(\cos \theta )$ are associated Legendre polynomials of degree 
$n$ and order 
$m$.
Spherical harmonics in Lamb's general solution for flow outside the droplet:
\begin{gather} p_{-n-1}=r^{-n-1}\sum_{m=0}^{n}\left(A_{-n-1,m}^{(q)}\cos (m\phi)+ \hat{A}_{-n-1,m}^{(q)}\sin (m\phi)\right)P_{n,m}(\cos \theta), \end{gather}
\begin{gather}\phi_{-n-1}=r^{-n-1}\sum_{m=0}^{n}\left(B_{-n-1,m}^{(q)}\cos (m\phi)+ \hat{B}_{-n-1,m}^{(q)}\sin (m\phi)\right)P_{n,m}(\cos \theta), \end{gather}
\begin{gather}\chi_{-n-1}=r^{-n-1}\sum_{m=0}^{n}\left(C_{-n-1,m}^{(q)}\cos (m\phi)+ \hat{C}_{-n-1,m}^{(q)}\sin (m\phi)\right)P_{n,m}(\cos \theta). \end{gather}Appendix B
Leading order velocity and pressure fields inside and outside the droplet:
\begin{align} \boldsymbol{u}_{i,r}^{(0)} &= \frac{a(r^{2}-1)}{8R^{2}}\left(-\frac{16a\cos \theta}{2+2Bq_{\kappa}+3\lambda}+ \frac{7ar^{2}(3\cos \theta+5\cos 3\theta)}{7+12Bq_{\kappa}+10Bq_{\mu}+7\lambda} \right.\nonumber\\ &\quad \left.-\frac{60br\cos \phi\sin 2\theta}{5+6Bq_{\kappa}+4Bq_{\mu}+5\lambda}\right). \end{align}
\begin{align} \boldsymbol{u}_{i,\theta}^{(0)} &= \frac{a}{8R^{2}}\left(8br\cos \phi-\frac{20br(-3+5r^{2})\cos 2\theta \cos \phi}{5+6Bq_{\kappa}+4Bq_{\mu}+5\lambda} \right.\nonumber\\ &\quad \left.+\frac{16a(-1+2r^{2})\sin \theta}{2+2Bq_{\kappa}+3\lambda}- \frac{7br^{2}(-2+3r^{2})(3+5\cos 2\theta)\sin \theta}{7+12Bq_{\kappa}+10Bq_{\mu}+7\lambda}\right). \end{align}
\begin{align} \boldsymbol{u}_{i,\phi}^{(0)} &= \frac{abr(-12Bq_{\kappa}-8Bq_{\mu}+5(-5+5r^{2}-2\lambda))\cos \theta \sin \phi}{2R^{2}(5+6Bq_{\kappa}+4Bq_{\mu}+5\lambda)}. \end{align}
\begin{align} & \boldsymbol{u}_{e,r}^{(0)} = \frac{a\cos \theta}{60R^{2}}\left(-\frac{16a(-3+2Bq_{\kappa}+3\lambda)}{r^{3} (2+2Bq_{\kappa}+3\lambda)}+40a\frac{(2Bq_{\kappa}+3\lambda)}{2+2Bq_{\kappa}+3\lambda} \right.\nonumber\\ &\quad -\frac{3a(12Bq_{\kappa}(-5+7r^{2})+10Bq_{\mu}(-5+7r^{2})-35\lambda+7r^{2}(2+7\lambda)) (-1+5\cos 2\theta)}{r^{5}(7+12Bq_{\kappa}+10Bq_{\mu}+7\lambda)} \nonumber\\ &\quad -120br\cos \phi\sin \theta \nonumber\\ &\quad +\frac{60b(-6(3Bq_{\kappa}+2Bq_{\mu})+10(-1+3Bq_{\kappa}+2Bq_{\mu})r^{2}+5\lambda(-3+5r^{2}))\cos \phi \sin \theta}{r^{4}(5+6Bq_{\kappa}+4Bq_{\mu}+5\lambda)} \nonumber\\ &\quad \left.\vphantom{\frac{60b(-6(3Bq_{\kappa}+2Bq_{\mu})+10(-1+3Bq_{\kappa}+2Bq_{\mu})r^{2}+5\lambda(-3+5r^{2}))\cos \phi \sin \theta}{r^{4}}}-60ar^{2}\sin ^{2}\theta\right). \end{align}
\begin{align} \boldsymbol{u}_{e,\theta}^{(0)} &= \frac{a}{48r^{5}R^{2}}\left(3\cos 2\theta \left(\frac{16br(6Bq_{\kappa}+4Bq_{\mu}+5\lambda)\cos \phi}{5+6Bq_{\kappa}+4Bq_{\mu}+5\lambda} \right.\right. \nonumber\\ &\quad \left.-\frac{a(-30(6Bq_{\kappa}+5Bq_{\mu})+14(1+6Bq_{\kappa}+5Bq_{\mu})r^{2}+7\lambda(-15+7r^{2}) \sin \theta)}{7+12Bq_{\kappa}+10Bq_{\mu}+7\lambda}\right) \nonumber\\ &\quad +\sin \theta\left(a\left(27-19r^{2}-32r^{5}+\frac{32(r^{2}+2r^{5})}{2+2Bq_{\kappa}+3\lambda}+ \frac{63(-3+r^{2})}{7+12Bq_{\kappa}+10Bq_{\mu}+7\lambda}\right) \right. \nonumber\\ &\quad +\left.\left.\vphantom{\frac{60}{70}}48r^{6}\sin \theta(2b\cos \phi+ar\sin \phi)\right)\right). \end{align}
\begin{align} \boldsymbol{u}_{e,\phi}^{(0)} = -\frac{ab(6Bq_{\kappa}+4Bq_{\mu}+5\lambda)\cos \theta \sin \phi}{r^{4}R^{2}(5+6Bq_{\kappa}+4Bq_{\mu}+5\lambda)}. \end{align}
\begin{align} p_{i}^{0} &= \frac{ar\lambda\cos \theta}{R^{2}}\left(-\frac{20a}{2+2Bq_{\kappa}+3\lambda}+ \frac{21ar^{2}(-3+5\cos ^{2}\theta)}{7+12Bq_{\kappa}+10Bq_{\mu}+7\lambda} \right.\nonumber\\ &\quad \left.-\frac{105br\cos \phi\sin \theta}{5+6Bq_{\kappa}+4Bq_{\mu}+5\lambda}\right). \end{align}
\begin{align} p_{e}^{0} &= \frac{a}{4r^{4}R^{2}}\left(-\frac{7a(2+12Bq_{\kappa}+10Bq_{\mu}+7\lambda) \cos \theta(-3+5\cos ^{2}\theta)}{7+12Bq_{\kappa}+10Bq_{\mu}+7\lambda} \right.\nonumber\\ &\quad \left.+\frac{20br(2+6Bq_{\kappa}+4Bq_{\mu}+5\lambda)\cos \phi \sin 2\theta}{5+6Bq_{\kappa}+4Bq_{\mu}+5\lambda}\right) \end{align}First-order velocity and pressure fields inside and outside the droplet:
\begin{align} \boldsymbol{u}_{i,r}^{(1)} &= \frac{aMa(r^{2}-1)}{8R^{2}}\left(\frac{16a\cos \theta}{(2+2Bq_{\kappa}+3\lambda)^{2}}- \frac{7ar^{2}(3\cos \theta+5\cos 3\theta)}{(7+12Bq_{\kappa}+10Bq_{\mu}+7\lambda)^{2}} \right.\nonumber\\ &\quad \left.+\frac{60br\cos \phi\sin 2\theta}{(5+6Bq_{\kappa}+4Bq_{\mu}+5\lambda)^{2}}\right). \end{align}
\begin{align} \boldsymbol{u}_{i,\theta}^{(1)} &= \frac{aMa}{8R^{2}}\left(\frac{20br(-3+5r^{2})\cos 2\theta \cos \phi}{(5+6Bq_{\kappa}+4Bq_{\mu}+5\lambda)^{2}}+\frac{16a(1-2r^{2})\sin \theta}{(2+2Bq_{\kappa}+3\lambda)^{2}} \right.\nonumber\\ &\quad \left.+\frac{7ar^{2}(-2+3r^{2})(3+5\cos 2\theta)\sin \theta}{(7+12Bq_{\kappa}+10Bq_{\mu}+7\lambda)^{2}}\right). \end{align}
\begin{equation} \boldsymbol{u}_{i,\phi}^{(1)} = -\frac{5abMar(-3+5r^{2})\cos \theta \sin \phi}{2R^{2}(5+6Bq_{\kappa}+4Bq_{\mu}+5\lambda)^{2}}. \end{equation}
\begin{align} \boldsymbol{u}_{e,r}^{(1)} &= \frac{aMa}{6R^{2}}\left(\frac{8a\cos \theta}{(2+2Bq_{\kappa}+3\lambda)^{2}}- \frac{8a\cos \theta}{r^{3}(2+2Bq_{\kappa}+3\lambda)^{2}} \right.\nonumber\\ &\quad \left. -\frac{21a(-1+r^{2})\cos \theta (-3+5\cos ^{2}\theta)}{r^{5}(7+12Bq_{\kappa}+10Bq_{\mu}+7\lambda)^{2}}+\frac{45b(-1+r^{2})\cos \phi \sin 2\theta}{r^{4}(5+6Bq_{\kappa}+4Bq_{\mu}+5\lambda)^{2}}\right). \end{align}
\begin{align} \boldsymbol{u}_{e,\theta}^{(1)} &= \frac{aMa}{48R^{2}}\left(-\frac{64a\sin \theta}{(2+2Bq_{\kappa}+3\lambda)^{2}}- \frac{32a\sin \theta}{r^{3}(2+2Bq_{\kappa}+3\lambda)^{2}} \right.\nonumber\\ &\quad \left.-\frac{21a(-3+r^{2})\sin \theta(3+5\cos ^{2}\theta)}{r^{5}(7+12Bq_{\kappa}+10Bq_{\mu}+7\lambda)^{2}}+ \frac{240b\cos \phi\cos 2\theta}{r^{4}(5+6Bq_{\kappa}+4Bq_{\mu}+5\lambda)^{2}}\right). \end{align}
\begin{equation} \boldsymbol{u}_{e,\phi}^{(1)}=-\frac{5abMa\cos \theta \sin \phi}{r^{4}R^{2}(5+6Bq_{\kappa}+4Bq_{\mu}+5\lambda)^{2}}. \end{equation}
\begin{align} p_{i}^{(1)} &= \frac{aMar\lambda\cos \theta}{R^{2}}\left(\frac{20a}{(2+2Bq_{\kappa}+3\lambda)^{2}}- \frac{21a(-r^{2})(-3+5\cos ^{2}\theta)}{(7+12Bq_{\kappa}+10Bq_{\mu}+7\lambda)^{2}} \right.\nonumber\\ &\quad \left.+\frac{105br\cos \phi\sin \theta}{(5+6Bq_{\kappa}+4Bq_{\mu}+5\lambda)^{2}}\right). \end{align}
\begin{equation} p_{e}^{(1)}=\frac{5aMa}{4r^{4}R^{2}}\left(-\frac{7a\cos \theta(-3+5\cos ^{2}\theta)}{(7+12Bq_{\kappa}+10Bq_{\mu}+7\lambda)^{2}}+ \frac{12br\cos \phi\sin 2\theta}{(5+6Bq_{\kappa}+4Bq_{\mu}+5\lambda)^{2}}\right). \end{equation}Second-order surfactant concentration distribution:
\begin{align} \Gamma^{(2)} &= P_{1}^{0}(\cos \theta) \left(\frac{2a^{2}Ma}{R^{2}(2+2Bq_{\kappa}+3\lambda)^{2}}\right) \nonumber\\ &\quad +P_{2}^{0}(\cos \theta) \frac{a^{2}}{378R^{4}} \left(\frac{45b^{2}(40+42Bq_{\kappa}+28Bq_{\mu}+35\lambda)}{(5+6Bq_{\kappa}+4Bq_{\mu}+5\lambda)^{2}}+ \frac{504a^{2}}{(2+2Bq_{\kappa}+3\lambda)^{2}} \right.\nonumber\\ &\quad \left.+\frac{49a^{2}}{(7+12Bq_{\kappa}+10Bq_{\mu}+7\lambda)^{2}}- \frac{378a^{2}}{(2+2Bq_{\kappa}+3\lambda)(7+12Bq_{\kappa}+10Bq_{\mu}+7\lambda)}\right) \nonumber\\ &\quad +P_{3}^{0}(\cos \theta)\left(-\frac{7a^{2}Ma}{6R^{2}(7+12Bq_{\kappa}+10Bq_{\mu}+7\lambda)^{2}}\right) \nonumber\\ &\quad +P_{4}^{0}(\cos \theta)\frac{a^{2}}{924R^{4}}\left(-\frac{1320b^{2}}{(5+6Bq_{\kappa}+4Bq_{\mu}+5\lambda)^{2}} \right.\nonumber\\ &\quad \left.-\frac{7a^{2}(1190+2070Bq_{\kappa}+1760Bq_{\mu}+1169\lambda)}{(2+2Bq_{\kappa}+3\lambda) (7+12Bq_{\kappa}+10Bq_{\mu}+7\lambda)^{2}}\right) \nonumber\\ &\quad +P_{6}^{0}(\cos \theta)\left(\frac{175a^{4}}{594R^{4}(7+12Bq_{\kappa}+10Bq_{\mu}+7\lambda)^{2}}\right) +P_{1}^{1}(\cos \theta)\cos \phi \nonumber\\ &\quad \times\left(\frac{a^{3}b(40+144Bq_{\kappa}^{2}) +80Bq_{\mu}^{2}+109\lambda+70\lambda^{2}}{2R^{4}(2+2Bq_{\kappa}+3\lambda) (5+6Bq_{\kappa}+4Bq_{\mu}+5\lambda)(7+12Bq_{\kappa}+10Bq_{\mu}+7\lambda)} \right.\nonumber\\ &\quad \left.+\frac{4Bq_{\mu}(29+39\lambda)+2Bq_{\kappa} (77+108Bq_{\mu}+102\lambda)}{R^{4}(2+2Bq_{\kappa}+3\lambda)(5+6Bq_{\kappa}+4Bq_{\mu}+5\lambda) (7+12Bq_{\kappa}+10Bq_{\mu}+7\lambda)}\right) \nonumber\\ &\quad +P_{2}^{1}(\cos \theta)\cos \phi \left(-\frac{5abMa}{3R^{2}(5+6Bq_{\kappa}+4Bq_{\mu}+5\lambda)^{2}}\right) +P_{3}^{1}(\cos \theta)\cos \phi \nonumber\\ &\quad \times\left(\frac{a^{3}b(2170+252Bq_{\kappa}^{2})+2457\lambda+315\lambda^{2}}{216R^{4}(2+2Bq_{\kappa}+3\lambda) (5+6Bq_{\kappa}+4Bq_{\mu}+5\lambda)(7+12Bq_{\kappa}+10Bq_{\mu}+7\lambda)} \right. \nonumber\\ &\quad \left.+\frac{12Bq_{\mu}(254+21\lambda)+ 2Bq_{\kappa}(1931+84Bq_{\mu}+294\lambda)}{216R^{4}(2+2Bq_{\kappa}+3\lambda)(5+6Bq_{\kappa}+4Bq_{\mu}+5\lambda) (7+12Bq_{\kappa}+10Bq_{\mu}+7\lambda)}\right) \nonumber\\ &\quad +P_{5}^{1}(\cos \theta)\cos \phi \left(-\frac{10a^{3}b}{27R^{4}(5+6Bq_{\kappa}+4Bq_{\mu}+5\lambda)(7+12Bq_{\kappa}+10Bq_{\mu}+7\lambda)}\right) \nonumber\\ &\quad +P_{2}^{2}(\cos \theta)\cos 2\phi \left(-\frac{5a^{2}b^{2}(20+42Bq_{\kappa}+28Bq_{\mu}+35\lambda)}{252R^{4} (5+6Bq_{\kappa}+4Bq_{\mu}+5\lambda)^{2}}\right) \nonumber\\ &\quad +P_{4}^{2}(\cos \theta)\cos 4\phi \left(\frac{5a^{2}b^{2}}{42R^{4}(5+6Bq_{\kappa}+4Bq_{\mu}+5\lambda)^{2}}\right). \end{align}
