Emergence of active turbulence in microswimmer suspensions due to active hydrodynamic stress and volume exclusion

Microswimmers exhibit an intriguing, highly-dynamic collective motion with large-scale swirling and streaming patterns, denoted as active turbulence – reminiscent of classical high-Reynolds-number hydrodynamic turbulence. Various experimental, numerical, and theoretical approaches have been applied to elucidate similarities and differences of inertial hydrodynamic and active turbulence. We use squirmers embedded in a mesoscale fluid, modeled by the multiparticle collision dynamics (MPC) approach, to explore the collective behavior of bacteria-type microswimmers. Our model includes the active hydrodynamic stress generated by propulsion, and a rotlet dipole characteristic for flagellated bacteria. We find emergent clusters, activity-induced phase separation, and swarming behavior, depending on density, active stress, and the rotlet dipole strength. The analysis of the squirmer dynamics in the swarming phase yields Kolomogorov-Kraichnan-type hydrodynamic turbulence and energy spectra for sufficiently high concentrations and a strong rotlet dipole. This emphasizes the paramount importance of the hydrodynamic flow field for swarming motility and bacterial turbulence. Microswimmers are ubiquitous in nature and present highly-dynamic collective motion. This study presents computational simulations of spheroidal model-microswimmers confined between two walls, revealing two collective dynamic states, giant motile clusters and active turbulence, depending on density and hydrodynamic interactions.

Fundamental insight into hydrodynamic turbulence is achieved via velocity correlation functions [47]. In particular, Kolmogorov predicted the universal power-law dependence for the energy spectrum E ∼ k −κ on the wavenumber k = |k|, with κ = 5/3 [41,47]. In fact, this relation applies for two-(2D) and three-dimensional (3D) systems [42]. Numerous studies on active systems reveal a wide spectrum of possible turbulent characteristics dependent on their constituents and the detailed (microscopic) interaction mechanisms, reflected in a wide range of exponents deviating from the Kolmogorov value, see Tab. I. Experiments on B. subtilis and E. coli bacteria [13,34] yield exponents significantly above and below the Kolmogorov value. Computer simulations employing various models have been performed and the energy spectrum has been calculated. Nonhydrodynamic particle-based simulations of an extension of the Vicsek model [48], accounting for short-range parallel and large-range antiparallel alignment, yield the same exponent [45] as in experiments on E. coli [13]. Simulations of self-propelled rodlike particles give a value close to the Kolmogorov value [13,44]. Lattice Boltzmann simulations of microswimmers represented by extended force dipoles (point particles) produce seemingly turbulent behavior for sufficiently large swimmer densities [46] (see Table I). For active nematics, the route to chaotic behavior has been studied experimentally and theoretically [23,49]. Their dynamics is characterized by an intrinsic I. Various aspects of experimentally, theoretically, and by simulations studied systems exhibiting features of active turbulence. The articles (Ref.) discuss microscale systems exhibiting a power-law energy spectrum E(k) ∼ k −κ for km < k < kc and E(k) ∼ kκ for k < km, with km and kc the wavenumber of the maximum in the energy spectrum and that of the microswimmer characteristic length, respectively. Note that in active nematic theory, km = 2π/la [18,39]. Cells comprise canine kidney, endothelial, myoblast, and fibroblast cells. Abbreviations: act. nem.: active nematics, SPR: self-propelled rod, LB: Lattice Boltzmann, MPC: multiparticle collision dynamics, HI: hydrodynamic interactions, φ: packing fraction. Symbols: " " aspect is present, "−" aspect is absent, "/" aspect has not been analyzed/considered . length scale l a , where l a is determined by the balance between the active and nematic elastic stress [18,39], and the creation and annihilation of topological defects. In addition, various theoretical studies have been performed with [39] and without [18] defects, where both yield similar energy spectra with distinct power-law exponents for length scales larger and smaller than l a (Tab. I). In contrast, we expect hydrodynamic interactions to dominate the chaotic and turbulent behavior in bacterial suspensions. Hence, it is a priori not evident that both types of chaotic dynamics exhibit the same kind of turbulent behavior, taken into account the disparity in the exponents κ andκ. There are two particular systems of mesoscopic active particles, namely spinners -short rodlike self-organized colloidal structures rotated by an external magnetic field [50] -and Marangoni surfers [28], where turbulent dynamics consistent with Kolmogorov scaling has been observed. Their Reynolds numbers Re ∼ O(10) are much smaller than that of classical inertial turbulence, but are much larger than those of microswimmer systems, where Re 1. As a major difference to hydrodynamic turbulence, various experimental and simulation studies of active turbulence suggest the presence of a characteristic upper length scale for the vortex size, only below which the energy spectrum decreases in a power-law manner with increasing wavenumber k [31,33,35]. This scale is typically on the order of ten microswimmer lengths. Theoretical studies based on a continuum approach [13,33,40,51], where the velocity field is described by the incompressible Toner-Tu equation [52,53] combined with a Swift-Hohenberg term [54] for pattern formation, support this observation. However, in contrast to high-Reynolds num-ber hydrodynamic turbulence, the internal stress due to self-propulsion and polar alignment interactions of the active agents is important, which, combined with the fluid dynamics, determines the vortex size [51].
The diversity of obtained energy spectra and characteristic power laws (Tab. I) indicates a strong dependence of the collective behavior on the detailed microswimmer interactions. Yet, it is not clear to which extent and under what circumstances hydrodynamic interactions are important.
By systematically varying the squirmer density, the active stress, and the rotlet dipole strength, our simulations provide insight into their influence on the collective dynamics of microswimmers. The combination of active stresses and a non-zero rotlet dipole suppresses phase separation and promotes swarming motility.
The analysis of the swarming phases reveals turbulentlike motion, where the energy spectrum displays powerlaw decays below the characteristic length scale discussed above, however, with an exponent depending on the squirmer concentration. Remarkably, we find the value κ = 5/3 for our largest density, strong active stress, and a non-zero rotlet dipole, consistent with the Kolmogorov prediction.

