Propulsion by stiff elastic filaments in viscous fluids

TL;DR

This paper develops mathematical tools to analyze how stiff elastic filaments propel themselves in viscous fluids, revealing how flexibility affects propulsion speed and filament deformation, with applications to biological and artificial microswimmers.

Contribution

It introduces a mathematical framework for studying filament-driven propulsion in the stiff limit, including effects of deformation on propulsion in viscous environments.

Findings

Filaments deform by stretching, compression, and bending due to fluid forces.

Deformation influences propulsion speed depending on key parameters.

New insights into the physical interpretation of filament flexibility effects.

Abstract

Flexible filaments moving in viscous fluids are ubiquitous in the natural microscopic world. For example, the swimming of bacteria and spermatozoa as well as important physiological functions at organ-level, such as the cilia-induced motion of mucus in the lungs, or individual cell-level, such as actin filaments or microtubules, all employ flexible filaments moving in viscous fluids. As a result of fluid-structure interactions, a variety of nonlinear phenomena may arise in the dynamics of such moving flexible filaments. In this paper we derive the mathematical tools required to study filament-driven propulsion in the asymptotic limit of stiff filaments. Motion in the rigid limit leads to hydrodynamic loads which deform the filament and impact the filament propulsion. We first derive the general mathematical formulation and then apply it to the case of a helical filament, a situation relevant for the swimming of flagellated bacteria and for the transport of artificial, magnetically actuated motors. We find that, as a result of flexibility, the helical filament is either stretched or compressed (conforming previous studies) and its axis also bends, a new result which we interpret physically. We then explore and interpret the dependence of the perturbed propulsion speed due to the deformation on the relevant dimensionless dynamic and geometric parameters.

Equations

$${\partial }_s\mathbf{r}={\mathbf{d}}_3.$$

$${\partial }_s{\mathbf{d}}_i$$

$$\mathbf{D}$$

$${\partial }_s\mathbf{F}+\mathbf{K}$$

$${\partial }_s\mathbf{M}+{\mathbf{d}}_3\wedge \mathbf{F}+\mathbf{N}$$

$$\mathbf{M}=E{I}^{(1)}(\delta D_1){\mathbf{d}}_1+E{I}^{(2)}(\delta D_2){\mathbf{d}}_2+{\mu }_{S}J(\delta D_3){\mathbf{d}}_3,$$

$$I^{(1)}=I^{(2)}=\pi r^4/4,J=\pi r^4/2.$$

$$\mathbf{K}(s)$$

$${\zeta }_{\perp }$$

$${\zeta }_{\parallel }$$

$${\mathbf{N}}_{visc}=-{\zeta }_r({\mathbf{d}}_3.\mathbf{\omega }){\mathbf{d}}_3,$$

$${\zeta }_r=4\pi \mu r^2.$$

$$\frac{|{\mathbf{N}}_{visc}|}{|\mathbf{r}\wedge \mathbf{K}|}\sim \frac{\mu \omega r^2}{\mu \omega R_0^2}\sim {\left(\frac{r}{R_0}\right)}^2.$$

$${\mathbf{V}}_{\mathbf{rel}}$$

$$\mathbf{F}(L)=0,\mathbf{M}(L)=0,$$

$$\mathbf{F}(s)$$

$$\mathbf{M}(s)$$

$$\mathbf{M}(s)=$$

$$\mathbf{M}(s)+\mathbf{r}(s)\wedge \mathbf{F}(s)$$

$$\delta {\mathbf{d}}_i(s)=\delta \mathbf{\varphi }(s)\wedge {{\mathbf{d}}_i}^{(0)}(s).$$

$$\delta ({\partial }_s{\mathbf{d}}_i)$$

$$[{\partial }_s(\delta \mathbf{\varphi })]\wedge {{\mathbf{d}}_i}^{(0)}+(\delta \mathbf{\varphi })\wedge ({\mathbf{D}}^{(0)}\wedge {{\mathbf{d}}_i}^{(0)})-(\delta \mathbf{D})\wedge {{\mathbf{d}}_i}^{(0)}-{\mathbf{D}}^{(0)}\wedge (\delta \mathbf{\varphi }\wedge {{\mathbf{d}}_i}^{(0)})=0.$$

$$[{\partial }_s(\delta \mathbf{\varphi })+(\delta \mathbf{\varphi })\wedge {\mathbf{D}}^{(0)}-\delta \mathbf{D}]\wedge {{\mathbf{d}}_i}^{(0)}$$

$${\partial }_s(\delta \mathbf{\varphi })+(\delta \mathbf{\varphi })\wedge {\mathbf{D}}^{(0)}-\delta \mathbf{D}=0.$$

$${\partial }_s(\delta \mathbf{\varphi })$$

$$=(\delta D_i){{\mathbf{d}}_i}^{(0)}+{D_i}^{(0)}(\delta {\mathbf{d}}_i)-{D_i}^{(0)}(\delta \mathbf{\varphi }\wedge {{\mathbf{d}}_i}^{(0)})$$

$$=(\delta D_i){{\mathbf{d}}_i}^{(0)},$$

$$(\delta \mathbf{\varphi })(s)$$

$$(\delta \mathbf{\varphi })(0)=\mathbf{\Phi }.$$

$$\mathbf{M}=EI\sum _i(\delta D_i){{\mathbf{d}}_i}^{(0)}.$$

$$(\delta \mathbf{\varphi })(s)$$

$$\delta {\mathbf{d}}_3(s)$$

$$\delta \mathbf{r}(s)=$$

Full text

I Introduction

Many biophysical processes and engineering problems exhibit rich nonlinear behaviour due to fluid-structure interactions DowellHallFluidStructureInteraction2001 ; HuangAlbenFSIarticle2016 . Examples include tall buildings responding to winds Braun2009 , the flapping of bird wings Kang2011flapping ; WangZhangFlapping2016 and the motion of aircrafts ThomasDowellHall2002 . Physiological flows HeilHazelFSIInternalPhysiologicalFlows2011 provide additional examples, for example the vibrating vocal folds Horacek2002 , heart valves opening and closing with blood flow Sacks2009 , lungs expanding and contracting with breathing Otis1950Breathingmechanics ; Dubois1956Lungs and pulse propagation in blood vessels PedleyBloodFlowArteriesVeins2000 . Going all the way down to the microscopic world, one comes across elastic structures being deformed due to hydrodynamic loads from flows which are, in turn, affected by the deformation. Examples including deforming cilia BrennenWinet77 ; Bray2000 , fluctuating actin filaments Howard2001 , polymerising microtubules Hawkins2010 , the waving flagella of spermatozoa and the rotating flagellar filaments of bacteria BergAnderson1973 ; LaugaPowers2009 .

A prevalent morphology in the microscopic natural world is that of a helix. DNA BustamanteBryantSmith2003DNAMech ; Stricketal1996 , spirochaetes Lietal2000Spirochete ; SFGoldsteinCharon1988spirocheteleptospira , spiroplasma YangWolgemuthHuber2009Spiroplasma , trypanosoma Trypanosoma2015Stark and bacterial flagellar filaments BergAnderson1973 , all thake the shape of helices and spirals, as illustrated in Fig. 1 (A-C). In particular, the chiral shape of a helix is able to couple hydrodynamically rotation to translation. As such, it is used as the propulsive machinery for bacteria swimming in viscous fluids. The helical flagellar filaments of bacteria such as Escherichia coli (E. coli) are several micrometers in length and $20~{}\rm{nm}$ in radius, and are rotated at a typical frequency of $100~{}\rm{Hz}$ by specialised rotary motors SamuelBerg1996 ; Lauga2016BactlHydrdyn .

In the micro-engineering world, different types of externally powered motors have been proposed, studied and constructed. Examples include rigid helical propellers often termed ‘artificial bacterial flagella’, Zhang2009 ; helixpropGhoshFischer2009 , and motors that use flexible filaments Dreyfus2005 ; Gaonanowire2010 ; Paknanowire2011 , both of which are illustrated in Fig. 1 (D-F). One of the aspiring application of such artificial microswimmers is non-invasive medicine Sitti2009voyage ; BJNelson2010 , in which swimmers are to access targeted locations in the body, such as tumours, in order to deliver drugs Gaocargotowing2012 or perform delicate surgical tasks eyesurg2013 . In view of the former application, one interesting possibility is that of the motor being propelled by the drug itself taking the shape of long strips that twist up into a helical shape upon rotation LiTanZhang2009 . Magnetically propelled microswimmers with a flexible helix can be also be manufactured Huangetal2016AddressableShapeControl ; Paknanowire2011 and display rich nonlinear dynamics such as velocity profiles peaked at certain operational frequencies Paknanowire2011 . In addition, nonlinearity can be exploited in achieving selective control of microswimmers in large numbers, a desirable feature for any practical application KatsambaLauga2016 ; KeiCheang20143bead ; Frutiger2010 .

