Modeling a spheroidal microswimmer and cooperative swimming in thin films

We propose a hydrodynamic model for a spheroidal microswimmer with two tangential surface velocity modes. This model is analytically solvable and reduces to Lighthill's and Blake's spherical squirmer model in the limit of equal major and minor semi-axes. Furthermore, we present an implementation of such a spheroidal squirmer by means of multiparticle collision dynamics simulations. We investigate its properties as well as the scattering of two spheroidal squirmers in a slit geometry. Thereby we find a stable fixed point, where two pullers swim cooperatively forming a wedge-like conformation with a small constant angle.


Introduction
Living matter exhibits a broad spectrum of unique phenomena which emerge as a consequence of its active constituents.Examples of such systems range from the macroscopic scale of flocks of birds and mammalian herds to the microscopic scale of bacterial suspensions. 1,2][10][11][12][13][14][15][16][17][18][19][20][21] The understanding of these collective phenomena requires the characterization of the underlying physical interaction mechanisms.][24] For micrometer-size biological unicellular swimmers, e.g., bacteria (E.][27][28][29][30] Generic models, which capture the essential swimming aspects, are crucial in theoretical studies of microswimmers.On the one hand, they help to unravel the relevant interaction mechanisms and, on the other hand, allow for the study of sufficiently large systems.A prominent example is the squirmer model introduced by Lighthill 31 and revised by Blake. 32Originally, it was intended as a model for ciliated microswimmers, such as Paramecia.Nowadays, it is considered as a generic model for a broad class of microswimmers, ranging from diffusiophoretic particles [33][34][35] to biological cells and has been applied to study collective effects in bulk, [36][37][38][39][40][41][42] at surfaces, 36,43,44 and in a narrow slit. 20n its simplest form, a squirmer is represented as a spherical rigid colloid with a prescribed surface velocity. 31,32,38Restricting the surface velocity to be tangential, the spherical squirmer is typically characterized by two modes accounting for its swimming velocity and its force-dipole.The latter distinguishes between pushers, pullers, and neutral squirmers.The assumption of a spherical shape is adequate for swimmers like Volvox, however, the shape of bacteria such as E. coli or the time-averaged shape of cells such as Chlamydomonas is nonspherical.Hence, an extension of the squirmer concept to spheroidal objects is desirable.In 1977, Keller and Wu proposed a generalization of the squirmer model to a prolate-spheroidal shape, which resembles lattice Boltzmann simulations, 36,41,50 the smoothed profile method, 42 and the force-coupling approach. 51In the following, we will apply the MPC method.MPC is a particle-based simulation technique which incorporates thermal fluctuations, [52][53][54] provides hydrodynamic correlations, 55,56 and is easily coupled with other simulation techniques such as molecular dynamics simulations for embedded particles. 53,54][59][60][61][62][63][64][65] Here, we implement our spheroidal squirmer model in MPC.More specifically, we study the resulting flow field and compare it with the theoretical prediction.Moreover, we present results for the cooperative swimming behavior of two spheroidal squirmers in a narrow slit.Two pullers exhibit a long-time stable configuration, where they swim together in a wedge-like conformation with a constant small angle due to the hydrodynamic interaction between the anisotropic squirmers as well as squirmers and walls.The cooperative and collective swimming motion of spheroidal squirmers in Stokes flow has been addressed in ref. 47 by an adopted boundary-element method.This approach neglects thermal fluctuations and tumbling of the squirmers completely; only hydrodynamic and excluded-volume interactions determine the squirmer motion.In contrast, our simulation approach includes thermal fluctuations, which affects the stability of the cooperative swimming motion due to the rotational diffusion of a spheroid.
2 Hydrodynamic model of a spheroidal squirmer

Spheroid geometry
We describe a nonspherical squirmer as a prolate spheroidal rigid body with a prescribed surface velocity u sq .In Cartesian coordinates (x, y, z), the surface equation of a spheroid, or ellipsoid of revolution, is with b z and b x the semi-major and semi-minor axis, respectively, and b z Z b x (cf.Fig. 2).We denote half of the focal length by , which yields the eccentricity e = c/b z .Furthermore, we define a swimmer diameter as s = 2b z .In terms of prolate (b z 4 b x ) spheroidal coordinates (z, t, j), the Cartesian coordinates are given by where À1 r z r 1, 1 r t r N, and 0 rj r 2p.All points with t = t 0 e À1 lie on the spheroid's surface.The intersection of the spheroid and a meridian plane, where j is constant, is an ellipse.The normal n and tangent s to this ellipse are given by the unit vectors e t and Àe z , respectively, which follow by partial derivative of eqn ( 2) with respect to the coordinates z and t.
For b x = b z , the spheroid becomes a sphere.The spherical coordinates (x, y, z) T = r(sin y cos j, sin y sin j, cos y) T (3) are obtained from eqn (2) for t -N, ct = r, and z = cos y.In this limit, the unit vectors turn into e te r and e z -Àe y (cf.Fig. 2).The Lame ´metric coefficients for prolate spheroidal coordinates , and

