# Activity-induced propulsion of a vesicle

Zhiwei Peng  Division of Chemistry and Chemical Engineering, California
Institute of Technology, Pasadena, California 91125, USA  0000-0002-9486-2837
Tingtao Zhou  Division of Engineering and Applied Science, California
Institute of Technology, Pasadena, California 91125, USA  0000-0002-1766-719X
John F. Brady Electronic mail<EMAIL_ADDRESS>Division of Chemistry and
Chemical Engineering, California Institute of Technology, Pasadena, California
91125, USA  Division of Engineering and Applied Science, California Institute
of Technology, Pasadena, California 91125, USA  0000-0001-5817-9128

###### Abstract

Modern biomedical applications such as targeted drug delivery require a
delivery system capable of enhanced transport beyond that of passive Brownian
diffusion. In this work an osmotic mechanism for the propulsion of a vesicle
immersed in a viscous fluid is proposed. By maintaining a steady-state solute
gradient inside the vesicle, a seepage flow of the solvent (e.g., water)
across the semipermeable membrane is generated which in turn propels the
vesicle. We develop a theoretical model for this vesicle-solute system in
which the seepage flow is described by a Darcy flow. Using the reciprocal
theorem for Stokes flow it is shown that the seepage velocity at the exterior
surface of the vesicle generates a thrust force which is balanced by the
hydrodynamic drag such that there is no net force on the vesicle. We
characterize the motility of the vesicle in relation to the concentration
distribution of the solute confined inside the vesicle. Any osmotic solute is
able to propel the vesicle so long as a concentration gradient is present. In
the present work, we propose active Brownian particles (ABPs) as a solute. To
maintain a symmetry-breaking concentration gradient, we consider ABPs with
spatially varying swim speed and ABPs with constant properties but under the
influence of an orienting field. In particular, it is shown that at high
activity the vesicle velocity is
$\bm{U}\sim[K_{\perp}/(\eta_{e}\ell_{m})]\int\Pi_{0}^{\mathrm{swim}}\bm{n}d\Omega$,
where $\Pi_{0}^{\mathrm{swim}}$ is the swim pressure just outside the thin
accumulation boundary layer on the vesicle interior surface, $\bm{n}$ is the
unit normal vector of the vesicle boundary, $K_{\perp}$ is the membrane
permeability, $\eta_{e}$ is the viscosity of the solvent, and $\ell_{m}$ is
the membrane thickness.

##  1 Introduction

Targeted drug delivery is an important goal of modern nanomedicine. Recent
advances in the design, manufacture and control of nanocarriers have enabled
the delivery of such cargoes into single cells for the purpose of imaging,
diagnostics and therapeutics [West and Halas, 2003, Gao et al., 2005, Rao et
al., 2007, Torchilin, 2012]. Commonly used pharmaceutical nanocarriers include
liposomes, micelles, nanoemulsions, polymeric nanoparticles and many others
[Torchilin, 2012]. In particular, liposomes have become an important class of
carriers for the encapsulation and transport of medical cargoes because of
several advantages including their biocompatibility with human cells, the
improved solubility of drugs and versatility for chemical targeting [Pattni et
al., 2015], among others.

A liposome is a vesicle that has an aqueous solution core encircled by a
hydrophobic membrane (lipid bilayer); hydrophilic solutes dissolved in the
core cannot readily pass through the membrane while lipophilic chemicals tend
to associate with the bilayer. As a result, a liposome can be loaded with
hydrophilic, lipophilic and/or amphiphilic cargoes in the context of drug
delivery. Recently, the Moderna vaccine developed to prevent coronavirus
disease 2019 (COVID-19) has utilized a lipid based nanovesicle to encapsulate
the mRNA vaccine that encodes the SARS-CoV-2 spike glycoprotein [Jackson et
al., 2020].

The liposome-encapsulated medical cargo is transported passively, either via
diffusion or advection due to local fluid flow, which limits its ability to
overcome biological barriers. To mitigate such limitations of passive drug
delivery, active drug delivery platforms using motile microrobots (or
microswimmers), either synthetic or biohybrid, have been proposed [Medina-
Sánchez et al., 2018, Erkoc et al., 2019, Singh et al., 2019, Bunea and
Taboryski, 2020]. By attaching nanoparticle cargoes to the surface of a motile
microswimmer, the delivery system can actively navigate, access regions that
are unreachable to passive drug delivery, and be directed to the desired site
using chemotaxis or an external magnetic field [Felfoul et al., 2016, Park et
al., 2017]. Due to self-propulsion of the microswimmer, the effective
dispersion of the attached cargo is greatly enhanced, sometimes by a few
orders of magnitude, compared to the long-time self diffusivity of the
passively-transported cargo [Singh et al., 2017].

Instead of attaching a cargo to the surface of a microswimmer, one can also
encapsulate both the cargo and the microswimmer inside the vesicle.
Encapsulated microswimmers have been studied by previous works. For example,
biological microswimmers and self-propelled Janus particles haven been
successfully encapsulated inside engineered giant unilamellar vesicles (GUVs)
[Trantidou et al., 2018, Takatori and Sahu, 2020, Vutukuri et al., 2020]. The
encapsulated microswimmer provides the vesicle with enhanced super-diffusive
motion mediated through hydrodynamic interactions between the microswimmer and
the vesicle provided that the fluid is allowed to pass through the membrane of
the vesicle [Marshall and Brady, 2021].

In the present work we consider a system that combines the benefits of the
vesicle for cargo encapsulation and the self-propulsion of microswimmers for
enhanced transport. We propose an alternate model system in which the vesicle
is propelled by an osmotic flow that is induced by an actively-maintained
concentration gradient of a solute inside the vesicle. This kind of osmotic
propulsion has been proposed as an alternate mechanism for tumor cells to
migrate under strong confinement, in which case other modes of motility such
as contractility is inhibited. Stroka et al. [2014] showed that through
physical and biochemical processes, the tumor cell establishes a spatial
gradient of solute (ions), which creates a net inflow of water at the cell
leading edge and a net outflow at the cell trailing edge. As a result, this
water permeation process enables the cell to migrate through narrow channels.
We are specifically interested in studying the motility of the vesicle as a
result of a prescribed concentration gradient of a solute that is confined
inside the vesicle. Because the solute particles are not allowed to pass
through the membrane, an osmotic flow of water is generated, which in turn
propels the vesicle immersed in water.

The main question we wish to address in this work is: What is the motility of
the vesicle system in relation to the concentration gradient of the solute?
More interestingly, does the vesicle move in the same or opposite direction of
the concentration gradient?

We show by explicit calculation that for a weakly permeable membrane the
translational velocity of a rigid spherical vesicle becomes

$\bm{U}=\frac{1}{4\pi}\frac{K_{\perp}}{\eta_{e}\ell_{m}}\int_{S^{2}}\Pi^{\mathrm{osmo}}_{0}\bm{n}d\Omega,$
(1.1)

where $\Pi_{0}^{\mathrm{osmo}}=n^{w}k_{B}T$ is the osmotic pressure of the
solute at the interior wall, $n^{w}$ is the local number density of the solute
in the absence of internal fluid flow, $k_{B}T$ is the thermal energy,
$K_{\perp}$ is the membrane permeability, $\eta_{e}$ is the viscosity of the
solvent (water) and $\ell_{m}$ is the thickness of the membrane. In equation
(1.1), $\bm{n}$ is the unit outward normal vector (see figure 1) and the
integration is over the solid angle in three dimensions (3D). In this limit,
the translational velocity of the vesicle is linearly proportional to the
driving force—the osmotic pressure. As expected, a number density at the
interior wall that breaks front-back symmetry is required in order to have a
nonzero translational velocity of the vesicle.

Equation (1.1) applies generally for any osmotic solute in the weak
permeability limit so that the interior fluid flow only slightly perturbs the
solute distribution. For example, a linear solute gradient,
$n_{0}=n_{0}(\bm{0})+\bm{x}\bm{\cdot}\bm{\nabla}n_{0}$, results in

$\bm{U}=\frac{1}{3}\frac{K_{\perp}}{\eta_{e}\ell_{m}}(R-\ell_{m})k_{B}T\bm{\nabla}n_{0},$
(1.2)

where $\bm{\nabla}n_{0}$ is a constant vector and $R$ is the exterior radius
of the vesicle. Therefore, for the simple prescribed linear-density gradient,
the vesicle translates in the same direction as the gradient in number
density.

The above discussion reveals that the vesicle is able to exhibit net motion
when an interior solute concentration gradient is given. A separate, but
important, question is: How can such a solute gradient be maintained? For a
biological cell, this is achieved by its internal physical and biochemical
processes [Stroka et al., 2014]. For a synthetic vesicle system for the
purpose of enhanced transport, alternate methods need to be implemented in
order to generate such a concentration gradient.

In this work, leveraging recent advances in the understanding of the dynamics
of active matter, we propose to use active Brownian particles (ABPs) as the
solute. In addition to normal thermal Brownian motion with translational
diffusivity $D_{T}$, ABPs self-propel with an intrinsic ‘swim’ speed $U_{s}$
in a direction $\bm{q}$. The orientation of the swimming direction $\bm{q}$
changes on a reorientation timescale $\tau_{R}$ that results from either
continuous random Brownian rotations or the often-observed discrete tumbling
events of bacteria. One important intrinsic length scale due to activity is
the run or persistence length $\ell=U_{s}\tau_{R}$. Previous works have shown
that a spatial variation in the swim speed leads to a spatial variation in the
concentration (or number density) of active particles [Schnitzer, 1993,
Tailleur and Cates, 2008, Row and Brady, 2020]. By tuning the swim speed
distribution of ABPs confined inside the vesicle, a spherically asymmetric
density distribution can emerge and lead to net motion of the vesicle.

For active particles with slow spatial variation in swim speed in 1D,
Schnitzer [1993] and later Tailleur and Cates [2008] showed that the local
number density $n$ is inversely proportional to the local swim speed $U_{s}$,
i.e., $nU_{s}=const$. This simple prediction has been validated experimentally
using bacteria that swim with an intensity-dependent speed when illuminated by
a spatial light pattern [Arlt et al., 2019]. Row and Brady [2020] generalized
this result and showed that the spatial variation in activity (e.g., swim
speed) can be utilized as a pump mechanism in which fluid flows from regions
of high concentration of particles to low. Employing this spatial variation,
we show that encapsulated ABPs with spatially varying activity can be used to
propel the vesicle.

In equations (1.1) and (1.2), the vesicle velocity appears to be linearly
proportional to $k_{B}T$. However, this does not imply that the driving force
is necessarily thermal in origin (in thermodynamic equilibrium no density
gradient is present). In the case of ABPs as solute, the active (non-
equilibrium) dynamics provides such a density gradient. Analogous to the
Stokes-Einstein-Sutherland relation $k_{B}T=\zeta D_{T}$, where $\zeta$ is the
Stokes drag coefficient, an active energy scale
$k_{s}T_{s}=\zeta\tilde{D}^{\mathrm{swim}}$ can be defined for active matter
systems [Takatori et al., 2014], where
$\tilde{D}^{\mathrm{swim}}=\tilde{U}_{s}^{2}\tau_{R}/6$ is the swim
diffusivity. We note that for ABPs with spatially varying swim speed a
characteristic swim speed $\tilde{U}_{s}$ is used in the definition of the
swim diffusivity; the local active energy $k_{s}T_{s}(\bm{x})$ can also be
defined by using the local swim speed $U_{s}(\bm{x})$ and/or local
reorientation time $\tau_{R}(\bm{x})$. An important parameter that quantifies
the activity of ABPs is the ratio
$k_{s}T_{s}/k_{B}T=\tilde{D}^{\mathrm{swim}}/D_{T}$. For many active matter
systems this ratio is very large, often exceeding $10^{3}$ [Takatori et al.,
2016]. In this high activity limit, the ABPs exhibit a thin accumulation
boundary layer at the interior surface of the vesicle. As we shall show in
section 3.2, the local density at the interior wall of the vesicle can be
related to the density just outside the boundary layer via the equation
$n^{w}k_{B}T=n^{0}k_{s}T_{s}(\bm{x})f=\Pi_{0}^{\mathrm{swim}}(\bm{x})f$, where
$\Pi_{0}^{\mathrm{swim}}$ is the swim pressure just outside the boundary layer
and $f$ is a factor that depends on the ratio of the run length to the size of
the vesicle. [ This factor is unity for the case of ABPs on one side of an
infinite planar wall [Yan and Brady, 2015].] For highly active ($k_{s}T_{s}\gg
k_{B}T$) ABPs, equation (1.1) becomes

$\bm{U}=\frac{1}{4\pi}\frac{K_{\perp}}{\eta_{e}\ell_{m}}\int_{S^{2}}n^{0}k_{s}T_{s}(\bm{x})f\bm{n}d\Omega=\frac{1}{4\pi}\frac{K_{\perp}}{\eta_{e}\ell_{m}}\int_{S^{2}}\Pi_{0}^{\mathrm{swim}}\bm{n}fd\Omega,$
(1.3)

