Clustering of microswimmers: interplay of shape and hydrodynamics

Abstract

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 and hydrodynamics. These direct interactions are typically anisotropic, and come in different “flavors”, such as spherical and elongated particle shapes, pusher and puller flow fields, 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. For spherical particles, the comparison between ABPs and hydrodynamic-squirmer ensembles reveals a suppression of MIPS due to hydrodynamic interactions. Yet, our analysis shows that dynamic clusters exist, with a broad size distribution. 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. The transition to a phase-separated state strongly depends on the nature of the swimmer's flow field—with an increased tendency toward MIPS for pullers, and a reduced tendency for pushers. Thus, hydrodynamic interactions show opposing effects on MIPS for spherical and elongated microswimmers, and details of the propulsion mechanism of biological microswimmers may be very important to determine their collective behavior.

---

1 Introduction

Motile bacteria at interfaces exhibit intriguing collective phenomena,^1,2 such as cluster formation, observed for Myxococcus xanthus³ or Thiovulum majus,⁴ as well as swarming, swirling, raft formation,^5–10 and the emergence of mesoscale turbulence,^8,11,12 observed for E. coli¹³ or Bacillus subtilis.¹⁴ Similarly, experiments on self-phoretic artificial spherical microswimmers, such as Janus particles, self-propelled liquid droplets, and photo-activated colloids, exhibit cluster formation and phase separation despite their isotropic shape and purely repulsive interactions.^15–21 The various active agents are propelled by different mechanisms and may exhibit different steric and propulsion-related interactions. Hence, it is a priori not evident which processes govern structure formation in the rather distinct systems. Unravelling the underlying universal features and discriminating them from specific aspects requires dedicated experimental and theoretical studies.

Computer simulations of model systems, such as spherical active Brownian particles (ABPs), yield motility-induced clustering (MIC) and motility-induced phase separation (MIPS)^20,22–27 in qualitative agreement with the above-mentioned experiments on synthetic self-phoretic 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.^21,28,29 However, in self-phoretic systems, the phoretic field has been shown to significantly contribute to clustering through field-induced attractive interactions.^30,31

Cluster formation, swarming, and mesoscale turbulence of anisotropic objects such as bacteria are 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,32–34 and even leads 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, a crucial aspect regarding MIC and MIPS of microswimmers is the role of hydrodynamics. At equilibrium, hydrodynamic interactions (HI) solely affect the dynamics and not the structural properties of a system. Since active systems are intrinsically out of equilibrium, this does not necessarily apply, and the steady-state—and hence the formation of clusters—is strongly affected by hydrodynamics.^28,35–38 Furthermore, the particular geometry of a microswimmer, such as the cell body and the flagellum or the flagella bundle of bacteria, may lead to a coupling of the orientational dynamics and fluid flow; such a coupling can lead to mesoscale turbulence.^39,40

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 are difficult to switch off, and therefore their contribution is hard to assess. In contrast, simulations can be performed with and without HI and the latter effect on MIPS can be analyzed, although with an increased numerical effort.

A frequently applied model for a self-propelled particle in the presence of a fluid is a squirmer^35,41–46—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).⁴⁶

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 (MPC) dynamics approach, a particle-based hydrodynamic simulation method to solve the Navier–Stokes equations.^47–49 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 ref. 35 and 28, which predict enhanced clustering of spherical squirmers or no MIPS in the presence of HI, respectively. On the other hand, we consider squirmers of prolate spheroidal shape. This is motivated by the different mechanisms 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 as compared to swimmers with elongated shapes.

As an important result, we find no MIPS for spherical squirmers in the quasi-2D system, which is in agreement with ref. 28 and 38, but in contrast to ref. 35. We attribute the discrepancy to peculiarities of the compressible MPC fluid employed in ref. 35, which requires a suitable choice of the parameters for simulations of low-compressibility fluids. Hence, HI 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 ratios, where we find a substantial enhancement of MIPS through HI. This applies for pushers, pullers, and neutral squirmers, as long as the force dipole is sufficiently weak. In the case of large force dipoles, pushers show the least tendency to phase separate, 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, and sometimes antagonistically.

The paper is structured as follows. Previous results of squirmer simulations of collective phenomena are briefly summarized in Section 2. Section 3 introduces the model for prolate squirmers. A comparison between MPC and ABP simulation results of spherical squirmers is presented in Section 4. Section 5 presents results for prolate spheroidal swimmers, and Section 6 discusses the various obtained aspects. Section 7 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), and the simulation method for ABPs (A.3). Appendix B collects definitions relevant to rigid-body dynamics of ABPs, and Appendix C discusses technical aspects of fluid compressibility in particle-based simulations of squirmers.

2 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 ref. 28 and 38, and ref. 35 and 50, respectively. In ref. 28, 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 HI, to highly flattened cylinders with 3D HI. Independent of the cylinder length (and dimensionality of HI), 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 HI compared to the configuration-independent rotational diffusion of ABPs. The studies of ref. 38 emphasize the importance of near-field HI in the suppression of MIPS. However, clustering of the squirmers is observed as well as the formation of polar order for nearly neutral squirmers. Note that in both ref. 28 and 38, thermal fluctuations are neglected. Hence, there is no thermal rotational diffusion leading to a decorrelation of the squirmer's orientational motion. In the simulations of ref. 35 and 50, squirmers are confined in a narrow slit (quasi-2D), slightly wider than 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/ρ² parallel to the confining surfaces in the slit center, where ρ is the radial distance from the squirmer.⁵¹ Most importantly, hydrodynamics is argued to enhance MIPS, in contrast to findings of ref. 28 and 38, and a phase diagram is presented in ref. 50.

