Anomalous rheology of puller-type microswimmer suspensions

Unique rheological properties, absent in passive suspensions, have been observed in active suspensions Hatwalne ; Sokolov ; Gachelin ; Lopez ; Liu ; PNAS2020 ; Rafai ; Mussler ; Review2 ; Marchetti ; Review3 . For instance, in suspensions of E. coli, categorized as pusher-type microswimmers, viscosity decreases sharply below the solvent viscosity at relatively low concentrations and shear rates Sokolov ; Gachelin ; Lopez ; Liu ; PNAS2020 , frequently resulting in a superfluid state with zero viscosity Lopez ; PNAS2020 . In contrast, in suspensions of Chlamydomonas reinhardtii, classified as puller-type microswimmers, a steeper increase in viscosity compared to non-motile cells has been experimentally observed Rafai ; Mussler . Such anomalies in viscosity were first predicted by a seminal phenomenological model proposed by Hatwalne et al. Hatwalne , suggesting that stress can vary due to the swimmers’ self-propulsive forces when swimmers have orientational order under shear flow.

For rod-like particles, thermal or athermal rotational diffusion combined with shear flow produces weak orientational order Hinch_Leal1 ; Hinch_Leal2 , with which self-propulsive forces can provide non-zero active stress Review3 . However, the substantial origin of rotational diffusion and orientational order in active suspensions remains elusive. One potential candidate is hydrodynamic interactions (HIs) among constituent swimmers, which may significantly contribute to the rotational diffusion of rod-like swimmers Ryan ; Gyrya . Recently, our hydrodynamic simulations Hayano using model pusher-type swimmers designed to replicate E. coli suspensions reveal that, under shear flow, weak orientational order systematically evolves via hydrodynamic scattering during travel between boundary walls. The simulations also show that swimmers near the walls move in the forward-flow direction, which is an important characteristic of rod-like pusher-type swimmers Berke ; Li-Tang ; Kaya . Furthermore, in Ref. Hayano , we argue that even at around 1$\%$ concentrations of E. coli, HIs can dominate over other effects, such as thermal fluctuations, contributing significantly to the unusual behavior of active suspensions.

Meanwhile, as noted above, it has been experimentally observed that the viscosity of suspensions of motile Chlamydomonas is significantly larger than that of non-motile ones Rafai ; Mussler , which contrasts with the case of pusher suspensions. However, it is not trivial, even qualitatively, to determine what steady-state swimming properties are responsible for such anomalous rheology distinct from those observed for pusher-type swimmers. Simulation studies by Jibuti et al. Jibuti , in which a model swimmer with a spherical body and two beating spherical flagella replicating Chlamydomonas is used, show that due to the beating motions of flagella, the rotational motions of the swimmers are notably different from those of rod-like particles, leading to a significant increase in suspension viscosity. Although the simulation results and their physical implications are compelling, the potential influence of many-body HIs, which can greatly impact steady-state behaviors, should be further investigated. Additionally, in actual experimental systems, boundary walls exist, inevitably affecting the dynamics of swimmers and the flow field near walls. It is unclear how these effects influence overall rheological properties and how they differ from those of pusher-type swimmers.

