Blaise Delmotte
email: blaise.delmotte@imft.fr

Eric E Keaveny
email: e.keaveny@imperial.ac.uk

Franck Plouraboué
email: franck.plouraboue@imft.fr

Eric Climent
email: eric.climent@imft.fr

Large-scale simulation of steady and time-dependent active suspensions with the force-coupling method

Keywords:

We present a new development of the force-coupling method (FCM) to address the accurate simulation of a large number of interacting micro-swimmers. Our approach is based on the squirmer model, which we adapt to the FCM framework, resulting in a method that is suitable for simulating semi-dilute squirmer suspensions. Other effects, such as steric interactions, are considered with our model. We test our method by comparing the velocity field around a single squirmer and the pairwise interactions between two squirmers with exact solutions to the Stokes equations and results given by other numerical methods. We also illustrate our method's ability to describe spheroidal swimmer shapes and biologically-relevant time-dependent swimming gaits. We detail the numerical algorithm used to compute the hydrodynamic coupling between a large collection (10 4 -10 5 ) of micro-swimmers. Using this methodology, we investigate the emergence of polar order in a suspension of squirmers and show that for large domains, both the steady-state polar order parameter and the growth rate of instability are independent of system size. These results demonstrate the effectiveness of our approach to achieve near continuum-level results, allowing for better comparison with experimental measurements while complementing and informing continuum models.

Introduction

Suspensions of active, self-propelled particles arise in both biological systems, such as populations of micro-organisms [START_REF] Dombrowski | Self-concentration and large-scale coherence in bacterial dynamics[END_REF][START_REF] Zhang | Collective motion and density fluctuations in bacterial colonies[END_REF][START_REF] Sokolov | Physical properties of collective motion in suspensions of bacteria[END_REF][START_REF] Ryan | Correlation properties of collective motion in bacterial suspensions[END_REF] and synthetic, colloidal systems [START_REF] Thutupalli | Swarming behavior of simple model squirmers[END_REF]. These suspensions can exhibit the formation of coherent structures and complex flow patterns which may lead to enhanced mixing of chemicals in the surrounding fluid, the alteration of suspension rheology, or, in the biological case, increased nutrient uptake by a population of micro-organisms. In addition to promising applications such as algae biofuels [START_REF] Croze | Dispersion of swimming algae in laminar and turbulent channel flows: consequences for photobioreactors[END_REF][START_REF] Bees | Mathematics for streamlined biofuel production from unicellular algae[END_REF], characterizing the collective dynamics found in these suspensions is of fundamental importance to understanding zooplankton dynamics [START_REF] Guasto | Fluid mechanics of planktonic microorganisms[END_REF][START_REF] Kiorboe | Flow disturbances generated by feeding and swimming zooplankton[END_REF] and mammal fertility [START_REF] Riedel | A self-organized vortex array of hydrodynamically entrained sperm cells[END_REF][START_REF] Creppy | Collective motility of sperm in confined cells[END_REF].

The mathematical modeling of active suspensions entails describing how individual swimmers move and interact in response to the flow fields that they generate [START_REF] Lauga | The hydrodynamics of swimming microorganisms[END_REF][START_REF] Ishikawa | Suspension biomechanics of swimming microbes[END_REF][START_REF] Saintillan | Active suspensions and their nonlinear models[END_REF]. It is particularly important for these models to be able to handle a large collection of swimmers in order to obtain suspension properties at the lab/in situ scale. The modeling of the collec-tive behavior of active matter has been a vibrant area of research during the last decade [START_REF] Ramaswamy | The mechanics and statistics of active matter[END_REF][START_REF] Koch | Collective hydrodynamics of swimming microorganisms: living fluids[END_REF][START_REF] Saintillan | Active suspensions and their nonlinear models[END_REF][START_REF] Marchetti | Hydrodynamics of soft active matter[END_REF], to cite only a few recent reviews. Generally speaking, the modeling approaches can be sorted into two categories: continuum theories and particle-based simulations. Most of the continuum models are generally valid for dilute suspensions where the hydrodynamic disturbances are given by a mean-field description of far-field hydrodynamic interactions [START_REF] Saintillan | Instabilities, pattern formation, and mixing in active suspensions[END_REF][START_REF] Baskaran | Statistical mechanics and hydrodynamics of bacterial suspensions[END_REF][START_REF] Koch | Collective hydrodynamics of swimming microorganisms: living fluids[END_REF]. Recent advances towards more concentrated suspensions include steric interactions [START_REF] Ezhilan | Instabilities and nonlinear dynamics of concentrated active suspensions[END_REF], but the inclusion of high-order singularities due to particle size remains outstanding. Despite this, these models are very attractive as they naturally provide a description of the dynamics at the population level and the resulting equations can be analyzed using a wide range of analytical and numerical techniques.

Particle-based simulations resolve the dynamics of each individual swimmer and from their positions and orientations, construct a picture of the dynamics of the suspension as a whole. As discussed in [START_REF] Koch | Collective hydrodynamics of swimming microorganisms: living fluids[END_REF], particle-based models provide opportunities to (i) test continuum theories, (ii) analyze finite-size effects resulting from a discrete number of swimmers, (iii) explore more complex interactions between swimmers and/or boundaries, and in some cases, (iv) reveal the effects of short-range hydrodynamic interactions and/or steric repulsion. Various models have been proposed in this context, each using different approximations to address the difficult problems of resolving the hydrodynamic interactions and incorporating the geometry of the swimmers. Some of the first such models used point force distributions to create dumbbell-shaped swimmers [START_REF] Hernandez-Ortiz | Transport and collective dynamics in suspensions of confined swimming particles[END_REF][START_REF] Hernandez-Ortiz | Fast computation of many-particle hydrodynamic and electrostatic interactions in a confined geometry[END_REF][START_REF] Underhill | Diffusion and spatial correlations in suspensions of swimming particles[END_REF], slender-body theory to model a slip velocity along the surfaces of rod-like swimmers [START_REF] Saintillan | Orientational order and instabilities in suspensions of self-locomoting rods[END_REF][START_REF] Saintillan | Emergence of coherent structures and large-scale flows in motile suspensions[END_REF], or the squirmer model [START_REF] Lighthill | On the squirming motion of nearly spherical deformable bodies through liquids at very small Reynolds numbers[END_REF][START_REF] Blake | A spherical envelope approach to ciliary propulsion[END_REF] to examine the interactions between spherical swimmers [START_REF] Ishikawa | Hydrodynamic interaction of two swimming model micro-organisms[END_REF]. These initial studies provided important fundamental results connecting the properties of the individual swimmers to the emergence of collective dynamics. Based on their success, these models have been more recently incorporated into a number of numerical approaches for suspension and fluid-structure interaction simulations including Stokesian dynamics [START_REF] Ishikawa | The rheology of a semi-dilute suspension of swimming model micro-organisms[END_REF][START_REF] Ishikawa | Development of coherent structures in concentrated suspensions of swimming model micro-organisms[END_REF][START_REF] Mehandia | The collective dynamics of self-propelled particles[END_REF], the immersed boundary method [START_REF] Lushi | Modeling and simulation of active suspensions containing large numbers of interacting micro-swimmers[END_REF][START_REF] Lambert | Active suspensions in thin films: nutrient uptake and swimmer motion[END_REF], Lattice Boltzmann methods [START_REF] Alarcon | Spontaneous aggregation and global polar ordering in squirmer suspensions[END_REF][START_REF] Pagonabarraga | The structure and rheology of sheared model swimmer suspensions[END_REF], and hybrid finite element/penalization schemes [START_REF] Decoene | Microscopic modelling of active bacterial suspensions[END_REF]. This has allowed for both increased swimmer numbers as well as the incorporation of other effects such as steric interactions, external boundaries, and aligning torques.

In this paper, we introduce an extension of the force coupling method (FCM) [START_REF] Maxey | Localized force representations for particles sedimenting in Stokes flow[END_REF][START_REF] Lomholt | Force-coupling method for particulate two-phase flow: Stokes flow[END_REF], an approach for the large-scale simulation of passive particles, to capture the many-body interactions between active particles. FCM relies on a regularized, rather than a singular, multipole expansion to account for the hydrodynamic interactions between the particles. It includes a higher-order correction due to particle rigidity by enforcing the constraint of zero-averaged strain rate in the vicinity of each particle. Since the force distributions have been regularized, the total particle force, including that associated with the constraints, can be projected onto a grid over which the fluid flow can be found numerically. This allows the hydrodynamic interactions for all particles to be resolved simultaneously. This mesh can be structured and simple such that an efficient parallel Stokes solver can be used to find the hydrodynamic interactions.

