Maximum in density heterogeneities of active swimmers

Suspensions of unicellular microswimmers such as flagellated bacteria or motile algae exhibit spontaneous density heterogeneities at large enough concentrations. Based on the relative location of the biological actuation appendages i.e. flagella or cilia) microswimmers' propulsion mechanism can be classified into two categories: (i) pushers, like \textit{E. coli} bacteria or spermatozoa, that generate thrust in their rear, push fluid away from them and propel themselves forward; (ii) pullers, like the microalgae \textit{Chlamydomonas reinhardtii}, that have two flagella attached to their front, pull the fluid in and thereby generate thrust in their front. 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\'eclet 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\'eclet 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][2][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 the spontaneous formation of spiral vortices [4], directed motion [5], swarming [6], bacterial turbulence [7], complex interaction with solid surfaces [8,9], and self-concentration [10].
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 R = avρ/η, where a 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., a ≈ 10 µm, v ≈ 30 µm/s, and for water ρ ≈ 10 3 kg/m 3 , η ≈ 10 −3 Pa s, which result in R ≈ 10 −5 . As noted by Purcell [11], 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 be-cause 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 [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29], the shape anisotropic raspberry swimmer [30][31][32], the force-counterforce model [33][34][35][36], the catalytic dimers [37], or other hydrodynamic models [38][39][40]. Experiments have confirmed that the flow field of flagellated bacteria like E. coli is to very good approximation modeled by a simple force dipole [41], whereas Chlamydomonas are modeled by three Stokeslets [42]. 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 microswimmmer 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 [43]. 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 [44], and active swimmers [19,23,25]. 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 show 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,[45][46][47]. 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 suppress 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 Sec. II we introduce the model for the microswimmer, the implementation of the fluid, and its coupling to the former. Section III describes the physical properties of the fluid and the microswimmer's flow field. In Sec. 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 suppressed by steric interactions at higher filling fractions. In Sec. V we present the analytical theory and show that we also find a maximum heterogeneity, which is mediated by the interplay of hydrodynamic and steric interactions. Finally, in Sec. VI we discuss our main results and summarize our conclusions.

II. MODEL
We employ a stroke-averaged model of biological microswimmers, similarly to [1,48], taking into account the asymmetric shape of biological microswimmers due to the cell's body and the flagella. Thus, the swimmer is modeled as an asymmetric dumbbell, as depicted in Fig. 1, that mimics a Chlamydomonas 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. We simulate the fluid surrounding the swimmer using MPCD (see Sec. II B). Precise measurements [41] show that pusher-type microswimmers are well modeled by a force dipole, represented by the red regions with embedded arrows in Fig. 1(b). Since the swimmer is shape-asymmetric the hydrodynamic center is shifted with respect to the center of mass and henceforth a symmetry breaking on the swimmer gives propulsion. Pullers, on the other hand, are well represented by a three-Stokeslet solution of the Stokes equation [42] [red regions in Fig. 1(a)].
The microswimmer is characterized by its mass M , center of mass position R and orientation q. The center of mass position in the swimmer's body frame and the inertia tensor I m are calculated in App. A. In the following we describe the equations governing the motion of the microswimmer, the fluid dynamics implemented through the MPCD, and their coupling.

A. 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 around an axis (the Mozzi axis) and a rotation around the same axis, as per the Mozzi-Chasles theorem [49]. 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, 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 P 3 (R) which admits a universal cover represented by the group of unit quaternions q = (q 0 , q 1 , q 2 , q 3 ) T , where the superscript T indicates the matrix transposition. In App. B we provide the basic definitions of quaternion algebra necessary for the following.
The equations of motion for the rigid body dynamics in three dimensions read [50] where Ω 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 Eq. (1), 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 given in App. B. The repulsive, steric interactions among swimmers are modeled using a Weeks-Chandler-Andersen potential [51] Φ(r ij,ab ) = 4 σ ab r ij,ab 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 [50], which was also used and discussed in detail in [24]. Given a vector in the laboratory frame f the transformation to the body frame vector f b is given by where the matrix D(q) is constructed from the quaternions and given in App. B. Thus, the orientation of the swimmer at any given time is found from D −1 (q(t))(0, 0, 1) T Note that all quantities that do not carry an index b are calculated in the laboratory frame.

