Squirmer rods as elongated microswimmers: flow fields and confinement

Abstract

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 different 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 different slab widths. Based on the hydrodynamic multipole expansion either for bulk or confinement between two parallel plates, we categorize the different 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.

---

1 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–3 Intriguing features such as active turbulence,^4,5 swarming,^6,7 polar patterns,⁸ 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.¹⁰ But also the elongated body shapes^4,7,8 and geometric confinement^11–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–21 Elongated and flexible active constituents are realized by polar biofilaments in motility assays⁸ or at a fluid–fluid interface,⁷ while bacteria provide a natural realization of active rods.⁴

Combining elements of the Toner–Tu²² and Swift–Hohenberg²³ theories, continuum models for bacterial microswimmers have been developed that are able to reproduce pattern formation such as vortices and active turbulence.^4,24–26 Alternative continuum theories for active matter based on liquid-crystal hydrodynamics reproduce the creation, annihilation, and motion of liquid-crystal defects²⁷ 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–31 Langevin dynamics simulations of active rods^4,32,33 or active filaments^34,35 suggest that many properties might already emerge from short-ranged 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–42 while with the method of multi-particle collision dynamics collective phases of spherical^43–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–54 and also biological organisms such as Volvox.⁵⁵ 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–45,60–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 large-scale 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).

Fig. 1 Perspective view of a squirmer rod moving in the midplane (in light gray) of a Hele-Shaw cell with slab width Δz. Periodic boundary conditions along the x and y directions are used.

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–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.

2 System and methods

The hydrodynamics of the squirmer rod at small Reynolds numbers is determined by the Stokes equations including the incompressibility condition,

η∇² u − ∇p + f  =  0 ∇·u  =  0. (1)

Here, η 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 multi-particle collision dynamics (MPCD) to simulate the flow fields together with the relevant parameters.

2.1 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 α = l_S/2R, where l_S is the rod length. We do not go beyond a maximum squirmer distance of d ≈ 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

v_s = B^s₁[(ê·x̂_s)x̂_s − ê] (2)

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₁ = B^s₁R³ as explained in the following section. The squirmer parameter B^s₁ also determines the swimming speed v₀, which amounts to v₀ = 2/3B^s₁ for a single squirmer.

Fig. 2 (a) Schematic of a single squirmer of radius R and with orientation given by the unit vector ê. The surface slip-velocity field of a neutral squirmer is indicated by blue arrows. (b) Schematic of the squirmer rod model. Always, n_sq = 10 spherical squirmers are placed on a straight line with distance d to form active rods. All squirmer orientations ê are aligned with the rod axis. (c) Implementation of a pusher type squirmer rod. The surface slip velocity on the rod surface is multiplied by the envelop function f(x_s*·ê) of eqn (3) such that the slip-velocity field is concentrated on the rear of the rod.

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.

2.2 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 force-free microswimmers, the flow field of the source dipole, and their higher-order derivatives. We introduce the polar angle θ relative to the squirmer-rod axis via cos θ = ê·r/r. Then, in spherical coordinates r, ϕ, and θ the respective radial and polar components of the velocity field can be written as^50,66Here, P_n(cos θ) are the ordinary Legendre polynomials andwhere ′ means derivative with respect to cos θ. Due to the rotational symmetry of a squirmer rods about its axis, the flow field is independent of the azimuthal angle ϕ and the azimuthal velocity component u_ϕ vanishes. The multipole coefficients A_n and B_n describe the strength of the nth-order force and source multipoles, where B₁ 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),using the completeness relation of the Legendre polynomials, which gives

u_r,n(r) = A_nr⁻ⁿ + B_nr⁻ⁿ⁻². (7)

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.

2.3 Fundamental solutions of Stokes flow in Hele-Shaw geometry

In a seminal paper Liron and Mochon¹² 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⁶⁵ one can generate the flow fields of higher-order force moments by applying the directional derivative ê·∇_p to the Stokeslet flow field.¹⁴ Here, ∇_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 ∇_p² 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.¹⁴ 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):†

