Swimming with a cage: Low-Reynolds-number locomotion inside a droplet

TL;DR

This study investigates how a spherical microswimmer inside a droplet can propel both itself and the droplet through hydrodynamic interactions, revealing stable co-swimming states and flow reversals relevant for targeted drug delivery.

Contribution

It provides an analytical and numerical analysis of droplet-encaged microswimmers, demonstrating stable co-movement, flow reversal effects, and stability scenarios based on swimmer gait.

Findings

Droplet can be propelled by the encaged swimmer.

Stable co-swimming states are achievable with specific actuation patterns.

Flow field around the swimmer is reversed by the droplet, affecting its pushers and pullers.

Abstract

Inspired by recent experiments using synthetic microswimmers to manipulate droplets, we investigate the low-Reynolds-number locomotion of a model swimmer (a spherical squirmer) encapsulated inside a droplet of comparable size in another viscous fluid. Meditated solely by hydrodynamic interactions, the encaged swimmer is seen to be able to propel the droplet, and in some situations both remain in a stable co-swimming state. The problem is tackled using both an exact analytical theory and a numerical implementation based on boundary element method, with a particular focus on the kinematics of the co-moving swimmer and droplet in a concentric configuration, and we obtain excellent quantitative agreement between the two. The droplet always moves slower than a swimmer which uses purely tangential surface actuation but when it uses a particular combination of tangential and normal actuations, the squirmer and droplet are able to attain a same velocity and stay concentric for all times. We next employ numerical simulations to examine the stability of their concentric co-movement, and highlight several stability scenarios depending on the particular gait adopted by the swimmer. Furthermore, we show that the droplet reverses the nature of the far-field flow induced by the swimmer: a droplet cage turns a pusher swimmer into a puller, and vice versa. Our work sheds light on the potential development of droplets as self-contained carriers of both chemical content and self-propelled devices for controllable and precise drug deliveries.

Tables

Table 1. Table 1: The constants for the fluid velocity field given in Eq.  13 and the power in Eq.  15 .

${\Delta }_n=\lambda \{{(2n+1)}^2({\chi }^{2n-1}-1)({\chi }^{2n+3}-1)$ $-(2n-1)(2n+3){({\chi }^{2n+1}-1)}^2\left\}\right(n\ge 2)$ $N_1=n\left(2n-1\right)\left\{\left({\chi }^{2n+1}-1\right)\lambda +1\right\}$ $-\left(n-2\right)\left\{\left(2n+1\right)\left({\chi }^{2n-1}-1\right)\lambda -2{\chi }^{2n-1}+2n+1\right\}$ $N_2=-2\left(2n-1\right)\left\{\left({\chi }^{2n+1}-1\right)\lambda +1\right\}$ $+2\left\{\left(2n+1\right)\left({\chi }^{2n-1}-1\right)\lambda -2{\chi }^{2n-1}+2n+1\right\}$ $N_3=\left(n-2\right)\left(2n+3\right)\left\{\left({\chi }^{2n+1}-1\right)\lambda +1\right\}$ $-n\left\{\left(2n+1\right)\left({\chi }^{2n+3}-1\right)\lambda +2{\chi }^{2n+3}+2n+1\right\}$ $N_4=-2\left(2n+3\right)\left\{\left({\chi }^{2n+1}-1\right)\lambda +1\right\}$ $+2\left\{\left(2n+1\right)\left({\chi }^{2n+3}-1\right)\lambda +2{\chi }^{2n+3}+2n+1\right\},$ $N_5={\chi }^{4n+2}[-(n+1)(2n+3)\{(1-{\chi }^{-2n-1})\lambda +1\}$ $+(n+3)\{(2n+1)(1-{\chi }^{-2n-3})\lambda +2{\chi }^{-2n-3}+2n+1\}]$ $N_6={\chi }^{4n+2}[-2(2n+3)\{(1-{\chi }^{-2n-1})\lambda +1\}$ $+2\{(2n+1)(1-{\chi }^{-2n-3})\lambda +2{\chi }^{-2n-3}+2n+1\}]$ $N_7={\chi }^{4n+2}[-(n+3)(2n-1)\{(1-{\chi }^{-2n-1})\lambda +1\}$ $+(n+1)\{(2n+1)(1-{\chi }^{-2n+1})\lambda -2{\chi }^{-2n+1}+2n+1\}]$ $N_8={\chi }^{4n+2}[-2(2n-1)\{(1-{\chi }^{-2n-1})\lambda +1\}$ $+2\{(2n+1)(1-{\chi }^{-2n+1})\lambda -2{\chi }^{-2n+1}+2n+1\}]$ $c_1=-\frac{n+3}{2}{N}_1{\chi }^{2n+3}-\frac{n+1}{2}{N}_3{\chi }^{2n+1}+\frac{n-2}{2}{N}_5{\chi }^2+\frac{n}{2}{N}_7$ $c_2=-\frac{n+3}{2}{N}_2{\chi }^{2n+3}-\frac{n+1}{2}{N}_4{\chi }^{2n+1}+\frac{n-2}{2}{N}_6{\chi }^2+\frac{n}{2}{N}_8$ $a_n=\left(\frac{2n+3}{n}-n-1\right)N_1-\left(n-1\right)N_3+\left(n+\frac{2n-1}{n+1}\right)N_5+\left(n+2\right)N_7$ $b_n=\left(\frac{2n+3}{n}-n-1\right)N_2-\left(n-1\right)N_4+\left(n+\frac{2n-1}{n+1}\right)N_6+\left(n+2\right)N_8$ $c_n=n\left(n+2\right)N_1+\left(n-1\right)\left(n+1\right)N_3+\left(n-1\right)\left(n+1\right)N_5+n\left(n+2\right)N_7$ $d_n=n\left(n+2\right)N_2+\left(n-1\right)\left(n+1\right)N_4+\left(n-1\right)\left(n+1\right)N_6+n\left(n+2\right)N_8$ $N_o=N_1+N_3+N_5+N_7$ $N_e=N_2+N_4+N_6+N_8$ ${\overline{N}}_o=-\frac{n+3}{2}{N}_1-\frac{n+1}{2}{N}_3+\frac{n-2}{2}{N}_5+\frac{n}{2}{N}_7$ ${\overline{N}}_e=-\frac{n+3}{2}{N}_2-\frac{n+1}{2}{N}_4+\frac{n-2}{2}{N}_6+\frac{n}{2}{N}_8$ $Z_1=2\left(2\lambda +3\right){\chi }^5-10\left(\lambda -1\right){\chi }^2+6\left(\lambda -1\right)$ $Z_2=\left(2\lambda +3\right){\chi }^5+10\left(\lambda -1\right){\chi }^2-12\left(\lambda -1\right)$ $Z_3=2\left\{\left(2\lambda +3\right){\chi }^5+3\left(\lambda -1\right)\right\}$

Equations

$$\nabla p^{\left(i\right)}={\mu }^{\left(i\right)}{\nabla }^2{\mathbf{v}}^{\left(i\right)},\nabla \cdot {\mathbf{v}}^{\left(i\right)}=0,$$