B. 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 [43] to simulate a fluid at the Navier-Stokes level of description. We include the Anderson thermostat and the con-servation of angular momentum into the MPCD dynamics; the resulting algorithm is usually denoted as MPC-AT+a [44,52,53]. The fluid is modeled using N fl pointlike 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 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 C(i) 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 [52] where the random velocity v ran i has components distributed according to a Gaussian distribution with zero mean and variance k B T /m, k B is the Boltzmann constant, N C(i) is the number of fluid and ghost particles (see Sec. II C) in cell C(i). The vector r j,c is the position of the neighboring particle j relative to the center of mass of the cell C(i). In Eq. (8), Π −1 is the inverse of the moment of inertia tensor Π ≡ j∈C(i) m [(r j · r j )I − r j ⊗ r j ] for the fluid particles in cell C(i), where I is the identity tensor, '·' is the scalar product and '⊗' the tensor product. Note that Π −1 is a dynamical quantity that has to be updated at every timestep, and it also includes the ghost particles within the swimmer (see Sec. II C).
To ensure Galilean invariance and avoid the build-up of spurious correlations in the velocities [54], 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].

C. 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.

Streaming step
To ensure the no-slip boundary condition the bounceback rule [55] 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 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 velocity reads 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 swimmers surface and then the bounce-back rule is applied. 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 (see [56] for external force) in the streaming step for each swimmer. This is done by modifying the streaming step inside the force regions to where the force in the lab frame reads and f ac is the active force discussed in the following. The flow fields are modeled by a 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. 2(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 swimmers orientation. The other two regions R pl,2 and R pl,2 are aligned 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.
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. 2(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 force Here, T is the distance between the MPCD particle and the center of the region r i ∈ R ps,1 or R ps,2 . As before the superscript b denotes the body frame, in which the swimmers orientation is aligned with the z axis. The constant f 0 gives the strength of the force that is applied and f 0 < 0 in R ps,1 and f 0 > 0 in R ps,2 .

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 [53]. The positions of the ghost particles r g i are uniformly distributed within the swimmer 1 , and are advected with the swimmer in each timestep. The ghost particles density is matched to the fluids density so as to make the swimmer neutrally buoyant. Before every collision step the ghost velocities v g i are updated according to where the components of v ran To conclude this section on our computational model, we note that this algorithm scales as O(N ), and thus is particularly prone to an efficient implementation with parallel programming. We therefore implemented the entire dynamics on graphics processing unit (GPU) cards.

D. 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 δt = 10 −2 ma 2 /(k B T ), whereas the MD timestep is δt MD = 5 × 10 −4 ma 2 /(k B T ). The resulting kinematic viscosity ν = η/ρ of the fluid for the MPC-AT+a algorithm (including both kinetic and collisional contribution) can be calculated exactly as ν = 3.88a k B T /m [44,53,57]. Simulations using a forced flow (for details see [56]) produced a viscosity of ν = 3.69a k B T /m.
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 Chlamydomonas a ratio d f /d b 2 is advisable [58]. 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 B T . 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 = 0.31.
To initialize the simulations, we distribute the swimmers homogeneously across a cubic box with periodic boundary conditions.

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) immersed in the MPCD fluid, that is, we carried out equilibrium simulations without activity f lab ac = 0. The equipartition theorem applied to the passive colloid for the translational and rotational motion predicts Hydrodynamics at low Reynolds numbers (relevant for micron-sized objects) allows a great simplification of the Navier-Stokes equations: the nonlinear, inertial effects can be neglected, and the governing equations are the Stokes equations where v(r) is the fluid velocity, p(r) the pressure, and f ext (r) is a body force acting on the fluid. Solving the Stokes equations means obtaining expressions for v and p that satisfy Eq. (20) 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 In free three-dimensional space, the Greens function is found by considering a point force f ext = f e δ(r) in an infinite, quiescent fluid, where e is the unit vector representing the direction of the force. A straightforward calculation [61] gives the Oseen tensor O(r) ≡ 1 8πηr (I +r ⊗r) wherer ≡ r/r, r = |r|, and the resulting flow field v(r) = f 8πηr [e + (r · e)r], which is termed a 'Stokeslet' and decays with distance as r −1 .
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 from the simulations. For pushers, the force in one of the two regions can be estimated by integrating Eq. (16), which yields f = f 0 5 48 πd 3 ps ρ. This takes into account the redirection of the force on the boundary of the swimmer and the density ρ of the fluid. For pullers the force in the region R pl,1 is f = f 0 1 6 πd 3 pl,1 ρ and in the regions R pl,2 , R pl,2 is f = f 0 1 6 πd 3 pl,2 ρ. 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 0 = 50k B T /a and carry out the full dynamics as described in Sec. II. Figure 3 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. 3(a)]. The situation is reversed in the pusher case [ Fig. 3(c)] where fluid is pushed out at the front and back of the swimmer. The theoretical predictions for both the puller [ Fig. 3(b)] and the pusher [ Fig. 3(d)] show the same characteristic behavior as the simulated flow fields.
The effective velocity v eff ≡ | e · U | of an isolated swimmer in the steady state depends linearly on the active force f 0 [1]. From our simulations, we calculate v eff for a pusher [see Fig. 4]. The linear fit has a slope of The analogous results for pullers are also shown in Fig. 4, where the slope of the linear fit is α = (4.4 ± 7 × 10 −2 ) × 10 −3 ma 2 k B T . For both pushers and pullers the maximal effective velocity we investigate is v eff ∼ 0.1 k B T /m, therefore the highest Reynolds number is Re = f0ασ ν ∼ 0.1, and thus we are in the low Reynolds number regime.

