# The role of body flexibility in stroke enhancements for finite-length undulatory swimmers in viscoelastic fluids

###### Abstract.

The role of passive body dynamics on the kinematics of swimming micro-organisms in complex fluids is investigated. Asymptotic analysis of small amplitude motions of a finite-length undulatory swimmer in a Stokes-Oldroyd-B fluid is used to predict shape changes that result as body elasticity and fluid elasticity are varied. Results from the analysis are compared with numerical simulations, and the numerically simulated shape changes agree with the analysis at both small and large amplitudes, even for strongly elastic flows. We compute a stroke-induced swimming speed that accounts for the shape changes, but not additional effects of fluid elasticity. Elasticity-induced shape changes lead to larger amplitude strokes for sufficiently soft swimmers in a viscoelastic fluid, and these stroke boosts can lead to swimming speed-ups. However, for the strokes we examine, we find that additional effects of fluid elasticity generically result in a slow-down. Our high amplitude strokes in strongly elastic flows lead to a qualitatively different regime in which highly concentrated elastic stresses accumulate near swimmer bodies and dramatic slow-downs are seen.

## 1. Introduction

There has been an intense effort over the past 10 years to understand the effect of fluid elasticity on micro-organism swimming. Experiments, analysis, and simulations of low-Reynolds number swimming of microorganisms in complex fluids, in particular viscoelastic fluids, has led to a variety of results – some complimentary, some apparently conflicting – on the effect of fluid elasticity on swimming speed. We know that gait, body stiffness, and nonlinear effects matter, but we still do not have a clear understanding of how they interact during locomotion.

Early work quantifying the effect of fluid elasticity on swimming using a linear constitutive law for the fluid and asymptotic analysis of small amplitude motions showed that elasticity had no effect on swimming speed but increased swimming efficiency [1, 2]. However, in [3] a full analysis of the classical Taylor swimming sheet for small amplitude undulatory motion showed that the nonlinearities of the viscoelastic fluid model must be included in a computation of swimming speed, and found that swimming speed is always hindered by fluid elasticity. Similar small amplitude asymptotic analysis was done for waving filaments and helices [4, 5] also predicting slow-downs due to fluid elasticity. [6] and [7] have demonstrated the importance of the details of the swimming gait in understanding the effect of fluid elasticity on swimming speed, indeed showing that elastic speed-ups are possible for some gaits, which is further highlighted by the analysis of three-sphere swimmers [8].

Biological swimmers have been shown to change their gait in response to changes in rheology [9, 10, 11], making it hard to interpret the mechanisms responsible for observed changes in swimming performance. In more controlled physical models of swimmers in different fluids a variety of results have shown that fluid elasticity can boost swimming speed [12, 13, 14] or retard swimming speed [15, 16]. Swimmers with large amplitude motions have been theoretically investigated using numerical simulations, and have added significant information about the response of swimmers to fluid elasticity with a variety of swimming gaits [17, 18, 19, 20, 21, 22, 23, 24, 25]. In addition, for recent reviews of swimming in complex fluids, see [26] for a theoretical view, and see [27] for an experimental view. These many studies have focused on different types of swimmers, in different fluid rheologies, and despite the wealth of results we still lack an understanding of the underlying principles of swimming in complex fluids.

To try to isolate physical mechanisms that are significant in a variety of biologically relevant problems, but simple enough to analyse, we focus here on undulatory swimmers in an Oldroyd-B fluid. Even in this more restrictive setting, we nevertheless still find apparently contradictory results and a lack of mechanistic explanations for those differences. Asymptotic analysis of infinitely long, prescribed shape, small-amplitude swimmers has shown that swimming is hindered by the addition of elastic stresses [3], although allowing for flexibility can lead to enhancements [28]. Biological experiments have shown a viscoelastic slow-down for C. elegans [9], while simulations of finite-length swimmers with large tail amplitudes [17, 22] give a non-monotonic boost as fluid elasticity is varied. In [22] we concluded that shape changes due to body flexibility and fluid elasticity are important, but those results did not explain the results from a physical experiment which showed monotonic speed-ups due to fluid elasticity in swimmers with large tail amplitudes [14]. Furthermore, recent numerical simulations [25] appear to contradict the speed-ups reported in [17, 22].

The relevance of body elasticity in viscoelastic speed enhancements was identified for small amplitude infinite length swimmers in [28], where the authors attribute the speed enhancements to a viscoelastic “suction” which results in an amplitude boost. However, their analysis does not extend to finite-length large amplitude swimmers where the role of elasticity-induced shape changes has not been addressed directly. The disparity of the results in [17, 22, 14, 25], all focusing on large amplitude, finite-length, undulatory swimmers in Oldroyd-B fluids, indicates that something is missing in our understanding of the problem.

There remains a gap between our understanding from analysis and what we see in computational, biological, and physical experiments. Here we combine analysis with numerical simulations of finite-length large amplitude swimmers to show how fluid elasticity induces shape changes in finite-length flexible swimmers and how those shape changes can lead to speed boosts. We show how shape changes depend on both body stiffness as well as fluid elasticity and analyse the effect that shape changes alone have on swimming speed.

## 2. Effect of passive body dynamics

### 2.1. Methodology

We follow the computational framework in [22, 29], where the swimmer is modeled as an inextensible flexible sheet of finite-length immersed in a 2D fluid. We describe the undulatory motion of the swimmer by a curvature of the form