A quasi-2D setup is certainly less favorable for cluster formation of ABPs compared to a strict 2D system,^24,52 since the swimmer's propulsion direction can point toward a wall, which reduces steric blocking and the lateral swim pressure.^53–55 Thus, compared to strictly 2D ABP systems, ABPs in quasi-2D require higher 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 HI predict a parallel alignment of the propulsion direction for pushers and a normal orientation for pullers at a no-slip surface.⁵⁶ 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 ref. 35 and 46.

For suspensions of athermal spherical squirmers in 3D, lattice-Boltzmann⁵⁷ and force-coupling³⁶ 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 large-scale MIPS is observed. In contrast, ABPs in 3D clearly exhibit MIPS^25–27 and large-scale collective motion.²⁵ Moreover, simulations of 2D, quasi-2D, and 3D systems of attractive squirmers show a substantially enhanced cluster formation compared to purely repulsive squirmers.^36,37 Also in such situations, hydrodynamics is found to reduce MIPS compared to ABPs.

So far, the influence of HI on MIPS has mainly been studied for spherical squirmers. In contrast, very little is known about the influence of HI on the structure and dynamics of elongated spheroidal squirmers confined in narrow slits. Here, the shape, squirmer–squirmer HI, and HI with the confining surfaces affect cluster formation and a possible MIC and MIPS.⁴⁶

3 Model and simulation approach

We model a nonspherical squirmer as a prolate spheroidal rigid body with the prescribed surface velocity (cf.Fig. 1)⁴⁶

u_sq = −B₁(e_ζ·e_z)(1 + βζ)e_ζ, (1)

which is an extension of and alternative to earlier approaches.^58–60 We employ spheroidal coordinates (ζ,τ,φ) (cf.Fig. 1), which are related to Cartesian coordinates viawhere , −1 ≤ ζ ≤ 1, 1 ≤ τ < ∞, and 0 ≤ φ ≤ 2π. The advantage of the choice (1) for the slip velocity is that an analytical solution of the flow field can be obtained.⁴⁶ Constant B₁ in eqn (1) determines the amplitude of the slip velocity and is related to the self-propulsion velocity U₀via⁴⁶

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

where τ₀ = b_z/c. The active stress is captured by coefficient β,^46,61 where, β < 0 corresponds to a pusher, β > 0 to a puller, and β = 0 to a neutral squirmer.

Fig. 1 Illustration of normal and tangent vectors of a spheroidal (left) and spherical (right) squirmer. Self-propulsion 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.⁴⁶

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 are 2b_x/L_y = 6/7, i.e., we focus on a quasi-2D geometry.

Fig. 2 Sketch of spheroidal microswimmers in a narrow slit. The ratio of the squirmer minor axis and the width of the slit is 2b_x/L_y = 6/7. The arrows indicate the swimming direction along a squirmer's major axis.

The squirmers are embedded in a fluid, which we simulate by multiparticle collision (MPC) dynamics, a mesoscale hydrodynamics simulation approach.^47–49 The details of the algorithm are described in Appendix A.1, and the implementation of the squirmer and the coupling with MPC in Appendix A.2.

For comparison, we perform simulations of active Brownian spheres and spheroids. The respective simulation approach is described in Appendix A.3.

4 Simulations—spherical squirmers

In the first step, we analyze the influence of HI on the phase behavior of spherical squirmers. The phase diagram for ABPs in two dimensions^23,26,52 shows a pronounced regime of MIPS above a certain Péclet number and a broadening range of densities with increasing Péclet numbers. Here, we consider two activities for ABPs and squirmers, respectively, and various densities within the phase separated regime. In order to provide a convincing and comprehensive characterization, we discuss various structural properties.

4.1 System setup and parameters

We consider two setups for the study 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 width of the slit, L_y = 7a, and the active-stress parameter, β ∈ {−1,0,1}.

Set 1 – Péclet number Pe = 115, packing fraction ϕ^2D = 0.6. The choice of the box length L = 96a = 32R corresponds to a 3D packing fraction of ϕ = 4πN_swR³/3L²L_y = 0.34, or the quasi-2D packing fractionWe employ the MPC time step and the mean number of particles in a cell 〈N_c〉 = 80, which yields a fluid viscosity of , as determined by independent simulations.⁶² 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₀, isfor the swim velocity U₀ = 2B₁/3, with , and from simulations. The latter is approximately 20% larger than the theoretical value. The chosen value of 〈N_c〉 = 80 is large compared to typical MPC studies, and results in an increased computational effort. However, it leads to a high viscosity, hence, a low D_R and a high Péclet number, a prerequisite to observe cluster formation in ABPs, as well as a low fluid compressibility. In ref. 63, a high Péclet number was achieved by a large B₁ and hence U₀. However, we find that such high values of B₁ lead to fluid-density inhomogeneities and consequently simulation artifacts (see Appendix C). Complementary to the squirmer simulations with MPC, we perform ABP simulations with the same values for U₀, D_R, R, and L_y, but consider additionally the larger system size L = 192a in order to reduce finite-size effects.

Set 2 – Péclet number Pe = 575, packing fractions ϕ^2D = 0.4, 0.5, and 0.6. The box lengths L/a = 96, 106, and 119 yield the packing fractions ϕ^2D = 0.6, 0.5, and 0.4, respectively. The MPC parameters and 〈N_c〉 = 480 result in a high viscosity and, with , the high Péclet number Pe = 575.