u(r) = u_SD(r) + u_FD(r) + u_FQ(r) + v_SO(r) + …. (8)

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 ũ(ρ,φ) in the x–y plane using polar coordinates ρ, φ. 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 givesHere, T_n(cos φ) = cos(nφ) are the Chebyshev polynomials of the first kind and [scr A, script letter A]_n, [scr B, script letter B]_n are the respective coefficients of the force and source multipole moments. They start with [scr A, script letter A]₁ for the force monopole, which does not exist for a force-free swimmer, and [scr B, script letter B]₁ 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 Δz. In Appendix A we motivate [scr A, script letter A]_n ∝ Δz and [scr B, script letter B]_n ∝ 1/Δz. 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 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 nφ = T_n(cos φ), we extract the multipole moments from eqn (9) by projecting ũ_ρ(ρ,φ) on the Chebyshev polynomials,

which yieldsIn 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 ([scr B, script letter B]₁), force dipole ([scr A, script letter A]₂), force quadrupole ([scr A, script letter A]₃), and source octupole ([scr B, script letter B]₃).

2.4 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).⁴⁴ 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₀, positions x_i(t), and velocities v_i(t) move ballistically during time Δt,

x_i(t + Δt) = x_i(t) + v_i(t)Δt. (12)

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₀ and centered around ξ. First, for the n_ξ particles in each cell with volume [scr V, script letter V]_ξ the mean velocity v_ξ and center-of-mass position x_ξ are determined. Then, in the center-of-mass frame the fluid particles are assigned new random velocities δv_i from a Boltzmann-distribution with temperature T₀. To restore overall momentum conservation, the total change in linear momentum, , has to be subtracted from the new velocities, while an additional term containing the difference of the angular momentum before the collision , and the change in angular momentum, , is added to preserve angular momentum. Here x_i,c denotes the particle position relative to the center-of-mass position x_ξ. Thus, the collision step can be summarized by

v^new_i = v_ξ + δv_i − Δv_ξ − x_i,c × I_ξ⁻¹(L_ξ − ΔL_ξ), (13)

Here, I_ξ is the moment-of-inertia tensor of the distribution of particles inside [scr V, script letter V]_ξ 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.⁶⁷ 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.

2.5 Geometry and fluid parameters

For the MPCD fluid we use a density of n₀ = 10/a₀³ fluid particles per cell and the same mass density ρ₀ = m₀n₀ for the immersed squirmer rod. With time step one obtains a dynamic fluid viscosity of .⁴⁴ Throughout this work, the radius and the squirner parameter of the single squirmer are chosen as R= 3a₀ and , respectively.

To determine how the hydrodynamic moments of the squirmer flow field in a bulk fluid varies with the aspect ratio α in Section 3.1, we simulate single rods in a cubic box of linear size L = 100a₀ 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₀ with periodic boundary conditions. We begin by simulating the system for a time 10⁴Δt to equilibrate the MPCD fluid flow and then average the flow fields over the time interval from 10⁴Δt to 5 × 10⁶Δt. To determine the radial component v_r(r,θ) 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₀ with no-slip walls and periodic boundary conditions along the x and y directions. For the slab width we investigate the three values Δz/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⁷Δt to compensate for the fact that the flow field cannot be averaged about the rod axis.

3 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 α. 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.

3.1 Neutral squirmer rods in bulk fluid

In the bulk fluid, rods of different aspect ratio α 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, α = 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 α = 1.75 the flow field already becomes stretched along the rod [cf.Fig. 3(b)]. At higher aspect ratio, α = 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 α.

Fig. 3 Simulated flow fields of an active squirmer rod for different aspect ratios α = 1 (a), α = 1.75 (b), and α = 4.0 (c) shown in the laboratory frame. The color-coded magnitude of the velocity field and streamlines are presented. In each case only a part of the simulated fluid volume is shown.

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 α = 4 and also include fitted polynomials in r⁻¹ as solid lines.

