Modeling a spheroidal microswimmer and cooperative swimming in a narrow slit

Abstract

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 particle-based mesoscale hydrodynamics simulations using the multiparticle collision dynamics approach. 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.

---

1 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 Specifically, active systems exhibit remarkable nonequilibrium phenomena and emergent behavior like swarming,^3–7 turbulence,⁶ and activity-induced clustering and phase transitions.^8–21 The understanding of these collective phenomena requires the characterization of the underlying physical interaction mechanisms. Experiments and simulations indicate that shape-induced interactions, such as inelastic collisions between elongated objects or of active particles with surfaces lead to clustering, collective motion, and surface-induced aggregation.^6,22–24 For micrometer-size biological unicellular swimmers, e.g., bacteria (E. coli), algae (Chlamydomonas), spermatozoa, or protozoa (Paramecium), hydrodynamic interactions are considered to be important for collective effects and determine their behavior adjacent to surfaces.^1,25–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³¹ and revised by Blake.³² Originally, 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–35 to biological cells and has been applied to study collective effects in bulk,^36–42 at surfaces,^36,43,44 and in a narrow slit.²⁰

In its simplest form, a squirmer is represented as a spherical rigid colloid with a prescribed surface velocity.^31,32,38 Restricting 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 real biological microswimmers such as Tetrahymenapyriformis, Spirostomum ambiguum, and Paramecium multimicronucleatum.⁴⁵ However, that squirmer model accounts for the swimming mode only and does not include a force-dipole mode. This is unfortunate, since the force-dipole mode determines swimmer–swimmer and swimmer–wall interactions.^25,37,39,46 A route to incorporate the force-dipole mode into the spheroidal squirmer model was proposed in ref. 44. However, to the best of our knowledge, the resulting hydrodynamic model is not solvable analytically. In this article, we propose an alternative model for a spheroidal squirmer, taking into account both, a swimming and a force-dipole mode. The major advantage of our approach is that the flow field can be determined analytically (cf.Fig. 1).

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

