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 Reference Schweitzer and Farmer2007; Morozov Reference Morozov2017). The interaction of many of these individual active particles can lead to complex collective dynamics (Ramaswamy Reference Ramaswamy2010). Natural examples include a flock of birds, a school of fish, or a suspension of bacteria (Toner, Tu & Ramaswamy Reference Toner, Tu and Ramaswamy2005), but active matter may also be composed of synthetic active particles (Bechinger et al. Reference Bechinger, Di Leonardo, Löwen, Reichhardt, Volpe and Volpe2016). 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. Reference Marchetti, Joanny, Ramaswamy, Liverpool, Prost, Rao and Simha2013; Saintillan Reference Saintillan2018).
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 Reference Happel and Brenner1965). 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 Reference Davis and Acrivos1985), eventually overcome by the pioneering work of Batchelor (Reference Batchelor1972), Jeffrey (Reference Jeffrey1974), Hinch (Reference Hinch1977), O'Brien (Reference O'Brien1979) 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 Reference Brady and Bossis1988). 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 Reference Durlofsky, Brady and Bossis1987). When the suspension is infinite or periodic, a modification of the method introduced by O'Brien (Reference O'Brien1979) 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. Reference Brady, Phillips, Lester and Bossis1988). 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 (Reference Taylor1951) through to several detailed reviews of microscale locomotion research (Lighthill Reference Lighthill1976; Brennen & Winet Reference Brennen and Winet1977; Lauga & Powers Reference Lauga and Powers2009). 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 Reference Saintillan2018). Interesting phenomena, such as particle clustering (motility-induced phase separation), have been observed for dense suspensions of active particles (Bechinger et al. Reference Bechinger, Di Leonardo, Löwen, Reichhardt, Volpe and Volpe2016), but very often numerical simulation of these suspensions is done with active Brownian particle models that neglect hydrodynamic interactions entirely (Cates & Tailleur Reference Cates and Tailleur2015). Others have used approaches for active suspensions that only approximate the Stokes equations, such as multiparticle collision dynamics (Zöttl & Stark Reference Zöttl and Stark2014) or lattice Boltzmann methods (Stenhammar et al. Reference Stenhammar, Nardini, Nash, Marenduzzo and Morozov2017) that still may not be accurate for very dense concentrations. Results for simplified swimmers in concentrated suspensions display qualitative differences (Ishikawa, Locsei & Pedley Reference Ishikawa, Locsei and Pedley2008; Evans et al. Reference Evans, Ishikawa, Yamaguchi and Lauga2011; Alarcón & Pagonabarraga Reference Alarcón and Pagonabarraga2013; Matas-Navarro et al. Reference Matas-Navarro, Golestanian, Liverpool and Fielding2014; Zöttl & Stark Reference Zöttl and Stark2014; Thutupalli et al. Reference Thutupalli, Geyer, Singh, Adhikari and Stone2018). Some argue that hydrodynamic interactions act to suppress phase separation in active matter (Matas-Navarro et al. Reference Matas-Navarro, Golestanian, Liverpool and Fielding2014), while others have shown that hydrodynamic interactions with boundaries can control phase separation (Thutupalli et al. Reference Thutupalli, Geyer, Singh, Adhikari and Stone2018). 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 Reference Burkholder and Brady2018). We believe that this mathematical structure provides an ideal formalism for theoretical analysis of active suspensions (Burkholder & Brady Reference Burkholder and Brady2019), much as it has for passive suspensions (Brady Reference Brady1993a,Reference Bradyb).
The Stokesian dynamics method was first adapted for use with self-propelled active particles by Mehandia & Nott (Reference Mehandia and Nott2008). 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. (Reference Ishikawa, Locsei and Pedley2008) adapted the Stokesian dynamics framework for use with spherical particles with a prescribed tangential slip velocity, so-called squirmer particles (Ishikawa et al. Reference Ishikawa, Locsei and Pedley2008). Using their own previous results for two-body hydrodynamic interaction between squirmer particles (Ishikawa, Simmonds & Pedley Reference Ishikawa, Simmonds and Pedley2006), Ishikawa et al. (Reference Ishikawa, Locsei and Pedley2008) 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 Reference Ishikawa and Pedley2007b), diffusion (Ishikawa & Pedley Reference Ishikawa and Pedley2007a) and coherent structures (Ishikawa & Pedley Reference Ishikawa and Pedley2008) 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. Reference Swan, Brady and Moore2011). 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. Reference Swan, Brady and Moore2011). 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 Reference Ichiki2002), was developed for suspensions of squirmer particles by using tensorial spherical harmonics (Singh, Ghose & Adhikari Reference Singh, Ghose and Adhikari2015). 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 Reference Maxey and Patel2001) 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 Reference Keaveny2014; Delmotte & Keaveny Reference Delmotte and Keaveny2015) and then to squirmer active particles, up to the stresslet level (Delmotte et al. Reference Delmotte, Keaveny, Plouraboué and Climent2015). Similarly, a fluctuating immersed boundary method (FIBM) was also developed to simulate suspensions of Brownian particles (Delong et al. Reference Delong, Balboa Usabiaga, Delgado-Buscalioni, Griffith and Donev2014). Like the FCM, the FIBM uses an explicit (fluctuating) solvent but uses kernels developed by Peskin (Reference Peskin2002) 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. Reference Balboa Usabiaga, Kallemov, Delmotte, Bhalla, Griffith and Donev2016). 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. Reference Balboa Usabiaga, Kallemov, Delmotte, Bhalla, Griffith and Donev2016; Sprinkle et al. Reference Sprinkle, Balboa Usabiaga, Patankar and Donev2017). 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. (Reference Singh, Ghose and Adhikari2015) and also shown here, a recent approach (for non-Brownian suspensions) uses instead vector spherical harmonics (Corona & Veerapaneni Reference Corona and Veerapaneni2018; Yan et al. Reference Yan, Corona, Malhotra, Veerapaneni and Shelley2020), 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
\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,
\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
\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,
\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 Reference Blake1971).
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
\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
\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:
\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),
\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,
\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,
\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
\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:
\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 Reference Claeys and Brady1993a,Reference Claeys and Bradyb,Reference Claeys and Bradyc; Nasouri & Elfring Reference Nasouri and Elfring2018).
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 Reference Hess2015). In this way, the deformation velocity of each active particle may be written as
\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
\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)}$,
$$\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}$$
$$\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}$$
$$\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. Reference Swan, Brady and Moore2011)
\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
\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):
\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:
\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 Reference Brady and Bossis1988); 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 Reference Anderson and Prieve1991; Stone & Samuel Reference Stone and Samuel1996; Elfring Reference Elfring2015). 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 (Reference Spagnolie and Lauga2010), 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. Reference Ishikawa, Simmonds and Pedley2006; Lauga & Michelin Reference Lauga and Michelin2016; Nasouri & Elfring Reference Nasouri and Elfring2018) of individual spherical active particles, where
\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 Reference Michelin and Lauga2014). This illustrates why it is particularly important to incorporate higher-order moments to capture accurately hydrodynamic interactions between active particles (Singh et al. Reference Singh, Ghose and Adhikari2015). 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. Reference Sanchez, Chen, DeCamp, Heymann and Dogic2012). Immotile active particles, sometimes called ‘shakers’ in contrast to ‘movers’ that are motile (Hatwalne et al. Reference Hatwalne, Ramaswamy, Rao and Simha2004), 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
\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
\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 Reference Happel and Brenner1965). 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 Reference Swan and Brady2007, Reference Swan and Brady2010). 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. Reference Swan, Brady and Moore2011).
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. (Reference Swan, Brady and Moore2011) and Nasouri & Elfring (Reference Nasouri and Elfring2018) 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 Reference Pedley2016). 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 Reference Lighthill1952; Blake Reference Blake1971), or as the outer solution of phoretic flow due to chemical concentrations confined to a thin layer near the sphere surface (Anderson Reference Anderson1989; Golestanian, Liverpool & Ajdari Reference Golestanian, Liverpool and Ajdari2005). The slip velocity is typically written as an expansion in Legendre polynomials:
\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 Reference Pak and Lauga2014). Recasting the slip velocity in terms of irreducible tensors of the surface normal such that
\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
$$\begin{gather} \boldsymbol{U}^s ={-}\tfrac{2}{3}B_1 \boldsymbol{p}, \end{gather}$$
$$\begin{gather}\boldsymbol{\varOmega}^s = \frac{1}{a^3}\,C_1\boldsymbol{p}, \end{gather}$$
$$\begin{gather}a\boldsymbol{E}^s ={-}\tfrac{3}{5}B_2 \overparen{\boldsymbol{p}\boldsymbol{p}}, \end{gather}$$
$$\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}$$
$$\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 Reference Pak and Lauga2014).
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:
$$\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}$$
$$\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
\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. (Reference Swan, Brady and Moore2011), 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
$$\begin{gather} \boldsymbol{U}_{\alpha} = \boldsymbol{U}_A+\boldsymbol{\varOmega}_A\times(\boldsymbol{x}_\alpha-\boldsymbol{x}_A), \end{gather}$$
$$\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. (Reference Swan, Brady and Moore2011), 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
\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:
\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
\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
\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
\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. (Reference Swan, Brady and Moore2011). 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
\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 Reference Purcell1977).
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
\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 Reference Batchelor1972). If the boundary integral equations are also expanded in irreducible moments (Durlofsky et al. Reference Durlofsky, Brady and Bossis1987), then we obtain a linear relationship between force and velocity moments:
\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. Reference Durlofsky, Brady and Bossis1987). 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. Reference Swan, Brady and Moore2011). 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. Reference Durlofsky, Brady and Bossis1987). In this manner, the hydrodynamic forces on the particles are decomposed,
\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 (Reference Sharifi-Mood, Mozaffari and Córdova-Figueroa2016) and Papavassiliou & Alexander (Reference Papavassiliou and Alexander2017) 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. (Reference Ishikawa, Simmonds and Pedley2006) 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
\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. Reference Swan, Brady and Moore2011).
More accurately, the exact two-body swim force contributions may be computed entirely separately, as done by Ishikawa et al. (Reference Ishikawa, Simmonds and Pedley2006), by integrating (3.4) directly. In this way, (3.7) is written as
\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
\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 Reference Batchelor1972). Brady et al. (Reference Brady, Phillips, Lester and Bossis1988) adapted the method of O'Brien (Reference O'Brien1979) 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. (Reference Brady, Phillips, Lester and Bossis1988) 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. Reference Ishikawa, Locsei and Pedley2008). 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. Reference Durlofsky, Brady and Bossis1987; Brady et al. Reference Brady, Phillips, Lester and Bossis1988), 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 Reference Sierou and Brady2001), an 
$O(N^{1.25} \ln N)$ Brownian accelerated Stokesian dynamics method (Banchio & Brady Reference Banchio and Brady2003), an 
$O(N \ln N)$ spectral Ewald accelerated Stokesian dynamics method (Wang & Brady Reference Wang and Brady2016), and recently an 
$O(N)$ fast Stokesian dynamics method (Fiore & Swan Reference Fiore and Swan2019). 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 Reference Delmotte and Keaveny2015; Sprinkle et al. Reference Sprinkle, Balboa Usabiaga, Patankar and Donev2017), 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.
Acknowledgements
G.J.E. acknowledges the hospitality of the Division of Chemistry and Chemical Engineering at the California Institute of Technology during a sabbatical stay that served as a formative period of this work. The authors dedicate this work to the memory of the late Professor J.W. Swan.
Funding
This work was supported the National Science Foundation (J.F.B., grant no. CBET 1803662), the Natural Sciences and Engineering Research Council of Canada (G.J.E., grant no. RGPIN-2020-04850) and by a UBC Killam Research Fellowship to G.J.E.
Declaration of interests
The authors report no conflict of interest.
Appendix A. A reciprocal theorem for active particles
We derive here an equation for the hydrodynamic forces on active particles as shown in (3.3) using the reciprocal theorem (see the excellent review by Masoud & Stone Reference Masoud and Stone2019). The presentation here largely follows that found in our other work (Elfring & Lauga Reference Elfring and Lauga2015; Elfring Reference Elfring2017) but is also found elsewhere (Papavassiliou & Alexander Reference Papavassiliou and Alexander2015).
Consider 
$N$ active free particles 
$\mathcal {B}_i$ with surfaces 
$\partial \mathcal {B}_i$, where 
$i\in [1,N]$, with boundary conditions 
$\boldsymbol {u}(\boldsymbol {x}\in \partial \mathcal {B}_i) = \boldsymbol {U}_i+\boldsymbol {\varOmega }_i\times \boldsymbol {r}_i+\boldsymbol {u}_i^s$, immersed in a background flow 
$\boldsymbol {u}^\infty$ (note that 
$\boldsymbol {u}^\infty$ describes the background flow without the presence of the particle). As an auxiliary problem, here denoted by a hat, consider 
$N$ bodies of the same instantaneous shape undergoing rigid-body motion 
$\hat {\boldsymbol {u}}(\boldsymbol {x}\in \partial \mathcal {B}_i) = \hat {\boldsymbol {U}}_i+\hat {\boldsymbol {\varOmega }}_i\times \boldsymbol {r}_i$ in a quiescent fluid (although not necessary, we take the fluids to have equal viscosity). All flow fields are incompressible, and we neglect inertia in the fluid so that we may write
\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\boldsymbol{\sigma}'\boldsymbol{\cdot}\hat{\boldsymbol{u}}-\hat{\boldsymbol{\sigma}}\boldsymbol{\cdot}\boldsymbol{u}'\right)=\boldsymbol{0}, \end{equation}
where we define disturbance flow 
$\boldsymbol {u}' = \boldsymbol {u} -\boldsymbol {u}^\infty$ and disturbance stress 
$\boldsymbol {\sigma }'=\boldsymbol {\sigma }-\boldsymbol {\sigma }^\infty$
We now integrate over a (sufficiently extended) volume of fluid exterior to 
$\mathcal {B}$ and apply the divergence theorem. Provided that the fields 
$\boldsymbol {u}'$ and 
$\boldsymbol {\sigma }'$ decay appropriately in the far field (Leal Reference Leal1980), we obtain
\begin{equation} \sum_i\int_{\partial\mathcal{B}_i}\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\sigma}'\boldsymbol{\cdot}\hat{\boldsymbol{u}}\,\text{d} S= \sum_i\int_{\partial\mathcal{B}_i}\boldsymbol{n}\boldsymbol{\cdot}\hat{\boldsymbol{\sigma}}\boldsymbol{\cdot}\boldsymbol{u}'\,\text{d} S. \end{equation}
Here, 
$\boldsymbol {n}$ is the normal to the surface 
$\partial \mathcal {B}_i$ pointing into the fluid. This is a statement of the equality of the virtual power of the motions of 
$\partial \mathcal {B}_i$ between 
$\boldsymbol {u}$ and 
$\hat {\boldsymbol {u}}$ (Happel & Brenner Reference Happel and Brenner1965).
Applying the boundary conditions for the rigid-body motion of the particles in the auxiliary problem, we obtain
\begin{equation} \sum_i\left[\boldsymbol{F}_i\boldsymbol{\cdot}\hat{\boldsymbol{U}}_i+\boldsymbol{L}_i\boldsymbol{\cdot}\hat{\boldsymbol{\varOmega}}_i\right] =\sum_i\int_{\partial\mathcal{B}_i}\boldsymbol{n}\boldsymbol{\cdot}\hat{\boldsymbol{\sigma}}\boldsymbol{\cdot}\boldsymbol{u}'\,\text{d} S, \end{equation}where the hydrodynamic force and torque on the active particles are, respectively,
$$\begin{gather} \boldsymbol{F}_i = \int_{\partial\mathcal{B}_i}\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\sigma}'\,\text{d} S, \end{gather}$$
$$\begin{gather}\boldsymbol{L}_i = \int_{\partial\mathcal{B}_i}\boldsymbol{r}_i\times(\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\sigma}')\,\text{d} S, \end{gather}$$where we drop the primes because the background flow is force- and torque-free.
Introducing a more compact notation, where
$$\begin{gather} \boldsymbol{\mathsf{F}} =\left[\boldsymbol{F}_1, \ \boldsymbol{L}_1, \ \boldsymbol{F}_2, \ \boldsymbol{L}_2, \ \ldots\right], \end{gather}$$
$$\begin{gather}\hat{\boldsymbol{\mathsf{U}}} =\left[\hat{\boldsymbol{U}}_1, \ \hat{\boldsymbol{\varOmega}}_1, \ \hat{\boldsymbol{U}}_2, \ \hat{\boldsymbol{\varOmega}}_2, \ \ldots\right], \end{gather}$$we obtain
\begin{equation} \boldsymbol{\mathsf{F}}\boldsymbol{\cdot}\hat{\boldsymbol{\mathsf{U}}}=\sum_i\int_{\partial\mathcal{B}_i}\boldsymbol{n}\boldsymbol{\cdot}\hat{\boldsymbol{\sigma}}\boldsymbol{\cdot}\boldsymbol{u}'\,\text{d} S.\end{equation}
Now, by linearity, we may write 
$\hat {\boldsymbol {\sigma }} = \boldsymbol{\mathsf{T}}_{\boldsymbol{\mathsf{U}}}\boldsymbol {\cdot }\hat {\boldsymbol{\mathsf{U}}}$, and substitution into (A8), upon discarding the arbitrary 
$\hat {\boldsymbol{\mathsf{U}}}$, leads to (3.3) for the hydrodynamic forces and torques on all 
$N$ active particles:
\begin{equation} \boldsymbol{\mathsf{F}}=\sum_i\int_{\partial\mathcal{B}}\boldsymbol{u}_i'\boldsymbol{\cdot}(\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\mathsf{T}}_{\boldsymbol{\mathsf{U}}})\,\text{d} S. \end{equation}This derivation extends naturally to higher-order force moments by taking the auxiliary problem to be rigid-body motion in an arbitrary background flow, represented as a series expansion (Elfring Reference Elfring2017; Nasouri & Elfring Reference Nasouri and Elfring2018).
Appendix B. Tensorial spherical harmonics
The (
$n$-adic) tensorial spherical harmonics are a set of tensors composed of irreducible products of the unit normal on a sphere (Brenner Reference Brenner1964a; Hess Reference Hess2015):
$$\begin{gather} \overparen{\boldsymbol{n}^2}=\overparen{\boldsymbol{n}\boldsymbol{n}}=\boldsymbol{n}\boldsymbol{n}-\tfrac{1}{3}\boldsymbol{I}, \end{gather}$$
$$\begin{gather}\overparen{\boldsymbol{n}^3}=\overparen{\boldsymbol{n}\boldsymbol{n}\boldsymbol{n}} = \boldsymbol{n}\boldsymbol{n}\boldsymbol{n} - \tfrac{1}{5}\left(\boldsymbol{I}\boldsymbol{n}+\boldsymbol{n}\boldsymbol{I}+ {}^{\rm T}(\boldsymbol{n}\boldsymbol{I})]\right), \end{gather}$$
\begin{gather} \left[\overparen{\boldsymbol{n}^4}\right]_{ijkl} = n_in_jn_kn_l - \tfrac{1}{7}\left(n_in_j\delta_{kl}+n_i\delta_{jk}n_l+\delta_{ij}n_kn_l+ \delta_{il}n_jn_k+\delta_{ik}n_jn_l+n_i\delta_{jl}n_k\right),\nonumber\\ \hspace{-10.5pc}+\tfrac{1}{35}\left(\delta_{ij}\delta_{kl}+\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}\right), \end{gather}
\begin{gather} \overparen{\boldsymbol{n}^n}=\frac{({-}1)^n}{(2n-1)!!}\,a^{n+1}\,\boldsymbol{\nabla}^{n}\left(\frac{1}{r}\right)_{r=a}, \end{gather}
where we use notation similar to Brenner (Reference Brenner1964b) such that 
$\boldsymbol {\nabla }^2 \equiv \boldsymbol {\nabla }\boldsymbol {\nabla }$ (distinct from 
$\nabla ^2=\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {\nabla }$), the overbracket indicates the irreducible (or fully symmetric and traceless) part of the tensor, and such that the transpose applies to the two adjacent indices (for example 
${}^{\rm T}\boldsymbol {a}\boldsymbol {b}\boldsymbol {c}=\boldsymbol {b}\boldsymbol {a}\boldsymbol {c}$ and 
$\boldsymbol {a}\boldsymbol {b}\boldsymbol {c}^{\rm T}=\boldsymbol {a}\boldsymbol {c}\boldsymbol {b}$).
The tensorial spherical harmonics are orthogonal, with the relationship
\begin{equation} \frac{1}{4{\rm \pi} a^2}\int w_n \overparen{\boldsymbol{n}^m} \overparen{\boldsymbol{n}^n} \,\text{d} S = \delta_{mn}\boldsymbol{\varDelta}^n, \end{equation}where the weight is
\begin{equation} w_n = \frac{(2n+1)!!}{n!}. \end{equation}
The isotropic tensor 
$\boldsymbol {\varDelta }^n$ is a 
$2n$-order tensor that projects an 
$n$-order tensor into its symmetric irreducible form (Hess Reference Hess2015); i.e. for the 
$n$-order tensor 
$\boldsymbol {A}$, 
$\boldsymbol {\varDelta }^{n}\odot \boldsymbol {A}=\overparen {\boldsymbol {A}}$, where 
$\odot$ is a complete tensor contraction. The first several symmetrizing tensors are
$$\begin{gather} \varDelta^0=1, \end{gather}$$
$$\begin{gather}\varDelta^1_{ii'}=\delta_{ii'}, \end{gather}$$
$$\begin{gather}\varDelta^2_{iji'j'}=\tfrac{1}{2}\left(\delta_{ii'}\delta_{jj'}+\delta_{ij'}\delta_{ji'}\right)- \tfrac{1}{3}\delta_{ij}\delta_{i'j'}, \end{gather}$$
\begin{align} \varDelta^3_{ijki'j'k'}&=\tfrac{1}{6}\left(\delta_{ii'}\delta_{jj'}\delta_{kk'}+ \delta_{ii'}\delta_{jk'}\delta_{kj'}+\delta_{ij'}\delta_{ji'}\delta_{kk'}+ \delta_{ij'}\delta_{jk'}\delta_{ki'}+\delta_{ik'}\delta_{ji'}\delta_{kj'}+ \delta_{ik'}\delta_{jj'}\delta_{ki'}\right)\nonumber\\ &\quad -\tfrac{1}{15}\left\{\left(\delta_{i'j'}\delta_{kk'}+\delta_{i'k'}\delta_{kj'}+ \delta_{j'k'}\delta_{ki'}\right)\delta_{ij}+\left(\delta_{i'j'}\delta_{ik'}+\delta_{i'k'} \delta_{ij'}+\delta_{j'k'}\delta_{ii'}\right)\delta_{jk}\right.\nonumber\\ &\quad\left. {}+\left(\delta_{i'j'}\delta_{jk'}+\delta_{i'k'}\delta_{jj'}+ \delta_{j'k'}\delta_{ji'}\right)\delta_{ik}\right\}, \end{align}where the primed indices are distinct from the unprimed ones.
Appendix C. Expansion of the boundary integral equation
We derive here the grand mobility relationship between velocity moments and force moments by means of a Galerkin projection onto tensorial spherical harmonics (Singh et al. Reference Singh, Ghose and Adhikari2015). Consider the boundary integral equation for a suspension of active particles
\begin{equation} \boldsymbol{u}(\boldsymbol{x}) -\boldsymbol{u}^\infty(\boldsymbol{x}) ={-}\sum_i\int_{\partial\mathcal{B}_i}\left[\boldsymbol{G}(\boldsymbol{x},\boldsymbol{y})\boldsymbol{\cdot} \boldsymbol{f}(\boldsymbol{y})+\boldsymbol{u}(\boldsymbol{y})\boldsymbol{\cdot}\boldsymbol{T}(\boldsymbol{x},\boldsymbol{y})\boldsymbol{\cdot}\boldsymbol{n}(\boldsymbol{y})\right]\,\text{d} S(\boldsymbol{y}), \end{equation}where
$$\begin{gather} \boldsymbol{G}(\boldsymbol{x},\boldsymbol{y}) = \frac{1}{8{\rm \pi}\eta}\left(\frac{\boldsymbol{I}}{\left|\boldsymbol{x}-\boldsymbol{y}\right|}+ \frac{(\boldsymbol{x}-\boldsymbol{y})(\boldsymbol{x}-\boldsymbol{y})}{\left|\boldsymbol{x}-\boldsymbol{y}\right|^3}\right), \end{gather}$$
$$\begin{gather}\boldsymbol{T}(\boldsymbol{x},\boldsymbol{y}) ={-}\frac{3}{4{\rm \pi}}\,\frac{(\boldsymbol{x}-\boldsymbol{y})(\boldsymbol{x}-\boldsymbol{y})(\boldsymbol{x}-\boldsymbol{y})}{\left|\boldsymbol{x}-\boldsymbol{y}\right|^5}. \end{gather}$$The velocities and tractions on each particle are now expressed in terms of expansions in tensorial spherical harmonics:
$$\begin{gather} \boldsymbol{u}(\boldsymbol{y} \in \partial \mathcal{B}_i)=\boldsymbol{U}_i+\boldsymbol{U}^s_i+(\boldsymbol{\varOmega}_i+\boldsymbol{\varOmega}_i^s)\times\boldsymbol{r}_i+\boldsymbol{r}_i\boldsymbol{\cdot}\boldsymbol{E}_i^s+\cdots, \end{gather}$$
$$\begin{gather}\boldsymbol{f}(\boldsymbol{y}\in\partial\mathcal{B}_i) = \frac{1}{4{\rm \pi} a_i^2}\,\boldsymbol{F}_i+\frac{3}{8{\rm \pi} a_i^3}\,\boldsymbol{L}_i\times\boldsymbol{n}+\frac{3}{4{\rm \pi} a_i^3}\,\boldsymbol{n}\boldsymbol{\cdot}\tilde{\boldsymbol{S}}_i+\cdots, \end{gather}$$where
$$\begin{gather} \boldsymbol{F}_i =\int_{\partial\mathcal{B}_i}\boldsymbol{f}\,\text{d} S, \end{gather}$$
$$\begin{gather}\boldsymbol{L}_i = \int_{\partial\mathcal{B}_i}\boldsymbol{r}_i\times\boldsymbol{f}\,\text{d} S, \end{gather}$$
$$\begin{gather}\tilde{\boldsymbol{S}}_i = \int_{\partial\mathcal{B}_i} \overparen{\boldsymbol{r}_i\boldsymbol{f}} \,\text{d} S. \end{gather}$$ Now taking moments of the flow over the surface of particle 
$\alpha$, 
$\boldsymbol {u}(\boldsymbol {x}\in \partial \mathcal {B}_\alpha )$, we systematically obtain mobility relationships for the 
$\alpha$ particle:
$$\begin{gather} \boldsymbol{U}_\alpha+\boldsymbol{U}_\alpha^s-\boldsymbol{U}_\alpha^\infty ={-}\boldsymbol{M}^{\alpha\alpha}_{UF}\boldsymbol{\cdot} \boldsymbol{F}_\alpha-\sum_{\beta\ne\alpha}[\boldsymbol{M}^{\alpha\beta}_{UF}\boldsymbol{\cdot}\boldsymbol{F}_\beta+ \boldsymbol{M}^{\alpha\beta}_{UL}\boldsymbol{\cdot}\boldsymbol{L}_\beta+\boldsymbol{M}^{\alpha\beta}_{US}:\boldsymbol{S}_\beta+\cdots], \end{gather}$$
$$\begin{gather}\boldsymbol{\varOmega}_\alpha+\boldsymbol{\varOmega}_\alpha^s-\boldsymbol{\varOmega}_\alpha^\infty ={-}\boldsymbol{M}^{\alpha\alpha}_{\varOmega L}\boldsymbol{\cdot} \boldsymbol{L}_\alpha-\sum_{\beta\ne\alpha}[\boldsymbol{M}^{\alpha\beta}_{\varOmega F}\boldsymbol{\cdot}\boldsymbol{F}_\beta+ \boldsymbol{M}^{\alpha\beta}_{\varOmega L}\boldsymbol{\cdot}\boldsymbol{L}_\beta+\boldsymbol{M}^{\alpha\beta}_{\varOmega S}:\boldsymbol{S}_\beta+\cdots], \end{gather}$$
$$\begin{gather}\boldsymbol{E}_\alpha^s-\boldsymbol{E}_\alpha^\infty ={-}\boldsymbol{M}^{\alpha\alpha}_{ES}\boldsymbol{\cdot}\boldsymbol{S}_\alpha-\sum_{\beta\ne\alpha} [\boldsymbol{M}^{\alpha\beta}_{EF}\boldsymbol{\cdot}\boldsymbol{F}_\beta+\boldsymbol{M}^{\alpha\beta}_{EL}\boldsymbol{\cdot} \boldsymbol{L}_\beta+\boldsymbol{M}^{\alpha\beta}_{ES}:\boldsymbol{S}_\beta+\cdots], \end{gather}$$where the stresslet for active particles includes a contribution from the double-layer kernel
\begin{equation} \boldsymbol{S}_\beta = \tilde{\boldsymbol{S}}_\beta -2\eta\,\frac{4{\rm \pi} a_\beta^3}{3}\,\boldsymbol{E}^s_\beta. \end{equation}The mobility tensors are identical to those for passive particles:
$$\begin{gather} \boldsymbol{M}^{\alpha\alpha}_{UF}=\frac{1}{6{\rm \pi} a_\alpha \eta}\,\boldsymbol{I}, \end{gather}$$
$$\begin{gather}\boldsymbol{M}^{\alpha\beta}_{UF}=\frac{1}{4{\rm \pi} a_\alpha^2}\int_{\partial\mathcal{B}_\alpha}\,\text{d} S(\boldsymbol{x})\, \frac{1}{4{\rm \pi} a_\beta^2}\int_{\partial\mathcal{B}_\beta}\boldsymbol{G}(\boldsymbol{x},\boldsymbol{y})\,\text{d} S(\boldsymbol{y}), \end{gather}$$
$$\begin{gather}\boldsymbol{M}^{\alpha\beta}_{UL}=\frac{1}{4{\rm \pi} a_\alpha^2}\int_{\partial\mathcal{B}_\alpha}\,\text{d} S(\boldsymbol{x})\, \frac{3}{8{\rm \pi} a_\beta^3}\int_{\partial\mathcal{B}_\beta}\boldsymbol{G}(\boldsymbol{x},\boldsymbol{y})\times\boldsymbol{n}(\boldsymbol{y})\,\text{d} S(\boldsymbol{y}), \end{gather}$$
$$\begin{gather}\boldsymbol{M}^{\alpha\beta}_{US}=\frac{1}{4{\rm \pi} a_\alpha^2}\int_{\partial\mathcal{B}_\alpha}\,\text{d} S(\boldsymbol{x})\, \frac{3}{8{\rm \pi} a_\beta^3}\int_{\partial\mathcal{B}_\beta}\boldsymbol{G}(\boldsymbol{x},\boldsymbol{y})\,(\boldsymbol{n}(\boldsymbol{y})+\boldsymbol{n}(\boldsymbol{y})^{\rm T})\,\text{d} S(\boldsymbol{y}), \end{gather}$$
$$\begin{gather}\boldsymbol{M}_{\varOmega L}^{\alpha\alpha}=\frac{1}{8{\rm \pi}\eta a_\alpha^3}\,\boldsymbol{I}, \end{gather}$$
$$\begin{gather}\boldsymbol{M}_{\varOmega F}^{\alpha\beta}=\frac{3}{8{\rm \pi} a_\alpha^3}\int_{\partial\mathcal{B}_\alpha}\,\text{d} S(\boldsymbol{x})\, \frac{1}{4{\rm \pi} a_\beta^2}\int_{\partial\mathcal{B}_\beta}\,\text{d} S(\boldsymbol{y})\,\boldsymbol{n}(\boldsymbol{x})\times\boldsymbol{G}(\boldsymbol{x},\boldsymbol{y}), \end{gather}$$
$$\begin{gather}\boldsymbol{M}_{\varOmega L}^{\alpha\beta}=\frac{3}{8{\rm \pi} a_\alpha^3}\int_{\partial\mathcal{B}_\alpha}\,\text{d} S(\boldsymbol{x})\, \frac{3}{8{\rm \pi} a_\beta^3}\int_{\partial\mathcal{B}_\beta}\,\text{d} S(\boldsymbol{y})\,\boldsymbol{n}(\boldsymbol{x})\times\boldsymbol{G}(\boldsymbol{x},\boldsymbol{y})\times\boldsymbol{n}(\boldsymbol{y}), \end{gather}$$
$$\begin{gather}\boldsymbol{M}_{\varOmega S}^{\alpha\beta}=\frac{3}{8{\rm \pi} a_\alpha^3}\int_{\partial\mathcal{B}_\alpha}\,\text{d} S(\boldsymbol{x})\, \frac{3}{8{\rm \pi} a_\beta^3}\int_{\partial\mathcal{B}_\beta}\,\text{d} S(\boldsymbol{y})\left(\boldsymbol{n}(\boldsymbol{x})\times\boldsymbol{G}(\boldsymbol{x},\boldsymbol{y})(\boldsymbol{n}(\boldsymbol{y})+\boldsymbol{n}(\boldsymbol{y})^{\rm T}\right), \end{gather}$$
$$\begin{gather}\boldsymbol{M}^{\alpha\alpha}_{ES}=\frac{3}{20{\rm \pi}\eta a_\alpha^3}\,\mathbb{I}, \end{gather}$$
$$\begin{gather}\boldsymbol{M}_{E F}^{\alpha\beta}=\frac{3}{8{\rm \pi} a_\alpha^3}\int_{\partial\mathcal{B}_\alpha}\,\text{d} S(\boldsymbol{x})\, \frac{1}{4{\rm \pi} a_\beta^2}\int_{\partial\mathcal{B}_\beta}\,\text{d} S(\boldsymbol{y})\left(\boldsymbol{n}(\boldsymbol{x})+{}^{\rm T}\boldsymbol{n}(\boldsymbol{x})\right)\boldsymbol{G}(\boldsymbol{x},\boldsymbol{y}), \end{gather}$$
$$\begin{gather}\boldsymbol{M}_{E L}^{\alpha\beta}=\frac{3}{8{\rm \pi} a_\alpha^3}\int_{\partial\mathcal{B}_\alpha}\,\text{d} S(\boldsymbol{x})\, \frac{3}{8{\rm \pi} a_\beta^3}\int_{\partial\mathcal{B}_\beta}\,\text{d} S(\boldsymbol{y})\left(\boldsymbol{n}(\boldsymbol{x})+{}^{\rm T}\boldsymbol{n}(\boldsymbol{x})\right)\boldsymbol{G}(\boldsymbol{x},\boldsymbol{y})\times\boldsymbol{n}(\boldsymbol{y}), \end{gather}$$
$$\begin{gather}\boldsymbol{M}_{E S}^{\alpha\beta}=\frac{3}{8{\rm \pi} a_\alpha^3}\int_{\partial\mathcal{B}_\alpha}\,\text{d} S(\boldsymbol{x})\,\frac{3}{8{\rm \pi} a_\beta^3}\int_{\partial\mathcal{B}_\beta}\,\text{d} S(\boldsymbol{y})\left((\boldsymbol{n}(\boldsymbol{x})+{}^{\rm T}\boldsymbol{n}(\boldsymbol{x})\right)\boldsymbol{G}(\boldsymbol{x},\boldsymbol{y})\left(\boldsymbol{n}(\boldsymbol{y})+\boldsymbol{n}(\boldsymbol{y})^{\rm T}\right), \end{gather}$$
where 
$\mathbb {I}$ is the fourth-order identity tensor. We give the mobilities here in integral form (Ichiki Reference Ichiki2002; Fiore & Swan Reference Fiore and Swan2018), but it is much more common to see them in the equivalent differential form, which may be found by Taylor expansion about the particle centres (Wajnryb et al. Reference Wajnryb, Mizerski, Zuk and Szymczak2013; Mizerski et al. Reference Mizerski, Wajnryb, Zuk and Szymczak2014; Fiore et al. Reference Fiore, Balboa Usabiaga, Donev and Swan2017).
For all 
$N$ particles, we write the mobility relationships between velocity moments and force moments in compact form as
$$\begin{gather} \boldsymbol{\mathsf{U}}+\boldsymbol{\mathsf{U}}^s-\boldsymbol{\mathsf{U}}^\infty ={-}\boldsymbol{\mathsf{M}}_{\boldsymbol{\mathsf{UF}}}\boldsymbol{\cdot} \boldsymbol{\mathsf{F}}-\boldsymbol{\mathsf{M}}_{\boldsymbol{\mathsf{US}}}:\boldsymbol{\mathsf{S}}+\cdots, \end{gather}$$
$$\begin{gather}\boldsymbol{\mathsf{E}}^s-\boldsymbol{\mathsf{E}}^\infty ={-}\boldsymbol{\mathsf{M}}_{\boldsymbol{\mathsf{EF}}}\boldsymbol{\cdot} \boldsymbol{\mathsf{F}}-\boldsymbol{\mathsf{M}}_{\boldsymbol{\mathsf{ES}}}:\boldsymbol{\mathsf{S}}+\cdots,\\ \vdots \nonumber \end{gather}$$
In this way, we form the grand mobility tensor (typically truncated at the 
$\boldsymbol{\mathsf{E}},\boldsymbol{\mathsf{S}}$ level shown above in Stokesian dynamics)
\begin{equation} \mathcal{U}' ={-}\mathcal{M}\boldsymbol{\cdot}\mathcal{F}. \end{equation}Appendix D. Dilute approximation
Dilute approximations are typically used throughout the literature in order to avoid the computational expense of inverting the far-field grand mobility tensor. As shown above in (C25), the far-field contribution to the swimming dynamics may be written in terms of mobilities. The stresslet is not prescribed but induced, so solving (C26) for 
$\boldsymbol{\mathsf{S}}$ and substituting into (C25) yields
\begin{equation} \boldsymbol{\mathsf{U}} = \boldsymbol{\mathsf{U}}^\infty-\boldsymbol{\mathsf{U}}^s -(\boldsymbol{\mathsf{M}}_{\boldsymbol{\mathsf{UF}}}- \boldsymbol{\mathsf{M}}_{\boldsymbol{\mathsf{US}}}:\boldsymbol{\mathsf{M}}_{\boldsymbol{\mathsf{ES}}}^{{-}1}: \boldsymbol{\mathsf{M}}_{\boldsymbol{\mathsf{EF}}})\boldsymbol{\cdot}\boldsymbol{\mathsf{F}}+\boldsymbol{\mathsf{M}}_{\boldsymbol{\mathsf{US}}}: \boldsymbol{\mathsf{M}}_{\boldsymbol{\mathsf{ES}}}^{{-}1}:(\boldsymbol{\mathsf{E}}^s-\boldsymbol{\mathsf{E}}^\infty)+\cdots, \end{equation}and for force-free particles, we have simply
\begin{equation} \boldsymbol{\mathsf{U}} = \boldsymbol{\mathsf{U}}^\infty-\boldsymbol{\mathsf{U}}^s+\boldsymbol{\mathsf{M}}_{\boldsymbol{\mathsf{US}}}: \boldsymbol{\mathsf{M}}_{\boldsymbol{\mathsf{ES}}}^{{-}1}:(\boldsymbol{\mathsf{E}}^s-\boldsymbol{\mathsf{E}}^\infty)+\cdots. \end{equation}
To leading order in a dilute approximation, only the single particle stresslet terms remain, 
$\boldsymbol{\mathsf{M}}_{\boldsymbol{\mathsf{ES}}}^{-1}$ is diagonal, and we have
$$\begin{gather} \boldsymbol{U}_\alpha =\boldsymbol{U}_\alpha^\infty-\boldsymbol{U}_\alpha^s -\sum_{\beta\ne\alpha}\boldsymbol{M}^{\alpha\beta}_{US}:\boldsymbol{S}_\beta+\cdots, \end{gather}$$
$$\begin{gather}\boldsymbol{\varOmega}_\alpha =\boldsymbol{\varOmega}_\alpha^\infty-\boldsymbol{\varOmega}_\alpha^s -\sum_{\beta\ne\alpha}\boldsymbol{M}^{\alpha\beta}_{\varOmega S}:\boldsymbol{S}_\beta+\cdots, \end{gather}$$where
\begin{equation} \boldsymbol{S}_\beta=(\boldsymbol{M}^{\beta\beta}_{ES})^{{-}1}:\left(\boldsymbol{E}_\beta^\infty-\boldsymbol{E}_\beta^s\right)= \frac{20{\rm \pi}\eta a_\beta^3}{3}\left(\boldsymbol{E}_\beta^\infty-\boldsymbol{E}_\beta^s\right). \end{equation}
Open access
