Activity-induced clustering in model dumbbell swimmers: The role of hydrodynamic interactions

Using a fluid-particle dynamics approach, we numerically study the effects of hydrodynamic interactions on the collective dynamics of active suspensions within a simple model for bacterial motility: each microorganism is modeled as a stroke-averaged dumbbell swimmer with prescribed dipolar force pairs. Using both simulations and qualitative arguments, we show that, when the separation between swimmers is comparable to their size, the swimmers' motions are strongly affected by activity-induced hydrodynamic forces. To further understand these effects, we investigate semidilute suspensions of swimmers in the presence of thermal fluctuations. A direct comparison between simulations with and without hydrodynamic interactions shows these to enhance the dynamic clustering at a relatively small volume fraction; with our chosen model the key ingredient for this clustering behavior is hydrodynamic trapping of one swimmer by another, induced by the active forces. Furthermore, the density dependence of the motility (of both the translational and rotational motions) exhibits distinctly different behaviors with and without hydrodynamic interactions; we argue that this is linked to the clustering tendency. Our study illustrates the fact that hydrodynamic interactions not only affect kinetic pathways in active suspensions, but also cause major changes in their steady state properties.


I. INTRODUCTION
It is well established that various microorganisms, such as bacteria and algae, propel themselves through a suspending medium using a nonreciprocal cyclic motion [1,2]. Despite the progress that has been made in elucidating the self-propulsion mechanisms of various microorganisms in isolation, an understanding of their collective dynamics is still elusive. Recent experimental and simulation studies have shown that complicated interactions arising from the microorganisms' activity produce diverse nonequilibrium cooperative phenomena, whose behaviors are notably different from those observed in passive systems. Such fascinating aspects of active systems are suggestive of new underlying principles; the search for these is the current focus of intense study among the soft matter physics community (see recent reviews [3][4][5] and references therein).
Hydrodynamic interactions (HIs), and their dynamic coupling to activity, are thought to be among the key elements that govern transport and rheological properties in suspensions of swimming microorganisms (see [6] for a review). The hydrodynamic interactions between a pair of microorganisms have been extensively studied and their general features have been well revealed in certain simple situations; for example, when two microorganisms are far apart compared to their sizes, a weak attractive interaction acts between them. This is because the dominant contribution is of a dipole-dipole character (which is attractive for particles that are free to rotate).
However, we are still far from a thorough understanding of the role of HIs in nondilute systems of active par-ticles. When the distances between microorganisms are comparable to their sizes, HIs in the near field become important, but the effects of such interactions are not easy to understand in general. This is because the details of both the swimming mechanisms and the particle shapes come into play with increased proximity (where not just dipolar terms but higher order multipoles contribute), so that the HIs among microorganisms at moderate or high density are more complicated and less universal than the dilute case. Thus, it is interesting and important to examine many-body HIs through simulations of minimal microorganism models. From such models it may be possible to extract some generic features of cooperative phenomena in active suspensions of various types.
It is well recognized that, in a wide variety of systems, activity produces marked dynamic interactions, resulting in collective fluctuations or structure formation. Examples include giant number fluctuations [7][8][9], clustering (or swarming) [10,11], and bulk phase separation [12][13][14][15][16][17][18][19]. These phenomena can occur without any direct potential interactions, and therefore have a purely dynamical origin. In systems manifesting such a dynamic cooperativity, a fundamental question is how HIs influence or change the activity-induced interactions.
Although certain collective effects characteristic of active systems are believed to be greatly influenced by HIs, there is still no general consensus on their roles in such phenomena. Indeed, for the reasons given above, this can depend on particle shape and/or swimming mechanism so that such questions must be asked within the context of well defined models. One popular model is that of squirmers: spherical particles with a prescribed surface velocity field [20,21]. For squirmers in two di-mensions it was found that the phase separation, otherwise caused by a collisional reduction in mean propulsion speed of swimmers at high density [19], is switched off by HIs [22]. Another popular model is based on dumbbell shaped swimmers subject to a set of discrete forces acting on the two solid particles and on the fluid [23][24][25]. As we explain further below, this model may need to be supplemented by an ad hoc rule to govern the case where a second swimmer is about to occupy the spatial position at which the active force from a first swimmer acts on the fluid. Even with such a rule, the model seems a better starting point than squirmers when modeling suspensions of rod-like motile bacteria. Specifically, squirmers can only exert active torques on each other, whereas rodlike swimmers can exert both active and passive torques. Note that torques are particularly important in active systems (e.g. [22]) since they rotate not just the particle, but also its swimming direction.
In this study, to further understand the role of HIs in collective dynamics for nonspherical swimmers, we numerically investigate the dynamics of non-dilute active suspensions of self-propelled dumbbells in three dimensions. We especially focus on activity-induced clustering at modest concentrations [10,11]. This can be seen as a precursor to activity-induced phase separation. However we do not address full phase separation, which can only be simulated accurately with enormous system sizes [26], such that a study of the role of near-field HIs in active phase separation remains out of reach with current computers. (Some far-field effects might instead be represented within a continuum model [27].) Our simulation method fully takes into account the solvent dynamics and thus is suitable for the investigation of HIs, as will be explained in the next section. In Sec. III, we demonstrate that the activity-induced HIs can strongly alter both the translational and rotational dynamics at modest densities. At a relatively small volume fraction, an activity-induced hydrodynamic attraction dominates the collective dynamics of the swimmers, resulting in a stronger slowing down, and hence enhanced clustering, compared to the case without HIs. The volume fraction dependence of the motility is found to be much stronger with HIs than without them; this also influences the clustering tendency, as we discuss in Sec. IV. Section V gives our conclusions and a further discussion.