We extend FCM to active particles by introducing the regularized singularities in the FCM multipole expansion that have a direct correspondence to the surface velocity modes of the squirmer model [START_REF] Blake | A spherical envelope approach to ciliary propulsion[END_REF]. With these terms included, we then rely on the usual FCM framework to resolve the hydrodynamic interactions in a very efficient manner. We show that by using the full capacity afforded by FCM, we are able to accurately simulate active particle suspensions in the semi-dilute limit with O (10 4 -10 5 ) swimmers. Using this method, we examine the influence of domain size on the steady-state polar order observed for squirmer suspensions. At the same time, we show that our method is quite versatile, being able to handle timedependent swimming gaits, ellipsoidal swimmer shapes, and steric interactions, each at a minimal additional computational cost. We explore in detail how to incorporate biologically-relevant, time-dependent swimming gaits by tuning our model to the recent measurements of the oscillatory flow around Chlamydomonas Rheinardtii [START_REF] Guasto | Oscillatory flows induced by microorganisms swimming in two dimensions[END_REF]. These experiments revealed that considering time-averaged flows for such micro-organisms may oversimplify the hydrodynamic interactions between neighbors. Time-dependency is also closely associated with the way zooplankton feed, mix the surrounding fluid, and interact with each other [START_REF] Croze | Dispersion of swimming algae in laminar and turbulent channel flows: consequences for photobioreactors[END_REF][START_REF] Kiorboe | Flow disturbances generated by feeding and swimming zooplankton[END_REF]. As stated in [START_REF] Saintillan | A quantitative look into microorganism hydrodynamics[END_REF], modeling micro-swimmers with a time-dependent swimming gait might be more realistic and should be included in mathematical models and computer simulations. We show that time-dependence can indeed affect the overall organization of the suspension.

We organize our paper as follows: In Section 2, we review FCM and present the theoretical background for its adaptation to active particles. Section 3 details the numerical method, its algorithmic implementation and how the computational work scales with the particle number. In Section 4, we validate the method and test its accuracy by comparing flow fields, trajectories, and pairwise interactions with previous results available in the literature. The effectiveness of our approach is demonstrated in Section 5 where we present results from large-scale simulations of active particle suspensions. Finally, extensions of FCM to more complex scenarios are introduced in Section 6. We simulate suspensions of spheroidal swimmers and demonstrate the new implementation of time-dependent swimming gaits. Here, we also present preliminary results showing the effect of time-dependence on suspension properties.

Squirmers using FCM

The force-coupling method (FCM) developed by Maxey and collaborators [START_REF] Maxey | Localized force representations for particles sedimenting in Stokes flow[END_REF][START_REF] Lomholt | Force-coupling method for particulate two-phase flow: Stokes flow[END_REF] is an effective approach for the largescale simulation of particulate suspensions, especially for moderately concentrated suspensions at low Reynolds number. In this context, it has been used to address a variety of problems in microfluidics [START_REF] Climent | Dynamics of self-assembled chaining in magnetorheological fluids[END_REF], biofluid dynamics [START_REF] Pivkin | Blood flow velocity effects and role of activation delay time on growth and form of platelet thrombi[END_REF], and micron-scale locomotion [START_REF] Keaveny | Spiral swimming of an artificial micro-swimmer[END_REF][START_REF] Keaveny | Interactions between comoving magnetic microswimmers[END_REF][START_REF] Majmudar | Experiments and theory of undulatory locomotion in a simple structured medium[END_REF]. FCM has also been extended to incorporate finite Reynolds number effects [START_REF] Xu | Numerical simulation of turbulent drag reduction using micro-bubbles[END_REF], thermal fluctuations [START_REF] Keaveny | Fluctuating force-coupling method for simulations of colloidal suspensions[END_REF], near contact lubrication hydrodynamics [START_REF] Dance | Incorporation of lubrication effects into the force-coupling method for particulate two-phase flow[END_REF][START_REF] Yeo | Simulation of concentrated suspensions using the force-coupling method[END_REF], and ellipsoidal particle shapes [START_REF] Liu | Force-coupling method for flows with ellipsoidal particles[END_REF]. With these additional features, FCM has been used to address questions in fundamental fluid dynamics in regimes where inertial effects are important and/or there is a high volume fraction of particles [START_REF] Yeo | Simulation of concentrated suspensions using the force-coupling method[END_REF][START_REF] Yeo | Dynamics of concentrated suspensions of non-colloidal particles in Couette flow[END_REF]. At the same time, FCM has been used to address problems of technological importance, such as micro-bubble drag reduction [START_REF] Xu | Numerical simulation of turbulent drag reduction using micro-bubbles[END_REF] and the dynamics of colloidal particles [START_REF] Keaveny | Fluctuating force-coupling method for simulations of colloidal suspensions[END_REF]. In this section, we expand on [START_REF] Keaveny | Dynamics of structures in active suspensions of paramagnetic particles and applications to artificial micro-swimmers[END_REF] and develop the theoretical underpinnings of FCM's further extension to active particle suspensions using the squirmer model proposed by Lighthill [START_REF] Lighthill | On the squirming motion of nearly spherical deformable bodies through liquids at very small Reynolds numbers[END_REF], advanced by Blake [START_REF] Blake | A spherical envelope approach to ciliary propulsion[END_REF], and employed by Ishikawa et al. [START_REF] Ishikawa | Hydrodynamic interaction of two swimming model micro-organisms[END_REF]. To begin this presentation, we give an overview of FCM, establishing also the notation that will be used throughout the paper.

FCM for passive particles

Consider a suspension of N p rigid spherical particles, each having radius a. Each particle n, (n = 1, . . . , N p ), is centered at Y n and subject to force F n and torque τ n . To determine their motion through the surrounding fluid, we first represent each particle by a low order, finite-force multipole expansion in the Stokes equations

∇p -η∇ 2 u = n F n n (x) + 1 2 τ n × ∇ n (x) + S n • ∇ n (x) ∇ • u = 0. ( 1
)
In Eq. ( 1), S n are the particle stresslets determined through a constraint on the local rate-of-strain as described below. Also in Eq. ( 1) are the two Gaussian envelopes,

n (x) = (2π σ 2 ) -3/2 e -|x-Y n | 2 /2σ 2 , n (x) = (2π σ 2 ) -3/2 e -|x-Y n | 2 /2σ 2 , ( 2
)
used to project the particle forces onto the fluid.

After solving Eq. ( 1), the velocity, V n , angular velocity, n , and local rate-of-strain, E n , of each particle n are found by volume averaging of the resulting fluid flow,

V n = u n (x)d 3 x (3) n = 1 2 [∇ × u] n (x)d 3 x, (4)
E n = 1 2 ∇u + (∇u) T n (x)d 3 x, (5)
where the integration is performed over the entire domain. In order for Eqs. ( 3)-( 5) to recover the correct mobility relations for a single, isolated sphere, namely that V = F/(6πaη) and = τ /(8πa 3 η), the envelope length scales need to be σ = a/ √ π and σ = a/ 6 √ π 1/3 . As the particles are rigid, the stresslets are found by enforcing the constraint that E n = 0 for each particle n [START_REF] Lomholt | Force-coupling method for particulate two-phase flow: Stokes flow[END_REF].

Squirmer model

In addition to undergoing rigid body motion in the absence of applied forces or torques, active and self-propelled particles are also characterized by the flows they generate. To model such particles, we will need to incorporate these flows into FCM. We accomplish this by adapting the axisymmetric squirmer model [START_REF] Lighthill | On the squirming motion of nearly spherical deformable bodies through liquids at very small Reynolds numbers[END_REF][START_REF] Blake | A spherical envelope approach to ciliary propulsion[END_REF][START_REF] Ishikawa | Hydrodynamic interaction of two swimming model micro-organisms[END_REF] to the FCM framework, though we note that a more general squirmer model [START_REF] Pak | Generalized squirming motion of a sphere[END_REF] with non-axysimmetric surface motion could also be considered.

The squirmer model consists of a spherically shaped, self-propelled particle that utilizes axisymmetric surface distortions to move through fluid with speed U in the direction p. If the amplitude of the distortions is small compared to the radius, a, of the squirmer, their effect can be represented by the surface velocity,

v(r = a) = v r r + v θ θ where v r = U cos θ + ∞ n=0 A n (t)P n (cos θ), (6)
v θ = -U sin θ - ∞ n=1 B n (t)V n (cos θ). (7)
Here, P n (x) are the Legendre polynomials,

V n (cos θ) = 2 n(n + 1) sin θ P ′ n (cos θ), (8)
the angle θ is measured with respect to the swimming direction p, and P ′ n (x) = d P n /dx. In order for the squirmer to be force-free, we have Following Ishikawa et al. [START_REF] Ishikawa | Hydrodynamic interaction of two swimming model micro-organisms[END_REF], we consider a reduced squirmer model where A n = 0 for all n and B n = 0 for all n > 2. We therefore only have the first two terms of the series. For this case, the resulting flow field in the frame moving with the swimmer is given by

U = 1 3 (2B 1 -A 1 ). (9)
u r (r, θ) = 2 3 B 1 a 3 r 3 P 1 (cos θ) + a 4 r 4 - a 2 r 2 B 2 P 2 (cos θ) (10)
u θ (r, θ) = 1 3 B 1 a 3 r 3 V 1 (cos θ) + a 4 r 4 B 2 V 2 (cos θ) (11)
where we have used U = 2B 1 /3 from Eq. ( 9). In terms of p, x, and r, this becomes

