Creeping motion of a solid particle inside a spherical elastic cavity. II. Asymmetric motion

TL;DR

This paper develops an analytical method to compute the hydrodynamic mobility of a small particle moving asymmetrically inside an elastic spherical cavity, revealing how membrane properties influence particle and cavity motion.

Contribution

It extends previous work by providing a general solution for arbitrary force directions and analyzing asymmetric motion within elastic cavities.

Findings

Particle mobility is higher in elastic cavities than in rigid ones.

Cavity translation is governed solely by membrane shear properties.

Analytical results agree with boundary integral simulations.

Abstract

An analytical method is proposed for computing the low-Reynolds-number hydrodynamic mobility function of a small colloidal particle asymmetrically moving inside a large spherical elastic cavity, the membrane of which is endowed with resistance toward shear and bending. In conjunction with the results obtained in the first part [Daddi-Moussa-Ider, L\"{o}wen, and Gekle, Eur. Phys. J. E 41, 104 (2018)], in which the axisymmetric motion normal to the surface of an elastic cavity is investigated, the general motion for an arbitrary force direction can be addressed. The elastohydrodynamic problem is formulated and solved using the classic method of images through expressing the hydrodynamic flow fields as a multipole expansion involving higher-order derivatives of the free-space Green's function. In the quasi-steady limit, we demonstrate that the particle self-mobility function of a particle moving tangent to the surface of the cavity is larger than that predicted inside a rigid stationary cavity of equal size. This difference is justified by the fact that a stationary rigid cavity introduces additional hindrance to the translational motion of the encapsulated particle, resulting in a reduction of its hydrodynamic mobility. Furthermore, the motion of the cavity is investigated, revealing that the translational pair (composite) mobility, which linearly couples the velocity of the elastic cavity to the force exerted on the solid particle, is solely determined by membrane shear properties. Our analytical predictions are favorably compared with fully-resolved computer simulations based on a completed-double-layer boundary integral method.

Equations

$$\eta {\nabla }^2{\mathbf{v}}^{(i)}-\nabla p^{(i)}+\mathbf{F}\delta \left(\mathbf{x}-{\mathbf{x}}_2\right)$$

$$\nabla \cdot {\mathbf{v}}^{(i)}$$

$$\eta {\nabla }^2{\mathbf{v}}^{(o)}-\nabla p^{(o)}$$

$$\nabla \cdot {\mathbf{v}}^{(o)}$$

$$\left|{\mathbf{v}}^{(i)}\right|<\infty$$

$${\mathbf{v}}^{(o)}\to 0$$

$${{\mathbf{v}}^{(o)}-{\mathbf{v}}^{(i)}|}_{r=a}$$

$${\left({\mathbf{\sigma }}^{(o)}-{\mathbf{\sigma }}^{(i)}\right)\cdot {\mathbf{e}}_r|}_{r=a}$$

$${\sigma }_{\theta r}$$

$${\sigma }_{\varphi r}$$

$${\sigma }_{rr}$$

$$E=\frac{{\kappa }_{\mathrm{S}}}{12}\left((I_1^2+2I_1-2I_2)+C{I}_2^2\right),$$

$$\Delta {\mathbf{f}}^{\mathrm{S}}={\Delta {\mathbf{f}}^{\mathrm{S}}|}_{\mathrm{Axi}}+{\Delta {\mathbf{f}}^{\mathrm{S}}|}_{\mathrm{Asy}}$$

$${\Delta f_{\theta }^{\mathrm{S}}|}_{\mathrm{Axi}}$$

$$\displaystyle\quad-u_{\theta}\left({\lambda}\cot^{2}\theta+{\lambda}-1\right)\bigg{)}\,,$$

$${\Delta f_{\varphi }^{\mathrm{S}}|}_{\mathrm{Axi}}$$

$${\Delta f_r^{\mathrm{S}}|}_{\mathrm{Axi}}$$

$${\Delta f_{\theta }^{\mathrm{S}}|}_{\mathrm{Asy}}$$

$${\Delta f_{\varphi }^{\mathrm{S}}|}_{\mathrm{Asy}}$$

$$\displaystyle\quad+\left(1-\cot^{2}\theta\right)\frac{u_{\phi}}{2}+\frac{u_{\phi,\theta}}{2}\cot\theta+\xi_{+}\frac{\cot\theta}{\sin\theta}\,u_{\theta,\phi}\bigg{)}\,,$$

$${\Delta f_r^{\mathrm{S}}|}_{\mathrm{Asy}}$$

$$\lambda :$$

$${\xi }_{\pm }$$

$$E_{\mathrm{B}}=2{\kappa }_{\mathrm{B}}{\left(H-H_0\right)}^2,$$

$$\Delta \mathbf{f}=-2{\kappa }_{\mathrm{B}}\left(2(H^2-K+H_{0}H)+{\Delta }_{\parallel }\right)(H-H_0)\mathbf{n},$$

$$\Delta f_r^{\mathrm{B}}={\Delta f_r^{\mathrm{B}}|}_{\mathrm{Axi}}+{\Delta f_r^{\mathrm{B}}|}_{\mathrm{Asy}},$$

$${\Delta f_r^{\mathrm{B}}|}_{\mathrm{Axi}}$$

$$\displaystyle\quad+\left(2-T^{2}\right)u_{r,\theta\theta}+2Tu_{r,\theta\theta\theta}+u_{r,\theta\theta\theta\theta}\bigg{)}\,,$$

$${\Delta f_r^{\mathrm{B}}|}_{\mathrm{Asy}}$$

$$\displaystyle\quad-2Tu_{r,\phi\phi\theta}+\left(1+T^{2}\right)u_{r,\phi\phi\phi\phi}\Big{)}\,,$$

$${\mathbf{v}|}_{r=a}=\frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t},$$

