Microrheology of active suspensions

Abstract

We study the microrheology of active suspensions through direct hydrodynamic simulations using model pusher-like microswimmers. We demonstrate that the friction coefficient of a probe particle is notably reduced by hydrodynamic interactions (HIs) among a moving probe and the swimmers. When a swimmer approaches a probe from the rear (front) side, the repulsive HIs between them are weakened (intensified), which results in a slight front-rear asymmetry in swimmer orientation distribution around the probe, creating a significant additional net driving force acting on the probe from the rear side. The present drag-reduction mechanism qualitatively differs from that of the viscosity-reduction observed in sheared bulk systems and depends on probing details. This study provides insights into our fundamental knowledge of hydrodynamic effects in active suspensions and serves as a practical example illuminating distinctions between micro- and macrorheology measurements.

---

1 Introduction

In active suspensions, the intrinsic activity of swimming particles leads to distinctive collective behaviors and transport/rheological properties deviating from those observed in passive particle suspensions.^1–3 A striking example is anomalous rheology.^4–23 In particular, in suspensions of pusher-like microswimmers, such as E. coli, the viscosity significantly decreases,^5–9 frequently establishing zero or negative viscosity states.^7,9,12 The underlying mechanism behind such anomalous rheology involves weak orientational order along the extension axis of the externally applied flow field, which could be attributed to rotational diffusivities and/or hydrodynamic interactions (HIs).^{11,13,16,18,22} The active dipolar forces with the orientational order intensify the mean flow, reducing the resistive stress required to drive the external flow and consequently diminishing the viscosity.⁴ The anomalous viscosity reduction in active suspensions contrasts with the viscosity behavior of dilute suspensions of passive particles, which is well described by the Einstein viscosity formula.²⁴

The local viscosity or viscoelastic properties of active suspensions are of great interest; probing rheological properties at the μm level can provide further insights into the underlying mechanisms and enable a more detailed characterization of the anomalous rheology. When probing smaller-scale rheological properties of complex fluids or soft materials, microrheology measurements are considered powerful tools (see recent reviews^25–27 and the references therein). By tracking the motions of small probes suspended in fluids, typically at a μm scale, microrheology allows for measuring viscoelastic properties across a wider range of temporal scales than conventional macrorheology techniques. This finer resolution in both time and length scales may offer a more comprehensive understanding of the material's rheological properties. However, microrheology still faces several unresolved issues. One such issue is whether the observations made using microrheology accurately reflect the actual local rheological properties of the material in the absence of probes. The interactions among probes and suspended constituents or inner structures can significantly influence the motion of the probes and potentially alter the local rheological properties around the probe. Understanding and accounting for these effects are crucial for accurately interpreting microrheology measurements and relating them to a material's intrinsic viscoelastic properties.

Microrheological investigations on active suspensions have provided various intriguing results.^28–40 A simulation study by Foffano et al.^32,33 demonstrated that a probe particle can experience a negative viscosity in active nematics composed of contractile-puller-type swimmers. More recent numerical studies^34–36,38 have predicted that the friction coefficient of a probe can be reduced in active Brownian particle baths without HIs, but not to a level smaller than that defined in the absence of active/passive Brownian particles, indicating that the measured viscosity is larger than the solvent viscosity. However, as described above, in bulk rheological experiments of extensile-pusher-type swimmers, it has been observed that the measured viscosity can be lower than the solvent viscosity and even approach zero. This raises a question of whether such a phenomenon can be replicated in microrheology measurements when HIs are taken into account. If so, this would prompt further investigation into whether the reduction of the friction coefficient of a probe shares a similar mechanism to the viscosity reduction observed in macrorheology.

2 Simulation methods

For the present purpose, we conduct direct hydrodynamic simulations using a model of active suspensions comprised of N rod-like dumbbell swimmers. Our model swimmer, shown in Fig. 1(a), is essentially the same as the one employed in our previous investigation of macrorheology simulations.²² In this model, each swimmer consists of a body and a flagellum: the body is treated as a rigid-body particle, while the flagellum is considered a massless “phantom” particle that simply follows the body's motions. This treatment maintains the relative position of the body and flagellum parts. For the α-th swimmer (α = 1, …, N), we assume that a force F_An̂_α acting on the (front) body is exerted by the (rear) flagellum and that the flagellum also exerts the force −F_An̂_α directly on the solvent fluid. Here, n̂_α is the direction of the α-th swimmer, and these forces compose a dipolar force of magnitude F_A[small script l]₀, with [small script l]₀ being the characteristic swimmer's length [see Fig. 1(a)]. Such a force-prescribed particle model emulates rod-like pusher-type microorganisms such as E. coli, as initially proposed in the ref. 41 and 42 and employed in subsequent studies.^{16,18,22,43–46}

Fig. 1 (a) Our model swimmer comprises body and flagellum parts with symmetric shapes. Each part is constituted by a superposition of three spheres with radius a. We assume that a force F_An̂_α is exerted on the body, while −F_An̂_α 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 of magnitude F_A[small script l]₀, with [small script l]₀ being the characteristic swimmer's length, which for the present model is given as the separation distance between the body and flagellum centers. In this study, to incorporate the present model swimmer into the SPM, the body and flagellum parts are represented through field variables, Ψ^(b)_α(r) and Ψ^(f)_α(r), respectively. In the bottom panel, we plot Ψ^(b)_α(r) + Ψ^(f)_α(r) in the xy-plane, where both R^(G)_α = (2.25a,0,0) and R^(CF)_α = (−2.25a,0,0) are included. (b) In our microrheology simulation, a single probe particle with radius R is immersed in a fluid. The probe particle is also described by the field variable Ψ_p(r). In the bottom panel, we plot Ψ_p(r) representing the probe particle in the xy-plane, where the probe center is set to R_p = (0,0,0). In (a) and (b), the discretized mesh size h is the same as that used in practical simulations (h = 0.3125a and ξ = 0.5h). Here, ξ is the interface thickness controlling the degree of smoothness of Ψ^(b)_α(r), Ψ^(f)_α(r), and Ψ_p(r). (c) A single probe particle immersed in a fluid is dragged by a constant force F_ex along the x-direction. The periodic boundary conditions are imposed in all directions with the linear dimension L.

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 a. The spheres composing the body are located at the positions R^(b)_i,α = R^(G)_α + (2 − i)an̂_α (i = 1, 2, 3), where R^(G)_α is the α-th swimmer's center-of-mass position. Similarly, the spheres composing the flagellum part are located at R^(f)_i,α = R^(CF)_α + (2 − i)an̂_α (i = 1, 2, 3), where R^(CF)_α = R^(G)_α − 4.5an̂_α = R^(G)_α − [small script l]₀n̂_α is the position of the center of the flagellum. The shape of the present model swimmer shows the head–tail 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 can obtain qualitatively the same results as long as these swimmers have rod-like forms with the prescribed force dipoles.

