2.1. The exterior flow

The exterior domain consists of solvent alone and its dynamics is governed by

(2.1a,b)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{{\sigma}}_f^{e} = \eta_e\,\nabla^{2} \boldsymbol{u}^{e} - \boldsymbol{\nabla} p_f^{e} =\boldsymbol{0}, \quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}^{e}=0. \end{equation}

Here, $\boldsymbol {{\sigma }}_f^{e}$ is the stress tensor, $\eta _e$ is the dynamic viscosity of the solvent, $p_f^{e}$ is the pressure field, and $\boldsymbol {u}^{e}$ is the velocity field. Far from the vesicle, the fluid is undisturbed and there is no background flow:

(2.2a,b)\begin{equation} p_f^{e} \to 0 \quad\mbox{and}\quad \boldsymbol{u}^{e} \to {\boldsymbol{0}}\quad\mbox{as}\ r \to \infty. \end{equation}

At the exterior surface of the vesicle, we have

(2.3)\begin{equation} \boldsymbol{u}^{e}(\boldsymbol{x} \in S_e) = \boldsymbol{U} +{\boldsymbol{\varOmega}} \times \boldsymbol{x} + \boldsymbol{u}^{s}(\boldsymbol{x}), \end{equation}

where $S_e$ denotes the exterior surface of the vesicle, $\boldsymbol {U}$ (${\boldsymbol {\varOmega }}$) is the rigid-body linear (angular) velocity of the vesicle, and $\boldsymbol {u}^{s}$ is the local seepage velocity at the exterior surface. The definition of $\boldsymbol {u}^{s}$ is deferred to § 2.4. We note that (2.3) is similar to the squirmer model where the closely packed cilia tips of a microorganism are modelled as a distribution of radial and tangential velocities on the cell body, often taken to be of spherical shape (Lighthill Reference Lighthill1952; Blake Reference Blake1971).

2.2. The interior suspension

The particles and solvent in the interior of the vesicle are treated as a continuum and governed by

(2.4a,b)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{{\sigma}}^{i} = \eta_i\,\nabla^{2} \boldsymbol{u}^{i} - \boldsymbol{\nabla} P =\boldsymbol{0} \quad\mbox{and}\quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}^{i} = 0, \end{equation}

where $\boldsymbol {{\sigma }}^{i}$ is the stress tensor, $\eta _i$ is the dynamic viscosity of the suspension, and $\boldsymbol {u}^{i}$ is the velocity field. Here, the total pressure is given by

(2.5)\begin{equation} P = p_f^{i} + n {k_BT}, \end{equation}

where $p_f^{i}$ is the fluid pressure, $n$ is the number density of the ABPs, and ${k_BT}$ is the thermal energy. In our model, the only contribution to the suspension stress from the ABPs is the osmotic pressure $n{k_BT}$.

Here, the swim pressure introduced by Takatori et al. (Reference Takatori, Yan and Brady2014) does not enter the analysis directly. Regardless of activity, the particle contribution to the stress is $\boldsymbol {{\sigma }}_p = - n k_BT \boldsymbol{\mathsf{I}}$. In the high-activity limit, however, as shown in (1.3), the vesicle motion results ultimately from the swim pressure variation. Furthermore, we note that additional stress contributions, such as the active hydrodynamic stresslet of ABPs (Saintillan & Shelley Reference Saintillan and Shelley2015), can be incorporated readily into our model. Since the osmotic pressure is present regardless of activity, in this paper we focus on the osmotic pressure and neglect additional stress contributions.

At the interior wall of the vesicle, we have

(2.6)\begin{equation} \boldsymbol{u}^{i}(\boldsymbol{x}\in S_i) = \boldsymbol{U} +{\boldsymbol{\varOmega}} \times \boldsymbol{x} + \boldsymbol{u}^{s}(\boldsymbol{x}), \end{equation}

where $S_i$ is the interior surface of the vesicle.

2.3. Dynamics of ABPs

The distribution of ABPs confined inside the vesicle is described by the probability density $\varPsi (\boldsymbol {x}, \boldsymbol {q}, t)$ as a function of space $\boldsymbol {x}$, orientation $\boldsymbol {q}$ ($|\boldsymbol {q}|=1$), and time $t$. The conservation of ABPs is governed by the Smoluchowski equation. At steady state, this is given by

(2.7)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{j}_T + \boldsymbol{\nabla}_R\boldsymbol{\cdot}\boldsymbol{j}_R =0, \end{equation}

where the translational and rotational fluxes are given by, respectively,

(2.8)$$\begin{gather} \boldsymbol{j}_T= \boldsymbol{u}^{i} \varPsi + U_s(\boldsymbol{x})\,\boldsymbol{q} \varPsi -D_T\,\boldsymbol{\nabla} \varPsi, \end{gather}$$

(2.9)$$\begin{gather}\boldsymbol{j}_R= \tfrac{1}{2}{\boldsymbol{\omega}}^{i} \varPsi - D_R\,\boldsymbol{\nabla}_R \varPsi. \end{gather}$$

Here, $D_T$ is the thermal diffusivity of ABPs, ${\boldsymbol {\omega }}^{i} = \boldsymbol {\nabla }\times \boldsymbol {u}^{i}$ is the vorticity vector, $D_R$ is the rotary diffusivity, $\boldsymbol {\nabla }_R = \boldsymbol {q}\times \boldsymbol {\nabla }_q$ is the rotary gradient operator, and $U_s(\boldsymbol {x})$ is the intrinsic swim speed of ABPs. The prescribed spatial variation of $U_s$ is the key ingredient of our model, and is responsible for the generation of a concentration gradient of ABPs inside the vesicle.

The conservation of ABPs requires that

(2.10)\begin{equation} \int_{V_i} n \,{\rm d}\boldsymbol{x} = N, \end{equation}

where $n = \int \varPsi \,{\rm d}\boldsymbol {q}$ is the number density, $N$ is the total number of ABPs, and $V_i$ is the volume of the interior of the vesicle. At the interior surface of the vesicle, the flux relative to the rigid-body motion must vanish. This no-flux condition can be written as

(2.11)\begin{equation} \boldsymbol{n} \boldsymbol{\cdot}\boldsymbol{j}_T = \boldsymbol{n}\boldsymbol{\cdot}\left(\boldsymbol{U}+ {\boldsymbol{\varOmega}}\times \boldsymbol{x} \right)\varPsi,\quad \boldsymbol{x} \in S_i, \end{equation}

where $\boldsymbol {n}$ is the unit normal vector as shown in figure 1. We note that as a model of active elements inside a cell, the rotary diffusivity $D_R$ is biological rather than thermal in origin. As a result, $D_R$ is independent of $D_T$ (which is assumed to be thermal in origin). The rotary diffusivity defines a reorientation time scale $\tau _R = 1/D_R$ that characterizes the relaxation of the swimming direction. The ABPs take a step of magnitude $\ell = U_s\tau _R$, which is often called the run (or persistence) length $\ell$, before its swimming direction changes significantly. Note that one might have a reorientation time $\tau _R(\boldsymbol {x})$ that is a function of position in addition to a spatially varying swim speed, as we show below that the important quantity is the run length $\ell (\boldsymbol {x})$.

In contrast to passive Brownian particles, the self-propulsion of active particles introduces a coupling between their rotational and translational dynamics via the swimming motion. That is, even for an isolated active Brownian sphere (which is geometrically isotropic), one must track both its orientation and position. One manifestation of such a coupling is the enhanced long-time self-diffusivity beyond the thermal diffusivity $D_T$, which for an ABP with constant properties in free space is $D^{eff} = D_T + D^{swim}$, where $D^{swim} = U_s^{2}\tau _R/6$ (in three dimensions) is the swim diffusivity. In the Smoluchowski equation (2.7), the orientation dynamics is described by the rotational flux – the active particle exhibits rotary Brownian motion and is rotated by the fluid vorticity.