$${\mathbf{v}|}_{r=a}=i\omega \mathbf{u}.$$

$${\mathbf{v}}^{(i)}={\mathbf{v}}^{\mathrm{S}}+{\mathbf{v}}^{\ast },$$

$${\mathbf{v}}^{\mathrm{S}}=\mathcal{G}\left(\mathbf{x}-{\mathbf{x}}_2\right)\cdot \mathbf{F}=\frac{1}{8\pi \eta }\left(\frac{1}{s}+\frac{\mathbf{s}\mathbf{s}}{s^3}\right)\cdot \mathbf{F},$$

$$\frac{1}{s}=\sum _{n=0}^{\infty }{R}^n{\phi }_n(r,\theta ).$$

$${\phi }_n(r,\theta ):=\frac{(\mathbf{d}\cdot \nabla {)}^n}{n!}\frac{1}{r}=\frac{1}{r^{n+1}}{P}_n(\cos \theta ),$$

$$\frac{\mathbf{s}\mathbf{s}}{s^3}=\mathbf{s}{\nabla }_2\left(\frac{1}{s}\right),$$

$${\nabla }_2{R}^n=-n{R}^{n-1}\mathbf{d},\left(\mathbf{d}\cdot {\nabla }_2\right)\mathbf{d}=0.$$

Full text

I Introduction

Many industrial and biological transport processes on the microscale predominantly occur under confinement, where hydrodynamic interactions with boundaries drastically alter the diffusive behavior of microparticles in viscous media. Prime examples include particle sorting in microfabricated fluidic devices stone04 ; fu99 ; lu04 ; huh05 ; schmid14 , membrane separation and purification in pharmaceutical industry darvishmanesh11 ; adamo13 ; gutmann15 , as well as intracellular drug delivery and targeting via multifunctional nanocarriers, which release therapeutic agents in disease regions such as tumor or inflammation sites colson12 ; hillaireau09 ; liu16 ; maeda13 ; naahidi13 ; rosenholm10 ; singh09 ; bareford07 . The uptake by cell membranes occurs via endocytosis or by direct penetration to reach target cellular compartments.

At small length scales, fluid flows are characterized by small Reynolds numbers, implying that viscous forces dominate inertial forces. In these situations, the fluid-mediated hydrodynamic interactions are fully encoded in the mobility tensor, which linearly couples the velocities of microparticles to the forces and torques exerted on them happel12 ; kim13 ; leal80 . Even for simple geometric confinements, finding closed analytical solutions of diverse flow problems can be challenging. Most theoretical approaches are based on the method of images, consisting of a set of (typically higher-order) singularities that are required to satisfy the prescribed boundary conditions at the confining boundaries blake71 . Using this approach, the solution of the Stokes equations in the presence of a point force singularity acting in a fluid domain bounded by a rigid spherical cavity has been obtained by Oseen oseen28 . Extensions of Oseen’s solution have further been proposed butler53 ; collins54 ; hasimoto56 ; hasimoto92 ; hasimoto97 ; shail87 ; shail88 ; sellier08 . A particularly more compliant solution that separately considers both axisymmetric and asymmetric Stokeslets has later been presented by Maul and Kim maul94 ; maul96 . Meanwhile, the hydrodynamic coupling and rotational mobilities have been calculated for point-like particles felderhof12b . In this context, the low-Reynolds-number swimming inside spherical containers has also attracted some attention tsemakh04 ; lavrenteva05 ; reigh17 ; zhu17 ; reigh17prfluids ; shaik18 .

In this manuscript, we examine the slow translational motion of a small colloidal particle moving inside a large spherical elastic cavity (that itself is floating in an infinitely-extended viscous fluid). This setup may be viewed as a relevant model system for transport processes within biological media, such as elastic cell membranes. The cavity membrane is modeled as a two-dimensional hyperelastic sheet, endowed with resistance toward shear elasticity and bending rigidity. This model has previously been employed to address the effect of elastic confinements on the diffusive behavior of colloidal particles moving close to planar daddi16 ; daddi18epje or curved elastic membranes daddi17b ; daddi17c ; daddi17d ; daddi17e .

The present article is a natural extension of a preceding paper daddi18cavity (hereinafter referred to as part I), where the axisymmetric motion was examined. The goal of the current study is to supplement and complement our previous results by quantifying the effect of the confining elastic cavity on the asymmetric motion of an encapsulated particle located at arbitrary position within the cavity. Our approach is based on the method of images employed by Fuentes and collaborators fuentes88 ; fuentes89 , who examined theoretically the hydrodynamic interactions between two unequally-sized spherical viscous drops at moderately small separations. Our analytical investigations proceed through the calculation of the Green’s functions associated with a point force acting inside a spherical elastic cavity. The problem treated here does not possess the symmetry properties of the simpler axisymmetric case considered in part I. This makes it necessary to employ an alternative mathematical framework to obtain the solution of the flow problem for the asymmetric case. The calculated hydrodynamic flow field is used to determine the frequency-dependent mobility functions for an enclosed point particle. This approximation is reasonable if the separation distance between the particle and the cavity surface is large compared to the particle size. Particularly, inside our deformable cavity, the mobility in the quasi-steady limit of vanishing frequency is shown to be always larger than the one predicted inside a rigid cavity with no-slip boundary condition. Our theoretical results favorably compare to numerical simulations.

The remainder of the paper is organized as follows. In sect. II, we use the multipole expansion method to find solutions of the elastohydrodynamic problem for the fluid inside and outside the cavity. We then provide in sect. III analytical expressions of the hydrodynamic self-mobility function for a particle moving tangent to the surface of the cavity. In sect. IV, we assess the motion of the large cavity and determine the deformation field induced by the motion of the particle. We provide in sect. V concluding remarks summarizing our findings. The appendix contains explicit expressions for the series coefficients arising from the multipole expansion.

