Higher-order force moments of active particles
Babak Nasouri Gwynn J. Elfring
physics.flu-dyn



Department of Mechanical Engineering, Institute of Applied Mathematics, University of British Columbia, Vancouver, British Columbia, Canada V6T 1Z4

gelfring@mech.ubc.ca
Department of Mechanical Engineering, Institute of Applied Mathematics, University of British Columbia, Vancouver, British Columbia, Canada V6T 1Z4






 
Active particles moving through fluids generate disturbance flows due to their activity. For simplicity, the induced flow field is often modeled by the leading terms in a far-field approximation of the Stokes equations, whose coefficients are the force, torque and stresslet (zeroth and first-order force moments) of the active particle. This level of approximation is quite useful, but may also fail to predict more complex behaviors that are observed experimentally. In this study, to provide a better approximation, we evaluate the contribution of the second-order force moments to the flow field and, by reciprocal theorem, present explicit formulas for the stresslet dipole, rotlet dipole and potential dipole for an arbitrarily-shaped active particle. As examples of this method, we derive modified Faxén laws for active spherical particles and resolve higher-order moments for active rod-like particles.






Higher-order force moments of active particles
    Gwynn J. Elfring
    December 30, 2023
==============================================



§ INTRODUCTION

Self-propulsion is ubiquitous in nature. Be it at the macroscopic scale of flying birds or the microscopic scale of swimming bacteria, the motion of active matter results from converting internal or ambient energy into mechanical work without any external input <cit.>. At sufficiently small scales in viscous fluids, inertia is irrelevant and viscous dissipation dominates the motion of the fluid and active particles within them <cit.>. In the absence of inertia, `reciprocal' body distortions are ineffective as a propulsion mechanism, and so active particles must propel themselves by other means in this realm <cit.>. There exist several techniques to achieve net locomotion in the low-Reynolds-number regime <cit.>. For instance, microorganisms such as Paramecium and Volvox use small appendages called cilia to facilitate motion <cit.>. Cilia generate thrust through a coordinated pattern of beating, which may arise from hydrodynamic <cit.> or basal <cit.> interactions. Propulsion can also be achieved synthetically by chemically-active particles with asymmetric non-uniform surface properties <cit.>. In both of these examples, the effect of surface activity is confined to a narrow region surrounding the particle and hence may be modelled using `apparent' slip velocities on the surface. This way, one can explicitly find the propulsion speed and thereby the disturbance flow field, in terms of prescribed (or measured) slip velocities <cit.>.