2.4. Transport in the membrane

We treat the fluid transport in the membrane using a macroscopic approach similar to Darcy's law; however, the porous region is modelled ultimately as a thin permeable interface. To this end, we first consider the membrane as having a network stress $\boldsymbol {{\sigma }}^{net}$ and a fluid stress $\boldsymbol {{\sigma }}_f^{m}$. The defining characteristic of the semipermeable membrane is that the fluid stress in the membrane balances the seepage velocity (Durlofsky & Brady Reference Durlofsky and Brady1987):

(2.12)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{{\sigma}}_f^{m} - \eta_e\boldsymbol{\mathsf{R}}_m \boldsymbol{\cdot}\boldsymbol{u}^{s}={\boldsymbol{0}}, \end{equation}

or $\boldsymbol {u}^{s} = \boldsymbol{\mathsf{K}} \boldsymbol {\cdot }\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {{\sigma }}_f^{m} /\eta _e$, where $\boldsymbol{\mathsf{K}} = \boldsymbol{\mathsf{R}}_m^{-1}$ is the permeability tensor, and $\boldsymbol{\mathsf{R}}_m$ is the membrane resistivity. The remaining network stress is responsible for maintaining the osmotic pressure difference across the membrane. That is, we have the force balance on the exterior and interior surfaces, respectively,

(2.13)$$\begin{gather} \boldsymbol{{\sigma}}_f^{e} \boldsymbol{\cdot} \boldsymbol{n} = \boldsymbol{{\sigma}}_f^{m} \boldsymbol{\cdot}\boldsymbol{n},\quad \boldsymbol{x} \in S_e, \end{gather}$$

(2.14)$$\begin{gather}\boldsymbol{{\sigma}}_f^{i} \boldsymbol{\cdot}\boldsymbol{n} = \boldsymbol{{\sigma}}_f^{m}\boldsymbol{\cdot}\boldsymbol{n},\quad \boldsymbol{x} \in S_i. \end{gather}$$

Note, importantly, that at the interior surface, $\boldsymbol {{\sigma }}_f^{i}$ is the interior fluid stress (it does not contain the osmotic pressure).

We model the membrane as a tangentially isotropic material with the permeability tensor

(2.15)\begin{equation} \boldsymbol{\mathsf{K}}(\boldsymbol{n}) = K_\perp \boldsymbol{n}\boldsymbol{n} + K_\parallel (\boldsymbol{\mathsf{I}} - \boldsymbol{n}\boldsymbol{n}), \end{equation}

where $K_\perp$ is the normal permeability, and $K_\parallel$ is the tangential one. For a thin membrane, the gradient in (2.12) can be approximated by a finite difference in the normal direction, which after applying the boundary conditions (2.13) and (2.14) leads to

(2.16)\begin{equation} \boldsymbol{u}^{s}(\boldsymbol{n}) = \frac{\boldsymbol{\mathsf{K}}}{\eta_e \ell_m}\boldsymbol{\cdot}\left( \boldsymbol{{\sigma}}_f^{e}|_{S_e} -\boldsymbol{{\sigma}}_f^{i}|_{S_i}\right)\boldsymbol{\cdot}\boldsymbol{n}. \end{equation}

Here, $\ell _m$ is the thickness of the membrane, and the thin membrane condition is $\ell _m \ll R$, with $R$ being the radius of the exterior surface. It is understood that in (2.16), $\boldsymbol {u}^{s}$ is a function of the local outward normal vector $\boldsymbol {n}$ (see figure 1). Equation (2.16) is a linear relation that specifies how a seepage velocity is generated in response to a jump in the fluid stress across the membrane.

In the absence of deviatoric stress, (2.16) reduces to

(2.17)\begin{equation} \boldsymbol{u}^{s} ={-} \frac{K_\perp}{\eta_e\ell_m}\left( p_f^{e}|_{S_e} -p_f^{i} |_{S_i}\right)\boldsymbol{n}, \end{equation}

which is the more familiar Darcy's law in terms of the fluid pressure difference. In general, the normal flow is driven by the fluid pressure difference as well as the shear stress.

We remark that different boundary conditions across membranes and macroscopic transport equations exist in the literature. For example, an empirical boundary condition was proposed by Beavers & Joseph (Reference Beavers and Joseph1967) and later rationalized by Saffman (Reference Saffman1971). This boundary condition was then generalized to a curved surface (Jones Reference Jones1973). Recently, using multiscale homogenization and matched asymptotic expansions between the near membrane and the far region, Zampogna & Gallaire (Reference Zampogna and Gallaire2020) developed a macroscopic condition to simulate the interaction between an incompressible fluid flow and a permeable thin membrane. For the purpose of the present work, (2.16) is sufficient.

Because the vesicle is rigid, the preservation of its volume dictates that

(2.18)\begin{equation} \int_{S_e} \boldsymbol{u}^{s} \boldsymbol{\cdot}\boldsymbol{n} \,{\rm d}S =0. \end{equation}

Henceforth, for simplicity we will assume that the membrane is not permeable in the tangential directions ($K_\parallel =0$), in which case the seepage velocity is normal to the vesicle surface.

In the above consideration, the vesicle membrane is treated as a rigid and thin porous region. To understand the material response of the vesicle, a proper treatment taking into consideration the constitutive law of the vesicle membrane is needed (Lebedev, Turitsyn & Vergeles Reference Lebedev, Turitsyn and Vergeles2007; Vlahovska & Gracia Reference Vlahovska and Gracia2007). In particular, the bending elasticity and local incompressibility give rise to a surface force density in the membrane, which is balanced by the jump in the traction from the fluid inside and the fluid outside the vesicle membrane. When such effects are included, the shape of the membrane is not known a priori and must be determined as part of the solution. If the departure from the spherical shape is small, then a perturbative approach can be adopted for both the membrane dynamics (Lebedev et al. Reference Lebedev, Turitsyn and Vergeles2007; Vlahovska & Gracia Reference Vlahovska and Gracia2007) and the fluid mechanics of a nearly spherical particle moving in a viscous fluid (Brenner Reference Brenner1964).

2.5. Dynamics of the vesicle

The rigid-body translational and rotational velocities of the vesicle are determined by the force/torque-free conditions given by

(2.19a,b)\begin{equation} \int _{S_e}\boldsymbol{{\sigma}}_f^{e} \boldsymbol{\cdot}\boldsymbol{n} \,{\rm d}S = {\boldsymbol{0}},\quad \int_{S_e}\boldsymbol{x} \times \boldsymbol{{\sigma}}_f^{e} \boldsymbol{\cdot}\boldsymbol{n} \,{\rm d} S = {\boldsymbol{0}}. \end{equation}

We can relate the rigid-body velocities $\boldsymbol {U}$ and ${\boldsymbol {\varOmega }}$ to the seepage velocity $\boldsymbol {u}^{s}$ at the exterior surface using the reciprocal theorem for Stokes flow (Masoud & Stone Reference Masoud and Stone2019). The formula for a general body shape is given in Elfring (Reference Elfring2015). For the case of a spherical particle, the rigid-body translational and rotational velocities are given by, respectively,

(2.20a,b)\begin{equation} \boldsymbol{U} ={-} \frac{1}{4 {\rm \pi}R^{2}}\int_{S_e}\boldsymbol{u}^{s} \,{\rm d} S,\quad {\boldsymbol{\varOmega}} ={-} \frac{3}{8{\rm \pi} R^{3}} \int_{S_e}\boldsymbol{n}\times \boldsymbol{u}^{s} \,{\rm d} S. \end{equation}