u(x) = - B 1 3 a 3 r 3 I -3 xx T r 2 p + a 4 r 4 - a 2 r 2 B 2 P 2 p • x r x r -3 a 4 r 4 B 2 p • x r I - xx T r 2 p = u B 1 + u B 2 (12)
Fig. 1 shows an example of a flow field given by Eq. ( 12), as well as the flows u B 1 and u B 2 related to the B 1 and B 2 contributions. While we have already seen that B 1 is related to the swimming speed, it can be shown [START_REF] Ishikawa | Hydrodynamic interaction of two swimming model micro-organisms[END_REF] that B 2 is directly related to the stresslet

G = 4 3 π ηa 2 (3pp -I) B 2 (13)
generated by the surface distortions. This term sets the leading-order flow field that decays like r -2 . We can introduce the parameter β = B 2 /B 1 which describes the relative stresslet strength. In addition, if β > 0, the squirmer behaves like a 'puller', bringing fluid in along p and expelling it laterally, whereas if β < 0, the squirmer is a 'pusher', expelling fluid along p and bringing it in laterally.

To adapt this model to the FCM framework, we first recognize that the flow given by Eq. ( 12) can be represented by the following singularity system in the Stokes equations

∇p -η∇ 2 u = G • ∇ δ(x) + a 2 6 ∇ 2 δ(x) + H∇ 2 δ(x) (14) ∇ • u = 0 (15)
where the degenerate quadrupole is related to B 1 through

H = - 4 3 π ηa 3 B 1 p (16)
and the stresslet G is given by Eq. [START_REF] Ishikawa | Suspension biomechanics of swimming microbes[END_REF]. We can draw a parallel between these singularities and the regularized singularities used with FCM. The stresslet term in Eq. ( 14) is the gradient of the singularity system for a single sphere subject to an applied force. Accordingly, the corresponding regularized singularity in FCM is ∇ (x), where (x) is given by Eq. [START_REF] Zhang | Collective motion and density fluctuations in bacterial colonies[END_REF]. It is important to note that even though we replace two singular force distributions with the one regularized FCM distribution, the particular choice of (x) will yield flows that are asymptotic to both singular flow fields [37]. For the degenerate quadrupole, however, there is not a corresponding natural choice for the regularized distribution. Following [START_REF] Keaveny | Modeling the magnetic interactions between paramagnetic beads in magnetorheological fluids[END_REF], we choose a Gaussian envelope with a length-scale small enough to yield an accurate representation of the singular flow, but not so small as to significantly increase the resolution needed in a numerical simulation (see Section 4.1). We therefore employ the FCM envelope for the force dipole and replace the singular distribution by ∇ 2 (x). Thus, for a single squirmer, the Stokes equations with the FCM squirmer force distribution are

∇p -η∇ 2 u = G • ∇ (x) + H∇ 2 (x) ∇ • u = 0. ( 17
)
As we show Section 4, the flow satisfying Eq. ( 17) closely matches that obtained using the original singularity distribution, Eq. ( 14).

Squirmer interactions and motion

Using FCM, the task of computing the interactions between squirmers is relatively straightforward. We now consider N p independent squirmers where each squirmer has swimming dipole G n and degenerate quadrupole H n . The squirmers may also be subject to external forces F n and torques τ n . Using the linearity of the Stokes equations, we combine Eqs. ( 1) and ( 17) to obtain

∇p -η∇ 2 u = n F n n (x) + 1 2 τ n × ∇ n (x) + S n • ∇ n (x) + G n • ∇ n (x) + H n ∇ 2 n (x) (18)
∇ • u = 0 ( 19
)
for the flow field generated by the suspension.

After finding the flow field, we determine the motion of the squirmers using Eqs. ( 3)-( 5) with two modifications. First, we need to add the swimming velocity, U p n , to Eq. ( 3). Second, we must subtract the artificial, self-induced velocity and the local rate-of-strain due to the squirming modes. The self-induced velocity is given by

W n = A • H n (x)d 3 x (20)
where the tensor A, in index notation, is [START_REF] Maxey | Localized force representations for particles sedimenting in Stokes flow[END_REF] A i j (x) =

1 4π ηr 3 δ i j - 3x i x j r 2 erf r σ √ 2 (21)
-

1 η δ i j - x i x j r 2 + δ i j - 3x i x j r 2 σ r 2 (x). (22)
The self-induced rate-of-strain is given by

K n = 1 2 ∇R • G n + (∇R • G n ) T (x)d 3 x (23)
where the expression for the third rank tensor R can be found in [START_REF] Lomholt | Force-coupling method for particulate two-phase flow: Stokes flow[END_REF]. Taking these self-induced effects into account, the motion of a squirmer n is given by

V n = U p n -W n + u n (x)d 3 x (24) n = 1 2 [∇ × u] n (x)d 3 x (25)
E n = -K n + 1 2 ∇u + (∇u) T n (x)d 3 x. ( 26
)
As before, the stresslets S n due to squirmer rigidity are obtained from the usual constraint on the local rate-of-strain, namely E n = 0 for all n. Squirmer positions Y n and orientations p n are then updated with the Lagrangian equations

dY n dt = V n , ( 27
)
dp n dt = n × p n . (28)
In Section 4, we show through a comparison with the boundary element simulations from [START_REF] Ishikawa | Hydrodynamic interaction of two swimming model micro-organisms[END_REF] that our FCM squirmer model recovers the velocities, angular velocities, stresslets (S n ) and trajectories for two interacting squirmers for a wide range of separations.

Including additional features

With FCM, additional effects can readily be incorporated into the squirmer model and in our subsequent simulations, we consider several of them to demonstrate the versatility of our approach. For example, forces due to steric interactions and external torques experienced by magnetotactic or gyrotactic organisms can be considered by including them in F n and τ n in Eq. ( 18). Additionally, we can extend the FCM squirmer model to ellipsoidal shapes by using the ellipsoidal FCM Gaussian distributions. Such a model can be carefully tuned by comparing with results from [START_REF] Kanevsky | Modeling simple locomotors in Stokes flow[END_REF] and [START_REF] Leshansky | A frictionless microswimmer[END_REF]. The effects of particle aspect ratio on suspension properties can then be explored systematically while still accounting for particle size effects such as Jeffery orbits. Finally, we are not limited to constant values for B 1 and B 2 . By allowing these parameters to be functions of time, the FCM squirmer model can be used to explore how the swimmers' strokes affect overall suspension dynamics. These extensions are addressed in Section 6.

Numerical methods

Fluid solver

The smoothness of the Gaussian force distributions allows FCM to be used with a variety of numerical methods to discretize the Stokes equations. It has been implemented with spectral and spectral element methods [START_REF] Pivkin | Blood flow velocity effects and role of activation delay time on growth and form of platelet thrombi[END_REF][START_REF] Liu | Force-coupling method for flows with ellipsoidal particles[END_REF][START_REF] Yeo | Simulation of concentrated suspensions using the force-coupling method[END_REF] and finite volume methods [START_REF] Loisel | The effect of neutrally buoyant finite-size particles on channel flows in the laminar-turbulent transition regime[END_REF][START_REF] Agbangla | Collective dynamics of flowing colloids during pore clogging[END_REF] in both simple and complex domain geometries.

Here, we use a Fourier spectral method with Fast Fourier Transforms (FFTs) to solve the Stokes equations ( 18)-( 19), for the fluid flow in a three-dimensional periodic domain. We set the zeroth Fourier mode of the velocity to zero, û(k = 0) = 0, to ensure no mean fluid flow through the unit cell. As the fluid encompasses the entire domain, this is equivalent to adding a uniform pressure gradient to balance the net force the particles exert on the fluid in order to maintain the system in equilibrium [START_REF] Hasimoto | On the periodic fundamental solutions of the Stokes equations and their application to viscous flow past a cubic array of spheres[END_REF][START_REF] Maxey | Localized force representations for particles sedimenting in Stokes flow[END_REF]. This also ensures we have a convergent evaluation of the stresslet interactions [START_REF] Brady | Dynamic simulation of hydrodynamically interacting suspensions[END_REF]. As the squirmers can move without exerting a force on the fluid, we do permit a net material flux of squirmers through the unit cell. This corresponds directly to the far-field Stokesian Dynamics computations performed in [START_REF] Ishikawa | Development of coherent structures in concentrated suspensions of swimming model micro-organisms[END_REF] for squirmer suspensions. The FFTs are parallelized with the MPI library P3DFFT. This library uses 2D decomposition of the 3D domain, introducing a better scalability than FFT libraries that implement a 1D decomposition. This decomposition has shown good scalability up to N c = 32,768 cores in Direct Numerical Simulation (DNS) of turbulence [START_REF] Pekurovsky | P3DFFT: a framework for parallel computations of Fourier transforms in three dimensions[END_REF].

Computational work

As explained in [START_REF] Yeo | Simulation of concentrated suspensions using the force-coupling method[END_REF], the number of floating point operations for FCM scales linearly with the number of particles, N p . We find the same scaling for our implementation of FCM. Fig. 2 shows the computational time per times-step for N p up to 80,000 particles with 384 3 ∼ 6 • 10 7 grid points and N c = 256 cores.

Steric interactions

