Orientational dynamics and rheology of active suspensions in weakly viscoelastic flows

Microswimmer suspensions in Newtonian fluids exhibit unusual macroscale properties, such as a superfluidic behavior, which can be harnessed to perform work at microscopic scales. Since most biological fluids are non-Newtonian, here we study the rheology of a microswimmer suspension in a weakly viscoelastic shear flow. At the individual level, we find that the viscoelastic stresses generated by activity substantially modify the Jeffery orbits well-known from Newtonian fluids. The orientational dynamics depends on the swimmer type; especially pushers can resist flow-induced rotation and align at an angle with the flow. To analyze its impact on bulk rheology, we study a dilute microswimmer suspension in the presence of random tumbling and rotational diffusion. Strikingly, swimmer activity and its elastic response in polymeric fluids alter the orientational distribution and substantially amplify the swimmer-induced viscosity. This suggests that pusher suspensions reach the superfluidic regime at lower volume fractions compared to a Newtonian fluid with identical viscosity.


INTRODUCTION
Systems of particulate matter suspended in fluids are prevalent in numerous natural and industrial processes.Active suspensions are particulate systems that are driven out of equilibrium by converting chemical energy or fuel into mechanical work to achieve selfpropulsion [1,2].Motile microorganisms are the prototypical example of active motion that is generated via metachronal actuation of hairlike cilia (Paramecium, Volvox ), by whipping cell-attached flagellar appendages (spermatozoa, algae), or by rotating a bundle of helical flagella (bacteria) [3].Microswimmers navigate through their environment by sensing or interacting with gradients in hydrodynamic, chemical, thermal, and light fields; a strategy known as taxis [4][5][6][7][8][9][10].A particular example is rheotaxis, where microswimmers experience hydrodynamic gradients and swim against the flow, which is relevant for biofilm formation and reproduction [2,8,11,12].
Pathogenic microswimmers when infiltrating human and animal bodies have to pass through mucus linings that lubricate and protect our respiratory tracks, eyes, urogenital and gastrointestinal systems [13].The presence of mucin fibres (typically 3-10 nm in length [14]) and DNA makes these mucus linings viscoelastic and shear thinning [15].The linings are often subjected to shearing motion, for instance, during blinking, coughing, reproduction, and continuos mucociliary clearance in respiratory systems.In such microbiological flows, the non-Newtonian fluid properties can alter the swimmer's rheotactic behaviour [16,17].
In this article we study the bulk rheology of microswimmers in viscoelastic flows.For Newtonian fluids, several rheological experiments on suspensions of exten-sile swimmers like E.coli have shown that activity can drastically reduce the effective viscosity [18][19][20], even down to the 'superfluidic' limit [21].The mechanism, first outlined by Hatwalne et al. [22] and elaborated further by Haines et al. [23] and Saintillan [24], is as follows.The orientations of elongated swimmers in shear flow follow periodic Jeffery orbits [25] and are also subject to thermal noise.This competition yields a mean orientation that points in the extensional quadrant of applied shear flow.With such an orientation extensile microswimmers (pushers) support the applied shear flow and thereby reduce the effective shear viscosity, whereas contractile swimmers like C. reinhardtii (pullers) resist it and thus increase viscosity [26].
Since almost all biological fluids are non-Newtonian, there has been a recent interest towards developing theoretical frameworks that capture the individual and collective dynamics of microswimmers in complex fluids [27][28][29].For example, experiments have shown reduced tumbling and increased persistence lengths of bacteria in polymeric fluids [30][31][32].Viscoelasticity can also initiate spatiotemporal order in active suspensions.For example, a recent study found that DNA polymers trigger oscillatory vortices in confined suspension of E.coli [33].However, determining the effective shear viscosity of an active suspension in a viscoelastic fluid has been an uncharted territory because of its complexity: the elastic relaxation of non-Newtonian fluids and their shearthinning/thickening property.
To gain better insights into biological and artificial motility in biological fluids, this communication investigates the role of non-linear polymeric stresses in non-Newtonian fluids.Specifically, this non-linearity allows us to directly couple a background shear flow to the active disturbance flow generated by a microswimmer.For spherical microswimmers in a Poiseuille flow, we already showed that due to such a coupling they experience a swimming lift force that depends on the swimmer type [17].Here, we will show that this nonlinear coupling influences the Jeffery orbit of an elongated swimmer in a shear flow, which thereby also fundamentally alters the bulk rheological response.Such activity-induced changes in the orientational dynamics cannot be observed in Newtonian fluids because the viscous stress tensor does not permit such a nonlinear coupling.De Corato and D'Avino [34] recently showed that a spherical microswimmer in the shear flow of a non-Newtonian fluid can exhibit a rotational velocity that is different from a Newtonian fluid.Although qualitative changes did not occur for weak viscoelasticity, strong viscoelasticity affected the orientational dynamics of pushers, pullers, and neutral swimmers distinctly.Here we consider an elongated particle and show that even weak viscoelasticity modifies the orientational dynamics, and ultimately, the rheology of a microswimmer suspension.
The current work uses the model of a second-order fluid with a single elastic relaxation time that captures the dynamics of polymeric fluids in the dilute limit (Boger fluids) [35,36].For weak elasticity, quantified by the Weissenberg number, we perform a perturbative analysis.Thereby we show that the inherent non-linearity in the polymeric stress tensor together with activity significantly alters the Jeffery orbits known from the Newtonian fluid and orbits of passive rods in a viscoelastic fluid [37,38].To evaluate the influence of thermal noise, we combine our results with the orientational Smoluchowski equation and find that the altered deterministic dynamics also substantially affects the orientational distribution, known as the suspension microstructure [39].It strongly differs between extensile (pushers) and contractile (pullers) microswimmers.The orientational distribution allows us to directly determine the effective viscosity of the active suspension from an orientational average over the swimmer stresslets.Our analysis shows that fluid elasticity reduces the effective viscosity of active suspensions for both pushers and pullers compared to Newtonian fluids.The reduction increases with activity and might allow to reach the superfluidic limit for smaller swimmer densities.

A. Setup and second-order fluid
We consider a dilute suspension of microswimmers in a shear flow of a viscoelastic fluid, for instance, consisting of polymers dissolved in a Newtonian fluid as shown in Fig. 1(a).Figure 1(b) depicts the coordinate system moving with the microswimmer that is modelled as an active prolate spheroid and swims with speed U s .The uniformly distributed polymers are much smaller than the microswimmers and hence modeled within a continuum description.The inertia-less hydrodynamics is governed by the mass and momentum conservation as ∇⋅V = 0 and ∇ ⋅ T = 0, respectively, where V is the velocity field and   T is the total stress tensor.It follows the second-order fluid (SOF) model: Here, E denotes the rate-of-strain tensor and S = 4E ⋅ E + 2δ ∆ E is the polymeric stress tensor, which is quadratic in E and contains the lower-convected time derivative.The SOF model not only allows us to capture the elastic effects pertinent to Boger fluids, but also to obtain analytical results for small Wi.Since we consider a steady shear rate in this work, we disregard the partial time derivative of E. In above governing equations, length, velocity, and pressure are already non-dimensionalized by the swimmer length (l = 2a), γl, and µ f γ, respectively, where µ f is the fluid shear viscosity.Also, the stress tensor is written in dimensionless units using characteristic numbers (Wi and δ).In particular, Wi = t relax γ is the shear based Weissenberg number that quantifies the importance of elasticity in the medium.It compares the shear rate γ or inverse shearing time to the polymer relaxation time t relax = (Ψ 1 + Ψ 2 )/µ f , where Ψ 1 and Ψ 2 are the normal stress coefficients.Furthermore, δ = −Ψ 1 /2(Ψ 1 + Ψ 2 ) is the viscometric parameter that typically varies between −0.7 to −0.5.In 'Methods: Hydrodynamic model' we explain in detail how we solve the governing equations using a systematic perturbation expansion in Wi in the limit of weak viscoelasticity.