$$v_r^{\left(1\right)}{|}_{r=a}$$

$$v_{\theta }^{\left(1\right)}{|}_{r=a}$$

$$v_r^{\left(1\right)}{|}_{r=b}=v_r^{\left(2\right)}{|}_{r=b}=U_D\cos \theta ,$$

$$v_{\theta }^{\left(1\right)}{|}_{r=b}=v_{\theta }^{\left(2\right)}{|}_{r=b},$$

$${\Pi }_{r\theta }^{\left(1\right)}{|}_{r=b}={\Pi }_{r\theta }^{\left(2\right)}{|}_{r=b},$$

$$U_0=\frac{2B_1-A_1}{3}\cdot$$

$$\displaystyle\bm{v}=\sum_{n=-\infty}^{\infty}\Big{[}\nabla\phi_{n}+\frac{n+3}{2\mu(n+1)(2n+3)}r^{2}\nabla p_{n}-\frac{n}{\mu(n+1)(2n+3)}\bm{r}p_{n}\Big{]},$$

$$p_n(r,\xi )={\tilde{p}}_n{r}^n{P}_n(\xi ),{\varphi }_n(r,\xi )={\tilde{\varphi }}_n{r}^n{P}_n(\xi ),$$

$$\displaystyle v_{r}=\sum_{n\geq 0}^{\infty}\Big{[}\bar{p}_{n}r^{n+1}+\bar{\phi}_{n}r^{n-1}+\bar{p}_{-(n+1)}r^{-n}+\bar{\phi}_{-(n+1)}\frac{1}{r^{n+2}}\Big{]}P_{n}(\xi),$$

$$\displaystyle v_{\theta}=\sum_{n\geq 1}^{\infty}\Big{[}-\frac{n+3}{2}\bar{p}_{n}r^{n+1}-\frac{n+1}{2}\bar{\phi}_{n}r^{n-1}+\frac{n-2}{2}\bar{p}_{-(n+1)}{r^{-n}}+\frac{n}{2}\bar{\phi}_{-(n+1)}r^{-\left(n+2\right)}\Big{]}V_{n}(\xi),$$

$${\bar{p}}_n=\frac{n}{2\mu (2n+3)}{\tilde{p}}_n,{\bar{\varphi }}_n=n{\tilde{\varphi }}_n.$$

$${\bar{p}}_{-1}^{\left(1\right)}+\frac{1}{a^2}{\bar{\varphi }}_{-1}^{\left(1\right)}=A_0,{\bar{p}}_{-1}^{\left(1\right)}+\frac{1}{b^2}{\bar{\varphi }}_{-1}^{\left(1\right)}=0,$$

$$a^2{\bar{p}}_1^{\left(1\right)}+{\bar{\varphi }}_1^{\left(1\right)}+\frac{1}{a}{\bar{p}}_{-2}^{\left(1\right)}+\frac{1}{a^3}{\bar{\varphi }}_{-2}^{\left(1\right)}=A_1+U_S,$$

$$-2a^2{\bar{p}}_1^{\left(1\right)}-{\bar{\varphi }}_1^{\left(1\right)}-\frac{1}{2a}{\bar{p}}_{-2}^{\left(1\right)}+\frac{1}{2a^3}{\bar{\varphi }}_{-2}^{\left(1\right)}=B_1-U_S,$$

$$b^2{\bar{p}}_1^{\left(1\right)}+{\bar{\varphi }}_1^{\left(1\right)}+\frac{1}{b}{\bar{p}}_{-2}^{\left(1\right)}+\frac{1}{b^3}{\bar{\varphi }}_{-2}^{\left(1\right)}=U_D,$$

$$(-2-\frac{1}{\lambda })b^2{\bar{p}}_1^{\left(1\right)}-{\bar{\varphi }}_1^{\left(1\right)}-\frac{1}{2b}{\bar{p}}_{-2}^{\left(1\right)}+(\frac{1}{2}-\frac{1}{\lambda })\frac{1}{b^3}{\bar{\varphi }}_{-2}^{\left(1\right)}=-\frac{1}{2}{U}_D.$$

$$\displaystyle\bar{p}_{-2}^{\left(2\right)}=\frac{1}{b^{2}}\bar{\phi}_{-2}^{\left(2\right)}+bU_{D},\hskip 15.0pt\bar{\phi}_{-2}^{\left(2\right)}=\frac{1}{\lambda}\Big{(}b^{5}\bar{p}_{1}^{\left(1\right)}+\bar{\phi}_{-2}^{\left(1\right)}\Big{)},$$

$$\mathbf{F}={\int }_{\hat{S}}{\mathbf{\Pi }}^{\left(1\right)}\cdot \hat{\mathbf{r}}dS=-4\pi \nabla \left[r^3{p}_{-2}^{\left(1\right)}\right]=0,$$

$$U_S=\frac{{\Xi }_1\lambda +{\Xi }_2}{\Delta },$$

$$U_D=\frac{10(A_1+B_1){\chi }^2}{\Delta },$$

$${\Xi }_1=2(2B_1-A_1){\chi }^5-10(A_1+B_1){\chi }^2+6(2A_1+B_1),$$

$${\Xi }_2=3(2B_1-A_1){\chi }^5+10(A_1+B_1){\chi }^2-6(2A_1+B_1),$$

$$\Delta =3\left[2({\chi }^5-1)\lambda +3{\chi }^5+2\right].$$

$$\displaystyle v_{r}^{\left(1\right)}=\frac{A_{0}}{\chi^{2}-1}\Big{\{}\chi^{2}\Big{(}\frac{a}{r}\Big{)}^{2}-1\Big{\}}P_{0}(\xi)+\frac{A_{1}+B_{1}}{\Delta}\Big{\{}6(\lambda-1)\Big{(}\frac{r}{a}\Big{)}^{2}-10(\lambda-1)\chi^{2}+2(2\lambda+3)\chi^{5}\Big{(}\frac{a}{r}\Big{)}^{3}\Big{\}}P_{1}(\xi)$$

$$\displaystyle\hskip 25.0pt+\sum_{n=2}^{\infty}\frac{1}{\Delta_{n}}\Big{\{}(N_{1}A_{n}+N_{2}B_{n})\Big{(}\frac{r}{a}\Big{)}^{n+1}+(N_{3}A_{n}+N_{4}B_{n})\Big{(}\frac{r}{a}\Big{)}^{n-1}+(N_{5}A_{n}+N_{6}B_{n})\Big{(}\frac{a}{r}\Big{)}^{n}$$

$$\displaystyle\hskip 90.0pt+(N_{7}A_{n}+N_{8}B_{n})\Big{(}\frac{a}{r}\Big{)}^{n+2}\Big{\}}P_{n}(\xi),$$

$$\displaystyle v_{\theta}^{\left(1\right)}=-\frac{A_{1}+B_{1}}{\Delta}\Big{\{}12(\lambda-1)\Big{(}\frac{r}{a}\Big{)}^{2}-10(\lambda-1)\chi^{2}-(2\lambda+3)\chi^{5}\Big{(}\frac{a}{r}\Big{)}^{3}\Big{\}}V_{1}(\xi)$$