(1) |

where is the body coordinate. Here is the curvature amplitude at the “tail” () and is the curvature amplitude at the “head” () of the swimmer.

We use the immersed boundary method to solve for the coupled motion of the fluid and the swimmer [30]. Both inextensibility and shape are imposed (approximately) by forces that are designed to penalize extension and deviation from the prescribed target curvature. These forces are derived from the variation of bending and extension (stretching) energy functionals. For example, the bending energy is

(2) |

where B is the bending stiffness, is the curvature of the sheet, and is the prescribed target curvature. One can interpret the model as an active sheet with bending stiffness B driven by an active body moment density . We scale forces relative to viscous forces so that for , the realized shape of the swimmer is very close to the prescribed shape. For , the realized shape is the result of fluid-structure interaction; i.e. passive body dynamics influence the resulting stroke.

The viscoelastic fluid is described by the Oldroyd-B model at zero Reynolds number [31], regularized by stress diffusion [32, 33]. The system of equations describing the fluid are

(3) | |||

(4) | |||

(5) |

where u is the fluid velocity, is the pressure, is the viscoelastic stress, is the rate of strain tensor, and is the elastic force density generated by the swimmer. Here is the polymer to solvent viscosity ratio, the Deborah number is the ratio of elastic relaxation time to stroke period, and is the stress diffusion coefficient.

### 2.2. Varying body stiffness

To understand the role of body elasticity, we use our simulations to calculate the Stokes-normalized swimming speed while varying B and De for a fixed period (). We use a stroke defined by equation (1) with and This gives a high-amplitude stroke like in [17, 22]. In figure 1(a) we plot normalized swimming speed as a function of De for three characteristic stiffness values of , which we refer to as very soft, moderately soft, and stiff, respectively. For very soft swimmers we see a monotonic boost in swimming speed, with a greater than 50 boost for high

This response is similar to what was reported in [14] using a physical model of a swimmer with a flexible tail (figure 1(b)). For moderately soft swimmers, we see a non-monotonic speed-up, including a smaller speed boost over the Newtonian speed, followed by a slow-down at larger This type of non-monotonic speed-up was first reported in [17] and again in [22] for a soft stroke with high amplitude (figure 1(c)). Finally, for stiff swimmers we see non-monotonic behavior but no boost over the Newtonian speed, again followed by a slow-down at larger This type of slow-down was reported in [22] for a stiff kicker (figure 1(d)).

In contrast to stiff, or rigid, swimmers, the dynamics of flexible swimmers involve an additional time scale. In a viscous fluid, rigid swimmers move with a velocity proportional to the beat frequency (the only time scale in the problem). The problem of a rigid swimmer in a viscoelastic fluid has two time scales, the beat frequency and the relaxation time, whose ratio is the dimensionless relaxation time The swimming speed of soft swimmers depends nonlinearly on the frequency because the shape changes with the frequency. Figure 2 shows the swimming speed in a Newtonian fluid for both a stiff swimmer, and a moderately soft swimmer over a range of beating frequencies. We see a linear response in the case of a stiff swimmer and a nonlinear response for the moderately soft swimmer. Inset swimmer shapes show how the shape changes as a function of the stiffness in the high frequency case.

To illustrate the significance of multiple time scales for flexible swimmers in viscoelastic fluids we compute the Stokes-normalized swimming speed as a function of De varied two ways: by varying the relaxation time for a fixed period, and by varying the period for a fixed relaxation time. Both simulations are performed with the same bending stiffness, , where passive body dynamics are significant. Results are shown in figure 3 (a) for a swimmer with the same stroke from figure 1, and the two curves show remarkable qualitative differences. For a rigid swimmer these would give equivalent results. Thus this third time-scale, arising from body flexibility, needs to be explicitly included in any discussion of swimming in elastic fluids. A more complete picture of how the swimming speed depends on both the relaxation time and period when the body is soft, is shown in figure 3 (b). Contours of constant are overlayed in black and the effect of body stiffness is clearly evident as you see the swimming speed vary significantly along any of the contours. The dashed lines denote the locations of the data in figure 3 (a).

## 3. Analysis of shape changes

The effect of body stiffness on swimming kinematics has been previously studied for viscous fluids [34, 35]. Shape changes in viscoelastic fluids have been examined [36], but the relationship between shape changes and swimming speed has not been examined for finite-length swimmers. Here we review the theory and compare it with numerical simulations.

### 3.1. Linear theory: Newtonian fluids

We begin by considering small amplitude displacements of a finite-length thin elastic rod in a Newtonian fluid driven by prescribed curvature, , (equivalently, prescribed moments) along the body with free ends. The theoretical analysis is similar in 2D and 3D, however we will make note of the differences when we compare with the numerical simulations in 2D. We proceed with the analysis in 3D for simplicity. The shape of the rod is determined by the balance between elastic forces and viscous drag. The vertical displacement, satisfies

(6) | |||

(7) |

Here is the perpendicular drag coefficient and is the bending stiffness of the rod. Note that the use of new notation is intended to distinguish this (dimensional) bending stiffness we use in the linear theory from our previously defined (non-dimensional) bending stiffness, B which we use in our numerical simulations.

Nondimensionalizing equation (6) using the body length and the period of the driving force results in the dimensionless parameter we call the body response time:

(8) |

