Axisymmetric spheroidal squirmers and self-diffusiophoretic particles

Self-propulsion of micro-organisms, such as bacteria, algae and protozoa, plays an important role in many aspects of nature. Whether a bacteria tries to reach a nutrient rich area or a sperm cell an unfertilized egg, motility often yields a substantial advantage over competitors. Due to their small size and velocity, the viscous friction experienced by micro-organisms swimming in water is very strong compared to inertial forces [1, 2]; consequently, they have developed swimming mechanisms adapted to these circumstances.

The motion of motile bacteria and other small organisms is typically induced by the beating of thin, thread-like appendages, the so called flagella; however, these species exhibit a broad morphology, in that they may posses either a single flagellum (e.g. monotrichous bacteria), several flagella (e.g. E. coli.), or a carpet of many small flagella, called cilia (e.g. Opalina). Focusing for the moment on the example of cilia covered micro-organisms, it was proposed by Lighthill that the emergence of motility can be understood without a detailed modeling of the complex, synchronous beating of the cilia [3]. Instead, the effect of this beating pattern, i.e. the time-averaged surface flow induced by the envelope of cilia tips, has been modeled as providing a prescribed active flow velocity ${\bf v}$ (actuation) at the surface of the particle (both within the tangent plane of the surface, ${\bf v}_s$, and in the direction normal to the surface, ${\bf v}_n$). The model, known by now as the 'squirmer' model, was subsequently corrected and extended by his student Blake [4].

In its simplest and most used form, the tangential slip velocity of a spherical squirmer is taken to be the superposition of the fore-aft asymmetric and fore-aft symmetric modes with the slowest decaying contributions to flow in the surrounding liquid, i.e. the leading order modes. Combined, they make for a squirmer endowed with motility (from the leading fore-aft asymmetric mode) and with a hydrodynamic stresslet (from the leading fore-aft symmetric mode) [5–9].

Considerable progress has already been made in understanding the behavior of spherical squirmers, and many interesting questions have been answered, such as: what does the flow field around a squirmer look like and how it compares with the ones produced by simple microorganisms [4, 10, 11]; what happens if the swimmer is not in free space, but rather disturbed by boundaries [5] or external flows [8]; and how do pairs or even swarms of these particles interact [6]. Although this model can also be used to better understand the motion of micro-organisms, e.g. Volvox [12], observed in experiments, its restriction to spherical swimmers limits a wider application. For example, Paramecia, one of the most studied ciliates, has an elongated body [13–15], which prevents a straightforward application of the traditional squirmer model. However, the aforementioned questions are also of interest for these and other elongated micro-organisms. A generalization of the squirmer model to non-spherical shapes is thus a natural, useful development.

Driven by advances in technology that allow increasingly sophisticated manufacturing capabilities, in the last decade significant efforts have been made towards the development of artificial swimmers [16–20]. The envisioned lab-on-a-chip devices and micro-cargo carriers [21, 22], e.g. for targeted drug deliveries and nanomachines focused on monitoring and dissolving harmful chemicals [23, 24], all need to perform precise motions on microscopic length scales. A better understanding of the general framework of microscale locomotion is required in order to optimally design and control such devices, and theoretical models facilitate new steps along this path. As one example, chemically active colloids achieve self-propulsion by harvesting local free energy. They catalyze a chemical reaction in the surrounding fluid and propel due to the ensuing chemical gradient. Spherical 'Janus' particles belong to this group, and it has been shown that the squirmer model can, in various circumstances, capture essential features of their motion [25]. However, there are many chemically active colloids with non-spherical shape, and rod-like particles are especially prevalent in experimental studies (see, e.g. [19, 26–28]).

It is well known that passive, non-spherical colloids exposed to an ambient flow exhibit significant qualitative differences—such as alignment with respect to the direction of the ambient flow [29], Jeffery orbits by ellipsoids in shear flow [30], nematic ordering arising from steric repulsion [31], or noise-induced migration away from confining surfaces [32]—in comparison to their spherical counterparts. It is thus reasonable to expect that endowing such objects with a self-propulsion mechanism will lead to rich, qualitatively novel dynamical behaviors, some of which may be advantageous for, while others may hinder, applications. Accordingly, it is important to develop an in-depth understanding of the shape-dependent behavior, and significant experimental and theoretical efforts have been made in this direction (see, e.g. [33–42] as well as the insightful reviews provided in [1, 21, 43, 44–49]).

Similar to the case of spherical swimmers, physical insight into the phenomenology exhibited by elongated swimmers can be gained from a corresponding squirmer model. For a model spheroidal swimmer moving by small deformations of its surface, [50] derived an analytic solution to the corresponding Stokes flow, from which the velocity of the swimmer could be calculated. In [51] it has been pointed out that an exact squirmer model for a spheroidal particle with a prescribed axi-symmetric, tangential slip velocity (active actuation of the fluid) on its surface can be written down by employing an available analytical solution for the axi-symmetric Stokes flow around a spheroidal object [52]. However, the approach has been used in the context of a somewhat restricted model particle, involving an additional fore-aft asymmetry of the surface slip velocity, because of the particular interest, for that work, in the question of 'hydrodynamically stealthy' microswimmers. While capturing the self-propulsion velocity of the swimmer, this removes certain characteristics of the flow, inter alia those that contribute to the corresponding stresslet [53] and would allow distinguishing between, e.g. 'puller' and 'pusher' type squirmers. Moreover, the particle stresslet is a key quantity connecting the microscopic dynamics of individual particles in a colloidal suspension with the continuum rheological properties of the suspension. For active particles like bacteria, the activity-induced stresslet can lead to novel material properties like 'superfluidity' and spontaneous flow [39].

We note that also a strongly truncated model, based on an ansatz for the slip velocity in the form of two terms which, in the corresponding limit of a sphere, reproduce the first two modes of the spherical squirmer model of Lighthill and Blake, has been discussed by [54]. This approach, however, is significantly affected by the fact that—as noticed in [50] and also discussed here (see section 4)—for spheroidal squirmers both their velocity and their stresslet involve significant contributions from the higher order slip modes, in contrast to the case of a spherical squirmer for which only the first two modes contribute to those observables.