$$\displaystyle\hskip 25.0pt+\sum_{n=2}^{\infty}\frac{1}{\Delta_{n}}\Big{\{}-\frac{n+3}{2}(N_{1}A_{n}+N_{2}B_{n})\Big{(}\frac{r}{a}\Big{)}^{n+1}-\frac{n+1}{2}(N_{3}A_{n}+N_{4}B_{n})\Big{(}\frac{r}{a}\Big{)}^{n-1}$$

$$\displaystyle\hskip 90.0pt+\frac{n-2}{2}(N_{5}A_{n}+N_{6}B_{n})\Big{(}\frac{a}{r}\Big{)}^{n}+\frac{n}{2}(N_{7}A_{n}+N_{8}B_{n})\Big{(}\frac{a}{r}\Big{)}^{n+2}\Big{\}}V_{n}(\xi),$$

$$\displaystyle v_{r}^{\left(2\right)}=\frac{10(A_{1}+B_{1})\chi^{5}}{\Delta}\Big{(}\frac{a}{r}\Big{)}^{3}P_{1}(\xi)-\sum_{n=2}^{\infty}\frac{c_{1}A_{n}+c_{2}B_{n}}{\Delta_{n}}\Big{\{}\frac{1}{\chi^{2}}\Big{(}\frac{a}{r}\Big{)}^{n}-\Big{(}\frac{a}{r}\Big{)}^{n+2}\Big{\}}P_{n}(\xi),$$

$$\displaystyle v_{\theta}^{\left(2\right)}=\frac{5(A_{1}+B_{1})\chi^{5}}{\Delta}\Big{(}\frac{a}{r}\Big{)}^{3}V_{1}(\xi)-\sum_{n=2}^{\infty}\frac{c_{1}A_{n}+c_{2}B_{n}}{2\Delta_{n}}\Big{\{}\frac{n-2}{\chi^{2}}\Big{(}\frac{a}{r}\Big{)}^{n}-n\Big{(}\frac{a}{r}\Big{)}^{n+2}\Big{\}}V_{n}(\xi),$$

$$\mathcal{P}=-{\int }_{\hat{S}}{\mathbf{v}}^{\left(1\right)}\cdot {\mathbf{\Pi }}^{\left(1\right)}\cdot {\mathbf{n}}_{\hat{S}}dS,$$

$$\displaystyle\frac{\mathcal{P}}{4\pi{\mu}_{1}a}=\bigg{[}2\frac{2\chi^{2}+1}{\chi^{2}-1}A_{0}^{2}+2(Z_{1}+Z_{2})Z_{3}\frac{(A_{1}+B_{1})^{2}}{\Delta^{2}}$$

$$\displaystyle+\sum_{n\geq 2}^{\infty}\frac{1}{(2n+1)\Delta_{n}^{2}}\bigg{\{}2\Big{(}a_{n}N_{o}A_{n}^{2}+b_{n}N_{e}B_{n}^{2}+(a_{n}N_{e}+b_{n}N_{o})A_{n}B_{n}\Big{)}$$

$$\displaystyle+\frac{4}{n(n+1)}\Big{(}c_{n}\bar{N}_{o}A_{n}^{2}+d_{n}\bar{N}_{e}B_{n}^{2}+(c_{n}\bar{N}_{e}+d_{n}\bar{N}_{o})A_{n}B_{n}\Big{)}\bigg{\}}\bigg{]}$$

$$+C_0{A}_0,$$

$$\mathbf{u}({\mathbf{x}}_0)=\frac{1}{2\pi }{\int }_{\tilde{S}}\kappa (\mathbf{x})\mathbf{n}(\mathbf{x})\cdot \mathbf{G}({\mathbf{x}}_0,\mathbf{x})dS(\mathbf{x})-\frac{1}{4\pi }{\int }_{\hat{S}}\mathbf{q}(\mathbf{x})\cdot \mathbf{G}({\mathbf{x}}_0,\mathbf{x})dS(\mathbf{x}),$$

$$\mathbf{G}({\mathbf{x}}_0,\mathbf{x})=\frac{\mathbf{\delta }}{r}+\frac{({\mathbf{x}}_0-\mathbf{x})({\mathbf{x}}_0-\mathbf{x})}{r^3},$$

$${\int }_{\tilde{S}}\kappa (\mathbf{x})\mathbf{n}(\mathbf{x})\cdot \mathbf{G}({\mathbf{x}}_0,\mathbf{x})d{S}_x={\int }_{\tilde{S}}[\kappa (\mathbf{x})-\kappa ({\mathbf{x}}_0)]\mathbf{n}(\mathbf{x})\cdot \mathbf{G}({\mathbf{x}}_0,\mathbf{x})d{S}_x,$$

Full text

Suggested keywords

I Introduction

Droplets have recently been used as small, isolated, aqueous compartments to encapsulate, incubate and manipulate cells for biological assays He et al. (2005). Such droplet-based cell encapsulation is commonly accomplished in microfluidic devices which are able to precisely produce and manipulate microdroplets of adjustable sizesKöster et al. (2008); Chabert and Viovy (2008). Current microfluidic technology allows a high-throughput and controllable analysis to be performed on individual cells in their own discrete microenvironments.

In related work, droplets have been used to cage motile organisms such as the nematode Caenorhabditis elegans (C. elegans)Clausell-Tormos et al. (2008); Wen et al. (2015) in order to carry out developmental work. In these studies, the size of an encaged adult C. elegans is comparable to the droplet radius. Despite their mobility, the worms failed to propel their liquid cages, because they were immobilized. In the work of Ref. Clausell-Tormos et al. (2008), the droplet was tightly squeezed inside a capillary tube, forming a plug thus immobilized hydrodynamically by the lubrication film while in the work of Ref. Wen et al. (2015), the droplet was anchored mechanically by a microfluidic trap.

Motivated by these droplet-based encapsulations of motile organisms, we raise in this paper a simple question: is it possible for a microswimmer encaged in a droplet to propel its viscous cage and co-swim with it? One could envision setups of this type of interest to the drug delivery community using droplets as small self-contained units propelled and steered by their internal synthetic swimmers.

Recently, microrobots propelled by a magnetically-rotated helical appendage mimicking the flagella of bacteria such as Escherichia coli (E. coli) were fabricatedZhang et al. (2009); Tottori et al. (2012), encapsulated and operated inside a water-in-oil droplet in microfluidic chipsDing et al. (2016). In this case, the droplets were not mobile, presumably for two reasons: the swimmer was much smaller than the droplet and the droplet was large compared to the height of the micro-fluidic chips so that it was tightly squeezed and thus anchored hydrodynamicallyClausell-Tormos et al. (2008). Excitingly, the same group managed however to use their microrobots to push a droplet of comparable size from the exterior when the droplet was unbounded or loosely bounded.