B. Active spheroid in viscoelastic shear flow
To determine the dynamics of an active spheroid, we first note that a swimmer disturbs the flow field both passively (due to its rigid body) and actively (due to self-propulsion).The disturbance fields are implemented using hydrodynamic multipoles [40,41].Flagellated microswimmers like E.coli and Chlamydomonas generate a force-dipole flow field and higher order disturbances: source dipole, rotlet dipole, and force quadrupole.We included all four of them and found that only the forcedipole flow field affects the swimmer dynamics in leading order in Wi.The force-dipole field around the active spheroid is σ r Here, σ is a non-dimensional parameter equal to the ratio of the force-dipole strength (σ * ) to the stresslet imposed by the shear flow (8πµ f γl 3 ), where σ > 0 represents pushers and σ < 0 pullers.For the current work, we focus on swimmers of 5µm size and shear rates of order 0.5−5s −1 , which corresponds to σ of typical wild-type E.coli being roughly 0.04 − 0.1 [42,43].Here, the Weissenberg number can be tuned either by varying the shear rate or relaxation time of the fluid.The latter for the current system is t relax ≲ 0.1s, which for instance, can be realised in PEO solutions of molecular weight ranging from 2 − 4 ×10 6 g mol -1 and concentration between 0.25 − 0.5 wt% [44].We vary Wi from 0.05 − 0.5 by fixing t relax and changing the shear rates between γ ∼ 0.5 − 5s −1 , following a protocol similar to that for varying σ.
In Newtonian fluids, the equation of motion for the orientation p (θ, ϕ) gives the Jeffery orbits.To formulate this equation for viscoelastic shear, we note that the polymeric stress tensor S is quadratic in the rate-of-strain tensor and vorticity [35].Thus, similar to Einarsson et al. [45] we determine all terms that by symmetry contribute to the rate of change ṗ up to first order in Wi.Neglecting small terms, we arrive at in non-dimensional form, where Ω is the angular velocity and the superscript ∞ denotes the quantities that belong to the prescribed shear flow.The shape factor Λ = −1+λ 2 1+λ 2 contains the aspect ratio λ = a/b, the ratio of major to minor axis.The shape factor approaches +1 and −1 for needle and disk-like particles, respectively.Since microswimmers are usually elongated, we focus on prolate spheroids of Λ > 0.9, which corresponds to λ > 4; it also helps in simplifying the calculations (see 'Methods: Hydrodynamic model').The terms with coefficients α i and β i represent the active and passive viscoelastic contributions, respectively.These coefficients are evaluated explicitly as outlined in 'Methods: Orientation dynamics' using the Lorentz reciprocal theorem, where we also derive Eq. (1).For our relevant parameters, we find α 1 ≫ |α 2 |.This along with Eq. (1) indicates that the modification arising from the activity predominantly depends on the extensional part (E ∞ ) of the viscoelastic shear flow rather than its rotational component.Since the coefficients vary only weakly with λ and δ, they are treated as constants in the following discussion of results (see also Supplementary Note 2A).
In the absence of polymers (Wi = 0), Eq. ( 1) reduces to the well-known Jeffery equation [25] that has an infinite number of neutrally stable solutions, i.e., the microswimmer's orientation traces periodic orbits that depend on initial conditions.During orbital time period T = 2π(λ + λ −1 )/ γ, they spent the majority of their time (∝ λ) aligned near the flow-vorticity plane.When polymers are present, the passive viscoelastic effect breaks the degeneracy of Jeffery orbits and the microswimmer slowly drifts towards an alignment with the vorticity axis, known as the "log-rolling" state [38,46,47].Figure 2a shows the curve traced by the orientation of a spheroid on a unit sphere while drifting towards the vorticity axis.This dynamics can be deduced from Eq. (1) for passive spheroids (σ = 0) when the terms with coefficients β 1 and β 2 are present.In Supplementary Note 2C we also compare our results of passive spheroids with Brunn [38].Now, we illustrate the influence of activity.For weak pushers (0 < σ ≲ 0.1), the orbits towards log rolling look similar to the ones of passive spheroids, albeit with a larger aspect ratio (see Fig. 2b) since they exhibit more skewed trajectories.This can directly be inferred from Eq. ( 1), where activity (σ) in the term with Wi σα 1 modifies the shape factor.Since α 1 > 0, a pusher (σ > 0) effectively increases Λ, which corresponds to a higher aspect ratio.The opposite occurs for weak pullers (σ < 0) as they behave like a passive spheroid with reduced aspect ratio.Here the trajectories are more circular as shown in Fig. 2d.
As the dipole strength of a pusher increases beyond a critical value σ c , the orientation drifts to the shear plane (θ = π/2) and aligns at an angle ϕ eq with the flow direction (see Fig. 2c).We quantify this transition in Fig. 3 and show that σ c (dots close to 0) decreases with increasing Wi while ϕ eq increases with both σ and Wi.This dynamical behaviour can again be discerned from Eq. ( 1), which essentially balances the effect of rotational (Ω ∞ ) and elongational (E ∞ ) flow on p. Now, activity together with viscoelasticity allows to control this balance.In particular, for our case of α 1 ≫ |α 2 |, we expect the elongational flow to dominate the dynamics for increasing activity.Indeed, when the effective shape factor Λ + Wi σα 1 exceeds unity at a critical value σ c , the dynamics of p transitions from the orbital to the alignment state as favored by the elongational flow.In 'Methods: Dynamical analysis 1', we show for a modified dynamical system that the relevant eigenvalue of the dynamical matrix becomes real at σ c , as expected for such a transition, and the corresponding eigenvector gives the alignment angle ϕ eq .As σ increases, ϕ eq grows from zero and approaches π/4 for large σ, since in this limit, the elongational flow (E ∞ ) with its principal axis along ϕ = π/4 completely dominates the dynamics of p.
For pullers, Fig. 2e shows that increasing the activity induces a transition from log rolling towards rotation in the shear plane.This orbit is also observed for passive particles with oblate shape (Λ < 0) in viscoelastic flows [38,46,47].It appears here and acts as an attracting limit cycle, when the effective shape factor Λ + Wiσα 1 in Eq. (1) becomes negative due to the puller activity σ < 0. We have performed a stability analysis for the shear-plane rotation in 'Methods: Dynamical analysis 2' and calculated the Lyapunov stability exponent L. It determines the exponential time variation of the disturbed limit cycle.In the inset of Fig. 3, we show the Lyapunov exponent as a function of σ.For weak and moderate pullers, L > 0 shows that shear-plane rotation is an unstable limit cycle, where the system drifts towards stable log rolling.As the dipole strength of the puller becomes more negative, L turns negative meaning that shear-plane rotation is stable.Upon analysing the orbital dynamics in the shear-plane rotation, we find that increasing |σ| grad-ually slows down the swimmer's rotation near the y-axis, along which the flow gradient is applied.Eventually, a transition to permanent alignment with ϕ eq = π/2 occurs, which can be similarly analyzed as the alignment transition of pushers.However, here it occurs at the effective shape factor Λ + Wi σ c α 1 = −1.With further increasing |σ| the active spheroid tilts against the flow, in contrast to pushers, and asymptotes at 3π/4, which is the direction of the second principal axis of the elongational flow (E ∞ ).This behavior is illustrated in Fig. 2f and Fig. 3. Finally, Fig. 3 also shows that as Wi increases, the regime of shear-plane rotation shrinks.

