Active Stokesian dynamics

Abstract

Since its development, Stokesian dynamics has been a leading approach for the dynamic simulation of suspensions of particles at arbitrary concentrations with full hydrodynamic interactions. Although developed originally for the simulation of passive particle suspensions, the Stokesian dynamics framework is equally well suited to the analysis and dynamic simulation of suspensions of active particles, as we elucidate here. We show how the reciprocal theorem can be used to formulate the exact dynamics for a suspension of arbitrary active particles, and then show how the Stokesian dynamics method provides a rigorous way to approximate and compute the dynamics of dense active suspensions where many-body hydrodynamic interactions are important.

1. Introduction

Active matter is a term used to describe matter that is composed of a large number of self-propelled active ‘particles’ that individually convert stored or ambient energy into systematic motion (Schweitzer & Farmer 2007; Morozov 2017). The interaction of many of these individual active particles can lead to complex collective dynamics (Ramaswamy 2010). Natural examples include a flock of birds, a school of fish, or a suspension of bacteria (Toner, Tu & Ramaswamy 2005), but active matter may also be composed of synthetic active particles (Bechinger et al. 2016). These out-of-equilibrium systems are most often in fluids, so understanding their dynamics and rheology involves a connection between fluid–body interactions and non-equilibrium statistical physics (Marchetti et al. 2013; Saintillan 2018).

The study of active matter at small scales is complicated by the fact that the Stokes equations, which govern momentum conservation of Newtonian fluids when inertia is negligible, feature a long-range decay of fluid disturbances (Happel & Brenner 1965). Because of this, active particles interact through the fluid over distances that are long relative to their individual size, and to properly capture the effect of the fluid in these systems, one may need to sum hydrodynamic interactions between all bodies, particularly at higher particle concentrations.

The difficulty of capturing accurately many-body hydrodynamic interactions is well known from the study of suspensions of passive particles, where early efforts to sum hydrodynamic interactions in infinite suspensions were plagued by problems of divergent sums (see, for example, the long literature on sedimenting particles; Davis & Acrivos 1985), eventually overcome by the pioneering work of Batchelor (1972), Jeffrey (1974), Hinch (1977), O'Brien (1979) and others. The Stokesian dynamics method that was developed soon after facilitated the efficient dynamic simulation of passive particle suspensions at arbitrary concentrations (Brady & Bossis 1988). The essential basis of the Stokesian dynamics method is a mixed asymptotic approach wherein hydrodynamic forces on particles due to interactions are computed distinctly when the particles are in close proximity versus widely separated. When the particles are widely separated, the method sums many-body hydrodynamic reflections between particles through inversion of a truncated grand mobility tensor, whereas when the particles are in close proximity, pairwise additive lubrication forces are used (Durlofsky, Brady & Bossis 1987). When the suspension is infinite or periodic, a modification of the method introduced by O'Brien (1979) is used to obtain absolutely convergent expressions for the hydrodynamic interactions among all particles, suitable for the numerical simulation of a wide range of problems from sedimentation to rheology (Brady et al. 1988). Since its inception, the Stokesian dynamics method has served as a foundational tool for the development of our understanding of suspension mechanics in the last several decades.

Unlike passive suspensions, in active suspensions each active particle in the fluid is endowed with non-trivial boundary conditions due to activity, and constantly injects energy into the fluid. Many advances have been made in understanding the dynamics of individual swimming microorganisms (biological and synthetic), from the pioneering work of Taylor (1951) through to several detailed reviews of microscale locomotion research (Lighthill 1976; Brennen & Winet 1977; Lauga & Powers 2009). However, in a fashion similar to the early development of the passive suspension literature, the majority of research on collective locomotion of many bodies and active suspensions has emphasized dilute suspensions where swimmer–swimmer interactions are greatly simplified, and far-field approximations are still valid (Saintillan 2018). Interesting phenomena, such as particle clustering (motility-induced phase separation), have been observed for dense suspensions of active particles (Bechinger et al. 2016), but very often numerical simulation of these suspensions is done with active Brownian particle models that neglect hydrodynamic interactions entirely (Cates & Tailleur 2015). Others have used approaches for active suspensions that only approximate the Stokes equations, such as multiparticle collision dynamics (Zöttl & Stark 2014) or lattice Boltzmann methods (Stenhammar et al. 2017) that still may not be accurate for very dense concentrations. Results for simplified swimmers in concentrated suspensions display qualitative differences (Ishikawa, Locsei & Pedley 2008; Evans et al. 2011; Alarcón & Pagonabarraga 2013; Matas-Navarro et al. 2014; Zöttl & Stark 2014; Thutupalli et al. 2018). Some argue that hydrodynamic interactions act to suppress phase separation in active matter (Matas-Navarro et al. 2014), while others have shown that hydrodynamic interactions with boundaries can control phase separation (Thutupalli et al. 2018). A complete understanding of the connection between individual particle activity, the hydrodynamic interactions between many particles that arise as a consequence of this activity, and the role this plays in the macroscopic dynamics of concentrated active suspensions has not been developed.

As we discuss in the following, the Stokesian dynamics methodology is easily adapted for the dynamic simulation of suspensions of active particles at any concentration, and as with passive suspensions, is particularly well suited for dense concentrations and periodic boundary conditions. The mathematical structure of the dynamical equations remains essentially unchanged between passive and active particles, so any implementation of the Stokesian dynamics method for passive particles may be modified simply and easily for use with active particles. By making the connection to passive suspensions and Stokesian dynamics, we obtain directly a long-developed framework for summing hydrodynamic interactions, adding non-hydrodynamic forces, including Brownian motion and constructing the equivalent Smoluchowski equations, for active suspensions with full hydrodynamic interactions (Burkholder & Brady 2018). We believe that this mathematical structure provides an ideal formalism for theoretical analysis of active suspensions (Burkholder & Brady 2019), much as it has for passive suspensions (Brady 1993a,b).

The Stokesian dynamics method was first adapted for use with self-propelled active particles by Mehandia & Nott (2008). In their work, they introduced spheres each with a prescribed virtual propulsive force, that interact through a prescribed stresslet whose magnitude sets the size of the virtual propulsive force, and an induced stresslet caused by particle rigidity in a bulk flow. The dynamics of these active spheres was then solved numerically using the Stokesian dynamics framework. The authors found that near-field interactions appeared important even at low concentrations, as particles tended to cluster, and they found qualitative differences in the dynamics between low- and high-volume fractions. Despite the novelty, the authors did not specify how the propulsive force arises from the surface boundary conditions, or how to generalize this approach. Shortly afterwards, Ishikawa et al. (2008) adapted the Stokesian dynamics framework for use with spherical particles with a prescribed tangential slip velocity, so-called squirmer particles (Ishikawa et al. 2008). Using their own previous results for two-body hydrodynamic interaction between squirmer particles (Ishikawa, Simmonds & Pedley 2006), Ishikawa et al. (2008) were able to incorporate both near-field interactions and many-body far-field interactions for the study of dense suspensions of (2-mode) squirmer particles. This framework was then used to study the rheology (Ishikawa & Pedley 2007b), diffusion (Ishikawa & Pedley 2007a) and coherent structures (Ishikawa & Pedley 2008) of these active suspensions. The Stokesian dynamics framework was then extended for use with passive and active spherical particles that had a fairly general surface velocity field (but were individually immotile), which could be linked together to form complex swimming assemblies (Swan et al. 2011). That machinery was then used to simulate a number of model swimming microorganisms, from pusher and puller swimmers to helical flagella, by using assemblies of spherical particles (Swan et al. 2011). Recently, a higher-order Stokesian-dynamics-like approach (without lubrication), namely constructing mobility tensors by a higher-order moment expansion of the boundary integral equations (Ichiki 2002), was developed for suspensions of squirmer particles by using tensorial spherical harmonics (Singh, Ghose & Adhikari 2015). Here, we show that this previous literature may all be synthesized into a fairly general theory for the dynamics of suspensions of arbitrary active particles. This work illustrates the similarity of hydrodynamic interactions in ‘passive’ colloidal suspensions versus active suspensions and how particle activity can be incorporated directly into the Stokesian dynamics approach.

