Flagellum pumping efficiency in shear-thinning viscoelastic fluids

Abstract

Microorganism motility often takes place within complex, viscoelastic fluid environments, e.g. sperm in cervicovaginal mucus and bacteria in biofilms. In such complex fluids, strains and stresses generated by the microorganism are stored and relax across a spectrum of length and time scales and the complex fluid can be driven out of its linear response regime. Phenomena not possible in viscous media thereby arise from feedback between the swimmer and the complex fluid, making swimming efficiency co-dependent on the propulsion mechanism and fluid properties. Here, we parameterize a flagellar motor and filament properties together with elastic relaxation and nonlinear shear-thinning properties of the fluid in a computational immersed boundary model. We then explore swimming efficiency, defined as a particular flow rate divided by the torque required to spin the motor, over this parameter space. Our findings indicate that motor efficiency (measured by the volumetric flow rate) can be boosted or degraded by relatively moderate or strong shear thinning of the viscoelastic environment.

1. Introduction

Microorganisms use a variety of techniques to swim in low Reynolds number environments (Pelczar, Gunsalus & Stanier 1960; Berg 2003). Organisms such as Escherichia coli use multiple rotating flagellar filaments that are connected to a rotary motor by a flexible short hook (Berg 2003). Spermatozoa have flagellar waveforms that are nearly planar (Simons, Fauci & Cortez 2015). Other organisms, such as Chlamydomonas reinhardtii, use an undulating motion to propel themselves forward (Lauga & Powers 2009). 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 1987). One such choice is the class of upper convected Maxwell-like constitutive models such as the Oldroyd-B model (Larson 1988; Beris & Edwards 1994; Morozov & Spagnolie 2015), 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 1988). 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. 1987).

The physics of locomotion of microorganisms has been extensively studied in Newtonian fluids (Lauga & Powers 2009), 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 (2021). For undulating swimmers, Lauga (2007) 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 2010; Thomases & Guy 2014). Fu, Wolgemuth & Powers (2009) 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 (2015) showed that the swimming speed can increase provided that kinematic waves travel in opposing directions along the sheet. Li et al. (2017) and Hyakutake, Sato & Sugita (2019) 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 2013; Gagnon, Keim & Arratia 2014), and these effects were deemed dominant over the hindrance from viscoelasticity, which is consistent with Hyakutake et al. In contrast, Dasgupta et al. (2013) 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. 2016; Li et al. 2021; Schuech, Cortez & Fauci 2022).

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. 2009; Li & Spagnolie 2015). However, for helical swimmers with a large pitch angle, the resulting swimming speed becomes nonlinear in the Deborah number. Liu, Powers & Breuer (2011) 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 (2013) 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. (2017) and Demir et al. (2020) 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 (2015) 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 (2020) 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. 2008). Goldstein et al. (2000) introduced a bistable energy formulation for Kirchhoff rod theory to permit two stable helical configurations. Darnton & Berg (2007) used Kirchhoff rod theory to estimate the bending rigidity of the flagellum by fitting a Kirchhoff rod model to experimental data. Lim & Peskin (2012) incorporated hydrodynamic interactions by creating a generalized immersed boundary (IB) method based on Kirchhoff rod theory. Ko et al. (2017) 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 (2012) and Griffith & Lim (2012) 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

(2.1)

(2.2)\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}\mathopen{}\left(\boldsymbol{x},t\right)\mathclose{} = 0, \end{gather}

(2.3)

(2.4)

(2.5)

(2.6)\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

(2.7)\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 2008).

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. 2024).

It is well established that the Oldroyd-B model exhibits instabilities near extensional points (Renardy & Thomases 2021), 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 2012). 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

(2.8)\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}

(2.9)\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. 2017). 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:

(2.10a,b)\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. (2017), we utilize a bistable energy formulation that allows for two different stable helices:

(2.11)\begin{gather} E = E_{bend} + E_{twist} + E_{shear} + E_{shear}, \end{gather}

(2.12)\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}

(2.13)\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}

(2.14)\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}

(2.15)\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

(2.16)\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}

(2.17)\begin{gather}N_1^{rod} = a_1\mathopen{}\left(\varOmega_1\mathopen{}\left(s,t\right)\mathclose{} - \kappa_1\right)\mathclose{}, \end{gather}

(2.18)\begin{gather}N_2^{rod} = a_2\mathopen{}\left(\varOmega_2\mathopen{}\left(s,t\right)\mathclose{} - \kappa_2\right)\mathclose{}, \end{gather}

(2.19)\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. 2008). 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

(2.20a,b)\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. (2017) 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:

(2.21)\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:

(2.22)\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}

(2.23)\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 2022). 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. 2008). 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

(2.24)\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

(2.25)\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

(2.26)\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}

(2.27)\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}

(2.28)\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 2013; Jabbarzadeh & Fu 2018). 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 2023), which is an MPI parallelized implementation of the IB method with adaptive mesh refinement (AMR). The software uses SAMRAI (Hornung & Kohn 2002; Hornung, Wissink & Kohn 2006) for AMR and PETSc (Balay et al. 1997, 2015, 2019) 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 1965). We use second-order finite differences to discretize the Laplacian, velocity gradients and stress divergence (Barrett 2019). The advective derivative is discretized with a second-order wave propagation algorithm (Ketcheson, Parsani & LeVeque 2013).

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 2009). The details of the spatial and temporal discretizations are discussed in previous work (Barrett 2019; Gruninger et al. 2024).

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. 2017). 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 1. Table of physical and computational parameters for both the model flagellum and the fluid.

Figure 1. 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}$.

Figure 2. 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.

Figure 3. 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 2008; Li et al. 2017). While the shape changes only slightly in viscoelastic fluids, shape changes can have drastic effects on flagellar bundling (Lee et al. 2018), 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

(3.1)\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).

Figure 4. 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 2014); 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.

Figure 5. 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.

Figure 6. 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

(3.2)\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}$.

Figure 7. 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. 2014; Gómez et al. 2017; Hyakutake et al. 2019; Qu & Breuer 2020) 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. 2010; Thomases & Guy 2014; Li et al. 2017). 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. 2017). 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.

Figure 8. 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.

Figure 9. 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. 2018).

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. 2000; Darnton et al. 2007; Ko et al. 2017), and this model has been successfully used to study flagellar bundling and locomotion in Newtonian fluids (Lim & Peskin 2012; Ko et al. 2017; Jabbarzadeh & Fu 2018). 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 (2017) 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 (2017).

The importance of shear-thinning behaviour for speed enhancement in helical swimmers has been demonstrated previously. Gómez et al. (2017) and Demir et al. (2020) 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 (2020), 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. 2017, 2021; Spagnolie et al. 2013), 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 (2020), 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 (2021), 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.