Squirmer rods as elongated microswimmers: flow fields and confinement

Microswimmers or active elements, such as bacteria and active filaments, have an elongated shape, which determines their individual and collective dynamics. There is still a need to identify what role long-range hydrodynamic interactions play in their fascinating dynamic structure formation. We construct rods of diﬀerent aspect ratios using several spherical squirmer model swimmers. With the help of the mesoscale simulation method of multi-particle collision dynamics we analyze the flow fields of these squirmer rods both in a bulk fluid and in Hele-Shaw geometries of diﬀerent slab widths. Based on the hydrodynamic multipole expansion either for bulk or confinement between two parallel plates, we categorize the diﬀerent multipole contributions of neutral as well as pusher-type squirmer rods. We demonstrate how confinement alters the radial decay of the flow fields for a given force or source multipole moment compared to the bulk fluid.


Introduction
Swimming at the micron scale and understanding how microorganisms overcome the constraints of low-Reynolds-number hydrodynamics have attracted a lot of attention among physicists. [1][2][3] Intriguing features such as active turbulence, 4,5 swarming, 6,7 polar patterns, 8 and vortex formation, 8,9 which arise from the active nature of the fluid, still pose challenges to scientists. The fluid environment and the flow fields initiated by microswimmers strongly determine both their isolated motion and their appealing collective patterns, which they develop in non-equilibrium. 10 But also the elongated body shapes 4,7,8 and geometric confinement [11][12][13][14][15] of microswimmers significantly contribute to their single and collective dynamics. However, it is still a matter of current debate to completely identify the individual contributions of direct steric and long-range hydrodynamic interactions. 2,16,17 By introducing an active squirmer rod in this article, we aim to contribute to the discussion with a detailed analysis of its flow field in bulk and Hele-Shaw geometry.
Artificial microscopic swimmers with various locomotion mechanisms have recently been constructed to study the principal properties of their individual and collective motion. Examples mostly include spherical microswimmers. 6,[18][19][20][21] Elongated and flexible active constituents are realized by polar biofilaments in motility assays 8 or at a fluid-fluid interface, 7 while bacteria provide a natural realization of active rods. 4 Combining elements of the Toner-Tu 22 and Swift-Hohenberg 23 theories, continuum models for bacterial microswimmers have been developed that are able to reproduce pattern formation such as vortices and active turbulence. 4,[24][25][26] Alternative continuum theories for active matter based on liquid-crystal hydrodynamics reproduce the creation, annihilation, and motion of liquid-crystal defects 27 or are able to describe the dynamics of the intracellular cytoskeleton gel. 2,28 In contrast, particle-based simulations are employed to study the collective dynamics of active rods. In the famous Vicsek model coarse-grained aligning interactions are able to reproduce dynamical states such as flocking and swarming. [29][30][31] Langevin dynamics simulations of active rods 4,32,33 or active filaments 34,35 suggest that many properties might already emerge from shortranged steric interactions. However, also implicit hydrodynamic pair interactions are included in such simulations to provide more realistic models with novel dynamic states. 36,37 An extension to single active filaments exists. 38,39 Explicit hydrodynamic simulation schemes have been applied, as well. For example, the lattice Boltzmann method was used to investigate the properties of microswimmers with various shapes [40][41][42] while with the method of multi-particle collision dynamics collective phases of spherical [43][44][45] and ellipsoidal microswimmers 46,47 were studied, as well as realizations of pusher and puller-type swimmers. 48,49 In this work, we introduce and characterize single squirmer rods with the idea to model elongated microswimmers in their fluid environment. They then can be used to disentangle the effect of direct steric and long-range hydrodynamic interactions in their collective behavior. The squirmer model swimmer 45,50,51 is a good model for spherical artificial microswimmers such as Janus particles [52][53][54] and also biological organisms such as Volvox. 55 Squirmers have in common that they propel themselves by an axisymmetric surface-velocity field, which acts on the surrounding fluid. 45,56,57 In biological systems this is realized by cilia, located all over the cell surface, which perform synchronized collective non-reciprocal motions. The squirmer has frequently been used in hydrodynamic simulations of the collective behavior of microswimmers 42,58,59 and also by our own group, [43][44][45][60][61][62] where we rely on the mesoscale method of multi-particle collision dynamics. 45,63,64 In this article we use squirmers to build rigid rod-shaped microswimmers of different aspect ratios and perform largescale simulations with multi-particle collision dynamics. We explore their hydrodynamic flow fields both in a bulk fluid and in a Hele-Shaw geometry with varying cell width, where we keep the squirmer rod in the midplane to mimic the fluid interface in the experiments of ref. 7 (cf. Fig. 1).
In our analysis, we rely on the hydrodynamic multipole expansion both in the bulk fluid 10,40,50,65 and between two parallel plates. [12][13][14] This will enable us to categorize the different multipoles of both neutral as well as pusher-type squirmer rods and to determine the different radial decays of their multipole flow fields.
The paper is organized as follows. In Section 2 we introduce the squirmer rod and the methods to generate and analyze its flow fields. We present the results of our analysis in Section 3 and conclude in Section 4.