Alternative approaches to the Stokesian dynamics methodology for the numerical simulation of passive and active suspensions have also been developed recently. An approach known as the force-coupling method (FCM) (Maxey & Patel 2001) instead uses regularized (Gaussian) force distributions in the fluid and integrals over the entire fluid domain (rather than just particle boundaries) to construct (far-field) mobilities. The FCM was extended to incorporate fluctuating hydrodynamics for passive suspensions (Keaveny 2014; Delmotte & Keaveny 2015) and then to squirmer active particles, up to the stresslet level (Delmotte et al. 2015). Similarly, a fluctuating immersed boundary method (FIBM) was also developed to simulate suspensions of Brownian particles (Delong et al. 2014). Like the FCM, the FIBM uses an explicit (fluctuating) solvent but uses kernels developed by Peskin (2002) for the immersed boundary method to mediate the fluid–particle interactions. This immersed boundary approach has also been extended to simulate rigid assemblies of particles and active particles (at the monopole level) (Balboa Usabiaga et al. 2016). A key benefit of these approaches is that exact Green's functions need not be found, which for complex boundaries may not be feasible at all. The method used to construct mobilities can be chosen depending on the boundary conditions, and this is the approach taken in the rigid multiblob method with exact Green's functions (at the Rotne–Prager level) for unbounded or half-space simulations, while using the FIBM for confined geometries (Balboa Usabiaga et al. 2016; Sprinkle et al. 2017). These methods are both fast and have well-supported code bases, and although they generally include only far-field hydrodynamic interactions, in principle, lubrication could also be added. Finally, rather than forming the mobilities by expanding in tensorial spherical harmonics as done by Singh et al. (2015) and also shown here, a recent approach (for non-Brownian suspensions) uses instead vector spherical harmonics (Corona & Veerapaneni 2018; Yan et al. 2020), which seems to lead to simpler formulas for higher-order terms, allowing simulation of the near field without resorting to a lubrication approximation.

We begin by developing a general kinematic description of an arbitrary active particle in § 2. Next, we show how the reciprocal theorem can be used to yield the exact dynamics for a suspension of $N$ arbitrary active particles in § 3. We then show how the Stokesian dynamics technique is used for the approximation and dynamic simulation of these exact equations in § 4.

2. Kinematics of an active particle

Consider an active particle identified with the region $\mathcal {B}$ as shown in figure 1. Changes in the spatial configuration of the active particle can be described by a map $\boldsymbol {\chi }$ from a reference configuration $\mathcal {B}_0$ such that $\boldsymbol {x}=\boldsymbol {\chi }(\boldsymbol {X},t)$ for $\boldsymbol {x}\in \mathcal {B}$ and $\boldsymbol {X}\in \mathcal {B}_0$. The motion of the body can be decomposed into shape change $\boldsymbol {\chi }_s$, which represents the swimming gait of the active particle, and rigid-body motion $\boldsymbol {\chi }_r$, which arises as a consequence of interaction with the fluid, so that

(2.1)\begin{equation} \boldsymbol{\chi}(\boldsymbol{X},t) = \boldsymbol{x}^c(t) + \boldsymbol{\varTheta}(t)\boldsymbol{\cdot}(\boldsymbol{\chi}_s(\boldsymbol{X},t)-\boldsymbol{\chi}_s(\boldsymbol{X}_0,t)), \end{equation}

where $\boldsymbol {x}^c$ is the translation, and $\boldsymbol {\varTheta }$ is the rotation (about $\boldsymbol {\chi }_s(\boldsymbol {X}_0,t)$) of the body under the action of $\boldsymbol {\chi }_r$. Upon differentiation, we obtain the velocity of the body,

(2.2)\begin{equation} \boldsymbol{u}(\boldsymbol{x}\in \mathcal{B}) = \boldsymbol{U}+\boldsymbol{\varOmega}\times\boldsymbol{r}+\boldsymbol{u}^s, \end{equation}

where the translational velocity is $\boldsymbol {U} = {\rm d}\boldsymbol {x}^c/{\rm d}t$, while the rotational velocity $\boldsymbol {\varOmega }$ is defined by ${\rm d}\boldsymbol {\varTheta }/{\rm d} t = \boldsymbol {\varOmega }\times \boldsymbol {\varTheta }$ and $\boldsymbol {r} = \boldsymbol {\varTheta }(t)\boldsymbol {\cdot }(\boldsymbol {\chi }_s(\boldsymbol {X},t)-\boldsymbol {\chi }_s(\boldsymbol {X}_0,t))$. The deformation velocity due to shape change is

(2.3)\begin{equation} \boldsymbol{u}^s=\boldsymbol{\varTheta}\boldsymbol{\cdot}\frac{{\rm d}(\boldsymbol{\chi}_s(\boldsymbol{X},t)-\boldsymbol{\chi}_s(\boldsymbol{X}_0,t))}{{\rm d} t}, \end{equation}

where the last term is the deformation velocity in the unoriented configuration,

(2.4)\begin{equation} \frac{{\rm d}(\boldsymbol{\chi}_s(\boldsymbol{X},t)-\boldsymbol{\chi}_s(\boldsymbol{X}_0,t))}{{\rm d}t}=\boldsymbol{\varTheta}^{{-}1}\boldsymbol{\cdot}\boldsymbol{u}^s\equiv\tilde{\boldsymbol{u}}^s. \end{equation}

In a purely kinematic description of the activity of the particle, we would consider the shape change $\boldsymbol {\chi }_s$ to be prescribed and then solve for the rigid-body translation and rotation of the active particle such that momentum (of the particle and of the fluid) is conserved. Often, in the process of modelling an active particle, the details of the deformation are coarse-grained away and one simply prescribes a surface velocity on a suitable reference configuration, $\tilde {\boldsymbol {u}}(\boldsymbol {X}\in \partial \mathcal {B}_0)$, such as ciliary beating on the surface of a microorganism that is represented as a slip velocity over a fixed geometry (Blake 1971).

Figure 1. Schematic of a deforming active particle.

3. Dynamics of active particles

Consider a suspension of $N$ particles, each labelled $\mathcal {B}_i$ where $i\in [1,N]$, immersed in an arbitrary background flow denoted $\boldsymbol {u}^\infty$. The disturbance velocity field generated by the particles is

(3.1)\begin{equation} \boldsymbol{u}' = \boldsymbol{u}-\boldsymbol{u}^\infty. \end{equation}

Neglecting the inertia of the active particles and of the Newtonian fluid in which they are immersed, the rigid-body dynamics of active particles is governed by an instantaneous force balance

(3.2)\begin{equation} \boldsymbol{\mathsf{F}}+\boldsymbol{\mathsf{F}}_{ext}=\boldsymbol{\mathsf{0}}, \end{equation}

where $\boldsymbol{\mathsf{F}}$ and $\boldsymbol{\mathsf{F}}_{ext}$ are respectively $6N$-dimensional vectors of hydrodynamic and external (or interparticle) forces/torques on all $N$ particles.

In general, the hydrodynamic forces may be easily shown, by the reciprocal theorem of low Reynolds number hydrodynamics, to be weighted integrals of the boundary traction during rigid-body motion:

(3.3)\begin{equation} \boldsymbol{\mathsf{F}} = \sum_i\int_{\partial\mathcal{B}_i}\boldsymbol{u}'\boldsymbol{\cdot}(\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\mathsf{T}}_{\boldsymbol{\mathsf{U}}})\,\text{d} S, \end{equation}

where $\boldsymbol {n}$ is the normal to the surface $\partial \mathcal {B}_i$ pointing into the fluid. The tensor field $\boldsymbol{\mathsf{T}}_{\boldsymbol{\mathsf{U}}}$ connects the rigid-body motion of $6N$ particles to stress (see Appendix A for a detailed derivation). Substitution of the boundary conditions of each active particle, (2.2), into (3.3) yields a decomposition of the hydrodynamic forces into three separate forces due to each aspect of the boundary motion: the hydrodynamic ‘swim’ force (or thrust),

