Maximum in density heterogeneities of active swimmers

Abstract

Suspensions of unicellular microswimmers such as flagellated bacteria or motile algae can exhibit spontaneous density heterogeneities at large enough concentrations. We introduce a novel model for biological microswimmers that creates the flow field of the corresponding microswimmers, and takes into account the shape anisotropy of the swimmer's body and stroke-averaged flagella. By employing multiparticle collision dynamics, we directly couple the swimmer's dynamics to the fluid's. We characterize the nonequilibrium phase diagram, as the filling fraction and Péclet number are varied, and find density heterogeneities in the distribution of both pullers and pushers, due to hydrodynamic instabilities. We find a maximum degree of clustering at intermediate filling fractions and at large Péclet numbers resulting from a competition of hydrodynamic and steric interactions between the swimmers. We develop an analytical theory that supports these results. This maximum might represent an optimum for the microorganisms' colonization of their environment.

---

I. Introduction

Physical interactions in suspensions of microswimmers consisting of bacteria or algae have been recognized to play an important role in the swimmers' collective behavior.^1–3 The nonequilibrium character of active suspensions, where the energy injection takes place at the scale of the microorganisms, produces myriad mesmerizing phenomena, such as complex interaction with solid surfaces,^4–6 the spontaneous formation of spiral vortices,⁷ directed motion,⁸ swarming,⁹ bacterial turbulence,¹⁰ and self-concentration.¹¹

Almost invariably, motile microorganisms move in an aqueous environment, where, because of their size, viscous forces dominate, and inertial forces are completely negligible. In fact, consideration of the Navier–Stokes equations identifies that the nature of the dynamics is dictated by the ratio of viscous to inertial forces, known as the Reynolds number [scr R, script letter R] = σvρ/η, where σ is the typical size of the microorganism, v its mean velocity, and ρ, η are the fluid's density and viscosity, respectively. For Escherichia coli, e.g., σ ≈ 10 μm, v ≈ 30 μm s⁻¹, and for water ρ ≈ 10³ kg m⁻³, η ≈ 10⁻³ Pa s, which result in [scr R, script letter R] ≈ 10⁻⁵. As noted by Purcell,¹² this means that if the propulsion of a swimmer were to suddenly disappear, it would only coast for 0.1 Å. Thus, the state of motion is only determined by the forces acting at that very moment, and inertia is negligible.

Due to the microswimmers' low Reynolds numbers, the sum of viscous drag and thrust balances out to zero, in most situations. A direct consequence of force-free motion is that the leading term of the solution of the Stokes equation for a microswimmer is a symmetric force dipole (or stresslet).

Biological microswimmers are complex systems because of the combination of biological, biochemical and physical processes all taking place at the same time. It is thus of great scientific value to develop theoretical models that isolate the relevant degrees of freedom and interactions. Considerable work has been done in recent years, and various models have been introduced, like the squirmer model,^13–31 the shape anisotropic raspberry swimmer,^32–34 the force–counterforce model,^35–38 the catalytic dimers,³⁹ or other hydrodynamic models.^40–43 Experiments have confirmed that the flow field of flagellated bacteria like E. coli is to very good approximation modeled by a simple force dipole,⁴⁴ whereas Chlamydomonas reinhardtii are modeled by three Stokeslets.⁴⁵ Furthermore, as cell shapes vary greatly in the natural world, and realistic steric interactions are important in dense suspensions, a model that allows for flexibility in the shape of a microswimmer is a highly desirable feature. In this article we fill this lacuna. We derive a model for a flexible-shape microswimmer that produces self-propulsion by means of a force dipole for pusher-like microswimmers, or three Stokeslets for puller-like microswimmers.

An efficient method to simulate fluids at mesoscopic scales, and their hydrodynamics is the multiparticle collision dynamics (MPCD) technique.⁴⁶ MPCD is a particle-based simulation method that correctly produces hydrodynamic modes. Due to its particle nature MPCD naturally includes thermal fluctuations, and can be easily coupled to molecular dynamics methods of solutes, colloids,⁴⁷ and active swimmers.^20,24,26 The MPCD technique in fact proves to be ideal for our purposes.

The nonequilibrium phase diagram of microswimmers has been subject to considerable interest, especially with regard to the emergence of density heterogeneities in the swimmers' distribution.² We explore the phase diagram of active swimmers and describe the presence of heterogeneities in the spatial distribution of both pushers and pullers. These heterogeneities arise due to the hydrodynamic interactions between the swimmers and relate to existing hydrodynamic theories.^2,48–50 Interestingly, we find a maximum in the heterogeneities as filling fraction and Péclet number are varied. By using both computer simulations and analytical theory, we demonstrate that this maximum results from a competition between hydrodynamic and steric interactions, where the latter temper the hydrodynamic instability at higher filling fractions. This optimum might have important biological implications on the ability of motile bacteria and algae to form colonies or biofilms.

The remainder of this article is organized as follows. In Section II we introduce the model for the microswimmer. Section III describes the physical properties of the fluid and the microswimmer's flow field. In Section IV we present the nonequilibrium phase diagram of our model microswimmers, and specifically we characterize the density heterogeneities emerging from their hydrodynamic interactions and show that these are tempered by steric interactions at higher filling fractions. In Section V we present an analytical theory and show that we also find a maximum in heterogeneity, which is mediated by the interplay of hydrodynamic and steric interactions. Finally, in Section VI we discuss our main results and summarize our conclusions.

II. Model

We employ a stroke-averaged model of biological microswimmers, similarly to,^1,10,51 taking into account the asymmetric shape of biological microswimmers due to the cell's body and the flagella. The swimmer is thus modeled as an asymmetric dumbbell, as depicted in Fig. 1, that mimics a C. reinhardtii or an E. coli cell. The smaller sphere models the swimmer's body and the larger sphere is a stroke average of the region spanned by flagellar motion. The rigid body dynamics of the dumbbells is simulated using Newton's equations and quaternion dynamics (for details see Appendix A1 and A2).

Fig. 1 Schematic representation in a perspective view of the active swimmer model for (a) a puller-type, and (b) a pusher-type microswimmer. In both panels, the small green spheres (1) represent the swimmer's body, and the larger transparent spheres (2) represent the stroke-averaged space spanned by the flagella. The red spheres (3) with embedded arrows represent the regions where the forces are applied. The golden arrows (4) represent the swimming direction. The lines with arrows (5) are a sketch of the hydrodynamic streamlines generated by the swimmers.