Flow field
The squirmer is immersed in an incompressible low-Reynoldsnumber fluid, which is described by the incompressible Stokes equations Here, v(r) is the fluid velocity field, p(r) the pressure field at the position r, and Z the viscosity.In an axisymmetric flow, The stream function itself satisfies the equation 66 with the operator 67 Each function C in the kernel of E 2 can be represented as 67 Cðt; zÞ ¼ with constants c i n and the functions Here, G n (x) and H n (x) are Gegenbauer functions of the first and second kind, respectively (see Appendix B).The velocity components follow from the stream function via 66 An important feature of a squirmer is the hydrodynamic boundary condition at its surface, which demands v(r) = u sq .For the squirming velocity u sq we propose Here, s is the tangent vector, e z = (0, 0, 1) T is the unit vector in z-direction, B 1 and B 2 are the two surface velocity modes, and b = B 2 /B 1 (cf.Fig. 2).In the swimmer's rest frame, and with eqn (12), the boundary value problem becomes Eqn ( 14) implies a constant background flow v = ÀU 0 e z infinitely far from the squirmer, eqn (15) guarantees v t = 0 at the spheroid surface, and eqn (16) demands v z = u sq (z)Áe z .Due to linearity of the Stokes stream function eqn (6), we can solve this boundary value problem for B 2 = 0 first, which yields the stream function C 1 .Subsequently we solve the problem Eqn ( 17) imposes a vanishing velocity field infinitely far from the squirmer, eqn (18) again guarantees v t = 0 at the spheroid surface, and eqn (19) demands v z = u sq (z,B 1 = 0)Áe z .We denote the solution of the problem eqn ( 17)-( 19) by C 2 .Finally, C = C 1 + C 2 solves the initial problem ( 14)-( 16) for arbitrary B 1 and B 2 .
The boundary value problem eqn ( 14)-( 16) for B 2 = 0 can be solved by the ansatz Here, the third term is found by the separation ansatz C(t,z) = g(t)(1 À z 2 ) for eqn (6).Eqn (14) directly yields a 1 = À2U 0 c 2 .The remaining coefficients a 2 and a 3 are determined by eqn (15) and (16), keeping in mind that B 2 = 0.This yields The boundary value problem eqn ( 17)-( 19) can be solved by the ansatz As before, the third term follows by a separation ansatz C(t,z) = g(t)z(1 À z 2 ) for eqn (6).Eqn (17) yields a 4 = 0.The coefficients a 5 and a 6 are determined by eqn (18) and (19) such that The total stream function C = C 1 + C 2 can be transformed to the laboratory frame (cf.Fig. 1) by adding the background flow v = U 0 e z , which yields The force on the spheroid by the fluid follows from a multipole expansion, 66,69 with a Stokeslet as the dominating contribution far away from a swimmer.Hence, for r -N the stream function C lab has to be equal to the stream function of a Stokeslet, 66 namely and, thus, the force on the spheroidal squirmer is given by 66 where As expected, C 2 does not contribute to the force, since it assumes a constant value at infinity.Since a swimmer must be force free, F z = 0, which implies a 3 = 0.Then, eqn (22) yields the swimming velocity of the squirmer (t 0 = 1/e) which was already found by Keller and Wu for the case B 2 = 0. 45 As a consequence, a 2 in eqn ( 22) simply becomes a 2 = 2B 1 c 2 t 0 (t 0 2 À 1).Examples of fluid velocity fields of a spheroidal squirmer are presented in Fig. 1 and 3.
Far field.The far field of a cylindrically symmetric microswimmer in terms of a multipole expansion is presented in ref.48 and the ESI of ref. 70.To obtain the far field expansion of our spheroidal squirmer, we expand the stream function, eqn (26), in powers of 1/t.Similarly, we determine the stream functions of the first few singularity solutions appearing in the multipole expansion (force dipole, force quadrupole, source dipole, rotlet dipole, etc.) and Taylor expand them in 1/t.Note that we can omit the non-axisymmetric singularity solutions of the multipole expansion like the rotlet dipole (v RD ), which is cylindrically symmetric but not axisymmetric (v RD Áe j a 0).Equating the coefficients of (1/t) n for n = 0, 1, 2, we find that the squirmer is well described in the far field by the flow fields of a force dipole, a source dipole, and a source quadrupole (30)   where which decay like r À2 , r À3 , and r À4 for large r, respectively. 48The multipole coefficients are with the coefficients a 5 and a 6 of eqn ( 24) and (25).The values of the multipole coefficients in the spherical limit (b zb x R, where R is the radius) follow from the above coefficients for t 0 -N, c -0, and ct 0 = R as k FD = ÀB 2 R 2 /2, k SD = ÀB 1 R 3 /3, and k SQ = B 2 R 4 /6 as expected for a spherical squirmer. 71