(3.4)\begin{equation} \boldsymbol{\mathsf{F}}_s = \sum_i\int_{\partial\mathcal{B}_i}\boldsymbol{u}^s_i\boldsymbol{\cdot}(\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\mathsf{T}}_{\boldsymbol{\mathsf{U}}})\,\text{d} S, \end{equation}

generated by each active particle as if held fixed in an otherwise quiescent fluid; the hydrodynamic drag force on each particle,

(3.5)\begin{equation} \boldsymbol{\mathsf{F}}_\infty ={-}\sum_i\int_{\partial\mathcal{B}_i}\boldsymbol{u}^\infty\boldsymbol{\cdot}(\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\mathsf{T}}_{\boldsymbol{\mathsf{U}}})\,\text{d} S, \end{equation}

as if inactive and held fixed in a background flow; and the hydrodynamic drag due to the rigid-body motion of each particle,

(3.6)\begin{equation} \boldsymbol{\mathsf{F}}_d=\sum_i\int_{\partial\mathcal{B}_i}(\boldsymbol{U}_i+\boldsymbol{\varOmega}_i\times\boldsymbol{r}_i)\boldsymbol{\cdot} (\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\mathsf{T}}_{\boldsymbol{\mathsf{U}}})\,\text{d} S={-}\boldsymbol{\mathsf{R}}_{\boldsymbol{\mathsf{FU}}}\boldsymbol{\cdot}\boldsymbol{\mathsf{U}}, \end{equation}

as if inactive (passive) in an otherwise quiescent fluid. The latter is written in terms of a ($6N\times 6N$) resistance tensor, which is the linear operator that gives hydrodynamic forces due to rigid-body translational/rotational velocities $\boldsymbol{\mathsf{U}}$ (another $6N$-dimensional vector).

Substitution of these forces into (3.2) and inversion of the resistance tensor gives

(3.7)\begin{equation} \boldsymbol{\mathsf{U}}=\boldsymbol{\mathsf{R}}_{\boldsymbol{\mathsf{FU}}}^{{-}1}\boldsymbol{\cdot}\left[\boldsymbol{\mathsf{F}}_{ext}+\boldsymbol{\mathsf{F}}_s+\boldsymbol{\mathsf{F}}_\infty\right]. \end{equation}

This relationship simply states that the rigid-body motion of active particles is linearly related to the forces exerted by or on those particles. The deterministic formula in (3.7) is exact and completely general; it governs the dynamics of a suspension of active (and passive) particles of arbitrary shape and activity in a general background flow. A stochastic Brownian force may also be included in the above force balance, with the associated thermal drift term that arises upon elimination of inertial degrees of freedom:

(3.8)\begin{equation} \boldsymbol{\mathsf{U}}=\boldsymbol{\mathsf{R}}_{\boldsymbol{\mathsf{FU}}}^{{-}1}\boldsymbol{\cdot}\left[\boldsymbol{\mathsf{F}}_{ext}+\boldsymbol{\mathsf{F}}_s+ \boldsymbol{\mathsf{F}}_\infty+\boldsymbol{\mathsf{F}}_B \right]+k_B T\,\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\mathsf{R}}_{\boldsymbol{\mathsf{FU}}}^{{-}1}. \end{equation}

The vector $\boldsymbol{\mathsf{U}}$ now represents discrete changes in position and orientation over an interval $\Delta t$, and the Brownian force is $\boldsymbol{\mathsf{F}}_B = \sqrt {2k_BT/\Delta t}\,\boldsymbol{\mathsf{R}}_{\boldsymbol{\mathsf{FU}}}^{1/2}\boldsymbol {\cdot } \boldsymbol {\varPsi }$, where $k_B$ is the Boltzmann constant, $T$ is the fluid temperature and $\boldsymbol {\varPsi }$ is a vector of standard Gaussian random variables.

Although (3.7) and (3.8) are exact, to compute the dynamics, the tensor field $\boldsymbol{\mathsf{T}}_{\boldsymbol{\mathsf{U}}}$ would need to be found at each instant, and this is prohibitively expensive for suspensions of large numbers of particles (and more so if they are changing shape). Instead, an approximate approach used in Stokesian dynamics is to evaluate a truncated set of moments of the traction operator $\boldsymbol {n}\boldsymbol {\cdot }\boldsymbol{\mathsf{T}}_{\boldsymbol{\mathsf{U}}}$ on the surfaces of the particles $\partial \mathcal {B}_i$. We outline this approach for spherical active bodies below, where the approach is particularly elegant and simplified, but the methodology can certainly be extended to anisotropic bodies (Claeys & Brady 1993a,b,c; Nasouri & Elfring 2018).

3.1. Spherical moments

In order to facilitate computation, we use the fact that one may write an arbitrary function on a sphere in terms of an expansion in irreducible tensors of the unit normal $\boldsymbol {n}$ (dimensionless tensorial spherical harmonics) (Hess 2015). In this way, the deformation velocity of each active particle may be written as

(3.9)\begin{equation} \boldsymbol{u}^s(\boldsymbol{x} \in \partial \mathcal{B}_i)=\boldsymbol{C}^{(1)}_i+\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{C}^{(2)}_i+{\overparen{\boldsymbol{n}\boldsymbol{n}}}:\boldsymbol{C}^{(3)}_i+\cdots, \end{equation}

where in this shorthand notation the superscript indicates the tensor order of the coefficients, and the overbracket means the irreducible (or fully symmetric and traceless) part of the tensor (see Appendix B for further details).

The coefficient tensors $\boldsymbol {C}^{(n)}$ may be obtained easily by appealing to the orthogonality of the tensorial spherical harmonics

(3.10)\begin{equation} \boldsymbol{C}_i^{(n+1)} = \frac{1}{4{\rm \pi} a^2}\int_{\partial\mathcal{B}_i} w_n \overparen{\boldsymbol{n}^n}\boldsymbol{u}^s\,\text{d} S, \end{equation}

where the weight is $w_n = (2n+1)!!/n!$. Hence we see that the coefficient tensors are (weighted) irreducible moments of $\boldsymbol {u}^s$. We can recast these coefficients in more familiar terms by separating symmetric and antisymmetric parts of $\boldsymbol {C}_i^{(2)}$,

(3.11)$$\begin{gather} \boldsymbol{U}^s_i =\frac{1}{4{\rm \pi} a_i^2}\int_{\partial\mathcal{B}_i}\boldsymbol{u}^s\,\text{d} S, \end{gather}$$

(3.12)$$\begin{gather}\boldsymbol{\varOmega}^s_i = \frac{3}{8{\rm \pi} a_i^4}\int_{\partial\mathcal{B}_i}\boldsymbol{r}_i\times\boldsymbol{u}^s\,\text{d} S, \end{gather}$$

(3.13)$$\begin{gather}\boldsymbol{E}^s_i = \frac{3}{4{\rm \pi} a_i^4}\int_{\partial\mathcal{B}_i} \left[\frac{1}{2}\,(\boldsymbol{r}_i\boldsymbol{u}^s+\boldsymbol{u}^s\boldsymbol{r}_i)\right]\,\text{d} S, \end{gather}$$

to rewrite the surface velocity in familiar form (Swan et al. 2011)

(3.14)\begin{equation} \boldsymbol{u}^s(\boldsymbol{x} \in \partial \mathcal{B}_i)=\boldsymbol{U}^s_i+\boldsymbol{\varOmega}_i^s\times\boldsymbol{r}_i+\boldsymbol{r}_i\boldsymbol{\cdot}\boldsymbol{E}_i^s+\cdots. \end{equation}

We can likewise express the background flow in terms of a moment expansion, and in this way write in a consistent fashion for the disturbance field