In the study of the rigid-body motion of micro-swimmers with prescribed kinematics (gaits), such as squirmers, the reciprocal theorem allows one to bypass the calculation of the unknown flow field, provided that one can solve the resistance/mobility problem for the swimmer shape. For the problem considered here, the seepage velocity of the vesicle is not known a priori; we need to determine the rigid-body motion, the exterior/interior flow fields and the distribution of ABPs simultaneously.

2.6. Non-dimensional equations for a spherical vesicle

For a spherical vesicle, the angular velocity vanishes (${\boldsymbol {\varOmega }} = {\boldsymbol {0}}$) and the torque balance is automatically satisfied. We define a characteristic swim speed $\tilde {U}_s$ such that

(2.21)\begin{equation} U_s(\boldsymbol{x}) = \tilde{U}_s\,\hat{U}_s(\boldsymbol{x}). \end{equation}

For a spatially homogeneous swim speed, $\hat {U}_s(\boldsymbol {x}) =1$. The average density of ABPs inside the vesicle is $\bar {n} = N/V_i$, where $V_i = 4{\rm \pi} (R-\ell _m)^{3}/3$ is the volume of the interior. We use this average density to scale the probability density such that

(2.22)\begin{equation} \varPsi = \bar{n} g, \end{equation}

where $g$ is the non-dimensional probability density. To render the governing equations non-dimensional, we scale pressures and stresses by $\bar {n}{k_sT_s}$, length by $R$, and fluid/vesicle velocities by $\bar {n}{k_sT_s}K_\perp /(\eta _e\ell _m)$. Recall that the activity is ${k_sT_s}=\zeta \tilde {U}_s^{2}\tau _R/6$.

Using the characteristic swim speed, we define the swim Péclet number as

(2.23)\begin{equation} {Pe_s} = \frac{\tilde{U}_s \tau_D}{R} = \frac{\tilde{U}_s R}{D_T}, \end{equation}

which compares the swim speed to the diffusive speed $R/\tau _D$, where $\tau _D = R^{2}/D_T$ is a diffusive time scale. Another dimensionless parameter for ABPs is defined as

(2.24)\begin{equation} \gamma = \sqrt{\frac{\tau_D}{\tau_R}} = \frac{R}{\delta}, \end{equation}

where $\delta = \sqrt {D_T \tau _R}$ is a microscopic length that quantifies the distance travelled by translational diffusion on the time scale of $\tau _R$. Alternative parameters, including $\ell /\delta$ and $\ell /R$, are often used in the literature. These parameters are direct comparisons between different length scales. We note that they are related to ${Pe_s}$ and $\gamma$ by ${Pe_s} = (\ell /\delta )^{2} (\ell /R)^{-1}$ and $\gamma = (\ell /R)^{-1}\ell /\delta$.

The non-dimensional exterior problem is given by

(2.25)$$\begin{gather} Da\, \nabla^{2} \boldsymbol{u}^{e} = \boldsymbol{\nabla} p_f^{e}, \end{gather}$$

(2.26)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}^{e}=0, \end{gather}$$

(2.27a,b)\begin{gather} \boldsymbol{u}^{e} \to {\boldsymbol{0}}\quad\mbox{and} \quad p_f^{e} \to 0\quad\mbox{as}\ r \to \infty, \end{gather}

(2.28)\begin{gather} \boldsymbol{u}^{e} = \boldsymbol{U}+ \boldsymbol{u}^{s}\quad\mbox{at}\ r = 1,\end{gather}

where

(2.29)\begin{equation} Da = \frac{K_\perp}{R \ell_m}, \end{equation}

is a Darcy number that compares the permeability of the membrane to its characteristic cross-sectional area.

In the interior, the rigid-body translation $\boldsymbol {U}$ has no effect on the fluid dynamics and we need to consider only the deviation $\boldsymbol {u}^{\prime } = \boldsymbol {u}^{i} - \boldsymbol {U}$. Thus the non-dimensional flow problem in the interior is governed by

(2.30)$$\begin{gather} \beta \,Da\, \nabla^{2} \boldsymbol{u}^{\prime} = \boldsymbol{\nabla} P, \end{gather}$$

(2.31)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}^{\prime}=0, \end{gather}$$

(2.32)$$\begin{gather}|\boldsymbol{u}^{\prime}|, P < \infty \quad \mbox{at}\ r=0, \end{gather}$$

(2.33)$$\begin{gather}\boldsymbol{u}^{\prime} = \boldsymbol{u}^{s} \quad \mbox{at}\ r = \varDelta. \end{gather}$$

Here,

(2.34)\begin{equation} \beta = \frac{\eta_i}{\eta_e} \end{equation}

is the interior-to-exterior viscosity ratio, and

(2.35)\begin{equation} \varDelta = \frac{R-\ell_m}{R} \end{equation}

is the radius ratio between the interior and the exterior surfaces of the membrane. For a thin membrane, $\ell _m/R\ll 1$, $\varDelta$ is $O(1)$. The non-dimensional total pressure is given by

(2.36)\begin{equation} P = p_f^{i} + \frac{k_BT}{k_sT_s}\,n= p_f^{i} + \frac{6\gamma^{2}}{{Pe_s}^{2}}\,n, \end{equation}

where we have used the relation $k_BT/(k_sT_s) = D_T / (\bar {U}_0^{2}\tau _R/6) = 6 \gamma ^{2}/{Pe_s}^{2}$.

The non-dimensional deviatoric stress tensors in the exterior and interior are, respectively,

(2.37a,b)\begin{equation} \boldsymbol{{\tau}}^{e} =Da \left[ \boldsymbol{\nabla} \boldsymbol{u}^{e} + \left(\boldsymbol{\nabla} \boldsymbol{u}^{e}\right)^{\rm T}\right],\quad \boldsymbol{{\tau}}^{i} = \beta \,Da\left[ \boldsymbol{\nabla} \boldsymbol{u}^{\prime} + \left(\boldsymbol{\nabla} \boldsymbol{u}^{\prime} \right)^{\rm T}\right]. \end{equation}

The seepage velocity is given by

(2.38)\begin{equation} \boldsymbol{u}^{s} = \boldsymbol{n}\boldsymbol{n}\boldsymbol{\cdot}\left( \boldsymbol{{\sigma}}_f^{e}|_{S_e} -\boldsymbol{{\sigma}}_f^{i} |_{S_i}\right)\boldsymbol{\cdot}\boldsymbol{n}, \end{equation}

where $\boldsymbol {{\sigma }}_f^{e} = -p_f^{e} \boldsymbol{\mathsf{I}} + \boldsymbol {{\tau }}^{e}$ and $\boldsymbol {{\sigma }}_f^{i} = -p_f^{i} \boldsymbol{\mathsf{I}} + \boldsymbol {{\tau }}^{i}$. The volume conservation of the vesicle is

(2.39)\begin{equation} \int_{S_e} \boldsymbol{u}^{s} \boldsymbol{\cdot} \boldsymbol{n} = 0. \end{equation}

The rigid-body translational velocity of the vesicle is then

(2.40)\begin{equation} \boldsymbol{U} ={-}\frac{1}{4{\rm \pi}}\int_{S^{2}} \boldsymbol{u}^{s} \,{\rm d}\varOmega. \end{equation}

The non-dimensional Smoluchowski equation and its fluxes, boundary condition and particle conservation are, respectively,

(2.41)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{j}_T + \boldsymbol{\nabla}_R\boldsymbol{\cdot}\boldsymbol{j}_R=0, \end{gather}$$