II Singularity solution

We examine the low-Reynolds-number motion of a small sphere of radius $b$ situated inside a large spherical elastic cavity of radius $a$. The fluid inside and outside the cavity is characterized by a constant dynamic viscosity $\eta$, and the flow is assumed to be incompressible. The center of the cavity at $\bm{x}_{1}$ coincides with the origin of the spherical coordinate system. The solid particle located at position $\bm{x}_{2}=R\bm{e}_{z}$ is moving under the action of an asymmetric external force $\bm{F}\perp\bm{e}_{z}$. An illustration of the system under consideration is shown in fig. 1.

The physical problem is thus equivalent to solving the forced Stokes equations inside the cavity happel12 ; kim13 ,

[TABLE]

and homogeneous (force-free) equations for the outer fluid,

[TABLE]

wherein $\bm{v}^{(i)}$ and $\bm{v}^{(o)}$ denote the flow velocity fields for the inner and outer fluids, respectively, and $p^{(i)}$ and $p^{(o)}$ are the corresponding pressure fields. Equations (1) and (2) are subject to the regularity conditions

[TABLE]

in addition to the standard boundary conditions of continuity of the velocity field and discontinuity of the hydrodynamic stresses at the cavity surface. In the present work, we assume that the cavity undergoes a small deformation only, so that the boundary conditions are evaluated at the undeformed surface of reference at $r=a$. Specifically,

[TABLE]

where $\bm{\sigma}=-p\bm{I}+2\eta\bm{E}$ is the viscous stress tensor. Here, $\bm{E}=\left(\bm{\nabla}\bm{v}+\bm{\nabla}\bm{v}^{\mathrm{T}}\right)/2$ denotes the rate-of-strain tensor, the components of which are given in spherical coordinates by

[TABLE]

where $\phi$ and $\theta$, respectively, denote the azimuthal and polar angles, such that $(\phi,\theta)\in[0,2\pi)\times[0,\pi]$ describes a point on the surface of the unit sphere. Furthermore, by convention, indices after a comma stand for the corresponding partial derivatives, e.g., $v_{r,r}=\partial v_{r}/\partial r$. Additionally, $\Delta\bm{f}^{\mathrm{S}}$ and $\Delta\bm{f}^{\mathrm{B}}$ denote the traction jumps stemming from shear and bending deformation modes, respectively. We further remark that, if the membrane cavity undergoes a large deformation, the boundary conditions should rather be evaluated at the displaced membrane positions, see, e.g., refs. sekimoto93, ; weekley06, ; salez15, ; saintyves16, ; rallabandi17, ; daddi18stone, ; rallabandi18, .

In this work, we model the elastic cavity as a spherical hyperelastic shell of vanishing thickness, the deformation of which is governed by the shear elasticity model proposed by Skalak skalak73 that is commonly employed when modeling, e.g., the membranes of red blood cells Freund_2014 ; krueger12 . Specifically, the areal strain energy density of the Skalak model is given by krueger11

[TABLE]

where $I_{1}$ and $I_{2}$ stand for the invariants of the right Cauchy-Green deformation tensor green60 ; zhu14 , and $C=\kappa_{\mathrm{A}}/\kappa_{\mathrm{S}}$ is the Skalak coefficient representing the ratio between the area dilatation modulus $\kappa_{\mathrm{A}}$ and shear modulus $\kappa_{\mathrm{S}}$ skalak73 . For $C=1$, the Skalak model is equivalent to the classical Neo-Hookean model for small membrane deformations lac04 .

Accordingly, the linearized traction jump due to shear is expressed in terms of the deformation field $\bm{u}$, and can be split into an axisymmetric and an asymmetric part as

[TABLE]

where

[TABLE]

and

[TABLE]

for the axisymmetric and asymmetric parts, respectively. Here, the asymmetric part includes all terms that depend on $u_{\phi}$ or involve derivatives with respect to $\phi$. Moreover, we have defined

[TABLE]

In addition, we introduce a resistance toward bending following the Helfrich model helfrich73 ; berndl90 ; seifert97 . The areal bending energy density thus is described by a curvature-elastic continuum model of a quadratic form given by GuckenbergerJPCM

[TABLE]

wherein $\kappa_{\mathrm{B}}$ denotes the bending modulus, $H_{0}$ stands for the spontaneous curvature (here taken as the corresponding value for the initial undeformed sphere), and $H:=b_{\alpha}^{\alpha}/2$ (summing over repeated indices) is the mean curvature, with $b_{\alpha}^{\beta}$ being the corresponding component of the curvature tensor kobayashi63 .

The traction jump equation across the membrane as derived from this model reads GuckenbergerJPCM

[TABLE]

where $\bm{n}$ is the outward-pointing unit normal vector to the spherical cavity, $K:=\mathrm{det\leavevmode\nobreak\ }b_{\alpha}^{\beta}$ stands for the Gaussian curvature, and $\Delta_{\parallel}$ denotes the Laplace-Beltrami operator deserno15 . Accordingly, bending introduces a traction jump along the normal direction which can be split into an axisymmetric and an asymmetric part as

[TABLE]

where

[TABLE]

and where we have used the shorthand notation $T:=\cot\theta$. We note that bending as derived from Helfrich’s model does not introduce discontinuities along the tangential directions. Accordingly, $\Delta f_{\theta}^{\mathrm{B}}=\Delta f_{\phi}^{\mathrm{B}}=\bm{0}$. These traction jumps reduce to the axisymmetric case considered in part I daddi18cavity for which $\left.\Delta\bm{f}^{\mathrm{S}}\right|_{\mathrm{Asy}}=\left.\Delta\bm{f}^{\mathrm{B}}\right|_{\mathrm{Asy}}=\bm{0}$. In this situation, $u_{\phi}=0$ and all derivatives with respect to $\phi$ drop out.