System setup
In our simulations, N sq prolate spheroidal squirmers with the semi-major, b z , and -minor, b x , axis are confined in a three-dimensional narrow slit between two parallel walls and periodic boundary conditions along the x and z direction (Fig. 1). The prescribed squirmer surface velocity yields swimming with the velocity v 0 , an active stress of strength β, and a rotlet dipole of strength λ (Sec. ). The embedding fluid is modeled via the multiparticle collision dynamics (MPC) method [55,56], applying the stochastic-rotation variant with angular momentum conservation (MPC-SRD+a) [74,75].

Structural properties
The simulation snapshots of Fig. 2 illustrate emergent structures for the various considered packing fractions, active stresses, and rotlet dipole strengths. Distinct motility patterns can be identified: (i) Motility-induced phase separation (A-MIPS) for |β| ≥ 1, λ = 0, φ 0.3. Since here the shape of the spheroids implies squirmer alignment and the formation of polar motile clusters, we use the notation A-MIPS, to distinguish it from the case of isotropic, non-aligning particles, which form immobile clusters (MIPS) [15,67,68]. (ii) Swarming motility for |β| > 1, λ = 4, φ 0.5, and (iii) gas of small clusters for φ 0.3. The clusters emerging by A-MIPS increase with increasing packing fraction and are system-spanning for φ 0.5, consistent with our previous studies [57]; they are denoted as local and global clusters in Ref. [16]. The clusters are rather dynamic and exhibit translational and rotational motion. In the dense swarming phase, clusters of squirmers migrate collectively, thereby forming dynamic swirling and streaming patterns [10,12,16,73]. A quantitative criterion for the classification into A-MIPS and swarming motility will be provided in terms of the cluster-size distribution function (Fig. 4). Some of the small clusters for φ 0.3 exhibit cooperative motion, where a few squirmers move together for some time. In general, the rotlet dipole enhances cluster formation, and squirmers align side by side, which is clearly visible for φ 0.4. The precise mechanism for this cooperative motion is unexplored, but could depend on squirmer wall interactions. In contrast, for larger packing fractions the rotlet dipole suppresses A-MIPS and enhances swarming.