C. Impact of noise on orientational dynamics
The deterministic orientational dynamics is disturbed by two types of stochastic reorientations of the swimming direction, which we now address with the help of the Smoluchowski equation.First, a bacterium tumbles, which is triggered when the rotation of one of its flagella reverses so that it leaves the flagellar bundle [48][49][50].For a wild type E.coli, tumbling occurs roughly every 1s, where it attains a new random orientation.Although tumbling is biased in the forward direction with a mean tumbling angle of roughly 68.5 ○ , we model it to be unbiased for computational ease because it only affects our results marginally as suggested by Nambiar et al. [51].Second, thermal rotational diffusion continuously reorients a microswimmer but its effect is small compared to tumbling.
We now evaluate the orientational distribution of an ensemble of non-interacting microswimmers in steady state, caused by the three mechanisms discussed so far: deterministic motion in background shear flow, tumbling, and rotational diffusion.The steady-state probability distribution function ψ(p) for the orientation vector of non-interacting microswimmers is governed by the Smoluchowski equation [52] Pe where ψ(p) is normalized to unity.The first term describes the orientational drift using ∇ p as the gradient operator on the unit sphere.To compare the strength of flow-induced reorientation to the mean time τ between two tumbling events, we introduce the flow Péclet number Pe f = γτ .Tumbling away from the swimming direction p is handled by the third term and rotational diffusion by the second term.For bacteria, the latter is typically small compared to tumbling since τ ∼ 1s [48] and D r ≲ 0.1s −1 , as an estimate from the Stokes-Einstein relation shows.In restricting ourselves to Eq. ( 2), we assume a spatially uniform system valid when hydrodynamic and steric interactions with bounding surfaces can be neglected [29,41].In particular, this means that the mean length of persistent swimming, U s τ , where U s is the swimming speed, is much smaller than the spatial extent of the system.Finally, we consider weak fluid elasticity with relaxation time t relax ≪ τ, D −1 r , so that the fluid relaxes faster than the time scales given by stochastic reorientations [53].Typical dilute and semidilute PEO solutions of molecular weight of ∼ 10 6 g mol -1 have t relax ≲ 0.1s [44], whereas the stochastic reorientation time scale is always 1s or larger.Thus, here we do not need to take into account any memory in the rotational noise.
In the following, we want to emulate a rheological experiment and vary the shear rate γ via Pe f , while the fluid and swimmer properties are kept constant.Therefore, in Eq. ( 1) for the rotational drift velocity ṗ, we rewrite Wi as Pe f De, where the Deborah number De = t relax /τ compares the fluid relaxation time to τ .Furthermore, the second relevant parameter, Wiσ, does not explicitly depend on γ but rather quantifies the strength of activity relative to fluid elasticity.To vary the activity independent of other parameters, we replace Wiσ by Pe a De/8π.Here, Pe a is the signed activity Peclet number that is positive for pushers and negative for pullers.Using U s ∼ σ * /(µ f l 2 ) for the swimming speed with the force dipole moment σ * [42], we can write it in the familiar form Pe a ∼ U s τ /l.Thus, its magnitude compares the persistence length to the body length l.A wild-type E.coli typically has Pe a ≲ 5 [48,54].We numerically solve Eq. ( 2) for arbitrary P e f by expanding ψ in spherical harmonics and taking into account the first 100 harmonics.The numerical solution is also verified analytically in the limit of Pe f ≪ 1 (see further details in 'Methods: Kinetic model').
Figure 4 shows the orientational probability distribution of passive and active particles for various cases.It quantifies the orientational 'microstructure' of a suspension of non-interacting and orientable particles subject to shear flow [51,55,56].We begin with discussing the case of weak shear rate (Pe f < 1) in Figs.4a-d.In the absence of shear flow, the microstructure is solely governed by rotational noise and is therefore isotropic.For weak shear rates the extensional part (E ∞ ) of the shear flow primarily distorts the microstructure, which peaks near the extensional axis as noted for Newtonian fluids by Hinch and Leal [57].We quantify this in 'Methods: Kinetic model 1' by solving Eq. ( 2) via a perturbation expansion in Pe f and obtain the first-order correction: In the absence of activity, there is no deviation from the Newtonian microstructure at this order.In Figs.4a-d we use a weak shear rate of Pe f = 0.5 and |Pe a | = 20, which is an activity close to mutated strains of tumbling E. coli (pusher) [58] and the algae C. reinhardtii (puller).The observed distributions follow the trend of Eq. ( 3).While pushers (Pe a > 0) enhance the alignment along the extensional axis of the applied shear flow, pullers (Pe a < 0) weaken it.We note that the peak of the orientational distributions in the shear plane (plot d) exhibits a small deviation from ϕ = π/4.This is due to the rotational flow (Ω ∞ ), which contributes to ψ(p) in second order in Pe f [56].
As the shear rate increases, the deterministic dynamics becomes more visible.We illustrate this in Figs.4e-h for Pe f = 5.Compared to Fig. 4a, the peak in ψ(p) for passive particles (Fig. 4e) shifts towards the flow axis and becomes more anisotropic along the flow-vorticity plane, as the deterministic velocity ṗ is smallest in this plane.Eventually, ṗ approaches zero in the deterministic log-rolling state observed in viscoelastic fluids (Fig. 2a).However, the deterministic log-rolling state is not completely reflected in the distribution because the slow relaxation towards the vorticity axis is always interrupted by the dominating stochastic reorientations (see Supplementary note 4, where we illustrate the orientation distribution for a larger weakly Brownian particle that resembles log-rolling).Note that due to the fluid elasticity, the peak of ψ(p) in Fig. 4h is slightly closer to the flow axis compared to the Newtonian case.Furthermore, the inset therein also illustrates a broader distribution in the flow-vorticity plane (ϕ = 0).Further increasing Pe f spreads the distribution even more in the flow-vorticity plane and, ultimately, for Pe f ≳ 100 the peak develops along the vorticity axis, which is entirely reminiscent of log rolling.However, this is a regime which we cannot strictly capture as our analysis requires De, Wi < 1, which means Pe f < 10.These findings are qualitatively consistent with earlier studies on passive fibre suspensions at moderate shear rates [53,59] (further illustrated in Supplementary Note 4).
Figure 4f demonstrates that the activity of pushers in conjunction with the higher shear rate (Pe f = 5) makes the distribution more anisotropic.According to Fig. 4h, pushers strongly focus at an angle ϕ = 18 ○ in comparison to a weaker peak of passive particles near to flow axis.This angle is close to the alignment angle ϕ eq = 21 ○ of the deterministic state with the parameters Wi = 0.5, σ = 0.16 in Fig. 3. Conversely to pushers, pullers reduce the anisotropy in the orientational distribution as shown in Fig. 4g.We can directly relate this observation to the deterministic velocity ṗ with parameters Wi = 0.5, σ = −0.16 that indicate again the log-rolling state; the resulting orbit is similar to that depicted in Fig. 2d.Compared to the passive case, the orbit is more circular, which we already related to a reduced effective aspect ratio when discussing the deterministic orientational dynamics.Therefore, the orientation is less constrained to the flow-vorticity plane.This is explained by the dynamics of the azimuthal angle ϕ in the inset of Fig. 4g that shows lesser time spent in the plane compared to the inset in Fig. 4e.As a result, the distribution is broader and less anisotropic [60].Consequently, the particle is less susceptible to the effect of stochastic reorientations and thus, its distribution has a peak closer to the flow axis.Hence, unlike active Newtonian suspensions [24], the microstructure of microswimmers in a viscoelastic fluid is clearly sensitive to their hydrodynamic signature being either pushers or pullers.
We also examined the orientational distributions of microswimmers that do not tumble and only experience rotational thermal noise.Therefore, in this case we eliminate the last term on the right-hand side of Eq. ( 2) and redefine our parameters Pe f , Pe a , and De replacing τ by D −1 r .For passive particles and activity Pe a ∼ O(10), the microstructure is qualitatively similar to the distributions in Figs.4a-h.However, microswimmers that do not tumble can have larger persistence lengths [61].In particular, one can genetically modify E.coli such that tumbling does not occur [42,58,62].Therefore, in Figs.4i and j we show the orientational distributions for Pe a = ±100.Compared to Fig. 4f, increasing the activity of a pusher moves the peak in the distribution function (at ϕ = 34 ○ ) even closer to the alignment angle (ϕ eq = 36 ○ ) of the deterministic alignment state.For pullers we observe an even more drastic change upon increasing the activity.The peak in the orientational distribution shifts to the quadrant where the compressional part of the shear flow occurs (see Fig. 4j).Interestingly, this is not due to the deterministic alignment of pullers as the parameters (Wi = 0.5, σ = −0.8)belong to the shear-plane rotation state (see Fig. 3).As already explained in the discussion of Fig. 3, the rotational velocity ṗ slows down near ϕ = π/2 when the effective shape factor Λ + Wi σ c α 1 becomes negative.Therefore, the active particle resembles a passive oblate spheroid [63].Indeed, the inset of Fig. 4j shows how the puller orientation spends more time near ϕ = π/2.This again shows that even in the presence of significant noise the microstructure is determined by the deterministic dynamics.