A closure of the problem is achieved by requiring a no-slip boundary condition at the undisplaced membrane. Accordingly, the velocity field at $r=a$ is assumed to be equal to that of the displaced material points of the elastic cavity, i.e.,

[TABLE]

which can be written in Fourier space as

[TABLE]

Our resolution methodology proceeds through writing the solution of the elastohydrodynamic problem inside the cavity as

[TABLE]

where $\bm{v}^{\mathrm{S}}=\bm{\mathcal{G}}\left(\bm{x}-\bm{x}_{2}\right)\cdot\bm{F}$ represents the velocity field induced by a point-force singularity acting at position $\bm{x}_{2}$ in an unbounded fluid — i.e., in the absence of the cavity — and $\bm{v}^{*}$ is the complementary term that is required to satisfy the imposed boundary conditions at the cavity. This type of complementary solution is often termed as the image system solution or sometimes known under the name of reflected flow field blake71 ; blake74 .

We now briefly outline the main steps in our resolution approach. First, we express the Stokeslet solution in terms of harmonics, which are then rewritten in terms of harmonics relative to the origin via the Legendre expansion zill11 . Second, the reflected flow field and the solution outside the cavity are expressed using Lamb’s general solution lamb32 with interior and exterior harmonics, respectively. This gives us a complete solution form involving a set of unknown series coefficients. These coefficients are determined from the underlying boundary conditions imposed at the cavity surface. Finally, the solution of the flow problem can then be employed to assess the effect of the confining cavity on the motion of the encapsulated spherical particle.

II.1 Stokeslet representation

For the remainder of this paper, we will scale all the lengths by the cavity radius $a$. In analogy with part I, we begin by writing the Stokeslet singularity located at position $\bm{x}_{2}$ as

[TABLE]

where we have defined $\bm{s}:=\bm{x}-\bm{x}_{2}$ and $s:=|\bm{s}|$. Here, $\bm{1}$ is the unit tensor. Using Legendre expansion, the harmonics located at $\bm{x}_{2}$ can conveniently be expressed in terms of harmonics centered at $\bm{x}_{1}$ as

[TABLE]

Here, $\varphi_{n}$ are harmonics of degree $n$, which are related to Legendre polynomials by abramowitz72

[TABLE]

where $\bm{d}:=(\bm{x}_{1}-\bm{x}_{2})/R$ is a unit vector, $\bm{r}=\bm{x}-\bm{x}_{1}$ is the position vector in the spherical coordinate system centered at the cavity center, and $r:=|\bm{r}|$. The dyadic product in eq. (15) can be written as

[TABLE]

with $\bm{\nabla}_{2}:={\partial}/{\partial\bm{x}_{2}}$. By making use of eq. (16), the derivatives with respect to $\bm{x}_{2}$ can readily be taken care of by noting that

[TABLE]

In the present work, we focus our attention on the asymmetric situation in which the force is purely tangent to the membrane surface and thus $\bm{F}\cdot\bm{d}=0$. By taking this into consideration, the Stokeslet stated in eq. (15) can therefore be expressed as

[TABLE]

Accordingly, the Stokeslet solution has now been expressed in terms of spherical harmonics positioned at the origin. By defining $\bm{t}=\bm{F}\times\bm{d}$, we have the recurrence relation

[TABLE]

In addition, imposing $\bm{F}\cdot\bm{d}=0$ yields

[TABLE]

where the harmonics $\psi_{n}$ and $\bm{\gamma}_{n}$ are, respectively, defined as

[TABLE]

These are related to each other via $\psi_{n}=\bm{\gamma}_{n}\cdot\bm{d}$.

In the following, the functions $\bm{\nabla}\psi_{n}$, $\bm{r}\psi_{n}$, and $\bm{\gamma}_{n}$ are chosen as vector basis functions to be used for expanding the velocity and pressure fields. Accordingly, the Stokeslet solution can be written in a final form as

[TABLE]

We next proceed to deriving analogous expansions for the flow fields inside and outside the spherical cavity.

II.2 The image system solution

The solution for the flow field in a spherical domain possesses a generic form known as Lamb’s general solution lamb32 ; kim13 . It involves three sets of unknown coefficients to be determined from the underlying boundary conditions, and can be written for an asymmetric situation as

[TABLE]

where we have defined

[TABLE]

Here, $a_{n}$, $b_{n}$, and $c_{n}$ are free parameters that will be determined from the boundary conditions. It is worth noting that the present solution involves three unknown coefficients for each $n$, while the simpler axisymmetric motion considered in part I only involves two sets of coefficients. Unfortunately, this also means that we are not able to to proceed as for the axisymmetric case, but have to derive the solutions using a notably different framework.

II.3 The exterior solution

The solution on the outside of the spherical cavity can be expressed in terms of exterior harmonics using Lamb’s general solution as

[TABLE]

The latter expression can be deduced from the solution for the inner fluid given by eq. (23) by making use of the substitution $n\leftarrow-(n+1)$.

The six unknown coefficients ($a_{n}$, $b_{n}$, and $c_{n}$ for the image system solution, and $A_{n}$, $B_{n}$, and $C_{n}$ for the exterior flow) can now be determined from the underlying boundary conditions of continuity of the flow velocity field and discontinuity of the hydrodynamic stress tensor across the membrane.

II.4 Velocity projections

Before proceeding with the determination of the unknown series coefficients, it is convenient to state explicitly the projected expressions of the velocity field along the radial and tangential directions.

II.4.1 Radial velocities

The radial projection of the three vector basis functions are given by

[TABLE]

In addition to that, since $\bm{e}_{r}$ and $\bm{r}$ are collinear, the scalar triple product $\bm{e}_{r}\cdot(\bm{t}\times\bm{r})\,\varphi_{n-1}$ vanishes. Moreover, the projection of eq. (20) onto the radial direction yields