In addition to the rigid body dynamics, we simulate the fluid surrounding the swimmers using MPCD (see Appendix A3), which is a mesoscopic, particle based method, that reproduces hydrodynamics at the Navier–Stokes level.⁴⁷ Precise measurements⁴⁴ show that the flow field of pusher-type microswimmers is well modeled by a force dipole. The pullers flow field, on the other hand, is well represented by a three-Stokeslet solution of the Stokes equation.⁴⁵ These two flow-field models are included in our simulations by adding force regions to the fluid, as depicted by the red regions (3) in Fig. 1. Furthermore, we couple the fluid and swimmers' dynamics by imposing a no-slip boundary condition on the swimmers' surface (for details see Appendix A4).

It should be noted that the dumbbell model introduced here has an anisotropic rigid shape that can be easily modified to more complex shapes. Squirmer models have so far been described for spherical or ellipsoidal shapes.²⁵ Our hydrodynamic flow field is very similar to the three-bead-spring model.⁵² However, the rigid body dynamics differ in the fact that in the three-bead-spring model the beads are connected by springs, whereas ours is a rigid model.

In the following, we express all physical quantities in terms of the MPCD particle mass m, the size of an MPCD grid cell a, and the temperature T of the fluid. We simulate (see Section IV) N = 300–1560 active swimmers in a cubic domain with side length of 100a, which is approximately 20 times the size of an individual swimmer, and with periodic boundary conditions. The resulting filling fraction ranges between ϕ = 0.05 and ϕ = 0.29. The average number of MPCD particles per cell is 〈N_C〉 = 20 such that the total number of MPCD particles in a simulation is 2 × 10⁷. The Reynolds number considered here ranges from [scr R, script letter R] ∼ 0.01 to [scr R, script letter R] ∼ 0.1 and the Péclet number reaches from 2.2 × 10² to 2.6 × 10³.

The full details of our numerical implementation are given in Appendix A.

III. Characterization of the fluid and active hydrodynamics

In the following we describe calculations aimed at characterizing the thermal (equilibrium) properties of our model in the passive case, and also the flow field generated by the active motion.

We first consider a passive colloid (with the same geometry described above, and in Fig. 1) immersed in the MPCD fluid, that is, we carried out equilibrium simulations without activity. The equipartition theorem applied to the passive colloid for the translational motion predicts that the average of each velocity component squared is 〈U_α²〉 = k_BT/M, where M is the colloid's mass, k_B is the Boltzmann constant, and α = x, y, z. For our system, we find a theoretical value of the translational motion 〈U_theory²〉 = 2.7 × 10⁻⁴k_BT/m and the simulations give 〈U_x²〉 = 〈U_y²〉 = 〈U_z²〉 = 2.3 × 10⁻⁴k_BT/m. In case of the rotational motion, the equipartition theorem predicts for the angular velocity 〈(Ω^b_α)²〉 = k_BT/I_mα, where I_mα is the moment of inertia tensor. Considering our swimmer, whose long axis is aligned with the z axis of a Cartesian reference frame, the theoretical prediction for the angular motion in x and y direction yields 〈(Ω^b_x,y)²〉 = 1.3 × 10⁻³k_BT/ma² while the simulations yield 〈(Ω^b_x)²〉 = 〈(Ω^b_y)²〉 = 1.1 × 10⁻³k_BT/ma². In the z direction, the theory predicts 〈(Ω^b_z)²〉 = 1.4 × 10⁻³k_BT/ma² and the simulations give 〈(Ω^b_z)²〉 = 1.2 × 10⁻³k_BT/ma².

We now consider the active motion at [scr R, script letter R] ≪ 1. Hydrodynamics at low Reynolds numbers (relevant for micron-sized objects) allow a great simplification of the Navier–Stokes equations: the nonlinear, inertial effects can be neglected, and the governing equations are the Stokes equations

η∇²u = ∇p − f^ext, ∇·u = 0, (1)

where u(r) is the fluid velocity, p(r) the pressure, f^ext(r) is a body force acting on the fluid at position r, and η is the viscosity of the fluid. Solving the Stokes equations means obtaining expressions for u and p that satisfy eqn (1) and the boundary conditions. From this knowledge, the stress tensor σ can be calculated. For a Newtonian fluid, σ depends linearly on the instantaneous values of the velocity gradient, so that one can write σ = −pI + η[∇ ⊗ u + (∇ ⊗ u)^T], where ⊗ indicates the tensor product, and I the identity tensor. Because the Stokes equations are linear, their solution can be formally written in terms of the convolution of a Green's function with the inhomogeneous term f^ext 53,54

In free, three-dimensional space, the Green's function is found by considering a point force f^ext = feδ(r) acting on an infinite, quiescent fluid, where e is the unit vector representing the direction of the force, and δ(r) is the Dirac distribution. A straightforward calculation⁵⁵ gives the Oseen tensor where r̂ ≡ r/r, r = |r|, and the resulting flow field , which is termed a ‘Stokeslet’ and decays with distance as r⁻¹.

A theoretical prediction for the puller flow field is constructed from three Stokeslets, and for the pusher we use two Stokeslets. The Stokeslets positions are placed at the midpoints of the respective force regions in the swimmer model (see Fig. 1).

We now consider the flow field generated by the active motion of our model microswimmer. We switch on the active motion with a force f₀ = 50k_BT/a and carry out the full dynamics as described in Section II (see also Appendix A for a full definition of f₀). Fig. 2 shows the flow fields of a pusher and a puller in the lab frame. As expected, the flow field of the puller is contractile, as fluid is drawn in from the front and the back, while fluid is pushed away normal to the swimming direction [Fig. 2(a)]. The situation is reversed in the pusher case [Fig. 2(c)] where fluid is pushed out at the front and back of the swimmer. The theoretical predictions for both the puller [Fig. 2(b)] and the pusher [Fig. 2(d)] show a good quantitative agreement with the simulated flow fields. For the Stokesian algebraic decay of u(r), the reader is referred to Section I in the ESI.†

Fig. 2 Time-averaged flow field generated by (a) our model puller, (b) theoretical puller, (c) our model pusher, and (d) theoretical pusher. We show cross-sections on the x–y plane at z = 0. The force strength is f₀ = 50k_BT/a. The large central white regions show the hard cores of the active swimmers. The thin lines with arrows mark the streamlines, while the color code shows the magnitude of the flow velocity normalized to the thermal velocity. The large black arrows indicate the direction of motion.

The effective velocity v_eff ≡ |〈e·U〉| of an isolated swimmer in the steady state depends linearly on the active force f₀.¹ Here, e is the swimmer's orientation and U is the swimmer's velocity. From our simulations, we calculate v_eff for a pusher [see Fig. 3]. The linear fit has a slope of . The analogous results for pullers are also shown in Fig. 3, where the slope of the linear fit is . The linear dependence of v_eff on f₀ shows that our simulations correspond to the Stokes flow regime.