The multitude of biological systems that involve the dynamics of an elastic helix in a viscous fluid environment and the number of micro-engineering applications are compelling evidence for the need to study the elasto-hydrodynamic coupling theoretically. In the context of bacterial flagellar filaments, there have been several experimental, numerical and analytical investigations to address this problem. Theoretical and experimental studies were combined to derive the relationship between the elongation of a flagellar filament (modelled as a chain of segments) and the flow that it is subjected to HoshikawaKamiya1985 . Experimental studies showed different speeds between forward and backwards swimming of the single polar-flagellated bacterium Vibrio alginolyticus Magariyama2001forwardbackwardspeeddiff . These were followed by a numerical analysis for the deformation of a flagellum rotating in a viscous fluid TakanoGoto2003Numer . Studies of the deformation according to the Kirchhoff rod model combined with the Calladine model of the detailed structure of the filament provided analytical expressions for the bending moment, curvature and torsion of deformed flagellar filaments of swimming bacteria Vibrio alginolyticus and Salmonella. The deformation was obtained numerically and a comparison with experimental data provided an estimate of the elastic bending coefficient of the flagellar filament on the order of $EI\approxeq 10~{}\rm{pN}(\mu\rm{m})^{2}$ Takanoetal2003HelicFlagel ; Takanoetal2005HelicFlagelSalmonella .

The extension or compression of an elastic helix by flow and external fields has been studied analytically in the limit of a long, stiff helix KimPowers2005HelixDeformation . Force-extension curves were also derived in studies of helical springs subjected to uniaxial tension WadaNetz2007 . The propulsive force from a rotating, flexible, helical rod in a viscous fluid and the onset of the buckling instability above a critical rotation velocity have been studied by means of experiments and simulations Vogel2012thesis ; VogelStarkBuckling ; Jawedetal2015FlexHelixinViscousFluid . In particular, simulations based on a model that uses Kirchhoff’s classical elasticity theory for curved rods were used to investigate the transitions between the polymorphic forms of the bacterial flagellum Vogel2012thesis ; VogelStarkBuckling .

In this work, we put forward the mathematical framework necessary to address the steady-state locomotion of stiff elastic, slender filaments in viscous fluids. The motion of the filament induces a hydrodynamic load that deforms it. This in turn affects the kinematics because the shape has changed. Implementing the overall force and torque balance at leading order involves integrating the hydrodynamic load with the unknown velocity and rotation rate over the new deformed shape as if it were rigid, and inverting the system to solve for the unknown velocity and rotation rate. In order to obtain the perturbation to the rigid kinematics, one needs in particular to perturb the hydrodynamic resistance matrices. In this paper we show how to do this for a long, slender filament of arbitrary shape. We next apply our analytical framework to study the setup of an elastic helix that is rotating and translating, or equivalently is in the presence of such an external flow of a viscous fluid. We study the limit where the helix is very stiff, so that any deformation is small, and very long, so that it rotates about its long axis and does not wobble YiLaugawobbling . This setup is relevant to bacterial flagellar filaments and a popular design for magnetically actuated artificial microswimmers. The latter consists of a flexible helix clamped onto a magnetic head on which an external magnetic torque is exerted. We calculate below the full three-dimensional deformation analytically and its feedback on the swimming speed of a bacterium during swimming in a straight line. Our results of the deformation agree with previous analytical results of the extension/compression in Ref. KimPowers2005HelixDeformation , and capture and explain the bending of the helix axis that was observed in the numerical results of Ref. Takanoetal2003HelicFlagel and whose origin has been unclear so far.

Our paper is organised as follows. In §II we outline the mathematics framework for the steady motion of a stiff-elastic, slender filament of any shape in a viscous fluid. We show how the Kirchhoff model for an elastic rod (§II.1) combined with resistive-force theory for the hydrodynamic load (§II.2) lead to the deformation (§II.3) and how this in turn perturbs the leading-order kinematics (§II.4). The latter is obtained by implementing the force and torque balance (§II.5). We formulate this dynamic balance for two specific setups: firstly that of a flexible filament actuated by a magnetic torque exerted on the head on which it is clamped (§II.6) and secondly the case of swimming bacterium rotating a flexible flagellar filament relative to its cell body (§II.7). Since the prevalent geometry for these two cases is a helical one, in §III we give the common details of the hydrodynamic load, bending moment and the deformation of the helix (§III.1). We next interpret the bending of the helix axis and proceed to investigate the feedback of the deformation on the kinematics (§III.2). We then calculate the perturbations to the resistance matrices due to the small deformation. Applying the appropriate forms for the force and torque balance, in both the artificial (§III.2.4) and biological (§III.2.5) setup, we derive the perturbation of the swimming velocity and discuss the physical interpretations of our results.

II Propulsion by elastic filaments: General framework

We first consider in this section the propulsion of elastic filaments of arbitrary shapes, and apply it to two cases: externally-actuated microswimmer propulsion driven by a rotating magnetic field and bacterial propulsion. In both cases, our setup involve a long elastic filament with one end clamped on a head or cell body and which rotates with angular velocity $\boldsymbol{\omega}$ (in the laboratory frame) and translates with velocity $\mathbf{U}$. The particular case of helical filaments will be considered in §III.

In the context of artificial microswimmers, the head is magnetised and actuated by an external magnetic field, $\mathbf{B}$, rotating about a fixed axis, say the $z$-axis. Assuming the head to have a constant dipole moment, $\mathbf{m}$, the field exerts a magnetic torque $\mathbf{M}_{mag}=\mu_{0}\mathbf{m}\wedge\mathbf{B}$ on the head which thus rotates and, with the proper filament shape, also translates. In the case of bacteria, there is instead a rotary motor embedded in the cell that produces a torque, $M_{mot}$, actuating the flagellar filament in rotation. The filament, which in many cases is stiff, is connected to the motor via a very flexible short hook BFhookSamateuetalNamba2004 .

Following the actuation of the magnetised head or the bacterial motor, the elastic filament attached onto it will rotate and (if of the right shape) also translate. Due to hydrodynamic loads, it will also start deforming, until it reaches shape equilibrium. In this work we aim at characterising the steady state equilibrium configuration of the filament obtained after all transients, where the shape no longer changes and for which elastic and hydrodynamic stresses balance. The swimming kinematics are then governed by the force and torque balance over the entire swimming organism (or device), which involve integrating hydrodynamic loads over the new deformed shape. The fluid-structure interactions manifests themselves therefore via the deformation induced by the hydrodynamic load and by the feedback of the deformation on the kinematics via dynamic balance.

In both applications, the magnetised head or cell body and the filament will translate at the same speed in steady state. In the case of the artificial motor, the head and tail will share the same rotation rate since the tail is clamped onto the head. In our model of bacterium, we do not include details of the hook, but instead assume that its high flexibility allows different rotation rates between the cell body and flagellar filament, the values of which are determined by torque balance.

II.1 Elastic Rod

Working in the frame of reference of the filament, let $\mathbf{r}(s)$ denote the location of its centreline and $\{\mathbf{d}_{i}(s)\}_{i=1,2,3}$ the local material frame (which we will take later to coincide with the Frenet-Serret frame) so that $\mathbf{d}_{3}$ is tangent to the centreline,

[TABLE]

The configuration of the material frame along the filament is then described by

[TABLE]

where $\mathbf{D}$ is the Darboux vector, measuring the strains in the rod. Its components in the material frame are the material curvatures $(D_{1}=\kappa^{(1)},D_{2}=\kappa^{(2)})$ and the material twist $(D_{3}=\tau)$ of the rod.

The elastic behaviour of the filament is governed by the classical Kirchhoff equations for a rod AudolyPomeau , which give the balance of forces and moments on a cross section

[TABLE]

where $\mathbf{F}$ is the internal force acting on a cross section of the rod, $\mathbf{M}$ is the bending moment and $\mathbf{K}$, $\mathbf{N}$ are the distributed (external) force and torque densities respectively. The constitutive equation for a Hookean material gives the bending moment as AudolyPomeau

[TABLE]

where $\delta D_{i}=D_{i}-{D_{i}}^{(0)}$ are the deviations of the material curvatures and twist in the deformed state from their values in the reference configuration. In the paper, we add a [math] superscript or subscript to indicate quantities pertaining to the reference configuration. As appropriate for a linearly elastic material, the deviations of the curvatures and twist are therefore linearly related to the components of the bending moment via Young’s modulus, $E$, and the shear modulus, $\mu_{S}$, of the material and depend on three geometrical coefficients: the principal moments of inertia, $I^{(1)}$ and $I^{(2)}$, and the twist rigidity, $J$, of the rod. For a rod with a circular cross-section of radius $r$ these take the values

[TABLE]

II.2 Hydrodynamic load

II.2.1 Resistive-force theory

The viscous tractions due to the motion of a slender filament in a viscous fluid at low Reynolds number are well captured by resistive-force theory gray55 ; Lighthill1976FlagellarHydrodynamics ; Rodenborn2012 ; LaugaPowers2009 . This technique integrates fundamental solutions of the Stokes equations along the centreline of a slender filament to give an expression for the local hydrodynamic force per unit length, $\mathbf{K}$, exerted on the filament due to its motion in a viscous fluid. At the position labelled by the contour-length parameter value $s$, the local hydrodynamic force per unit length exerted on the filament is then given

[TABLE]

where $\mathbf{V_{rel}}(s)$ is the local relative velocity between the filament and the fluid at that position and $\zeta_{\parallel},\zeta_{\perp}$ are the drag coefficients for motion parallel and perpendicular to the local tangent of the filament. There are many forms of the drag coefficients in various geometries Hancock1953 ; GrayHancock1955 ; Lighthill1976FlagellarHydrodynamics ; Rodenborn2012 and for the purpose of studying helical filaments in §III, we will use Lighthill’s coefficients given by