[TABLE]

By making use of eqs. (25) and (26) in the radial projection of eqs. (22), (24), and (23), the components of the fluid velocity fields along the radial direction can thus be expressed in terms of the harmonics $\psi_{n}$ as

[TABLE]

II.4.2 Tangential velocities

As for the tangential direction, we define the projection operator $\bm{\Pi}:=\bm{1}-\bm{e}_{r}\bm{e}_{r}$, which projects vectors on a plane tangent to the surface of the spherical cavity. By applying the projection operator to eq. (20), we readily obtain

[TABLE]

where we have defined the vector harmonics

[TABLE]

Additionally, the tangential projection of $\left(\bm{t}\times\bm{r}\right)\varphi_{n}$ can be taken care of by noting that

[TABLE]

Applying the projection relations stated by eqs. (28) and (29) to eqs. (22), (24), and (23), we finally obtain

[TABLE]

II.5 Determination of the series coefficients

To determine the unknown coefficients, we have to make recourse to the orthogonality properties of spherical harmonics edmonds96 . For this purpose, we introduce the following notation to describe the average of a given quantity $Q(\phi,\theta)$ over the surface of a sphere. Specifically, this means

[TABLE]

At the surface of the cavity, i.e., for $r=1$, the harmonics $\varphi_{n}$ and $\psi_{n}$ satisfy the orthogonality relations

[TABLE]

where $\delta_{mn}$ denotes the Kronecker symbol, i.e., the above terms vanish for $m\neq n$. Moreover, the vector harmonics $\bm{\Psi}_{n-1}$ and $\bm{\Gamma}_{n}$ satisfy at $r=1$ the orthogonality properties

[TABLE]

We further note that their derivatives with respect to $r$ (needed in the calculation of the stress jumps) satisfy the recurrence relations

[TABLE]

II.5.1 Pressure field

Knowing the velocity fields on both sides of the elastic cavity, the inner and outer pressure fields can readily be calculated from the fluid motion equations. The solution inside the spherical cavity, which comprises both contributions from the Stokeslet and the image system solution, can be expressed in terms of a multipole expansion as

[TABLE]

Outside the cavity, only the exterior harmonics that decay at larger distances should be accounted for, thus excluding contributions of the form $r^{2n+1}\psi_{n-1}$. After some algebra, we obtain

[TABLE]

II.5.2 Continuity of velocity

The projections of the fluid velocity field along the radial and tangential directions can be presented in a generic form as

[TABLE]

wherein $q=i$ holds for the fluid on the inside, and $q=o$ for the fluid on the outside. Moreover, $\rho_{n}^{(q)}$, $\alpha_{n}^{(q)}$, and $\beta_{n}^{(q)}$, for $q\in\{i,o\}$, are radially symmetric series functions that can readily be obtained by identification with eqs. (27) and (30) giving the radial and tangential velocities, respectively.

The unknown coefficients inside the cavity can conveniently be expressed in terms of those outside thanks to the natural continuity of the velocity field across the membrane. By making use of the orthogonality properties of the basis functions, we obtain

[TABLE]

II.5.3 Discontinuity of stress tensor

Expressions for the unknown coefficients $A_{n}$, $B_{n}$, and $C_{n}$ associated with the outer fluid can be obtained from the traction jump equations across the membrane. For the sake of clarity, and to make the calculations traceable, we will consider in the following the effects of shear and bending deformation modes separately.

Pure shear

The tangential traction jump equations due to shear can conveniently be cast in the form

[TABLE]

where we have defined

[TABLE]

Here, $\bm{F}_{n}$, $\bm{G}_{n}$, and $\bm{f}_{n}$ are known series vectors, the expressions of which can be obtained by identification with eq. (4b) upon substitution of the tangential velocity field from eq. (30). They satisfy the orthogonality relations

[TABLE]

with the basis vector harmonics $\bm{\Psi}_{m-1}$, and

[TABLE]

with $\bm{\Gamma}_{n}$, where we have defined

[TABLE]

with

[TABLE]

being the shear number. Combining these equations with the orthogonality relations given by eqs. (33) yields

[TABLE]

Using our representation, the normal traction jump due to shear reads

[TABLE]

which, upon using the orthogonality property of $\psi_{n-1}$, yields

[TABLE]

By combining eqs. (36), (39), and (41), the unknown series coefficients for the outer fluid can be obtained and cast in the form

[TABLE]

where $K_{1},\dots,K_{9}$ are rather complex functions of $\alpha,\lambda$ and $n$, the expressions of which are explicitly provided in the Appendix. In the limit $i\alpha\to\infty$, which physically corresponds to a cavity membrane with an infinite shear elasticity modulus (or equivalently to a vanishing actuation frequency), the expressions of the series coefficients inside the cavity reduce to

[TABLE]

In this limit, $A_{n}$, $B_{n}$, and $C_{n}$ vanish except for $n=1$, where $(A_{1},B_{1},C_{1})=(4,2/3,-R)$. It is worthwhile to note that the coefficients given by eqs. (43) correspond to the solution for an asymmetric point force acting inside a rigid cavity with no-slip boundary conditions.

Pure bending

We now use a similar resolution procedure to determine the unknown series coefficients for a cavity membrane with pure bending resistance, such as that of a fluid vesicle or a liposome used as a vehicle for pharmaceutical drugs rui98 ; torchilin05 ; zylberberg16 . Since the tangential components of the traction are continuous, we obtain

[TABLE]

which, after applying the orthogonality properties given by eqs. (33), leads to

[TABLE]

The normal traction jump due to bending as derived from the Helfrich model reads

[TABLE]

which, upon using the orthogonality relation

[TABLE]

leads to

[TABLE]

wherein

[TABLE]

denotes the bending number.

By combining eqs. (45) and (48) with eqs. (36), the unknown series coefficients for the fluid on the outside can be cast in the form

