Simulations of flow fields around microscopic objects typically require methods that both solve the Navier–Stokes equations and also include thermal fluctuations. One such method popular in the field of soft-matter physics is the particle-based simulation method of multi-particle collision dynamics (MPCD). However, in contrast to the typically incompressible real fluid, the fluid of the traditional MPCD methods obeys the ideal-gas equation of state. This can be problematic because most fluid properties strongly depend on the fluid density. In a recent article, we proposed an extended MPCD algorithm and derived its non-ideal equation of state and an expression for the viscosity. In the present work, we demonstrate its accuracy and efficiency for the simulations of the flow fields of single squirmers and of the collective dynamics of squirmer rods. We use two exemplary squirmer-rod systems for which we compare the outcome of the extended MPCD method to the well-established MPCD version with an Andersen thermostat. First, we explicitly demonstrate the reduced compressibility of the MPCD fluid in a cluster of squirmer rods. Second, for shorter rods, we show the interesting result that in simulations with the extended MPCD method, dynamic swarms are more pronounced and have a higher polar order. Finally, we present a thorough study of the state diagram of squirmer rods moving in the center plane of a Hele-Shaw geometry. From a small to large aspect ratio and density, we observe a disordered state, dynamic swarms, a single swarm, and a jammed cluster, which we characterize accordingly.

## I. INTRODUCTION

Simulations of hydrodynamic flows in the context of microscopic colloidal and active particles share a long history with the method of multi-particle collision dynamics (MPCD).^{1–3} Different modifications of the original algorithm^{4} have been applied in studies of colloidal suspensions,^{5–9} microswimmers,^{10–23} polymers,^{24–26} and even complex deformable bodies such as red blood cells^{27,28} or the African trypanosome,^{29,30} the causative agent of the sleeping sickness, to name a few examples. Our particular interest lies in the collective dynamics of elongated microswimmers and active filaments for which interesting states, such as swarming,^{21,31–33} vortex formation,^{34–37} active nematic patterns,^{38–40} and active turbulence,^{41–45} have been observed. Many of these features can already be described by continuum models^{46–49} that combine elements of the Toner–Tu^{50} and Swift–Hohenberg^{51} equations. However, it is still a matter of current debate to identify the importance of purely steric and hydrodynamic interactions on the scale of the constituting particles.^{52–57} In this context, we introduced elongated squirmer rods as a model for elongated bacteria using MPCD to simulate their flow fields.^{17} Our model thereby adds to the topmost layer of detail in a wide range of models ranging from coarse-grained interactions in the famous Vicsek model^{58–60} over Langevin dynamics simulation^{41,61–67} with implicit hydrodynamic interactions^{68,69} to explicit hydrodynamic simulations using the MPCD and lattice Boltzmann approach.^{70–73}

In this context, solving the Navier–Stokes equations at vanishing Reynolds numbers with the particle-based mesoscale method of MPCD has few advantages. It inherently includes the omnipresent thermal fluctuations,^{4,74} and complex geometries are straightforward to implement with this method.^{28,75–77} Furthermore, the collision rules of MPCD algorithms suit well for parallel computing on graphic cards so that also on desktop computers, one can perform extensive simulations.^{78,79}

However, a property that most MPCD algorithms share is that the MPCD fluid particles obey the equation of state of an ideal gas.^{4,80,81} This makes the model fluid rather compressible and results in a low speed of sound *c*_{s},^{5} which is in conflict with simulating the dynamics of incompressible solvents at a low Reynolds number. It is a matter of current debate how this affects the dynamics of immersed particles,^{20} especially for strongly inhomogeneous systems, such as clustered microswimmers.^{21} To render the MPCD fluid more incompressible, different strategies can be used. For the ideal-gas models, increasing the density of fluid particles *n*_{0} naturally decreases the compressibility.^{21} However, this comes with an immense increase in simulation time proportional to $n02$ if the system size is held constant. Alternatively, auxiliary field equations can be coupled to the dynamics of the MPCD particles to alter the equation of state of the model.^{81}

In a recent publication,^{82} we introduced an extended version of MPCD with a modified collision step that transfers momentum across a randomly oriented plane within the collision cell in order to generate a non-ideal equation of state. In addition, the collision probability depends on the local density and velocities. This is different from traditional MPCD algorithms, where collisions occur at a fixed rate. Our approach is motivated by the collision rate in the Boltzmann equation^{83} and has the advantage that it can easily be integrated into an existing MPCD program. A similar approach has been formulated for two-dimensional systems in Refs. 84 and 85. While we have already derived expressions for the equation of state and the viscosity in our previous publication,^{82} in the present work, we demonstrate its application to the dynamics of single squirmers and squirmer rods, but most importantly to explore their collective dynamics.^{17}

In Sec. II, we start by shortly reviewing the extended MPCD method^{82} and the squirmer-rod model for elongated microswimmers^{17} introduced by us in earlier publications. We then demonstrate the applicability of the extended MPCD method and, in particular, compare it to the traditional MPCD version. First, in Sec. III, we show that our new method accurately reproduces the leading hydrodynamic multipoles contributing to the flow field of a single squirmer in the bulk fluid and in confinement. Second, in Sec. IV, we explicitly demonstrate the lower compressibility of the extended MPCD method considering the specific example of a cluster of squirmer rods. Then, for shorter rods, we illustrate quantitative differences in their swarming behavior between both MPCD methods. Finally, in Sec. V, we employ the extended MPCD method to study in-depth the collective dynamics of squirmer rods. We present a state diagram of squirmer rods for variable aspect ratio and density, which includes dynamic swarms, single swarms, and jammed clusters, and we discuss their structural and dynamic properties. We finish with concluding remarks and an outlook in Sec. VI.