In this study, we numerically address these problems through direct hydrodynamic simulations using model puller-type microswimmers. Our model swimmer is composed of a body and two flagellum parts. The body part is treated as a rigid body, while the flagellum parts are represented as two spherical massless particles simply following the body’s motions. This treatment keeps the relative position of these body and flagellum parts unchanged. For the $\alpha$-th swimmer ($\alpha=1,\cdots,N$), with $N$ being the total number of swimmers, we assume that a force $F_{A}{\hat{\mbox{\boldmath$n$}}}_{\alpha}$ acting on the body is exerted by the flagellum parts and that each flagellum part also exerts the force $-(1/2)F_{A}{\hat{\mbox{\boldmath$n$}}}_{\alpha}$ directly on the solvent fluid. Here, ${\hat{\mbox{\boldmath$n$}}}_{\alpha}$ represents the direction of the $\alpha$-th swimmer. The present particle-based model, schematically shown in Fig. 1(a), is essentially the same as those proposed in Refs. Graham1 ; Graham2 and used in Refs. Haines ; Ryan ; Gyrya ; Decoene ; Furukawa_Marenduzzo_Cates ; Hayano . Notably, in most of these former studies, it was supposed that the swimmers were of the pusher-type, while puller ones were employed in the current simulations. In this study, to address the effects of the swimmer’s shape on rheological properties, we employed three different aspect ratios $r=1.0,1.5$, and $2.0$ with the same volume ${\mathcal{V}}_{b}$ for the body part, as illustrated in Fig. 1(a). On the other hand, for these three models, the flagellum parts are composed of two identical spheres of radius $R^{\rm(f)}(=R^{\rm(b3)})$. Periodic boundary conditions are imposed in the $x$- and $y$-directions with the linear dimension $L$, and the planar 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 of $V/2$ and $-V/2$, respectively, whereby the mean shear rate is $\dot{\gamma}=V/H$. In our simulations, HIs among the constituent swimmers are incorporated by adopting the smoothed profile method (SPM) SPM ; SPM2 ; SPM3 , which is one of the mesoscopic simulation techniques LB1 ; LB2 ; FPD ; MPC ; DPD . More details of the model and simulations are provided in the Supplemental Information (SI).

Next, the simulation results are presented. First, in Fig. 2, we show the viscosity $\eta$ under various conditions. In this study, the viscosity is defined as

where $\Sigma_{xz}(x,y,\pm H/2)$ is the $xz$ component of the stress tensor at the walls and $\langle\cdots\rangle$ hereafter denotes the time average at a steady-state. Here, the solvent viscosity is scaled to be 1. We find that $\eta$ increases with increasing $\phi$ for all aspect ratios $r$. For elongated swimmers ($r=2.0,1.5$), the viscosity is higher than the solvent viscosity for all $\phi$, aligning with experimental results for Chlamydomonas suspensions Rafai ; Mussler . In contrast, for spherical swimmers ($r=1.0$), the viscosity is lower than the solvent viscosity at $\phi\lesssim 0.04$. This behavior strongly depends on $r$ and the volume fraction of the swimmers, defined as $\phi=N{\mathcal{V}}^{(b)}/L^{2}H$.

The viscosity can be divided into three parts: $\eta=\eta_{s}+\Delta\eta_{p}+\Delta\eta_{a}$, where $\eta_{s}(=1$ in this study) is the solvent viscosity, $\Delta\eta_{p}$ and $\Delta\eta_{a}$ are the passive and active contributions, respectively (see SI for more details). In the present framework, $\Delta\eta_{a}$ is given as

where ${\hat{\mbox{\boldmath$n$}}}_{\alpha}(t)$ is the unit vector representing the $\alpha$-th swimmer’s orientation at time $t$, and ${\hat{n}}_{\alpha,x}$ and ${\hat{n}}_{\alpha,z}$ are its $x$ and $z$ components, respectively. Essentially identical expressions of Eq. (2) were previously derived (see Refs. Haines ; Saintillan1 ; Review3 for example). From Eq. (2), when swimmers tend to align along the extension direction of the flow field ($\langle\hat{n}_{\alpha,x}\hat{n}_{\alpha,z}\rangle>0$), $\Delta\eta_{a}>0$. It is noteworthy that for pusher-type swimmers, a minus sign is applied in Eq. (2); an opposite change takes place in the puller case.