System and methods
The hydrodynamics of the squirmer rod at small Reynolds numbers is determined by the Stokes equations including the incompressibility condition, Here, Z is the viscosity of the fluid, u(r) the fluid velocity, p(r) the pressure, and f(r) denotes a body force at point r. In the following we introduce the squirmer-rod model, the fundamental solutions of the Stokes equations both in bulk and a Hele-Shaw cell, in order to characterize its hydrodynamic multipole moments, and details of the method of multiparticle collision dynamics (MPCD) to simulate the flow fields together with the relevant parameters.

Squirmer rod model
To investigate the influence of shape anisotropy of microswimmers and their hydrodynamic interactions on collective motion, we introduce a new model of squirming active rods. To construct them, we arrange several spherical squirmers [cf. Fig. 2(a)] on a line to form a rigid body as shown in Fig. 2(b). On the one hand, this closely resembles the actual rod shape of many biological microswimmers such as bacteria, on the other hand the model can be extended to a flexible swimmer body by introducing bending rigidity in future works. In concrete, we always take N sq = 10 spherical squirmers and vary the squirmer distance d [cf. Fig. 2(b)] to form rods of different aspect ratios a = l S /2R, where l S is the rod length. We do not go beyond a maximum squirmer distance of d E 0.8R, which amounts to a maximum rod length l S = 9.5R, so that the surface is still smooth enough that swimmers can slide past each other. Each of the individual squirmers of a rod propels itself by an imposed axisymmetric surface slip-velocity field of a neutral squirmer, 45,56,57 that acts on the surrounding fluid. Here, x s is the unit vector along x s , which points from the squirmer center to its surface, and unit vector ê indicates the orientation and swimming direction of the squirmer. In the case of a single squirmer, the induced flow field in a bulk fluid agrees with the velocity field of a pure source-dipole singularity with strength B 1 = B s 1 R 3 as explained in the following section. The squirmer parameter B s 1 also determines the swimming speed v 0 , which amounts to v 0 = 2/3B s 1 for a single squirmer. Having the surface velocity field distributed over the whole squirmer rod [cf. Fig. 2(b)], our model resembles ciliated microorganisms such as Paramecium. In contrast, bacteria like  E. coli propel themselves with a bundle of rotating flagella that pushes fluid backwards while the cell body does not have any surface velocity field. To implement squirmer rods of such a pusher or also of puller type, we use the envelope function Here, x s * is a vector from the center of mass of the rod to any point on the surface [cf. Fig. 2(c)] and 10/l S is the step width of the envelope function. Multiplying the surface velocity field with f (x s *Áê), concentrates the slip velocity field on the rear half of the rod for the -sign in f and thereby a pusher squirmer rod is realized, while the + sign generates a puller.

Fundamental solutions of Stokes flow in bulk
To provide an understanding of the flow field u(r) that squirmer rods initiate in a bulk fluid and how they interact hydrodynamically, we review the multipole expansion of u(r). 50,66 It is composed of singular solutions of the Stokes equations including the Stokeslet as the flow field of the leading force monopole, which has to vanish for forcefree microswimmers, the flow field of the source dipole, and their higher-order derivatives. We introduce the polar angle y relative to the squirmer-rod axis via cos y = êÁr/r. Then, in spherical coordinates r, f, and y the respective radial and polar components of the velocity field can be written as 50,66 u r ðr; yÞ ¼ X 1 n¼1 A n r Àn þ B n r ÀnÀ2 Â Ã P n ðcos yÞ; Here, P n (cos y) are the ordinary Legendre polynomials and V n ðyÞ ¼ 2 sin y nðn þ 1Þ P n 0 ðcos yÞ; where 0 means derivative with respect to cos y. Due to the rotational symmetry of a squirmer rods about its axis, the flow field is independent of the azimuthal angle f and the azimuthal velocity component u f vanishes. The multipole coefficients A n and B n describe the strength of the nth-order force and source multipoles, where B 1 belongs to the source dipole. To determine them, we calculate the nth expansion or Legendre coefficient of the radial component of the flow field in eqn (4), u r;n ðrÞ 2n þ 1 2 ð p 0 u r ðr; yÞP n ðcos yÞ sin ydy; using the completeness relation of the Legendre polynomials, which gives u r,n (r) = A n r Àn + B n r ÀnÀ2 .
Thus, from the simulated flow fields of the squirmer rods, in particular, from the radial decay of u r,n (r), we can infer their leading force and source multipole moments. We will demonstrate this in Section 3.1. Note that the nth force and source monopole can be clearly distinguished by their radial decay.

Fundamental solutions of Stokes flow in Hele-Shaw geometry
In a seminal paper Liron and Mochon 12 determined the Stokeslet solution between two infinitely extended parallel plates with no-slip boundary condition (cf. Hele-Shaw geometry in Fig. 1) using an infinite series of image forces. For a force parallel to the plates (x À y plane), the x, y components of the resulting Stokeslet flow field are long-range and consist of a Poiseuille flow profile g(z) along the plate normal times the flow field of a two-dimensional source dipole in the plane (cf. Appendix A). As in the bulk fluid 65 one can generate the flow fields of higher-order force moments by applying the directional derivative êÁr p to the Stokeslet flow field. 14 Here, r p is the nabla operator acting on the position of the Stokeslet singularity and the unit vector ê gives the direction of the force multipole in the x À y plane, which we assume to be uniaxial for simplicity. Likewise, by applying the Laplace operator r p 2 to the Stokeslet flow field, one obtains the flow field of a source dipole, which here is up to a factor identical to the Stokeslet flow field. 14 Higher source multipole moments are again generated by directional derivatives. Relevant for the squirmer rod will be the flow fields of the source dipole (u SD ), force dipole (u FD ), force quadrupole (u SD ), and source octupole (u SD ): † To determine these contributions from the simulated flow field of the squirmer rod, we first average over the z coordinate and then treat the resulting velocity field ũ(r,j) in the x-y plane using polar coordinates r, j. As we demonstrate in Appendix A, for a microswimmer oriented along the x axis, ê = ê x , the multipole expansion for the radial component of the flow field givesũ r ðr; jÞ ¼ X 1