Fig. 3 Dependence of the pusher and puller velocity v_eff on the active force f₀. Lines are linear fits to the simulated data. The observed linear dependence is evidence of the Stokes flow regime.

Further simulations on the two particle interaction statistics and two particle flow fields can be found in Section I (ESI†).

IV. Density heterogeneities

In the following we will study the nonequilibrium phase behavior of our active swimmer model. The filling fraction we employ here is computed using the volume of the swimmer's body V_B, as well as the volume that is spanned by the flagella sphere V_F, while taking into account their overlap volume V_Olwhere V is the volume of the simulated system. Note, that in an experiment only the body of the cell would be taken into account, thus the filling fraction should then be rescaled by V_B/(V_B + V_F − V_Ol) = 7.5 × 10⁻². Furthermore, we vary the strength of the active force f₀, which changes the propulsion speed v_eff as well as the strength of the hydrodynamic interactions between the swimmers. The Péclet number captures the ratio of advection to diffusion, and can be computed usingwhere we used the linear relation (fitted slope) between active velocity and force dipole strength from Section III. Furthermore, σ = 5a is the typical length of the swimmer and for the diffusion constant we assume .

We analyze the system's density using a Voronoi tessellation⁵⁶ and compute the local volume for each swimmer. A global measure for the heterogeneity of a configuration of swimmers is given by the standard deviation of the distribution of local Voronoi volumes σ_loc. In order to remove trivial scaling factors, we compare σ_loc to the standard deviation σ_rnd of local Voronoi volumes for random homogeneous configurations of nonoverlapping, passive dumbbells (as in Fig. 1) with the same filling fraction. Fig. 4 shows the resulting phase diagram, with the dependence of σ_loc/σ_rnd on Péclet number and filling fraction. Here, positive values of [scr P, script letter P] correspond to pusher-type swimmers, whereas negative values are puller-type swimmers.

Fig. 4 Standard deviation of local Voronoi volume σ_loc compared to standard deviation σ_rnd of a homogeneous configuration. The Péclet number [scr P, script letter P] as well as the global filling fraction ϕ are varied. Positive Péclet numbers correspond to pusher-type and negative to puller-type swimmers.

The phase diagram shows that for both pullers and pushers, initially σ_loc/σ_rnd grows with [scr P, script letter P] and filling fraction, then it reaches a maximum and drops to lower values. The initial increase is related to an instability that is mediated by the hydrodynamic interactions of the microswimmers which has also been found in ref. 1, 2 and 48–50. Intuitively this can be understood in the following way: the hydrodynamic flow field from the swimmers creates heterogeneities in the fluid's velocity, which couple back to the swimmers and produce heterogeneities in the density. It should be noted that the Péclet number of the system has to be rather large, such that the hydrodynamic interaction between the swimmers is strong enough to produce heterogeneous structures. As the filling fraction increases, steric interactions grow in importance and compete with the hydrodynamic instability. Thus, we ascribe the presence of the maximum in σ_loc/σ_rnd to a tapering effect of the steric interactions on the hydrodynamic instability. This tapering effect becomes visible only when steric interactions are fully accounted for.

To test the hypothesis that steric interactions stabilize the hydrodynamic instability, we carry out two more types of simulations: first, we exclude the steric effects by setting the hard-core repulsion [eqn (A6)] to zero, thus only hydrodynamic effects are included. Second, we carry out Brownian dynamics simulations, which completely neglect hydrodynamic interactions. More details about the Brownian dynamics simulations are given in Section III of the ESI.†

In Fig. 5 we show our original simulations that fully account for hydrodynamic and steric interactions, the simulations without steric interactions, and the active Brownian simulations. The simulations with hydrodynamics alone give rise to a strong increase of σ_loc/σ_rnd, but no maximum occurs, while the active Brownian simulations show a monotonic increase of σ_loc/σ_rnd, which is much less pronounced. Our original simulations (with both hydrodynamics and steric effects) exhibit intermediate values of σ_loc/σ_rnd. Thus, we conclude that the maximum which we see is mediated by an interplay of the hydrodynamic interactions and the steric interactions, confirming our hypothesis.

Fig. 5 Dependence on the filling fraction of the standard deviation of local Voronoi volumes σ_loc compared to the standard deviation σ_rnd of a homogeneous configuration for pushers (a), and pullers (b). The Péclet number is fixed to [scr P, script letter P] = 2.6 × 10³ for pushers and [scr P, script letter P] = 2.4 × 10³ for pullers. The circles are the results for simulations with both hydrodynamic and steric interactions; the squares are the results for simulations including only hydrodynamic interactions; and the triangles are simulations including only steric interactions.

V. Theoretical analysis

To bolster our numerical results, we develop an analytical theory of microswimmers that explicitly includes hydrodynamic and steric interactions. As in our numerical model (see Section II), we consider the dynamics of an asymmetric dumbbell, which is described by the following effective Langevin equationswhere r_Li is the position of the front sphere of the swimmer i with radius a_L and r_Si the position of the respective back sphere with radius a_S. The front and back spheres of each swimmer are connected by an infinitely thin rigid rod of length l. The swimmer is coupled to the fluid velocity u, which is determined by the Stokes equation including a stochastic and an active force term

η∇²u = ∇p − f^active + f^noise. (7)

Here, the active force is given by a force dipolewhich points along the orientation of the swimmer e_i and has a force strength f. The orientation is defined as e_i = (r_Li − r_Si)/l, i.e., the unit vector connecting the back to the front sphere. Furthermore, fluctuations in the swimmers motion are added to the fluid viawhere ξ^L,S_i(t) are noise terms with 〈ξ^L,S_i(t)ξ^L,S_j(t′)〉 = 2Γ_L,SIk_BTδ_ijδ(t − t′), and Γ_L,S = 6πηa_L,S are the friction coefficients of the front and back sphere. Considering only the active term eqn (7) can be solvedwith the Oseen tensor , which is regularized by using O(|r| ≤ a_L,S) ≡ I/Γ_L,S. Given the solution (10), the flow velocity can be eliminated from eqn (5) and (6). Using the hydrodynamic center , we can simplify eqn (5) and (6) towhere v₀ is the propulsion speed and ζ_hy = 3πη(a_S + a_l) is the hydrodynamic friction coefficient. The leading order of the multipole expansion of the hydrodynamic force F_ij and torque τ_ij between swimmers i and j are given in Section IV (ESI†). The random forces ξ_i and ξ^R_i are Gaussian white noises with

〈ξ_i,α(t)ξ_j,β(t′)〉 = 2Dδ_ijδ_αβδ(t − t′), (13)

〈ξ^R_i,α(t)ξ^R_j,β(t′)〉 = 2D_Rδ_ijδ_αβδ(t − t′), (14)