Multiparticle collision dynamics
Multiparticle collision dynamics (MPC) is a stochastic, particlebased mesoscale hydrodynamic simulation method. 54Thereby, a fluid is modeled by N point particles with equal mass m, undergoing subsequent streaming and collision steps.In the streaming step, the particle positions r i , i = 1,. ..,N, are updated according to where v i are the particle velocities and h is denoted as collision time step.In the subsequent collision step, the particle velocities are changed by a stochastic process, which mimics internal fluid interactions.In order to define the local collision environment, particles are sorted into cells of a cubic lattice with lattice constant a. Different realizations for this stochastic process have been proposed. 52,72,73We employ the stochastic rotation dynamics (SRD) approach of MPC with angular momentum conservation (SRD+a), 74,75 which updates the particle velocities in a cell according to Here, r i,c = r i À r cm , where r cm is the center-of-mass position of the particles in the cell, and similarly, v i,c = v i À v cm , with the center-of-mass velocity v cm .R(a) is the rotation matrix, which describes a rotation around a randomly oriented axis by the angle a.The angle a is a constant, and the axis of rotation is chosen independently for each cell and time step.Finally, I is the moment-of-inertia tensor of the particles in the center-ofmass reference frame of the cell.Partition of the system into collision cells leads to a violation of Galilean invariance.To reestablish Galilean invariance, a random shift of the collisioncell lattice is introduced at every collision step. 76,77ince energy is not conserved in the collision step, we apply a cell level canonical thermostat at temperature T. 78,79 The latter ensures Maxwell-Boltzmann distributed velocities.The MPC algorithm is embarrassingly parallel.Hence, we implement it on a Graphics Processing Unit (GPU) for a high performance gain. 80he following simulations are performed with the mean number of particles per collision cell hN c i = 10, the rotation angle a = 1301, and the time step , which yields a fluid viscosity of Z ¼ 17:8 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi mk B T=a 4 p .

Implementation of a spheroidal squirmer in MPC
A spheroidal squirmer is a homogeneous rigid body characterized by its mass M, center-of-mass position C, orientation q, translational velocity U, and angular momentum l.Thereby, q = (q 0 ,q 1 ,q 2 ,q 3 ) is a rotation quaternion and can be related to the rotation matrix D, which transforms vectors from the laboratory frame to the bodyfixed frame 81 (see Appendix A (eqn (64))).We distinguish vectors in the laboratory frame and body-fixed frame by a superscript, i.e., v s is a vector in the laboratory (or space-fixed) frame while is the corresponding vector in the body-fixed frame.For vectors in the laboratory frame, we will frequently omit the superscript. .When needed, the angular velocity is calculated as For all simulations we choose a neutrally bouyant spheroid, i.e., M = r(4p/3)b z b x 2 , where r is the fluid mass density.

Streaming step
During the streaming step, a spheroid will collide with several MPC particles.Since the total change in (angular) momentum of a spheroid during one streaming step is small, we perform the collisions with MPC particles in a coarse-grained way. 82or the streaming step at time t, we determine the spheroid's position, velocity, orientation, and angular velocity at times t + h/2 and t + h, under the assumption that there is no interaction with MPC particles.However, steric interactions between spheroids, as well as spheroids and walls are taken into account as described in Section 4.3.
Subsequently, all MPC particles are streamed, i.e., their positions are updated according to r i (t + h) = r i (t) + hv i (t).Thereby, a certain fraction of MPC particles penetrates a spheroid.To detect those particles in an efficient way, possible collision cells intersected by the spheroid are identified first.For this purpose, we select all those cells, which are within a sphere of radius b z enclosing the spheroid instead of the spheroid itself, which is more efficient, since it avoids rotating candidate cells into the body-fixed frame during selection.A loop over all particles in respective collision cells identifies those particles, which are inside the spheroid and they are labeled with the spheroid index.Then, each particle i inside a spheroid at time t + h is moved back in time by half a time step and subsequently translated onto the spheroid's surface.The translation can be realized in different ways.One possibility is to construct a virtual spheroid with semiaxes b ˜z, b ˜x, b ˜z/b ˜x = b z /b x and r i (t + h/2) on its surface.The particle is then translated along the normal vector of the virtual spheroid until it is on the real spheroid's surface.Alternatively, the difference vector r i (t + h/2) À C(t + h/2) can be scaled such that the particle position lies on the spheroid' surface.We tried both approaches and found no significant difference.Once the MPC particle at time t + h/2 is located on the spheroid's surface, the momentum transfer due to a bounce-back collision at time t + h/2 is determined, taking into account the squirmer surface fluid velocity u sq of eqn (11). 83Thereby, a useful identity to determine s is given in eqn (8) of ref. 45, and z is given by The velocity of the MPC particle is updated according to v i 0 = v i À J i /m.Subsequently, the position r i (t + h) is obtained by streaming the MPC particle for the remaining time h/2 with velocity v i 0 , i.e., As a consequence of the elastic collisions, the center-of-mass velocity and rotation frequency of a spheroid are finally given by where J ¼ P i J i is total momentum transfer by the MPC fluid and Þ Â J i is the respective angular momentum transfer.

Collision step
In a first step, ghost particles are distributed inside each spheroid. 82,84The number density and mass are equal for ghost and fluid particles.The ghost particle positions r g i are uniformly distributed in the spheroid and their velocities are given by The Cartesian components of v R,i are Gaussian-distributed random numbers with zero mean and variance ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k B T=m p .The squirming velocity u sq,i is determined by eqn (11), with the ghost particle position projected onto the spheroid's surface (cf.Section 4.1).As a result of MPC collisions, a spheroid's linear and angular momenta change by where % v g i and v g i are the ghost particle's velocity after and before the MPC collision.Hence, the spheroid velocity and angular velocity become