Fig. 4 (a) Radial Legendre coefficients u_r,n(r) of the simulated flow field of a neutral squirmer rod with aspect ratio α = 4.0 calculated from eqn (6) (symbols) and with fitted power laws (solid lines). SD and FQ stand for source dipole and force quadrupole, respectively. (b) Hydrodynamic multipole moments in units of the squirmer parameter B^s₁R³ plotted versus aspect ratio α. (c) Swimming speed v_rod of a squirmer rod plotted versus α. v_rod is normalized by the swimming speed v₀ = 2/3B^s₁ of a single squirmer at α = 1. Simulation data (blue dots) and analytical values as given by eqn (15) (dashed black line).

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³ decay, indicating a zero force monopole moment A₁ while the source dipole moment B₁ 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₂ = 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³, which indicates an additional force quadrupole moment A₃. 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₁ ≠ 0), and a force quadrupole moment (A₃ ≠ 0), which both decay as 1/r³. In Fig. 4(b) we plot the two moments normalized by the squirmer parameter B^s₁versus the aspect ratio α. They both increase roughly linearly in α. While the source-dipole moment B₁ shows a modest increase in α starting from B₁ = B^s₁R³ for α = 1, the force quadrupole moment A₃ increases much more strongly from zero and clearly dominates beyond α = 2.

The additional force quadrupole does not break the head–tail 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₃(cos θ) compared to P₁(cos θ) 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³.

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, . 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₁) + F_hyd(y_Nsq). In the following, we refer the force of the other squirmers on this remaining force, δF_hyd(y_i) = F_hyd(y_i) − F_hyd,p/(N_sq − 2), so that they sum up to zero: . In Fig. 5(a) we plot the relative force δF_hyd(y_i), the schematic in Fig. 5(b) rationalizes the different signs of δF_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.

Fig. 5 (a) Local relative force δF_hyd(y_i), with which the squirmer rod acts on the fluid, plotted along the long axis y of the squirmer rod. The aspect ratio is α = 4. The dots correspond to the axial positions y_i of the constituent squirmers and the line is a guide to the eye. Note, in contrast to the rest of the article, N_sq = 20 squirmers are used to compose the rod in order to generate a smoother surface. (b) Sketch of δF_hyd(y) (blue arrows) to illustrate the force quadrupole. The grey arrow indicates the swimming velocity.

Finally, as we show in Fig. 4(c), the swimming speed v_rod of the rod increases with α 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 α 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.⁷⁰ 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₁ = 3v₀/2 on the straight region is assumed. For α = 1 we reproduce the swimming velocity v₀ 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.⁷⁰ This is used in the derivation of eqn (14) and obviously no longer valid for the passive squirmer rod.

3.2 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 (lower panel). For a squirmer rod of α = 4.0 and similar to the bulk fluid [cf.Fig. 2(c)], we observe how the streamlines at the sides of the rod buckle inwards (upper panel of 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 [scr B, script letter B]₁ is of longer range than the one of the force quadrupole with moment [scr A, script letter A]₃. 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.

Fig. 6 Simulated flow field of a squirmer rod with aspect ratio α = 4.0 in a Hele-Shaw geometry with slab width Δz = 6R = 18a₀ shown in the laboratory frame. The color-coded magnitude of the velocity field v and streamlines are presented. Upper panel: Flow field in the midplane; lower panel: flow field in the cross-sectional plane indicated by the red dashed line in the upper panel.

In the lower panel of Fig. 6, further away from the rod (|x/l_S| > 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| < 3/4, near-field 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 φ) using eqn (10). Since from the radial decay of the expansion coefficients ũ_ρ,n(ρ) 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 ũ_ρ,1(ρ) and ũ_ρ,3(ρ), plotted in Fig. 7 and 8, to a source dipole (SD) and force quadrupole (FQ), respectively.

Fig. 7 Multipole analysis of the simulated flow field of a neutral squirmer rod with aspect ratio α = 4 for different slab widths of the Hele-Shaw cell: (a) Δz = 2.7R, (b) Δz = 6.0R and (c) Δz = 9.3R. As in the bulk fluid the radial decays of the expansion coefficients ũ_ρ,n(ρ) 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 ρ⁻² and ρ⁻⁴, respectively. Grey vertical bar: the coefficients ũ_ρ,n(ρ) were determine for radial distances ρ ≥ l_S/2. The gray vertical dashed line indicates ρ = Δz.

Fig. 8 Multipole analysis of the simulated flow field of a neutral squirmer rod confined in a Hele-Shaw cell with slab width Δz = 2.7R for different aspect ratios: (a) α = 1.75, (b) α = 3.25, and (c) α = 4.0. Otherwise, the same description as in Fig. 7 is used.

In Fig. 7 we plot ũ_ρ,1/3(ρ) averaged over the cell height for different slab widths Δz. For a small width Δz = 2.7R = 0.34l_S the coefficients are in very good agreement with the expected power-law decay: 1/ρ² for the source dipole and 1/ρ⁴ for the force quadrupole. Both power laws fit well over the whole range of the radial distance. For slab width Δz = 6R = 0.75l_S the measured coefficients ũ_ρ,n(ρ) approach the predicted power-law decays at a radial distance of approximately ρ/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 Δz = 9.3R = 1.2l_S the measured coefficients ũ_ρ,n(ρ) show clear deviations from the expected power-law decay over the whole range of the radial distance ρ. All in all, when we vary the slab width, the force quadrupole moment stays the same, while the source dipole moment increases with decreasing Δz.

In Fig. 8 we choose the smallest slab width Δz/R = 2.7, where we obtained the best agreement in the previous figure and plot the coefficients ũ_ρ,1/3(ρ) for different aspect ratios α. 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).