In our simulations, we use the smoothed profile method (SPM)^47–49 to accommodate many-body hydrodynamic interactions (HIs) among the constituent swimmers. The SPM can accurately reproduce even the near-field HIs:⁴⁹ in the present study, the ratio of the simulation mesh size to the particle size a is on the order of 0.1, determining the spatial resolution of HIs. With this setting, the SPM can quantitatively replicate the near-field HIs (lubrication interactions) to a satisfactory degree up to a closer distance as h/a ∼ 0.1 (or even less),⁴⁹ where h represents the gap distance measured as the separation distance between the interfaces of adjacent particles. However, to fully reproduce the singular divergence of lubrication forces, which become crucial at much closer distances, finer resolutions are required. This situation is similar to other hybrid simulation methods.^50–55 However, such singular lubrication interactions act only on particles with extremely smooth and rigid surfaces, where the degree of surface roughness is significantly smaller than the particle size. Thus, for modeling HIs among microswimmers or those between the probe particle and microswimmers, reproducing singular lubrication interactions is not necessary.

Here, we describe a detailed scheme to simulate the present model swimmer system, essentially the same as that used in a previous study.²² The body and flagellum parts of a swimmer are represented through the field variables Ψ^(b)_α(r) and Ψ^(f)_α(r), respectively:

andIn this study, we adopt the following function to ψ aswhere μ = b, f and ξ is the interface thickness controlling the degree of smoothness. In Fig. 1(a), we show the cross section of the model swimmer described by Ψ^(b)_α(r) and Ψ^(f)_α(r), including both R^(G)_α and R^(CF)_α in the same plane.

In our microrheology simulation, a single probe particle (radius R) is immersed in a fluid, which is dragged by a constant force F_ex along the x-direction. Similarly to the swimmers, as shown in Fig. 1(b), the probe particle is also described by the field variable Ψ_p(r) as

Ψ_p(r) = ψ[r, R_p; R], (4)

where R_p is the position of the probe sphere's center and R is the probe radius.

The working equations for the velocity field  v(r,t) are given as

∇·v = 0. (7)

Eqn (5) is the usual Navier–Stokes equation.²⁴ Here, ρ is the solvent mass density, , given as eqn (6), is the viscous stress tensor with η_s being the solvent viscosity, and the hydrostatic pressure p is determined by the incompressibility condition, eqn (7). In addition, f_H is the body force required to satisfy the rigid-body condition for the swimmer's body and probe particle regions, and f^(f)_A is the active force directly exerted by the flagellum part to the fluid:where is the volume of the flagellum part. In addition, the volume of the body part is given as . In this study, because the shapes of the body and flagellum parts are assumed to be the same, .

For the model swimmers, the equations of motion for the center-of-mass velocity V^(G)_α and the angular velocity with respect to the center-of-mass Ω^(G)_α are

whereandare the mass and the moment of inertia of the α-th swimmer's body, respectively. Here, is the unit tensor and Δr_α = r − R^(G)_α. In this study, the swimmer's density is assumed to be the same as the solvent density. In eqn (9) and (10), F_α,H and N_α,H are the force and torque exerted on the α-th swimmer due to HIs, respectively. The explicit forms of F_α,H, N_α,H, and the body force f_H are given in the ESI.†F_α,int and N_α,int are the force and torque acting on the α-th swimmer's body, respectively, due to the particle–particle and particle–probe potential interactions:where i, j = 1, 2, 3 and μ, ν = b, f. Here, U^μν is the interaction potential between two spheres in which each comprise the body or the flagellum part of different swimmers, and W is the interaction potential between such a sphere and the probe sphere. The explicit forms of U^μν and W are provided below. The active force acting on the body part, F^(b)_α,A, is given as

F^(b)_α,A = F_An̂_α. (15)

Eqn (8) and (15) prescribe a force dipole F_A[small script l]₀n̂_α with [small script l]₀n̂_α = R^(G)_α − R^(CF)_α [see also eqn (34)].

Similarly, for the probe particle, the equations of motion for the center-of-mass velocity V_p and the angular velocity with respect to the center-of-mass Ω_p are

whereandare the mass and the moment of inertia of the probe particle, respectively. Here, Δr_p = r − R_p, is the volume of the probe particle, and the probe particle density is also assumed to be the same as the solvent density. In eqn (16), F_p,int is the force due to the interaction with swimmers:where i = 1, 2, 3 and μ = b, f.

We assume the following form of the interparticle potentials:

where is a positive energy constant, δ_μ,f is the Kronecker delta, and μ, ν = b, f. This form prevents the body part of a swimmer from overlapping on different swimmers but allows overlaps among the flagellum parts. The interaction potential W between the probe particle and the particles constituting a swimmer is introduced to prevent the penetration of swimmers through the probe boundary. In this study, W is assumed to be given aswhere we assume the same energy constant as that of U^μν.

For a direct comparison between simulations with and without HIs, we have made equivalent simulations without HIs, using the same parameters for the interactions. Furthermore, in eqn (9), (10), (16), and (17), the hydrodynamic forces and torques are replaced as

N_α,H → −ζ_RΩ^(G)_α, (24)

F_p,H → −6πη_sRV_p, (25)