[TABLE]

where $Q_{1},\dots,Q_{7}$ are complicated functions of $\alpha_{\mathrm{B}}$, $\lambda$, and $n$ which are given in the Appendix. In the limit $i\alpha_{\mathrm{B}}\to\infty$, corresponding to an infinite membrane bending modulus, or to a vanishing forcing frequency, the series coefficients are given by

[TABLE]

for the inner fluid, and

[TABLE]

for the outer fluid when $n\geq 2$. In addition, $(A_{1},B_{1},C_{1})=(4,2R^{2}/15,-R)$.

Combined shear and bending

An analogous resolution strategy can be adopted for the determination of the sum coefficients when the membrane is simultaneously endowed with both a resistance toward shear and bending. Analytical expressions of the coefficients can readily be obtained using computer algebra systems but these are not provided here due to their complexity and lengthiness. It is noteworthy that, in contrast to planar elastic membranes, a coupling between shear and bending deformation modes has been observed for curved membranes.

III Hydrodynamic mobility

The calculation of the flow field presented in the previous section can be utilized to assess the effect of the confining cavity on the motion of the encapsulated particle. This effect is quantified by the hydrodynamic self-mobility function $\mu$, which relates the translational velocity of a colloidal particle to the force exerted on its surface.

We now assume an arbitrary time-dependent external force $\bm{F}_{2}$ to be acting on the spherical particle positioned at $\bm{x}_{2}$. The zeroth-order solution for the translational velocity of the solid particle can readily be obtained from the Stokeslet solution as $\bm{V}_{2}^{(0)}=\mu_{0}\bm{F}_{2}$, where $\mu_{0}=1/(6\pi\eta b)$ is the usual Stokes mobility for a sphere moving in an unconfined viscous fluid. The leading-order correction to the hydrodynamic self-mobility can be calculated from the reflected flow field as

[TABLE]

The latter result is often denominated as the mobility correction in the point-particle approximation bickel06 ; bickel07 . Higher-order correction terms can be obtained by employing a combination of the multipole expansion and the Faxén theorem swan07 ; swan10 . However, we will show in the sequel that this approximation, despite its simplicity, can surprisingly lead to a good prediction of the mobility correction when comparing with fully-resolved computer simulations.

By making use of the relations

[TABLE]

in addition to inserting eq. (23) into eq. (53), we write the scaled mobility correction as

[TABLE]

where we have used the relation $P_{n}^{\prime}(1)=n(n+1)/2$ for the derivative at the end point. We further remark that $R\in[0,1)$ because all distances have been scaled by the cavity radius $a$. The general term in the latter series, which we denote by $f_{n}(\alpha,R)$, has an asymptotic behavior at infinity that does not depend on the shear and bending properties of the membrane. Specifically, we obtain as $n\to\infty$

[TABLE]

In particular, for $R=0$, the mobility correction simplifies to

[TABLE]

in full agreement with the result obtained in part I for a particle concentric with the elastic cavity. We recall that the shear number $\alpha$ has previously been defined by eq. (38), and the dimensionless parameter $\lambda$ associated with the Skalak ratio has been defined by eq. (8a).

In the quasi-steady limit of vanishing frequency, the scaled correction to the mobility reads

[TABLE]

wherein $\Delta\mu^{\mathrm{R}}/\mu_{0}$ is the scaled correction to the particle mobility associated with asymmetric motion inside a rigid spherical cavity. This correction can readily be obtained by substituting the series coefficients given by eq. (43) into eq. (55) to obtain

[TABLE]

in agreement with the results by Aponte-Rivera and Zia aponte16 ; aponte18 ; aponte-thesis , who provided the elements of the grand mobility tensor for general motion inside a rigid cavity. Interestingly, the particle mobility in the limit of infinite stiffness is found to be always larger than that inside a rigid cavity with no-slip velocity boundary condition on its interior surface. Mathematically, this behavior can be justified by the fact that the limit and sum operators cannot generally be swapped in every situation. In fact, using Fatou’s Lemma carothers00 , it can be shown that

[TABLE]

That is, evaluating the sum over $n$ before taking the limit $\alpha\to\infty$ (as for an elastic cavity) could lead, under some circumstances, to a larger value in magnitude compared to the case in which the sum is taken after taking the limit (as it is the case for a rigid cavity). This is explained by the fact that the dominated convergence theorem does not apply for the series function at hand billingsley13 .

We further mention that the same limit given by eq. (58) is obtained when the cavity membrane only possesses resistance toward shear. In the limit of infinite cavity radius, the classic result for motion parallel to a planar hard wall is recovered. Specifically,

[TABLE]

wherein $h=1-R$ denotes the distance between the center of the particle and the closest point of the cavity membrane.

Next, we consider an idealized cavity membrane with pure bending resistance and calculate the correction to the self-mobility function in the limit of $\alpha_{\mathrm{B}}\to\infty$, corresponding to an infinite bending modulus or to particle motion in the quasi-steady limit of vanishing frequency. After some algebra, we obtain

[TABLE]

wherein $\Delta\mu_{\mathrm{D}}/\mu_{0}$ is the scaled correction to the particle mobility for motion inside a spherical drop of infinite surface tension (with vanishing normal velocity on its surface), given by

[TABLE]

Again, the particle mobility in the vanishing-frequency limit for a membrane with pure bending is found to be always larger than that inside a spherical drop. Notably, the mobility correction vanishes in the concentric configuration corresponding to $R=0$ where the system behavior is solely determined by membrane shear properties. This is in agreement with the results of part I obtained by exactly solving the fluid motion equations for an extended particle of finite size concentric with an elastic cavity.

In the limit of infinite cavity radius, we recover the mobility correction near a planar fluid-fluid interface,

[TABLE]

in agreement with the result by Lee and Leal lee79 .