We note that could be called a body relaxation time. The same nondimensional group has appeared previously, but has been interpreted differently. In [37] a quantity similar to G was considered an “effective viscosity” of growing elastic filaments. The Sperm number () is the ratio of the body length to the viscous decay length [34, 4]. This interpretation is natural when considering filaments driven at one end rather than along the body as we do here.

We change variables from displacement to curvature deviation

to facilitate comparing with large amplitude simulations. For small displacements , and equation (6), (non-dimensionalized) becomes

(9) | |||

(10) |

For a given we use an orthogonal function expansion to solve the non-dimensional equations for . We let the driving curvature be given as

(11) |

and solve the eigenvalue problem,

for eigenvalues and eigenfunctions The expansion coefficients of the realized curvature, , are then

(12) |

From this solution we can see that as the rod is stiffened (), the resultant curvature tends to the prescribed curvature, We also see that for softer rods, i.e. smaller values of the body response time G, the amplitude of the curvature decreases and there is a phase lag relative to the prescribed shape.

As mentioned above, we use intrinsic coordinates and curvature deviations, to allow us to consider large prescribed curvatures. However we note that equation (9) lacks terms coming from geometric nonlinearities and inextensibility that may not be small when the prescribed curvature is large [38, 39]. In sections 3.3 and 4.2 we compare our simulations to theoretical analysis using (9), and in Appendix A we show that the influence of the additional terms is in fact small for the amplitudes we consider.

### 3.2. Linear theory: viscoelastic fluids

We can modify the linear theory for elastic rods to include fluid elasticity. This is similar to what was done in [40, 36]. In [40] modifications to linear rod theory to include linear viscoelastic fluid effects were presented, and the authors concluded that while fluid elasticity does not change swimming speed, it reduces total work and thus can boost efficiency. However, it was pointed out in [3] that it is essential to use a nonlinear elasticity model in these types of calculations because the swimming speed itself is second order in amplitude, where the nonlinear effects are relevant. We note that with these higher order terms [3] shows that swimming speed is always hindered by fluid elasticity for the case considered – infinite length low amplitude swimmer with sinusoidal undulations. In [36] the authors analyzed shape changes induced by fluid elasticity in a linearly elastic fluid. Unlike swimming speed, shape changes due to fluid elasticity come in to the asymptotic expansion at first order in amplitude, and hence it is reasonable to use a linearly elastic fluid to look at shape changes. [36] did not make conclusions about how these shape changes affect swimming speed. Here we perform a similar analysis as in [36], but by applying the analysis to deviations in curvature we are able to study shape changes in low and high amplitude finite length flexible rods. In Section 4 we discuss how these shape changes affect swimming speed.

As in equation (6) we can write a force balance relation between the force on a fluid and from the beam as

where the represents the normal force on the rod from the viscoelastic fluid. If we define the fluid force to be based on the total deviatoric stress then (upon linearization) using equation (5):

(13) |

where is the viscous drag force. Note that given the form of the system in equations (3)-(5), we have assumed a total viscosity of . The swimmer motion is time-periodic so we take the Fourier transform in time of equation (13) to solve for viscoelastic modifications to the fluid drag. This yields,

As in the viscous theory, we can solve for modifications to the curvature from body stiffness and use the modifications to the fluid drag to account for the fluid elasticity:

(14) |

The coefficients in equation (14) give an analytical expression for the modifications to the rod shapes relative to the prescribed shapes as fluid and body elasticity are varied.

### 3.3. Elastic shape changes: theory and numerical comparison

The analysis in the previous sections made use of resistive force theory which relates the drag force and velocity on a long thin cylindrical object. More generally, for small amplitude the vertical displacement satisfies

(15) |

where is the mobility operator and is the linearized bending force operator. Resistive force theory makes the approximation Our analysis of shapes (equation (14)) contains the quantity where G depends on the drag coefficient, (equation (8)). To use the more general linear theory in our analysis one can identify , where denotes the eigenvalue of the operator We relate to the dimensionless bending stiffness, B used in our simulations, through

To compare the linear analysis with our two-dimensional simulations we numerically approximate equation (15). For small deviations to the vertical displacement, is the integral operator which is the convolution of the vertical force with the fundamental solution to Stokes equations. We approximate using the method of regularized Stokeslets [41], which is a numerical method based on a regularized Greens function for the Stokes equations. We can also numerically approximate the eigenvalue of the bending force operator using a second-order finite difference, and we find that with point spacing the eigenvalues of are within 1% of the eigenvalues of the continuous operator. We give the eigenvalues for the first four nontrivial modes in table 1. Note that to compute we assume viscosity one. Also in table 1, we give the first four (mode dependent) drag coefficients computed as .

Mode () | |||
---|---|---|---|

1 | 32.55 | ||

2 | 42.94 | ||

3 | 54.69 | ||

4 | 66.84 |

In order to compare the predicted shape changes given by equation (14) with our numerical simulations we prescribe a curvature of the form

(16) |

in our model equations (3)–(1). The prescribed standing wave of constant curvature corresponds to a motion through circular arcs with peak curvature . By symmetry, this motion does not result in any horizontal translation of the body. We refer to these non-translating “swimmers” as flexors. We consider both low and high amplitude curvatures, and The shapes are shown inset in figure 4(a).