where D is the translational and D_R the rotational diffusion coefficient.

Following,¹ we derive the one-body Smoluchowski equation from eqn (11) and (12). As we have seen in Section IV steric interactions among swimmers may play a crucial role in their dynamics. Thus, we explicitly include the effect of steric interactions by means of the Ansatz^57,58

v(c) = v₀ − cζ. (15)

Here, c is the concentration and the constant ζ quantifies how much the swimmers are slowed down by the steric interactions (for more details see ref. 57 and 58). The Smoluchowski equation then readswhere p(r,e,t) is the one-particle probability distribution function of finding a swimmer at position r, with orientation e at time t. The first term on the right hand side of eqn (16) takes into account the active motion with a density dependent velocity, due to steric interactions; the second term accounts for the hydrodynamic forces and the third term for the hydrodynamic torques; the last two terms are responsible for translational and rotational diffusivity respectively. Here, F_hy are the hydrodynamic two-body forces, and τ_hy is the hydrodynamic torque between two particles (see also ESI,† Section V). To make progress with this equation, we follow the standard path: we consider a multipole expansion and compute moment equations for the concentration c, the polarization P and the nematic order tensor QThe full equations for c, P and Q are given in Section V (ESI†). We linearize these moment equations around the isotropic state, described by c = c₀ + δc, P = δP and Q = δQ and turn to Fourier space, with wave vector k, where the fields are denoted by δc̃(k), δP̃(k), and δQ̃(k). To first order in the fluctuations the equations governing the temporal evolution readwhere k = |k| is the absolute value of the wavevector,[A_ij]^ST is the symmetric traceless form of the tensor A_ij, and we now denote i, j = x, y, z. The terms involving the force dipole strength f stem from the hydrodynamic interactions and can cause instabilities in the system. To analyze the stability of the system we will first consider pullers (f < 0) and in the second step pushers (f > 0), where we systematically keep terms up to order k².

For pullers (f < 0) the concentration fluctuations δc̃ are dominant. To analyze the fluctuations in the concentration we use a large length-scale and long time-scale (D_Rt ≫ 1) approximation for the longitudinal polarization δP̃_‖ = k̂_iδP̃_i, k̂_i ≡ k_i/k, which reads

whereas the fluctuations δQ̃_ij are of higher order in k, when terms of order k² in eqn (20) are kept. Inserting the quasi-stationary solution eqn (24) into the eqn (20) yieldsFor pullers the term introduces an instability at low concentrations c₀, which are counteracted by the term , that stabilizes the system at higher concentrations. Since the first term is dominated by the hydrodynamic interactions and the second term by the steric interactions, we can draw the same conclusion as from the simulations: the hydrodynamic interactions cause heterogeneities in the system which are suppressed by the steric effects at larger c₀. Moreover, inspection of the fastest-growing eigenvalue reveals a maximum of instability as a function of c₀ and, hence, a maximum degree of heterogeneity. The position of the maximum of the filling fraction can be estimated from microscopic information extracted from scattering events between two swimmers; we find ϕ^max = 0.12, which is consistent with the simulation in Fig. 5(b) (for details see Appendix B).

We also find an instability in the splay fluctuations of the nematic tensor δQ̃_‖‖ = k̂·Q̃·k̂. Here, the approximation of the polarization fluctuations [eqn (21)] is different, since we need to consider eqn (22) for the counting of powers of k. In the large length and time scale limit we arrive at

Inserting eqn (26) into eqn (22) and projecting on the splay part yieldsCombing eqn (27) with eqn (22) givesFrom this term we also find a maximum instability, but the resulting instability is sub-dominant compared to the instability in the concentration fluctuations.

For pushers (f > 0) the bend component of the nematic tensor fluctuations δQ̃_‖⊥ = k̂·δQ̃·(I − k̂ ⊗ k̂) become unstable, since here k_i[scr M, script letter M]_ij(δ_jl − k_jk_l) < 0, ∀l = x, y, z. For the polarization we use the large length and time scale limit in eqn (26). Inserting eqn (26) into the term [k_iδP̃_j]^ST of eqn (22), and projecting onto the bend part gives

Therefore, we have a single equation for the bend nematic fluctuations δQ̃_‖⊥, which is given byHere, the term destabilizes the system at low concentrations c₀ through nematic fluctuations, which are counteracted by the term , that stabilizes the system for higher concentrations. Again, the first term, which destabilizes the system, is dominated by the hydrodynamic interactions, whereas the second, stabilizing term comes from the steric interactions. Additionally from microscopic information extracted from scattering events between two swimmers we can estimate the filling fraction of the maximum heterogeneities, with the result ϕ^max = 0.23, which is in accordance with the simulation from Fig. 5(a) (for details see Appendix B).

VI. Conclusions

We have presented a new model for biological microswimmers that is based on Stokeslets and the stroke averaged motion of their flagella. The Stokeslets were distributed to model the flow fields of C. reinhardtii or E. coli cells. Furthermore, our model takes into account the anisotropic shape of a microswimmer. Typical for this is the shape of a C. reinhardtii cell, which is well modeled by an asymmetric dumbbell.¹⁰ Self-propulsion is generated through a symmetry breaking due to the asymmetric shape and force-free motion. The fluid and the hydrodynamic interactions are explicitly included with MPCD.

We show that the flow fields produced in our simulations can be predicted using simple formulae from the literature. These formulae also correspond to the experimentally measured flow fields.^44,45 Additionally, we test the effective velocity of the microswimmer model, and find that it depends linearly on the applied force, in agreement with the Stokes flow regime.

We study the phase diagram in terms of filling fraction and Péclet number. We find that both pullers and pushers exhibit density heterogeneities. The density heterogeneities show a maximum at intermediate filling fractions and high Péclet number. To determine the mechanism underpinning this phenomenon, we perform additional simulations, in which either the steric interactions or the hydrodynamic interactions were switched off. Simulations with active Brownian particles showed a small linear increase in the density heterogeneities, while simulations without the steric interactions show strong density heterogeneities. This is an instability caused only by the hydrodynamic interactions, which is known from the literature.^1,2,48–50 However, no maximum arises in the simulations with only hydrodynamic or only steric interactions, which shows that the maximum in the density heterogeneities is mediated by an interplay of the hydrodynamic and the steric interactions. The hydrodynamic interactions destabilize the system, whereas the steric interactions stabilize the system as the filling fraction grows and thus a maximum in density heterogeneity arises.