n¼1
A n þ B n r nþ1 T n ðcos jÞ: Here, T n (cos j) = cos(nj) are the Chebyshev polynomials of the first kind and A n , B n are the respective coefficients of the force and source multipole moments. They start with A 1 for the force monopole, which does not exist for a force-free swimmer, and B 1 for the source dipole, and so on. It immediately becomes obvious that in contrast to the bulk fluid one cannot distinguish the flow fields of force and source multipoles with the same angular dependence (same order n) by the radial decay. We note that the coefficients scale differently with the slab width Dz. In Appendix A we motivate A n p Dz and B n p 1/Dz. However, our simulated data are not always sufficiently accurate to discriminate both cases. So, we assign the relevant multipole assuming that it is preserved from the bulk fluid. The multipole expansion captures the distribution † Although the source quadrupole shares the same radial dependence as the force dipole in the Hele-Shaw geometry, it is not considered in this list for two reasons. First, we solely observe a force dipole in the bulk fluid but no source quadrupole. Second, the strength of ũ r,2 (r) grows with the slab width Dz in our simulations, which identifies it as a force dipole rather than a source quadrupole. We will elaborate on this further below. of force and source multipoles in leading order. This distribution generated by the rod surface should remain the same in the Hele-Shaw geometry. However, since we deal with a voluminous rod in contrast to a point-like source, higher-order terms enter to fulfill the boundary conditions on both the rod surface and the bounding plates. The flow fields of these terms will decay faster than the leading bulk moments at some distance from the rod.
Finally, exploiting the orthogonality relations of cos nj = T n (cos j), we extract the multipole moments from eqn (9) by projecting ũ r (r,j) on the Chebyshev polynomials, 0ũ r ðr; jÞT n ðcos jÞdj; which yieldsũ r;n ðrÞ ¼ In both the bulk fluid and Hele-Shaw geometry we can restrict our focus on a particular set of Stokes flow singularities. These will be the source dipole (B 1 ), force dipole (A 2 ), force quadrupole (A 3 ), and source octupole (B 3 ).

MPCD fluid model
To model the fluid in our simulations, we apply the method of multi-particle collision dynamics together with the Andersen thermostat and angular momentum conservation (MPCD-AT+a). 44 The method solves the Navier-Stokes equations and thereby can also treat hydrodynamic interactions between the active squirmer rod and the confining walls of the Hele-Shaw geometry. 63,64 Furthermore, it also includes thermal fluctuations. The MPCD method considers point-like particles that represent the fluid and defines rules for (i) their motion and (ii) their collisions such that the resulting hydrodynamic flow field fulfills the Navier-Stokes equations. These rules are applied in alternating steps. We shortly introduce the essential features of the MPCD method and refer to ref. 45 for more details. During the streaming step (i) the point particles with masses m 0 , positions x i (t), and velocities v i (t) move ballistically during time Dt, During this step fluid particles collide with confining walls or the surface of the squirmer rod. By applying the so-called bounce-back rule 44,45 the collisions either enforce the no-slip boundary condition on confining walls or the slip-velocity field on the squirmer-rod surface. In addition, the collisions also transfer both linear and angular momentum, in particular, to the squirmer rod. This changes the center-of-mass and angular momentum of the rigid rod, which is calculated relative to the center-of-mass. In Appendix B we present formulas for the mass m rod and moment-of-inertia tensor I rod of the rigid squirmer rod.
For the collision step (ii), the simulation volume is divided by a cubic lattice and the fluid particles are grouped into the cubic unit cells of linear size a 0 and centered around x. First, for the n x particles in each cell with volume V x the mean velocity v x and center-of-mass position x x are determined. Then, in the center-of-mass frame the fluid particles are assigned new random velocities dv i from a Boltzmann-distribution with temperature T 0 . To restore overall momentum conservation, the total change in linear momentum, be subtracted from the new velocities, while an additional term containing the difference of the angular momentum before Here x i,c denotes the particle position relative to the center-of-mass position x x . Thus, the collision step can be summarized by Here, I x is the moment-of-inertia tensor of the distribution of particles inside V x calculated in the center-of-mass frame. During this step immersed boundaries are represented by so-called ''ghost'' particles. They are added to the collision cells divided by the boundaries and interact with the other fluid particles. These ghost particles are assigned the local velocity of the translating and rotating squirmer rod plus a random thermal velocity drawn from a Boltzmann distribution. The changes of linear and angular momentum of the ghost particles during the collision step are then added to the squirmer rod to ensure linear and angular momentum conservation. Finally, before performing each collision step, the lattice is randomly shifted to ensure Galilean invariance. Given the center-of-mass and angular momentum of the rigid squirmer rod, its location and orientation are updated 10 times during each collision step using a standard leapfrog algorithm. 67 The rod is treated as single rigid body the position and orientation vector of which is updated following ref. 67. To keep the squirmer rod in the midplane of the Hele-Shaw cell, we only use the x and y component of the fluid force acting on the rod to integrate its motion in time throughout the simulation. Due to the symmetry about the midplane only Brownian forces are acting on the squirmer normal to the midplane. Note due to this constraint we do not observe any oscillatory trajectories between the plates as observed in ref. 68.
We have implemented the MPCD model fluid together with the squirmer-rod model in C++ and CUDA to enable the use of graphic cards. Because the main performance bottleneck in MPCD is memory access, we integrate a sorting algorithm and lookup table following ref. 69. Additionally, we use variable size cooperative thread groups to add dynamic load balancing to our MPCD collision routines in CUDA.