In this paper we employ the available analytical solution for axi-symmetric Stokes flow around a spheroidal object [52] to study the velocity and the induced hydrodynamic flow field around a spheroidal squirmer with a tangential slip velocity possessing axial symmetry, but otherwise unconstrained. The model squirmer is introduced in section 2. The series representation of the incompressible, axi-symmetric, creeping flow field around a spheroid [52] is succinctly summarized in section 3. In section 4 we discuss the velocity and the flow field corresponding to the spheroidal squirmer, with particular emphasis on illustrating the contributions from the modes of various order. Additionally, in section 4.1 we discuss the straightforward extension of these results to deriving the flow field around a chemically active self-phoretic spheroid (for a similar mapping in the case of spherical particles see [55]). The final section 5 is devoted to a summary of the results and to the conclusions of the study.

The model system we consider is that of a spheroidal, rigid and impermeable particle immersed in an incompressible, unbounded, quiescent Newtonian liquid through which it moves due to a prescribed 'slip velocity' (active actuation) at its surface (see figure 1). The slip velocity ${\bf v}_s$, which is tangential to the surface $\Sigma$ of the particle (i.e. $\newcommand{\e}{{\rm e}} {\bf n} \cdot {\bf v}_s \equiv 0$ on $\Sigma$, with ${\bf n}$ denoting the outer (into the fluid) normal to $\Sigma$), is assumed to preserve the axial symmetry and to be constant in time, but it is otherwise arbitrary. The surface slip ${\bf v}_s$ is part of the model and thus it is a given function (or, alternatively as in, see, section 4.1, it is determined as the solution of a separate problem). There are no external forces or torques acting on either the particle or the liquid.

Zoom In Zoom Out Reset image size