In figure 4 we plot the theoretical predictions from equation (14) (solid lines) along with values computed from numerical simulations; low amplitude () are indicated by hollow markers, and high amplitude () are indicated with filled markers. In figure 4(a) we plot the normalized amplitude of the first mode () to see how the amplitude deviates from the prescribed amplitude as a function of bending stiffness B. We see that generically the amplitude of the flexor decreases as the flexor is softened for fixed De. For sufficiently soft flexors () viscoelasticity increases the amplitude monotonically with De, but for stiffer swimmers the amplitude changes nonmonotonically with fluid elasticity.

In figures 4(b) and (c) we renormalize the data by the amplitude in a viscous fluid to see the effects of viscoelasticity more clearly. Again we see that fluid elasticity can increase the amplitude significantly for a soft flexor, but that effect is lost as the flexor is stiffened. When we plot the amplitude as a function of De for the very soft, moderately soft and stiff cases we see again that three qualitatively different regimes emerge. For very soft flexors the amplitude is monotonically increased by elasticity, for moderately soft flexors the response is non-monotonic, and can decrease or increase the amplitude, and for stiff flexors there is little change in the amplitude due to fluid elasticity. It is notable that the linear theory does such a good job predicting shape changes for low and high amplitude and for low and high Deborah number. In Appendix A we derive the theory for both the limit of small amplitude and the limit of high stiffness. We see in figures 4 (a) and (b) that the largest differences are for moderate stiffnesses at high amplitude. We note that we are showing results only for the first mode. For higher modes the trends are similar but the transition from stiff to soft behavior occurs at lower values of B because the eigenvalues increase with .

## 4. Analysis of swimming speed

In a viscous fluid, increasing the stroke amplitude will increase the swimming speed, and we can infer from section 3 that soft swimmers in a viscoelastic fluid sometimes obtain an amplitude boost over the corresponding swimmer in a Newtonian fluid. However when comparing swimmers in a viscoelastic fluid to those in a viscous fluid, even with an amplitude boost the viscoelastic swimmer may not swim faster than the viscous swimmer due to additional fluid elastic forces that the swimmer will encounter. Thus the effect of elasticity-induced shape changes is difficult to decouple from the overall effect of fluid elasticity. Analytical expressions for swimming speed can be obtained in certain limits, or for specialized swimmers, but even in these cases we see that the effect of fluid elasticity depends on many factors. For example infinite-length small amplitude undulatory swimmers show that a slow-down is generically expected for stiff swimmers in a viscoelastic fluid [3], but allowing for body flexibility, shape changes can lead to speed boosts [28].

In regimes that are more challenging for analysis such as the large amplitude, finite length swimmers considered here, it is more difficult to attribute speed boosts or slow-downs to specific swimmer attributes. For large amplitude finite-length undulatory swimmers, it was conjectured [17] that speed boosts were related to large tail stresses, and in [22] stroke asymmetries were correlated with both slow-downs and speed-ups. Here we will compute a stroke-induced swimming speed that isolates the effect of fluid elasticity on shape changes, and how those shape changes affect swimming speed in a Newtonian fluid. We then compare that analysis with the full nonlinear numerical simulations where the effect of shape changes is coupled with the fluid elasticity.

### 4.1. Swimming speed: two-mode swimmer

To keep the analysis simple we define a gait whose swimming speed in a viscous fluid we can compute analytically. We define a “two-mode swimmer” given by the curvature:

(17) |

where the for are the first and second bending modes. The modulation of a single mode results in a standing wave and will not translate in a Newtonian fluid. We use a sum of the first two modes with a phase difference to generate a nonreciprocal motion. Shapes of the first, second, and sum of the first and second modes are plotted in figure 5 for both low and high amplitudes.

Using resistive force theory one can derive the (time-averaged) swimming speed for a given, small amplitude, motion:

(18) |

where is the swimming speed, is the vertical displacement of the swimmer, and and are the perpendicular and parallel drag coefficients, respectively, [34, 35].

For small amplitudes, the shape of the swimmer (up to translation and rotation) is given by integrating equation (17) twice in space to compute the swimming speed via equation (18). The swimming speed (in a viscous fluid) for the two-mode swimmer is proportional to the product of the amplitudes and the sine of the phase difference:

(19) |

With this expression and the theoretical prediction for shape changes, we define a stroke-induced swimming speed which is the swimming speed in a Newtonian fluid that depends on the shape changes due to fluid elasticity and body flexibility. Specifically, for a given De and we compute from equation (14) ( and ) and the stroke-induced swimming speed from equation (19). We parametrize the shape changes due to changes in De using a parameter we call the stroke-Deborah number [22], In other words, represents the value of De used in equation (14) to compute the stroke-induced swimming speed via equation (19). Our analytical expression for the shape changes is based on a linearly elastic fluid, but because the nonlinear elastic effects and swimming speed are both second order in amplitude, we do not expect the stroke-induced swimming speed to capture the true viscoelastic swimming speed. An analytical expression for the swimming speed in a nonlinear viscoelastic fluid, as was computed in [3], is not tractable in the finite length case, because translational invariance, which facilitates the calculation for infinite length swimmers, is lost.

We plot the stroke-induced swimming speed for the two-mode swimmer over a range of and as with the flexor, we see the emergence of three regimes dependent on the body stiffness, see figure 6 (a). Shape changes boost the stroke-induced swimming speed if the swimmer is very soft, a smaller boost is obtained for the moderately soft swimmer, and additionally there is a non-monotonic response to increasing elasticity including a regime where shape changes slow down the swimmer, and finally if the swimmer is stiff there is a negligible effect.