IV. DENSITY HETEROGENEITIES
We carry out simulations with N = 300 − 1560 swimmers and study the phase behavior of the active swimmers. The filling fraction we denote here is computed using the volume of the swimmers 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 Ol Here V Box is the volume of the simulation box. 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 Furthermore, we vary the strength of the active force f 0 , 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 using where we used the linear relation (fitted slope) between active velocity and force dipole strength from Sec. III. Furthermore, σ = 5a is the typical length of the swimmer and for the diffusion constant we assume D = k B T 6πησ . We analyze the systems 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 . We compare σ loc to the standard deviation of local Voronoi volumes for random homogeneous configurations σ rnd of the corresponding filling fraction. Figure 5 shows the resulting phase diagram, where the ratio σ loc /σ rnd as a function of Péclet number and filling fraction is shown. Here, positive values of the Péclet number correspond to pusher type swimmers, whereas negative values are puller swimmers. The phase diagram shows that for both pullers and pushers, initially σ loc /σ rnd grows with Péclet number 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 [1,2,[45][46][47]. 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 potential [Eq. (5)] 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 App. C.
In Fig. 6 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 an increase of σ loc /σ rnd , which is much less pronounced. The original simulations are at an intermediate regime. 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.

V. THEORETICAL ANALYSIS
To bolster our numerical results, we develop an analytical theory of microswimmers that explicitly includes hydrodynamic and steric interactions. As suggested in [1] swimmers are modeled by an asymmetric dumbbell, similar to our numerical model (see Sec. II), where the motion of the swimmer is described by the following effective Langevin equations where r Li is the position of the front sphere of the swimmer i and r Si the position of the respective back sphere. The front and back spheres of each swimmer are connected by an infinitely thin rigid rod. The motion is coupled to the fluid velocity v, which is determined by the Stokes equation including a stochastic and an active force term Here, the active force is given by a force dipole which points along the orientation of the swimmer e i and has a force strength f . Furthermore, fluctuations in the swimmers motion are added to the fluid via where ξ L,S i (t) are noise terms with ξ L,S i (t)ξ L,S β (t ) = 2Γ L,S Ik B T δ iβ δ(t − t ), and Γ L,S = 6πηa L,S are the friction coefficients of the front and back sphere. Following [1], we derive the one-body Smoluchowski equation from Eq. (24)- (28). As we have seen in Sec. 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 [62,63] v(c) = v 0 − cζ.
Here, v 0 is the bare swimming velocity and c is the concentration. The constant ζ quantifies how much the swimmers are slowed down by the steric interactions (for more details see [62,63] where 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 Eq. (30) 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 diffusivity, with translational diffusion constant D and rotational diffusion constant D R . To make progress with this equation, we follow the standard path and compute moment equations for the concentration c, the polarization P and the nematic order tensor Q c (r, t) = de p(r, e, t), the "splay" δQ =k · δQ ·k and the "bend" δQ ⊥ = k · δQ · δ −kk components of the nematic tensor. To first order in the fluctuations the equations governing the temporal evolution read Here, k = |k| is the absolute value of the wavevector.
The constants γ 1 f and γ 2 f (see App. D) stem from the hydrodynamic interactions between the swimmers, which are dominated by the force dipole strength f 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).
For pullers (f < 0) the bend fluctuations δQ ⊥ are stable. To analyze the fluctuations in the concentration we use a large length-scale and long time-scale approximation for the longitudinal polarization and the splay component of the nematic order tensor Inserting Eq. (39)-(40) into the Eq. (34) and keeping terms up to O(c 2 0 ) yields Since for pullers f < 0 the term γ1D 4D R f c 0 introduces an instability at low concentrations c 0 , which are counteracted by the term v0ζ D R c 2 0 , which stabilizes the system at higher concentrations. Since the first term stems from the hydrodynamic interactions and the second term from 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 0 . Moreover, we find a maximum of instability as a function of c 0 and, hence, a maximum heterogeneity.
For pushers (f > 0) the splay fluctuations are stable, whereas the bend fluctuations become unstable. In the same large length and time scale limit the transverse polarization becomes Inserting this into the Eq. (37) gives Here, the term −γ 2 f c 0 destabilizes the system at low concentrations c 0 through bend fluctuations, which are counteracted by the term k 2 1 10D R v 2 0 c 2 0 , that stabilizes the system for higher concentrations. Again, the first term, which destabilizes the system, stems from the hydrodynamic interactions, whereas the second, stabilizing term comes from the steric interactions. Considering only terms up to O(c 2 0 ), we also find a maximum of instability, which translates to a maximum heterogeneity in the system.

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 Chlamydomonas or E. coli cells. Furthermore, our model takes into account the anisotropic shape of a microswimmer. Typical for this is the shape of a Chlamydomonas cell, which is well modeled by a dumbbell [58]. Our model allows for great flexibility in the geometrical modeling of microswimmer shapes. Self-propulsion is generated through a symmetry breaking due to the asymmetric shape and force free motion. The fluid is 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 [41,42]. Additionally, we test the effective velocity of the microswimmer model, and find that it depends linearly on the applied force.
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,[45][46][47]. 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.
Additionally, we include steric as well as hydrodynamic interactions into an analytical theory, based on a Smoluchowski equation. We computed the hydrodynamic moments of this equation and performed a linear stability analysis of the moment equations around the homogeneous sate. For both puller and pusher-type swimmers, we found that at low concentration the system is destabilized by hydrodynamic interactions. At higher concentrations the instabilities are counteracted by the steric interaction. This gives rise to a maximum in the instability of the homogeneous state, and thus a maximum heterogeneity in the concentration of swimmers.
The 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 might have important biological implications: it means that there is an optimal filling fraction and Péclet number for the formation of heterogeneous structures. Bacteria or microalgae exhibiting these optimal parameters are more likely to form colonies or biofilms.
With the height h i of the spherical caps The moment of inertia for a spherical cap in the x-as well as y-direction is In the z-direction it is By using the moment of inertia of a sphere I Spi = (A7)