N_p,H → −8πη_sR³Ω_p. (26)

Here, the values of the friction coefficients ζ_‖, ζ_⊥, and ζ_R are numerically evaluated as those of an isolated swimmer with HIs; ζ_‖ = 9.6πη_sa, ζ_⊥ = 10.9πη_sa, and ζ_R = 38.7πη_sa³. This approach draws parallels with the conventional frameworks used in Brownian or relaxation dynamics simulations of spherical particles, for which the friction coefficient is usually defined as 6πηa, with η and a representing the solvent viscosity and the sphere radius, respectively.

In our simulations, as illustrated in Fig. 1(c), periodic boundary conditions are imposed in all directions with the linear dimension L. The probe dynamics are influenced by periodic image probes and swimmers; as elucidated below, the probe dynamics are significantly influenced by HIs between the probe and its nearby surrounding swimmers. Therefore, while HIs with distant image swimmers (located farther than the nearest neighbors) will have a small quantitative contribution, they should not alter the qualitative aspects of the simulation results. We make the equations dimensionless by measuring space and time in units of h and t₀, respectively. Here, h denotes the discretization mesh size used in solving eqn (5)–(7), and t₀ = ρh²/η_s 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₀ = η_s/(ρh), ρh²/t₀² = η_s²/(ρh²), ρh⁴/t₀² = η_s²/ρ, and ρh⁵/t₀² = η_s²h/ρ, respectively. Note that ρh⁴/t₀² = η_s²/ρ is the intrinsic force scale of a Newtonian fluid. In our simulations, we set and F_A = 20. The parameters determining the swimmer's shape are set to be a = 3.2, [small script l]₀ = |R^(G)_α − R^(CF)_α| = 4.5a and ξ = 0.5. In this study, the swimmer volume fraction is identified as that of the rigid body particles given by .

3 Results

3.1 Friction coefficients

Fig. 2(a), (c), and (d) present the scaled friction coefficient ζ/ζ₀ for various conditions. In this study, the friction coefficient of the probe particle ζ is defined aswhere V_p,x is the x-component of the probe velocity, 〈⋯〉 represents the time average in a steady state, and ζ₀ is the bare friction coefficient experienced by the probe particle suspended in a pure solvent.

Fig. 2 In (a), (c), and (d), we present the friction coefficients ζ scaled by ζ₀ against F_ex. For the cases with HIs, ζ₀ is the friction coefficient of the probe particle suspended in a pure solvent with F_ex = 10. On the contrary, in the cases without HIs, ζ₀ = 6πη_sR as set in eqn (25). In (a), ζ/ζ₀ is shown for three different conditions at a probe radius R = 20(= 6.25a) and ϕ = 0.032. With HIs, ζ is significantly smaller than ζ₀ for F_ex ≲ 10², contrasting with the behavior in the absence of HIs. We also present a passive case, where we use the same swimmer model but without the active force (F_A = 0). Note that for F_ex ≳ 10³, ζ for active and passive systems almost converges, simply because the influence of the self-propulsion of swimmers is relatively diminished.^34–36 In (b), we plot 〈V_p,x〉, with error bars, for F_ex ≤ 10² both with and without HIs, shown in the upper and lower panels, respectively. In (c), ζ/ζ₀ is shown for different ϕ at a probe radius R = 20(= 6.25a) with HIs. The main data are the results for L = 128, while the yellow-green closed circle represents the data for a larger system size L = 256 at R = 20, ϕ = 0.016, and F_ex = 10, almost corresponding to the case of L = 128 at the same R and ϕ. The upper-right pannel displays 〈δV_p,x²〉 with δV_p,x = V_p,x − 〈V_p,x〉. For relatively small F_ex, 〈δV_p,x²〉 is nearly constant and smaller for lower ϕ. In the lower-right panel, the normalized friction coefficient, (ζ − ζ₀)/(ζ₀ϕ), is plotted. For 3 ≲ F_ex ≲ 10², an almost linear dependence of ζ − ζ₀ on ϕ is shown. In (d), ζ/ζ₀versus F_ex is shown for several different probe sizes (R = 5, 10, and 20) at ϕ = 0.016. At R = 5 and 10, within the examined range of F_ex, the reduction of ζ is less pronounced than at R = 20.

In Fig. 2(a), ζ/ζ₀ is shown for three different conditions for a probe with radius R = 20(= 6.25a) and a volume fraction ϕ = 0.032. When HIs are considered, ζ is significantly smaller than ζ₀ for F_ex ≲ 10². To provide further insights into the role of HIs, equivalent simulations without HIs are performed, employing identical parameters for the particle–particle interactions. Without HIs, ζ is always larger than ζ₀, showing non-monotonic F_ex dependence⁵⁶ similar to that observed in earlier simulation study.³⁴ The observed distinction strongly suggests that HIs play a significant role in the reduction of ζ in active suspensions. Additionally, we present a passive case, where we use the same swimmer model but without the active force (F_A = 0). Notably, for F_ex ≳ 10³, the values of ζ for both active and passive systems with HIs tend to converge with a marked force thickening. This convergence is simply because the influence of the self-propulsion of swimmers is relatively diminished.^34–36 In the ref. 57 and 58, the thickening of microviscosity due to hydrodynamic lubrication is explored in the context of active microrheology in passive suspensions, and it is found that the thickening becomes notable when the attainable gap distance h can be approximately less than 0.1 of the particle size. In our simulation, a similar degree of the resolution of HIs is achieved, which may explain the observed increase in ζ at larger F_ex as reflecting this thickening effect.

Fig. 2(b) shows the average probe velocity 〈V_p,x〉 at a probe radius R = 20(= 6.25a) and ϕ = 0.032, both with and without HIs. For relatively smaller F_ex, strong fluctuations in V_p,x are found, reflecting significant back-and-forth motions due to frequent collisions with surrounding swimmers. However, each data point is derived from simulations conducted over extended periods to ensure the accuracy of 〈V_p,x〉; especially for smaller F_ex values like F_ex = 3 and 10, the simulation ran for a period where the probe particle traveled several hundred times its own size. For further details on related topics, please refer to Fig. 3 and the associated discussion.