[TABLE]

where $\mu$ is the dynamic viscosity of the fluid and $\Lambda,r$ are the contour wavelength and cross-sectional radius of a helical filament respectively. The ratio of the two coefficients is approximately equal to $1/2$,

[TABLE]

The fact that the perpendicular drag coefficient is approximately twice the parallel one is the extension to a curved filament of the fact that it is approximately twice as hard to pull a rod through a viscous fluid in a direction perpendicular to its length than lengthwise. This drag anisotropy is at the heart of locomotion of microorganisms and artificial microswimmers. For example, it is the crucial ingredient coupling rotation to translation for bacterial flagellar filaments Lauga2016BactlHydrdyn .

There are two sources of hydrodynamic moments acting on the filament. The first one is the moment due to the distribution of forces in Eq. 8. The second is the viscous torque per unit length, $\mathbf{N}_{visc}$, that resists the rotation of an element of the rod about $\mathbf{d}_{3}$ (i.e. its centreline) given by

[TABLE]

where $\zeta_{r}$ is the rotational drag coefficient

[TABLE]

This second viscous torque, $\mathbf{N}_{visc}$,

can be typically neglected for filaments with sufficiently-small cross-sectional radius $r$. Consider a helical filament with helical radius ${R}_{0}$ in its reference configuration and compare the magnitude of the moment (measured with respect to the helical axis) due the hydrodynamic force per unit length, $\mathbf{K}$, with that of $\mathbf{N}_{visc}$ in the expression for the bending moment, $\mathbf{M}$, in Eq. 18. Their ratio scales as KimPowers2005HelixDeformation

[TABLE]

This ratio is very small unless one is dealing with nearly straight filament configurations, which is typically not the case for the helical geometry of bacteria and artificial micromotors. The second viscous torque, $\mathbf{N}_{visc}$, can therefore be safely neglected.

Once steady state has been reached, the filament rotates with uniform angular velocity, $\boldsymbol{\omega}$, and translates with uniform velocity, $\boldsymbol{U}$. We may then consider the frame in which the filament is stationary. The relative velocity between the filament and the fluid is given by

[TABLE]

which allows access to the leading-order estimate for the force density, $\mathbf{K}$, in Eq. 8.

The next step is to calculate the internal force, $\mathbf{F}$, and bending moment, $\mathbf{M}$, by integrating the Kirchhoff equations, 4 and 5. We assume that the end-point at $s=L$ of the filament is free,

[TABLE]

allowing to obtain explicitly

[TABLE]

Using that $\mathbf{d}_{3}(s^{\prime})=\partial_{s^{\prime}}[\mathbf{r}(s^{\prime})]$ we can rewrite $\partial_{s^{\prime}}[\mathbf{r}(s^{\prime})]\wedge\mathbf{F}(s^{\prime})=\partial_{s^{\prime}}[\mathbf{r}(s^{\prime})\wedge\mathbf{F}(s^{\prime})]-\mathbf{r}(s^{\prime})\wedge\partial_{s^{\prime}}[\mathbf{F}(s^{\prime})]$, integrate by parts and use Eq. 4 to obtain

[TABLE]

Using the boundary conditions of Eq. 15, we then obtain the integral formula

[TABLE]

As we will see in §II.6, the second term in the left-hand side of Eq. 19 arises naturally when one considers the torque balance with respect to the origin, as the bending moment, $\mathbf{M}(s)$, is defined with respect to the centre of the cross-section at arclength position $s$ along the filament.

II.3 Deformation

As a result of the hydrodynamic forcing, the filament will deform. In the limit of stiff filaments, i.e. of high Young’s modulus, the filament will undergo small deformations and the material frame will be slightly perturbed. As the material frame needs to stay orthonormal, this can be represented by a set of three small local rotation vectors, $\delta\boldsymbol{\phi}(s)$, along the rod Audoly-Pomeau-Elasticity-and-geometry-2000 so that

[TABLE]

In order to relate the bending moment to the deformation, we need to relate the perturbations to the components of the Darboux vector, $\delta D_{i}$, that appear in the constitutive equation, Eq. 6, to the small rotations $\delta\boldsymbol{\phi}(s)$ of Eq. 20, following Ref. AudolyPomeau . This is done by considering the first-order perturbations to the two sides of Eq. 2. Perturbing the left-hand side of Eq. 2 and using Eq. 20 leads to

[TABLE]

Perturbing next the right-hand side of Eq. 2 leads to $(\delta\mathbf{D})\wedge{\mathbf{d}_{i}}^{(0)}+{\mathbf{D}}^{(0)}\wedge(\delta\mathbf{d}_{i})$ and equating the two expressions we obtain

[TABLE]

Using the Jacobi identity for the second and last terms of Eq. 22 leads to

[TABLE]

for all $i$. We thus have that

[TABLE]

Rearranging and expressing quantities in the reference configuration material frame gives

[TABLE]

where we used the summation convention for repeated indices and where we used $\delta\mathbf{D}=\delta D_{i}{\mathbf{d}_{i}}^{(0)}+{D_{i}}^{(0)}\delta{\mathbf{d}_{i}}^{(0)}$. Integrating Eq. (27) then one obtains the expression for the infinitesimal rotation $\delta\boldsymbol{\phi}(s)$

[TABLE]

Note that the lower limit of the integral makes use of the boundary condition

[TABLE]

A clamped boundary condition at $s=0$ means $\mathbf{\Phi}=\mathbf{0}$, however for a free end at $s=0$ one would have to solve for the value of $\mathbf{\Phi}$.

Using the classical assumption that $EI^{(1)}=EI^{(2)}=\mu J$ KimPowers2005HelixDeformation ; VogelStark2010 ; VogelStarkBuckling , we next note that the constitutive equation for the bending moment, Eq. 6, simplifies to

[TABLE]

Thus $\delta\boldsymbol{\phi}$ in Eq. 28 is obtained by integrating the bending moment along the filament as

[TABLE]

The perturbation to the tangent vector is given by Eq. 20, hence we obtain

[TABLE]

Integrating Eq. 32 next gives access to the deformation

[TABLE]

where the perturbation to the position of the $s=0$ end-point is denoted $\Delta\mathbf{r}_{0}$; it vanishes if the end-point has fixed position, otherwise for a free end-point, $\Delta\mathbf{r}_{0}$ would have to be solved for.

II.4 Feedback of the deformation on the kinematics

The swimming kinematics of the filament-body/head pair are determined by the force and torque balance. For large values of the filament Young’s modulus, the deformations are small and the kinematics are given by those of a rigid filament to leading order. In this section we quantity the impact of the deformation $\delta\mathbf{r}$ on the kinematics, i.e. calculate the next-order effect on both $\mathbf{U}$ and $\boldsymbol{\omega}$. The perturbation principle behind the calculation of the feedback of the deformation on the kinematics is illustrated in Fig. 2. Once the cell or artificial motor starts moving, the filament starts deforming due to the hydrodynamic load until it reaches steady state after which it undergoes rigid body motion. The kinematics will be slightly different from the rigid-body ones, because the hydrodynamic load is integrated along a slightly perturbed shape along which force and torque balance is enforced. We first express the force and torque balance for a general shape of the filament, leading to a linear system involving the hydrodynamic load integrated over the shape of the entire filament. Inputting the reference configuration shape gives the rigid filament kinematics while perturbing the system leads to a relationship linearly relating the perturbed kinematics to the perturbed resistance matrices of the filament, which themselves are linear functions of the deformation. Once inverted, this procedure gives access to the perturbation in swimming kinematics.

II.5 Relevant forces and torques

Using the term ‘head’ to mean either the cell body of the bacterium or the magnetised head of the artificial motor, the hydrodynamic force, $\mathbf{F}_{h,visc}$, and torque with respect to the origin, $\mathbf{M}_{h,visc}$, on the head translating with velocity $\mathbf{U}$ and rotating with rate $\boldsymbol{\omega}_{h}$ are given in a Stokes flow by the linear relationship KimKarrilaMicrohydrodynamicsBook

[TABLE]

where $\boldsymbol{\alpha}_{h},\boldsymbol{\beta}_{h},\boldsymbol{\gamma}_{h}$ and $\boldsymbol{\delta}_{h}$ are the resistance matrices for a spherical head of radius $a_{h}$. With respect to the centre of the head these are

[TABLE]

where $\mathbf{1}$ denotes the identity matrix.