## II. METHODS

### A. Extended MPCD method

As common to all MPCD algorithms, the fluid is modeled by point-like particles, the dynamics of which proceeds in two alternating steps: ballistic motion along the mean free path followed by a collision. Because the collision rule conserves linear and also angular momentum, the simulated fluid flow solves the Navier–Stokes equations.^{4} In MPCD, these steps occur simultaneously for all particles. During the streaming step (i), particles with masses *m*_{0}, positions **x**_{i}(*t*), and velocities **v**_{i}(*t*) move ballistically for the duration of the time step Δ*t*,

Moreover, the particles collide with bounding walls, colloids, or, in our case, squirmers^{3,13,14,16–18} and transfer linear and angular momenta. To enforce either the no-slip boundary condition or a predefined slip-velocity at a squirmer surface, the so-called bounce-back rule is applied.^{4,75}

Our extended MPCD method^{82} now modifies the collision step (ii). As usual, the simulation volume is divided by a cubic lattice with lattice constant *a*_{0} and all particles are grouped into cubic cells with volume $V\xi $ centered around position ** ξ**. Before each collision step, the lattice is randomly shifted against the previous lattice to ensure Galilean invariance.

^{86}The

*n*

_{ξ}particles inside a cell have the mean velocity

**v**

_{ξ}, center-of-mass position

**x**

_{ξ}, and moment-of-inertia tensor

**I**

_{ξ}defined relative to

**x**

_{ξ}.

Now, in our modified collision step, each cell is divided into two halves *A* and *B* by a plane $Pn\u0302,x\xi $ that contains the center of mass **x**_{ξ} and has the normal unit vector $n\u0302$, which, by definition, always points toward *A*. For each cell and time step, $n\u0302$ is randomly drawn from a discrete set of 13 possible orientations. The modified collision step transfers momentum across the randomly oriented plane and thereby generates a non-linear equation of state. We denote the respective number of particles on each side of the plane by *n*_{A} and *n*_{B} and their respective mean velocity components along $n\u0302$ taken relative to **v**_{ξ} by $v\u0304A$ and $v\u0304B$. We also introduce the relative velocity $\Delta u\u2261v\u0304B\u2212v\u0304A$ for further use. The collision step introduces for each particle velocity **v**_{i} in cell $V\xi $ a new value $vinew$ and is defined by

where **x**_{i,c} denotes the position vector of particle *i* relative to the center-of-mass position **x**_{ξ}. In contrast to other MPCD methods, the collision only occurs with a certain probability as explained below. To initiate a collision, the stochastic variable *χ*(Δ*u*) is set to one, otherwise it is zero. The term in square brackets in the first line sets the normal velocity component of particle *i* relative to $n\u0302\u22c5v\xi $ to a new random value *δv*_{i}, which we introduce below. It mainly transfers momentum across the collision plane $Pn\u0302,x\xi $, and all values add up to zero in order to preserve the total momentum. Finally, the term in the second line is added to preserve the angular momentum $L\xi =m0\u2211xi\u2208V\xi xi\xd7vi$ during collision.

Note that the momentum transfer across the randomly oriented plane is responsible for the decreased compressibility of the MPCD fluid. As we showed in Ref. 82, it generates a term proportional to $n\xi 2$ in the equation of state, which thereby deviates from the ideal-gas equation of the classical MPCD fluid. On average, a collision between particles in regions *A* and *B* occurs when the two clouds move toward each other meaning $\Delta u\u2261v\u0304B\u2212v\u0304A>0$. Following the collision term in the Boltzmann equation, we take for the collision rate *c* Δ*u n*_{A} *n*_{B}, where *c* is the scattering cross section.^{83} Thus, we set the probability for a collision to occur to

where Θ(Δ*u*) is the Heaviside step function, and we use the exponential function to guarantee *p*_{χ}(Δ*u*) ≤ 1.

Finally, we introduce the new velocity component along the normal $n\u0302$,