We have also developed an analytical theory based on a Smoluchowski equation which includes steric as well as hydrodynamic interactions. We computed the hydrodynamic moments of this equation and performed a linear stability analysis of the moment equations around the homogeneous state. For both puller and pusher-type swimmers, we found that at low concentration the system is destabilized by hydrodynamic interactions. At higher concentrations, however, the instabilities are counteracted by the steric interaction. This interplay gives rise to a maximum in the instability of the homogeneous state, and thus a maximum heterogeneity in the concentration of swimmers. The position of the maxima calculated from the analytics is in accordance with the simulations. Our continuum theory does not explicitly account for steric effects at the microscopic level, which induce short range correlations, and lubrication forces, and influence the short range hydrodynamics. These effects, however, are captured by the numerical simulations. The agreement between our theory and simulations gives us confidence that our assumptions effectively includes the dominant physical effects.

The physical pictures from both simulations and analytical theory fit together: both show that the homogeneous state is not stable and there is a maximum of instability. Also, both analyses show that the instability arises from hydrodynamic interactions and is suppressed by the steric interactions.

The maximum in the density heterogeneities might have important biological implications: it points to a possible, optimal filling fraction and Péclet number for the formation of heterogeneous distributions of motile microorganisms. Bacteria or microalgae exhibiting these optimal parameters are more likely to form colonies or biofilms.

Conflicts of interest

There are no conflicts of interest to declare.

Appendix A: model details

The microswimmer is characterized by its mass M, center of mass position R and orientation q. In the following we describe the equations governing the motion of the microswimmer, the fluid dynamics implemented through the MPCD, and their coupling.

1. Rigid body dynamics of swimmers

As we do not consider shape deformable swimmers, we are only concerned with rigid-body dynamics. The most general motion of a rigid body is the combination of a translation along an axis (the Mozzi axis) and a rotation around the same axis, as per the Mozzi–Chasles theorem.⁵⁹ Any orientation in space can be described using three numbers, that are commonly represented with the Euler angles, which correspond to three elementary rotations. However, there are a number of issues with the choice of the Euler angles. For instance, the composition of rotations with Euler angles or rotation matrices is rather complex, and involves trigonometric functions which lead to an accumulation of rounding-off errors. Eventually the matrices representing the rotations may become not orthogonal. Importantly, for some values of the Euler angles there are discontinuous jumps in the representation. More fundamentally, the Euler angles do not generate a covering map of the rotation group SO(3), that is, the map from Euler angles to SO(3) is not always a local homeomorphism. Fortunately, the topology of SO(3) is diffeomorphic to the real projective space [Doublestruck P]₃([Doublestruck R]) which admits a universal cover represented by the group of unit quaternions q = (q₀,q₁,q₂,q₃)^T, where the superscript T indicates the matrix transposition.

The equations of motion for the rigid body dynamics in three dimensions and in terms of quaternions read⁶⁰

where Ω^b is the angular velocity of the swimmer, I^b_m the moment of inertia tensor of the swimmer in the body frame, and the indices (α,β,γ) take on as values the cyclic permutations of (x,y,z). In eqn (A1), F = −∇Φ and T = R_F × F are the force and torque, respectively, acting on the swimmer due to steric interactions with the neighbor, where R_F is the vector connecting the center of mass of the swimmer to the point of contact with the neighbor, and the matrix W is (see also ref. 61)The repulsive, steric interactions among swimmers are modeled using a Weeks–Chandler–Andersen potential⁶²if r_ij,ab < 2^1/6σ_ab, and Φ(r_ij,ab) = 0 otherwise, where r_ij,ab ≡ |r_ia − r_jb| is the distance between sphere a of swimmer i and sphere b of swimmer j, ε is the energy scale and σ_ab is the sum of the radii of sphere a and sphere b. For the numerical integration we use the Verlet algorithm proposed in ref. 60, which was also used and discussed in detail in ref. 25.

Given a vector in the laboratory frame f the transformation to the body frame vector f^b is given by

f^b = Df, (A7)

where the matrix D(q) is constructed from the quaternions.

Quaternions are represented as q = q₀ + q₁i + q₂j + q₃k, with q₀,…,q₃ ∈ [Doublestruck R], and i² = j² = k² = ijk = −1.

The unitary matrix D that transforms vectors from the lab to the body frame is (see also ref. 61)

Thus, the orientation of the swimmer at any given time is found from D⁻¹(q(t))(0,0,1)^T.

Note that all quantities that do not carry an index b are calculated in the laboratory frame.

2. Center of mass and moment of inertia

Fig. 6 shows the detailed geometry and dimensions of our model microswimmers (both pullers and pushers). In the body frame, the swimmer is aligned with the z direction and the coordinates of the centers of the B and F spheres are z₁ and z₂, respectively. Given a homogeneous mass distribution the center of mass of the swimmer is given bywhere V_i are the volumes of the spheres and V_Sci are the volumes of their spherical caps, which are cut by the other sphere⁶³and h_i are the heights of the spherical capsThe moment of inertia for a spherical cap about the x- as well as y-direction isand about the z-direction isBy using the moment of inertia of a sphere and with the use of the parallel axis theorem, we compute the moments of inertia of the swimmer as

I_(x,y) = I_Sp1,(x,y) − I_Sc1,(x,y) + ρ(V₁ − V_Sc1)x₁², (A14)

I_z = I_Sp1,z + I_Sp2,z − I_Sc1,z − I_Sc2,z. (A15)

Fig. 6 Details on the geometry of our model puller-type swimmer (a) and pusher-type swimmer (b). Regions B (green) and F (empty circle) are the body with diameter d_b, and stroke-averaged flagella with diameter d_f, respectively, and they are separated by a distance l. The red arrow denotes the swimmers orientation and C is the center of mass. Black circles are the force poles acting on the fluid. For pullers (a) the region R_pl,1 has the diameter d_pl,1 and the regions R_pl,2, R_pl,2′ have the diameter d_pl,2. For pushers (b) the regions R_ps,1 and R_ps,2 have the diameter d_ps.

3. Multiparticle collision dynamics

To simulate a fluid at fixed density ρ and temperature T surrounding the swimmers, we use the MPCD algorithm, which is a mesoscopic, particle based method⁴⁶ to simulate a fluid at the Navier–Stokes level of description. We include the Andersen thermostat and the conservation of angular momentum into the MPCD dynamics; the resulting algorithm is usually denoted as MPC-AT+a.^47,64,65 The fluid is modeled using N_fl point-like particles of mass m, whose dynamics are executed through two steps: the streaming step and the collision step. In the streaming step the fluid particles' positions r_i, i ∈ [1,N_fl] are updated according to

r_i(t + δt) = r_i(t) + v_i(t)δt, (A16)

where v_i(t) is their velocity and δt is the MPCD timestep.

The collision step mediates the interactions between the particles. Here, the system is divided into N_c collision cells with a regular grid of lattice constant a. The center of mass velocity in each cell is calculated and remains constant during the collision step, whereas the fluctuating part of the velocity of every fluid particle i is randomized, which mimics the collision between particles. Hence, the velocity of particle i is updated as follows⁶⁴