Fig. 3 The mean square displacement (MSD) of the probe particle, normalized by R², in the transverse (y- and z-) directions, is plotted for F_ex = 3, 10, and 10² at ϕ = 0.032 and R = 20. The solid lines represent (〈V_p,x〉Δt/R)² with Δt being the elapsed time, for F_ex = 3, 10, and 10², from lower to higher, respectively. The crossover from ballistic to diffusive motions is observed at around Δt ∼ 10³. Notably, 〈Δy_p²〉/R² and 〈Δz_p²〉/R² roughly converge into a single curve, ensuring the axisymmetric behavior of the probe particle motions around the external force direction (x-direction). For F_ex ≲ 10, there exists a period where 〈V_p,x〉Δt is less than the average diffusive displacement of the probe particle in the transverse directions, which is consistent with that the magnitude of fluctuations in V_p,x is larger than 〈V_p,x〉 for the same F_ex range.

In Fig. 2(c), ζ/ζ₀ is shown for different ϕ at a probe radius R = 20(= 6.25a) with HIs. The reduction in ζ is enhanced as the swimmer volume fraction ϕ increases at least within the range examined in our simulations (ϕ ≲ 0.05). The upper-right pannel displays 〈δV_p,x²〉 with δV_p,x = V_p,x − 〈V_p,x〉. For relatively small F_ex, where the average probe velocity 〈V_p,x〉 is considerably smaller than the average speed of microswimmers, 〈δV_p,x²〉 is nearly constant and smaller for lower ϕ. This feature suggests that the fluctuations in probe velocities are determined by the details of probe-swimmer collisions and their statistics. In the lower-right panel, the normalized friction coefficient, (ζ − ζ₀)/(ζ₀ϕ), is plotted. It indicates a linear dependence of ζ − ζ₀ on ϕ for F_ex ≲ 10² within the present range of ϕ. However, this linear relationship does not extend to higher F_ex values. Furthermore, even for F_ex ≲ 3, the expected scaling also breaks. The precise origin of these deviations – whether they arise from statistical inaccuracies or more intrinsic phenomena – remains to be elucidated.

In macrorheology simulations using the same rod-like model swimmer, HIs induce a weak alignment of swimmer orientations along the elongation axis of the applied flow field, resulting in the acceleration of the mean flow and a subsequent reduction in the viscosity.²² Accordingly, we can expect that the decrease in ζ in the present microrheology simulations also reflects a similar mechanism involving some ordering of swimming directions in the bulk region. However, as shown below, this is not the case.

Fig. 2(d) shows ζ/ζ₀ against F_ex for three different probe sizes. For the F_ex ranges at R = 5 and 10, the magnitude of the induced velocity gradient in the bulk region (more than one swimmer size away from the probe surface) approximately corresponds to that observed for F_ex ∼ 10–10³ at R = 20, but the reduction in ζ is not as noticeable as for R = 20. Furthermore, in Fig. 2(c), the main data are the results for L = 128, while the yellow-green closed circle represents the data for a larger system size L = 256 at R = 20, ϕ = 0.016, and F_ex = 10. Although we examine only a single case, a significant system-size dependence is hardly observed, which contrasts with the macrorheology case, where the viscosity reduction shows a strong system-size dependence.²² These findings suggest that the drag-reduction mechanism depends on the local probing details and thus qualitatively differs from that of the viscosity reduction observed through macrorheology measurements.

3.2 Local swimmer states around the probe particle

To elucidate the mechanism behind the reduction in ζ, let us explore the local swimmer states around the probe and their impact on the probe dynamics. Before proceeding, it is important to verify the symmetry around the x-axis, which corresponds to the direction of the applied external force. In Fig. 3, we plot the probe particle's mean square displacement (MSD) along the y- and z-directions for various values of F_ex at ϕ = 0.032 and R = 20. The MSD in these transverse directions is proportional to Δt for larger Δt, where Δt represents the elapsed time, indicating clear diffusive behavior of the probe particle in these directions. Importantly, the diffusion coefficients in the y- and z-directions are nearly identical, and these diffusivities do not show a significant F_ex dependence. Therefore, this observation ensures the axisymmetric behavior of the probe particle around the external force direction (x-direction) and indicates that the probe motion along and across the x-direction does not interfere significantly with each other. The same behavior was also observed for different ϕ (not shown here). In the following analysis, we assume that axial symmetry about the x-axis approximately holds. Note that, with increasing ϕ, collisions with swimmers occur more frequently, resulting in the enhanced probe diffusivity.

In Fig. 4(a)–(f), we show the contour plots of the number density of the swimmers 〈ρ(x,y)〉. At R = 20(≳ [small script l]₀), as shown in Fig. 4(a), (b), (e), and (f), 〈ρ(x,y)〉 near the probe boundary is considerably larger than in the bulk region, irrespective of the presence of HIs. As indicated in the literature, both the steric⁵⁹ and hydrodynamic⁶⁰ effects cause the entrapment of microswimmers on the boundaries. A self-propulsive rod-like particle, upon colliding with the probe, cannot rebound immediately due to the steric/geometric constraints. Instead, it stays near the probe for a while, swimming along the surface. This behavior may result in a larger density near the probe surface. Note that HIs between the probe and the swimmers enhance the entrapment,⁶⁰ but they are not significantly dominant in the present condition. At R = 5(≲ [small script l]₀), as shown in Fig. 4(c) and (d), the peak of 〈ρ(x,y)〉 near the boundary decreases, both with and without HIs, suggesting a weaker entrapment effect for smaller R. Also, as shown in Fig. 4(g) and (h), without HIs, the increase in 〈ρ(x,y)〉 is slightly more pronounced on the front than the rear, but such a behavior diminishes with HIs.