Various mesoscale simulation techniques have been applied to study the dynamics of squirmers embedded in a fluid, comprising Stokesian dynamics,^39,40,43 the boundary-element method,^{38,44,46–48} the multiparticle collision dynamics (MPC) approach,^20,37,49 lattice Boltzmann simulations,^36,41,50 the smoothed profile method,⁴² and the force-coupling approach.⁵¹ In the following, we will apply the MPC method. MPC is a particle-based simulation technique which incorporates thermal fluctuations,^52–54 provides hydrodynamic correlations,^55,56 and is easily coupled with other simulation techniques such as molecular dynamics simulations for embedded particles.^53,54 The method has successfully been applied in various studies of active systems underlining the importance of hydrodynamic interactions for microswimmers.^{1,20,24,28,37,53,57–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

2.1 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

(x² + y²)/b_x² + z²/b_z² = 1, (1)

with b_z and b_x the semi-major and semi-minor axis, respectively, and b_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 σ = 2b_z. In terms of prolate (b_z > b_x) spheroidal coordinates (ζ, τ, φ), the Cartesian coordinates are given bywhere −1 ≤ ζ ≤ 1, 1 ≤ τ ≤ ∞, and 0 ≤φ ≤ 2π. All points with τ = τ₀ ≡ e⁻¹ lie on the spheroid's surface. The intersection of the spheroid and a meridian plane, where φ is constant, is an ellipse. The normal n and tangent s to this ellipse are given by the unit vectors e_τ and −e_ζ, respectively, which follow by partial derivative of eqn (2) with respect to the coordinates ζ and τ. For b_x = b_z, the spheroid becomes a sphere. The spherical coordinates

(x, y, z)^T = r(sin θ cos φ, sin θ sin φ, cos θ)^T (3)

are obtained from eqn (2) for τ → ∞, cτ = r, and ζ = cos θ. In this limit, the unit vectors turn into e_τ → e_r and e_ζ → −e_θ (cf.Fig. 2). The Lamé metric coefficients for prolate spheroidal coordinates are , , and .

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.

2.2 Flow field

The squirmer is immersed in an incompressible low-Reynolds-number fluid, which is described by the incompressible Stokes equations

ηΔv − ∇p = 0, ∇·v = 0. (4)

Here, v(r) is the fluid velocity field, p(r) the pressure field at the position r, and η the viscosity. In an axisymmetric flow, the velocity field can be expressed by the stream function Ψ as⁶⁶ The stream function itself satisfies the equation⁶⁶

E⁴Ψ = 0, (6)

with the operator⁶⁷ Each function Ψ in the kernel of E² can be represented as⁶⁷with constants cⁱ_n and the functions

Θ¹_n(τ,ζ) = G_n(τ)G_n(ζ), Θ²_n(τ,ζ) = G_n(τ)H_n(ζ),

Θ³_n(τ,ζ) = H_n(τ)G_n(ζ), Θ⁴_n(τ,ζ) = H_n(τ)H_n(ζ).

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⁶⁶

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

u_sq = −B₁(s·e_z)s − B₂ζ(s·e_z)s (11)

= −B₁(1 + βζ)(s·e_z)s (12)

Here, s is the tangent vector, e_z = (0, 0, 1)^T is the unit vector in z-direction, B₁ and B₂ are the two surface velocity modes, and β = B₂/B₁ (cf.Fig. 2). B₁ determines the swimming velocity, while the B₂ term introduces a force-dipole, or pusher (B₂ < 0) and puller (B₂ > 0) mode. Note that the spherical squirmer introduced by Lighthill and Blake with modes B₁ and B₂^31,32 is recovered for the spherical limit of a spheroid, where ζ → cos(θ) = n·e_z.

For B₂ = 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(ζ) = −B₁s·e_z(1 + βn·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.

In the swimmer's rest frame, and with eqn (12), the boundary value problem becomes

Ψ(τ₀,ζ) = 0 for all ζ, (15)

Eqn (14) implies a constant background flow v = −U₀e_z infinitely far from the squirmer, eqn (15) guarantees v_τ = 0 at the spheroid surface, and eqn (16) demands v_ζ = u_sq(ζ)·e_ζ. Due to linearity of the Stokes stream function eqn (6), we can solve this boundary value problem for B₂ = 0 first, which yields the stream function Ψ₁. Subsequently we solve the problem

Ψ(τ,ζ) converges for τ → ∞, (17)

Ψ(τ₀,ζ) = 0 for all ζ, (18)

Eqn (17) imposes a vanishing velocity field infinitely far from the squirmer, eqn (18) again guarantees v_τ = 0 at the spheroid surface, and eqn (19) demands v_ζ = u_sq(ζ,B₁ = 0)·e_ζ. We denote the solution of the problem eqn (17)–(19) by Ψ₂. Finally, Ψ = Ψ₁ + Ψ₂ solves the initial problem (14)–(16) for arbitrary B₁ and B₂.

The boundary value problem eqn (14)–(16) for B₂ = 0 can be solved by the ansatz

Ψ₁(τ,ζ) = α₁G₂(τ)G₂(ζ) + α₂H₂(τ)G₂(ζ) + α₃τ(1 − ζ²). (20)

Here, the third term is found by the separation ansatz Ψ(τ,ζ) = g(τ)(1 − ζ²) for eqn (6). Eqn (14) directly yields α₁ = −2U₀c². The remaining coefficients α₂ and α₃ are determined by eqn (15) and (16), keeping in mind that B₂ = 0. This yields The boundary value problem eqn (17)–(19) can be solved by the ansatz

Ψ₂(τ,ζ) = α₄G₃(τ)G₃(ζ) + α₅H₃(τ)G₃(ζ) + α₆ζ(1 − ζ²). (23)

As before, the third term follows by a separation ansatz Ψ(τ,ζ) = g(τ)ζ(1 − ζ²) for eqn (6). Eqn (17) yields α₄ = 0. The coefficients α₅ and α₆ are determined by eqn (18) and (19) such that The total stream function Ψ = Ψ₁ + Ψ₂ can be transformed to the laboratory frame (cf.Fig. 1) by adding the background flow v = U₀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 → ∞ the stream function Ψ^lab has to be equal to the stream function of a Stokeslet,⁶⁶ namely

and, thus, the force on the spheroidal squirmer is given by⁶⁶where and . As expected, Ψ₂ 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 α₃ = 0. Then, eqn (22) yields the swimming velocity of the squirmer (τ₀ = 1/e)

U₀ = B₁τ₀(τ₀ − (τ₀² − 1)coth⁻¹τ₀), (29)

which was already found by Keller and Wu for the case B₂ = 0.⁴⁵ As a consequence, α₂ in eqn (22) simply becomes α₂ = 2B₁c²τ₀(τ₀² − 1). Examples of fluid velocity fields of a spheroidal squirmer are presented in Fig. 1 and 3.

Fig. 3 Fluid velocity fields of a spheroidal squirmer in the laboratory frame for (a) B₁ = 1, B₂ = 0, and (b) B₁ = 0, B₂ = 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₁ = 0, B₂ = −1 is not shown, since it follows from that of the puller with B₁ = 0, B₂ = 1 by inverting the arrows.

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/τ. 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/τ. 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_φ ≠ 0). Equating the coefficients of (1/τ)ⁿ 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

v(r) = κ^FDv^FD(r) + κ^SDv^SD(r) + κ^SQv^SQ(r) + [scr O, script letter O](r⁻⁵), (30)

wherewhich decay like r⁻², r⁻³, and r⁻⁴ for large r, respectively.⁴⁸ The multipole coefficients are

κ^FD = −α₆, (34)

with the coefficients α₅ and α₆ of eqn (24) and (25). The values of the multipole coefficients in the spherical limit (b_z → b_x ≡ R, where R is the radius) follow from the above coefficients for τ₀ → ∞, c → 0, and cτ₀ = R as κ^FD = −B₂R²/2, κ^SD = −B₁R³/3, and κ^SQ = B₂R⁴/6 as expected for a spherical squirmer.⁷¹

3 Multiparticle collision dynamics

Multiparticle collision dynamics (MPC) is a stochastic, particle-based mesoscale hydrodynamic simulation method.⁵⁴ Thereby, 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

r_i(t + h) = r_i(t) + hv_i(t), (37)

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,73 We 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 toHere, 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(α) is the rotation matrix, which describes a rotation around a randomly oriented axis by the angle α. The angle α 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-of-mass 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 collision-cell lattice is introduced at every collision step.^76,77

Since 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.⁸⁰

The following simulations are performed with the mean number of particles per collision cell 〈N_c〉 = 10, the rotation angle α = 130°, and the time step , which yields a fluid viscosity of .

4 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₀,q₁,q₂,q₃) is a rotation quaternion and can be related to the rotation matrix D, which transforms vectors from the laboratory frame to the body-fixed frame⁸¹ (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

v^b = Dv^s (39)

is the corresponding vector in the body-fixed frame. For vectors in the laboratory frame, we will frequently omit the superscript. The orientation vector of a spheroid is e = D^Te^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² + b_z²) = I_y and I_z = (2M/5)b_x². When needed, the angular velocity is calculated as Ω^s = D^T(I^b)⁻¹Dl^s.