Geometry and fluid parameters
To determine how the hydrodynamic moments of the squirmer flow field in a bulk fluid varies with the aspect ratio a in Section 3.1, we simulate single rods in a cubic box of linear size L = 100a 0 using periodic boundary conditions in all three dimensions. For all other bulk simulations in Sections 3.1 and 3.3, we use box sizes of L = 180a 0 with periodic boundary conditions. We begin by simulating the system for a time 10 4 Dt to equilibrate the MPCD fluid flow and then average the flow fields over the time interval from 10 4 Dt to 5 Â 10 6 Dt. To determine the radial component v r (r,y) of the velocity field, we also exploit the rotational symmetry of the flow fields by averaging about the rod axis. Using eqn (6), we then obtain the expansion coefficients u r,n (r) for different multipole order n and extract the strengths of the hydrodynamic moments by fitting u r,n (r) from eqn (7) to the curves determined from the simulated flow fields.
For the simulations in Hele-Shaw geometry in Sections 3.2 and 3.3, we use box sizes of L = 200a 0 with no-slip walls and periodic boundary conditions along the x and y directions. For the slab width we investigate the three values Dz/R = 2.7, 6.0, and 9.3. When determining the flow fields in the Hele-Shaw geometry, the total simulation time is increased up to 10 7 Dt to compensate for the fact that the flow field cannot be averaged about the rod axis.

Results
In the following, we present our results. We first discuss the flow field of a neutral squirmer rod in the bulk fluid and how it depends on the aspect ratio a. We then study neutral squirmer rods in Hele-Shaw geometries of different slab widths and illustrate the different far-field behavior compared to the bulk fluid. Finally, we look at a squirmer rod of pusher type moving in both bulk fluid and the Hele-Shaw geometry.

Neutral squirmer rods in bulk fluid
In the bulk fluid, rods of different aspect ratio a create flow fields as shown in Fig. 3. In all panels (a), (b), and (c) only a part of the simulated fluid volume is shown. In the case of the spherical squirmer, a = 1, in Fig. 3(a) we observe a perfect match with the theoretical prediction for the flow field of a source dipole moment. For a short squirmer rod with a = 1.75 the flow field already becomes stretched along the rod [cf. Fig. 3(b)]. At higher aspect ratio, a = 4.0, the flow field shows a strong deviation from the one of a source dipole moment, with the streamlines buckling inwards at the sides of the rod. Thus, we expect higher-order moments to contribute strongly to squirmer rods with higher aspect ratios a.
To provide a quantitative analysis, we apply eqn (6) and decompose the flow field into its different angular contributions with the nth-order Legendre coefficients u r,n (r). As an example, we show in Fig. 4(a) the measured u r,n (r) as data points for a squirmer rod of a = 4 and also include fitted polynomials in r À1 as solid lines.
All radial velocity components u r,n (r) either show distinct power law behavior as for n = 1 and 3 or vanish. The first-order coefficient u r,1 (r) clearly shows a pure 1/r 3 decay, indicating a zero force monopole moment A 1 while the source dipole moment B 1 is present, as expected. This is consistent with the fact that we consider microswimmers free of external forces. Consistent with the head-tail symmetry of the forces, the squirmer rods exert on the fluid, we do not observe a force dipole moment, A 2 = 0, since the second-order coefficient u r,2 (r) vanishes. Furthermore, the third-order contribution u r,3 (r) is also present and decays with 1/r 3 , which indicates an additional force quadrupole moment A 3 . We elaborate on its origin further below. Both curves u r,1 (r) and u r,3 (r) fall off from the theoretically predicted power law at large r due to the finite size of the simulation box. They also show small deviations from the power law, which we attribute to the flow fields of the image swimmers introduced by the use of periodic boundary conditions. All other coefficients u r,n (r) are too small to be distinguished from noise and hence can be neglected in the following. Thus, the flow field of a neutral squirmer rod is a superposition of a source dipole moment (B 1 a 0), and a force quadrupole moment (A 3 a 0), which both decay as 1/r 3 . In Fig. 4(b) we plot the two moments normalized by the squirmer parameter B s 1 versus the aspect ratio a. They both increase roughly linearly in a. While the source-dipole moment B 1 shows a modest increase in a starting from B 1 = B s 1 R 3 for a = 1, the force quadrupole moment A 3 increases much more strongly from zero and clearly dominates beyond a = 2.
The additional force quadrupole does not break the headtail symmetry of the squirmer rod, which we therefore denote as neutral. However, the corresponding flow field clearly affects the shape of the overall flow field as we saw in Fig. 3. Since its angular dependence is governed by the third-order Legendre polynomial P 3 (cos y) compared to P 1 (cos y) of the source dipole, it will influence the hydrodynamic interactions with other squirmer rods, also because its flow field shows the same radial decay, 1/r 3 .
To shed some light on the origin of the force quadrupole, we determined the local force F hyd (y i ), with which the spherical squirmer component i placed at y i along the rod axis acts on the surrounding fluid. Of course, the rod is force free, thus all the forces F hyd (y i ) add up to zero, P Nsq i¼1 F hyd ðy i Þ ¼ 0. The forces from the terminal squirmers i = 1 and i = N sq are mainly determined by fluid pressure. They cancel each other up to some remaining force, F hyd,p = F hyd ( y 1 ) + F hyd (y N sq ). In the following, we refer the force of the other squirmers on this remaining force, dF hyd (y i ) = F hyd (y i ) À F hyd,p /(N sq À 2), so that they sum up to zero: P NsqÀ1 i¼2 dF hyd ðy i Þ ¼ 0. In Fig. 5(a) we plot the relative force dF hyd (y i ), the schematic in Fig. 5(b) rationalizes the different signs of dF hyd (y i ). In the center of the rod (y = 0), the relative local force is negative due to the surface velocity of the squirmer rod pushing fluid backwards, while at both ends of the rod the fluid is dragged with the rod. This generates the observed force quadrupole.
Finally, as we show in Fig. 4(c), the swimming speed v rod of the rod increases with a and then saturates. With increasing aspect ratio the swimming speed is more and more determined by the maximum surface velocity around the equator of the constituent squirmers and thus the speed increases. At large a the cap regions of the squirmer rod become irrelevant and the swimming speed saturates. To calculate an analytic expression for the swimming speed, we average the surface slip-velocity field v s of the squirmer rod over the whole surface with area S rod . 70 For the component along the squirmer axis ê we obtain The derivation is sketched in Appendix C. Note that the first two terms on the right-hand side of eqn (15) result when a perfect spherocylinder with the maximum slip velocity B s 1 = 3v 0 /2 on the straight region is assumed. For a = 1 we reproduce the swimming velocity v 0 of the spherical squirmer. Fig. 4(c) shows the analytic expression overestimates the swimming speed. The reason is that the surface average of eqn (14) to calculate the swimming speed only applies to spherical shapes since then the hydrodynamic surface stress for a passive sphere is constant along the surface. 70 This is used in the derivation of eqn (14) and obviously no longer valid for the passive squirmer rod.