Fig. 4 (a)–(f) Contour plots of the scaled number density 〈ρ(x,y)〉/ρ₀ for various conditions, where ρ₀ is the average density defined as ρ₀ = N/[L³ − 4πR³/3]. The left and right panels correspond to cases with and without HIs, respectively. The white semi-circular curves represent the boundary of the probe sphere. In (a), the swimmer is shown at the same scale as the probe to compare these sizes directly. We assume axial symmetry about the x-axis. At R = 20(≳ [small script l]₀), as shown in (a), (b), (e), and (f), 〈ρ(x,y)〉 near the probe boundary is considerably larger than ρ₀, irrespective of the presence of HIs. On the other hand, for R = 5(≲ [small script l]₀), as shown in (c) and (d), the peak of 〈ρ(x,y)〉 near the surface is less pronounced. In (g) and (h), 〈ρ(x,0)〉/ρ₀ is shown at R = 20, ϕ = 0.032, and F_ex = 10 with and without HIs, respectively. Here, the gray-colored regions indicate 0 ≤ |x| − R ≤ [small script l]₀. (i)–(l) The scaled density 〈ρ〉_s,θ/ρ₀ averaged over the range R < r < R + [small script l]₀ at the polar angle θ for F_ex = 10 and 100 with and without HIs. In (a), the dashed white semi-circular curve represents r = R + [small script l]₀.

The distinction between with and without HIs is illustrated in Fig. 4(i)–(l), where we plot the density, averaged over the range R < r < R + [small script l]₀ and scaled by ρ₀:

Here, 〈ρ(x,y)〉 is reinterpreted as 〈ρ(r,θ)〉 with x = r cos θ and y = r sin θ. As demonstrated, in the presence of HIs, the density at the rear is larger than at the front, whereas in their absence, the reverse is observed. This tendency becomes more pronounced with increasing F_ex.

Fig. 5 shows the orientation distribution of swimmers P([small straight theta, Greek, circumflex]; x, y) around the probe particle. It is defined as

where r = (x, y, 0), x̂ is the unit vector along the x-axis and t̂_α = n̂^‖_α/|n̂^‖_α| with n̂^‖_α denoting the projected swimming direction onto the x–y plane. The normalization factor c(x, y) specific to the present visualization in Fig. 5 is determined so that . The main text presents the results for two cases at F_ex = 10 and 100 with HIs. For cases without HIs, please refer to the ESI.† In Fig. 5, the panels (a)–(e) show P([small straight theta, Greek, circumflex]; x, y) near the probe boundary. The tilt angle of the swimmers (measured with respect to the tangential direction of the probe surface) is larger when they face the probe than when they face the bulk, but this distinction is subtle in (c) and (d) for F_ex = 10. Additionally, more swimmers face the probe for x < 0 than x > 0. These characteristics are more evident at F_ex = 100 than at F_ex = 10.

Fig. 5 The orientation distributions of swimmers P([small straight theta, Greek, circumflex]; x, y) averaged over a square region outlined by dashed lines at R = 20 and ϕ = 0.032 with HIs, with the block size being 2a = 0.32R(= 6.4). For the definition of P([small straight theta, Greek, circumflex]; x, y), please refer to the explanation presented in the main text. The probe-sphere's center is located at (x, y) = (0, 0). In (a)–(f), the left and right panels present the distributions at F_ex = 10 and 100, respectively. The violet and dark-green lines correspond to P([small straight theta, Greek, circumflex]; x, y) (x < 0) and P(180° − [small straight theta, Greek, circumflex]; x, y) (x > 0), respectively, calculated in the regions indicated by identical colored characters and normalized to make the total enclosed area equal 1. The dotted lines guide the normal and tangential directions along the probe sphere.

These observations are further illustrated in Fig. 6, where we plot various components of the following quantities:

Here, , , , and , where r = (x, y, 0) with x = r cos θ and y = r sin θ. Here, n̂_α,r and n̂_α,θ are r- and θ-components of the (unit vector) direction of the α-th swimmer, respectively. Eqn (30) and (31) represent r- and θ-components of the average orientation vector of the swimmers, respectively. On the other hand, eqn (32) and (33) represent rr- and rθ-components of the average nematic order parameter, respectively. These quantities are averaged values per swimmer within the range R < r < R + [small script l]₀ at the polar angle θ. They are evaluated under conditions where R = 20 and ϕ = 0.032, with HIs, for F_ex = 10 and 100. Fig. 6(a) and (b) demonstrate that the swimmers face the probe more likely than facing the bulk, and such a tendency is more enhanced at F_ex = 100 than F_ex = 10. As shown in Fig. 6(c) and (d), at F_ex = 100, more swimmers move along the probe surface from the rear- to the front-sides, while at F_ex = 10, such a tendency is not obvious. Furthermore, we can find that the rr-component of the nematic order is more enhanced in the rear than in the front, and this tendency is more enhanced at F_ex = 100 than at F_ex = 10. Such a front-rear asymmetry in 〈W_rθ〉_s,θ is weaker than that in 〈W_rr〉_s,θ. This difference is essential in the observed front-rear asymmetry in the active stress, as shown in Fig. 7.

Fig. 6 (a) and (b) show 〈m̂_r〉_s,θ defined in eqn (30), while (c) and (d) represent 〈m̂_θ〉_s,θ defined in eqn (31). Here, 〈m̂_r〉_s,θ and 〈m̂_θ〉_s,θ represent r- and θ-components of the average orientation vector, respectively. (e) and (f) show 〈W_rr〉_s,θ defined in eqn (32), while (g) and (h) represent 〈W_rθ〉_s,θ defined in eqn (33). Here, 〈W_rr〉_s,θ and 〈W_rθ〉_s,θ represent rr- and rθ-components of the average nematic order parameter, respectively. These quantities are averaged values per swimmer within the range R < r < R + [small script l]₀ at the polar angle θ. They are evaluated under conditions where R = 20 and ϕ = 0.032, with HIs, for F_ex = 10 and 100.