For all simulations we choose a neutrally bouyant spheroid, i.e., M = ρ(4π/3)b_zb_x², where ρ is the fluid mass density.

4.1 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.⁸²

For 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 semi-axes 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

J_i = 2m{v_i − U − Ω × (r_i − C) − D^Tu^b_sq[D(r_i − C)]} (40)

at time t + h/2 is determined, taking into account the squirmer surface fluid velocity u_sq of eqn (11).⁸³ Thereby, a useful identity to determine s is given in eqn (8) of ref. 45, and ζ is given by The velocity of the MPC particle is updated according to v_i′ = 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′, i.e., r_i(t + h) = r_i(t + h/2) + hv_i′/2.

As a consequence of the elastic collisions, the center-of-mass velocity and rotation frequency of a spheroid are finally given by

U(t + h)′ = U(t + h) + J/M, (42)

Ω(t + h)′ = Ω(t + h) + D^T(I^b)⁻¹DL, (43)

where is total momentum transfer by the MPC fluid and is the respective angular momentum transfer.

4.2 Collision step

In a first step, ghost particles are distributed inside each spheroid.^82,84 The 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

v^g_i = U + Ω × (r_i − C) + u_sq,i + v_R,i. (44)

The Cartesian components of v_R,i are Gaussian-distributed random numbers with zero mean and variance . 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 J^g_i = m(v̄^g_i − v^g_i) and L^g_i = (r^g_i − C) × J^g_i, 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

U′ = U + J^g/M, (45)

Ω′ = Ω + R^T(I^b)⁻¹RL^g. (46)

4.3 Rigid body dynamics for spheroids