II. NUMERICAL METHODS
Many-body hydrodynamic interactions in semidilute regimes are difficult to deal with even using simulations because of their intrinsically long-range and timedependent character; in principle we must solve the Navier-Stokes equation with moving boundary conditions at the solvent-particle interfaces. To confront these difficulties, in the past decade, several hybrid simulation techniques [28][29][30][31][32] as well as other mesoscopic methods [33,34] for the dynamics of complex colloidal suspen-sions have been developed. Of these methods, we here adopt the fluid-particle-dynamics (FPD) method [30,31] to incorporate HIs into the study of a model active suspension.
Within the framework of the FPD method, HIs can be approximately taken into account (i) by treating a rigid colloidal particle as a nondeformable but highly viscous fluid particle, and (ii) by replacing the particlefluid boundary with a smoothed boundary. This simple scheme considerably reduces the numerical cost and the mathematical complexity of the computation even while preserving the basic nature of HIs in many-body colloidal systems. The validity of this method has been examined by several authors [30,31,[35][36][37][38]. It is good at capturing both far-field and intermediate-field aspects of the HIs, though like many other methods it cannot (at reasonable computational cost) also resolve the divergence of lubrication forces that occurs when hard suspended particles come into very close contact [30,31]. In practice the model can be considered to describe colloids whose hard-core radius (set by the interparticle pair potential) slightly exceeds their hydrodynamic radius, so that thin lubrication films never dominate.
A. Model swimmer system Here, a swimming bacterium is represented by two spherical particles with radius aH = aT = a and a "phantom" spherical particle with radius aP = b. The center-to-center distance between the head and tail particles and that between the phantom and tail particles are fixed at ℓ0 and ℓ1, respectively. The back-to-back (a) and face-to-face (b) force configurations correspond respectively to the pusher (extensile) and puller (contractile) swimming mechanisms. The (green) solid and (purple) dashed arrows respectively indicate schematically the swimming direction and the flow field around the swimmer.
Our model of a swimming microorganism comprises a dumbbell with a prescribed dipolar force pair, which is essentially the same as the models used in Refs. [23,24] and later in Ref. [39]. Each dumbbell is composed of two real particles plus one "phantom" particle -so called be-cause (a) it merely follows the motion of the pair of real particles to which it is attached and (b) it can overlap with the other real particles in nearby swimmers. The phantom particle can be thought of as modeling the effect of a thin flagellar bundle, whereby a force is exerted on the fluid at a position displaced from the rod-like bacterial body (represented by the dimer of real particles).
The positions R m and the radii a m of the particles and the distances of separation between them determine the shape of the swimmer, where the subscripts m = H and T denote the head and tail particles of the dumbbell, respectively, and m = P denotes the phantom particle (see Fig. 1). The swimmer's orientation is characterized by a unit vectorn = ( The propulsion force −f actn is exerted on the fluid via the phantom particle; an equal and opposite force f actn is exerted on the tail particle, which ensures the dipolar character of the swimming mechanism. Hence there is no net external force on a box of fluid that fully encloses the swimmer. Face-to-face (f act > 0) and back-to-back (f act < 0) force configurations correspond to "pusher" (extensile) and "puller" (contractile) microorganisms respectively ( Fig. 1). Although, in this paper, we limit our study to the dynamics of pushers such as motile bacteria, it would be straightforward to perform a corresponding simulation for pullers of the same geometry.
A similar dumbbell model, and its hydrodynamic instabilities, were theoretically investigated in Ref. [40] by analyzing hydrodynamic equations derived by a coarsegraining approach. It is one of several models, used to address active suspensions, in which surface stresses or force configurations are prescribed [41][42][43][44][45]. An alternative approach is to specify a surface slip velocity, which is usually done for spherical particles resulting in the squirmer model mentioned previously [20,21,46]. Recently, several groups have investigated many-body HIs in squirmer suspensions [22,[47][48][49][50][51][52]. Although both force-prescribed (rodlike) and velocity-prescribed (spherical) particles exhibit some similar effects, such as an enhancement of the particle diffusion, their differences in shapes and swimming mechanisms may cause important changes in the transport properties [42]. For example, these models show distinct flow responses, largely owing to the difference in the shapes, resulting in different rheological properties [42,47]. The intrinsic difference in the models can cause strong variations in the nature of near-field HIs as well as in the interparticle torques at close approach, both of which will alter cooperative behaviors at finite density. We will draw attention to these differences occasionally in what follows, in relation to recent simulations of squirmers [22].

B. Basic equations
We now briefly explain how to adapt the FPD method to the case of self-propelled dumbbells. The two real particles in each dimer are represented by a smooth position-dependent viscosity so that where η 0 and η p are the viscosities of the solvent liquid and of the (nearly) rigid particles, respectively. The phantom particle (m = P) does not alter the local viscosity. The function φ α m (r) represents the smoothed profile of the m-th particle of the α-th swimmer, where ξ m is the interface thickness. The rigidity of the real particles is approximately sustained by the large viscosity difference, η p /η 0 ≫ 1.
With η(r) obeying Eq.(1), FPD describes the dynamics by simply solving the usual Navier-Stokes equation for the velocity field: is the viscous stress tensor, p is the pressure, and ↔ δ is the unit tensor. The hydrostatic pressure p is determined by the incompressibility condition In (3), ↔ σ R is the random stress tensor which, in three dimensions, satisfies the fluctuation-dissipation relation as follows where T is the temperature in units of the Boltzmann constant. In (3), F rev (r) is the (volumetric) reversible force density arising from the direct interaction potential between the real particles. With the total potential energy U , whose explicit form is provided below, F rev (r) is uniquely determined as [30,31,35] where Ω α m = drφ α m is the particle volume. Finally, F act (r) is the active force density represented as where δ m,l (m, l = T, P) is the Kronecker δ. Notice that the volume integral of F act (r) is equal to zero, which ensures the conservation of the total momentum. In FPD, the velocity of a particle is defined as the following average value: It is worth noting that without external forces, active forces, or thermal fluctuations, the time derivative of the total energy of the system, which is the sum of U and the kinetic energy To complete the model of a microswimmer, we must introduce suitable dynamical rules for the phantom particles. We assume that such a particle merely follows the motion of the dimer of real particles to which it is attached. In other words, the position R α P (t) obeys where ℓ 1 is the constant center-to-center distance between the tail and phantom particles. We assume that there is no interaction between the phantom particles which therefore can overlap each other. On the other hand, an overlap between the real and phantom particles leads to an unphysical effect: In our model, the real particle is regarded as a rigid particle, and the volumetric force acting on the rigid particle is treated as a homogeneous force density. However, any penetration of the phantom particle into the real particle results in an inhomogeneity of the force density within the rigid particle region, creating a contradiction between the required dynamics and the precepts of the FPD numerical scheme. There are at least two possible ways to avoid any such effect [53]. The first one is to introduce an additional dynamical rule, by which the active force is switched off, when the phantom particle of one swimmer overlaps on the body part of the other. The second one is to introduce additional repulsive interactions involving the phantom particles, which directly prevent such overlaps without losing the propulsive force (thus, in this case, the "phantom" particle is not exactly a phantom). In this study, we mainly investigate the first of these models, but the second will be discussed in the Appendix with some simulation results presented there for comparison. We find that in the present dumbbell model, the qualitative picture regarding the role of hydrodynamic interactions is relatively insensitive to such model details (although they do alter the quantitative results): the nonlocal and longrange nature of HIs allows them to dominate some important aspects of the collective dynamics.
Here, we detail the dynamical rule employed in the main text: when |R α m − R α ′ P | < a + b, where m = (H, T) and α = α ′ , the α ′ -th swimmer becomes passive. That is, if a swimmer's phantom particle overlaps with the realparticle 'body' of another swimmer, the propulsive force pair attached to the first swimmer becomes switched off until its phantom particle once again lies in a purely fluid region. Although somewhat ad hoc, this dynamical rule could reflect a physical situation in which the motility of a bacterium is weakened if its flagellar bundle encounters the body of a second bacterial body. This behavior (a slowing of propulsion in the neighborhood of other particles) will cause an effective attraction between swimming particles [12] and will therefore promote clustering. However, such an effective attraction is already present, even for dynamical rules in which the propulsive force is maintained throughout a collision, since the hard-core repulsions still cause swimming particles to slow down at high densities. The effects of the rule-and collision-induced slowing are both relevant for the clustering as shown in Sec. III C. The important aspect for our study is not the local details of our overlap rule but the fact that, at all times, momentum conservation is still satisfied. This will allow us below to make direct comparisons between the chosen dynamics with and without HIs.
The total potential energy of the system is written as where m =H,T. The first and second terms correspond to the intra-swimmer and inter-swimmer interactions, respectively, between the real particles. The head and tail particles in the same swimmer are stiffly connected by the following harmonic potential: where E 1 is a positive energy constant and ℓ 0 is the natural length of the dumbbell, which is regarded below as the swimmer size. We assume the following form of u: where E 2 > 0 is introduced to prevent the overlap of real particles on different swimmers. The numerical calculations are performed as follows. First, the (off-lattice) particle position R α m (t) at time t is given. Next, from Eqs. (1) and (2), we obtain the on-lattice fluid velocity field at time t + ∆t by solving the Navier-Stokes equation, Eq. (3). Finally, we update the particle position as follows: In implementing our simulations, we first make the equations dimensionless by measuring space and time respectively in units of λ, which is the discretization mesh size used when solving Eq.(3), and τ = ρλ 2 /η ℓ . (The latter represents the momentum diffusion time across the unit length.) A mass unit is then chosen to make λ/τ ,σ = ρ(λ 2 /τ 2 ) and ǫ =σλ 3 the units of velocity, stress and energy, respectively. We set E 2 = 10ǫ and E 1 = 3.2 × 10 5 ǫ. Because E 1 ≫ E 2 , ℓ 0 is nearly constant for the particle densities studied here. In addition, the temperature is, unless stated otherwise, set as T = 0.125E 2 . This value ensures that the effective hard-core particle radius set by Eq.(13) stays close to a and, in combination with subsequent parameter choices, sets the swimmers' Péclet number to a reasonable value (see Sec.III C below). To avoid cumbersome notation we use the above-defined units for the rest of the paper, so that the relevant symbols from now on represent scaled variables. Throughout our simulations we choose a H,T = a = 3.2, a P = 0.75a = 2.4, ξ H = ξ T = ξ P = 1, ℓ 0 = 2.5a, and ℓ 1 = 3a. The viscosity ratio is η p /η 0 = 50. The simulation boxes used are of size 64 3 and 128 3 .

A. Flow field induced by a single swimmer
Before proceeding to the many-body results, we investigate hydrodynamic effects in a few-body system, and show that the present model can reproduce the general features of the flow field induced by a single swimmer and by a pair of swimmers. The simulation box used here is L 3 = 128 3 . To assist this comparison we set T = 0 so that the dynamics is deterministic. Figure 2(a) shows the time evolution of the swim speed v(t) when the active force f act is suddenly applied at t = 0. In the inset of Fig. 2(a), the steady-state velocity v 0 is plotted against f act . Here, we define the Reynolds number of a swimmer as Re=ρv 0 ℓ 0 /η 0 , which varies from 0.0024 to 0.24 as f act changes from 0.1 to 10. These are unrealistically large Reynolds numbers for actual bacteria but, so long as they remain well below unity, the resulting physics should not be strongly affected (see [54] for a fuller discussion in the colloidal context).
The results for a passive dumbbell with an external force f ex of the same magnitude exerted on the tail particle (namely, without the phantom particle) are also shown in Fig. 2 (a). The active swimmer is found to reach the steady-state velocity faster than does the passive dumbbell. Because of the dipolar nature of the swimmer, the distorted region of the velocity field in the steady state is smaller than that caused by the passive dumbbell, and therefore, the time necessary to reach the steadystate velocity field via momentum diffusion is shorter. The scaled friction coefficient along the axis is defined as [55]. We findζ || ∼ = 7.8 for active and ∼ = 3.1 for passive particles.ζ || of the active dumbbell is larger than that of the passive one, which is because the phantom particle, behaving as a model flagellum, pulls the fluid surrounding the dumbbell backwards for f act > 0.  In the steady state, the force balances for the head and tail particles are given by Note that these expressions for F sp H and F sp T in the steady state do not depend on the method of decomposing the active force into the head and tail particles. In our model the active force acts only on the tail particle, and because of the front-back asymmetry of the streamlines, an imbalance in the hydrodynamic drag between the head and tail particles arises. That is, the hydrodynamic drag acting on the tail particle is stronger than that of the head particle, |F vis T | > |F vis H |, which is caused by the active force exerted through the phantom particle, resulting in |F sp T + F act | > |F sp H |. This is not the case for a passive dumbbell subject to an external force: there, because of the front-back symmetry of the streamlines, the viscous drag forces acting on the head and tail particles are then the same.

B. Flow field induced by a pair of swimmers
The hydrodynamic pair interactions of swimmers have been extensively studied by many authors in several dumbbell models [39,57,58] and in squirmers [46,48].
We here investigate the HIs between two swimmers in several situations relevant to our subsequent results on many-body suspensions.

Swimming side by side
At t = 0, two swimmers are placed in parallel in the (z = 0) plane, and active forces parallel to the x-axis are suddenly applied. In Fig. 3(a), we plot the trajectories of the head and tail particles of one swimmer for several different initial distances ∆ between the two swimmers. Although both the head and tail particles are initially attracted, they eventually change direction and repel one another. This behavior has been previously observed by several authors [39,48]. A snapshot of the resulting flow field is shown in Fig. 3(d). For a symmetric dumbbell swimmer [23,24], for which the streamlines exhibit front-back symmetry, the hydrodynamic attraction acts equally on both the head and tail particles in the situation considered here. The repulsion we observe is caused by the front-back asymmetry of the streamlines; the hydrodynamic attraction acting on the tail particle is stronger than that on the head particle, and this imbalance results in a torque, causing outward rotation of the swimming direction.
For ∆/ℓ 0 ≫ 1 the characteristic rotation time of the swimmer can be estimated using the usual Stokeslet point-force approximation as follows: Here, we suppose that the two swimmers are placed in the xy plane and swim in parallel in the x direction. The y component of the velocities of the head and tail particles of swimmer 1 created by swimmer 2 are given by where G yx (x, y) = xy/8πη 0 (x 2 + y 2 ) 3/2 is the yxcomponent of the Oseen tensor, and F sp m,2 (m =H,T) and F act m,2 (m =T,P) are the constraining spring and active forces acting on the swimmer 2, respectively. The rotation is given by which behaves as Here, v 0 ∼ f act /ζ || is the swim speed, ζ || is the friction coefficient along the swimming axis of the dumbbell estimated in the previous subsection, and the relations, ℓ 0 = 2.5a and ℓ 1 = 3a, are made use of. Thus, for ∆/ℓ 0 ≫ 1, initially parallel swimmers can keep moving in their original directions for a distance much greater than ℓ 0 . However, extrapolating this result to ∆/ℓ 0 ≃ 1 indicates that during the time τ a = ℓ 0 /v 0 , over which the swimmers travel distances of order their own size, their swimming directions are strongly altered. This is visible in the results for ∆/ℓ 0 = 1.6 shown in Fig. 3(a). In Figs. 3(b) and 3(c), we show the time dependence of the angular velocity around the center of mass and its maximum value, respectively. This tendency is maximized for ∆/ℓ 0 ≃ 1.6; down to this value, decreasing ∆/ℓ 0 enhances the rotation. Our rough estimation of the characteristic rotation ω ∼ (∆/ℓ 0 ) −4 is qualitatively consistent with the simulation results in this range. However, as ∆/ℓ 0 decreases further, this repelling tendency gets weaker, which is evident in the last panel (∆/ℓ 0 =1.0) of Fig. 3(a) and in Figs. 3(b) and 3(c). This happens for two reasons. First, near-field hydrodynamic interactions generically prevent two nearly contacting tail or head particles from changing their separations: this "squeeze-film" lubrication drag becomes stronger for smaller ∆/ℓ 0 . This results in slower rotations, and hence a slower onset of the resulting repulsions. Second, once the tail particles closely approach one another, their direct interparticle repulsion itself prevents further motion. Because of these two effects, when two initially parallel swimmers lie very close to one another, their rotational velocity becomes much smaller than the far-field estimate. As a result they stay close to each other for a time much larger than τ a = ℓ 0 /v 0 .

Passing each other
Suppose that at t = 0 two swimmers lie anti parallel in a staggered coplanar configuration; the active forces are now suddenly applied. Figure 4(a) shows the trajectories of the head and tail particles of two such swimmers for several different initial positions. At a large enough separation (∆/ℓ 0 ≫ 1), the two swimmers pass by each other without a strong distortion of their trajectories. However, with decreasing ∆/ℓ 0 , the HIs drastically influence their mutual dynamics. When the two swimmers are approaching, they initially repel each other, because the repulsive HIs act on them. Then, when they are passing each other, the head particle of one swimmer is strongly attracted to the streamline induced by the other swimmer, which is also evident in the behavior of the swim- Here, the initial position of the tail particle of the swimmer 1 is set to the origin; that is, ∆x = x − R 1 T,x (0) and ∆y = y − R 1 H,y (0). In (b) and (c), the time dependence of the angular velocity (angles measured in radians) around the center of mass and its maximum value are shown, respectively. For ∆/ℓ0 ≫ 1, the swimmers hardly change their swimming directions while traveling a distance of their own size (∼ ℓ0). However, with decreasing the separation between the swimmers, stronger distortions of the swimming directions are found: In the present case, the hydrodynamic attraction acting on the tail particle is stronger than that on the head particle, and this asymmetry is enhanced for smaller ∆, which leads to faster rotation. However, as ∆ decreases further (∆/ℓ0 < ∼ 1.6), this repelling tendency gets weaker, which may be ascribed to the near-field hydrodynamic and steric effects (see the text for discussion). (d) Snapshot of two swimmers swimming side by side for fact = 10 and ∆/ℓ0 = 1.6 at t ∼ = τa. At t = 0, they are parallel. We also show the velocity field v(x, y, z = 0), where the axes of the two swimmers lie in the plane. The color bar represents the magnitude of the fluid velocity scaled by the value of v0 at fact = 10 shown in the inset of Fig. 2(a). ming speed shown in Fig. 4(b). Therefore at small ∆/ℓ 0 , the speed of the swimmer first reduces and then accelerates. A snapshot of the swimmers at ∆/ℓ 0 = 1.6 is shown in Fig. 4(c).
Similarly to the analysis presented in Sec. III B 1, by using the far-field Stokeslet point-force approximation, for ∆/ℓ 0 ≫ 1, we can roughly estimate the rotation rate as (v 0 /ℓ 0 )(ℓ 0 /∆) 4 . The passing time is on the order of ℓ 0 /v 0 . Thus, for ∆/ℓ 0 ≫ 1 , the original swimming directions are not much changed during the encounter. On the contrary, an extrapolation of this far-field estimation to a near-field ∆/ℓ 0 ∼ 1 predicts a strong realignment of the swimming directions. This is certainly confirmed by our simulationsas shown in Fig. 4.

C. Dynamics in the semi-dilute regime:
Activity-induced clustering So far, we have investigated one and two-body systems and shown that when the distance between two swimmers is comparable to their sizes, the motion of the swimmer is strongly influenced by the flow field caused by the other swimmer. Here, we investigate how such near-field HIs affect the many-body dynamics by simulations in semidilute regimes with thermal fluctuations. The simulations presented in this subsection contain N = 320 swimmers (960 real and phantom particles) in a simulation box with a size of L 3 = 128 3 , and thus the volume fraction of real particles is Ψ = N (Ω H +Ω T )/L 3 =0.052. In the next section, we will discuss the density dependence of the motility in a smaller simulation box (L 3 = 64 3 ). In a randomly distributed state, the average distance between the swimmers is approximately (L 3 /N ) 1/3 ∼ 2ℓ 0 . For comparison we have made equivalent simulations without HIs, using the same parameters for the particle interactions and temperature; without HIs the friction coefficient of a particle is set to ζ || /2 (whose value was evaluated in Sec. III A). This gives the same value of the swim speed of an isolated swimmer both with and without HIs.
In the presence of thermal fluctuations, the active Péclet number is defined as Pe 0 = v 0 τ 0 R /ℓ 0 (so that Pe 0 ∼ ζ || v 0 ℓ 0 /T ∼ f act ℓ 0 /T ), taking values 8.4 and 16.8 for f act = 7 and 14, respectively, as used in the following simulations. Here, τ 0 R ∼ η 0 ℓ 3 0 /T is the rotational relaxation time of the swimming orientation measured in bulk. Although large Pe 0 suggests that active effects dominate over thermal fluctuations, it should be noted that in dilute solution rotational Brownian motion is essential to allow reorientation of the swimming direction. Another mechanism involves the "tumbling" of bacteria at random intervals, causing sudden non-Brownian reorientation; this is not in our model which therefore describes only "smooth swimming" bacteria. For these, the chosen Péclet numbers are within the achievable range [59]. Notice also that Pe 0 and v 0 are values for an isolated swimmer; in the semidilute regime, due to many-body interactions, a marked slowdown occurs, which will be discussed below.
The work done by the active force is transformed into the kinetic energy of the solvent and of the swimmers, which is eventually dissipated by viscosity. In the present simulation, the increase of the average kinetic energy density of the solvent due to the active force is at most approximately 1% of the prescribed value given by the equipartition law, in which the velocity of a thermal fluid obeys v i (r)v j (r ′ ) = T δ ij δ(r − r ′ )/ρ [60], where i and j denote Cartesian components of v.
Each simulation starts from a randomly distributed state, and finally reaches a steady state with nontrivial correlations. In Fig. 5, we show the radial distribution function in the steady state. Both with and without HIs, in an active suspension (f act = 0), g mm ′ (r) increases with increasing f act , suggesting that the swimmers transiently form clusters. This clustering is more evident from the structure factor given by which is shown in Fig. 6. However, there is a marked difference in the behavior with and without HIs; while with HIs the increase is most pronounced in g TT (r), without HIs it is more apparent in g HH (r). More significantly, the overall clustering tendency is significantly enhanced by the addition of HIs. These different behaviors should directly reflect the difference in the clustering mechanism with and without HIs. In order to further understand this difference, we consider the following pair distribution function, where cos −1 (n α ·n α ′ ) is the relative angle between α th and α ′ th swimmers, and the angular integral of g mm ′ (θ, ψ) gives g mm ′ (r). Figures 7(a) and 7(b) represent the steady-state distribution of g mm ′ (θ, ψ) around the first peak of g mm ′ (r) (1 < r/2a < 2), from which we find that neighboring swimmers have different configurations with and without HIs.
A key role in the clustering in the system without HIs can be ascribed to collision-induced stalling of the swimmers when the head of one particle bumps into another. From Fig. 6, we find that the present (on/off) dynamical rule further promotes clustering: once swimmers are inactive, the time needed to separate them from their neighbors becomes longer, which promotes further accumulation. Notice, however, that rule-induced inactivity by itself cannot lead to clustering, because Brownian dumbbells (f = 0) do not cluster without attractive interactions. Clustering (and indeed phase separation) can however result from a density dependence of the mean activity [12,13]. Though not a dominant effect at the present density of Ψ = 0.052, at higher volume fractions the proportion of the inactive swimmers becomes more significant, resulting in a steep reduction of the motility at what is still relatively modest volume fraction. This point will be discussed in the following section. For a perfectly head-to-head collision (θ = π) the stalling continues until Brownian motion changes the swimming directions. But such collisions are rare; the dominant effect is from those where the swimming directions are roughly at right angles, giving a larger cross section for collisions but only temporary stalling. The result is a sharp peak of g HH (r, θ) and g HT (r, θ) at r ∼ = 2a and θ ∼ π/2.
On the other hand, when HIs are included, these mediate noncontact interparticle forces. As shown previously, these create attractions at large distances, but at intermediate separations cause torques that cause swimmers to separate. At short distances, however, these torques are suppressed and a parallel swimmer in close contact will remain so for a long period. (For antiparallel configurations, there is also a long-lifetime state in which each phantom particle lies on top of the tail particle of the other swimmer, so propulsion is suppressed; unlike for parallel alignment, the existence of this state depends on the chosen propulsion rule.) Such mechanisms of HIinduced trapping should still be effective for the semidilute regime studied here. The observed broad peaks in g mm ′ (r, θ) at θ < ∼ π/2 corresponds to nearly aligned configurations in parallel. On the other hand, the peak in g TT (r, θ) at θ > ∼ π/2 reflects antiparallel configurations. These results confirm that the activity-induced clustering is driven by different mechanisms depending on whether in the steady state at fact = 14 with (a) and without (b) HIs. The purple, green, and yellow curves are the contours corresponding to 1/6, 3/6, and 5/6 of the maximum value of g mm ′ (r, θ), respectively.
HIs are present. Without them, particles spend brief but finite periods in close contact due to stalling of rectilinear trajectories; with HIs, they spend longer periods in contact because of trapping in parallel (and, with the chosen swimming rule, also antiparallel) configurations.
This difference in the clustering mechanism is responsible for the contrasting steady-state properties of clusters, such as their lifetime and the average size, which are directly visible from the sequential simulation snapshots of the swimmers as shown in Fig. 8. When identifying clusters, the α th and α ′ th swimmers are considered to be connected if for at least one combination of m and m ′ , where we set δR = 0.6a which is approximately the first-peak width of the radial distribution function. This choice of δR is rather arbitrary, but a small change of its value does not essentially change our conclusions. In Figs. 8(a) and 8(b), we show the typical time evolution of clusters with and without HIs, respectively. Without HIs transient clusters form due to the mixed effects of the collisionand rule-induced stalling, but with the present parameters, the persistence time of contacts is approximately the time τ a = ℓ 0 /v 0 for an isolated swimmer to move its own length. Clusters therefore remain weak and transient. On the other hand, in the system with HIs, once the swimmers are close enough they remain in boundstate configurations for a time much longer than (τ a ), resulting in stronger and longer-lived clusters. We emphasize that many-body hydrodynamic effects further reduce the average swimming speed of our dumbbells, defined as the instantaneous projection of their velocity onto the swimming direction. In our semidilute system (Ψ = 0.052), the average swim speed is evaluated as 0.41v 0 and 0.78v 0 with and without HIs, respectively. Without HIs, about 20% of the swimmers are inactive on average; in this case the reduction of the swim speed is largely owing to the dynamical rule. However, a further large reduction is induced by HIs: with these nonlocal effects switched on, about 35% of the swimmers are inactive on average. This difference in the slowing-down is also seen in Fig. 9, where the average time evolution of the separation distance of a pair of swimmers is plotted for various initial values. This clearly shows that with HIs, prolonged trapping of swimmers becomes more marked than without HIs. This hydrodynamic trapping stimulates a further accumulation of the swimmers and subsequently leads to an enhancement of the cluster growth, which is evident in Fig. 8(c), where we show the probability distribution for numbers in a cluster with and without HIs. Although clusters are longer-lived with HIs, they do finally break up: recall that the hydrodynamic torque tends to misalign the swimming directions and eventually a cluster will find a configuration where this effect dominates long enough to cause its disintegration.
The tendency of HIs to promote clustering is strongly linked to the near-field hydrodynamics, especially in its effects on rotation of the swimming velocity. These effects can be expected to vary with particle shape and swimming type, and indeed our observations are almost the opposite of what happens in (neutral) squirmers [22], where HIs act to suppress rather than enhance the formation of density inhomogeneities.

IV. DENSITY DEPENDENCE OF THE MOTILITY
In the previous section we showed that, for the swimmers at the chosen volume fraction Ψ = 0.052, there is a distinct difference in the swimmers' dynamics, and especially their clustering mechanism, with and without HIs. As the volume fraction increases, the effects of the many-body interactions (and also of the dynamical rule) on the dynamics should be more pronounced, and thus the dynamics is expected to be significantly slower at high densities. Here, we investigate the volume fraction dependence of the mean motility variables (swim speed and rotational relaxation time) and discuss the influence of these dependences on the clustering process. We use the same simulation parameters as in Sec. III C except for the box size, which is here L 3 = (64) 3 . Note that we checked that the essential results for Ψ = 0.052 shown in Sec. III C are barely affected by this change.
A key observation on increasing the volume fraction far beyond 0.052 is that the chosen dynamical rule (switching off activity on overlap of phantom and real particles) has an increasingly pronounced effect at high densities. This motivated the simulations on an alternative model (described in the Appendix) to clarify the role of HIs on the collective behavior. The comparison of these approaches suggests that in the present dumbbell model the details of the local interactions do not qualitatively alter the physical role played by the HIs. For while there are significant near-field hydrodynamic effects which the local interactions clearly do alter, the main role of HIs arises at intermediate and large distances where they can directly affect the collective dynamics.
In Fig. 10(a), we plot the dependence on the volume fraction Ψ of the average swim speed, which is here defined by where (V α H + V α T )/2 is the velocity of the center of mass of the α th swimmer, · · · represents the time average in the steady state. Without HIs, upon increasing Ψ to moderate values, v s (Ψ) decreases. This is largely because of the chosen dynamical rule, which causes an increase in the number of inactive swimmers at high density: at Ψ = 0.31, nearly 90% of phantom particles overlap with the real particles of other swimmers and thus become inactivated. (Note that the details of this decrease will depend on the swimmer's shape.) On the other hand, with HIs, v s (Ψ) drops sharply even at rather small Ψ, in a regime where there are much fewer inactivation events due to overlap. This reflects the long-range character of HIs, which allow each swimmer's motion to be strongly influenced by that of its rather distant neighbors, promoting collective motions even without direct collisions. Such many-body hydrodynamic effects seemingly increase the mean drag force, resulting in the steep decrease of v s (Ψ).
is the center of mass of the α th swimmer. The colors represent the scaled value of ∆R α,α ′ (t). With HIs, the initially adjacent pairs [∆R α,α ′ (∆t = 0) < ∼ ℓ0(= 2.5a)] stay close to each other for a longer time than τa, although without HIs, such swimmers separate themselves much quicker. A strong effect of HIs is also seen in the rotational dynamics, as characterized by the relaxation time of the orientational correlator which can be fit to an exponential form in the whole range of Ψ studied here. In Fig. 10(b), the Ψ-dependent rotational relaxation time, τ R (Ψ), is shown. Without HIs, τ R is nearly constant for low Ψ( < ∼ 0.2), but increases to exceed the Brownian relaxation time [τ R (Ψ ∼ = 0)] at higher Ψ( > ∼ 0.2). This increase is because steric hindrance to rotation increases with crowding at higher Ψ; even in equilibrium (f act = 0), deviations from the ideal-gas behavior due to the excluded volume effects become pronounced for Ψ > ∼ 0.2 (not shown here). On the other hand, upon adding HIs, a striking change in the behavior of τ R (Ψ) is found: first a sharp decrease as Ψ rises towards about 0.05, and then a clear upturn at higher Ψ. Both features can be ascribed to the cooperative nature of HIs. As was discussed in the previous section, HIs among the swimmers generate torques. At small Ψ( < ∼ 0.05), such hydrodynamic torques should enhance rotational relaxation which would otherwise rely solely on Brownian motion. This extra rotation results in the sharply decreasing trend in τ R (Ψ) for small Ψ. (Notice that this effect depends on both the rodlike shape and the front-back asymmetric streamlines around the swimmer within the present model.) However, with crowding, the mean separation between swimmers decreases, which should suppress reorientation while enhancing the hydrodynamic attractions. This presumably causes the reversal of the decreasing trend in τ R at around Ψ ∼ 0.05. Increasing the volume fraction further (Ψ > ∼ 0.15), one expects activity-induced flow fields to increasingly be suppressed, reducing motility and increasing τ R . Nonetheless, due to the incompressible nature of HIs, neighboring swimmers push and pull one another even without direct contact, causing the translational and rotational dynamics to become slower than in the system without HIs. These differences in the volume fraction dependence of the motility with and without HIs affect the physics of clustering. In Fig. 11, we show the structure factor for various Ψ with and without HIs. Without HIs, the clustering tendency increases with volume fraction until Ψ ∼ 0.25. On the other hand, with HIs, this increase continues only until Ψ ∼ 0.15. These clustering tendencies were also checked by direct visualization.
In recent theories for active Brownian particles (ABPs), it was predicted that, if the swim speed decreases steeply enough with density, local densification induces a further accumulation of particles, eventually leading to unstable growth of fluctuations and phase sep-aration [4,12,13]. This mechanism for activity-induced phase separation was confirmed by subsequent simulations on self-propelled spherical particles with rotational diffusion and no HIs [19]. Although the present dumbbell swimmers can be categorized as ABPs, we have shown that the physics of mutual slowing down depends on particle shape, so that even without HIs there is no exact relationship to these earlier simulations. However, the underlying theoretical analysis has some quite generic features [4,12,13] so it is useful to discuss the clustering behavior observed here in similar terms.
As shown in Fig. 10, without HIs and for Ψ < ∼ 0.2, v s (Ψ) decreases almost linearly with Ψ, while the rotational relaxation time τ R (Ψ) decreases only slightly from the value set by purely Brownian rotation. Within the framework of the theory [4,12,13], the effective compressibility S(0) should increase with Ψ, diverging at a spinodal set by the condition dv s /dΨ < −v s /Ψ. This condition holds for Ψ > ∼ 0.15 [evaluated from the fitting of v s (Ψ) in Fig. 10]. In this evaluation, dv s /dΨ + v s /Ψ is found to have a minimum (< 0) at Ψ ∼ = 0.23. The compressibility increase is certainly seen for Ψ < ∼ 0.25 as shown in Fig. 11(b), but beyond this point, rather than diverging, S(0) falls again even though the spinodal condition is satisfied by the observed v s (Ψ). This suggests that for the dumbbell system at this density, direct repulsion via excluded volume interactions (which are omitted from the theory) are strong enough to overcome the effective attraction caused by the density-dependent motility. However, we have performed additional exploratory simulations without HIs (in a simulation box of L 3 = 128 3 ) using a larger self-propulsion force; these suggest that at Ψ = 0.2 and 0.25 phase separation does occur at high enough values of f act .
A broadly similar scenario holds in the presence of HIs, as seen in Fig. 11(a). However the maximum compressibility is significantly greater, and furthermore arises at significantly lower density (Ψ ≃ 0.15). The peak coincides with the state of most pronounced clustering as observed by direct visualization. The data for v s (Ψ) reported in Fig. 11(a) again suggest that the spinodal condition dv s /dΨ < −v s /Ψ of [4,12,13] is satisfied for Ψ > ∼ 0.06, yet once more we observe a maximum, not a divergence, of S(0) near the state of maximal clustering at Ψ ∼ = 0.15. [Note that dv s /dΨ + v s /Ψ now has a minimum (< 0) at Ψ ∼ = 0.12; above this density, the activityinduced attraction itself becomes weakened.] The subsequent decrease might be attributable to steric repulsions, but their effects would have to be stronger than in the case without HIs to start dominating at this somewhat lower density. An additional factor could be the effect of HIs on reducing τ R , which is the mechanism for avoidance of phase separation discussed by Fielding [22]. We have performed additional exploratory simulations with HIs [at the density Ψ = 0.15 in a simulation box of L 3 = (128) 3 ] using a larger self-propulsion force; the structure factor S(k) at small k is found to be dramatically enhanced by formation of a space-spanning structure of the order of the system size, which suggests that bulk phase separation occurs for large enough f act . Notice, however, that the formed structure is not compact and exhibits large fluctuations, which may be due to HIs and the relatively small system size of the present simulation. To firmly assess whether bulk phase separation occurs, we need to perform simulations with larger system sizes.
Thus we have shown the observed clustering tendency of our dumbbells to be broadly consistent with a theory in which collisional interactions are represented at meanfield level by a density-dependent swim speed [4,12,13]. However, using the numerically obtained v s (Ψ) shown in Fig. 10, this theory predicts bulk phase separation in a regime where we see only clusters. (This applies both with and without HIs.) Although steric repulsions offer an obvious candidate mechanism for this discrepancy, the detailed conditions for bulk phase separation in the present dumbbell model are not well understood, and a further systematic study will be needed to clarify them.

V. SUMMARY AND DISCUSSION
In this study, we have numerically investigated the effects of HIs on the collective dynamics of active suspensions, modeling each micro-organism as a stroke-averaged dumbbell swimmer with a prescribed force dipole (and an overlap rule). Our results can be summarized as follows.
With HIs, when the separation distance between the swimmers is comparable to their sizes, their swimming motions are strongly influenced by one another. The activity-induced flow field creates both attractive and repulsive effects that depend on the relative positions and orientations of the swimmers. These effects are significantly more complex than could be expected for simpler models of spherical particles (such as squirmers) in which the nonhydrodynamic forces between particles are spherically symmetric. Although in the far field swimming dumbbells attract, at intermediate separations the HIinduced torques tend to disalign the swimming directions causing swimmers to move apart. Then, at closer distances, these repulsive effects are suppressed; hydrodynamic forces promote preferred structures, in which the swimmers are close to parallel alignment (or, with our overlap rule, antiparallel). These states can last much longer than the intrinsic characteristic time for an encounter, τ a = ℓ 0 /v 0 . Our simulation results strongly suggest that such activity-induced HIs should be significant at rather smaller volume fraction (Ψ ∼ 5%), where they cause strong enhancement of a clustering tendency that is already apparent, albeit in a weaker form and by a different mechanism, when HIs are absent.
Our simulation results for the density dependence of motility parameters show distinctly different behaviors with and without HIs. At a relatively small volume fraction, with HIs, the translational swimming motion becomes slower, while rotational diffusion becomes faster; these are natural consequences of the far-field and intermediate-range pair interactions just described, and are stronger than the corresponding effects without HIs. On the other hand, at large enough volume fractions, the translational and rotational motions are both strongly suppressed. With crowding, because less space is available around a swimmer, the ability of each force dipole to set up a coherent propulsive flow field is diminished, weakening both the motility and the HIs. However, due to the incompressible nature of HIs, neighboring swimmers push and pull each other around even without direct contact. This impedes both the translational and rotational motions much more than is the case without HIs, where only direct collisional forces act between the swimmers.
The different volume fraction dependences of the motility with and without HIs influence the tendency to form clusters. For our chosen interaction and propulsion parameters, clustering is maximal at an intermediate volume fraction in each case, but this volume fraction is lower, and clustering stronger, with HIs than without. The conditions under which bulk phase separation arises, rather than just enhanced density fluctuations, remain somewhat unclear; in particular steric or other interactions at volume fractions above Ψ ≃ 0.15 (with HIs) or Ψ ≃ 0.2 (without) apparently suppress phase separation under conditions where the observed dependence of the mean swim speed v s on Ψ might lead one to expect it [4,12,13]. Such corrections to the theory have also been reported for spherical self-propelled particles [19], but not at such low densities.
In considering the effects of HIs on collective dynamics, we note in addition important differences between alternative types of models. In a recent simulation of twodimensional squirming disks [22], it was found that the activity-induced phase separation is strongly suppressed by HIs. In that system, because the slip velocity is prescribed at their surfaces, when two squirmers come very close, they generally tend to strongly rotate and repel. Thus, crowding enhances the reorientation of the swimmers and reduces their tendency to trap one another. In contrast, in our model of force-prescribed dumbbells, although the reorientation is enhanced at low enough volume fractions, at higher ones both the translational and rotational motions become much slower than without HIs. Moreover, at a certain intermediate volume fraction, the hydrodynamic trapping effect (and indeed the clustering) is most pronounced. Further simulation studies and analysis will be required to clarify the detailed conditions for phase separation with and without HIs in our dumbbell model.
In previous works by Graham and co-workers [23][24][25], it was found that, in suspensions of force-prescribed dumbbell pushers, HIs among them cause coherent fluid motions with a correlation length larger than the swimmer size. Interestingly, such correlated motions are hardly seen in suspensions of pullers. Saintillan and Shelley [41][42][43] reported similar results by simulating a similar model, where the active particle is modeled as a stress-prescribed rigid rod. Furthermore, they also found that an isotropic or aligned homogeneous state becomes unstable due to HIs and some non-linear hydrodynamic effects lead to steady-state structure with a pronounced density inhomogeneity and local orientational order of swimmers [61]. Although these previous studies do not contradict our simulations, we found that, at certain intermediate densities, the clustering caused by activity-induced attractions (or hydrodynamic trapping) is greater than previously reported, which, as already discussed, can be attributed to the effects of near-field HIs. The advantage of the hybrid simulation methods such as that used in the present study is that, by directly taking the solvent dynamics into account, we can accurately treat HIs among finite-sized particles at moderate separations, and also capture some significant near-field effects. In the previous studies, on the other hand, active particles are essentially represented as pairs of point forces and HIs are simply evaluated by using the usual Stokeslet. However, such an approach based on far-field hydrodynamics is not able to fully capture the collective hydrodynamic effects in relatively dense active suspensions. Thus, our present study should complement the earlier works cited above.
Our dumbbell simulations include thermal fluctuations, which act like a repulsion to prevent the swimmers from clustering. Similar to an analysis performed in Sec. III B 1, we can use the Stokeslet point force approximation to estimate the energy scale of the dipolar hydrodynamic forces as f act ℓ 0 /10, where the interaction range is assumed to be from the core size (= 2a) to the mean separation distance of the swimmers (∼ 2ℓ 0 = 5a). Thus the ratio of this energy scale to the temperature is of order Pe 0 /10 where Pe 0 is the Péclet number introduced previously. By this rough estimate, thermal and hydrodynamic effects can compete in our simulations. (We have done further exploratory simulations without thermal fluctuations, and found that the clustering is certainly enhanced when these are absent.) Some quantitative aspects of our results depend on the nature of the local interactions between swimming dumbbells, in particular the rule that causes propulsion to be switched off when swimmers overlap. However, by varying this rule (see the Appendix) we have found a robust role for long-range hydrodynamic interactions in controlling the collective behavior. The HIs promote clustering through a collective slowing of propulsive motions at densities that are too low for the rules governing direct interparticle collisions to be their dominant cause. In contrast, without HIs the slowing is collisional and thus depends more strongly on the local interactions. Our finding that HIs have significant effects on clustering physics is far from obvious; due to the dipole-dipole nature of these interactions among self-propelled particles, they are much weaker than for particles subjected to external forces, and under many conditions their effects appear almost negligible [62]. Certainly one can argue that for situations where the main interaction between particles involves well-separated two-body collisions, the effects of HIs are indeed relatively weak; however this argument is dangerous when applied to the many-body collective behavior of swimmers at the fairly high densities actually arising within clusters.
We conclude with some further remarks on topics that we hope to study in future: (i) Collective hydrodynamic effects should strongly depend not only on the shapes of the swimmers but on their self-propulsion mechanism. We have done preliminary studies for pullers (f act < 0); these show quite different steady-state properties from the pushers studied in this paper. At an equal volume fraction of Ψ = 0.052, the near-field radial distribution functions (for separations r < ∼ ℓ 0 ) are significantly smaller than those of the equilibrium (f act = 0) dumbbells. This indicates that for pullers repulsive interactions dominate at these separations, in contrast to the case of pushers. Accordingly, a very different density dependence of the motility (and therefore different clustering and/or phase separation properties) from that of pushers is expected to arise, whose details remain to be explored.
(ii) In recent experiments on E. coli in the presence of additional attractive forces (created via a depletion potential due to polymer additives) it was shown experimentally and by simulation that activity produces a significant shift of the phase boundary compared to that of a passivated system with the same attractions [63]. However the configurations favored by such an attraction need not coincide with those stabilized by the activity-induced HIs. Therefore the case of active dumbbells with attraction requires separate investigation. One very recent study suggests a mechanism whereby the equilibrium phase separation caused by attractions is interrupted by activity-induced cluster breakup [64], but equilibrium attractions could equally interfere with the motility-induced phase separation mechanism discussed in Section IV.
where a H = a T = a and a P = b. Note thatŨ is given as the sum of U [Eq. (11)] and the repulsive interactions involving the phantom particles. Since R α P is not an independent variable but is uniquely determined by R α H and R α T according to Eq. (10), the force is given by where the second term of the last line is due to the additional repulsive interactions, whose explicit forms are Here, is the component of ∂Ũ /∂R α P perpendicular to the swimmer's axis, where ↔ δ is the unit tensor andn α is the unit vector along the swimming axis of the α th swimmer. The physical meaning of these additional forces ∆F α,m (m =H,T) is as follows. The tail and phantom particles interact via hypothetical stretching and bending potentials, which transmit a repulsive force acting on the phantom particle; in the tight-binding limit of these interactions, this force is immediately transformed into the translational and rotational forces acting on the head and tail particles (body part) according to Eqs. (A4) and (A5). A schematic of this situation is shown in Fig. 12. Note that the present treatment also conserves the total momentum. In this Appendix, similarly to the analysis performed in Sec. IV, we investigate the volume fraction dependence of the mean motility and its influence on the clustering behavior with the above-introduced additional interactions (and without the dynamical rule). Here, we use the same parameters and box size as those in Sec. IV.
In Figs. 13(a) and 13(b), we plot the dependence on the volume fraction Ψ of the average swimspeed, v s (Ψ), and the rotational relaxation time, τ R (Ψ), respectively. Similarly to the result shown in Sec. IV [ Fig. 10(a)], with HIs, v s (Ψ) exhibits a much faster decay than that without HIs. The rotational relaxation also exhibits a similar behavior to that shown in Fig.10(b). These observed differences in the (mean) motility for systems with and without HIs can be again attributed to the long-range and cooperative natures of HIs (see the discussion in Sec. IV). Notice that the reduction of τ R (Ψ) at lower Ψ is more significant than that observed in Fig. 10(b). This may be because of both the addition of the repulsions and the removal of the dynamical rule, which enhance the collision-induced transfer of the angular momentum.
As is discussed in Sec. IV, these differences in the volume fraction dependence of the motility with and without HIs should affect the clustering behavior. In Fig. 14, we show the structure factor for various Ψ with and without HIs. Without HIs, the weak clustering occurs at smaller volume fraction (Ψ < ∼ 0.2) but is highly suppressed at higher Ψ. From the data of v s (Ψ), the condi- The red and green curves are for the cases with and without HIs, respectively. These curves exhibit similar behaviors to those with the dynamical rule shown in Fig. 10 but with an overall shift to higher Ψ (due to the absence of the rule-induced slowing down).
FIG. 14: (Color online) The structure factor with (a) and without (b) HIs for various Ψ at fact = 14. With HIs, within the range of Ψ investigated here, the structure factor at small k(= L/2π) increases with increasing Ψ. tion dv s /dΨ < −v s /Ψ holds for Ψ > ∼ 0.3, so that according to the theory [4,12,13], the effective compressibility is expected to be negative, leading to the instability of the homogeneous state. However, as noticed in Sec. IV, for the present dumbbell system, at such relatively large volume fractions, direct repulsions (which are omitted in the theory) give a dominant contribution to the compressibility (or pressure). This stabilizes the homogeneous state. Although the structure factor is found to grow significantly on increasing the magnitude of the active force, further numerical investigations are necessary to clarify the detailed conditions for bulk phase separation, which remain the subject of a future study. On the other hand, with HIs, as seen in Fig. 14(a), the clustering becomes enhanced with increasing Ψ. The fitting of v s (Ψ) represented in Fig. 13(a) suggests that the spinodal condition dv s /dΨ < −v s /Ψ of Refs. [4,12,13] is satisfied for Ψ > ∼ 0.1. Notice that the clustering tendency at relatively larger volume fractions is much more enhanced than that exhibited by the model with the switch-off rule [for comparison, see Fig. 11(a)]. In the model used here, the active forces of the swimmers are not reduced at close proximity, and therefore the activity-induced hydrodynamic attractions are strong enough to overcome repulsive interactions even for larger Ψ.
In this Appendix, by performing an analysis similar to that in Sec. IV, we have explored the clustering behavior of the dumbbell model with a different treatment of the phantom particle. Without HIs, some marked differences arise between the two treatments; that is, the collective behaviors are strongly dependent on the local rules or interactions. On the other hand, because the long-range and non-local natures of HIs govern the collective dynamics, the global picture regarding the hydrodynamic effects is not much altered. This robustness to the local interaction rules is an important result of the present study. The above simulation results, along with those shown in Sec. IV, show that theory can qualitatively explain the characteristic features of the effective attraction caused by the density dependent motility, although there are some discrepancies between the theory [4,12,13] and the present simulation, especially at higher volume fractions, where some effects omitted from the theory (for example, steric repulsions due to the excluded volume) can come to dominate. Although the effects of long-range HIs are largely separable from near-field and collisional effects, the latter are separately important, at least at the quantitative level. Here one expects strong dependencies both on particle shape and on whether the propulsion is described by a force-prescribed or a velocity-prescribed mechanism.
of the particle velocity [Eq. (9)]. In the Brownian simulations, the active force is multiplied by this factor to match the swim-speed of an isolated swimmer with and without HIs.