Local packing fraction
Clustering and A-MIPS of the squirmers are analyzed quantitatively by a Voronoi tessellation of the accessible volume [57,70,76,77]. Figure 3 provides examples of density distributions for the average packing fractions φ = 0.4 and 0.6. The pronounced peak at the local packing fraction φ loc ≈ 0.75 for φ = 0.4, β = −1, and λ = 0 indicates A-MIPS, with a dense phase in contact with a dilute phase, consistent with the snapshots of Fig. 2. Results for large |β| imply a disintegration of the large aggregate and ultimately, for β < −3, P φ displays a maximum at the average packing fraction, which indicates the absence of phase separation. Similarly, at φ = 0.6, the peaks in Fig. 3(b) for λ = 0 indicate phase separation, even for β as negative as β = −5. The rotlet dipole prevents formation of large clusters, but even for β = −5 and λ = 4 a broad range of cluster sizes exists.

Cluster-size distribution
The cluster-size distribution function represents the fraction of squirmers belonging to a cluster of size n, where p(n) is the number of clusters of size n. The distribution is normalized such that Nsq n=1 N (n) = 1. We use a distance and an orientation criterion to define a cluster: a squirmer belongs to a cluster, when its closest distance to another squirmer of the cluster is d s < 1.8(2 1/6 − 1)σ s and the angle between the orientations of the two squirmers is < π/6 (see Methods section for the definition). The latter allows us to identify different clusters even at high packing fractions.
The cluster-size distribution function is a useful quantity to characterize the motility pattern of a microswimmer system [16,78]. In the homogeneous phase, the distribution function decays exponentially, whereas a second peak (bimodal distribution) indicates the formation  of giant clusters (A-MIPS). At the percolation transition, N becomes scale free and decays by a power law, N ∼ x −γ [78]. The swarming phase is characterized by a power-law decay with an exponential cut-off and a characteristic scale determined by an average vortex size [16]. The distribution functions presented in Fig. 4 confirm our above conclusions on the emergent phases and motility patterns.
The probability distribution functions of the local packing fraction (Fig. 3) and cluster-size distribution functions (Figs. 4) clearly reveal a marked effect of the rotlet dipole on the collective behavior of the squirmers. In particular, A-MIPS is suppressed, but formation of highly dynamic clusters prevails, with a rather broad distribution of cluster sizes for high squirmer densities.

Rotational diffusion
An individual squirmer in the slit exhibits rotational diffusion around a minor body axis. Interactions between squirmers, either steric or by their flow fields, change their diffusive behavior substantially [57,60]. Figure 5(a) displays the time dependence of the autocorrelation function e(t)·e(0) of the propulsion direction of the squirmers. The various curves reflect a marked dependence of the rotational dynamics on the active stress and the rotlet dipole strength. The correlation function of the systems for (β, λ) = (−1, 0), (−3, 0), (−1, 4) exhibit a non-singleexponential decay. Steric interactions between squirmers with a preference to cluster formation as well as between finite-size clusters lead to a rotation of whole clusters, which implies a faster decay of the rotational correlation compared to thermal fluctuations alone (cf. movie M4) [80].
We characterize the rotational motion by fitting the initial decay of the correlation function with the expo-  nential as displayed in Fig. 5(a). The factor C 0 R ≈ 1.03 is included to account for a non-exponential decay for very short times. Squirmers with large active stresses and a rotlet dipole ((β, λ) = (−5, 0), (−3, 4), (−5, 4)) exhibit an exponentially decaying correlation function of C R over more than a order of magnitude. The extracted rotational diffusion coefficients D R obey D R /D 0 R > 1 (Fig. 5(b)), which reveals an accelerated rotational motion by shape-induced steric interactions and hydrodynamic flow fields. Note that D 0 R in a dilute system is independent of β. The diffusion coefficient D R increases with increasing squirmer concentration, reaches a packing fraction-dependent maximum and decreases again for  larger φ. An increasing number of squirmer contacts with increasing φ (φ 0.5) leads to a faster reorientation. However, at larger φ, clusters are formed, which move collectively and more persistently, which reduces D R . The larger D R values for larger |β| demonstrate the substantial contribution of active stress to the reorientation of the squirmers. At smaller φ and β < −1, the presence of a rotlet dipole with λ = 4 evidently reduces D R compared to that for λ = 0, which is associated with the appearance of small clusters of side-by-side swimming squirmers (cf. Fig. 2 and movie M4). In contrast, at high packing fractions, a rotlet dipole implies a larger D R as a consequence of an enhanced orientational motion of smaller clusters, specifically at large |β| = 5.