Due to the surface actuation, a hydrodynamic flow around the particle is induced and the particle sets in motion; we assume that the linear size L of the particle, the viscosity $\newcommand{\m}[1]{\mathrm{#1}} \mu$ and the density $\newcommand{\re}{{\rm Re}} \renewcommand{\rho}{\varrho} \rho$ of the liquid, and the magnitude $|{\bf v}_s|$ of the slip velocity are such that the Reynolds number $\newcommand{\m}[1]{\mathrm{#1}} \newcommand{\re}{{\rm Re}} \renewcommand{\rho}{\varrho} \mathrm{Re} := \rho |{\bf v}_s| L/\mu$ is very small. Accordingly, after a short transient time a steady state hydrodynamic flow is induced around the particle and the particle translates steadily with velocity ${\bf U}$ (with respect to a fixed system of coordinates, the 'laboratory frame'). (Owing to the axial symmetry, in the absence of thermal fluctuations, which are neglected in this work, there is no rigid-body rotation of the particle in this model.)

The analysis is more conveniently carried out in a system of coordinates attached to the particle (co-moving system). This is chosen with the origin at the center of the particle and such that the z-axis is along the axis of symmetry of the particle (thus ${\bf U} = U {\bf e}_z$). The semi-axes of the particle are accordingly denoted by b_z and b_x (see figure 1); their ratio $r_e = \frac{b_x}{b_z}$ is called the aspect ratio and determines the slenderness of the particle. The values r_e   <  1, r_e   =  1, and r_e   >  1 correspond to prolate, spherical, and oblate shapes, respectively.

In the co-moving system of coordinates, the flow ${\bf v}$ and the velocity ${\bf U}$ are determined as the solution of the Stokes equations

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:Stokes} \nabla \cdot \boldsymbol{\sigma} = 0,~\nabla \cdot {\bf v} = 0, \nonumber \end{align} \tag{ 1 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\m}[1]{\mathrm{#1}} \displaystyle \label{eq:stress_tens} \boldsymbol{\sigma} := - p\, \mathcal{I} + \mu [\nabla {\bf v} + (\nabla {\bf v})^T], \nonumber \end{align} \tag{ 2 }$$

is the Newtonian stress tensor, with p  denoting the pressure (enforcing incompressibility), $\newcommand{\m}[1]{\mathrm{#1}} \mathcal{I}$ denoting the unit tensor, and ()^T denoting the transpose, subject to:

• the boundary conditions (BC)
  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:BC_surf} {\bf v}|_\Sigma = {\bf v}_s, \nonumber \end{align} \tag{ 3a }$$
  and
  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:BC_infty} {\bf v}(|{\bf r}| \to \infty) \rightarrow -{\bf U}. \nonumber \end{align} \tag{ 3b }$$
• the force balance condition (overdamped motion of the particle in the absence of external forces)
  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:F_bal} \int_\Sigma \boldsymbol{\sigma} \cdot {\bf n} \, {\rm d}\Sigma = 0. \nonumber \end{align} \tag{ 4 }$$

The hydrodynamic flow $\newcommand{\m}[1]{\mathrm{#1}} {\bf v}_\mathrm{lab}$ in the laboratory frame, if desired, is then obtained as $\newcommand{\m}[1]{\mathrm{#1}} {\bf v}_\mathrm{lab}(\mathrm{P}) = {\bf v}(\mathrm{P}) + {\bf U}$, with $\newcommand{\m}[1]{\mathrm{#1}} \mathrm{P}$ denoting the observation point.

The boundary-value problem above can be straightforwardly solved numerically by using, e.g. the boundary element method (BEM) (see, e.g. [8, 56–58]) as well as analytically. In the following, we shall focus on the analytical solution of this problem; the corresponding numerical solutions obtained by the BEM are presented in appendix E and used as a means of testing the convergence of the series representation of the analytical solution. Before proceeding, we note that we will focus the discussion on the case of prolate shapes, i.e. r_e   <  1 (see figure 1(a)); the case of an oblate squirmer (r_e   >  1) can be obtained from the results for the prolate shapes via a certain mapping (see, e.g. [51] and appendix A), while the case of a sphere is obtained from the results corresponding to a prolate spheroid by taking the limit $r_e \to 1^{-}$.

The calculation of the hydrodynamic flow ${\bf v}$ induced by the squirmer and of the velocity ${\bf U}$ of the squirmer are most conveniently carried out by employing the modified prolate spheroidal coordinates $\newcommand{\p}{\partial} (1\leqslant\tau\leqslant\infty,-1\leqslant\xi\leqslant1,0\leqslant\phi\leqslant2\pi)$ as in [52]. These are defined in the co-moving system of reference, which has the origin in the center of the particle, and are related to the Cartesian coordinates via [52]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:prol_coord} \tau = \frac{1}{2c}\Big(\sqrt{x^2+y^2+(z+c)^2}+ &\sqrt{x^2+y^2+(z-c)^2}\large{\Big)},\nonumber \\ \zeta= \frac{1}{2c}\Big(\sqrt{x^2+y^2+(z+c)^2}-&\sqrt{x^2+y^2+(z-c)^2}\Big),\nonumber \\ \varphi = \arctan & \left(\frac{y}{x}\right),\nonumber \end{align} \tag{ 5 }$$

where $c=\sqrt{b_z^2-b_x^2}$ is a purely geometric quantity which ensures smooth convergence into spherical coordinates in the limit $b_x \rightarrow b_z$. The unit vectors ${\bf e}_\tau$ and ${\bf e}_\zeta$ (see figure 1) are related to the ones of the Cartesian coordinates via

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:unit_vec} {\bf e}_\tau &= \left(\frac{2\tau \cdot \sqrt{1-\zeta^2}}{\sqrt{\tau^2-1}} {\bf e}_x + \zeta {\bf e}_z \right) \cdot \frac{\sqrt{\tau^2-1}}{\sqrt{\tau^2-\zeta^2}}, \nonumber \\ {\bf e}_\zeta &= \left (\frac{2\zeta \cdot \sqrt{\tau^2-1}} {\sqrt{1-\zeta^2}} {\bf e}_x + \tau {\bf e}_z \right) \cdot \frac{\sqrt{1-\zeta^2}}{\sqrt{\tau^2-\zeta^2}}, \nonumber \end{align} \tag{ 6 }$$

and the corresponding Lamé metric coefficients are given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:metric_coef} & h_\zeta =c\frac{\sqrt{\tau^2-\zeta^2}}{\sqrt{1-\zeta^2}}, \nonumber \\ & h_\tau=c\frac{\sqrt{\tau^2-\zeta^2}}{\sqrt{\tau^2-1}}, \nonumber \\ & h_\varphi = c \sqrt{\tau^2-1}\sqrt{1-\zeta^2}.\nonumber \end{align} \tag{ 7 }$$

The iso-surfaces $\tau =$ const are confocal prolate spheroids with common center $\newcommand{\m}[1]{\mathrm{#1}} \mathrm{O}$; the surface of the particle corresponds to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:tau_part} \tau_0 = \frac{b_z}{c} = \frac{1}{\sqrt{1-r_e^2}} > 1, \nonumber \end{align} \tag{ 8 }$$

and the values $\tau > \tau_0$ ($1 \leqslant \tau < \tau_0$) correspond to the exterior (interior) of the prolate ellipsoid, respectively. The coordinate $\zeta$ takes the values $\newcommand{\p}{\partial} \zeta = \pm 1$ at the $\newcommand{\p}{\partial} z = \pm b_z$ apexes, respectively, and $\zeta = 0$ on the equatorial (x  −  y  plane cut) circle. Note that, as shown in figure 1(a), ${\bf e}_\tau$ coincides with the normal ${\bf n}$ while ${\bf e}_\zeta$ and $\newcommand{\p}{\partial} {\bf e}_\phi$ span the tangent plane.

3.1. Velocity of the squirmer

The velocity of the squirmer can be determined from equations (1)–(4) as a linear functional of the axi-symmetric slip velocity ${\bf v}_s = v_s(\zeta) {\bf e}_\zeta$ without explicitly solving for the flow ${\bf v}$. This follows via a straightforward application of the Lorentz reciprocal theorem [59, chapter 3–5]. For the case of the prolate spheroid with ${\bf U}=U{\bf e}_z$, this renders (for details see, e.g. [51, 53, 60])

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle U = -\frac{\tau_0}{2}\int^1_{-1} {\rm d}\zeta \, \sqrt{\frac{1-\zeta^2}{\tau_0^2-\zeta^2}} \, v_s(\zeta). \label{eq:Ulorentz} \nonumber \end{align} \tag{ 9 }$$

Therefore, ${\bf U}$ can be considered as known; accordingly, the boundary value problem defined by equations (1)–(4) is specified and the solution ${\bf v}$ can be determined.

3.2. Stokes stream function and hydrodynamic flow

Since the boundary value problem defined by equations (1)–(3) has axial symmetry, one searches for an axisymmetric solution ${\bf v}({\bf r})$ expressed in terms of a Stokes stream function $\newcommand{\p}{\partial} \psi({\bf r})$ as $\newcommand{\p}{\partial} {\bf v}({\bf r}) = \nabla \times \frac{\psi({\bf r})}{h_\varphi}{\bf e}_\varphi$ [59, chapter 4]; this renders

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\p}{\partial} \displaystyle v_\tau(\tau,\zeta) = \frac{1}{h_\zeta h_\varphi}\frac{\partial \psi}{\partial \zeta},\label{velos} \nonumber \end{align} \tag{ 10 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\p}{\partial} \displaystyle v_\zeta(\tau,\zeta) = - \frac{1}{h_\tau h_\varphi}\frac{\partial \psi}{\partial \tau}. \label{velo1} \nonumber \end{align} \tag{ 11 }$$

The general 'semiseparable' solution for the stream function in prolate coordinates has been derived by [52] in the form of a series representation

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\p}{\partial} \displaystyle \psi(\tau,\zeta)& = g_0(\tau)G_0(\zeta)+g_1(\tau)G_1(\zeta) \nonumber \\ &\quad +\sum_{n=2}^{\infty}[g_n(\tau)G_n(\zeta)+h_n(\tau)H_n(\zeta)] \label{Ansatz1} \nonumber \end{align} \tag{ 12 }$$

where G_n and H_n are the Gegenbauer functions of the first and second kind, respectively (see [61]). The functions $g_n(\tau)$ and $h_n(\tau)$ are given by certain linear combinations of $G_k(\tau)$ and $H_k(\tau)$, the coefficients of which are fixed by the corresponding boundary conditions (see below).

The general solution above is applied to our particular system as follows. Noting that for our system the solution ${\bf v}({\bf r})$ should be bounded (i.e. $|{\bf v}({\bf r})| < \infty$ everywhere), none of the terms involving the Gegenbauer functions of the second kind $H_n(\zeta)$, which are divergent along the z-axis (i.e. $\newcommand{\p}{\partial} \zeta = \pm 1$) [59, chapter 4–23], can be present in the solution; accordingly, $\newcommand{\e}{{\rm e}} h_n(\tau) \equiv 0$ for all $n \geqslant 2$. Furthermore, due to the same requirement of bounded magnitude of the flow, the terms involving the Gegenbauer functions of the first kind $G_0(\zeta) = 1$ and $G_1(\zeta) = - \zeta$ also cannot be present in the solution because they lead to divergences of $v_\zeta$ at $\newcommand{\p}{\partial} \zeta = \pm 1$ (see equations (11) and (7)); accordingly, $\newcommand{\e}{{\rm e}} g_0(\tau) \equiv 0$ and $\newcommand{\e}{{\rm e}} g_1(\tau) \equiv 0$. Therefore, for our system only the functions $g_{n \geqslant 2}$ will be of interest; these are given by [52]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle g_2(\tau)=&\frac{C_2}{6}G_1(\tau)+C_4H_4(\tau)+D_2H_2(\tau)\nonumber \\ &+F_2 G_2(\tau)+E_4 G_4(\tau)\nonumber \\ g_3(\tau)=&-\frac{C_3}{90}G_0(\tau)+C_5H_5(\tau)+D_3H_3(\tau)\label{coefficients} \nonumber \\ &+F_3 G_3(\tau)+E_5 G_5(\tau) \nonumber \\ g_{n\geqslant 4 }(\tau)=& C_{n+2}H_{n+2}(\tau) + C_{n}H_{n-2}(\tau)+D_nH_n(\tau)\nonumber \\ &+F_n G_n(\tau)+E_{n+2} G_{n+2}(\tau)+E_{n} G_{n-2}(\tau),\nonumber \end{align} \tag{ 13 }$$

where the constants $\{C_n, D_n\}_{n \geqslant 2}$, $\{E_n\}_{n \geqslant 4}$, and $\{F_n\}_{n \geqslant 2}$, are fixed, as noted above, by requiring that the solution satisfies the BCs, equation (3), and the force balance condition, equation (4).

Since the hydrodynamic force is proportional to the coefficient C₂ (see $g_2(\tau)$) [51], the force balance (equation (4)) implies

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:forc_bal_psi} C_2 = 0. \nonumber \end{align} \tag{ 14a }$$