If the centre of the head is at position $\mathbf{r}_{h}$ with respect to the origin, then the hydrodynamic torque on the head has an extra term given by $\mathbf{r}_{h}\wedge\mathbf{F}_{h,visc}$. This can be seen by writing the torque with respect to the origin, $\int_{S}\mathbf{r}\wedge(\boldsymbol{\sigma}\cdot\mathbf{n})\;\mathrm{d}S$, where $\boldsymbol{\sigma}$ is the viscous stress tensor and $\mathbf{n}$ is the normal to the surface $S$ of the head pointing outwards, as $\int_{S}(\mathbf{r}-\mathbf{r}_{h})\wedge(\boldsymbol{\sigma}\cdot\mathbf{n})\;\mathrm{d}S+\int_{S}\mathbf{r}_{h}\wedge(\boldsymbol{\sigma}\cdot\mathbf{n})\;\mathrm{d}S$, with the first integral being the torque with respect to the centre of the head $\mathbf{r}_{h}$, and the second integral simplifying to $\mathbf{r}_{h}\wedge\int_{S}(\boldsymbol{\sigma}\cdot\mathbf{n})\;\mathrm{d}S=\mathbf{r}_{h}\wedge\mathbf{F}_{h,visc}$. This adds the extra terms $\tilde{\boldsymbol{\gamma}}_{h}=\left(\underline{\underline{\underline{\boldsymbol{\epsilon}}}}\cdot\mathbf{r}_{h}\right)\cdot\boldsymbol{\alpha}_{h}$ and $\tilde{\boldsymbol{\delta}}_{h}=\left(\underline{\underline{\underline{\boldsymbol{\epsilon}}}}\cdot\mathbf{r}_{h}\right)\cdot\boldsymbol{\beta}_{h}$ to $\boldsymbol{\gamma}_{h}$ and $\boldsymbol{\delta}_{h}$ respectively, where the third rank tensor $\underline{\underline{\underline{\boldsymbol{\epsilon}}}}$ is the Levi-Civita tensor with components $\epsilon_{ijk}$ such that $\left(\underline{\underline{\underline{\boldsymbol{\epsilon}}}}\cdot\mathbf{r}\right)\cdot\boldsymbol{\omega}=\mathbf{r}\wedge\boldsymbol{\omega}$. In the case of a spherical head, these extra terms become $-6\pi\mu a_{h}\left(\underline{\underline{\underline{\boldsymbol{\epsilon}}}}\cdot\mathbf{r}_{h}\right)$ and $-8\pi\mu a_{h}^{3}\left(\underline{\underline{\underline{\boldsymbol{\epsilon}}}}\cdot\mathbf{r}_{h}\right)$ respectively. We will later focus on the case of a spherical head whose centre lies on the axis of a helical filament that is rotating and translating about its axis, in which case the cross products of $\mathbf{r}_{h}\wedge\mathbf{U}$ and $\mathbf{r}_{h}\wedge\boldsymbol{\omega}_{h}$ vanish, and hence also the components of the extra terms $\tilde{\boldsymbol{\gamma}}_{h},\tilde{\boldsymbol{\delta}}_{h}$ along the relevant axis also vanish.

We can write similar expressions for the force and torque (with respect to the origin) acting on the filament using Eqs. 16 and 19,

[TABLE]

where $\boldsymbol{\mathcal{A}}_{fil},\boldsymbol{\mathcal{B}}_{fil},\boldsymbol{\mathcal{C}}_{fil},\boldsymbol{\mathcal{D}}_{fil}$ are the resistance matrices of the filament

[TABLE]

and the third rank tensor $\underline{\underline{\underline{\boldsymbol{\epsilon}}}}$ in Eqs. 42 and 43 is the Levi-Civita tensor with components $\epsilon_{ijk}$ such that $\left(\underline{\underline{\underline{\boldsymbol{\epsilon}}}}\cdot\mathbf{r}\right)\cdot\boldsymbol{\omega}=\mathbf{r}\wedge\boldsymbol{\omega}$. Notice that $\boldsymbol{\mathcal{A}}=\boldsymbol{\mathcal{A}}^{\intercal},~{}\boldsymbol{\mathcal{D}}=\boldsymbol{\mathcal{D}}^{\intercal},~{}\boldsymbol{\mathcal{B}}=\boldsymbol{\mathcal{C}}^{\intercal}$, as expected for Stokes flows KimKarrilaMicrohydrodynamicsBook .

In the case of artificial propellers, the head and filament rotate at the same rate, so we may define overall resistance matrices for the head and filament together which are just the sums of the corresponding matrices

[TABLE]

In the case of a bacterium, the rotation rate of the head (i.e. the cell body) is different from that of the filament; in fact the head rotates in the opposite direction in order to satisfy the overall torque balance. The difference in actuation of the filament in the artificial and biological applications is summarised in Fig. 3 and we now consider the force and torque balance for each case separately.

II.6 Artificial propellers with a filament clamped on a magnetised head

We first focus on the artificial propeller actuated by a rotating magnetic field exerting a torque on the magnetised head on which the filament is clamped (Fig. 3, left). The internal force acting through the cross-section at $s=0$, transmitted from the filament to the head, needs to balance the viscous hydrodynamic force, $\mathbf{F}_{h,visc}$, resisting the motion of the head. Although the same is true for the axial torque balance on the head, care must be taken because the bending moment, $\mathbf{M}(0)$, transmitted through the cross-section at $s=0$ from the filament to the head is calculated with respect to the centre of the cross-section of the filament while viscous moments on the head and measured by its centre. The torque balance includes therefore the viscous hydrodynamic torque resisting the motion of the head, $\mathbf{M}_{h,visc}$, the external magnetic moment, both defined with respect to the origin, as well as the bending moment with respect to the origin, $\mathbf{M}(0)+\mathbf{r}(0)\wedge\mathbf{F}(0)$.

The dynamic equations are therefore written as

[TABLE]

Using the results of Sec. II.5, we see that the system to be solved in order to find $\mathbf{U},\boldsymbol{\omega}$ in matrix form is

[TABLE]

The leading-order kinematics are given by inverting the linear system

[TABLE]

For the first-order perturbation one needs to evaluate the perturbations to the resistance matrices

[TABLE]

The next order perturbation to the force and torque balance in Eq. 52 is

[TABLE]

For a more concise notation let us use $\boldsymbol{\mathcal{R}}$ to denote the large resistance matrix and $\boldsymbol{\mathcal{M}}$ to denote its inverse, the large motility matrix

[TABLE]

Inverting Eqs. 52 and 59 then leads to

[TABLE]

For analytically tractable calculations, we focus on long chiral filaments that have the $z$-axis as their long axis of rotation. As a result we can assume translation along and rotation about the $z$-axis, i.e. $\mathbf{U}=U\mathbf{e}_{z},~{}\boldsymbol{\omega}=\omega\mathbf{e}_{z}$, and only consider the axial components of the force and torque balances. For example, the $z$-component of $\boldsymbol{\mathcal{A}}\cdot\mathbf{U}$ in Eq. 51 reduces to just $\mathcal{A}_{zz}U$. In this one-dimensional limit, we will drop the $zz$ indices for notation convenience and use $\mathcal{A}$ to mean $\mathcal{A}_{zz}$, and similarly for other matrix components.

The leading-order result of Eq. 61 describes a rigid filament in its reference configuration with both translation and rotation proportional to the external moment, $M_{mag}$, as

[TABLE]

while for the next-order correction we obtain

[TABLE]

II.7 Biological locomotion with a filament rotated by a motor

In the case of swimming bacteria (Fig. 3, right), a motor embedded in the cell wall applies a constant torque, $M_{mot}$, via the short hook in order to rotate a long filament. The rotation rate of the filament in the laboratory frame is $\boldsymbol{\omega}_{fil}$ while the head rotates at a different rate denoted by $\boldsymbol{\omega}_{h}$. Both the head and the filament translate at the same velocity, $\mathbf{U}$.

In this case, there are three dynamic balances to consider: the overall force balance as well as the torque balances on the filament and the head

[TABLE]

Substituting in the terms of the resistance matrices we obtain

[TABLE]

Inverting Eq. 73 for $\boldsymbol{\omega}_{h}$ and substituting in Eq. 71 leads to

[TABLE]

The system to be solved can thus be written in matrix form as

[TABLE]

which in the special case of a spherical head further reduces to

[TABLE]

Notice how the system in Eq. 76 is mathematically similar to that for the artificial motor in Eq. 51 if one makes the substitution $\boldsymbol{\mathcal{A}}\to\boldsymbol{\mathcal{A}}$, $\boldsymbol{\mathcal{B}}\to\boldsymbol{\mathcal{B}}_{fil}$, $\boldsymbol{\mathcal{C}}\to\boldsymbol{\mathcal{C}}_{fil}$, $\boldsymbol{\mathcal{D}}\to\boldsymbol{\mathcal{D}}_{fil}$, and $\mathbf{M}_{mag}\to\mathbf{M}_{mot}$. Making these substitutions in the projected one-dimensional versions of the artificial motor equations, Eqs. 64-67, allows to obtain the corresponding ones for a bacterium translating along, and rotating about, the axis of its long chiral filament. The rigid result is

[TABLE]

while the first correction is given by

[TABLE]

III Propulsion by elastic helical filaments

Up to this point in the paper, we have formulated everything in terms of an arbitrary shape of the filament in the rigid limit. In this section we focus on the case of filaments with helical geometry, as this is the most common shape in both cases of interest discussed thus far, and we carry out the calculations of the expressions we formulated in §II. As done in §II.6 and II.7, in order to keep the analytical calculations tractable we will only study translation along, and rotation about, the $z$-axis of the helical shape, and not the full three-dimensional (3D) motion. This is an appropriate limit to consider if the filament is long enough to not wobble YiLaugawobbling .

III.1 Calculating the deformation

III.1.1 Reference configuration: Centreline geometry