D. Shear viscosity of active suspensions
Finally, we evaluate the effective shear viscosity of a dilute active suspension from the total stress tensor, Σ = Σ f + Σ p , where subscripts f and p refer to fluid and particle contributions, respectively.The deviatoric component of Σ provides the shear viscosity, µ = µ f + φµ p , where φµ p = Σ p xy / γ is due to the suspended microswimmers and φ = na 3 is proportional to their volume fraction with n being the particle density.Note that the shear viscosity µ f of the second-order fluid does not depend on the shear rate, while µ p is dependent on it.To evaluate µ p , we need to average over the stresslets Π(p) generated by the suspended particles.Using the orientational distribution, we obtain where Π(p) is the sum of three contributions [24], which give the following stress tensors: ] due to thermal reorientations [64], and Σ F p due to the resistance of passive particles to shear.For Newtonian fluids, the latter was first derived by Einstein [65] for spherical particles, and then later generalized to elongated particles by Hinch and Leal [39,55].For the second-order fluid, we follow Férec et al. [66] and approximate it using the Geisekus form, Σ F p = −µ f nl 3 A γ ∇ ⟨pp⟩/2, where A = π/6 ln(2λ) is a shape factor to access the friction of a long slender body and ∇ ⟨pp⟩ is the upper-convected time derivative of the second moment of ψ(p).The expression for Π F was originally derived for dumbbells suspended in a Newtonian fluid [67].As in Ref. [66], we use it here as an approximation for the stress response of particles in a weakly viscoelastic fluid, as there are no closed form expressions available.The consequences of the secondorder fluid come in through the orientational dynamics ṗ in Eq. ( 1) and further terms in the upper-convected derivative (see 'Methods: Rheology').The expressions for the three contributions to the particle stress tensor are pertinent to studies on slender rods and fibres (i.e., Λ → 1), and henceforth, we will be focusing on spheroids with larger aspect ratios (λ = 10).Using Eq. ( 4) in the definition of the particle-induced contribution to shear viscosity, φµ p = Σ p xy / γ, together with our characteristic parameters, we obtain (5) We use numerical solutions of Eq. ( 2) for the orientational distribution ψ(p) to evaluate µ p for varying Pe f and Pe a .We also match the numerical values of µ p with an expression, which we obtain using the analytic form of ψ(p) from Eq. (3) in the limit of small Pe f (the derivation is provided in 'Methods: Rheology').In Supplementary Figure 4, we also show that our results agree with Saintillan [24] in the Newtonian limit.
Figure S4a shows the particle-induced contribution to shear viscosity normalized by the bare fluid viscosity µ f , which is plotted versus shear strength Pe f .We begin with discussing suspensions in a Newtonian fluid (dashed lines) [24].When a suspension of passive rods (black) is sheared weakly (Pe f ≪ 1), the orientational microstructure is governed by Eq. ( 3) with De = 0 and resembles Fig. 4a.Here, the rods are aligned along the extensional axis of the applied shear flow and resist shearing (µ p > 0), which enhances viscosity.For a suspension of pushers (brown dashed) in the weak-shear regime, the microstructure still resembles Fig. 4a, but now the extensile force dipoles support the elongational part of the shear flow, which results in µ p < 0 so that the total shear viscosity is smaller than µ f .Conversely, pullers (green dashed) with their contractile force dipoles additionally resist shear flow and thereby enhance viscosity.Note that, since the microstructure of active suspensions in a Newtonian fluid is unchanged by the hydrodynamic signature of microswimmers, the magnitude of the activity contribution to µ p is identical for pushers and pullers.As the shear rate increases, the magnitude of µ p for both          passive and active rods reduces, which resembles a typical shear-thinning/thickening behaviour.For passive particles the microstructure is similar to Fig. 4e.The alignment close to the flow axis explains why the resistance of particles to the applied shear flow is reduced.For active particles, the activity-induced flow becomes less and less important with increasing shear rate as quantified by the inverse proportionality to Pe f in Eq. ( 5).Now, we discuss the results for the second-order fluid (SOF) in Fig. S4a (solid lines).For the suspension of passive rods (black), we do not find an appreciable deviation from the Newtonian case until Pe f ≈ 10, similar to Refs.[53,59].However, for active suspensions the modifications are substantial and qualitatively different for pushers (brown) and pullers (green).For pushers, the reduction in viscosity is further enhanced compared to Newtonian fluids, whereas the viscosity enhancement of pullers is reduced.For a fixed viscosity of µ f and a volume fraction of φ = 0.1, a pusher suspension of Pe a = predicts a 10% viscosity reduction, whereas in a Newtonian fluid the reduction is 6%.We understand this viscosity reduction by using the orientational distributions in Fig. 4: pushers are even more aligned in the extensional quadrant of the applied shear flow (Fig. 4b,d,f), which explains their further increased support of shear flow.In contrast, pullers are less aligned compared to the Newtonian case (Fig. 4c,d,g) and, therefore, oppose the shear flow to a reduced extent.Figure S4b for pushers clearly shows that the active stress (Σ A ) dominates and hence is primarily responsible for the observed viscous response.
In the limit of shear rate much smaller than tumbling rate, Pe f ≪ 1 , we can calculate µ p /µ f analytically using Eq. ( 3) (see 'Methods: Rheology') and obtain: We elaborate the result here since the influence of passive and active rods on the shear viscosity is the strongest in this limit.In Fig. S4c we plot µ p /µ f vs. Pe a for tumbling microswimmers (solid lines).As Eq. ( 6) shows, elasticity in the fluid adds a quadratic dependence in Pe a , while in a Newtonian fluid the dependence is only linear [24].For all Pe a , except a small region close to Pe a = 0, the values of µ p /µ f remain below those obtained in the Newtonian limit.Thus, fluid elasticity reduces the shear viscosity for both pusher and puller suspensions.
Note that, as we increase the activity of pullers (Pe a ≲ −50), they too can reduce the shear viscosity (µ p /µ f < 0) like pushers.This is a consequence of the orientational distribution shown in Fig. 4j; pullers with large activity align preferentially in the compressional quadrant of the shear flow and thereby also support the shearing fluid similar to pushers aligned in the extensional quadrant.
As described earlier, we can treat the case of nontumbling microswimmers by replacing τ by D −1 r in the definitions of Pe f , Pe a , and De.With these parameters, Eq. ( 6) can be formulated in the limit τ → ∞ to show that the viscosity contribution of non-tumbling active rods also follows a quadratic dependence in Pe a (Fig. S4c, dashed lines).Non-tumbling microswimmers can exhibit high persistence lengths.Thus, in Fig. S4d we show for |Pe a | = 100 that their contribution to viscosity is negative over a wide range of shear rates for both pushers and pullers in a second-order fluid.Hence, we find that the elasticity of a fluid always results in a reduced total viscosity as compared to suspensions in a Newtonian fluid of identical base viscosity (µ f ).