During the streaming step, the spheroids move according to rigid-body dynamics, governed by⁸⁵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 (α, β, γ) = (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⁸¹ 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 τ is performed as follows:

(i) Update C and q according to (cf.eqn (49) and (50))

The parameter [small lambda, Greek, tilde] is introduced to guarantee q² = 1.

(ii) Calculate forces and torques F^s(t + τ) and T^s(t + τ).

(iii) Update U and l^s according to

5 Simulations – thermal properties and flow field

5.1 Passive colloid

For the passive spheroidal colloid (B₁ = B₂ = 0), we perform equilibrium simulations and determine 〈U_α²〉 as well as 〈(Ω^b_α)²〉 for α ∈ {x, y, z}. Due to the equipartition of energy, we expect We fix the aspect ratio b_z/b_x = 2 and vary b_x in the range b_x ∈ [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 σ_r = (〈x_theo²〉 − 〈x_sim²〉)/〈x_theo²〉 is larger for Ω^b_α than for U_α. We find the largest relative error for 〈(Ω^b_z)²〉, namely σ_r = 9.5%, 5.3%, and 3.1% for b_x = 2a, 3a, and 4a. Hence, we choose the minor axis b_x ≥ 3a in the following.

In addition, we determine the orientation correlation function 〈e(t)·e(0)〉. The theory of rotational Brownian motion⁸⁶ predicts

〈e(t)·e(0)〉 = exp(−2D^⊥_Rt), (58)

where D_R = (2D^⊥_R + D^∥_R)/3, D^∥_R = k_BT/ξ^∥, D^⊥_R = k_BT/ξ^⊥, and ξ^∥ and ξ^⊥ are the parallel and perpendicular rotational friction coefficients of a prolate spheroid with respect to the major semi-axis; explicitly⁶⁹ 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.

Fig. 4 Orientation correlation functions 〈e(t)·e(0)〉 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).

5.2 Squirmer

We determine the steady state swimming velocity of a squirmer via 〈e·U〉, which should be equal to U₀ (cf.eqn (3)). Results for various eccentricities are displayed in Fig. 5. The velocity U₀ increases with increasing eccentricity e in close agreement with the theoretical prediction of eqn (3). We confirm that the force-dipole parameter β does not affect the velocity of the squirmer, as long as the Reynolds number Re is low, i.e., Re = ρU₀b_z/η ≲ 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).

Fig. 5 Mean swimming velocity as function of the eccentricity e for a spheroidal squirmer with and B₂ = 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.

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 two-dimensional flow field of the MPC fluid, averaged over the rotation angle φ, 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₁ 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 , which we attribute to a non-vanishing Reynolds number (Re ≈ 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 squirmer–squirmer 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.

Fig. 6 Fluid flow fields of a spheroidal squirmer in the laboratory frame with b_x = 3a, b_z = 6a, , and β = 3 (a–c), and with (d–f). The magnitude of the velocity field (in units of ) 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 Δv_α = |v^theo_α − v^sim_α|/[(|v^theo_α| + |v^sim_α|)/2]. Note, due to the discrete representation of the velocity field, some streamlines end abruptly.

6 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.⁸⁴ Steric interactions between two squirmers and between a squirmer and a wall are taken into account by the procedure described in Appendix C. The initial positions and orientations of the two squirmers (i = 1, 2) are

e_1/2 = (±cos(α₀), 0, sin(α₀))^T. (63)

Here, d_cm is the initial center-of-mass distance and α₀ = (π − θ₀)/2, where θ₀ is the initial angle between e₁ and e₂. The swimming mode is chosen as and the force dipole mode β ∈ {−4,0,4}. We choose d_cm such that the squirmers are well separated and vary θ₀. 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 〈d_s〉 and the mean alignment 〈e₁·e₂〉 = 〈cos θ〉 are shown in Fig. 7 for pushers, pullers, and neutral swimmers with an initial angle θ₀ = 3π/8. Due to the setup, the squirmers initially approach each other and collide at tU₀/σ ≈ 0.5. The (persistence) Péclet number Pe = U₀/(2D^⊥_Rσ) ≈ 60 is sufficiently high, such that the squirmer orientation has hardly changed before collision. When the neutral swimmers collide, they initially align parallel (cos θ ≈ 1 at tU₀/σ ≈ 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,³⁷ but at tU₀/σ ≈ 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,³⁷ swim cooperatively and reach a stable orientation with 〈cos(θ)〉 ≈ 0.77 shortly after they collided (at tU₀/σ ≈ 1). Thereby, their cooperative swimming velocity is about 0.8U₀. 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⁸⁸ 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 initial conditions θ₀ ∈ (0, π/2). Only pullers that are nearly parallel initially (θ₀ = π/8, cos θ₀ ≈ 0.92 in Fig. 9), repel each other such that they will not reach the fixed point. For Péclet numbers Pe < 60, the fixed point remains at 〈cos(θ)〉 ≈ 0.77. However, it becomes more likely for the swimmers to escape (or never reach) the fixed point.

Fig. 7 Average surface-to-surface distance d_s and orientation of squirmers, where cos(θ) = e₁·e₂, as function of time. The solid blue, dashed black, and dotted red lines correspond to pullers β = 4, neutrals β = 0 and pushers β = −4. The standard deviation of the blue line (β = 4) is indicated by the cyan shaded region.

Fig. 8 Flow field of two cooperatively swimming pullers in the laboratory frame. The magnitude of the velocity field (in units of ) 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 Time dependence of the average alignment 〈e₁·e₂〉 = 〈cos θ〉 of two pullers with β = 4, b_z/b_x = 2 in a slit of height L_y = 7a, and various initial angles θ₀ ∈ (0,π/2) (see ESI,† Movie beta 4).

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 fact 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⁸⁹ are essential for the swimmers' dynamics.

Fig. 10 Average surface-to-surface distance d_s between two pullers with β = 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).

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 pair-swimming by the authors—is unstable to perturbations that displace one swimmer out of the swimming plane.⁹¹ Since 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 θ ∈ (0, π/2) by nearest-neighbor two-body interactions, where all angles between 0 and π/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 β, 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 β = 4. With decreasing β, the stable alignment disappears, i.e., the pullers' distance increases after collision (see ESI,† Movies beta 3 and beta 4). For increasing β the fixed point remains, but the value of cos θ decreases, i.e., the squirmers form a larger angle. With increasing wall separation, the fixed-point value of cos θ 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 〈cos θ〉 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 ≥ 2, the minimal value of β 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 β ≈ 3.5.

7 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³ to 10⁴ swimmers) can be addressed by the implementation of MPC and the squirmer dynamics on current GPUs.