(3.15)\begin{align} \boldsymbol{u}'(\boldsymbol{x} \in \partial \mathcal{B}_i)=\boldsymbol{U}_i+\boldsymbol{U}_i^s-\boldsymbol{U}_i^\infty+ (\boldsymbol{\varOmega}_i+\boldsymbol{\varOmega}_i^s-\boldsymbol{\varOmega}_i^\infty)\times\boldsymbol{r}_i+\boldsymbol{r}_i\boldsymbol{\cdot}(\boldsymbol{E}_i^s-\boldsymbol{E}_i^\infty)+\cdots, \end{align}

where $\boldsymbol {U}_i^\infty$, $\boldsymbol {\varOmega }_i^\infty$ and $\boldsymbol {E}_i^\infty$ are defined as in (3.11)–(3.13) in terms of moments of the background flow $\boldsymbol {u}^\infty$. If the background flow is linear, then $\boldsymbol {\varOmega }_i^\infty =\boldsymbol {\varOmega }^\infty$ and $\boldsymbol {E}_i^\infty =\boldsymbol {E}^\infty$ are constants everywhere in the flow (but still may be arbitrary functions of time).

Using the expansion (3.15) in (3.3), one may write the hydrodynamic forces in terms of a set of moments of the traction operator $\boldsymbol {n}\boldsymbol {\cdot }\boldsymbol{\mathsf{T}}_{\boldsymbol{\mathsf{U}}}$ on the surfaces of the particles $\partial \mathcal {B}_i$ (forming resistance tensors):

(3.16)\begin{equation} \boldsymbol{\mathsf{F}} ={-}\boldsymbol{\mathsf{R}}_{\boldsymbol{\mathsf{FU}}}\boldsymbol{\cdot}(\boldsymbol{\mathsf{U}}+\boldsymbol{\mathsf{U}}^s-\boldsymbol{\mathsf{U}}^\infty)- \boldsymbol{\mathsf{R}}_{\boldsymbol{\mathsf{FE}}}:(\boldsymbol{\mathsf{E}}^s-\boldsymbol{\mathsf{E}}^\infty)+\cdots. \end{equation}

Now using the above expression for the hydrodynamic forces, together with Newton's second law (3.2), we obtain the translational and rotational velocities of the spherical active particles:

(3.17)\begin{equation} \boldsymbol{\mathsf{U}}={-}\boldsymbol{\mathsf{U}}^s+\boldsymbol{\mathsf{U}}^\infty+\boldsymbol{\mathsf{R}}_{\boldsymbol{\mathsf{FU}}}^{{-}1}\boldsymbol{\cdot} \left[\boldsymbol{\mathsf{F}}_{ext}-\boldsymbol{\mathsf{R}}_{\boldsymbol{\mathsf{FE}}}:(\boldsymbol{\mathsf{E}}^s-\boldsymbol{\mathsf{E}}^\infty)+\cdots\right]. \end{equation}

These equations have the same functional form as the governing equations of motion for passive particles, except with active particles one takes the difference in moments of the background flow to surface velocity, e.g. $\boldsymbol{\mathsf{E}}^\infty \rightarrow \boldsymbol{\mathsf{E}}^\infty -\boldsymbol{\mathsf{E}}^s$. If the particles are passive, $\boldsymbol {u}^s=\boldsymbol {0}$, then we recover equations of motion for passive particles in a background flow (Brady & Bossis 1988); however, computing the far-field hydrodynamic interactions of active spherical particles is no more difficult than for passive spherical particles, assuming that the velocities on the boundaries of the active particle, $\boldsymbol {u}^s_i$, are prescribed.

If hydrodynamic interactions are completely neglected, then the particles all move with their respective single-particle velocities $\boldsymbol{\mathsf{U}} =-\boldsymbol{\mathsf{U}}^s+\boldsymbol{\mathsf{U}}^\infty$ (higher-order moments do not contribute to self-propulsion for isolated spherical particles by symmetry), and we recover the classic result for single active spheres (Anderson & Prieve 1991; Stone & Samuel 1996; Elfring 2015). It is technically possible to devise a perfect stealth swimmer that does not disturb the surrounding fluid by setting $\boldsymbol {u}^s = \boldsymbol {U}^s$, for example by the jetting mechanism proposed by Spagnolie & Lauga (2010), but this is an extreme case, and in general, higher-order moments lead to hydrodynamic interactions. We emphasize that hydrodynamic interactions due to moments of the particle activity enter in exactly equivalent form to interactions due to moments of the background flow. For example, the leading-order change in the dynamics of the particles due to hydrodynamic interactions is given by $\boldsymbol{\mathsf{U}}+\boldsymbol{\mathsf{U}}^s-\boldsymbol{\mathsf{U}}^\infty =-\boldsymbol{\mathsf{R}}_{\boldsymbol{\mathsf{FU}}}^{-1}\boldsymbol {\cdot } \boldsymbol{\mathsf{R}}_{\boldsymbol{\mathsf{FE}}}:(\boldsymbol{\mathsf{E}}^s-\boldsymbol{\mathsf{E}}^\infty )$, where the resistance tensors $\boldsymbol{\mathsf{R}}_{\boldsymbol{\mathsf{FU}}}^{-1}\boldsymbol {\cdot }\boldsymbol{\mathsf{R}}_{\boldsymbol{\mathsf{FE}}}$ act to couple the particles in precisely the same fashion for active particles as passive particles. In Appendix D, we give the leading-order hydrodynamic interactions (a dilute approximation) in the mobility formulation more commonly employed in the literature.

We see that the leading-order hydrodynamic interactions due to activity are given by the symmetric first moment of activity, $\boldsymbol {E}^s_i$, of each active particle. This is not a surprise as the term $\boldsymbol {E}^s$ sets the active component of the stresslet $\boldsymbol {S}$ (Ishikawa et al. 2006; Lauga & Michelin 2016; Nasouri & Elfring 2018) of individual spherical active particles, where

(3.18)\begin{equation} \boldsymbol{S} = \frac{1}{2}\int_{\partial\mathcal{B}}[\boldsymbol{r}\boldsymbol{\sigma}\boldsymbol{\cdot}\boldsymbol{n}+\boldsymbol{\sigma}\boldsymbol{\cdot} \boldsymbol{n}\boldsymbol{r} -2\eta(\boldsymbol{u}\boldsymbol{n}+\boldsymbol{n}\boldsymbol{u})] \,\text{d} S = \frac{20{\rm \pi}\eta a^3}{3}\left(\boldsymbol{E}^\infty-\boldsymbol{E}^s\right). \end{equation}

The ‘active strain rate’ $\boldsymbol {E}^s$ can be zero, but then the leading-order term will generally arise at the second-moment level for self-motile active particles (that is, ones with non-zero surface averaged velocity). This is the case for so-called neutral squirmers (see § 3.2 below) or symmetric phoretic particles (Michelin & Lauga 2014). This illustrates why it is particularly important to incorporate higher-order moments to capture accurately hydrodynamic interactions between active particles (Singh et al. 2015). Active particles may also be immotile (not self propelling), meaning that $\boldsymbol {U}^s=\boldsymbol {0}$ but higher-order moments are non-zero. A canonical example of immotile active particles is extensile microtubule bundles that are driven by kinesin motors (Sanchez et al. 2012). Immotile active particles, sometimes called ‘shakers’ in contrast to ‘movers’ that are motile (Hatwalne et al. 2004), still interact hydrodynamically through higher-order active moments in much the same way as motile active particles.

These equations, above all, simply reflect the linear relationship between velocity and force moments. Using still more compact notation for all disturbance velocity moments $\mathcal {U}' = [\boldsymbol{\mathsf{U}}+\boldsymbol{\mathsf{U}}^s-\boldsymbol{\mathsf{U}}^\infty, \boldsymbol{\mathsf{E}}^s-\boldsymbol{\mathsf{E}}^\infty,\ldots ]^{\rm T}$ and hydrodynamic force moments $\mathcal {F} = [\boldsymbol{\mathsf{F}}, \boldsymbol{\mathsf{S}},\ldots ]^{\rm T}$, we may write, more generally, the linear relationship