We now compute the steady-state perturbation results where we assume the actuation from the rotating magnetic field or bacterial motor is weak enough, or the helix stiff enough, that the helix deforms only slightly. The reference configuration is a long, uniform, stress-free rod in the shape of a helix of uniform pitch ${P}_{0}$, radius ${R}_{0}$, helix angle ${\alpha}_{0}$, such that $\tan{{\alpha}_{0}}=(2\pi{R}_{0}/{P}_{0})$ and its axis is aligned with the $z$ axis. We define the chirality index $h$ which takes the value $\pm 1$ according to whether the helix is right-handed (RH, $h=1$) or left-handed (LH, $h=-1$). In its reference configuration, the centreline of the helix is therefore given by

[TABLE]

where ${\Lambda}_{0}=\sqrt{{P}_{0}^{2}+4\pi^{2}{R}_{0}^{2}}~{}$ is the helix wavelength measured along the arclength $s$. This geometry is illustrated in Fig. 4.

Taking the material frame $\{{\mathbf{d}_{1}}^{(0)},{\mathbf{d}_{2}}^{(0)},{\mathbf{d}_{3}}^{(0)}\}$ to coincide with the Serret-Frenet frame $(\mathbf{n},\mathbf{b},\mathbf{t})$ similarly to Ref. KimPowers2005HelixDeformation , we use $\partial_{s}{\mathbf{r}}_{0}(s)=\mathbf{t}={\mathbf{d}_{3}}^{(0)}$, $\partial_{s}\mathbf{t}={\kappa}_{0}\mathbf{n}$ and $\mathbf{b}=\mathbf{t}\wedge\mathbf{n}={\mathbf{d}_{2}}^{(0)}$ to obtain

[TABLE]

with

[TABLE]

and, using $\partial_{s}\mathbf{b}=-{\tau}_{0}\mathbf{n}$,

[TABLE]

Identifying the Serret-Frenet equations

[TABLE]

with $\partial_{s}{\mathbf{d}_{i}}^{(0)}={\mathbf{D}}^{(0)}\wedge{\mathbf{d}_{i}}^{(0)}$, gives the components of the Darboux vector

[TABLE]

The material frame at various positions along the filament is illustrated in Fig. 5, with the tangent, normal and binormal vectors shown in dash-dotted red, solid green, and dashed blue respectively.

III.1.2 Hydrodynamic load

Actuated by the magnetic field, once steady state has been reached the helix rotates with uniform angular velocity $\omega_{fil}\mathbf{e_{z}}$ and translates at uniform velocity $U\mathbf{e_{z}}$. While it is deformed from its reference configuration shape due to the forces and torques acting on it, it does undergo rigid body motion as its shape is no longer changing. We consider the frame in which the helix is stationary. The relative velocity between the helix and the fluid is given by

[TABLE]

where we approximated the shape to the reference helical one and used the rigid helix kinematics. This allows us to calculate the leading-order force density using Eq. 8 and noting that ${\mathbf{d}_{3}}^{(0)}\cdot{\mathbf{V}_{rel}}^{(0)}={\omega_{fil}}^{(0)}h{R}_{0}\sin{\alpha}_{0}+{U}^{(0)}\cos{\alpha}_{0}$,

[TABLE]

Rearranging gives

[TABLE]

where

[TABLE]

Integrating the first Kirchhoff equation, $\partial_{s}\mathbf{F}+\mathbf{K}=0$, with the boundary condition that $\mathbf{F}(L)={\bf 0}$ at the free end gives access to the distribution of internal force as

[TABLE]

Similarly, integrating the second Kirchhoff equation, $\partial_{s}\mathbf{M}+\mathbf{d}_{3}\wedge\mathbf{F}+\mathbf{N}_{visc}=\mathbf{0}$, with the free-end boundary condition, $\mathbf{M}(L)=\bf 0$, allows to compute the bending moment as

[TABLE]

III.1.3 Leading-order kinematics

For a stiff elastic filament, the leading-order kinematics are given by the rigid limit. The $z$-components of the expressions for the total force and torque give the $zz$-components of the resistance matrices for a rigid helical filament,

[TABLE]

The $zz$-components of the overall resistance matrices for both the helical filament and the spherical head are

[TABLE]

III.1.3.1 Unified approach.

We now make use of the similarity in the expressions for the kinematics and their perturbation for the artificial helical motors and the swimming bacterium, noting that the formulae for the artificial motors can be changed into those for bacteria by replacing $\mathcal{D}$ by $\mathcal{D}_{fil}$, and ${M}_{mag}$ by $M_{mot}$. We thus write

[TABLE]

where $M$ is to be substituted by either $M_{mag}$ or $M_{mot}$ for artificial motors and bacteria respectively, and write

[TABLE]

where we use the indicator function $\mathbbm{1}$ to take the value $1$ for the case of artificial motors and [math] for bacteria. The coefficients $\mathcal{A}$ and $\mathcal{B}$ are defined as in Eqs. 100 and 101.

III.1.3.2 Artificial motors.

The leading-order kinematics, given by Eqs. 103 and 104, are obtained as

[TABLE]

where the denominator is given explicitly by

[TABLE]

We note that the factor $\sin{\alpha}_{0}\cos{\alpha}_{0}$ present in the numerator of Eq. 107 gives vanishing speeds for helix angles [math] or $\pi/2$; this is expected by symmetry since chirality is lost in these two limits. In particular, we note that taking the limit of the helix angle to zero while keeping the total contour length $L$ and the number of turns $n$ fixed means that the helical radius ${R}_{0}=L\sin{\alpha}_{0}/(2\pi n)$ is also shrinking to zero. This limit gives vanishing speeds since the expression for the denominator given in Eq. 109 is non-zero due to the presence of the terms involving the head and the viscous rotational torque coefficient ($\zeta_{r}$).

The leading-order results in the limit of a filament which is slender ($r\ll L$) and long compared to the size of the magnetised head $(a_{h}\ll L)$, and for a non-vanishing helix angle are

[TABLE]

with relative errors of order ${\mathcal{O}}\left(r^{2}/L^{2},a_{h}/L\right)$.

III.1.3.3 Bacteria.

In the case of a swimming bacterium, the leading-order kinematics, given by Eqs. 103 and 104, are

[TABLE]

where the denominator is given by

[TABLE]

Here also the speed vanishes in the limiting cases of helix angles [math] or $\pi/2$ as expected. When the cell body is small compared to the length of the slender flagellar filament, $a_{h}\ll L$, in the case of a non-vanishing helix angle we can further simplify Eq. 112 as

[TABLE]

with relative errors of order ${\mathcal{O}}\left(r^{2}/L^{2},a_{h}/L\right)$.

III.1.4 Deformation for clamped Helix

The full details of the deformation of the helix are shown in Appendix A. Under the assumption that $EI^{(1)}=EI^{(2)}=\mu J$, we obtain at location $s/{R}_{0}\gg 1$

[TABLE]

The $z$-component of the deformation in Eq. 117 giving the extension/compression is in agreement with the results obtained in Ref. KimPowers2005HelixDeformation . In addition, the non-zero $x,y$ components of the deformation lead to bending of the helical axis. Such bending of the axis was observed in the numerical results of Ref. Takanoetal2003HelicFlagel , but its origin remained unclear. We now discuss and interpret these two aspects of the deformation.

III.1.5 Compression or extension?

The total amount of extension / compression is given by the $z$-component of the deformation evaluated at $s=L$,

[TABLE]

This recovers exactly the result of Ref. KimPowers2005HelixDeformation which investigated the compression/extension of a clamped helical filament subject to a uniform translating flow or a rotating flow. The contribution from the axial flow is the motion of the microswimmer propelling with a positive axial velocity, ${U}^{(0)}>0$, which is equivalent to the swimmer being fixed and subject to a flow $-{U}^{(0)}\mathbf{e}_{z}$. This leads to compression, regardless of the helix handedness, as expected from intuition. In contrast, the part due to the rotational flow predicts that a helix rotating in the positive sense about the $z$-axis $({\omega}_{0}>0)$ will be extended if it is right-handed (RH, $h>0$) and compressed if left-handed (LH, $h<0$).

In the cases we are investigating, the kinematics are set by the force and torque balances and ${U}^{(0)}$ and ${\omega_{fil}}^{(0)}$ are not independent. Instead both are proportional to the actuating torque, as per Eqs. 103 and 104. A helix actuated by a positive torque, $M>0$, will rotate in the positive sense (${\omega_{fil}}^{(0)}>0$) regardless of its handedness, but will translate with a positive velocity along the $z$-axis if RH and negative if LH. The effect of the translation is then that a RH helix ($h>0$) actuated with a positive magnetic torque will have a positive velocity, ${U}^{(0)}>0$, and hence be compressed, whereas a LH helix will move in the opposite direction, ${U}^{(0)}<0$, and will thus be extended. In contrast, the effect of rotation is the opposite to that of translation. Helices of both handedness will rotate in the positive sense, and a RH helix is extended as it is experiencing a rotating flow in the direction that ‘uncoils’ it, whereas a LH helix is compressed. As a result, translation and rotation have opposite effects in terms of extension/compression.

Substituting Eqs. 103 and 103 in Eq. 120 gives

[TABLE]

where

[TABLE]

and $(\mathcal{A}\mathcal{D}_{*}-\mathcal{B}^{2})$ in the denominator is given in Eq. 106. Under the long and slender approximations and for small head, this simplifies to

[TABLE]

As a result, the effect of rotation dominates and a RH helix ($h>0$) is always extended whereas a LH helix is compressed.

III.1.6 Bending of the helix axis