Fig. 7 Contour plots of the rr- and rθ-components of the active stresses, denoted as 〈Σ^act_rr(x,y)〉 (a) and 〈Σ^act_rθ(x,y)〉 (b), respectively, for F_ex = 10 and ϕ = 0.032. The white semi-circular curves represent the boundary of the probe sphere. In (a), the swimmer is shown at the same scale as the probe to compare these sizes directly. In (a), the inset shows 〈Σ^act_rr(x,y) − Σ^act_rr(−x,y)〉, indicating apparent front-rear asymmetry in 〈Σ^act_rr(x,y)〉. This asymmetry means that the swimmers exert more force from the rear. It is noteworthy to mention the following. If we make the active stress tensor traceless, as is often done in the literature, 〈Σ^act_rr〉 approaches zero with increasing distance from the probe. However, in our definition of the active stress eqn (34), we leave as , resulting in a finite, though nearly constant 〈Σ^act_rr〉 at distances away from the probe. In (b), 〈Σ^act_rθ(x,y)〉 shows more significant fluctuations than 〈Σ^act_rr(x,y)〉.

Hydrodynamic interactions (HIs) cause swimmers to be repelled more significantly when approaching from the front than from the rear. When approaching from the front, HIs between the swimmer and probe tend to push them away from each other. On the other hand, when approaching from the rear, such a tendency is diminished, and for sufficiently large F_ex, the swimmer can be pulled towards the probe depending on factors such as their relative velocity and angle of approach. In this situation, swimmers suffer from significant scatterings at the front side, while scatterings are comparatively weaker at the rear. This difference in scattering causes to the observed asymmetry between the front and rear in both the orientation and the nematic orders of the swimmers. Additional supporting simulations for this observation are provided in the ESI.†

3.3 Driving force exerted on the probe particle by swimmers

The observed local swimmers’ behavior significantly affects the probe motions. To illustrate this, we display the contour plots of the rr- and rθ-components of the active stresses, denoted as Σ^act_rr(x,y) and Σ^act_rθ(x,y), respectively, for F_ex = 10 and ϕ = 0.032 in Fig. 7.

We here derive an expression for the local active stress tensor defined in a small subsystem denoted as K, with a volume V_K and a representative position of (x,y,z). Referring to a similar procedure given in the previous studies,^22,48 the active stress can be derived as

where α ∈ K indicates that the position R_α is within the region K. Essentially, identical expressions of the local active stress were previously derived (see the ref. 61 and 62, for example). More detailed derivation of eqn (34) is presented in the ESI.†

By using the conversion formula to transform the tensor from Cartesian to polar coordinates, we obtain the following expressions

Here, we assume axial symmetry about the x-axis of the local swimmers' properties in steady states. In our calculations, we evaluate the active stress in a small region with a linear dimension of a, which is not large enough to define an instantaneous stress tensor. However, this does not pose any practical problem when investigating the average properties of the active stress in steady states.

The rr-component 〈Σ^act_rr(x,y)〉 reflects the degree of the extent to which the swimmers face the probe sphere, either from the tail or the head. As shown in Fig. 7(a), around the probe, the magnitude of 〈Σ^act_rr(x,y)〉 is slightly less pronounced at the front side than at the rear, reflecting the front-rear asymmetry in P([small straight theta, Greek, circumflex]; x, y). This behavior intuitively suggests that the swimmers exert more force on the probe from the rear than from the front, which is more clearly seen in the inset plotting 〈Σ^act_rr(x,y) − Σ^act_rr(−x,y)〉. On the other hand, as shown in Fig. 7(b), 〈Σ^act_rθ(x,y)〉 exhibits a notable negative value around the regions indicated as (c)–(e) for x < 0 in Fig. 5. Note, however, that despite averaging over a period in which the probe particle moves several hundred times more than R, 〈Σ^act_rθ(x,y)〉 still exhibits significant fluctuations.

Then, we phenomenologically argue the contribution from the active stress to the reduction of the probe friction as follows. Integrating the net stress tensor over an arbitrary closed surface encapsulating the region including the probe particle, denoted as S_p, gives the friction force acting on that region, which is equal to −F_exx̂ balancing with the external force F_exx̂:

where  n̂ represents the unit vector that is normal to and pointing outward from the enclosing surface S_p, and dS denotes the surface area differential on S_p. The net stress comes from two different sources; one is the stress acting through the solvent, , and the other is the active stress . First, let us adopt S_p to the probe surface. In this case, on S_p because the swimmers do not exist due to steric repulsions. Therefore, the solvent solely accounts for the stress exerted on S_p. As shown in Fig. 2, the average probe speed 〈V_p,x〉 = F_ex/ζ can be increased compared with the case when moving in a pure solvent. This results in a smaller friction coefficient ζ than ζ₀, meaning that the average velocity gradient at the probe surface is decreased to the same extent as the increase in 〈V_p,x〉. We then ask: what causes such an increase in 〈V_p,x〉? To further argue this question, we adopt S_p to the surface slightly exterior to the probe surface, wrapping the peak of 〈ρ(x,y)〉. In this case, the contribution from the active stress is estimated aswhere cos θ =  x̂·n̂ and we assume axial symmetry around the x-axis and that the non-zero contribution from the active stress remains along the x-axis after taking the time average in a steady state. Eqn (38) represents the force exerted by the swimmers on the region surrounded by S_p and, hereafter, is referred to as F^act. As shown in Fig. 7, the magnitude of 〈Σ^act_rθ(x,y)〉 is approximately 0.1 of 〈Σ^act_rr(x,y)〉 on S_p, and F^act mostly depends on 〈Σ^act_rr(x,y)〉. Noticing that 〈Σ^act_rr(x,y)〉 − 〈Σ^act_rr(−x,y)〉 ∼ −10⁻³ for x < 0 on S_p, as indicated in the inset of Fig. 7(a), and the area of S_p is approximately 10⁴, F^act can be comparable with F_ex. This F^act is positive and thus can be considered a negative friction force, effectively acting as a driving force. On the other hand, the force from the solvent viscous stress acting on the same surface may be approximated as −ζ₀〈V_p,x〉, whose magnitude surpasses F_ex(= ζ〈V_p,x〉). The smaller friction coefficient ζ than ζ₀ indicates that the average flow velocity becomes larger than those observed when moving in a pure solvent environment. Essentially, velocity gradients are larger in the regions slightly outside the probe, while those are smaller in the thin region directly adjacent to the probe surface. The active stress serves to connect these varying behaviors across the different regions. With these two contributions of the active and solvent stresses, we may approximately express the force balance in a steady state as: −F_ex = −ζ〈V_p,x〉 ∼ −ζ₀〈V_p,x〉 + F^act, which is further rewritten as: (ζ₀ − ζ) ∼ F^act/〈V_p,x〉.