(2.42)$$\begin{gather}\boldsymbol{j}_T = \alpha \,Da\, \boldsymbol{u}^{\prime} g + {Pe_s}\,\hat{U}_s(\boldsymbol{x})\,\boldsymbol{q} g - \boldsymbol{\nabla} g, \end{gather}$$

(2.43)$$\begin{gather}\boldsymbol{j}_R = \tfrac{1}{2}\alpha \,Da\, {\boldsymbol{\omega}}^{\prime} g - \gamma^{2}\,\boldsymbol{\nabla}_R g, \end{gather}$$

(2.44)$$\begin{gather}\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{j}_T = 0\quad\mbox{at}\ r=\varDelta, \end{gather}$$

(2.45)$$\begin{gather}\int g \,{\rm d} \boldsymbol{q} \,{\rm d}\boldsymbol{x} = \frac{4{\rm \pi}}{3}\varDelta^{3}, \end{gather}$$

where we have introduced three non-dimensional parameters, $\alpha$, ${Pe_s}$ and $\gamma$. The first parameter is a reduced osmotic pressure given by

(2.46)\begin{equation} \alpha = \frac{\bar{n} {k_sT_s} \tau_D}{\eta_e}. \end{equation}

Physically, this is a comparison between the active driving pressure ($\bar {n}{k_sT_s}$) and a viscous resistive ‘pressure’ ($\eta _e/\tau _D$ ) on the time scale $\tau _D$.

In the equations above, variables $\{\boldsymbol {u}^{e}, p_f^{e}, \boldsymbol {x}, r, \boldsymbol {U}, \boldsymbol {u}^{s}, P, \boldsymbol {u}^{\prime }\}$ and gradient operators are non-dimensional even though the same symbols as their dimensional counterparts are used. This is to avoid inconvenience in notation, and henceforth we will work with non-dimensional quantities unless otherwise noted. In table 1, we summarize the independent non-dimensional parameters and their physical implications.

Table 1.

Independent non-dimensional parameters.

It is convenient to consider the orientational moments of the probability density function. The zeroth-order moment, or the number density, is given by

(2.47)\begin{equation} n(\boldsymbol{x}) = \int_{S^{2}}g \,{\rm d}\boldsymbol{q}, \end{equation}

where $S^{2}$ is the surface of the unit sphere in $\mathbb {R}^{3}$, which represents all possible orientations that $\boldsymbol {q}$ takes. Integrating the Smoluchowski equation over all orientations, we obtain a conservation equation for the number density:

(2.48a)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{j}_{n}=0, \end{gather}$$

(2.48b)$$\begin{gather}\boldsymbol{j}_n = \alpha \,Da\, \boldsymbol{u}^{\prime} n + {Pe_s}\,\hat{U}_s(\boldsymbol{x})\,\boldsymbol{m} - \boldsymbol{\nabla} n. \end{gather}$$

This equation is coupled to the first moment, or polar order,

(2.49)\begin{equation} \boldsymbol{m}(\boldsymbol{x}) = \int_{S^{2}} \boldsymbol{q} g \,{\rm d}\boldsymbol{q}. \end{equation}

The no-flux condition (2.11) becomes $\boldsymbol {n}\boldsymbol {\cdot }\boldsymbol {j}_n =0$ for $\boldsymbol {x}\in S_i$. Multiplying the Smoluchowski equation by $\boldsymbol {q}$ and integrating over $S^{2}$, we obtain a governing equation for the polar order:

(2.50a)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{j}_{m} - \tfrac{1}{2}\alpha \,Da\, {\boldsymbol{\omega}}^{i}\times \boldsymbol{m} + 2 \gamma^{2} \boldsymbol{m}={\boldsymbol{0}}, \end{gather}$$

(2.50b)$$\begin{gather}\boldsymbol{j}_{m} = \alpha \,Da\, \boldsymbol{u}^{\prime} \boldsymbol{m} + {Pe_s}\,\hat{U}_s(\boldsymbol{x})\left( \boldsymbol{\mathsf{Q}} + \tfrac{1}{3}n \boldsymbol{\mathsf{I}}\right) -\boldsymbol{\nabla}\boldsymbol{m}, \end{gather}$$

where

(2.51)\begin{equation} \boldsymbol{\mathsf{Q}} = \int_{S^{2}} \left(\boldsymbol{q}\boldsymbol{q}-\tfrac{1}{3}\boldsymbol{\mathsf{I}}\right)g \,{\rm d}\boldsymbol{q} \end{equation}

is the trace-free nematic order tensor, and $\boldsymbol{\mathsf{I}}$ is the identity tensor of rank 2. The no-flux condition at the interior surface for the polar order becomes $\boldsymbol {n}\boldsymbol {\cdot }\boldsymbol {j}_m = {\boldsymbol {0}}$. Different from the conservation of the total number of ABPs, the polar order is not conserved, as indicated by the presence of the sink term $2\gamma ^{2}\boldsymbol {m}$ in (2.50a) even in the absence of flow. This sink term describes the randomization, due to rotary diffusion, of any polar order.

As can be inferred from the above discussion, there is an infinite hierarchical structure to the moment equations. To truncate this infinite set of equations, a closure model such as $\boldsymbol{\mathsf{Q}} ={\boldsymbol {0}}$ is often considered in the literature (Saintillan & Shelley Reference Saintillan and Shelley2015; Yan & Brady Reference Yan and Brady2015). A closure leads to a set of closed equations that can be solved as an approximation to the Smoluchowski equation. We note that a closure approximation is often not uniformly accurate across different regimes of physical parameters or different spatial/time domains, and care must be taken when interpreting results obtained from such methods (Burkholder & Brady Reference Burkholder and Brady2020; Dulaney & Brady Reference Dulaney and Brady2020; Peng & Brady Reference Peng and Brady2020). A systematic approach to derive low-order closure models that are able to approximate the full solution of the Smoluchowski equation is still lacking.

In the context of active nematic (apolar) suspensions, the Bingham closure (Chaubal & Leal Reference Chaubal and Leal1998) has been shown to agree well with the full kinetic theory, and recently, a numerical scheme has been developed to evaluate efficiently the Bingham closure (Weady, Shelley & Stein Reference Weady, Shelley and Stein2022). With this closure, simulations with high spatial resolution are performed for active nematics. As note by Weady et al. (Reference Weady, Shelley and Stein2022), their closure is formulated for apolar suspensions, and the generalization to polar active matter remains (e.g. ABPs) to be considered. Furthermore, a comparison of the accuracy of different closure models for ABPs is largely unexplored.

The mechanism for an induced concentration gradient from a prescribed activity gradient in the absence of flow has been studied in previous works (Schnitzer Reference Schnitzer1993; Tailleur & Cates Reference Tailleur and Cates2008; Row & Brady Reference Row and Brady2020). To illustrate this mechanism and motivate later discussions, we summarize the simple one-dimensional result here. In the absence of external linear or angular velocities, such as due to flow or orienting field, the governing equation in one dimension for highly active ABPs is $\boldsymbol {\nabla }\boldsymbol {\cdot }(\hat {U}_s\boldsymbol {m})=0$, where the diffusive term is neglected. The solution in one dimension is simply $m=0$. Then (2.50a) reduces to $n \hat {U}_s = {\rm const}$. Further, Row & Brady (Reference Row and Brady2020) showed that this spatial variation of activity and concentration can drive a reverse osmotic flow, i.e. fluid flow from regions of high concentration to low. In this work, we exploit this spatial variation to propel a vesicle that is able to maintain an activity gradient in the swim speed of ABPs confined inside.

3.1. Governing equations

The behaviour of the system in this small-Péclet limit can be derived systematically by considering a weakly permeable membrane, $Da \ll 1$.