In this paper, we conduct a combined theoretical and numerical study of a three-dimensional ($3$D) model microswimmer encapsulated in a droplet in free space. The size of the swimmer is of the same order as the radius of the droplet and we attempt to answer the following fundamental questions: Will the droplet co-swim with the swimmer? What is the swimming velocity of the droplet compared to that of the swimmer? How are the kinematics and energetics of the microswimmer affected by the confinement due to the presence of the droplet? How stable is the co-movement of the concentric pair of swimmer and droplet?

II Problem description

We consider, in the creeping-flow regime, the locomotion of a $3$D microswimmer encapsulated in a droplet. Due to hydrodynamic interactions, the motion of the swimmer is influenced by the presence of the droplet interface. The geometrical setup is shown in Fig. 1a. We use a spherical, axisymmetric squirmerLighthill (1952); Blake (1971) as our model swimmer. It achieves locomotion by squirming, i.e. by generating tangential and/or normal velocities on its fixed spherical surface. This is a classical model for physical actuation of microorganisms continuously deforming their bodies or beating their densely-packed cilia, and has been employed in the past to address a variety of biophysical aspects of locomotionMagar et al. (2003); Ishikawa et al. (2006); Michelin and Lauga (2010); Doostmohammadi et al. (2012); Zöttl and Stark (2012); Pak and Lauga (2014); Datt et al. (2015); Delfau et al. (2016). The shape of the droplet is maintained as spherical by maintaining a sufficiently large surface tension $\gamma$ on its interface, i.e. we assume to remain in the low-Capillary number limit. The radius of the squirmer is denoted by $a$ while that of the droplet is $b>a$, respectively, and $\chi=b/a>1$ is the size ratio. The fluid phases inside and outside the droplet are marked as phase $1$ and $2$. Both are Newtonian, with dynamic viscosities of ${\mu}^{\left(1\right)}$ and ${\mu}^{\left(2\right)}$, and $\lambda=\mu^{\left(2\right)}/\mu^{\left(1\right)}$ denotes the viscosity ratio. Both Cartesian ($x,y,z$) and spherical ($r,\theta,\phi$) coordinate systems are used, shown in Fig. 1b.

We then solve the steady Stokes equations for fluid phase $1$ and $2$,

[TABLE]

where $p^{\left(i\right)}$ is the dynamic pressure and $\bm{v}^{\left(i\right)}$ the fluid velocity in phase $(i)$, where $i=1$ or $2$. Following classical workLighthill (1952); Blake (1971), we impose normal and/or tangential squirming velocities on the surface of the swimmer $r=a$ to represent its effective swimming motion. These velocities are assumed to be time-independent and axisymmetric about its swimming direction, i.e., the $z$ axis passing through the centers of the squirmer and droplet. The squirmer drives the droplet to co-swim in the same direction, and hence the problem is fully axisymmetric about the $z$ axis. The velocity of the swimmer and droplet are denoted by $U_{S}$ and $U_{D}$ respectively.

In the laboratory frame of reference, the fluid velocity components $\bm{v}^{\left(1\right)}=(v_{r}^{\left(1\right)},v_{\theta}^{\left(1\right)})$ on the swimmer surface, $r=a$, are given by

[TABLE]

where $A_{n}$ (respectively $B_{n}$) indicates the $n$-th mode of the normal (respectively tangential) squirming velocities, $P_{n}$ are the Legendre polynomial, $\xi\equiv\cos{\theta}$, $V_{n}=-2P_{n}^{1}(\xi)/(n^{2}+n)$, and $P_{n}^{1}$ is the associated Legendre function of the first kind of order 1. In Eq. (2), $U_{S}$ is the value of the unknown swimming velocity of the swimmer, and $U_{D}$ is the unknown swimming speed of the droplet.

On the droplet interface $r=b$, the normal velocities in the droplet frame vanishes because the droplet does not deform. In addition, the tangential velocities and tangential stresses are continuous across the interface. These boundary conditions formulated in the laboratory frame are written as

[TABLE]

where $\bm{\Pi}^{\left(i\right)}=-p^{\left(i\right)}\mathbf{I}+\mu^{\left(i\right)}\left[\nabla\bm{v}^{\left(i\right)}+\left(\nabla\bm{v}^{\left(i\right)}\right)^{\mathrm{T}}\right]$ is the stress tensor for fluid $i$. Furthermore, the velocity $\bm{v}^{\left(2\right)}$ decays to zero in the far field $r\gg b$.

Finally, the total hydrodynamic forces exerted on both the swimmer and on the droplet interface are zero, which will be used to determine the values of both swimming velocities, $U_{S}$ and $U_{D}$. For an unbounded squirmer in a single-phase fluid, the velocity $U_{S}\equiv U_{0}$ is given byLighthill (1952); Blake (1971)

[TABLE]

III Analytical theory

We first solve the problem analytically. The methodology is based on Lamb’s general solution of the Stokes equations in spherical coordinatesLambs (1932); Happel and Brenner (1973). For a single-phase fluid with viscosity $\mu$, the fluid velocity field $\bm{v}$ can be expanded in spherical harmonics as

[TABLE]

where $p_{n}$ and $\phi_{n}$ are solid spherical harmonics satisfying $\nabla^{2}p_{n}=0$ and $\nabla^{2}\phi_{n}=0$, respectively. In axisymmetric flow, $p_{n}$ and $\phi_{n}$ are expressed by a series of Legendre functions as

[TABLE]

where $\tilde{p}_{n}$ and $\tilde{\phi}_{n}$ are constants independent of $r$ and $\xi$.

The radial and tangential velocity components $v_{r}$ and $v_{\theta}$ are then obtained as

[TABLE]

where

[TABLE]

Note that the solution for the flow in region 1 may contain all terms in the brackets of Eq. (III) while those in region 2 only contain the last two terms due to the boundary condition at infinity.

Applying this framework to our case, we use Eq. III for both the inner and outer fluid, solving for the unknown constants $\bar{p}_{n}^{\left(i\right)}$, $\bar{\phi}_{n}^{\left(i\right)}$, $\bar{p}^{(i)}_{-(n+1)}$ and $\bar{\phi}^{(i)}_{-(n+1)}$ ($i=1,2$) using the boundary conditions, Eqs. 2-3, together with the condition at infinity. Taking the $n=0,1$ terms in the series expansion of Eq. (III) with the use of Eqs. 2-3 leads to the system for the inner fluid

[TABLE]

Hence, the constants $\bar{p}^{\left(1\right)}_{n}$ and $\bar{\phi}^{\left(1\right)}_{n}$ ($n=-2,-1\textrm{\; and \;}1$) are obtained explicitly in terms of both $U_{S}$ and $U_{D}$. The constants in the outer fluid are then given by

[TABLE]

and the condition at infinity leads trivially to $\bar{p}_{-1}^{\left(2\right)}=0$ and $\bar{\phi}_{-1}^{\left(2\right)}=0$.

Applying the force-free condition for the swimmer, we have

[TABLE]