Including steric repulsion is straightforward with FCM. These forces are introduced to both prevent particles from overlapping during the finite time-step and to account for contact forces. For spherical particles, we use the steric barrier described in [START_REF] Dance | Collision barrier effects on the bulk flow in a random suspension[END_REF]. For particles n and m, let r nm = Y m -Y n and r nm = r nm . The repulsive force experienced by n due to steric interactions is

F b n =      - F ref 2a R 2 ref -r 2 nm R 2 ref -4a 2 2γ r nm , for r nm < R ref , 0, otherwise. ( 29
)
where F ref is the magnitude of the force, the cut-off distance R ref sets the distance over which the force acts, and the exponent γ can be adjusted to control the stiffness of the force. Unless specified, all the simulations are run with

F ref /6π ηaU = 4, R ref = 2.
2a and γ = 2. From Newton's third law, we obtain

F b m = -F b n .
It is important to note that this force barrier is generic and it is not intended to model any particular physics. While we could instead use an exponential DLVO-like potential [START_REF] Agbangla | Collective dynamics of flowing colloids during pore clogging[END_REF], the force barrier can be evaluated more rapidly than an exponential potential and by adjusting its parameters we can incorporate repulsion of a similar strength and length scale as DLVO forces. A thorough study on the effect of this force barrier on the dynamics of particulate suspensions is provided in [START_REF] Dance | Collision barrier effects on the bulk flow in a random suspension[END_REF]. Steric repulsion for spheroidal particles requires different modeling. In our study, steric forces and torques are introduced by using a soft repulsive potential with the surface-to-surface distance approximated by the Berne-Pechukas range parameter [START_REF] Allen | Expressions for forces and torques in molecular simulations using rigid bodies[END_REF].

Using a direct pairwise calculation, the evaluation of steric interactions between all particle pairs at each time step would require O (N 2 p /(2N c )) computations per core. This cost is much greater than the O (N p ) cost of the hydrodynamic aspects of FCM. Therefore, instead of a direct calculation, we use the linked-list algorithm described in [START_REF] Allen | Computer Simulation of Liquids[END_REF]. This method divides the computational domain into smaller sub-domains into which the particles are sorted. The edge-length of each sub-domain is slightly larger than R ref .

These sub-domains are distributed over the cores where the steric interactions are evaluated. For a homogeneous suspension, this results in a per core cost for steric interactions that is O 14N 2 p /(2N c N s ) , where N s is the total number of sub-domains. Since for large system sizes, we have N s ≈ L 3 /R 3 ref ≫ 14 and N c ≫ 1, the linked-list algorithm for steric interactions is much more efficient than the direct computation.

Algorithm

We summarize the overall procedure to simulate large populations of micro-swimmers in Stokes flow with the FCM:

• Initialize particle positions Y (0)
n and orientations p

• Start time loop:

for k = 1, . . . , N it 1. Compute Gaussians (k) n (x) and (k) n (x), Eq. (2), 2. Update swimming multipoles G (k)
n , Eq. ( 13), and

H (k)
n , Eq. ( 16), which both depend on p (k) n , 3. Compute steric interactions, Eq. ( 29), with the linked-list algorithm, 4. Add additional forcing if any (gyrotactic torques, magnetic dipoles, . . . ), 5. Project the Gaussian distributions onto the grid (RHS of Eq. ( 19)), 6. Solve Stokes equations, Eqs. ( 18)- (19), to obtain the fluid velocity field u (k) (x), n following [START_REF] Yeo | Simulation of concentrated suspensions using the force-coupling method[END_REF], a) Project all the multipoles onto the grid (RHS of Eq. ( 19)), b) Solve for Stokes equations, Eqs. ( 18)- (19), to obtain the fluid velocity field u (k) (x),

Validation

Comparison with Blake's solution

As explained in Section 2.3, in FCM the velocity field around an isolated squirmer n in an infinite, quiescent fluid is given by

u FCM = A • H n + R : G n , ( 30
)
where A is a second rank tensor given by Eq. ( 21) and the third rank tensor R can be found in [START_REF] Lomholt | Force-coupling method for particulate two-phase flow: Stokes flow[END_REF]. Fig. 3 shows Blake's solution Eq. ( 12) and the velocity field provided by Eq. [START_REF] Ishikawa | Development of coherent structures in concentrated suspensions of swimming model micro-organisms[END_REF] with two different values for the degenerate quadrupole Gaussian envelope size, σ and σ /2, that appears in the tensor A. The streamlines are identical, except for a near-field recirculat- ing region that appears when the width of the degenerate quadrupole envelope is σ (Fig. 3(b)). In this region, however, the magnitude of the velocity is small compared to the swimming speed. This region should not significantly impact the squirmer-squirmer hydrodynamic interactions that we are aiming to resolve.

A quantitative comparison of the velocity field is provided in Fig. 4. The agreement with Blake's solution is very good for r/a > 1.25 when using σ for the width of the degenerate quadrupole envelope. As shown in Fig. 4(b), the smaller envelope size (σ /2) matches Blake's solution more closely for r/a < 1.2, with clear improvement at the front and rear of the squirmer. For this envelope size, the velocity field induced by the degenerate quadrupole u FCM B 1 matches exactly the analytical solution u B 1 in Eq. ( 12) (not shown here). Below this width no quantitative improvement is observed as the remaining error comes from the dipolar contribution, u B 2 .

While this quantitative comparison provides a nice way to choose the degenerate quadrupole envelope size, we must also keep in mind the computational cost associated with decreasing this length scale. Even though it would yield a flow field slightly more in register with Blake's solution, resolving the length scale σ /2 in a 3D simulation would require a grid with 8 times as many points as that needed for σ , the smallest length-scale already in FCM. This would increase computation times by at least an order of magnitude. In addition, the FCM volume averaging used to determine the squirmer translational and angular velocities, Eqs. ( 3)-( 5), will reduce the contribution of the localized velocity field discrepancies to the squirmer-squirmer interactions. In our subsequent simulations, we therefore utilize σ for the degenerate quadrupole envelope size since reducing this length-scale would significantly increase the computational cost, but only provide a minimal improvement.

Pairwise interactions of squirmers

In [START_REF] Ishikawa | Hydrodynamic interaction of two swimming model micro-organisms[END_REF], the authors performed a variety of simulations using the Boundary Element Method (BEM) to compute to high accuracy the pairwise interactions between a squirmer and an inert sphere, and between two squirmers. Here, we consider the same scenarios as those authors and compare results from our FCM simulations with their BEM results.

Interactions between a squirmer and an inert sphere

We first consider the interactions between an inert sphere (labeled "2") located at a point r from the center of a puller squirmer (labeled "1") with β = 5. The direction r/r forms the angle θ with the swimming direction p of the squirmer. The problem setup is depicted in Fig. 5. Fig. 6 compares the velocity of the sphere obtained using FCM simulations with the BEM results and far-field analytical solutions from [START_REF] Ishikawa | Hydrodynamic interaction of two swimming model micro-organisms[END_REF]. We computed the far-field velocity in Fig. 6(a) using the expressions provided in [START_REF] Ishikawa | Hydrodynamic interaction of two swimming model micro-organisms[END_REF] since it was not plotted in Fig. 6(a) of [START_REF] Ishikawa | Hydrodynamic interaction of two swimming model micro-organisms[END_REF] for r/a < 2.5. We suspect that the BEM results were also not plotted for this range. As FCM also resolves the mutually induced particle stresslets, it provides a more accurate estimation than far-field approximation of [START_REF] Ishikawa | Hydrodynamic interaction of two swimming model micro-organisms[END_REF]. As a result, we see that the FCM results very closely match BEM, even in the range where r/a < 3. Similar trends are observed for the angular velocity of the inert sphere (Fig. 6(b)) and the stresslet components (Figs. 6(c),6(d)). These comparisons illustrate the accuracy of the results that can be obtained using FCM, which closely matches BEM but incurs a fraction of the computational cost.

Trajectories of two interacting squirmers

We compute the trajectories of two puller squirmers (β = 5) and compare the results with the BEM simulations of [START_REF] Ishikawa | Hydrodynamic interaction of two swimming model micro-organisms[END_REF]. One squirmer initially swims in the x-direction, p 1 = e x , and the other one in the opposite direction, p 2 = -e x . They are placed with initial separation distance δ y = 1a, . . . , 10a in the transverse direction and δx = 10a in the x-direction. The problem set-up is depicted in Fig. 7. Since the squirmers may collide, we also include steric interactions provided by the force barrier (29).

As shown in Fig. 8, the trajectories match very well with the BEM results for δ y ≥ 2a. When δ y = 1a, the collision barrier and near-field hydrodynamic interactions play an important role in determining the overall squirmer trajectories. Fig. 9 shows the effect of the steric repulsion parameters F ref and R ref on the squirmer trajectories. The specific values of these parameters are provided in Appendix B. We see that by varying the barrier parameters one can obtain trajectories that closely match the results of [START_REF] Ishikawa | Hydrodynamic interaction of two swimming model micro-organisms[END_REF].

Simulation results

Large suspensions of swimming micro-organisms