3.3 Pusher-type squirmer rods

As explained in Section 2.1 we also implemented a pusher-type squirmer rod with aspect ratio α = 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² starting from r ≈ 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³ decay (as indicated by the solid line) to a decay of roughly 1/r². 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 ≈ 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⁵ decay of a source octupole with positive moment B₃ at distances r/l_S < 1, while at r/l_S >1 we indicate the decay of a force quadrupole with negative moment A₃. Again, the latter is probably strongly disturbed by the periodic images, resulting in a decay of roughly 1/r², 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.

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 α = 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.

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, ũ_ρ,2(ρ) ∝ 1/ρ³, as described by theory. At distances ρ/l_S < 1.3 the strong but shorter ranged force-dipole flow field dominates, while for ρ/l_S > 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 Δz. Thus, the dipole-force moment increases with Δz in qualitative agreement with [scr A, script letter A]_n ∝ Δz (cf. Section 2.3), while the source-dipole moment decreases with Δz again in qualitative agreement with [scr B, script letter B]_n ∝ 1/Δz.

Fig. 10 Multipole analysis of the simulated flow field of a pusher-type squirmer rod with aspect ratio α = 4 for different slab widths of the Hele-Shaw cell: (a) Δz = 2.7R, (b) Δz = 6.0R and (c) Δz = 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.

The third-order coefficient ũ_ρ,3(ρ) ∝ 1/ρ⁴ does not change very significantly in one direction with increasing Δz. We take this as an indication that both the force quadrupole and source octupole contribute to the flow field. A more detailed analysis of [scr A, script letter A]₃ and [scr B, script letter B]₃ 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.