If the vesicle is non-permeable ($Da = 0$), then no external or internal flows can be generated, and the vesicle remains stationary despite the non-uniform density distribution and accumulation of the ABPs at the boundary. Due to the scaling of the dimensional velocities by the permeability, the leading-order non-dimensional velocities are $O(1)$ as $Da \to 0$. To study the motion of the vesicle in the $Da \ll 1$ limit, we pose regular expansions for all fields:

(3.1)$$\begin{gather} \boldsymbol{u}^{e} = \boldsymbol{u}_0^{e} + Da\, \boldsymbol{u}_1^{e} +\cdots, \end{gather}$$

(3.2)$$\begin{gather}p_f^{e} = p_{f,0}^{e} +Da\, p_{f,1}^{e} +\cdots, \end{gather}$$

(3.3)$$\begin{gather}\boldsymbol{u}^{\prime} = \boldsymbol{u}_0^{\prime} + Da\, \boldsymbol{u}_1^{\prime} +\cdots, \end{gather}$$

(3.4)$$\begin{gather}P = P_0 +Da\, P_1 +\cdots, \end{gather}$$

(3.5)$$\begin{gather}g= g_0 +Da\, g_1 +\cdots. \end{gather}$$

The dimensionless number density is given by $n=\int g \,{\rm d}\boldsymbol {q} = n_0 +Da\, n_1+\cdots$. Similarly, the expansions for the translational and seepage velocities are, respectively,

(3.6)$$\begin{gather} \boldsymbol{U} = \boldsymbol{U}_0 + Da\, \boldsymbol{U}_1 +\cdots, \end{gather}$$

(3.7)$$\begin{gather}\boldsymbol{u}^{s}= \boldsymbol{u}_0^{s} +Da\, \boldsymbol{u}_1^{s} +\cdots. \end{gather}$$

From (2.37a,b), we know that the leading-order deviatoric stresses are $O(Da)$, which does not contribute to the $O(1)$ seepage velocity. As a result, the seepage velocity at leading order is driven by the fluid pressure difference across the membrane,

(3.8)\begin{equation} \boldsymbol{u}_0^{s} = \left(p_f^{i} |_{S_i} - p_f^{e}|_{S_e} \right)\boldsymbol{n}. \end{equation}

Inserting these expansions into the exterior Stokes equations (2.25) and (2.26) gives, to leading order,

(3.9a,b)\begin{equation} \boldsymbol{\nabla} p_{f,0}^{e} = 0,\quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_0^{e}=0. \end{equation}

The kinematic boundary condition at the exterior surface is $\boldsymbol {u}_0^{e}(r=1) = \boldsymbol {U}_0 + \boldsymbol {u}_0^{s}$. Due to the linearity of Stokes flow, we need only solve (3.9a,b) using the seepage velocity condition ($\boldsymbol {u}_0^{e}(r=1) = \boldsymbol {u}_0^{s}$); the rigid-body translation is determined from the reciprocal theorem given by (2.40). Because $\boldsymbol {u}_0^{s}$ is in the radial direction, the exterior flow is radial and given by

(3.10a,b)\begin{equation} p_{f,0}^{e} = 0,\quad \boldsymbol{u}_0^{e} = \frac{\boldsymbol{u}_0^{s}}{r^{2}}. \end{equation}

Similarly, the leading-order equation governing the interior flow is given by

(3.11a,b)\begin{equation} \boldsymbol{\nabla} P_{0} = 0, \quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_0^{\prime}=0. \end{equation}

At the interior surface, the flow field satisfies the condition $\boldsymbol {u}^{\prime }(r=\varDelta ) = \boldsymbol {u}^{s}$. We note that the interior flow field is not analytically tractable, but it is not required in order to determine the vesicle motion. The total pressure at leading order is a constant, consisting of spatially varying fluid pressure and osmotic pressure:

(3.12)\begin{equation} p_{f,0}^{i} + 6\gamma^{2} n_0/{Pe_s}^{2}=P_0 ={\rm const}. \end{equation}

Inserting the expansions into the Smoluchowski equations (2.41)–(2.45), we obtain, at leading order,

(3.13)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot} \left({Pe_s}\,\hat{U}_s(\boldsymbol{x})\,\boldsymbol{q} g_0 - \boldsymbol{\nabla} g_0 \right)- \gamma^{2}\,\boldsymbol{\nabla}_R^{2} g_0=0, \end{gather}$$

(3.14)$$\begin{gather}\boldsymbol{n}\boldsymbol{\cdot}\left({Pe_s}\,\hat{U}_s(\boldsymbol{x})\,\boldsymbol{q} g_0 - \boldsymbol{\nabla} g_0 \right)=0\quad\mbox{at}\ r=\varDelta, \end{gather}$$

(3.15)$$\begin{gather}\int g_0 \,{\rm d} \boldsymbol{q} \,{\rm d}\boldsymbol{x} = \frac{4{\rm \pi}}{3}\varDelta^{3}. \end{gather}$$

Using (2.40), (3.8) and (3.12), we obtain

(3.16)\begin{equation} \boldsymbol{U}_0 = \frac{3\gamma^{2}}{2{\rm \pi}\,{Pe_s}^{2}} \int_{S^{2}} n_0(r=\varDelta) \boldsymbol{n} \,{\rm d} \varOmega. \end{equation}

It is more intuitive to examine the above expression in its dimensional form

(3.17)\begin{equation} \boldsymbol{U}_0 = \frac{1}{4{\rm \pi}}\,\frac{K_\perp}{\eta_e \ell_m} \int_{S^{2}} \varPi^{osmo}_0 \boldsymbol{n} \,{\rm d} \varOmega, \end{equation}

where $\varPi _0^{osmo} = n^{w} k_BT$ is the dimensional osmotic pressure of ABPs in the absence of flow.

To sum up, one needs to solve (3.13)–(3.15) to obtain the density distribution of ABPs in the absence of flow, and then use (3.16) to calculate the vesicle motion. In the remainder of § 3, the subscript ‘0’ (e.g. $g_0$, $\boldsymbol {U}_0$) will be dropped for notational convenience.

In general, one can represent the number density distribution at the spherical interior wall by the complete spherical harmonic expansion

(3.18)\begin{equation} n_0(\varDelta, \theta, \phi) = \sum_{l=0}^{\infty}\sum_{m={-}l}^{m=l}C_{l,m}\,Y_l^{m}(\theta, \phi), \end{equation}

where $Y_l^{m} \!=\! \sqrt {(2l+1)(l-m)!/[4{\rm \pi} (l+m)!]}\,P_l^{m}(\cos \theta )\exp ({\rm i} m \phi )$, and $P_l^{m}$ is the associated Legendre polynomial of degree $l$ and order $m$. Using (3.16), a direct integration shows that only the $l=1$ modes contribute to the translational velocity of the vesicle. This is similar to the tangential spherical squirmer model in which only the ‘$B_1$’ mode – the coefficient of $P_1^{1}(\cos \theta )$ – contributes to the velocity of the squirmer.

3.2. High activity

We now explore the limit of high activity, ${k_sT_s}/({k_BT}) = \tilde {D}^{swim}/D_T = \ell ^{2}/(6\delta ^{2})\gg 1$, which is often observed in active matter systems (Takatori et al. Reference Takatori, De Dier, Vermant and Brady2016). Equivalently, we define $\epsilon = 1/\gamma ^{2}$ (note that ${Pe_s} = \gamma ^{2} \ell /R$) and consider the limit $\epsilon \to 0$. Expanding the probability density function $g = g^{(0)}+\epsilon g^{(1)}+\cdots$, we obtain, at leading order,