(3.19)\begin{equation} \mathcal{F}={-}\mathcal{R}\boldsymbol{\cdot}\mathcal{U}', \end{equation}

where $\mathcal {R}$ is the grand resistance tensor, an (unbounded) linear operator that maps velocity moments to force moments. In this notation, the hydrodynamic force is written compactly as

(3.20)\begin{equation} \boldsymbol{\mathsf{F}} ={-}\mathcal{R}_{\boldsymbol{\mathsf{F}}\mathcal{U}}\boldsymbol{\cdot}\mathcal{U}'. \end{equation}

In order to capture the dynamics of active particles, we seek an effective and efficient way to form $\mathcal {R}_{\boldsymbol{\mathsf{F}}\mathcal {U}}$. The grand resistance tensor is a purely geometric operator, depending only on the position (and orientation if they are anisotropic) of each active particle (Happel & Brenner 1965). Perhaps less obvious is that the grand resistance tensor does not depend on the prescribed surface activity of the particles, and is thus identical to the case when they are passive. This also applies to particles in a bounded geometry – $\mathcal {R}$ is a function of geometry only (Swan & Brady 2007, 2010). We have assumed here that the surface activity of the particles is prescribed; however, the surface activity may depend on the traction on the boundary, as it would, for example, for biological particles that have power-limited surface actuation, but the linear relationship between force moments and velocity moments makes it straightforward to prescribe force moments (Swan et al. 2011).

We focus here on spherical particles, as is common in the literature for colloidal suspensions; however, the method described above can be generalized to other geometries. We formed moments of forces and velocities by projection onto tensorial spherical harmonics, but for other geometries, a more suitable basis for the vector fields on the particle surfaces $\partial \mathcal {B}_i$ would be used. Alternatively, and more generally, one may perform Taylor series expansion of the boundary integral equations about the centre of each particle, which naturally projects tractions onto force moments for particles of arbitrary geometry (see recent work by Swan et al. (2011) and Nasouri & Elfring (2018) for details of this method applied to active particles). A problem with this approach is that the particle activity $\boldsymbol {u}^s$ might be defined only on the particle surfaces (in the form of surface slip as in § 3.2), but this difficulty can be ameliorated by lifting $\boldsymbol {u}^s$ to a suitably continuous function defined in $\mathbb {R}^3$. Despite this complication, fundamentally, the linear relationship between velocity and force moments remains, regardless of geometry.

3.2. Squirmers

A squirmer is a spherical particle whose surface slip velocity is tangential to the surface (Pedley 2016). Most often, the slip velocity is taken to be axisymmetric; here, the direction of the axis of symmetry of the particle is denoted by $\boldsymbol {p}$ (the particle director). A purely tangential slip velocity is of course an idealization, but one that arises quite naturally, for example in the limit of small-amplitude deformations that are projected onto a time-averaged spherical manifold (Lighthill 1952; Blake 1971), or as the outer solution of phoretic flow due to chemical concentrations confined to a thin layer near the sphere surface (Anderson 1989; Golestanian, Liverpool & Ajdari 2005). The slip velocity is typically written as an expansion in Legendre polynomials:

(3.21)\begin{equation} \boldsymbol{u}^s = \sum_n \frac{2}{n(n+1)}P_n'\, (\boldsymbol{p}\boldsymbol{\cdot}\boldsymbol{n})\boldsymbol{p}\boldsymbol{\cdot}\left[B_n(\boldsymbol{I}-\boldsymbol{n}\boldsymbol{n})+C_n\boldsymbol{I}\times\boldsymbol{n}\right], \end{equation}

where $P_n$ is the Legendre polynomial of degree $n$, and $P_n'(x)=({{\rm d}}/{{\rm d}\kern 0.06em x})P_n(x)$. The polar slip coefficients $B_n$ are often called ‘squirming’ modes, while the azimuthal slip (or ‘swirl’), with coefficients $C_n$, is not often considered but can lead to particle spin, for instance (change $C_n$ to $C_n {n(n+1)}/{2a^{n+1}}$ for the coefficients used by Pak & Lauga 2014). Recasting the slip velocity in terms of irreducible tensors of the surface normal such that

(3.22)\begin{equation} \boldsymbol{u}^s(\boldsymbol{x} \in \partial \mathcal{B})=\boldsymbol{U}^s+\boldsymbol{\varOmega}^s\times\boldsymbol{r}+\boldsymbol{r}\boldsymbol{\cdot}\boldsymbol{E}^s+ \overparen{\boldsymbol{r}^2}:\boldsymbol{B}^s+\overparen{\boldsymbol{r}^3}\odot\,\boldsymbol{C}^s+\cdots, \end{equation}

we obtain

(3.23)$$\begin{gather} \boldsymbol{U}^s ={-}\tfrac{2}{3}B_1 \boldsymbol{p}, \end{gather}$$

(3.24)$$\begin{gather}\boldsymbol{\varOmega}^s = \frac{1}{a^3}\,C_1\boldsymbol{p}, \end{gather}$$

(3.25)$$\begin{gather}a\boldsymbol{E}^s ={-}\tfrac{3}{5}B_2 \overparen{\boldsymbol{p}\boldsymbol{p}}, \end{gather}$$

(3.26)$$\begin{gather}a^2\boldsymbol{B}^s = B_1\boldsymbol{\varDelta}^{2}\boldsymbol{\cdot}\boldsymbol{p}-\tfrac{5}{7}B_3\overparen{\boldsymbol{p}\boldsymbol{p}\boldsymbol{p}}-C_2(\boldsymbol{\varDelta}^2\boldsymbol{\cdot}\boldsymbol{p})\times\boldsymbol{p}, \end{gather}$$

(3.27)$$\begin{gather}a^3\boldsymbol{C}^s = B_2\boldsymbol{\varDelta}^{3}:\overparen{\boldsymbol{p}\boldsymbol{p}}-\tfrac{35}{36}B_4\overparen{\boldsymbol{p}\boldsymbol{p}\boldsymbol{p}\boldsymbol{p}}- \tfrac{5}{4}C_3(\boldsymbol{\varDelta}^3:\overparen{\boldsymbol{p}\boldsymbol{p}})\times\boldsymbol{p}, \end{gather}$$

where $\boldsymbol {\varDelta }^{n}$ is an isotropic $2n$-order tensor that when applied on a tensor of rank $n$, projects onto the symmetric traceless part of that tensor (see Appendix B for further details). By symmetry, each coefficient is necessarily composed only of products of the particle director $\boldsymbol {p}$. We see that the swimming speed is given by the first squirming mode $B_1$ as $\boldsymbol {U}=-\boldsymbol {U}^s=(2/3)B_1\boldsymbol {p}$ for an isolated squirmer, but note that the first mode also contributes a higher-order term to hydrodynamic interactions between particles embedded in $\boldsymbol {B}^s$. The stresslet due to surface activity of a particle is given by the second squirming mode, $\boldsymbol {S} = 4{\rm \pi} \eta a^2 B_2\overparen {\boldsymbol {p}\boldsymbol {p}}$. This determines if a squirmer is a pusher or a puller, but can easily be zero – a so-called neutral squirmer – and in that case, the leading-order term contributing to hydrodynamic interactions is necessarily given by $B_1$ (and $B_3$ if non-zero). Azimuthal slip leads naturally to rotation given by the $C_1$ mode, $\boldsymbol {\varOmega }=-\boldsymbol {\varOmega }^s = -(C_1/a^3)\boldsymbol {p}$ for an isolated squirmer, while the $C_2$ mode leads to a rotlet dipole contribution in the far field (Pak & Lauga 2014).

As an example of the framework developed here, consider an active squirmer particle, labelled $\mathcal {B}_1$, in the presence of a freely suspended passive sphere, labelled $\mathcal {B}_2$, as shown in figure 2. Using (3.17), we obtain the velocities of the two particles in terms of moments of the surface activity of the active particle:

(3.28)$$\begin{gather} \boldsymbol{\mathsf{U}}_1 ={-}\boldsymbol{\mathsf{U}}_1^s - (\boldsymbol{\mathsf{M}}_{\boldsymbol{\mathsf{UF}}}^{11}\boldsymbol{\cdot} \boldsymbol{\mathsf{R}}_{\boldsymbol{\mathsf{FE}}}^{11}+ \boldsymbol{\mathsf{M}}_{\boldsymbol{\mathsf{UF}}}^{12}\boldsymbol{\cdot} \boldsymbol{\mathsf{R}}_{\boldsymbol{\mathsf{FE}}}^{21}):\boldsymbol{E}_1^s+\cdots, \end{gather}$$

(3.29)$$\begin{gather}\boldsymbol{\mathsf{U}}_2 ={-} (\boldsymbol{\mathsf{M}}_{\boldsymbol{\mathsf{UF}}}^{21}\boldsymbol{\cdot} \boldsymbol{\mathsf{R}}_{\boldsymbol{\mathsf{FE}}}^{11} + \boldsymbol{\mathsf{M}}_{\boldsymbol{\mathsf{UF}}}^{22}\boldsymbol{\cdot} \boldsymbol{\mathsf{R}}_{\boldsymbol{\mathsf{FE}}}^{21}):\boldsymbol{E}_1^s+\cdots, \end{gather}$$

where the superscripts, for example $\boldsymbol{\mathsf{R}}_{\boldsymbol{\mathsf{FE}}}^{\alpha \beta }$, indicate the linear relationship between particle $\alpha$ and particle $\beta$, while $\boldsymbol{\mathsf{M}}_{\boldsymbol{\mathsf{UF}}} = \boldsymbol{\mathsf{R}}_{\boldsymbol{\mathsf{FU}}}^{-1}$. The first term on the right-hand side of (3.28) represents the self-propulsion of the active particle, while the second term represents the change in the velocity due to hydrodynamic interactions induced by the surface strain rate of the active particle $\boldsymbol {E}^s_1$ (and higher-order moments). Hydrodynamic interactions also induce the motion of the passive particle. In essence, the moments of $\boldsymbol {u}^s$ on the active particle result in a ‘swim’ force on both particles, which must then be balanced by drag due to rigid-body motion. The trajectories of both active and passive particles are illustrated in figure 2.

Figure 2. Trajectory of a (pusher) active particle (labelled $\mathcal {B}_1$) in the presence of a passive particle (labelled $\mathcal {B}_2$).

As another example, consider a suspension of immotile ‘shaker’ particles. Using (3.17), the velocities of the particles are

(3.30)\begin{equation} \boldsymbol{\mathsf{U}}_i ={-} \sum_j\sum_k(\boldsymbol{\mathsf{M}}_{\boldsymbol{\mathsf{UF}}}^{ij}\boldsymbol{\cdot} \boldsymbol{\mathsf{R}}_{\boldsymbol{\mathsf{FE}}}^{jk}):\boldsymbol{E}_k^s+\cdots. \end{equation}

Here, there is no self-propulsion, only the effects of hydrodynamic interactions induced by the surface strain rate of the active particles $\boldsymbol {E}^s_k$ (and higher-order moments), as illustrated in figure 3.

Figure 3. An illustration of dynamics mediated by hydrodynamic interactions between immotile active particles (shakers). (a) 512 shaker particles are initially randomly oriented on a cubic lattice. Hydrodynamic interactions alone drive particle dynamics and mixing (b,c). Images are snapshots in time from left to right.

3.3. Assemblies

As detailed by Swan et al. (2011), assemblies of active (or passive) particles can be dealt with easily within the Stokesian dynamics framework, and in what follows we outline the presentation given in that work.

A set of particles may be constrained to move as a rigid body, that is, particle $\alpha$ in a rigid assembly $A$ will move as

(3.31)$$\begin{gather} \boldsymbol{U}_{\alpha} = \boldsymbol{U}_A+\boldsymbol{\varOmega}_A\times(\boldsymbol{x}_\alpha-\boldsymbol{x}_A), \end{gather}$$

(3.32)$$\begin{gather}\boldsymbol{\varOmega}_{\alpha} = \boldsymbol{\varOmega}_A, \end{gather}$$

where $\boldsymbol {x}_A$ is a convenient point on the assembly. Following the notation in Swan et al. (2011), this may be written compactly in terms of six-dimensional vectors as $\boldsymbol{\mathsf{U}}_\alpha = \boldsymbol {\varSigma }_{\alpha A}^{\rm T}\boldsymbol {\cdot }\boldsymbol{\mathsf{U}}_A$, where $\boldsymbol {\varSigma }_{\alpha A}^{\rm T}$ projects the translational and rotational velocity of the assembly onto particle $\alpha$. The rigid-body translational and rotational velocities of all $N$ particles in an assembly may then be written in terms of $6N$-dimensional vectors and tensors as

(3.33) \begin{equation} \boldsymbol{\mathsf{U}} = \boldsymbol{\varSigma}_{A}^{{\rm T}}\boldsymbol{\cdot}\boldsymbol{\mathsf{U}}_A. \end{equation}

The forces and torques that enforce the rigid constraints on the assembly, $\boldsymbol{\mathsf{F}}_c$, must be included in the sum of forces on the particles:

(3.34)\begin{equation} \boldsymbol{\mathsf{F}}+\boldsymbol{\mathsf{F}}_c+\boldsymbol{\mathsf{F}}_{ext}=\boldsymbol{\mathsf{0}}. \end{equation}

These constraint forces are internal forces, and as such exert no net force or torque on the assembly. This may be written as

(3.35)\begin{equation} \boldsymbol{\varSigma}_{A}\boldsymbol{\cdot}\boldsymbol{\mathsf{F}}_c = \boldsymbol{\mathsf{0}}, \end{equation}

where the operator $\boldsymbol {\varSigma }_{A}$, the transpose of the projection above, sums forces and torques (about $\boldsymbol {x}_A$) on the assembly. In this way, the force balance on the assembly is

(3.36)\begin{equation} \boldsymbol{\varSigma}_{A}\boldsymbol{\cdot}\boldsymbol{\mathsf{F}}+\boldsymbol{\varSigma}_{A}\boldsymbol{\cdot}\boldsymbol{\mathsf{F}}_{ext}=\boldsymbol{\mathsf{0}}. \end{equation}

Substitution of the relevant hydrodynamic forces and the kinematic constraint in (3.33) into this force balance leads to the rigid-body motion of the assembly given by

(3.37)\begin{equation} \boldsymbol{\mathsf{U}}_A=\left[\boldsymbol{\varSigma}_{A}\boldsymbol{\cdot}\boldsymbol{\mathsf{R}}_{\boldsymbol{\mathsf{FU}}}\boldsymbol{\cdot} \boldsymbol{\varSigma}_A^{\rm T}\right]^{{-}1}\boldsymbol{\cdot}\boldsymbol{\varSigma}_{A}\boldsymbol{\cdot}\left(\boldsymbol{\mathsf{F}}_{ext}+\boldsymbol{\mathsf{F}}_s+\boldsymbol{\mathsf{F}}_\infty\right), \end{equation}

where $\boldsymbol {\varSigma }_{A}\boldsymbol {\cdot }\boldsymbol{\mathsf{R}}_{\boldsymbol{\mathsf{FU}}}\boldsymbol {\cdot }\boldsymbol {\varSigma }_A^{\rm T}=\boldsymbol{\mathsf{R}}_{\boldsymbol{\mathsf{FU}}}^A$ is the hydrodynamic resistance of the assembly. Equation (3.37) is an exact description of the dynamics of an assembly of active or passive particles; no approximation has yet been made. In particular, we note that while (3.37) yields the instantaneous rigid-body motion of the assembly, it does not mean that the assembly cannot deform. Indeed, through the prescription of the activity of each particle, by way of $\boldsymbol {u}^s$, we may construct an assembly of virtually any shape and kinematics. This approach is also extended straightforwardly to multiple assemblies through an extended operator $\boldsymbol {\varSigma }$ that sums forces on each assembly as shown by Swan et al. (2011). As discussed above, a natural method of solution is to use Stokesian dynamics to resolve hydrodynamic forces as a truncated set of moments.