DISCUSSION
In this article we show how activity influences the dynamics of a sheared microswimmer suspension in a viscoelastic fluid at an individual level and in the bulk.At the individual level, the orientational dynamics of passive rods is well known from Jeffery's orbits in Newtonian shear [25] and 'log-rolling' orbits in viscoelastic shear flow [37,38].Our analytical result [Eq.( 1)], derived for elongated active spheroids in a second-order fluid, demonstrates how the active flow field of a swimmer modifies the orientational dynamics.Extensile swimmers like E.coli (pushers) drift to the shear plane and align at an angle to the flow direction.For contractile swimmers like C. reinhardtii (pullers), activity effectively lowers their aspect ratio.With increasing dipole strength, pullers show log rolling, transition to shear-plane rotation, and ultimately align at an angle against the flow direction.The latter occurs for strong pullers at dipole strengths relevant to artificial swimmers with large propulsion speed [68,69].
To demonstrate the impact of the individual swimmer dynamics on the bulk rheological behaviour, we employ the Smoluchowski equation to evaluate the orientational probability distribution of an ensemble of active spheroids called the suspension microstructure [39].Accounting for tumbling and rotary diffusion, we find that the microstructure is sensitive to the hydrodynamic signature of the swimmer unlike suspensions in Newtonian fluids [22,24,51].Pushers align more strongly in the extensional quadrant of the applied shear flow, while the alignment of pullers is significantly weaker.This activityspecific microstructure significantly modifies the effective shear viscosity compared to Newtonian fluids [21].In particular, the viscosity reduction of a dilute pusher suspension is more pronounced under viscoelastic shear flow, while the viscosity enhancement of pullers is weaker.Thus, the activity of a microswimmer contributes in two ways to the rheology of swimmer suspensions in a viscoelastic fluid; directly through its active stresses and indirectly by coupling to the elasticity of the fluid and thereby influencing the orientational microstructure.In particular, in the weak-shear limit the particle-induced contribution to viscosity scales quadratically with the activity Pe a .In total, we find that the elasticity of a second-order fluid always reduces the total viscosity of the microswimmer suspension, as compared to a Newtonian fluid of identical base viscosity.Especially for pushers, this might help to reach the regime of superfluidity at lower volume fractions compared to Newtonian fluids.We note that superfluidity is reported to be also associated with onset of collective motion [70,71], and requires further analysis at higher microswimmer densities.
Most biological fluids are viscoelastic.We presented a first systematic study of the individual dynamics and the bulk rheology of microswimmers suspended in such fluids.Earlier investigations on active suspensions in quiescent [72] and vortical viscoelastic flows [73] assumed conventional Newtonian Jeffery orbits.In the light of our results, these studies on the collective behaviour might need to be revisited.Our investigations makes several assumptions to address the complexity of microbial flows, which offers ample opportunities for future developments.Biological fluids can be more complex and exhibit shear-thinning/-thickening properties or several characteristic relaxation times, which requires more evolved modeling using, for example, the FENE-P or Giesekus model.Furthermore, mucus besides being significantly shear-thinning has relaxation times larger than 1 second, which is of the order of the bacterial mean free tumble time [15,28,29].Thus, noise becomes non-Markovian and memory needs to be incorporated in the stochastic description, for example, within a generalized Langevin equation [74].

A. Hydrodynamic Model
The hydrodynamics around the spheroidal particle is governed by the mass and momentum conservation as ∇ ⋅ V = 0, ∇ ⋅ T = 0.Here T is the total stress tensor defined as T = −P I+2 E+Wi S [35], where S = 4E⋅E+2δ ∆ E. The lower-convected derivative is defined as Since we consider a steady shear flow and a steadily active swimmer, we can disregard the temporal derivative in the above stress tensor.In particle frame of reference, the boundary condition on its surface is the no-slip condition V = Ω p × r, where the rotational velocity Ω p is currently unknown and will be evaluated by using the torque-free condition.We split the velocity field into the disturbance field v and background flow field The equations governing the disturbance flow field are Here e is the rate of strain tensor for the disturbance flow (∇v + ∇v † )/2, whereas w, ∆ e and ∆ w are the different constituents of the polymeric tensor (s) due to the disturbance flow field.Specifically, w is the interaction tensor defined as E ∞ ⋅ e + e ⋅ E ∞ , arising from the interaction between background flow and disturbance field; and ∆ w is the lower-convected derivative of e and E ∞ with respect to V ∞ and v, respectively: The boundary conditions of the disturbance field at the particle surface and far away are respectively.Here ∆Ω is the difference in particle angular velocity and background vorticity (Ω p − Ω ∞ ).
For small values of Wi, the disturbance field variables are expanded as f = f (0) +Wi f (1) +⋯.Here, f is a generic field variable that represents velocity (v), pressure (p), and angular velocity (Ω p ).We substitute this expansion in the Eq. ( 8) and obtain the O(1) Stokes problem as with the boundary condition at particle surface: and a decaying condition at infinity (v (0) → 0).Using the finite multipole expansion around the spheroid [40,45,75], we obtain the disturbance velocity for a passive spheroid at O(1) as Here Q represent the spheroidal multipoles i.e. integral representation of fundamental singularities spread on a line extending from one foci to another (−c to c, where c = √ λ 2 − 1/λ); Q R , Q S represent rotlet and stresslet around the spheroid, respectively.Here χ = 1 8(λ 2 −1) represents the strength of higher order quadrupolar field (Q Q ij,llk ).Since we focus on elongated particles, we exclude this component for computational ease (λ = {5, 10} correspond to χ = {0.005,0.001}).Following Einarsson et al. [45], we simplify these multipoles as Here I and J represent various integrals defined as In Eq. ( 13), the fourth order orientation tensors (p A , p B , p C ) and other coefficients (A, B, C) are described in the Supplementary Note 1.
At O(1), we add the activity using the far-field descriptions, which consists of a force-dipole (FD), source-dipole (SD), rotlet-dipole (RD), and force-quadrupole (FQ) [41]: At O(Wi), the governing equation is ), where ), (17) The boundary condition at the particle surface being with a decaying condition at infinity (v (1) → 0).Conventionally, a solution to Eqs. ( 17) and (18) (i.e., v (1) ) is sought, which, on the implementation of torque-free condition, reveals the modification of Jeffery orbits (Ω p ).We employ the reciprocal theorem [75,76] to avoid solving for the first order flow field and directly obtain the rotation velocity as Here, superscript t refers to the test problem where the particle rotates at a unit velocity Ω t = e i , where i corresponds to either of the three Cartesian coordinates.The test torque is: We solve Eq. ( 19) for three test field torques (along the three Cartesian coordinates) and obtain the following system of equations: R † ⋅ Ω p = {I 1 , I 2 , I 3 }, where † represents transpose and I i is the solution of volume integral (19) for the test field in i th unit vector.We further simplify the above relation by noting that R is a symmetric matrix and ṗ(1) = Ω (1) p × p: Eq. ( 21) can be evaluated over a discretized angular grid, where each point requires solving the volume integral in Eq. (19).The polymeric stress therein (s (0) ) is given by Eq. ( 8) and is entirely dependent on O(1) disturbance flow fields i.e., passive disturbance in Eq. ( 13) and active disturbance in Eq. ( 16).We find that out of all active fields, only force-dipole (∼ r −2 ) provides a modification at the current order of approximation; the rest decays in odd powers of distance (∼ r −3 ), and due to antisymmetry, give zero contribution to the volume integral at O(W i).

B. Orientation dynamics
We evaluate ṗ(1) analytically by noting that it stems from the leading order polymeric stress [specifically s (0)  in Eq. ( 19)]; its form in Eq. ( 8) suggests that the modification for a passive particle will be quadratic in the flow gradient tensor, which can be decomposed in symmetric E ∞ and antisymmetric O ∞ (rotation-rate tensor) components.Following Einarsson et al. [45], this can be written in the general form as: where superscript P denotes passive contribution.The coefficients of the fifth-order tensor K are composed of all possible permutations of the orientation vector p with E ∞ and O ∞ : (β [1]  n p n1 p n2 p n3 p n4 p n5 + β [2]  n p n1 p n2 p n3 δ n4n5 +β [3]  n p n1 δ n2n3 δ n4n5 ) .
Here β represents the unique coefficients for all the 5! terms.
When activity is added in the hydrodynamics, the form of polymeric stress tensor in Eq. ( 8) suggests that the gradients of disturbance velocity (directed with p) will multiply with gradients of passive disturbance (which originate from E ∞ and O ∞ ) to yield further modification to ṗ(1) .Thus, the contribution resulting from interaction of active disturbances (p) with background flow (E ∞ and O ∞ ) takes the following form at O(Wi): where coefficients of the third-order tensor L are composed of the following combinations (α [1]  n p n1 p n2 p n3 + α [2]  n p n1 δ n2n3 ) .
Here α represent the unique coefficients for all the 3! terms.We simplify Eq. ( 22) and ( 23) by using the symmetry arguments and noting that magnitude of p is always unity.We obtain the following six irreducible terms: We calculate the coefficients (α i and β i ) by evaluating ṗ(1) numerically for six independent orientations (p) and solve the system of equations from Eq. ( 24) to extract the coefficients.We find that β 3 and β 4 are nearly zero (≲ 10 −5 ) and thus neglect them.We obtain the analytical approximation as We use the above equation in the main text as Eq. ( 1), where we write O ∞ ⋅p as Ω ∞ ×p.In Supplementary Note 2A, we show the weak variation of these coefficients with particle aspect ratio λ and the viscometric coefficient δ.
C. Dynamical analysis