The squirmer model has been widely used to investigate both the behavior of single swimmers [START_REF] Doostmohammadi | Low-Reynolds-number swimming at pycnoclines[END_REF][START_REF] Wang | Inertial squirmer[END_REF][START_REF] Michelin | Unsteady feeding and optimal strokes of model ciliates[END_REF][START_REF] Zhu | Low-Reynolds-number swimming in a capillary tube[END_REF], as well as their collective dynamics and interactions [START_REF] Ishikawa | The rheology of a semi-dilute suspension of swimming model micro-organisms[END_REF][START_REF] Ishikawa | Development of coherent structures in concentrated suspensions of swimming model micro-organisms[END_REF][START_REF] Evans | Orientational order in concentrated suspensions of spherical microswimmers[END_REF][START_REF] Alarcon | Spontaneous aggregation and global polar ordering in squirmer suspensions[END_REF][START_REF] Lambert | Active suspensions in thin films: nutrient uptake and swimmer motion[END_REF]. Simulations of suspensions revealed that the overall population dynamics depend strongly on the squirming parameter β. In particular, when |β| is small, the isotropic state for a periodic suspension has been shown to be unstable, and the suspension evolves to a polar steady-state with a non-zero value of the polar order parameter P (t) = : barrier used in Fig. 8. Fig. 10. Snapshots of the orientational state in a semi-dilute suspension (φ v = 0.1) containing N p = 37,659 swimmers with β = 1. p is the mean steadystate orientation vector on the unit sphere defined in Eq. [START_REF] Lushi | Modeling and simulation of active suspensions containing large numbers of interacting micro-swimmers[END_REF]. p n • p quantifies the degree of alignment with the mean direction p of swimmer n.

This instability has been studied numerically by [START_REF] Ishikawa | Development of coherent structures in concentrated suspensions of swimming model micro-organisms[END_REF] and [START_REF] Evans | Orientational order in concentrated suspensions of spherical microswimmers[END_REF] using Stokesian dynamics with N p = 64 swimmers for volume fractions φ v = 0.01 → 0.5. Using the Lattice-Boltzmann method, [START_REF] Alarcon | Spontaneous aggregation and global polar ordering in squirmer suspensions[END_REF] observed the same behavior for N p = 2000 and φ v = 0.1.

Polar order parameter

Using our FCM model, we study this instability and the resulting polar order of a squirmer suspension. In particular, we examine how the domain size affects both the growth rate of the instability and the final steady-state. The influence of domain size has not been addressed previously for squirmer suspensions though it has been observed in simulations of rod-like swimmers [START_REF] Saintillan | Emergence of coherent structures and large-scale flows in motile suspensions[END_REF]. We perform simulations of semi-dilute suspensions (φ v = 0.1) of puller squirmers (β = 1) in triply-periodic square domains with edge lengths ranging from L/a = 14 to L/a = 116. As the volume fraction φ v is fixed, varying L/a increases of the number of swimmers from N p = 64 (as in [START_REF] Ishikawa | Development of coherent structures in concentrated suspensions of swimming model micro-organisms[END_REF][START_REF] Evans | Orientational order in concentrated suspensions of spherical microswimmers[END_REF]) to N p = 37,659. We initialize a homogeneous, isotropic suspension by distributing the swimmer positions uniformly in the domain and the swimming directions uniformly over the unit sphere. Depending on domain size, the simulations are run to final time t f = 1000 -1500a/U with a time step t = 0.005a/U . Thus, each simulation requires between 2 × 10 5 -3 × 10 5 time-steps. Fig. 10 illustrates the polar ordered state for a simulations with L/a = 116, N p = 37,659. We quantify the degree of alignment of each swimmer n, p n • p , with the mean steady-state orientation p given by

p = lim t→∞ 1 N p N p n=1 p n (t) P (t) , (32)
where p lies on the unit sphere S 2 . At t = 0, there is no clear mean orientation, while at t = 1000a/U , a significant proportion of particles are aligned with the mean direction p . Fig. 11a shows P (t), Eq. ( 31), for the simulations with different domain sizes. For each domain size, we see that the suspension evolves from the initial isotropic state to one that has polar order. We observe, however, that the final value, P ∞ , of the polar order parameter depends on the domain size. We find that it decreases as L/a (and N p ) increases. The data also shows that as L/a → ∞, P ∞ decays like (L/a) -3/2 (or

N -1/2 p
) and reaches an asymptotic value of P ∞ → 0.452.

These results would indicate that polar order should also arise in an unbounded suspension. We also observe that in the polar ordered state, the average swimming speed is 4% less than its value in isotropic state which itself is nearly the free swimming value. Indeed, we remind that in our simulations, a zero net flux of the suspension is imposed, therefore, when the swimmers are aligned, they are slowed down by the backflow.

From our simulations, we may also analyze how L/a affects the time evolution of the instability. In Fig. 11(a), we clearly see that the time it takes to reach the final polar state increases with domain size. We analyze this data in more detail in Fig. 12, now plotting it in semilogarithmic scale. For each case, we find that after an initial transient state, the instability grows exponentially and with a growth rate that is independent of the system size (see inset figure). It would certainly be interesting to investigate if this result could be reproduced by a linear stability analysis of continuum model for a squirmer suspension.

Orientational distribution

We can examine the polar order in more detail by computing the orientational distribution (θ, φ, t) defined over the unit sphere [START_REF] Saintillan | Extensional rheology of active suspensions[END_REF]. Here, θ = cos -1 (p z ) corresponds to the elevation angle while φ = tan -1 (p y /p x ) gives the azimuthal angle. Fig. 13 shows (θ, φ, t) normalized by the isotropic distribution 0 = 1/(4π ) at times before and after the transition to polar order for the case where N p = 11, 158 and L = 78a. As expected, before the transition, the distribution is nearly uniform over the surface of the sphere. After the polar state is reached, we see that (θ, φ) is narrowly distributed around the mean direction p . We note that the steady state mean direction depends both on the random initial seeding of swimmers and the domain geometry. We find that system first aligns along an arbitrary direction due to the random initial seeding of swimmers, but as time goes on, the mean director tends to align with the normal to one of the periodic boundaries. As a consequence, the steady state mean direction is biased by the domain shape that breaks radial symmetry. We would like to stress that, however, as demonstrated here and in [START_REF] Ishikawa | Development of coherent structures in concentrated suspensions of swimming model micro-organisms[END_REF], the polar ordering itself is a consequence of squirmer hydrodynamic and steric interactions rather than the domain shape. Fig. 14 shows the steady-state averaged distribution (θ, φ)| p in the frame where the mean direction is given by θ = 0 and φ = 0. We see that the resulting distribution is axisymmetric in that it does not depend on φ. This is not surprising as the flow field induced by a squirmer is axisymmetric. Orientational distributions could also be obtained from continuum models, as was done for swimmers in external flow fields [START_REF] Saintillan | Extensional rheology of active suspensions[END_REF]. Our results could be compared with continuum models of squirmer suspensions, though, to the best of our knowledge, the continuum models that are currently available in the literature predict a stable isotropic state for spherical swimmer suspensions.

Spatial distribution

Along with the orientational distribution, we also examine the spatial distribution of squirmers by computing the Voronoi tessellation for the set of points corresponding to the squirmers' centers. We have used the C++ library Voro++ [START_REF] Rycroft | Voro++: a three-dimensional Voronoi cell library in c++[END_REF] to perform the computation and determine the volume, V V , of each Voronoi cell. Fig. 15 shows the evolution of the standard deviation, σ V , of the V V distribution for the case where N p = 11,158 and L = 78a. We see sudden increase in σ V as the suspension transitions to polar order, corresponding to a widening of the distribution. The inset in Fig. 15 shows the time-averaged distribution before and after the transition. Before polar ordering occurs, the mode of the distribution is V V ≈ 34a 3 , which is slightly less than the average value V V = L 3 /N p = 42a 3 that one might expect. During the polar steady-state, we find that the mode decreases to V V ≈ 30a 3 , which is indicative of a slight clustering of the particles.

To further examine the spatial distribution, we compute the steady-state pair distribution function g(r, θ) which gives the probability of finding a squirmer at a distance r = |r| and with elevation angle θ = cos -1 (r • p/r) from another squirmer that has swimming direction p. Fig. 16(a) shows g(r, θ) for the case where N p = 11,158 and L = 78a. We see clearly that g(r, θ) depends not only on r, but on θ as well. At steady-state, for a given squirmer there is a significantly higher probability of finding another squirmer in front of it (θ = 0) rather than behind it. If we integrate g(r, θ) over θ , we obtain the radial distribution function shown in Fig. 16(b). We find that g(r) exhibits a peak at two radii from the swimmer surface. This peak suggests the existence of particle clusters whose sizes are greater than two individuals.

Orientational correlations