Neutral squirmer rods in Hele-Shaw geometry
To discuss the flow fields of a squirmer rod in the confining Hele-Shaw geometry, we start with Fig. 6 where we show the color-coded strength of the flow field and streamlines in the midplane of the cell (upper panel) and the cross-sectional plane   Fig. 6) due to the flow field of the force quadrupole moment. However, as explained in Section 2.3 and detailed in Appendix A, in the Hele-Shaw geometry the flow field of the source dipole with moment B 1 is of longer range than the one of the force quadrupole with moment A 3 . Indeed, already at a distance l S from the rod the buckling of the flow lines is much weaker compared to the bulk fluid [cf. Fig. 3(c)] and the flow field assumes the shape of a pure source dipole.
In the lower panel of Fig. 6, further away from the rod (|x/l S | 4 1) the streamlines are approximately parallel to the bounding plates in agreement with the Poiseuille flow profile of the single force and source multipoles as detailed in Appendix A. However, in the direct vicinity of the rod, |x/l S | o 3/4, nearfield flow along the z direction occurs, which is due to terms with exponential decay in the full hydrodynamic solution in a slab geometry.
We now present a more quantitative analysis of the simulated flow fields in Fig. 7 and 8. As for the bulk fluid, we decompose the flow fields into the different angular contributions given by the Chebyshev polynomials T n (cos j) using eqn (10). Since from the radial decay of the expansion coefficients ũ r,n (r) we cannot distinguish between the force multipole of nth order and the source multipole of n + 1th order, we are guided by the analysis in the bulk fluid from Section 3.1 and attribute the expansion coefficients ũ r,1 (r) and ũ r,3 (r), plotted in Fig. 7 and 8, to a source dipole (SD) and force quadrupole (FQ), respectively.
In Fig. 7 we plot ũ r,1/3 (r) averaged over the cell height for different slab widths Dz. For a small width Dz = 2.7R = 0.34l S the coefficients are in very good agreement with the expected power-law decay: 1/r 2 for the source dipole and 1/r 4 for the force quadrupole. Both power laws fit well over the whole range of the radial distance. For slab width Dz = 6R = 0.75l S the measured coefficients ũ r,n (r) approach the predicted power-law decays at a radial distance of approximately r/l S = 1. This is larger than the slab width, where we expect the far-field solutions of the different multipoles to become valid. At this distance the corresponding streamlines shown in the lower panel of Fig. 6 are parallel to the bounding plates as predicted by the multipole flow fields. Finally, increasing the slab width further to Dz = 9.3R = 1.2l S the measured coefficients ũ r,n (r) show clear deviations from the expected power-law decay over the whole range of the radial distance r. All in all, when we vary the slab width, the force quadrupole moment stays the same, while the source dipole moment increases with decreasing Dz.
In Fig. 8 we choose the smallest slab width Dz/R = 2.7, where we obtained the best agreement in the previous figure and plot the coefficients ũ r,1/3 (r) for different aspect ratios a. Again there is very good agreement between the coefficients and the expected power law decay. Only at radial distance smaller than l S one realizes deviations, as expected. This is, in particular, visible in plot (a).