Mean square displacement
The mean-square displacement of the squirmers at high packing fractions (φ ≥ 0.6, Fig. 6) exhibit the typical ballistic motion for short times and a crossover to a diffusive motion for long times tD 0 R 0.1 [15,67], at least for systems with λ = 4. (The resolution of the longtime behavior of the phase separated systems for λ = 0 requires longer simulations.) There is only a slight difference in the swimming speed of the various squirmers at short times. The presence of a rotlet dipole causes an earlier deviation from a strict ballistic motion toward a ballistic-like motion with an exponent somewhat smaller than 2 as time increases compared to squirmers without such a dipole. Most remarkable, the systems with (β, λ) = (−5, 0), (−1, 4), (−3, 4), (−5, 4) exhibit a crossover from a ballistic or near ballistic to a diffusive motion at a displacement roughly corresponding to 12b z , i.e., 6 squirmer lengths. We may consider this as a characteristic length scale in the system, separating the scale of persistent motion from that of diffusive motion.

Velocity distribution function
Thermal and active fluctuations imply strongly varying instantaneous squirmer velocities, with magnitudes exceeding the swimming velocities by far. Hence, for the calculation of the velocity distribution function, we determine a swimming velocity by the difference quotient During the selected time interval ∆t = 10 3 ma 2 /(k B T ), a squirmer moves at most the distance 2b z /3. The distribution function P v of the Cartesian in-plane velocities ∆v = (v x/z −v x/z ) -the two spatial dimensions are equivalent -, wherev x/z are the average velocities along the Cartesian directions x and z, of an single squirmer (dilute system) in the slit, is Gaussian due to the thermal noise of the fluid. (The velocitiesv x/z are typically very small and non-zero only due to finitesize effects and statistical inaccuracy.) Collective effects modify the distribution function and P v deviates from a Gaussian in general, as shown in Fig. 7. Even for a pronounced active stress, β = −5, the distribution functions for packing fractions φ < 0.6 deviate from a Gaussian ( Fig. 7(a)) independent of λ -the curves are flattened at the maximum and are wider or narrower in the tails. Similar, at φ = 0.68 ( Fig. 7(b)), P v is broadened for all systems with λ = 0, as well as for β = −1 and λ = 4, although the distribution function are close to a Gaussian. Remarkably, the velocity distribution function for the system with φ = 0.68 and (β, λ) = (−5, 4) is very well described by a Gaussian despite pronounced collective swimming. Evidently, steric and flow-field interactions induce sufficient randomness to yield isotropic two-dimensional Gaussian distributed velocities. This aspect is particularly relevant for active turbulence, because velocities in high-Reynolds-number turbulent flows are Gaussian distributed [13,33] Active turbulence The characteristic features of the swimmer flow fields at higher densities are illustrated in Fig. 8. The clusters depicted in Fig. 8(a) exhibit a chaotic collective motion with regions of low and high velocity ( Fig. 8(b)) and vorticity (Fig. 8(c)) (see movies M3, M5, and M6 for the packing fraction φ = 0.68). The patterns are similar to those observed in experiments on bacteria [13,16,31,33,35], previous simulations [13,45], and continuum theory [13,32,33].