which leads to $\bar{p}^{\left(1\right)}_{-2}=0$. Applying the same condition for the droplet, we obtain $\bar{p}^{\left(2\right)}_{-2}=0$. Plugging the two constants into Eq. 7, we obtain the values of all underdetermined constants together with the velocity of the swimmer, $U_{S}$, and that of the droplet, $U_{D}$, as,

[TABLE]

and

[TABLE]

where

[TABLE]

Similarly to case of an unbounded squirmer (see Eq. 4), the swimming velocities $U_{S}$ and $U_{D}$ are seen to be independent of the squirming modes $A_{n}$ or $B_{n}$ for $n\geq 2$, but depend only on $A_{1}$ and $B_{1}$.

In order to complete the calculation and charactarize the flow in both fluids, we need to calculate the values of the constants $\bar{p}_{n}$, $\bar{\phi}_{n}$, $\bar{p}_{-(n+1)}$ and $\bar{\phi}_{-(n+1)}$ for $n\geq 2$ in the series expansion from Eq. (III). The velocities inside and outside the droplet in the laboratory frame are then obtained to be

[TABLE]

where the values of all undefined constants are provided in Appendix A.

We can finally calculate the power consumption of the squirmer, $\mathcal{P}$, which is equal to rate of working done by the squirmer on the fluid,

[TABLE]

where $\bm{n}_{\hat{S}}$ denotes the normal vector on $\hat{S}$ pointing towards the fluid. We obtain

[TABLE]

where $C_{0}$ is given in terms of the surface tension of the droplet, $\gamma$, as $C_{0}=\{\gamma-\mu^{\left(2\right)}A_{0}(2\chi^{2}+1)/(\chi^{2}-1)\}/(\pi\mu^{(1)}a^{2}\chi)$ based on the condition $\Pi^{\left(2\right)}_{rr}-\Pi^{\left(1\right)}_{rr}=2\gamma/b$. Again, all undefined constants are given in Appendix A.

IV Numerical simulations

In parallel with our theoretical approach, we use numerical simulations based on a $3$D boundary element method. By choosing the characteristic length, velocity, and stress as $b$, $\lambda\gamma/\{\mu^{(2)}(1+\lambda)\}$, and $\gamma/b$ respectively, the nondimensional boundary integral formulation for the matching-viscosity case ($\lambda=1$) can be obtained. The nondimensional velocity $\bm{u}(\bm{x}_{0})$ at position $\bm{x}_{0}$ everywhere in the domain is classically written as

[TABLE]

where $\tilde{S}$ and $\hat{S}$ denote the surface of the droplet and swimmer respectively, $\bm{n}$ the normal vector on $\tilde{S}$ towards the outer fluid, $\kappa=-\frac{1}{2}\nabla_{s}\cdot\bm{n}$ the mean curvature of $\tilde{S}$, and $\bm{q}$ the density of the single-layer potential on $\hat{S}$. The tensor $\bm{G}$ is the free-space Green’s function, also known as the Stokeslet or the Oseen-Burgers tensor,

[TABLE]

where $\bm{\delta}$ is identity tensor and $r=|\bm{x}_{0}-\bm{x}|$. As shown in Eq. 16, only single-layer integration is performed, which is sufficient for the rigid body motion of the swimmer and the dynamics of a matching-viscosity dropletPozrikidis (1992).

The surfaces of the swimmer and droplet are discretized using zero-order flat quadrilateral and second-order curved triangular elements respectively. For the spherical swimmer, a six-patch structured meshHigdon and Muldowney (1995); Zhu et al. (2013) consisting of $600$ (before mesh refinement) elements is constructed. The number of elements on the droplet interface is around $2500$ ($\sim 5000$ discretized points). Gauss-Legendre quadrature is applied on the quadrilateral elements to compute nonsingular integrations; on triangular elements, we compute the integrations using a symmetric Gaussian quadrature ruleDunavant (1985). When $\bm{x}_{0}$ is on the surfaces $\tilde{S}$ or $\hat{S}$, the surface integrals become singular and different desingularization strategies are chosen: on the droplet interface $\tilde{S}$, the well-known integral identity for $\bm{G}$ is exploited and hence the first integral in Eq. 16 becomes

[TABLE]

where the $\mathcal{O}(r^{-1})$ singularity of the original integrand is removed; on the squirmer surface $\tilde{S}$, each quadrilateral element is divided into four triangular sub-elements, where polar coordinates transformationPozrikidis (2002) with Gauss-Legendre quadrature is adopted to desingularize the integral. Both integrals in Eq. 16 tend to be nearly singular when the distance between the two surfaces $\tilde{S}$ and $\hat{S}$ is too small. Desingularizing measures are hence taken for them: on the droplet interface $\tilde{S}$, a high-order near-singularity subtraction is implemented by following Ref. Zinchenko and Davis (2006) and on the swimmer surface $\hat{S}$, adaptive mesh refinement is utilized. Figure 2 presents a schematic view of the adaptively-refined mesh.

A crucial numerical difficulty arising from droplet/bubble simulations based on Lagrangian interface representation is to maintain the quality of the mesh of the interface. In order to guarantee the smoothness and orthogonality of the triangle mesh over a long time evolution, we implement a so-called ‘passive mesh stabilization’ scheme Zinchenko et al. (1997); Zinchenko and Davis (2013). At each time step, the scheme searches the optimal tangential field that is added to the normal velocity to update the Lagrangian points, minimizing a global kinetic-energy-like norm that quantifies the clustering and distortion of the mesh. This scheme significantly slows down mesh degradation. Its effectiveness was proved in the previous study on a squeezed pancake droplet in a microfluidic chip based on an accelerated boundary integral implementationZhu and Gallaire (2016).

In contrast to the infinite-surface tension assumed in the theory, a large but finite surface tension is adopted in the simulations and hence the numerical droplet is not strictly spherical but slightly deformable. The strength of the typical ratio of viscous stresses to surface tension forces is measured by the capillary number, $\text{Ca}\equiv\mu^{(2)}B_{1}/\gamma$, and $\text{Ca}=0$ corresponds to the theoretical limit of infinite surface tension. We vary $\mathrm{Ca}$ numerically from $10^{-3}$ to $10^{-2}$, without detecting significant changes in the kinematics of the swimmer. We hence use $\mathrm{Ca}=10^{-3}$ throughout our study, and are able to approximate well the $\mathrm{Ca}=0$ limit from our a posteriori comparison with the theory.

V Results

V.1 Squirming with purely tangential velocities