### 4.2. Swimming speed: theory and numerical comparison

We simulate a two-mode swimmer of both low and high amplitude by prescribing a curvature of the form given in equation (17) with , , , for (low), and (high). These values come from projections of the stroke used to generate figure 1. The Stokes-normalized swimming speeds for a very soft, moderately soft, and stiff swimmer at both low (hollow markers) and high (filled markers) amplitude are plotted in figure 6 (b). We note that the three regimes seen in figure 6 (a) still emerge from these simulations, but for this two-mode swimmer the simulation swimming speeds are always slower than the stroke-induced swimming speeds. The ratio of swimming speed to stroke-induced swimming speed is shown in figure 6 (c). This quantity can be interpreted as the effect of fluid elasticity that is not related to shape changes. It is notable that these curves collapse onto a single curve for the low amplitude swimmers at all stiffnesses as well as the high-amplitude swimmer in the very soft regime. This additional elastic fluid effect on swimming speed is likely to be highly stroke dependent.

The additional effects of fluid elasticity are fundamentally different for the large amplitude, large De regime. In figure 6 (d) we plot the ratio of swimming speeds for the high-to-low amplitude strokes, and see that for (for sufficiently stiff swimmers) a significant difference in swimming speed arises. This difference is not related to shape changes because, like the flexor, the elasticity-induced shape changes predicted by the theory for the two-mode swimmer agree very well with the simulation results, for all at low and high amplitudes; see figure 7. At low amplitude (not shown) the relative error between theoretical and numerical predictions is less than 1 and for high amplitude the error at low De is at most 5 and at high De the error is at most 9. These results indicate that the theoretically predicted shape changes and their isolated effects on swimming speed can be well approximated by the analytical results for the amplitudes simulated and the range of De considered. A mechanistic understanding is lacking to explain what causes the dramatic slow-downs of swimmers in the high amplitude, high De regime.

We conjecture that the slow-downs in the high amplitude, high De regime must be attributed in part to the large localized stresses that accumulate near the body [17, 22]. To explore this conjecture, in figure 8 (a) we plot the average elastic to viscous stress ratio, computed as the time average over one period of where is the Frobenius norm, over a range of body stiffness B for both low (hollow markers) and high amplitude (filled markers) strokes for There is a notable transition in the stress ratio in the high De, high amplitude swimmer as the body is stiffened, while this stress ratio is flat for both low amplitude and low De swimmers. Stiff swimmer shapes along with the elastic-viscous stress ratio are plotted on a log-scale in figure 8 (b). The low amplitude strokes are surrounded by elastic stresses that are at least two orders of magnitude smaller than the high amplitude strokes, but even at large amplitude the low De swimmer still has relatively low stress near the body. Lastly, we plot the tail amplitude, as one measure of the swimmer stroke, in figure 8 (c). We see that for sufficiently soft swimmers the “high amplitude” stroke has a lower amplitude, which explains why in figure 6 the very soft high amplitude swimmer behaves like the low amplitude swimmers. For the high amplitude strokes it is only in the large amplitude and large De regime where significant stress accumulates near the swimmer.

## 5. Conclusions

In [22] we showed that stroke related speed-ups depend on body stiffness, and the analysis from this paper shows explicitly how the stroke changes depend on body stiffness and fluid elasticity through two dimensionless “relaxation times”: the fluid relaxation time, De and the body relaxation time, When we look at apparently contradictory results from the literature, we see that calculating G will determine which regime the swimmer falls into.

In
[14] the Sperm number is reported to be between
but even with the awareness that these are soft swimmers
the authors “conjecture that the effect [due to shape changes] is not
significant”. We use their reported
parameters^{1}^{1}1mm, cross-sectional radius m,
Young’s modulus GPa, viscosity Pa-s. For moment of
inertia we get
and thus and a characteristic
frequency of s and find We
cannot directly conclude that this value lies in the very soft regime () due
to differences between 2D and 3D as well as the way that the micro-swimmer is driven (it is a
flexible tail with a magnetically driven head). In Appendix B we repeat the calculation from section 4.1
to compute an equivalent stroke-induced swimming speed for a flexible filament driven at one end. This calculation shows that speed boosts still arise for sufficiently soft swimmers, despite the difference in driving mechanism. The most significant boost in speed from viscoelastic shape changes occurs for but for it would be reasonable
to conclude that the significant speed-ups observed in the experiment are related to shape changes.

In [25] the parameter reported for what they consider to be a soft swimmer is . However their swimmer length is mm (with characteristic length 1 mm) hence an equivalent dimensionless body response time G must be multiplied by . This pushes their “soft” simulations into the stiff regime where there are no speed-ups from shape changes, also agreeing with their results. Furthermore, in [25] it is conjectured that stress diffusion, used to regularize the simulations in [17, 22], is the source of the speed-ups, but the speed-ups we see are theoretically predicted, and realized in our simulations, even in the low amplitude regime where no regularization is necessary.

In our analysis we quantify the effect of body and fluid elasticity-induced shape changes on swimming speed. We see that the shape change analysis holds for the amplitudes simulated and the range of De considered, and in this case we see an additional elastic slow-down that is reminiscent of the type of slow-down predicted by asymptotic analysis of infinite-length small amplitude undulatory swimmers [3]. It may be tractable to apply asymptotic analysis [28, 6] to determine the form of the elastic slow-down for low amplitude finite-length swimmers. A fundamentally different regime arises for large amplitude swimmers in highly elastic fluids. A different approach is needed to understand the mechanisms that cause large localized stresses and their effect on swimming.

