Clustering of microswimmers: Interplay of shape and hydrodynamics

The spatiotemporal dynamics in systems of active self-propelled particles is controlled by the propulsion mechanism in combination with various direct interactions, such as steric repulsion, hydrodynamics, and chemical fields. Yet, these direct interactions are typically anisotropic, and come in different 'flavors', such as spherical and elongated particle shapes for steric repulsion, pusher and puller flow fields for hydrodynamics, etc. The combination of the various aspects is expected to lead to new emergent behavior. However, it is a priori not evident whether shape and hydrodynamics act synergistically or antagonistically to generate motility-induced clustering (MIC) and phase separation (MIPS). We employ a model of prolate spheroidal microswimmers - called squirmers - in quasi-two-dimensional confinement to address this issue by mesoscale hydrodynamic simulations. For comparison, non-hydrodynamic active Brownian particles (ABPs) are considered to elucidate the contribution of hydrodynamic interactions on MIC and MIPS. For spherical particles, the comparison between ABP and hydrodynamic-squirmer ensembles reveals a suppression of MIPS due to hydrodynamic interactions. The fundamental difference between ABPs and squirmers is attributed to an increased reorientation of squirmers by hydrodynamic torques during their collisions. In contrast, for elongated squirmers, hydrodynamics interactions enhance MIPS. Thus, hydrodynamic interactions show opposing effects on MIPS for spherical and elongated microswimmers.

