Hydrodynamic interactions in anomalous rheology of active suspensions

We explore a mechanism of the anomalous rheology of active suspensions by hydrodynamic simulations using model pusher swimmers. Our simulations demonstrate that hydrodynamic interactions under shear flow systematically orient swimmers along the extension direction, which is responsible for determining the global swimming states and the resulting significant viscosity reduction. The present results indicate the essential role of hydrodynamic interactions in the elementary processes controlling the rheological properties in active suspensions. Furthermore, such processes may be the substance of the previously proposed scenario for anomalous rheology based on the interplay between the rotational diffusivities and the external shear flow.


I. INTRODUCTION
Anomalous rheology observed in the broad class of active suspensions is one of the most typical phenomena highlighting distinctive differences from passive systems [1][2][3][4][5][6][7][8][9][10] .In particular, for rod-like extensile pusher microswimmers (such as E. coli), a significant viscosity reduction has been experimentally observed at lower shear rates and volume fractions [2][3][4][5][6], which frequently leads to a superfluid state with zero viscosity [4,6].A seminal study by Hatwalne et al. [1] predicted that if an orientational order along the extension axis of the applied flow is somehow realized, the active dipolar forces intensify the mean flow, reducing the resistive stress required to drive the external flow and thus the viscosity.Following that, many theoretical attempts have been made to predict or explain the anomalous rheology in active suspensions (see papers [11][12][13][14][15][16][17][18][19] and the references therein).
In dilute suspensions of rod-like particles, the orientational distribution of particles under shear flow are known to be enhanced along the extension axis when (thermal or athermal) random rotational diffusion processes exist [20,21].By taking such fluctuation effects into account, the viscosity reduction was successfully modeled within the framework of continuum kinetic theory [13,15].
In dilute/semidilute active suspensions, hydrodynamic interactions (HIs) are expected to play a crucial role in couplings among constituents [22,23].In Ref. [16], it is theoretically demonstrated that the long-range HIs induce marked density fluctuations that provide additional sources of the effective rotational noise, resulting in a decrease in the viscosity.Indeed, recent experiments [6] indicate a close link between collective many-body properties and anomalous rheology.Nevertheless, due to the highly nonlinear and nonequilibrium nature of HIs, our understanding of the extent to which interactions among swimmers are involved in the rheological properties is still lacking beyond the effective one-body theory.
In this study, we investigate the mechanism of the anomalous viscosity reduction observed in active sus-pensions by revisiting the role of HIs.Our analysis, along with a phenomenological explanation, elucidates that swimming along the extension axis of the applied flow is hydrodynamically favorable, resulting in a significant reduction of the viscosity.Furthermore, we argue that in usual swimming bacteria, such as E. coli, the self-propulsive forces are strong enough that the induced HIs can compete or dominate other effects like thermal fluctuations even in dilute suspensions.

