1. Introduction
Micron and submicron scale active particles have been investigated for decades. At these length scales, the active media or ‘swimmers’ are typically in the low-Reynolds-number regime where inertial effects of the encompassing fluid are negligible. Unlike most simply periodic mechanical motions that are adopted by macroscopic active matter for self-propulsion, the swimming mechanisms governing the propulsion of such swimmers are often more complex due to the restrictions that arise with the absence of inertia (Purcell's ‘scallop theorem’ (Purcell Reference Purcell1977)). Apart from the natural mechanisms found in small-scale living organisms, self-propulsion of artificially designed active matter or ‘colloidal motors’ is also of great interest (Sundararajan et al. Reference Sundararajan, Lammert, Zudans, Crespi and Sen2008; Solovev et al. Reference Solovev, Xi, Gracias, Harazim, Deneke, Sanchez and Schmidt2012; Wang & Gao Reference Wang and Gao2012; Abdelmohsen et al. Reference Abdelmohsen, Peng, Tu and Wilson2014) and has been achieved through mechanisms including ‘bubble propulsion’ (Gibbs & Zhao Reference Gibbs and Zhao2009), ‘neutral diffusiophoresis’ (Golestanian, Liverpool & Ajdari Reference Golestanian, Liverpool and Ajdari2005) and ‘self-electrophoresis’ (Wang et al. Reference Wang, Hernandez, Bartlett, Bingham, Kline, Sen and Mallouk2006). In particular, diffusiophoresis is a commonly cited mechanism wherein a particle is placed in a gradient of another solute species which then interacts through attractive or repulsive forces with the particle leading to a net imbalance in the forces acting on the particle (Anderson Reference Anderson1989; Sharifi-Mood, Koplik & Maldarelli Reference Sharifi-Mood, Koplik and Maldarelli2013). The particle is thus forced to translate by this force which is then balanced by the hydrodynamic drag from the surrounding fluid. Janus particles have been favoured in the past to design motors based on this mechanism due to their non-homogeneous surface properties that give them the ability to interact with the environment in an asymmetric manner. Janus particles can be designed with varying surface properties so that a portion of the particle (often labelled the ‘active side’) can itself produce or selectively catalyse the production of the surrounding solute species thereby creating and maintaining the required asymmetric concentration distribution of the solute which in turn drives the Janus particle (Golestanian et al. Reference Golestanian, Liverpool and Ajdari2005; Pawar & Kretzschmar Reference Pawar and Kretzschmar2010; Popescu, Uspal & Dietrich Reference Popescu, Uspal and Dietrich2016). Realizations of this mechanism have revealed that the resulting propulsion speeds are generally quite low with magnitudes of the order of 1–10 body lengths per second which typically translates to speeds of the order of 
$10~{\rm \mu}$m s
$^{-1}$ (Howse et al. Reference Howse, Jones, Ryan, Gough, Vafabakhsh and Golestanian2007; Ebbens & Howse Reference Ebbens and Howse2010; Wang & Wu Reference Wang and Wu2014). This is because the forces involved are mostly weak van der Waals forces (for neutral solutes). This is one of the key limitations of this form of colloidal motor.
Surface tension forces can often play a dominant role in the dynamics of fluids when liquid interfaces are involved. Self-propulsion using capillary forces has also been studied. For example, Lord Rayleigh (Strutt Reference Strutt1890) attempted to experimentally quantify the sensitivity of grease in arresting the spontaneous motion of camphor scrapings floating atop a bath of water. He describes some of the earliest efforts that helped cement the understanding of some of the relevant physics, including the importance of surface tension differences in driving these scrapings.
The term ‘Marangoni surfer’ is commonly used to refer to particles located at or near the junction of two fluids that are propelled along the fluid–fluid interface by generating a gradient in the surface tension forces surrounding them. Camphor boats are a classic example of such surfers where the boat is some particle partially immersed in a liquid interface loaded with camphor (a surface-active agent). As the camphor dissolves, it reduces the local tension at the fluid interface. The particle can generate a surface tension imbalance around it by releasing the camphor asymmetrically which causes the propulsion of the boat. This simple demonstration of the concept of Marangoni surfing using everyday items has evoked much academic inquiry in the subject matter. There has been much work devoted to understanding the single particle and collective dynamics of particles coated with camphor or other such surfactants at a liquid interface (Nagayama et al. (Reference Nagayama, Nakata, Doi and Hayashima2004) and references therein).
Janus particles have also been used to design propulsion mechanisms involving surface tension gradients or Marangoni forces. Metallic Janus particles can also be preferentially heated to produce a temperature gradient around the particle (Würger Reference Würger2014; Dietrich et al. Reference Dietrich, Jaensson, Buttinoni, Volpe and Isa2020). This results in a similar reduction in the interfacial tension and the resulting imbalance of forces (thermocapillary forces) is also a possible driving force for such surfers. Such surface tension-based propulsion mechanisms have proved to be quite effective with the resulting velocities reaching scales of 
${\sim }O$(cm s
$^{-1}$) which translates to 1000–10 000 body lengths per second for micron-scale colloidal motors (Sur, Masoud & Rothstein Reference Sur, Masoud and Rothstein2019; Dietrich et al. Reference Dietrich, Jaensson, Buttinoni, Volpe and Isa2020)!
Early modelling efforts used a mix of analytical and semiempirical methods in an attempt to capture the interesting dynamics of the particles without fully resolving the flow field surrounding them (Nakata et al. Reference Nakata, Iguchi, Ose, Kuboyama, Ishii and Yoshikawa1997; Nagayama et al. Reference Nagayama, Nakata, Doi and Hayashima2004; Heisler et al. Reference Heisler, Suematsu, Awazu and Nishimori2012), etc. More recently, Lauga & Davis (Reference Lauga and Davis2012) approached the problem by solving the Stokes equations for the edgewise translation of a thin cylindrical disk with a prescribed, asymmetric concentration of some surfactant along the circular edge of the disk. They solved the surfactant transport equations at the diffusive limit with negligible surface Péclet number 
$(Pe^s \to 0)$ while assuming that the surfactant is insoluble in the bulk volume of the fluid. However, commonly studied Marangoni propulsion problems require that the propelled body, such as the Janus particle, either be at or near a fluid–fluid interface (Domínguez et al. Reference Domínguez, Malgaretti, Popescu and Dietrich2016; Gidituri, Panchagnula & Pototsky Reference Gidituri, Panchagnula and Pototsky2019) which make it difficult for this mechanism to be adapted for applications involving propulsion within the bulk of the fluid.
Another well studied class of low-Reynolds-number propulsion problems involve the motion of individual bubbles and drops or a collection of them in the bulk of the fluid. Here the earlier issue is sidestepped as the species being propelled carries with it a fluid interface. One of the earliest fully analytical solutions in this class of problems is the well-known solution to a spherical drop translating with a constant velocity under an applied temperature gradient (Young, Goldstein & Block Reference Young, Goldstein and Block1959) where the heat transfer is assumed to be purely diffusive (low Péclet number). The solutions also assume a low capillary number regime to retain the spherical shape of the drop. Other interesting developments include solutions to pairs of drops interacting with an externally applied temperature gradient (Wang, Mauri & Acrivos Reference Wang, Mauri and Acrivos1994), drops translating under a temperature/concentration gradient that are affected by the adsorption of surfactants at the interface (Kim & Subramanian Reference Kim and Subramanian1989; Chen & Stebe Reference Chen and Stebe1997) and more. While these types of solutions have been put to practical use in various engineering applications, they have not been used to design active particles.
In the above examples, it can be noted that the propulsion of the drop or particles depends on its ability to sample a temperature or solute gradient around it. This is because, in the quasistatic diffusion dominated regime of heat/mass transport, it would not be possible to spontaneously generate thrust forces if the particle were to only interact symmetrically with its surroundings. When considering particles and drops that are able to generate a symmetric distribution of some interacting species such as surfactant molecules or other solutes that can interact with the particle or drop surface via, for example, van der Waals forces, the net forces generated from these interactions typically cancel and the particle or drop remains quiescent. As mentioned above, a simple way to generate a gradient around the motor is to use an active Janus particle or a similar anisotropic active particle (Schnitzer & Yariv Reference Schnitzer and Yariv2015). In fact, if the geometry around the active particle is non-isotropic, it allows for asymmetric distributions of the solute or heat around the particle due to diffusion. Thus, in such cases self-propulsion can still be achieved and even guided by the boundaries around the active particle. (Popescu, Dietrich & Oshanin Reference Popescu, Dietrich and Oshanin2009; Michelin & Lauga Reference Michelin and Lauga2015; Domínguez et al. Reference Domínguez, Malgaretti, Popescu and Dietrich2016; Mozaffari et al. Reference Mozaffari, Sharifi-Mood, Koplik and Maldarelli2016; Uspal et al. Reference Uspal, Popescu, Dietrich and Tasinkevych2016; Yariv Reference Yariv2016; Popescu et al. Reference Popescu, Uspal, Dominguez and Dietrich2018).
However, spontaneous motion can be achieved if the Péclet number for the mass transfer around the motor is sufficiently high. It is well known that camphor particles that are symmetric in shape can also propel themselves at an interface under the appropriate conditions (Kitahata et al. Reference Kitahata, Hiromatsu, Doi, Nakata and Islam2004; Nagayama et al. Reference Nagayama, Nakata, Doi and Hayashima2004). The mechanism behind this spontaneous propulsion has also been theoretically investigated in analogous systems such as a particle or drop in the bulk of the fluid exhibiting spontaneous propulsion driven by phoretic or Marangoni forces (Michelin, Lauga & Bartolo Reference Michelin, Lauga and Bartolo2013; Izri et al. Reference Izri, Van Der Linden, Michelin and Dauchot2014; Michelin & Lauga Reference Michelin and Lauga2014). The advective coupling of the concentration field with the velocity field implies that small perturbations to the velocity field can lead to a deviation from the symmetric distribution of the concentration of interacting species along the surface of the particle or drop. These perturbations in turn generate additional thrusts on the particle or drop that are unbalanced leading to its propulsion. As a result, the system can either restore itself to the earlier state of quiescence or deviate further from this state by further propelling the particle. A stability analysis may help determine the right conditions to potentially induce spontaneous propulsion in the system.
To retain the practical uses of a solid particle and to satisfy the need for an interface at the vicinity of the particle, a suitable design combining drops with particles is then sought. In particular, particles or drops encapsulated in liquid drops (‘compound drops or particles’ (Johnson & Sadhal Reference Johnson and Sadhal1985; Chaithanya & Thampi Reference Chaithanya and Thampi2019)) can be used as an alternative to an isolated drop or particle. Their ability to carry the core particle in the presence of gravitational forces is one of the earliest such analytical solutions in the literature (Sadhal & Oguz Reference Sadhal and Oguz1985). This study was further extended to the case of thermocapillary motion under an externally applied temperature gradient (Morton, Subramanian & Balasubramaniam Reference Morton, Subramanian and Balasubramaniam1990). Other related literature can be found in the book by Subramanian, Balasubramaniam & Clark (Reference Subramanian, Balasubramaniam and Clark2002). However, the need for an externally applied gradient falls short of our requirement that the motor be self-propelled. For compound drops, the inner drop can be loaded with surfactants that are released asymmetrically with respect to the outer drop (Tsemakh, Lavrenteva & Nir Reference Tsemakh, Lavrenteva and Nir2004). For compound particles, it is necessary to incorporate one of the earlier mechanisms by which the motor can produce and maintain its own temperature or surfactant gradient.
The focus of this article is then to design a novel active species using the familiar Janus particle at its core and a liquid interface around the particle so that the combined compound drop/Janus particle system can behave as an active motor which employs the more significant Marangoni propulsion force to propel itself in the bulk of the fluid while retaining any benefits of a solid colloidal motor. In recent literature, active compound particles have been investigated where the drop is driven by the actuation of the fluid inside using different mechanisms including for example an anisotropic distribution of internal forces, for example an active encapsulated squirmer (Chaithanya & Thampi Reference Chaithanya and Thampi2020; Sprenger et al. Reference Sprenger, Shaik, Ardekani, Lisicki, Mathijssen, Guzmán-Lastra, Löwen, Menzel and Daddi-Moussa-Ider2020; Kree, Rueckert & Zippelius Reference Kree, Rueckert and Zippelius2021; Kree & Zippelius Reference Kree and Zippelius2022) Although these studies also involve the propulsion of a compound particle/drop, they differ fundamentally from the current analysis in that the drop interface is not necessary to generate the propulsion forces unlike in the case of our Marangoni-driven mechanism.
2. Problem statement
Consider a spherical Janus particle of radius 
${{R}_{1}}$ encapsulated by a spherical drop of radius 
${{R}_{2}}$ of fluid 1 in an infinite medium of fluid 2 as depicted in figure 1. The inert portion of the Janus particle is denoted by the surface 
${{S}_{p}}$ and the active spherical cap is denoted by the surface 
${{S}_{c}}$. The interface of the two fluids is labelled 
${{S}_{I}}$. Subscripts 1 and 2 will be used to specify the individual variables in either the drop phase or the infinite medium, respectively. The Janus cap is active in the sense that it generates a flux of some solute species which exhibits surface activity at the interface of the two fluids. This can be achieved by loading the cap with a surfactant, or by allowing the cap to react with the surrounding fluid to produce the desired species as in catalytic active motors. Equivalently, the Janus particle can be preferentially heated using laser-based mechanisms to produce a heat flux from the (typically metallic) conductive cap. Since the surface tension forces are sensitive to the species concentration and temperature gradients at the interface, it is necessary to solve the coupled scalar species or energy and momentum conservation equation in both fluids to completely determine the corresponding concentration (or) temperature and velocity fields.
Figure 1. Illustration of the compound drop/Janus motor generating a flux of a surface-active species.
2.1. Concentration field
2.1.1. Governing equations
Following the approach of (Lauga & Davis Reference Lauga and Davis2012; Masoud & Stone Reference Masoud and Stone2014), we neglect the contributions of advection and restrict ourselves to a regime of purely diffusive transport of surfactants in either fluid 
$( P{{e}_{i}}\rightarrow 0 )$. In doing so, the governing equations reduces to the simple, linear Laplace equation for the species concentration in either fluid 
${{C}_{i}}( x )$ (or analogously the temperature field for heat transport) (Leal Reference Leal2007),
\begin{equation} {{\nabla}^{2}}{{C}_{i}}=0.\end{equation}
The fluid motion is still driven by the gradients of the surface tension 
$\gamma$ across the drop interface which requires a one-way coupling of the equations at the minimum.
2.1.2. Boundary conditions
We assume that the Janus particle produces a constant flux 
$J$ of the surface-active species from the active cap region to approximate the transfer of the solute species from the Janus particle surface into the surrounding liquid. This assumption is valid if the variations in the bulk concentration of the solute near the Janus particle surface have negligible consequences on the solute efflux rate (i.e. when the diffusion of the solute from the particle surface is faster than the rate of production of the solute (Córdova-Figueroa, Brady & Shklyaev Reference Córdova-Figueroa, Brady and Shklyaev2013; Shklyaev, Brady & Córdova-Figueroa Reference Shklyaev, Brady and Córdova-Figueroa2014)). Similar assumptions can be found in the literature for Janus particles driven via diffusiophoresis (Sharifi-Mood et al. Reference Sharifi-Mood, Koplik and Maldarelli2013). It can sometimes be instructive to assume a constant value for the concentration/temperature at the cap surface. The latter case is particularly useful in analysing temperature gradients near hot Janus particles with highly conductive metallic caps (Würger Reference Würger2014).
Such a limit is not considered in this study. At the drop interface, we restrict the analysis to solutes which exhibit faster adsorption and desorption kinetics relative to the surface transport of the species. Therefore, the surface excess concentration 
$\varGamma ( x )$ of the solute equilibrates with the bulk concentration instantaneously and the disturbance due to surface transport is negligible. Under dilute conditions, this equilibrium value is linearly proportional to the corresponding bulk concentrations in either fluid. Consequently, we assume that the surface concentrations of the solute species are purely determined by the instantaneous bulk concentration at every point along the interface via the use of partition coefficients 
$( {{k_1},{k_2}} )$ (Berg Reference Berg2010). For a heat-driven compound drop/Janus particle system, the temperature field is continuous at the interface 
$( {{k_1}={k_2}}=1 )$.
The boundary conditions for the concentration field satisfying (2.1) are given as
\begin{gather} - {D_1}\hat{\boldsymbol n} \boldsymbol{\cdot}\boldsymbol{\nabla} {C_1}( x ) = J \{ x \in {S_c}\} , \end{gather}
\begin{gather}\hat{\boldsymbol n} \boldsymbol{\cdot} \boldsymbol{\nabla} {C_1}( x ) = 0 \{ x \in {S_p}\} , \end{gather}
\begin{gather}\varGamma ( x ) = {k_1}{C_1}( x ) = {k_2}{C_2}( x ) \{ x \in {S_I}\} , \end{gather}
\begin{gather}- {D_1}\hat{\boldsymbol n} \boldsymbol{\cdot} \boldsymbol{\nabla} {C_1}( x ) ={-} {D_2}\hat{\boldsymbol n} \boldsymbol{\cdot} \boldsymbol{\nabla} {C_2}( x )\quad \{ x \in {S_I}\}. \end{gather}
Here 
$\hat{\boldsymbol n}$ refers to the normal pointing into the fluid 1 at the prescribed surface and 
${{D}_{i}}$ is the diffusion coefficient of the solute in the respective fluid i.
2.2. Velocity field
We shall limit the analysis to axisymmetric configurations of the compound drop/Janus particle system wherein the instantaneous position and orientation of the Janus particle relative to the drop is such that the axis of the spherical cap is along the line joining the centres of the two spheres. In such configurations the entire compound drop/Janus particle system can be described as a body of revolution. By arguments of symmetry, the motion is purely translatory (parallel to the axis of revolution).
2.2.1. Governing equations
We neglect any effects of shear and dilatational surface viscosities due to surfactant at the surface. The small length scale typical in problems involving microscale colloidal motors means the governing equations for the velocity and pressure fields 
$({{u}_{i}},{{p}_{i}} )$ are well approximated by the Stokes equations 
$({{Re}_{i}}\sim 0 )$ given by (assuming no external body forces or other singularities)
\begin{gather} -\boldsymbol{\nabla} {{p}_{i}}+{{\mu }_{i}}{{\nabla}^{2}}{\boldsymbol{{u}_{i}}}=0, \end{gather}
\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{{u}_{i}}}=0, \end{gather}
where 
${{\mu }_{i}}$ is the dynamic viscosity in the fluid i.
2.2.2. Boundary conditions
Assuming the Janus particle and drop exhibit negligible inertia, the hydrodynamic drag will have to balance out the propulsive forces on either body. Under these force-free conditions, the drop and the Janus particle translate with velocities U, V, respectively. The boundary conditions for the fluid are no slip at the Janus particle surface and the continuity of velocities at the fluid–fluid interface. As previously mentioned, the drop is assumed to retain its spherical shape. This is a fair assumption if the surface tension forces dominate over the viscous forces to restore the drop to the energetically favourable spherical shape 
$( Ca\rightarrow 0 )$. Thus the kinematic boundary condition reduces to the simpler no-penetration condition which restrains the normal velocity at the drop surface. Finally, the tangential stress jump at the interface is balanced by the additional stress generated due to the surface tension gradients (Marangoni stresses) (Leal Reference Leal2007),
\begin{gather} {\boldsymbol{{u}_{1}}}( x )=V{{\hat{\boldsymbol{e}}_{z}}} \{x\in {{S}_{p}}, {{S}_{c}}\}, \end{gather}
\begin{gather}{\boldsymbol{{u}_{1}}}( x )={\boldsymbol{{u}_{2}}}( x ) \{x\in {{S}_{I}}\}, \end{gather}
\begin{gather}{\boldsymbol{{u}_{1}}}( x )\boldsymbol{\cdot} \boldsymbol{\hat{n}}=U{{\hat{\boldsymbol{e}}_{z}}}\boldsymbol{\cdot}\boldsymbol{\hat{n}} \{x\in {{S}_{I}}\}, \end{gather}
\begin{gather}\boldsymbol{\hat{n}}\boldsymbol{\cdot} ( \boldsymbol{{{\sigma }_{2}}}( x )-\boldsymbol{{{\sigma }_{1}}}( x ) )\boldsymbol{\cdot}{{\hat{\boldsymbol{t}}}}={{\hat{\boldsymbol{t}}}}\boldsymbol{\cdot} \boldsymbol{\nabla_s} \gamma ( x ) \{x\in {{S}_{I}}\}, \end{gather}
where 
$( \boldsymbol {\nabla _s}=\boldsymbol {\nabla } -( \boldsymbol {\hat {n}}\boldsymbol {{\cdot }}\boldsymbol {\nabla } )\boldsymbol {\hat {n}} )$ and 
$( {{{\hat {\boldsymbol {t}}}}} )$ represents each of the pair of orthogonal tangential unit vectors at any given point on the surface. Together 
$( \boldsymbol {\hat {n}},{{{\hat {\boldsymbol {t}}}}})$ are a mutually orthogonal set of unit vectors at the prescribed surface and 
$\boldsymbol {{{\sigma }_{i}}}$ is the Cauchy stress tensor for fluid i.
A linear relationship between the surface tension and the concentration field is adopted,
\begin{equation} \gamma ={{\gamma }_{0}}( 1-\beta \varGamma ),\end{equation}
where 
${{\gamma }_{0}}$ is the clean interface surface tension and the coefficient 
$\beta =-({1}/{{{\gamma }_{0}}})({\partial \gamma }/{\partial \varGamma })$ is assumed to have a prescribed constant value. This assumption has been shown to be valid for sufficiently dilute concentrations of the surfactant (Manikantan & Squires Reference Manikantan and Squires2020).
The forces on the drop and the Janus particle directed along the axis of symmetry (say 
${{\hat {\boldsymbol {e}}_{z}}}$) are given by
\begin{gather} \mathop{{\int\!\!\!\!\!\int}\mkern-21mu \bigcirc}_{{S_I}} { {\boldsymbol{{{\hat e}_z}} \boldsymbol{\cdot} {\boldsymbol{\sigma _2}} \boldsymbol{\cdot} \hat{\boldsymbol n}} } \,{\rm d} s = 0, \end{gather}
\begin{gather}\mathop{{\int\!\!\!\!\!\int}\mkern-21mu \bigcirc}_{{S_p} \cup {S_c}} { {\boldsymbol{{{\hat e}_z}} \boldsymbol{\cdot} {\boldsymbol{\sigma _1}} \boldsymbol{\cdot} \hat{\boldsymbol n}} } \,{\rm d} s = 0. \end{gather}
For convenience, we can denote the relevant ratios of the solute partition and diffusion coefficients and the fluid viscosities as 
${{{k_2}}/{{k_1}}} = k, {{{\mu _2}}/{{\mu _1}}} = \mu, {{{D_2}}/{{D_1}}} = D$. Since the governing equations for both the fluid mechanics and the concentration field are linear and quasisteady, it is possible to solve the problem analytically for various realizations of the compound drop/Janus particle system involving different relative sizes, positions, cap sizes, etc. The relevant non-dimensional groups are defined as 
$P{e_i} = {{U{R_1}}/{{D_i}}}$, 
$Ca = {{{\mu _1}U}/{{\gamma _0}}}$, 
${{{{Re}}_i} = {{{\rho _i}U{R_1}}/{{\mu _i}}}}$.
3. Analytical solutions
3.1. Concentric configuration
In this analysis, the drop and the particle are assumed to be concentric. It is convenient to use spherical coordinates with the origin at the common centre as shown in figure 2. By symmetry, both the Janus particle and the surrounding drop experience a translational motion along the axis of revolution denoted as 
${{\hat {e}}_{z}}$. To proceed with the solution, it is convenient to make the governing equations dimensionless using the scaling factors 
${{R}_{1}},( J{{R}_{1}}/{{D}_{1}} )$ and 
$U$ for length, concentration and velocity.
Figure 2. Spherical basis for concentric configuration depicting the active cap facing the 
${{\hat {e}}_{z}}$ direction.
3.1.1. Concentration field
The general axisymmetric solution to (2.1) in spherical coordinates is given as (Subramanian et al. Reference Subramanian, Balasubramaniam and Clark2002)
\begin{equation} C = \sum_{n = 0}^\infty {( {{X_n}{r^n} + {Y_n}{r^{ - ( {n + 1} )}}} )} {p_n}( \vartheta ), \end{equation}
where 
$\vartheta =\cos \theta$ and 
${{p}_{n}}( \vartheta )$ are the 
$n$th-order Legendre polynomials. Since the concentration decays to 0 far field, i.e. 
${{C}_{2}}\to 0$ as 
$r\to \infty$, the outer and inner concentration fields are given as
\begin{gather} {C_1} = \sum_{n = 0}^\infty {( {{X_n}{r^n} + {Y_n}{r^{ - ( {n + 1} )}}} )} {p_n}( \vartheta ), \end{gather}
\begin{gather}{C_2} = \sum_{n = 0}^\infty {( {{y_n}{r^{ - ( {n + 1} )}}} )} {p_n}( \vartheta ). \end{gather}The boundary conditions (2.2)–(2.5) can be written in spherical coordinates as follows:
\begin{gather} \frac{{\partial {C_1}}}{{\partial r}} = D\frac{{\partial {C_2}}}{{\partial r}}\{ r = {\chi ^{ - 1}}\} , \end{gather}
\begin{gather}{C_1} = k {C_2}\{ r = {\chi ^{ - 1}}\} , \end{gather}
\begin{gather}\frac{{\partial {C_1}}}{{\partial r}} ={-} 1 \{ {r = 1,0 \le \theta \le \alpha } \}, \end{gather}
\begin{gather}\frac{{\partial {C_1}}}{{\partial r}} = 0\{ {r = 1,\alpha < \theta \le {\rm \pi}} \}, \end{gather}
where 
$\chi ={{R}_{1}}/{{R}_{2}}$ and 
$\alpha$ is the cap angle. Using orthogonality relations of Legendre polynomials
\begin{equation} \int_{{-}1}^1 {{p_n}} ( \vartheta ){p_m}( \vartheta )\, {\rm d} \vartheta = \left\{ \begin{array}{@{}cc} {\dfrac{2}{{2n + 1}}} & {\{ {n = m} \}}\\ 0 & {\{ {n \ne m} \}} \end{array} \right., \end{equation}
the boundary conditions can be converted into a system of equations for the coefficients 
$(X,Y)$. As the equations for individual orders of the Legendre polynomial are decoupled, they allow for an explicit solution for each coefficient and this is given in the Supplementary material available at https://doi.org/10.1017/jfm.2023.5.
3.1.2. Velocity field
It is convenient to solve the problem in the moving reference frame attached to the drop. As the problem is completely axisymmetric, it is possible to simplify the vector equations into a single equation for the Stokes stream function 
$\psi$. The corresponding equation for the stream function for Stokes flow in spherical coordinate system is given by (Leal Reference Leal2007)
\begin{equation} {E^4}( \psi ) = 0,\end{equation}
where 
${{E}^{2}}= {{{\partial }^{2}}}/{\partial {{r}^{2}}}+({( 1-{{\vartheta }^{2}} )}/{{{r}^{2}}})({{{\partial }^{2}}}/{\partial {{\vartheta }^{2}}})$. The velocity fields can be obtained from the stream function using the relation 
$\boldsymbol {u}=\boldsymbol {\nabla } \times ( ({\psi }/{\omega }){{{\hat {\boldsymbol {e}}}_{\phi }}} )$ where 
$\omega$ is the perpendicular distance from the axis of rotation,
\begin{gather} \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{{\hat e_r}} ={-} \frac{1}{{{r^2}}}\frac{{\partial \psi }}{{\partial \vartheta }}, \end{gather}
\begin{gather}\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{{\hat e_\theta }} ={-} \frac{1}{{r\sqrt {1 - {\vartheta ^2}} }}\frac{{\partial \psi }}{{\partial r}}. \end{gather}
The stream function is scaled by 
$U{{R}_{1}}^{2}$. The general solutions to the stream function are given as
\begin{equation} \psi = \sum_{n = 1}^\infty {( {{A_n}{r^{n + 3}} + {B_n}{r^{n + 1}} + {C_n}{r^{2 - n}} + {D_n}{r^{ - n}}} )} {g_n}( \vartheta ),\end{equation}
where 
${{g}_{n}}( \vartheta )$ are the modified Gegenbauer polynomials of order n (Leal Reference Leal2007). These are related to the traditional Gegenbauer polynomials 
$\mathbb {C}_{n}^{m}$ as 
${{g}_{n}}( \vartheta )=-\mathbb {C}_{n+1}^{-1/2}( \vartheta )$. The derivative of these polynomials can also be related to the Legendre polynomials of corresponding order,
\begin{equation} \frac{{{\rm d}{g_n}( \vartheta )}}{{{\rm d}\vartheta }} = {p_n}( \vartheta ).\end{equation}The modified Gegenbauer polynomials obey the following orthogonality relations:
\begin{equation} \int_{ - 1}^1 {\frac{{{g_n}( \vartheta ){g_m}( \vartheta )}}{{( {1 - {\vartheta ^2}} )}}\, {\rm d} \vartheta } = \left\{ \begin{array}{@{}cc} \dfrac{2}{{( n )( {n + 1} )( {2n + 1} )}} & {\{ {n = m} \}}\\ 0 & {\{ {n \ne m} \}} \end{array} \right.. \end{equation}
As in the previous analysis, the stream functions cannot contain singular terms within the domain of their applicability. Also, since the flow far from the drop is a uniform streaming flow 
$-U{{\boldsymbol {\hat {e}}_{z}}}$ in the reference frame attached to the drop, the stream function tends to
\begin{equation} {\boldsymbol{u_2}} \to - 1{\hat{\boldsymbol e_z}}\ {\rm or} \ {\psi _2} \to \frac{{{r^2}}}{2}( {{\vartheta ^2} - 1} )\{ {r \to \infty } \}. \end{equation}Correspondingly, we seek general solutions to the stream function in each fluid of the form
\begin{gather} {\psi _1} = \sum_{n = 1}^\infty {( {{A_n}{r^{n + 3}} + {B_n}{r^{n + 1}} + {C_n}{r^{2 - n}} + {D_n}{r^{ - n}}} )} {g_n}( \vartheta ), \end{gather}
\begin{gather}{\psi _2} = {r^2}{g_1}( \vartheta ) + \sum_{n = 1}^\infty {( {{c_n}{r^{2 - n}} + {d_n}{r^{ - n}}} )} {g_n}( \vartheta ). \end{gather}Substituting (3.10), (3.11) in the boundary conditions (2.8)–(2.11), we get
\begin{gather} {\psi _1} = {\psi _2} = 0\{ {r = {\chi ^{ - 1}}} \}, \end{gather}
\begin{gather}\frac{{\partial {\psi _1}}}{{\partial r}} = \frac{{\partial {\psi _2}}}{{\partial r}}\{ {r = {\chi ^{ - 1}}} \}, \end{gather}
\begin{gather}
\frac{1}{r}\left( {\mu \frac{{{\partial ^2}{\psi
_2}}}{{\partial {r^2}}} - \frac{{{\partial ^2}{\psi
_1}}}{{\partial {r^2}}}} \right) - \frac{2}{{{r^2}}}( {\mu
- 1} )\frac{{\partial {\psi _1}}}{{\partial r}}\nonumber\\
\quad ={-} \left( {\frac{{\beta {\gamma _0}{k_2}}}{{U{\mu
_1}}}\left( {\frac{{J{R_1}}}{{{D_1}}}} \right)}
\right)\sum_{n = 0}^\infty {\frac{{n( {n + 1}
){b_n}}}{{{{\left( {\dfrac{1}{\chi }} \right)}^{n + 2}}}}}
{g_n}( \vartheta )\{ {r = {\chi ^{ - 1}}} \},
\end{gather}
\begin{gather} - \frac{1}{{{r^2}}}\frac{{\partial {\psi _1}}}{{\partial \vartheta }} = \left( {\frac{V}{U}} - 1 \right)\vartheta \{ {r = 1} \}, \end{gather}
\begin{gather}- \frac{1}{r}\frac{{\partial {\psi _1}}}{{\partial r}} = \left( {\frac{V}{U}} - 1 \right)( {{\vartheta ^2} - 1} )\{ {r = 1} \}. \end{gather}
Upon substituting (3.16) and (3.17) into the above equations, a set of linear equations is obtained. These equations can be simultaneously solved to evaluate the coefficients and thereby the stream function for either fluid (see Supplementary material). The solution for the two velocities can then be evaluated by using the force-free conditions (2.13), (2.14) which reduces to setting the coefficients of 
${{g}_{1}}$ linear in the radius for the two stream functions to 0, i.e. 
$( {{C}_{1}},{{c}_{1}}=0 )$ (see Leal Reference Leal2007, p. 462). The explicit solutions for the terminal velocity are thus obtained to be
\begin{gather} \frac{U}{{ {\dfrac{{\beta {\gamma _0}{k_2}{J{R_1}}}}{{{\mu _1}{D_1}}} } }} = \frac{{ - 3( { - 1 + {\chi ^{ - 5}}} ) {{{\sin }^2}\alpha }}}{{4\left( {2 - 2\mu + \dfrac{{3 + 2\mu }}{{{\chi ^5}}}} \right)\left( {k - D - \dfrac{{k + 2D}}{{{\chi ^3}}}} \right)\chi }}, \end{gather}
\begin{gather}V = \left( {1 + \frac{{( { - 1 + {{\chi ^{ - 1}}}} )( {2 + 4 {{\chi ^{ - 1}}} + 6 {{\chi ^{ - 2}}} + 3 {{\chi ^{ - 3}}} } )}}{{2( {1 + {{\chi ^{ - 1}}} + {{\chi ^{ - 2}}} + {{\chi ^{ - 3}}} + {{\chi ^{ - 4}}} } )}}} \right)U. \end{gather}
Evidently, the natural scale for the Marangoni-driven motor velocity is 
${{ {{{\beta {\gamma _0}{k_2}{J{R_1}}}}/{{{\mu _1}{D_1}}} } }}$.
3.2. Eccentric configuration
Since both the drop and the Janus particle move along the common axis of revolution, it is evident that the Janus particle will translate along the line joining the centres of the two spheres at least initially. Further, since such eccentric configurations are also axisymmetric, there will be no torque on the Janus particle or the drop due to symmetry arguments. Thus, we shall consider all eccentric configurations where the Janus particle is displaced along the line joining their centres and the entire compound drop/Janus particle system can be described as a body of revolution about the axis passing through both their centres. However, it must be noted that the Janus particle and the drop need not remain axisymmetric under the action of perturbations. A stability analysis must be conducted to determine the operational regime wherein these eccentric configurations are stable to outside disturbance. Such an analysis is left for a future study.
In order to solve this system, it is convenient to formulate the problem in bispherical coordinates 
$( \xi,\eta,\phi )$ as the constant coordinate surfaces of 
$\xi$ can be used to construct the eccentric spheres making it convenient to satisfy the boundary conditions on either surface (figure 3). The relevant conversions between the bispherical and the more familiar cylindrical coordinates 
$( \xi,\eta,\phi )\leftrightarrow ( z,\omega,\phi )$ and corresponding unit vectors are given as (Subramanian et al. Reference Subramanian, Balasubramaniam and Clark2002)
\begin{equation} z = \frac{{p\sinh \xi }}{{\cosh \xi - \cos \eta }},\quad \omega = \frac{{p\sin \eta }}{{\cosh \xi - \cos \eta }}. \end{equation}From Happel & Brenner (Reference Happel and Brenner2012), the metric coefficients h (in our notation) corresponding to the coordinates q are given as
\begin{equation} {h_1} = {h_2} = \frac{{( {\cosh \xi - \cos \eta } )}}{p}, \quad {h_3} = \frac{{( {\cosh \xi - \cos \eta } )}}{{p\sin \eta }}, \end{equation}
where 
$( {{q_1},{q_2},{q_3}} ) = ( {\eta,\xi,\phi } )$ and 
$p=|{{R}_{1}}\sinh {{\xi }_{1}}|$. The unit vector along the axis of assumed motion then becomes
\begin{equation} {\hat{\boldsymbol e_z}} = {\hat{\boldsymbol e_\xi} }( {{h_2}} )\frac{{\delta z}}{{\delta \xi }} + {\hat{\boldsymbol e_\eta} }( {{h_1}} )\frac{{\delta z}}{{\delta \eta }}. \end{equation}
It is convenient to use 
$p,\,Jp/{{D}_{1}},U,U{{p}^{2}}$ as the scaling factors for the length, concentration, velocity and stream function and 
$\zeta =\cos \eta$. Thus,
\begin{equation} {\hat{\boldsymbol e_z}} = {\hat{\boldsymbol e_\xi} }\frac{{( {1 - \zeta \cosh \xi } )}}{{( {\cosh \xi - \zeta } )}} - {\hat{\boldsymbol e_\eta} }\frac{{\sinh \xi \sqrt {( {1 - {\zeta ^2}} )} }}{{( {\cosh \xi - \zeta } )}}. \end{equation}
The range of relative positions of the two spheres can be classified into two cases (see figure 3). Case 1 corresponds to the configurations where 
${{\xi }_{1}}>{{\xi }_{2}}, \{ {{\xi }_{i}}>0 \}$ and case 2 refers to the remaining configurations 
${{\xi }_{1}}<{{\xi }_{2}}, \{ {{\xi }_{i}}<0 \}$. Further discussion of the bispherical coordinate system can be undertaken simultaneously for both cases by choosing the upper sign for case 1 and the lower sign for case 2 wherever applicable. The coordinate system is then completely defined by choosing 
$\chi ={{R}_{1}}/{{R}_{2}}$, 
$\upsilon =d/({{R}_{2}}-{{R}_{1}})$ (where d is the centre-to-centre distance between the Janus particle and drop) and 
$\lambda =\mp \upsilon$ which represents the non-dimensional position of the Janus particle relative to the drop centre along the axis joining the two spheres (Subramanian et al. Reference Subramanian, Balasubramaniam and Clark2002):
Figure 3. Bispherical basis for eccentric configuration(s) when the position vector of the particle relative to the drop centre is along (a) 
${{\hat {\boldsymbol {e}}_{z}}}$ (case 1) or along (b) 
$-{{\hat {\boldsymbol {e}}_{z}}}$ (case 2).
\begin{gather} {\xi _1} ={\pm} \,{\rm arcosh}\left[ {\frac{{{\upsilon ^2}( {\chi - 1} ) + ( {1 + \chi } )}}{{2\upsilon \chi }}} \right], \end{gather}
\begin{gather}{\xi _2} ={\pm} \,{\rm arcosh}\left[ {\frac{{{\upsilon ^2}( {\chi - 1} ) - ( {1 + \chi } )}}{{2\upsilon }}} \right]. \end{gather}
As a convention, we shall take positive and negative values of 
$\lambda$ to mean displacement along and against the direction faced by the Janus cap.
3.2.1. Concentration field
The Laplace equation for an axisymmetric scalar field can be written in bispherical coordinates as (see Leal Reference Leal2007, p. 448)
\begin{equation} {\nabla ^2} = {( {\cosh \xi - \zeta } )^3}\left( {\frac{\partial }{{\partial \xi }}\left( {\frac{1}{{\left( {\cosh \xi - \zeta } \right)}}\frac{\partial }{{\partial \xi }}} \right) + \frac{\partial }{{\partial \zeta }}\left( {\left( {\frac{{1 - {\zeta ^2}}}{{\cosh \xi - \zeta }}} \right)\frac{\partial }{{\partial \zeta }}} \right)} \right), \end{equation}and the general axisymmetric solution to (2.1) is (Subramanian et al. Reference Subramanian, Balasubramaniam and Clark2002)
\begin{equation} C = \sqrt {\cosh \xi - \zeta } \sum_{n = 0}^\infty {[ {{X_n}\,{{\rm e}^{(n + 1/2)\xi }} + {Y_n}\,{{\rm e}^{ - (n + 1/2)\xi }}} ]} {p_n}( \zeta). \end{equation}
As shown in figure 3, the coordinate 
$\xi$ tends to infinity at a finite distance from the drop surface. To avoid the resulting singularities in the field variables, the general solution for the concentration field in both fluids can be written in the following form convenient for applying the appropriate boundary conditions:
\begin{gather} {C_1} = \sqrt {\cosh \xi - \zeta } \sum_{n = 0}^\infty {[ {{X_n}\,{{\rm e}^{( {n + 1/2} )( {\xi - {\xi _2}} )}} + {Y_n}\,{{\rm e}^{ - ( {n + 1/2} )( {\xi - {\xi _2}} )}}} ]} {p_n}( \zeta ), \end{gather}
\begin{gather}{C_2} = \sqrt {\cosh \xi - \zeta } \sum_{n = 0}^\infty {[ {{x_n}\,{{\rm e}^{ {\pm} ( {n + 1/2} )( {\xi - {\xi _2}} )}}} ]} {p_n}( \zeta ), \end{gather}
since 
$\xi \to \mp \infty$ is within the domain of the outer fluid. Applying (2.4), (2.5) at 
$\xi ={{\xi }_{2}}$ gives
\begin{gather} {C_1} = k {C_2}, \end{gather}
\begin{gather}\frac{{\partial {C_1}}}{{\partial \xi }} = D \frac{{\partial {C_2}}}{{\partial \xi }}. \end{gather}
Applying (2.2), (2.3) at 
$\xi ={{\xi }_{1}}$ we get
\begin{equation} \frac{{\partial {C_1}}}{{\partial \xi }} ={\pm} \frac{F}{{( {\cosh {\xi _1} - \zeta } )}} , \quad \left\{ F = \begin{array}{c@{}} {1, \eta \in \{ {0,{\eta ^*}} \}}\\ {0, \eta \in \{ {{\eta ^*},{\rm \pi} } \}} \end{array} \right\}. \end{equation}
Here F is the non-dimensional flux along 
$\hat {n}=\mp \hat {\xi }$ and the range of values of 
$\eta$ between 
$(0,\eta ^*)$ depicts the active cap. Unlike the simpler case with the concentric configuration, these equations do not easily yield a simple explicit solution, but upon substituting the general solution for both the fluids into these equations, a coupled system of linear equations for the coefficients of the general solution is obtained (see Supplementary material). The corresponding matrix for the infinite system of linear equations can be truncated as required to obtain a sufficiently converged solution to the coefficients of the concentration field for the inner and outer fluid.
3.2.2. Velocity field
Due to the axisymmetric nature of the problem, we adopt a similar approach as used in the earlier section for the concentric configuration and seek solutions for the velocity fields using the Stokes stream functions. However, we shall solve for the velocity fields in the stationary reference frame as it proves to be more convenient while applying the far field boundary conditions,
\begin{equation} ( {{u_\xi },{u_\eta }} ) = \left( {\frac{{\cosh \xi - \zeta }}{{\sin \eta }}} \right)\left( {\frac{{\partial \psi }}{{\partial \eta }}, - \frac{{\partial \psi }}{{\partial \xi }}} \right), \end{equation}
where the stream function satisfies (3.9) (Subramanian et al. Reference Subramanian, Balasubramaniam and Clark2002) with 
${{E}^{2}}=( \cosh \xi -\zeta )( ({\partial }/{\partial \xi })( ( \cosh \xi -\zeta )({\partial }/{\partial \xi }) )+( 1-{{\zeta }^{2}} )({\partial }/{\partial \zeta })( ( \cosh \xi -\zeta )({\partial }/{\partial \zeta }) ) )$. Based on similar considerations as mentioned for the concentration field, the general solutions for the stream functions of the two fluids are written as
\begin{gather} {\psi _1} = {( {\cosh \xi - \zeta } )^{ - 3/2}}\sum_{n = 1}^\infty {\left[ \begin{aligned} & {A_n}\,{{\rm e}^{( {n - 1/2} )( {\xi - {\xi _2}} )}} + {B_n}\,{{\rm e}^{ - ( {n - 1/2} )( {\xi - {\xi _2}} )}}\\ & \quad + {C_n}\,{{\rm e}^{( {n + 3/2} )( {\xi - {\xi _2}} )}} + {D_n}\,{{\rm e}^{ - ( {n + 3/2} )( {\xi - {\xi _2}} )}} \end{aligned} \right]} {g_n}( \zeta ), \end{gather}
\begin{gather}{\psi _2} = {( {\cosh \xi - \zeta } )^{ - 3/2}}\sum_{n = 1}^\infty {[ {{a_n}\,{{\rm e}^{ {\pm} ( {n - 1/2} )( {\xi - {\xi _2}} )}} + {b_n}\,{{\rm e}^{ {\pm} ( {n + 3/2} )( {\xi - {\xi _2}} )}}} ]} {g_n}( \zeta ). \end{gather}
Applying (2.9)–(2.11) at 
$\xi ={{\xi }_{2}}$ gives
\begin{gather} {\boldsymbol{u_1}} \boldsymbol{\cdot} {\hat{\boldsymbol e_\xi} } = {\boldsymbol{u_2}} \boldsymbol{\cdot} {\hat{\boldsymbol e_\xi} } = 1 {{\boldsymbol{{\hat e}_\xi} } \boldsymbol{\cdot} {\boldsymbol{{\hat e}_z}}} , \end{gather}
\begin{gather}{\boldsymbol{u_1}} \boldsymbol{\cdot} {\hat{\boldsymbol e_\eta} } = {\boldsymbol{u_2}} \boldsymbol{\cdot} {\hat{\boldsymbol e_\eta} }, \end{gather}
\begin{gather}{\hat{\boldsymbol e_\xi} } \boldsymbol{\cdot} ( {{\boldsymbol{\sigma _1}} - {\boldsymbol{\sigma _2}}} ) \boldsymbol{\cdot} {\hat{\boldsymbol e_\eta} } ={\mp} \frac{1}{{{\mu _1}U}}( {{\boldsymbol{\nabla} _s}\gamma } ) \boldsymbol{\cdot} {\hat{\boldsymbol e_\eta} }, \end{gather}where
\begin{equation} {\hat{\boldsymbol e_\xi} } \boldsymbol{\cdot} \boldsymbol{\sigma_i} \boldsymbol{\cdot} {\hat{\boldsymbol e_\eta} } = {\sigma_{i\xi \eta }} =\left(\frac{\mu _i}{\mu _1}\right) \frac{{{{( {\cosh \xi - \zeta } )}^{3/2}}}}{{\sqrt {1 - {\zeta ^2}} }} \left[ \begin{aligned} & ( {1 - {\zeta ^2}} )\left[ {( {\cosh \xi - \zeta } )\frac{{{\partial ^2}\psi_i }}{{\partial {\zeta ^2}}} - 3\frac{{\partial \psi_i }}{{\partial \zeta }}} \right]\\ & \quad - \left[ {( {\cosh \xi - \zeta } )\frac{{{\partial ^2}\psi_i }}{{\partial {\xi ^2}}} + 3\sinh \xi \frac{{\partial \psi_i }}{{\partial \xi }}} \right] \end{aligned}\right].\end{equation}
Using the partition coefficient, the surface concentration can be related to either bulk concentration (here for the outer fluid) 
$\varGamma ={{k}_{2}}{{C}_{2}}$. Thus
\begin{equation} {\boldsymbol{\nabla}_s}\gamma ={-} \left( {\frac{{{\beta {\gamma _0}{k_2}}Jp}}{{{D_1}}}} \right)( {\cosh \xi - \zeta } )\frac{{\partial {C_2}}}{{\partial \eta }}.\end{equation}
Similarly, applying boundary condition (2.8) at 
$\xi ={{\xi }_{1}}$ gives
\begin{gather} {\boldsymbol{u_1}} \boldsymbol{\cdot} {\hat{\boldsymbol e_\xi} } = \left( \frac{V}{U} \right){\boldsymbol{{{\hat e}_\xi }} \boldsymbol{\cdot} \boldsymbol{{{\hat e}_z}}} , \end{gather}
\begin{gather}{\boldsymbol{u_1}} \boldsymbol{\cdot} {\hat{\boldsymbol e_\eta} } = \left( \frac{V}{U} \right) {\boldsymbol{{{\hat e}_\eta } }\boldsymbol{\cdot} \boldsymbol{{{\hat e}_z}}}. \end{gather}As in the case of the concentration field, the boundary conditions can be used to set up a system of linear equations for the coefficients of the general solution (see Supplementary material). We can then truncate the resulting matrix to obtain an approximate solution for the unknown coefficients. Once the velocity fields are known, the terminal velocities U and V can be found using the force-free relations (2.13), (2.14) on the drop and the Janus particle, respectively.
4. Results and discussion
4.1. Concentric configuration
The terminal velocities of the Janus particle and the drop as a function of the dimensionless thickness of the annulus of liquid between the Janus particle and the drop interface 
$( \tau =({{{R}_{2}}-{{R}_{1}}})/{{{R}_{1}}} )$ when they are concentrically placed is plotted in figure 4. These and further results are computed using the first 20 terms of the series expansions for both the stream function and the concentration field and the velocities of the particle and the drop are reported as dimensionless quantities scaled by the Marangoni scaling factor 
${\beta {{\gamma }_{0}}{{k}_{2}}J{{R}_{1}} }/{{{\mu }_{1}}{{D}_{1}}}$ in all subsequent plots. The terminal velocities and the solutions for the relevant concentration fields and stream functions are depicted using the solutions from § 2.1.
Figure 4. (a) Plot of the (dimensionless) velocities of the a(i) drop and the a(ii) Janus particle as a function of the thickness for different cap angles 
$( \alpha ={{\rm \pi} }/{3},{{\rm \pi} }/{2},{2{\rm \pi} }/{3} )$ keeping 
$( k=D=\mu =1 )$; (b) heat map of the solute concentration for cap angles 
$( \alpha ={{\rm \pi} }/{2}\pm {{\rm \pi} }/{6} )$ keeping 
$( k=D=\mu =1 )$ and 
$( \tau =1 )$.
It is evident from the analytical solutions that the magnitude of the velocity of the Janus particle relative to that of the drop is always greater than 0. Consequently, the drop and the Janus particle cannot remain concentric. However, since the Stokes equations are quasisteady, the solutions for the concentric case are correct in the instant that the Janus particle migrates to the drop centre. Another interesting feature of these solutions is that the terminal velocities of both the drop and the particle are functions of 
$\sin {\alpha }$, where 
$\alpha$ is the cap angle. Thus, the particles with cap angles 
$( \alpha ={{\rm \pi} }/{2}\pm {{\alpha }_{0}} )$ behave identically. Figure 4(a) depicts this phenomenon for 
${{\alpha }_{0}}={{\rm \pi} }/{6}$. This result can be understood by comparing the gradients of the concentration field generated by these cap sizes. Figure 4(b) is a heat map of the concentration field. The propulsion forces generated from the concentration gradients can be similar for large and small cap sizes even though the larger cap generates more solute and thereby can have a larger surfactant concentration near the interface.
Finally, figure 4(a) also indicates that there exists an optimum drop size at which the particle and drop velocities attain a maximum value. This can be understood as a consequence of the competing effects of increasing the drop size. On one hand, a larger drop has greater interfacial area and therefore can generate larger propulsion forces and on the other, it can increase the hydrodynamic drag from the surrounding fluid.
4.2. Eccentric configuration
The motion of the Janus particle and the drop were analysed using the formulations in § 2 for various eccentric positions. In figure 5, the dimensionless velocities of the drop and the Janus particle are plotted as a function of 
$\lambda$ for three different Janus angles. The points along the 
${x} = 0$ axis correspond to the situation where the drop and the particle are concentric and are the results derived from the analysis using the concentric configuration. The remaining data points are computed using the eccentric analysis from § 2.2. We can numerically verify that the bispherical solutions agree with the concentric solutions as 
$\lambda \to 0$.
Figure 5. Dimensionless velocities of the (a–c) drop and the (d–f) Janus particle for different cap angles 
$( \alpha ={{\rm \pi} }/{3},{{\rm \pi} }/{2},{2{\rm \pi} }/{3})$ keeping 
$( k=D=\mu =1 )$.
We can define two directions ‘forward’ and ‘backward’ to represent the direction along the Janus particle/drop axis penetrating through the Janus cap and the reverse direction (i.e. 
$+{{\hat {e}}_{z}}$ and 
$-{{\hat {e}}_{z}}$, respectively, in our formulations). When the Janus particle is located in the region 
$0\le \lambda <1$, the Janus cap generates a region of higher concentration of surface active solute in its immediate vicinity which in turn adsorbs onto the drop surface creating a region of low surface tension. Since the location on the drop surface closest to the Janus cap is the point of intersection of the drop surface and the axis of the Janus particle, this region will have the lowest tension. From the symmetry in the problem, it is then obvious that the surface tension monotonically increases until it reaches its highest value at the point on the drop surface along the same axis in the backward direction. The gradients thus created drive the Janus particle/drop in the forward direction. Therefore, it is conceivable that the closer the particle is to the drop surface, the stronger these gradients get and the stronger the resulting propulsive forces. On the other hand, when the Janus particle is in the region 
$-1<\lambda <0$, the concentration field produced by the Janus cap exhibits a maximum value along the drop interface at a location that is not necessarily along the axis of the drop/Janus particle. Thus, the symmetry of the problem requires the setting up of two opposing directions of concentration gradients (see figure 6a). The resulting competing forces can propel the Janus particle and the drop either forward or backward depending on various other factors.
Figure 6. (a) Illustration of the opposing directions of the excess surface concentration gradient generated due to the location maximum in the surface excess being off the line of centres of the drop and the Janus particle 
$( \theta =0 )$ for positions of the Janus particle 
$\lambda \in ( -1,0 )$; (b) illustration of a stationary position where the drop and the Janus particle (dimensionless) velocities are equal for 
$( \tau =1,\alpha ={2{\rm \pi} }/{3})$ keeping 
$( k=D=\mu =1 )$; (c) surface excess concentration distribution for the configuration in panel (b) for the interpolated stationary position of 
$\lambda = -0.48$ (from inset) with 
${\varGamma _0} = {k_2}J{R_1}/{D_1}$ (
$\theta$ is defined anticlockwise about the drop centre from the 
$z$ axis as shown in panel (a)).
The presence of positions where the drop and the Janus particle may move together is of interest for practical applications. While the concentric concentration solutions explicitly disallow equal velocities, the system of equations for the bispherical formulation (in this case) have been numerically inverted and further analysis is necessary to make any quantitative statements regarding the existence and prediction of such positions. However, it is worth noting that any such position that lies in the second or third quadrants of the plots in figure 5 is inherently unstable in the sense that any slight displacement to the position of the particle can lead to the system moving away from that configuration. This also implies that the Janus particle and the drop do not naturally attain this position unless artificially positioned that way. For example, figure 6(b) depicts one such stationary position for a particular choice of the problem parameters. Figure 6(c) depicts the surface concentration of the solute for the stationary position depicted in figure 6(b). It is evident that the surface concentration has a maximum at an arbitrary location along the drop surface and will therefore exhibit surface concentration gradients along two directions. From (2.11) and (2.12), it is then evident that the surface concentration gradient will affect the direction of the induced surface stresses along the drop surface leading to the stationary position.
If the terminal velocity of the Janus particle always differs from that of the drop, the Janus particle will approach the limiting eccentric positions of 
$\lambda \to \pm 1$ eventually. While the bispherical solutions are valid for the entire range of 
$\lambda$, the solutions do not converge sufficiently fast as the Janus particle comes close to the drop surface. This is due to the inherent mathematical limitation of the bispherical coordinate system which is not well defined in the limiting case of the Janus particle touching the drop surface. Furthermore, as the Janus particle approaches the drop surface, the need to force fluid out of the thin region between the drop and the particle will generate large restoring lubrication pressures that reduce the velocity of the Janus particle. To study these limits, we would have to employ numerical techniques or resort to an asymptotic matching of the velocity fields using a lubrication model for the inner flow in the thin gap and a solution based on tangent sphere basis (Moon & Spencer Reference Moon and Spencer2012) for the outer flow.
Even well within the drop, the trends displayed by the drop and Janus particle velocities at different eccentric locations are widely different and we can attempt to explain this by comparing these solutions to those from asymptotic formulations based on the plausible dominant physics for different drop thicknesses. In particular, the reduced velocities of the Janus particle and the drop for low values of 
$\tau$ can be a consequence of lubrication forces which would also explain why the terminal velocities are not nearly as monotonic as the solutions for larger 
$\tau$. Thus, a simplified model of the problem highlighting the effects of lubrication should be capable of reproducing these trends and validating our reasoning. Similarly, the limit of large 
$\tau$ produces nearly symmetric behaviour about the two directions of eccentricity. As the particle becomes sufficiently small, its size becomes increasingly negligible in determining the dynamics of the two spheres. It is then conceivable that the dominant physics can be reproduced using an asymptotic model that takes advantage of the existence of this small parameter (namely, the relative size of the Janus particle to the drop). In the subsequent sections, the effects of the limiting values of 
$\tau$ are reproduced using such simplified analytical models.
4.2.1. Lubrication limit 
$( \tau \to 0 )$
When the Janus particle radius approaches that of the drop, it is more convenient to think of the drop as a thin liquid layer on the spherical Janus particle. To avoid the issue of multiple scales, we restrict ourselves to the configurations wherein the eccentricity in the Janus particle location relative to the drop is also of the same order as the thickness of the drop annulus. The problem then reduces to a lubrication analysis for the fluid dynamics inside the drop along with a leading-order analysis of the flow in the unbounded outer fluid. The velocities for the lubrication regime are
\begin{gather} U \sim V \sim \tau {U'} + O( {{\tau ^2}} ), \end{gather}
\begin{gather}\frac{{U'}}{{\dfrac{{{k_2}{\gamma _0}\beta J{R_1}}}{{{\mu _1}{D_1}}}}} = \frac{1}{8}\left( {6B'\left( {\frac{{( {\lambda + ( {{\lambda ^2} - 1} )\,{\rm arctanh}\,\lambda } )}}{{{\lambda ^3}}}} \right) - A'}\right). \end{gather}
See Appendix A for the complete analysis. The coefficients 
${A}'$ and 
${B}'$ are defined in the appendix and are not repeated here for the sake of brevity. The Marangoni scaling for the terminal velocities
$( {{{k}_{2}}{{\gamma }_{0}}\beta J{{R}_{1}}}/{{{\mu }_{1}}{{D}_{1}}} )$ is retained in the lubrication limit. These results can be validated by plotting them against the bispherical solutions for small values of 
$\tau$ (figure 7).
Figure 7. Validation of the lubrication analysis: comparison of the (dimensionless) velocities from the lubrication solution (4.2) with the bispherical solutions for different cap angles 
$( \alpha ={{\rm \pi} }/{3}, {{\rm \pi} }/{2}, {2{\rm \pi} }/{3})$ keeping 
$( k=D=\mu =1 )$ and 
$( \tau =0.01 )$.
From this analysis, we conclude that the drop and Janus particle velocities are equal to leading order and that they scale as 
$O( \tau )$. This result is unsurprising as the analytical solutions for the terminal velocity from the concentric configuration (which are a good indicator for the scale of the velocity) can be expanded in 
$\tau$ to obtain the same information. Nevertheless, having an analytical formula for all eccentric configurations gives us further insights. In this limit, the velocities of the drop and the Janus particle are essentially the same and can by plotted for different cap angles for comparison (figure 8). Interestingly, there seems to exist an antisymmetry in the dynamics of the drop along the directions of eccentricity for Janus particles with supplementary cap angles. Further, the maximum velocity of the drop/Janus particle is achieved when the Janus particle is concentrically located with its cap spanning exactly half its surface area (i.e. 
$\alpha ={\rm \pi} /2$). The fact that the relative velocity between the drop and the Janus particle scales as at least 
$O( {{\tau }^{2}} )$ enables us to design colloidal motors in this regime which can retain their relative configuration for extended periods of time.
Figure 8. Leading-order (dimensionless) velocities of both drop and Janus particle in the lubrication limit 
$( \tau =0.01 )$ for different cap angles 
$( \alpha ={{\rm \pi} }/{6}, {{\rm \pi} }/{4}, {{\rm \pi} }/{3}, {{\rm \pi} }/{2}, {2{\rm \pi} }/{3}, {3{\rm \pi} }/{4}, {5{\rm \pi} }/{6})$ keeping 
$( k=D=\mu =1 )$.
4.2.2. Point particle limit $( \tau \to \infty )$
(a) Point source formulation.
Due to the diffusive nature of the solute transport, the relative orientation of the particle (direction of the Janus cap) becomes unimportant and when the Janus particle is sufficiently small, its position with respect to the centre of the drop becomes the dominating geometric factor in determining its dynamics. In this limit, we can consider the particle to be represented as a point source of the solute. The symmetry of the point source ensures that the system is still axisymmetric and can admit solutions in the spherical coordinate system centred at the drop. Correspondingly, we can modify the governing equation of the concentration field in the inner fluid to
\begin{equation} D_{1}{\nabla^2}{C_1} = 4{\rm \pi} {N_p}\delta ( {r{\boldsymbol{{\hat e}_r}} - {R_p}{\boldsymbol{{\hat e}_z}}} ), \end{equation}
where the point source is located at 
${{R}_{p}}{{\hat {\boldsymbol {e}}_{z}}}$ and has the strength 
${{N}_{p}}=({1}/{4{\rm \pi} })\int _{{{S}_{I}}}{J\, {\rm d} s}=J{{R}_{1}}^{2}({( 1-\cos \alpha )}/{2})$.
It is evident from the force-free nature of the problem that the monopole hydrodynamic contribution of the sphere (i.e. the point force contribution) to the velocity field is zero. Hence, in the point particle limit, the problem reduces to that of finding the velocity of the drop translating under the gradient generated by a point source of the solute/heat at the prescribed eccentric location of the Janus particle. Since the flow is governed by Stokes equations, the force-free Janus particle in the point size limit naturally assumes the velocity of the drop at the point.
The general solution to the homogeneous component (i.e. the Laplace equation) is a series expansion in Legendre polynomials much like the solutions in § 3. We can simply append a particular solution to account for the singular inhomogeneity. The dimensionless inner solution to the concentration field then takes the form
\begin{equation} {C_1} = \sum_{n = 0}^\infty { {{X_n}{r^n}} } {p_n}( \vartheta ) + \frac{1}{{4{\rm \pi} {{( {{r^2} + {R^2} - 2Rr\cos \theta } )}^{1/2}}}}. \end{equation}
Here 
$R={{{R}_{p}}}/{{{R}_{2}}}$ and the natural scales are 
${{R}_{2}}$ (length), 
$U$(velocity) and 
${4{\rm \pi} {{N}_{p}}}/{{{D}_{1}}{{R}_{2}}}$ (concentration). The outer concentration, boundary conditions at the interface, the stream function formulations and the force-free conditions for the drop surface are identical to those from § 3. The absence of the boundary conditions at the Janus particle surface is compensated by the reduced number of permitted (non-singular) eigenfunctions from the solutions to the field variables in the inner fluid. The point source concentration solution (and thus the concentration gradient along the interface) are geometrically determined by the size of the drop and the position of the point source. Therefore, as expected the terminal velocity plots exhibit mirror symmetry about the concentric configuration and are necessarily zero when the Janus particle and the drop become concentric due to the lack of propulsion (figure 9).
(b) Method of reflections formulation.
While the point source formulation is able to capture the physics of the compound drop/Janus particle system at the limit of 
$\tau \to \infty$, it is not an accurate representation of the system as the eccentricity approaches the limiting cases of 
$\lambda \to \pm 1$ as the particle size becomes relatively important. Furthermore, the asymmetry in the concentration field due to the oriented Janus cap is necessary to explain the velocity of the system when the particle approaches the origin since the point source solutions invariably lead to a symmetric concentration field which cannot propel the drop/Janus particle.
Figure 9. Leading-order (dimensionless) velocities of both (a) drop and (b) Janus particle in the point source limit 
$( \tau =99 )$ for different cap angles 
$( \alpha ={{\rm \pi} }/{6}, {{\rm \pi} }/{4}, {{\rm \pi} }/{3}, {{\rm \pi} }/{2}, {2{\rm \pi} }/{3}, {3{\rm \pi} }/{4}, {5{\rm \pi} }/{6})$ keeping 
$( k=D= \mu =1 )$.
Assuming that the particle is sufficiently small relative to the drop, the distance from the particle surface to the drop is of the scale of the drop radius in most regions except very near the drop surface. Therefore, we use the method of reflections (Happel & Brenner Reference Happel and Brenner2012) to generate an expansion for the terminal velocities of the drop and the Janus particle in increasing powers of 
$\chi$ (see Appendix B for full details including the number of reflections required to obtain the first two terms of the expansion depicted here)
\begin{gather} \frac{U}{{\left( {\dfrac{{{k_2}\beta {\gamma _0}J{R_1}}}{{{\mu _1}{D_1}}}} \right)}} = \frac{2}{{( {9 + 6\mu } )}}( {\chi {f_1}( 1 ) + {{\chi ^2}} {f_2}( 1 )} ) , \end{gather}
\begin{gather}{\tiny
\frac{V}{{\left( {\dfrac{{{k_2}\beta {\gamma
_0}J{R_1}}}{{{\mu _1}{D_1}}}} \right)}} = \frac{{{f_1}( 1
)}}{{( {9 + 6\mu } )}}\left( {\left( \begin{aligned} &
\left( {5 - 3{{\left( {\frac{d}{{{R_2}}}} \right)}^2}}
\right)\left( {\chi + \left( {\frac{{{f_2}( 1 )}}{{{f_1}( 1
)}}} \right) {{\chi ^2}} } \right)\\ & \quad + \sum_{n =
2}^\infty {\left( \begin{aligned} & \frac{{n\left( {n + 1}
\right)\left( {9 + 6\mu } \right)}}{{2\left( {1 + 2n}
\right)\left( {1 + \mu } \right)}}\left( {\frac{{{f_1}( n
)}}{{{f_1}( 1 )}}} \right)\\ & \quad \times \left(
{{{\left( {\frac{d}{{{R_2}}}} \right)}^{n - 1}} - {{\left(
{\frac{d}{{{R_2}}}} \right)}^{n + 1}}} \right)\left( {\chi
+ \left( {\frac{{{f_2}( n )}}{{{f_1}( n )}}} \right){{\chi
^2}}} \right) \end{aligned}\right)}
\end{aligned}\right)} \right)} ,
\end{gather}
\begin{gather}
{f_1}( n ) = \frac{{2n + 1}}{{2( {nk + D( {n + 1} )} )}}\nonumber\\
\quad \times \int_{ - 1}^1 {\left( {\frac{{{R_2}(
{{d^2}n + ( {1 + n} ){R_2}^2 - d( {1 + 2n} ){R_2}x} )(
{\cos \alpha - 1} )}}{{2{{( {{d^2} + {R_2}^2 - 2d{R_2}x}
)}^{3/2}}}}} \right){p_n}( x )\,{{\rm d} x}},
\end{gather}
\begin{gather}{f_2}( n ) = \frac{{2n \!+\! 1}}{{2( {nk \!+\! D( {n + 1} )} )}}\int_{ - 1}^1 {\left( {\frac{{3{R_2}^2\left( \begin{array}{l} - {d^3}n + {d^2}( {2 + 3n} ){R_2}x\\ + ( {2 + n} ){R_2}^3x\\ - d{R_2}^2( {3 + n + {x^2} + 2n{x^2}} ) \end{array} \right)\sin {\alpha ^2}}}{{8{{( {{d^2} + {R_2}^2 - 2d{R_2}x} )}^{5/2}}}}} \right){p_n}( x )\,{{\rm d} x}}. \end{gather}
From the above equations, it is evident that the leading velocities for the Janus particle and the drop are both at most 
$O( \chi )$ due to the scaling. This can also be seen from the plots of the bispherical solution which illustrate the similarity between the magnitude of the velocities for low and high values of 
$\tau$. We note that the leading-order solutions from the method of reflections are numerically identical to the point source solutions. The higher-order correction thus presents the leading-order effects of the size of the Janus particle in the dynamics of the compound drop/Janus particle system. We validate these solutions by comparing them with the bispherical solutions in figure 10.
Figure 10. Validation of the method of reflections: comparison of the (dimensionless) velocities from the method of reflection solutions (4.5) and (4.6) with the bispherical solutions for different cap angles (a) 
$\alpha = {{\rm \pi} }/{3}$, (b) 
$\alpha = {{\rm \pi} }/{2}$, (c) 
$\alpha = {2{\rm \pi} }/{3}$ keeping 
$( k=D=\mu =1 )$ and 
$( \tau =9 )$.
These solutions also allow us to illustrate the importance of the size of the particle in generating the competing asymmetry to the particle location. As shown in figure 10, the method of reflections solutions presented above can predict the position where the drop or the Janus particle becomes stationary due to the previously mentioned competition between the contributions to the concentration gradient generated by the direction of the Janus cap and the Janus particle position relative to the centre of the drop. The solutions for the case of 
$\lambda =0$ are derived to from the Taylor series expansion of the solutions to the concentric configuration to 
$O( {{\chi }^{2}} )$.
4.2.3. Effect of transport properties
The effect of the various transport properties including the fluid viscosities, the diffusivity of the solute in the two fluids and the partition coefficients of the concentration field at the drop interface on the terminal velocities are illustrated in figure 11. These results are plotted using the bispherical and spherical coordinate formulations. However, the asymptotic formulations from the earlier discussion can also reproduce these results when 
$\tau \to (0,\infty )$. The effect of varying 
$\tau$ has already been established in the previous sections and is seen to be valid in all sampled combinations of 
$( k,D,\mu )$. While changing the relative magnitude of the transport properties affects the magnitude of the terminal velocities, the effect is seen to be consistent over the sampled values of 
$\tau$. For example, increasing 
$\mu$ reduces the terminal velocities irrespective of the value of 
$\tau$. Thus, we can discuss the effect of the transport properties independently. Although, the quantitative effect of the different properties may not be fully discernible from the numerical results plotted above, some qualitative features can be understood by examining the effect of these properties on the explicit analytical solutions (3.23), (3.24) for the concentric configuration. The equations show that the drop velocity can decrease monotonically when any one of the ratios 
$( k,D,\mu )$ is increased while maintaining fixed values for the others and 
$\tau$. Correspondingly, this causes a monotonic variation in the velocity of the Janus particle. Increasing 
$\mu$ increases the drag on the drop thereby decreasing the velocity of the drop. On the other hand, increasing 
$( k,D )$ decreases the propulsive Marangoni force thereby also decreasing the velocity of the compound drop/Janus particle system. It can be noted that increments in 
$D$ reduce the amount of the solute at the interface at steady state as the solute species diffuses more easily into the unbounded fluid. This reduces the magnitude of the forces generated by the Marangoni tractions. However, upon increasing 
$k$, it is likely that the surface active solute accumulates within the drop thereby reducing the tangential gradients along the fluid interface and consequently decreases the propulsion forces on the drop.
Figure 11. Dimensionless velocities of the drop and the Janus particle for low and high values of the ratio (a–c) 
$D$, (d–f) 
$k$ and (g–i) 
$\mu$ for 
$(\tau =0.01,1,99)$ keeping the other corresponding ratios as 1 and 
$\alpha ={{\rm \pi} }/{2}$.
5. Conclusion
The compound drop/Janus particle design was studied with the explicit intention of taking advantage of the dominant capillary forces at length scales where inertia and advective forces are negligible. The terminal velocities achieved by the drop and the Janus particle can be predicted by analysis of the coupled concentration (or) heat and momentum transport equations using a bispherical coordinate system or using various asymptotic models (where explicit analytical models are derivable) as explained in the previous section. We showed that either limit of the relative thickness of the annulus of liquid within the drop with respect to the size of the particle 
$\tau \to ( 0,\infty )$ produces velocities of smaller magnitudes. In the limit when 
$\tau \to 0$ , the drop size is approximately the Janus particle size and the resulting thin film dynamics within the drop create large lubrication resistances. On the contrary, the latter limit of a small Janus particle (
$\tau \to \infty$), the compound motor faces large hydrodynamic resistances from the relatively large drop size. Furthermore, the effects of the various transport properties and the idea of competing asymmetries (between the effects of the orientation and the position of the Janus particle) resulting in flipping the propulsion direction were also discussed.
If the Janus particle were undergoing neutral diffusiophoretic propulsion in fluid 1 in the absence of an interface and a second fluid, the velocity scale would typically be given by 
${{U}_{d}}\sim {{b}_{d}}J/{{D}_{1}}$ where 
${{b}_{d}}$ is the slip coefficient that relates the surface slip velocity 
${{u}_{s}}$ to the concentration gradient of the solute at the particle surface 
$( {\boldsymbol {{u}_{s}}}={{b}_{d}}{{\boldsymbol {\nabla } }_{s}}C )$. We can compare the potency of the Marangoni mechanism prescribed in this paper with neutral diffusiophoresis by analysing the magnitude of the ratio of the two velocities 
${U}/{{{U}_{d}}}\sim {{{k}_{2}}\beta {{\gamma }_{0}}{{R}_{1}}}/{{{\mu }_{1}}{{b}_{d}}}$. From Sharifi-Mood et al. (Reference Sharifi-Mood, Koplik and Maldarelli2013), the slip coefficient at a given temperature 
$T$ is given by
\begin{equation} {b_d} = \frac{{{k_b}T}}{{{\mu _1}}}{\ell _s}^2\int_0^\infty {y( {{{\rm e}^{{\varPhi _1}( y )/{k_b}T}} - 1} )}\, {{\rm d} y}, \end{equation}
where 
${{k}_{b}}$ is the Boltzmann constant and 
${{\ell }_{S}}$ is the characteristic length scale of the inner region where the effect of the interaction potential 
${{\varPhi }_{1}}$ is significant. Using the Gibb's formula for the surface excess of the solute (Peng et al. Reference Peng, Duignan, Zhao and Nguyen2020), the partition coefficient as defined in (2.4) can be written as
\begin{equation} {k_2} = {\ell _I}\int_0^\infty {( {{{\rm e}^{{\varPhi _2}( y )/{k_b}T}} - 1} )}\, {{\rm d} y}, \end{equation}
where 
${{\ell }_{I}}$ is the characteristic length scale of the interfacial adsorption zone of the solute in the fluid and 
${{\varPhi }_{2}}$ is the interaction potential between the solute and the interface. Finally, we can assume a Langmuir-type adsorption isotherm for the surface pressure 
$\varPi$ to get
\begin{equation} \frac{{\partial \varPi }}{{\partial \varGamma }} = \beta {\gamma _0} = {{{k_b}T}/{( {1 - {\varGamma/{{\varGamma _\infty }}}} )}},\end{equation}
where 
${{\varGamma }_{\infty }}$ is the maximum surface excess due to steric considerations. Substituting (5.1), (5.2) and (5.3) into the ratio of velocities gives
\begin{equation} \frac{U}{{{U_d}}} \sim \left( {\frac{{\left( {\dfrac{{{R_1}{\ell _I}}}{{{\ell _s}^2}}} \right)}}{{( {1 - {\varGamma/{{\varGamma _\infty }}}} )}}} \right)\left( {\frac{{\displaystyle \int_0^\infty {( {{{\rm e}^{{\varPhi _2}( y )/{k_b}T}} - 1} )}\, {{\rm d} y}}}{{\displaystyle \int_0^\infty {y( {{{\rm e}^{{\varPhi _1}( y )/{k_b}T}} - 1} )}\, {{\rm d} y}}}} \right). \end{equation}
Here 
${{\ell }_{I}}$ and 
${{\ell }_{S}}$ both have similar origins and are typically of the order of a few molecular lengths (Anderson, Lowell & Prieve Reference Anderson, Lowell and Prieve1982). We restrain our analysis to the dilute limit of the solute surface excess where 
${\varGamma }/{{{\varGamma }_{\infty }}} \to 0$ which is in line with our earlier assumption that the partition coefficient is independent of the bulk concentration. Further, Anderson et al. (Reference Anderson, Lowell and Prieve1982) showed that the ratio of the integrals in (5.4) is of 
$O( 1 )$ if we assume that the (typically attractive) potentials 
${{\varPhi }_{1}},{{\varPhi }_{2}}$ are similar in magnitude. The ratio of velocities is then directly proportional to the ratio of the particle size 
$\sim O( {{10}^{-6}}\ {\rm m})$ and the adsorption zone length scale which is typically 
$\sim O( {{10}^{-9}}\ {\rm m})$ giving 
${U}/{{{U}_{d}}} \sim O( {{10}^{3}} )$. Thus, it is evident that a Marangoni-based propulsion mechanism is a suitable alternative to self-diffusiophoresis as a mechanism to drive a colloidal motor based on Janus particles. It can be noted that while our calculations (see figures 4–11) do not suggest that the velocities go as high as 
$O$(cm s
$^{-1}$), surface-tension-driven systems can indeed realize such high velocities. In these cases, it is necessary to study the effect of the advection on the solute transport as the corresponding Péclet numbers (here defined in § 2.2.2) have been shown to exhibit non-negligible values. The inclusion of fluid inertia and/or advection can significantly alter the dynamics of the self-propelled particle as demonstrated in recent literature involving Marangoni surfers (Kang et al. Reference Kang, Sur, Rothstein and Masoud2020; Ender et al. Reference Ender, Froin, Rehage and Kierfeld2021).
Further research will be necessary to determine the dynamics of the system if the particle and the drop come into contact. As the particle approaches the drop interface (i.e. as 
$\lambda \rightarrow \pm 1$ ), a few different possibilities can occur. If the particle is strongly wet by the drop phase, then the interactions between the particle surface and the interface in the form of a disjoining pressure (De Gennes Reference De Gennes1985) may become important impeding the particle from breaching the interface. However, if the particle exhibits a finite contact angle at the drop interface, the particle is likely to straddle the interface. In the latter case, the solutions using the bispherical basis will no longer be valid and the analysis will need to be conducted using tangent sphere and/or toroidal coordinates (Danov, Dimova & Pouligny Reference Danov, Dimova and Pouligny2000; Tsemakh et al. Reference Tsemakh, Lavrenteva and Nir2004). In addition, it may be worthwhile to analyse the stability of the axisymmetric configuration studied here to non-axisymmetric perturbations. Further, the stability might be easier to examine for the lubrication limit which can be more interesting for applications as it can maintain its configuration for extended periods (due to the small relative velocity of the particle with respect to the drop) and has simpler explicit formulae for the leading-order compound drop/Janus particle velocity which allow for better control and predictability of the system.
Finally, the kinetics of the solute exchange at the interface have been assumed to fast relative to the diffusion in the bulk. (i.e. Equation (2.4)). For very small droplets, that is when the diffusion length scale is sufficiently small, the kinetic transport from and to the interface might exhibit time scales that are comparable to diffusion and hence might need to also be considered (see for example Alvarez, Walker & Anna (Reference Alvarez, Walker and Anna2010)).
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2023.5.
Funding
J.K. and C.M. were supported in part by NSF grants CBET-1805554 and CBET-1929502.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Lubrication analysis in the limit of $\tau \to 0$
Consider the compound drop/Janus particle system in the limiting case wherein the drop and the Janus particle are nearly identical in radius. In this limit, there exists a small parameter; namely, the dimensionless thickness of the inner liquid layer. The Janus particle can adopt eccentric positions such that the relative distance between the centres of the two spheres is as depicted in figure 12(a). Adopting a similar terminology to the lubrication analysis of concentric cylinders in Leal (Reference Leal2007), the eccentricity is given, with positive and negative values for defined in the same way as in the bispherical solutions in the paper. The problem is best solved using the spherical coordinate system with the origin at the centre of the Janus particle and moving with the velocity of the drop as shown in figure 12(b) (This is possible as the solutions are quasi steady.) The dimensional radial coordinate of the liquid interface is then given by 
$R_{I}(\theta)$ and the equation for the interface can be stated simply as
\begin{equation} H = r - {R_I}( \theta ) = 0.\end{equation}
Figure 12. (a) Graphical depiction of lubrication limit; (b) schematic diagram for lubrication analysis.
Using this, the normal and tangential unit vectors can be written as
\begin{gather} \hat{\boldsymbol n} = \frac{{ - \boldsymbol{\nabla} H}}{{|\boldsymbol{\nabla} H|}} ={-} \left( {\frac{1}{{\sqrt {1 + {{\left( {\left( {\dfrac{1}{{1 + h}}} \right)\dfrac{{\partial h}}{{\partial \theta }}} \right)}^2}} }}{{{\hat{\boldsymbol e}_r}} } - \frac{{\left( {\left( {\dfrac{1}{{1 + h}}} \right)\dfrac{{\partial h}}{{\partial \theta }}} \right)}}{{\sqrt {1 + {{\left( {\left( {\dfrac{1}{{1 + h}}} \right)\dfrac{{\partial h}}{{\partial \theta }}} \right)}^2}} }}{{\hat{\boldsymbol e}_\theta} }} \right), \end{gather}
\begin{gather}\hat{\boldsymbol t} = \left( {\frac{{\left( {\left( {\dfrac{1}{{1 + h}}} \right)\dfrac{{\partial h}}{{\partial \theta }}} \right)}}{{\sqrt {1 + {{\left( {\left( {\dfrac{1}{{1 + h}}} \right)\dfrac{{\partial h}}{{\partial \theta }}} \right)}^2}} }}{{{\hat{\boldsymbol e}_r}} } + \frac{1}{{\sqrt {1 + {{\left( {\left( {\dfrac{1}{{1 + h}}} \right)\dfrac{{\partial h}}{{\partial \theta }}} \right)}^2}} }}{{\hat{\boldsymbol e}_\theta} }} \right). \end{gather}
Note that a negative sign is used to be consistent with our earlier definition of the unit normal (pointing into fluid 1). Scaling length with 
${{R}_{1}}$, velocity by U and pressure by a characteristic scale 
${{P}_{0}}$ (yet to be determined)
\begin{equation} h = \tau ( {1 - \lambda \cos \theta } ).\end{equation}
Here the drop is assumed to retain its spherical shape. If the drop were allowed to deform, then we would have to solve the kinematic boundary conditions to determine 
$h( \theta )$. The normal vectors have both 
$O( 1 )$ and 
$O( \tau )$ terms,
\begin{align}
\hat{\boldsymbol n} &=- \left( {\frac{1}{{\sqrt {1 +
{{\left( {\dfrac{{\tau \lambda \sin \theta }}{{1 + \tau (
{1 - \lambda \cos \theta } )}}} \right)}^2}}
}}{{\hat{\boldsymbol e}_r}} + \frac{{\left( {\dfrac{{\tau
\lambda \sin \theta }}{{1 + \tau ( {1 - \lambda \cos \theta
} )}}} \right)}}{{\sqrt {1 + {{\left( {\dfrac{{\tau \lambda
\sin \theta }}{{1 + \tau ( {1 - \lambda \cos \theta } )}}}
\right)}^2}} }}{{\hat{\boldsymbol e}_\theta} }}
\right)\nonumber\\
& =- {\hat e_r} + O( \tau ), \end{align}
\begin{align}\hat{\boldsymbol
t} &= \left( { - \frac{{\left( {\dfrac{{\tau \lambda \sin
\theta }}{{1 + \tau ( {1 - \lambda \cos \theta } )}}}
\right)}}{{\sqrt {1 + {{\left( {\dfrac{{\tau \lambda \sin
\theta }}{{1 + \tau ( {1 - \lambda \cos \theta } )}}}
\right)}^2}} }}{{{\hat{\boldsymbol e}_r}} } +
\frac{1}{{\sqrt {1 + {{\left( {\dfrac{{\tau \lambda \sin
\theta }}{{1 + \tau ( {1 - \lambda \cos \theta } )}}}
\right)}^2}} }}{{\hat{\boldsymbol e}_\theta} }} \right)
\nonumber\\
&= {{\hat{\boldsymbol e}_\theta} } + O(
\tau ),
\end{align}
taking the scaled velocities 
${\boldsymbol {{u}_{i}}}={{u}_{ri}}{{\hat {\boldsymbol {e}}_{r}}}+{{u}_{\theta i}}{{\hat {\boldsymbol {e}}_{\theta }}} , \{ i\in 1,2 \}$. The governing equations of the continuity equation and the individual components of the momentum balance equations for either fluid are given by
\begin{gather} \frac{1}{{{r^2}}}\frac{\partial }{{\partial r}}( {{r^2}{u_{ri}}} ) + \frac{1}{{r\sin \theta }}\frac{\partial }{{\partial \theta }}( {\sin \theta {u_{\theta i}}} ) = 0, \end{gather}
\begin{gather}\left( {\frac{{{P_0}{R_1}}}{{{\mu _1}U}}} \right)\left( { - \frac{{\partial {p_i}}}{{\partial r}}} \right) + \left( \begin{aligned} & \frac{1}{{{r^2}}}\frac{\partial }{{\partial r}}\left( {{r^2}\frac{{\partial {u_{ri}}}}{{\partial r}}} \right) + \frac{1}{{{r^2}\sin \theta }}\frac{\partial }{{\partial \theta }}\left( {\sin \theta \frac{{\partial {u_{ri}}}}{{\partial \theta }}} \right)\\ & \quad - \frac{{2{u_{ri}}}}{{{r^2}}} - \frac{2}{{{r^2}}}\frac{{\partial {u_{\theta i}}}}{{\partial \theta }} - \frac{{2{u_{\theta i}}\cot \theta }}{{{r^2}}}\end{aligned}\right) = 0, \end{gather}
\begin{gather}\left( {\frac{{{P_0}{R_1}}}{{{\mu _1}U}}} \right)\left( { - \frac{{\partial {p_i}}}{{\partial r}}} \right) + \left( \begin{aligned} & \frac{1}{{{r^2}}}\frac{\partial }{{\partial r}}\left( {{r^2}\frac{{\partial {u_{\theta i}}}}{{\partial r}}} \right) + \frac{1}{{{r^2}\sin \theta }}\frac{\partial }{{\partial \theta }}\left( {\sin \theta \frac{{\partial {u_{\theta i}}}}{{\partial \theta }}} \right)\\ & \quad + \frac{2}{{{r^2}}}\frac{{\partial {u_{ri}}}}{{\partial \theta }} - \frac{{{u_{\theta i}}}}{{{r^2}{{\sin }^2}\theta }} \end{aligned} \right) = 0. \end{gather}
The boundary conditions in this moving reference frame can be written as follows: at 
$r =1$,
\begin{gather} {u_{r1}} = \left( {\frac{V}{U} - 1} \right)( {{\boldsymbol{{\hat e}_z}} \boldsymbol{\cdot} {\boldsymbol{{\hat e}_r}}} ), \end{gather}
\begin{gather}{u_{\theta 1}} = \left( {\frac{V}{U} - 1} \right)( {{\boldsymbol{{\hat e}_z}} \boldsymbol{\cdot} {\boldsymbol{{\hat e}_\theta} }} ); \end{gather}
at 
$r=1+h( \theta )$,
\begin{gather} {u_{r1}} = {u_{r2}}, \end{gather}
\begin{gather}{u_{\theta 1}} = {u_{\theta 2}}, \end{gather}
\begin{gather}{\boldsymbol{u_{1}}} \boldsymbol{\cdot} \hat{\boldsymbol n} = 0, \end{gather}
\begin{gather}\mu \left( {r\frac{\partial }{{\partial r}}\left( {\frac{{{u_{\theta 2}}}}{r}} \right) + \frac{1}{r}\frac{{\partial {u_{r2}}}}{{\partial \theta }}} \right) - \left( {r\frac{\partial }{{\partial r}}\left( {\frac{{{u_{\theta 1}}}}{r}} \right) + \frac{1}{r}\frac{{\partial {u_{r1}}}}{{\partial \theta }}} \right) = \left( {\frac{{{k_2}{\gamma _0}\beta }}{{{\mu _1}U}}\frac{{J{R_1}}}{{{D_1}}}} \right)\left( {\frac{1}{r}} \right)\frac{{\partial {C_2}}}{{\partial \theta }}. \end{gather}
Following the approach in Kang, Nadim & Chugunova (Reference Kang, Nadim and Chugunova2016) for a thin viscous coating over a sphere, we define a secondary length scale y as 
$r=1+\tau y$ for the inner fluid. This length scale is the appropriate length scale for the radial direction in the inner fluid as 
$\tau \to 0$. The effects of the lubrication forces as are dominant only in the inner fluid as it is confined in between the spherical drop surface and the Janus particle. Further, to simplify the problem analytically, we can restrict 
$( k,D,\mu )\sim O( 1 )$. Therefore, the solutes are free to traverse through the interface and does not accumulate in the inner fluid irrespective of the drop thickness. Thus, the lubrication effects are only present in the fluid mechanics of the problem. In balancing the order of the terms in the continuity equation, we get 
${{u}_{r1}}=O( \tau )$. Rewriting the velocity explicitly in terms of the thickness
\begin{equation} {u_{r1}} = \tau {u'_{r1}}. \end{equation}
The 
$\theta$ component of the momentum balance also reveals that 
${{P}_{0}}={{{\mu }_{1}}U}/{{{R}_{1}}{{\tau }^{2}}}$. The governing equations reduce to
\begin{gather} \frac{{{\rm d}{{u'}_{r1}}}}{{{\rm d} y}} + \frac{1}{{\sin \theta }}\frac{{\rm d}}{{{\rm d}\theta }}( {{u_{\theta 1}}\sin \theta }) = 0, \end{gather}
\begin{gather}\frac{{\partial {p_1}}}{{\partial \theta }} = \frac{{{{\rm d}^2}{u_{\theta 1}}}}{{{\rm d}{y^2}}}, \end{gather}
\begin{gather}\frac{{\partial {p_1}}}{{\partial y}} = 0. \end{gather}
These are similar to those derived by Kang et al. (Reference Kang, Nadim and Chugunova2016). Before discussing the boundary conditions, it is necessary to establish the scaling in the outer fluid. In this region, the fluid is not confined by the interface and will not exhibit lubrication effects. Thus, the natural length scale for the radial and angular directions remains 
$O( 1 )$. As dictated by mass continuity condition, the velocity fields in these directions are also scaled by a common uniform scaling which as of yet remains unknown. However, this also means that the scaled equations are identical to the dimensionless Stokes equations. Thus, the general solution for the axisymmetric stream function in spherical coordinates (3.17) can be adopted to satisfy the boundary conditions. This is possible despite the origin of the basis being located at the centre of the Janus particle as the drop surface is still a constant coordinate surface to leading order.
The boundary conditions to leading order at the Janus particle surface are
\begin{gather} {u'_{r1}} = V'\cos \theta , \end{gather}
\begin{gather}{u_{\theta 1}} = 0. \end{gather}
Here 
${V}'=({1}/{\tau })({V}/{U}-1 )$. At 
$y=1-\lambda \cos \theta$, the boundary conditions to leading order along the aforementioned normal and tangential unit vectors (A5), (A6) are
\begin{gather} {u'_{r1}} = ( {\lambda \sin \theta } ){u_{\theta 1}}, \end{gather}
\begin{gather}{u_{\theta 1}} = {u_{\theta 2}}, \end{gather}
\begin{gather}\frac{{\partial {u_{\theta 1}}}}{{\partial y}} ={-} \left( {\frac{{{k_2}{\gamma _0}\beta J{R_1}}}{{{\mu _1}U'{D_1}}}} \right)\left( {\frac{{\partial {C_2}}}{{\partial \theta }}} \right) ={-} T( \theta ). \end{gather}
The angular component of the velocity field in the outer fluid has the same scaling as that of the inner fluid as necessitated by its continuity at the liquid interface. Since the velocity components in the outer fluid were shown to have a uniform scale, the scaled boundary condition at infinity remains 
${\boldsymbol {{u}_{2}}}\to -1 {\boldsymbol {{e}_{z}}}$ as necessitated by the moving reference frame adopted from the very beginning. Further, the pressure is scaled by the natural scaling 
$( {{{\mu }_{2}}U}/{{{R}_{1}}})$. Notice that the tangential stress balance contains only the derivative of the velocity field from the inner fluid. This comes out naturally from the earlier ideas of scaling. Intuitively, this is expected as the gradients of the velocity field normal to the interface in the inner fluid must be balanced by the surface tension forces as the absence of a confining surface disallows such strong gradients in the outer fluid. The equations derived here are like those derived in an earlier study (Kang, Nadim & Chugunova Reference Kang, Nadim and Chugunova2017) with an additional one way coupling between the angular component of the outer and inner fluid velocities. However, since the inner fluid is completely defined even without this condition, the solution procedure is fairly straightforward. Integrating the governing equation(s) along with the mass continuity equation and applying the appropriate boundary conditions give
\begin{equation} {u_{\theta 1}} = \frac{{{\rm d}{p_1}}}{{{\rm d}\theta }}\left( {\frac{{{y^2}}}{2} - ( {1 - \lambda \cos \theta } )y} \right) - T( \theta )y, \end{equation}
\begin{align} {{u'}_{r1}}&= \int { - \frac{1}{{\sin \theta }}\frac{{\rm d}}{{{\rm d} \theta }}({{u_{\theta 1}}\sin \theta } )\,{{\rm d} y}} \nonumber\\ &= \left( { - \frac{1}{{\sin \theta }}\frac{{\rm d}}{{{\rm d} \theta }}\left( {\sin \theta \left( {\frac{{{\rm d}{p_1}}}{{{\rm d} \theta }}\left( {\frac{{{y^3}}}{6} - ( {1 - \lambda \cos \theta } )\frac{{{y^2}}}{2}} \right) - T( \theta )\frac{{{y^2}}}{2}} \right)} \right)} \right) + V'\cos \theta. \end{align}
It is convenient to define the following functions: 
$( 1-{{\vartheta }^{2}} )({{\rm d}{{p}_{1}}}{{\rm d}\vartheta })=K,( 1- {{\vartheta }^{2}} )^{1/2}T( \theta ) =Q$, where 
$\vartheta =\cos \theta$ as defined in the earlier analysis of the concentric configuration. It follows from the axisymmetric nature of the problem that 
$( K,Q=0 ), \vartheta =( -1,1 )$. Applying the remaining boundary conditions, we get
\begin{equation} \frac{{{{( {1 - \lambda \vartheta } )}^3}}}{3}\frac{{{\rm d}K}}{{{\rm d}\vartheta }} - \lambda K{( {1 - \lambda \vartheta } )^2} + V'\vartheta = \frac{{\rm d}}{{{\rm d} \vartheta }}\left( {\frac{{Q{{( {1 - \lambda \vartheta } )}^2}}}{2}} \right).\end{equation}The complete solution for the coefficients of the stream function for the outer fluid is straightforward once the boundary conditions are written in terms of the stream function (see (3.18)–(3.20)). However, for the purposes of our calculation, it turns out that only the coefficients proportional to the first modified Gegenbauer polynomial are necessary. Applying (3.14) and the force-free condition (2.13) we can eliminate the coefficients in the continuity of the angular component of velocity boundary condition to get (to leading order)
\begin{align} &\int_{ - 1}^1 {\left( {\frac{{{\rm d}{p_1}}}{{{\rm d}\vartheta }}\frac{{{{( {1 - \lambda \vartheta } )}^2}}}{2}( {1 - {\vartheta ^2}} ) - T( {1 - \lambda \vartheta } ){{( {1 - {\vartheta ^2}} )}^{1/2}}} \right)\, {\rm d} \vartheta } \nonumber\\ &\quad = \int_{{-}1}^1 {\left( {K\frac{{{{( {1 - \lambda \vartheta } )}^2}}}{2} - Q( {1 - \lambda \vartheta } )} \right)\, {\rm d} \vartheta } = 2. \end{align}Finally, the force-free condition for the Janus particle reduces to
\begin{equation} \boldsymbol{F} \boldsymbol{\cdot} {\hat{\boldsymbol e_z}} = \int_0^{2{\rm \pi} } {\int_0^{\rm \pi} {( {{\boldsymbol{{\hat e}_r}} \boldsymbol{\cdot} ( { - {p_1}{\boldsymbol{I} } + {\boldsymbol{\sigma _1}}}) \boldsymbol{\cdot} {\boldsymbol{{\hat e}_z}}} )\sin \theta } \,{\rm d}\theta \,} {\rm d}\phi = 0,\end{equation}
where 
$\boldsymbol {I}$ is the identity tensor. Here upon scaling, the force assumes a natural scale from the product of the pressure scaling and length squared 
$( {{{\mu }_{1}}U{{R}_{1}}}/{{{\tau }^{2}}} )$ and the viscous stress does not contribute to the leading order in 
$\tau$. The force-free condition can be simplified using the axisymmetry conditions to the integral 
$\int _{-1}^{1}{K}\,{\rm d}\vartheta =0$. Integrating and simplifying (A29), we get the following relations for the terminal velocity after a bit of algebra:
\begin{gather} \frac{{U'}}{{\dfrac{{{k_2}{\gamma _0}\beta J{R_1}}}{{{\mu _1}{D_1}}}}} = \frac{1}{8}\left( {6B'\left( {\frac{{( {\lambda + ( {{\lambda ^2} - 1} )\,{\rm arctanh}\,\lambda } )}}{{{\lambda ^3}}}} \right) - A'} \right), \end{gather}
\begin{gather}V' = \left( {\frac{{{k_2}{\gamma _0}\beta J{R_1}}}{{{\mu _1}{D_1}U'}}} \right)B'. \end{gather}With the constants defined as
\begin{gather} A' ={-} \left( {\sum_{n = 1}^\infty {\left( {\left( {\frac{{n\left( {2n + 1} \right)}}{2}\left( {\frac{1}{D}} \right)\int_1^{\cos \alpha } {{p_n}( \vartheta )\, {\rm d} \vartheta } } \right)\int_{ - 1}^1 {{g_n}( \eta )( {1 - \lambda \eta } )\, {\rm d}\eta } } \right)} } \right), \end{gather}
\begin{gather}B' ={-} \left( \begin{aligned} & \left( {\frac{{{\lambda ^3}}}{{2\left( {\dfrac{\lambda }{{ - 1 + {\lambda ^2}}} + {\rm Tanh}^{ - 1}\lambda } \right)}}} \right)\\ & \quad \times \sum_{n = 1}^\infty {\left( {\left( {\frac{{n( {2n + 1} )}}{2}\left( {\frac{1}{D}} \right)\int_1^{\cos \alpha } {{p_n}( \vartheta )\, {\rm d} \vartheta } } \right)\int_{ - 1}^1 {\frac{{{g_n}( \eta )}}{{( {1 - \lambda \eta } )}}\, {\rm d}\eta } } \right)} \end{aligned}\right), \end{gather}
where have used the Legendre polynomial expansion for the outer fluid concentration to leading order. Since 
$( {V}/{U}-1 )\sim O( {{\tau }^{2}} )$ from definition of 
${V}'$, the dimensional leading-order solution to the terminal velocities for both the Janus particle and the drop is given by
\begin{equation} U \sim V \sim {\tau U' + O( {{\tau ^2}} )}.\end{equation}Thus the Janus particle and the drop can move together with a common terminal velocity to leading order.
Appendix B. Method of reflection analysis for the limit of $\tau \to \infty$ while retaining the leading-order effect of the size of the Janus particle
Since the problems of species transport and the momentum transport are completely linear, the solutions can be written as a sum of infinitely many reflections 
$({}^{i}{{M}_{j}} ), ( M\in \{ C,\psi \}, i\in \{ 0,\infty \},j\in \{ 1,2 \} )$ where the 
$M$ represents the field variable, the superscript represents the index of the reflection and the subscript represents the fluid. The expression 
${}^{n}{{C}_{1}}$is the 
$n$th reflection of the concentration field in the inner fluid. In this notation, the complete solutions to the concentration and stream functions can be written as
\begin{gather} {C_i} = \sum_{j = 0}^\infty {{}^j{C_i}} , \end{gather}
\begin{gather}{\psi _i} = \sum_{j = 0}^\infty {{}^j{\psi _i}} , \end{gather}
where each successive reflection is solved in spherical coordinates with the same 
${{\hat {\boldsymbol {e}}}_{\phi }}$ and the origin alternatingly at the Janus particle 
$( {{r}_{p}},{{\theta }_{p}} )$ or the drop centres 
$( {{r}_{d}},{{\theta }_{d}} )$. The conversion between the two can be done using the following relations:
\begin{align}
\left.\begin{array}{c} {r_d} = \sqrt {{{( {\tau
{R_1}\lambda } )}^2} + 2{r_p}( {\tau {R_1}\lambda } )\cos
{\theta _p} + {r_p}^2} , \quad \cos {\theta _d} = \dfrac{{(
{\tau {R_1}\lambda } ) + {r_p}\cos {\theta _p}}}{{\sqrt
{{{( {\tau {R_1}\lambda } )}^2} + 2{r_p} ( {\tau
{R_1}\lambda } )\cos {\theta _p} + {r_p}^2} }}\\ {r_p} =
\sqrt {{{( {\tau {R_1}\lambda } )}^2} - 2{r_d}( {\tau
{R_1}\lambda } )\cos {\theta _d} + {r_d}^2} ,\quad \cos
{\theta _p} = \dfrac{{{r_d}\cos {\theta _d} - ( {\tau
{R_1}\lambda } )}}{{\sqrt {{{( {\tau {R_1}\lambda } )}^2} -
2{r_d} ( {\tau {R_1}\lambda } )\cos {\theta _d} + {r_d}^2}
}} \end{array}
\right\}.\end{align}
We reiterate that the eccentric location of the Janus particle is given by the expression 
$( \tau {{R}_{1}}\lambda )$, but with the majority of the configurations, the centre-to-centre distance scales with the drop radius. Thus, we can define a scaled variable 
${{d}_{p}}=( \tau \lambda \chi )$ which remains finite in the limits of 
$\tau \to \infty,\chi \to 0$. Therefore, it is convenient to work with 
$( {{R}_{2}}{{d}_{p}} )$ instead of 
$( \tau {{R}_{1}}\lambda )$. The goal of this calculation is to determine the leading-order correction to the earlier computed terminal velocities from the point source formulation.
B.1. Concentration field
The linearity of the equations governing the concentration field and their boundary conditions enables us to independently satisfy the governing equation with each reflection,
\begin{equation} {\nabla^2}( {{}^i{C_1}} ) = 0,\quad i \in \{ {0,\infty } \}.\end{equation}It is convenient to begin with the unbounded particle problem which is scaled with the length scale with the following boundary condition(s):
at 
${{r}_{p}}=1$,
\begin{equation} - \frac{{\partial ( {{}^0{C_1}} )}}{{\partial {r_p}}} = \left\{ \begin{array}{cc} 1, & { ( {{\theta _p} \le \alpha } )}\\ 0, & {( {{\theta _p} > \alpha } )} \end{array} \right . \end{equation}
The concentration field is scaled as 
$JR_1/D_1$. Solving this system using an expansion in Legendre polynomials, we get the solution
\begin{equation} {}^0{C_1} = \sum_{n = 0}^\infty {\left( {\frac{{{}^0{Y_n}}}{{{r_p}^{n + 1}}}} \right)} {p_n}( {\cos {\theta _p}} ), \end{equation}where the solutions to the coefficients can be obtained using the orthogonality properties of the Legendre polynomials (3.8). Clearly, the reflections emanating from the Janus particle surface do not have a corollary in the outer fluid and are defined only for the inner fluid. However, they do contribute to the concentration field in the outer fluid as the inhomogeneity from these reflections enable nontrivial solutions to the corresponding reflections from the drop surface which are indeed defined in both the outer and inner fluids.
The first reflection is evaluated at the drop surface with the aid of the general solutions as defined in § 3 by satisfying the necessary boundary conditions in the presence of the contributions from the original solution. The corresponding boundary conditions for the reflection from the drop surface without the particle are as follows:
at 
${{r}_{d}}=1$,
\begin{gather} ( {{}^1{C_1} + {}^0{C_1}} ) = k( {{}^1{C_2}} ), \end{gather}
\begin{gather}\frac{{\partial ( {{}^1{C_1} + {}^0{C_1}} )}}{{\partial {r_d}}} = D\frac{{\partial ( {{}^1{C_2}} )}}{{\partial {r_d}}}, \end{gather}
where we need to rescale the problem with the length scale 
${{R}_{2}}$ in order to apply the boundary conditions at the drop surface. The concentration field retains the earlier scale. The solutions from 
${}^{0}{{C}_{1}}$ for the inner fluid can be expanded in the small parameter 
$\varepsilon =\chi$ as
\begin{equation} {}^0{C_1} = \sum_{n = 0}^\infty {( {{\varepsilon ^{n + 1}}} )\left( {\frac{{{}^0{Y_1}}}{{{{( {{r_d}^2 \!+\! {d_p}^2 \!+\! 2{r_d}{d_p}\cos {\theta _d}} )}^{({{n + 1}})/{2}}}}}} \right)} {p_n}\left( {\frac{{{r_d}\cos {\theta _d} \!-\! {d_p}}}{{\sqrt {{d_p}^2 \!-\! 2{r_d}{d_p}\cos {\theta _d} + {r_d}^2} }}} \right). \end{equation}The corresponding solutions for the concentration field are given by
\begin{gather} {}^1{C_1} = \sum_{n = 0}^\infty {( {{}^1{X_n}{r_d}^n} )} {p_n}( {\cos {\theta _d}} ), \end{gather}
\begin{gather}{}^1{C_2} = \sum_{n = 0}^\infty {\left( {\frac{{{}^1{y_n}}}{{{r_d}^{n + 1}}}} \right)} {p_n}( {\cos {\theta _d}} ). \end{gather}
It is convenient to further expand the coefficients of the concentration field into a power series in the small parameter, i.e. 
$( {}^{j}{{M}_{n}} )=\sum _{k=0}^{\infty }{( {{\varepsilon }^{k}} )}( {}^{j}{{M}_{n}}^{k} ), ( M\in \{ X,Y,x,y \} )$. The solutions to the equations are then calculable using the above equations by matching orders of 
$\varepsilon$ (note, the subscript for the coefficient denotes the order of the Legendre polynomial to which it is a coefficient and is not an index representing the inner/outer fluid like in the case of the other field variables).
It is evident that the solution from the particle surface is at least an order of magnitude smaller at the surface of the drop. Since 
${}^{0}{{C}_{1}}$ is the only inhomogeneity in the equations for the first reflection, the solutions to 
${}^{1}{{C}_{1}},{}^{1}{{C}_{2}}$ are to leading order at most 
$O( \varepsilon )$. Similarly, it can by shown that the contributions from the further reflections onto the drop surface are either constant at the surface or at most 
$O( {{\varepsilon }^{3}} )$. Since the velocity fields are driven by the tangential gradients of the concentration field at the drop surface, it can be shown that no further reflections need to be computed to get the required first-order correction to the terminal velocities. The coefficients of the concentration field at the drop surface for the outer fluid can then be calculated to within a constant to the order of 
$O({\varepsilon }^{2})$,
\begin{equation} {C_2}( {{\theta _d}} ){|_{{r_d} = 1}} = \sum_{n = 0}^\infty {( {( \varepsilon )( {{}^1{y_n}^1} ) + ( {{\varepsilon ^2}} )( {{}^1{y_n}^2} )} )} {p_n}( {\cos {\theta _d}} ) + O( {{\varepsilon ^3}} ) + {{\rm const.}} \end{equation}B.2. Velocity field
The problem is solved in a stationary reference frame. The velocity field is then calculated using the axisymmetric stream functions in the same manner as in the concentric case solutions discussed in the paper,
\begin{equation} {\boldsymbol{u_i}} = \boldsymbol{\nabla} \times \left( {\frac{1}{\omega }\sum_{j = 0}^\infty {{}^j{\psi _i}}{{\hat{\boldsymbol{e}}}_{\phi }} } \right), \end{equation}
with 
$\omega ={{r}_{d}}\sin {{\theta }_{d}}={{r}_{p}}\sin {{\theta }_{p}}$. To allow for successive corrections to the terminal velocities, the unknown velocities must be expanded in the small parameter,
\begin{gather} \frac{U}{{{U_0}}} = \left( {1 + \sum_{k = 1}^\infty {{\varepsilon ^k}{U_k}} } \right), \end{gather}
\begin{gather}\frac{V}{{{U_0}}} = \left( {\sum_{k = 0}^\infty {{\varepsilon ^k}{V_k}} } \right). \end{gather}
The natural scale for the velocity in all reflections is taken to be 
${{U}_{0}}$ which automatically implies the varying scales 
$( {{U}_{0}}{{R}_{1}}^{2},{{U}_{0}}{{R}_{2}}^{2} )$ for the stream functions in reflections originating from the Janus particle and the drop surface, respectively. It is convenient to begin the calculations with an original solution for either fluid 
$( {}^{0}{{\psi }_{i}} )$ from the drop surface and then to take successive reflections from the Janus particle and the drop alternatingly. It is to be noted that when taking the reflection from the surface of the Janus particle, the inhomogeneity from the drop surface gets rescaled from 
${{R}_{2}}$ to 
${{R}_{1}}$. In this process, it is convenient to expand out the coordinates in the small parameter as follows:
\begin{gather} \frac{{{r_d}}}{{{R_2}}} = \varepsilon \frac{{{r_d}}}{{{R_1}}} = \sqrt {{d_p}^2 + 2\varepsilon \left( {\frac{{{r_p}}}{{{R_1}}}} \right){d_p}\cos {\theta _p} + {\varepsilon ^2}{{\left( {\frac{{{r_p}}}{{{R_1}}}} \right)}^2}} \approx {d_p} + \varepsilon \left( {\frac{{{r_p}}}{{{R_1}}}} \right)\cos {\theta _p} + O( {{\varepsilon ^2}} ), \end{gather}
\begin{gather}\cos {\theta _d} = \frac{{{d_p} + \varepsilon \left( {\frac{{{r_p}}}{{{R_1}}}} \right)\cos {\theta _p}}}{{\sqrt {{d_p}^2 + 2\varepsilon \left( {\frac{{{r_p}}}{{{R_1}}}} \right){d_p}\cos {\theta _p} + {\varepsilon ^2}{{\left( {\frac{{{r_p}}}{{{R_1}}}} \right)}^2}} }} \approx 1 - \frac{{{\varepsilon ^2}}}{2}\frac{{( {{{\sin }^2}{\theta _p}} )}}{{{d_p}^2}}{\left( {\frac{{{r_p}}}{{{R_1}}}} \right)^2} + O( {{\varepsilon ^3}} ). \end{gather}
Since the inhomogeneity is a stream function (or a concentration field), it will be convenient to express the Gegenbauer polynomials as a Taylor series about zero in the small parameter 
$\varepsilon$, as follows:
\begin{equation} {g_n}( {\cos {\theta _d}} ) = {g_1}( {\cos {\theta _p}} ){\left( {\frac{{{r_p}}}{{{d_p}{R_1}}}} \right)^2}( {{\varepsilon ^2}} ) - 2{g_1}( {\cos {\theta _p}} ){\left( {\frac{{{r_p}}}{{{d_p}{R_1}}}} \right)^3}( {{\varepsilon ^3}} ) + O( {{\varepsilon ^4}}), \end{equation}
where we have used 
${{g}_{n}}( 1 )=0,{{p}_{n}}( 1 )=1$ along with (3.13).
As usual, the tangential stress balance at the drop surface couples the velocity field to the earlier evaluated concentration field. Since the concentration field is at most of 
$O( \varepsilon )$ (see (B12)), the scale of the velocity must be in itself 
$O( \varepsilon )$, i.e. 
${{U}_{0}}=\varepsilon {{{U}'}_{0}}$ with 
${{{U}'}_{0}}\sim ( {{{k}_{2}}\beta {{\gamma }_{0}}J{{R}_{1}}}/{{{\mu }_{1}}{{D}_{1}}} )\sim O( 1 )$.
The general solutions to the stream functions assume the form (3.12). The coefficients corresponding to the 
${{n}^{th}}$ Gegenbauer polynomials of the eigenfunction expansions for the stream function 
$( {}^{j}A_{n}^{k},{}^{j}B_{n}^{k},{}^{j}C_{n}^{k},{}^{j}D_{n}^{k},{}^{j}a_{n}^{k},{}^{j}b_{n}^{k},{}^{j}c_{n}^{k},{}^{j}d_{n}^{k} )$ are further expanded in terms of the small parameter like in the case of the concentration field and they can be evaluated by satisfying the boundary conditions pertinent to each reflection. The force-free conditions can be employed at each order of the small parameter to obtain the corresponding correction to the terminal velocities. In particular, we can show that the force-free conditions reduce to the following conditions:
\begin{gather} O( 1 ):{}^1C_1^0 = {}^0c_1^0 = 0, \end{gather}
\begin{gather}O( \varepsilon ):{}^1C_1^1 = {}^0c_1^1 = 0. \end{gather}Using these conditions, we can estimate the dimensional terminal velocities to within the first correction for the particle size giving
\begin{gather} U = \varepsilon {U'_0} + {\varepsilon ^2}{U'_0}{U_1}, \end{gather}
\begin{gather}V = \varepsilon {U'_0}{V_0} + {\varepsilon ^2}{U'_0}{V_1}. \end{gather} The calculations reveal that under force-free conditions, the contributions of the reflections 
$( {}^{0}{{\psi }_{1}},{}^{0}{{\psi }_{2}},{}^{1}{{\psi }_{1}} )$ are sufficient to evaluate the required 
$O( \varepsilon )$ correction to the terminal velocities of the drop and the Janus particle. Further reflections of the concentration field and the stream function need to be evaluated to get higher-order corrections to the terminal velocities. The actual steps of the calculation are identical in spirit to that of the concentration field and are being omitted herewith in the interest of brevity.