The hydrodynamic and steric interactions between the squirmers also lead to correlations in their orientations. We compute the steady-state orientational correlation function [START_REF] Ishikawa | Development of coherent structures in concentrated suspensions of swimming model micro-organisms[END_REF] where the brackets • denote the ensemble average that we compute by averaging over both time and squirmer pairs. Fig. 17(a) shows the correlations in the frame of a squirmer located at the origin and with p = ẑ and the definitions of r and θ are identical to those in Section 5.1.3. We find that the highest correlations occur near contact (r ≈ 2) at the angles θ ≈ π/3 and θ ≈ 3π /4. Despite the overall polar order, we can identify two distinct regions around θ = 0 and θ = π where the orientations are uncorrelated. We do, however, find a positive value for the correlations, even far away from the origin. This is most evident in the θ -averaged correlation, I p (r), shown in Fig. 17(b) which approaches a finite value I p = 0.22 as r increases. This is consistent with the observed long-range polar order of the suspension. We also find that I p (r) > 0 for all r. This again may be a result of the strong polar ordering of the suspension.

I p (r) = p(x) • p(x + r) , (33)

Extensions to ellipsoidal swimmers and time-dependent swimming gaits

In this section, we demonstrate the extension of our approach to both spheroidal swimmer shapes and time-dependent swimming gaits. This illustrates the versatility of FCM while preserving good computational scalability and an accurate treatment of hydrodynamic interactions.

FCM for ellipsoidal particles

To extend FCM to passive ellipsoidal particles, one simply modifies the Gaussian envelopes in Eq. ( 2) [START_REF] Liu | Force-coupling method for flows with ellipsoidal particles[END_REF]. For example, taking the orthonormal vectors ê1 , ê2 , and ê3 to be aligned with the ellipsoid semi-axes having lengths a 1 , a 2 , and a 3 respectively, the Gaussian envelope corresponding to n (x)

is ell n (x) = (2π ) -3/2 (σ ;1 σ ;2 σ ;3 ) -1 exp - 1 2 (x -Y n ) T Q T Q(x -Y n ) (34)
where

σ ;i = a i / √ π for i = 1, 2, 3, Q = (ê 1 ê2 ê3 ), and =    σ -2 ;1 0 0 0 σ -2 ;2 0 0 0 σ -2 ;3    . (35)
A similar expression is used for ell n (x) with σ ;i = a i / 6 √ π 1/3 . Beyond this, the underlying algorithm of projecting the particle forces onto the fluid and volume averaging the resulting fluid flow remains unchanged. The constraint that E n = 0 for each n is still used to find the stresslets. Thus, using FCM to compute the motion and hydrodynamic interactions of ellipsoids does not require any additional steps in the algorithm described in Section 3.4. As explained in Section 3.3, the modeling of steric interactions differs from that for spherical particles and our simple force barrier could not be used. An extensive validation of FCM for ellipsoidal particles is presented in [START_REF] Liu | Force-coupling method for flows with ellipsoidal particles[END_REF] where many of the classical results for ellipsoidal particles in Stokes flow, e.g. Jeffery's orbits, are shown to be recovered exactly with FCM. We note that for very slender particles, FCM would not be appropriate and computations based on slender-body theories [START_REF] Tornberg | A numerical method for simulations of rigid fiber suspensions[END_REF] should be performed. To extend FCM to active ellipsoidal particles requires adding the stresslet and possible potential dipole terms to the multipole expansion. As for spherical particles (see Section 2.3, Eqs. ( 20) and ( 23)), these additional multipoles will lead to artificial, self-induced velocities and local rates-of-strain. These effects must be subtracted away using the formula derived in Appendix A. It is worth reiterating that these are the only steps that need to be added to the FCM algorithm for passive particles to simulate active ones.

Spheroidal swimmer simulations

Previous particle-based simulations [START_REF] Saintillan | Emergence of coherent structures and large-scale flows in motile suspensions[END_REF][START_REF] Lushi | Modeling and simulation of active suspensions containing large numbers of interacting micro-swimmers[END_REF] and continuum theory results [START_REF] Saintillan | Instabilities, pattern formation, and mixing in active suspensions[END_REF][START_REF] Baskaran | Statistical mechanics and hydrodynamics of bacterial suspensions[END_REF][START_REF] Ezhilan | Instabilities and nonlinear dynamics of concentrated active suspensions[END_REF] predict an unstable isotropic state for suspensions of prolate spheroidal pushers. Using FCM, we simulate a dilute suspension, φ v = 0.05, of N p = 1500 spheroidal pushers with aspect ratio a 1 a 2 = a 1 a 3 = 3. For these simulations, the stresslet parameter is B 2 = -1.5 and the degenerate quadrupole is set to zero, B 1 = 0. Steric interactions are included by using a soft repulsive potential with a ∼ r -5 decay and surface-to-surface distance approximated by the Berne-Pechukas range parameter [START_REF] Allen | Expressions for forces and torques in molecular simulations using rigid bodies[END_REF]. Fig. 18(a) shows a snapshot of the suspension where clusters of swimmers have velocities nearly 1.6 times larger than the isolated swimming speed. As in [START_REF] Saintillan | Emergence of coherent structures and large-scale flows in motile suspensions[END_REF][START_REF] Lushi | Modeling and simulation of active suspensions containing large numbers of interacting micro-swimmers[END_REF], we see that the isotropic state is not stable and the polar order parameter increases with time (Fig. 18(b)). We do not see as large an increase in P (t) as in [START_REF] Saintillan | Emergence of coherent structures and large-scale flows in motile suspensions[END_REF][START_REF] Lushi | Modeling and simulation of active suspensions containing large numbers of interacting micro-swimmers[END_REF] due to the relatively low aspect ratio of our swimmers. The tumbling effect due to Jeffery orbits and the steric torques that tend to align adjacent swimmers are much lower than they are for rod-like swimmers.

Suspension dynamics of time-dependent swimmers

Here, we show how to incorporate time-dependence into our model. Using the procedure outlined in [START_REF] Ghose | Irreducible representations of oscillatory and swirling flows in active soft matter[END_REF], we determine the time-dependent multipole coefficients B 1 (t) and B 2 (t) from the experimental data provided in [START_REF] Guasto | Oscillatory flows induced by microorganisms swimming in two dimensions[END_REF]. We then perform suspension simulations that show time-dependence at the level of individual swimmers can affect the overall suspension properties.

A single time-dependent swimmer

Recent experiments [START_REF] Guasto | Oscillatory flows induced by microorganisms swimming in two dimensions[END_REF] quantified the periodic swimming gait and resulting flow field of the algae cell Chlamydomonas Rheinardtii. They extracted the swimming speed, the induced velocity field, and the power dissipation, showing also that all can be represented as periodic functions of time. A recent theoretical investigation [START_REF] Ghose | Irreducible representations of oscillatory and swirling flows in active soft matter[END_REF] showed that these quantities could be reproduced using a multipole-based model. In their study, they considered three time-dependent multipoles: a stresslet, a degenerate quadrupole (or potential dipole) and a "septlet". The stresslet decays as r -2 whereas the degenerate quadrupole and the "septlet" decay like r -3 . Here, we utilize only the stresslet and degenerate quadrupole terms and find that they are sufficient to reproduce the features measured by [START_REF] Guasto | Oscillatory flows induced by microorganisms swimming in two dimensions[END_REF].

The measured swimming speed from [START_REF] Guasto | Oscillatory flows induced by microorganisms swimming in two dimensions[END_REF], U (t), can be represented using the truncated Fourier series in time (see [START_REF] Ghose | Irreducible representations of oscillatory and swirling flows in active soft matter[END_REF]): U (t) = a 0 + a 1 cos (ωt) + a 2 cos (2ωt) + b 1 sin (ωt) + b 2 sin (2ωt) [START_REF] Decoene | Microscopic modelling of active bacterial suspensions[END_REF] where ω is the frequency of the swimming gait. The mean swimming speed over one beat period, T = 2π /ω, is given by a 0 = 49.54a s -1 , where a = 2.5 µm is the radius of the microorganism. The values for the remaining coefficients are : FCM, : results from [START_REF] Guasto | Oscillatory flows induced by microorganisms swimming in two dimensions[END_REF].

provided in the supplementary information of [START_REF] Ghose | Irreducible representations of oscillatory and swirling flows in active soft matter[END_REF]. From Eq. ( 9), we can immediately determine the time-dependent degenerate quadrupole strength

B 1 (t) = 3 2 U (t) (37)
in order to preserve the instantaneous force-free condition. Unlike this term, there is more than one way to calibrate the stresslet strength, B 2 (t). For example, one could determine B 2 (t) using the power dissipation measurements from [START_REF] Guasto | Oscillatory flows induced by microorganisms swimming in two dimensions[END_REF] and Eq. ( 9) from [START_REF] Blake | A model for the micro-structure in ciliated organisms[END_REF]

d (t) = 2 3 πηa 8B 1 (t) 2 + 4B 2 (t) 2 (38)
for the power dissipated by a squirmer. Using this approach, we found that our resulting flow field did not match the experimental results from [START_REF] Guasto | Oscillatory flows induced by microorganisms swimming in two dimensions[END_REF]. We instead determine B 2 (t) directly from the experimental flow field by fitting to the location of the moving stagnation point. This is the same approach adopted by [START_REF] Ghose | Irreducible representations of oscillatory and swirling flows in active soft matter[END_REF]. We utilize three Fourier modes to describe the time evolution of B 2