The boundary condition at infinity (equation (3b )) implies the asymptotic behavior of the stream function

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\p}{\partial} \displaystyle \psi(\tau\to \infty,\zeta) \propto \frac{1}{2} U c^2 (\tau^2-1)(1-\zeta^2); \nonumber \end{align*}$$

by noting the asymptotic behaviors $\newcommand{\p}{\partial} G_n(\tau \gg 1) \propto \tau^n$ and $\newcommand{\p}{\partial} H_n(\tau \gg 1) \propto \tau^{-(n+1)}$ [61], this implies that the functions $g_{n\geqslant 2}(\tau)$ cannot have contributions from terms involving the Gegenbauer functions $G_n(\tau)$ of index n  >  2, i.e.

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\m}[1]{\mathrm{#1}} \displaystyle \label{eq:far_field_psi_1} E_n = 0,~n = 4,5,\dots~\mathrm{and}~F_n = 0,n = 3,4,\dots \nonumber \end{align} \tag{ 14b }$$

Furthermore, in order to exactly match the asymptotic $(\tau^2-1)(1-\zeta^2)$ form and its prefactor, it is necessary that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle F_2 = -2 c^2 U. \label{eq:velo_bc} \nonumber \end{align} \tag{ 14c }$$

The impenetrability of the surface and equation (10) imply that $\newcommand{\p}{\partial} \psi(\tau_0,\zeta)$ is a constant; since $\newcommand{\p}{\partial} \psi(\tau_0,\pm 1) = 0$ as $\newcommand{\p}{\partial} G_n(\zeta=\pm 1)=0$, one concludes that $\newcommand{\p}{\partial} \newcommand{\e}{{\rm e}} \psi(\tau_0,\zeta) \equiv 0$ and thus

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:BC_norm_psi} g_n(\tau_0) = 0,~n = 2,3,\dots \nonumber \end{align} \tag{ 14d }$$

Finally, the tangential slip velocity condition at the surface of the particle (equation (3a )) together with: the expression for $v_\zeta$ in equation (11), the orthogonality of the Gegenbauer functions $G_{n\geqslant 2}(\zeta)$ (see appendix B), and the relation between the Gegenbauer functions $G_{n\geqslant 2}(\zeta)$ and the associated Legendre polynomials $P_{n\geqslant 1}^1(\zeta)$ (see the appendix B) leads to

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\m}[1]{\mathrm{#1}} \displaystyle \label{eq:BC_tang_psi} \left.\frac{{\rm d} g_n}{{\rm d}\tau}\right|_{\tau = \tau_0} =& c^2(n-\frac{1}{2}) \int_{-1}^{+1} \mathrm{d} \zeta (\tau_0^2-\zeta^2)^{1/2} v_s(\zeta) P_{n-1}^1(\zeta),\nonumber \\ &~n = 2,3,\dots. \nonumber \end{align} \tag{ 14e }$$

For given $v_s(\zeta)$, and with U determined by equation (9), the relations (14d ) and (14e ) provide a system of linear equations determining all the remaining unknown coefficients $\{C_{n \geqslant 3},D_{n \geqslant 2}\}$, as detailed in the appendix C.

The form of equation (14e ) suggests (for the expansion of the slip velocity ${\bf v}_s$) the functions

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:def_Vn} V_n(\zeta) := (\tau_0^2-\zeta^2)^{-1/2} P_n^1(\zeta), \nonumber \end{align} \tag{ 15 }$$

as a suitable basis over the space of square integrable functions $f(\zeta)$ satisfying $\newcommand{\p}{\partial} f(\zeta = \pm 1) = 0$ (note that $V_n$ has a parametric dependence on $\tau_0 > 1$). Defining, in this space, the weighted scalar product

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\m}[1]{\mathrm{#1}} \displaystyle \label{eq:scalar_prod} \langle\,f_1(\zeta), f_2(\zeta) \rangle_{w(\zeta)} := \int_{-1}^{+1} \mathrm{d} \zeta w(\zeta)\,f_1(\zeta)\,f_2(\zeta), \nonumber \end{align} \tag{ 16 }$$

and choosing the weight as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:def_weight} w(\zeta) = \tau_0^2-\zeta^2 > 0, \nonumber \end{align} \tag{ 17 }$$

one infers (form the known properties of the associated Legendre polynomials) that indeed the set $\{V_n\}_{n \geqslant 1}$ is an orthogonal basis,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:orthog} \langle V_n(\zeta), V_m(\zeta) \rangle_{w(\zeta)} = \frac{n (n+1)}{n+1/2} \delta_{n,m}. \nonumber \end{align} \tag{ 18 }$$