4 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 becomes dominant beyond α = 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 Δz we find a good match with our simulations. The flow field of the source dipole now decay as 1/ρ² and dominates at radial distances ρ > Δz for all α and Δz over the field of the force quadrupole, which now decays as 1/ρ⁴.

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². However, in the Hele-Shaw geometry the radial decay changes to 1/ρ³ 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 strain-rate tensor⁷¹ 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,¹³ 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 [capital G, script] for solving the Stokes equations in a Hele-Shaw cell of width Δz (cf.Fig. 1) using an infinite series of image point forces.¹² 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₃(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 two-dimensional hydrodynamics source dipole, which appears in the second line,Here, the indices i, j ∈ {1, 2} belong to the x–y plane, ρ 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₃ = 0 and Δz, and an equivalent term with the z coordinate of the swimmer, r^sw₃, appears.

In the Stokes flow regime any flow field surrounding a solid body can be written as the sum of hydrodynamic multipoles or 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 ê·∇_pn times on [capital G, script]_ijê_j, where ∇_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 Δz from the Stokeslet to the force-multipole moments, so that we have [scr A, script letter A]_n ∝ Δz. Again following ref. 14 and 65, we derive the flow field of a source dipole by the operation, u_i,SD(r) ∝ ∇_p²[capital G, script]_ijê_j. While the second derivative with respect to r^sw_s removes the term with r^sw_s and generates a factor 1/Δz², the remaining two-dimensional Laplace operator acting on the source dipole term gives zero. Thus, one arrives atwhich is nearly identical to the Stokeslet. Flow fields of higher source multipoles are generated by the directional derivate ê·∇_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/Δz to the source-multipole moments, so that we have [scr B, script letter B]_n ∝ 1/Δz. 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 (ρ,φ,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 atThe resulting radial velocity component can also be decomposed into cell-height averaged multipole moments ũ_ρ,n. For the relevant multipole moments we have

ũ_ρ(ρ,φ) =  ũ_ρ,FD(ρ,φ) + ũ_ρ,FQ(ρ,φ) + ũ_ρ,SD(ρ,φ) +ũ_ρ,SO(ρ,φ) + …. (23)

which we calculate using eqn (18)–(21). For the force dipole and quadrupole, we obtainandrespectively. The source dipole and octupole giveandrespectively, where the T_n(cos φ) = cos(nφ) 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 [scr A, script letter A]₂, [scr A, script letter A]₃,₁, and [scr B, script letter B]₃ by a numerical fit of the power law ρ⁻ⁿ to ũ_ρ,n(ρ).

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 (ρ,φ,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 ρ₀, 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 ρ₀, , we first calculate the I_zz components of the constituent parts. Integrating ρ₀ρ² over the enclosed volume yields

For the component I_xx the integrand becomes , which yieldsWe 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′ 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′ − mz_com². The distance z_com between z = 0 and the center of mass follows by integrating zρ₀/m over the enclosed volume,For the middle segments the mass distribution is symmetric about z = 0 and hence I_xx = I_xx′.

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. (32)

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.⁷⁰ If (v,σ) are the velocity and stress tensor that solve the Stokes equation for the force-free swimmer and (v′,σ′) the corresponding fields for a passive object of the same shape pulled now by a force F′, thenHere 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′ = −6πηRU′, where U′ is the translational velocity, and the surface force density is constant and proportional to U′. 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′ 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 and θ_max = π − θ_min. We obtain S_rod  =  4πR²α. Similarly, we integrate the scalar product of the surface velocity v_s and the orientation vector ê

v_s·ê = B^s₁(cos² θ − 1) (37)

over the surface of one complete sphere and n_s − 1 middle segments,and obtainThe swimming velocity of the squirmer rod is ultimately given by eqn (15).

Acknowledgements

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

References

1. E. Lauga and T. R. Powers, Phys. Rev. Lett., 2009, 72, 096601.
2. M. C. Marchetti, J.-F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao and R. A. Simha, Rev. Mod. Phys., 2013, 85, 1143.
3. J. Elgeti, R. G. Winkler and G. Gompper, Rep. Prog. Phys., 2015, 78, 056601.
4. H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher, R. E. Goldstein, H. Löwen and J. M. Yeomans, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 14308.
5. D. Nishiguchi, I. S. Aranson, A. Snezhko and A. Sokolov, Nat. Commun., 2018, 9, 1–8.
6. S. Thutupalli, R. Seemann and S. Herminghaus, New J. Phys., 2011, 13, 073021.
7. T. Sánchez, D. T. N. Chen, S. J. DeCamp, M. Heymann and Z. Dogic, Nature, 2012, 491, 431.
8. V. Schaller, C. Weber, C. Semmrich, E. Frey and A. R. Bausch, Nature, 2010, 467, 73–77.
9. H. Wioland, F. G. Woodhouse, J. Dunkel, J. O. Kessler and R. E. Goldstein, Phys. Rev. Lett., 2013, 110, 268102.
10. A. Zöttl and H. Stark, J. Phys.: Condens. Matter, 2016, 28, 253001.
11. J. Blake and A. Chwang, J. Eng. Mech., 1974, 8, 23–29.
12. N. Liron and S. Mochon, J. Eng. Mech., 1976, 10, 287–303.
13. T. Brotto, J.-B. Caussin, E. Lauga and D. Bartolo, Phys. Rev. Lett., 2013, 110, 038101.
14. A. J. Mathijssen, A. Doostmohammadi, J. M. Yeomans and T. N. Shendruk, J. Fluid Mech., 2016, 806, 35–70.
15. S. Thutupalli, D. Geyer, R. Singh, R. Adhikari and H. A. Stone, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 5403–5408.
16. D. L. Koch and G. Subramanian, Annu. Rev. Fluid Mech., 2011, 43, 637–659.
17. M. Bär, R. Großmann, S. Heidenreich and F. Peruani, Annu. Rev. Condens. Matter Phys., 2019, 11, 441–466.
18. I. Theurkauff, C. Cottin-Bizonne, J. Palacci, C. Ybert and L. Bocquet, Phys. Rev. Lett., 2012, 108, 268303.
19. J. Palacci, S. Sacanna, A. P. Steinberg, D. J. Pine and P. M. Chaikin, Science, 2013, 339, 936–940.
20. D. Nishiguchi and M. Sano, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 92, 052309.
21. C. C. Maass, C. Krüger, S. Herminghaus and C. Bahr, Annu. Rev. Condens. Matter Phys., 2016, 7, 171–193.
22. J. Toner and Y. Tu, Phys. Rev. Lett., 1995, 75, 4326.
23. J. Swift and P. C. Hohenberg, Phys. Rev. A: At., Mol., Opt. Phys., 1977, 15, 319.
24. J. Dunkel, S. Heidenreich, K. Drescher, H. H. Wensink, M. Bär and R. E. Goldstein, Phys. Rev. Lett., 2013, 110, 228102.
25. J. Dunkel, S. Heidenreich, M. Bär and R. E. Goldstein, New J. Phys., 2013, 15, 045016.
26. S. Heidenreich, J. Dunkel, S. H. Klapp and M. Bär, Phys. Rev. E, 2016, 94, 020601.
27. A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans and F. Sagués, Nat. Commun., 2018, 9, 1–13.
28. F. Jülicher, S. W. Grill and G. Salbreux, Rep. Prog. Phys., 2018, 81, 076601.
29. T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen and O. Shochet, Phys. Rev. Lett., 1995, 75, 1226.
30. G. Grégoire and H. Chaté, Phys. Rev. Lett., 2004, 92, 025702.
31. H. Chaté, F. Ginelli, G. Grégoire, F. Peruani and F. Raynaud, Eur. Phys. J. B, 2008, 64, 451–456.
32. H. Wensink and H. Löwen, J. Phys.: Condens. Matter, 2012, 24, 464130.
33. M. C. Bott, F. Winterhalter, M. Marechal, A. Sharma, J. M. Brader and R. Wittmann, Phys. Rev. E, 2018, 98, 012601.
34. Ö. Duman, R. E. Isele-Holder, J. Elgeti and G. Gompper, Soft Matter, 2018, 14, 4483–4494.
35. K. Prathyusha, S. Henkes and R. Sknepnek, Phys. Rev. E, 2018, 97, 022606.
36. M. Hennes, K. Wolff and H. Stark, Phys. Rev. Lett., 2014, 112, 238104.
37. H. Jeckel, E. Jelli, R. Hartmann, P. K. Singh, R. Mok, J. F. Totz, L. Vidakovic, B. Eckhardt, J. Dunkel and K. Drescher, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 1489–1494.
38. A. Laskar, R. Singh, S. Ghose, G. Jayaraman, P. S. Kumar and R. Adhikari, Sci. Rep., 2013, 3, 1–5.
39. A. Laskar and R. Adhikari, Soft Matter, 2015, 11, 9073–9085.
40. J. de Graaf, H. Menke, A. J. Mathijssen, M. Fabritius, C. Holm and T. N. Shendruk, J. Chem. Phys., 2016, 144, 134106.
41. A. Pandey, P. S. Kumar and R. Adhikari, Soft Matter, 2016, 12, 9068–9076.
42. M. Kuron, P. Stärk, C. Burkard, J. De Graaf and C. Holm, J. Chem. Phys., 2019, 150, 144110.
43. A. Zöttl and H. Stark, Phys. Rev. Lett., 2014, 112, 118101.
44. J. Blaschke, M. Maurer, K. Menon, A. Zöttl and H. Stark, Soft Matter, 2016, 12, 9821–9831.
45. A. Zöttl and H. Stark, Eur. Phys. J. E: Soft Matter Biol. Phys., 2018, 41, 61.
46. M. Theers, E. Westphal, G. Gompper and R. G. Winkler, Soft Matter, 2016, 12, 7372–7385.
47. M. Theers, E. Westphal, K. Qi, R. G. Winkler and G. Gompper, Soft Matter, 2018, 14, 8590–8603.
48. F. J. Schwarzendahl and M. G. Mazza, Soft Matter, 2018, 14, 4666–4678.
49. F. J. Schwarzendahl and M. G. Mazza, J. Chem. Phys., 2019, 150, 184902.
50. M. Lighthill, Commun. Pure Appl. Math., 1952, 5, 109–118.
51. J. R. Blake, J. Fluid Mech., 1971, 46, 199–208.
52. F. Sciortino, A. Giacometti and G. Pastore, Phys. Rev. Lett., 2009, 103, 237801.
53. S. Jiang, Q. Chen, M. Tripathy, E. Luijten, K. S. Schweizer and S. Granick, Adv. Mater., 2010, 22, 1060–1071.
54. I. Buttinoni, J. Bialké, F. Kümmel, H. Löwen, C. Bechinger and T. Speck, Phys. Rev. Lett., 2013, 110, 238301.
55. T. J. Pedley, D. R. Brumley and R. E. Goldstein, J. Fluid Mech., 2016, 798, 165–186.
56. T. Ishikawa, M. Simmonds and T. Pedley, J. Fluid Mech., 2006, 568, 119–160.
57. M. T. Downton and H. Stark, J. Phys.: Condens. Matter, 2009, 21, 204101.
58. I. O. Götze and G. Gompper, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 82, 041921.
59. G.-J. Li and A. M. Ardekani, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 013010.
60. J.-T. Kuhr, J. Blaschke, F. Rühle and H. Stark, Soft Matter, 2017, 13, 7458–7555.
61. F. Rühle, J. Blaschke, J.-T. Kuhr and H. Stark, New J. Phys., 2018, 20, 025003.
62. J.-T. Kuhr, F. Rühle and H. Stark, Soft Matter, 2019, 15, 5685–5694.
63. A. Malevanets and R. Kapral, J. Chem. Phys., 1999, 110, 8605–8613.
64. G. Gompper, T. Ihle, D. Kroll and R. Winkler, Advanced computer simulation approaches for soft matter sciences III, Springer, Berlin, 2009, pp. 1–87.
65. S. E. Spagnolie and E. Lauga, J. Fluid Mech., 2012, 700, 105–147.
66. J. Blake, Math. Proc. Cambridge Philos. Soc., 1971, 70, 303–310.
67. M. P. Allen and D. J. Tildesley, Computer simulation of liquids, Oxford University Press, Oxford, 2017.
68. J. de Graaf, A. J. Mathijssen, M. Fabritius, H. Menke, C. Holm and T. N. Shendruk, Soft Matter, 2016, 12, 4704–4708.
69. M. P. Howard, A. Z. Panagiotopoulos and A. Nikoubashman, Comput. Phys. Commun., 2018, 230, 10–20.
70. H. A. Stone and A. D. Samuel, Phys. Rev. Lett., 1996, 77, 4102.
71. G. B. Jeffery, Proc. R. Soc. London, Ser. A, 1922, 102, 161–179.

---

Footnote

† 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 ũ_ρ,2(ρ) grows with the slab width Δz in our simulations, which identifies it as a force dipole rather than a source quadrupole. We will elaborate on this further below.

---