Computer simulations of model systems, such as spherical active Brownian particles (ABPs), yield motilityinduced clustering (MIC) and motility-induced phase separation (MIPS) [19,[21][22][23][24][25][26] in qualitative agreement with the above-mentioned experiments on synthetic selfphoretic particles.The intuitive explanation for the emergence of MIC and MIPS is a positive feedback between blocking of persistent particle motion by steric interactions, and an enhanced probability of collisions with further particles at sufficiently large concentrations and activities [20,27,28].However, in self-phoretic systems, the phoretic field has been show to significantly contribute to clustering by an field-induced attractive interaction [29,30].
Cluster formation, swarming, and mesoscale turbulence of anisotropic objects such as bacteria is explained by steric interaction-induced alignment.Computer simulations of two-dimensional assemblies of rods suggest that the interplay of rod geometry, self-propulsion, and steric interaction suffices to facilitate aggregation into clusters [3,[31][32][33] and even lead to swarming and turbulence [12,13].
Aside from steric, other interaction mechanisms between active particles are present, evidently and most prominently fluid-mediated interactions, since motile bacteria and phoretic microswimmers propel themselves via the embedding fluid.Hence, an crucial aspect regarding MIC and MIPS of microswimmers is the role of hydrodynamics.At equilibrium, hydrodynamic interactions (HIs) solely affect the dynamics and not the structural properties of a system.Since active systems are intrinsically out of equilibrium, this does not necessary apply, and the steady-state-hence the formation of clusters-is strongly affected by hydrodynamics [27,[34][35][36][37].
The phase behavior of active particles in the presence of hydrodynamic interactions has received much less attention than ABPs.In simulations, this is largely due to the substantial computational challenges posed by hydrodynamics.Specifically, the long-ranged flow field created by the swimmers has to be accounted for adequately.In experiments, HI is difficult to switch off, and therefore its contribution hard to assess.In contrast, simulations can be performed with and without HI and the latter effect on MIPS be analysed, although with an increased numerical effort.
A frequently applied model for a self-propelled particles in the presence of a fluid is a squirmer [34,[38][39][40][41][42][43]-a spheroidal colloid with a prescribed slip velocity on its surface.It was originally introduced to describe ciliated microswimmers such as Paramecia and Volvox.Nowadays, it is applied to a broad class of microswimmers, both, biological as well as synthetic ones.A squirmer is typically characterized by two modes accounting for its swimming velocity and its active stress.The latter distinguishes between pushers (e.g., E. coli, sperm), pullers (e.g., Chlamydomonas), and neutral squirmers (e.g., Paramecium) [43].
In order to shed light onto the importance of hydrodynamics for structure formation in active systems, we perform mesoscale hydrodynamic simulations applying the multiparticle collision dynamics approach (MPC), a particle-based hydrodynamic simulation method to solve the Navier-Stokes equations [44][45][46].Specifically, we consider squirmers confined in a thin slit, i.e., in a quasi-2D geometry.Here, we want to address and clarify several issues.On the one hand, we want to resolve the contradiction between Refs.[34] and [27], which predict enhanced clustering of spherical squirmers or no MIPS in the presence of HIs, respectively.On the other hand, we consider squirmers of prolate spheroidal shape.This is motivated by the different mechanism of cluster formation for spherical and elongated particles (blocking vs. alignment).Thus, the effect of hydrodynamic interactions on the collective behavior of spherical swimmers can be qualitatively different than for swimmers with elongated shapes.
As an important result, we find no MIPS for spherical squirmers in quasi-2D, in agreement with Ref. [27] and [37], but in contrast to Ref. [34].We attribute the discrepancy to peculiarities of the compressible MPC fluid employed in Ref. [34], which requires a suitable choice of the parameters for simulations of low-compressibility fluids.Hence, HIs suppress cluster formation and MIPS of self-propelled spherical squirmers compared to ABPs.Interestingly, the opposite is true for spheroidal squirmers with sufficiently large aspect ratio, where we find a substantial enhancement of MIPS by HIs.This applies for pushers, pullers, and neutral squirmers, as long as the force dipole is sufficiently weak.In case of large force dipoles, pushers show the least tendency to phase separa-tion, with even ABPs exhibiting stronger cluster formation.A density-aspect-ratio phase diagram for moderate force dipole strength shows most pronounced phase separation for pullers, followed by neutral squirmers, pushers, and finally ABPs.Thus, the effects of shape and hydrodynamics are highly non-additive, and act sometimes synergistically, sometimes antagonistically.
The paper is structured as follows.Previous results of squirmer simulations of collective phenomena are briefly summarized in Sec.II.Section III introduces the model for prolate squirmers.A comparison between MPC and ABP simulation results of spherical squirmers is presented in Sec.IV.Section V presents results for prolate spheroidal swimmers, and Sec.VI discusses the various obtained aspects.Section VII summarizes our results.Appendix A presents the simulation approach for the fluid (MPC) (A 1), the coupling of a squirmer with the fluid (A 2), as well as the simulation method for ABPs (A 3).Appendix B collects definitions relevant for rigidbody dynamic of ABPs, and App.C discusses technical aspects of fluid compressibility in particle-based simulations of squirmers.

II. COLLECTIVE BEHAVIOR OF SQUIRMERS-BRIEF SUMMARY OF PREVIOUS SIMULATION STUDIES
The dynamics of spheroidal squirmers in two dimensions (2D) and in a three-dimensional slit geometry (quasi-2D) has been studied and compared with ABPs in Refs.[27,37] and [34,47], respectively.In Ref. [27], the squirmer dynamics is strictly two dimensional, but embedded in a 3D fluid.The shape of the squirmers is changed from infinitely long cylinders (perpendicular to the 2D plane), corresponding to 2D HIs, to highly flattened cylinders with 3D HIs.Independent of the cylinder length (and dimensionality of HIs), the simulation studies show no evidence of phase separation, in contrast to comparable simulations of ABPs.The qualitatively different behavior is attributed to a faster decorrelation of the squirmer's swimming direction due to HIs compared to the configuration-independent rotational diffusion of ABPs.The studies of Ref. [37] emphasizes the importance of near-field HIs in the suppression of MIPS.However, clustering of the squirmers is observed and the formation of polar order for nearly neutral squirmers.Note that in both Refs.[27,37], thermal fluctuations are neglected.Hence, there is no thermal rotational diffusion leading to a decorrelation of the squirmers orientational motion.In the simulations of Refs.[34,47], squirmers are confined in a narrow slit (quasi-2D), slightly wider then their diameter, which captures the geometry of experiments, where glass plates are used for confinement, and thermal fluctuations are taken into account.Here, the swimming direction is free to orient in three dimensions, and hydrodynamics is screened and decays as 1/ 2 parallel to the confining surfaces in the slit center, where is the radial distance from the squirmer [48].Most importantly, hydrodynamics is argued to enhances MIPS, in contrast to findings of Refs.[27,37], and a phase diagram is presented in Ref. [47].
A quasi-2D setup is certainly less favorable for cluster formation of ABPs compared to a strict 2D system, since the swimmer's propulsion direction can point toward a wall, which reduces steric blocking and the lateral swim pressure [49][50][51].Thus, compared to strictly 2D ABP systems, ABPs in quasi-2D require larger Péclet numbers and packing fractions to form clusters, and those clusters exhibit a finite life time.Surface hydrodynamic interactions may lead to preferred squirmer orientations.Here, considerations based on far-field HIs predict a parallel alignment of the propulsion direction for pushers and a normal orientation for pullers at a no-slip surface [52].Hence, quasi-2D confinement may result in a different behavior, and near-field effects (specifically with the confining walls) might be important for cluster formation as indicated by the studies of Refs.[34,43].
Suspensions of athermal spherical squirmers in 3D, lattice-Boltzmann [53] and force-coupling [35] simulations have also found a pronounced polar order for pullers and neutral squirmers and clustering over a certain range of force-dipole strengths.However, no largescale MIPS is observed.In contrast, ABPs in 3D clearly exhibit MIPS [24][25][26] and large-scale collective motion [24].Moreover, simulations of 2D, quasi-2D, and 3D systems of attractive squirmers show a substantially enhanced cluster formation compared to purely repulsive squirmers [35,36].Also in such situations, hydrodynamics is found to reduce MIPS compared to ABPs.
So far, the influence of HIs on MIPS has mainly been studied for spherical squirmers.In contrast, very little is known about the influence of HIs on the structure and dynamics of elongated spheroidal squirmers confined in narrow slits.Here, the shape, squirmer-squirmer HIs, and HIs with the confining surfaces affect cluster formation and a possible MIC and MIPS [43].

III. MODEL AND SIMULATION APPROACH
We model a nonspherical squirmer as a prolate spheroidal rigid body with the prescribed surface velocity (cf.Fig. 1) [43,[54][55][56] Here, we employ spheroidal coordinates (ζ, τ, ϕ) (cf.Fig. 1), which are related to Cartesian coordinates via and 0 ≤ ϕ ≤ 2π, which is an extension of and alternative to Illustration of normal and tangent vectors of a spheroidal (left) and spherical (right) squirmer.Selfpropulsion along the body-fixed orientation vector e, here e = (0, 0, 1) T , is achieved by an axisymmetric prescribed surface velocity in the direction of the tangent vector s [43].

FIG. 2.
Sketch of spheroidal microswimmers in a narrow slit.The ratio of the squirmer minor axis and the width of the slit is 2bx/Ly = 6/7.The arrows indicate the swimming direction along a squirmers major axis.
earlier approaches [54,56].The advantage of the choice (1) for the slip velocity is that an analytical solution of the flow field can be obtained [43].The constant B 1 in Eq. ( 1) determines the amplitude of the slip velocity and is related to the self-propulsion velocity U 0 via [43] where τ 0 = b z /c.The active stress is captured by the coefficient β [43,57].Thereby, β < 0 corresponds to a pusher, β > 0 to a puller, and β = 0 to a neutral squirmer.
We consider the collective dynamics of squirmers confined in a narrow slit of dimensions L×L y ×L (cf.Fig. 2), where the ratio of the squirmer minor axis and the width of the slit is 2b x /L y = 6/7, i.e., we focus on a quasi-2D geometry.
The squirmers are embedded in a fluid, which we simulate by multiparticle collision dynamics (MPC), a mesoscale hydrodynamics simulation approach [44][45][46].The details of the algorithm are described in App.A 1, and the implementation of the squirmer and the coupling with MPC in App.A 2.
For comparison, we perform simulations of active Brownian spheres and spheroids.The respective simulation approach is described in App.A 3.

IV. SIMULATIONS-SPHERICAL SQUIRMERS A. System Setup and Parameters
We consider two setups for the studies of spherical swimmers, with different packing fractions and propulsion velocities.Common to both cases is the number of swimmers, N sw = 196, the radius of the squirmer, R = 3a, the with of the slit, L y = 7a, and the active-stress parameter, β ∈ {−1, 0, 1}.
Set 1 -High packing fraction, small Péclet number: The choice of the box length L = 96a = 32R corresponds to a 3D packing fraction of φ = 4πN sw R 3 /3L 2 L y = 0.34, or the quasi-2D packing fraction We employ the time step h = 0.02 ma 2 /(k B T ) and the mean number of particles in a cell N c = 80, which yields a fluid viscosity of η = 178 mk B T /a 4 , as determined by independent simulations [58].The resulting Péclet number, which compares the time scale for rotational diffusion, 2D R (D R is the rotational diffusion coefficient), to the swimming time scale, 2R/U 0 , is ABP system exhibits strong clustering and hexagonal order, whereas no large clusters and no pronounced order emerges for neutral squirmers.Hence, hydrodynamics strongly affects the aggregation behavior.
In order to characterize the clustering behavior more quantitatively, we analyse the density distribution by performing a Voronoi decomposition of the fluid volume with the squirmer's centers as generator points [24,60].A local packing fraction φ loc of squirmers is then introduced as the ratio of a squirmer volume to the volume of its Voronoi cell.Note that the Voroni construction yields in general a distribution function which is different from the distribution function of the local density.In dilute regions, φ loc is low, since the Voronoi volume is large, and P (φ loc ) is low too, because of the low particle density.Figure 4(a) shows distributions of local packing fractions for squirmers and ABPs at P e = 115.We consider two system sizes for the ABPs and pullers (β = 1), namely L = 96a and L = 192a, to provide clear evidence of a possible MIPS.The ABP system exhibits pronounced peaks at the packing fractions φ 2D loc ≈ 0.8 (L/a = 96, solid) and φ 2D loc ≈ 0.85 (L/a = 192, dashed), respectively, close to the maximum value of φ 2D = 0.907 of 2D hexagonal packing, and lower probabilities for dilute regions.The two curves reflect a strong system-size dependence, as is well-known for first-order phase transitions, because the system is affected by interfaces [61].Accordingly, and in agreement with previous simulations [21][22][23], our APB system phase separates into a dense Cluster-size distribution function for squirmers (pusher: β = −1, puller: β = 1, neutral: β = 0) and ABPs for the packing fraction φ 2D = 0.6 and the Péclet number P e = 575.The solid and dashed lines correspond to the system sizes L = 96a and L = 192a, respectively.The lines (light blue) indicate the power-law decay P (n) ∼ n −δ , with δ = 1 (pusher: β = −1), 1 (neutral: β = 0), 1.4 (puller: β = 1), and δ = 2, 2.1 for the small and large ABP system, respectively.
giant cluster comprising the vast majority of the ABPs (cf.Fig. 5) and a dilute gas-like region.In contrast, squirmers show a far less pronounced high-density region.Most remarkable, the distribution function P (φ loc ) is independent of the system size, as indicated in Fig. 4(a) for pullers, which we consider as indication that squirmers exhibit no MIPS.The force dipole evidently strongly affects cluster formation as is reflected by the differences of the distribution functions.Thereby, pulling (β = 1) promotes formation of larger clusters, whereas pushing (β = −1)is disadvantage for cluster formation.At the lower packing fraction φ 2D = 0.4 and the high Péclet number P e = 575 (Fig. 4(b)), ABPs exhibit a two-peak distribution, which becomes more pronounced with increasing system size, and peaks shift to smaller and larger concentrations, respectively, as expected for a first order phase transition, indicating a MIPS in the limit L → ∞.For the squirmers, the maximum of the local density for pullers and neutral squirmers coincides with the global density φ 2D = 0.4.However, for pullers a long tail toward higher local densities indicates formation of large temporary clusters (MIC).Considering the higher packing fraction φ 2D = 0.6 for P e = 575, we find little differences to the density distributions of Fig. 4(a).
Another quantity to characterize the structure is the cluster-size distribution function N (n), which is defined as where N (n) is the average fraction of particles that are members of a cluster of size n, and p(n) is the number of clusters of size n.The distribution is normalized We use a distance criterion to define a cluster: a swimmer belongs to a cluster, when its center-to-center distance to another swimmer of the cluster is within 1.02σ.As shown in Fig. 5, predominately smaller clusters are formed for squirmers.However, ABPs exhibit a peak for clusters comprising almost all swimmers, consisted with MIPS.This is particularly evident for the larger system size.With increasing system size, the probability for intermediate-size clusters drops significantly and N (n) is dominated by small and giant clusters.However, the initial power-law decay of N (n) changes only slightly from N (n) ABP ∼ n −2.0 for L = 96a to N (n) ABP ∼ n −2.1 for L = 192a.For larger systems, we expect an even more pronounced giant cluster consistent with the appearance of MIPS.Compared to ABPs, squirmers exhibit a higher probability for medium-size clusters in the range 2 < n 150.The cluster-size distributions decay as a power-law for n 30.Thereby, pusher exhibit the slowest decay with N (n) ∼ n −1.0 , and pullers the somewhat fastest decay N (n) ∼ n −1.4 .This power-law decrease of N (n) is independent of the system size.The qualitative and quantitative behavior changes when the concentration is reduced (not shown).For φ 2D = 0.4, all systems exhibit an initial power-law decay N (n) ∼ n −1 and no system-spanning cluster appears.However, pushers and neutral squirmers show a stronger, non-power-law decay for large cluster sizes.Hence, there is a crossover from a power-law to a faster decay of the cluster-size distribution function with decreasing concentration for β 0. A similar behavior has been observed in Ref. [36].
A pronounced long-range translational order for ABPs is clearly visible in the pair-correlation function g(r) [62,63] presented in Fig. 6.Even on the length scale of half of the system size L, peaks are visible for ABPs.Squirmers exhibit a less pronounced translation order for φ 2D = 0.6.Thereby, as already indicated in the local packing fraction, pullers (β = 1) exhibit the most pronounced order, but the order decays faster with distance than for ABPs.Pushers (β = −1) exhibit the lowest tendency to long-rang order, and g(r) assumes the idealgas value essentially on the scale of three to four neighbor distances.A decreasing packing fraction leads to decreasing long-range order.For φ 2D = 0.5, APBs still show the most pronounced order, however, the correlation function reaches unity for r/σ 7.At φ 2D = 0.4, only short-range order is present for r/σ 4. Interestingly, pullers exhibit the most pronounced order in this case.
The peak positions of Fig. 6 are consistent with a hexagonal packing of the swimmers.The preference of such an order is confirmed by the hexagonal order parameter |q 6 | 2 defined as [23,34,64] where N k 6 is the set of the six nearest neighbors of swimmer k, and ϑ kj is the angle between the vector R k − R j connecting the centers of particles j and k and the xaxis.For a perfect hexagonal lattice |q 6 | 2 = 1. Figure 7 displays probability distribution functions of the order parameter for squirmers and ABPs at P e = 115.Distributions for the larger Péclet number P e = 575 closely agree with those of Fig. 7. Consistent with the high local packing, ABPs exhibit an pronounced probability for hexagonal order.The respective values for squirmers are smaller.Thereby, the tendency to form hexagonal order is most distinct for pullers.In contrast, pushers exhibit hardly any local hexagonal order.
Strictly 2D ABP systems exhibit stronger aggregation and a more pronounced phase separation [23,34] compared to similar systems in quasi-2D [34].The reason is that the freedom of the propulsion vector e to indepen- dently change in a diffusive manner in 3D or quasi-2D systems reduces the lateral swim pressure when e is oriented normal to the confining wall, i.e., the driving force for cluster formation is reduced.Interestingly, in our systems the swim direction of the squirmers is preferentially parallel to the surfaces, as shown in Fig. 8. Hence, despite a preferential 2D swimming motion, squirmers are less likely to form clusters in the quasi-2D geometry compared to the strict 2D case.Moderate fluctuations in the propulsion direction seem to be sufficient to modify the onset of cluster formation.
The suppression of MIC and MIPS is attributed to a strong reorientation (scattering) of squirmers by hydrodynamic torques during their collisions [27].This conclusion is quantitatively confirmed by our simulations in Fig. 9, where the orientation autocorrelation function for squirmers and ABPs is shown.It decays exponentially as exp(−t/τ sq ) with a characteristic time τ sq .For ABPs, the relaxation time is given by 1/2D R , in agreement with theoretical expectations.However, for squirmers, τ sq is approximately an order of magnitude smaller, and rather similar for the various force-dipole strengths.This can be interpreted as a reduced effective Péclet number P e = τ sq U 0 /R ≈ 10, a value, where no phase transition is expected even for ABPs.However, the differences in the tendency of the squirmers to form clusters (cf.Fig. 4(a)) indicate the importance of hydrodynamic interactions, which are not captured in the Péclet number calculate on the basis of the rotational diffusion coefficient at infinite dilution.
The relaxation time τ sq depends on the squirmer concentration, as shown in Fig. 10(a) for P e = 575.Compared to ABPs, τ sq of squirmers is reduced by orders of magnitude (factor 50 − 80) and decreases almost linearly in the considered range of concentrations (for P e = 115, 2D R τ sq ≈ 0.1 (cf.Fig. 9).)Naturally, τ sq will increase in a nonlinear manner for φ 2D → 0, because at φ 2D = 0, 2D R τ sq = 1.In addition, the average swimming velocity |U | decreases with increasing concentration, also for ABPs (cf.Fig. 10(b)).Thereby, the velocity for ABPs seems to decrease linearly, whereas those of the various squirmers decline nonlinearly.A decreasing swimming velocity and relaxation time with increasing areal fraction have also been found in Ref. [27] for strictly 2D systems.

V. SIMULATIONS-SPHEROIDAL SQUIRMERS A. System Setup and Parameters
We consider spheroidal squirmers with the minor axis b x = 3a and the aspect ratios b z /b x ∈ {1, 1.5, 2, 3, 4} in the quasi-2D geometry described above (see Fig. 2).We employ the time step h = 0.02 ma 2 /(k B T ), rotation angle α = 130 • , and the mean number of particles N c = 10, which yields the fluid viscosity η = 17.8 mk B T /a 4 .
To investigate the collective behavior for different packing fractions φ = 4πN sw b 2 x b z /3L 2 L y , or equivalently area packing fractions φ 2D = N sw πb x b z /L 2 , the length L of the simulation box is varied, while the number of swimmers is fixed at N sw = 196 for b z /b x = 1 and N sw = 200 for all other aspect ratios.We choose the values of the swimming mode B 1 / k B T /m = 0.01, 0.007, 0.005, 0.003, and 0.002, which yields, with Eq. ( 3), the identical Péclet number for all aspect ratios.Here, D ⊥ R = k B T /ξ ⊥ is the rotational diffusion coefficient around the minor axis.Note that this Péclet number is much smaller than those of the previous section for our studies on MIPS of spheres.As before, β ∈ {−1, 0, 1} and simulations of ABPs are performed for comparison.

B. Results
Figure 11 illustrates structure formation for various aspect ratios.Although the Péclet number P e = 12 is quite small, spheroidal squirmers clearly reveal clusters for all aspect ratios b x /b z 2. Thereby, a larger aspect ratio is beneficial for cluster formation.For spherical squirmers, no phase separation occurs at this Péclet number (Fig. 11 (a)).The anisotropic nature of the spheroids leads to shape-induced jamming and alignment, where alignment increases with increasing aspect ratio as is evident from the comparison of Figs.11(b) and 11(d).
Again, we employ Voronoi tessellation to quantify clustering and phase separation.However, we cannot use spheroid centers as generator points.For spheroids, the distance from a point in space to a swimmer is mea-sured as the Euclidian distance to the nearest point on the swimmer's surface.For an efficient calculation of the (approximate) Voronoi volume, we employ a vertex-mesh on the surface of every spheroid [65].Those vertices serve as generator points and the Voronoi cell of a swimmer is defined as the union of the cells of its vertices.
The probability distribution P (φ loc ) of the local packing fraction is displayed in Fig. 12.No phase separation occurs for ABPs and pushers (β = −1) at the concentration φ 2D loc = 0.3 (Fig. 12(a)); however, clusters appear for neutral squirmers (β = 0) and pullers (β = 1).At the higher concentration φ 2D loc = 0.5 (Fig. 12(b)), all squirmers with β = 0, ±1 exhibit dense clusters.In contrast, spheroidal APBs exhibit only a weak tendency to form clusters.This points toward a major role of hydrodynamics in cluster formation of elongated squirmers.Simulations of pushers with a strong active stress (β = −5) emphasized this aspect.Here, we find a clear maximum in P (φ loc ) at the average density and, hence, no MIPS.
By varying the aspect ratio (b z /b x ∈ {1, 1.5, 2, 3}) and the squirmer concentration (φ 2D ∈ {0.1, 0.2, 0.3, 0.4, 0.5}), we establish the state diagram displayed in Fig. 13.The various colored areas indicate the coexistence of dense cluster region and a dilute region.The respective phase borders should be viewed as approximations to the infinite-systems-size phase diagram.An accurate calculation of the phase boundaries requires a substantial simulation effort due to finite-size effects.The considered system sizes are certainly too small to ensure absence of finite-size effects.The blue area indicates that pullers (β = 1) already phase separate for relatively small aspect ratios and at low packing fractions.For decreasing β (black area β = 0, red area β = −1), the phase-separation line shifts to the top right, i.e, higher aspect ratios and packing fractions are required for phase separation to occur.The fact that the area for phase separation of ABPs is farthest to the top-right indicates that hydrodynamics enhances phase separation for elongated squirmers.Hence, the influence of hydrodynamic interactions is reversed compared to spherical swimmers, where hydrodynamics suppresses cluster formation.However, the qualitative effect of the active stress remains unchanged-pulling enhances cluster formation compared to pushing.

VI. DISCUSSION
Our simulation results for spherical squirmers agree with certain aspects of some previous studies, but disagree with others.In particular, the cluster-size distribution function has been analyzed in Ref. [36] for the same squirmer model with the Lattice Boltzmann approach, but a somewhat different simulation setup.Significantly larger systems with many more squirmers have been considered, and three-dimensional hydrodynamic interactions, but squirmers confined in 2D at the density The blue, black, red, and purple areas indicate giant clusters for pullers (β = 1), neutrals (β = 0), pushers (β = 1), and active Brownian particles, respectively.Thereby clusters of pullers appear in all colored areas, clusters of neutral squirmers in the black, red, and purple area, etc. See Supplemental Material at ??? for individual state diagrams of the various swimmers.
φ 2D = 0.1.We like to emphasize that our studies yield probability distribution functions for the local packing fraction of squirmers, which are independent of the system size.In agreement with our Fig. 5, in Ref. [36] a power-law decay (p(n) ∼ n −2.0 ) for small cluster sizes (n 10 3 ) is obtained for pullers with 0 < β 1.For other β (β > 2 and β < 0), a stronger decay appears for n 10 in [36].This non-power-law decay is stronger than that observed in our simulations.Hence, at smaller concentrations we expect a stronger drop of the cluster distribution for β < 0 also in our quasi-2D systems.
Our simulation results of spherical squirmers disagree qualitatively and quantitatively with results of Refs.[34,47], where comparable systems have been considered.We find that hydrodynamic interactions suppress cluster formation of squirmers, in contrast to Ref. [34].According to Ref. [47], the Péclet numbers of our neutral squirmers are deep in the phase-separated region of the phase diagram (Fig. 3 of Ref. [47]) [66].However, we do not see any signature for phase separation, not even for neutral squirmer.Moreover, the distribution function of the bond order parameter differs significantly.For squirmers, we find the most pronounced local order for pullers, followed by neutral squirmers and weakest order for pushers, in agreement with other studies [35][36][37].The studies of Ref. [34] show strongest cluster formation and local order for neutral squirmers, comparable to strictly 2D ABPs and more pronounced compared to ABPs in quasi-2D.
We suspect that the origin of the discrepancy between our results and those of Refs.[34,47] is the different compressibilities of the employed MPC fluids.It is shown in App.C that the squirmer-induced fluid transport can cause strong MPC-fluid density inhomogeneities when squirmers are blocked, either by other squirmers or by a surface, at too strong propulsion velocities.A significant fluid-density modulation leads to a reduced propulsion efficiency and a reduced active pressure inside the cluster.As a consequence, propulsion of squirmers confined inside a cluster is effectively reduced, which leads to more stable clusters.By our choice of simulation parameters, we ensured that critical fluid density variations are avoided.Thus, the results of the simulation studies of Refs.[34,47] are incompatible with ours; however, if (strong) compressibility effects would be considered as relevant, like fluids near a critical point, these studies are applicable.
In the case of spheroidal squirmers, two aspects contribute to the hydrodynamic enhancement of cluster formation, namely near-field HIs between squirmers, and between a squirmer and the confining walls [43].The effect of the active stress β is evident from the state diagram of Fig. 13.Note that the force (FD) and source dipole (SD) flow-field contributions of a squirmer [43] with e y = 0 decay as −3 and −2 in the center between the walls, respectively, with the distance from a squirmer center parallel to the walls [48].Our simulation studies of Ref. [43] on the cooperative motion of two squirmers in a slit geometry already reveal a long-time stable configuration of two pullers, in which they swim together in a wedge-like conformation with a constant small angle.A similar configuration of neutral squirmers or pushers is far less stable and the squirmers scatter more strongly, which emphasizes the relevance of specificities in fluid interactions of the squirmers due to their close proximity, and particular with the walls.Applying three-dimensional hydrodynamic interactions rather than no-slip conditions to a still confined pair of pullers yields short cooperative motion only, emphasizing that surface hydrodynamics is extremely important [43].Hence, the seed for the cluster formation is the cooperative motion of two squirmers.Blockage by additional encountering squirmers further stabilizes the emerging cluster.The latter is reflected in the density dependence of the cluster formation.

VII. SUMMARY AND CONCLUSIONS
We have studied the collective dynamics of spherical and spheroidal squirmers confined in a narrow slit by mesoscale hydrodynamic simulations (MPC).To elucidate the role of hydrodynamic interactions, we have complemented the studies by simulations of spherical and spheroidal active Brownian particles.
We find that hydrodynamics suppresses phase separation for spherical microswimmers, in contrast to Ref. [34], but in accordance with Ref. [27].We attribute the contradiction with Ref. [34] to compressibility effects of the MPC fluid.In agreement with Ref. [27], we confirm that a hydrodynamically enhanced reorientation of squirmers during swimmer-swimmer collisions implies an effectively smaller Péclet number and, hence, a suppression of phase separation.However, the squirmers form shortlived clusters over a broad range of cluster sizes.As already pointed out in Ref. [37], near-field hydrodynamic interactions play a major role in cluster formation.This is reflected in the clear differences between pullers, neutral squirmers, and pushers [37].A strong near-field effect can be expected, since the multipole contributions to a swimmer's flow field decay quickly with distance parallel to the confining walls in a narrow channel [48].
For anisotropic swimmers, alignment due to steric interactions is a key mechanism for MIC, while isotropic swimmers form clusters due to jamming and blocking.Our studies of elongated, prolate spheroidal squirmers yield hydrodynamically enhanced MIC compared to ABPs.This result is surprising at first glance.Based on our studies on the cooperative motion of pairs of confined elongated squirmers, we attribute this enhancement to near-field hydrodynamics and, most importantly, swimmer-surface hydrodynamic interactions.Yet, the anisotropic shape leads to shape-induced phase separation for ABPs at considerably smaller Péclet numbers than for spheres.This is related to differences in the clustering mechanisms of isotropic and anisotropic swimmers.
For both, spherical and spheroidal squirmers, the active stress parameter β correlates positively with cluster formation, i.e., pullers (β > 0) show the strongest tendency to form clusters, followed by neutral squirmer (β = 0), and pushers (β < 0).In fact, this is well consistent with our previous study of two squirmers in a narrow slit, which shows that cooperativity is most pronounced for pullers [43].
Simulations of dumbbell swimmers [67,68] in three dimensions have also been performed and hydrodynamically enhanced MIC has been observed [69].Since dumbbell swimmers are anisotropic, this finding is consitent with our results for spheroidal squirmers.However, we have to be careful with too general conclusions.As we demonstrated, MIC depends decisively on the hydrody-namic flow field, i.e., the details of the flow field matters.This is reflected by our results for very strong pushers (β = −5), which show a suppression of phase separation.
The observed clustering of spheroidal squirmers is in contrast to the behavior of motile bacteria, e.g., E. coli, in suspensions, which display active turbulence [12].At the moment, there is no detailed understanding of the different behaviors and underlying mechanisms.Could a dependence on the Péclet number explain the differences?The rotational diffusion coefficient D e R = 0.057/s of non-tumbling E. coli has been measured in Ref. [12].With the swimming velocity U e 0 = 30µ/s and the length 2b z ≈ 10µ, we obtain a Péclet number P e e ≈ 50.Our simulations show that such P e even further enhance clustering.There seems to be an additional relevant interaction in bacteria suspensions.Of course, it could be related to the structure of bacteria with the cell body and the flagella bundle [70], which form a rather flexible structure.Depending on the body-bundle orientation, a bacterium exhibits a wobbling motion, which perturbs alignment and leads to a fast decay of orientational correlations of nearby swimmers [71].Furthermore, hydrodynamic interactions could be responsible-the propulsion of the bacterium with a rotating flagella bundle and a counterrotating cell body leads to a rotlet dipole [72].Such a flow field leads to circular motion on a planar surface [72][73][74][75] and may also enhance the decorrelation of the orientational order.
Squirmers are highly idealized models of biological microswimmers, not only in terms of shape, but also in terms of the stationary, time-independent flow fields.The strong interplay between shape and hydrodynamics, which is revealed by our simulations, indicates that the details of the cell shape and beat pattern may play a much more important role in the collective dynamics of biological microswimmers than previously expected.
In MPC [44][45][46], a fluid is represented by N point particles of mass m with continuous positions r i and velocities v i (i ∈ {1, . . ., N }).The particle dynamics proceeds in two steps-streaming and collision.During the streaming step of duration h (collision time), the particles move ballistically, i.e., their positions are updated according to (A1) In the collision step, the particles interact locally by an instantaneous stochastic process, for which we apply the stochastic rotation dynamics (SRD) variant of the MPC method [44][45][46]76] with angular momentum conservation (MPC-SRD+a) [58,77,78].For this purpose, the simulation volume is partitioned into cubic collision cells of length a.In each of the cells the particle velocities are then updated as [58,77] Here, R(α) is the rotation matrix for the rotation of the relative velocity v i,c = v i − v cm , with respect to centerof-mass velocity v cm of the cell, by an angle α around a randomly oriented axis.r i,c = r i − r cm is the position relative to the center of mass of the cell r cm , and I is the moment of inertia tensor of the particles in the cell's center-of-mass frame.The grid of collision cells is shifted randomly before each collision step to ensure Galilean invariance [79].Thermal fluctuations are intrinsic to the MPC method [58,80,81].Hence, a cell-level canonical thermostat (MBS thermostat) is applied after every collision step, which maintains the temperature T [80,81].The algorithm conserves mass, linear, and angular momentum on the collision cell level, which gives rise to hydrodynamics on large length and long time scales [45,82,83].We measure lengths in units of the collision cell size a and time units of ma 2 /k B T , where k B is the Boltzmann factor, which corresponds to the choice a = m = k B T = 1.The rotation angle is set to α = 130 • .The viscosity η of the fluid can be tuned by changing the mean number of particles in a cell N c or the time step h [58,77,81].

Boundary Conditions and Squirmer Implementation
No-slip boundary conditions are implemented at the confining walls (cf.Fig. 2) via the bounce-back rule v i (t) = −v i (t), where v i and v i are the velocity before and after the particle's collision with the wall, respectively.Additionally, spatially uniformly distributed phantom particles with Gaussian distributed velocity components are inserted into the walls and interact with fluid particles during the collision step [46,84].Parallel to the walls, periodic boundary conditions are applied.
During the streaming step, a squirmer moves like a passive colloid according to rigid-body equations of motion.The equations of motion of the rotational degrees of freedom are solve by introducing quaternions [43,63] (see also App.A 3). Interacting fluid particles experience the modified bounced back rule v i = −v i + 2v S , where v S is the squirmer's surface velocity at the point of contact.This surface velocity includes its translational and rotational velocity as well as its slip velocity u sq (1).The change of a fluid particle's linear and angular momentum is transferred to the squirmer, such that the total momenta are conserved.For the collision step, as for walls, spatially uniformly distributed phantom particles are generated inside the squirmer with Gaussian distributed velocity components, where their mean value is given by the rigid-body translational and rotational velocity plus the squirming velocity, and the variance is equal to k B T /m [42,43].After the collision step, the linear and angular momentum changes of all phantom particles are transferred to the squirmer.A more detailed description is provided in Ref. [43].Repulsive forces and torques due to steric interactions between squirmers, and between a squirmer and a wall, during a streaming step are implemented as outlined in the appendix of Ref. [43] (see also Ref. [85]).

Simulation Algorithm for Spheroidal Active Brownian Particles
For active Brownian particles, as for the hydrodynamic case, we solve the rigid-body equations of motion, but now in the presence of translational and rotational noise.Thereby, the squirmer orientation is parameterized by a unit quaternion q = (q 0 , q 1 , q 2 , q 3 ) [63].The integration scheme for the overdamped dynamics of the center-ofmass R(t) and the quaternion q(t) of the spheroid is [86] where h is the time step, F = F st + F th and T = T st + T th are the force and torque on the ABP due to steric interactions and thermal noise.The 4×4 matrix Q comprised of the quaternions is defined in Appendix B. As before, U 0 is the propulsion velocity and e ≡ D T e z points along the major axis of the squirmer and in the direction of the desired propulsion.The rotation matrix D that rotates vectors from the laboratory frame into the body-fixed frame is given in Appendix B. The Lagrangian multiplier λ q is introduced to ensure normalization of q.Hence, we perform the quaternion update without the term λ q q(t) to obtain some q and then determine λ q such that q + λ q q is normalized [86].
The translational and rotational friction coefficients of a spheroid are given by γ = 6πηb z   MPC is a compressible fluid and correspondingly strong density inhomogeneities can emerge for certain choices of MPC parameters.Figure 14 provides an example of the fluid-particle density distribution for neutral spherical squirmers confined in a narrow slit (cf.Fig. 2).The squirmer density is φ 2D = 0.6 and the Péclet number P e = 115.However, the Péclet number is realized by different combinations of MPC fluid and squirmer parameters.In Fig. 14 Figure 14(a) shows large density inhomogeneities with dense (orange) and rarefied regions (purple).Thereby, the rarefied regions coincides with the area occupied by the large squirmer cluster (cf.inset of Fig. 14(a)).In contrast, the MPC fluid distribution in Fig. 14(b) is rather homogeneous aside from thermal fluctuations and no link to the squirmer distribution is evident.It seems that a pressure gradient is build up by active clustered squirmers by expelling fluid form the cluster-an effect related to the compressibility of MPC.For density inhomogeneities as large as those in Fig. 14(a), MPC simulations do not describe incompressible fluids, since the viscosity of the MPC fluid depends linearly on density and would, hence, depend strongly on position.
The mechanism in a MPC fluid, which opposes the diminished density in the gaps between squirmers is fluidparticle diffusion.The ratio of the time t dif f necessary for a MPC particle to diffuse over a swimmer diameter σ = 2R, compared to the time required for a MPC particle to be advected by the surface activity of a squirmer with the velocity U 0 over the same distance, yields the dimensionless pumping number where D f = 0.013 √ a 2 k B T m is the MPC-fluid diffusion coefficient.For small pumping number, P u, we expect the density inhomogeneities to disappear, i.e., MPC fluidparticle diffusion is faster than activity-induced advection.The above choices of parameters yield the pumping numbers P u = 5 (Fig. 14(a)) and P u = 0.5 (Fig. 14(b)).Consistent with our expectation, the difference in the pumping number explains the appearing density modulations.
In terms of the Péclet number, P u is given by P u = σ 2 D R P e/6D f .Since P e 1 for active systems, P u < 1 requires σ 2 D R /6D f 1.Hence, the MPC parameters have to be chosen such that D R 1.Because D R ∼ 1/η, this is achieved for N c 1.Such an increase in N c leaves D f virtually unchanged, however, implies an increase in computational effort.

FIG. 1 .
FIG. 1.Illustration of normal and tangent vectors of a spheroidal (left) and spherical (right) squirmer.Selfpropulsion along the body-fixed orientation vector e, here e = (0, 0, 1) T , is achieved by an axisymmetric prescribed surface velocity in the direction of the tangent vector s[43].
Figure 3 displays snapshots of structures of spherical ABPs and squirmers for the parameter Set 1.The

FIG. 8 .FIG. 9 .
FIG.8.Probability distribution functions of the propulsiondirection component ey normal to the confining walls for ABPs and squirmers.The squirmers are preferentially orientation parallel to the walls.The Péclet number is P e = 575 and the concentration φ 2D = 0.6.The distribution functions for P e = 115 and the same density are nearly identical.

FIG. 10 .
FIG. 10.(a) Decay time τsq of the orientational correlation function Ce(t) and (b) average swimming velocity as function of the squirmer concentration for neutral squirmers (β = 0), pushers (β = −1), pullers (β = 1), and active Brownian particles.The decay time is scaled by the rotational relaxation time 1/2DR and the swimming velocity by that of the set value U0.The Péclet number is P e = 575.The lines are guides to the eye (either linear or parabolic functions).

2 −
FIG. 14. Fluid particle density ρ in a slit with spherical neutral squirmers at P e = 115 and φ 2D = 0.6.The total fluid density is show, comprising the MPC particle density ρ f l and that of the phantom particles ρ ph , i.e., ρ = m Nc /a 3 = ρ f l + ρ ph .ρ0 denotes the average density.In (a), the choice B1 = 0.1 kBT /m and Nc = 10 yields the pumping number P u = 5, and in (b), B1 = 0.01 kBT /m and Nc = 80 yields P u = 0.5; the other parameters are the same.The insets show the corresponding squirmer locations.Please, note the differences in the color bars in (a) and (b).
Appendix C: Avoiding MPC Simulation Artifacts Related to Density Inhomogeneities (a), we set B 1 = 0.1 k B T /m and N c = 10, and in (b) B 1 = 0.01 k B T /m and N c = 80.The latter gives a ten times larger viscosity.The other parameters are the same as defined in App.A 1 and IV A.