showing that the velocity of the vesicle is proportional to the swim pressure.
More precisely, it is the variation of the swim pressure [due to the variation
in swim speed or run length $\ell(\bm{x}$)] that gives rise to net motion.

Instead of using ABPs with spatially varying swim speed or run length, one can
also consider using an external field that orients constant-property ABPs
towards a certain direction. External fields such as chemical gradients or
magnetic fields can affect the swimming behavior of microorganisms to
facilitate their movement towards a favorable region. In the laboratory, an
externally applied magnetic field has been used to guide nanocarriers for the
purpose of targeted drug delivery [Felfoul et al., 2016, Pattni et al., 2015].
In the presence of an external orienting field, even for ABPs with constant
properties, the front-back symmetry is broken, and net motion of the vesicle
is generated. The balance of the strength of the orienting field and the
random reorientation due to rotary diffusion is characterized by the Langevin
parameter, $\chi_{R}=\Omega_{c}\tau_{R}$, where $\Omega_{c}$ is the strength
of the angular velocity induced by the field [Takatori and Brady, 2014].
Noting that the force exerted by the active particles on the wall
$\bm{F}^{w}=k_{B}T\int n^{w}\bm{n}dS$ [Yan and Brady, 2015], we rewrite
equation (1.1) as $\bm{U}=K_{\perp}\bm{F}^{W}/(4\pi R^{2}\eta_{e}\ell_{m})$.
In other words, we need to know the net force the active particles exert on
the wall to determine the net vesicle motion. The force on the wall scales as
$N^{w}\zeta U_{s}$, where $N^{w}$ is the total number of particles at the wall
and each particle pushes against the wall with at most its swim force $\zeta
U_{s}$. The balance of this force due to the ABPs with the drag force of the
porous vesicle moving through an external viscous fluid gives the net motion.
Of particular interest is the strong-field limit, where the number of
particles on the wall is on the same order as the total number of particles,
$N^{w}/N=O(1)$, and the net speed of the vesicle is the largest, $U\sim
K_{\perp}N\zeta U_{s}/(R^{2}\eta_{e}\ell_{m})$.

This last example where we argued that the vesicle motion can be deduced from
the net swim force of the ABPs balancing the drag of the vesicle also applies
to the so-called ‘dry’ active matter [Marchetti et al., 2013]. Dry active
matter describes bacteria (or other organisms) that crawl (or even walk) on a
surface of a medium of resistivity $\zeta$. Active particles confined to a
‘container’ that is able to slide along the surface in response to a lateral
force will be able to push the container via their ‘swim’ force if there is an
asymmetric distribution of ABPs. The net swim force would scale as $N^{w}\zeta
U_{s}$, and the container would translate with the speed $U_{c}\sim N^{w}\zeta
U_{s}/\zeta_{c}$, where $\zeta_{c}$ is the resistivity for sliding the
container along the surface. For dry active matter there is no fluid and thus
one does not have the notion of a semipermeable membrane nor a seepage
velocity driven by an osmotic pressure difference. Nevertheless, the mechanics
are the same: like the seepage velocity, the substrate surface must move
across the container boundary as it slides along the surface, and the ABPs
achieve their propulsive ‘crawling’ force by pushing off the substrate just
like swimmers push off the fluid. Thus, at least at high activity, the results
derived here apply equally well to dry active matter with an appropriate
change in notation.

In the case of a spherical vesicle, its net motion is induced by an asymmetric
number density distribution on the vesicle interior surface. An alternate
route for the generation of net motion is to use a vesicle with an asymmetric
shape. Because the accumulation of ABPs at the interior surface depends on the
local curvature of the boundary, a vesicle that has a front-back asymmetry in
its shape is able to exhibit net motion. Indeed, the exterior version of the
problem where a passive object is immersed in a bath of active particles has
been studied. It has been shown in experiments and simulations that for an
object with shape asymmetry, net motion can be achieved [Sokolov et al., 2010,
Kaiser et al., 2014, Yan and Brady, 2018].

To obtain the results for the vesicle motility, in section 2 we describe the
model and derive a theoretical formulation that governs the dynamics of the
vesicle, the interior solute suspension and the exterior fluid flow. A Darcy-
like constitutive law that models the response of the fluid seepage velocity
in relation to the fluid stress differences across the membrane is used. This
formulation is at the continuum level, where the vesicle is large compared to
the size of the ABPs so that the interior (fluid and ABPs) is treated as a
suspension; the suspension stress includes the fluid stress and the osmotic
pressure of the ABPs. The exterior flow field satisfies the boundary condition
that the fluid velocity at the exterior surface of the vesicle consists of the
rigid body motion and a seepage velocity. Because the vesicle is force- and
torque-free, we can relate the rigid body motion to the seepage velocity
distribution at the exterior surface using the reciprocal theorem. This
approach is similar to treatments of the swimming of microorganisms using the
squirmer model [Stone and Samuel, 1996] where the boundary velocity at the
surface of the swimmer is decomposed into rigid-body motion and the slip
velocity distribution.

In situations relevant for the vesicle model as we consider here, the interior
fluid flow is often weak compared to the active self-propulsion. In section 3,
by neglecting the interior fluid flow we show that the total (fluid and
osmotic) pressure inside the vesicle is constant and the leading-order
translational velocity of the vesicle is driven by the difference in the fluid
pressure across the membrane. As a result, one only needs to compute the
distribution of ABPs in the absence of flow and the resulting number density
distribution at the interior wall is used to obtain the translational
velocity. The effect of an external orienting field on the dynamics of
confined ABPs and the motion of the vesicle is considered in section 3.4. The
behavior of ABPs with slow spatial variation in their swim speed where fluid
motion is explicitly considered is discussed in section 4. Finally, we
conclude in section 5 with a discussion of the limitations and extensions of
this vesicle-ABPs propulsion system.

##  2 Problem formulation

Consider a rigid vesicle or cell consisting of a thin membrane and a solution
core immersed in an otherwise quiescent viscous fluid (see figure 1). The
interior of the vesicle is a suspension of potentially active elements, which
we model as active Brownian particles. The boundary or membrane of the vesicle
is permeable to the solvent (i.e. water) but not to the solute (ABPs). In
other words, the membrane is an osmotic membrane and serves as a confining
boundary for the ABPs. Relative to the vesicle, the fluid domain is
partitioned into interior, exterior and the thin porous (in the membrane)
regions. The solvent in all regions is identical.

The ABPs encapsulated inside the vesicle swim with a prescribed spatially
varying swim speed, which is the driving mechanism for a spatially varying
number density.

At small scales relevant to the vesicle-ABP system proposed here, the inertia
of the fluid, the ABPs and the vesicle are negligible. In particular, for
motile bacteria such as E. coli, which has a characteristic size of $\sim$1µm
and a swim speed of $\sim$30µm/s, the Reynolds number in water is $3\times
10^{-5}$. The resulting speed of the vesicle and the Reynolds number based on
the size of the vesicle and its speed are also small. In this low Reynolds
number limit, the dynamics of the fluid is governed by the Stokes equations
and there is no external force/torque on the vesicle.

Figure 1: Left: A rigid spherical vesicle with a semipermeable membrane
immersed in an otherwise quiescent viscous fluid. Active Brownian particles
are confined inside the vesicle. Right: Schematic of the semipermeable
membrane with a permeability tensor $\bm{K}$ and thickness $\ell_{m}$. The
seepage velocity in the membrane is $\bm{u}^{s}$, which in general depends on
the local position vector.

###  2.1 The exterior flow

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

$\bm{\nabla}\bm{\cdot}\bm{\sigma}_{f}^{e}=\eta_{e}\nabla^{2}\bm{u}^{e}-\bm{\nabla}p_{f}^{e}=\bm{0},\quad\bm{\nabla}\bm{\cdot}\bm{u}^{e}=0.$
(2.1)

Here, $\bm{\sigma}_{f}^{e}$ is the stress tensor, $\eta_{e}$ is the dynamic
viscosity of the solvent, $p_{f}^{e}$ is the pressure field and $\bm{u}^{e}$
is the velocity field. Far from the vesicle, the fluid is undisturbed and
there is no background flow:

$p_{f}^{e}\to 0\quad\mbox{and}\quad\bm{u}^{e}\to\bm{0}\quad\mbox{as}\quad
r\to\infty.$ (2.2)

At the exterior surface of the vesicle, we have

$\bm{u}^{e}\left(\bm{x}\in
S_{e}\right)=\bm{U}+\bm{\varOmega}\times\bm{x}+\bm{u}^{s}(\bm{x}),$ (2.3)

where $S_{e}$ denotes the exterior surface of the vesicle, $\bm{U}$
($\bm{\varOmega}$) is the rigid-body linear (angular) velocity of the vesicle
and $\bm{u}^{s}$ is the local seepage velocity at the exterior surface. The
definition of $\bm{u}^{s}$ is deferred to Section 2.4. We note that equation
(2.3) is similar to the squirmer model where the closely packed cilia tips of
a microorganism are modeled as a distribution of radial and tangential
velocities on the cell body, often taken to be of spherical shape [Lighthill,
1952, Blake, 1971].

###  2.2 The interior suspension

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

$\bm{\nabla}\bm{\cdot}\bm{\sigma}^{i}=\eta_{i}\nabla^{2}\bm{u}^{i}-\bm{\nabla}P=\bm{0}\quad\mbox{and}\quad\bm{\nabla}\bm{\cdot}\bm{u}^{i}=0,$
(2.4)

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

$P=p_{f}^{i}+nk_{B}T,$ (2.5)

where $p_{f}^{i}$ is the fluid pressure, $n$ the number density of the ABPs
and $k_{B}T$ is the thermal energy. In our model, the only contribution to the
suspension stress from the ABPs is the osmotic pressure $nk_{B}T$.

Here, the swim pressure introduced by Takatori et al. [2014] does not directly
enter the analysis. Regardless of activity, the particle contribution to the
stress is $\bm{\sigma}_{p}=-nk_{B}T\bm{I}$. In the high activity limit,
however, as shown in equation (1.3), the vesicle motion ultimately results
from the swim pressure variation. Furthermore, we note that additional stress
contributions such as the active hydrodynamic stresslet of ABPs [Saintillan
and Shelley, 2015] can be readily incorporated 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

$\bm{u}^{i}\left(\bm{x}\in
S_{i}\right)=\bm{U}+\bm{\varOmega}\times\bm{x}+\bm{u}^{s}(\bm{x}),$ (2.6)

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(\bm{x},\bm{q},t)$ as a function of space
$\bm{x}$, orientation $\bm{q}$ ($|\bm{q}|=1$) and time $t$. The conservation
of ABPs is governed by the Smoluchowski equation. At steady state, this is
given by

$\bm{\nabla}\bm{\cdot}\bm{j}_{T}+\bm{\nabla}_{R}\bm{\cdot}\bm{j}_{R}=0,$ (2.7)

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

$\displaystyle\bm{j}_{T}$ $\displaystyle=$
$\displaystyle\bm{u}^{i}\varPsi+U_{s}(\bm{x})\bm{q}\varPsi-
D_{T}\bm{\nabla}\varPsi,$ (2.8) $\displaystyle\bm{j}_{R}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\bm{\omega}^{i}\varPsi-D_{R}\bm{\nabla}_{R}\varPsi.$
(2.9)