Pusher-type squirmer rods
As explained in Section 2.1 we also implemented a pusher-type squirmer rod with aspect ratio a = 4.0 and show our analysis in Fig. 9. In the bulk fluid, the force dipole now dominates the flow field as the new coefficient u r,2 (r) demonstrates, which decays with 1/r 2 starting from r E l S . Compared to the neutral squirmer rod and the dominating force dipole, the flow field of the source dipole is roughly an order of magnitude smaller. We also observe a deviation from the expected 1/r 3 decay (as indicated by the solid line) to a decay of roughly 1/r 2 . One explanation might be the hydrodynamic interaction of the squirmer rod with its images due to the periodic boundary conditions, which especially affects the weaker multipole moments. Compared to the neutral rod, the force-dipole field is more long-ranged and therefore hydrodynamic interactions between the images are stronger, which especially influences the field of higher-order multipoles. Very different from the neutral squirmer rod, the third-order coefficient u r,3 (r) changes sign at a radial distance of r E l S from positive to negative with increasing r. This is why we plot the magnitude of the coefficients u r,n (r). It is tempting to fit the 1/r 5 decay of a source octupole with positive moment B 3 at distances r/l S o 1, while at r/l S 41 we indicate the decay of a force quadrupole with negative moment A 3 . Again, the latter is probably strongly disturbed by the periodic images, resulting in a decay of roughly 1/r 2 , similar to our observations before. In conclusion, since the flow fields of both the source dipole and force quadrupole are one order of magnitude weaker compared to the force dipole, the overall flow field is well described by the latter.
In Fig. 10 we analyze the flow field of the pusher-type squirmer rod in the Hele-Shaw geometry for different slab widths. In addition to the findings for the neutral squirmer rod (cf. Fig. 7), we also observe the force dipole. Compared to the bulk fluid, it now has a stronger radial decay, ũ r,2 (r) p 1/r 3 , as described by theory. At distances r/l S o 1.3 the strong but shorter ranged force-dipole flow field dominates, while for r/l S 4 1.3 the slower decay of the source-dipole field takes over. Increasing the slab width, the flow field of the force dipole always dominates in the given radial range [cf. Fig. 10(b) and (c)]. The reason is the green curve shifts upwards while the blue curve shifts downwards with increasing Dz. Thus, the dipole-force moment increases with Dz in qualitative agreement with A n p Dz (cf. Section 2.3), while the source-dipole moment decreases with Dz again in qualitative agreement with B n p 1/Dz. The third-order coefficient ũ r,3 (r) p 1/r 4 does not change very significantly in one direction with increasing Dz. We take this as an indication that both the force quadrupole and source octupole contribute to the flow field. A more detailed analysis of A 3 and B 3 is not feasible with the hydrodynamic MPCD method since due to thermal fluctuations it requires long averaging in order to obtain smooth flow lines at large distances. Therefore, larger system sizes would require an immense computational effort.

Summary and conclusion
In this paper we introduced the squirmer rod as a new model for elongated microswimmers. By varying aspect ratio and surface slip-velocity field, it is able to describe artificial and biological microswimmers of different shape and propulsion type, such as pushers, pullers, and neutral swimmers. To quantify the generated hydrodynamic flow fields of the squirmer rod, we determined the hydrodynamic multipole moments both in the bulk fluid and the Hele-Shaw geometry, by projecting the simulated flow fields on Legendre or Chebyshev polynomials, respectively. The corresponding expansion coefficients showed the expected radial decay.
The flow field of the neutral squirmer rod in the bulk fluid shows the expected source dipole, while a force quadrupole moment develops linearly with increasing aspect ratio and Fig. 7 Multipole analysis of the simulated flow field of a neutral squirmer rod with aspect ratio a = 4 for different slab widths of the Hele-Shaw cell: (a) Dz = 2.7R, (b) Dz = 6.0R and (c) Dz = 9.3R. As in the bulk fluid the radial decays of the expansion coefficients ũ r,n (r) feature a source dipole (n = 1, blue dots) and force quadrupole (n = 3, red pluses). The solid lines show fits with the corresponding power laws of r À2 and r À4 , respectively. Grey vertical bar: the coefficients ũ r,n (r) were determine for radial distances r Z l S /2. The gray vertical dashed line indicates r = Dz.  Fig. 7 is used. Fig. 9 Absolute values of the radial Legendre coefficients |u r,n (r)| of the simulated flow field of a pusher-type squirmer rod with aspect ratio a = 4.0 calculated from eqn (6) (symbols) and with fitted power laws (solid lines). FD, FQ stand for force dipole and quadrupole, while SD, SO mean source dipole and octupole, respectively.
becomes dominant beyond a = 2. It is due to a non-uniform distribution of the force, with which the rod acts on the fluid. By taking an average of the surface slip-velocity field over the rod surface, the actual swimming velocity is overestimated. In the Hele-Shaw geometry the radial decay of the multipole flow fields changes as predicted by theory. Especially at low slab width Dz we find a good match with our simulations. The flow field of the source dipole now decay as 1/r 2 and dominates at radial distances r 4 Dz for all a and Dz over the field of the force quadrupole, which now decays as 1/r 4 .
For the pusher-type squirmer rod with noticeable elongation we observe that the flow field is composed of four hydrodynamic moments: force dipole, source dipole, force quadrupole, and source octupole. In bulk the force dipole completely dominates the flow field and determines the radial decay with 1/r 2 . However, in the Hele-Shaw geometry the radial decay changes to 1/r 3 and is less long-ranged. Nevertheless, for larger slab widths and the recorded radial distances it dominates the flow field, while for small slab widths we see a cross over to the longer-ranged source-dipole field. This is in qualitative agreement with the expected behavior of the strength of the multipole moments with increasing slab width: the force dipole becomes stronger while the source dipole weakens. Finally, in the Hele-Shaw geometry the flow fields of force quadrupole and source octupole have the same radial decay. Varying slab width suggests that they both contribute.
Our work shows how elongated microswimmers generate additional hydrodynamic multipole moments compared to swimmers of spherical shape, which leads to a more complex appearance of the generated flow field. Since rods experience additional torques in a non-uniform flow field via the strainrate tensor 71 this will generate different behavior in suspensions of squirmer rods compared to their pure steric interactions. Furthermore, the source dipole is predicted to be the dominant hydrodynamic moment in a Hele-Shaw geometry, 13 which we confirm for neutral squirmer rods for varying slab width. However, for pusher and puller rods the dominance of the force dipole depends on the slab width and radial distance.
Therefore, in the continuation of this work, we plan to investigate the collective dynamic behavior of the squirmer rods introduced in this work. Thereby we will gain an understanding how hydrodynamic interactions through the self-generated flow fields contribute to different types of dynamic behavior such as swarming in active nematics or active turbulent phases. As a further extension of the presented model, we plan to introduce bending rigidity between the different components of the squirmer rod to model active flexible filaments used, for example, in ref. 7 and 8.

