Signatures of elastoviscous buckling in the dilute rheology of stiff polymers

As a stiff polymer tumbles in shear flow, it experiences compressive viscous forces that can cause it to buckle and undergo a sequence of morphological transitions with increasing flow strength. We use numerical simulations to uncover the effects of these transitions on the steady shear rheology of a dilute suspension of stiff polymers. Our results agree with classic scalings for Brownian rods in relatively weak flows but depart from them above the buckling threshold. Signatures of elastoviscous buckling include enhanced shear thinning and an increase in the magnitude of normal stress differences. We discuss our findings in the light of past work on rigid Brownian rods and non-Brownian elastic fibers and highlight the subtle role of thermal fluctuations in triggering instabilities.


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 2005;Schroeder 2018). 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 2005). The persistence length 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 † Email address for correspondence: dstn@ucsd.edu  Table 1. Asymptotic scalings for the relative viscosity η 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 1972;Brenner 1974). The first and second normal stress differences N 1 > 0 and N 2 < 0 scale similarly as the diagonal stress components.
conformations and viscous forces tending to stretch the molecules, as quantified by the Weissenberg number Wi ≡γ τ r or product of the applied strain rateγ with the longest polymer relaxation time τ r . Along with coarse-grained simulations (Jendrejack, de Pablo & Graham 2002) and kinetic theories (Winkler 2006), 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 2000;Schroeder et al. 2005), and lead to the 'coil-stretch' transition in extensional flow (Schroeder et al. 2003).
In the other limit, the rheology of rigid Brownian rod-like suspensions is also well understood (Doi & Edwards 1988). Thermal fluctuations, in this case, result in orientational diffusion with rotational diffusivity d r = 3k B T ln(2r)/πμL 3 , where k B T is the thermal energy, μ 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 ≡γ d −1 r , where d −1 r 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 (1971), Hinch & Leal (1972) and Brenner (1974) (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.
We address this problem here with a focus on stiff polymers characterized by p 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 2001;Manikantan & Saintillan 2015;Liu et al. 2018;du Roure et al. 2019;Chakrabarti et al. 2020), 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 2001;Tornberg & Shelley 2004) 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.

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 x(s, t), where s is the arc length (Liu et al. 2018). Hydrodynamics is captured by the local slender body theory for Stokes flow, in which the centreline position evolves as (2.1) Here, u ∞ = (γ y, 0, 0) is the background shear flow with constant shear rateγ . The force per unit length exerted by the filament on the fluid is modelled as where B is the bending rigidity, σ is a Lagrange multiplier that enforces inextensibility of the filament and can be interpreted as line tension, and f b is the Brownian force density obeying the fluctuation-dissipation theorem. The local mobility operator Λ accounts for drag anisotropy and is given by where c = ln( 2 e) < 0 is an asymptotic geometric parameter and = 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- ( and is a suitable description as long as the local curvature κ a −1 . In the case of slender filaments used here for which a 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. 2018). 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 τ r = 8πμL 4 /B, elastic forces by the bending force scale B/L 2 and Brownian forces by L/ p B/L 2 . The dimensionless equation of motion then reads where ζ is a Gaussian random vector with zero mean and unit variance. Two dimensionless groups appear: (i) the elastoviscous number 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/ p captures the importance of thermal shape fluctuations. The limit of rigid Brownian rods formally corresponds toμ → 0 (no flow-induced deformations) and L/ p → 0 (no thermal shape fluctuations). We note thatμ and L/ p are related to the rotary Péclet number: 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 2004), and we provide only a brief outline here. In (2.1), the line tension σ (s, t), which acts as a Lagrange multiplier, is an unknown. To solve for it, we make use of the inextensibility constraint x s · 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 σ as explained in detail by Tornberg & Shelley (2004), which is subsequently solved with the boundary condition σ = 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 x ss = 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 Δt ∼ 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 (2013) and Liu et al. (2018). 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 1996).

Measures of stress
A calculation of the extra stress in a dilute suspension of force-and torque-free particles was provided by Batchelor (1970). 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 1950). For our polymer model, the expression for the extra stress is where n is the number density in the suspension, f is the dimensionless force density exerted on the fluid with contributions from both elastic deformations and Brownian fluctuations, and · denotes the ensemble average. This expression is extremely convenient for non-Brownian fibres (Becker & Shelley 2001). However, in simulations of Brownian polymers, fluctuations have contributions of O(Δt −1/2 ), where Δt is the integration time step. These contributions also enter the Lagrange multiplier σ (s, t) that enforces inextensibility, which results in a poor convergence of the ensemble average as first pointed out by Doyle, Shaqfeh & Gast (1997). Because our interest is in the steady-state extra stress, we instead use the Giesekus stress expression commonly used for polymers (Öttinger 1996;Doyle et al. 1997), where R = Λ −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 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.