In the following, we assess the appropriateness and validity of our analytical calculations by direct comparison with computer simulations based on a completed-double-layer boundary integral method pozrikidis01 . The method is perfectly suited for solving numerically diverse flow problems in the Stokes regime involving both rigid and elastic boundaries. For technical details regarding the computational method and its numerical implementation, we refer the reader to refs. daddi16b, and guckenberger16, .

To probe the effect of the confining elastic cavity on the motion of an encapsulated particle, we present in fig. 2 the variations of the scaled correction to the self-mobility as a function of the forcing frequency, for a cavity membrane possessing only shear (green), only bending (red), or both shear and bending deformation modes (black). Here, the particle of radius $b=1/10$ is positioned at $R=4/5$ from the cavity center. We observe that the real (reactive) part of the mobility correction (shown as dashed lines) is a monotonically increasing function with frequency and approaches zero for larger forcing frequencies. In contrast to that, the imaginary (dissipative) part (shown as solid lines) exhibits the typical bell-shaped profile which peaks at around $\omega\sim\kappa_{\mathrm{S}}/(\eta a)$. In the low-frequency regime, the mobility correction approaches the plateau values predicted by eqs. (58) and (62) for a cavity membrane with only shear elasticity or pure bending, respectively. Overall, there is strong quantitative agreement between the full numerical solutions (symbols) and the theoretical predictions. The small observed discrepancy notably for the real part in the low-frequency regime is most probably due to the finite size effect, because the analytical predictions are based on the point-particle approximation, whereas the numerical simulations necessarily account for the finite radius of the solid particle.

IV Cavity motion and membrane deformation

IV.1 Pair (composite) mobility

The hydrodynamic self-mobility discussed in sect. III represents the particle response function to an external force. In this regard, one can also define an analogous response function for the whole elastic cavity and its interior, that relates the translational velocity $\bm{V}_{1}$ of the cavity centroid to the force $\bm{F}_{2}$ exerted on the encapsulated particle via $\bm{V}_{1}=\bm{\mu}^{\mathrm{P}}\cdot\bm{F}_{2}$. In accordance to part I, we call the tensor $\bm{\mu}^{\mathrm{P}}$ the pair (composite) mobility. By symmetry, $\bm{V}_{1}\parallel\bm{F}_{2}$ holds, so that the components of $\bm{\mu}^{\mathrm{P}}$ reduce to a single entry $\mu^{\mathrm{P}}$ connecting the corresponding magnitudes via $V_{1}=\mu^{\mathrm{P}}F_{2}$.

Without loss of generality, we assume in the following that $\bm{F}_{2}$ is exerted along the $x$-direction. Accordingly, the translational velocity of the elastic cavity can be calculated by integration over the fluid domain inside the cavity as felderhof14

[TABLE]

where $\Omega=4\pi/3$ is the scaled volume of the undeformed cavity, and

[TABLE]

The resulting frequency-dependent pair mobility function is obtained as

[TABLE]

so that only the term corresponding to $n=1$ remains after volume integration. Upon simplification and rearrangement, the result can be presented in a scaled form as

[TABLE]

where the parameters $\alpha$ and $\lambda$ are defined by eq. (38) and eq. (8a), respectively.

Consequently, $\mu^{\mathrm{P}}$ depends only on the membrane shear properties and can be described by a simple Debye model with a single relaxation time $\tau/\tau_{\mathrm{S}}=15/\left(2(2\lambda-1)\right)$, where $\tau_{\mathrm{S}}=a\eta/\kappa_{\mathrm{S}}$ is a characteristic time scale for shear. Remarkably, the pair mobility can also become independent of frequency for $R=\sqrt{30}/6\approx 9/10$, a value for which $6\pi\eta\mu^{\mathrm{P}}=1$. Nevertheless, as $R\sim 1$, it becomes essential to ensure that the inequality $R+b\ll 1$ remains satisfied, for the point-particle approximation employed here to be applicable.

In fig. 3, we show the variations of the pair mobility (scaled by $6\pi\eta$) as a function of the scaled frequency. Results are shown for a cavity membrane with pure shear (green), pure bending (red), and both shear and bending (black). The pair mobility for a bending-only membrane remains unchanged upon varying the frequency and amounts to $3/2-3R^{2}/5$. In contrast to that, the reactive part for a membrane possessing a shear resistance shows a logistic sigmoid curve varying between 1 (when $\alpha\to\infty$) and $3/2-3R^{2}/5$ (when $\alpha=0$), whereas the dissipative part exhibits a Gaussian-like – or – bell-shaped profile. In all cases, there is strong agreement between the series-expansion theory (solid lines) and the full numerical solutions (symbols), confirming our theoretical predictions that the pair mobility is solely dependent on membrane shear properties (and independent of bending properties).

In analogy to the above discussion of the translational motion of the cavity in the presence of a force acting on the enclosed particle, one can also consider the corresponding rotational response. The angular velocity $\bm{\Omega}$ of the cavity is (due to symmetry) of the form $\bm{\Omega}=\Omega\bm{e}_{y}$ and has to fulfill $\bm{v}^{(q)}(\bm{r})=\bm{\Omega}\times\bm{r}$ at the surface of the cavity ($r=1$), with $q\in\{i,o\}$. After some algebra, one obtains

[TABLE]

upon inserting our solution for the flow field, with angular brackets again denoting the surface average defined in eq. (31). Here, we find the same value of $\Omega$ for the different series coefficients obtained for pure shear, pure bending, as well as combined shear and bending. Additionally, we note that only the term $\propto\bm{\Gamma}_{0}$ in eqs. (35) contributes to the rotation of the cavity, while all other terms lead to vanishing contributions. Accordingly, the angular velocity effectively stems only from the Stokeslet solution and does not depend on membrane shear and bending properties.

IV.2 Membrane deformation