For the inertialess motion of sufficiently small particles in viscous fluids, the flow field is often approximated by far-field singularity solutions of the Stokes equations. To leading order, the flow field decays linearly by distance (∼1/|𝐱|) and, at this level of approximation, the particle is replaced by a point force (i.e., zeroth-order force moment) that leads to flow <cit.>. The next-order correction to the flow field, which decays quadratically (∼1/|𝐱|^2), can be expressed using a force-dipole (i.e., first-order force moment), which is decomposed into a torque (the antisymmetric part) and a stresslet (the symmetric part) <cit.>. In the absence of an external force, the over-damped motion of the particle has no net hydrodynamic force or torque and so the stresslet governs the leading-order flow field. The importance of the stresslet in characterizing the interactions of active particles <cit.>, the rheology <cit.> and stability <cit.> of active suspensions, and the collective locomotion of bacteria <cit.> is well documented. However, the stresslet term alone fails to explain behaviors such as the `dancing' of two Volvox colonies when they are in proximity of one another <cit.>, or the vortices induced due to the motion of C. reinhardtii <cit.>. An emerging picture is that modeling the motion using only terms up to the stresslet may limit understanding of how active particles interact with their environment and motivates investigation of higher-order force moments. In a recent study, <cit.> showed that the swirling motion of an active spherical particle only appears in the flow field decaying as ∼1/|𝐱|^3 and ∼1/|𝐱|^4 and derived expressions for higher-order force moments of a sphere. In this study, we generalize their results by investigating the effects of higher-order force moments on an arbitrarily-shaped active particle, thereby extending recent general results for the stresslet term by <cit.>. Using the boundary integral equations, we express the flow field around an active particle through a multipole expansion up to the contribution of the second-order force moments. We then provide explicit formulas for these force moments by exploiting the reciprocal theorem using a framework developed in <cit.>.


The reciprocal theorem for low-Reynolds-number hydrodynamics has long been an avenue to simplify calculations in Stokes flow <cit.>. Its application has ranged from the inertialess jet propulsion <cit.>, boundary-driven channel flow <cit.>, to Marangoni motion of a droplet covered with bulk-insoluble surfactants in a Poiseuille flow <cit.>. In particular, <cit.> showed that the kinematics of an active particle can be determined explicitly, using the flow field induced by the rigid-body motion of a passive particle of the same instantaneous shape. Subsequently, the reciprocal theorem has been widely used to determine the kinematics of active particles both in Newtonian <cit.> and non-Newtonian fluids <cit.>. This approach was recently extended to determine the stresslet of active particles <cit.>. More recently, a general framework has been developed for finding the force moments (of any order) of an active particle in a Newtonian (or non-Newtonian) fluid <cit.>. Following that approach, in this study, we provide formulas for calculating the force moments up to the second order, for any arbitrarily-shaped active particle. 


The paper is organized as follows. In Sec. <ref>, we employ the boundary integral equation to describe the disturbance flow field caused by an active particle. Using an asymptotic expansion of the far-field flow, we show how the force moments contribute to the disturbance flow field. We then in Sec. <ref>, use the reciprocal theorem to find general expressions for these force moments and, as examples, evaluate them explicitly for a spherical active particle, a generalized squirmer and an active slender rod.





§ MULTIPOLE EXPANSION
 

We consider a particle with boundary ∂ℬ in an otherwise unbounded Newtonian fluid of viscosity μ and background flow field 𝐮^∞, as shown in Fig. <ref>. Using the boundary integral equations, the disturbance flow field 𝐮^'=𝐮-𝐮^∞ can be expressed as a summation of single-layer and double-layer potentials <cit.>

    𝐮^'(𝐱)=-1/8πμ∫_∂ℬ𝐟^'(𝐲)·𝐉(𝐱-𝐲) dS(𝐲)-1/8π∫_∂ℬ𝐮^'(𝐲)𝐧(𝐲):𝐊(𝐱-𝐲) dS(𝐲),

where 𝐟^'=𝐧·σ^' is the traction of disturbance stress tensor σ^', 𝐲 is the position that is integrated over the particle surface and 𝐧 is the surface normal pointing into the fluid. Here 𝐉(𝐱)=𝐈/|𝐱|+𝐱𝐱/|𝐱|^3 is the Green's function of Stokes equations (or the Oseen tensor) and 𝐊(𝐱)=-6𝐱𝐱𝐱/|𝐱|^5 is its associated stress tensor. 

Expanding in 𝐲 about a convenient point in the body 𝐱_0 (for example the center of mass),

    𝐉(𝐱-𝐲)   =𝐉(𝐱-𝐱_0) - 𝐫·∇𝐉(𝐱-𝐱_0) +1/2𝐫𝐫: ∇∇𝐉(𝐱-𝐱_0) +⋯,
    𝐊(𝐱-𝐲)   =𝐊(𝐱-𝐱_0) - 𝐫·∇𝐊(𝐱-𝐱_0)+⋯,

where 𝐫 = 𝐲-𝐱_0. Equation (<ref>) then takes the form

    𝐮^'(𝐱)=   -1/8πμ[ <𝐟^'>·𝐉(𝐱-𝐱_0) - <𝐟^'𝐫>:∇𝐉(𝐱-𝐱_0) +1/2<𝐟^'𝐫𝐫>⊙∇∇𝐉(𝐱-𝐱_0)+⋯]
       -1/8π[<𝐮^'𝐧>:𝐊(𝐱-𝐱_0) - <𝐮^'𝐧𝐫>⊙∇𝐊(𝐱-𝐱_0)+⋯].

For convenience, in this work, we denote the surface integral by ∫_∂ℬ⋯d S≡<⋯>.We also use ⊙ to denote a k-fold contraction where k=min{a,b} and a and b are the tensorial orders of the contracted tensors, e.g.,  [ <𝐮^'𝐧𝐫>⊙∇𝐊]_i=<u^'_jn_k𝐫_m>∇_m K_kji. We define 𝐮^'(𝐱)=∑_i=1𝐮^'(i)(𝐱), where 𝐮^'(i)(𝐱) is the flow field that decays as |𝐱|^-i. At leading order, one can recognize the net hydrodynamic force 𝐅=<𝐟^'> as the zeroth-order force moment. The presence of particle at this order is represented by a point force of strength -𝐅 and the flow field is simply governed by a Stokeslet, 

    𝐮^'(1)(𝐱)=-1/8πμ𝐅·𝐉(𝐱-𝐱_0).


To find the flow field decaying as |𝐱|^-2, it is useful to decompose 𝐟^'𝐫 and 𝐮^'𝐧, to their symmetric and antisymmetric parts as 𝐟^'𝐫=𝐟^'𝐫+𝐫𝐟^'/2+𝐟^'𝐫-𝐫𝐟^'/2 and 𝐮^'𝐧=𝐮^'𝐧+𝐧𝐮^'/2+𝐮^'𝐧-𝐧𝐮^'/2. Noting the symmetry of 𝐊(𝐱) and also recalling

    ∇_k J_ij=x_i δ_jk + x_j δ_ik - x_k δ_ij/(x_l x_l)^3/2-3x_i x_j x_k/(x_l x_l)^5/2,

we have

    𝐮^'(2)(𝐱)   =-1/8πμ[<𝐫×𝐟^'>·𝐂(𝐱-𝐱_0)-1/2<𝐟^'𝐫+𝐫𝐟^'/2-1/3(𝐟^'·𝐫)𝐈-μ(𝐮^'𝐧+𝐧𝐮^')>:𝐊(𝐱-𝐱_0)],
       =-1/8πμ[𝐋·𝐂(𝐱-𝐱_0) -1/2𝐒:𝐊(𝐱-𝐱_0)],

where 𝐋=<𝐫×𝐟^'> is the antisymmetric first-order force moment, i.e. torque and 𝐂(𝐱)=𝐈×𝐱/|𝐱|^3 is the associated rotlet (or couplet) tensor <cit.>. In using the cross product, we follow the convention [𝐫×𝐟^']_i=ϵ_ijkr_jf^'_k and [𝐈×𝐫]_ij=ϵ_jskδ_isr_k, where ϵ is the third-order permutation tensor. The symmetric and deviatoric first-order force moment (i.e. 𝐟^'𝐫) along with the contribution of the double-layer potential lead to 𝐒=<𝐟^'𝐫-2μ𝐮^'𝐧>, namely the stresslet <cit.>. The over-bracket denotes the fully-symmetric and deviatoric part of a tensor which are defined

    [ ]_ij   =(1/2)([ ]_ij+[ ]_ji)-(1/3)[ ]_ssδ_ij,
    [ ]_ijk   =(1/6)([ ]_ijk+[ ]_ikj+[ ]_jik+[ ]_jki+[ ]_kij+[ ]_kji)
       -(1/15){([ ]_ssi+[ ]_sis+[ ]_iss)δ_kj+([ ]_ssj+[ ]_sjs+[ ]_jss)δ_ik+([ ]_ssk+[ ]_sks+[ ]_kss)δ_ij},

for the second and third-order tensors, respectively. 

To determine the flow field decaying as |𝐱|^-3, we decompose the third-order tensors 𝐟^'𝐫𝐫 and 𝐮^'𝐧𝐫 to their irreducible parts (see <cit.> for the decomposition technique). Now by taking the second gradient of the Oseen tensor

    ∇_m∇_k J_ij   =δ_imδ_jk+δ_jmδ_ik-δ_kmδ_ij/(x_lx_l)^3/2+15x_ix_jx_kx_m/(x_lx_l)^7/2   -3(x_m x_iδ_jk+x_j x_mδ_ik-x_k x_mδ_ij+x_j x_kδ_im+x_i x_kδ_jm+x_i x_jδ_km)/(x_lx_l)^5/2,

we obtain the next order correction for the flow field

    𝐮^'(3)(𝐱)=   -1/32πμ[   <𝐟^'𝐫𝐫-4μ𝐮^'𝐧𝐫>⊙∇𝐊(𝐱-𝐱_0)        ]
    -1/24πμ[<𝐫(𝐫×𝐟^')>:∇[∇×𝐉(𝐱-𝐱_0)] ]
       -1/80πμ[<2|𝐫|^2𝐟^'-(𝐫·𝐟^')𝐫 +3μ[4(𝐮^'·𝐧)𝐫 - (𝐮^'·𝐫)𝐧-(𝐫·𝐧)𝐮^']>·∇^2 𝐉(𝐱-𝐱_0)].

We may then identify stresslet dipole 𝐒_𝒟=<𝐟^'𝐫𝐫-4μ𝐮^'𝐧𝐫>, rotlet dipole 𝐂_𝒟=<𝐫(𝐫×𝐟^')> and potential dipole 𝐝=<2|𝐫|^2𝐟^'-(𝐫·𝐟^')𝐫 +3μ[4(𝐮^'·𝐧)𝐫 - (𝐮^'·𝐫)𝐧-(𝐫·𝐧)𝐮^']>. Finally, we find the flow field around the particle that decays slower than |𝐱|^-4 as

    𝐮^'(𝐱)=-1/8πμ[𝐅·𝐉(𝐱-𝐱_0)+𝐋·𝐂(𝐱-𝐱_0) -1/2𝐒:𝐊(𝐱-𝐱_0)+1/4𝐒_𝒟⊙∇𝐊(𝐱-𝐱_0)-𝐂_𝒟:Υ(𝐱-𝐱_0) +1/10𝐝·∇^2𝐉(𝐱-𝐱_0) ],

where

    Υ_ijk= ϵ_iksx_sx_j+ϵ_jksx_sx_i/(x_lx_l)^5/2,

is the tensor associated with the rotlet dipole. 

Introducing a more compact notation through 

    𝖲   =[𝐅,𝐋,𝐒,𝐒_𝒟,𝐂_𝒟,𝐝,⋯],
       =[<𝐟^'>,<𝐫×𝐟^'>,<𝐟^'𝐫-2μ𝐮^'𝐧>,<𝐟^'𝐫𝐫-4μ𝐮^'𝐧𝐫>,<𝐫(𝐫×𝐟^')>,<2|𝐫|^2𝐟^'-(𝐫·𝐟^')𝐫 +3μ[4(𝐮^'·𝐧)𝐫 - (𝐮^'·𝐫)𝐧-(𝐫·𝐧)𝐮^']>,⋯],

and

    𝖩   =[𝐉,𝐂,(-1/2)𝐊,(1/4)∇𝐊,-Υ,(1/10) ∇^2𝐉, ⋯],

we obtain

    𝐮^'(𝐱)=-1/8πμ𝖲⊙𝖩.

We note in particular that 

    𝖲 = 𝖥^'+μ<𝖣^'>,

represents the strengths of the multipoles of 𝖩 which contains force moments from the single-layer integral,

    𝖥^'=[<𝐟^'>,<𝐫×𝐟^'>,<𝐫𝐟^'>,<𝐟^'𝐫𝐫>,<𝐫(𝐫×𝐟^')>,<2|𝐫|^2𝐟^'-(𝐫·𝐟^')𝐫>, ⋯],

and terms due to surface disturbance velocity from the double-layer,

    𝖣^'=[0,0,-2𝐧𝐮^',-4𝐮^'𝐧𝐫,0,12(𝐮^'·𝐧)𝐫 - 3(𝐮^'·𝐫)𝐧-3(𝐫·𝐧)𝐮^', ⋯].

Finally, we note 𝖲=𝖲^' because all terms in 𝖲^∞ are zero as the boundary integral equation vanishes identically in that case <cit.>.

To summarize, the second-order force moment is decomposed into a rotlet dipole, a stresslet dipole and a potential dipole, with contributions from the double-layer potentials. To better understand the physical interpretation of these force moments, consider an active particle that propels itself forward using flagella (e.g., E. coli). The force exerted by the flagella is completely balanced by the drag force on the body and since the distribution of these two forces is separated in position (e.g., tail and head), the induced flow field is captured by a stresslet. However, the length scales of the cell body and the flagella differ, often by orders of magnitude, and so this asymmetry gives rise to a stresslet dipole. For a case wherein flagella use rotation to generate thrust, the cell body must counter-rotate to maintain the torque-free motion thereby generating a rotlet dipole. Finally, the finite size of the cell body accounts for the presence of a potential dipole (further discussion may be found in Refs. <cit.>).

We should emphasize that the multipole expansion given in (<ref>) is valid for an active (or passive) particle of any arbitrary shape, in a viscous fluid. Given that 𝖩 is generic for any unbounded single particle, determining the flow field is reduced to finding the strengths 𝖲. Thus, both traction 𝐟^' and disturbance surface velocity 𝐮^'(𝐱∈∂ℬ) are needed. Although the latter may be explicitly prescribed (e.g. through slip velocity), finding the traction generally requires solving the flow field in full. However, one can avoid such calculations by using the Lorentz reciprocal theorem for the Stokes flow <cit.> to calculate moments of 𝐟^'. In the following, we employ the general framework given in <cit.> to find the force moments and hence the multipole strengths 𝖲. Recovering the expressions for force, torque and stresslet, we report explicit formulas for the stresslet dipole, rotlet dipole and potential dipole for an arbitrarily-shaped active particle.





§ EVALUATING THE FORCE MOMENTS OF AN ACTIVE PARTICLE


We are interested in the motion of an active particle with boundary conditions

    𝐮(𝐱∈∂ℬ)=𝐔 + Ω×𝐫+𝐮^s,

where 𝐔 and Ω are the rigid-body translation and rotation of the particle while 𝐮^s is a velocity due to surface activity, such as diffusiophoretic slip or a swimming gait <cit.>. As a dual (or auxiliary) problem, here denoted by a hat, we take the passive motion of a rigid body of the same instantaneous shape,

    𝐮̂(𝐱∈∂ℬ)=𝐔̂ + Ω̂×𝐫.


The reciprocal theorem indicates that the virtual power of the motion of these two bodies is equal <cit.>, namely

    μ̂<𝐧·σ^'·𝐮̂^'>=μ<𝐧·σ̂^'·𝐮^'>,

where 𝐮̂^'=𝐮̂-𝐮̂^∞ and σ̂^'=σ̂-σ̂^∞ are the disturbance flow and stress field for the auxiliary problem, respectively. As we will show below, by using (operators of) the dual problem, the force moments of the active particle may be obtained without resolution of the disturbance field 𝐮^' nor traction 𝐟^'. 

Following <cit.>, we expand the background flow of the auxiliary problem around a point in the body 𝐱_0 as

    𝐮̂^∞(𝐱∈∂ℬ)= 𝐔̂^∞(𝐱_0)+𝐫·∇𝐮̂^∞(𝐱_0)+1/2𝐫𝐫:∇∇𝐮̂^∞(𝐱_0)+⋯,

which by decomposing ∇𝐮̂^∞ and ∇∇𝐮̂^∞ to their irreducible parts, can be rewritten

    𝐮̂^∞(𝐱∈∂ℬ)= 𝐔̂^∞+Ω̂^∞×𝐫+𝐫·𝐄̂^∞+𝐫𝐫:Γ̂^∞+(ϵ·𝐫)𝐫:Λ̂^∞+(2|𝐫|^2𝐈-𝐫𝐫)·𝐞̂^∞+⋯.

Here, 𝐔̂^∞ and Ω̂^∞ indicate the translation and rotation of the background flow of auxiliary problem at 𝐱_0. 𝐄̂^∞=∇𝐮̂^∞ and Γ̂^∞=(1/2)∇∇𝐮̂^∞ are the fully symmetric and deviatoric (Ê_ii^∞=0; Γ̂_iij^∞=Γ̂_iji^∞=Γ̂_jii^∞=0) second and third-order tensors, respectively. Λ̂^∞=∇(∇×𝐮̂^∞) is a second-order symmetric tensor and 𝐞̂^∞=(1/10)∇^2𝐮̂^∞. Using this expansion and relying on the linearity of the Stokes equations, we can write

    p̂^'   =μ̂𝖯̂⊙𝖴̂^',
    𝐮̂^'   =𝖦̂⊙𝖴̂^',
    σ̂^'   =μ̂𝖳̂⊙𝖴̂^',
    𝖥̂^'   =-𝖱̂⊙𝖴̂^',

wherein the velocity gradients

    𝖴̂^'   =[𝐔̂-𝐔̂^∞,Ω̂-Ω̂^∞,-𝐄̂^∞,-Γ̂^∞,-Λ̂^∞,-𝐞̂^∞,⋯],

are linearly mapped to the disturbance pressure, velocity, and stress fields by

    𝖯̂   =[𝐏̂_U,𝐏̂_Ω,𝐏̂_E,𝐏̂_Γ,𝐏̂_Λ,𝐏̂_e,⋯],
    𝖦̂   =[𝐆̂_U,𝐆̂_Ω,𝐆̂_E,𝐆̂_Γ,𝐆̂_Λ,𝐆̂_e,⋯],
    𝖳̂   =[𝐓̂_U,𝐓̂_Ω,𝐓̂_E,𝐓̂_Γ,𝐓̂_Λ,𝐓̂_e,⋯],

which are functions of the position in space and the geometry of the particle and each term maintains the symmetry of the term, against which it operates. 𝖱̂ is the grand resistance tensor that linearly maps velocity moments to force moments in the dual problem. The symmetry of the stress tensor implies that all components of 𝖳̂ are symmetric in their first two (or non-contracted) indices. Under these definitions, the reciprocal theorem can be rewritten as

    <(𝐧·σ^')·𝖦̂>⊙𝖴̂^'=μ<𝐮^'·(𝐧·𝖳̂)>⊙𝖴̂^'.

Importantly, given the boundary condition of the dual problem and the decomposition of the background field, we know the terms in the operator 𝖦̂(𝐱∈∂ℬ) on the boundary of the particle, which may be defined as

    Ĝ_U,ij   =δ_ij,
    Ĝ_Ω,ij   =ϵ_ijsr_s,
    Ĝ_E,ijk   =δ_ijr_k^jk,
    Ĝ_Γ,ijkm   =δ_ijr_kr_m^jkm,
    Ĝ_Λ,ijk   =ϵ_ijsr_sr_k^jk,
    Ĝ_e,ijkm   =2δ_ijr_sr_s-r_ir_j,
         ⋮

where the over-brackets are identical to (<ref>) but only operate over the specified indices. In this way, the operator acts precisely to map the traction 𝐟'=𝐧·σ^' to the force moments

    𝖥^'=<𝐟'·𝖦̂>=[<𝐟^'>,<𝐫×𝐟^'>,<𝐫𝐟^'>,<𝐟^'𝐫𝐫>,<𝐫(𝐫×𝐟^')>,<2|𝐫|^2𝐟^'-(𝐫·𝐟^')𝐫>, ⋯].


Now, given that 𝖴̂^' is arbitrarily chosen, we may discard it from both sides of equation (<ref>). Applying the boundary conditions on 𝐮^', and also expressing the rigid body translation and rotation as 𝖴=[𝐔,Ω,0,⋯], equation (<ref>) can be reduced to <cit.>

    𝖥^'   =-μ/μ̂𝖱̂⊙𝖴+μ<(𝐮^s-𝐮^∞)·(𝐧·𝖳̂)>.



Then by way of (<ref>), one obtains the multipole strengths as

    𝖲=-μ/μ̂𝖱̂⊙𝖴+μ<(𝐮^s-𝐮^∞)·(𝐧·𝖳̂)+𝖣^'>,
 
where the double-layer potential contribution,

    𝖣^'=[0,0,-2𝐧(𝐮^s-𝐮^∞),-4(𝐮^s-𝐮^∞)𝐧𝐫,0,12[(𝐮^s-𝐮^∞)·𝐧]𝐫 - 3[(𝐮^s-𝐮^∞)·𝐫]𝐧-3(𝐫·𝐧)(𝐮^s-𝐮^∞), ⋯],

is simplified as the terms associated with rigid-body motion integrate to zero <cit.>. 

Equation (<ref>) provides the tensorial relationship between the boundary motion and the strength of multipoles for any arbitrarily-shaped active particle in Stokes flow. We note that the multipole strengths are split into terms arising from the rigid-body motion of the particle, 𝖴, and those associated with the (disturbance) surface activity 𝐮^s-𝐮^∞. 
Using this equation, one can derive explicit formulas for 𝖲 provided 𝖳̂ is known, as we illustrate in the following.



We begin with 𝐅 and 𝐋. In self-propulsion, in the absence of any external force and torque, the net force and torque on the particle are strictly zero. However, to find the net translational and rotational velocity (which are unknown at this point), we may use the reciprocal theorem for the force and torque. Upon setting 𝐅=0 and 𝐋=0 in equation (<ref>) we may solve for 𝖴 directly

    [ [ 𝐔; Ω;   ]]
      =μ̂/μ([ 𝐑̂_FU 𝐑̂_LU; 𝐑̂_FΩ 𝐑̂_LΩ;       ])
      ^-1·[ [ 𝐅^s; 𝐋^s;     ]],

where 𝐑̂_FU=-μ̂<𝐧·𝐓̂_U>, 𝐑̂_FΩ=-μ̂<𝐧·𝐓̂_Ω>, 𝐑̂_LU=-μ̂<𝐫×(𝐧·𝐓̂_U)> and 𝐑̂_LΩ=-μ̂<𝐫×(𝐧·𝐓̂_Ω)> are the components of the grand resistance tensor 𝖱̂ associated with rigid-body motion. Here 𝐅^s=μ<(𝐮^s-𝐮^∞)·(𝐧·𝐓̂_U)> is the hydrodynamic force arising solely from the surface activities (often referred to as the thrust) and 𝐋^s=μ<(𝐮^s-𝐮^∞)·(𝐧·𝐓̂_Ω)> is surface activity driven torque. Equation (<ref>)  simply illustrates the balance between the force and torque generated by the surface activities and the hydrodynamic drag.

We may now determine other components of 𝖲 by using (<ref>) at higher orders. We obtain

    𝐒   =-μ/μ̂(𝐑̂_SU·𝐔+𝐑̂_SΩ·Ω)+𝐒^s,
    𝐒_𝒟   =-μ/μ̂(𝐑̂_S_𝒟U·𝐔 +𝐑̂_S_𝒟Ω·Ω)+𝐒_𝒟^s, 
    𝐂_𝒟   =-μ/μ̂(𝐑̂_C_𝒟U·𝐔 +𝐑̂_C_𝒟Ω·Ω)+𝐂_𝒟^s, 
    𝐝   = -μ/μ̂(𝐑̂_dU·𝐔 +𝐑̂_dΩ·Ω)+𝐝^s,

where the resistance tensors may be written in terms of 𝖳̂ as follows

    R̂_SU,ijk   =-μ̂<n_sT̂_U,sikr_j^ij>,
    R̂_SΩ,ijk   =-μ̂<n_sT̂_Ω,sikr_j^ij>,
    R̂_S_𝒟U,ijkm   =-μ̂<n_sT̂_U,simr_j r_k^ijk>,
    R̂_S_𝒟Ω,ijkm   =-μ̂<n_sT̂_Ω,simr_j r_k^ijk>,
    R̂_C_𝒟U,ijk   =-μ̂<n_sT̂_U,snkϵ_jmnr_ir_m^ij>,
    R̂_C_𝒟Ω,ijk   =-μ̂<n_sT̂_Ω,snkϵ_jmnr_ir_m^ij>,
    R̂_dU,ij   =-μ̂<2n_sT̂_U,sijr_l r_l -n_sT̂_U,smjr_mr_i>,
    R̂_dΩ,ij   =-μ̂<2n_sT̂_Ω,sijr_l r_l -n_sT̂_Ω,smjr_mr_i>,

and the contributions of the surface activities are likewise

    𝐒^s   =μ<(𝐮^s-𝐮^∞)𝐧:𝐓̂_E-2(𝐮^s-𝐮^∞)𝐧>,
    𝐒_𝒟^s   =μ<(𝐮^s-𝐮^∞)𝐧:𝐓̂_Γ-4(𝐮^s-𝐮^∞)𝐧𝐫>,
    𝐂_𝒟^s   =μ<(𝐮^s-𝐮^∞)𝐧:𝐓̂_Λ>,
    𝐝^s   =μ<(𝐮^s-𝐮^∞)𝐧:𝐓̂_e+12[𝐧·(𝐮^s-𝐮^∞)]𝐫-3[𝐫·(𝐮^s-𝐮^∞)]𝐧-3(𝐫·𝐧)(𝐮^s-𝐮^∞)>.
 
We should emphasize that all components of 𝖳̂ are unique for a given particle geometry. Therefore by finding them once, we can determine the force moments for any prescribed surface activity provided the shape does not change. 







 §.§ Sphere

We now resolve the force moments of an active spherical particle, using the expressions reported in the previous section. We take 𝐱_0 to be the center of the sphere 𝐫=a𝐧, where a is the radius. Details of the auxiliary flow field and stress field corresponding to each force moment (i.e., 𝖯̂, 𝖦̂ and 𝖳̂) are reported in the appendix.

Having 𝐓̂_U and 𝐓̂_Ω at hand, the rigid-body resistance tensors can be evaluated. We find 𝐑̂_FU=6πμ̂ a𝐈, 𝐑̂_LΩ=8πμ̂ a^3𝐈, 𝐑̂_dU=10πμ̂ a^3𝐈 and [𝐑̂_FΩ,𝐑̂_LU,𝐑̂_SU,𝐑̂_SΩ, 𝐑̂_S_𝒟U,𝐑̂_S_𝒟Ω, 𝐑̂_C_𝒟U,𝐑̂_C_𝒟Ω,𝐑̂_dΩ]=0. 

The force and torque are respectively

    𝐅   = -6π aμ𝐔-3μ/2a<𝐮^s-𝐮^∞> = -6π aμ(𝐔-(1+a^2/6∇^2)𝐔^∞)-3μ/2a<𝐮^s>,
    𝐋   = -8π a^3μΩ+3μ<(𝐮^s-𝐮^∞)×𝐧>=-8π a^3μ(Ω-Ω^∞)+3μ<𝐮^s×𝐧>,

where 𝐔^∞ and Ω^∞ are the velocity and rotation rate of the background flow at the center of the sphere. If the particle is passive, 𝐮^s=0, we recover Faxén's first and second laws as expected. In the absence of an external force and torque, the rigid-body translation and rotation of spherical active particle with surface velocities 𝐮^s are given by

    𝐔   = -1/4π a^2<𝐮^s-𝐮^∞>,
    Ω   = -3/8π a^4<𝐫×(𝐮^s-𝐮^∞)>,

as first shown by <cit.> and later generalized <cit.>.  
Using the expression for stresslet given in equation (<ref>), we find 

    𝐒   =-5μ<(𝐮^s-𝐮^∞)𝐧>-2/3μ<(𝐮^s-𝐮^∞)·𝐧>𝐈,
       =20πμ a^3 /3(1+a^2/10∇^2)𝐄^∞-5μ<𝐮^s𝐧>-2/3μ<𝐮^s·𝐧>𝐈,

where 𝐄^∞=∇𝐮^∞(𝐱_0). We note that this expression for the stresslet amends a typographical error in the results of <cit.> (equation (10) in their reference). When the sphere is passive, equation (<ref>) recovers Faxén's third law.
By using 𝐓̂_Γ, we determine the stresslet dipole 

    𝐒_𝒟   =-35/4μ a<(𝐮^s-𝐮^∞)𝐧𝐧>,
       =14μπ a^5/3(1+a^2/14∇^2)Γ^∞-35/4μ a<𝐮^s𝐧𝐧>,

where Γ^∞=(1/2)∇∇𝐮^∞(𝐱_0). 
The rotlet dipole is then similarly found

    𝐂_𝒟   =4μ a <[(𝐮^s-𝐮^∞)×𝐧]𝐧>,
       =16μπ a^5/15Λ^∞+4μ a <[𝐮^s×𝐧]𝐧>,

where Λ^∞=∇(∇×𝐮^∞)(𝐱_0). 
Finally, for the potential dipole, we arrive at

    𝐝   =-10πμ a^3𝐔+15aμ/2<2[(𝐮^s-𝐮^∞)·𝐧]𝐧-𝐮^s+𝐮^∞>,
       =-10πμ a^3(𝐔-𝐔^∞)+30μπ a^5 𝐝^∞+15aμ/2<2(𝐮^s·𝐧)𝐧-𝐮^s>,

with 𝐝^∞=(1/10)∇^2𝐮^∞(𝐱_0). In total, for a spherical active particle, we have

    𝖲=[   -6π aμ𝐔-3μ/2a<𝐮^s-𝐮^∞>,-8π a^3μΩ+3μ<(𝐮^s-𝐮^∞)×𝐧>,-5μ<(𝐮^s-𝐮^∞)𝐧>-2/3μ<(𝐮^s-𝐮^∞)·𝐧>𝐈,
       -35/4μ a<(𝐮^s-𝐮^∞)𝐧𝐧>,4μ a <[(𝐮^s-𝐮^∞)×𝐧]𝐧>,-10πμ a^3𝐔+15aμ/2<2[(𝐮^s-𝐮^∞)·𝐧]𝐧-𝐮^s+𝐮^∞>,⋯].






 §.§ Generalized squirmer

We now examine the expressions obtained above for the specific case of a sphere with purely tangential surface activity, i.e., a squirmer <cit.>. One may then express 𝐮^s=u^s_θ𝐞_θ+u^s_ϕ𝐞_ϕ in spherical coordinates (r,θ,ϕ) as <cit.>

    u^s_θ   =∑_n=1^∞∑_m=0^n[-2sinθP_n^m^'(ξ)/na^n+2(B_mncos mϕ + B̃_mnsin mϕ)+m P_n^m(ξ)/a^n+1sinθ(C̃_mncos mϕ - C_mnsin mϕ)],
    
    u^s_ϕ   =∑_n=1^∞∑_m=0^n[sinθP_n^m^'(ξ)/a^n+1(C_mncos mϕ + C̃_mnsin mϕ)+2m P_n^m(ξ)/n a^n+2sinθ(B̃_mncos mϕ - B_mnsin mϕ)],

where P_n^m(ξ) is a Legendre function of order m and degree n, and the prime in P_n^m^'(ξ) indicates differentiation with respect to ξ=cosθ. Here, B_mn,B̃_mn,C_mn and C̃_mn are constant coefficients representing different modes of the surface activity. We find the net translational and rotational velocities in terms of these coefficients and Cartesian unit vectors 𝐞_x, 𝐞_y and 𝐞_z as

    𝐔   =-4/3a^3(B_01𝐞_z -B_11𝐞_x-B̃_11𝐞_y),
    Ω   =-1/a^3(C_01𝐞_z -C_11𝐞_x-C̃_11𝐞_y).

Stresslet, stresslet dipole, rotlet dipole and potential dipole can be similarly determined

    𝐒=   -12πμ/a^2[B_02𝐞_z𝐞_z-2B_12𝐞_x𝐞_z-2B̃_12𝐞_y𝐞_z+2B_22(𝐞_x𝐞_x-𝐞_y𝐞_y)+4B̃_02𝐞_x𝐞_y],
    𝐒_𝒟=   -8πμ/a^2[5/3B_03𝐞_z𝐞_z𝐞_z+B_13(𝐞_x𝐞_y𝐞_y-4𝐞_x𝐞_z𝐞_z+𝐞_x𝐞_x𝐞_x)+B̃_13(𝐞_y𝐞_x𝐞_x-4𝐞_y𝐞_z𝐞_z+𝐞_y𝐞_y𝐞_y)
       +10B_23(𝐞_x𝐞_x𝐞_z-𝐞_y𝐞_y𝐞_z)+20B̃_23𝐞_x𝐞_y𝐞_z+10B_33(3𝐞_x𝐞_y𝐞_y-𝐞_x𝐞_x𝐞_x)-10B̃_33(3𝐞_y𝐞_x𝐞_x-𝐞_y𝐞_y𝐞_y)],
    𝐂_𝒟=   -48πμ/5[C_02𝐞_z𝐞_z-2C_12𝐞_x𝐞_z-2C̃_12𝐞_y𝐞_z+2C_22(𝐞_x𝐞_x-𝐞_y𝐞_y)+4C̃_22𝐞_x𝐞_y],
    𝐝=   -80πμ/3(B_01𝐞_z -B_11𝐞_x-B̃_11𝐞_y).

By only keeping B_0n and C_0n terms in equations (<ref>) and (<ref>) and setting the other coefficients to zero, the solution reduces to the axisymmetric motion of a squirmer <cit.>. In this case, by symmetry, all the force moments generated by the surface activity are invariant by rotation with respect to 𝐞_z. Thus, as one can see from (<ref>) to (<ref>), they must be of form 𝐞_z, 𝐞_z𝐞_z, 𝐞_z𝐞_z𝐞_z, ⋯ which are the irreducible traceless rotation-invariant tensors of 𝐞_z. The contribution of non-zero force moments of an axisymmetric squirmer to the flow field is illustrated in Fig. <ref>. Note that, to express the axisymmetric solutions in terms of tangential squirming modes B_n used in <cit.> and <cit.>, one can simply set C_0n=0 and substitute B_0n by -a^n+2/n+1B_n in equations (<ref>) to (<ref>).










 §.§ Axisymmetric slender rod

Let us now consider a slender rod, whose orientation is given by a unit vector 𝐩, with an axisymmetric swimming gait along its length 𝐮^s=α(s)𝐩 in an otherwise quiescent fluid. Here s parameterizes (by arclength) the centerline of the rod, e.g. s∈[-l/2,l/2] where l is the length. 

To find the force moments, we first decompose a surface integral into an integration around the perimeter (denoted by R) in a plane with normal 𝐩 and one over the length of the rod (denoted by s) so that <⋯>=<<⋯>_R>_s. Using the resistive force theory for slender rods <cit.>, we may approximate the force density per unit length

    <𝐧·σ̂^'>_R=-[ζ̂_∥𝐩𝐩+ζ̂_⊥(𝐈-𝐩𝐩)]·𝐮̂^',

where ζ̂_∥ and ζ̂_⊥ are the parallel and perpendicular drag coefficients. Under this approximation, finding 𝖳̂ does not require details of the auxiliary flow field as we illustrate in the following. Recalling that 𝐮̂^'=𝖦̂⊙𝖴̂^' and σ̂^'=μ̂𝖳̂⊙𝖴̂^', one can write

    μ̂<𝐧·𝖳̂>=-<[ζ̂_∥𝐩𝐩+ζ̂_⊥(𝐈-𝐩𝐩)]·𝖦̂>_s.

We note that 𝖦̂(𝐱∈∂ℬ) is known from equations (<ref>) to (<ref>). 

To find the resistance tensors, from equation (<ref>) we have

    μ̂<𝐧·𝐓̂_U>=-<ζ̂_∥𝐩𝐩+ζ̂_⊥(𝐈-𝐩𝐩)>_s,

thus 

    𝐑̂_FU=-μ̂<𝐧·𝐓̂_U>=<ζ̂_∥𝐩𝐩+ζ̂_⊥(𝐈-𝐩𝐩)>_s=[ζ̂_∥𝐩𝐩+ζ̂_⊥(𝐈-𝐩𝐩)]l.

Similarly, we find [𝐑̂_FΩ,𝐑̂_LU,𝐑̂_SU,𝐑̂_S_𝒟Ω,𝐑̂_C_𝒟Ω,𝐑̂_dΩ]=0, 𝐑̂_LΩ=-l^3/12ζ̂_⊥(𝐈-𝐩𝐩), 𝐑̂_dU=l^3/12[ζ̂_∥𝐩𝐩+2ζ̂_⊥(𝐈-𝐩𝐩)] and

    R̂_SΩ,ijk=   l^3/24ζ̂_⊥(ϵ_iksp_j+ϵ_jksp_i)p_s,
    R̂_C_𝒟U,ijk=   l^3/24ζ̂_⊥(p_iϵ_kjm+p_jϵ_kim)p_m,
    R̂_S_𝒟U,ijkm=   l^3/12[ζ̂_∥p_ip_jp_kp_m+ζ̂_⊥(δ_imp_jp_k^ijk-
    p_ip_jp_kp_m)].





Now to find 𝖲, we substitute (<ref>) in (<ref>). Noting that <𝖣>=0 since <𝐧>_R=0, we arrive at

    𝖲   =-μ/μ̂𝖱̂⊙𝖴+μ<𝐮^s·(𝐧·𝖳̂)>,
       =-μ/μ̂𝖱̂⊙𝖴-<α𝐩·[ζ_∥𝐩𝐩+ζ_⊥(𝐈-𝐩𝐩)]·𝖦̂>_s,
       =-μ/μ̂𝖱̂⊙𝖴-ζ_∥𝐩·<α𝖦̂>_s.

Note that ζ_∥=(μ/μ̂)ζ̂_∥ and ζ_⊥=(μ/μ̂)ζ̂_⊥. From this equation, we find 𝐋^s=0, 𝐂^s_𝒟=0 and

    𝐅^s   =-ζ_∥𝐩·<α𝐆̂_U>_s=-ζ_∥<α>_s𝐩,
    𝐒^s   =-ζ_∥𝐩·<α𝐆̂_E>_s=-ζ_∥<sα>_s𝐩𝐩,
    𝐒^s_𝒟   =-ζ_∥𝐩·<α𝐆̂_Γ>_s=-ζ_∥<s^2α>_s𝐩𝐩𝐩,
    𝐝^s   =-ζ_∥𝐩·<α𝐆̂_e>_s=-ζ_∥<s^2α>_s𝐩,

which are in the form of irreducible traceless rotation-invariant tensors with regard to symmetry axis 𝐩, as expected. 


With no external force or torque acting on the rod, we can then determine the translational velocity 

    𝐔=-(1/l)<𝐮^s>_s=-(1/l)<α>_s𝐩,

as also shown by <cit.>. Recalling that 𝐋^s=0, we find Ω=0.
 
 

§ CONCLUSION

In this paper, we investigated the effects of higher-order force moments on the flow field induced by an active particle. Using the boundary integral equations, we expressed the flow as a multipole expansion and decomposed the contribution of second-order force moments into a stresslet dipole, rotlet dipole and a potential dipole. Then, via the reciprocal theorem, we derived explicit formulas for these force moments which are valid for an active particle of arbitrary shape and then evaluated them for a spherical particle, a squirmer and an axisymmetric slender rod. We believe that by providing simple and explicit formulas for more accurate approximations of the flow-fields generated by active particles, we may enhance our understanding of how these particles interact with their surroundings. Given the generality of the employed framework, our results can be extended to capture the effect of third (or higher) order force moments and also can be adapted to study the hydrodynamic interactions between two <cit.> or many <cit.> active particles or active particles near boundaries <cit.>.






§ ACKNOWLEDGEMENT

The authors acknowledge funding from the NSERC Grant No. RGPIN-2014-06577.



*


§ 

Here we present expressions for the auxiliary flow field and stress field considered for an active spherical particle. For simplicity, we define the flow associated with rigid-body translation as 𝐮̂^'_U=𝐆̂_U·(𝐔̂-𝐔̂^∞), which from equation (<ref>) gives 𝐮̂^'_U(𝐱∈∂ℬ)=𝐔̂-𝐔̂^∞. Thus, for a sphere, 𝐏̂_U, 𝐆̂_U and 𝐓̂_U  can be simply found

    P̂_U,i   =3 a/2(x_lx_l)^3/2x_i,
    Ĝ_U,ij   =(3a/4(x_lx_l)^1/2+a^3/4(x_lx_l)^3/2)δ_ij+(3a/4(x_lx_l)^3/2-3a^3/4(x_lx_l)^5/2)x_ix_j,
    T̂_U,ijk   =(-9a/2(x_l x_l)^5/2+15 a^3/2(x_l x_l)^7/2)x_ix_jx_k-3a^3/2(x_l x_l)^5/2(x_i δ_jk+x_j δ_ik+x_k δ_ij).

Similarly, 𝐮̂^'_Ω(𝐱∈∂ℬ)=(Ω̂-Ω̂^∞)×𝐫 leads to

    P̂_Ω,i   = 0,
    Ĝ_Ω,ij   =a^3/(x_lx_l)^3/2ϵ_ijkx_k,
    T̂_Ω,ijk   =3a^3/(x_l x_l)^5/2[ϵ_kisx_j +ϵ_kjs x_i ]x_s.

Relevant to the stresslet calculations, we impose 𝐮̂^'_E(𝐱∈∂ℬ)=𝐫·(-𝐄̂^∞) in which 𝐄̂^∞ is a symmetric and deviatoric second-order tensor and so

    P̂_E,ij   = 5 a^3/(x_lx_l)^5/2x_ix_j,
    Ĝ_E,ijk   =a^5/(x_lx_l)^5/2δ_ijx_k^jk+5/2(a^3/(x_lx_l)^5/2-a^5/(x_lx_l)^7/2) x_ix_jx_k,
    T̂_E,ijkm   =a^5/(x_lx_l)^5/2(δ_kjδ_im + δ_kiδ_jm)+5(7a^5/(x_lx_l)^9/2-5a^3/(x_lx_l)^7/2)x_ix_jx_kx_m
       +5a^3/2(x_lx_l)^5/2(δ_kj x_i x_m+δ_ki x_j x_m+δ_mj x_i x_k+δ_mi x_j x_k)
       -5a^5/(x_lx_l)^7/2(δ_im x_k x_j+δ_jm x_k x_i+δ_ij x_k x_m+δ_jk x_i x_m+δ_ik x_jx_m ).

In determining the stresslet dipole, we have 𝐮̂_Γ^'(𝐱∈∂ℬ)=𝐫𝐫:(-Γ̂^∞). Noting that Γ̂^∞ is a fully symmetric, deviatoric third-order tensor, we find

    P̂_Γ,ijk   =35/4a^5/(x_lx_l)^7/2x_ix_jx_k,
    Ĝ_Γ,ijkm   =1/8(15a^7/(x_lx_l)^7/2-7a^5/(x_lx_l)^5/2) δ_imx_jx_k^jkm+35/8(a^5/(x_lx_l)^7/2-a^7/(x_lx_l)^9/2) x_ix_jx_kx_m,
    T̂_Γ,ijkmn   =15/8a^7{δ_ni[δ_jk x_m +  δ_mjx_k/(x_lx_l)^7/2 - 7x_jx_kx_m/(x_lx_l)^9/2] +δ_nj[δ_ik x_m +  δ_imx_k/(x_lx_l)^7/2 - 7x_ix_k x_m/(x_lx_l)^9/2]  }
       -7/8a^5{δ_ni[δ_jk x_m +  δ_mjx_k/(x_lx_l)^5/2 - 5x_jx_kx_m/(x_lx_l)^7/2] +δ_nj[δ_ik x_m +  δ_imx_k/(x_lx_l)^5/2 - 5x_ix_k x_m/(x_lx_l)^7/2]  }
       +35/4a^5[-7x_i x_j x_k x_m x_n/(x_lx_l)^9/2]-35/4a^7[δ_ij x_m x_n x_k/(x_lx_l)^9/2-9x_i x_j x_m x_n x_k/(x_lx_l)^11/2]
       +35/8(a^5/(x_lx_l)^7/2 - a^7/(x_lx_l)^9/2)[δ_jkx_i x_m x_n+δ_jmx_i x_k x_n+δ_jnx_i x_k x_m+δ_ikx_j x_m x_n + δ_imx_j x_k x_n + δ_inx_j x_k x_m  ].

For a sphere with boundary condition 𝐮̂^'_Λ(𝐱∈∂ℬ)=(ϵ·𝐫)𝐫:(-Λ̂^∞), wherein Λ̂^∞ is a second-order symmetric and deviatoric tensor, we have

    P̂_Λ,ij   =0,
    Ĝ_Λ,ijk   =a^5ϵ_ijmx_mx_k^jk/(x_lx_l)^5/2,
    T̂_Λ,ijkm   =a^5/(x_lx_l)^5/2(ϵ_iksδ_mjx_s+ϵ_jksδ_mix_s+ϵ_kmix_j+ϵ_kmjx_i)
       -5a^5/(x_lx_l)^7/2x_mx_s(  ϵ_iksx_j+   ϵ_jmsx_i).

Finally, 𝐮̂_e^'(𝐱∈∂ℬ)= (2|𝐫|^2𝐈-𝐫𝐫)·(-𝐞̂^∞) yields

    P̂_e,i   =  5 a^3/2(x_lx_l)^3/2x_i,
    Ĝ_e,ij   =(3a^5/4(x_lx_l)^3/2+5a^3/4(x_lx_l)^1/2) δ_ij-(9a^5/4(x_lx_l)^5/2-5a^3/4(x_lx_l)^3/2) x_ix_j,
    T̂_e,ijk   =15/2(3a^5/(x_lx_l)^7/2-a^3/(x_lx_l)^5/2)x_ix_jx_k-9a^5/2(x_lx_l)^5/2(δ_ijx_k+δ_ikx_j+δ_jkx_i).