To further investigate the swimming states involved in the anomalous viscosity, the following steady-state quantities were examined: ${\rho}(z)=\sum_{\alpha=1}^{N}\langle\delta[{\mbox{\boldmath$r$}}-{\mbox{\boldmath$R$}_{\alpha}}(t)]\rangle$, ${\mbox{\boldmath$p$}}(z)=\sum_{\alpha=1}^{N}\langle{\hat{\mbox{\boldmath$n$}}}_{\alpha}(t)\delta[{\mbox{\boldmath$r$}}-{\mbox{\boldmath$R$}_{\alpha}}(t)]\rangle$, and ${\stackrel{{\scriptstyle\leftrightarrow}}{{\mbox{\boldmath$W$}}}}(z)=\sum_{\alpha=1}^{N}\langle[{\hat{\mbox{\boldmath$n$}}}_{\alpha}(t){\hat{\mbox{\boldmath$n$}}}_{\alpha}(t)-{\stackrel{{\scriptstyle\leftrightarrow}}{{\mbox{\boldmath$\delta$}}}}/3]\delta[{\mbox{\boldmath$r$}}-{\mbox{\boldmath$R$}_{\alpha}}(t)]\rangle$. Here, ${\mbox{\boldmath$R$}_{\alpha}}(t)$ is the position corresponding to the black circles drawn in Fig. 1(a) (see also the caption), $\rho(z)$ is the density, and ${\mbox{\boldmath$p$}}(z)/\rho(z)$ and ${\stackrel{{\scriptstyle\leftrightarrow}}{{\mbox{\boldmath$W$}}}}(z)/\rho(z)$ represent the polarization vector and the nematic order parameter tensor, respectively Review2 . These quantities are shown in Figs. 3(a)-(h) for various conditions. Here, the following points should be noted. The present puller model used here tends to produce large collective motions due to HIs among the swimmers, with strong polar order and fluctuations of swimmer density, in particular for the case of $r=1.0$. Therefore, although these quantities should ideally be symmetric with respect to $z=0$, achieving this in practical simulations requires an extremely large number of samples and long-time averages. In fact, the larger the number of samples is, the closer the simulation is to the ideal situation. These points are discussed in the SI.

In Figs. 3(a) and (b), $\rho(z)$ has significant peaks near the boundary walls, and otherwise, it is almost constant (except at higher $\phi$), indicating that the walls effectively attract swimmers. Such behaviors have been previously reported and discussed in the literature (for example, Refs. Berke ; Li-Tang ; Review1 ; LaugaB ; Figueroa-Morales ; Bianchi ; Denissenko ; Ezhilan ; Ezhilan-Saintillan ; Yan-Brady ). Note that for $r=1.0$, the symmetry of the $z$-axis is slightly disrupted, which is discussed in the SI.

Figures 3(c)-(f) display ${\mbox{\boldmath$p$}}(z)/{\rho}(z)$. $p_{z}(z)/\rho(z)$ exhibits similar behavior regardless of the aspect ratio $r$ and concentration $\phi$, indicating that swimmers temporarily trapped at the walls tend to tilt their ”heads” toward the walls. This behavior has been previously reported for rod-like pusher swimmers (see, for example, Refs. Vigeant ; Spagnolie_Lauga ; Sipos ). On the other hand, $p_{x}(z)/\rho(z)$ shows significant $r$-dependence: For swimmers with a rod-like body ($r=2.0$), $p_{x}(z)/\rho(z)$ is positive for $z>0$, and negative for $z<0$, suggesting that the swimmers move along the flow direction on average, which becomes most noticeable around the walls. In contrast, swimmers with a spherical body ($r=1.0$) move in directions opposite to the mean flow direction on average, which is also most significant around the walls. At $r=1.5$, a crossover between these two states takes place; a nearly flat $z$-dependence of $p_{x}(z)/\rho(z)$ is observed in the bulk region, while the swimmers around the walls still show a tendency similar to that at $r=1.0$.

In terms of viscosity enhancement or reduction, among Figs. 3(a)-(h), of particular interests are (g) and (h), exhibiting $W_{xz}(z)/\rho(z)$. For greater clarity, Eq. (2) is rewritten as