Eigenvalue analysis of alignment dynamics
Following Bretherton [77], we explain the alignment of a particle in the shear plane by extracting a fundamental matrix solution of Eq. (1).For this, we first consider the non-dimensional Jeffery equation in the expanded form as which can be represented in an unconserved form that facilitates a fundamental matrix solution [78]: Here q(t) |q(t)| = p(t) and |q(t represents the exponential elongation of q.We study the angular dynamics of q, as its normalization yields back the orientation vector p.The solution to Eq. ( 27) is The eigensystem of the above exponential matrix determines the orbital dynamics of spheroid.For the twodimensional shear flow, the eigenvalues of the matrix in the exponent (ΛE ∞ + O ∞ ) are: 0, ± √ Λ 2 − 1/2.The nonzero eigenvalues are purely imaginary as 0 < Λ < 1, where Λ = 1 for an infinitely slender particle.As a consequence of this imaginary pair, we observe degenerate infinite solutions of orbits in Newtonian fluids (i.e., Jeffery orbits).For the case of pure elongational flow (O ∞ ≡ 0), the eigenvalues are real: 0, ±Λ/2, which reveals the absence of periodic orbits.The normalized eigenvector corresponding to the positive eigenvalue determines the equilibrium orientation: {1/ √ 2, 1/ √ 2, 0} i.e., 45°or 225°in the two extension quadrants.
We now use Eq. ( 1) for the second-order fluid and consider only the active viscoelastic component because the passive effects do not contribute to the alignment [this assumption is later verified by matching the results with numerical integration of complete Eq. ( 1)].In this case, the equation of motion for the unconserved vector q is which yields the solution q(t) We obtain the eigenvalues of the fundamental matrix as 0, ± 1 2 These eigenvalues turn real when σ > For instance, for a particle of aspect ratio λ = 5, the negative and positive critical values are σ crit = −0.448/Wi;0.016/Wi, which matches with Fig. 3 generated from the integration of full Eq.( 1) (see Supplementary Note 2B).In the alignment regimes of Fig. 3, we further find that the normalized eigenvector corresponding to the positive eigenvalue matches the angle of alignment (ϕ eq ) obtained via numerical integration.In addition to validating our results of Fig. 3, Eq. ( 30) provides a key insight that the contribution from active disturbances is to effectively alter the strength of elongation (prefactor of E ∞ ) and the rotation rate (prefactor of O ∞ ).The first prefactor is the more decisive as α 1 ≫ |α 2 |.It suggests that a pusher (σ > 0) disturbs the local shear flow in order to increase the weight of elongation.Once the activity exceeds the critical limit, the local shear flow transforms into an effective elongation flow whose axis of elongation points in the direction of the normalized eigenvector associated with the positive eigenvalue.As pusher's activity further increases, the locally elongated flow asymptotically approaches the pure elongation state where the impact of O ∞ is negligible in comparison (similar to our aforementioned case of O ∞ ≡ 0 in a previous paragraph).

Stability analysis of orbits
Here we find the stability exponent of the complete non-linear Eq. ( 1) to analyze the onset of shear-plane rotation state as shown in Fig. 2,3.For Newtonian fluids, θ = π/2 is one of the infinite neutral Jeffery orbits.For SOF, the polymeric stress in the fluid lifts this degeneracy.To quantify this, we write down the angular dynamics as and perform a Taylor expansion near the shear-plane: θ(t) = π/2+ϵ(t).Here ϵ represents the deviation from Jeffery orbit and and we determine its growth i.e., whether it grows or shrinks with time.At O(ϵ) we obtain which can be simplified and integrated over an orbit to obtain Here the integration is over −2π because the particle's rotation due to shear is in −ϕ direction.Eq. ( 32) can be expressed in the form of Lyapunov exponent (ϵ = ϵ 0 exp[LT ]) [79] as where T is the time period of a Jeffery orbit 2π(λ+λ −1 )/ γ.We show the solution to Eq. ( 33) in the inset of Fig. 3.
First, we detail the evaluation of microstructure near equilibrium (Pe f = γτ ≪ 1) i.e., when stochastic reorientation dominates the shear-induced reorientation.In this limit, we expand the orientation distribution as We substitute this in Eq. ( 2) and collect the zeroth and first order terms in Pe f .At O(1), we get ψ (0) = 1/4π as the isotropic microstructure.At O(Pe f ), we have Here we denote the equation of motion (Eq. 1) as having three parts: ṗ = ṗ(0) +DePe f ṗ(1) P +DePe a ṗ(1) A , where subscripts represent passive and active components.Upon substituting this in Eq. ( 35), we find that the ṗ(1) P term only contributes at O(Pe 2 f ), and is thus neglected.Upon simplification, we find that the microstructure at O(Pe f ) is governed by Following [51,56], we solve the above inhomogeneous linear differential equation using the Green's function G(p|p ′ ), which is governed by where δ represents the Dirac-delta function.Following [80], we expand this G into spherical harmonic series as Here C n,m are the series coefficients and Y m n (p) represent the spherical harmonics.Similarly, the Dirac-delta function can also be related to the spherical harmonics as , where Y m n represents the corresponding complex conjugate spherical harmonic [80].Eq. ( 37) is now expressible as Using the properties of spherical harmonics [80], we substitute ∇ 2 p Y m n = −n(n+1)Y m n , take the inner product with respect to Y M N (p) on both sides of Eq. ( 39), and use the orthogonality property to obtain where the normalization condition yields C 0,0 = 1/ √ 4π.We finally use the Green's function solution to find ψ (1) : In the limit of weak shear, this solution matches with the numerically obtained ψ for arbitrary Pe f , whose evaluation is discussed next.

Numerical solution valid for arbitrary Pe f
Now we detail the numerical solution of Eq. ( 2), which uses the decomposition of ψ as We substitute this expansion in Eq. ( 2) and use the properties of spherical harmonics to obtain Here H(Y m n ) is a collection of expressions obtained by simplifying ∇ p ⋅ ( ṗψ) .These are represented in terms of angular momentum operators for computational convenience [56] (detailed in Supplementary Note 3).We take an inner product with Y j i on both sides of Eq. (S8) and use the orthogonality property to obtain Since the n = 0 mode is already known from normalization, the above equations can be recasted as the following linear system of equations where C i,j (for i ≥ 1) is unknown: To find C i,j coefficients, the above equations are solved using the Clebsch-Gordon coefficient formulation, where first 100 modes in n are used.We note that although the numerical approach is valid for arbitrary Pe f , there is an indirect upper bound of weak viscoelasticity that we must adhere to, as we are using the SOF model.In nondimensional terms this bound is: De Pe f < 1.Thus, for De = 0.1, we explore the results within an upper bound of Pe f = 10.

E. Rheology
First we show the expanded version of the flow-induced component of the particle stress where a i is a shorthand notation for the ensemble of orientation moment of i th order: ⟨p ⊗i ⟩.Expanding the upper-convected derivative, we obtain This Σ F p matches with the Newtonian bulk stress for De = 0 [55].For near-equilibrium results (Pe f ≪ 1), the above expression is substituted in Eq. ( 5) and orientation distribution from Eq. ( 41) is used to perform the ensemble integral.Since the passive viscoelastic contribution is O(Pe 2 f ), it does not contribute at the O(Pe f ) calculations.The viscosity relates to the deviatoric particleinduced stress as φµ p = Σ p xy / γ, where φ = na 3 is the volume fraction.This viscometeric relation can be simplified to obtain the near-equilibrium viscosity ratio The last term in the above expression is negligible and is generated from the non-Newtonian component of Σ F p .Thus, the major modifications due to fluid's viscoelasticity scale as Pe 2 a , as also depicted in Fig. S4d.

C. Passive spheroid in weakly viscoelastic shear flow
The non-dimensional equation of motion for the orbits is This can be written in the trigonometric form as To quantify the log-rolling and compare it with Brunn [38], we use the following transformation [60]: Here C represents the orbit constant (C = 0 and C = ∞ being the log-rolling and shear-plane rotation orbits, respectively) and τ is essentially the time coordinate that is scaled with T /2π, where T is the Jeffery orbit.Using the above relation, we obtain the evolution equation of the orbit constant as We now compare our result with Brunn [38] who, for a rigid tri-dumbbell, showed that the particle slowly drifts towards the log-rolling orbit.We plot Eq. (S6) and Eq.(4.17) from Brunn in Supplementary Figure S3 and find a good match.Note that over one period the particle's orbit declines in a wiggly manner.Interestingly, these wiggles have a physical significance and can be characterized as orbit reduction with two distinct plateaus.Supplementary Figure S3b shows that the long plateau (around t/T = 1/4) occurs first and is the state where particle is in the vicinity of flow-vorticity plane (where it slows down and spends most of its time).Once the particle exits the alignment state, it flips quickly to the other side (−x axis).This flip generates the short plateau around t/T = 1/2. .

III. DETAILS OF KINETIC MODEL
Here we provide full expressions and further details involved in the computational steps of kinetic model.The spectral decomposition of ψ is carried out as Substituting this in the Smoluchowski equation of main text and using the properties of spherical harmonics yields Here H(Y m n ) is a term obtained by simplifying ∇ p ⋅ ( ṗψ) and is represented in terms of angular momentum operators for computational convenience [56].Using such transformations, Eq. (S9) yields The angular momentum operators are further simplified using the following relations [80] L y = L y |Y m n ⟩ = We substitute above relations in Eq. (S12) to simplify H(Y m n ) and use that in Eq. (S8).Next, we take an inner product with Y j i on both sides of Eq. (S8) and use the orthogonality property to obtain The first term is only non-zero for the zeroth mode (i = 0), which yields −1/ √ 4π.Since the C 0,0 = 1/ √ 4π is already known from normalization, the first term cancels the i = 0 mode of second term.Using these identities, Eq. ( 14) can be recasted as the following linear system of equations in terms of C i,j : where the right hand side is known.The two integral terms in the above set of equations contain the inner product of three spherical harmonics.These are represented in terms of Clebsch-Gordon coefficients that exploit the orthogonal relations between harmonics for computational ease [80]: The bracket terms here are the Clebsch-Gordan coefficients inbuilt in Mathematica 13.
Below we show the match between our results (obtained from the above kinetic model) in the limit De = 0 and those obtained by Saintillan [24,   Comparing the shear viscosity as a function of Pe f from our analysis in the Newtonian fluid limit with that of Ref. [24] for |Pea| 8πA = 1, Λ = 1, τ = 1, Dr = 0.1.Lines (blue is puller, red is pusher) indicate our viscosity ratio and dots correspond to Ref. [24], whose viscosity is defined here with respect to the particle volume fraction φ = na 3 .Ref. [24] used φ = nl 3 A.

IV. LOG-ROLLING STATE IN THE MICROSTRUCTURE
In Supplementary Figure S5, we illustrate how log-rolling state manifests in the microstructure for large (weakly Brownian) passive particles.We consider a larger particle of length 30µm, which corresponds to D r = 2 × 10 −4 s −1 .As discussed in the main text, Supplementary Figure S5a shows that when shear is weak, the isotropic distribution of rods is modified primarily by extensional component (E ∞ ) of the shear flow.Hence, the microstructure peaks near the extensional axis.As the shear rate increases, the rotational flow (Ω ∞ ) contributes and the peak shifts towards the flow axis.Supplementary Figure S5b shows that particle almost aligns with the flow direction (positive and negative x-axis).Further increasing Pe f spreads the distribution in the flow-vorticity plane (Supplementary Figure S5c) and, ultimately for Pe f > 10 3 , a peak develops along the vorticity axis, which is entirely reminiscent of log rolling (Supplementary Figure S5d,e).A further increase concentrates the orientation distribution near the vorticity axis (Supplementary Figure S5f).
t e x i t s h a 1 _ b a s e 6 4 = " b h J + f 5 3 R 3 x T L G r 2 y C Q C c O I j d l p 4 = " > A A A C J H i c b V D L S s Q w F E 1 9 O 7 6 q L t 0 E B 8 H V 0 I q o I M K g G 5 c K z g O m d U g z 6 R h M 0 p L c i q X 0 Y 9 z 4 K 2 5 c + M C F G 7 / F z G M x P g 6 E H M 6 5 N 7 n 3 R K n g B j z v 0 5 m a n p m d m 1 9 Y r C w t r 6 y u u e s b T Z N k m r IG T U S i 2 x E x T H D F G s B B s H a q G Z G R Y K 3 o 9 m z g t + 6 Y N j x R V 5 C n L J S k r 3 j M K Q E r d d 3 j I E p E z + T S X k W z v C 4 Cr m L I S 3 y C g 1 4 C R d A n U p I S 5 3 i y k J X d 4 r 7 s u l W v 5 g 2 B / x J / T K p o j I u u + 2 b f p J l k C q g g x n R 8 L 4 W w I B o 4 F a y s B J l h K a G 3 p M 8 6 l i o i m Q m L 4 Z I l 3 r F K D 8 e J t k c B H q q T H Q W R Z j C e r Z Q E b s x v b y D + 5 3 U y i I / C g q s 0 A 6 b o 6 K M 4 E x g S P E g M 9 7 h m F E R u C a G a 2 1 k x v S G a U L C 5 V m w I / u + V / 5 L m X s 0 / q O 1 f 7 l f r p + M 4 F t A W 2 k a 7 y E e H q I 7 O 0 Q V q I I o e 0 B N 6 Q a / O o / P s v D s f o 9 I p Z 9 y z i X 7 A + f o G p Z + m r g = = < / l a t e x i t > V 1 = ˙ ye x < l a t e x i t s h a 1 _ b a s e 6 4 = " K j P x 5 y U r P V M O 9 8 K Q D 3 i O 9 b c r 4 K 0 = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m k q M e i F 4 8 V T V t o Q 9 l s J + 3 S z S b s b o R S + h O 8 e F D E q 7 / I m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 M B V c G 9 f 9 d g p r 6 x u b W 8 X t 0 s 7 u 3 v 5 9 E i b I l D Z m r v y c m N N Z 6 H I e 2 M 6 Z m q J e 9 m f i f 1 8 l M d B 1 M u E w z g 5 I t F k W Z I C Y h s 7 9 J n y t k R o w t o U x x e y t h Q 6 o o M z a d k g 3 B W 3 5 5 l T Q v q t 5 l t X Z f q 9 R v 8 j i K c A K n c A 4 e X E E d 7 q A B P j A Y w D O 8 w p s j n B f n 3 f l Y t B a c f O Y Y / s D 5 / A E 9 V I 3 I < / l a t e x i t > U s < l a t e x i t s h a 1 _ b a s e 6 4 = " E M z Z G w J t Y 7 o V I h 0 Y K P v w w m A Z e 6 E = " > A A A B 9 X i c b V D L S g M x F L 1 T X 7 W + q i 7 d B I v g q s y I q M u i G 5 c V 7 A P a s W Q y m T Y 0 k w x J R i l D / 8 O N C 0 X c + i / u / B s z 7 S y 0 9 U D I 4 Z x 7 y c k J E s 6 0 c d 1 v p 7 S y u r a + U d 6 s b

Fig. 1 .
Fig. 1.Microswimmers in viscoelastic shear flow.Schematics showing (a) the dilute suspension of microswimmers in an external shear flow V ∞ of a viscoelastic fluid, (b) coordinate system moving with an individual swimmer that is modelled as an active prolate spheroid of major axis a and minor axis b.Here p and Us denote the orientation and swimming speed, respectively.

Fig. 3 .
Fig. 3. State diagram showing the different states of the orientational dynamics for three Weissenberg numbers.With increasing dipole strength σ from left to right, we observe alignment (Align.),shear-plane rotation (SPR), log rolling (LR) and again alignment.In the alignment state, the alignment angle ϕeq is plotted versus dimensionless activity σ as solid lines.The dots on the horizontal axis indicate the critical values σc at which alignment (Align.)occurs either to the left of the dashed-dotted lines (pullers, σ < 0) or to the right (pushers, σ > 0).The dashed lines separate SPR and LR from each other at σ = −1.93,−0.49, −0.19 for increasing Wi.The inset shows the variation of Lyapunov exponent (L) with σ.Other parameters are chosen as in Fig. 2.
t e x i t s h a 1 _ b a s e 6 4 = " 9 P H / J r 7 + H n Y t O o o Q H n 0 f O k V 0 9 x g = " > A A A B 8 n i c b V B N S 8 N A E N 3 U r 1 q / q h 6 9 L B b B U 0 m k q M e i F 4 8 V 7 A e 0 o W y 2 k 3 b p Z h N 2 J 2 I J / R l e P C j i 1 V / j z X / j t s 1 B W x 8 M P N 6 b Y W Z e k E h h 0 H W / n c L a + s b m V n G 7 t L O 7 t 3 9 Q 3 k l b w 4 6 L 8 6 7 8 7 F o L T j 5 z D H 5 A + f z B 6 r Q k Y Q = < / l a t e x i t > Pe f < l a t e x i t s h a 1 _ b a s e 6 4 = " X w + 5 D n J qI Y k n P 7 k n M T G D Y n L T / q Q = " > A A A B + H i c b V B N S 8 N A E N 3 U r 1 o / G v X o J V g E T y U pR b 0 I R S 8 e K 9 g P a E P Y b K f t 0 s 0 H u x O x h v 4 S L x 4 U 8 e p P 8 e a / c d P m o K 0 P B h 7 v z T A z z 4 8 F V 2 j b 3 0 Z h b X 1 j c 6 u 4 X d r Z 3 d s v m w e H b R U l k k G L R S K S X Z 8 q E D y E F n I U 0 I 0 l 0 M A X 0 P E n N 5 n f e Q C p e B T e 4 z Q G N 6 C j k A 8 5 o 6 g l z y z 3 E R 4 x b c L M o 1 c 1 8 / g D C o p K B < / l a t e x i t > Pe a = 20 (a)  ! " < l a t e x i t s h a 1 _ b a s e 6 4 = " 9 P H / J r 7 + H n Y t O o o Q H n 0 f O k V 0 9 x g = " > A A A B 8 n i c b V B N S 8 N A E N 3 U r 1 q / q h 6 9 L B b B U 0 m k q M e i F 4 8 V 7 A e 0 o W y 2 k 3 b p Z h N 2 J 2 I J / R l e P C j i 1 V / j z X / j t s 1 B W x 8 M P N 6 b Y W Z e k E h h 0 H W / n c L a + s b m V n G 7 t L O 7 t 3 9 Q 3 k l b w 4 6 L 8 6 7 8 7 F o L T j 5 z D H 5 A + f z B 6 r Q k Y Q = < / l a t e x i t > t e x i t s h a 1 _ b a s e 6 4 = " b TZ z + W w K w B 7 1 6 A C C b R Z G x B W + y j 0 = " > A A A C C n i c b V B N S 8 N A E N 3 4 W e t X 1 a O X a B G 8 W J J S 1 I t Q 9 O K x g v 2 A p o T N d t I u 3 W z C 7 k Q s o W c v / h U vH h T x 6 i / w 5 r 8 x a X u o r Q 8 G H u / N M D P P i w T X a F k / x t L y y u r a e m 4 j v 7 m 1 v b N b 2 N t v 6 D B W D O o s F K F q e V S D 4 B L q y F F A K 1 J A A 0 9 A 0 x v c Z H 7 z A Z T m o b z H Y Q S d g P Y k 9 z m j m E p u 4 c h B e M S k B i O X X p U t x 8 n P C m d l K H w s T Q z P L x e x y B Q z F M C W U K Z 7 e a r I + V Z R h m l 4 W g j 3 / 8 i J p l E v 2 e a l y V y l W r 6 d x 5 M g h O S a n x C Y X p E p u S Y 3 U C S N P 5 I W 8 k X f j 2 X g 1 P o z P S e u S M Z 0 5 I H 9 g f P 0 C P W u Z U Q = = < / l a t e x i t > Pea = 20 Pea = 20 < l a t e x i t s h a 1 _ b a s e 6 4 = " b T Z z + W w K w B 7 1 6 A C C b R Z G x B W + y j 0 = " > A A A C C n i c b V B N S 8 N A E N 3 4 W e t X 1 a O X a B G 8 W J J S 1 I t Q 9 O K x g v 2 A p o T N d t I u 3 W z C 7 k Q s o W c v / h U v H h T x 6 i / w 5 r 8 x a X u o r Q 8 G H u / N M D P P i w T X a F k / x t L y y u r a e m 4 j v 7 m 1 v b N b 2 N t v 6 D B W D O o s F K F q e V S D 4 B L q y F F A K 1 J A A 0 9 A 0 x v c Z H 7 z A Z T m o b z H Y Q S d g P Y k 9 z m j m E p u 4 c h B e M S k B i O X X p U t x 8 n P C m d l K H w s T Q z P L x e x y B Q z F M C W U K Z 7 e a r I + V Z R h m l 4 W g j 3 / 8 i J p l E v 2 e a l y V y l W r 6 d x 5 M g h O S a n x C Y X p E p u S Y 3 U C S N P 5 I W 8 k X f j 2 X g 1 P o z P S e u S M Z 0 5 I H 9 g f P 0 C P W u Z U Q = = < / l a t e x i t > t e x i t s h a 1 _ b a s e 6 4 = " T r Z q 5 r k j D d I k n M T k m b y S N w e d F + f d + V i 0 F p x 8 5 p j 8 g f P 5 A 6 M 8 k X 8 = < / l a t e x i t > Pe a < l a t e x i t s h a 1 _ b a s e 6 4 = "

20 Pe
I H M U O 3 k P 9 1 u x J M A Q u S S a d 1 y 7 B j b K V M o u I R R z k 0 0 x I z f s h 6 0 D A 1 Z A L q d T t 4 Z 0 a J R u t S P l K k Q 6 U T 9 P Z G y Q O t B 4 J n O g G F f z 3 t j 8 T + v l a B / 1 k 5 F G C c I I Z 8 u 8 h N J M a L j b G h X K O A o B 4 Y w r o S 5 l f I + U 4 y j S T B n Q n D m X 1 4 k 9 e O S c 1 I q X 5 c L l Y t Z H F l y Q A 7 J E X H I K a m Q K 1 I l N c L J A 3 k i L + T V e r S e r T f r f d q a s W Y z + + Q P r I 9 v 8 o O b Z Q = = < / l a t e x i t > |Pe a | = t e x i t s h a 1 _ b a s e 6 4 = " 9 P H / J r 7 + H n Y t O o o Q H n 0 f O k V 0 9 x g = " > A A A B 8 n i c b V B N S 8 N A E N 3 U r 1 q / q h 6 9 L B b B U 0 m k q M e i F 4 8 V 7 A e 0 o W y 2 k 3 b p Z h N 2 J 2 I J / R l e P C j i 1 V / j z X / j t s 1 B W x 8 M P N 6 b Y W Z e k E h h 0 H W / n c L a + s b m V n G 7 t L O 7 t 3 9 Q 3 k l b w 4 6 L 8 6 7 8 7 F o L T j 5 z D H 5 A + f z B 6 r Q k Y Q = < / l a t e x i t > Pe a =-20 Pe a =20 Puller < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 e D U 1 r J e W 3 a H k 2 v e 6 F u P / J M e s / g = " > A A A B + n i c b V B N S 8 N A E N 3 U r 1 q / U j 1 6 C R b B U 0 l E 1 I t Q 9 O K x g v 2 A N o T N d t I u 3 X y w O 1 F L 2 p / i x Y M i X v 0 l 3 v w 3 b t s c t P X B w O O 9 G W b m + Y n g C m 3 7 2 y i s r K 6 t b x Q 3 S 1 v b O 7 t 7 Z n m / q e J U M m i w W M S y 7
Fig.S4.Comparing the shear viscosity as a function of Pe f from our analysis in the Newtonian fluid limit with that of Ref.[24] for |Pea| 8πA = 1, Λ = 1, τ = 1, Dr = 0.1.Lines (blue is puller, red is pusher) indicate our viscosity ratio and dots correspond to Ref.[24], whose viscosity is defined here with respect to the particle volume fraction φ = na 3 .Ref.[24] used φ = nl 3 A.