where the random velocity v^ran_i has components sampled from a Gaussian distribution with zero mean and variance , is the number of fluid and ghost particles (see Section A4) in cell . The vector r_j,c is the position of the neighboring particle j relative to the center of mass of the cell . In eqn (A17), Π⁻¹ is the inverse of the moment of inertia tensor for the fluid particles in cell . Note that Π⁻¹ is a dynamical quantity that has to be updated at every timestep, and it also includes the ghost particles within the swimmer (see Section A4).

To ensure Galilean invariance and avoid the build-up of spurious correlations in the velocities,⁶⁶ the usual grid shift is performed at each timestep, that is, the grid is shifted by a random vector, whose components are uniformly distributed in the interval [−a/2,a/2].

4. Coupling of the swimmer's and fluid's dynamics

No velocity field is prescribed in our model of microswimmers. Locomotion is achieved by obeying the conservation of momentum in the collisions between the fluid particles and the swimmers; the shape asymmetry then induces self-propulsion. Two physical effects need to be included: we impose no-slip boundary conditions on the model swimmer's surface, and the force poles are explicitly included (see Fig. 1). Both effects induce modifications of the streaming and collision steps of the MPCD algorithm that we explain in the following.

a. Streaming step. To ensure the no-slip boundary condition the bounce-back rule⁶⁷ is applied to the MPCD particles that hit the spheres during the streaming step. The velocity of the fluid particle is reversed and the change in momentum is given by

J_i = 2m[v_i − U − Ω × (r̃_i − R)], (A18)

where R is the center-of-mass position of the swimmer colliding with the fluid particle, U and Ω are the linear and angular velocity of the swimmer, r̃_i is the position of the fluid particle upon collision with the sphere. The updated fluid particle velocity reads

v_i′ = v_i − J_i/m. (A19)

In addition, the fluid particles are reflected back along the direction of their initial velocity. For this, we use an exact ray tracing method to detect the collision of the MPCD particle onto the swimmer's surface. If a collision is detected the MPCD particle is propagated back onto the swimmer's surface and then the bounce-back rule is applied. If the particle travels a fraction λ, 0 < λ < 1 of the timestep towards the swimmer, then after the collision with the swimmer, it will travel for the time (1 − λ)δt away from the swimmer's surface. Furthermore, we allow for multiple collisions within the same timestep; this has been shown to prevent spurious depletion forces among colloids.⁶⁸ The new linear and angular velocities of the swimmer after the collision with the fluid particles read

The force poles are added as external force regions⁶⁹ in the streaming step for each swimmer. This is done by modifying the streaming step inside the force regions to

v_i(t + δt) = v_i(t) + f^lab_acδt, (A23)

where the force in the lab frame reads

f^lab_ac ≡ f_ac − [U + Ω × (r_i − R)]/δt, (A24)

and f_ac is the active force discussed in the following. The flow fields are modeled by force poles. While mathematically such force poles are point forces, any numerical implementation must mollify this requirement.

Pullers. The flow field is modeled by three Stokeslets, and the active force f^lab_ac is applied to R_pl,1, R_pl,2 and R_pl,2′ [see Fig. 6(a)]. The region R_pl,1 with diameter d_pl,1 is located at the rear of the swimmer and its force points into the direction of the swimmer's orientation. The other two regions R_pl,2 and R_pl,2′ are placed on the side of the swimmer and have the opposite orientation. The angle α_pl between the orientation of the puller and the line connecting the center of mass C and the midpoint of the region R_pl,2 (or R_pl,2′) defines their position on the boundary of the swimmer. The diameter of both R_pl,2 and R_pl,2′ is d_pl,2 = d_pl,1/(2)^1/3, such that they have half the volume of R_pl,1, making the fluid force free. The total force in the region R_pl,1 is and in the regions R_pl,2, R_pl,2′ is .
Pushers. The flow field is modeled by a force dipole. We apply the force f^lab_ac to all fluid particles located within spherical regions [see Fig. 6(b)]. The regions R_ps,1 and R_ps,2 where f^lab_ac is applied are equally sized spheres with diameter d_ps and the two forces f_ac are equal and opposite, to ensure that the fluid is overall force free. To generate a smooth flow on the boundary of the swimmer the direction of the applied force in regions R_ps,1 and R_ps,2 is modeled as follows. For fluid particles r_i ∈ R_ps,1 or R_ps,2 we apply the forceHere, s^b = (s^b_x,s^b_y,s^b_z)^T is the distance between the MPCD particle and the center of the region R_ps,1 or R_ps,2. The small, black arrows in Fig. 6(b) give a schematic representation of the flow field arising from eqn (A25). As before the superscript b denotes the body frame, in which the swimmer's orientation is aligned with the z axis. The constant f₀ gives the strength of the force that is applied. The total force in one of the two regions can be estimated by integrating eqn (A25), which yields . This takes into account the redirection of the force on the boundary of the swimmer and the density ρ of the fluid.

b. Collision step. To guarantee the no-slip boundary conditions on the surface of the swimmers, it is necessary to fill each swimmer with ghost particles, such that the collision step can be properly executed.⁶⁵ The positions of the ghost particles r^g_i are uniformly distributed within the swimmer,‡ and are advected with the swimmer in each timestep. The ghost particles density is matched to the fluid's density so as to make the swimmer neutrally buoyant. Before every collision step the ghost velocities v^g_i are updated according to

v^g_i = U + Ω × (r^g_i − R) + v^ran_i, (A26)

where the components of v^ran_i are sampled from a Gaussian distribution. The ghost particles then (together with the fluid particles) take part in the collision step [see eqn (A17)], and their velocities are updated to . The resulting change in linear momentum due to the ghost particles is and the change in angular momentum is L^g_i = (r^g_i − R) × J^g_i. These changes are then transferred to the swimmer¹⁷

5. Algorithm implementation

In this section we explain how the present algorithm is implemented. First, note that the MPCD algorithm scales as [scr O, script letter O](N), and thus is particularly prone to an efficient implementation with parallel programming. We therefore implemented the entire dynamics (both swimmer and MPCD) on graphics processing unit (GPU) cards using the Nvidia CUDA language. Parallelization of the MPCD algorithm is rather straightforward. The streaming step is performed for each fluid particle independently in a separate CUDA kernel, whereas ghost particles are simply translated together with the corresponding swimmer. Then, a kernel for each particle is started to carry out the bounce-back rule and afterwards a kernel for each particle is started to apply the periodic boundary conditions. The collision step eqn (A17) is implemented with the following kernels:

1. the cell of each particle is found and the center of mass for each cell is found by starting a kernel for each particle;