The authors would like to thank Henry Fu and Roberto Zenit for interesting discussions and suggestions on this work, and Michael Shelley for suggesting the term “flexors”. The authors would also like to thank the anonymous referees for suggesting useful modifications to our original manuscript. The work of RDG was partially supported by NSF grants DMS-1160438 and DMS-1226386.

## Appendix A Derivation of PDE for rod dynamics

In this appendix we give the derivation for the equation of motion for a thin filament in a viscous fluid which includes terms arising from inextensibility and geometric nonlinearties. For more details on similar calculations see for example [38, 39].

### a.1. Geometric Relations

Consider an inextensible thin rod whose centerline position is given by , where is arclength coordinate. We suppose that the deformation of the rod is planar. Let and be the tangent and normal vectors, be the tangent angle, and be the curvature. We have the following relationships between these quantities:

(20) |

### a.2. Equation of Motion

Let be the force density along the rod. Using resistive force theory, the motion of the rod is given by

(21) |

where and are the normal and tangential drag coefficients, respectively. Taking the derivative of this equation with respect to arclength gives

(22) |

Because , the left side of the above equation can be expressed as , and thus the normal terms give and the tangential terms must be zero:

(23) | |||

(24) |

Equation (24) represents a constraint on the forces that must be satisfied to maintain inextensibility. The evolution equation for the curvature is obtained by differentiating equation (23) with respect to arclength to obtain

(25) |

### a.3. Expression for Elastic Forces

The elastic forces are obtained from the variation of an elastic energy functional. The total elastic energy is the sum of a bending term from equation (2) with an energy associated with inextensibility:

(26) |

where is a tension used to enforce inextensibility. The force comes from the variation of the energy

(27) |

Using the natural free boundary conditions

(28) |

the force is

(29) |

We define the total tension in the rod as

(30) |

and the expression for the force is

(31) | ||||

(32) |

The normal and tangential force densities on the rod are thus

(33) | ||||

(34) |

With these forces, the evolution equation for the curvature (25) and the inextensibility constraint which determines the tension (24) are

(35) | |||

(36) |

These equations together with the boundary conditions at the ends of the rod

(37) |

determine the motion.

### a.4. Asymptotic Expansions

We examine the leading order behavior in two different limits: (1) small amplitude motion in which , and (2) high stiffness, , in which . Note that for soft bodies with large the realized amplitude, , is in fact small. Thus we expect the largest discrepancy between the asymptotic solutions at large at intermediate stiffness. In fact this is what we see in figure 4 for the flexor and figure 7 for the two-mode swimmer.

#### a.4.1. Small Amplitude

In the limit of small curvature, we see from equation (36) that the size of the tension is like the square of the curvature. At leading order the tension is zero, and the equation for the curvature is

(38) |

Changing to curvature deviation, , and nondimensionalizing gives equation (9) analyzed in the main text.

#### a.4.2. Large Stiffness

In the limit , . We introduce the variable

(39) |

to denote the deviation from the prescribed curvature. For increasing stiffness, . Changing variables from to and linearizing about small gives

(40) | |||

(41) |

and the boundary conditions are

(42) |

These equations contain many terms that are absent for small curvatures. However, as demonstrated in the text by comparing with numerical results, the low curvature equations appear to give a reasonable approximation at the amplitudes tested. Below we show why the two approximations are similar.

For simplicity we consider the flexor in which is only a function of time. With this simplification, equations (40)-(41) become

(43) | |||

(44) |

These equations can be solved by orthogonal function expansion. First we eliminate the tension by solving the second equation. We express as the series

(45) |

We then write the function as

(46) |

where represents the sequence of coefficients and is the orthogonal operator which maps these coefficients to , i.e. is like the Fourier transform operator. Using this expansion, the operator applied to in equation (44) diagonalizes, and the solution is

(47) |

where is a diagonal matrix with elements on the diagonal. After using this expression to eliminate in equation (43), after some simplification, we get the equation for the curvature deviation as

(48) |

where is a diagonal matrix with elements on the diagonal

(49) |

Because , the values of can be bounded as .

Equation (48) contains one additional term involving the second derivative that is not present in the corresponding low amplitude equation (9). Below we argue that the additional term is small even when itself is not small. Our analysis in the main text relied on performing an eigenfunction expansion using the eigenfunctions of the beam equation. We use the same expansion here

(50) |

and we relate to its expansion coefficients by

(51) |

where is the orthogonal operator that maps the expansion coefficients to . We can transform equation (48) into a system of differential equation for the expansion coefficients as

(52) |

where we define so that , and is a diagonal matrix containing the eigenvalues of the beam equation. That is, the diagonal entry of , , is related to the eigenvalue, , by . One expects that the contribution of to the equation to scale like . As argued above the norm of is about 1, and so we expect the additional terms relative to the bending terms to contribute . The smallest eigenvalue is about , and thus we expect these additional terms to be small.

In figure 9 we compare the expansion coefficients of the first mode of the small amplitude expansion and the high stiffness expansion for the high-amplitude flexor () for a range of stiffnesses. In agreement with our numerical results, the qualitative behavior of the two expansions is the same, and at high and low stiffnesses the quantitative behavior is the same. The largest difference occurs for moderately soft bodies where the difference is less than 25%.