Accordingly, the slip velocity function ${\bf v}_s = v_s(\zeta) {\bf e}_\zeta$ has a unique representation in terms of a series in the functions $\{V_n\}_{n \geqslant 1}$,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{squirmmodel} v_s(\zeta) = \tau_0 \sum\limits_{n \geqslant 1} B_n V_n(\zeta). \nonumber \end{align} \tag{ 19 }$$

The coefficients of the expansion are written in the form above to ensure that in the limit of a sphere ($b_x \to b_z$, i.e. $\tau_0 \to \infty$) one arrives at the usual form employed for a spherical squirmer, i.e. that of an expansion in associated Legendre polynomials $\newcommand{\re}{{\rm Re}} \renewcommand{\theta}{\vartheta} P_n^1(\cos\theta)$ (see, e.g. [4, 62]). This can be seen by noticing that in the limit $\tau_0 \to \infty$ one has $\tau_0 (\tau_0^2-\zeta^2){}^{-1/2} P_n^1(\zeta) \to P_n^1(\zeta)$ and that $b_x \to b_z$ implies $\newcommand{\re}{{\rm Re}} \renewcommand{\theta}{\vartheta} \zeta \to \cos(\theta)$; accordingly, it follows that, in the limit of the shape approaching that of a sphere, $\newcommand{\re}{{\rm Re}} \renewcommand{\theta}{\vartheta} \tau_0 V_n(\zeta) {\bf e}_\zeta \to - P_n^1(\cos\theta) {\bf e}_\theta$, and the expansion in [4] is matched identically upon changing $B_n \rightarrow \frac{2}{n(n+1)} B_n$.

With the representation of the slip velocity in terms of the function $V_n$, equation (19), upon exploiting the orthogonality relation (equation (18)) the boundary condition in equation (14e ) renders the simple relations

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\p}{\partial} \displaystyle \left.\frac{\partial g_{n}(\tau)}{\partial\tau}\right|_{\tau=\tau_0} =\tau_0 c^2 n (n-1)B_{n-1} ,~n = 2,3,\dots \label{final2} \nonumber \end{align} \tag{ 20 }$$

For a given slip velocity, thus a given set of amplitudes B_n of the slip modes $V_n(\zeta)$, equations (14d ) and (20) evaluated for $n \geqslant 1$ provide a system of coupled linear equations for the last unknown coefficients in the stream function, $C_{n\geqslant 3}$ and $D_{n \geqslant 2}$. Inspection of this system reveals that it splits into two subsystems of coupled equations, one involving only the coefficients of even index, n  =  2k, and the other one involving only the coefficients of odd index n  =  2k  +  1 (see appendix C). Furthermore, from equation (20) it can be inferred that in the case that the slip velocity is given by a pure slip mode $V_{n_0}(\zeta)$, i.e. $B_n = \delta_{n,n_0}$ for $n \geqslant 1$, then the parity of n₀ selects one of the two subsystems. If n₀ is even, then all the coefficients C_k and D_k of even index k vanish and the stream function, equation (12), involves only the functions g_k of odd index k, while for odd n₀ all the coefficients C_k and D_k of odd index k vanish and the stream function, equation (12), involves only the functions g_k of even index k. Finally, we note that a pure slip mode $V_{n_0}(\zeta)$ selects, as discussed above, either the odd or even number terms in the series representation of the stream function, equation (12), but not only a single term in the stream function expansion, as is the case for the spherical squirmer. Accordingly, even simple distributions of active slip velocities on the surface of the particle can give rise to quite complex hydrodynamic flows around the squirmer.

With these general results, we can proceed to the discussion of prolate squirmers. We will focus on the usual quantities employed to characterize an axisymmetric squirmer [3, 4, 50], i.e. the velocity U of the squirmer, the magnitude S of the stresslet associated with the squirmer (which determines the far-field hydrodynamic flow of the squirmer), and the general characteristics of the hydrodynamic flow ${\bf v}({\bf r})$ around the squirmer.