Spatial velocity correlation function
Quantitative insight into the turbulent dynamics of the squirmers is obtained by their spatial velocity correlation function, a concept well established in classic hydrodynamic turbulence [13,41,43,47]. For the discrete particle system, we define the spatial velocity correlation function as [44,70,81] C where r i is the center-of-mass position of squirmer i. Moreover, we introduce a normalized velocity correlation function as C 0 v (R) = C v (R)/c 0 , with c 0 = i v 2 i /N sq . (For an homogeneous and isotropic system, C v (R) is a function of R = |R| only.) Results of C 0 v for the packing fractions φ = 0.4 and 0.6 are presented in Fig. 9.  5, 4)). The case (i) corresponds to long-range correlations over the entire simulation box, consistent with A-MIPS and the appearance of a large cluster (Fig. 4). As shown in Fig. 9(b), such C 0 v (R) can be fitted by the function Specifically for φ = 0.6, we obtain the parameters of Tab. III. The respective velocity correlation functions decay approximately exponentially, with characteristic lengths scales between 2.3 and 5.6 swimmer lengths. The smaller value ξ/(2b z ) = 2.3 for λ = 4 indicates that a non-zero rotlet dipole implies weaker spatial correlation and, hence, smaller clusters. The distinct decay patterns support our conclusion on the motility mode as discussed in relation the cluster-size distribution functions (Fig. 4). However, a clear-cut separation of swarming and cluster dynamics is difficult to establish based on C v (R). An important feature of bacterial turbulence is a finite vortex size, which marks a characteristic length scale in the system and is reflected in a minimum of the velocity correlation function [13,30,33,81]. Our simulations yield such a minimum, e.g., for φ = 0.6, 0.68, (β, λ) = (−5, 4). Hence, we expect such squirmer system to exhibit active turbulence. A characteristic length scale can also exist for lower densities, e.g., for φ = 0.4, (β, λ) = (−5, 0), (−5, 4), where only small clusters are present. We would not denote the dynamics of such systems as turbulent. Insight into the turbulent behavior is gained by the energy spectrum which is obtained as Fourier transform of the spatial velocity correlation function (5) [47], and manifests the distribution of kinetic energy over different length scales.
In the calculation of E(k), we apply a left-shift of the correlation function C v (R) (Fig. 9) such that the decay starts at R = 0 in order to avoid artifacts in the Fourier transformation by a truncated correlation function. As for bacterial suspensions, the energy injection scale is the length scale of a microswimmer (2b z ), which yields the characteristic (maximum) wavenumber k c = π/b z for our squirmers. Figure 10 displays the energy spectrum for (β, λ) = (−5, 4) and the two packing fractions φ = 0.6, 0.68, and various system sizes. The simulations show two powerlaw regimes for a given density, namely E(k) ∼ kκ for k < k m and E(k) ∼ k −κ for k m < k < k c , with k m corresponding to the peak position of E(k). Such a maximum in E(k) is a feature of microswimmer active turbulence, and reflects a characteristic vortex size [13,33,35]. Our simulations yield approximate vortex sizes of 5 (10b z )  and 10 squirmer lengths (20b z ) for φ = 0.6 and 0.68, respectively. They are roughly consistent with the patterns of Fig. 8, the crossover from ballistic to diffusive motion in the mean-square displacement of Fig. 6, and the minimum of the correlation function of Fig. 9(b). Vortex sizes on the order of 5 − 10 microswimmer lengths are also found in experiments [13,34,35].
For k m < k < k c , corresponding to R > 2b x , our simulations yield turbulent flow patterns (Fig. 8). The exponent of the scaling regime depends on the squirmer density, with the values κ = −2 for φ = 0.6 and κ = −5/3 for φ = 0.68. The latter is consistent with the Kolmogorov-Kraichnan prediction for classical 2D turbulence [42]. This is remarkable, considering the wide scatter of exponents found in simulations and experiments (cf. Tab. I). Density seems to play an important role for the observed turbulent behavior. The squirmers of both densities exhibit swarming, namely, collective motion with largescale swirling and streaming patterns. However, only the dynamics in the higher density system exhibits the exponent κ = 5/3.
In the small k-value regime, we obtain the exponents  10. Energy spectrum. Energy spectra of systems with β = −5 and λ = 4 for the packing fractions φ = 0.6 (red) and 0.68 (blue). Various system sizes (see legend) have been explored in order to verify absence of finite-size effects. The dashed lines indicate power-laws in the respective regimes. κ = 1 for φ = 0.68 andκ = 5/3 for φ = 0.6, which reflect an increase of the energy with increasing k. The dependence k 5/3 is consistent with that observed theoretically and experimentally in Ref. [13], as well as in simulations [45]. However, other studies yield rather different dependencies (Tab. I). Theoretical models suggest that the small-k slope is governed by finite-system-size effects, i.e., depends in the boundary condition and physical parameters [40]. The curves in Fig. 10 reflect a weak dependence on the system size.