As an illustrative example of a deforming assembly, consider a simple reciprocal two-sphere (or dumbbell) swimmer (see figure 4a). In this model swimmer, two spheres labelled $\mathcal {B}_1$ and $\mathcal {B}_2$ (where $\mathcal {B}_A = \mathcal {B}_1\cup \mathcal {B}_2$), of radii $a$ and $\lambda a$, respectively, have a prescribed distance between their centres, $L(t)$, that changes periodically in time. We describe the shape change of this swimmer as the motion of sphere $\mathcal {B}_1$ relative to sphere $\mathcal {B}_2$; in this way, $\boldsymbol {u}^s$ is non-zero only on $\mathcal {B}_1$. Written in terms of an expansion in moments as in (3.22), $\boldsymbol {u}^s(\boldsymbol {x}\in \mathcal {B}_1)=\boldsymbol {U}_1^s=\dot {L}\boldsymbol {p}$, with all other terms exactly zero, while $\boldsymbol {u}^s(\boldsymbol {x}\in \mathcal {B}_2)=\boldsymbol {0}$. The total velocity of $\mathcal {B}_2$ is then due solely to the rigid-body motion of the assembly $\boldsymbol {u}(\boldsymbol {x}\in \mathcal {B}_2)=\boldsymbol {U}_A$, while $\mathcal {B}_1$ has an additional component due to shape change, $\boldsymbol {u}(\boldsymbol {x}\in \mathcal {B}_1)=\boldsymbol {U}_A+\boldsymbol {U}_1^s$. By symmetry, this swimmer does not rotate, $\boldsymbol {\varOmega }_A=\boldsymbol {0}$. The choice of reference is not unique and affects what is delineated as rigid-body motion versus shape change at any particular instant; however, typically we are concerned with the time-averaged motion of the body, which is invariant to the choice of reference for periodic gaits, and the flexibility allows one to take advantage of simplifications implied by a particular choice.

Figure 4. (a) Schematic of a reciprocal dumbbell swimmer. (b) Schematic of a dumbbell squirmer swimmer.

Due to the lack or rotation or torque, only the force–velocity resistance tensor of the assembly, $\boldsymbol {R}_{FU}^A$, and the linear operator that gives stress due to translation, $\boldsymbol{\mathsf{T}}_U$, are required. Substitution of the swim force (3.4) into (3.37) and simplification leads to

(3.38)\begin{equation} \boldsymbol{U}_A ={-}(\boldsymbol{R}_{FU}^{A})^{{-}1}\boldsymbol{\cdot}\left(\boldsymbol{R}_{FU}^{11}+\boldsymbol{R}_{FU}^{21}\right)\boldsymbol{\cdot}\boldsymbol{U}_1^s, \end{equation}

where the resistance tensors $\boldsymbol {R}_{FU}^{11}+\boldsymbol {R}_{FU}^{21}$ and $\boldsymbol {R}_{FU}^{A}=\boldsymbol {R}_{FU}^{11}+\boldsymbol {R}_{FU}^{12}+\boldsymbol {R}_{FU}^{21}+\boldsymbol {R}_{FU}^{22}$ are functions of the length $L(t)$, and hence depend on time. We may further simplify by noting that the propulsive force and velocity will be collinear with the axis of symmetry, so only a scalar coefficient for each resistance is required. A symmetric swimmer with $\lambda =1$ has $\boldsymbol {U}_A =-\frac {1}{2}\boldsymbol {U}_1^s= -\frac {1}{2}\dot {L}\boldsymbol {p}$, so the dumbbell moves opposite to the deformation with half the speed, as expected. This reciprocal motion clearly leads to zero net displacement over a period when $L(t)$ is periodic. Less obvious, but also true, is that this holds for any $\lambda$, by the scallop theorem (Purcell 1977).

The previous example was particularly straightforward because $\boldsymbol {u}^s$ was uniform, hence only the zeroth moment, $\boldsymbol {U}^s$, was non-zero. In contrast, consider a dumbbell swimmer with a fixed length $L=const.$, but where sphere $\mathcal {B}_1$ is a squirmer particle (see figure 4b), namely $\boldsymbol {u}^s(\boldsymbol {x}\in \mathcal {B}_1) = \boldsymbol {U}^s_1+\boldsymbol {r}\boldsymbol {\cdot }\boldsymbol {E}_1^s+\cdots$, with the moments of the surface velocity given by the squirming modes. In this case, the velocity of the assembly is given by

(3.39)\begin{equation} \boldsymbol{U}_A ={-}(\boldsymbol{R}_{FU}^{A})^{{-}1}\boldsymbol{\cdot}\left[\left(\boldsymbol{R}_{FU}^{11}+\boldsymbol{R}_{FU}^{21}\right)\boldsymbol{\cdot} \boldsymbol{U}_1^s+\left(\boldsymbol{R}_{FE}^{11}+\boldsymbol{R}_{FE}^{21}\right):\boldsymbol{E}_1^s+\cdots\right], \end{equation}

and when the spheres are equal in size, $\lambda = 1$, we have simply $\boldsymbol {U}_A = -\frac {1}{2}\boldsymbol {U}_1^s-(\boldsymbol {R}_{FU}^{A})^{-1}\boldsymbol {\cdot }$ $\left [\boldsymbol {R}_{FE}^{1}:\boldsymbol {E}_1^s+\cdots \right ]$. Note that this swimmer can self-propel even when $\boldsymbol {U}_1^s=\boldsymbol {0}$ due to hydrodynamic interactions with the second sphere.

4. Stokesian dynamics

The configuration-dependent $N$-body resistance tensors may be formed indirectly by first constructing the grand mobility tensor $\mathcal {M}=\mathcal {R}^{-1}$. In the Stokesian dynamics approach, one takes irreducible moments of the velocity field, as given by the boundary integral equation, over the surfaces of all the particles, yielding Faxén's laws for the velocity moments of the active particles (Batchelor 1972). If the boundary integral equations are also expanded in irreducible moments (Durlofsky et al. 1987), then we obtain a linear relationship between force and velocity moments:

(4.1)\begin{equation} \mathcal{U}' ={-}\mathcal{M}\boldsymbol{\cdot}\mathcal{F}. \end{equation}

For active particles, the force moments contain contributions from the double-layer kernel due to the surface activity (see Appendix C for details). The grand mobility tensor is then inverted to obtain the grand resistance tensor $\mathcal {R}=\mathcal {M}^{-1}$, thereby summing many-body hydrodynamic interactions among the particles (Durlofsky et al. 1987). In principle, to capture near-field lubrication effects, the entire unbounded set of moments would need to be computed, and in practice this is unfeasible. The coupling between the $m$th moment of velocity and $n$th moment of force scales as $r^{-(1+m+n)}$, so higher-order moments decay quite quickly with separation distance $r$ between two particles, and a reasonable and common far-field approximation of the mobility is to truncate at the first moment level – we label this truncated mobility $\mathcal {M}^{ff}$ (Swan et al. 2011). This level of approximation is inappropriate for particles that are nearly touching, and the compromise used within Stokesian dynamics is to use a mixed asymptotic approach wherein close interactions are computed separately using pairwise exact solutions (Durlofsky et al. 1987). In this manner, the hydrodynamic forces on the particles are decomposed,

(4.2)\begin{equation} \boldsymbol{\mathsf{F}} = \boldsymbol{\mathsf{F}}^{ff}+\boldsymbol{\mathsf{F}}^{2B,exact}-\boldsymbol{\mathsf{F}}^{2B,ff}, \end{equation}