Far from the clamped end, i.e. for $s\sim L\gg{R}_{0}$, the deformation is given by

[TABLE]

with relative errors of order ${\mathcal{O}}\left({{R}_{0}}/{L},{r}/{L}\right)$. We can see from Eq. 124 that the helix axis bends parallel to the direction $\left(\cos\left(\frac{2\pi L}{{\Lambda}_{0}}\right),h\sin\left(\frac{2\pi L}{{\Lambda}_{0}}\right),0\right)$, i.e. the $xy$ projection of the position vector of the free end-point ${\mathbf{r}}_{0}(s=L)$ in the reference configuration. In order to physically interpret this result we evaluate the forces and torques acting at the clamped end of the filament

[TABLE]

and

[TABLE]

Hence for ${R}_{0}/L\ll 1$,

[TABLE]

We can now explain the bending of the axis of the helix using an ‘effective’ rod analogy, as has been used in a variety of settings BiezenoGrammel1939 ; Haringx1942 ; Haringx1948 ; Vogel2012thesis . Consider an effective rod around which the helix is wind, following the shape of the helix axis. This is a straight rod clamped at the origin and parallel with the $z$-axis in its reference configuration, so that the material frame of the effective rod is

[TABLE]

A bending torque, $\tilde{\mathbf{M}}=\mathcal{M}{\tilde{\mathbf{d}_{1}}}^{(0)}$ with $\mathcal{M}=(A_{x}{R}_{0})\cot{\alpha}_{0}L$, exerted at the clamped end of the rod will lead to a perturbation to the Darboux vector of the rod with only non-zero component

[TABLE]

The resulting infinitesimal rotation is given by

[TABLE]

Calculating the perturbations to the material frame of the rod as

[TABLE]

we find

[TABLE]

The perturbed tangent vector is given by

[TABLE]

which gives the deformed shape as

[TABLE]

This has the same quadratic bending away from the $z$-axis as the helix axis in Eq. 124 and captures the bending of the axis in the direction perpendicular to the bending torque exerted at its clamped end and also the quadratic terms in Eq. 124 for s such that $({R}_{0}\ll s\ll L)$. This bending of the helix axis is illustrated in Fig. 6.

III.2 Feedback of the deformation on the kinematics

III.2.1 Scalings

Having computed the deformation of the helical filament, we now proceed to investigate its feedback on the locomotion kinematics. Since we assumed that the helix was long enough so as not to wobble, we may take the velocity and rotation rate to be axial. We also assume that the perturbed velocity and rotation rate are along the long axis. Although there are non-axial deformations due to the bending of the axis, they average out to zero upon a complete rotation of the swimmer.

In order to find the expected order of magnitudes for the perturbed quantities, we consider the following scaling arguments. From Eqs. 93 and 94 we obtain the scalings $A_{x},A_{z}\sim\zeta_{\perp}{U}^{(0)}\sim\zeta_{\perp}{\omega_{fil}}^{(0)}{R}_{0}$. The leading-order kinematics in Eq. 103 give ${U}^{(0)}\sim hM/\zeta_{\perp}{R}_{0}L$ where $M$ stands $M_{mag}$ or $M_{mot}$ for artificial motors and bacteria respectively. The small rotations of the material frame due to its deformation defined by Eq. 20, scale as $\delta\phi\sim\mathcal{P}$ where $\mathcal{P}$ is a small dimensionless number

[TABLE]

Such scaling can be seen in Eq. 168 of Appendix A for example. Thus we have $\delta\mathbf{d}_{3}\sim\mathcal{P}$ and the deformation scales as $\delta\mathbf{r}\sim\mathcal{P}L$. We note that the bending of the axis is quantified by an angle $\beta$ between the end of the bent axis and the $z$ axis, where $\tan\beta\sim|\delta\mathbf{r}_{xy~{}proj}|/L\cos{\alpha}_{0}$, hence $\beta\sim\mathcal{P}$. The effect of the bending of the axis to the kinematics will therefore be an ${\mathcal{O}}(\mathcal{P})$ velocity component perpendicular to the $z$ axis that, due to rotation of the swimmer, will average to 0 at ${\mathcal{O}}(\mathcal{P})$, and similarly for ${\omega_{fil}}^{(0)}$. It is thus appropriate to project the force and torque balances along the $z$-axis and assume that $\delta\mathbf{U}$ and $\delta\boldsymbol{\omega}$ are aligned with the $z$-axis. As $\delta\mathbf{d}_{3}\sim\mathcal{P}$, from Eqs. 66, 67, 79 and 80 we therefore expect $\delta U$ and $\delta\omega$ to be of order ${\mathcal{O}}(\mathcal{P})$.

III.2.2 Perturbing the resistance matrices

The leading-order kinematics in $\mathcal{P}$ are given in § III.1.3. In this section we perform an asymptotic analysis for the next-order kinematics, i.e. the feedback of the deformation on the kinematics. We start by calculating the projections along the $z$-axis of the perturbations to the resistance matrices

[TABLE]

We may write these results as linear combinations of four integrals, $J_{i}$ ($1\leq i\leq 4$),

[TABLE]

where the four integrals are given by

[TABLE]

which depend on the deformation integrated along the entire length of the helix. The details of this long calculation are given in Appendices B and C where the cases of either clamped or free $s=0$ ends are addressed, as well as the effects of the viscous rotational torque.

In the case of a clamped end at $s=0$ ($\boldsymbol{\Phi},\Delta\mathbf{r}_{0}=\mathbf{0}$), the perturbations to the resistance matrices projected along the $z$ axis are

[TABLE]

Notably, the viscous rotational torque, $\mathbf{N}_{visc}$, does not contribute to any of the perturbations of the resistance matrices.

III.2.3 Unified Approach

We once again invoke the unified approach of section §III.1.3. The leading-order kinematics are given by Eqs. 103 and 104 and the perturbation to the swimming velocity is given by

[TABLE]

where $M$ stands for either $M_{mag}$ or $M_{mot}$ for artificial motors and bacteria respectively, and where $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{D}_{*}$ are given by Eqs. 100, 101 and 105 respectively. The denominator $\mathcal{A}\mathcal{D}_{*}-\mathcal{B}^{2}$ is given by Eq. 106. After substituting in the values of the resistance matrices, the quantities $A_{x},A_{z}$ become

[TABLE]

Considering only the clamped contribution, the perturbations to the resistance matrices in Eqs. 148-151 which are linear in $A_{x}$ and $A_{z}$ can be calculated (with details shown in Appendix D). The perturbation $\delta U$ is finally obtained as

[TABLE]

where $\mathcal{A}\mathcal{D}_{*}-\mathcal{B}^{2}$ is given by Eq. 106 and

[TABLE]

with relative errors of order ${\mathcal{O}}\left({R}_{0}/L\right)$ and where we have used the shorthand notation $s\equiv\sin{\alpha}_{0}$ and $c\equiv\cos{\alpha}_{0}$. The total velocity, $U_{tot}$, is simply $U_{tot}={U}^{(0)}+\delta U$. The details of the calculation are given in Appendix D.

We now proceed by considering the applications of this unified approach to artificial motors and bacteria separately.

III.2.4 Artificial motors

A stiff helical microswimmer.

Typical profiles of the dimensionless velocities ${U}^{(0)}/\hat{U}$ in the reference configuration and $U_{tot}/\hat{U}$ after the deformation (left) and the perturbation $\delta U/\hat{U}$ due to the deformation (right) are shown in Fig. 7 as a function of the reference helix angle ${\alpha}_{0}$ (the choice of speed scale $\hat{U}$ is given below). These results are obtained using Eqs. 103 and 152 respectively, using the full expressions for the resistance matrices (within the limit ${R}_{0}/L\ll 1$), for a RH artificial bacterial flagellum with the choice $\mathcal{P}\equiv M_{mag}L/EI=0.035$. For the non-dimensionalisation we use the length scale $L$, time scale $T=\mu L^{3}/M_{mag}$ and force scale $\mu L^{2}/T=M_{mag}/L$ (or, equivalently, mass scale $\mu LT=\mu^{2}L^{4}/M_{mag}$), so that the scale for the speed is $\hat{U}=M_{mag}/\mu L^{2}$. We use values for the geometrical parameters similar in order of magnitude to those of the rigid ‘artificial bacterial flagellum’ of Ref. Zhang2009 with the choice $n=8$ for the number of turns so that we can satisfy the requirement that ${R}_{0}/L\ll 1$. As such, our model microswimmer has $a_{h}=2~{}\mu\rm{m}$, $r=50~{}\rm{nm}$, $\Lambda=5~{}\mu\rm{m}$, $n=8$, $L=40~{}\mu\rm{m}$, $M_{mag}=4.3~{}\times 10^{-17}~{}\rm{Nm}$ and is moving in water with dynamic viscosity $\mu=10^{-3}~{}\rm{Pas}$. As a result, the scale for speed, $\hat{U}$, takes the value of $27~{}\mu\rm{m~{}s}^{-1}$. We set the Young’s modulus to be $E=10^{10}~{}\rm{Pa}$, as an indicative example of a stiff helical microswimmer, according to the ‘equivalent’ Young’s modulus of nanowire filaments of existing microswimmer designs Paknanowire2011 . This value is also relevant for ZnO nanohelices GaoMaiWang2006 .