|E|/(ak B T /m)
The presences of a small-distance cut-off, where energy input by the squirmers occurs, and the peak in E(k), corresponding to a characteristic vortex size, limits the krange over which the energy spectrum decays in a powerlaw manner. This is in stark contrast to classical high-Reynolds-number turbulence, where the energy cascade extents over many orders of magnitude.

CONCLUSIONS
We have performed large-scale mesoscale hydrodynamics simulations of spheroidal squirmers in a narrow slit in order to analyze the emerging structures, motility patterns, and turbulent behavior for various packing fractions, active stresses, and rotlet-dipole strengths.
Our studies reveal a strong dependence of the motility pattern on the microswimmer concentration and their propulsion-induced flow field. The classification of the distinct motion pattern into the various categories -swimming and collective motion of very small clusters (cluster gas), phase separation by activity and anisotropic swimmer shape (A-MIPS), and swarmingis accomplished by visual inspection of snapshots (Fig. 2) and the characteristic features of the cluster-size distri-bution function (Fig. 4). A-MIPS appears for small active stresses, |β| 3, and all packing fractions φ > 0.2. Squirmers with stronger forces dipoles, |β| 3, at concentrations φ < 0.4 exhibit small clusters and strong cooperative effects for λ = 4. At higher packing fractions, φ ≥ 0.5, a swarming phase appears for λ = 4, where clusters of squirmers move collectively, and even exhibit active turbulence for high packing fractions (φ = 0.6, 0.68) and sufficiently large |β| (Fig. 2). Importantly, the rotlet dipole suppresses A-MIPS.
Our simulations clearly reveal the difficulty to characterize turbulence in active systems. Even more fundamental is the question, which criteria should be applied to classify a mesoscale system as turbulent, Considering microswimmer systems, chaotic flow patterns are evidently not sufficient. Inspired by classical hydrodynamic turbulence, we propose the following "minimal" criteria: • Reynolds numbers Re < 1 • presence of chaotic flow patterns with large-scale collective behavior • characteristic vortex size and a negative velocity correlation function • Gaussian velocity distribution function of the microswimmer's Cartesian velocity components • energy spectrum with power-law decay E(k) ∼ k −κ , κ > 0, on length scales below the characteristic vortex size.
The presence of small and large length-scale cut-offs by the microswimmer and vortex size implies a universal, scale-free behavior only over a limited range of length scales.
Analyzing the swarming motion of the squirmers, we find non-Gaussian distribution functions for the velocities parallel to the confining walls for φ < 0.6. According to our criteria, we classify such systems as nonturbulent. Yet, we obtain a Gaussian velocity distribution for φ = 0.68 and (β, λ) = (−5, 4) (Fig. 7). The energy spectrum of that system exhibits a power-law decay with the exponent κ = 5/3, characteristic for Kolmogorov-Kraichnan-type turbulence in the inertial range. Hence, this systems fulfills all the above criteria, and we consider it as fully turbulent.
The slope of the power-law regime depends on the squirmer density. At the smaller packing fraction φ = 0.6 and (β, λ) = (−5, 4), the energy spectrum decreases faster, with the exponent κ = 2. At the same time, the velocity distribution function is non-Gaussian. Hence, the system is not showing active turbulence in the above sense, yet, exhibiting swarming motility. This suggests a tight link between the energy spectrum and the velocity distribution function, a relation which needs further considerations.
As typically observed in turbulent bacterial suspensions [13,33,35], we also obtain a maximum in the energy spectrum at 5 − 10 squirmer lengths, as well as a negative spatial velocity correlation function, in agreement with the presence of a characteristic vortex size.
Inertia of the collective active motion could play an important role, since the crossover from the active ballistic motion -equivalent to inertia of a passive systemto active diffusion appears on the length scale of approximately 6 squirmers lengths, which is comparable to the characteristic vortex size. Yet, the Reynolds number on the scale of a vortex (approximately 10 microswimmer lengths) is still smaller than unity. Here, more detailed theoretical studies of a suitable model are required to assess the relevance of the various interactions on active turbulence.
Despite the similarities of our squirmer systems with bacterial suspensions, there is one major difference, namely, the swimming speed of bacteria increasing in the swarming phase, whereas it decreases in our case [82]. This may point toward a particular role of bacterial flagella in the propulsion of the dense bacterial system.
We like to emphasize that hydrodynamic interactions are paramount for microswimmer swarming and active turbulence, specifically the active stress and the rotlet dipole determine their swarming behavior. However, for Kolmogorov-Kraichnan-type characteristics to merge, in addition, density plays a major role, and ensures an isotropic and homogeneous dynamics on lengths scales larger than approximately a squirmer length. Our simulations provide a benchmark for further theoretical and simulation studies on bacterial turbulence to elucidate the interplay between hydrodynamic stress -specifically a rotlet dipole -, alignment interactions by anisotropic swimmer shapes, and volume exclusion.