The elastic deformation of the membrane can be assessed by calculating the displacement field $\bm{u}(\theta,\phi,\omega)$ resulting from the external force acting on the particle. This field quantifies the motion of the material points of the cavity membrane relative to their initial positions in the undeformed state. In the small deformation regime, the displacement field can readily be obtained from the no-slip boundary condition given by eq. (13), to obtain

[TABLE]

We now define the reaction tensor $\bm{R}$, a frequency-dependent tensorial quantity relating the membrane displacement field of the cavity to the asymmetric point force as bickel07

[TABLE]

By considering a harmonic oscillation of the form $\bm{F}=\bm{K}e^{i\omega_{0}t}$, of amplitude $\bm{K}$ and frequency $\omega_{0}$, the membrane displacement in real space can readily be obtained from inverse Fourier transform as bracewell99

[TABLE]

An exemplary displacement field is displayed in fig. 4 as a function of the polar angle for three different forcing frequencies. The azimuthal angle $\phi$ is chosen to represent the planes of maximum deformation for the respective components, as described in the figure caption. Here, the cavity membrane is endowed with both shear and bending rigidities. We observe that the radial component $u_{r}$ vanishes at the upper pole and shows a peak around $\theta/\pi\approx 1/8$, before decaying quasi linearly to zero upon increasing $\theta$. The in-plane displacements $u_{\theta}$ and $u_{\phi}$ display a maximum value at the upper pole, and monotonically decay as $\theta$ increases. Our analytical predictions are in good agreement with numerical simulations. Notably, we observe a small deviation in the plot for $u_{r}$ shown in panel $(a)$ which is most probably due to a finite-sized effect. In contrast to the axisymmetric case discussed in part I, the deformation here is (in general) largest in the tangential direction.

In typical biological situations, the forces that could be exerted by optical tweezers on particles are of the order of 1 pN cipparrone10 . The spherical cavity may have a radius of 10^{-6}\leavevmode\nobreak\m and a shear modulus of \kappa_{\mathrm{S}}=5\times 10^{-6}\leavevmode\nobreak\N/m Freund_2014 . For a scaled frequency $(3\eta a\omega)/(2\kappa_{\mathrm{S}})=4$, the membrane cavity is expected to undergo a maximal deformation of only about 1 % of its initial undeformed radius. Consequently, cavity deformations and deviations from the spherical shape are notably small.

V Conclusions

In summary, we have presented an analytical theory to describe the low-Reynolds-number motion of a spherical particle moving inside a spherical membrane cavity endowed with both shear elasticity and bending rigidity. Here, we have focused on the situation in which the force exerted on the particle is directed tangent to the surface of the cavity. Together with the axisymmetric results obtained in an earlier paper daddi18cavity , the solution of the elastohydrodynamic problem for a point force acting inside a spherical elastic cavity is thus obtained.

We have expressed the solution of the flow problem using the method of images. For this purpose, the hydrodynamic flow field is represented by a multipole expansion, summing over modes in terms of spherical harmonics, in analogy with familiar methods in electrostatics. The unknown series coefficients associated with each mode have been determined analytically from the prescribed boundary conditions of continuity of the fluid velocity field at the membrane cavity and discontinuity of hydrodynamic stresses as derived from Skalak and Helfrich elasticity models, associated with shear and bending deformation modes, respectively.

We have then explored the role of confinement on the motion of the encapsulated particle by calculating the frequency-dependent mobility functions. The latter linearly couple the translational velocity of the particle to the external force exerted on it. In the quasi-steady limit of vanishing actuation frequency, we have demonstrated that the hydrodynamic mobility inside a spherical elastic cavity is always larger than that predicted inside a rigid cavity of equal size with no-slip surface conditions. In addition, we have quantified the translational and rotational motion of the confining cavity, finding that the translational pair (composite) mobility is uniquely determined by membrane shear elasticity and that bending does not play a role in the dynamics of the cavity. We have further assessed the membrane deformation caused by the motion of the particle, showing that the cavity membrane primarily experiences deformation along the tangential direction.

Finally, we have assessed the appropriateness and applicability of our theoretical approach by supplementing our analytical calculations with fully-resolved computer simulations of truly-extended particles using the boundary integral method. Good agreement is obtained between theoretical predictions and numerical simulations over the full range of applied forcing frequencies. The developed method may find applications in the simulation of hydrodynamically interacting microparticles confined by a spherical elastic cavity, or medical capsules that are directed to a requested site by magnetic forces acting on incorporated magnetic particles.

Acknowledgements.

We would like to thank Stephan Gekle and the Biofluid Simulation and Modeling group at the University of Bayreuth, where the boundary integral code used in this work has been developed. The authors gratefully acknowledge support from the DFG (Deutsche Forschungsgemeinschaft) through the projects LO 418/17-2 (H.L.), ME 3571/2-2 (A.M.M.), and DA 2107/1-1 (A.D.M.I).

Author contribution statement

A.D.M.I. conceived the study and performed the numerical simulations. C.H. and A.D.M.I. carried out the analytical calculations and drafted the manuscript. C.H., H.L., A.M.M., and A.D.M.I. discussed and interpreted the results, edited the text, and finalized the manuscript.

Appendix A Expression of the coefficients

In this Appendix, we provide explicit expressions for the coefficients stated in eqs. (42) and (50) of the main body of the paper. For an idealized membrane with pure shear, the coefficients are given by

[TABLE]

where we recall that $\alpha=2\kappa_{\mathrm{S}}/(3\eta i\omega)$ is the shear number, and $\lambda=C+1$ is the dimensionless parameter associated with the Skalak ratio. For an idealized membrane with pure bending resistance, the corresponding coefficients read

[TABLE]

wherein $\alpha_{\mathrm{B}}=\kappa_{\mathrm{B}}/(\eta i\omega)$ denotes the bending number.