(3.19)\begin{equation} \frac{\ell}{R}\,\boldsymbol{\nabla}\boldsymbol{\cdot} \left[ \hat{U}_s \boldsymbol{q} g^{(0)}\right] - \frac{1}{\hat{\tau}_R}\,\nabla^{2}_R g^{(0)}=0, \end{equation}

where we have included the spatial variation of $\tau _R(\boldsymbol {x})$ and defined $\tau _R = \tilde {\tau }_R\hat {\tau }_R$, similar to the case of spatially varying swim speed. Integrating over the orientation space leads to an equation for the polar order:

(3.20)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\left( \hat{U}_s \boldsymbol{m}^{(0)} \right)=0. \end{equation}

Equation (3.19) is incompatible with the no-flux boundary condition and thus is valid only in the bulk of the interior. At the interior membrane surface, the swimming flux is balanced by the diffusive flux, which implies the existence of an accumulation boundary layer of thickness $O(\epsilon )$. In this high-activity limit, the number of particles in the boundary layer is still finite, which suggests that the probability density is $O(1/\epsilon )$ as $\epsilon \to 0$. Therefore, the probability density in the boundary layer admits an expansion of the form $g(y, \theta, \phi, \boldsymbol {q}) = g^{(-1)}/\epsilon + g^{(0)} +\cdots$. Defining a stretched boundary-layer coordinate in the radial direction $y = (\varDelta -r)/\epsilon$, the Smoluchowski equation, to leading order, is

(3.21)$$\begin{gather} \frac{\ell}{R}\,\hat{U}_s\big\rvert_{S_i} \boldsymbol{q}\boldsymbol{\cdot}\boldsymbol{e}_r\,\frac{\partial g^{({-}1)}}{\partial y} + \frac{\partial^{2} g^{({-}1)}}{\partial y^{2}}=0, \end{gather}$$

(3.22)$$\begin{gather}\frac{\ell}{R}\,\hat{U}_s\big\rvert_{S_i} \boldsymbol{q}\boldsymbol{\cdot}\boldsymbol{e}_rg^{({-}1)} + \frac{\partial g^{({-}1)}}{\partial y}=0\quad\mbox{at}\ y=0, \end{gather}$$

(3.23)$$\begin{gather}g^{({-}1)} \to 0 \quad\mbox{as}\ y\to +\infty. \end{gather}$$

Here, the Taylor expansion $\hat {U}_s(r, \theta, \phi ) = \hat {U}_s\rvert _{S_i} -\epsilon y ({{\rm d}\hat {U}_s}/{{\rm d}r})\rvert _{S_i}+\cdots$ is used. The solution is obtained readily as

(3.24)\begin{equation} g^{({-}1)} =\begin{cases} A_1(\theta, \phi,\boldsymbol{q}) \exp\left(-\dfrac{\ell}{R}\,\hat{U}_s|_{S_i} \boldsymbol{q}\boldsymbol{\cdot}\boldsymbol{e}_r y\right), & \boldsymbol{q}\boldsymbol{\cdot}\boldsymbol{e}_r >0,\\ 0, & \mathrm{otherwise}. \end{cases} \end{equation}

This singular accumulation occurs only for particles with orientation pointing towards the wall ($\boldsymbol {q}\boldsymbol {\cdot }\boldsymbol {e}_r >0$) because otherwise they would swim away. In (3.24), $A_1$ is an unknown function that can be determined only from the next-order solution. The boundary-layer solution $g^{(0)}(y, \theta,\phi,\boldsymbol {q})$ in the limit $y\to \infty$ needs to be matched with the solution in the bulk as $r\to \varDelta$.

At the interior surface of the vesicle ($y=0$), the leading-order density is large and given by $\gamma ^{2} \int _{\boldsymbol {q}\boldsymbol {\cdot } \boldsymbol {e}_r >0} A_1 \,{\rm d}\boldsymbol {q}$. Just outside the boundary layer (i.e. $y\to \infty$), the density is $O(1)$ as $\gamma ^{2} \to \infty$. This boundary-layer structure allows us to relate the osmotic pressure at the interior surface of the vesicle to the swim pressure outside the boundary layer. To this end, we consider the ratio $n^{w} {k_BT}/(n^{0}{k_sT_s})$, where all quantities are dimensional. The density at the wall $n^{w}$ and the density outside the boundary layer $n^{0}$ are defined locally along the interior surface and are functions of the local surface normal vector $\boldsymbol {n}$. From the above analysis, we have

(3.25)\begin{equation} \frac{n^{w}{k_BT}}{n^{0}{k_sT_s}} = \frac{\gamma^{2} \int_{\boldsymbol{q}\boldsymbol{\cdot} \boldsymbol{e}_r >0} A_1 \,{\rm d}\boldsymbol{q} }{\int g^{(0)}(y\to \infty, \theta,\phi,\boldsymbol{q})\,{\rm d}\boldsymbol{q}}\,\frac{{k_BT}}{{k_sT_s}} = f(\ell/R, \varDelta), \end{equation}

where $\gamma ^{2}{k_BT}/({k_sT_s})=6R^{2}/\ell ^{2}$ is not a function of the thermal diffusivity $D_T$ (or $\ell /\delta$). Because in general $A_1$ is not tractable analytically, the factor $f(\ell /R, \varDelta )$ in the preceding equation cannot be obtained explicitly. Nevertheless, (3.25) reveals the important fact that at high activity,

(3.26)\begin{equation} \varPi^{osmo} = n^{w}{k_BT} = \varPi^{swim}_0 f(\ell/R, \varDelta), \end{equation}

where $\varPi ^{swim}_0 = n^{0}{k_sT_s}$. In other words, the osmotic pressure at the wall is equal to the swim pressure in the bulk of the interior just outside the boundary layer, but modified by a scale factor that is a function of $\ell /R$ and $\varDelta$. We emphasize that in (3.26), all quantities are defined locally along the interior surface of the vesicle. This is a generalization of the result of Yan & Brady (Reference Yan and Brady2015) for ABPs outside an infinite planar wall, where $n^{w}{k_BT} = n^{0}{k_sT_s}$ in the limit $\gamma ^{2} \to \infty$ because of the absence of curvature of the geometry.

Equation (3.26) allows us to obtain the dimensional speed of the vesicle:

(3.27)\begin{equation} \boldsymbol{U} = \frac{1}{4{\rm \pi}}\,\frac{K_\perp}{\eta_e\ell_m} \int_{S^{2}} \varPi_0^{swim}\,f(\ell/R, \varDelta)\,\boldsymbol{n} \,{\rm d} \varOmega. \end{equation}

We note that this relation holds for ABPs with spatially varying swim speed or reorientation time.

To understand the dependence of the motion of the vesicle on $\ell /R$, we approach the problem from a micromechanical perspective using Brownian dynamics simulations that resolve the Langevin equations of motion governing the stochastic dynamics of an ABP in its physical and orientation space. The details of the simulation method are given in Appendix A. The ABPs are treated as point particles, and their hard-particle interaction with the vesicle interior boundary is implemented using the potential-free algorithm (Heyes & Melrose Reference Heyes and Melrose1993). In this approach, the force exerted on the wall due to the collision with ABPs is obtained readily. Consider a simulation of $N$ ABPs that interact only with the boundary independently, but not among themselves. After a time step ${\rm \Delta} t$, some particles might have moved outside the interior wall. For particle $i$ that is now outside, we add a displacement ${\rm \Delta} \boldsymbol {x}_i$ to the particle such that after the move, the particle is in contact with the boundary. The total force exerted on the wall is then $\boldsymbol {F}^{w}=-\zeta \sum _{i \in \mathcal {I}}{\rm \Delta} \boldsymbol {x}_i/{\rm \Delta} t$, where $\mathcal {I}$ is the set of all particles that are outside the boundary before the hard-sphere move. As seen in (3.17), the net speed of the vesicle is proportional to the net force $\boldsymbol {F}^{w}$.