The above argument explains the notable reduction in ζ for smaller F_ex. However, as F_ex increases, while the front-rear asymmetry in the active stress is enhanced, F^act may not linearly increase in F_ex, resulting in less noticeable contributions to the reduction in ζ for larger F_ex.

4 Concluding remarks

Let us evaluate the applicability of our results to realistic situations by considering a typical experimental setup of dilute E. coli suspensions with a volume fraction similar to our simulations (∼0.01). Here, we assume that the swimming speed of E. coli is v_s ∼ 10 μm s⁻¹, the magnitude of the active force is F_A ∼ 1 pN, and the cell size is [small script l]₀ ∼ 1 μm. In a typical condition of magnetic-bead microrheology, the force exerted on the probe is in the pN range, and magnetic colloidal particles have diameters ranging from 1 to 10 μm (= 2R). In our simulation conditions, where the drag reduction is observed, the force ratio is F_ex/F_A ≲ 3 and the size ratio is 2R/[small script l]₀ ∼ 3, and these can be met in real experiments.²⁶ Therefore, we can expect that a substantial drag reduction occurs, if similar density and orientational distributions to those obtained in our simulations are realized in practical experiments.

In summary, using a pusher-type swimmer model, the present study investigated the microrheology of active suspensions through direct hydrodynamic simulations. We revealed that, with HIs, the friction coefficient of the probe particle can be significantly smaller than when immersed in a pure solvent for relatively small drag forces. The local swimmer states near the probe boundary predominantly influence the observed drag reduction. That is, in active suspensions, microrheology measurements are strongly influenced by local probing details and thus can be qualitatively different from macrorheology measurements. Hydrodynamic interactions induce the front-rear asymmetry in the swimmer orientation distributions. This asymmetry further causes a force/stress imbalance, resulting in an additional driving force acting on the probe particle through solvent-mediated interactions (HIs). However, without HIs, the local states of swimmers outside the (short) range of the direct interaction potential do not contribute to the friction force. Our simulation studies illuminate the fundamental importance of HIs in active suspensions and suggest potential microfluidics applications. These issues will be further explored in subsequent studies.

Author contributions

Takahiro Kanazawa: investigation, formal analysis, visualization, writing – review and editing. Akira Furukawa: conceptualisation, methodology, software, supervision, funding acquisition, investigation, formal analysis, visualization, project administration, writing – original draft preparation, review, and editing.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by KAKENHI (Grants No. 20H05619) from the Japan Society for the Promotion of Science (JSPS). We are also grateful for UROP (Undergraduate Research Opportunity Program) of the Institute of Industrial Science, the University of Tokyo.

Notes and references