B 2 (t) = c 0 + c 1 cos (ωt + ϕ c1 ) + c 2 cos (2ωt + ϕ c2 ) + c 3 cos (3ωt + ϕ c3 ) + s 1 sin (ωt + ϕ s1 ) + s 2 sin (2ωt + ϕ s2 ) + s 3 sin (3ωt + ϕ s3 ) , (39)
where the amplitudes c 0 , c 1 , s 1 , . . . and phases ϕ c1 , ϕ s1 , . . . are fitted manually. Appendix C gives the resulting values of these parameters. The phase diagram in Fig. 19(a) shows the value of B 1 versus B 2 and it is similar to that found by [START_REF] Ghose | Irreducible representations of oscillatory and swirling flows in active soft matter[END_REF]. We extract the average value of

β(t) = B 2 (t)/B 1 (t) over one beat cycle β = 1 T T 0 B 2 (t)/B 1 (t)dt = 0.1 (40)
which corresponds to a puller squirmer with a relatively small stresslet magnitude. Fig. 19(b) shows the resulting power dissipation as determined from Eq. [START_REF] Lomholt | Force-coupling method for particulate two-phase flow: Stokes flow[END_REF]. It reaches a peak value at t/T ≈ 0.3, which coincides with the time at which the swimming speed reaches its maximum value. We note that unlike [START_REF] Guasto | Oscillatory flows induced by microorganisms swimming in two dimensions[END_REF] our swimmers generate axisymmetric flow fields and we are considering a 3D periodic domain. Despite this, we achieve a qualitatively similar power dissipation profile with slightly greater values during the first half of the beat cycle. Fig. 20 shows the flow field around our model of C. Rheinardtii at six different times during its beat cycle. These time points are chosen to correspond to those in Fig. 3 of [START_REF] Guasto | Oscillatory flows induced by microorganisms swimming in two dimensions[END_REF]. We achieve very similar streamlines and, by construction, the position of the stagnation point matches the experimental data very well. We also note that our flow field is similar to that given by the multipole model found in [START_REF] Ghose | Irreducible representations of oscillatory and swirling flows in active soft matter[END_REF] even though we do not include the rapidly decaying "septlet" term in our model. These results illustrate that our properly tuned, time-dependent squirmer model can yield flow fields very similar to those of real organisms.

Interactions between two model C. Rheinardtii

We first consider pairwise interactions between two model C. Rheinardtii for two initial configurations, δ y = 1a and 2a (cf. Section 4.2.2). We introduce a phase shift, ϕ, between their swimming cycles. When ϕ = 0, the swimmers are synchronized, whereas for ϕ = π they are completely out of phase. Fig. 21 shows the effect of this phase shift on the trajectories. For each case, we show only the trajectory of the swimmer labeled "2" in Fig. 7. We also provide the trajectories of steady swimmers with β = 0 and β = 0.1 for comparison. Recall that β = 0.1 corresponds to the average value of B 2 (t)/B 1 (t) in our model. We see that the trajectories, including the final positions and orientations, do depend on ϕ. We find, however, that these variations are small compared to the swimmer separation distance (see the inset of Fig. 21 for the relative trajectories in the frame of one swimmer). The trajectories are also close to those for steady squirmers with β = 0, 0.1.

Collective dynamics of time-dependent swimmers

In this section, we present results from suspension simulations using our unsteady model for C. Rheinardtii. To the best of our knowledge, similar results have not appeared in the literature. We aim to show in this initial study that timedependence of individual swimmers can affect their overall distribution. Specifically, we show that the distribution of beat phases affects the steady-state polar order of the suspension. Fig. 22 shows a sketch of the simulations with time-dependent swimmers.

We consider a suspension of N p = 1395 time-dependent swimmers distributed uniformly in a triply periodic domain. The volume fraction is φ v = 0.1. We take the initial swimming directions to be distributed uniformly over the unit-sphere.

For swimmer n, its gait is characterized by B 1 (t + ϕ n ) and B 2 (t + ϕ n ), where ϕ n is the value of its beat phase. We consider two different distributions of the beat phases. The phases are either uniformly distributed with ϕ n ∈ [0; 2π ], or ϕ n = 0 for all n such that all swimmers are synchronized. We run the simulations for 3700 dimensionless time units, corresponding to 4000 beat cycles. As mentioned in Section 6.3.1, the mean velocity of C. Rheinardtii is U = 49.54a s -1 . Since we take the time-step t = 0.0025a/U , the total time for our simulations correspond to 1.5 × 10 6 t.

Fig. 23(a) shows the evolution of the polar order parameter for these two distributions of beat phase. Also shown is the polar order parameter for suspensions of steady swimmers with swimming speeds equal to the mean swimming speed for the time-dependent case and either β = 0 or β = 0.1. We find that when there is a uniform distribution of the beat phase, we achieve results that match the steady case with β = 0.1. On the other hand, if the swimmers are synchronized, we find that the polar order parameter matches that for the steady case with β = 0. Fig. 23(b) shows the steady-state orientational distributions (θ)| p about the mean director. Again, we see that the distribution for case of a uniform distribution of beat phases matches that for the steady case with β = 0.1, while the synchronized suspension has the same distribution as the β = 0 case. These new results show that the distribution of beat phases can strongly affect the orientational order of the suspension. : uniform distribution 0 (θ) = 1/π .

Conclusions

In this study, we presented an extension of FCM to compute the hydrodynamic interactions between a large number (10 4 -10 5 ) of active self-propelled particles in semi-dilute suspensions. Our approach builds from the spherical squirmer model developed by Blake [START_REF] Blake | A spherical envelope approach to ciliary propulsion[END_REF] by including in the usual FCM force distribution the regularized singularities that correspond to the surface squirming modes. We have shown that our model readily allows for different swimming gaits, as well as spheroidal swimmer shapes, and can also account for effects such as steric repulsion. We demonstrated the accuracy of our model by comparing velocity fields, pairwise interactions, and trajectories with the analytical or numerical results given by the full squirmer model. In our squirmer suspension simulations, we recovered the polar order found in other studies, but also quantified the effects of domain size on the final steady-state. Specifically, we showed that for puller suspensions with β = 1, the final value of the polar order parameter decreases and the time needed to reach steady state increases as the system gets larger. We did, however, find that these quantities appeared to converge to an asymptotic value and that polar order should be present even for an infinitely large domain. The scale of our simulations allowed us to compile robust statistics and examine in detail the orientational distribution and local micro-structure of the suspension. By extending the model to time-dependent swimming gaits, we also illustrated for the first time, to best of our knowledge, that time-dependence at the level of individual swimmers can affect the final steady-state distribution of the suspension.

We see that the squirmer model within the FCM framework provides an effective computational approach to simulate active suspensions, allowing for particle numbers that begin to approach continuum levels. As a result, we feel that our computational model can both be used to complement experimental research, as well as inform the development of continuum models. For example, we have shown that the time-dependency of swimming gaits can be readily included in our model by tuning it to available experimental data [START_REF] Guasto | Oscillatory flows induced by microorganisms swimming in two dimensions[END_REF]. By tuning our model to similar experimental data, but for a wider zoology of microorganisms [START_REF] Kiorboe | Flow disturbances generated by feeding and swimming zooplankton[END_REF], we could assess the possible differences in collective dynamics exhibited by different species, or even look into how one species might interact with another. In addition, our model could be an effective approach to address questions regarding the mixing of background scalar fields [START_REF] Lambert | Active suspensions in thin films: nutrient uptake and swimmer motion[END_REF] or passive tracers dynamics [START_REF] Leptos | Dynamics of enhanced tracer diffusion in suspensions of swimming eukaryotic microorganisms[END_REF][START_REF] Lin | Stirring by squirmers[END_REF] in active, time-dependent suspensions. The effects of tracer Brownian motion [START_REF] Keaveny | Fluctuating force-coupling method for simulations of colloidal suspensions[END_REF] can be included with squirmer model and the competition between tracer diffusion and advection due to the flows produced by the swimmers can be assessed.

Fig. 1 .

 1 Fig. 1. Decomposition of squirmer velocity field for β = 1.

Fig. 2 .

 2 Fig. 2. Scaling of the FCM with the number of squirmers N p . Computational time per time-step versus N p , or equivalently, versus volumetric fraction φ v . N c = 256 cores work in parallel for a cubic domain with 384 3 grid points.

7 .

 7 Compute particle rate of strains E (k) n , Eq. (26), 8. If E (k) n > ε, compute stresslets S(k)

9 .

 9 Compute particle velocities V (k) n , Eq. (24), and rotations (k) n , Eq. (25), 10. Integrate Eqs. (27) and (28) using the fourth-order Adams-Bashforth scheme to obtain Y

Fig. 3 .

 3 Fig. 3. Velocity field u/U around a puller squirmer (β = 1) swimming to the right. (a) Blake's solution; (b) FCM solution with σ for the degenerate quadrupole envelope. (c) FCM solution with σ /2 for the degenerate quadrupole envelope.