In figure 2, we show the dimensionless net force exerted on the interior vesicle surface by the ABPs, $\boldsymbol {F}^{w}/(4{\rm \pi} R_i^{2}\bar {n}{k_sT_s})$, as a function of $\ell /R_i$ for ABPs with no $D_T$ (infinitely active, $\ell /\delta =\infty$) and a spatially varying swim speed. The swim speed profile is a step function given by

(3.28)\begin{equation} \hat{U}_s = \begin{cases} 1, & x < 0,\\ 1/2, & x >0. \end{cases} \end{equation}

The net force points to the side with a larger swim speed, and only the force magnitude is shown in figure 2. As $\ell /R_i$ increases, the net force decreases. For large $\ell /R_i$, the ABPs spend most of their time pushing against and sliding along the interior vesicle surface until rotary Brownian motion reorients them towards the bulk of the interior. In this limit, the number of particles pushing against the interior surface on the side of slow speed is comparable to the side of high speed.

Figure 2.

The magnitude of the dimensionless net force on the interior vesicle surface $\boldsymbol {F}^{w}/(4{\rm \pi} R_i^{2}\bar {n}{k_sT_s})$ as a function of $\ell /R_i$ for ABPs with spatially varying swim speed. The speed profile is a step function where the swim speed in one of the hemispheres is half of that in the other. The reorientation time $\tau _R$ is a constant. The net force points towards the side with a higher swim speed.

As discussed earlier, in one dimension the relation $n U_s={\rm const}$ holds for ABPs with spatially varying properties. In the interior of a vesicle, this relation is still useful for the qualitative understanding of the distribution of ABPs and the motion of the vesicle. Taking the step function given by (3.28) as an example, $n^{0} U_s={\rm const}$ means that in the bulk of the interior the density on the right side ($x>0$) is higher than that on the left ($x<0$), $n^{0}_R > n^{0}_L$. Because $n^{w} \sim n^{0} {k_sT_s}/({k_BT}) \sim n^{0} U_s \zeta \ell /({k_BT})$ and $n^{0}U_s={\rm const}$, we have $n^{w} \sim U_s$ for ABPs with constant $\tau _R$. Therefore, the density at the interior vesicle surface on the right side is lower than that on the left ($n^{w}_R < n^{w}_L$), which is opposite to the behaviour of the bulk density. Because only the ABPs at the interior surface contribute to the net force, and they can push only against the boundary, this leads to the fact that the net force is in the negative $x$-direction (to the left). If one had only observations of the number density in the bulk, one would conclude that the vesicle moves in the direction of a lower concentration – a ‘reverse’ osmotic propulsion (cf. (1.2)).

The number density profile in the bulk and the boundary layer is sketched in figure 3 (red line) for a general swim speed profile that decreases from the left to the right. The variation of the swim speed leads to a gradient in the number density in the bulk of the interior. Two thin accumulation boundary layers are established at the left and right sides of the interior vesicle surface. Because the density at the wall on the right is smaller than that on the left, $n^{w}_R < n^{w}_L$, the dimensional version of (3.12) then leads to a larger fluid pressure on the low-density side (right), $p_{f,R}^{i} > p_{f,L}^{i}$. Since the fluid pressure in the exterior is homogeneous, the fluid is pushed out of the vesicle from the right and drawn in from the left by conservation of mass. For the vesicle–ABPs system as a whole, it effectively moves by way of jet propulsion. This kind of non-inertial jet propulsion has been proposed and studied in detail by Spagnolie & Lauga (Reference Spagnolie and Lauga2010) as an alternative mechanism for the locomotion of microswimmers. In their paper, the jetting velocity distribution of a microswimmer ($\boldsymbol {u}^{s}$) is prescribed, and then the swim speed is determined from the reciprocal theorem.

Figure 3.

Schematic of the number density profile (red) and the flow direction (blue) in the high-activity limit for a swim speed profile that decreases from the left to the right. A weak density gradient is present in the bulk of the interior due to the variation of the swim speed. Two accumulation boundary layers are established at the left and right sides of the interior wall, with the density at the wall on the left larger than that on the right. The vesicle–ABPs system as a whole moves by way of jet propulsion.

Using the approximation $n^{0} U_s={\rm const}$ and the relation $\varPi _0^{swim} = n^{0}{k_sT_s} =n^{0} U_s\zeta \ell /6$, we see that it is the variation of run length $\ell (\boldsymbol {x})$ that is responsible for the net force on the vesicle interior surface and ultimately the vesicle motion. Using (3.27), a Taylor series expansion about the centre of the vesicle leads to the scaling relation $\boldsymbol {U} \sim K_\perp R \zeta n^{0} U_s \boldsymbol {\nabla } \ell /(\eta _e\ell _m)$, where $\boldsymbol {\nabla }\ell$ is the gradient of the run length at the centre of the vesicle.

3.3. A large vesicle

When the vesicle is large, the confinement is weak, $\ell /R \ll 1$, and ABPs exhibit a thin accumulation boundary layer at the wall and a uniform distribution in the bulk of the interior to leading order. To study this large-vesicle limit of $\ell /R \ll 1$, we first write (3.13) equivalently as

(3.29)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\left[\frac{\ell}{R}\,\hat{U}_s(\boldsymbol{x})\,\boldsymbol{q} g - \left(\frac{\ell}{R}\right)^{2} \left(\frac{\ell}{\delta}\right)^{{-}2} \boldsymbol{\nabla} g \right] - \nabla_R^{2} g=0. \end{equation}

In this subsection, we use the definition $\epsilon = \ell /R$ and consider the limit as $\epsilon \to 0$. In the bulk of the interior, we have the expansion $g = g^{(0)} + \epsilon g^{(1)}+\cdots$ and the leading-order equation $\nabla _R^{2} g^{(0)}=0$. The solution in the bulk is then $g^{(0)}(\boldsymbol {x}, \boldsymbol {q}) = n^{(0)}(\boldsymbol {x})/(4{\rm \pi} )$. The boundary-layer thickness is determined by a balance between the swimming and diffusive fluxes, which leads to the leading-order equation

(3.30)$$\begin{gather} - \frac{\partial}{\partial \rho} \left(\frac{\ell}{\delta}\, \hat{U}_s|_{S^{i}}\,\boldsymbol{q}\boldsymbol{\cdot}\boldsymbol{e}_r g^{(0)} + \frac{\partial }{\partial \rho}g^{(0)}\right) - \nabla_R^{2} g^{(0)}=0, \end{gather}$$

(3.31)$$\begin{gather}\frac{\ell}{\delta}\,\hat{U}_s|_{S^{i}}\,\boldsymbol{q}\boldsymbol{\cdot}\boldsymbol{e}_r g^{(0)} + \frac{\partial }{\partial \rho}g^{(0)} =0\quad\mbox{at}\ \rho=0. \end{gather}$$

Here, we have used the stretched coordinate $\rho = (\varDelta -r)/\epsilon$. Since $\ell \ll R$, curvature of the domain has no effect at $O(1)$, and the boundary-layer equation is similar to that in a planar domain. The $O(1)$ probability density in the boundary layer does not contribute to the $O(1)$ conservation because the boundary-layer thickness is $O(\epsilon )$. This means that the total conservation is given by the density outside the boundary layer alone, $\int n^{(0)}(\boldsymbol {x})\,{\rm d}\boldsymbol {x} = 4{\rm \pi} \varDelta ^{3}/3$. In the absence of curvature terms, just like the problem of ABPs on one side of an infinite planar wall (Yan & Brady Reference Yan and Brady2015), the number density at the interior wall of the vesicle at $O(1)$ can be determined analytically; the result is given by