Importantly, the perturbation in velocity, $\delta U$, changes sign at an angle ${\alpha}_{0}^{*}\approx 0.60~{}\rm{rad}~{}(34.4^{\circ})$, being negative for smaller values of ${\alpha}_{0}$ and positive for larger ones. Note that inside the resistance matrices used in Eq. 152 we have kept all terms involving $a_{h}$ and $\zeta_{r}$, and have not imposed the ratio $a_{h}/L$ to be small.

Approximation.

We proceed by further approximating the expression for $\mathcal{Q}$ in Eq. 157, assuming that both the helical radius and the head are negligible in size compared to the total contour length of the helix, such that ${R}_{0}$ and $a_{h}$ satisfy the limits

[TABLE]

The resulting expression for $\delta U$ in Eq. 156 simplifies to

[TABLE]

where $\mathcal{A}\mathcal{D}-\mathcal{B}^{2}$ is given by Eq. 109 and

[TABLE]

When reducing the ratios ${R}_{0}/L$ and $a_{h}/L$, for non-vanishing helix angles for which the ratio ${R}_{0}/L$ is kept small, the profile of $\delta U$ vs. ${\alpha}_{0}$ that we obtain from the full expression of Eq. 152 converges to the theoretical approximation of Eq. 159. Indeed, if we keep all parameter values as above and only decrease the size of the head, the angle $\alpha_{0}^{*}$ at which $\delta U$ vanishes converges to $33.1^{\circ}$, which is well approximated by the root of $\Gamma({\alpha}_{0})$ in Eq. 160 equal to $32.8^{\circ}$.

For non-vanishing ${\alpha}_{0},{R}_{0}$,

[TABLE]

with a relative error of order of magnitude ${\mathcal{O}}\left({R}_{0}/L,~{}a_{h}/L,~{}a_{h}^{3}/{R}_{0}^{3}L,~{}a_{h}^{4}/{R}_{0}^{2}L^{2}\right)$, so that the corrected velocity is

[TABLE]

with the dimensionless parameter $\epsilon$ defined as

[TABLE]

where $\mathcal{P}$ is the parameter defined in §III.2.1. We therefore see that a stiff elastic helix has a quadratic perturbation in $M_{mag}$ to the classical linear relation between $U$ and $M_{mag}$ for a rigid helix. This is of course valid for a weak enough actuation, or a stiff enough helix, that $\epsilon$ is kept small.

In our numerical results, $\zeta_{\perp}$ and $\rho$ are given by Eqs. 9 and 10. For the purpose of simplicity and to allow physical interpretation, we can further assume an integer number of turns for the helix and take the approximate value $\rho=0.5$ for the ratio of drag coefficients, so that the expressions for $\delta U$ and the corrected velocity $U_{tot}$ further simplify to

[TABLE]

Feedback to the speed of propulsion.

We may now investigate whether the small amount of elasticity speeds up or slows down a stiff elastic helical filament compared to its rigid kinematics, which is governed by the sign of $\delta U$. We first note that the sign of ${U}^{(0)}$ is opposite for right-handed (RH) and left-handed (LH) helices, as shown by the presence of the chirality index $h$ in its expression. In contrast, the expression for $\delta U$ in Eqs. 156 and 157 does not involve the value of $h$. Thus, if all other geometrical parameters are kept constant, the magnitude of the speed for a stiff helix, $|{U}^{(0)}+\delta U|$, will be increased or decreased according to the handedness of the helix. In addition, as discussed above, the sign of $\delta U$ changes at a critical angle ${\alpha}_{0}^{*}$ (with a value which depends on the geometry). The case of our typical stiff helical swimmer above had ${\alpha}_{0}^{*}=34.4^{\circ}$. For a microswimmer with a smaller head but otherwise same geometrical parameters, we obtained instead ${\alpha}_{0}^{*}=33.1^{\circ}$. Thus the perturbation to the speed changes sign according to a subtle interplay between the handedness and the reference helix angle.

We can interpret the change in the sign of $\delta U$ by approximating the new, steady-state, perturbed shape by an effective helix over which the extension or compression has been uniformly distributed. As illustrated in the top panel of Fig. 8, the extension (or compression) of a RH (or LH) helix can be interpreted as a decrease (or increase) in the helix angle. This provides us with an intuitive reasoning for the sign of the perturbation to the speed.

Consider the rigid speed profile as a function of the reference helix angle. As illustrated in the bottom panel of Fig. 8, if we start at an angle greater than the optimum, a RH helix that extends will reduce its effective helix angle, thereby moving towards the maximum of the speed profile and as a result will speed up. In contrast, a LH helix with a reference helix angle greater than the optimum will move away from the optimum, as its effective helix angle increases due to the compression. If we start at an angle less than the optimum, a RH helix that extends will reduce its effective helix angle, thereby moving away from the maximum of the speed profile, and will thus slow down. In contrast, a LH helix with a reference helix angle less than the optimum will move towards the maximum and speed up, as its effective helix angle increases due to the compression. We think that the quantitive discrepancy between the angle for which $\delta U$ vanishes and the optimum value of the uniform rigid helix velocity profile is likely due to the nonuniformity of the extension/compression along the filament and the bending of the helix axis.

III.2.5 Bacteria

Swimming bacteria will have, in principle, a contribution in their speed from the values of $\boldsymbol{\Phi},\Delta\mathbf{r}_{0}$ (in Eq. 34) as the end of the flagellar filament connected to the hook. This is due to the flexibility of the hook at the base of the flagellum. In fact, for some single-flagellated bacteria such as V. alginolyticus, the buckling of the hook after the transition from a backward to a forward swimming period is a mechanism to change the angle between the cell body and the flagellar filament axis SonGuastoStockeronHook2013 . Our study focuses however on the swimming along straight lines and hence for our purposes we will neglect buckling of the hook and will assume that both $\boldsymbol{\Phi}$ and $\Delta\mathbf{r}_{0}$ are zero.

In the case of bacteria, we can make use of Eqs. 156 and 157, where the indicator function now takes the value [math] and thus the terms multiplied by it are removed. Typical profiles of the dimensionless velocities in the reference configuration, ${U}^{(0)}/\hat{U}$, and after the deformation, $U_{tot}/\hat{U}$, are shown in Fig. 9 (left) as a function of the reference helix angle, ${\alpha}_{0}$, for a RH bacterial flagellar filament with $\mathcal{P}\equiv M_{mot}L/EI=0.19$ and using the full expressions for the resistance matrices (within the limit ${R}_{0}/L\ll 1$) (Eqs. 112). We also plot in Fig. 9 (right) the perturbation in swimming speed due to the deformation, $\delta U/\hat{U}$, as a function of ${\alpha}_{0}$. Here the non-dimensionalisation used the velocity scale $\hat{U}=M_{mot}/\mu L^{2}$. The parameter values are chosen as the typical values for Salmonella bacteria given in Refs. NambaVonderviszt1997 ; KamiyaAsakura1976 (and references therein), so that we have $a_{h}=1~{}\mu\rm{m}$, $r=11.5~{}\rm{nm}$, $\Lambda=2.3~{}\mu\rm{m}$, $n=6$, $L=13.8~{}\mu\rm{m}$, $M_{mot}=2~{}\times 10^{-18}~{}\rm{Nm}$ and take the viscosity to be that of water, $\mu=10^{-3}~{}\rm{Pas}$. The scale for speed, $\hat{U}$, takes as a result the value of $10.5~{}\mu\rm{m~{}s}^{-1}$. Further we take the Young’s modulus to be $E=1.04\times 10^{10}~{}\rm{Pa}$ as in Ref. TrachtenbergHammel1992 so that the bending stiffness is $EI=1.43\times 10^{-22}~{}\rm{Nm^{2}}$. We obtain as a result that the perturbation in swimming speed, $\delta\hat{U}$, transitions in sign at the angle ${\alpha}_{0}\approx 0.69~{}\rm{rad}~{}(39.5^{\circ})$, being negative for smaller values of ${\alpha}_{0}$ and positive for larger ones.

Importantly, the reported values for the elasticity of flagellar filaments in the literature show large variations, which naturally leads to a wide range of values for the bending stiffness, $EI$ KimPowers2005HelixDeformation . Trachtenberg and Hammel reported values of the Young’s modulus ranging from $E=1.04\times 10^{10}~{}\rm{Pa}$ to $E=1.77\times 10^{11}\rm{Pa}$ for a variety of filament types TrachtenbergHammel1992 , which correspond to bending stiffness of $EI=10^{-22}~{}\rm{Nm^{2}}$ and $2\times 10^{-21}~{}\rm{Nm^{2}}$, and with the parameter values above, the dimensionless number $\mathcal{P}$ would range from $\mathcal{P}=0.19$ to $\mathcal{P}=0.011$. Other reported values of the filament elasticity include $(2-4)\times 10^{-24}~{}\rm{Nm^{2}}$ for the bending stiffness $EI$ of Salmonella typhimurium filaments FujimeMaruyamaAsakura1972 and $10^{10}~{}\rm{N/m^{2}}$ for the shear modulus of S. typhimurium flagellar filament which gives the value $\mu_{S}J=10^{-22}~{}\rm{Nm^{2}}$ for the twist modulus of a filament of radius of $10~{}\rm{nm}$ HoshikawaKamiya1985 . Takano et al. estimated $EI=10^{-24}~{}\rm{Nm^{2}}$ for Vibrio alginolyticus Takanoetal2003HelicFlagel whereas Kim and Powers estimated $EI=10^{-24}~{}\rm{Nm^{2}}$ for S. typhimurium KimPowers2005HelixDeformation by reinterpreting the data of Hoshikawa and Kamiya HoshikawaKamiya1985 .