which shows that $W_{xz}(z)$ is directly related to the viscosity enhancement/reduction. Distinct behaviors are observed between elongated swimmers and spherical swimmers. For $r=2.0$, $W_{xz}(z)/\rho(z)$ is positive near the walls and slightly negative in the bulk region, and this trend becomes more pronounced with increasing $\phi$. On the contrary, for $r=1.0$, $W_{xz}(z)/\rho(z)$ is negative near the walls and positive in the bulk region. These results along with those in Fig. 2 suggest that the behaviors of $W_{xz}(z)$ around the walls predominantly contribute to the viscosity enhancement (reduction) in the elongated (spherical) swimmers. For $r=1.5$, an intermediate behavior between that at $r=2.0$ and $r=1.0$ is observed.

In Fig. 4(a), the $x$-component of the steady-state velocity field averaged over the $xy$-plane, $\langle v_{x}\rangle(z)$, is shown for $r=1.0$, $1.5$, and $2.0$ at $\phi=0.022$. Additionally, ${\rm d}\langle v_{x}\rangle(z)/{\rm d}z$, which is identical to the $z$-dependent solvent stress, is illustrated in Fig. 4(b). In the bulk region, the velocity gradient is enhanced and reduced at $r=2.0$ and $1.0$, respectively. On the other hand, near the walls, significantly negative velocity gradients and downward peaks in the solvent stress are observed for $r=2.0$, while the opposite occurs for $r=1.0$. Note that, in both the bulk and wall regions, at $r=1.5$, the observed change in ${\rm d}\langle v_{x}\rangle(z)/{\rm d}z$ is much smaller than those at $r=2.0$ and $1.0$. The viscosity observed at the walls results purely from solvent velocity gradients because the swimmers are not located in the narrow region near the walls due to the repulsive forces from the walls. Thus, this enhanced (reduced) viscosity immediately indicates the accelerated (decelerated) shear flow at the walls. In Figs. 5(a) and (b), we schematically display how swimming state affects solvent flow for the case of $r=2.0$. The following points are worth noting. In the present puller systems, significant orientational order is more apparent around the walls than in the bulk regions. Consequently, steady-state swimming properties near the walls predominantly determine the resultant viscosity [see Eq. (3)]. However, this tendency is completely different from the pusher case Hayano , where the swimming properties in the bulk region control rheology for similar situations. This difference can be attributed to the difference in swimming mechanisms, specifically to the induced velocity fields, between pusher and puller swimmers.

Finally, we discuss how the orientational and polar orders, shown in Fig. 3, are realized in the present model active suspension, and we attempt to deduce a possible description of global steady-states. To address this issue in more detail, we introduce the following conditionally averaged density and polarization vector:

where $H(x)$ is Heaviside’s step function. In Eqs. (4) and (5), $+$ and $-$ denote that the summation is limited to swimmers with $\hat{n}_{z}>0$ and $\hat{n}_{z}<0$, respectively.

Figures 6(a)-(d) show these quantities for $r=2.0$ and $1.0$ at $\phi=0.022$, $\dot{\gamma}=10^{-3}$ and $H=128$. As shown in Figs. 6(a) and (b), $\rho^{+}(z)$ has large and small peaks near the top and bottom walls, respectively, while for $\rho^{-}(z)$, the opposite is observed. In the bulk region, these values are nearly constant. This behavior indicates that most swimmers near the walls tilt their ”heads” toward the walls, while the swimmers tilting their ”heads” toward the bulk region quickly escape from the wall regions. For $r=2.0$, the differences between $\rho^{+}(z)$ and $\rho^{-}(z)$ near the walls are smaller than for $r=1.0$, suggesting that rod-like ($r=2.0$) swimmers are more easily peeled off from the walls than spherical ($r=1.0$) ones. In Fig. 6(c), $p_{x}^{+}(z)/\rho^{+}(z)$ and $p_{x}^{-}(z)/\rho^{-}(z)$ are plotted. At $r=2.0$, the average swimming directions are oriented along the mean flow direction both for swimmers moving towards and away from the walls. Contrary to this, at $r=1.0$, swimmers move against the mean flow direction. Additionally, Fig. 6(d) shows $p_{z}^{+}(z)/\rho^{+}(z)$ and $p_{z}^{-}(z)/\rho^{-}(z)$, demonstrating that the swimmers moving towards the walls tend to orient vertically more than those moving away from the walls. This tendency is slightly enhanced for $r=2.0$ compared with $r=1.0$.