Summary of rigid rod rheology
A slender non-Brownian rod in shear flow undergoes a periodic tumbling motion known as a Jeffery orbit (Jeffery 1922). 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. 2019), which results in an anisotropic orientational probability distribution ψ(p) at steady state, where p is a unit vector that identifies the orientation of the rod (Chen & Jiang 1999). 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 nL 3 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 1970;Leal & Hinch 1971;Hinch & Leal 1972;Brenner 1974): where E ∞ is the rate-of-strain tensor of the applied flow u ∞ . The extra stress Σ = Σ f + Σ b is thus entirely determined from the second and fourth moments, pp and pppp , of rod orientations. These moments can be computed from the steady-state orientation distribution function ψ( p), which is set by the balance of the advective rotational flux arising from the flow and of the Brownian diffusive flux. Hinch & Leal (1972) 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 η, and the first and second normal stress differences N 1 and N 2 ,

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 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 η, and first and second normal stress differences N 1,2 as functions of Péclet number Pe (or, equivalently, elastoviscous numberμ). In weak flows (Pe 1, pink region), η 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 η ∼ Pe −1/3 and N 1 ∼ 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 ).
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μ (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 2001). 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. 2013). As uncovered in our past work (Liu et al. 2018), the transition to this new mode of transport occurs at μ (2) = 1.8 × 10 3 . Both of these thresholds are indicated by vertical lines in figure 1(a-c) and, for the chosen values of p /L and r, fall within the intermediate scaling range of 1 Pe 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 μ (1) . This result is in contrast with the two-dimensional (2-D) rheology of non-Brownian elastic fibres studied by Becker & Shelley (2001) and Tornberg & Shelley (2004), 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π/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 ψ( p) and, therefore, the probability of a buckling event is negligible atμ =μ (1) . Asμ is increased beyondμ (1) , buckling becomes increasingly more likely, and indeed η 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 η ∼ Pe −0.51±0.02 and N 1 ∼ Pe 1±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 Σ xx , while Σ yy and Σ 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 1, the values of Σ yy and Σ 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 whose dominant eigenvector is used to define the mean filament orientation θ(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 1), θ asymptotes to the value of π/4 for an isotropic orientation distribution and, correspondingly, the diagonal components of G , which describe the variance of the polymer mass distribution along the coordinate directions, all tend to the same value of 1/4π, as shown in figure 2(b). Alignment of the filament with the flow is accompanied by a decrease in the mean orientation angle θ with increasing shear rate, and is also indicated by the growth and saturation of G xx while G yy rapidly decays. Increasing flow strength also forces the filament toward the shear plane leading to a decrease in G zz . In strong flows, the initiation of U turns leads to a sharper decrease in both G yy and G zz , as the emergent folded conformations remain increasingly aligned with the flow direction, as seen in the third row of figure 1(f ).
It is interesting to note that the scaling law for the viscosity is altered at a slightly lower value ofμ 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 1f ), 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.

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 p /L 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 p /L while keeping rotational diffusion important, with a smaller 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 η and normal stress difference N 1 as functions of shear rate for two combinations of ( 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. 2018), 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 η 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 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 p /L, because Pe and p /L are related by (2.5) and the buckling threshold occurs at a fixedμ =μ (1) . Beyond the threshold, close agreement is found between the 2-D and 3-D results, and identical scalings are obtained for η and N 1 for both values of 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. 2013;Liu et al. 2018), is outside the theoretical and numerical frameworks developed here.
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 2001), 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 (2004) and Becker & Shelley (2001) 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 η 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 wherē μ >μ (2) and where the dynamics is dominated by U turns with occasional S-shaped modes (Liu et al. 2018). In this experiment, we perform a standard Brownian simulation with p /L = 50 but artificially switch-off Brownian forces 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 Σ xy and N 1 in figure 4. These stresses are plotted versusμ, which is more relevant than Pe in the case of non-Brownian fibres that do not experience orientational diffusion. Recall thatμ and Pe are proportional to each other in (2.5), while η and Σ xy are related throughγ in (2.9a-c).
The scalings obtained for Σ xy and N 1 , with respect toμ, translate into η ∼ Σ xy /γ ∼ Pe −0.56±0.03 and N 1 ∼ Pe 0.87±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 (2001) 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 p 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. 2014;Kirchenbuechler et al. 2014;Lang et al. 2019). Extensions of the present work may consider the case of semi-flexible polymers with L ∼ p , which is more challenging numerically, as well as the role of hydrodynamic interactions in semi-dilute suspensions.