Appendix

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

B Gegenbauer functions

For n ≥ 2 and x ∈ [Doublestruck 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 For n = 0, 1, they are defined as

G₀(x) = −H₁(x) = 1, G₁(x) = H₀(x) = −x. (67)

For the reader's convenience, we give the formula for the Gegenbauer functions of the first kind for n = 2, 3, and x ∈ [Doublestruck R] Furthermore, the Gegenbauer functions of the second kind for n = 2, 3, and x > 1 are given by Here, we used coth⁻¹(x) = ln([x + 1]/[x − 1])/2.

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

1 = [scr A, script letter A](x) ≡ (x − C)^TA(x − C), (72)

where the orientation matrix A can be expressed as

A = (1 − ee^T)/b_x² + ee^T/b_z². (73)

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.

C.1 Interaction between spheroids. We introduce a repulsive interaction potential between spheroids to prevent their overlap. The potential is given byHere, σ₀ and ε₀ correspond to a length and energy scale, respectively. We choose ε₀ = k_BT and σ₀ = 2d_v. The directional contact distance d_R between two spheroids, with orientation matrices A₁, A₂ and center positions C₁, C₂, is an approximation to their true distance of closest approach and is defined by

d_R = R(1 − F(A₁,A₂)^−1/2) (75)

Here, R = C₂ − C₁, R = |R|, and F(A₁,A₂) is the elliptic contact function, defined as⁹⁴ Minimization with respect to x demands ∇[scr S, script letter S](x,λ) = 0, and hence,

x(λ) = {λA₁ + (1 − λ)A₂}⁻¹{λA₁C₁ + (1 − λ)A₂C₂}. (78)

The critical value λ = λ_c that maximizes [scr S, script letter S](x(λ),λ) can be found by the root finding problem

[scr A, script letter A]₁(x(λ)) − [scr A, script letter A]₂(x(λ)) = 0. (79)

We implement Brent's root finding approach.⁹⁵ The forces and torques arising from the potential (74) can be calculated analytically and are given by‡ ⁹⁴and

× F^−3/2(x_c − C) × X_c (82)

for the first spheroid, where X_c = 2λ_cA₁(x_c − C₁). The force and torque on the second spheroid follow by Newton's action-reaction law, namely

F₂ = −F₁, (83)

T₂ = −T₁ + R × F₁. (84)

We restrict ourselves to short-range repulsive interactions by setting the potential U to a constant value for , which implies that F₁ and T₁ 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 if . This inequality is checked before a numerical calculation of d_R is employed.

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 [scr A, script letter A](x) = 1. Using the method of Lagrange multipliers, we have to minimize Λ(x,λ) = h(x) + λ([scr A, script letter A](x) − 1). The necessary condition for a minimum ∂Λ/∂x = 0 yields

e_y + λ∇[scr A, script letter A](x) = e_y + 2λA(x − C) = 0, (85)

and hence,

x = C − A⁻¹e_y/(2λ). (86)

Substitution of eqn (86) into [scr A, script letter A](x) = 1 yields Finally, we obtain the point closest to the wall asHere, the minus sign has to be chosen, which can be visualized by the example of a sphere of radius R, for which A = R⁻²1. This finally yields the height We employ the Lennard-Jones potentialfor a repulsive wall, and U_w assumes a constant value for all . We can derive the force F_α = −∂U_w/∂C_α and torque T_α = −∂U_w/∂ψ_α acting on the spheroid analytically. For the force, we findand for the torquewithHere, we used the relationwhich 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

x = C + A⁻¹e_y[(A⁻¹)_yy]^−1/2. (96)

The formulas for torque and force do not change, except that we have to insert 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 two-dimensional packing fraction N_sqπb_xb_z/(L_xL_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.⁸¹ The 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.

Table 1 Memory requirement to store fluid particle, collision cell, and spheroid properties

Memory usage per cell:
Particles per cell	4 bytes
Random seed	8 bytes
Cell energy/thermostat	4 bytes
Center of mass velocity	24 bytes
Rotation matrix	24 bytes
Center of mass	24 bytes
Inertia tensor	48 bytes
Angular momentum	24 bytes
Total per cell	160 bytes
Memory usage per fluid particle:
Particle position (single precision)	12 bytes
Particle velocity (double precision)	24 bytes
Cell index	4 bytes
Temporary memory for performance optimisation	12 bytes
Total per particle	52 bytes
Memory usage per spheroid:
Parameters describing current state	296 bytes
Parameters describing state after half time step	296 bytes
Total per spheroid	592 bytes

---

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 〈N_c〉 = 10 particles per cell, the required GPU memory is 680 bytes per collision cell. Therefore, approximately 1.5 × 10⁶ 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⁶–2 × 10⁷ 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³–2.4 × 10⁴. 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 (≈10³ squirmers) for acceptable total simulation times.

Acknowledgements

We thank A. Wysocki, A. Varghese, B. U. Felderhof, and A. J. T. M. Mathijssen for helpful discussions. Support of this work by the DFG priority program SPP 1726 on “Microswimmers – from Single Particle Motion to Collective Behaviour” is gratefully acknowledged.

References

1. J. Elgeti, R. G. Winkler and G. Gompper, Rep. Prog. Phys., 2015, 78, 056601.
2. T. Vicsek and A. Zafeiris, Phys. Rep., 2012, 517, 71.
3. M. F. Copeland and D. B. Weibel, Soft Matter, 2009, 5, 1174.
4. N. C. Darnton, L. Turner, S. Rojevsky and H. C. Berg, Biophys. J., 2010, 98, 2082.
5. D. B. Kearns, Nat. Rev. Microbiol., 2010, 8, 634–644.
6. H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher, R. E. Goldstein, H. Löwen and J. M. Yeomans, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 14308.
7. J. D. Partridge and R. M. Harshey, J. Bacteriol., 2013, 195, 909.
8. J. Bialké, T. Speck and H. Löwen, Phys. Rev. Lett., 2012, 108, 168301.
9. I. Buttinoni, J. Bialké, F. Kümmel, H. Löwen, C. Bechinger and T. Speck, Phys. Rev. Lett., 2013, 110, 238301.
10. B. M. Mognetti, A. Šarić, S. Angioletti-Uberti, A. Cacciuto, C. Valeriani and D. Frenkel, Phys. Rev. Lett., 2013, 111, 245702.
11. I. Theurkauff, C. Cottin-Bizonne, J. Palacci, C. Ybert and L. Bocquet, Phys. Rev. Lett., 2012, 108, 268303.
12. Y. Fily, S. Henkes and M. C. Marchetti, Soft Matter, 2014, 10, 2132.
13. X. Yang, M. L. Manning and M. C. Marchetti, Soft Matter, 2014, 10, 6477.
14. J. Stenhammar, D. Marenduzzo, R. J. Allen and M. E. Cates, Soft Matter, 2014, 10, 1489.
15. Y. Fily, A. Baskaran and M. F. Hagan, Soft Matter, 2014, 10, 5609.
16. G. S. Redner, M. F. Hagan and A. Baskaran, Phys. Rev. Lett., 2013, 110, 055701.
17. Y. Fily and M. C. Marchetti, Phys. Rev. Lett., 2012, 108, 235702.
18. R. Großmann, L. Schimansky-Geier and P. Romanczuk, New J. Phys., 2012, 14, 073033.
19. V. Lobaskin and M. Romenskyy, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 87, 052135.
20. A. Zöttl and H. Stark, Phys. Rev. Lett., 2014, 112, 118101.
21. A. Wysocki, R. G. Winkler and G. Gompper, EPL, 2014, 105, 48004.
22. F. Peruani, A. Deutsch and M. Bär, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 74, 030904.
23. Y. Yang, V. Marceau and G. Gompper, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 82, 031904.
24. J. Elgeti and G. Gompper, EPL, 2009, 85, 38002.
25. A. P. Berke, L. Turner, H. C. Berg and E. Lauga, Phys. Rev. Lett., 2008, 101, 038102.
26. E. Lauga and T. R. Powers, Rep. Prog. Phys., 2009, 72, 096601.
27. E. Lauga, W. R. DiLuzio, G. M. Whitesides and H. A. Stone, Biophys. J., 2006, 90, 400.
28. J. Hu, A. Wysocki, R. G. Winkler and G. Gompper, Sci. Rep., 2015, 5, 9586.
29. R. Di Leonardo, D. Dell'Arciprete, L. Angelani and V. Iebba, Phys. Rev. Lett., 2011, 106, 038101.
30. L. Lemelle, J.-F. Palierne, E. Chatre, C. Vaillant and C. Place, Soft Matter, 2013, 9, 9759.
31. M. J. Lighthill, Commun. Pure Appl. Math., 1952, 5, 109–118.
32. J. R. Blake, J. Fluid Mech., 1971, 46, 199–208.
33. J. R. Howse, R. A. L. Jones, A. J. Ryan, T. Gough, R. Vafabakhsh and R. Golestanian, Phys. Rev. Lett., 2007, 99, 048102.
34. A. Erbe, M. Zientara, L. Baraban, C. Kreidler and P. Leiderer, J. Phys.: Condens. Matter, 2008, 20, 404215.
35. G. Volpe, I. Buttinoni, D. Vogt, H. J. Kümmerer and C. Bechinger, Soft Matter, 2011, 7, 8810.
36. I. Llopis and I. Pagonabarraga, J. Non-Newtonian Fluid Mech., 2010, 165, 946.
37. I. O. Götze and G. Gompper, EPL, 2010, 92, 64003.
38. T. Ishikawa and M. Hota, J. Exp. Biol., 2006, 209, 4452–4463.
39. T. Ishikawa and T. J. Pedley, J. Fluid Mech., 2007, 588, 399–435.
40. A. A. Evans, T. Ishikawa, T. Yamaguchi and E. Lauga, Phys. Fluids, 2011, 23, 111702.
41. F. Alarcón and I. Pagonabarraga, J. Mol. Liq., 2013, 185, 56–61.
42. J. J. Molina, Y. Nakayama and R. Yamamoto, Soft Matter, 2013, 9, 4923.
43. T. Ishikawa and T. J. Pedley, Phys. Rev. Lett., 2008, 100, 088103.
44. K. Ishimoto and E. A. Gaffney, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 88, 062702.
45. S. R. Keller and T. Y. Wu, J. Fluid Mech., 1977, 80, 259–278.
46. G.-J. Li and A. M. Ardekani, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 013010.
47. K. Kyoya, D. Matsunaga, Y. Imai, T. Omori and T. Ishikawa, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 92, 063027.
48. S. E. Spagnolie and E. Lauga, J. Fluid Mech., 2012, 700, 105–147.
49. A. Zöttl and H. Stark, Phys. Rev. Lett., 2012, 108, 218104.
50. I. Pagonabarraga and I. Llopis, Soft Matter, 2013, 9, 7174–7184.
51. B. Delmotte, E. E. Keaveny, F. Plouraboué and E. Climent, J. Comput. Phys., 2015, 302, 524–547.
52. A. Malevanets and R. Kapral, J. Chem. Phys., 1999, 110, 8605–8613.
53. R. Kapral, Adv. Chem. Phys., 2008, 140, 89–146.
54. G. Gompper, T. Ihle, D. M. Kroll and R. G. Winkler, in Advanced Computer Simulation Approaches for Soft Matter Sciences III, ed. P. C. Holm and P. K. Kremer, Springer, Berlin, Heidelberg, 2009, pp. 1–87.
55. E. Tüzel, T. Ihle and D. M. Kroll, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 74, 056702.
56. C.-C. Huang, G. Gompper and R. G. Winkler, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 86, 056711.
57. S. Y. Reigh, R. G. Winkler and G. Gompper, Soft Matter, 2012, 8, 4363–4372.
58. S. Y. Reigh, R. G. Winkler and G. Gompper, PLoS One, 2013, 8, e70868.
59. G. Rückner and R. Kapral, Phys. Rev. Lett., 2007, 98, 150603.
60. M. Yang and M. Ripoll, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 84, 061401.
61. D. J. Earl, C. M. Pooley, J. F. Ryder, I. Bredberg and J. M. Yeomans, J. Chem. Phys., 2007, 126, 064703.
62. J. Elgeti, U. B. Kaupp and G. Gompper, Biophys. J., 2010, 99, 1018–1026.
63. J. Elgeti and G. Gompper, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 4470.
64. M. Theers and R. G. Winkler, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 88, 023012.
65. J. Hu, M. Yang, G. Gompper and R. G. Winkler, Soft Matter, 2015, 11, 7843.
66. J. Happel and H. Brenner, Low Reynolds number hydrodynamics: with special applications to particulate media, Springer Science & Business Media, 1983.
67. G. Dassios and P. Vafeas, Phys. Res. Int., 2008, 2008, e135289.
68. A. M. Leshansky, O. Kenneth, O. Gat and J. E. Avron, New J. Phys., 2007, 9, 145.
69. S. Kim and S. J. Karrila, Microhydrodynamics: Principles and Selected Applications, Butterworth-Heinemann, 2013.
70. A. J. T. M. Mathijssen, D. O. Pushkin and J. M. Yeomans, J. Fluid Mech., 2015, 773, 498–519.
71. O. S. Pak and E. Lauga, J. Eng. Math., 2014, 88, 1–28.
72. E. Allahyarov and G. Gompper, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 66, 036702.
73. H. Noguchi, N. Kikuchi and G. Gompper, EPL, 2007, 78, 10005.
74. H. Noguchi and G. Gompper, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 016706.
75. M. Theers and R. G. Winkler, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91, 033309.
76. T. Ihle and D. M. Kroll, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2001, 63, 020201.
77. T. Ihle and D. M. Kroll, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 67, 066705.
78. C. C. Huang, A. Chatterji, G. Sutmann, G. Gompper and R. G. Winkler, J. Comput. Phys., 2010, 229, 168–177.
79. C.-C. Huang, A. Varghese, G. Gompper and R. G. Winkler, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91, 013310.
80. E. Westphal, S. P. Singh, C. C. Huang, G. Gompper and R. G. Winkler, Comput. Phys. Commun., 2014, 185, 495–503.
81. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1987.
82. J. T. Padding, A. Wysocki, H. Löwen and A. A. Louis, J. Phys.: Condens. Matter, 2005, 17, S3393.
83. M. T. Downton and H. Stark, J. Phys.: Condens. Matter, 2009, 21, 204101.
84. A. Lamura, G. Gompper, T. Ihle and D. M. Kroll, EPL, 2001, 56, 319.
85. I. P. Omelyan, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1998, 58, 1169–1172.
86. L. D. Favro, Phys. Rev., 1960, 119, 53–62.
87. A. J. T. M. Mathijssen, A. Doostmohammadi, J. M. Yeomans and T. N. Shendruk, 2015, arXiv:1511.01859 [cond-mat, physics:physics].
88. K. Drescher, J. Dunkel, L. H. Cisneros, S. Ganguly and R. E. Goldstein, Proc. Natl. Acad. Sci. U. S. A., 2011, 10940, 108.
89. J. d. Graaf, A. J. T. M. Mathijssen, M. Fabritius, H. Menke, C. Holm and T. N. Shendruk, Soft Matter, 2016, 12, 4704–4708.
90. I. O. Götze and G. Gompper, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 82, 041921.
91. T. Ishikawa, M. P. Simmonds and T. J. Pedley, J. Fluid Mech., 2006, 568, 119–160.
92. B. U. Felderhof, 2016, arXiv:1603.08574 [physics].
93. Z. X. Wang and D. R. Guo, Special Functions, World Scientific, 1989.
94. L. Paramonov and S. N. Yaliraki, J. Chem. Phys., 2005, 123, 194111.
95. W. H. Press, Numerical Recipes 3rd Edition: The Art of Scientific Computing, Cambridge University Press, 2007.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c6sm01424k

‡ Note that eqn (54) of ref. 94 contains a typographical error. The factor 24 needs to be replaced by 12.

---