As schematically shown in Fig. 7, rod-like swimmers ($r=2.0$), whether moving toward or away from the walls, tend to align with the flow direction. In contrast, swimmers with spherical bodies ($r=1.0$) tend to move against the flow direction and remain near the walls longer than those with $r=2.0$, resulting in a higher swimmer density around the walls. The key factor in understanding these behaviors is the competing effects of torques from the shear flow in the $xy$- and $xz$-planes. The torque in the $xz$-plane induces clockwise rotation in the swimmers. For swimmers moving opposite to the flow direction around the walls, further rotation is hindered by wall repulsion, acting as a counter torque to the shear flow. However, for swimmers moving along the flow direction, this steric hindrance is absent, allowing the $xz$-plane torque to orient them toward the bulk region, enabling escape from the walls. Meanwhile, the torque in the $xy$-plane aligns swimmers tilting toward the walls with the flow direction, a phenomenon known as the ”weathercock effect.” For rod-like swimmers, the influence of the $xy$-plane torque may be more significant than that due to the $xz$-plane torque, whereas for spherical swimmers, the latter one is weaker. This difference in the torque effect may explain why rod-like and spherical swimmers exhibit opposite tendencies in the moving directions. However, it is worth mentioning the following: the present puller swimmers tend to line up vertically due to HIs, frequently generating collective motions with polar orders (shown in SI). Such accompanying vertical motions seem to be geometrically more stable for shorter-body swimmers and persist for longer durations. These collective effects may significantly influence the determination of swimming direction more than the simple torque effects described above.

In this study, we investigated the steady-state swimming properties and associated viscosity changes for model puller-type microswimmers under various conditions. Although our models were designed to mimic Chlamydomonas by adjusting the shapes and relative positions of the swimmer body and flagellum parts, the simulation results show greater variation than the experimental observations, depending on the aspect ratio of the swimmer’s bodies. Swimming direction near walls depends significantly on the aspect ratio of the swimmer’s bodies. Rod-like swimmers tend to swim toward the flow direction, resulting in an increase in suspension viscosity. On the other hand, swimmers with spherical bodies swim against the flow direction, leading to a decrease in suspension viscosity. The present results contrast with our previous simulations of suspensions of pusher-type swimmers, where weak orientation order in bulk regions determines steady-state properties Hayano . Additionally, for the present puller systems, polar and orientational orders are more pronounced than those for pusher swimmers under similar conditions (compare Fig. 3 with Fig. 3 in Ref. Hayano ), which may occur because the puller-type model swimmers tend to line up vertically due to HIs, bringing their ”heads” and ”tails” closer. Such polar order has been reported for squirmers in previous studies Evans ; Alarcon . These findings strongly indicate that the rheological properties of active suspensions vary widely depending on the swimmers’ properties, such as their shapes and swimming mechanisms.

Finally, we address the discrepancy with experimental results. As noted above, the present model puller swimmers were designed to replicate actual Chlamydomonas, which is rather spherical. The experiments Rafai ; Mussler report an anomalous increase in viscosity, while our simulations show the opposite tendency for spherical puller swimmers. This discrepancy suggests the potential importance of the flagellum’s beating motions, which generate self-propulsive forces. However, in our simulations, active forces are represented by constant prescribed ones. Differences in several conditions between the experiments Rafai ; Mussler and our simulations may also explain this discrepancy. In the experiments, the applied shear rate $\dot{\gamma}$ is comparable to or larger than the intrinsic shear rate given by the swimmer’s velocity divided by the swimmer’s size. In the present simulations, on the other hand, $\dot{\gamma}(=10^{-3})$ is nearly ten times smaller than the intrinsic shear rate. Moreover, the volume fraction in our simulations is at most 0.056, while it is higher in the experiments. Other factors differing from experimental conditions, such as the system size, may also contribute. We will investigate these issues in future studies.