Fig. 4 .

 4 Fig. 4. Comparison with Blake's solution, Eq. (12). (a) Normalized difference, u Blakeu FCM /U , between FCM and Blake's solution for a puller squirmer (β = 1). The half-width for the degenerate quadrupole envelope is σ . : 10% iso-value. (b) Velocity profile along the swimmer axis for Blake's solution and the FCM approximation for the two different degenerate quadrupole envelope sizes.

Fig. 5 .

 5 Fig. 5. Sketch showing the set-up of our computations of the interactions between squirmer "1" and inert sphere "2".

Fig. 6 .

 6 Fig. 6. The (a) radial velocity |U r,2 |, (b) angular velocity | z,2 |, (c) stresslet component |S xx,2 |, and (d) stresslet component |S xy,2 | for the inert sphere "2" at a distance r from a puller squirmer (β = 5).

Fig. 7 .

 7 Fig. 7. Initial configuration of the squirmers in the trajectory simulations.

Fig. 8 .

 8 Fig.8. Trajectories of two squirmers swimming in opposite directions with transverse initial distance δ y = 1a, . . . , 10a. Lines: data from[START_REF] Ishikawa | Hydrodynamic interaction of two swimming model micro-organisms[END_REF]. Symbols: FCM results.

Fig. 9 .

 9 Fig. 9. Influence of the collision barrier parameters on the squirmer trajectories for the case where initially δ y = 1a. (a) Trajectories. Crosses: data from [28]. (b) Force barrier profiles (see also Appendix B). The line styles correspond those showing the trajectories in (a).: barrier used in Fig.8.

Fig. 11 .Fig. 12 .

 1112 Fig. 11. Polar order P (t) in a semi-dilute suspension (φ v = 0.1) of squirmer pullers (β = 1). (a) Time evolution of polar order depending on the number of swimmers. : L/a = 14, N p = 64; : L/a = 19, N p = 174; : L/a = 38, N p = 1395; : L/a = 58, N p = 4707; : L/a = 77, N p = 11,158; : L/a = 116, N p = 37,659. (b) Steady state value P ∞ depending on N p ∼ (L/a) 3 . : fit linear with 1/ N p . (Inset): Dependence

Fig. 13 .

 13 Fig. 13. Orientational distribution on the unit sphere (φ, θ). (a) Distribution before the transition to polar order: tU /a = 90. (b) Distribution once the polar ordered state is reached: tU /a = 1110.

Fig. 14 .

 14 Fig. 14. Time-averaged steady-state orientational distribution in the frame of the mean orientation vector (θ, φ)| p . (a) Distribution over the unit sphere. (b)

Fig. 15 .

 15 Fig. 15. Distribution of Voronoi cell volumes. (Main figure): Time dependence of the standard deviation, σ V , of the Voronoi volume distribution normalized by the mean value V V . (Inset): Time average of the Voronoi volume distribution before () and after ( ) the polar order transition.

Fig. 16 .

 16 Fig. 16. Steady state pair distribution function (a) g(r, θ). (b) g(r).

Fig. 17 .

 17 Fig. 17. Correlations in squirmer orientation at steady state. (a) I p (r, θ). (b) I p (r).

Fig. 18 .

 18 Fig. 18. Simulation of a dilute suspensions, φ v = 0.05, of prolate spheroidal pushers, B 2 = -1.5, B 1 = 0, with aspect ratio a1 a2 = a1 a3 = 3. (a) Snapshot at t = 30U /a 1 . Colors indicate the velocity magnitude normalized by the individual swimming speed. (b) Time evolution of polar order P (t).

:

  isotropic value of the polar order parameter 1/ N p = 0.026.

Fig. 19 .

 19 Fig. 19. (a) (B 1 (t), B 2 (t)) phase diagram for one beat cycle. (b) d (t) over one beat cycle.

Fig. 20 .

 20 Fig. 20. Snapshots of the time-dependent flow field around a model Chlamydomonas. Background grey levels represent the natural logarithm of the norm of the velocity field in rad s -1 . : position of the stagnation point given by the FCM model. : position of the stagnation point measured by [39]. (Insets) Swimming speed along the beat cycle.

Fig. 21 .

 21 Fig. 21. Relative trajectories of two model C. Rheinardtii swimming in opposite directions with stroke phase shift ϕ and initial transverse distance δ y = 1a and 2a. The grey levels lighten as ϕ increases from ϕ = 0 → ϕ = 7π /4.: trajectory of two steady squirmers with β = 0.

Fig. 22 .

 22 Fig. 22. Sketch of the simulations with time-dependent swimmers. Each swimmer n generates time-dependent flow disturbances with the same amplitude and frequency, but the phase ϕ n may be different for each swimmer.

Fig. 23 .

 23 Fig. 23. Polar order parameter and orientational distribution for suspensions of steady and time-dependent swimmers. : time-dependent swimmers with random phase; : synchronized time-dependent swimmers; : steady swimmers with β = 0; : steady swimmers with β = 0.1. (a) Time evolution of polar order. (b) Steady state orientation distribution around the mean director (θ)| p .: uniform distribution 0 (θ) = 1/π .

Table C . 1

 C1 Magnitudes of the Fourier modes used to describe B 2 (t) in rad s -1 .Phases of the Fourier modes used to describe B 2 (t) in rad.

	c 0	c 1	c 2	c 3	s 1	s 2	s 3
	4.5347	64.053	-84.7192	-10.4545	91.4529	-92.6420	-5.9849
	Table C.2
	ϕ c1	ϕ c2	ϕ c3	ϕ s1	ϕ s2		ϕ s3
	1.7373	3.5761	-0.9154	0.1666	2.0054	-1.7125

Acknowledgements

The authors would like to thank Prof. Martin Maxey for many insightful comments and acknowledge his instrumental role in bringing the authors together to carry out this work. We would also like to acknowledge Prof. Raymond Goldstein, Dr. Ottavio Croze, Navish Wadhwa and Prof. Pierre Degond for interesting discussions. We are also grateful to Annaig Pedrono and Dr. Pascal Fede for providing helpful support on numerical issues. This work is developed within the MOTIMO ANR Project. Simulations were performed on the Calmip supercomputing mesocenter. We thank INPT for funding the international collaboration between IMFT and Imperial College (grant no. SMI 2014). We would also like to acknowledge travel support from EPSRC Mathematics Platform grant EP/I019111/1.

Appendix A. Self-induced effects for spheroidal swimmers

In this appendix, we show how to compute the artificial self-induced effects due to the squirming modes (see Section 2.3) for the case of spheroidal swimmers described in Section 6.1.

In the following calculations, we consider an ellipsoidal squirmer located at the origin Y = 0 with swimming direction p = ê1 .

A.1. Self-induced velocity W due to the degenerate quadrupole H

Using the FCM envelopes for ellipsoidal particles, the force-distribution corresponding to the degenerate quadrupole generated by the squirmer will be given by f H = H∇ 2 ell (x), (A.1) where we also have σ = σ ;1 , σ ;2 , σ ;3 , the width of the Gaussian envelope ell (x) for each of the semi-axes. The resulting fluid flow can be obtained in Fourier space as

The self-induced velocity, W, will be given by volume averaging this fluid velocity against the monopole Gaussian envelope ell (x) for ellipsoidal particles. Accordingly, the expression for W is

where we write the width of the Gaussian envelope ell (x) for each semi-axes as σ = σ ;1 , σ ;2 , σ ;3 . This expression can be viewed as

is the FCM self-mobility matrix that relates the particle velocity to the degenerate quadrupole coefficient. For a spheroidal particle (σ ;2 = σ ;3 and σ ;2 = σ ;3 ) and H in the direction ê1 , we need only to consider M HU 11 to obtain the self-induced effects. This mobility matrix entry can be written as

and dz dt = r 1t 2 3/2 and Eq. (A.6) becomes

After integration by parts with respect to r, we have

and eventually

1 + γ 2 -1 t 2 3/2 dt.

(A.9)

From the integrals in the Appendix of [START_REF] Liu | Force-coupling method for flows with ellipsoidal particles[END_REF], this maybe be expressed compactly as

where C = 32π 3/2 µσ 3

and I 0 and I 1 are coefficients whose expressions are provided in the Appendix of [START_REF] Liu | Force-coupling method for flows with ellipsoidal particles[END_REF].

The self-induced velocity W is then simply [START_REF] Creppy | Collective motility of sperm in confined cells[END_REF] and its artificial effect can be subtracted away from the total volume average velocity as done in Eq. [START_REF] Saintillan | Orientational order and instabilities in suspensions of self-locomoting rods[END_REF].

A.2. Self-induced rate of strain K due to the squirming dipole G

In the frame of the swimmer, the non-zero entries of the self-induced rate of strain K ′ are given by

where C is given above and I 2 is a coefficient whose expression is detailed in the Appendix of [START_REF] Liu | Force-coupling method for flows with ellipsoidal particles[END_REF]. To obtain the selfinduced rate of strain in the lab frame, K, one just needs to apply a rotation operator onto K ′ :

where Q = ê1 ê2 ê3 is the rotation matrix of the swimmer.

Appendix B. Repulsive force parameters

Table B.1 lists the parameter values used in our squirmer trajectory calculations.