## Appendix B Rod driven at one end

In this appendix we give the derivation for the shape-induced swimming speed of a thin filament in a viscous or viscoelastic fluid which is driven by oscillations at one end. This type of motion is akin to the experiments of [14], and we perform the calculation to compute the swimming speed as a function of the dimensionless body response time G and the fluid relaxation time We also show the range of G where a viscoelastic speed-up is theoretically predicted for this type of motion.

### b.1. Shape of the swimmer

The problem we consider is a flexible filament with one free end and one clamped end. The motion is driven by prescribing sinusoidal oscillations in the angle at the clamped end. The shape of the filament satisfies

(53) | |||

(54) | |||

(55) | |||

(56) | |||

(57) |

We now change variables using

(58) |

so that represents deviations from the infinitely stiff case of a straight rod. Letting be complex valued, the equation for is then

(59) |

with homogeneous boundary conditions. This equation can be solved using an expansion of eigenfunctions, , which satisfy

We express the function using an eigenfunction expansion

(60) |

and look for a solution to the PDE of the form

(61) |

The transformation of the PDE yields

(62) |

which gives

(63) |

We can then write the shape as

(64) | ||||

(65) |

Notice that the factor multiplying above is exactly the same as the one that appears in the shape analysis of swimmers driven by active moments in (12). We can express the shape as

(66) |

where is defined by (12). As in the main body of the paper, to add viscoelastic effects, we simply use (14) in place of (12) to define the expansion coefficients. Although these expressions are the same, we note that the eigenfunctions and eigenvalues are different for this problem.

### b.2. Expression for Swimming Speed

The expansion for the shape of the swimmer can be written as

(67) |

where

(68) | |||

(69) |

To compute the swimming speed, we average in space and time the product

(70) |

Averaging the above expression by integration gives the swimming speed

(71) |

We use this expression with the first six modes to compute the shape-induced Stokes-normalized swimming speed; see figure 10. As with the problem from the paper, for sufficiently soft bodies, we see an almost monotonic speed-up from the shape changes. For sufficiently stiff swimmers () we see a monotonic slow down. There is a transition range around . The location of this transition is evident in figure 10 (b) where we show the swimming speed as a function of G for low and high . The qualitative results from the paper do not change in the sufficiently soft regime, but this problem has a different driving mechanism and hence there is a different effect in the stiff regime. As we see in figure 10 (c) as the body is stiffened the swimming speed goes to zero and viscoelasticity always slows the swimmer.

## References