Microswimmer model: prolate squirmer
The prescribed surface velocity of the prolate spheroidal squirmer, a homogeneous colloidal particle of mass M , is given by the [58][59][60][61] in terms of spheroidal coordinates τ, ζ, ϕ (1 ≤ τ < ∞, −1 ≤ ζ ≤ 1, 0 ≤ ϕ < 2π) ( Fig. 1(a)) [60,83]. For a squirmer with propulsion direction e = (0, 0, 1), the Cartesian coordinates of a point on the spheroid surface r s = (x s , y s , z s ) T are withr s = x 2 s + y 2 s , r s = |r s |, and τ = τ 0 = b z / b 2 z − b 2 x and the lengths b z and b x along the semimajor and -minor axis ( Fig. 1(a)). The terms with the coefficients B 1 and β (β < 0, pusher) account for swimming in the direction e and an active stress, respectively [60,74,83]. The rotlet-dipole term of strength λ accounts for the torque-free nature of swimming bacteria with a counter rotating cell body compared to the rotating flagellar bundle [63]. The swimming velocity of a squirmer is related to B 1 as To insure quasi-two-dimensional motion between the walls (Fig. 1), a strong repulsive interaction between squirmers and walls is implemented by the truncated and shifted Lennard-Jones potential for y < 2 1/6 σ w and zero else, where y is the closest distance between a wall and the surface of a squirmer. Here, σ w and w determine to the length and energy scale, respectively. Squirmer volume-exclusion interactions are described by a separation-shifted Lennard-Jones potential with parameters σ s and s , where y → d s + σ s in Eq. (10), and d s is the distance between the two closest points on the surfaces of two interacting spheroids [60,83].
The solid-body equations of motion of the squirmers -the center-of-mass translational motion and the rotational motion described by quaternions -are solved by the velocity-Verlet algorithm [60,83].
We apply the MPC approach with angular momentum conservation (MPC-SRD+a) [74,75]. The algorithm proceeds in two steps -streaming and collision. In the streaming step, the MPC point particles of mass m propagate ballistically over a time interval h, denoted as collision time. In the collision step, fluid particles are sorted into the cells of a cubic lattice of lattice constant a defining the collision environment, and their relative velocities, with respect to the center-of-mass velocity of the collision cell, are rotated around a randomly oriented axes by a fixed angle α. The algorithm conserves mass, linear, and angular momentum on the collision-cell level, which implies hydrodynamics on large length and long time scales [55,84]. A random shift of the collision cell lattice is applied at every collision step to ensure Galilean invariance [96]. Thermal fluctuations are intrinsic to the MPC method. A cell-level canonical thermostat (Maxwell-Boltzmann scaling (MBS) thermostat) is applied after every collision step, which maintains the temperature at the desired value [97]. The MPC method is highly parallel and is efficiently implemented on a graphics processing unit (GPU) for a high-performance gain [98].
Squirmer-fluid interactions appear during streaming and collision. While streaming squirmers and fluid particles, fluid particles are reflected at a squirmer's surface by application of the bounce-back rule and addition of the surface velocity u s (8). To minimize slip, phantom particles are added inside of the squirmers, which contribute when collision cells penetrate squirmers. In all cases, the total linear and angular momenta are included in the squirmer dynamics. More details are described in Ref. [60] and the supplementary material of Ref. [83].