II. MODEL SWIMMER SYSTEM
For the present purpose, we perform hydrodynamic simulations of model active suspensions composed of N rod-like dumbbell swimmers with a prescribed force dipole.Our model swimmer, schematically shown in Fig. 1(a), is composed of body and flagellum parts.The body part is treated as a rigid body, while the flagellum part is regarded as a massless "phantom" particle simply following the body's motions.This treatment always keeps the relative position of these two parts unchanged.For the α-th swimmer (α = 1, • • • , N ), it is assumed that the force F A nα acting on the (front) body is exerted by the (rear) flagellum and that the flagellum also exerts the force −F A nα directly on the solvent fluid.Here, nα is the direction of the α-th swimmer, and these forces compose a dipolar force (please refer to Appendix A for details).The present particle-base model is essentially the same as those proposed in Refs.[24,25] and used in Refs.[13,16,[27][28][29].Continuum kinetic models of hydrodynamically interacting rod-like swimmers with prescribed stresses or forces were also developed [13,16,[30][31][32][33].In Refs.[30][31][32], it was demonstrated that nonlinear hydrodynamic effects can lead to larger-scale correlated motions with marked density fluctuations.
As illustrated in Fig. 1(a), the body and flagellum parts are assumed to have the same shape and are each described by a superposition of three spheres with a common radius R. The spheres composing the body are located at the positions R α − ℓ 0 nα is the position of the center of the flagellum.The shape of the present model swimmer shows the headtail symmetry, and the mid point is thus given by R α = (R (G)  α + R (CF ) α )/2.Although arbitrary shapes of swimmers with an imposed head-tail asymmetry can be composed, we may obtain qualitatively the same results as long as these swimmers have rod-like forms with the prescribed force dipoles.
Periodic boundary conditions are imposed in the xand y-directions with the linear dimension L, and the planner top and bottom walls are placed at z = H/2 and −H/2, respectively.The shear flow is imposed by moving the top and bottom walls in the x-direction at constant velocities V /2 and −V /2, respectively, whereby the mean shear rate is γ = V /H.This situation is illustrated in Fig. 1(b).Hydrodynamic interactions among the swimmers are incorporated by adopting the smoothed profile method (SPM) [34][35][36], which is one of the mesoscopic simulation techniques [37][38][39][40][41][42] .In the SPM [34][35][36], the dynamics of rigid particles and a host fluid can be considered simultaneously with vastly reducing numerical costs by replacing particle-fluid boundaries with smoothed ones and by taking particle rigidity into account through the body force term in the Navier-Stokes equation.The details of the simulation methods are presented in Appendix A and B.
body flagellum FIG. 1: (Color online) (a) Our model swimmer comprises "body" and "flagellum" parts with symmetric shapes.Each parts are constituted by a superposition of three spheres with radius R. We assume that the force FA nα is exerted on the body, while −FA nα is directly exerted on the solvent through the flagellum part, with nα being the orientation of the α-th swimmer.These forces constitute a force dipole with the magnitude FAℓ0.Here, ℓ0 is the characteristic swimmer's length and, for the present model, is given as the separation distance between the body and flagellum centers.(b) The periodic boundary conditions are imposed in the x-and y-directions with the linear dimension L. The shear flow is imposed by moving the top and bottom walls in the x-direction at constant velocities V /2 and −V /2, respectively.These two walls are separated in the z-direction by the distance H. First, in Fig. 2, we show the viscosity η for various conditions.In this study, the viscosity is defined as where Σ xz (x, y, ±H/2) is the xz component of the stress tensor at the walls and ⟨• • • ⟩ hereafter denotes taking the time average in a steady state.Here, the solvent viscosity is scaled to be 1.At a relatively low shear rate γ, we find that η takes lower values than the solvent viscosity, which qualitatively agrees with the experimental results [2][3][4][5][6].This behavior strongly depends on H and the volume fraction of the swimmers defined as ϕ = N V (b) /L 2 H with V (b) being the volume of the body part.The viscosity can be divided into three parts: η = η s + ∆η p + ∆η a , where η s (= 1 in this study) is the solvent viscosity, ∆η p and ∆η a are the passive and active contributions, respectively (see Appendix A for more details).In the present framework, ∆η a is given as where nα (t) is the unit vector representing the α-th swimmer's orientation at time t, and nα,x and nα,z are its x and z components, respectively.Essentially identical expressions of Eq. ( 2) were previously derived (see Refs. [10,13,15] for example).From Eq. ( 2), when swimmers tend to align along the extension direction of the flow field (⟨n α,x nα,z ⟩ > 0), ∆η a < 0. Since the contribution of ∆η p to η is positive in general, a significant decrease in the viscosity occurs from the negative ∆η a .
To further explore what swimming states are involved in the viscosity reduction, we investigate the following steady-state quantities: Here, R α (t) is the center of the force dipole (̸ = the center-of-mass position), ρ(z) is the denisty, and p(z)/ρ(z) and ↔ W (z)/ρ(z) represent the polarization vector and the nematic order parameter tensor, respectively [8].These quantities, which depend only on z at steady state, are shown in Figs.3(a)-(h) for various conditions.In Figs.3(a) and (b), ρ(z) has significant peaks near the boundary walls, and otherwise, it is almost constant, indicating that the walls attract swimmers.Such behaviors were already reported and discussed in the literature (for example, Refs.[22,23,[43][44][45][46][47][48][49][50]).In the present model, without thermal fluctuations, when placing one swimmer near the wall, it continues to swim along the wall, which suggests that the force-dipole prescribed to the swimmer contributes to the wall attraction [43].However, in the many-swimmer case, significant disturbances are induced by interactions among the swimmers.Such disturbances produce an outgoing flux from the wall to the bulk region.Meanwhile, self-propulsive motions give incoming flux to the wall from the bulk.Competition between these two flux terms should determine the amount of accumulation of swimmers at the walls [49,50].
Figures 3(c)-(f) show p(z)/ρ(z).Due to the flow and geometrical symmetries, p y (z) = 0 for all z.For ϕ ≳ 0.03, swimmers trapped at the walls tilt their "heads" to the walls [51][52][53].Moreover, the tilting angle is greater for larger ϕ, which may be caused by HIs among the swimmers on the wall.These issues will be studied elsewhere.
In terms of the viscosity reduction, among Figs.3(a)-(h), of particular interest are (g) and (h), exhibiting W xz (z)/ρ(z).At ϕ = 0.01 and H = 128, where the viscosity reduction is absent (see Fig. 2), W xz (z) ≲ 0 as a whole.For larger ϕ and H, in contrast, W xz (z) > 0 for all z.Within the present range of ϕ and H, by increasing these parameters, the upward convex form of W xz (z)/ρ(z) tends to grow.Equation ( 2) is rewritten as through which W xz (z) is directly related to the viscosity reduction.
The reduced viscosity immediately indicates the reduced shear rate at the walls.Here, we briefly review this behavior.For a swimmer with nα,x nα,z > 0, the active-force reinforces the applied flow.More specifically, in a small region including the swimmer, the velocity gradient is intensified, while in outer regions, the opposite happens.A superposition of such contributions gives the net effects on the mean flow, and we observe a lower shear  rate at the walls, γw , in exchange for a greater shear rate in the interior region.These situations are schematically illustrated in Figs. 4 (a) and (b).Consequently, as shown in Fig. 4(c), the shear stress required to maintain the applied shear rate γ = V /H is reduced [1], and the viscosity is given by which is smaller than η s for γw < γ.Such a modulation of the velocity field accompanying with the viscosity reduction was certainly observed in experiments of E. coli suspensions [6].Notably, in contrast to active (pusher) suspensions, dispersed particles in the usual passive system suppress the velocity gradient in the interior region, and the observed viscosity is increased.

B. Key role of hydrodynamic interactions in determining the global swimming states
Then, we investigate how such observed global steady states are realized.To this end, taking the flow and geometrical symmetries into account, it is useful to classify one-swimmer states into the following four states (see Fig. 5): state 1, (n α,x > 0, nα,z > 0), state 2, (n α,x < 0, nα,z > 0), state 3, (n α,x < 0, nα,z < 0), and state 4, (n α,x > 0, nα,z < 0).States 1(2) and 3(4) are equivalent; that is, they can be converted to each other by simply rotating the coordinate frame about the y-axis by π.Figures 3(g) and (h) with Eq. ( 3) indicate that, when ∆η a < 0, states 1 and 3 with nα,x nα,z > 0 are realized more favorably than states 2 and 4 with nα,x nα,z < 0. Thus, we may further classify states 1 and 3 into the favorable F -state and states 2 and 4 into the unfavorable U F -state.Now, the question is why states 1 and 3, which contribute to W xz > 0, are more favorable than states 2 and 4. In our simulations, thermal effects are absent, and excluded volume effects are almost irrelevant because the suspensions mainly considered here are dilute.Instead, hydrodynamic effects are expected to determine the overall swimming states.We expect swimmer's motions to be largely disturbed by HIs even without approaching the contact distance to other swimmers; we regard such events as hydrodynamic collisions.We support this per- spective by analyzing the transition probabilities between the swimming states: we pick up a pair of swimmers whose separation at time t is less than a certain close distance d 0 ; then, the transition probabilities are determined by comparing their states at t − ∆t and t + ∆t.
In this study, we set d 0 = 0.7ℓ 0 and ∆t = 1.25τ 0 .Here, τ 0 = ℓ 0 /v s is the time to travel the distance of the swimmer size (∼ ℓ 0 ), with v s being the average swimming speed.In the present range of ϕ, the average distance between neighboring swimmers, l N = (ϕ/V (b) ) −1/3 is 2 ∼ 3 times larger than d 0 .Although quantitative evaluations of the transition probabilities significantly depend on d 0 and ∆t, the qualitative discussion presented below is not affected as far as d 0 and ∆t are sufficiently smaller than l N and l N /v s , respectively.In Appendix C, we discuss how the present definition can capture hydrodynamic collisions with the settings of d 0 and ∆t.Table I shows the numerically obtained probability of the µ-state, P µ , and the transition probability from the µ-to ν-states, W µ|ν , (µ, ν=F, U F ) for various conditions.Here, P µ and W µ|ν are calculated for swimmers in the region −0.3 ≤ z/H ≤ 0.3.At γ = 0, because the F -and U F -states are not distinguished, W F |U F ∼ = W U F |F , and P F ∼ = P U F .However, for γ ̸ = 0, we find that W F |U F is significantly smaller than W U F |F .As H and ϕ increase (in the dilute regime), the population of the Fstate swimmers with nα,x nα,z > 0 increases, indicating that an increase in the collision frequency or time further promotes transitions.
For swimmers trapped at the walls, the hydrodynamic torques arising from the applied flow weakly align them along the flow direction because their heads are slightly tilted against the walls.Thus, for trapped swimmers, the population of the F -state is slightly larger than that of the U F -state.After longer-term traps, swimmers leave from the bottom (top) wall by raising (dropping) their heads, which changes their states (F ⇄ U F ). Reflecting such conditions, for swimmers just after leaving the walls, the population of the U F -state is slightly larger than that of the F -state (not shown here).As swimmers move inward from the boundary walls, transitions from the U F -to F -states are gradually promoted by collisions.Due to the geometrical symmetry, the population of the F -state is maximized at z = 0, leading to the upward convex form of W xz (z)/ρ(z).In an ideal bulk system or a system with periodic boundary conditions without walls, a detailed balance between the F -and U F -states, P F W F |U F = P U F W U F |F , should be realized.Such a detailed balance may nearly hold at larger H and ϕ in the present system, but that was not investigated in detail.Tables II and III show the numerically obtained transition probabilities of a swimmer in states 1 and 2 before a collision, respectively, at ϕ = 0.01, H = 384, and γ = 10 −3 .Here, w (ν) λ|µ represents the transition probability from states λ to µ through a collision with another swimmer in state ν.Note that similar results are obtained at different parameters where negative ∆η a is obtained.We find significant differences between w (ν) 1|µ and w (ν) 2|µ .For both cases the majorities are w (ν) µ|µ , whereas a swimmer in state 1 is more likely to retain its state unchanged by a collision than one in state 2.
We can understand the role of HIs in the elementary processes of these state transitions through the following phenomenological arguments, which are separately provided for different cases.Let us first consider that two swimmers in the F -and U F -states are approaching.There are essentially two different cases: case (A) where swimmers in states 1 and 2 (equivalently, 3 and 4) are approaching, and case (B) where swimmers in states 3 and 2 (equivalently, 1 and 4) are approaching.
For case (A), as schematically shown in Fig. 6(a), HIs tend to rotate the swimmers in opposite directions [23], while the externally applied shear flow rotates them in the same direction.As the swimming directions become parallel to each other and perpendicular to the flow direction, the torques due to HIs grow weaker, but those arising from the shear flow grow stronger.Furthermore, once two swimmers move nearly side by side (Fig. 6(a)), a hydrodynamic attraction acts on them, which may make a swimmer in state 1 drag one in state 2 into eventually moving in the same direction in the collision process.These hydrodynamic effects are expected to promote the transition from states 2 to 1, responsible for w In contrast, in case (B), due to similar asymmetry in the net torques, w    2|3 is less notable than that between w 2|1 : namely, hydrodynamic collisions of case (A) predominantly contribute to the transition from the U F to F states, whereas those of case (B) are marginal.

Hydrodynamic collisions between two swimmers in the same states
Here, we consider the following two cases: case (A') where two swimmers are both in state 1 (equivalently, both in state 3), and case (B') where those are both in state 2 (equivalently, both in state 4).For these cases, schematics are shown in Figs.7(a) and (b).
For both cases (A') and (B'), the torques caused by HIs rotate the swimmers in opposite directions and grow weaker as the swimmers become parallel to each other.On the other hand, for the torques caused by the shear flow, in case (A'), as the collision proceeds, the torque resisting the transition to state 2 grows stronger, whereas the other torque, which helps the transition to state 4, grows weaker.In (B'), the opposite occurs: one torque due to the shear flow promoting the transition to state 1 grows stronger, while the other one resisting the transition to state 3 grows weaker.The difference in how the torques contribute to the transition is expected to be responsible for the measured difference in the transition probabilities.That is, as shown in Tables II and III, w

Hydrodynamic collisions between two swimmers in states 1 and 3, and 2 and 4
When two approaching swimmers are in states 1 and 3, there are essentially two different cases, (A"1) and (A"2), which are illustrated in Figs.8(a) and (c), respectively.As torques arising from HIs rotate the two swimmers in opposite directions, the swimming directions become perpendicular and parallel to the flow in cases (A"1) and (A"2), respectively.In case (A"1), the torques due to the shear flow, which prevent the swimmers from changing from the F -to U F -states, grow stronger.On the other hand, in case (A"2), such torques promoting the changes to the U F -state grow weaker.When swimmers in states 2 and 4 are approaching, there are two different cases, (B"1) and (B"2), as schematically shown in Figs.8(b) and (d), respectively.With similar arguments, in case (B"1), the torques due to the shear flow, which promote the transition to the F -state, grow stronger, while, in case (B"2), those preventing the changes to the F -state grow weaker.

These differences may result in the difference between w
(3) 1|µ and w (4) 2|µ .That is, when two swimmers in states 1 and 3 are approaching each other, the swimmers tend to remain their states unchanged more than when they are in states 2 and 4: w 2|2 .Note that hydrodynamic collisions of cases (A"1) and (B"1) predominantly contribute to the transition from the U F to F states, whereas those of cases (A"2) and (B"2) are marginal.
The realistic collision processes are more complicated; thus, the present arguments are oversimplified.However, they qualitatively explains why HIs systematically promote the transition from the U F to F states.

IV. DISCUSSION AND CONCLUDING REMARKS
It has been known that rod-like particles in a shear flow tend to orientate to the extension direction due to an interplay between flow and rotational diffusivities [20,21]: for a rod-like particle, although the torque due to shear flow becomes unidirectional and stronger as its orientation becomes perpendicular to the flow direction, the torque due to thermal rotational diffusivities is bidirectional and does not depend on the rod orientation.By considering such an effect, an explanation for the viscosity reduction was provided [13,15].Moreover, it was proposed that the activity-induced HIs provide a source of random rotations in addition to thermal fluctuations and tumbling [16].The present study further illuminates the role of even starting from a random state, our results suggest that steady global states where the swimmers are weakly aligned along the extension axis may form as self-organization by repeated hydrodynamic collisions.A study of this issue would be an interesting task for future studies.
In typical microorganisms systems, the propulsive forces are sufficiently strong that hydrodynamic effects may dominate over thermal fluctuations.Below, we validate this condition by considering a typical experimental situation [6]: an E. coli suspension at a volume fraction ϕ(= V (b) /l 3 N ) =0.01 at room temperature (∼ 300K), for which the average separation distance is l N ∼ 5µm and the thermal rotational diffusion coefficient is D T ≲1s −1 .Hereafter, we assume that the swimming speed is v s ∼ 10µm/s, the magnitude of the force dipole is P ∼ 10 −18 N•m, the cell size is ℓ 0 ∼ 1µm, the cell volume is V (b) ∼ 1µm 3 , and the solvent viscosity is η s ∼ 10 −3 Pa•s.For a duration ∼ 1/D T ≳ 1s, a swimmer may at least once approach another swimmer closer than l c ∼ 2µm, estimated by πl 2 c v s × (1/D T ) × (1/l N ) 3 ∼ 1.The magnitude of the rotational flow field, ω, induced at a distance r from a swimmer is approximately given as ω ∼ 0.1P/(η s r 3 ) [23].Therefore, at r = l N , ω ∼ 1s −1 , while at r = l c , ω ∼ 10s −1 .By such a hydrodynamic "collision" process, which lasts for approximately ℓ 0 /v s ∼ 0.1s, swimming motions can be largely affected more than by thermal fluctuations.In other words, reorientation due to HIs may be a faster process than thermal rotational diffusion.
In this study, we have explored a mechanism of the anomalous rheology of active suspensions, focusing on the role of HIs.Before closing, we present the following remarks.(1) Our pusher model is transformed into a puller model by simply changing the sign of the active forces.Our preliminary results shown in Fig. 2(b) suggest that the viscosity of the puller model is increased more than that in the passive systems, which agrees with experimental observations for motile and immotile puller bacterial suspensions [7].(2) In Ref. [5], under Poiseuille flow, lower viscosity is observed for smaller separation between the walls.This contrasts with the present result, where the viscosity reduction is enhanced by increasing H under simple shear.In addition, the viscosity reduction occurs at larger shear rates than those of the experiments of Refs.[4,6].These differences may be attributed to the difference in the flow geometry.We plan to investigate these issues further elsewhere.
where i, j = 1, 2, 3 and µ, ν = b, f .Here, U µν is the interaction potential between two spheres which each comprise the body or the flagellum part of different swimmers, and W µ is the interaction potential between such a sphere and the planar wall.The explicit forms of U µν and W µ are provided below.In Eqs.(A10) and (A11), F α,ext and N α,ext are the force and torque exerted on the α-th swimmer due to the external field, which are absent in the present study.The active force acting on the body part, F α,A , is given as Eqs. (A9) and (A16) prescribe a force dipole [see also Eq. ( A19)].Finally, F α,H and N α,H are the force and torque exerted on the α-th swimmer due to HIs.The explicit forms of F α,H , N α,H , and the body force f H can be given in the discretized equations of motion as Eqs.(B6), (B7), and (B9), respectively in the next section.
We assume the following form of the interparticle potential: where ϵ is a positive energy constant and δ µ,f is the Kronecker delta.This form prevents the body part of a swimmer from overlapping on different swimmers but allows overlaps among the flagellum parts.The wall-particle interaction potential W µ is introduced to prevent the penetration of particles through the boundary walls and is assumed to be given as where we assume the same energy constant as that of U µν .In Eqs.(A17) and (A18), µ, ν = b, f .In our simulations, we make the equations dimensionless by measuring space and time in units of h, which is the discretization mesh size used when solving Eqs.(A6)-(A8), and t 0 = ρh 2 /η s , which is the momentum diffusion time across the unit length.Accordingly, the scaled solvent viscosity is 1, and the units of velocity, stress, force, and energy are chosen to be h/t 0 , ρh 2 /t 2 0 , ρh 4 /t 2 0 and ρh 5 /t 2 0 , respectively.In our simulations, we set ϵ = 30 and F A = 20.The parameters determining the swimmer's shape are set to be R = 3.2, 5R and ξ = 0.5.In this study, the swimmers' volume fraction is identified as that of the rigid body particles given by ϕ = N V (b) α /HL 2 .In Ref. [35], the general scheme deriving the volumeaverage stress tensor, ↔ S, in the framework of the SPM is provided: ↔ S is divided into the three parts, ↔ s s , ↔ s p , and ↔ s a , due to the solvent, passive, and active contributions, respectively.The passive part is further divided into two parts arising from HIs and the potentials (U µν and W µ ).Such sources to the stress tensor exist without active forces, so we call ↔ s p the "passive" stress.In this Appendix, according to a similar procedure given in Ref. [35], we derive an expression of the active stress as follows. ↔ As usual, we may redefine the active stress by making S are divided into three parts, and then, we have η = η s + ⟨s p,xz ⟩/ γ + ⟨s a,xz ⟩/ γ.The passive part < s p,xz > / γ positively contributes to the viscosity, and therefore, the viscosity reduction is entirely due to (weak) alignment of the force dipoles along with the extension direction, ⟨n α,x nα,z ⟩ > 0. In the main text, we also discuss how such an alignment of the swimmers in the interior regions reduces the velocity gradient at the boundaries, which is observed as a viscosity reduction.
We may define the swimmer's Reynolds number as Re = ρv s ℓ 0 /η s with v s being the average swimming speed.In the present simulations, v s ∼ 0.1, giving Re ∼ 1, which is unrealistically large, but may not be problematic for the following reason.It is known that unless Re is much greater than unity, the induced flow patterns around a moving particle do not change enough for the inertia effect to significantly change the resultant dynamics and transport properties.In Ref. [38], an excellent discussion is presented for colloidal simulations.In practice, by using a smaller value of ρ, it is possible to reduce Re as much as possible.For this treatment, a smaller time increment is required to ensure the stability of the time integration of the equations of motion, whereas it leads to longer simulation times.In our preliminary simulations at Re ∼ 0.1, which is still unrealistically large, we confirm that the main results obtained in the present study remain almost unchanged.

Appendix B: Explicit time-integration algorithm
Here, following Refs.[34,35], we describe an explicit time-integration scheme for solving model equations as follows.
The set of physical variables is assumed to be clearly defined at the discrete time step t n = n∆t.
First, we solve Eq. (A6) without including f H as

FIG. 2 :
FIG. 2: (Color online) (a) The γ-dependence of the viscosity η for various ϕ at H = 128.For smaller γ and larger ϕ, η becomes smaller.The inset shows H-dependence of η for two different ϕ at γ = 10 −3 .In the present range of H, η is significantly smaller for larger H.(b) The ϕ-dependent viscosity for active and passive rod-like particles at H = 128 and γ = 10 −3 .Our preliminary results for passive (FA = 0) and active-puller (FA = −20) cases show that η increases with increasing ϕ; the viscosity enhancement is much stronger in the active-puller case.In both (a) and (b), the dashed lines indicate the solvent viscosity ηs(= 1).

FIG. 4 :
FIG.4: (Color online) (a) For a swimmer with nα,x nα,z > 0, the surrounding velocity gradient is intensified, while away from the swimmer, the opposite occurs.(b) The net flow is determined by a superposition of the individual swimmers' contributions.(c) The x component of the velocity field averaged over the steady state at ϕ = 0.032 and H = 256.At the walls, the shear rate is lower than γ as γw = (η/ηs) γ < γ, in exchange for a greater shear rate in the interior region.(d) Wxz(z)/ρ(z) at the same condition as (c).

FFIG. 5 :
FIG. 5: (Color online) Fourfold classification of (oneswimmer) swimming states.These four states are further classified into the favorable F -and the unfavorable U F -states.

1 .
Hydrodynamic collisions between two swimmers in states 1 and 2, and 2 and 3

FIG. 6 :
FIG. 6: (Color online) Schematic illustrations for two approaching swimmers in the F -and U F -states: case (A) where the two swimmers are in states 1 and 2 (a), and case (B) where those are in states 2 and 3 (b).The torques due to HIs and the shear flow are described by the black and green curved arrows, respectively.
this case, the torques both due to HIs and the shear flow are reduced as their swimming directions become parallel along with the flow direction (see Fig.6(b) for a schematic), and therefore, the difference between w

FIG. 7 :
FIG. 7: (Color online) Schematic illustrations of the two cases: case (A') where the two swimmers are both in state 1 (a) and case (B') where those are both in state 2 (b).The torques due to HIs and the shear flow are described by the black and green curved arrows, respectively.

2 FIG. 8 :
FIG. 8: (Color online) Schematic illustrations of cases (A"1) and (A"2), shown in (a) and (c), respectively, where the two swimmers are in states 1 and 3, and those of cases (B"1) and (B"2), shown in (b) and (d), respectively, where the two swimmers are in states 2 and 4. The torques due to HIs and the shear flow are described by the black and green curved arrows, respectively.

TABLE I :
State probabilities Pµ and the transition probabilities W µ|ν for various conditions(10 2 ϕ, H, γ) PF PUF W F |F W F |U F W U F |F W U F |U F

TABLE III :
Transition probability from states 2 to µ through a collision with another swimmer in state ν, w