(3.32)\begin{equation} \frac{n^{w}}{n^{0} } = 1 + \frac{1}{6}\left( \frac{\ell}{\delta}\right)^{2} \hat{U}^{2}_0|_{S^{i}}. \end{equation}

In dimensional terms, this means that the osmotic pressure at the wall is $\varPi _0^{osmo}=n^{w}{k_BT} = n^{0} {k_BT} + n^{0} {k_sT_s} \hat {U}_s^{2}$, where $n^{0}$ is the density outside the boundary layer. To determine $n^{0}$, one needs to solve (3.30) and then match the boundary-layer solution to that in the bulk.

The dimensional translational velocity in the large-vesicle limit is written as

(3.33)\begin{equation} \boldsymbol{U} = \frac{1}{4{\rm \pi}}\,\frac{K_\perp}{\eta_e\ell_m}\int_{S^{2}} \left[n^{0} {k_BT} + n^{0} {k_sT_s}\, \hat{U}_s^{2}|_{S^{i}} \right]\boldsymbol{n} \,{\rm d} \varOmega. \end{equation}

For a large vesicle, the accumulation boundary layer has a structure similar to that obtained in the high-activity limit. Even for weakly active ABPs, this accumulation boundary layer exists so long as $\ell /R \ll 1$. As expected, (3.33) reduces to a form of (3.27) if the activity is high.

3.4. Vesicle motion due to an external orienting field

Another way to achieve motion is to apply an external orienting field, which affects the orientational dynamics but not the swim speed of the ABPs. Takatori & Brady (Reference Takatori and Brady2014) showed that net directed motion of ABPs in free space can be achieved due to the fact that the external field can orient particles to move in the same direction. Instead of having ABPs with spatially varying swim speed, we consider the same orienting field as in Takatori & Brady (Reference Takatori and Brady2014), but now with ABPs confined inside the vesicle. The only change to the orientational dynamics is that the orienting field exerts an external torque that depends on the orientation of the particle relative to the field direction; the dimensional rotary flux now becomes $\boldsymbol {j}_R = \varOmega _c \boldsymbol {q}\times \hat {\boldsymbol {H}} g - D_R\,\boldsymbol {\nabla }_R g$, where $\varOmega _c$ characterizes the rate of reorientation due to the field, and $\hat {\boldsymbol {H}}$ is the direction of the field. When an ABP is aligned with the field direction ($\boldsymbol {q} \parallel \hat {\boldsymbol {H}}$), the external torque vanishes. The Smoluchowski equation (3.13) for ABPs with constant properties in the presence of an orienting field is then

(3.34)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot} \left({Pe_s}\,\boldsymbol{q} g - \boldsymbol{\nabla} g \right)+\gamma^{2}\, \boldsymbol{\nabla}_R\boldsymbol{\cdot}\left(\chi_R \boldsymbol{q} \times\hat{\boldsymbol{H}} g - \boldsymbol{\nabla}_R g\right) =0, \end{equation}

while the no-flux boundary condition (3.14) and the total conservation (3.15) remain unchanged. Here, we have defined the Langevin parameter $\chi _R=\varOmega _c\tau _R$, which measures the strength of the orienting field compared to rotary diffusion.

In the high-activity limit, an accumulation boundary layer is established at the interior wall. The boundary-layer structure is identical to that obtained for ABPs with spatially varying swim speed. At leading order, the probability density in the bulk of the interior is governed by

(3.35)\begin{equation} \frac{\ell}{R}\,\boldsymbol{q}\boldsymbol{\cdot}\boldsymbol{\nabla} g^{(0)} +\boldsymbol{\nabla}\boldsymbol{\cdot} \left(\chi_R \boldsymbol{q} \times \hat{\boldsymbol{H}}g^{(0)} - \boldsymbol{\nabla}_R g^{(0)} \right) =0. \end{equation}

Compared to (3.19) for spatial variation, the preceding equation has a constant swim speed, and the orientational dynamics is affected by the orienting field. In the boundary layer, the leading-order equation is identical to (3.21), and the density at the wall is large.

Because (3.34) together with its no-flux boundary condition is not analytically tractable, we again make use of Brownian dynamics simulations. In figure 4(a), we show the dimensionless net force exerted on the interior wall by the ABPs, $\boldsymbol {F}^{w}/(4{\rm \pi} R_i^{2}\bar {n}{k_sT_s})$, as a function of the field strength for different values of $\ell /R$. We note that the net force is in the field direction $\hat {\boldsymbol {H}}$. In figure 4(b), the same data are plotted but with the dimensionless net force multiplied by $\ell /R_i$. This rescaling allows us to collapse all data onto a single curve. In the linear response regime, the net force is proportional to $\chi _R$. On the other hand, the net force asymptotes to a finite value in the strong-field limit. This is due to the fact that at most all $N$ particles are aligned with $\hat {\boldsymbol {H}}$ and are pushing against the vesicle; further increasing of the field strength beyond this limit has no effect.

Figure 4.

(a) The magnitude of the dimensionless net force on the interior wall $\boldsymbol {F}^{w}/(4{\rm \pi} R_i^{2}\bar {n}{k_sT_s})$ as a function of the field strength $\chi _R$ for different values of $\ell /R_i$. (b) The rescaled net force $\boldsymbol {F}^{w} \ell /(4{\rm \pi} R_i^{3}\bar {n}{k_sT_s})$ as a function of $\chi _R$ for different values of $\ell /R$. All data collapse into one curve in (b). The values of $\ell /R_i$ in both panels are the same and are thus shown only in (b). In both panels, the translational diffusion is absent, i.e. $D_T\equiv 0$. In the weak-field limit $\chi _R \ll 1$, the net force is linearly proportional to $\chi _R$, as shown by the dashed line.

In ‘wet’ active matter systems such as the vesicle problem, the fluid mechanics is ultimately responsible for the motion of the vesicle and needs to be treated properly. Nevertheless, the perspective offered by the dry active matter force balance as discussed in § 1 gives the right answer for the speed of the vesicle. In particular, consider the case in which the vesicle is driven by an orienting field. The ratio $N^{w}/N$ is a function of the field strength $\chi _R$, $N^{w}/N = f(\chi _R)$. As a result, we have the qualitative scaling relation $F^{w}\sim N \zeta U_s\, f(\chi _R)$. Noting that $\bar {n} \sim N/R_i^{3}$ and ${k_sT_s}\sim \zeta U_s^{2}\tau _R$, we have

(3.36)\begin{equation} \frac{F^{w}}{4{\rm \pi} R_i^{2} \bar{n}{k_sT_s}}\sim \frac{ N \zeta U_s}{R_i^{2} \bar{n}{k_sT_s}}\,f(\chi_R) \sim \frac{R_i}{\ell}\,f(\chi_R). \end{equation}

In the weak-field limit, $f(\chi _R)\sim \chi _R$. For large $\chi _R$, $f(\chi _R) \sim 1$ (independent of $\chi _R$). The above scaling argument also explains the collapse of the data as shown in figure 4(b). The maximum that $F^{w}$ may achieve is $N \zeta U_s$, which gives the result that $F^{w}\ell /(4{\rm \pi} R_i^{3} \bar {n}{k_sT_s} ) = 2$, which is plotted as a horizontal dashed line in figure 4(b).

We note that in figure 4, the translational diffusion is absent ($D_T\equiv 0$), which allows the system to achieve the maximum in the net force on the wall. For finite thermal diffusion, the net force is reduced and so is the speed of the vesicle.