In this section, we start by investigating the instantaneous dynamics of a droplet encapsulating a squirmer using solely tangential surface velocities, i.e. with $A_{n}=0$. If one further sets the $B_{n}$ ($n\geq 3$) modes to zero, as is classically done for the squirmer model Ishikawa et al. (2006), the swimming gait consists of only $B_{1}$ and $B_{2}$ modes: the $B_{1}$ mode determines the swimming velocity while the $B_{2}$ mode captures the leading order disturbance flow induced by the swimmer, namely a stresslet (or force dipole). We define $\beta\equiv B_{2}/B_{1}$ to measure the relative strength of the stresslet. The squirmer is said to be neutral when $\beta=0$, while it is a pusher (respectively a puller) when $\beta$ is negative (respectively positive). Varying the value of $\beta$ allows to model the majority of swimming microorganisms and synthetic microswimmers: pushers model flagellated bacteria such as E. coli Berg (2004) while biflagellated green algae such as C. reinhardtiiHernandez-Ortiz et al. (2009); Spagnolie and Lauga (2012); Goldstein (2015) are pullers. Neutral swimmers may be considered as special cases of synthetic swimmers such as Janus particles self-propelling owing to various phoretic mechanisms or some active droplets driven by Marangoni stressesYoshinaga et al. (2012); Schmitt and Stark (2013); Herminghaus et al. (2014); Maass et al. (2016); Golestanian et al. (2005); Anderson (1989); Wang et al. (2013); Colberg et al. (2014).

V.1.1 Velocity of the swimmer and droplet: theory

For a tangential squirmer, the velocities $U_{S}$ and $U_{D}$ are given analytically by

[TABLE]

where $\Delta$ is defined in Eq. III and we use the swimming velocity $U_{0}$ of an unbounded squirmer as the reference scale, $U_{0}=2B_{1}/3$. Both velocities are functions solely of the size ratio, $\chi$, and the viscosity ratio, $\lambda$.

We plot in Fig. 3a the dependence of the swimmer velocity $U_{S}$ on $\chi$ and $\lambda$. The velocity decreases monotonically with $\lambda$. When the outer and inner phase have matching viscosities ($\lambda=1$), $U_{S}$ is not affected by the presence of the droplet, and is thus equal to the unbounded velocity $U_{0}$ for all values of $\chi$. The squirmer swims faster than the unbounded one when the outer phase is less viscous than the inner ($\lambda<1$), and swims slower in the opposite limit, $\lambda>1$. When $\lambda\neq 1$, the velocity $U_{S}$ varies with the size ratio $\chi$ non-monotonically, reaching its maximum value for $\lambda<1$ when the swimmer is tightly confined, $\chi\approx 1.1\sim 1.2$, namely when the droplet is slightly larger than the swimmer. The result is similar when $\lambda>1$ and the minimum is reached. For any viscosity ratios, $U_{S}=U_{0}$ in the limit of $\chi=1$ and $\chi\rightarrow\infty$. The former corresponds to the situation when the droplet exactly encompasses the swimmer and the latter to when the droplet is much larger than the swimmer. In Fig. 3b, we further show that the velocity $U_{D}$ of the droplet decreases monotonically with $\chi$, as well as with $\lambda$. The inset of Fig. 3b presents the ratio $U_{D}/U_{S}$ of the droplet velocity over the swimmer velocity as a function of $\chi$ in a log-log form, this ratio decays as $\chi^{-3}$ for large $\chi$. It is important to note that for any values of $\chi$ or $\lambda$ the swimmer is always faster than the droplet, $U_{S}>U_{D}$. The concentric configuration is thus not a steady state if the swimmer only applies tangential forcing.

V.1.2 Comparisons between theory and

simulations

Here we consider the dynamics of a neutral swimmer ($\beta=0$), a pusher with $\beta=-5$ and a puller with $\beta=5$ encapsulated inside a same-viscosity droplet ($\lambda=1$). For simplicity we further take $B_{n}=0$ for $n\geq 3$ and $A_{n}=0$ for all $n$. Since the velocities $U_{D}$ and $U_{S}$ are independent of $\beta$, the ratio $U_{D}/U_{S}$ only depends on the value of $\chi$. This functional dependence is plotted in Fig. 4, showing an excellent agreement between the theory (green lines) and numerical data (red squares).

Next in Fig. 5a-c, we plot the flow velocity field, $\bm{v}/B_{1}$, in the laboratory frame for the pusher (a), neutral (b) and puller (c) swimmers respectively. The size ratio is $\chi=2$. Theoretical results are shown on the left panel and numerical data on the right. The numerical predictions show good agreement with the theoretical data in most of the flow domain except very close to the droplet interface where numerical errors arise from the nearly-singular integration.

For the neutral swimmer, note that the velocity field is not affected at all by the presence of the droplet. This is corroborated by the fact that neither the swimming velocity nor the power are impacted by the droplet, as implied by Eq. 10 and Eq. 15. This results from the vanishing radial velocity in the droplet frame, such that the spherical droplet interface introduces no perturbation and hence does not influence the swimming dynamics.

For the pusher in a drop, similarly to a pusher in free space, fluid is locally pushed away from the anterior ($\theta=0$) and posterior ($\theta=\pi$) parts of the swimmer and comes to the lateral directions ($\theta=\pi/2$). Due to the non-penetrating nature of the droplet interface, two counter-rotating toroidal vortices form inside. Outside the droplet the fluid is drawn towards its poles and expelled away on the equatorial plane. Interestingly the flow signature of a local pusher turns therefore into a puller in a far field. More quantitatively one can show that the velocity fields of a puller with $\beta>0$ and a pusher with $-\beta$ satisfy the relation

[TABLE]

which indicates that the mirror symmetry about the equatorial plane $\theta=\pi/2$ of the flow field of the pusher with $-\beta$ is equivalent to the reversed flow field of the puller with $\beta$.

We next investigate the spatial variation of $\bm{v}\left(z\right)/B_{1}$ along the $z$ axis in Fig. 5d, e and f for the three swimmers. Here again, numerical data (empty red circles) agree very well with the theoretical predictions (solid bule line). The velocity magnitude, $|\bm{v}|/B_{1}$, decays in the far-field from the swimmer center as $r^{-2}$ for the pusher/puller and $r^{-3}$ for the neutral swimmer. The velocity distribution $\bm{v}\left(z\right)$ over $z$ for the pusher and that for the puller are symmetric about $z=0$, as implied by Eq. V.1.2. For both swimmers, two stagnation points appear near the droplet interface $r=b$, one close to the frontal interface and other close to the rear. They can be observed in Fig. 5.

It is worth emphasizing the result that the presence of droplet reverses the direction of the far-field flow with respect to that of a pusher/puller in free space (Fig. 5a and c). This can be made more precise by an analysis of the theoretical predictions in Eq. 13. With only $B_{1}$ and $B_{2}$ modes, the leading-order contribution to the radial velocity $v_{r}^{\left(2\right)}$ in the outer phase is

[TABLE]

and that to the radial velocity of an unbounded pusher/puller is given by Ref. Blake (1971) as

[TABLE]

Their ratio is $v_{r}^{\left(2\right)}|_{\mathrm{leading}}/v_{r}|_{\mathrm{leading}}=c_{2}/\left(\Delta_{2}\chi^{2}\right)$, which is negative for any size ratio $\chi>1$ hence rationalizing the velocity inversion.

V.1.3 Power consumption

When the viscosities inside and outside the droplet are equal ($\lambda=1$) and the swimmer uses tangential surface actuations alone, the power consumption $\mathcal{P}$ based on Eq. 15 is simplified to

[TABLE]

where

[TABLE]

Restricting then our attention to the simplest squirmer with $B_{n}=0$ for $n\geq 3$, the power becomes