Parameters
Multiple squirmers with the semi-major axis b z = 6a and semi-minor axis b x = 2a are distributed in a narrow slit of width L y = 8a, where a is the length of the MPC fluid collision cell. Parallel to the walls, periodic boundary conditions are applied. We set σ w = 1.8a and w = 18k B T . Squirmer propulsion requires fluid particles adjacent to its surface. To avoid MPC particle depletion when two squirmers approach each other, we introduce a safety layer of thickness d v = 0.25a around every squirmer, corresponding to the effective squirmer semi-axes b z +d v and b x +d v , respectively. The squirmersquirmer Lennard-Jones parameters are set to σ s = 0.5a, s = 5k B T . d s (see microswimmer model) is now the distance between two closest points on the surfaces of the two interacting squirmers with effective (larger) semiaxes [57,60,83].
To avoid MPC-particle depletion [57], we employ a high average particle number N c = 60 in a collision cell. Furthermore, we choose a small collision-time step h = 0.02 ma 2 /(k B T ) and the large rotation angle α = 130 • . This results in the fluid viscosity η = 127.8 mk B T /a 4 and the 2D rotational diffusion coefficient around a minor axis D 0 R = 5.2 × 10 −6 k B T a 2 /m. This is in close agreement with the theoretical value of a spheroid D 0 R = 5.5 × 10 −6 k B T a 2 /m. In the far-field, the microswimmer flow field is dominated by the force-dipole term of strength [6,15,62,99] where P = f D l D is the magnitude of the force dipole of force f D and length l D . The latter parameters can be determined from experiments [62] and simulations [63]. The far-field expansion of the flow field of a spheroidal squirmer provides the relation between χ and the active stress parameter β [60]: With the approximation of the bacteria cell body by a spheroid, Eq. (12) provides an estimation of β for a given χ.
In both cases, the viscosity of water is used. These β values approximately fall into the range of active stresses considered in our simulations.

ACKNOWLEDGMENTS
This work has been supported by the DFG priority program SPP 1726 "Microswimmers -from Single Particle Motion to Collective Behaviour". The authors gratefully acknowledge the computing time granted through JARA-HPC on the supercomputer JURECA at Forschungszentrum Jülich.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.

AUTHOR CONTRIBUTIONS
R.G.W. and G.G. designed the study. K.Q. and E.W. wrote the simulation code and K.Q. performed the simulations. K.Q., R.G.W, and G.G.. analyzed and discussed the results. R.G.W, G.G., and K.Q. wrote the paper.