Here, the first term transfers momentum between regions *A* and *B*. The first index applies to particles *i* in region *A*. They take over the momentum $m0v\u0304B$ from region *B* and vice versa. The prefactor *n*_{B/A}/*n*_{A/B} guarantees that the overall momentum is conserved. Essentially, the total momenta from regions A and B are just swapped. The second term assigns each particle a thermal velocity $\delta viMB$ drawn from a Maxwell–Boltzmann distribution at temperature *T*, which thereby serves as a thermostat for the fluid. The last term subtracts the mean thermal velocities in each region (*A* and *B*),

in order to conserve momentum.

As in all MPCD algorithms, the so-called “ghost” particles are introduced to participate in the collision step. They are added where overlap with boundaries leads to partially empty collision cells, thereby improving the validity of the no-slip boundary condition.^{75} After the collision step, the changes in the linear and angular momentum of the ghost particles are assigned to movable objects, which in our case are squirmers. This ensures that the linear and angular momenta are conserved.

### B. Squirmer and squirmer rods

To investigate the flow fields of microswimmers and their hydrodynamic interactions in collective motion, we use spherical squirmers and build squirmer rods from them. Spherical squirmers propel by an imposed axisymmetric slip-velocity field on their surfaces,^{3,87,88}

which initiates flow in the surrounding fluid. Here, the unit vector $x\u0302s$ points from the squirmer center toward the squirmer surface and $e\u0302$ indicates the orientation and swimming direction of the squirmer. The resulting flow field in the bulk fluid is a superposition of multipole solutions of the Stokes equations. In particular, it contains a source dipole singularity with strength $B1=B1sR3$ and a force dipole singularity with strength $A2=\beta B1sR2$. The latter is controlled by the parameter *β*, which determines the squirmer type, and *R* denotes the radius of the squirmer. For different values of *β*, the squirmer represents either a pusher (*β* < 0), a neutral swimmer (*β* = 0), or a puller (*β* > 0). The squirmer parameter $B1s$ also determines the swimming speed $v0=2/3B1s$.

To model the actual rod shape of many biological microswimmers, we use the squirmer rod. It is constructed by arranging several neutral spherical squirmers on a line such that they overlap and form a single rigid body. By varying the squirmer distance *d*, rods of different aspect ratios *α* = *l*_{S}/2*R* are realized where *l*_{S} is the rod length. In this work, we always use *N*_{sq} = 10 spherical squirmers and a maximum squirmer distance of *d* ≈ 0.8*R* so that the surface is still smooth enough that swimmers can slide past each other. This amounts to a maximum aspect ratio of *α* = 4.75. To model steric interactions between pairs of squirmer rods, we employ the Weeks–Chandler–Andersen potential^{89} that acts between two squirmers from different rods. The force constant is chosen strong enough (*ε*_{WCA} ≈ 10^{4}*k*_{B}*T*) to ensure that there is no significant overlap between two rods.

The squirmer rod model described so far resembles ciliated micro-organisms such as *Paramecium*. To model bacteria, such as *E. Coli*, which are driven by a bundle of rotating flagella extending from its rear, one can concentrate the surface velocity on the back of the rod by multiplying the slip velocity with an envelope or step function.^{17}

### C. Parameters

We choose the parameters for the extended MPCD method such that the total viscosity $\eta =\eta str+\eta coll+A\u224816m0kBT/a02$ (see Ref. 82) given in MPCD units matches previous studies where other MPCD methods are employed.^{13,14,17,21} In particular, the collision parameter *c* = 1/100 and the average number density *n*_{0} = 20 are chosen such that the average collision rate ⟨*χ*(Δ*u*)⟩ ≈ 0.14 is low enough to be in the quadratic regime of the equation of state.^{82} Together with a time step of $\Delta t=0.005a0kBT/m0$, the total viscosity takes the prescribed value. For the squirmers, we use a radius of *R* = 3*a*_{0} and a source dipole strength of $B1s=0.1$, which results in a Péclet number of Pe = *Rv*_{0}/*D*_{T} ≈ 354.6, where *D*_{T} is the translational diffusion coefficient of the spherical squirmer.^{14}

To initialize a simulation, we distribute the fluid particles randomly and choose thermal velocities from the Maxwell–Boltzmann distribution. In addition, the squirmer rods are randomly distributed in the midplane of the simulation box for area fractions *ϕ* below 0.6. For *ϕ* ≥ 0.6, we initiate the simulation by arranging the swimmers in a rectangular pattern and by choosing their orientations at random along the longer side of the rectangular unit cell.

We use the following geometries and simulation times in our computer studies. To simulate the flow fields of pusher-type squirmers, we use periodic boundary conditions for the bulk fluid with an edge length of *L* = 150*a*_{0} of the cubic simulation box. In the Hele-Shaw geometry, the confining parallel plates have a separation of Δ*z* = 8*a*_{0} and an edge length of *L* = 300*a*_{0}, and the squirmers are constrained to move in the midplane of the cell (see Fig. 1) by setting the normal component of the velocity *v*_{z} and the components *ω*_{x}, *ω*_{y} of the angular velocities of the squirmer rods to zero after each time step. The flow fields are first equilibrated during 50 000 time steps Δ*t*, and 20 independent simulation runs are used to average the flow fields over a total of 10^{7}Δ*t*. To simulate the collective motion of squirmer rods, we choose a Hele-Shaw geometry with a plate separation of Δ*z* = 18*a*_{0}, otherwise the edge length *L* = 300*a*_{0} is the same. The system is equilibrated during 10^{7} time steps Δ*t*, and data are sampled during another 10^{7}Δ*t*.

## III. FLOW FIELDS IN THE EXTENDED MPCD METHOD

In this section, we present a first test for the extended MPCD method by demonstrating the accuracy with which it simulates flow fields of microswimmers. We compare the simulation results of the flow field of a pusher-type squirmer with the analytic multipole expression both in the bulk fluid and in confinement between two parallel plates.

We briefly outline the underlying theory. In the Stokes regime, the flow field generated by an object moving through the bulk fluid can be written as a series of hydrodynamic multipoles. This series contains at the lowest order the Stokeslet or hydrodynamic force monopole, which is not allowed for autonomous swimmers, the source dipole, and their higher-order derivatives. For swimmers with radial symmetry about their swimming direction $e\u0302$, the flow fields of these multipoles can be expanded into a Legendre polynomial *P*_{n}(cos *θ*). In spherical coordinates *r*, *ϕ*, and *θ* with $z\u0302=e\u0302$, the radial velocity component of the multipole expansion reads

We use the orthogonality of the Legendre polynomials to extract the coefficients *u*_{r,n}(*r*) = *A*_{n}*r*^{−n} + *B*_{n}*r*^{−n−2} from the simulated flow field with the radial component *u*_{r}(*r*, *θ*) and then determine the coefficients of the force multipoles *A*_{n} and source multipoles *B*_{n} by a polynomial fit in *r*^{−n} and *r*^{−n−2}.

When the object moves in a fluid confined between two parallel plates with distance Δ*z* (Hele-Shaw geometry), a multipole expansion for the flow field averaged along the plate normal exists.^{17,90} Using polar coordinates *ρ* and *φ* in the plane, the radial velocity component reads

where *φ* is measured against the swimming direction $e\u0302$. Along the normal or *z* direction, the flow field is parabolic, which is integrated out and included into the coefficients $An$ and $Bn$. Similar to the bulk case, we use the orthogonality of the basis functions cos(*nφ*) to extract the coefficients $u\u0303\rho ,n(\rho )=(An+Bn)\rho \u2212n\u22121$. The relative strength of the amplitudes $An$ and $Bn$ can be deduced either by extrapolating from the bulk case or by their different scaling behavior $An\u221d\Delta z$ and $Bn\u221d1/\Delta z$. Most interesting about Eq. (8) is the radial decay of the different multipoles compared to the bulk case. For a force-free microswimmer ($A1=A1=0$), the source dipole moment $B1$ becomes the dominant term in the Hele-Shaw geometry, while in the bulk fluid the force dipole *A*_{2} dominates.

In Fig. 2, we demonstrate the multipole analysis of the flow field of the pusher squirmer simulated with the extended MPCD method. Figure 2(a) shows all non-vanishing radial velocity coefficients *u*_{r,1}(*r*) and *u*_{r,2}(*r*) in the bulk fluid for the squirmer type with *β* = −1. In agreement with the analytic expression, each follows a distinctive power law specific for the multipole moment. The radial decay of *u*_{r,1}(*r*) clearly has the exponent −3 corresponding to the source dipole moment *B*_{1}, and *u*_{r,2}(*r*) decays with an exponent −2 in the farfield as expected for the force dipole moment *A*_{2}. However, for *u*_{r,2}(*r*) close to the squirmer, we also recognize the importance of the short-ranged flow field of the source quadrupole moment *B*_{2}, which has the exponent −4. Thus, the polynomial fit deviates from the straight line. The additional source quadrupole moment with a negative value is necessary to balance the force dipole *A*_{2} such that the correct surface velocity field of the pusher squirmer is implemented. In contrast, the flow field of the source dipole moment *B*_{1} obeys the boundary condition by itself so that *u*_{r,1}(*r*) follows a pure power law.

In the Hele-Shaw geometry, we choose *β* = −3 to get a stronger signal from the force dipole. We observe the corresponding radial coefficients $u\u0303\rho ,1(\rho )$ and $u\u0303\rho ,2(\rho )$, which perfectly show the expected radial decay [see Fig. 2(b)]. However, note that due to the confinement, the exponents of the moments of force and source dipole are swapped, and hence, the source dipole ultimately dominates the far field. Note that, unlike in the bulk fluid, force and source multipoles in $u\u0303\rho ,n(\rho )$ decay with the same power law in a Hele-Shaw geometry. Here, the deviation of the force dipole field from theory at *ρ* < Δ*z* (gray dashed line) is due to short-ranged contributions to the flow field that are neglected in Eq. (8).^{90,91}

## IV. COMPARISON OF TRADITIONAL AND EXTENDED MPCD

### A. Fluid density in squirmer-rod systems

In this section, we demonstrate how the low compressibility of the extended MPCD algorithm results in a much more uniform density of the MPCD particles compared to traditional MPCD methods in systems where the fluid is strongly forced. In particular, we consider a dense suspension of squirmer rods that are constrained to move in the center plane between two confining parallel plates, which we will further study in Sec. V. We choose squirmer rods with an aspect ratio of *α* = 4.0 at an area fraction of *ϕ* = 0.6.

As a reference, we also perform simulations with the standard MPCD-AT+A algorithm with angular momentum conservation and an Andersen thermostat.^{3} We choose the simulation parameters such that the relevant physical quantities characterizing the system, the fluid viscosity and the Péclet number, match in both MPCD algorithms. Hence, for the standard MPCD method, we use *n*_{0}′ = 10 and $\Delta t\u2032=0.02a0kBT/m0$, which results in the same viscosity $\eta \u224816m0kBT/a02$.^{14,82} Because the two MPCD algorithms use different time steps, we choose 20 × 10^{6} and 5 × 10^{6} time steps, respectively, to simulate the same period with the extended and conventional MPCD method.

With both algorithms, the initial random distribution of squirmer rods evolves into a stable compact swarm, which moves through the system [insets of Figs. 3(a) and 3(b)]. Above and below the swarm, the fluid experiences strong forcing due to the slip velocity fields on the rod surfaces. This generates a strong pressure gradient between the region of the cluster and the rest of the system. However, the MPCD fluid densities simulated with the traditional and extended MPCD algorithms behave differently, as shown in the main parts of Figs. 3(a) and 3(b). While the relative fluid densities $\u27e8\rho /\rho 0\u27e9z$ and $\u27e8\rho /\rho 0\u2032\u27e9z$, averaged over the *z* direction, become non-uniform in both simulations, the deviations from the mean density when using the extended MPCD algorithm are much smaller.

For the traditional MPCD algorithm [cf. Fig. 3(a)], the pressure gradient creates regions sparse of fluid particles above and below the swarming cluster as well as between the rods. These fluid particles are pressed into regions outside of the cluster, where the density increases. Note that the density field also includes the ghost particles of the squirmer rods, which are clearly visible in white because their density is held constant at *n*_{0} by definition. In the vicinity of the swarming cluster, the relative density drops to values of about 60% and it rises to 140% in regions outside the cluster. In simulations with stronger confinement, i.e., smaller Δ*z*, this effect becomes even stronger.

In contrast, for the extended MPCD algorithm, the deviations in relative density are significantly lower with around ±6% [cf. Fig. 3(b)]. Thus, the density is only mildly decreased in the vicinity of the cluster and increased in the remaining regions of the system. Note that the fluid density directly between squirmer rods in the midplane follows the trend of the height-averaged density in Figs. 3(a) and 3(b), but the deviations are stronger. With the traditional MPCD method, the relative density drops to 30%, while with the extended MPCD method, the minimum relative density is 90%. All this is in accordance with the low compressibility of the fluid in the extended MPCD method, which is about five times lower compared to the simulations with traditional MPCD.^{82} However, note that this comes with an increase in simulation time.

Thus, the extended method better satisfies the incompressibility condition ∇ · **v** = 0. Furthermore, also the viscosity, which depends on the local density, is more uniform in the extended MPCD method and the simulation thus closer to the Stokes limit. This is particularly important for simulating microswimmers, where typically a constant Péclet number Pe ∝ *η* is assumed for the entire system.

To be more quantitative, we show in Fig. 4 the pressure $\u27e8p/p0\u27e9z$ and viscosity fields $\u27e8\eta /\eta 0\u27e9z$ derived from the density fields shown in Fig. 3 using the respective analytical formulas of the traditional method and our extended MPCD method.^{2,82} Note that, since both pressure and viscosity depend on density *ρ*, they are not independent of each other. While the viscosity of the traditional MPCD method is linear in the density and varies accordingly by ±40% [cf. Fig. 4(c)], variations are much smaller for the extended method. Here, the relation between density and viscosity is more complex and only results in a variation of ±10% [cf. Fig. 4(d)]. The pressure fields show a similar behavior. A closer inspection shows that the deviations from the mean value are marginally smaller compared to the viscosity variations, as one can see in Fig. 4(a) for the traditional method and Fig. 4(b) for the extended method.

Now, with the significant increase in the local Péclet number Pe ∝ *η* outside of the cluster for the traditional MPCD method, we would assume enhanced clustering compared to the much less compressible fluid, which has been reported for spherical squirmers.^{21} However, clustering in the two simulations reported here is very similar, suggesting that steric interactions between the elongated rods dominate. In Subsection IV B, we will therefore compare two simulations with shorter rods such that the alignment due to the elongated shape is less pronounced and hydrodynamic interactions between the rods become more important.

### B. Influence on hydrodynamic interactions

To investigate how the lower compressibility of the extended MPCD method affects the hydrodynamic interactions between the rods and ultimately their clustering, we choose a lower aspect ratio of *α* = 2.5 and an area fraction of *ϕ* = 0.4 to reduce the effect of steric interactions. Otherwise, we use the same system parameters as in Sec. IV A and perform simulations with the two different MPCD algorithms.

In the temporal evolution of both systems, small clusters form due to direct collisions. However, only with the extended MPCD method, they begin to form swarms that grow to a significant size, while with the traditional MPCD method the collision clusters remain small and dissolve rather quickly. The snapshots of Figs. 5(a) and 5(b) as well as videos V1 and V2 illustrate this difference. To quantify it, we identify individual rod clusters such that all spherical squirmers that are separated by a center-to-center distance less than Δ*x* < 1.02*R* belong to one cluster. The cluster size *N*_{c} then is the number of rods within a cluster determined by the ratio of squirmer number in the cluster to squirmers per rod. Additionally, we define the cluster velocity as the instantaneous average of rod velocities over the cluster.

Figure 5(c) shows the mean cluster velocity ⟨*v*_{c}⟩ relative to the squirmer velocity *v*_{0}, and Fig. 5(d) shows the probability density for a rod to lie in a cluster with size *N*_{c}. Note that the swimming speed of single squirmer rods is ∼1.1*v*_{0}.^{17} The difference between both methods is immediately visible. For small groups *N* < 10, we observe similar velocities *v*_{c} ≈ 1.1*v*_{0} in Fig. 5(c) close to the value of a single rod in the bulk fluid. These rods form a dilute background, which exists equally in both systems, as the identical probability distributions of cluster sizes with *N*_{c} < 10 in Fig. 5(d) show.

Large clusters also exist in both systems, but they behave differently. Visual inspection of video V1 shows that clusters in simulations with traditional MPCD are formed rather by collisions between rods or other clusters and dissolve quickly. As a result, the mean cluster velocity *v*_{c} in Fig. 5(c) decreases drastically with the cluster size *N*_{c}. Furthermore, Fig. 5(d) shows that the clusters do not exceed a size of roughly *N*_{c} = 200.

In contrast, the formation of stable swarms in simulations with the extended MPCD method leads to the plateau value of the mean cluster velocity at *v*_{c} ≈ 0.6*v*_{0} for *N*_{c} > 30 [Fig. 5(c)]. In addition, the distribution *P*(*N*_{c}) in Fig. 5(d) shows the existence of larger clusters. Both results indicate that these swarms are more ordered compared to the simulations with the traditional MPCD method. They accumulate new members during the collective swimming so that the polar order remains constant. Visual inspection of video V2 shows a dynamic swarming behavior where swarms still merge and dissolve, but overall they are more stable.

This is an interesting result because it shows that the extended MPCD method with its lower compressibility stabilizes clusters in the current system. Moreover, we conclude that the extended method mediates hydrodynamic interactions better such that rods align and form swarms.

## V. COLLECTIVE DYNAMICS OF SQUIRMER RODS

We now present a thorough investigation of the collective dynamics of squirmer rods varying aspect ratio *α* and 2D area fraction *ϕ*. Overall, we observe four different dynamic states, which we briefly describe here based on visual inspections before turning to a detailed analysis in Subsections V A and V B. The state diagram together with exemplary snapshots is given in Fig. 6, and videos V3–V6 in the supplementary material illustrate the different dynamic states and correspond to the respective snapshots [Figs. 6(b)–6(e)].

At small densities and also at larger densities below a critical value of the aspect ratio, *α*_{min} ≈ 2, we observe a disordered state with little variations in the local density. Thereby, we support the conclusion of Ref. 21 that the low fluid compressibility reduces clustering. Increasing *α* above 2 and also for sufficiently large density, swarms of squirmer rods form. They consist of clusters with inherent polar order such that the squirmer rods move collectively in one direction. However, as described already in Sec. IV B, the swarms are dynamic, meaning that they are not stable in time but frequently break apart and reform creating an interesting dynamics reminiscent of dynamic clustering found for spherical active particles in the presence of a chemical field.^{92,93} Crossing the dashed black line in the state diagram by increasing density and/or aspect ratio further, the dynamic swarms merge into one large swarm, which is stable in time. Depending on the individual history of such a swarm, it can show either continuous translational or rotational motion. Further increasing area fraction *ϕ*, the global swarms grow bigger until they span over the whole system. In this situation, the cluster is jammed, where all movements inside the cluster are blocked and only some drift velocity of the whole cluster remains. Our findings generally show parallels to the surface dynamics of *B. subtilis* reported in Ref. 94. In particular, the dynamic swarms, single swarm, and the jammed cluster are observed. For systems of elongated spheroidal squirmers, the disordered and single swarm states have been reported as well.^{21}

For a detailed quantitative analysis of these dynamic states, we first study in Sec. V A the distribution of the local area fraction *ϕ*_{loc} as well as the cluster size and velocity. To analyze the alignment of the squirmer rods in the dilute and also cluster states, we use the spatial orientational pair-correlation function in Sec. V B.

### A. Distributions of local area fraction, cluster size, and velocity

The observed swarms are clustered squirmer rods with additional polar order. Thus, to differentiate swarms from dilute regions in the system and to characterize the different states, we determine the distribution of local area fraction as well as cluster sizes and their velocities. To assign each squirmer rod a local area fraction *ϕ*_{loc}, we use a Voronoi tessellation following our earlier work.^{14} The concept is extended to the squirmer rods as follows and illustrated in Fig. 7(b). First, the ordinary Voronoi tessellation is calculated for all spherical squirmers. Here, periodic images are added around the central simulation box to achieve correct results across the periodic boundary conditions. Then, the Voronoi cells of the squirmers that belong to the same rod form the Voronoi cell of the rod, which are no longer convex polygons. Finally, the local area fraction follows from the ratio of the squirmer rod’s cross section and the area of the Voronoi cell. In addition, we identify individual rod clusters and determine their velocity as described in Sec. IV B. To calculate the distribution of local area fractions *P*(*ϕ*_{loc}) and to identify clusters, we use all snapshots of the system after an equilibration time of *t*_{eq} = 10 × 10^{6}Δ*t* in steps of 5 × 10^{4}Δ*t*.

We now characterize the different dynamic states using Fig. 7(a), where we plot the normalized distribution functions *P*(*ϕ*_{loc}) for different *global* or mean area fractions *ϕ*, always for squirmer rods with an aspect ratio of *α* = 3.25. In addition, Fig. 7(c) shows individual clusters with their size *N*_{c} and velocity *v*_{c} so that for each area fraction a dynamic state is represented by a cloud of points. A dashed line is added for each area fraction *ϕ*, indicating the maximum possible cluster size equal to the total number of rods. Note that the distribution in the velocity is naturally broader for small clusters, which give rise to larger fluctuations.

Increasing *ϕ* from small to large values, we observe the following characteristics. At *ϕ* = 0.1, the distribution *P*(*ϕ*_{loc}) is narrow and strongly peaked at the global value, which indicates the absence of pronounced clustering in the system. The distribution of cluster sizes in Fig. 7(c) reflects this behavior indicating only small groups of colliding squirmer rods. This is where we locate the disordered state, which is homogeneous in the local area fraction. At the larger area fraction of *ϕ* = 0.2, *P*(*ϕ*_{loc}) shows the same trend but now the peak at low *ϕ*_{loc} is broader and a pronounced tail toward higher values of *ϕ*_{loc} is visible. Most of the system is still dilute, but small denser clusters are present, as shown in Fig. 7(c). This indicates that we are at the transition to the state of dynamic swarms. However, Fig. 7(c) indicates that these clusters are still mainly caused by collisions and not by longer lasting swarms with intrinsic rod alignment because their velocities are low.

At *ϕ* = 0.4, we are right in the subsequent state of dynamic swarms. Here, the distribution *P*(*ϕ*_{loc}) becomes very broad and, in particular, bimodal. The bimodal distribution of area fractions has also been found for swarming *B. Subtilis*.^{94} This signature is caused by parts of the system forming a dilute disordered background, while also dense dynamic swarms are visible. Their preferred area fraction indicated by the right peak of *P*(*ϕ*_{loc}) is close to the maximum *ϕ*_{loc} = 1. Thus, the clusters in the swarms are compact. Note that the peak for the dilute disordered background becomes more pronounced when decreasing the global area fraction and/or the aspect ratio toward the disordered state. The very broad distribution of cluster sizes and velocities in Fig. 7(c) also reflects this dynamic state of swarms. However, the size of the clusters always remains well below the total number of squirmer rods.

At the global aspect ratio of *ϕ* = 0.6, we are in the state where the squirmer rods form a single moving swarm. Here, the distribution *P*(*ϕ*_{loc}) only shows one peak close to *ϕ*_{loc} ≈ 1 and the dilute background is nearly absent. The distribution of cluster sizes in Fig. 7(c) is narrow, meaning that this cluster is highly stable and its size almost reaches the total number of rods in the system. The jammed cluster at the highest aspect ratio *ϕ* = 0.77 shows almost the same trend. However, the maximum of *P*(*ϕ*_{loc}) slightly moved to the left and the distribution is broader compared to *ϕ* = 0.6. This is in agreement with visual inspections, which show that the jammed clusters are less dense. When they form, the density is so high that they cannot all align properly as in the swarms, which we will confirm in Sec. V B. As shown in Fig. 7(c), the size of the cluster now reaches the total number of squirmer rods and the velocity of the jammed cluster is significantly smaller compared to the single swarm state.

For squirmer rods with other aspect ratios *α*, we observe the same behavior, with differences only in the quantitative description. At smaller aspect ratios *α* < 3.25, the distributions *P*(*ϕ*_{loc}) are broader and the bimodality in the swarming state is less pronounced. For higher aspect ratios *α* > 3.25, the trend is inverted and also the peaks of *P*(*ϕ*_{loc}) are more pronounced. Equally, the distribution of cluster sizes becomes more compact. This is a clear consequence of the fact that the rods tend to align better at larger aspect ratios *α*, which stabilizes the swarms.

### B. Orientational correlations

As a second measure to characterize the different states of our system, we use the orientational pair-correlation function $\u27e8e\u0302i\u22c5e\u0302j\u27e9$, where $e\u0302i$ and $e\u0302j$ are the orientation vectors of two squirmer rods with a distance Δ*x*. It helps us to distinguish the different swarming and jamming states by quantifying the degree of alignment in the clusters. Figure 8(a) presents the orientational pair-correlation function for different global area fractions for squirmer rods with an aspect ratio of *α* = 3.25. The plot uses the same parameters as Figs. 7(a) and 7(c) in Sec. V A and supplements the results reported there.

For dilute rod suspensions with an area fraction of *ϕ* = 0.1, we observe that the squirmer rods are indeed disordered. Although the orientational pair correlations have a pronounced peak at the side-by-side contact distance of two squirmer rods, they quickly decay to almost zero at two rod lengths. The strong correlations at a short distance make sense since the squirmer rods collide with each other even in the disordered state. Because of their active movement, pairs of squirmer rods are more likely to approach and collide if their swimming directions are already nearly aligned to each other. Moreover, a side-by-side configuration is more persistent, which causes the first peak with $\u27e8e\u0302i\u22c5e\u0302j\u27e9\u22481$. At Δ*x* = *l*_{S}, we observe a second smaller peak in the orientational pair correlations, which is due to two squirmer rods swimming behind each other. This configuration persists for longer times than a head on collision, and therefore, the correlation function is positive. At *ϕ* = 0.2, we are close to the dynamic-swarm state. Still, the pair correlations are similar to *ϕ* = 0.1, but overall they are increased. While the height of the first peak is almost the same, in particular, the second peak is more pronounced. This indicates the presence of small clusters with squirmer rods swimming side by side and also behind each other.

At *ϕ* = 0.4, the orientational pair-correlations clearly change. They are of longer range, and the system is in the dynamic-swarm state. For example, the correlation only decays slowly to $\u27e8e\u0302i\u22c5e\u0302j\u27e9\u22480.25$ at 4*l*_{S}. The peak at Δ*x* = *l*_{S} for two squirmer rods swimming behind each other increases, and a careful inspection of the curve also shows an additional peak at Δ*x* = 0.62, twice the side-by-side distance. This agrees with visual inspections, which show that squirmer rods tend to swim in successive rows within a swarm. Furthermore, the very weak peaks at 2*l*_{S} and 3*l*_{S} agree with the approximate extension of *R*_{S} ≈ 3*l*_{S} of the dynamic swarms. Note that the alignment is not perfect within clusters, and hence, the correlation reduces with increasing distance.

For the higher area fraction *ϕ* = 0.6, the squirmer rods now form a single big swarm. Here, we see a further increase in the strength and range of the orientational correlations. Interestingly, Fig. 8(b) shows that the shape of the correlation function depends on the initial conditions. As described in Sec. II C, for *ϕ* ≥ 0.6, we only change the orientation randomly along one direction. As a result, in each simulation run, the cluster assumes a different shape and behaves differently. In particular, in the first simulation run the cluster of squirmer rods rotates rather than performing a straight motion as observed in simulation run 2. Thus, the rods at opposite sides of the cluster are anti-parallel to each other, causing the negative correlations $\u27e8e\u0302i\u22c5e\u0302j\u27e9$ at distances Δ*x* > 6*l*_{S}. In simulation runs 3 and 4, the configuration of the cluster creates both straight motion and rotation so that it moves on a circle. Thus, the orientational pair-correlation function has an intermediate characteristic.

Increasing the area fraction further to *ϕ* = 0.77, the cluster is jammed. Compared to the swarming cluster, the orientational correlations are weaker and of shorter range. This indicates that both local and global alignments are reduced and swarming does not occur. The system is too dense to allow alignment during swimming, which initiated the swarm formation at lower densities. For *α* = 3.25, the jammed cluster is not entirely stable and can partially break up and reform again, as observed by visual inspection. As a result, the orientational correlation function looks quite similar to that for the swarm state at *ϕ* = 0.4. For the larger aspect ratios *α* = 4.0 and *α* = 4.75, the jammed cluster is entirely stable and the orientational correlation function decays to zero at even shorter distances.

For aspect ratios different from *α* = 3.25, the results are qualitatively similar. For smaller aspect ratios, the orientational correlations are expectably lower and decay faster than for higher aspect ratios. For *α* > 3.25, the orientational correlations are stronger at short distances Δ*x* < 2*l*_{S}, especially in the dynamic-swarm state. However, in the single-swarm state, the squirmer rods cannot align as much as before, which results in lower cluster velocities.

## VI. CONCLUSIONS AND OUTLOOK

In this article, we have tested the extended MPCD method introduced in our previous article^{82} against known cases and also used it to explore the collective dynamics of squirmer rods. As we have demonstrated in Sec. III, the extended MPCD method accurately solves the Navier–Stokes equations in the Stokes limit. In particular, the hydrodynamic multipoles that occur due to the prescribed surface flow field of the pusher-type squirmer are correctly represented in the simulated velocity field. In both, the bulk fluid and Hele-Shaw geometry, the hydrodynamic source dipole and force dipole decay with the correct power law. In the simulations of the bulk fluid, we could also identify the source quadrupole, which becomes noticeable close to the squirmer.

Due to the non-linear equation-of-state, the MPCD fluid of the extended method is less compressible and, therefore, the density, viscosity, and pressure are significantly more uniform compared to traditional MPCD methods. Small variations in the fluid density remain since also in the extended model the fluid density encodes fluid pressure. Nevertheless, squirmer rods with a high aspect ratio *α* = 4 show very similar clustering in both methods, which suggests that steric interactions due to the strongly elongated shape are dominant. To reduce the tendency to align and thereby increase the influence of hydrodynamic interactions, we lowered the aspect ratio to *α* = 2.5. Interestingly, the observed swarming of dynamic clusters is more pronounced with the extended MPCD method. They have higher polar order, are more stable, and form larger clusters. Our understanding is that in a less compressible fluid where the fluid density is not reduced between nearby rods as in the traditional MPCD, the hydrodynamic flow reduces the occurrence of hard collisions, and thereby, the rods have the possibility to better align along each other. This favors the formation of swarms due to the prolonged duration of polar encounters.

Finally, we systematically investigated the collective dynamics of squirmer rods and found four different dynamic states with increasing aspect ratio *α* and area fraction *ϕ*, namely, a disordered state, dynamic swarms, a single swarm, and a jammed cluster. Similar states are identified for *B. Subtilis* in the experiments of Ref. 94. The disordered state and a global swarm also exist for spheroidal squirmers.^{21} To characterize the dynamic states quantitatively, we determined the distribution of local area fraction *P*(*ϕ*_{loc}) as a measure for the degree of clustering in the system. Additionally, we identified clusters of squirmer rods and determined their mean velocities, while the intrinsic structure of the swarms was analyzed with the orientational pair-correlation function.

For low density *ϕ* and aspect ratio *α*, we observe a disordered state with homogeneous density in which squirmer rods only interact during frequent collisions. Above a critical value of *α* ≈ 2 and at larger densities, the squirmer rods form dynamic swarms, where the system clearly separates into a dilute background and very dynamic polar clusters with a broad size distribution. When the density is increased further, the dynamic swarms turn into a single stable swarm. The dilute background of squirmer rods disappears, and the swarm keeps a high polarization. Depending on the initial condition, the swarm can exhibit directed motion, rotation, or a combination of both, which results in a circular trajectory. At the highest density *ϕ* = 0.77 in our simulations, the system is so crowded that the formation of highly ordered swarms is inhibited. The rod cluster becomes jammed and moves with a significantly lower velocity compared to the swarms. Our findings show similarities with previous studies^{64–67} of dry active rods with different penetrability. In particular, the disordered and also motile cluster or swarming state is found in all systems. These states appear to originate from the rod shape and frequent collisions rather than the particular details of the steric interactions. Additionally, Refs. 64–67 found states of giant clusters of different nature. If particles cannot cross each other, also immobile giant clusters are observed.^{65,67} In contrast, the giant clusters in our study are more polar and therefore more motile, which we attribute to the hydrodynamic interactions present in our work. In the case when rods are allowed to cross each other, the tendency to jam is reduced.^{66} However, the mechanism for forming giant clusters is now different. Reference 65 adds an interesting study of hydrodynamically interacting flagella of pusher type. Thus, their collective dynamics is less comparable to our system.

In future investigations, we plan to extend our study to the collective dynamics of pusher- and puller-type squirmer rods as introduced in Ref. 17. In addition, it is straightforward to extend the squirmer-rod model to active flexible filaments by adding bending rigidity. This will provide valuable insights how hydrodynamic interactions influence the collective dynamics of active filaments and bridge between dry models^{33,95} and models with implicit hydrodynamic interactions,^{96,97} on the one hand, and continuum descriptions,^{98,99} on the other hand. Clearly, the extended MPCD method will prove beneficial for simulations of dense suspensions of active filaments by providing a less compressible fluid environment in the regime of vanishing Reynolds numbers.

## SUPPLEMENTARY MATERIAL

See the supplementary material for additional video files V1–V6, which show the dynamic evolution of the systems shown in Figs. 5(a), 5(b), and 6(b)–6(e), in the order of their occurrence.

## ACKNOWLEDGMENTS

We thank Josua Grawitter for helpful discussions on the topic of the manuscript. We also acknowledge financial support from the Collaborative Research Center 910 funded by the Deutsche Forschungsgemeinschaft.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

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