- [1] T K Chaudhury. On swimming in a visco-elastic liquid. Journal of Fluid Mechanics, 95(01):189–197, 1979.
- [2] LD Sturges. Motion induced by a waving plate. Journal of Non-Newtonian Fluid Mechanics, 8(3-4):357–364, 1981.
- [3] Eric Lauga. Propulsion in a viscoelastic fluid. Phys. Fluids, 19(8):83104, 2007.
- [4] Henry C Fu, Thomas R Powers, and Charles W Wolgemuth. Theory of swimming filaments in viscoelastic media. Physical review letters, 99(25):258101, 2007.
- [5] Henry C Fu, Charles W Wolgemuth, and Thomas R Powers. Swimming speeds of filaments in nonlinearly viscoelastic fluids. Physics of Fluids (1994-present), 21(3):33102, 2009.
- [6] Gwynn J Elfring and Gaurav Goyal. The effect of gait on swimming in viscoelastic fluids. Journal of Non-Newtonian Fluid Mechanics, 234:8–14, 2016.
- [7] Emily E Riley and Eric Lauga. Small-amplitude swimmers can self-propel faster in viscoelastic fluids. Journal of theoretical biology, 382:345–355, 2015.
- [8] Mark P Curtis and Eamonn A Gaffney. Three-sphere swimmer in a nonlinear viscoelastic medium. Physical Review E, 87(4):043006, 2013.
- [9] X N Shen and Paulo E Arratia. Undulatory swimming in viscoelastic fluids. Physical review letters, 106(20):208101, 2011.
- [10] DA Gagnon, NC Keim, and PE Arratia. Undulatory swimming in shear-thinning fluids: experiments with caenorhabditis elegans. Journal of Fluid Mechanics, 758:R3, 2014.
- [11] B Qin, A Gopinath, J Yang, Jerry P Gollub, and PE Arratia. Flagellar kinematics and swimming of algal cells in viscoelastic fluids. Scientific reports, 5, 2015.
- [12] Bin Liu, Thomas R Powers, and Kenneth S Breuer. Force-free swimming of a model helical flagellum in viscoelastic fluids. Proceedings of the National Academy of Sciences, 108(49):19516–19520, 2011.
- [13] Nathan C. Keim, Mike Garcia, and Paulo E. Arratia. Fluid elasticity can enable propulsion at low Reynolds number. Physics of Fluids (1994-present), 24(8):81703, 2012.
- [14] Julian Espinosa-Garcia, Eric Lauga, Roberto Zenit, San Diego, and La Jolla. Fluid elasticity increases the locomotion of flexible swimmers. Physics of Fluids (1994-present), 25(3):31701, 2013.
- [15] Moumita Dasgupta, Bin Liu, Henry C Fu, Michael Berhanu, Kenneth S Breuer, Thomas R Powers, and Arshad Kudrolli. Speed of a swimming sheet in Newtonian and viscoelastic fluids. Physical Review E, 87(1):13015, 2013.
- [16] Francisco A Godínez, Lyndon Koens, Thomas D Montenegro-Johnson, Roberto Zenit, and Eric Lauga. Complex fluids affect low-reynolds number locomotion in a kinematic-dependent manner. Experiments in Fluids, 56(5):1–10, 2015.
- [17] Joseph Teran, Lisa Fauci, and Michael Shelley. Viscoelastic fluid response can increase the speed and efficiency of a free swimmer. Physical review letters, 104(3):38101, January 2010.
- [18] NJ Balmforth, D Coombs, and S Pachmann. Microelastohydrodynamics of swimming organisms near solid boundaries in complex fluids. QJ Mech. Appl. Math, 63(3):267–294, 2010.
- [19] Lailai Zhu, Eric Lauga, and Luca Brandt. Self-propulsion in viscoelastic fluids: Pushers vs. pullers. Physics of Fluids (1994-present), 24(5):51902, 2012.
- [20] Saverio E Spagnolie, Bin Liu, and Thomas R Powers. Locomotion of helical bodies in viscoelastic fluids: enhanced swimming at large helical amplitudes. Physical review letters, 111(6):68101, 2013.
- [21] Thomas D. Montenegro-Johnson, David J. Smith, and Daniel Loghin. Physics of rheologically enhanced propulsion: Different strokes in generalized Stokes. Physics of Fluids, 25(8):081903, 2013.
- [22] Becca Thomases and Robert D. Guy. Mechanisms of Elastic Enhancement and Hindrance for Finite-Length Undulatory Swimmers in Viscoelastic Fluids. Physical Review Letters, 113(9):098102, August 2014.
- [23] G-J Li, A Karimi, and AM Ardekani. Effect of solid boundaries on swimming dynamics of microorganisms in a viscoelastic fluid. Rheologica acta, 53(12):911–926, 2014.
- [24] Gaojin Li and Arezoo M Ardekani. Undulatory swimming in non-newtonian fluids. Journal of Fluid Mechanics, 784:R4, 2015.
- [25] Daniel Salazar, Alexandre M Roma, and Hector D Ceniceros. Numerical study of an inextensible, finite swimmer in stokesian viscoelastic flow. Physics of Fluids (1994-present), 28(6):063101, 2016.
- [26] Gwynn J Elfring and Eric Lauga. Theory of locomotion through complex fluids. In Complex fluids in biological systems, pages 283–317. Springer, 2015.
- [27] JosuÚ Sznitman and Paulo E Arratia. Locomotion through complex fluids: an experimental view. In Complex Fluids in Biological Systems, pages 245–281. Springer, 2015.
- [28] Emily E Riley and Eric Lauga. Enhanced active swimming in viscoelastic fluids. EPL (Europhysics Letters), 108(3):34003, 2014.
- [29] Robert D Guy and Becca Thomases. Computational challenges for simulating strongly elastic flows in biology. In Complex Fluids in Biological Systems, pages 359–397. Springer, 2015.
- [30] Lisa J Fauci and Charles S Peskin. A computational model of aquatic animal locomotion. Journal of Computational Physics, 77(1):85–108, 1988.
- [31] R. B. Bird, O. Hassager, R. C. Armstrong, and C. F. Curtiss. Dynamics of Polymeric Liquids, Vol. 2: Kinetic Theory. John Wiley and Sons, 1980.
- [32] R Sureshkumar and Antony N Beris. Effect of artificial stress diffusivity on the stability of numerical calculations and the flow dynamics of time-dependent viscoelastic flows. J. Non-Newton. Fluid Mech., 60(1):53–80, 1995.
- [33] Becca Thomases. An analysis of the effect of stress diffusion on the dynamics of creeping viscoelastic flow. Journal of Non-Newtonian Fluid Mechanics, 166(21-22):1221–1228, 2011.
- [34] Chris H Wiggins and Raymond E Goldstein. Flexive and propulsive dynamics of elastica at low reynolds number. Physical Review Letters, 80(17):3879, 1998.
- [35] Eric Lauga. Floppy swimming: Viscous locomotion of actuated elastica. Physical Review E, 75(4):041916, 2007.
- [36] Henry C Fu, Charles W Wolgemuth, and Thomas R Powers. Beating patterns of filaments in viscoelastic fluids. Physical Review E, 78(4):041913, 2008.
- [37] Michael J Shelley and Tetsuji Ueda. The stokesian hydrodynamics of flexing, stretching filaments. Physica D: Nonlinear Phenomena, 146(1):221–245, 2000.
- [38] Raymond E. Goldstein and Stephen A. Langer. Nonlinear dynamics of stiff polymers. Phys. Rev. Lett., 75:1094–1097, Aug 1995.
- [39] Sébastien Camalet and Frank Jülicher. Generic aspects of axonemal beating. New Journal of Physics, 2(1):24, 2000.
- [40] Glenn R Fulford, David F Katz, and Robert L Powell. Swimming of spermatozoa in a linear viscoelastic fluid. Biorheology, 35(4):295–309, 1998.
- [41] Ricardo Cortez. The method of regularized stokeslets. SIAM Journal on Scientific Computing, 23(4):1204–1225, 2001.