Our analysis mainly focuses on the asymptotics and feedback of the deformation to the kinematics, and to first order, the perturbation $\delta U$ scales linearly with $1/EI$. In particular, we note that within the ‘stiff’ filament regime and first-order analysis presented here, the value of the reference angle ${\alpha}_{0}$ for which the sign of $\delta U$ changes is independent of the Young’s modulus of the filament.

The resulting expression for $\delta U$ in Eq. 156 simplifies to

[TABLE]

where $\mathcal{A}\mathcal{D}_{fil}-\mathcal{B}_{fil}^{2}$ and $\Gamma({\alpha}_{0})$ are given by Eqs. 114 and 160 respectively. For non-vanishing ${\alpha}_{0},{R}_{0}$, in the limit of small cell body and helical radius compared to the length of the flagellar filament, the expression for $\delta U$ reduces to

[TABLE]

with a relative error of order of magnitude ${\mathcal{O}}\left({R}_{0}/L,~{}a_{h}/L,~{}a_{h}^{3}/{R}_{0}^{3}L,~{}a_{h}^{4}/{R}_{0}^{2}L^{2}\right)$, similarly to the result for artificial motors in Eq. 164. Thus, similar conclusions can be drawn as for the case of artificial motors.

IV Discussion

In this paper, motivated by the locomotion of bacteria and artificial microswimmers, we have developed the mathematical framework to study the steady-state motion of individual, stiff elastic filaments attached on a rigid body and propelling in a viscous fluid. The richness of this fluid-structure interaction problem comes from the hydrodynamic loads that induce deformation and whose integrated effect along the whole deformed shape leads to perturbed swimming kinematics.

Our analytical approach consisted in first integrating the hydrodynamic load along the initial reference configuration geometry and obtaining the induced deformation via the linear constitutive equation for the bending moment given by the classical Kirchhoff equations. For a helical geometry, the axis bends and the helix is extended or compressed if it is right- or left-handed, respectively. Our mathematical expressions describe the full three-dimensional deformation and an effective rod analogy provides an explanation for the bending of the helix axis whose origin was previously unclear. Imposing the force and torque balance along the new, deformed shape linearly relates the velocity perturbation to the deformation integrated along the entire length of the filament. As a result, the propulsion speed acquires a quadratic correction in the actuation.

Our analytical expression for the correction to the speed for stiff elastic helical filaments changes sign according to a subtle interplay between the handedness and reference helix angle. Approximating the new, steady-state, perturbed shape as an ‘effective’ uniform helix, the extension (or compression) of a right- (or left-) handed helix can be rationalised as due to a decrease (or increase) in the helix angle. This interpretation provides us with an intuitive reasoning for the sign of the speed perturbation. If the value of the original helix angle in the reference configuration is greater than the optimal angle maximising the speed profile of a rigid helix as a function of the helix angle, then the shift in the effective helix angle due to extension (or compression) will give rise to increased (or decreased) speeds depending on the handedness of the helix.

We note that our mathematical formulation allows to tackle any slender filament shape, and need not be limited to helices. Furthermore, one could include in it a model for the activity of the filament. A natural extension of this type would address magnetised tails such as the recently-manufactured hydrogels with embedded magnetic nanoparticles which not only guide self-folding during fabrication Huangetal2016Nature but could also display interesting non-linear relationship between the external actuation and the propulsion speed.

Acknowledgements

This project has received funding from the EPSRC (PK) and the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement 682754 to EL).

Appendix A Calculation of the Deformation for helical filaments

The small rotations of the material frame due to its deformation (defined by Eq. 20) are calculated as the integral of the bending moment (given in Eq. 96), in accordance to Eq. 31,

[TABLE]

The perturbation to the tangent vector is given by Eq. 20 as $\delta\mathbf{d}_{3}=\delta\boldsymbol{\phi}\wedge{\mathbf{d}_{3}}^{(0)}$, which gives

[TABLE]

The deformation of the helix is given be Eq. 33, as the integral of the perturbation to the tangent vector.

Integrating the expression for $\delta\mathbf{d}_{3}(s)$ in Eq. 169 gives

[TABLE]

Appendix B Clamped $s=0$ end-point contribution

In this appendix, we calculate the leading-order estimates of the quantities $J_{1}^{clamp}-J_{4}^{clamp}$ for our helical geometry, clamped at one end.

B.1 Calculation of $J_{1}^{clamp}-J_{4}^{clamp}$

We first calculate the contributions to $J_{1}^{clamp}-J_{4}^{clamp}$, without the terms due to the viscous rotational torque, as those will be calculated separately in the next subsection. Straightforward integrations give

[TABLE]

We can rewrite $J_{4}^{clamp}$ and $J_{2}^{clamp}$ as

[TABLE]

We thus need to calculate the two integrals $I_{1}=\int_{0}^{L}\left[\cos\left(2\pi s/\Lambda_{0}\right)\delta x+h\sin\left(2\pi s/\Lambda_{0}\right)\delta y\right]\;\mathrm{d}s$ and $I_{2}=\int_{0}^{L}\left[\sin\left(2\pi s/\Lambda_{0}\right)(\delta\mathbf{d}_{3})_{x}-h\cos\left(2\pi s/\Lambda_{0}\right)(\delta\mathbf{d}_{3})_{y}\right]\;\mathrm{d}s$. Straightforward integrations lead to

[TABLE]

Using integration by parts, we obtain

[TABLE]

The boundary terms are

[TABLE]

Thus

[TABLE]

and

[TABLE]

Putting it all together we obtain

[TABLE]

with the terms arising from the viscous rotational torque to be calculated in the next subsection.

B.2 Contribution from the viscous rotational torque

In this subsection we calculate the contributions $J_{1}^{rot}-J_{4}^{rot}$ arising from the viscous rotational torque to $J_{1}^{clamp}-J_{4}^{clamp}$.

[TABLE]

where

[TABLE]

and

[TABLE]

The calculation to obtain the result of Eq. 191 is the following:

[TABLE]

Hence

[TABLE]

Integrating by parts gives

[TABLE]

Hence

[TABLE]

with the last equality arising by substitution of the formulae for $\delta x|_{L},\delta y|_{L}$ terms arising from the viscous rotational torque.

[TABLE]

In summary, $J_{1}^{rot}=0$, $J_{2}^{rot}=0$ and $J_{4}^{rot}=\cos{\alpha}_{0}J_{3}^{rot}$, with

[TABLE]

These values of the contributions to the $J_{1}-J_{4}$ integrals, $J_{1}^{rot}=0,~{}J_{2}^{rot}=0,~{}J_{4}^{rot}=\cos{\alpha}_{0}J_{3}^{rot}$, mean that the viscous rotational torque does not contribute to the perturbations to the resistance matrices,

[TABLE]

We summarize the resulting $J_{1}-J_{4}$ for the clamped case, ($\boldsymbol{\Phi},\Delta\mathbf{r}_{0}=\mathbf{0}$), including the terms arising from the viscous rotational torque

[TABLE]

The resulting perturbations to the resistance matrices (projected along the z axis) are

[TABLE]

where we note again that the viscous rotational torque, $\mathbf{N}_{visc}$, does not contribute to any of the perturbations of the resistance matrices.

Appendix C Free $s=0$ end-point contribution

In this Appendix we calculate the contributions of $\boldsymbol{\Phi},\Delta\mathbf{r}_{0}$ to the perturbations to the resistance matrices. We first find the contributions to $J_{1}^{free}-J_{4}^{free}$ from $\boldsymbol{\Phi},\Delta\mathbf{r}_{0}$ for the case of free $s=0$ end-point, i.e. with $\boldsymbol{\Phi},\Delta\mathbf{r}_{0}\neq\mathbf{0}$. The calculation is as follows:

[TABLE]

[TABLE]

[TABLE]

With the contributions $J_{1}^{free}-J_{4}^{free}$ of $\Phi,\Delta{\mathbf{r}}_{0}$ to $J_{1}-J_{4}$ calculated, we now proceed to find the contribution of $\Phi,\Delta{\mathbf{r}}_{0}$ to the perturbations to the resistance martices.

[TABLE]

[TABLE]

In summary, we have

[TABLE]

and the resulting contributions to the perturbations of the resistance matrices from $\boldsymbol{\Phi},\Delta\mathbf{r}_{0}$ are

[TABLE]

Appendix D Feedback to the kinematics

We note that throughout this Appendix, we will use the shorthand notation $s=\sin{\alpha}_{0}$ and $c=\cos{\alpha}_{0}$. Using the rigid, leading order kinematics of Eqs. 103-104, and Eqs. 93-94, we can express $A_{x}$ and $A_{z}$ as

[TABLE]

Substituting into Eqs. 148-151 (and using the simplified notation $\delta\mathcal{A}$ instead of $\delta\mathcal{A}_{zz}^{clamped}$ etc) we obtain

[TABLE]

and substituting into Eq. 152 gives the result of Eqs. 156,157.