2. the two sums over the velocities in eqn (A17) are computed, where a kernel for each velocity component of each particle is started;

3. a kernel for each particle is started to compute the number of particles in each cell;

4. a kernel for each cell is started to normalize the velocities and compute the center of mass;

5. a kernel for each particle is started to compute its position with respect to the center of mass;

6. a kernel for each particle is started to compute the cross product in eqn (A17) and the contributions in each cell are summed;

7. six kernels for each particle are started to compute the contribution to the six components of the inertia tensor in the corresponding cell;

8. a kernel for each cell is started to compute the inverse of the respective inertia tensors;

9. for each cell a kernel is started to compute the curly bracket in eqn (A17);

10. a kernel for each particle is started to add all contributions of eqn (A17) and finish the collision step, while considering the rule for ghost particles eqn (A26);

11. a kernel for each ghost particle is started, which computes the momentum and angular momentum transfer [eqn (A27) and (A28)] to the swimmers.

Some of these computations could be combined into single kernels, but it is computationally more efficient to start many kernels with small computations, which we have opted for in this algorithm. After the MPCD steps the rigid body dynamics part (see also Section A1) of the code is executed. It consists of the following kernels, one for each swimmer:

1. the velocity and angular velocity stemming from the fluid interaction are added to the swimmers;

2. the positions and quaternions are updated;

3. periodic boundary conditions are applied;

4. the neighbor list between swimmers is updated (here we use a linked list);

5. the steric forces and torques are computed;

6. the velocities and angular velocities are updated.

This concludes the algorithm, which now goes back to the streaming step of the MPCD part.

6. Computational details

We carried out three-dimensional simulations with an average of 〈N_C〉 = 20 fluid particles per cell. The timestep of the MPCD algorithm is fixed to , whereas the MD timestep is . The resulting kinematic viscosity ν = η/ρ of the fluid for the MPC-AT+a algorithm (including both kinetic and collisional contribution) can be calculated exactly as .^47,65,70 Simulations using a forced flow (for details see ref. 69) produced a viscosity of .

The large sphere F associated to the stroke-averaged flagella of the swimmer has a diameter of d_f = 7a, while the small sphere B associated to the body of the swimmer has d_b = 3a, and the distance between the spheres centers is l = 7a. The choice of the geometrical parameters is dictated by a combination of factors. First, it is computationally convenient to make the swimmers' linear size a few times the grid spacing a. Second, inspired by the geometric properties of C. reinhardtii a ratio d_f/d_b ≳ 2 is advisable.¹⁰ For the sake of clarity in the comparison of our results, we maintain the same geometry also for pushers. The energy scale of the steric interactions is set to ε = 10k_BT. For pushers we fix the diameter of the force dipole regions R_ps,1 and R_ps,2 to d_ps = 3a. The region R_pl,1 of the pullers has the same diameter d_pl,1 = 3a and accordingly the regions R_pl,2 and R_pl,2′ have the diameter d_pl,2 = 3a/(2)^1/3. The angle between the swimmers orientation and the line connecting the center of mass of the pullers C to the midpoint of the regions R_pl,2 and R_pl,2′ is α_pl = 107°.

To initialize the simulations, we distribute the swimmers homogeneously (and without overlaps) across a cubic box with periodic boundary conditions.

Appendix B: estimate of maximum from theoretical prediction

To estimate the maximum predicted from the analytical treatment, we computed the maximum of the eigenvalues in eqn (30) for pushers and eqn (25) for pullers. As a first step, we have to find estimates for both the wavenumber k and the constant quantifying the steric interactions ζ. Since we expect the hydrodynamic interactions to be relevant on the size of the swimmer, we chose the distance between the force poles l to determine the characteristic wavenumber k = 2π/l. The constant ζ quantifying the steric interactions can be estimated from ref. 58

ζ = v₀²σ_sτ_c (A29)

where v₀ is the propulsion speed, σ_s the geometrical cross section, and τ_c is the collision time. The collision time can be estimated from the center-of-mass distance of two colliding swimmers, where we considered collisions as shown in Fig. S5 and S6 (ESI†). The resulting collision time for pullers is and for pushers . The propulsion speed is extracted from Fig. 3 and the geometrical cross section is σ_s = 4πσ². Here, σ = l/2 = 7.5a is the effective steric radius of the swimmer also estimated from the distance of the force poles. To finally obtain a filling fraction we use ϕ = c₀(V_B + V_F − V_Ol) [see eqn (3)].

Acknowledgements

We gratefully acknowledge insightful conversations with Johannes Blaschke, Jens ElgetiStephan Herminghaus, Sebastian Mair, and Kuang-Wu Lee. We also gratefully acknowledge support from the Deutsche Forschungsgemeinschaft (SFB 937, project A20). Open Access funding provided by the Max Planck Society.

References