Conflicts of interest
There are no conflicts to declare.

A Force and source multipole moments in Hele-Shaw geometry
In their seminal paper Liron and Mochon derived Green's function G for solving the Stokes equations in a Hele-Shaw cell of width Dz (cf. Fig. 1) using an infinite series of image point forces. 12 For a point force along the z axis normal to the bounding plates, the Stokeslet decays exponentially and the same applies to the flow field component u 3 (r) along the z direction. In the following we concentrate on the case, where the point force is directed along unit vector ê in the x-y plane and only flow in this plane is monitored. Then the main feature of the Stokeslet is the long-range flow field of a twodimensional hydrodynamics source dipole, which appears in the second line, Here, the indices i, j A {1, 2} belong to the x-y plane, r is the corresponding radial distance from the point force, and r the radial unit vector. In the first line on the right hand side a Poiseuille flow profile is visible, which vanishes at the plate locations r 3 = 0 and Dz, and an equivalent term with the z coordinate of the swimmer, r sw 3 , appears. In the Stokes flow regime any flow field surrounding a solid body can be written as the sum of hydrodynamic multipoles or Fig. 10 Multipole analysis of the simulated flow field of a pusher-type squirmer rod with aspect ratio a = 4 for different slab widths of the Hele-Shaw cell: (a) Dz = 2.7R, (b) Dz = 6.0R and (c) Dz = 9.3R. In addition to the source dipole (n = 1, blue dots) and force quadrupole/source octupole (n = 3, red pluses), the force dipole (n = 2, green crosses) is observed. Otherwise, the same description as in Fig. 7 is used.
singularities, which can be derived from the Stokeslet in eqn (16). For the flow field of the squirmer rod we list the relevant terms: u(r) = u FD (r) + u FQ (r) + u SD (r) + u SO (r). . .. (17) Since the squirmer rod is force free, there only appears the flow field of a force dipole [u FD (r)] and, in addition, the flow fields of a force quadrupole [u FQ (r)], a source dipole [u SD (r)], and a source octupole [u SO (r)]. 14,65 Following ref. 14 and 65 we derive the flow fields of force multipoles with uniaxial symmetry from the Stokeslet in eqn (16) by applying the directional derivative êÁr p n times on G ij ê j , where r p acts on the location of the point force. For n = 1 and 2 we thus obtain the respective flow fields of a hydrodynamic force dipole, and force quadrupole, Note that we allocated a factor Dz from the Stokeslet to the forcemultipole moments, so that we have A n p Dz. Again following ref. 14 and 65, we derive the flow field of a source dipole by the operation, u i,SD (r) p r p 2 G ij ê j . While the second derivative with respect to r sw s removes the term with r sw s and generates a factor 1/Dz 2 , the remaining two-dimensional Laplace operator acting on the source dipole term gives zero. Thus, one arrives at v i;SD ðrÞ ¼ À6 which is nearly identical to the Stokeslet. Flow fields of higher source multipoles are generated by the directional derivate êÁr p acting on v SD (r). Of relevance is the second derivative, which gives the flow field of a source octupole: Again, we have allocated a factor 1/Dz to the source-multipole moments, so that we have B n p 1/Dz. In the following, we restrict ourselves to a squirmer rod oriented along the x-axis and swimming in the midplane of the Hele-Shaw cell, i.e., ê = ê x and To convert the flow fields in the Hele-Shaw geometry to a form similar to the one used in free Stokes flow, we switch to cylindrical coordinates (r,j,z). We eliminate the z-dependence of the multipole flow fields by integrating over the Poiseuille profile along z, project the resulting velocity field on the radial direction, and arrive atũ r ðr; jÞ ê r Dz Á ð Dz 0 uðr; j; zÞdz: The resulting radial velocity component can also be decomposed into cell-height averaged multipole moments ũ r,n . For the relevant multipole moments we have ũ r (r,j) = ũ r,FD (r,j) + ũ r,FQ (r,j) + ũ r,SD (r,j) +ũ r,SO (r,j) + . . ..
which we calculate using eqn (18)- (21). For the force dipole and quadrupole, we obtaiñ u r;FD ðr; jÞ ¼ ÀA 2 T 2 ðcos jÞ r 3 respectively, where the T n (cos j) = cos(nj) are Chebyshev polynomials of the first kind. These formulas agree with the general result given in eqn (9) in the main text. Using eqn (10) and (11), we can then determine the force and source moments A 2 , A 3 , 1 , and B 3 by a numerical fit of the power law r Àn to ũ r,n (r).