[TABLE]

of a similar form to that of an unbounded squirmerBlake (1971)

[TABLE]

Theoretical and numerical values of $\mathcal{P}$ show excellent agreement, as shown in Fig. 6. The power of an encapsulated squirmer, $\mathcal{P}$, always exceeds that of an unbounded one, $\mathcal{P}_{0}$. From a practical standpoint, $\mathcal{P}$ approximately doubles when the radius of the droplet is $50\%$ larger than that of the swimmer. We further observe that $\mathcal{P}$ is negatively correlated to $\chi$, and the swimmer expends more energy due to a stronger confinement. The inset log-log plot indicates that scaled excessive power $\mathcal{P}/\mathcal{P}_{0}-1$ decreases with the size ratio as $\chi^{-3}$.

V.2 Co-swimming by combining tangential and normal

squirming

We have shown in the previous sections that a swimmer employing solely tangential squirming modes, $B_{n}$, is always faster than the droplet, i.e. $U_{S}>U_{D}$. Thus, the swimmer and droplet cannot remain concentric. With the idea of using artificial swimmers encapsulated in a droplet for controllable cargo delivery, it is attempting to try and tune the squirming gait such that the swimmer and droplet co-move with a same velocity $U_{S}=U_{D}$ and maintain a concentric configuration. We find that a squirmer combining both tangential and normal velocities is able to accomplish this, as shown below.

The results in Eq. 10 and 11 imply that the swimming velocities $U_{S}$ and $U_{D}$ only depend on the first modes, $A_{1}$ and $B_{1}$. We define $\alpha\equiv A_{1}/B_{1}$ to indicate the relative strength of the modes. By comparing Eq. 10 and 11, we find that a particular value of $\alpha$, denoted by $\alpha^{\mathrm{co}}$, allows to obtain equal velocities, namely

[TABLE]

leading to a co-swimming squirmer and droplet velocity, $U_{SD}^{\mathrm{co}}$, given by

[TABLE]

For any size ratio $\chi>1$, $\alpha^{\mathrm{co}}>0$ and thus a positive $A_{1}$ mode, which contributes to the swimming velocity negatively and therefore enables the squirmer to co-swim with the droplet.

The influence of confinement $\chi$ and viscosity ratio $\lambda$ on the resulting co-swimming speed is depicted in Fig. 7 by plotting the scaled co-moving speed $U_{SD}^{\mathrm{co}}/U_{0}$, where $U_{0}=2B_{1}/3$ is the velocity of an unbounded squirmer with pure tangential modes. Even for small viscosity ratio ($\lambda=0.1$), the co-moving velocity $U_{SD}^{\mathrm{co}}$ of the pair remains below $0.7U_{0}$. Simulations have been performed to determine the values of $\alpha^{\mathrm{co}}$ and $U_{SD}^{\mathrm{co}}$ for the $\lambda=1$ case, and here again the numerical results show excellent agreement with the theory (not shown).

The relation between the mode strength $\alpha$ and the size ratio $\chi$ required to achieve concentric co-swimming is given by Eq. 27 for arbitrary viscosity ratio $\lambda$. When $\lambda$ is fixed, the particular value $\alpha^{\mathrm{co}}$ ensuring co-swimming is easily chosen as a function of $\chi$. Conversely, one may determine a particular size ratio $\chi^{\mathrm{co}}$ as a function of $\alpha$ by solving the quintic equation. In the case of $\lambda=1$, the required size ratio $\chi^{\mathrm{co}}$ is simply given by

[TABLE]

It implies that for a given swimmer with fixed modes one may select a particular size of droplet transportable by the swimmer in a co-swimming state. This encouraging result points to a practical route toward building self-propelled chemical droplets.

V.3 Stability of co-swimming state: axisymmetric configuration

While the analysis above shows that co-swimming is possible, it is not clear a priori if such configuration would be stable. In order to address the stability of swimmers, we perform numerical simulations for a swimmer-droplet pair which are initially off-center but axisymmetric. The stability problem depends on many parameters including the size ratio $\chi$, the viscosity ratio $\lambda$, the value of the mode ratio $\alpha^{\mathrm{co}}$, the stresslet strength $\beta$, and the initial offset distance $z_{\textrm{off}}$. In order to make the problem tractable, we restrict the parameter values as $\chi=0.5$, $\lambda=1$, $\alpha^{\mathrm{co}}=1.4$ and $\beta=-5,0,5$. We use $z_{\textrm{off}}=z_{\textrm{sq}}-z_{\textrm{dp}}$ to denote the offset distance in the axial direction, where $z_{\textrm{sq}}$ and $z_{\textrm{dp}}$ are the axial positions of the swimmer and droplet respectively and all simulations start with $z_{\textrm{off}}(t=0)=\pm 0.2a$.

Figure 8 (top row) displays the time evolution of $z_{\textrm{off}}$ for a swimmer which starts initially ahead (blue dot-dashed lines) or behind (red solid lines) using a tangential squirming of $\beta=-5$ (pusher, a), $\beta=0$ (neutral, b) and $\beta=5$ (puller, c). The physical characteristic time $T=b/B_{1}$ is used to scale the time $t$. For the co-moving pusher as shown in Fig. 8a, the offset $z_{\textrm{off}}(0)$ decays to zero regardless of its sign: the concentric co-moving state is recovered and remains stable. The influence of $z_{\textrm{off}}(0)$ for the co-moving neutral swimmer is shown in Fig. 8b. The concentric co-movement is seen to be stable if the swimmer is initially ahead of the droplet, but it is unstable and yields a finite-time collision between the swimmer and the droplet interface, when the swimmer is initially behind. In contrast, for the puller illustrated in Fig. 8c, the swimmer eventually touches the rear interface indicating instability when $z_{\textrm{off}}(0)<0$, while when $z_{\textrm{off}}(0)>0$, the pair reaches an eccentric co-moving state that is asymptotically stable. In the later case, the swimmer is close to the front droplet interface but separated by a thin lubrication film which acts to stabilize their co-movement via hydrodynamic interactions. The asymptotically steady thickness of the film is about $0.08a$.

The stability properties of the co-moving state seen in Fig. 8 may be interpreted physically by examining the disturbance flow field induced by the swimmer. We plot in Fig. 8 (middle row) the disturbance flow patterns corresponding to the co-moving swimming gaits which consist of normal squirming $\alpha^{\mathrm{co}}$ (dashed magenta lines) and tangential squirming $\beta$ (solid black lines). The disturbance flow of the pusher and puller is characterized by a stresslet oriented in the swimming direction, decaying as $1/r^{2}$; that of the neutral swimmer resembles a source dipole along the same direction, decaying faster as $1/r^{3}$. The analysis of Ref. Blake (1971) shows that the flow induced by the $A_{1}$ mode squirming is equivalent to that by a neutral swimmer with $B_{1}=A_{1}$. The details of this disturbance flow dictate hydrodynamic interactions between the swimmer and its environment. As can be seen in Fig. 8d, a body located in front of or behind a pusher tends to be repelled by it while it will tend to be attracted for a puller. In contrast for a neutral swimmer with $A_{1}>0$, ahead of the swimmer will be repulsive while it will tend to be attractive behind it.