1. S. Ramaswamy, Annu. Rev. Condens. Matter Phys., 2010, 1, 323–345.
2. M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao and R. A. Simha, Rev. Mod. Phys., 2013, 85, 1143.
3. C. Bechinger, R. Di Leonardo, H. Loẅen, C. Reichhardt, G. Volpe and G. Volpe, Rev. Mod. Phys., 2016, 88, 045006.
4. Y. Hatwalne, S. Ramaswamy, M. Rao and R. A. Simha, Phys. Rev. Lett., 2004, 92, 118101.
5. A. Sokolov and I. S. Aranson, Phys. Rev. Lett., 2009, 103, 148101.
6. J. Gachelin, G. Miño, H. Berthet, A. Lindner, A. Rousselet and E. Clément, Phys. Rev. Lett., 2013, 110, 268103.
7. H. M. López, J. Gachelin, C. Douarche, H. Auradou and E. Clément, Phys. Rev. Lett., 2015, 115, 028301.
8. A. Liu, K. Zhang and X. Cheng, Rheol. Acta, 2019, 58, 439–451.
9. V. A. Martinez, E. Clément, J. Arlt, C. Douarche, A. Dawson, J. Schwarz-Linek, A. K. Creppy, V. Škultéty, A. N. Morozov, H. Auradou and W. C. K. Poon, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 2326–2331.
10. S. Rafaï, L. Jibuti and P. Peyla, Phys. Rev. Lett., 2010, 104, 098102.
11. D. Saintillan, Exp. Mech., 2010, 50, 1275–1281.
12. M. C. Marchetti, Nature, 2015, 525, 37–39.
13. D. Saintillan, Annu. Rev. Fluid Mech., 2018, 50, 563–592.
14. T. Ishikawa and T. J. Pedray, J. Fluid Mech., 2007, 588, 399–435.
15. M. E. Cates, S. M. Fielding, D. Marenduzzo, E. Orlandini and J. M. Yeomans, Phys. Rev. Lett., 2008, 101, 068102.
16. B. M. Haines, A. Sokolov, I. S. Aranson, L. Berlyand and D. A. Karpeev, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 80, 041922.
17. L. Giomi, T. B. Liverpool and M. C. Marchetti, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81, 051908.
18. S. D. Ryan, B. M. Haines, L. Berlyand, F. Ziebert and I. S. Aranson, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 83, 050904(R).
19. M. Moradi and A. Najafi, EPL, 2015, 109, 24001.
20. T. M. Bechtel and A. S. Khair, Rheol. Acta, 2017, 56, 149–160.
21. S. C. Takatori and J. F. Brady, Phys. Rev. Lett., 2017, 118, 018003.
22. H. Hayano and A. Furukawa, Phys. Rev. Res., 2022, 4, 043091.
23. A. G. Bayram, F. Jan Schwarzendahl, H. Loẅen and L. Biancofiore, Soft Matter, 2023, 19, 4571–4578.
24. L. D. Landau and E. M. Lifshitz, Fluid Mechanics, Pergamon Press, New York, 1959.
25. T. M. Squires and T. G. Mason, Annu. Rev. Fluid Mech., 2010, 42, 413–438.
26. E. M. Furst and T. M. Squires, Magnetic bead microrheology, Microrheology, Oxford University Press, New York, 2017, ch. 8.
27. R. N. Zia, Annu. Rev. Fluid Mech., 2018, 50, 371–405.
28. X.-L. Wu and A. Libchaber, Phys. Rev. Lett., 2000, 84, 3017–3020.
29. D. T. N. Chen, A. W. Lau, L. A. Hough, M. F. Islam, M. Goulian, T. C. Lubensky and A. G. Yodh, Phys. Rev. Lett., 2007, 99, 148302.
30. K. C. Leptos, J. S. Guasto, J. P. Gollub, A. I. Pesci and R. E. Goldstein, Phys. Rev. Lett., 2009, 103, 198103.
31. J. J. Molina, Y. Nakayama and R. Yamamoto, Soft Matter, 2013, 9, 4923–4936.
32. G. Foffano, J. S. Lintuvuori, A. N. Morozov, K. Stratford, M. E. Cates and D. Marenduzzo, Eur. Phys. J. E: Soft Matter Biol. Phys., 2012, 35, 98.
33. G. Foffano, J. S. Lintuvuori, K. Stratford, M. E. Cates and D. Marenduzzo, Phys. Rev. Lett., 2012, 109, 028103.
34. E. W. Burkholder and J. F. Brady, J. Chem. Phys., 2019, 150, 184901.
35. E. W. Burkholder and J. F. Brady, Soft Matter, 2020, 16, 1034–1046.
36. Z. Peng and J. F. Brady, J. Rheol., 2002, 66, 955–972.
37. M. Knežević and H. Stark, New J. Phys., 2020, 22, 113025.
38. M. Knežević, L. E. A. Podgurski and H. Stark, Sci. Rep., 2021, 11, 22706.
39. A. Loisy, J. Eggers and T. B. Liverpool, Phys. Rev. Lett., 2018, 121, 018001.
40. C. Reichhardt and C. J. O. Reichhardt, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91, 032313.
41. J. P. Hernandez-Ortiz, C. G. Stoltz and M. D. Graham, Phys. Rev. Lett., 2005, 95, 204501.
42. P. T. Underhill, J. P. Hernandez-Ortiz and M. D. Graham, Phys. Rev. Lett., 2008, 100, 248101.
43. D. Saintillan and M. Shelley, Phys. Rev. Lett., 2007, 99, 058102.
44. V. Gyrya, I. S. Aranson, L. V. Berlyand and D. Karpeev, Bull. Math. Biol., 2010, 72, 148–183.
45. A. Decoene, S. Martin and B. Maury, Math. Modell. Nat. Phenom., 2011, 6, 98–129.
46. A. Furukawa, D. Marenduzzo and M. E. Cates, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 022303.
47. Y. Nakayama and R. Yamamoto, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 71, 036707.
48. J. J. Molina, K. Otomura, H. Shiba, H. Kobayashi, M. Sano and R. Yamamoto, J. Fluid Mech., 2016, 792, 590–619.
49. R. Yamamoto, J. J. Molina and Y. Nakayama, Soft Matter, 2021, 17, 4226–4253.
50. A. J. C. Ladd, Phys. Rev. Lett., 1993, 70, 1339–1342.
51. M. E. Cates, K. Stratford, R. Adhikari, P. Stainsell, J.-C. Desplat, I. Pagonabarraga and A. J. Wagner, J. Phys.: Condens. Matter, 2004, 16, S3903.
52. H. Tanaka and T. Araki, Phys. Rev. Lett., 2000, 85, 1338–1341.
53. A. Furukawa, M. Tateno and H. Tanaka, Soft Matter, 2018, 14, 3738–3747.
54. P. Español and P. Warren, Europhys. Lett., 1995, 30, 191–196.
55. S. Ji, R. Jiang, R. G. Winkler and G. Gompper, J. Chem. Phys., 2011, 135, 134116.
56. We speculate on the origin of significant non-monotonic behavior in ζ for the case without HIs, which is shown in Fig. 2(a). For small enough F_ex, the probe particle experiences collisions with the swimmers almost equally from both the front and the rear sides. However, as F_ex increases, the swimmers colliding from the front side are scattered more significantly than those from the rear. That is, more swimmers face the probe at the rear than at the front. Consequently, the force exerted from the front becomes gradually weaker as F_ex increases, which may cause the first thinning behavior exhibited in Fig. 2(a). Then, as F_ex increases more, very few swimmers catch up with the probe from the rear, making the impact of the swimmer colliding from the front more dominant. Notably, the average velocities of the probe and the swimmer are comparable at around the peak of ζ. As a result, the resistance increases with increasing F_ex, resulting in the force-thickening observed at intermediate values of F_ex. As F_ex is further increased, the active force becomes negligibly small compared to F_ex, leading to a decrease in ζ, which is observed as the second thinning. Even with HIs, a slight non-monotonicity in ζ is also found at smaller F_ex, to which a similar explanation may be applicable. A detailed investigation of these behaviors will be presented elsewhere.
57. A. S. Khair and J. F. Brady, J. Fluid Mech., 2006, 557, 73–117.
58. J. W. Swan and R. N. Zia, Phys. Fluids, 2013, 25, 083303.
59. G. Li and J. X. Tang, Phys. Rev. Lett., 2009, 103, 078101.
60. A. P. Berke, L. Turner, H. C. Berg and E. Lauga, Phys. Rev. Lett., 2008, 101, 038102.
61. R. A. Simha and S. Ramaswarmy, Phys. Rev. Lett., 2002, 89, 058101.
62. D. Saintillan and M. Shelley, Phys. Fluids, 2008, 20, 123304.

---

Footnote

† Electronic supplementary information (ESI) available: Explicit time-integration algorithm in the present simulations and supplementary simulation results. See DOI: https://doi.org/10.1039/d4sm00408f

---