Rigid body dynamics for spheroids
During the streaming step, the spheroids move according to rigid-body dynamics, governed by 85 MC ¨= F, Here, Q(q) is defined in Appendix A (eqn (68)) and F and T are the force and torque acting on the spheroid.Forces and torques are derived from steric interaction potentials as presented in Appendix C. Eqn (50) are Euler's equations for rigid body dynamics and hold for (a, b, g) = (x, y, z), (y, z, x), and (z, x, y).Whenever necessary, body-fixed and laboratory-frame quantities can be related by the rotation matrix D which is given in terms of the quaternion q in Appendix A (eqn (64)).
For the numerical integration of the equations of motion, the widely applied leap-frog method 81 is not useful, since velocity, angular momentum, position, and orientation are required at the same point in time for the coupling to the MPC method.Hence, we employ the Verlet algorithm for rigid-body rotational motion proposed in ref. 85.Integration for a time step t is performed as follows: (i) Update C and q according to (cf.eqn ( 49) and ( 50) The parameter l is introduced to guarantee q 2 = 1.
(ii) Calculate forces and torques F s (t + t) and T s (t + t).
(iii) Update U and l s according to 5 Simulations -thermal properties and flow field

Passive colloid
For the passive spheroidal colloid (B 1 = B 2 = 0), we perform equilibrium simulations and determine hU a 2 i as well as h(O b a ) 2 i for a A {x, y, z}.Due to the equipartition of energy, we expect i is larger for O b a than for U a .We find the largest relative error for h(O b z ) 2 i, namely s r = 9.5%, 5.3%, and 3.1% for b x = 2a, 3a, and 4a.Hence, we choose the minor axis b x Z 3a in the following.
In addition, we determine the orientation correlation function he(t)Áe(0)i.The theory of rotational Brownian motion 86 predicts where x > , and x J and x > are the parallel and perpendicular rotational friction coefficients of a prolate spheroid with respect to the major semi-axis; explicitly 69 Simulation results for the orientational auto-correlation function are shown in Fig. 4 for two spheroids of different eccentricity.The correlation functions decay exponentially.However, for the spheroid with the smaller eccentricity, we find a somewhat faster decay than predicted by theory, whereas good agreement is found for the larger spheroid.We attribute the difference to finite-size effects related to the discreteness of the collision lattice.For larger objects, discretization effects become smaller.
This journal is © The Royal Society of Chemistry 2016

Squirmer
We determine the steady state swimming velocity of a squirmer via heÁUi, which should be equal to U 0 (cf.eqn (3)).Results for various eccentricities are displayed in Fig. 5.The velocity U 0 increases with increasing eccentricity e in close agreement with the theoretical prediction of eqn (3).We confirm that the forcedipole parameter b does not affect the velocity of the squirmer, as long as the Reynolds number Re is low, i.e., Re = rU 0 b z /Z t 0.1.We also determine the orientational correlation function and find that a squirmer exhibits the same orientational decorrelation as the corresponding passive particle (cf.Fig. 4).Moreover, we calculate the flow field from the simulation data and compare it with the theoretical prediction.As shown in Fig. 6, the two fields are in close agreement.The twodimensional flow field of the MPC fluid, averaged over the rotation angle j, is determined at the vertices of a fine resolution mesh.The velocities at these vertices include averages over time of an individual realization as well as ensemble averages over various realizations.By the latter, we determine an estimate for the error of the mean velocity.The median (over vertices) of this error is approximately 5% for the parameters of Fig. 6(b) and 10% for that of Fig. 6(e).Note that we choose a smaller swimming mode B 1 for the puller (Fig. 6(b)) than for the neutral squirmer (Fig. 6(e)).The reason is that the agreement with theory was not satisfactory for the puller with B 1 ¼ 0:05 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k B T=m p , which we attribute to a non-vanishing Reynolds number (Re E 0.1) in the simulation.In Fig. 6(c) and (f), we observe lines of high relative errors (yellow in the color code).They appear because theory predicts v % r = 0 or v z = 0 for these lines, which is difficult to achieve in simulations.Hence, the overall agreement between simulations and theory is very satisfactory, and the implementation is very valuable for the simulation of squirmersquirmer and squirmer-wall interactions, where the details of the flow field matter.For a benchmark of the code on a current GPU see Appendix D.

Cooperative swimming in a narrow slit
We simulate the cooperative swimming behavior of two squirmers in a slit geometry.The slit is formed by two parallel no-slip walls located at y = 0 and y = L y .The no-slip boundary condition is implemented by applying the bounce-back rule and ghost particles of zero mean velocity in the walls. 84Steric interactions between two squirmers and between a squirmer and a wall are taken into account by the procedure described in Appendix C. The plot shows the simulation data (blue and black solid lines), an exponential fit to that data (red dashed), and the theoretical prediction according to eqn (58) (green dotted).The initial positions and orientations of the two squirmers (i = 1, 2) are e 1/2 = (AEcos(a 0 ), 0, sin(a 0 )) T .( 63) Here, d cm is the initial center-of-mass distance and a 0 = (p À y 0 )/2, where y 0 is the initial angle between e 1 and e 2 .The swimming mode is chosen as B 1 ¼ 0:05 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k B T=m p and the force dipole mode b A {À4,0,4}.We choose d cm such that the squirmers are well separated and vary y 0 .The squirmers major and minor axes are b x = 3a and b z = 6a, respectively, and the simulation box size is L x = L z = 15b z , and L y = 7a.Note that L y \ b x which keeps the swimming orientation essentially in the x-z plane.
Results for the mean surface-to-surface distance between squirmers hd s i and the mean alignment he 1 Áe 2 i = hcos yi are shown in Fig. 7 for pushers, pullers, and neutral swimmers with an initial angle y 0 = 3p/8.Due to the setup, the squirmers initially approach each other and collide at tU 0 /s E 0.5.The (persistence) Pe ´clet number Pe = U 0 /(2D > R s) E 60 is sufficiently high, such that the squirmer orientation has hardly changed before collision.When the neutral swimmers collide, they initially align parallel (cos y E 1 at tU 0 /s E 1 in Fig. 7), but their trajectories start to diverge immediately thereafter.Pushers remain parallel for an extended time window, which is expected as pushers are known to attract each other, 37 but at tU 0 /s E 3 (cf.Fig. 7) their trajectories diverge as well.This is probably due to noise, since we observe several realizations where pushers remain parallel.Interestingly, pullers, which are known to repel each other when swimming in parallel, 37 swim cooperatively and reach a stable orientation with hcos(y)i E 0.77 shortly after they collided (at tU 0 /s E 1).Thereby, their cooperative swimming velocity is about 0.8U 0 .The flow field of this stable state, determined by MPC simulations, is shown in Fig. 8.Note that the velocity field in the swimming plane is left-right symmetric, and that there is a stagnation point in the center behind the swimmers.Fig. 8 reveals that this point actually corresponds to a line normal to the walls.In ref. 87, it was shown that the flow field of a force dipole in a narrow slit exhibits a recirculating pattern with loops in a plane parallel to the walls 88 and a parabolic flow profile perpendicular to them.Both features can be observed in Fig. 8. Fig. 9 shows that the fixed point of cooperatively swimming pullers is reached for nearly all simulated  This journal is © The Royal Society of Chemistry 2016 initial conditions y 0 A (0, p/2).Only pullers that are nearly parallel initially (y 0 = p/8, cos y 0 E 0.92 in Fig. 9), repel each other such that they will not reach the fixed point.For Pe ´clet numbers Pe o 60, the fixed point remains at hcos(y)i E 0.77.However, it becomes more likely for the swimmers to escape (or never reach) the fixed point.
A detailed study reveals that the fixed point vanishes, when the walls are replaced by periodic boundary conditions (we use a cubic simulation box of length 10b z ).This is even true when we apply three-dimensional periodic boundaries, but keep the wall potential implemented, i.e., the squirmers are still confined in a narrow slit of height 7a as before (the fluid simulation box is cubic with length 10b z ).These observations are quantified in Fig. 10.Note that the slow increase of the average squirmer surface-to-surface distance in the case of a periodic system with wall potential is due to the that the pullers still swim together in several realizations.However, in these realizations they no longer swim in a plane as for the slit geometry, but rather show three-dimensional configurations with tilted major axes.Evidently, hydrodynamic interactions with the walls 89 are essential for the swimmers' dynamics.
In addition, we studied the swim behavior of spherical squirmers.Here, we observe diverging trajectories for all squirmer types, i.e., pushers, neutral squirmers, and pullers.Such diverging trajectories have already been reported in ref. 90 for spherical squirmers in bulk.In contrast, in ref. 91 a cooperative swimming mode for spherical squirmers has been observed (see Fig. 22(c) of ref. 91).However, this cooperative swimming-termed pairswimming by the authors-is unstable to perturbations that displace one swimmer out of the swimming plane. 91Since our simulations and those of ref. 90 include thermal fluctuations, we consequently do not observe the cooperative swimming mode of ref. 91.
Hence, the stable close-by cooperative swimming of pullers is governed by the squirmer anisotropy, by the hydrodynamic interactions between them and, importantly, between pullers and confining surfaces.
This conclusion is in contrast to results presented in ref. 47, where a monolayer of spheroidal squirmers is considered, with their centers and orientation vectors fixed in the same plane, however, without confining walls.The study reports a stable cooperative motion for pullers with angles y A (0, p/2) by nearestneighbor two-body interactions, where all angles between 0 and p/2 are stable.The difference to our study is that in ref. 47 cooperative features were extracted from a simulation of many swimmers, whereas we explicitly studied two swimmers.Furthermore our study explicitly models no-slip walls and includes thermal noise.
To shed light on the stability of the cooperative puller motion, we varied the puller strength b, the aspect ratio b z /b x , and the width of the slit L y .Thereby, we started from our basic parameter set b x = 3a, b z = 6a, L y = 7a, and b = 4.With decreasing b, the stable alignment disappears, i.e., the pullers' distance increases after collision (see ESI, † Movies beta 3 and beta 4).For increasing b the fixed point remains, but the value of cos y decreases, i.e., the squirmers form a larger angle.With increasing wall separation, the fixed-point value of cos y decreases, i.e., the angle between the swimmers increases.For L y /b z \ 11/6, the mean distance between swimmers increases rapidly after collision (see Fig. 10).An increase of the aspect ratio b z /b x from 2 to 3 and 4 increases the fixed-point value of hcos yi from 0.77 to 0.84 and 0.88.The more elongated shape leads to a more parallel alignment of the squirmers.For b z /b x Z 2, the minimal value of b required to achieve cooperative motion depends weakly on the aspect ratio.In particular, for b z /b x = 2, 3, and 4, we find the minimal value b E 3.5.

Summary and conclusions
We have introduced a spheroidal squirmer model, which comprises the swimming and force-dipole modes.It is a variation of previously proposed squirmer models.On the one hand, it includes the force-dipole mode as an extension to the model of ref. 45.On the other hand, it is an alternative approach compared to ref. 44 and 47, with the major advantage that our model allows for the analytical calculation of the flow field.In the present calculations we employed the Stokes stream function equation.Very recently a full set of solutions to Stokes' equations in spheroidal coordinates were given in ref. 92, which opens an alternative approach to derive the flow field for our choice of boundary conditions.Furthermore, we have presented an implementation of our spheroidal squirmer in a MPC fluid.In contrast to other frequently employed simulation approaches, MPC includes thermal fluctuations.The comparison between the fluid flow profile of a squirmer extracted from the simulation data with the theoretical prediction yields very good agreement.As a consequence of the MPC approach with its discrete collision cells, the minor axis of the spheroid has to be larger than a few collision cells to avoid discretization effects.The analysis of the squirmer orientation correlation function shows that very good agreement between theory and simulations is already obtained for b z = 9a (major axis) and b x = 3a (minor axis).
To shed light on the cooperative swimming motion and on near-field hydrodynamic interactions, we investigated the collision of two spheroidal squirmers in a slit geometry.We found a stable stationary state of close-by swimming for spheroidal pullers, which is determined by hydrodynamic interactions between the anisotropic squirmers, and, even more important, by squirmers and surfaces.This stationary state disappears for low puller strengths and low eccentricities.We expect the stable close-by swimming of pullers to strongly enhance clustering in puller suspensions in narrow slits.
Our studies confirm that spheroidal squirmers can accurately be simulated by the MPC method.The proposed implementation opens an avenue to study collective and non-equilibrium effects in systems of anisotropic microswimmers.Even large-scale systems (10 3 to 10 4 swimmers) can be addressed by the implementation of MPC and the squirmer dynamics on current GPUs.

A Quaternion matrices
The rotation matrix D introduced in eqn ( 39) is given in terms of the rotation quaternion q as The matrix Q(q) in eqn ( 49) is given by QðqÞ ¼ q 0 Àq 1 Àq 2 Àq 3 q 1 q 0 Àq 3 q 2 q 2 q 3 q 0 Àq 1

B Gegenbauer functions
For n Z 2 and x A R the Gegenbauer functions of the first and second kind G n and H n , are defined in terms of the Legendre functions of the first and second kind P n and Q n as 67,93 G n ðxÞ ¼ P nÀ2 ðxÞ À P n ðxÞ 2n À 1 For n = 0, 1, they are defined as For the reader's convenience, we give the formula for the Gegenbauer functions of the first kind for n = 2, 3, and Furthermore, the Gegenbauer functions of the second kind for n = 2, 3, and x 4 1 are given by Here, we used coth

C Steric interactions
Here, we illustrate our implementation of the excluded-volume interactions between spheroids and walls following the approach provided in ref. 94.The spheroid's surface in the laboratory frame is given by the quadratic form where the orientation matrix A can be expressed as C.1 Interaction between spheroids.We introduce a repulsive interaction potential between spheroids to prevent their overlap.The potential is given by Here, s 0 and e 0 correspond to a length and energy scale, respectively.We choose e 0 = k B T and s 0 = 2d v .The directional contact distance d R between two spheroids, with orientation matrices A 1 , A 2 and center positions C 1 , C 2 , is an approximation to their true distance of closest approach and is defined by This journal is © The Royal Society of Chemistry 2016 Here, R = C 2 À C 1 , R = |R|, and F(A 1 ,A 2 ) is the elliptic contact function, defined as 94 Minimization with respect to x demands =S(x,l) = 0, and hence, The critical value l = l c that maximizes S(x(l),l) can be found by the root finding problem We implement Brent's root finding approach. 95The forces and torques arising from the potential (74) can be calculated analytically and are given by ‡ 94 and for the first spheroid, where X c = 2l c A 1 (x c À C 1 ).The force and torque on the second spheroid follow by Newton's actionreaction law, namely We restrict ourselves to short-range repulsive interactions by setting the potential U to a constant value for d R 4 ffiffi ffi 2 6 p À 1 À Á s 0 , implies that F 1 and T 1 are zero for this range of d R values.Note that an upper bound to d R is R À 2b z , which means that two spheroids will not interact C.2 Interaction between a spheroid and a wall.We assume that two parallel walls are positioned at y = 0, L y , which-taking into account the safety distance d v -results in the effective wall positions y = d v and L y À d v .We propose an interaction between a spheroid and a wall in the style of the spheroid-spheroid interaction presented in ref. 94.First, we find the point x on the spheroid's surface that is closest to a wall.For the wall at y = d v , this is achieved by minimizing the height h(x) = e y Áx À d v under the constraint A(x) = 1.Using the method of Lagrange multipliers, we have to minimize L(x,l) = h(x) + l(A(x) À 1).The necessary condition for a minimum @L/@x = 0 yields e y + l=A(x) = e y + 2lA(x À C) = 0, (85)   and hence, Substitution of eqn (86) into A(x) = 1 yields Finally, we obtain the point closest to the wall as Here, the minus sign has to be chosen, which can be visualized by the example of a sphere of radius R, for which A = R À2 1.This finally yields the height We employ the Lennard-Jones potential for a repulsive wall, and U w assumes a constant value for all h !ffiffi ffi 2 6 p À 1 À Á s 0 .We can derive the force F a = À@U w /@C a and torque T a = À@U w /@c a acting on the spheroid analytically.For the force, we find F ¼ À @U w @h @h @C y e y (91) and for the torque with Here, we used the relation which holds for an invertible matrix B = B(t) depending on a scalar parameter t, and eqn (C9) from ref. 94.
For the wall at y = L y À d v , we have to minimize h(x) = L y À d v À e y Áx, with x on the spheroid's surface.This yields The formulas for torque and force do not change, except that we have to insert h ¼ L y À d v À C y À ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi A À1 ð Þ yy q and need to change the sign of the force.

D Benchmark simulation
Computation time.We perform benchmark simulations on an NVIDIA K80 GPU.First we determine the computation time for one MPC step in a periodic system, normalized by the number of fluid particles.We obtain approximately 4.5 nanoseconds per time step and particle, which is in accordance with the value of about 3 nanoseconds reported in ref. 80, since the extension to angular momentum conserving MPC roughly doubles the computation time, while the NVIDIA K80 is approximately twice as fast as the NVIDIA GTX580 employed in ref. 80.
Next, we perform a benchmark MPC simulation of many spheroidal squirmers with b x = 3a, b z = 6a in a narrow slit.The slit height and length are chosen as L y = 7a and L x = L z = 600a, respectively.We vary the number N sq of squirmers, i.e., the twodimensional packing fraction N sq pb x b z /(L x L z ), and determine the computation time per MPC time step, normalized by the number of fluid and ghost particles.A simulation without squirmers takes approximately 4.9 nanoseconds per time step and particle.Hence, the addition of no-slip walls leads only to a minor increase in simulation time by approximately 10%.When spheroidal squirmers are introduced, the computation time increases linearly with the number of squirmers up to a value of 13.2 nanoseconds for the two-dimensional packing fraction 0.6, which corresponds to 3825 squirmers in our set-up.Similar values are found for smaller (L x = L y = 300a) and larger systems (L x = L y = 900a).The increase in computation time is mostly related to the squirmers' ghost particles.Unlike the fluid particles and the ghost particles for the no-slip walls, they are not spatially ordered in memory, which leads to a significantly lower processing speed (cf.Fig. 3

of ref. 80).
The steric interactions are computed sequentially on the CPU exploiting a cell-linked list. 81The computation time spent for the steric interactions increases linearly with N sq .However, it only contributes by about 4% to the total simulation time for the considered set-up.
Memory limitations.The size of the systems that can be studied with our code is limited by the available memory on the used GPU.Explicitly, the required memory to store fluid particle, collision cell, and spheroid properties are listed in Table 1.
Compared to the number of fluid particles, the number of spheroidal squirmers is several orders of magnitude smaller and their impact on GPU memory usage is negligible.With a typical density of hN c i = 10 particles per cell, the required GPU memory is 680 bytes per collision cell.Therefore, approximately 1.5 Â 10 6 collision cells can be studied per 1 GB of GPU memory.Hence, the memory capacity of 4-12 GB of recent GPUs corresponds approximately to 6 Â 10 6 -2 Â 10 7 collision cells.The volume of a spheroidal squirmer with b x = 3a, b z = 6a is equivalent to 226 cells.Assuming a 30% volume fraction of squirmers, memory limits the total number of squirmers to 8 Â 10 3 -2.4Â 10 4 .Note that for systems of many squirmers, the time scales of interest can be very long.Hence, in applications, the considered system size has often to be smaller (E10 3 squirmers) for acceptable total simulation times.