1. A. Baskaran and M. C. Marchetti, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 15567–15572.
2. B. Ezhilan, M. J. Shelley and D. Saintillan, Phys. Fluids, 2013, 25, 070607.
3. J. Elgeti, R. G. Winkler and G. Gompper, Rep. Prog. Phys., 2015, 78, 056601.
4. E. Lauga, W. R. DiLuzio, G. M. Whitesides and H. A. Stone, Biophys. J., 2006, 90, 400–412.
5. A. P. Berke, L. Turner, H. C. Berg and E. Lauga, Phys. Rev. Lett., 2008, 101, 038102.
6. T. Ostapenko, F. J. Schwarzendahl, T. J. Böddeker, C. T. Kreis, J. Cammann, M. G. Mazza and O. Bäumchen, Phys. Rev. Lett., 2018, 120, 068002.
7. H. Wioland, F. G. Woodhouse, J. Dunkel, J. O. Kessler and R. E. Goldstein, Phys. Rev. Lett., 2013, 110, 268102.
8. H. Wioland, E. Lushi and R. E. Goldstein, New J. Phys., 2016, 18, 075002.
9. M. F. Copeland and D. B. Weibel, Soft Matter, 2009, 5, 1174–1187.
10. 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–14313.
11. C. Dombrowski, L. Cisneros, S. Chatkaew, R. E. Goldstein and J. O. Kessler, Phys. Rev. Lett., 2004, 93, 098103.
12. E. M. Purcell, Am. J. Phys., 1977, 45, 3–11.
13. M. J. Lighthill, Commun. Pure Appl. Math., 1952, 5, 109–118.
14. J. R. Blake, J. Fluid Mech., 1971, 46, 199–208.
15. T. Ishikawa and T. J. Pedley, Phys. Rev. Lett., 2008, 100, 088103.
16. M. T. Downton and H. Stark, J. Phys.: Condens. Matter, 2009, 21, 204101.
17. I. O. Götze and G. Gompper, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 82, 041921.
18. I. Llopis and I. Pagonabarraga, J. Non-Newtonian Fluid Mech., 2010, 165, 946–952.
19. A. A. Evans, T. Ishikawa, T. Yamaguchi and E. Lauga, Phys. Fluids, 2011, 23, 111702.
20. A. Zöttl and H. Stark, Phys. Rev. Lett., 2012, 108, 218104.
21. F. Alarcón and I. Pagonabarraga, J. Mol. Liq., 2013, 185, 56–61.
22. K. Ishimoto and E. A. Gaffney, Phys. Rev. E, 2013, 88, 062702.
23. J. J. Molina, Y. Nakayama and R. Yamamoto, Soft Matter, 2013, 9, 4923–4936.
24. A. Zöttl and H. Stark, Phys. Rev. Lett., 2014, 112, 118101.
25. M. Theers, E. Westphal, G. Gompper and R. G. Winkler, Soft Matter, 2016, 12, 7372–7385.
26. J. Blaschke, M. Maurer, K. Menon, A. Zöttl and H. Stark, Soft Matter, 2016, 12, 9821–9831.
27. F. Alarcon, C. Valeriani and I. Pagonabarraga, Soft Matter, 2017, 13, 814–826.
28. T. Ishikawa, J. T. Locsei and T. J. Pedley, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 82, 021408.
29. G.-J. Li and A. M. Ardekani, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 013010.
30. I. Pagonabarraga and I. Llopis, Soft Matter, 2013, 9, 7174–7184.
31. R. Matas-Navarro, R. Golestanian, T. B. Liverpool and S. M. Fielding, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 032304.
32. J. de Graaf, H. Menke, A. J. T. M. Mathijssen, M. Fabritius, C. Holm and T. N. Shendruk, J. Chem. Phys., 2016, 144, 134106.
33. L. P. Fischer, T. Peter, C. Holm and J. de Graaf, J. Chem. Phys., 2015, 143, 084107.
34. J. de Graaf, T. Peter, L. P. Fischer and C. Holm, J. Chem. Phys., 2015, 143, 084108.
35. R. W. Nash, R. Adhikari, J. Tailleur and M. E. Cates, Phys. Rev. Lett., 2010, 104, 258101.
36. R. W. Nash, R. Adhikari and M. E. Cates, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 77, 026709.
37. J. P. Hernandez-Ortiz, C. G. Stoltz and M. D. Graham, Phys. Rev. Lett., 2005, 95, 204501.
38. J. Stenhammar, C. Nardini, R. W. Nash, D. Marenduzzo and A. Morozov, Phys. Rev. Lett., 2017, 119, 028005.
39. L. F. Valadares, Y.-G. Tao, N. S. Zacharia, V. Kitaev, F. Galembeck, R. Kapral and G. A. Ozin, Small, 2010, 6, 565–572.
40. D. Saintillan and M. J. Shelley, Phys. Rev. Lett., 2007, 99, 058102.
41. R. Singh, S. Ghose and R. Adhikari, J. Stat. Mech.: Theory Exp., 2015, 2015, P06017.
42. J. W. Swan, J. F. Brady, R. S. Moore and C. 174, Phys. Fluids, 2011, 23, 071901.
43. D. Saintillan and M. J. Shelley, J. R. Soc., Interface, 2011 DOI:10.1098/rsif.2011.0355.
44. K. Drescher, J. Dunkel, L. H. Cisneros, S. Ganguly and R. E. Goldstein, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 10940–10945.
45. K. Drescher, R. E. Goldstein, N. Michel, M. Polin and I. Tuval, Phys. Rev. Lett., 2010, 105, 168101.
46. A. Malevanets and R. Kapral, J. Chem. Phys., 1999, 110, 8605–8613.
47. G. Gompper, T. Ihle, D. M. Kroll and R. G. Winkler, Multi-particle collision dynamics: a particle-based mesoscale simulation approach to the hydrodynamics of complex fluids, Springer, 2009, vol. Advanced computer simulation approaches for soft matter sciences III, pp. 1–87.
48. T. Ishikawa, J. T. Locsei and T. J. Pedley, J. Fluid Mech., 2008, 615, 401–431.
49. D. Saintillan and M. J. Shelley, Phys. Rev. Lett., 2008, 100, 178103.
50. P. T. Underhill, J. P. Hernandez-Ortiz and M. D. Graham, Phys. Rev. Lett., 2008, 100, 248101.
51. A. Wysocki, J. Elgeti and G. Gompper, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91, 050302.
52. E. Lushi, V. Kantsler and R. E. Goldstein, Phys. Rev. E, 2017, 96, 023102.
53. S. E. Spagnolie and E. Lauga, J. Fluid Mech., 2012, 700, 105–147.
54. A. Zöttl and H. Stark, J. Phys.: Condens. Matter, 2016, 28, 253001.
55. J. K. G. Dhont, An introduction to dynamics of colloids, Elsevier, 1996.
56. C. H. Rycroft, Chaos, 2009, 19, 041111.
57. J. Bialké, H. Löwen and T. Speck, EPL, 2013, 103, 30008.
58. J. Stenhammar, A. Tiribocchi, R. J. Allen, D. Marenduzzo and M. E. Cates, Phys. Rev. Lett., 2013, 111, 145702.
59. H. Goldstein, C. P. Poole and J. L. Safko, Classical Mechanics, Addison-Wesley, 3rd edn, 2001.
60. I. P. Omelyan, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1998, 58, 1169–1172.
61. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford University Press, 1989.
62. J. D. Weeks, D. Chandler and H. C. Andersen, J. Chem. Phys., 1971, 54, 5237–5247.
63. I. N. Bronstein and K. A. Semendjajew, Taschenbuch der Mathematik, BSB B. G. Teubner Verlagsgesellschaft, Nauka-Verlag, Leipzig, Moskau, 19th edn, 1979.
64. H. Noguchi, N. Kikuchi and G. Gompper, Europhys. Lett., 2007, 78, 10005.
65. I. O. Götze, H. Noguchi and G. Gompper, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 046705.
66. T. Ihle and D. M. Kroll, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2001, 63, 020201.
67. A. Lamura, G. Gompper, T. Ihle and D. M. Kroll, Europhys. Lett., 2001, 56, 319.
68. J. T. Padding and A. A. Louis, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 74, 031402.
69. D. S. Bolintineanu, J. B. Lechman, S. J. Plimpton and G. S. Grest, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 86, 066703.
70. H. Noguchi and G. Gompper, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 016706.

---

Footnotes

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

‡ We recommend to fill the swimmer with multiple ghost particles, rather than using a single ghost particle of large mass, as this would result in a wrong value of the torque.

---