By combining equations (14e ), (19), (18) and (15), and noting that the polynomials $P_n^1(x)$ are even (odd) functions of x for n odd (even), one arrives at

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\m}[1]{\mathrm{#1}} \displaystyle \label{eq:U_terms_Bn} U(\tau_0) &= \frac{\tau_0^2}{2} \sum_{n \geqslant 1, n~\mathrm{odd}} B_n \int_{-1}^{1} {\rm d}x \frac{P_1^1(x) P_n^1(x)}{\tau_0^2 - x^2}\nonumber \\ &= \sum_{n \geqslant 1, n~\mathrm{odd}} B_n U_n(\tau_0).\nonumber \end{align} \tag{ 21 }$$

Accordingly, it follows that (a) a squirmer may exhibit self-motility (i.e. $U \neq 0$) only if the slip velocity involves at least an odd index n slip mode $V_n$; (b) in contrast to the case of a spherical squirmer, for which the velocity is determined solely by the slip mode n  =  1 irrespective of the details of the slip velocity ${\bf v}_s$, for a spheroidal squirmer all the slip modes of odd index contribute to the velocity (see also figure 2); consequently, (c) spheroidal squirmers with B₁  =  0 can be self-motile (due to contributions from other odd index slip modes), and spheroidal squirmers with $B_1 \neq 0$ can yet be non-motile if the contributions from other slip modes of odd index n precisely balance the contribution of the mode B₁ (which clearly pinpoints the shortcomings of a model with only two slip modes as in [54]).

Zoom In Zoom Out Reset image size

In terms of the dependence on the slenderness parameter r_e (which determines the value of $\tau_0$, see equation (8)), there are two findings. (d) for every slip mode n, the contribution |U_n | (in absolute value) is a decreasing function of r_e (see figure 2); second, (e) while at low values of the aspect ratio the contributions |U_n | of the n  >  1 slip modes are significant, the contributions from the modes $n \geqslant 3$ decay steeply with increasing r_e and become negligible, compared to |U₁| (which remains non-zero), as the aspect ratio r_e of the spheroid approaches that of a sphere ($r_e \to 1^-$). This ensures a smooth transition into the spherical case, where, as previously mentioned, higher mode squirmers are not motile and $\newcommand{\p}{\partial} U \propto B_1$. Finally, we note that, by comparison with the expression in equation (14c ), one infers that the series in the last line of equation (21) is proportional to the coefficient F₂ in the expansion of the stream function.

In what concerns the stresslet $S(\tau)$, which (similarly to the case of a spherical squirmer) allows classification into pullers (positive stresslet, S  >  0), pushers (negative stresslet, S  <  0), and neutral swimmers (vanishing stresslet, S  =  0), we note that it can also be expressed in terms of the slip velocity $u_s(\zeta)$ as [53]:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\m}[1]{\mathrm{#1}} \displaystyle \label{eq:stresslet} S &=-\frac{2A\mu}{F(r_e^{-1})J(r_e^{-1})} \int^1_{-1}u_s(\zeta)\zeta\sqrt{\frac{r_e^{-2}(1-\zeta^2)}{\zeta^2+r_e^{-2} (1-\zeta^2)}}{\rm d}\zeta\nonumber \\ &= \sum_{n \geqslant 1, n~\mathrm{even}} \mu B_n S_n(\tau_0)\nonumber \end{align} \tag{ 22 }$$

where A denotes the surface area of the spheroid and

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle F(x)=\frac{1}{(x^2-1)^2}\left(-3x^2+\frac{x(1+2x^2)}{\sqrt{1-x^2}}\cos^{-1}(x)\right), \nonumber \end{align*}$$

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle J(x)=1+\frac{x^2}{\sqrt{x^2-1}}\cos^{-1}\left(\frac{1}{x}\right). \nonumber \end{align*}$$

As in the case of the velocity, there are certain significant differences from the case of a spherical squirmer (see figure 3): (a) a necessary condition for a prolate squirmer to exhibit a nonvanishing stresslet is that at least one even (k  =  2n) mode $V_k$ contributes to the slip velocity; (b) as for the velocity U (where more than a single mode contributes), the stresslet depends on all even squirmer modes; hence, (c) the stresslet contribution of $B_2 \neq 0$ can be offset or even inverted by other even, active modes $B_k \neq 0$; (d) At a given aspect ratio r_e , the contribution (in absolute value) |S_2n | is a decreasing function of n; and (e) while for elongated shapes (small values r_e ) the contributions |S_n | from the slip modes $n \geqslant 4$ are significant, these contributions are steeply decreasing towards zero with increasing r_e towards the value 1⁻. In contrast, |S₂| is increasing with r_e . As in the case of the velocity, this behavior ensures the smooth transition to the case of a spherical shape, where $\newcommand{\p}{\partial} S \propto B_2$.

Zoom In Zoom Out Reset image size

Since the stresslet is, by definition, the amplitude of the r⁻² far-field term in the flow field (in the lab frame) of the squirmer, the series in the last line of equation (22) can be connected with one of the coefficients in the expansion of the stream function as follows. In the laboratory frame, which is related to the one ($\newcommand{\p}{\partial} \psi_{\rm particle}$) in the co-moving frame via

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\p}{\partial} \displaystyle \psi_{\rm lab}=\psi_{\rm particle}+\frac{1}{2} U_0(1-r_e^2)(\tau^2-1)(1-\zeta^2), \label{labframe} \nonumber \end{align} \tag{ 23 }$$

the slowest decaying term with the distance $\tau$ from the squirmer is $-\frac{C_3}{90}{\rm G}_0(\tau){\rm G}_3(\zeta)$; by equation (10), this term leads to a contribution $\sim C_3/\tau^2$ to the flow. Accordingly, $\newcommand{\p}{\partial} S \propto C_3$ and the pusher or puller squirmers ($S \neq 0$) indeed exhibit the expected far-field hydrodynamics, while for the neutral squirmers (S  =  0) the far-field flow necessarily decays at least as  ∼$1/\tau^3$.

Turning now to the flow field around the prolate squirmer, we will discuss separately the flow fields generated by the first few pure slip modes, i.e. the cases $B_n = \delta_{n,n_0}$ with n₀  =  1,2,3,4; these flows are shown, in the laboratory frame, in figure 4. From the discussion above, we know that only a subset (either the odd index ones, if n₀ is even, or vice versa) of the terms in the series representation, equation (12), of the stream functions contributes to the flow. Since the metric factors are even functions of $\zeta$, while $G_n(\zeta)$ is an odd (even) function of $\zeta$ when n is odd (even), the flow has the following fore-aft symmetries. For n₀ odd, the stream function involves the functions G_k of index k an even number, and thus $\newcommand{\p}{\partial} \psi(\tau,\zeta) = \psi(\tau,-\zeta)$; this implies $v_\tau(\tau,\zeta) = - v_\tau(\tau,-\zeta)$ and $v_\zeta(\tau,\zeta) = v_\zeta(\tau,-\zeta)$ (see figures 5 and 6); i.e. the z– (x–) flow components are fore-aft (anti)symmetric (see figure 4), and, accordingly, it contributes to the motility because it provides a 'fore to back' streaming. Vice versa, for n₀ an even number one has $v_\tau(\tau,\zeta) = v_\tau(\tau,-\zeta)$ and $v_\zeta(\tau,\zeta) = - v_\zeta(\tau,-\zeta)$, i.e. the x– (z–) flow component are fore-aft symmetric (see figure 4); consequently, they cannot be associated to a motile particle.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, we note that the similar analysis and results for the case of an oblate squirmer can be obtained from the available ones for a prolate squirmer via a simple transformation of the coordinates system and of the stream function (i.e. a mapping), as discussed in the appendix A. As an illustration of using this mapping, the results shown in figure 4, corresponding to a prolate squirmer, have been used to determine the flows around the corresponding (i.e. of slenderness parameters $r'_e = 1/r_e$) oblate squirmers induced by the pure slip modes $B_n =\delta_{n,n_0}$ with n₀  =  1,2,3,4; these are shown in figure 7. The analysis of the velocity U and stresslet S (see figures 8 and 9) for oblate spheroids leads to conclusions that are similar with those drawn in the case of prolate shapes. The only significant difference is that now for oblate squirmers the contributions from the higher orders n slip modes decay to zero with decreasing aspect ratio $r_e \to 1^+$; again, this ensures a smooth transition to the case of a spherical shape, where only B₁ or B₂ are relevant.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.1. Self-phoretic particle

Squirmers with a wide range of active slip modes occur naturally in the context of model self-phoretic particles. One of the often employed realizations of such systems consists of micrometer-sized silica or polystyrene spherical particles partially coated with a Pt layer and immersed in an aqueous peroxide solution [58, 64, 65]. The catalytic decomposition of the peroxide at the Pt side creates gradients in the chemical composition of the suspension; as in the case of classic phoresis [66], these gradients, in conjunction with the interaction between the colloid and the various molecular species in solution, give rise to self-phoretic motility [67]. The mechanism of steady-state motility can be intuitively understood in terms of the creation of a so-called phoretic slip velocity tangential to the surface of the particle [20, 66, 68]. For such chemically active, axi-symmetric, spherical particles in unbounded solutions, the approximation of phoretic slip velocity leads to a straightforward mapping [55] onto a squirmer model; accordingly, the translational velocity of the particle and the hydrodynamic flow around the particle in terms of the phoretic slip can be directly inferred from the corresponding results in [4].

A similar mapping can be developed from a spheroidal, self-phoretic colloid to a spheroidal squirmer model for spheroidal, chemically active colloids. Following [60] the slip velocity of such particles can be written as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle v_{\rm S}({\bf r}_{\rm P})=\sum_{l\geqslant0}\frac{b} {c} c_l(\chi) {\rm Q}_l(\tau_0) {\rm P}^1_l(\zeta) (\tau_0^2-\zeta^2)^{-\frac{1}{2}}, \label{vschem} \nonumber \end{align} \tag{ 24 }$$

where b is the so-called phoretic mobility (for simplicity, here assumed to be a constant) over the surface of the particle, and ${\rm Q}_l(\tau)$ denotes the Legendre polynomial of the second kind [61]. The coverage $\chi = h-1$ is defined in terms of the height of the active cap (measured, from the bottom apex, in units of b_z ), i.e. $\chi = -1$ corresponds to a chemically inactive spheroid and $\chi = +1$ corresponds to all the surface being active. The coverage dependent coefficients $c_l(\chi)$ (see [60]) describe the decoration of the surface of the particle by the chemically active element (e.g. Pt in the example discussed above) as an expansion in Legendre polynomials ${\rm P}_l(\zeta)$. Knowledge of these parameters allows us to formulate a slip velocity boundary condition similar to (20), i.e.

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\p}{\partial} \displaystyle \frac{\partial g_{n}(\tau)}{\partial\tau}\Big|_{\tau=\tau_0} = \tau_0c^2n(n-1)\tilde{B}_{n-1}; \label{final3} \nonumber \end{align} \tag{ 25 }$$

by identifying the effective squirmer modes $\tilde{B}_l =\frac{b}{c\cdot\tau_0}c_l{\rm Q}_l(\tau_0)$, the desired mapping is achieved. As an illustration of this mapping, we show in figure 10 (left column) the flow fields for particles with parameters r_e   =  0.3, 0.5 and $\chi=0,-0.3$, respectively, and compare with the corresponding results obtained by direct numerical solutions obtained using BEM [63]. (Note that, if necessary, the accuracy of the analytical estimate can be systematically improved simply by increasing the order of the truncation in the series expansion of the stream function (see the appendix D)).

Zoom In Zoom Out Reset image size

We have studied in detail the most general axi-symmetric spheroidal squirmer and we have shown that, in analogy with the situation for spherical shapes [55], model chemically active, self-phoretic colloids can be mapped onto squirmers. By using the semiseparable ansatz derived in [52] for the stream function, and representing the active slip (squirming) velocity of the squirmer in a suitable basis (equation (19), chosen such that in the limit of a spherical shape it smoothly transforms into the usually employed expansion for the classical spherical squirmer [4], the velocity of the squirmer, the stresslet of the squirmer, and the hydrodynamic flow around the squirmer have been determined analytically (section 4). The corresponding series representations have been validated by cross-checking against direct numerical calculations, obtained by using the BEM, of the flow around the squirmer (appendix E).

The main conclusion emerging from the study is that for spheroidal squirmers (or self-phoretic particles) the squirming modes beyond the second are, in general, as important as the first two ones in what concerns the contributions to the velocity and stresslet of the particle (and, implicitly, to the flow field around the particle, even in the far-field). The velocity is contributed by all the odd-index components (but none of the even-index ones) of the slip velocity; accordingly, in contrast with the case of spherical squirmers, it is possible to have spheroidal squirmers with a non-zero first mode and yet not motile, as well as ones missing the first slip mode and yet motile. Similarly, the stresslet value is contributed by all the even-index slip modes; thus, distinctly from the case of spherical squirmers, one can have pushers/pullers even if the second slip mode is vanishing, as well as neutral squirmers in spite of a non-vanishing second slip mode. Finally, even a single slip mode leads to a large number of non-vanishing terms in the series expansion of the stream function, and thus spheroidal squirmers with simple distributions of slip on their surface can lead to very complex flows around them (see figures 4 and 7).

This raises the interesting speculative question as whether the spheroidal shape is providing an evolutionary advantage; i.e. with small modifications of the squirming pattern, e.g. switching from a sole B₁ mode to a sole B₃ mode, a microrganism could maintain its velocity unchanged but dramatically alter the topology of the flow around it (compare the first and third rows in figure 4). In other words, compared to a squirmer with spherical shape, for which multiple modes must be simultaneously activated in order to change the structure of the flow, a spheroidal squirmer possesses simple means for acting in hydrodynamic disguise, which can be advantageous as either predator or prey.

The similarity between oblate and prolate spheroids allows us to obtain the flow field around an oblate microswimmer of aspect ratio $r_e' > 1$ as a series in the oblate coordinates by using a mapping from the results, in prolate coordinates, corresponding to a prolate microswimmer with aspect ratio $r_e = 1/r_e' < 1$ (and vice versa).

The oblate spheroidal coordinate system is defined by

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \lambda = \frac{1}{2\bar{c}}\Big(\sqrt{x^2+y^2+(z-{\rm i}\bar{c})^2}+ &\sqrt{x^2+y^2+(z+{\rm i}\bar{c})^2}\large{\Big)}, \nonumber \\ \zeta= \frac{1}{-2{\rm i}\bar{c}}\Big(\sqrt{x^2+y^2+(z-{\rm i}\bar{c})^2}-&\sqrt{x^2+y^2+(z+{\rm i}\bar{c})^2}\Big), \nonumber \\ \varphi = \arctan & \left(\frac{y}{x}\right), \nonumber \end{align*}$$

with $0\leqslant\lambda\leqslant\infty$ and $\bar{c}=\sqrt{r_e'^2-1}$, and the corresponding Lamé metric coefficients are given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:metric_coef_oblate} & h_\zeta =\bar{c}\frac{\sqrt{\lambda^2+\zeta^2}}{\sqrt{1-\zeta^2}}, ~~h_\lambda=\bar{c}\frac{\sqrt{\lambda^2+\zeta^2}}{\sqrt{\lambda^2+1}}, \nonumber \\ & h_\varphi = \bar{c} \sqrt{\lambda^2+1}\sqrt{1-\zeta^2}.\nonumber \end{align} \tag{ A.1 }$$

Noting that the oblate coordinates and the metric factors can be obtained from the expressions of the corresponding prolate coordinates via the transformations [52]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau\rightarrow {\rm i}\lambda \qquad c \rightarrow -{\rm i}\bar{c} , \nonumber \end{align} \tag{ A.2 }$$

one concludes that the equation and boundary conditions obeyed by the stream function in oblate coordinates in the domain outside the oblate of aspect ratio $r_e' = 1/r_e$ can be obtained, by using the same transformation, from the ones in prolate coordinates outside a prolate of aspect ratio r_e . Accordingly, it follows that $\newcommand{\p}{\partial} \psi_{{\rm obl},r_e'=1/r_e}(\lambda,\zeta) = \psi_{{\rm pro},r_e}(\tau = {\rm i}\lambda,\zeta)$. The flow field then follows as the curl of the stream function, i.e.

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\p}{\partial} \displaystyle {\bf v}_{{\rm obl},r_e'}(\lambda,\xi,\varphi)= \nabla \times \left(\frac{\psi_{{\rm obl},r_e'} (\lambda,\xi)}{h_\varphi} \,\,{\bf e}_\varphi \right) . \nonumber \end{align} \tag{ A.3 }$$