We then link in Fig. 8 (bottom row) the disturbance flow of the swimmer and its relative movement with respect to the droplet, where solid/dashed circles denote the swimmer’s initial/final location (the dot-dashed circles denotes an intermediate position). As seen in Fig. 8g, for a co-moving pusher initially ahead of the droplet center, the repulsive flow in front of the swimmer, consisting of both repulsive flows from tangential squirming of $\beta=-5$ and normal squirming of $\alpha=1.4$, is stronger than its rear counterpart and brings the swimmer back to the center (stable). For the same swimmer but initially closer to the rear of the droplet as depicted in Fig. 8j, the rear flows dominate. While the flows induced by the two squirming modes are of opposite sign, the repulsive flow arising from tangential squirming is likely to overcome the attractive one of the normal squirming due to the faster-decaying and shorter-ranged disturbance flow of the latter ($1/r^{3}$ vs. $1/r^{2}$).

The behavior of the co-moving neutral swimmer can be understood along the same vein, as illustrated in Fig. 8h and k, and similarly for the puller when its is initially located behind that of the droplet (Fig. 8l). The only non-intuitive result is the asymptotically-stable eccentric location of the co-moving puller that is originally closer to the droplet front as illustrated in Fig. 8i. Initially, the gap between the swimmer and interface is relatively large, therefore the longer-ranged attractive flow from the tangential squirming will outweigh the shorter-ranged repulsive one from the normal squirming, and the swimmer will be attracted towards the interface. As the gap width decreases, the repulsive short-range flow becomes stronger, eventually dominating and preventing the swimmer from further approaching the interface. This explains, at least qualitatively, why hydrodynamic interactions lead in this situation to a stable eccentric configuration.

Additional simulations were then performed with $1/\chi$ ranging from 0.3 to 0.7 and $\beta$ ranging from $-5$ to $5$. These simulations show that the stability properties of the co-moving state is independent of the size ratio $\chi$ and depend only on $\beta$. As shown in Fig. 9, when $\beta\leq 0$, the concentric co-movement state is stable regardless of the sign of the initial offset $z_{\textrm{off}}$. When $\beta\geq 1$, the eccentric co-moving state is stable if the swimmer is initially ahead ($z_{\textrm{off}}>0$) while no stable co-moving configuration is observed otherwise ($z_{\textrm{off}}<0$).

V.4 Stability of co-swimming state: non-axisymmetric configuration

We next address the issue of stability when the initial position of the swimmer center is not aligned with the droplet along the $z$ axis. Since the system is not axisymmetric in this case, we employ numerical simulations allowing the swimmer to display rotational motion. We track the two offset distance in $x$ and $z$ directions with $x_{\textrm{off}}=x_{\textrm{sq}}-x_{\textrm{dp}}$ and $z_{\textrm{off}}=z_{\textrm{sq}}-z_{\textrm{dp}}$. When $\chi=2$ and $\alpha^{\mathrm{co}}=1.4$, we consider three types of swimmers, namely a pusher with stresslet strength $\beta=-5$, a neutral swimmer with $\beta=0$, and a puller with $\beta=5$.

We first plot in Fig. 10a the trajectories of pullers in the laboratory frame with an initial offset $(x_{\textrm{off}},z_{\textrm{off}})=(0.2a,0.2a)$. Initially the system is not axisymmetric but after a slight rotation the swimmer settles in an axisymmetric configuration. Although the rotational motion is small, it occurs early in the dynamics, in particular before the swimmer closely approaches the droplet. After that, the system becomes equivalent to the axisymmetric situation considered in Fig 8c and the swimmer reaches a stable state maintaining a thin gap with the droplet.

Next we show in Fig. 10b the trajectories of pushers with an initial offset $(x_{\textrm{off}},z_{\textrm{off}})=(0.2a,-0.2a)$. The swimmer slightly rotates but in this case does not align with the droplet axisymmetrically. Instead, due to the attractive flows in the lateral directions, the pusher approaches the droplet and eventually collides with it. Other cases with the initial offset $(x_{\textrm{off}},z_{\textrm{off}})=(0.2a,0.2a)$ or $(0.2a,0)$ exhibit similar behaviors as in Fig. 10b with no stable configurations. Also pullers with the initial offset $(x_{\textrm{off}},z_{\textrm{off}})=(0.2a,-0.2a)$ or $(0.2a,0)$ and neutral swimmers with $(x_{\textrm{off}},z_{\textrm{off}})=(0.2a,\pm 0.2a)$ or $(0.2a,0)$ do not settle a stable configuration. Additional simulations by changing the size ratio and stresslet strength leads to similar results.

VI Conclusion

In this paper, we have studied in the creeping flow regime the dynamics of a spherical squirmer encapsulated in an undeformable droplet using both theory and computations. The incompressible Stokes equations were first solved analytically, and when the swimmer and droplet are concentric, we obtained exact solutions of the swimmer and droplet velocities, the flow velocity fields and its dissipated power. Along with this analytic approach, numerical simulations based on a boundary element method were performed and the numerical results agreed well with the theoretical results.

The analytical solutions provide a useful physical picture of the instantaneous dynamics for the concentric configuration of the squirmer and droplet. For a squirmer using pure tangential surface actuations, although their movement are doomed to be transient, the theoretical results state that the swimmer is always faster than the droplet. When the normal surface velocities are incorporated on top of tangential modes, the squirmer and droplet are able to co-swim with a same velocity and thus to remain concentric.

When the swimmers are slightly displaced from the concentric position, we found that they would either return to the center (stable), deviate further and eventually touch the droplet interface (unstable), or reach an eccentric steady-state position (stable). Such final states depend on swimming gaits or relative locations of swimmers.

The ultimate goal of encaging swimmers is to help transport and deliver small chemical payloads, and thus a lot of future work lies ahead for swimmer-droplet complexes. Questions including swimming near complex boundaries or near walls, or non-axisymmetrically, will have to be tackled. Surfactants, which are commonly used in droplet-based microfluidics to prevent coalescence, could perhaps be used here to prevent collision between swimmers and interface, with interesting physical consequences. Finally, if heterogeneous fluid mixtures are to be transported in the droplet, it will important to quantify their mixing and chemical fate as they move along with the swimmer.

VII Acknowledgements

Gioele Balestra is acknowledged for his helpful suggestions on making the $3$D schematic plot. The computer time is provided by the Swiss National Supercomputing Centre (CSCS) under project ID s603 and by SNIC (Swedish National Infrastructure for Computing). A VR International Postdoc Grant from Swedish Research Council (L.Z.), an ERC starting grant ’SimCoMiCs 280117’ (F.G.), a Marie Curie CIG Grant (E.L.) and an ERC Consolidator grant (E.L.) are gratefully acknowledged.

Appendix A Constants in flow solution

The undefined constants for the fluid velocity fields in Eq. 13 and the power calculations in Eq. 15 are given in Table 1.