4.2 Results

Fig. 3 displays snapshots of structures of spherical ABPs and squirmers for the parameter Set 1. The ABP system exhibits strong clustering and a hexagonal order, whereas no large clusters and no pronounced order emerge for neutral squirmers. Hence, hydrodynamics strongly affects the aggregation behavior.

Fig. 3 Snapshots of spherical squirmers at Pe = 115, the packing fraction ϕ^2D = 0.6, and the simulation box length L = 192a. (a) Active Brownian particles exhibit MIPS and a local hexagonal order. (b) Neutral squirmers (β = 0) exhibit no long-range order and no MIPS. Note that the cluster in (a) partly dissolves in the course of time and a new system-spanning cluster is formed subsequently. Similarly, the clusters in (b) are unstable. See the ESI† for movies of the two cases (M1 and M2).

4.2.1 Density distribution. In order to characterize the clustering behavior more quantitatively, we analyze the spatial density distribution by performing a Voronoi decomposition of the fluid volume with the squirmers’ centers as generator points.^25,64 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 Voronoi 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.

Fig. 4(a) shows distributions of local packing fractions for squirmers and ABPs at Pe = 115 (parameter set 1). 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.⁶⁵ Accordingly, and in agreement with previous simulations,^22–24,52 our ABP system phase separates into a dense 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 remarkably, the distribution function P(ϕ_loc) is independent of the system size, as indicated in Fig. 4(a) for pullers, which we consider as an indication that squirmers exhibit no MIPS. The force dipole evidently strongly affects cluster formation as is reflected by the differences of the distribution functions. Therefore, pulling (β = 1) promotes the formation of larger clusters, whereas pushing (β = −1) is a disadvantage for cluster formation.

Fig. 4 Probability distribution of local packing fractions of spherical squirmers and ABPs (purple solid and dashed lines) for (a) the Péclet number Pe = 115 and the global area packing fraction ϕ^2D = 0.6, and (b) Pe = 575 and ϕ^2D = 0.4. Results for pushers (β = −1, red), pullers (β = 1, blue), and neutral squirmers (β = 0, black) are displayed. The solid and dashed blue and purple lines correspond to the system sizes L = 96a and L = 192a, respectively.

Fig. 5 Cluster-size distribution function for squirmers (pusher: β = −1, puller: β = 1, and neutral: β = 0) and ABPs for the packing fraction ϕ^2D = 0.6 and the Péclet number Pe = 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.

At the lower packing fraction ϕ^2D = 0.4 and the high Péclet number Pe = 575 (parameter set 2) (Fig. 4(b)), ABPs exhibit a two-peak distribution, which becomes more pronounced with an increasing system size, and the 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 the formation of large temporary clusters (MIC). Considering the higher packing fraction ϕ^2D = 0.6 for Pe = 575, we find little differences in the density distributions of Fig. 4(a).