B Moment-of-inertia tensor of the squirmer rod
The moment-of-inertia tensor I rod of the squirmer rod can be calculated applying Steiner's theorem to the constituent parts of the squirmer rod. First, we calculate the masses m mid and m cap as well as the moment-of-inertia tensors I mid and I cap of the middle segments and end caps relative to their respective centers of mass. Due to the rotational symmetry around the rod axis, I rod , I mid and I cap have two degenerate eigenvalues in the xÀy plane perpendicular to the rod axis and a different eigenvalue along the rod axis ê = ẑ.
In each of the constituent parts of the squirmer rods, we use cylindrical coordinates (r,j,z), where z = 0 is the center of each spherical squirmer. The constituent parts are sections of spheres of radius R that extend from z = h min to z = h max . Integration over the enclosed volume with constant mass density r 0 , gives the mass For the middle segments of the rod, Àh min = h max = d/2, and for the end caps h min = Àd/2 and h max = R. Based on the known formula for the moment-of-inertia tensor for a solid body with uniform mass density r 0 , I ij ¼ r 0 Ð ðx 2 À x i x j ÞdV, we first calculate the I zz components of the constituent parts. Integrating r 0 r 2 over the enclosed volume yields For the component I xx the integrand becomes r 0 ðz 2 þ r 2 sin 2 jÞ, which yields We attach a dash here since for the end caps the mass distribution is not symmetric about z = 0 and therefore the center of mass of the caps does not coincide with the reference point x = 0, for which I xx 0 was calculated. For further use below, we need the moment of inertia I xx for a rotation axis perpendicular to the z axis through the center of mass and calculate it using Steiner's theorem: I xx = I xx 0 À mz com 2 . The distance z com between z = 0 and the center of mass follows by integrating zr 0 /m over the enclosed volume, For the middle segments the mass distribution is symmetric about z = 0 and hence I xx = I xx 0 .
We now use the previous results to calculate the relevant moments of inertia I rod,ij for the squirmer rod. For the moment I rod,zz the relevant rotational axis along ẑ goes through all the centers of mass of the n s segments. Thus, we can just add up their moments of inertia: I rod,zz = 2I cap,zz + (n s À 2)I mid,zz .
For the remaining non-zero moments, I rod,xx = I rod,yy , we have to shift the respective moments of the middle segments and end caps, I mid,xx and I cap,xx , using their distances from the rod's center of mass and Steiner's theorem: for an even number of squirmers per rod n s .

C Swimming velocity
In principle, following ref. 70 the swimming velocity of a squirmer of arbitrary shape can be calculated analytically using the reciprocal theorem. 70 If (v,s) are the velocity and stress tensor that solve the Stokes equation for the force-free swimmer and (v 0 ,s 0 ) the corresponding fields for a passive object of the same shape pulled now by a force F 0 , then Here U is the swimming velocity, n the surface normal, v s the slip velocity at the swimmer surface and S the surface of the swimmer's body. For a sphere of radius R the acting force is F 0 = À6pZRU 0 , where U 0 is the translational velocity, and the surface force densityn Á s 0 ¼ À 3Z 2R U 0 is constant and proportional to U 0 . Thus, for a sphere the swimming velocity U is simply the negative average of the surface slip velocity v s To provide an estimate for the swimming velocity of the squirmer rod, we apply the approximation that the surface force density is constant along the surface of the passive rod, when pulled by force F 0 along its long axis. The estimate for the swimming velocity is thus the average of the surface slip velocity v s as given in eqn (14).
To determine the swimming velocity, we calculate the surface of the squirmer rod, which adds up the surface of one complete sphere and n s À 1 middle segments of a sphere. The integration limits for the polar angles of the middle segments are y min ¼ arccos a À 1 n s À 1 and y max = p À y min . We obtain S rod = 4pR 2 a. Similarly, we integrate the scalar product of the surface velocity v s and the orientation vector ê v s Áê = B s 1 (cos 2 y À 1) over the surface of one complete sphere and n s À 1 middle segments, ð and obtain ð S rod v s ÁêdS ¼ 4pR 2 B s 1 a À 1 3 À n s À 1 3 a À 1 n À 1 3 " # : The swimming velocity of the squirmer rod is ultimately given by eqn (15).