The Gegenbauer functions of first kind G_n are also known as the Gegenbauer polynomials $C^\alpha _n$ with parameter $\alpha = -1/2$ [61]; they are defined in terms of the Legendre polynomials P_n as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle G_n(x)=\frac{1}{2n-1}\left(P_{n-2}(x)-P_n(x)\right), \nonumber \end{align} \tag{ B.1 }$$

and fulfill the orthogonality relation

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \int_{-1}^1\frac{G_n(x)G_m(x)}{1-x^2}{\rm d}x=\frac{2}{n(n-1)(2n-1)}\delta_{n,m}\quad n,m\geqslant2. \nonumber \end{align} \tag{ B.2 }$$

The Gegenbauer functions can be related to the associated Legendre polynomials ${\rm P}^1_n$, which are defined in terms of the Legendre polynomials P_n as

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle (1-x^2)^{1/2}{\rm P}^1_l = -(1-x^2)\frac{{\rm d}}{{\rm d}x}P_l(x). \nonumber \end{align*}$$

Using the relations $\frac{x^2-1}{n}\frac{{\rm d}}{{\rm d}x}P_n(x)=x{\rm P}_n(x)-{\rm P }_{n-1} (x)$ and $(n+1){\rm P}_{n+1}(x)=$ $(2n+1)x{\rm P}_n(x)-n{\rm P}_{n-1}(x)$, one arrives at the following relation between the lth associated Legendre polynomial and the l  +  1th Gegenbauer function.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle (1-x^2)^{1/2}{\rm P}^1_l =-(l^2+l){\rm G}_{l+1}(x). \nonumber \end{align} \tag{ B.3 }$$