Fig. 1
Fig. 1 Flow field of a spheroidal puller with b = 3, (a) in the laboratory frame, and (b) in the body-fixed frame.The magnitude of the velocity field is color coded logarithmically.

Fig. 2
Fig. 2 Sketch of normal and tangent vectors of a spheroidal (left) and spherical (right) squirmer.In the squirmer model, self-propulsion (in z-direction) is achieved by a prescribed tangential surface velocity in direction of the tangent vector s.
B 1 determines the swimming velocity, while the B 2 term introduces a force-dipole, or pusher (B 2 o 0) and puller (B 2 4 0) mode.Note that the spherical squirmer introduced by Lighthill and Blake with modes B 1 and B 2 31,32 is recovered for the spherical limit of a spheroid, where zcos(y) = nÁe z .For B 2 = 0, this model of a spheroidal squirmer was already introduced and analysed in ref. 45 and 68.An additional force-dipole mode has been introduced in ref. 44 and 47 as u sq (z) = ÀB 1 sÁe z (1 + bnÁe z )s.However, we prefer the squirming velocity introduced in eqn (12), since it yields an analytically solvable boundary value problem for the Stokes equation.The two approaches provide a somewhat different flow field in the vicinity of the squirmer, but both yield the model of Lighthill and Blake in the limit of zero eccentricity.