Here, $D_{T}$ is the thermal diffusivity of ABPs,
$\bm{\omega}^{i}=\bm{\nabla}\times\bm{u}^{i}$ is the vorticity vector, $D_{R}$
is the rotary diffusivity, $\bm{\nabla}_{R}=\bm{q}\times\bm{\nabla}_{q}$ is
the rotary gradient operator and $U_{s}(\bm{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

$\int_{V_{i}}nd\bm{x}=N,$ (2.10)

where $n=\int\Psi d\bm{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

$\bm{n}\bm{\cdot}\bm{j}_{T}=\bm{n}\bm{\cdot}\left(\bm{U}+\bm{\varOmega}\times\bm{x}\right)\varPsi,\quad\bm{x}\in
S_{i},$ (2.11)

where $\bm{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 timescale $\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}(\bm{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(\bm{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^{\text{eff}}=D_{T}+D^{\text{swim}}$, where
$D^{\text{swim}}=U_{s}^{2}\tau_{R}/6$ (in 3D) 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 ultimately modelled as a
thin permeable interface. To this end, we first consider the membrane as
having a network stress $\bm{\sigma}^{\mathrm{net}}$ and a fluid stress
$\bm{\sigma}_{f}^{m}$. The defining characteristic of the semi-permeable
membrane is that the fluid stress in the membrane balances the seepage
velocity [Durlofsky and Brady, 1987]:

$\bm{\nabla}\bm{\cdot}\bm{\sigma}_{f}^{m}-\eta_{e}\bm{R}_{m}\bm{\cdot}\bm{u}^{s}=\bm{0},$
(2.12)

or
$\bm{u}^{s}=\bm{K}\bm{\cdot}\bm{\nabla}\bm{\cdot}\bm{\sigma}_{f}^{m}/\eta_{e}$,
where $\bm{K}=\bm{R}_{m}^{-1}$ is the permeability tensor and $\bm{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,

$\displaystyle\bm{\sigma}_{f}^{e}\bm{\cdot}\bm{n}$ $\displaystyle=$
$\displaystyle\bm{\sigma}_{f}^{m}\bm{\cdot}\bm{n},\quad\bm{x}\in S_{e}$ (2.13)
$\displaystyle\bm{\sigma}_{f}^{i}\bm{\cdot}\bm{n}$ $\displaystyle=$
$\displaystyle\bm{\sigma}_{f}^{m}\bm{\cdot}\bm{n},\quad\bm{x}\in S_{i}.$
(2.14)

Note, importantly, that at the interior surface, $\bm{\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

$\bm{K}(\bm{n})=K_{\perp}\bm{n}\bm{n}+K_{\parallel}(\bm{I}-\bm{n}\bm{n}),$
(2.15)

where $K_{\perp}$ is the normal permeability and $K_{\parallel}$ is the
tangential one. For a thin membrane the gradient in equation (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

$\displaystyle\bm{u}^{s}(\bm{n})=\frac{\bm{K}}{\eta_{e}\ell_{m}}\bm{\cdot}\left(\bm{\sigma}_{f}^{e}\big{\rvert}_{S_{e}}-\bm{\sigma}_{f}^{i}\big{\rvert}_{S_{i}}\right)\bm{\cdot}\bm{n}.$
(2.16)

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 equation (2.16) $\bm{u}^{s}$ is a function
of the local outward normal vector $\bm{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, equation (2.16) reduces to

$\bm{u}^{s}=-\frac{K_{\perp}}{\eta_{e}\ell_{m}}\left(p_{f}^{e}\rvert_{S_{e}}-p_{f}^{i}\big{\rvert}_{S_{i}}\right)\bm{n},$
(2.17)

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 and Joseph [1967] and later
rationalized by Saffman [1971]. This boundary condition was then generalized
to a curved surface [Jones, 1973]. Recently, using multiscale homogenization
and matched asymptotic expansions between the near membrane and the far
region, Zampogna and Gallaire [2020] 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, equation (2.16) is
sufficient.

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

$\displaystyle\int_{S_{e}}\bm{u}^{s}\bm{\cdot}\bm{n}dS=0.$ (2.18)

Henceforth, for simplicity we shall 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 et al., 2007, Vlahovska and Gracia, 2007]. 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, a perturbative approach can be adopted for both
the membrane dynamics [Lebedev et al., 2007, Vlahovska and Gracia, 2007] and
the fluid mechanics of a nearly spherical particle moving in a viscous fluid
[Brenner, 1964].

###  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

$\displaystyle\int_{S_{e}}\bm{\sigma}_{f}^{e}\bm{\cdot}\bm{n}dS=\bm{0},\quad\int_{S_{e}}\bm{x}\times\bm{\sigma}_{f}^{e}\bm{\cdot}\bm{n}dS=\bm{0}.$
(2.19)

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

$\displaystyle\bm{U}=-\frac{1}{4\pi
R^{2}}\int_{S_{e}}\bm{u}^{s}dS,\quad\bm{\varOmega}=-\frac{3}{8\pi
R^{3}}\int_{S_{e}}\bm{n}\times\bm{u}^{s}dS.$ (2.20)

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 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
($\bm{\varOmega}=\bm{0}$) and the torque balance is automatically satisfied.
We define a characteristic swim speed $\tilde{U}_{s}$ such that

$U_{s}(\bm{x})=\tilde{U}_{s}\hat{U}_{s}(\bm{x}).$ (2.21)

For a spatially homogeneous swim speed, $\hat{U}_{s}(\bm{x})=1$. The average
density of ABPs inside the vesicle is $\overline{n}=N/V_{i}$, where
$V_{i}=4\pi(R-\ell_{m})^{3}/3$ is the volume of the interior. We use this
average density to scale the probability density such that

$\varPsi=\overline{n}g,$ (2.22)

where $g$ is the non-dimensional probability density. To render the governing
equations non-dimensional, we scale pressures and stresses by
$\overline{n}k_{s}T_{s}$, length by $R$ and fluid/vesicle velocities by
$\overline{n}k_{s}T_{s}K_{\perp}/(\eta_{e}\ell_{m})$. Recall that the activity
$k_{s}T_{s}=\zeta\tilde{U}_{s}^{2}\tau_{R}/6$.

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

$Pe_{s}=\frac{\tilde{U}_{s}\tau_{D}}{R}=\frac{\tilde{U}_{s}R}{D_{T}}$ (2.23)

that compares the swim speed to the diffusive speed $R/\tau_{D}$, where
$\tau_{D}=R^{2}/D_{T}$ is a diffusive timescale. Another dimensionless
parameter for ABPs is defined as

$\gamma=\sqrt{\frac{\tau_{D}}{\tau_{R}}}=\frac{R}{\delta},$ (2.24)

where $\delta=\sqrt{D_{T}\tau_{R}}$ is a microscopic length that quantifies
the distance traveled by translational diffusion on the timescale of
$\tau_{R}$. Alternate 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

$\displaystyle Da\nabla^{2}\bm{u}^{e}$ $\displaystyle=$
$\displaystyle\bm{\nabla}p_{f}^{e},$ (2.25)
$\displaystyle\bm{\nabla}\bm{\cdot}\bm{u}^{e}$ $\displaystyle=$ $\displaystyle
0,$ (2.26) $\displaystyle\bm{u}^{e}$ $\displaystyle\to$
$\displaystyle\bm{0}\quad\mbox{and}\quad p_{f}^{e}\to 0\quad\mbox{as}\quad
r\to\infty,$ (2.27) $\displaystyle\bm{u}^{e}$ $\displaystyle=$
$\displaystyle\bm{U}+\bm{u}^{s}\quad\mbox{at}\quad r=1.$ (2.28)

where

$Da=\frac{K_{\perp}}{R\ell_{m}},$ (2.29)

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

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

$\displaystyle\beta Da\nabla^{2}\bm{u}^{\prime}$ $\displaystyle=$
$\displaystyle\bm{\nabla}P,$ (2.30)
$\displaystyle\bm{\nabla}\bm{\cdot}\bm{u}^{\prime}$ $\displaystyle=$
$\displaystyle 0,$ (2.31) $\displaystyle|\bm{u}^{\prime}|,P$ $\displaystyle<$
$\displaystyle\infty\quad\mbox{at}\quad r=0,$ (2.32)
$\displaystyle\bm{u}^{\prime}$ $\displaystyle=$
$\displaystyle\bm{u}^{s}\quad\mbox{at}\quad r=\Delta.$ (2.33)

Here,

$\beta=\frac{\eta_{i}}{\eta_{e}}$ (2.34)

is the interior-to-exterior viscosity ratio and

$\Delta=\frac{R-\ell_{m}}{R}$ (2.35)

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

$P=p_{f}^{i}+\frac{k_{B}T}{k_{s}T_{s}}n=p_{f}^{i}+\frac{6\gamma^{2}}{Pe_{s}^{2}}n,$
(2.36)

where we have used the relation
$k_{B}T/k_{s}T_{s}=D_{T}/(\overline{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,

$\displaystyle\bm{\tau}^{e}=Da\left[\bm{\nabla}\bm{u}^{e}+\left(\bm{\nabla}\bm{u}^{e}\right)^{T}\right],\quad\bm{\tau}^{i}=\beta
Da\left[\bm{\nabla}\bm{u}^{\prime}+\left(\bm{\nabla}\bm{u}^{\prime}\right)^{T}\right].$
(2.37)

The seepage velocity is given by

$\bm{u}^{s}=\bm{n}\bm{n}\bm{\cdot}\left(\bm{\sigma}_{f}^{e}\big{\rvert}_{S_{e}}-\bm{\sigma}_{f}^{i}\big{\rvert}_{S_{i}}\right)\bm{\cdot}\bm{n},$
(2.38)

where $\bm{\sigma}_{f}^{e}=-p_{f}^{e}\bm{I}+\bm{\tau}^{e}$ and
$\bm{\sigma}_{f}^{i}=-p_{f}^{i}\bm{I}+\bm{\tau}^{i}$. The volume conservation
of the vesicle is

$\int_{S_{e}}\bm{u}^{s}\bm{\cdot}\bm{n}=0.$ (2.39)

The rigid-body translational velocity of the vesicle is then

$\bm{U}=-\frac{1}{4\pi}\int_{S^{2}}\bm{u}^{s}d\Omega.$ (2.40)

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

$\displaystyle\bm{\nabla}\bm{\cdot}\bm{j}_{T}+\bm{\nabla}_{R}\bm{\cdot}\bm{j}_{R}=0,$
(2.41) $\displaystyle\bm{j}_{T}=\alpha
Da\bm{u}^{\prime}g+Pe_{s}\hat{U}_{s}(\bm{x})\bm{q}g-\bm{\nabla}g,$ (2.42)
$\displaystyle\bm{j}_{R}=\frac{1}{2}\alpha
Da\bm{\omega}^{\prime}g-\gamma^{2}\bm{\nabla}_{R}g,$ (2.43)
$\displaystyle\bm{n}\bm{\cdot}\bm{j}_{T}=0\quad\mbox{at}\quad r=\Delta,$
(2.44) $\displaystyle\int gd\bm{q}d\bm{x}=\frac{4\pi}{3}\Delta^{3},$ (2.45)

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

$\alpha=\frac{\overline{n}k_{s}T_{s}\tau_{D}}{\eta_{e}}.$ (2.46)

Physically, this is a comparison between the active driving pressure
($\overline{n}k_{s}T_{s}$) and a viscous resistive ‘pressure’
($\eta_{e}/\tau_{D}$) on the timescale $\tau_{D}$.

In the equations above, variables
$\\{\bm{u}^{e},p_{f}^{e},\bm{x},r,\bm{U},\bm{u}^{s},P,\bm{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 shall work with non-dimensional quantities unless otherwise
noted.

Non-dimensional parameter | Mathematical definition | Physical description
---|---|---
$\alpha$ | $\overline{n}k_{s}T_{s}\tau_{D}/\eta_{e}$ | Reduced osmotic pressure
$\beta$ | $\eta_{i}/\eta_{e}$ | Viscosity ratio
$\gamma$ | $R/\delta$ | Comparison of $R$ and $\delta$
$Da$ | $K_{\perp}/(R\ell_{m})$ | Darcy number
$Pe_{s}$ | $\tilde{U}_{s}\tau_{D}/R$ | Swim Péclet number
$\Delta$ | $(R-\ell_{m})/R$ | Radius ratio
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

$n(\bm{x})=\int_{S^{2}}gd\bm{q},$ (2.47)

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

$\bm{\nabla}\bm{\cdot}\bm{j}_{n}=0,$ (2.48a) $\bm{j}_{n}=\alpha
Da\bm{u}^{\prime}n+Pe_{s}\hat{U}_{s}(\bm{x})\bm{m}-\bm{\nabla}n.$ (2.48b)

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

$\bm{m}(\bm{x})=\int_{S^{2}}\bm{q}gd\bm{q}.$ (2.49)

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

$\bm{\nabla}\bm{\cdot}\bm{j}_{m}-\frac{1}{2}\alpha
Da\bm{\omega}^{i}\times\bm{m}+2\gamma^{2}\bm{m}=\bm{0},$ (2.50a)
$\bm{j}_{m}=\alpha
Da\bm{u}^{\prime}\bm{m}+Pe_{s}\hat{U}_{s}(\bm{x})\left(\bm{Q}+\frac{1}{3}n\bm{I}\right)-\bm{\nabla}\bm{m}$
(2.50b)

where

$\bm{Q}=\int_{S^{2}}\left(\bm{q}\bm{q}-\frac{1}{3}\bm{I}\right)gd\bm{q},$
(2.51)

is the trace-free nematic order tensor and $\bm{I}$ is the identity tensor of
rank two. The no-flux condition at the interior surface for the polar order
becomes $\bm{n}\bm{\cdot}\bm{j}_{m}=\bm{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}\bm{m}$ in equation (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 $\bm{Q}=\bm{0}$ is often considered in
the literature [Saintillan and Shelley, 2015, Yan and Brady, 2015]. 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 [Dulaney and Brady, 2020, Burkholder and
Brady, 2020, Peng and Brady, 2020]. 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 and Leal, 1998] has been shown to agree well with the full kinetic
theory and recently a numerical scheme has been developed to efficiently
evaluate the Bingham closure [Weady et al., 2022]. With this closure,
simulations with high spatial resolution are performed for active nematics. As
note by Weady et al. [2022], 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,
1993, Tailleur and Cates, 2008, Row and Brady, 2020]. To illustrate this
mechanism and motivate later discussions, we summarize the simple one-
dimensional (1D) result here. In the absence of external linear or angular
velocities, such as due to flow or orienting field, the governing equation in
1D for highly active ABPs is $\bm{\nabla}\bm{\cdot}(\hat{U}_{s}\bm{m})=0$,
where the diffusive term is neglected. The solution in 1D is simply
$\bm{m}=0$. Then, equation (2.50a) reduces to $n\hat{U}_{s}=const$. Further,
Row and Brady [2020] 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 Vesicle motion in the limit of weak interior flow

In many situations, the advection due to the interior fluid flow is much
weaker compared to the self-propulsion of the ABPs or its active swim
diffusion (small Péclet number), and we may neglect the effect of the fluid
velocity disturbance on the distribution of ABPs.

###  3.1 Governing equations

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

If the vesicle is non-permeable ($Da=0$), no external or internal flows can be
generated, and the vesicle remains stationary despite the nonuniform 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:

$\displaystyle\bm{u}^{e}$ $\displaystyle=$
$\displaystyle\bm{u}_{0}^{e}+Da\bm{u}_{1}^{e}+\cdot\cdot\cdot,$ (3.1)
$\displaystyle p_{f}^{e}$ $\displaystyle=$ $\displaystyle
p_{f,0}^{e}+Dap_{f,1}^{e}+\cdot\cdot\cdot,$ (3.2)
$\displaystyle\bm{u}^{\prime}$ $\displaystyle=$
$\displaystyle\bm{u}_{0}^{\prime}+Da\bm{u}_{1}^{\prime}+\cdot\cdot\cdot,$
(3.3) $\displaystyle P$ $\displaystyle=$ $\displaystyle
P_{0}+DaP_{1}+\cdot\cdot\cdot,$ (3.4) $\displaystyle g$ $\displaystyle=$
$\displaystyle g_{0}+Dag_{1}+\cdot\cdot\cdot.$ (3.5)

The dimensionless number density is given by $n=\int
gd\bm{q}=n_{0}+Dan_{1}+\cdot\cdot\cdot.$ Similarly, the expansions for the
translational and the seepage velocities are, respectively,

$\displaystyle\bm{U}$ $\displaystyle=$
$\displaystyle\bm{U}_{0}+Da\bm{U}_{1}+\cdot\cdot\cdot,$ (3.6)
$\displaystyle\bm{u}^{s}$ $\displaystyle=$
$\displaystyle\bm{u}_{0}^{s}+Da\bm{u}_{1}^{s}+\cdot\cdot\cdot.$ (3.7)

From equation (2.37), 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,

$\displaystyle\bm{u}_{0}^{s}=\left(p_{f}^{i}\big{\rvert}_{S_{i}}-p_{f}^{e}\big{\rvert}_{S_{e}}\right)\bm{n}.$
(3.8)

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

$\displaystyle\bm{\nabla}p_{f,0}^{e}=0,\quad\bm{\nabla}\bm{\cdot}\bm{u}_{0}^{e}=0.$
(3.9)

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

$p_{f,0}^{e}=0,\quad\bm{u}_{0}^{e}=\frac{\bm{u}_{0}^{s}}{r^{2}}.$ (3.10)

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

$\displaystyle\bm{\nabla}P_{0}=0,\quad\bm{\nabla}\bm{\cdot}\bm{u}_{0}^{\prime}=0.$
(3.11)

At the interior surface, the flow field satisfies the condition
$\bm{u}^{\prime}(r=\Delta)=\bm{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,

$p_{f,0}^{i}+6\gamma^{2}n_{0}/Pe_{s}^{2}=P_{0}=const.$ (3.12)

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

$\displaystyle\bm{\nabla}\bm{\cdot}\left(Pe_{s}\hat{U}_{s}(\bm{x})\bm{q}g_{0}-\bm{\nabla}g_{0}\right)-\gamma^{2}\bm{\nabla}_{R}^{2}g_{0}=0,$
(3.13)
$\displaystyle\bm{n}\bm{\cdot}\left(Pe_{s}\hat{U}_{s}(\bm{x})\bm{q}g_{0}-\bm{\nabla}g_{0}\right)=0\quad\mbox{at}\quad
r=\Delta,$ (3.14) $\displaystyle\int
g_{0}d\bm{q}d\bm{x}=\frac{4\pi}{3}\Delta^{3}.$ (3.15)

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

$\bm{U}_{0}=\frac{3\gamma^{2}}{2\pi
Pe_{s}^{2}}\int_{S^{2}}n_{0}(r=\Delta)\bm{n}d\Omega.$ (3.16)

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

$\bm{U}_{0}=\frac{1}{4\pi}\frac{K_{\perp}}{\eta_{e}\ell_{m}}\int_{S^{2}}\Pi^{\mathrm{osmo}}_{0}\bm{n}d\Omega,$
(3.17)

where $\Pi_{0}^{\mathrm{osmo}}=n^{w}k_{B}T$ is the dimensional osmotic
pressure of ABPs in the absence of flow.

To sum up, one needs to solve equations (3.13)-(3.15) to obtain the density
distribution of ABPs in the absence of flow, and then using equation (3.16) to
calculate the vesicle motion. In the remainder of section 3, the subscript ‘0’
(e.g., $g_{0}$, $\bm{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

$n_{0}(\Delta,\theta,\phi)=\sum_{l=0}^{\infty}\sum_{m=-l}^{m=l}C_{l,m}Y_{l}^{m}(\theta,\phi),$
(3.18)

where
$Y_{l}^{m}=\sqrt{(2l+1)(l-m)!/[4\pi(l+m)!]}P_{l}^{m}(\cos\theta)\exp(im\phi)$
and $P_{l}^{m}$ is the associated Legendre polynomial of degree $l$ and order
$m$. Using equation (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_{s}T_{s}/k_{B}T=\tilde{D}^{\mathrm{swim}}/D_{T}=\ell^{2}/(6\delta^{2})\gg
1$, which is often observed in active matter systems [Takatori et al., 2016].
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)}+\cdot\cdot\cdot$,
we obtain at leading order

$\frac{\ell}{R}\bm{\nabla}\bm{\cdot}\left[\hat{U}_{s}\bm{q}g^{(0)}\right]-\frac{1}{\hat{\tau}_{R}}\nabla^{2}_{R}g^{(0)}=0,$
(3.19)

where we have included the spatial variation of $\tau_{R}(\bm{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

$\bm{\nabla}\bm{\cdot}\left(\hat{U}_{s}\bm{m}^{(0)}\right)=0.$ (3.20)

Equation (3.19) is incompatible with the no-flux boundary condition and thus
is only valid 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,\bm{q})=g^{(-1)}/\epsilon+g^{(0)}+\cdot\cdot\cdot$. Defining
a stretched boundary-layer coordinate in the radial direction
$y=(\Delta-r)/\epsilon$, the Smoluchowski equation to leading order is

$\displaystyle\frac{\ell}{R}\hat{U}_{s}\big{\rvert}_{S_{i}}\bm{q}\bm{\cdot}\bm{e}_{r}\frac{\partial
g^{(-1)}}{\partial y}+\frac{\partial^{2}g^{(-1)}}{\partial y^{2}}$
$\displaystyle=$ $\displaystyle 0,$ (3.21)
$\displaystyle\frac{\ell}{R}\hat{U}_{s}\big{\rvert}_{S_{i}}\bm{q}\bm{\cdot}\bm{e}_{r}g^{(-1)}+\frac{\partial
g^{(-1)}}{\partial y}$ $\displaystyle=$ $\displaystyle 0\quad\mbox{at}\quad
y=0,$ (3.22) $\displaystyle g^{(-1)}\to 0\quad\mbox{as}\quad y\to+\infty.$
(3.23)

Here, the Taylor expansion
$\hat{U}_{s}(r,\theta,\phi)=\hat{U}_{s}\big{\rvert}_{S_{i}}-\epsilon
y\frac{d\hat{U}_{s}}{dr}\big{\rvert}_{S_{i}}+\cdot\cdot\cdot$ is used. The
solution is readily obtained

$g^{(-1)}=\begin{cases}A_{1}(\theta,\phi,\bm{q})\exp\left(-\frac{\ell}{R}\hat{U}_{s}\big{\rvert}_{S_{i}}\bm{q}\bm{\cdot}\bm{e}_{r}y\right),&\bm{q}\bm{\cdot}\bm{e}_{r}>0,\\\
0,&\mathrm{otherwise}.\end{cases}$ (3.24)

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

At the interior surface of the vesicle ($y=0$), the leading-order density is
large and given by
$\gamma^{2}\int_{\bm{q}\bm{\cdot}\bm{e}_{r}>0}A_{1}d\bm{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_{B}T/\left(n^{0}k_{s}T_{s}\right)$, 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 $\bm{n}$. From the above
analysis, we have

$\frac{n^{w}k_{B}T}{n^{0}k_{s}T_{s}}=\frac{\gamma^{2}\int_{\bm{q}\bm{\cdot}\bm{e}_{r}>0}A_{1}d\bm{q}}{\int
g^{(0)}(y\to\infty,\theta,\phi,\bm{q})d\bm{q}}\frac{k_{B}T}{k_{s}T_{s}}=f(\ell/R,\Delta),$
(3.25)

where $\gamma^{2}k_{B}T/k_{s}T_{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 analytically tractable, the factor $f(\ell/R,\Delta)$ in the preceding
equation cannot be explicitly obtained. Nevertheless, equation (3.25) reveals
the important fact that at high activity

$\Pi^{\mathrm{osmo}}=n^{w}k_{B}T=\Pi^{\mathrm{swim}}_{0}f(\ell/R,\Delta),$
(3.26)

where $\Pi^{\mathrm{swim}}_{0}=n^{0}k_{s}T_{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 $\Delta$. We emphasize that in equation (3.26), all
quantities are defined locally along the interior surface of the vesicle. This
is a generalization of the result of Yan and Brady [2015] for ABPs outside an
infinite planar wall, where $n^{w}k_{B}T=n^{0}k_{s}T_{s}$ in the limit
$\gamma^{2}\to\infty$ because of the absence of curvature of the geometry.

Figure 2: The magnitude of the dimensionless net force on the interior vesicle
surface $\bm{F}^{w}/(4\pi R_{i}^{2}\overline{n}k_{s}T_{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 hemisphere 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.

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

$\bm{U}=\frac{1}{4\pi}\frac{K_{\perp}}{\eta_{e}\ell_{m}}\int_{S^{2}}\Pi_{0}^{\mathrm{swim}}f(\ell/R,\Delta)\bm{n}d\Omega.$
(3.27)

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 is given in section 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 and
Melrose, 1993]. In this approach, the force exerted on the wall due to the
collision with ABPs is readily obtained. Consider a simulation of $N$ ABPs
that only interact with the boundary independently but not among themselves.
After a time step $\Delta t$, some particles might have moved outside the
interior wall. For particle $i$ that is now outside, we add a displacement
$\Delta\bm{x}_{i}$ to the particle such that after the move the particle is at
contact with the boundary. The total force exerted on the wall is then
$\bm{F}^{w}=-\zeta\sum_{i\in\mathcal{I}}\Delta\bm{x}_{i}/\Delta t$ where
$\mathcal{I}$ is the set of all particles that are outside the boundary before
the hard-sphere move. As seen in equation (3.17), the net speed of the vesicle
is proportional to the net force $\bm{F}^{w}$.

In figure 2, we show the dimensionless net force exerted on the interior
vesicle surface by the ABPs, $\bm{F}^{w}/(4\pi
R_{i}^{2}\overline{n}k_{s}T_{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

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

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.

As discussed earlier, in 1D the relation $nU_{s}=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 equation (3.28)
as an example, $n^{0}U_{s}=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_{s}T_{s}/k_{B}T\sim
n^{0}U_{s}\zeta\ell/k_{B}T$ and $n^{0}U_{s}=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 behavior of the bulk
density. Because only the ABPs at the interior surface contribute to the net
force, and they can only push against the boundary, this leads to the fact
that the net force is in the negative $x$ direction (to the left). If one only
had 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. equation (1.2)].

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.

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 equation (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 noninertial
jet propulsion has been proposed and studied in detail by Spagnolie and Lauga
[2010] as an alternate mechanism for the locomotion of microswimmers. In their
paper, the jetting velocity distribution of a microswimmer ($\bm{u}^{s}$) is
prescribed, and then the swim speed is determined from the reciprocal theorem.

Using the approximation $n^{0}U_{s}=const$ and the relation
$\Pi_{0}^{\mathrm{swim}}=n^{0}k_{s}T_{s}=n^{0}U_{s}\zeta\ell/6$, we see that
it is the variation of run length $\ell(\bm{x})$ that is responsible for the
net force on the vesicle interior surface and ultimately the vesicle motion.
Using equation (3.27), a Taylor series expansion about the center of the
vesicle leads to the scaling relation $\bm{U}\sim K_{\perp}R\zeta
n^{0}U_{s}\bm{\nabla}\ell/(\eta_{e}\ell_{m})$, where $\nabla\ell$ is the
gradient of the run length at the center of the vesicle.

###  3.3 A large vesicle

When the vesicle is large, the confinement is weak, $\ell/R\ll 1$, 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 equation (3.13)
equivalently as

$\bm{\nabla}\cdot\left[\frac{\ell}{R}\hat{U}_{s}(\bm{x})\bm{q}g-\left(\frac{\ell}{R}\right)^{2}\left(\frac{\ell}{\delta}\right)^{-2}\bm{\nabla}g\right]-\nabla_{R}^{2}g=0.$
(3.29)

In this section, 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)}+\cdot\cdot\cdot$ and the leading order equation
$\nabla_{R}^{2}g^{(0)}=0$. The solution in the bulk is then
$g^{(0)}(\bm{x},\bm{q})=n^{(0)}(\bm{x})/(4\pi)$. The boundary-layer thickness
is determined by a balance between the swimming and the diffusive fluxes,
which leads to the leading-order equation

$\displaystyle-\frac{\partial}{\partial\rho}\left(\frac{\ell}{\delta}\hat{U}_{s}\big{\rvert}_{S^{i}}\bm{q}\bm{\cdot}\bm{e}_{r}g^{(0)}+\frac{\partial}{\partial\rho}g^{(0)}\right)-\nabla_{R}^{2}g^{(0)}$
$\displaystyle=$ $\displaystyle 0,$ (3.30)
$\displaystyle\frac{\ell}{\delta}\hat{U}_{s}\big{\rvert}_{S^{i}}\bm{q}\bm{\cdot}\bm{e}_{r}g^{(0)}+\frac{\partial}{\partial\rho}g^{(0)}$
$\displaystyle=$ $\displaystyle 0\quad\mbox{at}\quad\rho=0.$ (3.31)

Here, we have used the stretched coordinate $\rho=(\Delta-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)}(\bm{x})d\bm{x}=4\pi\Delta^{3}/3$. In the absence of curvature
terms, just like the problem of ABPs on one side of an infinite planar wall
[Yan and Brady, 2015], the number density at the interior wall of the vesicle
at $O(1)$ can be determined analytically; the result is given by

$\frac{n^{w}}{n^{0}}=1+\frac{1}{6}\left(\frac{\ell}{\delta}\right)^{2}\hat{U}^{2}_{0}\big{\rvert}_{S^{i}}.$
(3.32)

In dimensional terms, this means that the osmotic pressure at the wall
$\Pi_{0}^{\mathrm{osmo}}=n^{w}k_{B}T=n^{0}k_{B}T+n^{0}k_{s}T_{s}\hat{U}_{s}^{2}$
where $n^{0}$ is the density outside the boundary layer. To determine $n^{0}$,
one needs to solve equation (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

$\bm{U}=\frac{1}{4\pi}\frac{K_{\perp}}{\eta_{e}\ell_{m}}\int_{S^{2}}\left[n^{0}k_{B}T+n^{0}k_{s}T_{s}\hat{U}_{s}^{2}\big{\rvert}_{S^{i}}\right]\bm{n}d\Omega.$
(3.33)

For a large vesicle, the accumulation boundary layer has a similar structure
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,
equation (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 and Brady [2014] 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 and
Brady [2014] 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
$\bm{j}_{R}=\Omega_{c}\bm{q}\times\hat{\bm{H}}g-D_{R}\bm{\nabla}_{R}g$, where
$\Omega_{c}$ characterizes the rate of reorientation due to the field and
$\hat{\bm{H}}$ is the direction of the field. When an ABP is aligned with the
field direction ($\bm{q}\parallel\hat{\bm{H}}$), the external torque vanishes.
The Smoluchowski equation (3.13) for ABPs with constant properties in the
presence of an orienting field is then

$\bm{\nabla}\bm{\cdot}\left(Pe_{s}\bm{q}g-\bm{\nabla}g\right)+\gamma^{2}\bm{\nabla}_{R}\bm{\cdot}\left(\chi_{R}\bm{q}\times\hat{\bm{H}}g-\bm{\nabla}_{R}g\right)=0,$
(3.34)

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}=\Omega_{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

$\frac{\ell}{R}\bm{q}\bm{\cdot}\bm{\nabla}g^{(0)}+\bm{\nabla}\bm{\cdot}\left(\chi_{R}\bm{q}\times\hat{\bm{H}}g^{(0)}-\bm{\nabla}_{R}g^{(0)}\right)=0.$
(3.35)

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.

Figure 4: (a): The magnitude of the dimensionless net force on the interior
wall $\bm{F}^{w}/(4\pi R_{i}^{2}\overline{n}k_{s}T_{s})$ as a function of the
field strength $\chi_{R}$ for different values of $\ell/R_{i}$. (b): The
rescaled net force, $\bm{F}^{w}\ell/(4\pi R_{i}^{3}\overline{n}k_{s}T_{s})$,
as a function of $\chi_{R}$ for different values of $\ell/R$. All data
collapse into one curve in panel (b). The values of $\ell/R_{i}$ in both
panels are the same and are thus only shown in (b). In both panels, the
translational diffusion is absent, $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.

Because equation (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, $\bm{F}^{w}/(4\pi R_{i}^{2}\overline{n}k_{s}T_{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{\bm{H}}$. In figure 4(b), the
same data is 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{\bm{H}}$ and are pushing against the vesicle; further
increasing of the field strength beyond this limit has no effect.

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 section 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
$\overline{n}\sim N/R_{i}^{3}$ and $k_{s}T_{s}\sim\zeta U_{s}^{2}\tau_{R}$, we
have

$\frac{F^{w}}{4\pi R_{i}^{2}\overline{n}k_{s}T_{s}}\sim\frac{N\zeta
U_{s}}{R_{i}^{2}\overline{n}k_{s}T_{s}}f(\chi_{R})\sim\frac{R_{i}}{\ell}f(\chi_{R}).$
(3.36)

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\pi R_{i}^{3}\overline{n}k_{s}T_{s})=2$, this 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 does the
speed of the vesicle.

##  4 Slow variation in activity

In the previous section, the dynamics of the vesicle is determined by the
distribution of ABPs in the absence of flow. To understand the effect of
interior fluid flow on the distribution of ABPs and the dynamics of the
vesicle, we consider the case of slow variation in activity. When the activity
gradient is small, any smooth variation of the swim speed can be approximated
by a Taylor series expansion about the origin. Here, we consider the first
effect of a small gradient by keeping the linear term only. The non-
dimensional swim speed can be written as

$\hat{U}_{s}(\bm{x})=1+\epsilon\bm{e}\bm{\cdot}\bm{x},$ (4.1)

where $\epsilon=|\bm{\nabla}U_{s}|R/U_{s}\ll 1$ and
$\bm{e}=\bm{\nabla}U_{s}/|\bm{\nabla}U_{s}|$ is a constant unit vector in the
direction of the gradient. If $\epsilon$ is identically zero, we have a
spatially homogeneous swim speed and there is no vesicle motion due to
spherical symmetry (see discussion in section 5). In this case of
$\epsilon\equiv 0$, the solution is
$\bm{u}_{0}^{\prime}=\bm{u}_{0}^{e}=\bm{u}_{0}^{s}=\bm{U}_{0}=\bm{0}$,
$p_{f,0}^{e}=0$ and $P_{0}=const$. The distribution of ABPs is governed by
equations (3.13), (3.14) and (3.15) but with $\hat{U}_{s}=1$, i.e., this
problem reduces to that of ABPs confined inside a fixed spherical domain. This
spherical symmetry means that the number density is a function of the radial
coordinate only, $n_{0}(\bm{x})=n_{0}(r)$. As shown by Yan and Brady [2015],
the number density is a monotonically increasing function that obtains its
maximum at the interior wall. Because the total pressure $P_{0}$ is a
constant, this variation of number density (osmotic pressure) maintains a
fluid pressure gradient with its maximum at the center of the interior domain.
The fluid pressure across the membrane is constant, and no seepage velocity is
generated.

To probe the first effect of a small linear gradient, we pose regular
expansions for all fields and the translational velocity:

$\displaystyle g$ $\displaystyle=$ $\displaystyle g_{0}+\epsilon
g_{1}+\cdot\cdot\cdot,$ (4.2)
$\displaystyle\left(P,p_{f}^{e},p_{f}^{i}\right)$ $\displaystyle=$
$\displaystyle(P_{0},0,0)+\epsilon\left(P_{1},p_{1}^{e},p_{1}^{\prime}\right)+\cdot\cdot\cdot,$
(4.3) $\displaystyle\left(\bm{u}^{\prime},\bm{u}^{e},\bm{u}^{s},\bm{U}\right)$
$\displaystyle=$
$\displaystyle\bm{0}+\epsilon\left(\bm{u}_{1}^{\prime},\bm{u}_{1}^{e},\bm{u}_{1}^{s},\bm{U}_{1}\right)+\cdot\cdot\cdot.$
(4.4)

At $O(\epsilon)$, the exterior fluid and the interior suspension are still
governed by equations (2.25)–(2.28) and (2.30)–(2.33). Similarly, the seepage
velocity is related to the jump in the fluid stress across the membrane given
by equation (2.38). The disturbance to the distribution of ABPs at this order
is governed by the inhomogeneous equation

$\displaystyle\bm{\nabla}\bm{\cdot}\left(Pe_{s}\bm{q}g_{1}-\bm{\nabla}g_{1}\right)-\gamma^{2}\nabla_{R}^{2}g_{1}$
$\displaystyle=$ $\displaystyle-\bm{\nabla}\bm{\cdot}\left(\alpha
Da\bm{u}_{1}^{\prime}g_{0}+Pe_{s}\bm{e}\bm{\cdot}\bm{x}\bm{q}g_{0}\right)$
(4.5) $\displaystyle-\frac{1}{2}\alpha
Da\bm{\nabla}_{R}\bm{\cdot}\left(\bm{\omega}_{1}^{\prime}g_{0}\right),$

with the boundary condition

$\bm{n}\bm{\cdot}\left(Pe_{s}\bm{q}g_{1}-\bm{\nabla}g_{1}\right)=-\alpha
Da\bm{n}\bm{\cdot}\bm{u}_{1}^{\prime}g_{0}-Pe_{s}\bm{e}\bm{\cdot}\bm{x}\bm{n}\bm{\cdot}\bm{q}g_{0}\quad\mbox{at}\quad
r=\Delta.$ (4.6)

The net disturbance is zero, $\int g_{1}d\bm{x}d\bm{q}=0$. As can be seen from
equation (4.5), the disturbance fields must be linear to the vector $\bm{e}$,
which allows us to write the number density in the form

$n_{1}=\bm{e}\bm{\cdot}\bm{x}h_{1}(r),$ (4.7)

where $h_{1}(r)$ is a scalar function of the radial coordinate only.

Due to linearity of the Stokes equations, the interior flow problem at
$O(\epsilon)$ admits a solution of the form

$\displaystyle P_{1}$ $\displaystyle=$ $\displaystyle
A_{1}\bm{e}\bm{\cdot}\bm{x},$ (4.8) $\displaystyle\bm{u}^{\prime}_{1}$
$\displaystyle=$ $\displaystyle
A_{2}\bm{e}+A_{3}\bm{e}\bm{\cdot}\left(\bm{x}\bm{x}-\frac{1}{3}r^{2}\bm{I}\right)+\frac{1}{2\beta
Da}P_{1}\bm{x}.$ (4.9)

Here, the momentum equation (2.30) is solved using a linear combination of the
growing tensor harmonic functions [Leal, 2007]. The continuity equation (2.31)
gives a constraint

$5A_{3}+\frac{3A_{1}}{\beta Da}=0.$ (4.10)

We can solve the external flow problem by considering two separate problems
with different boundary conditions: (1) $\bm{u}^{e}_{1}=\bm{u}_{1}^{s}$ and
(2) $\bm{u}^{e}_{1}=\bm{U}_{1}$ at $r=1$. Instead of solving the flow field
due to the second boundary condition in terms of the yet unknown velocity
$\bm{U}_{1}$, it will be determined from the reciprocal theorem (2.40). As a
result, one only needs to compute the exterior flow field due to the seepage
velocity $\bm{u}_{1}^{s}$. The exterior flow problem with the first boundary
condition has a solution of the form

$\displaystyle p_{1}^{e}$ $\displaystyle=$ $\displaystyle
A_{4}\bm{e}\bm{\cdot}\frac{\bm{x}}{r^{3}},$ (4.11)
$\displaystyle\bm{u}^{e}_{1}$ $\displaystyle=$ $\displaystyle
A_{5}\bm{e}\frac{1}{r}+A_{6}\bm{e}\bm{\cdot}\left(\frac{\bm{I}}{r^{3}}-3\frac{\bm{x}\bm{x}}{r^{5}}\right)+\frac{1}{2Da}p_{1}^{e}\bm{x},$
(4.12)

where the decaying tensor harmonic functions are used. To satisfy the
continuity equation (2.26), we must have

$A_{4}=2DaA_{5}.$ (4.13)

The seepage velocity connects the interior and exterior flow field via

$\bm{u}_{1}^{\prime}(\bm{x}=\Delta\bm{e}_{r})=\bm{u}_{1}^{s}=\bm{u}_{1}^{e}(\bm{x}=\bm{e}_{r}),$
(4.14)

which reduces to

$A_{2}-\frac{1}{3}\Delta^{2}A_{3}=A_{5}+A_{6}\quad\mbox{and}\quad
A_{3}\Delta^{2}+\frac{A_{1}\Delta^{2}}{2\beta Da}=-3A_{6}+\frac{A_{4}}{2Da}.$
(4.15)

The volume conservation (2.39) is satisfied. The velocity of the vesicle is
obtained from the reciprocal theorem, which gives

$\bm{U}_{1}=-\frac{1}{4\pi}\int_{S^{2}}\bm{u}_{1}^{s}d\Omega=-\left(A_{2}+\frac{A_{1}\Delta^{2}}{6\beta
Da}\right)\bm{e}.$ (4.16)

Finally, to solve equation (2.38) at this order, we need to compute the fluid
stress at the interior and the exterior wall. At the interior wall, we have

$\displaystyle\bm{\sigma}_{f,1}^{i}\bm{\cdot}\bm{e}_{r}=$
$\displaystyle-\Delta\left(A_{1}-h_{1}(\Delta)\frac{k_{B}T}{k_{s}T_{s}}\right)\bm{e}\bm{\cdot}\bm{e}_{r}\bm{e}_{r}$
(4.17) $\displaystyle+\Delta\left(\frac{7}{3}A_{3}\beta
Da+\frac{3}{2}A_{1}\right)\bm{e}\bm{\cdot}\bm{e}_{r}\bm{e}_{r}+\Delta\left(\frac{1}{3}A_{3}\beta
Da+\frac{1}{2}A_{1}\right)\bm{e}.$

The traction at the exterior wall has two contributions. The first is due to
the vesicle translating at a constant speed $\bm{U}_{1}$, which is given by
[Guazzelli and Morris, 2011, pp. 44]

$\bm{\sigma}_{U_{1}}^{e}\bm{\cdot}\bm{e}_{r}=-\frac{3}{2}Da\bm{U}_{1}.$ (4.18)

The second contribution is from the seepage velocity boundary condition
$\bm{u}_{1}^{s}$, which is given by

$-A_{4}\bm{e}\bm{\cdot}\bm{e}_{r}\bm{e}_{r}+Da\left(-A_{5}-6A_{6}+\frac{A_{4}}{2Da}\right)\bm{e}+Da\left(-A_{5}+18A_{6}-\frac{3A_{4}}{2Da}\right)\bm{e}\bm{\cdot}\bm{e}_{r}\bm{e}_{r}.$
(4.19)

Using equations (4.16)–(4.19) we can obtain the jump in the fluid stress
across the membrane, which then allows us to calculate the seepage velocity
using equation (2.38). Equating this result with the seepage velocity obtained
from equation (4.12) by setting $r=1$, we arrive at the following equations
for the coefficients:

$A_{5}+A_{6}=0,$ (4.20)

and

$\displaystyle\frac{A_{4}}{2Da}-3A_{6}=$ $\displaystyle
A_{1}\Delta\left(-1+\frac{\Delta}{4\beta}\right)+\frac{3}{2}DaA_{2}-\frac{8}{3}A_{3}\beta\Delta
Da-2A_{4}$ (4.21) $\displaystyle-2A_{5}Da+12A_{6}Da-\Delta
h_{1}(\Delta)\frac{k_{B}T}{k_{s}T_{s}}.$

Equation (4.20) implies that $\bm{u}_{1}^{s}$ is proportional to
$\bm{e}\bm{\cdot}\bm{e}_{r}\bm{e}_{r}$ and the component proportional to
$\bm{e}$ is zero, which is consistent with the fact that the seepage velocity
is in the normal ($\bm{e}_{r}$) direction. At this stage, we have obtained 6
equations for the 6 unknown coefficients $A_{i}$ ($i=1\cdot\cdot\cdot 6$),
which are given by equations (4.10), (4.13), (4.15), (4.20) and (4.21). Using
these equations, one could express $A_{i}$ in terms of the boundary value of
$h_{1}$ at the interior wall, i.e., $h_{1}(\Delta)$. These relations are
obtained as

$A_{6}=\frac{\Delta^{2}}{4}\frac{k_{B}T}{k_{s}T_{s}}\frac{h_{1}(\Delta)}{\Delta+Da(6\beta+4\Delta)},$
(4.22)

and

$\displaystyle A_{1}$ $\displaystyle=$
$\displaystyle\frac{40Da\beta}{\Delta^{2}}A_{6},\quad A_{2}=-8A_{6},\quad
A_{3}=-\frac{24}{\Delta^{2}}A_{6},$ (4.23) $\displaystyle A_{4}$
$\displaystyle=$ $\displaystyle-2DaA_{6},\quad A_{5}=-A_{6}.$ (4.24)

From equation (4.16), we have the net motion of the vesicle

$\bm{U}_{1}=\frac{4}{3}A_{6}\bm{e}=\frac{\Delta^{2}}{3}\frac{k_{B}T}{k_{s}T_{s}}\frac{h_{1}(\Delta)}{\Delta+Da(6\beta+4\Delta)}\bm{e}.$
(4.25)

Equation (4.25) is the main result of this section. In obtaining (4.25) the
only assumption made is the small linear gradient in the swim speed;
therefore, it applies generally for all ranges of the parameters $\alpha$,
$\beta$, $Da$, $Pe_{s}$ and $\gamma$. In particular, no restriction on the
activity of the ABPs (e.g., $\ell/\delta$) is made. We note that
$h_{1}(\Delta)$ depends parametrically on all the above parameters.

To obtain $h_{1}(r)$, we need to solve equation (4.5) that governs the
disturbance probability density distribution of the ABPs. As an approximation,
we consider the general solution using the $\bm{Q}=\bm{0}$ closure. At $O(1)$,
the spherical symmetry allows us to write the number density and polar order
in the form

$\displaystyle n_{0}(\bm{x})$ $\displaystyle=$ $\displaystyle n_{0}(r),$
(4.26) $\displaystyle\bm{m}_{0}(\bm{x})$ $\displaystyle=$
$\displaystyle\bm{x}f(r),$ (4.27)

which, when inserted into equations (2.48a) and (2.50a), leads to a couple of
ordinary differential equations (ODEs) for $n_{0}(r)$ and $f(r)$. The
solutions to $n_{0}$ and $\bm{m}_{0}$ under this assumption are obtained by
Yan and Brady [2015].

Next, we consider the disturbance distribution of ABPs at $O(\epsilon)$. At
this order, the number density distribution is governed by

$\bm{\nabla}\bm{\cdot}\bm{j}_{n,1}=0\quad\mbox{and}\quad\bm{j}_{n,1}=\alpha
Da\bm{u}_{1}^{\prime}n_{0}+Pe_{s}\bm{m}_{1}+Pe_{s}\bm{e}\bm{\cdot}\bm{x}\bm{m}_{0}-\bm{\nabla}n_{1}.$
(4.28)

The no-flux boundary condition is $\bm{n}\bm{\cdot}\bm{j}_{n,1}=0$ at
$r=\Delta$. Similarly, the governing equation for the polar order (assuming
$\bm{Q}_{1}=\bm{0}$) is

$\bm{\nabla}\bm{\cdot}\bm{j}_{m,1}+2\gamma^{2}\bm{m}_{1}-\frac{1}{2}\alpha
Da\bm{\omega}_{1}^{\prime}\times\bm{m}_{0}=0,$ (4.29)

and

$\bm{j}_{m,1}=\alpha
Da\bm{u}_{1}^{\prime}\bm{m}_{0}+\frac{1}{3}Pe_{s}\left(n_{0}\bm{e}\bm{\cdot}\bm{x}+n_{1}\right)\bm{I}-\bm{\nabla}\bm{m}_{1}.$
(4.30)

No-flux at $r=\Delta$ is $\bm{n}\bm{\cdot}\bm{j}_{m,1}=\bm{0}$. Similar to
equation (4.7), linearity and symmetry allows us to write the solution to the
polar order in the form

$\bm{m}_{1}=\bm{e}h_{2}(r)+\bm{e}\bm{\cdot}\bm{x}\bm{x}h_{3}(r),$ (4.31)

where $h_{2}(r)$ and $h_{3}(r)$ are functions of the radial coordinate only
and satisfy a coupled set of ODEs that can be found in Appendix B.

Figure 5: The dimensionless speed of the vesicle $U_{1}$ as a function of
$\ell/\delta$ for different fixed values of $\ell/R$. All other parameters are
fixed: $\Delta=0.98,\alpha=1,Da=0.1$ and $\beta=1.0$.

In figure 5 we show the dimensionless speed of the vesicle ($U_{1}$) as a
function of $\ell/\delta$ for $\ell/R=\\{0.1,1\\}$. With other dimensionless
parameters fixed, the increase of $\ell/\delta$ means the decrease of the
translational diffusivity and thus the increase of activity. The speed of the
vesicle vanishes as the activity approaches zero, $\ell/\delta\to 0$. As
$\ell/\delta$ increases, the speed of the vesicle increases and asymptotes to
a finite value for large $\ell/\delta$. The speed is larger for a smaller
$\ell/R$ because a thin boundary layer near the interior wall develops that
enhances the front-back asymmetry of the density distribution.

##  5 Concluding remarks

In this paper we have proposed a composite low-Reynolds-number propulsion
system made up of active Brownian particles encapsulated in a vesicle for the
purpose of enhanced transport beyond that of passive Brownian diffusion.
Instead of using the self-propulsion of a microswimmer directly, such as by
attaching a cargo to its surface, we considered an alternate mechanism in
which the vesicle is propelled by a fluid seepage velocity generated by a
concentration gradient of these encapsulated particles. In the present work,
we considered the cases in which the concentration gradient is generated by
either a prescribed activity gradient in the swim speed of these ABPs or an
external orienting field. By tuning the spatial pattern of variation in the
swim speed, one could obtain a concentration profile that in turn propels the
vesicle with a certain speed or in a desired direction. Alternatively, the
application of an external orienting field can push the ABPs against the wall
and generate net thrust for the vesicle. We provided a continuum formulation
governing the dynamics of the vesicle-ABPs system and explicitly analyzed its
behavior in the limits of weak interior flow and small activity gradient. For
the composite system as a whole, it moves by jet propulsion at low Reynolds
number, i.e., fluid is drawn in from one side of the vesicle and expels from
the other. The encapsulation of ABPs only provides a mechanism to generate
such a seepage flow.

We emphasize that in the present model it is the concentration gradient rather
than the species of the solute particles that is ultimately responsible for
vesicle locomotion. Any osmotic solute, not necessarily active, is able to
propel the vesicle so long as a concentration gradient is maintained. For a
passive solute, one can maintain a concentration gradient using chemical
reaction, e.g., by placing a distribution of sources and sinks. In this paper,
we analyzed how such a concentration gradient may be generated by an activity
gradient or by the application of an external orienting field. For
magnetotactic bacteria or synthetic active particles, an aligning magnetic
field can be used to control the direction of the concentration gradient and
therefore the direction of motion of the vesicle.

In an experimental setting, a spatial variation of the swim speed of
photokinetic bacteria can be achieved by exposing the bacteria to external
light intensity gradients. These light-powered bacteria exhibit a larger swim
speed in regions of higher light intensity. Under spatially patterned light
fields, light-responsive bacteria can self-assemble into reconfigurable
structures—‘painting’ with bacteria [Arlt et al., 2018, Frangipane et al.,
2018]. Another possible mechanism for inducing a spatially varying swim speed
could be the spatial modulation of ‘fuel’ (food sources).

For magnetotactic bacteria, instead of spatial modulation of swim speed one
can use an external static magnetic field that tends to align the bacteria in
a certain direction. For static or slowly-varying magnetic fields, the
magnitude of the induced electric field in this low frequency limit ($\ll 100$
kHz) is small so that its effect on the membrane dynamics is negligible [Ye
and Curcuru, 2015].

In obtaining the results we assumed that the ABPs can be treated as a
continuum and only contribute to the suspension stress via the osmotic
pressure. We note that additional constitutive models at the continuum level
for the suspension stress can be readily incorporated into our model. The
hydrodynamic interactions of the active particles with each other or the
confining vesicle boundary is neglected. These effects can be studied using a
colloidal approach by considering the detailed interactions among the active
particles and with the boundary. For example, this is considered in the study
of a single squirmer encapsulated in a porous container by Marshall and Brady
[2021] and for the case of a collection of squirmers inside a droplet that is
immersed in another fluid by Huang et al. [2020].

To achieve net motion of the spherical vesicle, a number density distribution
at the vesicle interior wall that breaks the front-back symmetry is required.
Instead of maintaining an asymmetric density distribution in a spherical
vesicle using ABPs with spatially-varying properties or ABPs with constant
properties but in an orienting field, one can also consider an asymmetric
vesicle. For ABPs with constant properties confined in an asymmetric
container, a symmetry-breaking density distribution will emerge because the
accumulation of ABPs at the wall depends on the local curvature. The effect of
vesicle shape on its net motion is left for a future study.

The enhancement of transport revealed by our study may be useful for the
development of synthetic microscale propelling systems for the purpose of
delivery of therapeutic payloads, penetrating complex media, or clearing
clogged arteries. We hope that our proposed theoretical designs can prompt new
experimental implementations.

## Funding

This work is supported by the National Science Foundation under Grant No. CBET
1803662.

##  Appendix A Brownian dynamics simulations

The dynamics of ABPs confined inside the vesicle in an external orienting
field can be resolved using Brownian dynamics (BD) simulations. Each ABP
follows the Langevin equations of motion given by

$\bm{0}=-\zeta\left(\bm{U}-U_{s}\bm{q}\right)+\bm{F}^{B}+\bm{F}^{w}\quad\mathrm{and}\quad\bm{0}=-\zeta_{R}\bm{\varOmega}+\bm{L}^{B}+\bm{L}^{e},$
(A.1)

where $\bm{U}$ ($\bm{\varOmega}$) is the instantaneous linear (angular)
velocity, $\bm{F}^{B}$ is the Brownian force, $\bm{F}^{w}$ is the hard-sphere
force due to collisions with the interior wall, $\zeta_{R}$ is the rotary
Stokes drag coefficient, $\bm{L}^{B}$ is the Brownian torque and $\bm{L}^{e}$
is the external torque due to the field.

The Brownian force and torque satisfy the white noise statistics:
$\overline{\bm{F}^{B}}=\bm{0},\overline{\bm{F}^{B}(0)\bm{F}^{B}(t)}=2k_{B}T\zeta\delta(t)\bm{I}$,
and
$\overline{\bm{L}^{B}}=\bm{0},\overline{\bm{L}^{B}(0)\bm{L}^{B}(t)}=2\zeta_{R}^{2}\delta(t)\bm{I}/\tau_{R}$.
Here, $\delta(t)$ is the delta function. In the BD simulations, the particle
orientations are represented using unit quaternions. At each time step, the
instantaneous particle velocities are computed and then used to update the
positions and orientations. The kinematic equation relating the angular
velocity and the rate-of-change of the quaternion is given by Delong et al.
[2015].

In figure 4, all data points are obtained by averaging over the long-time
behavior of the system. In each simulation, $10^{5}$ noninteracting ABPs are
used, and the system is evolved for a sufficiently long time such that the
steady state is reached.

##  Appendix B Equations for $h_{1},h_{2}$ and $h_{3}$

In this appendix we provide the detail on the derivation of the ODEs for
$h_{1},h_{2}$ and $h_{3}$. Note that the conservation

$\int_{|\bm{x}|\leq\Delta}n_{1}d\bm{x}=0$ (B.1)

is satisfied.

Note that

$\bm{\nabla}f(r)=\bm{x}\frac{1}{r}f^{\prime},$ (B.2)

and

$\bm{\nabla}(\bm{e}\bm{\cdot}\bm{x}f)=\bm{e}f+\bm{e}\bm{\cdot}\bm{x}\bm{x}\frac{1}{r}f^{\prime}.$
(B.3)

Using the identity

$\bm{\nabla}\bm{\cdot}(\underbrace{\bm{x}\bm{x}\cdot\cdot\cdot\bm{x}}_{k}f(r))=\left[(d+k-1)f+rf^{\prime}\right]\underbrace{\bm{x}\bm{x}\cdot\cdot\cdot\bm{x}}_{k-1},$
(B.4)

we can obtain

$\displaystyle\bm{\nabla}\bm{\cdot}[\bm{e}\bm{\cdot}\bm{x}\bm{x}f(r)]$
$\displaystyle=$
$\displaystyle\bm{e}\bm{\cdot}[\bm{\nabla}\bm{\cdot}(\bm{x}\bm{x}f(r))]=\bm{e}\bm{\cdot}\bm{x}(4f+rf^{\prime}),$
(B.5)
$\displaystyle\bm{\nabla}\bm{\cdot}[\bm{e}\bm{\cdot}\bm{x}\bm{x}\bm{x}f(r)]$
$\displaystyle=$
$\displaystyle\bm{e}\bm{\cdot}[\bm{\nabla}\bm{\cdot}(\bm{x}\bm{x}\bm{x}f(r))]=\bm{e}\bm{\cdot}\bm{x}\bm{x}(5f+rf^{\prime}).$
(B.6)

Similarly, we have

$\nabla^{2}f=\frac{2}{r}f^{\prime}+f^{\prime\prime},$ (B.8)
$\nabla^{2}[\bm{e}\bm{\cdot}\bm{x}f]=\bm{e}\bm{\cdot}\bm{x}\left(\frac{4f^{\prime}}{r}+f^{\prime\prime}\right),$
(B.9)
$\nabla^{2}(\bm{e}\bm{\cdot}\bm{x}\bm{x}f(r))=2\bm{e}f+\bm{e}\bm{\cdot}\bm{x}\bm{x}\left(\frac{6}{r}f^{\prime}+f^{\prime\prime}\right).$
(B.10)

The equation for $h_{1}$ is given by

$\displaystyle\alpha
Da\frac{dn_{0}}{dr}\left(\frac{1}{r}A_{2}+\frac{2}{3}rA_{3}+\frac{r}{2\beta
Da}A_{1}\right)+Pe_{s}\left(\frac{1}{r}\frac{dh_{2}}{dr}+4h_{3}+r\frac{dh_{3}}{dr}\right)$
$\displaystyle+Pe_{s}\left(4f+r\frac{df}{dr}\right)-\frac{4}{r}\frac{dh_{1}}{dr}-\frac{d^{2}h_{1}}{dr^{2}}=0.$
(B.11)

The no-flux condition is given by

$Pe_{s}(r^{2}f+h_{2}+r^{2}h_{3})-h_{1}-r\frac{dh_{1}}{dr}+A_{2}\alpha
Dan_{0}+\frac{\alpha}{6\beta}r^{2}n_{0}(3A_{1}+4A_{3}Da\beta)=0,$ (B.12)

evaluated at $r=\Delta$. The governing equation for $h_{2}$ is

$\displaystyle\alpha
Da\left(A_{2}-\frac{1}{3}r^{2}A_{3}\right)f+\frac{1}{3}Pe_{s}\left(n_{0}+h_{1}\right)-\frac{2}{r}\frac{dh_{2}}{dr}-\frac{d^{2}h_{2}}{dr^{2}}-2h_{3}$
$\displaystyle+2\gamma^{2}h_{2}+\frac{1}{2}\alpha
Da\left(\frac{5}{3}A_{3}+\frac{A_{2}}{2\beta Da}\right)r^{2}f=0.$ (B.13)

The no-flux condition at $r=\Delta$ is

$\frac{dh_{2}}{dr}=0.$ (B.14)

The governing equation for $h_{3}$ is

$\displaystyle\alpha Da\left(A_{3}+\frac{A_{1}}{2\beta Da}\right)f+\alpha
Da\frac{1}{r}\frac{df}{dr}\left(A_{2}+\frac{2}{3}r^{2}A_{3}+\frac{A_{1}r^{2}}{2\beta
Da}\right)$
$\displaystyle+\frac{1}{3}Pe_{s}\frac{1}{r}\left(\frac{dn_{0}}{dr}+\frac{dh_{1}}{dr}\right)-\frac{6}{r}\frac{dh_{3}}{dr}-\frac{d^{2}h_{3}}{dr^{2}}$
$\displaystyle+2\gamma^{2}h_{3}-\frac{1}{2}\alpha
Da\left(\frac{5}{3}A_{3}+\frac{A_{2}}{2\beta Da}\right)f=0.$ (B.15)

The no-flux condition is

$\alpha Darf\left(A_{2}+\frac{2}{3}r^{2}A_{3}+\frac{r^{2}A_{1}}{2\beta
Da}\right)+\frac{1}{3}Pe_{s}r(n_{0}+h_{1})-2rh_{3}-r^{2}\frac{dh_{3}}{dr}=0,$
(B.16)

evaluated at $r=\Delta$. We solve these equations in MATLAB using a Chebyshev
collocation method [Trefethen, 2000].

## References

  * Arlt et al. [2018] Jochen Arlt, Vincent A Martinez, Angela Dawson, Teuta Pilizota, and Wilson CK Poon.  Painting with light-powered bacteria.  _Nat. Commun._ , 9(1):1–7, 2018.  doi: 10.1038/s41467-018-03161-8.
  * Arlt et al. [2019] Jochen Arlt, Vincent A Martinez, Angela Dawson, Teuta Pilizota, and Wilson CK Poon.  Dynamics-dependent density distribution in active suspensions.  _Nat. Commun._ , 10(1):1–7, 2019.  doi: https://doi.org/10.1038/s41467-019-10283-0.
  * Beavers and Joseph [1967] Gordon S. Beavers and Daniel D. Joseph.  Boundary conditions at a naturally permeable wall.  _J. Fluid Mech._ , 30(1):197–207, 1967.  doi: 10.1017/S0022112067001375.
  * Blake [1971] J. R. Blake.  A spherical envelope approach to ciliary propulsion.  _J. Fluid Mech._ , 46(1):199–208, 1971.  doi: 10.1017/S002211207100048X.
  * Brenner [1964] H. Brenner.  The stokes resistance of a slightly deformed sphere.  _Chem. Eng. Sci._ , 19(8):519–539, 1964.  ISSN 0009-2509.  doi: 10.1016/0009-2509(64)85045-4.
  * Bunea and Taboryski [2020] Ada-Ioana Bunea and Rafael Taboryski.  Recent advances in microswimmers for biomedical applications.  _Micromachines_ , 11(12), 2020.  ISSN 2072-666X.  doi: 10.3390/mi11121048.
  * Burkholder and Brady [2020] Eric W. Burkholder and John F. Brady.  Nonlinear microrheology of active brownian suspensions.  _Soft Matter_ , 16:1034–1046, 2020.  doi: 10.1039/C9SM01713E.
  * Chaubal and Leal [1998] Charu V. Chaubal and L. Gary Leal.  A closure approximation for liquid-crystalline polymer models based on parametric density estimation.  _J. Rheol._ , 42(1):177–201, 1998.  doi: 10.1122/1.550887.
  * Delong et al. [2015] Steven Delong, Florencio Balboa Usabiaga, and Aleksandar Donev.  Brownian dynamics of confined rigid bodies.  _J. Chem. Phys._ , 143(14):144107, 2015.  doi: 10.1063/1.4932062.
  * Dulaney and Brady [2020] A. R. Dulaney and J. F. Brady.  Waves in active matter: The transition from ballistic to diffusive behavior.  _Phys. Rev. E_ , 101:052609, May 2020.  doi: 10.1103/PhysRevE.101.052609.
  * Durlofsky and Brady [1987] L. Durlofsky and J. F. Brady.  Analysis of the Brinkman equation as a model for flow in porous media.  _Phys. Fluids_ , 30(11):3329–3341, 1987.  doi: 10.1063/1.866465.
  * Elfring [2015] Gwynn J. Elfring.  A note on the reciprocal theorem for the swimming of simple bodies.  _Phys. Fluids_ , 27(2):023101, 2015.  doi: 10.1063/1.4906993.
  * Erkoc et al. [2019] Pelin Erkoc, Immihan C. Yasa, Hakan Ceylan, Oncay Yasa, Yunus Alapan, and Metin Sitti.  Mobile microrobots for active therapeutic delivery.  _Adv. Ther._ , 2(1):1800064, 2019.  doi: https://doi.org/10.1002/adtp.201800064.
  * Felfoul et al. [2016] Ouajdi Felfoul, Mahmood Mohammadi, Samira Taherkhani, Dominic De Lanauze, Yong Zhong Xu, Dumitru Loghin, Sherief Essa, Sylwia Jancik, Daniel Houle, Michel Lafleur, et al.  Magneto-aerotactic bacteria deliver drug-containing nanoliposomes to tumour hypoxic regions.  _Nat. Nanotechnol._ , 11(11):941–947, 2016.  doi: 10.1038/nnano.2016.137.
  * Frangipane et al. [2018] Giacomo Frangipane, Dario Dell’Arciprete, Serena Petracchini, Claudio Maggi, Filippo Saglimbeni, Silvio Bianchi, Gaszton Vizsnyiczai, Maria Lina Bernardini, and Roberto Di Leonardo.  Dynamic density shaping of photokinetic E. coli.  _eLife_ , 7:e36608, aug 2018.  ISSN 2050-084X.  doi: 10.7554/eLife.36608.
  * Gao et al. [2005] Xiaohu Gao, Lily Yang, John A Petros, Fray F Marshall, Jonathan W Simons, and Shuming Nie.  In vivo molecular and cellular imaging with quantum dots.  _Curr. Opin. Biotech._ , 16(1):63–72, 2005.  ISSN 0958-1669.  doi: https://doi.org/10.1016/j.copbio.2004.11.003.
  * Guazzelli and Morris [2011] Elisabeth Guazzelli and Jeffrey F Morris.  _A physical introduction to suspension dynamics_ , volume 45.  Cambridge University Press, 2011.
  * Heyes and Melrose [1993] D.M. Heyes and J.R. Melrose.  Brownian dynamics simulations of model hard-sphere suspensions.  _J. Non-Newtonian Fluid Mech._ , 46(1):1–28, 1993\.  ISSN 0377-0257.  doi: https://doi.org/10.1016/0377-0257(93)80001-R.
  * Huang et al. [2020] Zhihan Huang, Toshihiro Omori, and Takuji Ishikawa.  Active droplet driven by a collective motion of enclosed microswimmers.  _Phys. Rev. E_ , 102:022603, Aug 2020.  doi: 10.1103/PhysRevE.102.022603.
  * Jackson et al. [2020] Lisa A. Jackson, Evan J. Anderson, Nadine G. Rouphael, Paul C. Roberts, Mamodikoe Makhene, Rhea N. Coler, Michele P. McCullough, James D. Chappell, Mark R. Denison, Laura J. Stevens, Andrea J. Pruijssers, Adrian McDermott, Britta Flach, Nicole A. Doria-Rose, Kizzmekia S. Corbett, Kaitlyn M. Morabito, Sijy O’Dell, Stephen D. Schmidt, Phillip A. Swanson, Marcelino Padilla, John R. Mascola, Kathleen M. Neuzil, Hamilton Bennett, Wellington Sun, Etza Peters, Mat Makowski, Jim Albert, Kaitlyn Cross, Wendy Buchanan, Rhonda Pikaart-Tautges, Julie E. Ledgerwood, Barney S. Graham, and John H. Beigel.  An mRNA vaccine against SARS-CoV-2 — preliminary report.  _N. Engl. J. Med._ , 383(20):1920–1931, 2020\.  doi: 10.1056/NEJMoa2022483.
  * Jones [1973] I. P. Jones.  Low reynolds number flow past a porous spherical shell.  _Proc. Camb. Phil. Soc._ , 73(1):231–238, 1973\.  doi: 10.1017/S0305004100047642.
  * Kaiser et al. [2014] Andreas Kaiser, Anton Peshkov, Andrey Sokolov, Borge ten Hagen, Hartmut Löwen, and Igor S. Aranson.  Transport powered by bacterial turbulence.  _Phys. Rev. Lett._ , 112:158101, Apr 2014.  doi: 10.1103/PhysRevLett.112.158101.
  * Leal [2007] L Gary Leal.  _Advanced transport phenomena: fluid mechanics and convective transport processes_ , volume 7.  Cambridge University Press, 2007.
  * Lebedev et al. [2007] V. V. Lebedev, K. S. Turitsyn, and S. S. Vergeles.  Dynamics of nearly spherical vesicles in an external flow.  _Phys. Rev. Lett._ , 99:218101, Nov 2007.  doi: 10.1103/PhysRevLett.99.218101.
  * Lighthill [1952] MJ Lighthill.  On the squirming motion of nearly spherical deformable bodies through liquids at very small Reynolds numbers.  _Commun. Pure Appl. Math._ , 5(2):109–118, 1952\.  doi: 10.1002/cpa.3160050201.
  * Marchetti et al. [2013] M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, Madan Rao, and R. Aditi Simha.  Hydrodynamics of soft active matter.  _Rev. Mod. Phys._ , 85:1143–1189, Jul 2013.  doi: 10.1103/RevModPhys.85.1143.
  * Marshall and Brady [2021] Kevin J. Marshall and John F. Brady.  The hydrodynamics of an active squirming particle inside of a porous container.  _J. Fluid Mech._ , 919:A31, 2021.  doi: 10.1017/jfm.2021.276.
  * Masoud and Stone [2019] Hassan Masoud and Howard A. Stone.  The reciprocal theorem in fluid dynamics and transport phenomena.  _J. Fluid Mech._ , 879:P1, 2019.  doi: 10.1017/jfm.2019.553.
  * Medina-Sánchez et al. [2018] Mariana Medina-Sánchez, Haifeng Xu, and Oliver G Schmidt.  Micro- and nano-motors: the new generation of drug carriers.  _Ther. Deliv._ , 9(4):303–316, 2018.  doi: 10.4155/tde-2017-0113.
  * Park et al. [2017] Byung-Wook Park, Jiang Zhuang, Oncay Yasa, and Metin Sitti.  Multifunctional bacteria-driven microswimmers for targeted active drug delivery.  _ACS Nano_ , 11(9):8910–8923, 2017.  doi: 10.1021/acsnano.7b03207.
  * Pattni et al. [2015] Bhushan S. Pattni, Vladimir V. Chupin, and Vladimir P. Torchilin.  New developments in liposomal drug delivery.  _Chem. Rev._ , 115(19):10938–10966, 2015.  doi: 10.1021/acs.chemrev.5b00046.
  * Peng and Brady [2020] Zhiwei Peng and John F. Brady.  Upstream swimming and Taylor dispersion of active Brownian particles.  _Phys. Rev. Fluids_ , 5:073102, Jul 2020.  doi: 10.1103/PhysRevFluids.5.073102.
  * Rao et al. [2007] Jianghong Rao, Anca Dragulescu-Andrasi, and Hequan Yao.  Fluorescence imaging in vivo: recent advances.  _Curr. Opin. Biotech._ , 18(1):17–25, 2007.  ISSN 0958-1669.  doi: https://doi.org/10.1016/j.copbio.2007.01.003.
  * Row and Brady [2020] Hyeongjoo Row and John F. Brady.  Reverse osmotic effect in active matter.  _Phys. Rev. E_ , 101:062604, Jun 2020.  doi: 10.1103/PhysRevE.101.062604.
  * Saffman [1971] P. G. Saffman.  On the boundary condition at the surface of a porous medium.  _Stud. Appl. Math._ , 50(2):93–101, 1971.  doi: https://doi.org/10.1002/sapm197150293.
  * Saintillan and Shelley [2015] David Saintillan and Michael J. Shelley.  _Theory of Active Suspensions_ , pages 319–355.  Springer New York, New York, NY, 2015.  ISBN 978-1-4939-2065-5.  doi: 10.1007/978-1-4939-2065-5_9.
  * Schnitzer [1993] Mark J. Schnitzer.  Theory of continuum random walks and application to chemotaxis.  _Phys. Rev. E_ , 48:2553–2568, Oct 1993.  doi: 10.1103/PhysRevE.48.2553.
  * Singh et al. [2017] Ajay Vikram Singh, Zeinab Hosseinidoust, Byung-Wook Park, Oncay Yasa, and Metin Sitti.  Microemulsion-based soft bacteria-driven microswimmers for active cargo delivery.  _ACS Nano_ , 11(10):9759–9769, 2017.  doi: 10.1021/acsnano.7b02082.
  * Singh et al. [2019] Ajay Vikram Singh, Mohammad Hasan Dad Ansari, Peter Laux, and Andreas Luch.  Micro-nanorobots: important considerations when developing novel drug delivery platforms.  _Expert Opin. Drug Deliv._ , 16(11):1259–1275, 2019.  doi: 10.1080/17425247.2019.1676228.
  * Sokolov et al. [2010] Andrey Sokolov, Mario M. Apodaca, Bartosz A. Grzybowski, and Igor S. Aranson.  Swimming bacteria power microscopic gears.  _Proc. Natl. Acad. Sci. U.S.A._ , 107(3):969–974, 2010.  ISSN 0027-8424.  doi: 10.1073/pnas.0913015107.
  * Spagnolie and Lauga [2010] Saverio E. Spagnolie and Eric Lauga.  Jet propulsion without inertia.  _Phys. Fluids_ , 22(8):081902, 2010.  doi: 10.1063/1.3469786.
  * Stone and Samuel [1996] Howard A. Stone and Aravinthan D. T. Samuel.  Propulsion of microorganisms by surface distortions.  _Phys. Rev. Lett._ , 77:4102–4104, Nov 1996.  doi: 10.1103/PhysRevLett.77.4102.
  * Stroka et al. [2014] Kimberly M. Stroka, Hongyuan Jiang, Shih-Hsun Chen, Ziqiu Tong, Denis Wirtz, Sean X. Sun, and Konstantinos Konstantopoulos.  Water permeation drives tumor cell migration in confined microenvironments.  _Cell_ , 157(3):611–623, 2014.  ISSN 0092-8674.  doi: https://doi.org/10.1016/j.cell.2014.02.052.
  * Tailleur and Cates [2008] J. Tailleur and M. E. Cates.  Statistical mechanics of interacting run-and-tumble bacteria.  _Phys. Rev. Lett._ , 100:218103, May 2008.  doi: 10.1103/PhysRevLett.100.218103.
  * Takatori et al. [2014] S. C. Takatori, W. Yan, and J. F. Brady.  Swim pressure: Stress generation in active matter.  _Phys. Rev. Lett._ , 113:028103, Jul 2014.  doi: 10.1103/PhysRevLett.113.028103.
  * Takatori and Brady [2014] Sho C. Takatori and John F. Brady.  Swim stress, motion, and deformation of active matter: effect of an external field.  _Soft Matter_ , 10:9433–9445, 2014.  doi: 10.1039/C4SM01409J.
  * Takatori and Sahu [2020] Sho C. Takatori and Amaresh Sahu.  Active contact forces drive nonequilibrium fluctuations in membrane vesicles.  _Phys. Rev. Lett._ , 124:158102, Apr 2020.  doi: 10.1103/PhysRevLett.124.158102.
  * Takatori et al. [2016] Sho C Takatori, Raf De Dier, Jan Vermant, and John F Brady.  Acoustic trapping of active matter.  _Nat. Commun._ , 7(1):1–7, 2016.  doi: 10.1038/ncomms10694.
  * Torchilin [2012] Vladimir P. Torchilin.  Multifunctional nanocarriers.  _Adv. Drug Deliv. Rev._ , 64:302–315, 2012.  ISSN 0169-409X.  doi: https://doi.org/10.1016/j.addr.2012.09.031.
  * Trantidou et al. [2018] Tatiana Trantidou, Linda Dekker, Karen Polizzi, Oscar Ces, and Yuval Elani.  Functionalizing cell-mimetic giant vesicles with encapsulated bacterial biosensors.  _Interface Focus_ , 8(5):20180024, 2018.  doi: 10.1098/rsfs.2018.0024.
  * Trefethen [2000] Lloyd N. Trefethen.  _Spectral Methods in MATLAB_.  Society for Industrial and Applied Mathematics, 2000.  doi: 10.1137/1.9780898719598.
  * Vlahovska and Gracia [2007] Petia M. Vlahovska and Ruben Serral Gracia.  Dynamics of a viscous vesicle in linear flows.  _Phys. Rev. E_ , 75:016313, Jan 2007.  doi: 10.1103/PhysRevE.75.016313.
  * Vutukuri et al. [2020] Hanumantha Rao Vutukuri, Masoud Hoore, Clara Abaurrea-Velasco, Lennard van Buren, Alessandro Dutto, Thorsten Auth, Dmitry A Fedosov, Gerhard Gompper, and Jan Vermant.  Active particles induce large shape deformations in giant lipid vesicles.  _Nature_ , 586(7827):52–56, 2020.  doi: https://doi.org/10.1038/s41586-020-2730-x.
  * Weady et al. [2022] Scott Weady, Michael J. Shelley, and David B. Stein.  A fast chebyshev method for the bingham closure with application to active nematic suspensions.  _J. Comput. Phys._ , 457:110937, 2022.  ISSN 0021-9991.  doi: 10.1016/j.jcp.2021.110937.
  * West and Halas [2003] Jennifer L. West and Naomi J. Halas.  Engineered nanomaterials for biophotonics applications: Improving sensing, imaging, and therapeutics.  _Annu. Rev. Biomed. Eng._ , 5(1):285–292, 2003\.  doi: 10.1146/annurev.bioeng.5.011303.120723.
  * Yan and Brady [2015] Wen Yan and John F. Brady.  The force on a boundary in active matter.  _J. Fluid Mech._ , 785:R1, 2015.  doi: 10.1017/jfm.2015.621.
  * Yan and Brady [2018] Wen Yan and John F. Brady.  The curved kinetic boundary layer of active matter.  _Soft Matter_ , 14:279–290, 2018.  doi: 10.1039/C7SM01643C.
  * Ye and Curcuru [2015] Hui Ye and Austen Curcuru.  Vesicle biomechanics in a time-varying magnetic field.  _BMC Biophys._ , 8(1):1–14, 2015.
  * Zampogna and Gallaire [2020] Giuseppe A. Zampogna and François Gallaire.  Effective stress jump across membranes.  _J. Fluid Mech._ , 892:A9, 2020.  doi: 10.1017/jfm.2020.144.