into far-field interactions between many bodies $\boldsymbol{\mathsf{F}}^{ff}$, and two-body interactions computed exactly for nearby bodies $\boldsymbol{\mathsf{F}}^{2B,exact}$. Note that the last term arises because the far-field interactions must be removed between any two bodies where the interactions are computed exactly, to avoid double counting.

To render exact solutions for active two-body hydrodynamic interactions, one needs to obtain the tensor field $\boldsymbol{\mathsf{T}}_{\boldsymbol{\mathsf{U}}}$ from the two-particle rigid-body motion problem. The general passive two-sphere problem for arbitrary separations may be constructed from a basis of four simplified two-sphere problems that have all been solved in the literature and are nicely summarized by Sharifi-Mood, Mozaffari & Córdova-Figueroa (2016) and Papavassiliou & Alexander (2017) in the context of two-body interactions between diffusiophoretic Janus particles and spherical squirmers, respectively. Asymptotic solutions for lubrication interactions, which are valid strictly only when the particles are very close, may be used alternatively and are given by Ishikawa et al. (2006) for spherical squirmers.

It is important to note that unlike passive particles in a linear background flow, active particles can have higher-order velocity moments due to surface activity. In § 3.2, we showed that even two-mode squirmer particles contribute third- and fourth-order velocity moments. For far-field interactions, these higher-order moments may not be significant due to the decay of the associated flow disturbances. However, for near-field interactions there is no rationale, other than the convergence of the series of tensorial spherical harmonics, to discard the contributions of higher-order velocity moments in the swim force. If higher-order moments are nonetheless discarded, then the approach for active particles is virtually identical to that of passive particles in Stokesian dynamics; the dynamics are given by (3.17), and the resistance tensor used is modified to include both exact two-body interactions $\mathcal {R}^{2B,exact}$ and a truncation of the moment expansion valid for far-field interactions, $(\mathcal {M}^{ff})^{-1}$, such that

(4.3)\begin{equation} \mathcal{R} = (\mathcal{M}^{ff})^{{-}1}+\mathcal{R}^{2B,exact}-\mathcal{R}^{2B,ff}, \end{equation}

where the two-body interactions that are captured by the near-field approach must be subtracted in the far-field solution to avoid double counting (Swan et al. 2011).

More accurately, the exact two-body swim force contributions may be computed entirely separately, as done by Ishikawa et al. (2006), by integrating (3.4) directly. In this way, (3.7) is written as

(4.4)\begin{align} \boldsymbol{\mathsf{U}}&=\boldsymbol{\mathsf{U}}^\infty+\boldsymbol{\mathsf{R}}_{\boldsymbol{\mathsf{FU}}}^{{-}1}\boldsymbol{\cdot} \left[\boldsymbol{\mathsf{F}}_{ext}+\boldsymbol{\mathsf{R}}_{\boldsymbol{\mathsf{FE}}}:\boldsymbol{\mathsf{E}}^\infty+\cdots\right]\nonumber\\ &\quad+\boldsymbol{\mathsf{R}}_{\boldsymbol{\mathsf{FU}}}^{{-}1}\boldsymbol{\cdot}\left[\boldsymbol{\mathsf{F}}_{s}^{2B,exact}- \boldsymbol{\mathsf{R}}'_{\boldsymbol{\mathsf{FU}}}\boldsymbol{\cdot}\boldsymbol{\mathsf{U}}^s-\boldsymbol{\mathsf{R}}'_{\boldsymbol{\mathsf{FE}}}:\boldsymbol{\mathsf{E}}^s+\cdots\right]. \end{align}

Here, the terms on the first line represent the dynamics of passive spheres exactly as in conventional Stokesian dynamics, while the second line accounts for the contribution of activity, separated into two-body and far-field contributions. The primed resistance tensors contribute only far-field interactions

(4.5)\begin{equation} \mathcal{R}' = (\mathcal{M}^{ff})^{{-}1}-\mathcal{R}^{2B,ff}. \end{equation}

If Brownian motion is included, then it is $\boldsymbol{\mathsf{R}}_{\boldsymbol{\mathsf{FU}}}$ from the total resistance (4.3) that sets the magnitude of the Brownian force.

4.1. Infinite suspensions

The method described above was for a finite system of $N$ active particles for which the fluid can be assumed to decay in the far field. For an infinite or periodic suspension of particles (active or passive), no such assumption can be made, and indeed naively extending $N\rightarrow \infty$ leads to divergent integrals, a problem that plagued the earlier suspension literature (Batchelor 1972). Brady et al. (1988) adapted the method of O'Brien (1979) wherein the fluid domain for a set of particles is bounded by a large macroscopic surface over which suspension averages can be performed. Specifically, $\boldsymbol {U}^\infty$, $\boldsymbol {\varOmega }^\infty$ and $\boldsymbol {E}^\infty$ become the average values of the suspension – particle plus fluid. Individual particle motion is then relative to the volume-averaged quantities. Suspension-averaged terms serve to regularize the formulas leading to absolutely convergent expressions for fluid and particle velocities. Periodic boundary conditions may then be employed easily, and as Brady et al. (1988) showed, the far-field mobility matrix $\mathcal {M}^{ff}$ may be simply replaced by the appropriated Ewald-summed mobility matrix $\mathcal {M}^{ff*}$. As discussed above, the mobility matrix is unchanged if particles are active or passive; only the force and velocity moments are altered by activity, and mobility is unchanged whether or not the suspension averaged quantities are non-zero. Therefore, the Ewald-summed mobility matrix used for periodic passive suspensions is unchanged for active suspensions (Ishikawa et al. 2008). It is important to note that for self-propulsion there is no net volume displacement of material as the body moves: as the body advances, an equal volume of fluid moves in the opposite direction. In contrast, a body moving in response to an external force drags fluid along with it, and to have no net flux of mass, an external pressure gradient must be imposed.

4.2. Accelerated methods

Since the original development of the Stokesian dynamics method (Durlofsky et al. 1987; Brady et al. 1988), which naively requires $O(N^3)$ computations due to the inversion of the far-field mobility matrix, there have been several numerical implementations that have improved the algorithmic efficiency of the method. These include an $O(N \ln N)$ deterministic accelerated Stokesian dynamics method (Sierou & Brady 2001), an $O(N^{1.25} \ln N)$ Brownian accelerated Stokesian dynamics method (Banchio & Brady 2003), an $O(N \ln N)$ spectral Ewald accelerated Stokesian dynamics method (Wang & Brady 2016), and recently an $O(N)$ fast Stokesian dynamics method (Fiore & Swan 2019). In principle, because the mathematical structure shown in (3.17) remains essentially unchanged between passive and active particles, any of these approaches may be used to simulate active suspensions with minor modification, and we do not suggest any particular numerical implementation here. Indeed, the recent fast Stokesian dynamics method utilizes an imposed Brownian ‘slip’ velocity in order to obtain the stochastic rigid-body motion of passive Brownian particles as shown in previous work (Delmotte & Keaveny 2015; Sprinkle et al. 2017), and we have adapted the fast Stokesian dynamics approach to include particle activity (see figure 5) but will discuss algorithmic details in a future work.

Figure 5. Simulation of a suspension of 4000 identical active particles using the fast Stokesian dynamics method (warmer colours indicate faster speeds).

5. Conclusions

In this work, we have given a detailed exact theoretical description of the dynamics suspensions of active particles in fluids in the absence of inertia, including full hydrodynamic interactions among particles. We argue that, as is done for passive particles, hydrodynamic interactions are ideally separated into near-field and far-field forces, with the latter expanded in a truncated set of moments. The resulting mathematical structure of the dynamical equations remains virtually unchanged between passive and active particles save for the addition of velocity moments due to particle activity. Because of this, any implementation of the Stokesian dynamics method for passive particles may be modified simply and easily for use with active particles. Moreover, we believe that this mathematical structure provides an ideal formalism for theoretical analysis of hydrodynamic interactions in active matter, much as it has for passive suspensions.