Fig. 3
Fig. 3 Fluid velocity fields of a spheroidal squirmer in the laboratory frame for (a) B 1 = 1, B 2 = 0, and (b)B 1 = 0, B 2 = 1.The corresponding stream function is given by eqn(26).The magnitude of the velocity field is color coded logarithmically.Note that the pusher velocity field with B 1 = 0, B 2 = À1 is not shown, since it follows from that of the puller with B 1 = 0, B 2 = 1 by inverting the arrows.

b z 2 )
The orientation vector of a spheroid is e = D T e b = D T (0,0,1) T .The moment of inertia tensor in the body-fixed frame I b is a constant diagonal matrix with diagonal elements I x = (M/5)(b x 2 + = I y and I z = (2M/5)b x 2

2 i À hx sim 2 i)/hx theo 2
) We fix the aspect ratio b z /b x = 2 and vary b x in the range b x A [2a, 4a].The simulation results agree very well with the theoretical values (56) and(57).As expected, the deviations from theory decrease with increasing spheroid size, due to a better resolution in terms of collision cells.In general, the relative error s r = (hx theo

Fig. 4
Fig. 4 Orientation correlation functions he(t)Áe(0)i for passive spheroids with b z = 6a, b x = 3a (bottom blue line) and b z = 9a, b x = 3a (top black line).The plot shows the simulation data (blue and black solid lines), an exponential fit to that data (red dashed), and the theoretical prediction according to eqn (58) (green dotted).

Fig. 5
Fig. 5 Mean swimming velocity as function of the eccentricity e for a spheroidal squirmer with B 1 ¼ 0:05 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k B T=m p and B 2 = 0.The solid line shows the theoretical prediction of eqn (3).Black dots are simulation results.The eccentricity was varied by changing b z and keeping b x = 3a constant.For the red triangle, we simulated a larger spheroid with b x = 6a, which shows a better agreement with theory.

Fig. 6
Fig. 6 Fluid flow fields of a spheroidal squirmer in the laboratory frame with b x = 3a, b z = 6a, B 1 ¼ 0:01 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k B T=m p , and b = 3 (a-c), and with B 1 ¼ 0:05 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k B T=m p ; b ¼ 0 (d-f).The magnitude of the velocity field (in units of ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k B T=m p ) is color coded logarithmically.The plots (a) and (d) show theoretical results, (b) and (e) simulation results, and (c) and (f) relative errors.The relative error of the flow field is defined as Dv a = |v theo a À v sim a |/[(|v theo a | + |v sim a |)/2].Note, due to the discrete representation of the velocity field, some streamlines end abruptly.

Fig. 7
Fig. 7 Average surface-to-surface distance d s and orientation of squirmers, where cos(y) = e 1 Áe 2 , as function of time.The solid blue, dashed black, and dotted red lines correspond to pullers b = 4, neutrals b = 0 and pushers b = À4.The standard deviation of the blue line (b = 4) is indicated by the cyan shaded region.

Fig. 8
Fig. 8 Flow field of two cooperatively swimming pullers in the laboratory frame.The magnitude of the velocity field (in units of ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k B T=m p ) is color coded logarithmically.We denote the direction normal to the walls by y, the cooperative swimming direction by z, the remaining Cartesian axis by x, and choose the swimmers' center of mass as origin.Panel (a) shows the flow field at x = 0 in the zy-plane, while panel (b) shows the flow field at y = 0 in the xz-plane.The black elliptical shapes (''transparent'', solid) indicate the projection of the swimmers onto the considered plane.

Fig. 9
Fig. 9 Time dependence of the average alignment he 1 Áe 2 i = hcos yi of two pullers with b = 4, b z /b x = 2 in a slit of height L y = 7a, and various initial angles y 0 A (0,p/2) (see ESI, † Movie beta 4).

Fig. 10
Fig. 10 Average surface-to-surface distance d s between two pullers with b = 4, b z = 6a, and b x = 3a as a function of time for various geometries.The distance d s increases rapidly after collision for squirmers in a cubic simulation box of side length 10b z with periodic boundary conditions (black dotted line).Trajectories still diverge, when a confining potential for the squirmers is present (red dashed line); for details see text.Squirmers confined in a slit of width L y = 10a with no-slip walls swim together (blue solid line).However, their trajectories diverge for wide slits (L y = 11a for the green dash-dotted line and L y = 12a for the cyan dash-dot-dotted line).

)
For the steric interactions, we introduce a virtual safety distance d v , which is small compared to b x and b z .When computing steric interactions, we replace b x and b z by b x + d v and b z + d v , respectively.In this paper we used d v = 0.05a for all simulations.