1. Introduction
Understanding and relating bulk rheological properties of complex fluids to the orientations, deformations and interactions of their microscopic constituents has been a long-standing challenge in fluid mechanics (Larson Reference Larson2005; Schroeder Reference Schroeder2018). With the ability to visualize single molecules using fluorescence microscopy, over the last two decades, a large body of research using deoxyribonucleic acid (DNA) has focused on deciphering the dynamics and rheological properties of dilute long-chain polymer solutions in simple flows (Shaqfeh Reference Shaqfeh2005). The persistence length 
$\ell _p$ of DNA molecules is much shorter than their typical contour length 
$L$, and in this limit conformational properties are governed by a competition between entropic forces favouring coiled conformations and viscous forces tending to stretch the molecules, as quantified by the Weissenberg number 
$Wi\equiv \dot {\gamma }\tau _r$ or product of the applied strain rate 
$\dot {\gamma }$ with the longest polymer relaxation time 
$\tau _r$. Along with coarse-grained simulations (Jendrejack, de Pablo & Graham Reference Jendrejack, de Pablo and Graham2002) and kinetic theories (Winkler Reference Winkler2006), these molecular rheology studies have highlighted how microscopic conformational dynamics give rise to macroscopic rheological properties such as shear thinning and normal stress differences in shear flow (Hur, Shaqfeh & Larson Reference Hur, Shaqfeh and Larson2000; Schroeder et al. Reference Schroeder, Teixeira, Shaqfeh and Chu2005), and lead to the ‘coil-stretch’ transition in extensional flow (Schroeder et al. Reference Schroeder, Babcock, Shaqfeh and Chu2003).
In the other limit, the rheology of rigid Brownian rod-like suspensions is also well understood (Doi & Edwards Reference Doi and Edwards1988). Thermal fluctuations, in this case, result in orientational diffusion with rotational diffusivity 
$d_r = {3 k_B T \ln (2 r)}/{{\rm \pi} \mu L^3}$, where 
$k_B T$ is the thermal energy, 
$\mu$ is the shear viscosity of the solvent and 
$r = L/a$ is the aspect ratio of the rods with characteristic radius 
$a$. The competition between orientational diffusion favouring an isotropic distribution and the background shear that tries to align the rods is characterized by the rotary Péclet number 
$Pe \equiv \dot {\gamma } d_r^{-1}$, where 
$d_r^{-1}$ describes the orientational relaxation time. With increasing 
$Pe$, the preferential alignment of these rod-like polymers near the flow axis reduces viscous dissipation, which results in shear thinning and also yields non-zero normal stress differences. Three distinct scaling regimes for the shear viscosity and normal stress differences as functions of 
$Pe$ have been identified and characterized by the foundational theoretical analyses of Leal & Hinch (Reference Leal and Hinch1971), Hinch & Leal (Reference Hinch and Leal1972) and Brenner (Reference Brenner1974) (see § 2.4 and table 1 for a summary). However, so far, only rigid Brownian rods have been considered in detail, and the role of flow-induced deformations on the rheology of these suspensions remains largely unexplored.
Table 1. Asymptotic scalings for the relative viscosity 
$\eta$ and the individual stress components in dilute suspensions of rigid Brownian rods in various regimes of the rotational Péclet number 
$Pe$ and aspect ratio 
$r$ (Hinch & Leal Reference Hinch and Leal1972; Brenner Reference Brenner1974). The first and second normal stress differences 
$N_1>0$ and 
$N_2<0$ scale similarly as the diagonal stress components.
We address this problem here with a focus on stiff polymers characterized by 
$\ell _p \gg L$, the opposite limit compared with DNA. While in weak flows these filaments behave as rigid rods, they are known to undergo various buckling instabilities in stronger flows (Becker & Shelley Reference Becker and Shelley2001; Manikantan & Saintillan Reference Manikantan and Saintillan2015; Liu et al. Reference Liu, Chakrabarti, Saintillan, Lindner and du Roure2018; du Roure et al. Reference du Roure, Lindner, Nazockdast and Shelley2019; Chakrabarti et al. Reference Chakrabarti, Liu, LaGrone, Cortez, Fauci, du Roure, Saintillan and Lindner2020), yet clear insight into the specific role of these instabilities in the rheology of dilute suspensions is lacking. Here, we use numerical simulations to relate the morphological transitions of stiff Brownian filaments in simple shear flow to the rheology of their dilute suspensions. We also contrast our predictions with known results for non-Brownian deformable fibres (Becker & Shelley Reference Becker and Shelley2001; Tornberg & Shelley Reference Tornberg and Shelley2004) and uncover how they are altered by shape fluctuations and orientational diffusion.
The paper is organized as follows. In § 2, we provide details of the polymer model, measures of the extra stress and a brief summary of the scaling laws for Brownian rigid rods. We present numerical results for the rheology in both two and three dimensions in § 3, and discuss our predictions in the context of the rheology of rigid Brownian rods as well as non-Brownian elastic fibres.
2. Problem description and methodology
2.1. Governing equations
In the dilute limit, we simulate the dynamics of a single polymer modelled as a fluctuating, inextensible Euler elastica with centreline parametrized as 
${\boldsymbol x}(s,t)$, where 
$s$ is the arc length (Liu et al. Reference Liu, Chakrabarti, Saintillan, Lindner and du Roure2018). Hydrodynamics is captured by the local slender body theory for Stokes flow, in which the centreline position evolves as
\begin{equation} 8 {\rm \pi}\mu [\partial_t{\boldsymbol x}(s,t) - {\boldsymbol u}_\infty] ={-}\boldsymbol{\varLambda} \boldsymbol{\cdot} {\boldsymbol f}(s,t). \end{equation}
Here, 
${\boldsymbol u}_\infty = (\dot {\gamma } y, 0, 0)$ is the background shear flow with constant shear rate 
$\dot {\gamma }$. The force per unit length exerted by the filament on the fluid is modelled as 
${\boldsymbol f} = B {\boldsymbol x}_{ssss} - (\sigma {\boldsymbol x}_s)_s + {\boldsymbol f}^{{b}}$, where 
$B$ is the bending rigidity, 
$\sigma$ is a Lagrange multiplier that enforces inextensibility of the filament and can be interpreted as line tension, and 
${\boldsymbol f}^{{b}}$ is the Brownian force density obeying the fluctuation-dissipation theorem. The local mobility operator 
$\boldsymbol {\varLambda }$ accounts for drag anisotropy and is given by
\begin{equation} \boldsymbol{\varLambda}\boldsymbol{\cdot}{\boldsymbol f}= [(2-c) \boldsymbol{\mathsf{I}} - (c+2) {\boldsymbol x}_s {\boldsymbol x}_s] \boldsymbol{\cdot} {\boldsymbol f}, \end{equation}
where 
$c = \ln (\epsilon ^2 e) < 0$ is an asymptotic geometric parameter and 
$\epsilon = r^{-1}$ is the inverse aspect ratio. This geometrically nonlinear description of the centreline elasticity has been extensively used to successfully describe various elastohydrodynamic problems in both low- (Shelley & Ueda Reference Shelley and Ueda2000; Tornberg & Shelley Reference Tornberg and Shelley2004; Young & Shelley Reference Young and Shelley2007; Lim et al. Reference Lim, Ferent, Wang and Peskin2008; Manikantan & Saintillan Reference Manikantan and Saintillan2013; Chakrabarti et al. Reference Chakrabarti, Liu, LaGrone, Cortez, Fauci, du Roure, Saintillan and Lindner2020) and high-Reynolds-number flows (Allende, Henry & Bec Reference Allende, Henry and Bec2018; Banaei, Rosti & Brandt Reference Banaei, Rosti and Brandt2020). The present model is identical to the classical planar Euler elastica problem (Singh & Hanna Reference Singh and Hanna2019; Audoly & Pomeau Reference Audoly and Pomeau2000) and is a suitable description as long as the local curvature 
$\kappa \ll a^{-1}$. In the case of slender filaments used here for which 
$a\ll L$, such extreme deformations do not occur over the range of flow strengths considered, thus justifying the use of this model. Note that the mobility operator of (2.2) does not account for non-local hydrodynamic interactions between distant parts of the filaments. Including these interactions can be achieved using the non-local slender-body operator as done in our past work (Liu et al. Reference Liu, Chakrabarti, Saintillan, Lindner and du Roure2018). This results, however, in an increased computational cost that is prohibitive for the present study as it requires averaging over very long times. To test the consequences of this approximation, we have also performed a few select simulations with full hydrodynamics, where we observed only slight quantitative differences with the local drag model in terms of the magnitude of stresses, but no difference in the scalings of the various rheological quantities with respect to flow strength.
We scale lengths by 
$L$, time by the characteristic relaxation time of bending deformations 
$\tau _r = 8 {\rm \pi}\mu L^4/B$, elastic forces by the bending force scale 
$B/L^2$ and Brownian forces by 
$\sqrt {L/\ell _p} B/L^2$. The dimensionless equation of motion then reads
\begin{equation} \partial_t{\boldsymbol x}(s,t)= \bar{\mu} {\boldsymbol u}_\infty - \boldsymbol{\varLambda} \boldsymbol{\cdot} [{\boldsymbol x}_{ssss} - (\sigma {\boldsymbol x}_s)_s + \sqrt{L/\ell_p} \boldsymbol{\zeta} ], \end{equation}
where 
$\boldsymbol {\zeta }$ is a Gaussian random vector with zero mean and unit variance. Two dimensionless groups appear: (i) the elastoviscous number
\begin{equation} \bar{\mu} \equiv \frac{\dot{\gamma}\tau_r}{c}=\frac{8 {\rm \pi}\mu \dot{\gamma} L^4}{Bc} , \end{equation}
serves as the measure of hydrodynamic forcing against internal elasticity and plays a role analogous to the Weissenberg number for flexible polymers, and (ii) 
$L/\ell _p$ captures the importance of thermal shape fluctuations. The limit of rigid Brownian rods formally corresponds to 
$\bar {\mu }\rightarrow 0$ (no flow-induced deformations) and 
$L/\ell _p\rightarrow 0$ (no thermal shape fluctuations). We note that 
$\bar {\mu }$ and 
$L/\ell _p$ are related to the rotary Péclet number:
\begin{equation} Pe = \frac{\bar{\mu} c}{24 \ln (2 r)} \frac{\ell_p}{L},\end{equation}which will facilitate comparisons of our simulations of Brownian polymers including finite bending resistance with analytical predictions for rigid Brownian rods.
2.2. Numerical methods
Numerical schemes to solve (2.1) in the absence of Brownian motion have been described previously in great detail (Tornberg & Shelley Reference Tornberg and Shelley2004), and we provide only a brief outline here. In (2.1), the line tension 
$\sigma (s,t)$, which acts as a Lagrange multiplier, is an unknown. To solve for it, we make use of the inextensibility constraint 
${\boldsymbol x}_s \boldsymbol {\cdot } {\boldsymbol x}_s = 1$. Differentiating this constraint with respect to time and using the slender-body-theory equation provides a second-order linear ordinary differential equation for 
$\sigma$ as explained in detail by Tornberg & Shelley (Reference Tornberg and Shelley2004), which is subsequently solved with the boundary condition 
$\sigma = 0$ at 
$s = 0,1$. The time marching of (2.1) is performed using an implicit–explicit second-order accurate backward finite difference scheme, where the stiff linear terms arising from bending are treated implicitly while the nonlinear terms are handled explicitly. The boundary conditions for time marching are the force- and moment-free conditions for the filament ends, which translate to 
${\boldsymbol x}_{ss} = {\boldsymbol x}_{sss} = 0$ at 
$s = 0,1$. We used 
$N=64-128$ points to discretize the filament centreline and typical time steps were in the range of 
$\varDelta t \sim 10^{-6} - 10^{-9}$.
Treatment of the spatially and temporally uncorrelated Brownian forces in (2.1) requires special attention, and has been described previously by Manikantan & Saintillan (Reference Manikantan and Saintillan2013) and Liu et al. (Reference Liu, Chakrabarti, Saintillan, Lindner and du Roure2018). Specifically, we apply a low-pass filter to smooth out the noise along the centreline and typically remove 50 % of the high-frequency components in the process. The algorithm was benchmarked against standard equilibrium properties of semiflexible polymers (Wilhelm & Frey Reference Wilhelm and Frey1996).
2.3. Measures of stress
A calculation of the extra stress in a dilute suspension of force- and torque-free particles was provided by Batchelor (Reference Batchelor1970). The single-particle contribution to the bulk stress tensor in the dilute limit is given by the stresslet, which generalizes the Kirkwood formula commonly used for molecular systems (Irving & Kirkwood Reference Irving and Kirkwood1950). For our polymer model, the expression for the extra stress is
\begin{equation} \boldsymbol{\varSigma} ={-}n \left\langle\int_{0}^{L} \left[\tfrac{1}{2}({\boldsymbol x}{\boldsymbol f} +{\boldsymbol f}{\boldsymbol x})-\tfrac{1}{3} {\boldsymbol{\mathsf{I}}} ({\boldsymbol x} \boldsymbol{\cdot}{\boldsymbol f}) \right]\,{\mathrm{d}} s \right\rangle, \end{equation}
where 
$n$ is the number density in the suspension, 
${\boldsymbol f}$ is the dimensionless force density exerted on the fluid with contributions from both elastic deformations and Brownian fluctuations, and 
$\langle \cdot \rangle$ denotes the ensemble average. This expression is extremely convenient for non-Brownian fibres (Becker & Shelley Reference Becker and Shelley2001). However, in simulations of Brownian polymers, fluctuations have contributions of 
$O(\varDelta t^{-1/2})$, where 
$\varDelta t$ is the integration time step. These contributions also enter the Lagrange multiplier 
$\sigma (s,t)$ that enforces inextensibility, which results in a poor convergence of the ensemble average as first pointed out by Doyle, Shaqfeh & Gast (Reference Doyle, Shaqfeh and Gast1997). Because our interest is in the steady-state extra stress, we instead use the Giesekus stress expression commonly used for polymers (Öttinger Reference Öttinger1996; Doyle et al. Reference Doyle, Shaqfeh and Gast1997),
\begin{equation} \boldsymbol{\varSigma} ={-}n \left\langle\int_{0}^{L} \left[\tfrac{1}{2}({\boldsymbol x} {\boldsymbol{\mathsf{R}}} \boldsymbol{\cdot} {\boldsymbol u}_\infty +{\boldsymbol{\mathsf{R}}} \boldsymbol{\cdot} {\boldsymbol u}_\infty {\boldsymbol x})-\tfrac{1}{3} {\boldsymbol{\mathsf{I}}} ({\boldsymbol x} \boldsymbol{\cdot}{\boldsymbol{\mathsf{R}}} \boldsymbol{\cdot} {\boldsymbol u}_\infty) \right] \, {\mathrm{d}} s \right\rangle,\end{equation}
where 
${\boldsymbol{\mathsf{R}}} = \boldsymbol {\varLambda }^{-1}$ is the local resistance tensor along the centreline. Results obtained with (2.7) were tested against (2.6) in various regimes of 
$Pe$. In weak flows (
$Pe\ll 1$), the Giesekus stress was found to slightly underestimate the magnitude of the viscosity as it omits Brownian contributions, but excellent agreement was found in stronger flows and identical scalings in terms of 
$Pe$ were obtained by both methods across all regimes of 
$Pe$. In the following, we present results based on (2.7), which is computationally more efficient than the Kirkwood expression.
2.4. Summary of rigid rod rheology
A slender non-Brownian rod in shear flow undergoes a periodic tumbling motion known as a Jeffery orbit (Jeffery Reference Jeffery1922). During this periodic tumbling, the particle spends most of its time aligned with the flow direction and equal amounts of time in the extensional and compressional quadrants of the flow. This dynamics is fundamentally altered in the presence of rotational diffusion. While the shear flow still results in quasi-periodic tumbling, the Brownian rod is now able to stochastically sample different orbits (Zöttl et al. Reference Zöttl, Klop, Balin, Gao, Yeomans and Aarts2019), which results in an anisotropic orientational probability distribution 
$\psi ({\boldsymbol p})$ at steady state, where 
${\boldsymbol p}$ is a unit vector that identifies the orientation of the rod (Chen & Jiang Reference Chen and Jiang1999). This distribution leads to a mean orientation of the rod in the extensional quadrant, which gives rise to a contractile stresslet as the inextensible rod resists stretching by the flow. This stresslet in turn alters the effective viscosity of the system. In the dilute limit of 
$n L^3 \ll 1$, computing the extra stress reduces to obtaining the steady-state orientation distribution of a single Brownian rod. Contributions to the stresslet arise from the external flow and from Brownian diffusion, and can be computed using the slender body theory (Batchelor Reference Batchelor1970; Leal & Hinch Reference Leal and Hinch1971; Hinch & Leal Reference Hinch and Leal1972; Brenner Reference Brenner1974):
\begin{equation}
\boldsymbol{\varSigma}^{f} = \tfrac{{\rm \pi} \mu n L^3}{6 \ln(2
r)} [\langle\, {\boldsymbol p} {\boldsymbol p}
{\boldsymbol p} {\boldsymbol p} \rangle -
\tfrac{1}{3}{\boldsymbol{\mathsf{I}}}\langle {\boldsymbol
p} {\boldsymbol p} \rangle ] \boldsymbol{:}
{\boldsymbol{\mathsf{E}}}_\infty, \quad
\boldsymbol{\varSigma}^{b} = 3 n k_B T
[\langle\,{\boldsymbol p} {\boldsymbol p} \rangle -
\tfrac{1}{3}{\boldsymbol{\mathsf{I}}} ],
\end{equation}
where 
${\boldsymbol{\mathsf{E}}}_{\infty }$ is the rate-of-strain tensor of the applied flow 
${\boldsymbol u}_\infty$. The extra stress 
$\boldsymbol {\varSigma }=\boldsymbol {\varSigma }^{f}+\boldsymbol {\varSigma }^{b}$ is thus entirely determined from the second and fourth moments, 
$\langle\,{\boldsymbol p} {\boldsymbol p} \rangle$ and 
$\langle\,{\boldsymbol p} {\boldsymbol p} {\boldsymbol p} {\boldsymbol p} \rangle$, of rod orientations. These moments can be computed from the steady-state orientation distribution function 
$\psi (\,\boldsymbol {p})$, which is set by the balance of the advective rotational flux arising from the flow and of the Brownian diffusive flux. Hinch & Leal (Reference Hinch and Leal1972) solved for the distribution function and associated particle stress in three distinct asymptotic regimes of the rotary Péclet number 
$Pe$. These regimes and corresponding scalings are summarized in table 1.
The three rheological measures of primary interest to us are the relative polymer viscosity 
$\eta$, and the first and second normal stress differences 
$N_1$ and 
$N_2$,
\begin{equation} \eta = \frac{\varSigma_{xy}}{\mu \dot{\gamma} nL^3} , \quad N_1 = \frac{\varSigma_{xx}-\varSigma_{yy}}{n k_B T}, \quad N_2 = \frac{\varSigma_{yy}-\varSigma_{zz}}{n k_B T}. \end{equation}3. Numerical results and discussion
3.1. Three-dimensional rheology of stiff polymers
We present numerical results on the rheology in three dimensions, with a focus on the case of stiff slender polymers with 
$\ell _p/L = 1000$ and aspect ratio 
$r=220$. In this limit, the main effect of Brownian motion is to cause orientational diffusion with negligible shape fluctuations, and any deformations are thus the result of elastoviscous buckling. Figure 1(
$a$–
$c$) shows the relative polymer viscosity 
$\eta$, and first and second normal stress differences 
$N_{1,2}$ as functions of Péclet number 
$Pe$ (or, equivalently, elastoviscous number 
$\bar {\mu }$). In weak flows (
$Pe\ll 1$, pink region), 
$\eta$ exhibits a plateau whereas 
$N_{1,2}$ both grow from zero as 
$Pe^2$, with 
$N_1>0$ and 
$N_2<0$. With increasing flow strength, shear thinning takes place as the polymers start aligning with the flow, and a second regime begins with scalings of 
${\eta \sim Pe^{-1/3}}$ and 
${N_{1} \sim Pe^{2/3}}$ in perfect agreement with the theoretical predictions of table 1. The data for the second normal stress difference are very noisy in this range of 
$Pe$ and fail to capture the expected scaling of 
${Pe^{2/3}}$ for reasons we explain later. In this regime, the filament remains straight and tumbles quasi-periodically as shown in the first row of figure 1(
$f$).
Figure 1. (
$a$) Polymer viscosity 
$\eta$, (
$b$) first normal stress difference 
$N_1$, (
$c$) negative second normal stress difference 
$-N_2$, (
$d$) shear stress and (
$e$) normal stress components as functions of 
$Pe$ (bottom axis) and 
$\bar {\mu }$ (top axis) for polymers with 
$\ell _p/L=1000$ and 
$r=220$ in three dimensions. The low- and intermediate-
$Pe$ scalings are theoretical predictions for rigid rods, for which they are valid inside the pink (
$Pe\ll 1$) and green (
$1\ll Pe\ll r^3+r^{-3}$) regions. The vertical dashed lines show the onsets of 
$C$ buckling and 
$U$ turns (Liu et al. Reference Liu, Chakrabarti, Saintillan, Lindner and du Roure2018). (
$\,f$) Typical sequences of conformations during tumbling, 
$C$ buckling and a 
$U$ turn. Corresponding regimes and values of 
$Pe$ are labelled in (
$a$). In all figures, marked scalings before the buckling transition are rigid rod predictions (solid line), whereas scalings past the transition are numerical observations (dotted line).
As the filament performs a tumble, it rotates across the compressional quadrant of the flow, where it experiences compressive viscous stresses. Above a critical value of the elastoviscous number of 
$\bar {\mu }^{(1)} = 306.6$, these stresses can overcome bending rigidity and drive an Euler buckling instability leading to deformed configurations reminiscent of a 
$C$ shape, typical of the first mode of buckling (Becker & Shelley Reference Becker and Shelley2001). The filament then rotates as a 
$C$ and stretches out again as it enters and sweeps through the extensional quadrant. With increasing flow strength, the filament becomes more likely to buckle while remaining nearly aligned with the flow direction, and this ultimately gives rise to distinctive folded 
$U$-shaped conformations that perform tank-treading motions while maintaining a mean orientation close to the flow axis (Harasim et al. Reference Harasim, Wunderlich, Peleg, Kröger and Bausch2013). As uncovered in our past work (Liu et al. Reference Liu, Chakrabarti, Saintillan, Lindner and du Roure2018), the transition to this new mode of transport occurs at 
$\bar {\mu }^{(2)} = 1.8 \times 10^3$. Both of these thresholds are indicated by vertical lines in figure 1(
$a$–
$c$) and, for the chosen values of 
$\ell _p/L$ and 
$r$, fall within the intermediate scaling range of 
$1\ll Pe \ll r^3+r^{-3}$ (green region). Typical conformations from tumbling, 
$C$ buckling and 
$U$ turns are shown in figure 1(
$f$).
Quite remarkably, we find that the onset of buckling has no immediate signature on the rheology, with the Brownian rod-like scaling laws persisting past the critical value of 
$\bar {\mu }^{(1)}$. This result is in contrast with the two-dimensional (2-D) rheology of non-Brownian elastic fibres studied by Becker & Shelley (Reference Becker and Shelley2001) and Tornberg & Shelley (Reference Tornberg and Shelley2004), where buckling is responsible for shear thinning and non-zero normal stress differences. This discrepancy is attributed to the presence of three-dimensional (3-D) rotational diffusion in our simulations. In shear flow, the viscous compressive force experienced by the filament is a function of its orientation and reaches a maximum at an angle of 
$3 {\rm \pi}/4$ with the direction of flow in the plane of shear. In the presence of rotational noise, the filament orientation is not restricted to the shear plane and the maximum compression experienced is reduced. This translates to a set of measure zero for the probability density function 
$\psi (\,\boldsymbol {p})$ and, therefore, the probability of a buckling event is negligible at 
$\bar {\mu }=\bar {\mu }^{(1)}$. As 
$\bar {\mu }$ is increased beyond 
$\bar {\mu }^{(1)}$, buckling becomes increasingly more likely, and indeed 
$\eta$ and 
$N_{1,2}$ start to depart from the intermediate scalings before the onset of tank treading, as shown in figure 1. This departure marks a transition to new scalings of 
${\eta \sim {Pe}^{{-0.51\pm 0.02}}}$ and 
$N_1 \sim Pe^{{1 \pm 0.07}}$, and is accompanied by a sharp increase in 
$N_2$. These rheological changes are clear signatures of elastoviscous buckling, as they occur within the range of validity of the intermediate rigid rod scalings.
A more complete picture is provided in figure 1(
$d$,
$e$), which shows the shear and diagonal components of the extra stress tensor. In particular, we find that normal stresses in figure 1(
$b$) are dominated by 
$\varSigma _{xx}$, while 
$\varSigma _{yy}$ and 
$\varSigma _{zz}$, which are negative, have smaller magnitudes. All three components follow the same scaling, with the intermediate scaling of 
$Pe^{2/3}$ giving way to a nearly linear scaling into the buckling regime. Note that for 
$Pe\gg 1$, the values of 
$\varSigma _{yy}$ and 
$\varSigma _{zz}$ are almost identical, which explains the strong noise in the data for 
$N_2$ in figure 1(
$c$), especially in the intermediate scaling regime.
To relate these findings to filament conformations, we introduce the gyration tensor
\begin{equation} {\boldsymbol{\mathsf{G}}}(t) = \int_0^1 ({\boldsymbol x} - {\boldsymbol x}_c ) ({\boldsymbol x} - {\boldsymbol x}_c ) \, \mathrm{d}s, \end{equation}
whose dominant eigenvector is used to define the mean filament orientation 
$\theta (t)$ with respect to the flow direction. Its ensemble average is shown as a function of Péclet number in figure 2(
$a$). In weak flows (
$Pe\ll 1$), 
$\langle \theta \rangle$ asymptotes to the value of 
${\rm \pi} /4$ for an isotropic orientation distribution and, correspondingly, the diagonal components of 
$\langle {\boldsymbol{\mathsf{G}}}\rangle$, which describe the variance of the polymer mass distribution along the coordinate directions, all tend to the same value of 
$1/4{\rm \pi}$, as shown in figure 2(
$b$). Alignment of the filament with the flow is accompanied by a decrease in the mean orientation angle 
$\langle \theta \rangle$ with increasing shear rate, and is also indicated by the growth and saturation of 
$\langle {G}_{xx}\rangle$ while 
$\langle G_{yy}\rangle$ rapidly decays. Increasing flow strength also forces the filament toward the shear plane leading to a decrease in 
$\langle {G}_{zz}\rangle$. In strong flows, the initiation of 
$U$ turns leads to a sharper decrease in both 
$\langle G_{yy}\rangle$ and 
$\langle G_{zz}\rangle$, as the emergent folded conformations remain increasingly aligned with the flow direction, as seen in the third row of figure 1(
$f$).
Figure 2. (
$a$) Mean polymer orientation 
$\langle \theta \rangle$, defined as the angle made by the dominant eigenvector of the gyration tensor with the flow direction, as a function of 
$Pe$ and 
$\bar {\mu }$. (
$b$) Diagonal components 
$\langle G_{xx}\rangle$, 
$\langle G_{yy}\rangle$ and 
$\langle G_{zz}\rangle$ of the mean gyration tensor.
It is interesting to note that the scaling law for the viscosity is altered at a slightly lower value of 
$\bar {\mu }$ compared with the gyration tensor or mean orientation angle. This solidifies the idea that the signature of deformations observed in the stress and viscosity stems from the gradually more frequent occurrence of buckling events. During 
$C$ buckling (second row of figure 1
$f$), the deformed filament still tumbles as a whole, which leads to negligible changes in the behaviour of the gyration tensor and mean orientation angle, even though the stress components are affected.
3.2. Discussion
Above, we have discussed the case of stiff polymers and have compared our results with known scaling laws for rigid rods in three dimensions. Stiff polymers and rods both experience strong rotational diffusion. Shape fluctuations are absent in the case of rigid rods and remain weak for the stiff polymers in our simulations performed in the limit of 
$\ell _p/L\gg 1$. The main difference between the two systems is thus the occurrence of buckling instabilities above a given threshold for the stiff polymers. Numerical results show a clear signature of such elastoviscous buckling on the rheology, with enhanced shear thinning and normal stress differences compared with the case of rigid rods.
We now address the role of shape fluctuations on the rheology. Remaining in the limit of weak fluctuations for stiff polymers, we can tune fluctuations by varying 
$\ell _p/L$ while keeping rotational diffusion important, with a smaller 
$\ell _p/L$ corresponding to stronger shape fluctuations. For simplicity, we present here results obtained in two dimensions, and show in figure 3 the variation of the shear viscosity 
$\eta$ and normal stress difference 
$N_1$ as functions of shear rate for two combinations of (
$\ell _p/L$, 
$r$). Note that the theoretical threshold for buckling is identical in two dimensions and three dimensions, because the first instance of buckling occurs for a filament lying in the shear plane. The theoretical onset of 
$U$ turns was predicted in our past work using a 2-D reduced-order model (Liu et al. Reference Liu, Chakrabarti, Saintillan, Lindner and du Roure2018), but also faithfully describes the transition in three dimensions because the dynamics is primarily 2-D in strong flows, as previously shown in figure 2(
$b$). Scaling laws have been determined for both 
$\eta$ and 
$N_1$ before and after the buckling threshold, with a transition region where the scaling is changing continuously. Before the threshold, the measured scalings are in agreement with 3-D rigid rod predictions within error bars, and quantitatively identical results are obtained for both values of 
$\ell _p/L$ because elasticity plays a negligible role in that regime. As the Péclet number increases, buckling occurs first for the smaller value of 
$\ell _p/L$, because 
$Pe$ and 
$\ell _p/L$ are related by (2.5) and the buckling threshold occurs at a fixed 
$\bar {\mu }=\bar {\mu }^{(1)}$. Beyond the threshold, close agreement is found between the 2-D and 3-D results, and identical scalings are obtained for 
$\eta$ and 
$N_1$ for both values of 
$\ell _p/L$, which indicates that varying the importance of shape fluctuations, while remaining in the limit of weak fluctuations, does not alter the observed rheology noticeably. The limit of strong shape fluctuations, relevant to a number of experimental systems and to previous work (Harasim et al. Reference Harasim, Wunderlich, Peleg, Kröger and Bausch2013; Liu et al. Reference Liu, Chakrabarti, Saintillan, Lindner and du Roure2018), is outside the theoretical and numerical frameworks developed here.
Figure 3. (
$a$) Polymer viscosity 
$\eta$, and (
$b$) normal stress difference 
$N_1$ as functions of 
$Pe$ for 
$(\ell _p/L,r)=(1000,220)$ and 
$(50,100)$ in 2-D simulations. Vertical dashed lines show the two buckling thresholds, which correspond to the same value of 
${\bar {\mu }=\bar {\mu }^{(1)}}$.
Another observation can be made from figure 3. In three dimensions in figure 1, the changes in the scaling laws of viscosity and normal stress differences do not occur right at the theoretical onset of buckling, but are delayed owing to the presence of rotational diffusion, as the probability for the polymer to align perfectly with the direction of maximum compression is negligible right at the buckling threshold. On the contrary, in two dimensions, any polymer performing a tumble is required to sweep through the entire compressional quadrant and is therefore more likely to buckle. This interpretation is confirmed by our 2-D results shown in figure 3, in which the change in scaling owing to deformations now occurs slightly closer to the theoretical buckling threshold. Nevertheless, and in contrast to fully non-Brownian systems (Becker & Shelley Reference Becker and Shelley2001), the transition is not abrupt, likely owing to the presence of Brownian noise. The fact that the slope changes occur closer to the buckling threshold confirms again that the occurrence of 
$C$ buckling, rather than 
$U$ turns, is responsible for the change in behaviour. In addition, we find that the emergence of 
$U$ turns is not accompanied by another clear change in scaling. Rather, both buckling events and 
$U$ turns have the same signature on the rheology, in spite of their distinct morphologies and dynamics.
Our results can also be discussed in comparison to previous observations made for non-Brownian elastic fibres in two dimensions. For non-Brownian fibres, shape fluctuations as well as rotational diffusion are absent. Tornberg & Shelley (Reference Tornberg and Shelley2004) and Becker & Shelley (Reference Becker and Shelley2001) predicted the onset of normal stress differences and shear thinning above the buckling threshold from 2-D simulations. In contrast to the case of Brownian rods, normal stress differences are zero and the shear viscosity is constant in the absence of buckling. No scaling laws for 
$\eta$ and 
$N_1$ exist for non-Brownian fibres in the dilute limit.
Simulations performed with non-Brownian fibres require the use of an initial perturbation for the buckling instability to be triggered, and the properties chosen might influence the obtained results. Random fluctuations are naturally present in Brownian systems, as investigated here, and thus no initial perturbation needs to be imposed. To connect our Brownian simulations more directly to the non-Brownian case, we perform a numerical experiment, which focuses for the sake of illustration on a regime where 
$\bar {\mu } > \bar {\mu }^{(2)}$ and where the dynamics is dominated by 
$U$ turns with occasional 
$S$-shaped modes (Liu et al. Reference Liu, Chakrabarti, Saintillan, Lindner and du Roure2018). In this experiment, we perform a standard Brownian simulation with 
$\ell _p/L=50$ but artificially switch-off Brownian forces 
$\boldsymbol {f}^b$ at the initiation of any buckling event, detected by the threshold of 
$R_{ee}/L < 0.98$ on the end-to-end distance. In this way, we produce initial conditions set by Brownian noise and thus identical to those present in our Brownian simulations, but can investigate the arising stresses for a situation without Brownian noise. Stresses are only evaluated once the noise has been switched off. This artificial situation does not allow a meaningful viscosity to be calculated and thus we show 
$\varSigma _{xy}$ and 
$N_1$ in figure 4. These stresses are plotted versus 
$\bar {\mu }$, which is more relevant than 
$Pe$ in the case of non-Brownian fibres that do not experience orientational diffusion. Recall that 
$\bar {\mu }$ and 
$Pe$ are proportional to each other in (2.5), while 
$\eta$ and 
$\varSigma _{xy}$ are related through 
$\dot {\gamma }$ in (2.9a–c).
Figure 4. Numerical experiment illustrating the respective roles of elastic instabilities and thermal fluctuations. Perturbations for these simulations are obtained from a Brownian filament with 
$\ell _p/L=50$. Variation of (
$a$) shear stress and (
$b$) first normal stress difference versus 
$\bar {\mu }$.
The scalings obtained for 
$\varSigma _{xy}$ and 
$N_1$, with respect to 
$\bar {\mu }$, translate into 
$\eta \sim \varSigma _{xy}/\dot {\gamma } \sim Pe^{-0.56 \pm 0.03}$ and 
$N_1 \sim Pe^{0.87 \pm 0.06}$. Surprisingly, these exponents differ only slightly from those observed in the fully Brownian simulations in figure 3. As this numerical experiment does not account for the Brownian stress that arises from shape fluctuations and orientational diffusion during buckling, our results suggest that the leading contribution to the stress and its scaling with flow strength is set by elastic instabilities, with fluctuations mainly serving to trigger polymer tumbles, as in the rigid rod case, while also exciting the dominant buckling modes. Note that Becker & Shelley (Reference Becker and Shelley2001) found 
$N_1<0$ for non-Brownian fibres in sufficiently strong flows, in disagreement with our findings. We believe this discrepancy may arise from the particular way the fibre backbones were perturbed in their simulations.
Our numerical simulations have shed light on the role of elastohydrodynamic instabilities on the rheology of dilute suspensions of stiff polymers in the limit of 
$\ell _p\gg L$. The leading effect of buckling was shown to enhance shear thinning in strong flows while also driving an increase in normal stress differences in comparison to the case of rigid Brownian rods. Detailed rheological measurements in the dilute regime and in monodisperse systems are a challenge and have yet to be performed, but would be of great use to confirm our numerical predictions and connect them with past observations in more concentrated systems (Huber et al. Reference Huber, Harasim, Wunderlich, Kröger and Bausch2014; Kirchenbuechler et al. Reference Kirchenbuechler, Guu, Kuraniawan, Koenderink and Lettinga2014; Lang et al. Reference Lang, Kohlbrecher, Porcar, Radulescu, Sellinghoff, Dhont and Lettinga2019). Extensions of the present work may consider the case of semi-flexible polymers with 
$L\sim \ell _p$, which is more challenging numerically, as well as the role of hydrodynamic interactions in semi-dilute suspensions.
Acknowledgements
The authors thank M. Shelley for useful discussions.
Funding
A.L., B.C. and Y.L. acknowledge funding from the ERC Consolidator Grant PaDyFlow (Agreement 682367). D.S. acknowledges funding from a Paris Sciences Chair from ESPCI Paris.
Declaration of interest
The authors report no conflict of interest.
Open access
