1. Introduction
Microorganisms use a variety of techniques to swim in low Reynolds number environments (Pelczar, Gunsalus & Stanier Reference Pelczar, Gunsalus and Stanier1960; Berg Reference Berg2003). Organisms such as Escherichia coli use multiple rotating flagellar filaments that are connected to a rotary motor by a flexible short hook (Berg Reference Berg2003). Spermatozoa have flagellar waveforms that are nearly planar (Simons, Fauci & Cortez Reference Simons, Fauci and Cortez2015). Other organisms, such as Chlamydomonas reinhardtii, use an undulating motion to propel themselves forward (Lauga & Powers Reference Lauga and Powers2009). Although the efficacy of these motions is well understood in Newtonian fluids, in nature, these organisms frequently navigate viscoelastic environments composed of biomolecular proteins and polymeric networks solvated in a viscous Newtonian solvent. Viscoelastic dispersions induce elastic and viscous, frequency- and amplitude-dependent responses to imposed stresses and strains. This is a completely different fluid–structure interaction system than swimming in viscous fluids in which the fluid response is instantaneous: stored elastic stress and strain in the viscoelastic fluid are released over time on spatial and temporal scales dictated by the underlying structure and organization of the macromolecular species. To avoid the expense of modelling the polymeric network at the molecular scale, the ‘extra stress’ stored by the fluid's polymeric network is often modelled using continuum constitutive equations (Bird, Armstrong & Hassager Reference Bird, Armstrong and Hassager1987). One such choice is the class of upper convected Maxwell-like constitutive models such as the Oldroyd-B model (Larson Reference Larson1988; Beris & Edwards Reference Beris and Edwards1994; Morozov & Spagnolie Reference Morozov and Spagnolie2015), which has served as an idealized testbed for numerical method development, due to the simplicity of the material model. The Oldroyd-B model captures some features of polymeric fluids, such as the generation of normal stresses in pure shear, but it fails to capture the prevalent viscoelastic property of shear thinning. Generalizations of the Oldroyd-B model, such as the Giesekus model, were developed at the scale of polymeric kinetic theory and then coarse grained via moment closure analysis to arrive at a constitutive law for the second moment of the configurational probability distribution, the so-called extra-stress tensor (Larson Reference Larson1988). A more sophisticated polymeric kinetic theory is the Rouse linear entangled polymers (Rolie–Poly) model for entangled polymer solutions, which along with the Giesekus model and other constitutive equations, resolves shear thinning as well as other nonlinear responses (Bird et al. Reference Bird, Armstrong and Hassager1987).
The physics of locomotion of microorganisms has been extensively studied in Newtonian fluids (Lauga & Powers Reference Lauga and Powers2009), but a complete understanding of the motion of microorganisms in a viscoelastic fluid is still lacking. There is a deep literature for non-Newtonian swimming, with results that heavily depend on both the fluid model and swimmer model. A comprehensive review of motility in non-Newtonian fluids was performed by Li, Lauga & Ardekani (Reference Li, Lauga and Ardekani2021). For undulating swimmers, Lauga (Reference Lauga2007) showed that an infinite undulating sheet is always hindered by viscoelasticity. In contrast, a finite sheet can achieve a speed boost if its body shape and body elasticity is tuned to the fluid elasticity (Teran, Fauci & Shelley Reference Teran, Fauci and Shelley2010; Thomases & Guy Reference Thomases and Guy2014). Fu, Wolgemuth & Powers (Reference Fu, Wolgemuth and Powers2009) extended the analysis to three spatial dimensions and found that, generally, nonlinear viscoelasticity decreases the swimming velocity for small-amplitude waves. In contrast, Riley & Lauga (Reference Riley and Lauga2015) showed that the swimming speed can increase provided that kinematic waves travel in opposing directions along the sheet. Li et al. (Reference Li, Qin, Gopinath, Arratia, Thomases and Guy2017) and Hyakutake, Sato & Sugita (Reference Hyakutake, Sato and Sugita2019) found that organisms adopt a different gait in non-Newtonian fluids to enhance their swimming speed and efficiency, although Hyakutake et al. suggested the shear thinning provides a greater boost than viscoelasticity. For shear-thinning fluids, undulatory sheets can move more efficiently than in Newtonian counterparts (Vélez-Cordero & Lauga Reference Vélez-Cordero and Lauga2013; Gagnon, Keim & Arratia Reference Gagnon, Keim and Arratia2014), and these effects were deemed dominant over the hindrance from viscoelasticity, which is consistent with Hyakutake et al. In contrast, Dasgupta et al. (Reference Dasgupta, Liu, Fu, Berhanu, Breuer, Powers and Kudrolli2013) showed that there are some regimes in which viscoelasticity and shear thinning can enhance organism motility cooperatively, depending on the fluid. Additionally, there is also evidence that continuum models do not capture the complete story and that discrete spring-dashpot models of viscoelasticity are needed to accurately understand the dynamics of swimming in these regimes (Wróbel et al. Reference Wróbel, Lynch, Barrett, Fauci and Cortez2016; Li et al. Reference Li, Lauga and Ardekani2021; Schuech, Cortez & Fauci Reference Schuech, Cortez and Fauci2022).
The literature for helical swimmers is similarly complicated. For small pitch angle infinite swimmers, viscoelasticity always results in a decrease in swimming speed (Fu et al. Reference Fu, Wolgemuth and Powers2009; Li & Spagnolie Reference Li and Spagnolie2015). However, for helical swimmers with a large pitch angle, the resulting swimming speed becomes nonlinear in the Deborah number. Liu, Powers & Breuer (Reference Liu, Powers and Breuer2011) found that, experimentally, helical swimmers in viscoelastic fluids have an enhanced swimming speed in fluids for which the relaxation time matches that of the rotation time of the helix. Spagnolie, Liu & Powers (Reference Spagnolie, Liu and Powers2013) determined that enhanced swimming speed depends on many complex factors, such as the helical geometry, the material properties of the fluid and the rotation rate. When considering shear-thinning fluids, Gómez et al. (Reference Gómez, Godínez, Lauga and Zenit2017) and Demir et al. (Reference Demir, Lordi, Ding and Pak2020) suggested that shear thinning results in enhanced helical swimming speeds, and that the only argument consistent with their findings is a confinement-like effect on the swimmer. Further, their results suggest that shear thinning is the dominant effect in speed enhancement because of the magnitude of speed gains over viscoelasticity. Li & Spagnolie (Reference Li and Spagnolie2015) found that the introduction of a confining cylindrical tube around a swimmer significantly enhances the swimming speed, although the details depend on the helical pitch. Qu & Breuer (Reference Qu and Breuer2020) determined experimentally that shear-thinning behaviour is dominant over elasticity, at least in weakly elastic fluids in which the solvent shear thins.
Herein, we determine conditions for increased vs decreased efficiency of a helical flagellum rotating in shear-thinning vs non-shear-thinning fluids. We show that, for high elasticity and low shear thinning, the swimmer can more effectively pump fluid than in a comparable Newtonian fluid. Further, we contrast our results with prior work and suggest that the underlying model for shear thinning is important. To understand the efficacy of a swimmer, one must characterize the material properties of the fluid in which they are swimming and determine an appropriate model for that particular fluid. We additionally quantify the geometric properties the swimmer takes on in these different fluids and demonstrate that the viscoelastic fluid has only mild effects on the geometry of the swimmer.
Elastic bacterial flagella can be effectively modelled using Kirchhoff rod theory because the flagellum is long (
${\approx }10\,\mathrm {\mu }{\rm m}$) but very thin (
${\approx }20\,{\rm nm}$). Kirchhoff rod theory describes the forces and torques generated by an elastic rod in terms of the position of the centreline and an orthonormal set of director vectors attached to the centreline (Lim et al. Reference Lim, Ferent, Wang and Peskin2008). Goldstein et al. (Reference Goldstein, Goriely, Huber and Wolgemuth2000) introduced a bistable energy formulation for Kirchhoff rod theory to permit two stable helical configurations. Darnton & Berg (Reference Darnton and Berg2007) used Kirchhoff rod theory to estimate the bending rigidity of the flagellum by fitting a Kirchhoff rod model to experimental data. Lim & Peskin (Reference Lim and Peskin2012) incorporated hydrodynamic interactions by creating a generalized immersed boundary (IB) method based on Kirchhoff rod theory. Ko et al. (Reference Ko, Lim, Lee, Kim, Berg and Peskin2017) incorporated a bistable energy formulation into the generalized immersed boundary method to model polymorphic transformation.
In this study, we use the generalized IB method developed by Lim & Peskin (Reference Lim and Peskin2012) and Griffith & Lim (Reference Griffith and Lim2012) to simulate a flagellum in a shear-thinning viscoelastic fluid. We assess the pumping efficiency of the flagellum and determine the shape assumed by the flagellum in different fluids.
2. Mathematical model
We consider the motion of a single helical flagellum immersed in an incompressible Giesekus fluid. The fluid is described by the Cauchy stress tensor 
, which consists of stress from the Newtonian solvent 
and stress from the embedded polymers 
. The fluid equations of motion are
\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}\mathopen{}\left(\boldsymbol{x},t\right)\mathclose{} = 0, \end{gather}
\begin{gather}\mathop{\mathbb{C}}\limits^{\triangledown}\mathopen{}\left(\boldsymbol{x},t\right)\mathclose{} = -\frac{1}{\lambda}\mathopen{}\left(\mathbb{C}\mathopen{}\left(\boldsymbol{x},t\right)\mathclose{}-\mathbb{I}\right)\mathclose{} - \frac{\alpha}{\lambda}\mathopen{}\left(\mathbb{C}\mathopen{}\left(\boldsymbol{x},t\right)\mathclose{}-\mathbb{I}\right)\mathclose{}^2 + D\Delta\mathopen{}\left(\mathbb{C}\mathopen{}\left(\boldsymbol{x},t\right)\mathclose{} - \mathbb{I}\right)\mathclose{}, \end{gather}
in which 
$\mathbb {D} = \frac {1}{2}\mathopen {}\left (\boldsymbol {\nabla }\boldsymbol {u}\mathopen {}\left (\boldsymbol {x},t\right )\mathclose {}+\boldsymbol {\nabla }\boldsymbol {u}\mathopen {}\left (\boldsymbol {x},t\right )\mathclose {}^{\mkern -1.5mu\mathsf {T}}\right )\mathclose {}$ is the rate of strain tensor, 
$\boldsymbol {u}\mathopen {}\left (\boldsymbol {x},t\right )\mathclose {}$ is the Eulerian fluid velocity, 
$p\mathopen {}\left (\boldsymbol {x},t\right )\mathclose {}$ is the pressure, 
$\boldsymbol {f}\mathopen {}\left (\boldsymbol {x},t\right )\mathclose {}$ is the external body force density acting on the fluid, 
$\mu _{n}$ and 
$\mu _{p}$ are the Newtonian solvent and polymeric contributions to the viscosity, respectively, 
$\lambda$ is the relaxation time of the fluid, 
$D$ is the diffusion coefficient of the stress and 
$\mathop {\mathbb {C}}\limits ^{\triangledown }\mathopen {}\left (\boldsymbol {x},t\right )\mathclose {}$ is the upper convected derivative defined by
\begin{equation} \mathop{\mathbb{C}}\limits^{\triangledown}\mathopen{}\left(\boldsymbol{x},t\right)\mathclose{} = \frac{\partial \mathbb{C}\mathopen{}\left(\boldsymbol{x},t\right)\mathclose{}}{\partial t} + \boldsymbol{u}\mathopen{}\left(\boldsymbol{x},t\right)\mathclose{}\boldsymbol{\cdot}\boldsymbol{\nabla}\mathbb{C}\mathopen{}\left(\boldsymbol{x},t\right)\mathclose{} - \mathopen{}\left(\mathbb{C}\mathopen{}\left(\boldsymbol{x},t\right)\mathclose{}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}\mathopen{}\left(\boldsymbol{x},t\right)\mathclose{}^{\mkern-1.5mu\mathsf{T}}+\boldsymbol{\nabla}\boldsymbol{u}\mathopen{}\left(\boldsymbol{x},t\right)\mathclose{}\boldsymbol{\cdot}\mathbb{C}\mathopen{}\left(\boldsymbol{x},t\right)\mathclose{}\right)\mathclose{}. \end{equation}
The parameter 
$\alpha$ governs the strength of the nonlinear anisotropic drag term, which we refer to as the nonlinear parameter in the remainder of this work. This parameter determines the degree to which the fluid experiences shear thinning. Large values of 
$\alpha$ correspond to an enhanced capacity for shear thinning. The Giesekus model reduces to the Oldroyd-B model in the limit of 
$\alpha \rightarrow 0$. In the continuum equations of motion, the conformation tensor 
$\mathbb {C}\mathopen {}\left (\boldsymbol {x},t\right )\mathclose {}$ should remain positive definite at all times during the computation. If, during the course of the simulation, the conformation tensor loses positive definiteness, we project the conformation tensor onto the nearest non-negative-definite tensor (Guy & Fogelson Reference Guy and Fogelson2008).
We note the use of the unsteady Stokes equations in (2.1) as opposed to the steady Stokes equations. In our simulations, we utilize adaptive mesh refinement around the flagellum, which is precisely where we want to capture greater accuracy. However, our discretization of the Stokes operator near coarse–fine interfaces limits the solver to non-zero Reynolds numbers (Gruninger et al. Reference Gruninger, Barrett, Fang, Gregory Forest and Griffith2024).
It is well established that the Oldroyd-B model exhibits instabilities near extensional points (Renardy & Thomases Reference Renardy and Thomases2021), such as those found around the flagellum motor. For simulations using the Oldroyd-B model, we use a small diffusion coefficient 
$D$ that is proportional to the grid spacing. While it is unclear the degree to which the Giesekus model exhibits similar instabilities, our numerical tests indicate that stress diffusion is not needed for the non-zero values of 
$\alpha$ used in this study.
The flagellar element is described by a version of Kirchhoff rod theory, in which the flagellum is modelled as a thin rod and stresses are applied on cross-sections along the rod (Lim & Peskin Reference Lim and Peskin2012). The configuration of the flagellum is described by the current physical configuration of its centreline, 
$\boldsymbol {\chi }\mathopen {}\left (s,t\right )\mathclose {}$, and an orthonormal set of unit director vectors attached to the centreline, 
$\{\boldsymbol {D}^1\mathopen {}\left (s,t\right )\mathclose {},\boldsymbol {D}^2\mathopen {}\left (s,t\right )\mathclose {},\boldsymbol {D}^3\mathopen {}\left (s,t\right )\mathclose {}\}$, in which 
$s$ is a Lagrangian parameter with 
$0\leq s\leq L$, where 
$L$ is the length of the flagellum and 
$t$ is the time.
To describe the force and torque balance along the filament, let 
$\boldsymbol {F}^{rod}\mathopen {}\left (s,t\right )\mathclose {}$ and 
$\boldsymbol {N}^{rod}\mathopen {}\left (s,t\right )\mathclose {}$ be the force and moment, respectively, that are generated across a cross-section of the rod at point 
$s$ along the centreline. Let 
$\boldsymbol {F}\mathopen {}\left (s,t\right )\mathclose {}$ and 
$\boldsymbol {N}\mathopen {}\left (s,t\right )\mathclose {}$ be the applied force and torque densities of the fluid on the filament. Then the momentum and angular momentum balance equations are
\begin{gather} \boldsymbol{0} = \boldsymbol{F}\mathopen{}\left(s,t\right)\mathclose{} + \frac{\partial \boldsymbol{F}^{rod}\mathopen{}\left(s,t\right)\mathclose{}}{\partial s} , \end{gather}
\begin{gather}\boldsymbol{0} = \boldsymbol{N}\mathopen{}\left(s,t\right)\mathclose{} + \frac{\partial \boldsymbol{N}^{rod}\mathopen{}\left(s,t\right)\mathclose{}}{\partial s} + \mathopen{}\left(\frac{\partial \boldsymbol{\chi}\mathopen{}\left(s,t\right)\mathclose{}}{\partial s}\times \boldsymbol{F}^{rod}\mathopen{}\left(s,t\right)\mathclose{}\right)\mathclose{}. \end{gather}
The constitutive relations for the force and moment are derived from a variational argument using the elastic energy potential. In particular, we use a bistable energy formulation that allows for two different stable helices (Ko et al. Reference Ko, Lim, Lee, Kim, Berg and Peskin2017). We expand 
$\boldsymbol {F}^{rod}\mathopen {}\left (s,t\right )\mathclose {}$ and 
$\boldsymbol {N}^{rod}\mathopen {}\left (s,t\right )\mathclose {}$ in the basis of the local director vectors:
\begin{equation} \boldsymbol{F}^{rod}\mathopen{}\left(s,t\right)\mathclose{} = \sum_{i = 1}^3 F_i^{rod}\boldsymbol{D}^i\mathopen{}\left(s,t\right)\mathclose{} \quad\text{and}\quad \boldsymbol{N}^{rod}\mathopen{}\left(s,t\right)\mathclose{} = \sum_{i = 1}^3 N_i^{rod}\boldsymbol{D}^i\mathopen{}\left(s,t\right)\mathclose{}. \end{equation}Following Ko et al. (Reference Ko, Lim, Lee, Kim, Berg and Peskin2017), we utilize a bistable energy formulation that allows for two different stable helices:
\begin{gather} E = E_{bend} + E_{twist} + E_{shear} + E_{shear}, \end{gather}
\begin{gather}E_{bend} = \frac{1}{2}\int_0^L \mathopen{}\left(a_1\left(\varOmega_1\mathopen{}\left(s,t\right)\mathclose{}-\kappa_1\right)^2 + a_2\left(\varOmega_2\mathopen{}\left(s,t\right)\mathclose{}-\kappa_2\right)^2\right)\mathclose{}\,\text{d}s, \end{gather}
\begin{gather}E_{twist} = \int_0^L \mathopen{}\left(\frac{a_3}{4}\left(\varOmega_3\mathopen{}\left(s,t\right)\mathclose{}-\tau_1\right)^2\left(\varOmega_3\mathopen{}\left(s,t\right)\mathclose{}-\tau_2\right)^2 + \frac{\gamma^2}{2}\left(\frac{\partial\varOmega_3\mathopen{}\left(s,t\right)\mathclose{}}{\partial s}\right)^2\right)\mathclose{}\,\text{d}s, \end{gather}
\begin{gather}E_{shear} = \frac{1}{2}\int_0^L\mathopen{}\left(b_1\left(\frac{\partial\boldsymbol{\chi}\mathopen{}\left(s,t\right)\mathclose{}}{\partial s}\boldsymbol{\cdot}\boldsymbol{D}^1\mathopen{}\left(s,t\right)\mathclose{}\right)^2 + b_2\left(\frac{\partial\boldsymbol{\chi}\mathopen{}\left(s,t\right)\mathclose{}}{\partial s}\boldsymbol{\cdot}\boldsymbol{D}^2\mathopen{}\left(s,t\right)\mathclose{}\right)^2\right)\mathclose{}\,\text{d}s, \end{gather}
\begin{gather}E_{stretch} = \frac{1}{2}\int_0^L b_3\left(\frac{\partial\boldsymbol{\chi}\mathopen{}\left(s,t\right)\mathclose{}}{\partial s}\boldsymbol{\cdot}\boldsymbol{D}^3\mathopen{}\left(s,t\right)\mathclose{}-1\right)^2\,\text{d}s, \end{gather}which results in the forces
\begin{gather} F_i^{rod} = b_i\mathopen{}\left(\frac{\partial \boldsymbol{\chi}\mathopen{}\left(s,t\right)\mathclose{}}{\partial s}\boldsymbol{\cdot}\boldsymbol{D}^i\mathopen{}\left(s,t\right)\mathclose{} - \delta_{3i}\right)\mathclose{},\quad i=1,2,3, \end{gather}
\begin{gather}N_1^{rod} = a_1\mathopen{}\left(\varOmega_1\mathopen{}\left(s,t\right)\mathclose{} - \kappa_1\right)\mathclose{}, \end{gather}
\begin{gather}N_2^{rod} = a_2\mathopen{}\left(\varOmega_2\mathopen{}\left(s,t\right)\mathclose{} - \kappa_2\right)\mathclose{}, \end{gather}
\begin{gather}N_3^{rod} = a_3\mathopen{}\left(\varOmega_3\mathopen{}\left(s,t\right)\mathclose{} - \tau_1\right)\mathclose{}\mathopen{}\left(\varOmega_3\mathopen{}\left(s,t\right)\mathclose{} - \tau_2\right)\mathclose{}\mathopen{}\left(\varOmega_3\mathopen{}\left(s,t\right)\mathclose{} - \frac{\tau_1+\tau_2}{2}\right)\mathclose{} - \gamma^2\frac{\partial^2\varOmega_3\mathopen{}\left(s,t\right)\mathclose{}}{\partial s^2}, \end{gather}
in which 
$\delta _{3i}$ is the Kronecker delta function and 
$\gamma$ is the twist-gradient coefficient. The coefficients 
$a_1$ and 
$a_2$ are the bending moduli and 
$a_3$ is the twist modulus. In standard Kirchhoff rod theory, 
$\boldsymbol {D}^3$ is constrained to be the tangential vector of the structure. In the present work, (2.16) with 
$i = 3$ instead provides a penalty force that approximately enforces this constraint (Lim et al. Reference Lim, Ferent, Wang and Peskin2008). The coefficients 
$b_1$ and 
$b_2$ are the shear moduli and 
$b_3$ is the stretching modulus. The strain twist vector 
$\mathopen {}\left (\varOmega _1, \varOmega _2, \varOmega _3\right )\mathclose {}$, in which 
$\varOmega _i = {\partial \boldsymbol {D}^j\mathopen {}\left (s,t\right )\mathclose {}}/{\partial s}\boldsymbol {\cdot }\boldsymbol {D}^k\mathopen {}\left (s,t\right )\mathclose {}$ and 
$(i,j,k)$ is a cyclic permutation of 
$(1,2,3)$, determines the geometric properties of the helical flagellum. The parameter 
$\kappa = \sqrt {\kappa _1^2+\kappa _2^2}$ is the intrinsic curvature, and 
$\tau _1$ and 
$\tau _2$ are the intrinsic twists, which we assume are constant throughout the simulation. Given the curvature 
$\hat {\kappa } = \sqrt {\varOmega _1^2 + \varOmega _2^2}$ and twist 
$\hat {\tau } = \varOmega _3$ of the flagellum, the resulting helix has the radius 
$R$ and pitch 
$P$ with
\begin{equation} R = \frac{\hat{\kappa}}{\hat{\kappa}^2+\hat{\tau}^2} \quad\text{and}\quad P = \frac{2{\rm \pi}\hat{\tau}}{\hat{\kappa}^2+\hat{\tau}^2}. \end{equation}
Although we do not denote them as such, the values of 
$\kappa _1$, 
$\kappa _2$, 
$\tau _1$ and 
$\tau _2$ need not be constant along the flagellum, which we exploit to set different properties to the flexible hook that separates the motor from the flagellar filament. Although this work does not consider transition between helical shapes, we retain the bistable energy formulation used in by Ko et al. (Reference Ko, Lim, Lee, Kim, Berg and Peskin2017) so that future studies can study the effect of polymorphic transition in viscoelastic fluids.
The force density 
$\boldsymbol {f}\mathopen {}\left (\boldsymbol {x},t\right )\mathclose {}$ is generated by the deformation of the rotating elastic flagellum:
\begin{equation} \boldsymbol{f}\mathopen{}\left(\boldsymbol{x},t\right)\mathclose{} = \int_0^L -\boldsymbol{F}\mathopen{}\left(s,t\right)\mathclose{}\,\delta_w\mathopen{}\left(\boldsymbol{x}-\boldsymbol{\chi}\mathopen{}\left(s,t\right)\mathclose{}\right)\mathclose{}\,\text{d}s + \frac{1}{2}\boldsymbol{\nabla}\times\int_0^L -\boldsymbol{N}\mathopen{}\left(s,t\right)\mathclose{}\,\delta_w\mathopen{}\left(\boldsymbol{x}-\boldsymbol{\chi}\mathopen{}\left(s,t\right)\mathclose{}\right)\mathclose{}\,\text{d}s, \end{equation}
in which 
$\boldsymbol {F}\mathopen {}\left (s,t\right )\mathclose {}$ and 
$\boldsymbol {N}\mathopen {}\left (s,t\right )\mathclose {}$ are the force and torque applied by the fluid on the flagellum and 
$\delta _w(\boldsymbol {x})$ is a smooth, compactly supported kernel function that mediates coupling between the Lagrangian and Eulerian variables. The linear and angular velocities at the filament are computed by interpolating the Eulerian velocity onto the filament:
\begin{gather} \frac{\partial \boldsymbol{\chi}\mathopen{}\left(s,t\right)\mathclose{}}{\partial t} = \int_{\mathcal{B}} \boldsymbol{u}\mathopen{}\left(\boldsymbol{x},t\right)\mathclose{}\,\delta_w\mathopen{}\left(\boldsymbol{\chi}\mathopen{}\left(s,t\right)\mathclose{}-\boldsymbol{x}\right)\mathclose{}\,\text{d}\kern0.06em \boldsymbol{x}, \end{gather}
\begin{gather}\frac{\partial \boldsymbol{D}^i\mathopen{}\left(s,t\right)\mathclose{}}{\partial t} = \frac{1}{2}\int_{\mathcal{B}} \boldsymbol{\nabla}\times\boldsymbol{u}\mathopen{}\left(\boldsymbol{x},t\right)\mathclose{}\,\delta_w\mathopen{}\left(\boldsymbol{\chi}\mathopen{}\left(s,t\right)\mathclose{}-\boldsymbol{x}\right)\mathclose{}\,\text{d}\kern0.06em \boldsymbol{x},\quad \text{for } i=1,2,3, \end{gather}
in which 
$\mathcal {B}$ is the fixed Eulerian domain. In this work we use a delta function based on the three-point B-spline kernel (Lee & Griffith Reference Lee and Griffith2022). We remark that the smooth kernel function 
$\delta _w\mathopen {}\left (\boldsymbol {x}\right )\mathclose {}$ appears in both the continuum equations as well as the discretized equations. This is in contrast to traditional IB formulations, in which a smooth kernel function like 
$\delta _w\mathopen {}\left (\boldsymbol {x}\right )\mathclose {}$ appears in the discrete equations of motion but not in the continuum equations (Lim et al. Reference Lim, Ferent, Wang and Peskin2008). In conventional IB methods, the width of the smoothed delta is a numerical parameter proportional to the Eulerian mesh width that controls the accuracy of the regularized approximation to the singular delta function. In this model, the width of the delta function is a physical parameter of the model which can be viewed as controlling the effective thickness of the rod. Kernels with larger width yields rods with larger effective thickness.
The initial radius of the helical filament is
\begin{equation} R\mathopen{}\left(s\right)\mathclose{} = \left\{\begin{array}{ll} 0, & 0\leq s\leq L_h, \\ R_0\mathopen{}\left(1-e^{{-}c\mathopen{}\left(s-L_h\right)\mathclose{}^2}\right)\mathclose{}, & L_h\leq s\leq L_h+L_f,\end{array}\right. \end{equation}
in which 
$c = 2\,\mathrm {\mu }{\rm m}^{-2}$. The helical filament is then initialized as
\begin{equation} \boldsymbol{X}\mathopen{}\left(s,0\right)\mathclose{} = \left[R(s)\cos(2{\rm \pi} s / \bar{\omega}), R(s)\sin(2{\rm \pi} s / \bar{\omega}), s\right]^{\mkern-1.5mu\mathsf{T}} \end{equation}
in which 
$\bar {\omega } = 0.468\,\mathrm {\mu }{\rm m}$ is the wavenumber of the helix. The helical radius is 
$0\,\mathrm {\mu }{\rm m}$ for the hook and gradually increases to a radius of 
$R_0$. The vector 
$\boldsymbol {D}^3\mathopen {}\left (s,t\right )\mathclose {}$ is initially set as the unit tangent vector to the flagellum, and 
$\boldsymbol {D}^1\mathopen {}\left (s,t\right )\mathclose {}$ and 
$\boldsymbol {D}^2\mathopen {}\left (s,t\right )\mathclose {}$ are the normal and binormal unit vectors. The flagellum is driven by a rotary motor. We fix the point at 
$s = 0$ in space and specify the rotation of the director vectors as
\begin{gather} \boldsymbol{D}^1\mathopen{}\left(0,t\right)\mathclose{} = \mathopen{}\left(\cos\mathopen{}\left(2{\rm \pi}\omega t\right)\mathclose{}, -\sin\mathopen{}\left(2{\rm \pi}\omega t\right)\mathclose{}, 0\right)\mathclose{}, \end{gather}
\begin{gather}\boldsymbol{D}^2\mathopen{}\left(0,t\right)\mathclose{} = \mathopen{}\left(\sin\mathopen{}\left(2{\rm \pi}\omega t\right)\mathclose{}, \cos\mathopen{}\left(2{\rm \pi}\omega t\right)\mathclose{}, 0\right)\mathclose{}, \end{gather}
\begin{gather}\boldsymbol{D}^3\mathopen{}\left(0,t\right)\mathclose{} = \mathopen{}\left(0,0,1\right)\mathclose{}, \end{gather}
in which 
$\omega$ is the specified rotation rate. The sign of 
$\omega$ determines the direction of the rotation. The rotation of the motor generates a torque that is then transmitted to the flagellum through the compliant hook. The length of the hook for an E. coli bacterium ranges from 
$50$ to 
$80\,{\rm nm}$. Here, we specify the hook's length to be 
$L_{h} = 80\,{\rm nm}$. To make the hook flexible, we specify its bending modulus to be two orders of magnitude smaller than that of the filament (Son, Guasto & Stocker Reference Son, Guasto and Stocker2013; Jabbarzadeh & Fu Reference Jabbarzadeh and Fu2018). We also assume the hook to be intrinsically straight, so that 
$\tau = \kappa = 0\,\mathrm {\mu }{\rm m}^{-1}$.
The model is implemented in the open source software IBAMR (Griffith Reference Griffith2023), which is an MPI parallelized implementation of the IB method with adaptive mesh refinement (AMR). The software uses SAMRAI (Hornung & Kohn Reference Hornung and Kohn2002; Hornung, Wissink & Kohn Reference Hornung, Wissink and Kohn2006) for AMR and PETSc (Balay et al. Reference Balay, Gropp, McInnes and Smith1997, Reference Balay2015, Reference Balay2019) for iterative linear solvers along with custom preconditioners for AMR discretizations of the incompressible Navier–Stokes equations.
The equations are discretized on a staggered grid in which the normal components of velocities are stored on cell sides and the pressures and components of the conformation tensor are stored on cell centres (Harlow & Welch Reference Harlow and Welch1965). We use second-order finite differences to discretize the Laplacian, velocity gradients and stress divergence (Barrett Reference Barrett2019). The advective derivative is discretized with a second-order wave propagation algorithm (Ketcheson, Parsani & LeVeque Reference Ketcheson, Parsani and LeVeque2013).
We discretize in time using the implicit trapezoidal rule for the viscous terms and an explicit predictor–corrector method for the conformation tensor. The resulting linear system is solved using a GMRES algorithm with a projection method as a preconditioner (Griffith Reference Griffith2009). The details of the spatial and temporal discretizations are discussed in previous work (Barrett Reference Barrett2019; Gruninger et al. Reference Gruninger, Barrett, Fang, Gregory Forest and Griffith2024).
The filament is placed in a computational domain 
$\mathcal {B}$ which is a cube of length 
$H = \frac {10}{3}L = 20\,\mathrm {\mu }{\rm m}$. At the physical boundaries, we specify zero velocity and use linear extrapolation for the stress. The computational domain is discretized such that 
$N = 512$ points in each direction would fill the finest level of the AMR grid. We discretize the Lagrangian structure so that the structure contains approximately one point per Eulerian grid cell.
3. Results
The flagellum is placed in a viscoelastic fluid, and the motor rotates at a prescribed rate of 
$\omega = 2{\rm \pi} 100\,{\rm rad}\,{\rm s}^{-1}$. We can classify the elasticity of the fluid by using the non-dimensional Deborah number, which is the ratio of the polymeric relaxation time 
$\lambda$ to a characteristic time scale of the structure 
$T_{s}$. We define 
$T_{s}$ to be the time for the motor to complete one turn, so that 
$T_{s} = ({2{\rm \pi} }/{\omega })\,{\rm s} = 0.01\,{\rm s}$. This gives a Deborah number of 
${\textit {De}} = {\lambda }/({2{\rm \pi} /\omega })$. The total viscosity of the fluid is given by 
$\mu = \mu _{n} + \mu _{p}$. The viscosities for the viscoelastic fluids are 
$\mu _{p} = \frac {1}{3}\times 10^{-6}\,{\rm g}\,(\mathrm {\mu }{\rm m}\,{\rm s})^{-1}$ and 
$\mu _{n} = \frac {2}{3}\times 10^{-6}\,{\rm g}\,(\mathrm {\mu }{\rm m}\,{\rm s})^{-1}$, giving a viscosity ratio of 
$\beta = {\mu _{p}}/({\mu _{n} + \mu _{p}}) = \frac {1}{3}$. We also simulate a flagellum in two Newtonian fluids with 
$\mu _{p} = 0$: one with 
$\mu _{n} = 1\times 10^{-6}\,{\rm g}\,(\mathrm {\mu }{\rm m}\,{\rm s})^{-1}$, which we refer to as 
$N^{high}$, and one with 
$\mu _{n} = \frac {2}{3}\times 10^{-6}\,{\rm g}\,(\mathrm {\mu }{\rm m}\,{\rm s})^{-1}$, which we refer to as 
$N^{low}$. The 
$N^{high}$ fluid has viscosity matching the total viscosity of the viscoelastic fluids, whereas the 
$N^{low}$ fluid has viscosity matching only the Newtonian solvent viscosity of the viscoelastic fluids. We define the Reynolds number as 
${\textit {Re}} = {\rho U R_0}/{\mu }$, in which 
$U = R_0 / T_{s}$ is the characteristic velocity of the system. This yields a Reynolds number of approximately 
${\textit {Re}} \approx 4.27\times 10^{-6}$ for 
$N^{high}$ and 
${\textit {Re}} \approx 6.41\times 10^{-6}$ for 
$N^{low}$. Note that the Reynolds number for the Giesekus fluid is not well defined, as the effective viscosity of the fluid is shear-rate dependent. If we use the zero shear-rate viscosity, the Reynolds number for all Giesekus fluids is 
${\textit {Re}}\approx 4.27\times 10^{-6}$, that of the 
$N^{high}$ fluid. The remaining parameters are given in table 1. For the simulations that follow, we fix the rotation rate and vary the Deborah number 
${\textit {De}}$ by the parameter 
$\lambda$ and the nonlinear parameter 
$\alpha$ of the polymeric fluid model. Non-zero values of 
$\alpha$ tune the quadratic nonlinearity of extra stress in the Giesekus model, which controls the degree of shear thinning, and which we demonstrate has profound influences, particularly on the pumping efficiency of the flagellum. In the results that follow, a periodic steady state is achieved after several rotations of the motor. For plots in which average values are shown, we compute the averages after a periodic steady state has been achieved. For fluids with high Deborah numbers and low 
$\alpha$ values, more cycles are required to achieve a periodic steady state, as the fluid's elastic energy storage capacity is increased. Figure 1(a) shows the shape of superimposed flagella rotating in fluids of varying 
$\alpha$ values. Figure 1(b) shows the trace of the conformation tensor around the flagellum. The trace of the conformation tensor is proportional to the elastic strain energy density of the polymeric material (Li et al. Reference Li, Qin, Gopinath, Arratia, Thomases and Guy2017). We observe trace values that are approximately 25 times larger for the Oldroyd-B fluid (
$\alpha = 0.0$) than for the Giesekus fluid with 
$\alpha = 0.3$. In addition, we observe that large regions of substantial viscoelastic stress exist through the entire path the flagellum traces for 
$\alpha = 0.0$. For 
$\alpha = 0.3$, the viscoelastic stresses are localized only near the current location of the flagellum.
Table of physical and computational parameters for both the model flagellum and the fluid.
Panel (a) shows three superimposed flagella from the side and top perspectives. The black curve corresponds to a flagellum in the fluid 
$N^{high}$, the green curve corresponds to a flagellum in a viscoelastic fluid with 
${\textit {De}} = 1.0$ and 
$\alpha = 0.3$ and the red curve corresponds to a flagellum in a viscoelastic fluid with 
${\textit {De}}=1.0$ and 
$\alpha =0.0$. The shapes between the green and black curves are closer than the shape of the red curve. We also show the motor triads at the beginning of the flagellum. The physical location is fixed in place, and the rotation rate of the triads is prescribed. Panel (b) shows the trace of the conformation tensor 
$\mathbb {C}$ along a slice parallel to the flagellum with 
${\textit {De}} = 1.0$ and 
$\alpha = 0.0$ (b i) and 
$\alpha = 0.3$ (b ii). Note that, for 
$\alpha = 0.0$, we observe trace values roughly 25 times larger than for 
$\alpha = 0.3$. Additionally, large regions of stress are found throughout the entire path the flagellum traces for 
$\alpha = 0.0$, while for 
$\alpha = 0.3$, the stress quickly dissipates away from the flagellum. Brighter and darker regions of the flagellum highlight the flagellum being in front of or behind the plane.
3.1. Domain dependence
Because the model is at very low Reynolds number, the domain can have a profound impact on the results of the simulation. We test the effect that the size of the domain has on the flow by systematically increasing the size of the domain while keeping the grid spacing fixed. We measure the total flux through a circular disk 
$\mathcal {D}$ above the helix, as will be described in a later section. The flux is shown in figure 2 for a Newtonian fluid. We observe an increase in flow rates as we increase the size of the domain. We note that due to the hyperbolic nature of the viscoelastic constitutive law, there is little transport of viscoelastic stress towards the boundary, so we expect the boundary effects from viscoelasticity to be comparatively minor. The simulations that follow use a domain size of 
$L = 10\,\mathrm {\mu }{\rm m}$.
The average flux along with maximum and minimum values through the region 
$\mathcal {D}$ as a function of the domain size, normalized by the length of the flagellum 
$L$. We observe an increase in flow rates as the domain size increases. The error bars signify the maximum and minimum flow rates.
3.2. Flagellum shape
We assess the shape of the flagellum by computing the radius and pitch of the helix using the strain twist vector in (2.20a,b). The twist of the helix is computed as 
$\hat {\tau } = \varOmega _3$, and the curvature is computed as 
$\hat {\kappa } = \sqrt {\varOmega _1^2 + \varOmega _2^2}$. Figure 3 shows the time-average, maximum and minimum values of the radius and pitch along the middle (
$s = L/2$) and end (
$s = L$) of the flagellum as we vary both the Deborah number 
${\textit {De}}$ and the nonlinear parameter 
$\alpha$. The reported radii and pitches are normalized with respect to the time-average radius and pitch for the 
$N^{high}$ fluid. Across various Deborah numbers, we find that fluids characterized by values of 
$\alpha$, which indicate a reduced shear-thinning capacity, tend to result in smaller measurements for both the pitch and radius of the flagellum. Additionally, we have found that significant variations in the shape of the flagella primarily result from lower 
$\alpha$ values and higher Deborah numbers. Therefore, we infer that major changes in the flagellar helix shape stem from increased elastic responses within the fluid. Conversely, the presence of shear thinning seems to mitigate these effects.
The time-averaged and maximum and minimum radii (a,b) and pitch (c,d) of the flagellum along the middle (a,c) and end (b,d), as calculated by (2.20a,b). The radius and pitch have been normalized by those in 
$N^{high}$. For a small value of the nonlinear parameter 
$\alpha = 0.01$, we observe initial increases in pitch and radius as we increase the Deborah number 
${\textit {De}}$, followed by larger decreases, eventually becoming smaller than that of 
$N^{high}$. For all other fluids, we observe increases in pitch and radius as we increase 
${\textit {De}}$ approaching those of the 
$N^{low}$ fluid. Note that, in all cases, the changes are no more than 10 % of those in 
$N^{high}$ and frequently are only a few percentage points.
Overall, we find that the shape of the helix remains relatively consistent in the viscoelastic fluids when compared with Newtonian fluids. The shape changes for all parameters tested are less than 10 % different than in a Newtonian fluid, and typically only vary within a few percentage points. This is in contrast to many undulatory swimmers, which swim with modified shapes in viscoelastic fluids (Fu, Wolgemuth & Powers Reference Fu, Wolgemuth and Powers2008; Li et al. Reference Li, Qin, Gopinath, Arratia, Thomases and Guy2017). While the shape changes only slightly in viscoelastic fluids, shape changes can have drastic effects on flagellar bundling (Lee et al. Reference Lee, Kim, Griffith and Lim2018), which are vital to bacterial locomotion. This effect remains an important area of future study. We also note that the apparent large shape differences in figure 1(a) are a result of the varying pitch. If we project the helix onto the 
$y\unicode{x2013}z$ plane, we obtain a sinusoidal curve, whose amplitude is proportional to the radius of the helix and frequency is proportional to the pitch. Slight changes in the pitch can therefore result in a helix appearing out of phase at different time points.
3.3. Pumping efficiency
To characterize pumping efficiency, we measure the volumetric flow rate through a disk 
$\mathcal {D}$ of radius 
$R = 1\,\mathrm {\mu }{\rm m}$ at the plane 
$z = 5.5\,\mathrm {\mu }{\rm m}$. Although the flagellum's length is 
$6\,\mathrm {\mu }{\rm m}$, the flagellum's coiled length is approximately 
$5.25\,\mathrm {\mu }{\rm m}$. Consequently, the tip of the flagellum never passes through the disk 
$\mathcal {D}$. The volumetric flow rate is computed as
\begin{equation} Q = \int_{\mathcal{D}} \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{n} \,\text{d}A, \end{equation}
in which 
$\boldsymbol {n}$ is the upward unit normal of the disk 
$\mathcal {D}$. We note that the choice of the size of the disk 
$\mathcal {D}$ determines the flux measure, because in the limit as the radius of the disk tends to infinity, the volumetric flow rate 
$Q$ tends to zero as we encapsulate more regurgitation flow. Figure 4 shows the region through which the flow rate is calculated. Because we specify the rotation rate of the motor, we can also calculate the torque required to turn the motor at the specified frequency. We calculate the total torque acting on the hook from the motor using (2.9).
The disk 
$\mathcal {D}$ through which the volumetric flow rate is measured is shown in red. The velocity vectors shown on the region 
$\mathcal {D}$ correspond to the fluid with 
$\alpha = 0.3$ and 
${\textit {De}} = 1.0$ and have been normalized by the characteristic velocity 
$U = R_0/T_{s} \approx 126\,\mathrm {\mu }{\rm m}\,{\rm s}^{-1}$.
Figure 5 shows the flow rate 
$Q$ through the region 
$\mathcal {D}$ as well as the total torque acting on the hook for two fixed Deborah numbers as we vary the nonlinear parameter 
$\alpha$. We find that a steady state is quickly reached after the first few rotations of the motor. The exception is for larger Deborah numbers and small 
$\alpha$ values, which require additional rotations to reach a steady state. This is indicative of the instabilities inherent in the Oldroyd-B model for very high Deborah numbers (Thomases & Guy Reference Thomases and Guy2014); although, because of the inclusion of stress diffusion, we do anticipate that steady state values will eventually be achieved for all Deborah numbers. Figure 6 shows the time-averaged, maximum and minimum flow rates and time-averaged torque for each fluid, normalized with respect to the time-averaged values obtained for the 
$N^{high}$ fluid. For a small nonlinear parameter value of 
$\alpha = 0.01$, we observe a non-monotonic required torque to spin the motor as a function of the Deborah number. The required torque initially decreases, achieving a minimum at approximately 
${\textit {De}} = 0.3$, followed by increases as we increase the Deborah number further. For Deborah number 
${\textit {De}} = 2.0$, the torque required to spin the motor is nearly 1.4 times that of the required torque in the 
$N^{high}$ fluid. For all other values of 
$\alpha$, we observe slow decreases in the required torque as we increase the Deborah number, with the lowest required torque found with 
$\alpha = 0.3$. For all fluids tested, the required torques were higher than that of flagellum in the 
$N^{low}$ fluid. This suggests that, for large enough 
$\alpha$ values, the shear-thinning behaviour of the viscoelastic fluid reduces the impact of its elasticity, lowering the effective viscosity near the motor below that of the 
$N^{high}$ fluid and towards that of the 
$N^{low}$ fluid.
Panels (a,c) show the torque generated by the motor for varying nonlinear parameter 
$\alpha$ and fixed Deborah number 
${\textit {De}}$ normalized by that in the 
$N^{high}$ fluid. We observe a decrease in required torque as we increase the nonlinear parameter 
$\alpha$. Panels (b,d) show the volumetric flow rate through the fixed disk 
$\mathcal {D}$ normalized by that in the 
$N^{high}$ fluid. Panels (a,b) use a fixed Deborah number of 
${\textit {De}} = 0.5$. Panels (c,d) use 
${\textit {De}} = 1.0$. We observe substantial decreases in the flow rates as the nonlinear parameter increases.
The time-averaged torque and flow rate for various viscoelastic fluids normalized by the torque and flow rate for the 
$N^{high}$ fluid. Also shown by the error bars is the maximum and minimum flow rates across a motor rotation. Error bars are not shown for the torque, as the torque is relatively constant in time, see figure 5. We observe that, for small 
$\alpha$ values, the torque and flow rates increase as we increase the Deborah number. For larger 
$\alpha$ values, the torque decreases as we increase the Deborah number, and the flow rates marginally increase.
The flow rate for every fluid considered is larger than that of the 
$N^{high}$ fluid. Generally, we observe increases in flow rates as we increase the Deborah number or decrease the nonlinear parameter 
$\alpha$ with the largest flow rates observed corresponding to the largest 
$De$ and smallest 
$\alpha$ values tested. Similar to the shape of the flagellum, the flow rate varies across a complete rotation of the motor. Further, congruent with our findings regarding the shape variation of the flagellum, we find that flow variations increase as the Deborah number is increased and the nonlinear parameter 
$\alpha$ is decreased.
To assess the performance of the motor, we compute its instantaneous efficiency
\begin{equation} E_i(t) = \frac{Q_i(t) / \tau_i(t)}{Q_{N^{high}}(t)/\tau_{N^{high}}(t)}, \end{equation}
in which 
$\tau _{N^{high}}(t)$ and 
$Q_{N^{high}}(t)$ are the torque and flow rate from the 
$N^{high}$ fluid and 
$\tau _i(t)$ and 
$Q_i(t)$ are the torque and flow rate from fluid 
$i$. We then average the efficiency 
$E_i(t)$ over the last five rotations of the motor. Efficiency values greater than one imply that the motor in fluid 
$i$ is able to pump fluid more efficiently than the motor in fluid 
$N^{high}$. Figure 7 shows the efficiency for all fluids tested. We perform linear interpolation of the efficiency value for combinations of 
${\textit {De}}$ and 
$\alpha$ that were not simulated. In all simulations, the efficiency of the motor is greater than the efficiency in the 
$N^{high}$ fluid. The solid black curve in figure 7 shows the contour on which the efficiency is equal to that of 
$N^{low}$. We observe that, for small 
$\alpha$ values and large 
${\textit {De}}$ values, the efficiency of the motor is greater than that for 
$N^{low}$.
The efficiency of the motor as computed by (3.2) as we vary the nonlinear parameter 
$\alpha$ and the Deborah number 
${\textit {De}}$. The black circles are data points from the simulation, and linear interpolation is used to compute the remaining efficiencies. The black curve is the contour on which the efficiency is equal to that of 
$N^{low}$. For every simulation, we observe a greater efficiency than in the 
$N^{high}$ fluid. For small 
$\alpha$ and large 
${\textit {De}}$ values, we observe a greater efficiency than in the 
$N^{low}$ fluid.
Our findings reveal a nuanced interplay between the viscoelastic fluid's elasticity and shear-thinning properties and their impact on motor efficiency. Increased elasticity necessitates a higher torque to spin the motor, yet it boosts the flow rate to a greater extent, resulting in enhanced motor efficiencies at higher Deborah numbers. Conversely, while shear thinning reduces the required torque to spin the motor, it substantially reduces the flow rate, thereby diminishing motor efficiency at larger values of 
$\alpha$. Additionally, the enhancement in motor efficiency driven by increasing Deborah number values is more pronounced at lower values of 
$\alpha$.
These results appear to be inconsistent with recent results on locomotion in complex fluids. In generalized Newtonian fluids that exhibit shear thinning with negligible elasticity, several groups (Gagnon et al. Reference Gagnon, Keim and Arratia2014; Gómez et al. Reference Gómez, Godínez, Lauga and Zenit2017; Hyakutake et al. Reference Hyakutake, Sato and Sugita2019; Qu & Breuer Reference Qu and Breuer2020) have demonstrated that organisms swim faster than in comparable Newtonian fluids. The common explanation is that swimmers benefit from a soft confinement effect. The fluid exhibits a lower viscosity in the immediate vicinity of the swimmer, allowing the swimmer to ‘push off’ the more viscous fluid layer. For fluids with significant elastic properties, a prevailing theory is that swimmers gain a speed boost from the stored energy in the fluid (Teran et al. Reference Teran, Fauci and Shelley2010; Thomases & Guy Reference Thomases and Guy2014; Li et al. Reference Li, Qin, Gopinath, Arratia, Thomases and Guy2017). For Deborah numbers 
${\textit {De}} \geq 1$, the relaxation time of the fluid exceeds the rotation period of the flagellum, meaning the flagellum continues to return to a volume of stored elastic stress. This provides the flagellum with elastic resistance allowing the flagellum to swim faster, analogous to a pusher.
For the viscoelastic model used here, however, these two stories are at odds with each other. Shear thinning is induced in a neighbourhood around the flagella, providing less viscous resistance to the flagellum's motion. Increased shear thinning requires less torque to spin the motor and therefore weaker elastic resistance for pushing. The consequence of this can be observed in figure 7, which shows that the efficiency of the motor decreases as we increase the shear-thinning capacity of the fluid. The elastic energy of the fluid is given by the trace of the conformation tensor (Li et al. Reference Li, Qin, Gopinath, Arratia, Thomases and Guy2017). Figure 8 shows a plot of the trace of the conformation tensor along the middle of the flagellum as we vary 
$\alpha$ and 
${\textit {De}}$. As expected, for 
${\textit {De}} \geq 1.0$, the elastic energy of the fluid remains high throughout the path that the flagellum traces out in the fluid. However, the magnitude of the trace drastically decreases as we increase 
$\alpha$. For large 
$\alpha$ values, the flagellum is returning to areas with lower elastic energies.
The trace of the conformation tensor 
$\mathbb {C}\mathopen {}\left (\boldsymbol {x},t\right )\mathclose {}$ as we vary the Deborah number and the nonlinear parameter. The trace is shown on a slice along the midpoint of the flagellum. We observe that as the Deborah number increases, the flagellum begins to return to regions of stored elastic energy, which enhances the pumping capacity of the flagellum. As we increase the nonlinear parameter, the magnitude of the stored stress sharply decreases, which reduces the effect of fluid elasticity on the flagellum. Note the colour bar uses a logarithmic scaling.
We also do not see a clear confinement effect around the flagellum. Figure 9 shows the magnitude of the velocity as we vary 
${\textit {De}}$ and 
$\alpha$. For a fixed small nonlinear parameter of 
$\alpha = 0.01$, we actually observe that the region that experiences non-zero flow increases in size as we increase the Deborah number: the opposite of a confinement effect. As we increase the shear-thinning capacity of the fluid, the size of this region decreases, approaching a size that is comparable to that of the 
$N^{high}$ and 
$N^{low}$ fluid. Because the shear thinning is limited to the polymeric stress, we do not observe a confinement effect overall, as the effect from the Newtonian solvent becomes dominant.
The velocity magnitude as we vary the Deborah number 
${\textit {De}}$ and the nonlinear parameter 
$\alpha$. The velocity magnitude is shown on a slice along the midpoint of the flagellum and has been normalized by the non-dimensional velocity 
$U = R_0 / T_{s}$. As we increase the Deborah number, we observe that the size of the region of non-zero velocities increases. However, as we increase 
$\alpha$, the size of the region decreases. Also shown are the velocity magnitudes for the Newtonian 
$N^{low}$ and 
$N^{high}$ fluids. As we increase 
$\alpha$, the size of the region approaches that of the Newtonian fluids, suggesting marginal confinement effects.
4. Discussion and concluding remarks
We have studied both the shape of the flagellum and the motor efficiency in viscoelastic fluids with varying capacities of shear thinning. We compared these flagella with those in Newtonian fluids, one with viscosity equalling the total of the polymeric and solvent viscosities 
$N^{high}$, and one with viscosity equalling that of only the solvent 
$N^{low}$. We have observed different helical shapes assumed by the flagellum, which for the range of parameters considered in this study, differed by less than 10 % when compared with a Newtonian fluid. While it is unclear how much effect this slight change has on the pumping efficiency of the flagellum, the bundling ability of multiple flagella could be significantly affected and is worth further investigation (Lee et al. Reference Lee, Kim, Griffith and Lim2018).
We note that the flagellum elasticity model utilized in this study is based on physical characterizations of the flagella of E. Coli (Goldstein et al. Reference Goldstein, Goriely, Huber and Wolgemuth2000; Darnton et al. Reference Darnton, Turner, Rojevsky and Berg2007; Ko et al. Reference Ko, Lim, Lee, Kim, Berg and Peskin2017), and this model has been successfully used to study flagellar bundling and locomotion in Newtonian fluids (Lim & Peskin Reference Lim and Peskin2012; Ko et al. Reference Ko, Lim, Lee, Kim, Berg and Peskin2017; Jabbarzadeh & Fu Reference Jabbarzadeh and Fu2018). The elasticity model has several time scales related to stretching, twisting and bending, and the interplay between the time scales of the flagellum and the elastic time scale of the fluid is not fully explored in this manuscript. The time scales of the flagellum are derived in the supplementary material available at https://doi.org/10.1017/jfm.2024.666. For undulatory swimmers, Thomases & Guy (Reference Thomases and Guy2017) found that swimmers can achieve greater speed boosts if their gait can adjust to the elasticity of the fluid. To the authors’ knowledge, a similar analysis has not been performed for helical swimmers. For the results discussed herein, when compared with the rotation period of the motor, we find the time scales of the flagellum are much faster than that of the motor, making the flagellum a stiff body. The time scales of the flagellum are also faster than that of the elasticity of the fluid, and the flagellum's shape is not substantially changed by fluid elasticity. However, shear thinning of the fluid reduces the elastic energy stored by the fluid, which complicates the conclusion reached by Thomases & Guy (Reference Thomases and Guy2017).
The importance of shear-thinning behaviour for speed enhancement in helical swimmers has been demonstrated previously. Gómez et al. (Reference Gómez, Godínez, Lauga and Zenit2017) and Demir et al. (Reference Demir, Lordi, Ding and Pak2020) both studied helical swimming in shear-thinning fluids with negligible viscoelasticity. They both demonstrated that shear thinning can greatly enhance swimming speed for helical swimmers, and they concluded that confinement effects are the primary cause of speed enhancement. Demir et al. further hypothesized that the shear-thinning effects dominate any viscoelastic effects, due to the magnitude of the swimming enhancement observed in their studies. Qu & Breuer (Reference Qu and Breuer2020), in their studies of helical swimming in viscoelastic fluids with shear thinning, also found that the magnitude of shear-thinning effects are dominant over the viscoelastic effects.
We find that the motor efficiency in all viscoelastic fluids is greater than the efficiency in the 
$N^{high}$ fluid. For fluids with large Deborah numbers and small values of the nonlinear parameter 
$\alpha$, the efficiency is greater than the 
$N^{low}$ fluid. As we increase 
$\alpha$, the efficiency of the motor decreases. This decrease in motor efficiency can be explained by the complex interplay between the elastic and shear-thinning responses of the viscoelastic fluid. Although greater elasticity enhances the storage of elastic energy in the fluid, facilitating higher swimming speeds and thus increased flow rates (Li et al. Reference Li, Qin, Gopinath, Arratia, Thomases and Guy2017, Reference Li, Lauga and Ardekani2021; Spagnolie et al. Reference Spagnolie, Liu and Powers2013), shear thinning has a contrasting effect. It diminishes the torque required for the flagella to rotate, which results in a decreased amount of elastic energy stored in the surrounding fluid and decreased flow rates. Further, in our results, we do not observe any confinement effect, and the velocity profile approaches that of a purely Newtonian fluid as 
$\alpha$ is increased. In fact, for small values of 
$\alpha$, we observe that the region of substantial flow around the flagellum increases in size as we increase the Deborah number. This finding is in contrast to Qu & Breuer (Reference Qu and Breuer2020), and this difference in results can be partially attributed to the low Deborah numbers utilized in Qu & Breuer's experiments. Our findings suggest that the Deborah number needs to be larger than unity to fully observe speed enhancements.
Our results suggest that too much shear thinning in the Giesekus model actually reduces the efficiency of the swimmer. We note that this is not in contradiction with other results. Instead, we claim that the mechanism for shear thinning is important. For the Giesekus model, the solvent behaves as a Newtonian fluid. If the viscoelastic contribution to the stress experiences excessive shear thinning, the viscoelastic stresses become subdominant to the Newtonian solvent's stresses, making the fluid appear Newtonian. Further, when making comparisons with Newtonian fluids, it is important to carefully select the appropriate Newtonian fluid with which to compare. The Reynolds number measures the relative importance of inertial and shear forces. The classical definition relies on a choice of viscosity, which for shear-thinning fluids is shear-rate dependent, which varies in both space and time in our context. In the results considered here, we compare with Newtonian fluids whose Reynolds numbers match with the zero shear-rate viscosity of the shear-thinning fluids. As discussed by Thompson & Oishi (Reference Thompson and Oishi2021), more appropriate comparisons would need to determine the characteristic viscosity near the flagellum. An important step to fully elucidating the mechanisms behind non-Newtonian swimming is to determine the detailed rheology of the fluids in which bacteria swim, and from this, determine an accurate model that fully captures the physical mechanisms of shear thinning and viscoelasticity. This can be a complicated undertaking, however, because biological fluids typically have multiple relaxation times as well as a multitude of other properties, such as polymer chain entanglements of mucus in the respiratory, gut and cervicovaginal tracts.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2024.666.
Acknowledgements
The authors thank R. Guy for helpful conversations about motility in viscoelastic fluids.
Funding
We acknowledge support from NSF awards DMS 1929298, DMS 1853591, DMS 1410873, DMS 1664645, OAC 1450327, OAC 1652541 and OAC 1931516 and NIH awards U01HL143336 and R01HL157631. M.G.F. acknowledges support by the Alfred P. Sloan Foundation award G-2021-14197. C.G. is grateful for support from the Department of Defense (DoD) through the National Defense Science and Engineering Graduate (NDSEG) Fellowship Program. S.L. acknowledges support from the Charles Phelps Taft Research Center, University of Cincinnati.
Declaration of interests
The authors report no conflict of interest.
Open access