The coefficients C_n and D_n appear in the functions $g_k(\tau)$, entering the series expansion of the stream function, at various indexes k (e.g. $C_{n \geqslant 4}$ appears at both k  =  n  −  2 and k  =  n); thus the infinite system of linear equations is strongly coupled. However, since the even and odd terms are not entering the same equations, the system splits into two decoupled subsystems, which are solved by using different methods.

The first subsystem consists of the even numbered terms in the series expansion of $\newcommand{\p}{\partial} \psi(\tau,\zeta)$ and the set of conditions (equations (14a ), (14c ), (14d ) and (20)). Because each of the equations (14a ) and (14c ) fix one of the even index coefficients (i.e. C₂ and F₂), the equations (14d ) and (20) evaluated at n  =  2 involve only two unknowns, C₄ and D₂, and thus can be solved as a sub-subsystem. With C₄ and D₂ known, the rest of the coefficients, up to the order $N_{\rm max}$ at which the system is truncated (i.e. B_k is set to zero for $k > N_{\rm max}$), are solved iteratively by noting that equations (14d ) and (20) evaluated at k  =  n involve only two unknowns, the rest of the coefficients being already determined in the previous iterations up to k  =  n  −  2.

The second subsystem includes all odd-index terms in the series expansion of $\newcommand{\p}{\partial} \psi(\tau,\zeta)$ and equations (14d ) and (20). Complementing the previous case, only the even squirmer modes contribute. However, unlike in the case of the first subsystem, there are no lower level decouplings; accordingly, this subsystem is solved in the standard manner by truncation at the cut-off $N_{\rm max}$ (above which all the coefficients are set to zero) and inversion of the resulting finite system of linear equations.

The choice of $N_{\rm max}$ is done by varying the value of the cut-off and testing the changes in the coefficients. For the cases we analyzed in this work, a value $N_{\rm max} = 16$ was found to be sufficient.

To show the importance of the higher orders in the case of a self-phoretic swimmer, we compare the analytic results for keeping different numbers of squirmer modes B_n in figure C1. Even for an aspect ratio r_e   =  0.8 (i.e. close to a spherical shape), the higher orders have significant influence on the flow field around the particle: e.g. compare the occurrence of a region of high magnitude flow near the point $(0,-1.5)$ (behind the particle) instead of the correct location at near the point $(0,1.5)$ (in front of the particle).

Zoom In Zoom Out Reset image size

In figures E1–E3 we present a more detailed comparison of the results obtained by using the series representation (top row) with numerical results obtained by using the BEM to directly solve the governing equations (second row). In addition, the relative error, defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta v= \frac{\sqrt{(v_{x,{\rm ana}}-v_{x,{\rm num}})^2+ (v_{z,{\rm ana}}-v_{z,{\rm num}})^2}}{\frac{1}{2}(\sqrt{v_{x,{\rm ana}}^2+v_{z,{\rm ana}}^2}+\sqrt{v_{x, {\rm num}}^2+v_{z,{\rm num}}^2})} \label{error} \nonumber \end{align} \tag{ E.1 }$$

is shown in the bottom row. The comparison confirms the expected quantitative agreement, with regions of significant relative error ($\Delta v > 10^{-1}$) corresponding precisely to the regions where the flow is anyway very weak ($|\frac{({\bf v}/{\bf r})}{U}|<10^{-3}$).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size