4.2.2 Cluster size distribution. Another quantity to characterize the structure is the cluster-size distribution function [scr N, script letter N](n), which is defined aswhere [scr N, script letter 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 such that . 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, consistent with MIPS. This is particularly evident for the larger system size. With increasing system size, the probability for intermediate-size clusters drops significantly and [scr N, script letter N](n) is dominated by small and giant clusters. However, the initial power-law decay of [scr N, script letter N](n) changes only slightly from [scr N, script letter N](n)_ABP ∼ n^−2.0 for L = 96a to [scr N, script letter 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 in a power-law manner for n ≲ 30, where, pushers exhibit the slowest decay with [scr N, script letter N](n) ∼ n^−1.0, and pullers the somewhat fastest decay [scr N, script letter N](n) ∼ n^−1.4. This power-law decrease of [scr N, script letter 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 [scr N, script letter N](n) ∼ n⁻¹ 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. 37.

4.2.3 Pair correlation function. A pronounced long-range translational order for ABPs is clearly visible in the pair-correlation function g(r)^66,67 presented in Fig. 6 for ϕ^2D = 0.6. Even on the length scale of half of the system size, peaks are visible for ABPs. Squirmers exhibit a less pronounced translation order. Here, 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-range order, and g(r) assumes the ideal-gas value essentially on the scale of three to four neighbor distances. A decreasing packing fraction leads to a decreasing long-range order. For ϕ^2D = 0.5, ABPs 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.

Fig. 6 Two-dimensional pair-correlation functions for squirmers (pushers: β = −1, pullers: β = 1, neutral squirmers: β = 0) and ABPs for the packing fraction ϕ^2D = 0.6 and the Péclet number Pe = 575. For ABPs, the system sizes are L = 96a (solid) and L = 192a (dashed). Please note the logarithmic scale of the y-axis.

4.2.4 Distribution of hexagonal order parameters. 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₆|² defined as^24,35,68where N^k₆ is the set of the six nearest neighbors of swimmer k, and ϑ_kj is the angle between vector R_k − R_j connecting the centers of particles j and k and the x-axis. For a perfect hexagonal lattice |q₆|² = 1. Fig. 7 displays probability distribution functions of the order parameter for squirmers and ABPs at Pe = 115. Distributions for the higher Péclet number Pe = 575 closely agree with those shown in Fig. 7. Consistent with the high local packing, ABPs exhibit a pronounced probability for hexagonal order. The respective values for squirmers are smaller. Therefore, the tendency to form hexagonal order is most distinct for pullers. In contrast, pushers exhibit hardly any local hexagonal order.

Fig. 7 Probability distribution of the hexagonal order parameter |q₆|² of spherical ABPs (purple; L = 96a (solid), L = 192a (dashed)) and squirmers (pullers: β = −1, red; neutral: β = 0, black; pullers: β = 1, blue) for Pe = 115 and ϕ^2D = 0.6.

4.2.5 Distribution of propulsion direction. Strictly 2D ABP systems exhibit stronger aggregation and a more pronounced phase separation^24,35,52 compared to similar systems in quasi-2D.³⁵ The reason is that the freedom of the propulsion vector e to independently 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.

Fig. 8 Probability distribution functions of the propulsion-direction component e_y normal to the confining walls for ABPs and squirmers. The squirmers are preferentially orientated parallel to the walls. The Péclet number is Pe = 575 and the concentration ϕ^2D = 0.6. The distribution functions for Pe = 115 and the same density are nearly identical.

4.2.6 Orientational correlations. The suppression of MIC and MIPS is attributed to a strong reorientation (scattering) of squirmers by hydrodynamic torques during their collisions.²⁸ This conclusion is quantitatively confirmed by our simulations, as shown in Fig. 9, using the average orientation autocorrelation functionof squirmers and ABPs. 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 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 , 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 calculated on the basis of the rotational diffusion coefficient at infinite dilution.

Fig. 9 Orientation correlation function C_e(t) averaged over all swimmers for a system of spherical squirmers and active Brownian particles at ϕ^2D = 0.6 and Pe = 115. The dashed line is the theoretical expectation.

The relaxation time τ_sq depends on the squirmer concentration, as shown in Fig. 10(a) for Pe = 575. Compared to ABPs, the τ_sq of squirmers is reduced by orders of magnitude (factor 50–80) and decreases almost linearly in the considered range of concentrations (for Pe = 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)). 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. 28 for strictly 2D systems.

Fig. 10 (a) Decay time τ_sq of the orientational correlation function C_e(t) and (b) average swimming velocity as a 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/2D_R and the swimming velocity by U₀. The Péclet number is Pe = 575. The lines are guides to the eye (either linear or parabolic functions).

5 Simulations—spheroidal squirmers

To elucidate the influence of HI on the phase behavior of spheroidal squirmers, we focus on the Péclet number Pe = 12 and packing fractions 0.4 ≤ ϕ^2D ≤ 0.6. For this Péclet number and densities, spherical ABPs are far from a phase transition and in the fluid-like regime.⁵² The phase behavior of spheroidal ABPs depends on the aspect ratio; no systematic studies to obtain a phase diagram have been performed so far.

5.1 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 MPC time step , rotation angle α = 130°, and the mean number of particles 〈N_c〉 = 10, which yields the fluid viscosity .

To investigate the collective behavior for different packing fractions ϕ = 4πN_swb_x²b_z/3L²L_y, or equivalently area packing fractions ϕ^2D = N_swπb_xb_z/L², 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 , which yields, with eqn (3), the Péclet number

Pe = U₀/(2D^⊥_Rb_z) = 12 (9)

for all aspect ratios. Here, D^⊥_R = k_BT/ξ^⊥ is the rotational diffusion coefficient around the minor axis, and ξ^⊥ is the rotational friction coefficient provided in Appendix B. Note that this Péclet number is much lower than those of the previous section for our studies on MIPS of spheres. As discussed before, β ∈ {−1,0,1} and simulations of ABPs are performed for comparison.

5.2 Results

Fig. 11 illustrates structure formation for various aspect ratios. Although the Péclet number Pe = 12 is quite low, spheroidal squirmers clearly reveal clusters for all aspect ratios b_x/b_z ≥ 2. Therefore, a larger aspect ratio is beneficial for cluster formation. For spherical squirmers, no phase separation occurs at this Péclet number (Fig. 11(a)), consistent with the results of ref. 52. 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 Fig. 11(b) and (d).

Fig. 11 Snapshots of the structure of neutral (β = 0) spheroidal squirmers for Pe = 12, the area packing fraction ϕ^2D = 0.5, and the aspect ratios (a) b_z/b_x = 1, (b) 2, (c) 3, and (d) 4. See the ESI† for movies of the various cases (M3–M6).

5.2.1 Density distribution. Again, we employ Voronoi tessellations 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 measured as the Euclidean 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.⁶⁹ These 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 ABPs 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) emphasize this aspect. Here, we find a clear maximum in P(ϕ_loc) at the average density and, hence, no MIPS. Near-field HI by the stronger flow field lead to a faster reorientation dynamics and, hence, a less pronounced order.

Fig. 12 (a) Probability distribution of the local packing fraction for spheroidal squirmers with the aspect ratio b_z/b_x = 2 at the average area packing fraction ϕ^2D = 0.3, i.e., ϕ = 0.17. (b) Corresponding results for ϕ^2D = 0.5, i.e., ϕ = 0.29. The Péclet number is Pe = 12.

5.2.2 Phase diagram. 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 a dense cluster region and a dilute region. The respective phase borders should be viewed as approximations to the infinite-system-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 the 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.

Fig. 13 State diagram for squirmers and ABPs at Pe = 12. The blue, black, red, and purple areas indicate giant clusters for pullers (β = 1), neutrals (β = 0), pushers (β = 1), and active Brownian particles, respectively. Therefore clusters of pullers appear in all colored areas, clusters of neutral squirmers in the black, red, and purple area, etc. See the ESI† for individual state diagrams of the various swimmers.

6 Discussion

For spherical squirmers, the various structural quantities—distribution of local packing fraction, cluster-size distribution, pair correlation function, and probability distribution of hexagonal order parameter—reveal suppression of MIPS by hydrodynamic interactions over a wide-range of densities and Péclet numbers deep inside the MIPS parameter regime of ABPs. The distinct flow fields of pushers, pullers, and neutral squirmers influence the structural properties of the fluid, and MIC emerges, most pronounced for pullers and least for pushers. However, no phase separation occurs for packing fractions ϕ^2D ≲ 0.6. We did not explore the high-density part of the phase diagram, where further phases and transitions might be expected. Interestingly, HI reduce finite-system-size effects as distribution functions for the local packing fraction and the cluster size are essentially independent of the system size.

Our simulation results for spherical squirmers agree with certain aspects of previous studies, but disagree with others. In particular, the cluster-size distribution function has been analyzed in ref. 37 for the same squirmer model using the Lattice Boltzmann approach, but using a somewhat different simulation setup. Significantly larger systems with many more squirmers have been considered, and three-dimensional hydrodynamic interactions, but with squirmers confined in 2D at the density ϕ^2D = 0.1. In agreement with our Fig. 5, in ref. 37 a power-law decay (p(n) ∼ n^−2.0) for small cluster sizes (n ≲ 10³) is obtained for pullers with 0 < β ≲ 1. For other β (β > 2 and β < 0), a stronger decay appears for n ≳ 10. 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 system.

Our simulation results of spherical squirmers disagree qualitatively and quantitatively with the results of ref. 35 and 50, where comparable systems have been considered. We find that hydrodynamic interactions suppress cluster formation of squirmers, in contrast to ref. 35. According to ref. 50, the Péclet numbers of our neutral squirmers are deep in the phase-separated region of the phase diagram (Fig. 3 of ref. 50).⁷⁰ However, we do not see any signature for phase separation, not even for neutral squirmers. 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 the weakest order for pushers, in agreement with other studies.^36–38 The studies of ref. 35 show the 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 ref. 35 and 50 is the different compressibilities of the employed MPC fluids. It is shown in Appendix 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 ref. 35 and 50 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, HI strongly enhance phase separation even for Péclet numbers as low as Pe = 12 and densities, where spheroidal ABPs are deep inside the fluid state. Two aspects contribute to the hydrodynamic enhancement of cluster formation, namely near-field HI between squirmers, and between a squirmer and the confining walls.⁴⁶ The effect of the active stress, β, is evident from the state diagram of Fig. 13. Our simulation studies of ref. 46 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 in particular with the walls. Applying three-dimensional hydrodynamic interactions rather than no-slip boundary conditions to a still confined pair of pullers yields short cooperative motion only, emphasizing that surface hydrodynamics is extremely important.⁴⁶ 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.

7 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. 35, but in accordance with ref. 28. We attribute the contradiction with ref. 35 to compressibility effects of the MPC fluid. In agreement with ref. 28, we confirm that a hydrodynamically enhanced reorientation of squirmers during swimmer–swimmer collisions implies an effectively lower Péclet number and, hence, a suppression of phase separation. However, the squirmers form short-lived clusters, MIC, over a broad range of cluster sizes. As already pointed out in ref. 38, near-field hydrodynamic interactions play a major role in cluster formation. This is reflected in the clear differences between pullers, neutral squirmers, and pushers.³⁸ 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.⁵¹

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 lower 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 squirmers (β = 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.⁴⁶

Simulations of dumbbell swimmers^71,72 in three dimensions have also been performed and hydrodynamically enhanced MIC has been observed.⁷³ Since dumbbell swimmers are anisotropic, this finding is consistent with our results for spheroidal squirmers. However, we have to be careful with too general conclusions. As we demonstrated, MIC depends decisively on the hydrodynamic flow field, i.e., the details of the flow field matter. This is reflected by our results for very strong pushers (β = −5), which show a suppression of the 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.¹³ At the moment, there is no unambiguous understanding of the different behaviors and underlying mechanisms. Could the 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. 74. With the swimming velocity U^e₀ = 30 μs and the length 2b_z ≈ 10μ, we obtain a Péclet number Pe^e ≈ 50. Our simulations show that such a Pe even further enhances 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,⁷⁵ which forms 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.⁷⁶ 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.⁷⁷ Such a flow field induces a circular motion on a planar surface^77–80 and may also enhance the decorrelation of the orientational order. Another mechanism has been proposed in ref. 39 and 40. The particular geometry of a microswimmer, such as the inhomogeneous, elongated shape of bacteria, may lead to a coupling of the orientational dynamics and the fluid flow. This could be a mechanism for an enhanced rotational motion, and the emergence of mesoscale turbulence. Naturally, further mechanisms may exist, an example is a fore-aft asymmetry,⁸¹ which lead to the required enhanced reorientation dynamics, necessary for the suppression of MIPS and the maintenance of a dynamic state. Here, systematic studies are desirable to elucidate the relevance and importance of the various mechanisms and interactions.

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 the shape and hydrodynamics, which is revealed by our simulations, indicates that the details of the cell shape and the propulsion pattern may play a much more important role in the collective dynamics of biological microswimmers than previously expected.

Conflicts of interest

There are no conflicts to declare.

Appendix A: simulation approaches

A.1 Multiparticle collision dynamics (MPC)

In MPC,^47–49 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

r_i(t + h) = r_i(t) + hv_i(t). (10)

In the collision step, the particles interact locally through an instantaneous stochastic process, for which we apply the stochastic rotation dynamics (SRD) variant of the MPC method^47–49,82 with angular momentum conservation (MPC-SRD + a).^62,83,84 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^62,83Here, R(α) is the rotation matrix for the rotation of the relative velocity v_i,c = v_i − v_cm, with respect to center-of-mass velocity v_cm of the cell, by angle α around a randomly oriented axis. r_i,c = r_i − r_cm is the position relative to the center of mass r_cm of the cell, 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.⁸⁵ Thermal fluctuations are intrinsic to the MPC method.^62,86,87 Hence, a cell-level canonical thermostat (MBS thermostat) is applied after every collision step, which maintains the temperature T.^86,87 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.^48,88,89

We measure lengths in units of the collision cell size a and time in units of , where k_B is the Boltzmann factor, which corresponds to the choice a = m = k_BT = 1. The rotation angle is set to α = 130°. The viscosity η of the fluid can be tuned by changing the mean number of particles 〈N_c〉 in a cell or the time step h.^62,83,87

A.2 Boundary conditions and squirmer implementation

No-slip boundary conditions are implemented at the confining walls (cf.Fig. 2) using 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.^49,90 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 solved by introducing quaternions^46,67 (see also Appendix 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 (eqn (1)). The change of the 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_BT/m.^45,46 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. 46. 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. 46 (see also ref. 91).

A.3 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. The squirmer orientation is parameterized by a unit quaternion q = (q₀,q₁,q₂,q₃).⁶⁷ The integration scheme for the overdamped dynamics of the center-of-mass R(t) and quaternion q(t) of the spheroid is⁹²

R(t + h) = R(t) + μFh + U₀eh, (12)

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 discussed before, U₀ is the propulsion velocity and e ≡ D^Te_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 λ_qq(t) to obtain q̃ and then determine λ_q such that q̃ + λ_qq is normalized.⁹²

The mobility matrix μ relates the applied force F to the resulting ABP velocity UviaU = μF, while the rotational mobility matrix μ^R relates the applied torque T to the resulting angular velocity ΩviaΩ = μ^RT. In the body-fixed frame, the respective mobility matrices are diagonal, with diagonal elements given by the inverse translational and rotational friction coefficients γ_x⁻¹, γ_y⁻¹, γ_z⁻¹ and ξ_x⁻¹, ξ_y⁻¹, ξ_z⁻¹, respectively. For a prolate spheroid γ_x = γ_y = γ^⊥, γ_z = γ^‖, ξ_x = ξ_y = ξ^⊥, and ξ_z = ξ^‖. The parallel and perpendicular friction coefficients are provided in Appendix B. For zero eccentricity, i.e., a sphere with b_x = b_z = R, we obtain γ^‖ = γ^⊥ = 6πηR and ξ^‖ = ξ^⊥ = 8πηR³.

In the body-fixed frame, the Cartesian components of the thermal force F^th are Gaussian distributed random numbers of zero mean and variance 2k_BTγ_α/h (α ∈ {x,y,z}). The thermal torque T^th is generated similarly, but with variance 2k_BTξ_α/h. The implementation of the repulsive forces and torques due to steric interactions between ABPs and between an ABP and a wall is described in the appendix of ref. 46.

Appendix B: rigid-body dynamics: quaternions

The rotation matrix D and the matrix Q introduced in eqn (13) are given byin terms of the rotation quaternion q = (q₀,q₁,q₂,q₃).⁶⁷

The translational and rotational friction coefficients of a spheroid are given by

where η is the fluid viscosity, ê = 1/τ₀ = c/b_z is the eccentricity, and L = log((1 + ê)/(1 − ê)). For ê = 0, i.e., a sphere with b_x = b_z = R, follows γ^‖ = γ^⊥ = 6πηR and ξ^‖ = ξ^⊥ = 8πηR³.

Appendix C: avoiding MPC simulation artifacts related to density inhomogeneities

MPC is a compressible fluid and correspondingly strong density inhomogeneities can emerge for certain choices of MPC parameters. Fig. 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 Pe = 115. However, the Péclet number is realized by different combinations of the MPC fluid and squirmer parameters. In Fig. 14(a), we set and 〈N_c〉 = 10, and in (b) and 〈N_c〉 = 80. The latter gives a ten times larger viscosity. The other parameters are the same as defined in Appendix A.1 and Section 4.1.

Fig. 14 Fluid particle density ρ in a slit with spherical neutral squirmers at Pe = 115 and ϕ^2D = 0.6. The total fluid density is shown, comprising the MPC particle density ρ_fl and that of the phantom particles ρ_ph, i.e., ρ = m〈N_c〉/a³ = ρ_fl + ρ_ph. ρ₀ denotes the average density. In (a), the choice and 〈N_c〉 = 10 yields the pumping number Pu = 5, and in (b), and 〈N_c〉 = 80 yields Pu = 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).

Fig. 14(a) shows large density inhomogeneities with dense (orange) and rarefied regions (purple). Therefore, the rarefied regions coincide 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 built up by active clustered squirmers by expelling the fluid from 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, hence, depends strongly on the location.

The mechanism in a MPC fluid, which opposes the diminished density in the gaps between squirmers is fluid-particle diffusion. The ratio of time t_diff 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 velocity U₀ over the same distance, yields the dimensionless pumping number

where is the MPC-fluid diffusion coefficient. For small pumping number, Pu, we expect the density inhomogeneities to disappear, i.e., MPC fluid-particle diffusion is faster than activity-induced advection. The above choices of parameters yield the pumping numbers Pu = 5 (Fig. 14(a)) and Pu = 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, Pu is given by Pu = σ²D_RPe/6D_f. Since Pe ≫ 1 for active systems, Pu < 1 requires σ²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.

Acknowledgements

We thank A. Wysocki and R. Hornung for helpful discussions. The support of this work from the DFG priority program SPP 1726 “Microswimmers – from Single Particle Motion to Collective Behaviour” is gratefully acknowledged. We also acknowledge the partial financial support from the Bavarian Ministry of Economic Affairs, Energy and Technology within joint projects between Helmholtz institutes. In addition, the authors gratefully acknowledge the computing time granted for JARA-HPC on supercomputer JURECA at Forschungszentrum Jülich.

References

1. J. Elgeti, R. G. Winkler and G. Gompper, Rep. Prog. Phys., 2015, 78, 056601.
2. M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao and R. A. Simha, Rev. Mod. Phys., 2013, 85, 1143.
3. F. Peruani, J. Starruß, V. Jakovljevic, L. Søgaard-Andersen, A. Deutsch and M. Bär, Phys. Rev. Lett., 2012, 108, 098102.
4. A. P. Petroff, X.-L. Wu and A. Libchaber, Phys. Rev. Lett., 2015, 114, 158102.
5. D. B. Kearns, Nat. Rev. Microbiol., 2010, 8, 634.
6. N. C. Darnton, L. Turner, S. Rojevsky and H. C. Berg, Biophys. J., 2010, 98, 2082.
7. X.-L. Wu and A. Libchaber, Phys. Rev. Lett., 2000, 84, 3017.
8. C. Dombrowski, L. Cisneros, S. Chatkaew, R. E. Goldstein and J. O. Kessler, Phys. Rev. Lett., 2004, 93, 098103.
9. M. F. Copeland and D. B. Weibel, Soft Matter, 2009, 5, 1174.
10. R. Zhang, L. Turner and H. C. Berg, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 288–290.
11. N. H. Mendelson, A. Bourque, K. Wilkening, K. R. Anderson and J. C. Watkins, J. Bacteriol., 1999, 181, 600.
12. H. H. Wensink and H. Löwen, J. Phys.: Condens. Matter, 2012, 24, 460130.
13. 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.
14. A. Sokolov and I. S. Aranson, Phys. Rev. Lett., 2012, 109, 248109.
15. R. Golestanian, T. B. Liverpool and A. Ajdari, Phys. Rev. Lett., 2005, 94, 220801.
16. I. Theurkauff, C. Cottin-Bizonne, J. Palacci, C. Ybert and L. Bocquet, Phys. Rev. Lett., 2012, 108, 268303.
17. S. Thutupalli, R. Seemann and S. Herminghaus, New J. Phys., 2011, 13, 073021.
18. J. Palacci, S. Sacanna, A. P. Steinberg, D. J. Pine and P. M. Chaikin, Science, 2013, 339, 936.
19. I. Buttinoni, J. Bialké, F. Kümmel, H. Löwen, C. Bechinger and T. Speck, Phys. Rev. Lett., 2013, 110, 238301.
20. M. C. Marchetti, Y. Fily, S. Henkes, A. Patch and D. Yllanes, Curr. Opin. Colloid Interface Sci., 2016, 21, 34.
21. C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe and G. Volpe, Rev. Mod. Phys., 2016, 88, 045006.
22. Y. Fily and M. C. Marchetti, Phys. Rev. Lett., 2012, 108, 235702.
23. G. S. Redner, M. F. Hagan and A. Baskaran, Phys. Rev. Lett., 2013, 110, 055701.
24. J. Bialké, T. Speck and H. Löwen, Phys. Rev. Lett., 2012, 108, 168301.
25. A. Wysocki, R. G. Winkler and G. Gompper, EPL, 2014, 105, 48004.
26. J. Stenhammar, D. Marenduzzo, R. J. Allen and M. E. Cates, Soft Matter, 2014, 10, 1489.
27. A. Wysocki, R. G. Winkler and G. Gompper, New J. Phys., 2016, 18, 123030.
28. R. Matas-Navarro, R. Golestanian, T. B. Liverpool and S. M. Fielding, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 032304.
29. M. E. Cates and J. Tailleur, Annu. Rev. Condens. Matter Phys., 2015, 6, 219.
30. R. Golestanian, Phys. Rev. Lett., 2012, 108, 038303.
31. O. Pohl and H. Stark, Phys. Rev. Lett., 2014, 112, 238303.
32. F. Peruani, A. Deutsch and M. Bär, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 74, 030904.
33. Y. Yang, V. Marceau and G. Gompper, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 82, 031904.
34. M. Abkenar, K. Marx, T. Auth and G. Gompper, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 88, 062314.
35. A. Zöttl and H. Stark, Phys. Rev. Lett., 2014, 112, 118101.
36. B. Delmotte, E. E. Keaveny, F. Plouraboué and E. Climent, J. Comput. Phys., 2015, 302, 524.
37. F. Alarcón, C. Valeriani and I. Pagonabarraga, Soft Matter, 2017, 13, 814.
38. N. Yoshinaga and T. B. Liverpool, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2017, 96, 020603.
39. S. Heidenreich, J. Dunkel, S. H. L. Klapp and M. Bär, Phys. Rev. E, 2016, 94, 020601.
40. H. Reinken, S. H. L. Klapp, M. Bär and S. Heidenreich, Phys. Rev. E, 2018, 97, 022613.
41. M. J. Lighthill, Comm. Pure Appl. Math., 1952, 5, 109.
42. J. R. Blake, J. Fluid Mech., 1971, 46, 199.
43. T. Ishikawa and M. Hota, J. Exp. Biol., 2006, 209, 4452.
44. I. Pagonabarraga and I. Llopis, Soft Matter, 2013, 9, 7174.
45. I. O. Götze and G. Gompper, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 82, 041921.
46. M. Theers, E. Westphal, G. Gompper and R. G. Winkler, Soft Matter, 2016, 12, 7372.
47. A. Malevanets and R. Kapral, J. Chem. Phys., 1999, 110, 8605.
48. R. Kapral, Adv. Chem. Phys., 2008, 140, 89.
49. G. Gompper, T. Ihle, D. M. Kroll and R. G. Winkler, Adv. Polym. Sci., 2009, 221, 1.
50. J. Blaschke, M. Maurer, K. Menon, A. Zöttl and H. Stark, Soft Matter, 2016, 12, 9821.
51. A. J. T. M. Mathijssen, A. Doostmohammadi, J. M. Yeomans and T. N. Shendruk, J. Fluid Mech., 2016, 806, 35.
52. P. Digregorio, D. Levis, A. Suma, L. F. Cugliandolo, G. Gonnella and I. Pagonabarraga, Phys. Rev. Lett., 2018, 121, 098003.
53. S. C. Takatori, W. Yan and J. F. Brady, Phys. Rev. Lett., 2014, 113, 028103.
54. S. C. Takatori and J. F. Brady, Soft Matter, 2015, 11, 7920.
55. R. G. Winkler, A. Wysocki and G. Gompper, Soft Matter, 2015, 11, 6680.
56. S. E. Spagnolie and E. Lauga, J. Fluid Mech., 2012, 700, 105.
57. F. Alarcón and I. Pagonabarraga, J. Mol. Liq., 2013, 185, 56.
58. S. R. Keller and T. Y. Wu, J. Fluid Mech., 1977, 80, 259.
59. A. M. Leshansky, O. Kenneth, O. Gat and J. E. Avron, New J. Phys., 2007, 9, 145.
60. K. Ishimoto and E. A. Gaffney, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 88, 062702.
61. T. Ishikawa, M. P. Simmonds and T. J. Pedley, J. Fluid Mech., 2006, 568, 119.
62. M. Theers and R. G. Winkler, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91, 033309.
63. A. Zöttl and H. Stark, Phys. Rev. Lett., 2012, 108, 218104.
64. C. H. Rycroft, Chaos, 2009, 19, 041111.
65. M. S. S. Challa, D. P. Landau and K. Binder, Phys. Rev. B: Condens. Matter Mater. Phys., 1986, 34, 1841.
66. J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids, Academic Press, London, 1986.
67. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1987.
68. P. J. Steinhardt, D. R. Nelson and M. Ronchetti, Phys. Rev. B: Condens. Matter Mater. Phys., 1983, 28, 784.
69. P. Persson and G. Strang, SIAM Rev., 2004, 46, 329.
70. Note that the definition of the Péclet number of ref. 50 is different; it is three times ours.
71. A. Suma, G. Gonnella, G. Laghezza, A. Lamura, A. Mossa and L. F. Cugliandolo, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 052130.
72. J. T. Siebert, J. Letz, T. Speck and P. Virnau, Soft Matter, 2017, 13, 1020.
73. A. Furukawa, D. Marenduzzo and M. E. Cates, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 022303.
74. K. Drescher, J. Dunkel, L. H. Cisneros, S. Ganguly and R. E. Goldstein, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 10940.
75. H. C. Berg, E. Coli in Motion, Springer, New York, 2004.
76. T. Eisenstecken, J. Hu and R. G. Winkler, Soft Matter, 2016, 12, 8316.
77. J. Hu, A. Wysocki, R. G. Winkler and G. Gompper, Sci. Rep., 2015, 5, 9586.
78. E. Lauga, W. R. DiLuzio, G. M. Whitesides and H. A. Stone, Biophys. J., 2006, 90, 400.
79. R. Di Leonardo, D. Dell'Arciprete, L. Angelani and V. Iebba, Phys. Rev. Lett., 2011, 106, 038101.
80. L. Lemelle, J.-F. Palierne, E. Chatre, C. Vaillant and C. Place, Soft Matter, 2013, 9, 9759.
81. Q.-S. Chen, A. Patelli, H. Chaté, Y.-Q. Ma and X.-Q. Shi, Phys. Rev. E, 2017, 96, 020601.
82. T. Ihle and D. M. Kroll, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2001, 63, 020201(R).
83. H. Noguchi and G. Gompper, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 016706.
84. M. Yang, M. Theers, J. Hu, G. Gompper, R. G. Winkler and M. Ripoll, Phys. Rev. E, 2015, 92, 013301.
85. T. Ihle and D. M. Kroll, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 67, 066705.
86. C.-C. Huang, R. G. Winkler, G. Sutmann and G. Gompper, Macromolecules, 2010, 43, 10107.
87. C.-C. Huang, A. Varghese, G. Gompper and R. G. Winkler, Phys. Rev. E, 2015, 91, 013310.
88. T. Ihle, Phys. Chem. Chem. Phys., 2009, 11, 9667.
89. C.-C. Huang, G. Gompper and R. G. Winkler, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 86, 056711.
90. A. Lamura, G. Gompper, T. Ihle and D. M. Kroll, Europhys. Lett., 2001, 56, 319–325.
91. L. Paramonov and S. N. Yaliraki, J. Chem. Phys., 2005, 123, 194111.
92. I. M. Ilie, W. J. Briels and W. K. D. Otter, J. Chem. Phys., 2015, 142, 114103.

---

Footnote

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

---