Home | About | Journals | Submit | Contact Us | Français |

**|**Proc Math Phys Eng Sci**|**PMC5014114

Formats

Article sections

- Abstract
- 1. Introduction
- 2. Problem formulation
- 3. Perturbation expansion of the dynamics
- 4. Approximate expressions for Lighthill's efficiency
- 5. Conclusion
- Supplementary Material
- References

Authors

Related links

Proc Math Phys Eng Sci. 2016 August; 472(2192): 20160425.

PMCID: PMC5014114

Faculty of Mechanical Engineering, Technion – Israel Institute of Technology, Haifa 32000, Israel

e-mail: li.ca.noinhcet@izi

Received 2016 June 2; Accepted 2016 July 15.

Copyright © 2016 The Author(s)

Published by the Royal Society. All rights reserved.

This work studies the motion of Purcell's three-link microswimmer in viscous flow, by using perturbation expansion of its dynamics under small-amplitude strokes. Explicit leading-order expressions and next-order correction terms for the displacement of the swimmer are obtained for the cases of a square or circular gait in the plane of joint angles. The correction terms demonstrate the reversal in movement direction for large stroke amplitudes, which has previously only been shown numerically. In addition, asymptotic expressions for Lighthill's energetic efficiency are obtained for both gaits. These approximations enable calculating optimal stroke amplitudes and swimmer's geometry (i.e. ratio of links’ lengths) for maximizing either net displacement or Lighthill's efficiency.

The study of micrometre-size swimmers dynamics has in recent years become a highly active research area and has possible implications on understanding the motion of swimming microorganisms and biological infections [1,2]. The considerable advances made in the field of micro- and nanotechnology have promoted the possibility of manufacturing miniature robotic devices that operate in these small scales [3–5]. Such mini-robots may have many applications in medicine, performing medical procedures in a minimally invasive way and delivering drugs with high precision [6–8]. All this requires an understanding of swimming dynamics at the low Reynolds number regime.

The Reynolds number represents the ratio of the inertial forces to the viscous ones. When dealing in microfluidics, where the characteristic lengths are extremely small, hydrodynamic forces are typically governed by very low Reynolds numbers (*Re*1) [9]. The result of this is that the strategy of motion in this regime needs to be drastically different from the familiar motion of larger organisms, such as fish, that rely on imparting momentum to the surrounding fluid. In his famous lecture [10], Purcell introduced the ‘Scallop theorem’, which states that a swimmer that changes its shape and then changes back to the original shape by reversing the same sequence will return to the point it started its motion with no regard to the speed in which any part of the motion is made. Any reciprocal motion of the swimmer will result in zero net translation. There are several ways of overcoming the scallop theorem and generating net motion in a low Reynolds regime. One of them is continuously performing a unidirectional rigid body motion such as a rotating corkscrew. In this way, there is no reciprocal motion and the result is generation of net motion of the swimmer. This is precisely the method used by the *Escherichia coli* bacteria to propel itself [11]. Another way for overcoming the scallop theorem is by making a non-reciprocal periodic shape change, which will be henceforth called a ‘gait’. The simplest version of such a swimmer, known as ‘Purcell's three-link swimmer’ (figure 1), was suggested by Purcell in [10] and can be seen as a simplified version of the travelling wave ‘Taylor sheet’ [12] in two dimensions which is discretized to have only 2 degrees of freedom. Purcell's swimmer comprises three rigid links connected by two rotary joints (figure 1). Purcell indicated that this swimmer could propel itself along a straight line by alternately rotating its front and back links in a non-reciprocal way (square gait, figure 2*a*). Through symmetry considerations alone, it can be shown that the three-link swimmer will move along the *x*-axis when using the shape changes suggested by Purcell. Purcell claimed that determining the direction of net motion (i.e. forward or backward) is trivial and left it as ‘an exercise for the student’. Only 26 years later, Becker *et al.* [13] obtained an explicit formulation for the dynamics of this microswimmer, and surprisingly found that the direction of net motion actually depends on the stroke amplitude of joint angles. Specifically, for small amplitudes the swimmer will move in one direction but for larger amplitudes the swimmer will move in the opposite direction. Additionally, Becker *et al.* [13] studied Lighthill's energetic efficiency, which is roughly equivalent to the maximal mean speed achievable under a given mean expenditure of mechanical power. It is worth noting that both [10,13] considered only the case where the joint angles rotate alternately, creating a square gait and did not study other possible periodic shape change, such as the circular gait (figure 2*b*), at which the two joint angles oscillate sinusoidally with a quarter-period phase shift.

Gait representation in joint angle plane. (*a*) Square gait and (*b*) circular gait. (Online version in colour.)

Some previous works have examined Purcell's swimmer from different perspectives. In [14], Gutman and Or analyse the symmetries of Purcell's swimmer, derive conditions on gaits that result in movement along the swimmer's principal directions and present motion experiments of a macro-scale robotic swimmer in a highly viscous fluid. Avron & Raz [15] and Hatton and colleagues [16,17] use tools of differential geometry in order to obtain geometric visualization of performance measures in shape space (plane). By observing curvature maps of the swimmer's dynamics, those works can qualitatively explain the reversal in direction of motion for large-amplitude strokes, and also provide guidelines that help in obtaining optimal gaits. Tam & Hosoi [18], through numerical computation only, found gaits that achieve optimal translation during a period as well as energetically optimal gaits. They also found the *optimal geometry* of the swimmer, i.e. ratio between the swimmer's links lengths, for these two optimality criteria, again only through numeric calculations.

Aside from numerical computations, another useful approach is obtaining closed-form expressions under some simplifying scaling assumptions, by using asymptotic analysis. In microswimmers, this approach dates back to the analysis of Taylor [12], who showed that the swimming speed of Taylor's sheet scales quadratically with the wave's amplitude at leading order. In [19], this concept was employed on analysis of self-propulsion of spherical swimmers performing a squirming motion. For the shape-changing three-link swimmer, one has to use the method of perturbation expansion [20] under the assumption of small stroke amplitudes. This method was used in [13] for obtaining leading-order expressions for the motion of Purcell's swimmer and finding its optimal geometry. In [21], the leading-order expressions were derived using a similar approach and then used to find optimal geometry of the swimmer for the square gait and for derivation of an optimal polygonal gait under small joint angles. Interestingly, there are differences in the results found in [13,18,21] in terms of the swimmer's optimal geometry. Recently, perturbation expansion has become a more common tool for analysis of simple microswimmers motion in various cases such as a three-link swimmer with a passive elastic joint [22], a two-link magnetically actuated microswimmer [23,24] and wobbling of a magnetically actuated helix [25].

The goal of this work is to introduce a systematic method for analysing the dynamics of Purcell's three-link swimmer using perturbation expansion and to exploit this method to perform optimizations on the swimmer geometry and stroke amplitude. The main new contributions of this work compared with previous literature are (i) formulation of leading-order expression for the displacement under the circular gait (figure 2*b*); (ii) derivation of the next-order correction terms for the displacement under both square and circular gaits, which explicitly show the reversal in movement direction first found in [13], as well as obtaining an approximation for the optimal stroke amplitude; (iii) obtaining leading-order expressions for Lighthill's energetic efficiency of the two presented gaits and also a next-order term for the square gait, as well as using these expressions to perform optimization of the swimmer's geometry and stroke amplitude for maximal efficiency. (iv) Resolving the previously unexplained differences between the findings in [13,18,21] regarding the optimal geometry.

The paper is organized as follows. The next section presents Purcell's three-link swimmer model and its equations of motion, the two gaits in question and the concept of Lighthill's energetic efficiency. Section 3 focuses on perturbation expansion and presents a systematic way for deriving the net displacement of the swimmer in the form of a power series for any gait, and the first two terms are found explicitly for the square and circular ones. The optimal geometry and stroke amplitude for maximal displacement are also found. This section also includes the use of perturbation expansion to find the period time under constant joint torque for the square gait, and under constant power for both square and circular gaits. In §4, the previous results are used in order to obtain the leading-order terms of the energetic efficiency of the two gaits, as well as next-order term for the square gait. Optimal geometry and amplitude are then found for maximal efficiency. Finally, §5 offers our conclusions, indicating possible future extensions and consequences.

In this section, we introduce Purcell's swimmer model, formulate its dynamic equations of motion, define the square and circular gaits and, finally, review the concept of Lighthill's energetic efficiency. The swimmer in question consists of three thin rigid links with lengths *l*_{0},*l*_{1},*l*_{2}, with *l*_{1}=*l*_{2}, *l*=*l*_{0}+*l*_{1}+*l*_{2} and *η*=*l*_{0}/*l*. The links are connected by two rotary joints whose angles are denoted by *ϕ*_{1} and *ϕ*_{2} (figure 1). The shape of the swimmer will be described by these two angles ** ϕ**=(

The kinematic relation between body velocity, joints velocities and links velocities is given by

$$\begin{array}{rl}{\mathbf{\text{v}}}_{i}& ={\mathbf{\text{v}}}_{b}+{\omega}_{b}\mathbf{\text{z}}\times ({\mathbf{\text{r}}}_{i}-{\mathbf{\text{r}}}_{0})+\sum _{j}{I}_{ij}{\dot{\varphi}}_{j}\mathbf{\text{z}}\times ({\mathbf{\text{r}}}_{i}-{\mathbf{\text{b}}}_{j})\\ \text{and}\phantom{\rule{2em}{0ex}}{\omega}_{i}& ={\omega}_{b}+\sum _{j}{I}_{ij}{\dot{\varphi}}_{j},\end{array}\}$$

2.1

where **r**_{i} is the position of the centre of the *i*th link and **b**_{j} is the position of the *j*th joint. In matrix form, the velocity **V**_{i} is related to the body velocity **V**_{b} and shape velocity $\dot{\mathit{\varphi}}$ through

$${\mathbf{\text{V}}}_{i}={\mathbf{\text{T}}}_{i}(\mathbf{\text{q}},\mathit{\varphi})\dot{\mathbf{\text{q}}}+{\mathbf{\text{E}}}_{i}(\mathbf{\text{q}},\mathit{\varphi})\dot{\mathit{\varphi}}.$$

2.2

The matrices **T**_{i}(**q**,** ϕ**) and

$$\begin{array}{rl}{\mathbf{\text{T}}}_{0}& =\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right],\phantom{\rule{1em}{0ex}}{\mathbf{\text{E}}}_{0}=\left[\begin{array}{cc}0& 0\\ 0& 0\\ 0& 0\end{array}\right]\\ {\mathbf{\text{T}}}_{1}& =\left[\begin{array}{ccc}1& 0& -0.5{l}_{0}\mathrm{sin}{\alpha}_{0}-0.5{l}_{1}\mathrm{sin}{\alpha}_{1}\\ 0& 1& 0.5{l}_{0}\mathrm{cos}{\alpha}_{0}+0.5{l}_{1}\mathrm{cos}{\alpha}_{1}\\ 0& 0& 1\end{array}\right],\phantom{\rule{1em}{0ex}}{\mathbf{\text{E}}}_{1}=\left[\begin{array}{cc}-0.5{l}_{1}\mathrm{sin}{\alpha}_{1}& 0\\ 0.5{l}_{1}\mathrm{cos}{\alpha}_{1}& 0\\ 1& 0\end{array}\right]\\ \text{and}\phantom{\rule{2em}{0ex}}{\mathbf{\text{T}}}_{2}& =\left[\begin{array}{ccc}1& 0& 0.5{l}_{0}\mathrm{sin}{\alpha}_{0}+0.5{l}_{2}\mathrm{sin}{\alpha}_{2}\\ 0& 1& -0.5{l}_{0}\mathrm{cos}{\alpha}_{0}-0.5{l}_{2}\mathrm{cos}{\alpha}_{2}\\ 0& 0& 1\end{array}\right],\phantom{\rule{1em}{0ex}}{\mathbf{\text{E}}}_{2}=\left[\begin{array}{cc}0& -0.5{l}_{2}\mathrm{sin}{\alpha}_{2}\\ 0& 0.5{l}_{2}\mathrm{cos}{\alpha}_{2}\\ 0& -1\end{array}\right],\end{array}\}$$

2.3

where *α*_{0}=*θ*, *α*_{1}=*ϕ*_{1}−*θ* and *α*_{2}=*θ*+*ϕ*_{2} are the absolute orientation angles of each link. Next, we invoke *resistive force theory* [26,27], which states that the viscous drag force **f**_{i} and torque *m*_{i} on the *i*th slender link under planar motion are proportional to its linear and angular velocities. Thus, we can write the expression for the drag force and torque exerted on each link:

$$\begin{array}{rl}{\mathbf{\text{f}}}_{i}& =-{c}_{\mathrm{t}}{l}_{i}({\mathbf{\text{v}}}_{i}\cdot {\mathbf{\text{t}}}_{i}){\mathbf{\text{t}}}_{i}-{c}_{\mathrm{n}}{l}_{i}({\mathbf{\text{v}}}_{i}\cdot {\mathbf{\text{n}}}_{i}){\mathbf{\text{n}}}_{i}\\ \text{and}\phantom{\rule{2em}{0ex}}{m}_{i}& =-\frac{1}{12}{c}_{\mathrm{n}}{l}_{i}^{3}{\omega}_{i},\end{array}\}$$

2.4

where ${\mathbf{\text{t}}}_{i}={(\mathrm{cos}{\alpha}_{i},\mathrm{sin}{\alpha}_{i})}^{\mathrm{T}}$ is a unit vector in the axial direction of the *i*th link and ${\mathbf{\text{n}}}_{i}={(-\mathrm{sin}{\alpha}_{i},\mathrm{cos}{\alpha}_{i})}^{\mathrm{T}}$ is a unit vector in the normal direction. The resistance coefficients for the normal and axial directions are ${c}_{\mathrm{n}}=2{\mathrm{c}}_{\mathrm{t}}=4\pi \mu /\mathrm{log}({l}_{\mathrm{c}}/a),$ where *μ* is the dynamic viscosity of the fluid, *a* is the radius of the links and *l*_{c} is a characteristic length. As the ratio of the link's length to its radius is assumed to be very large, the difference in the resistance coefficients between the links is very small and will be neglected. It is also assumed that effects of hydrodynamic interaction between the links are negligible. Denoting the vector of forces and torques on the *i*th link as **F**_{i}=(**f**_{i},*m*_{i}), the relation (2.4) can be written in matrix form

$${\mathbf{\text{F}}}_{i}=-{\mathcal{R}}_{i}(\mathbf{\text{q}},\mathit{\varphi}){\mathbf{\text{V}}}_{i},$$

2.5

where

$${\mathcal{R}}_{i}(\mathbf{\text{q}},\mathit{\varphi})={\mathrm{c}}_{\mathrm{t}}^{(i)}{l}_{i}\left[\begin{array}{ccc}1+{\mathrm{sin}}^{2}{\alpha}_{i}& -\mathrm{cos}{\alpha}_{i}\mathrm{sin}{\alpha}_{i}& 0\\ -\mathrm{cos}{\alpha}_{i}\mathrm{sin}{\alpha}_{i}& 1+{\mathrm{cos}}^{2}{\alpha}_{i}& 0\\ 0& 0& \frac{1}{6}{l}_{i}^{2}\end{array}\right]$$

2.6

is called the resistance tensor. The total drag forces and torques acting on the swimmer's body are given by

$${\mathbf{\text{f}}}_{b}=\sum _{i=0}^{2}{\mathbf{\text{f}}}_{i},\phantom{\rule{1em}{0ex}}{m}_{b}=\sum _{i=0}^{2}({m}_{i}+(({\mathbf{\text{r}}}_{b}-{\mathbf{\text{r}}}_{i})\times {\mathbf{\text{f}}}_{i})\cdot \mathbf{\text{z}}).$$

2.7

Using the matrices **T**_{i} from the kinematic relation (2.2) and augmenting in a vector **F**_{b}=(**f**_{b},*m*_{b}), (2.7) is written in matrix form as

$${\mathbf{\text{F}}}_{b}=\sum _{i=0}^{2}{\mathbf{\text{T}}}_{i}^{\mathrm{T}}{\mathbf{\text{F}}}_{i}=-\sum _{i=0}^{2}{\mathbf{\text{T}}}_{i}^{\mathrm{T}}{\mathcal{R}}_{i}({\mathbf{\text{T}}}_{i}{\mathbf{\text{V}}}_{b}+{\mathbf{\text{E}}}_{i}\dot{\mathit{\varphi}}).$$

2.8

Denoting

$${\mathcal{R}}_{bb}=\sum _{i=0}^{2}{\mathbf{\text{T}}}_{i}^{\mathrm{T}}{\mathcal{R}}_{i}{\mathbf{\text{T}}}_{i}\phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{1em}{0ex}}{\mathcal{R}}_{bu}=\sum _{i=0}^{2}{\mathbf{\text{T}}}_{i}^{\mathrm{T}}{\mathcal{R}}_{i}{\mathbf{\text{E}}}_{i},$$

2.9

equation (2.8) becomes

$${\mathbf{\text{F}}}_{b}={\mathcal{R}}_{bb}{\mathbf{\text{V}}}_{b}+{\mathcal{R}}_{bu}\dot{\mathit{\varphi}}.$$

2.10

Our choice of coordinates of body location **q** implies that the body velocity satisfies $\dot{\mathbf{\text{q}}}={\mathbf{\text{V}}}_{b}$. Assuming quasi-static motion, the swimmer is in static equilibrium **F**_{b}=0. Substituting this into (2.10) then yields the nonlinear differential equations that govern the swimmer's dynamics:

$$\dot{\mathbf{\text{q}}}=\mathcal{G}(\mathbf{\text{q}},\mathit{\varphi})\dot{\mathit{\varphi}},$$

2.11

where

$$\mathcal{G}(\mathbf{\text{q}},\mathit{\varphi})=-{\mathcal{R}}_{bb}^{-1}{\mathcal{R}}_{bu}.$$

2.12

As there are no external boundary conditions in the unbounded fluid domain, the swimmer's velocity expressed in body-fixed reference frame is independent of the position **q** and thus (2.11) can be written as (cf. [14])

$$\begin{array}{rl}\dot{\mathbf{\text{q}}}& =\mathbf{\text{D}}(\theta )\mathbf{\text{G}}(\mathit{\varphi})\dot{\mathit{\varphi}},\\ \text{where}\phantom{\rule{2em}{0ex}}\mathbf{\text{D}}(\theta )& =\left[\begin{array}{ccc}\mathrm{cos}\theta & -\mathrm{sin}\theta & 0\\ \mathrm{sin}\theta & \mathrm{cos}\theta & 0\\ 0& 0& 1\end{array}\right].\end{array}\}$$

2.13

Note that (2.13) implies that the angular velocity of the swimmer $\dot{\theta}$ is independent of *θ*. The matrix **G**(** ϕ**) obeys some symmetry relations due to the swimmer's structure (see [14] for details):

$$\begin{array}{rl}\mathbf{\text{G}}(-\mathit{\varphi})& =\left[\begin{array}{ccc}-1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]\mathbf{\text{G}}(\mathit{\varphi}))\\ \text{and}\phantom{\rule{2em}{0ex}}\mathbf{\text{G}}(\mathbf{\text{S}}\mathit{\varphi})& =\left[\begin{array}{ccc}-1& 0& 0\\ 0& 1& 0\\ 0& 0& -1\end{array}\right]\mathbf{\text{G}}(\mathit{\varphi})\mathbf{\text{S}},\end{array}\}$$

2.14

where $\mathbf{\text{S}}=\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]$ represents interchanging between the joint angles *ϕ*_{1} and *ϕ*_{2}.

Next, we compute the actuation torques acting at the joints, which are denoted by the vector ** τ**=(

$${\tau}_{j}=-\sum _{i}{I}_{ij}({m}_{i}+(({\mathbf{\text{r}}}_{i}-{\mathbf{\text{b}}}_{j})\times {\mathbf{\text{f}}}_{i})\cdot \mathbf{\text{z}}).$$

2.15

This relation can be written in matrix form using the previously defined matrices **E**_{i}:

$$\mathit{\tau}=-\sum _{i=0}^{2}{\mathbf{\text{E}}}_{i}^{\mathrm{T}}{\mathbf{\text{F}}}_{i}.$$

2.16

Substituting (2.2), (2.5) and (2.11) in (2.16) gives

$$\begin{array}{rl}\mathit{\tau}& =\sum _{i=0}^{2}{\mathbf{\text{E}}}_{i}^{\mathrm{T}}{\mathcal{R}}_{i}{\mathbf{\text{V}}}_{i}\\ & =\sum _{i=0}^{2}{\mathbf{\text{E}}}_{i}^{\mathrm{T}}{\mathcal{R}}_{i}({\mathbf{\text{T}}}_{i}{\dot{\mathbf{\text{q}}}}_{i}+{\mathbf{\text{E}}}_{i}{\dot{\mathit{\varphi}}}_{i})\\ & =(\sum _{i=0}^{2}{\mathbf{\text{E}}}_{i}^{\mathrm{T}}{\mathcal{R}}_{i}{\mathbf{\text{E}}}_{i}-\sum _{i=0}^{2}{\mathbf{\text{E}}}_{i}^{\mathrm{T}}{\mathcal{R}}_{i}{\mathbf{\text{T}}}_{i}{\left(\sum _{i=0}^{2}{\mathbf{\text{T}}}_{i}^{\mathrm{T}}{\mathcal{R}}_{i}{\mathbf{\text{T}}}_{i}\right)}^{-1}\sum _{i=0}^{2}{\mathbf{\text{T}}}_{i}^{\mathrm{T}}{\mathcal{R}}_{i}{\mathbf{\text{E}}}_{i})\dot{\mathit{\varphi}}\end{array}$$

2.17

denoting ${\mathcal{R}}_{uu}=\sum _{i=0}^{2}{\mathbf{\text{E}}}_{i}^{\mathrm{T}}{\mathcal{R}}_{i}{\mathbf{\text{E}}}_{i}$; this equation can be written as

$$\mathit{\tau}=({\mathcal{R}}_{uu}-{\mathcal{R}}_{bu}^{\mathrm{T}}{\mathcal{R}}_{bb}^{-1}{\mathcal{R}}_{bu})\dot{\mathit{\varphi}}$$

2.18

with ${\mathcal{R}}_{bb}$ and ${\mathcal{R}}_{bu}$ as defined earlier in (2.9). Denoting

$$\mathbf{\text{W}}(\mathit{\varphi})={\mathcal{R}}_{uu}-{\mathcal{R}}_{bu}^{\mathrm{T}}{\mathcal{R}}_{bb}^{-1}{\mathcal{R}}_{bu},$$

2.19

we have

$$\mathit{\tau}=\mathbf{\text{W}}(\mathit{\varphi})\dot{\mathit{\varphi}}.$$

2.20

The resulting dynamic equations in (2.11) and (2.20) are strongly nonlinear. Nevertheless, they can be expanded by assuming small-amplitude changes of the angles about ** ϕ**=0, as explained in §3.

This work considers time-periodic inputs of shape changes ** ϕ**(

$$\mathrm{square}:\hspace{0.17em}{s}_{1}(t)=\{\begin{array}{ll}1,& t\in [0,2]\\ (3-t),& t\in [2,4]\\ -1,& t\in [4,6]\\ (t-7),& t\in [6,8]\end{array},\phantom{\rule{1em}{0ex}}{s}_{2}(t)=\{\begin{array}{ll}(1-t),& t\in [0,2]\\ -1,& t\in [2,4]\\ (t-5),& t\in [4,6]\\ 1,& t\in [6,8]\end{array}$$

2.21

and

$$\mathrm{circular}:\hspace{0.17em}{s}_{1}(t)=\mathrm{sin}(t+\frac{\pi}{4}),\phantom{\rule{1em}{0ex}}{s}_{2}(t)=\mathrm{cos}(t+\frac{\pi}{4}),\hspace{0.17em}t\in [0,2\pi ].$$

2.22

These equations describe circular and square-shaped trajectories with stroke amplitude of 1, which are then scaled to stroke of *ε* by setting *ϕ*_{i}(*t*)=*εs*_{i}(*t*). As the equation of motion (2.11) is time invariant, the net motion is independent of time parametrization of the gait. Here, time parametrization of the gaits was chosen arbitrarily so that the period times are *T*=8 for the square gait and *T*=2*π* for the circular one. The net displacement in the *x*-direction of one full period will be denoted *X* and is calculated through $X={\int}_{0}^{T}\dot{x}\hspace{0.17em}\mathrm{d}t$. The mean swimming speed is denoted by *V* =*X*/*T*.

The energetic efficiency of a stroke can be defined in a several different ways, cf. [28]. First, we write an expression for the power exerted by the swimmer. The mechanical power dissipated by the fluid's viscous drag forces and torques on all three links is

$$P=-\sum _{i=0}^{2}{\mathbf{\text{F}}}_{i}^{\mathrm{T}}{\mathbf{\text{V}}}_{i}=\sum _{i=0}^{2}{\mathbf{\text{V}}}_{i}^{\mathrm{T}}{\mathcal{R}}_{i}{\mathbf{\text{V}}}_{i}.$$

2.23

On the other hand, the mechanical power expended internally by the actuation torques is

$$P={\mathit{\tau}}^{\mathrm{T}}\dot{\mathit{\varphi}}={\dot{\mathit{\varphi}}}^{\mathrm{T}}\mathbf{\text{W}}\dot{\mathit{\varphi}}.$$

2.24

These last two expressions are equivalent, which can be proved using the relations (2.19) and (2.20). The total work can be calculated from here by $W={\int}_{0}^{T}P(t)\hspace{0.17em}\mathrm{d}t$.

Owing to time-invariance of the swimmer's dynamics, a well-known observation (cf. [13,18]) states that the energy expenditure under a given gait trajectory can be made arbitrarily small by pacing along the trajectory in a sufficiently slow rate. Thus, energy per unit distance cannot serve as a reasonable performance measure. Following Becker *et al.* [13] and Tam & Hosoi [18], we use an energetic efficiency criterion similar to that defined by Lighthill in [29] which is the ratio of average power $\overline{P}=W/T$ exerted by the swimmer to the power needed to drag the swimmer as a rigid body at the same mean speed:

$$\stackrel{~}{\xi}=\frac{{c}_{\mathrm{t}}l{V}^{2}}{\overline{P}}=\frac{{c}_{\mathrm{t}}l{X}^{2}}{\overline{P}{T}^{2}}.$$

2.25

This definition is non-unique for a given gait trajectory ** ϕ**(

The period time *T* under a constant mechanical power *P*=*P*_{o} can be found using the following derivation. Consider the gait ** ϕ**(

$$P(\sigma (t))=({\mathit{\varphi}}^{\prime}{(\sigma )}^{\mathrm{T}}\mathbf{\text{W}}(\sigma ){\mathit{\varphi}}^{\prime}(\sigma )){\dot{\sigma}}^{2},\phantom{\rule{1em}{0ex}}\text{where}\hspace{0.17em}{\mathit{\varphi}}^{\prime}(\sigma )=\frac{\mathrm{\partial}\mathit{\varphi}}{\mathrm{\partial}\sigma}$$

2.26

denoting *F*(*σ*)=** ϕ**′(

$$\dot{\sigma}=\frac{\mathrm{d}\sigma}{\mathrm{d}t}=\sqrt{\frac{P}{F(\sigma )}}\to \mathrm{d}t=\sqrt{\frac{F(\sigma )}{P}}\hspace{0.17em}\mathrm{d}\sigma .$$

2.27

This transformation is well defined for any non-degenerate trajectory and time parametrizations such that ** ϕ**′(

$${T}_{{p}_{o}}=\frac{1}{\sqrt{{P}_{o}}}{\int}_{{\sigma}_{0}}^{{\sigma}_{1}}\sqrt{F(\sigma )}\hspace{0.17em}\mathrm{d}\sigma .$$

2.28

Substituting (2.28) into the expression for Lighthill's efficiency (2.25) under a constant power $\overline{P}={P}_{o}$, it is clear that the *P*_{o} cancels out and so the calculations can be done for *P*_{o}=1 without loss of generality. Thus, Lighthill's efficiency of a gait ** ϕ**(

$$\xi =\frac{{X}^{2}}{{T}_{p}^{2}},$$

2.29

where *T*_{p} now denotes the period time under constant power of *P*_{o}=1, as given in (2.28), and the constants *c*_{t}*l* from (2.25) are dropped.

In order to find the leading-order expressions for the swimmer's motion using perturbation expansion [20], the dynamics of the swimmer are expanded as power series of the stroke amplitude *ε*. First, expansion of the swimmer's orientation angle *θ*(*t*) is obtained, followed by expansions of the instantaneous body velocity and of the resulting net displacement *X* under each specific gait.

The position and orientation of the swimmer are now expanded into a power series in *ε* as

$$\mathbf{\text{q}}(t)=\epsilon {\mathbf{\text{q}}}^{(1)}(t)+{\epsilon}^{2}{\mathbf{\text{q}}}^{(2)}(t)+{\epsilon}^{3}{\mathbf{\text{q}}}^{(3)}(t)+\cdots .$$

3.1

And for each coordinate

$$\begin{array}{rl}x(t)& =\epsilon {x}^{(1)}(t)+{\epsilon}^{2}{x}^{(2)}(t)+{\epsilon}^{3}\hspace{0.17em}{x}^{(3)}(t)+\cdots ,\end{array}$$

3.2

$$\begin{array}{rl}y(t)& =\epsilon {y}^{(1)}(t)+{\epsilon}^{2}{y}^{(2)}(t)+{\epsilon}^{3}\hspace{0.17em}{y}^{(3)}(t)+\cdots \end{array}$$

3.3

$$\begin{array}{rl}\text{and}\phantom{\rule{2em}{0ex}}\theta (t)& =\epsilon {\theta}^{(1)}(t)+{\epsilon}^{2}{\theta}^{(2)}(t)+{\epsilon}^{3}{\theta}^{(3)}(t)+\cdots .\end{array}$$

3.4

From equation (2.13), it can be observed that $\dot{\theta}$ is independent of *θ*, and thus the ODE for *θ* can be written as

$$\dot{\theta}=\sum _{j=1}^{2}{\mathbf{\text{G}}}_{3j}(\mathit{\varphi}){\dot{\varphi}}_{j}=\epsilon \sum _{j=1}^{2}{\mathbf{\text{G}}}_{3j}(\mathit{\varphi}){\dot{s}}_{j}.$$

3.5

Using Taylor expansion we get

$$\begin{array}{rl}\dot{\theta}& =\sum _{j=1}^{2}[{\mathbf{\text{G}}}_{3j}(0)+({\varphi}_{1}{\frac{\mathrm{\partial}}{\mathrm{\partial}{\varphi}_{1}}|}_{(0,0)}+{\varphi}_{2}{\frac{\mathrm{\partial}}{\mathrm{\partial}{\varphi}_{2}}|}_{(0,0)}){\mathbf{\text{G}}}_{3j}\\ & \phantom{\rule{1em}{0ex}}+\hspace{0.17em}\frac{1}{2!}{({\varphi}_{1}{\frac{\mathrm{\partial}}{\mathrm{\partial}{\varphi}_{1}}|}_{(0,0)}+{\varphi}_{2}{\frac{\mathrm{\partial}}{\mathrm{\partial}{\varphi}_{2}}|}_{(0,0)})}^{2}{\mathbf{\text{G}}}_{3j}+\cdots ]{\dot{\varphi}}_{j}.\end{array}$$

3.6

Substituting the expansion for *θ* and *ϕ*_{i}(*t*)=*εs*_{i}(*t*):

$$\begin{array}{rl}\epsilon {\dot{\theta}}^{(1)}+{\epsilon}^{2}{\dot{\theta}}^{(2)}+{\epsilon}^{3}{\dot{\theta}}^{(3)}+\cdots & =\sum _{j=1}^{2}[{\mathbf{\text{G}}}_{3j}(0)+\epsilon ({s}_{1}{\frac{\mathrm{\partial}}{\mathrm{\partial}{\varphi}_{1}}|}_{(0,0)}+{s}_{2}{\frac{\mathrm{\partial}}{\mathrm{\partial}{\varphi}_{2}}|}_{(0,0)}){\mathbf{\text{G}}}_{3j}\\ & \phantom{\rule{1em}{0ex}}+\hspace{0.17em}{\epsilon}^{2}\frac{1}{2!}{({s}_{1}{\frac{\mathrm{\partial}}{\mathrm{\partial}{\varphi}_{1}}|}_{(0,0)}+{s}_{2}{\frac{\mathrm{\partial}}{\mathrm{\partial}{\varphi}_{2}}|}_{(0,0)})}^{2}{\mathbf{\text{G}}}_{3j}+\cdots ]\epsilon {\dot{s}}_{j}.\end{array}$$

3.7

Now we can collect terms of different orders in *ε*.

$$\begin{array}{rl}{\dot{\theta}}^{(1)}(t)& =\sum _{j=1}^{2}{\mathbf{\text{G}}}_{3j}(0){\dot{s}}_{j}(t),\end{array}$$

3.8

$$\begin{array}{rl}{\dot{\theta}}^{(2)}(t)& =\sum _{j=1}^{2}(\frac{\mathrm{\partial}{\mathbf{\text{G}}}_{3j}}{\mathrm{\partial}{\varphi}_{1}}(0){s}_{1}+\frac{\mathrm{\partial}{\mathbf{\text{G}}}_{3j}}{\mathrm{\partial}{\varphi}_{2}}(0){s}_{2}){\dot{s}}_{j}(t)\end{array}$$

3.9

$$\begin{array}{rl}\text{and}\phantom{\rule{2em}{0ex}}{\dot{\theta}}^{(3)}(t)& =\sum _{j=1}^{2}(\frac{1}{2}\frac{{\mathrm{\partial}}^{2}{\mathbf{\text{G}}}_{3j}}{\mathrm{\partial}{\varphi}_{1}^{2}}(0){s}_{1}^{2}+\frac{1}{2}\frac{{\mathrm{\partial}}^{2}{\mathbf{\text{G}}}_{3j}}{\mathrm{\partial}{\varphi}_{2}^{2}}(0){s}_{2}^{2}+\frac{{\mathrm{\partial}}^{2}{\mathbf{\text{G}}}_{3j}}{\mathrm{\partial}{\varphi}_{1}\mathrm{\partial}{\varphi}_{2}}(0){s}_{1}{s}_{2}){\dot{s}}_{j}(t).\end{array}$$

3.10

Explicit expressions for the derivatives are given in the electronic supplementary material. The first-order derivatives at the origin are zero as **G**_{31}, **G**_{32} are even functions in (*ϕ*_{1},*ϕ*_{2}) and so ${\dot{\theta}}^{(2)}=0$ in (3.9), implying that *θ*^{(2)}(*t*)=0 is zero (under zero initial conditions). Now, using the expression for *θ*(*t*) we can calculate expansions of the net motion **q**(*t*). First, we expand the matrices in (2.13) as

$$\begin{array}{rl}\mathbf{\text{D}}(\theta )& =\mathbf{\text{I}}+\theta \mathbf{\text{J}}+\frac{1}{2!}{\theta}^{2}{\mathbf{\text{J}}}^{2}+\frac{1}{3!}{\theta}^{3}{\mathbf{\text{J}}}^{3}+\cdots \\ \text{and}\phantom{\rule{2em}{0ex}}\mathbf{\text{J}}& =\left[\begin{array}{ccc}0& -1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right]\end{array}\}$$

3.11

and

$$\begin{array}{r}\mathbf{\text{G}}(\mathit{\varphi})=\mathbf{\text{G}}(0)+\epsilon ({s}_{1}{\frac{\mathrm{\partial}}{\mathrm{\partial}{\varphi}_{1}}|}_{(0,0)}+{s}_{2}{\frac{\mathrm{\partial}}{\mathrm{\partial}{\varphi}_{2}}|}_{(0,0)})\mathbf{\text{G}}+{\epsilon}^{2}\frac{1}{2!}{({s}_{1}{\frac{\mathrm{\partial}}{\mathrm{\partial}{\varphi}_{1}}|}_{(0,0)}+{s}_{2}{\frac{\mathrm{\partial}}{\mathrm{\partial}{\varphi}_{2}}|}_{(0,0)})}^{2}\mathbf{\text{G}}+\cdots .\end{array}$$

3.12

This, of course, is the same process done for the last row of **G** in (3.6). Expanding equation (2.13) by substituting expansions for *θ*(*t*), **D**(*θ*), **G**(** ϕ**) and also of

$$\begin{array}{rl}{\epsilon}^{1}:\phantom{\rule{1em}{0ex}}{\dot{\mathbf{\text{q}}}}^{(1)}& =\mathbf{\text{G}}(0)\dot{\mathbf{\text{s}}}(t)\\ {\epsilon}^{2}:\phantom{\rule{1em}{0ex}}{\dot{\mathbf{\text{q}}}}^{(2)}& =[{\theta}^{(1)}\mathbf{\text{J}}\mathbf{\text{G}}(0)+({s}_{1}{\frac{\mathrm{\partial}}{\mathrm{\partial}{\varphi}_{1}}|}_{(0,0)}+{s}_{2}{\frac{\mathrm{\partial}}{\mathrm{\partial}{\varphi}_{2}}|}_{(0,0)})\mathbf{\text{G}}]\dot{\mathbf{\text{s}}}(t)\\ {\epsilon}^{3}:\phantom{\rule{1em}{0ex}}{\dot{\mathbf{\text{q}}}}^{(3)}& =[\frac{1}{2}{\theta}^{(1)2}{\mathbf{\text{J}}}^{2}\mathbf{\text{G}}(0)+{\theta}^{(1)}\mathbf{\text{J}}({s}_{1}{\frac{\mathrm{\partial}}{\mathrm{\partial}{\varphi}_{1}}|}_{(0,0)}+{s}_{2}{\frac{\mathrm{\partial}}{\mathrm{\partial}{\varphi}_{2}}|}_{(0,0)})\mathbf{\text{G}}\\ & \phantom{\rule{1em}{0ex}}+\frac{1}{2}{({s}_{1}{\frac{\mathrm{\partial}}{\mathrm{\partial}{\varphi}_{1}}|}_{(0,0)}+{s}_{2}{\frac{\mathrm{\partial}}{\mathrm{\partial}{\varphi}_{2}}|}_{(0,0)})}^{2}\mathbf{\text{G}}]\dot{\mathbf{\text{s}}}(t)\\ {\epsilon}^{4}:\phantom{\rule{1em}{0ex}}{\dot{\mathbf{\text{q}}}}^{(4)}& =[({\theta}^{(3)}\mathbf{\text{J}}+\frac{1}{3!}{\theta}^{(1)3}{\mathbf{\text{J}}}^{3})\mathbf{\text{G}}(0)+\frac{1}{2}{\theta}^{(1)2}{\mathbf{\text{J}}}^{2}({s}_{1}{\frac{\mathrm{\partial}}{\mathrm{\partial}{\varphi}_{1}}|}_{(0,0)}+{s}_{2}{\frac{\mathrm{\partial}}{\mathrm{\partial}{\varphi}_{2}}|}_{(0,0)})\mathbf{\text{G}}\\ & \phantom{\rule{1em}{0ex}}+\frac{1}{2}{\theta}^{(1)}\mathbf{\text{J}}{({s}_{1}{\frac{\mathrm{\partial}}{\mathrm{\partial}{\varphi}_{1}}|}_{(0,0)}+{s}_{2}{\frac{\mathrm{\partial}}{\mathrm{\partial}{\varphi}_{2}}|}_{(0,0)})}^{2}\\ & \phantom{\rule{1em}{0ex}}+\frac{1}{6}{({s}_{1}{\frac{\mathrm{\partial}}{\mathrm{\partial}{\varphi}_{1}}|}_{(0,0)}+{s}_{2}{\frac{\mathrm{\partial}}{\mathrm{\partial}{\varphi}_{2}}|}_{(0,0)})}^{3}\mathbf{\text{G}}]\dot{\mathbf{\text{s}}}(t).\end{array}\}$$

3.13

Explicit expressions for the derivatives of elements of **G** are given in the electronic supplementary material. Owing to symmetry of the swimmer and of the gaits, it is known that there will only be net translation along *x*-direction, which is the axis of the centre link, while rotation and translation in *y*-direction are cancelled out [14]. Thus, expansions of net displacement *X* for specific gaits are calculated next.

Once the series expansion for the instantaneous velocities of the swimmer is derived, it is possible to obtain the net displacement *X* over a period for a specific gait via integration. The following proposition summarizes the results for the square and circular gaits.

*For a symmetric, three linked ‘Purcell swimmer’ performing a square or circular gait with amplitude*
*ε*, *the leading-order term and next-order correction for the displacement X over one full stroke in the direction of*

$$X={f}_{2}(\eta ){\epsilon}^{2}-{f}_{4}(\eta ){\epsilon}^{4}+O({\epsilon}^{6}).$$

3.14

*With*

$$\begin{array}{rl}{f}_{2}(\eta )& ={C}_{2}\eta l{(1-\eta )}^{3}(\eta +3)\\ \text{and}\phantom{\rule{2em}{0ex}}{f}_{4}(\eta )& ={C}_{4}\eta l{(1-\eta )}^{3}({\eta}^{7}+3{\eta}^{6}-10{\eta}^{5}-22{\eta}^{4}+29{\eta}^{3}+95{\eta}^{2}+44\eta +20),\end{array}\}$$

3.15

*where for the square gait given in* (2.21), *we have*
${C}_{2}=\frac{1}{4},\hspace{0.17em}{C}_{4}=\frac{1}{192}$
*and for the circular gait given in* (2.22) *C*_{2}=*π*/16, *C*_{4}=*π*/1024.

The proof of this proposition is given in the electronic supplementary material. It can easily be seen that the leading-order term of *X* is quadratic in *ε* and that the next-order corrections are of opposite sign to the leading-order terms. This implies that the displacement grows monotonically with *ε* for small stroke amplitudes, whereas for larger amplitudes with *ε*>1 the *O*(*ε*^{4})-term causes reversal in the direction of net motion. This, in turn, indicates the existence of a locally optimal value of *ε* that achieves maximal displacement. Note that, a leading-order term for displacement under the square gait is found in [21] using Lie brackets, and the results are identical to those given here.

As a demonstration of the utility of proposition 3.1, figure 3 shows plots of the displacement *X* versus stroke amplitude *ε* with link ratio of $\eta =\frac{1}{3}$ for both gaits (figure 3*a*, square, and and33*b*, circular). The solid curves are obtained from numerical integration of the nonlinear equations of motion (2.11), whereas the dashed and dash-dotted curves are *O*(*ε*^{2}) and *O*(*ε*^{4}) approximations, respectively. While the *O*(*ε*^{2})-approximation is monotonic in *ε* and works only for small stroke amplitudes, the *O*(*ε*^{4})-approximation with next-order correction captures the reversal in the direction of the displacement for intermediate amplitudes, though there is an increasing deviation from the numerical values for larger amplitudes, which is of *O*(*ξ*^{6}).

First, we discuss the dependence of net displacement *X* on the swimmer's geometric ratio *η* and derive its locally optimal values for both gaits, based on equations (3.14) and (3.15). As a demonstration, plots of *X* versus *η* under the square gait are shown in figure 4*a*,*b* for stroke amplitudes of *ε*=*π*/4 and *ε*=*π*/2, respectively. It can be seen that the *O*(*ε*^{2}) approximation (dashed) has a larger deviation from the exact value obtained by numerical integration (solid curves), compared with that of the *O*(*ε*^{4}) approximation (dash-dotted), where the deviations are further increased for larger amplitudes *ε*. Nevertheless, in all cases it is obvious that the displacement *X* vanishes at extreme cases of *η*→0 or *η*→1, where either the middle link or side links vanish, and that an optimal value of *η* exists that maximizes the displacement *X*.

Displacement over one cycle *X* as a function of *η* under the square gait with stroke amplitudes of (*a*) *ε*=*π*/4 and (*b*) *ε*=*π*/2. Solid curves, numerical integration. Dashed curves, *O*(*ε*^{2}) approximation. Dash-dotted **...**

We now obtain approximations of the optimal value of *η* based on the leading-order *O*(*ε*^{2}) terms in (3.14). Interestingly, it can be seen from (3.15) that the leading-order expressions for both square and circular gaits are identical up to multiplication by a constant. An important observation is that this relation is fairly general. That is, for any small-amplitude gait and any swimmer's dynamics with the same structure (i.e. not necessarily assuming resistive force theory), the leading-order term of *X* can always be decomposed into a function of swimmer's geometry multiplied by a function of the gait's unscaled trajectory. This key statement is summarized in the following theorem.

*Consider a swimmer with two shape variables whose dynamics can be written in the form of equation (2.13), under a small-amplitude gait*
*ϕ**(t)=ε***s***(t). The leading-order approximation of the swimmer's displacement X over one period can be written in the form: X*^{(2)}*=C(***s***)f(η), where C(***s***) depends on the shape of the unscaled trajectory and f(η) depends on the swimmer's geometric structure.*

The theorem, whose proof is detailed in the electronic supplementary material, implies that based on leading order approximation, the optimal geometry is independent of the gait's shape. For our particular swimmer model, the geometry-dependent function in (3.14) is *ηl*(1−*η*)^{3}(*η*+3). As expected, for the cases of *η*=0 and *η*=1, the net displacement vanishes, and the optimal value of *η* is easily obtained as ${\eta}^{\ast}=0.4\sqrt{10}-1=0.2649$, which results in maximal displacement of *X**=0.0859*ε*^{2} for the square gait and *X**=0.0675*ε*^{2} for the circular gait. The value of *η** is in agreement with the optimal geometry obtained in [21], but not with the one in [13]. This disagreement, which originates from differences in definitions and scaling, is discussed in the sequel.

Next, we use both *O*(*ε*^{2}) and *O*(*ε*^{4}) terms in (3.14) in order to obtain both optimal ratio *η* and amplitude *ε* for the two gaits. The polynomials *f*_{2}(*η*) and *f*_{4}(*η*) in (3.15) are the same for both gaits up to multiplication by a scalar. However, this is not true for all gaits. For example, it can be shown that for an elliptical gait trajectory, *f*_{4}(*η*) is a different polynomial. For a given value of *η*, (3.14) implies that the locally optimal amplitude *ε** for maximal displacement is given by

$${\epsilon}^{\ast}=\sqrt{\frac{{f}_{2}(\eta )}{2\hspace{0.17em}{f}_{4}(\eta )}}.$$

3.16

Substituting into (3.14) one obtains the optimal displacement as

$${X}^{\ast}=\frac{{f}_{2}(\eta )}{4\hspace{0.17em}{f}_{4}(\eta )}.$$

3.17

This implies the existence of an optimal combination of link ratio *η* and amplitude *ε* that achieve a local maximum of the displacement. Differentiating (3.17) with respect to *η*, the optimal value is obtained by finding the roots of a ninth order polynomial in *η*. As the expression for *X** in (3.17) is the same for both square and circular gaits up to a multiplicative constant, the optimal geometry for both gaits is obtained as *η**=0.1784. Substituting into (3.16), optimal amplitudes *ε* and displacements *X* are obtained for each gait and are presented in table 1. Numerically calculating the optimal geometry *η* and stroke amplitude *ε* using Matlab's function **fmincon** gives the results also presented in table 1.

The locally optimal amplitude *ε** based on *O*(*ε*^{4}) approximation is quite close to the exact value obtained numerically. On the other hand, the large difference in *η** can be explained by the large deviation in *X* for *ε*≈*π*/2, as seen in figure 4*b*. The figure also shows that the *O*(*ε*^{4}) expression for *X* overestimates the exact value. Finally, it is important to note that all optimization results discussed here are limited to finding local rather than global maxima. In fact, continuing to higher order terms in the expansions of (3.14) reveals a minimum of *X*<0 for *ε*≈3*rad* with larger absolute value than the first positive maximum, and even higher extremum points for larger non-physical values of *ε* [30].

As mentioned above, the work of Becker *et al.* [13] obtained a different value of optimal geometry *η*. The explanation for this difference is twofold: first, Becker *et al.* [13] considered maximization of the mean forward velocity under a constant torque at the joint rather than net displacement. Second, they [13] scaled swimming distance by the side link's length *l*_{1} rather than the total length *l*.

We now show that by adopting the definitions of [13] and using our leading-order expressions in proposition 3.1, one obtains optimal geometry that agrees with [13]. First, we consider the case where a constant torque is applied on the rotating joint in each part of the square gait. The relation between the joint torques and the joint velocities is given in (2.16). For the first quarter of the square gait we have that ${\dot{\varphi}}_{1}=0$ and so the second row of (2.16) is reduced to ${\tau}_{2}={W}_{22}{\dot{\varphi}}_{2}$, with *W*_{ij} the elements of **W**(** ϕ**) defined in (2.19). Assuming a constant torque of

$$\begin{array}{rl}{T}_{\tau}& =\int \mathrm{d}t=4{\int}_{\epsilon}^{-\epsilon}\frac{\mathrm{d}t}{\mathrm{d}{\varphi}_{2}}\hspace{0.17em}\mathrm{d}{\varphi}_{2}=4{\int}_{\epsilon}^{-\epsilon}\frac{{W}_{22}(\mathit{\varphi})}{{\tau}_{o}}\hspace{0.17em}\mathrm{d}{\varphi}_{2}\\ & =4{\int}_{\epsilon}^{-\epsilon}\frac{{W}_{22}(0)}{{\tau}_{o}}\hspace{0.17em}\mathrm{d}{\varphi}_{2}+O({\epsilon}^{3})=\frac{1}{12}{l}^{3}{(1-{\eta}^{2})}^{3}\frac{{\mathrm{c}}_{\mathrm{t}}}{{\tau}_{o}}\epsilon +O({\epsilon}^{3}).\end{array}$$

3.18

The mean forward velocity *V*
_{τ}=*X*/*T*_{τ} can now be computed to leading order as

$${V}_{\tau}^{(1)}=\frac{{X}^{(2)}}{{T}_{\tau}^{(1)}}=\epsilon \frac{3\eta (\eta +3)}{{l}^{2}{(\eta +1)}^{3}}\frac{{\tau}_{o}}{{\mathrm{c}}_{\mathrm{t}}}.$$

3.19

From (3.19), the optimal geometry that maximizes the leading-order mean velocity *V*
^{(1)} is obtained as *η**=0.646, which is fundamentally different from the optimal geometry for maximal displacement *X* obtained above.

Becker *et al.* [13] also used a different definition of geometric ratio, denoted here as *η*_{b}=*l*_{0}/*l*_{1}. That is, the length is normalized by that of the side links. The relation between the two definitions of *η* is given by

$${\eta}_{b}=\frac{{l}_{0}}{{l}_{1}}=\frac{2\eta}{1-\eta}.$$

3.20

Leading-order expressions for the time of quarter period under constant torque *T*_{τ} and the mean velocity *V* are given in [13] as (after multiplying *T*_{τ} by four and adapting some notation)

$${T}_{\tau}=[\frac{16}{3}\frac{{({\eta}_{b}+1)}^{3}{l}_{1}^{3}}{({\eta}_{b}+4){({\eta}_{b}+1)}^{2}+3{\eta}_{b}+4}+O({\epsilon}^{2})]\left(\frac{{\mathrm{c}}_{\mathrm{t}}{l}_{0}^{3}}{{\tau}_{o}}\right)$$

3.21

and

$${V}_{\tau}=[\frac{3}{4}\frac{{\eta}_{b}(2{\eta}_{b}+3)}{{({\eta}_{b}+1)}^{3}({\eta}_{b}+2){l}_{1}^{2}}+O({\epsilon}^{2})]\left(\frac{{\tau}_{o}}{{\mathrm{c}}_{\mathrm{t}}{l}_{0}^{2}}\right).$$

3.22

Substituting the transformation (3.20), it can be shown that these expressions agree with our derivations in (3.18), (3.19), and also with our leading order expression for net displacement *X*=*V*
_{τ}*T*_{τ} as given in (3.14). From (3.22), optimal geometry for maximizing *V*
_{τ} has been obtained in [13] as ${\eta}_{b}^{\ast}=0.54$, which corresponds to *η*=0.213. The difference from the optimal value of *η*=0.646 obtained by maximizing *V*
_{τ} in (3.19) is now due to the fact that (3.19) is scaled by *l* while (3.22) is scaled by *l*_{1}.

We now compute series expansions for the period time *T*_{p} under constant power *P*_{0}=1 as formulated in (2.28), for both square and circular gaits. These expressions are used in the next section for obtaining approximate expressions of Lighthill's energetic efficiency according to equation (2.29). The period time is expanded into a power series as

$${T}_{\mathrm{p}}=\epsilon {T}_{\mathrm{p}}^{(1)}+{\epsilon}^{2}{T}_{\mathrm{p}}^{(2)}+{\epsilon}^{3}{T}_{\mathrm{p}}^{(3)}+\cdots .$$

3.23

Explicit expressions for the different terms ${T}_{\mathrm{p}}^{(i)}$ depend on the choice of the gait. For the square and circular gaits, these expressions are summarized in the following proposition.

*For a symmetric, three linked ‘Purcell swimmer’ performing a square or circular gait with amplitude*
*ε*, *the period times of one full stroke with constant power of*
*P*_{o}=1 *exerted by the joints are expanded as follows*.

*Square gait*:

$$\begin{array}{r}{T}_{\mathrm{p}}=\frac{\sqrt{6}}{3}\sqrt{{\mathrm{c}}_{\mathrm{t}}{l}^{3}{({1-\eta}^{2})}^{3}}\epsilon +\frac{\sqrt{6{\mathrm{c}}_{\mathrm{t}}{l}^{3}}{(1-\eta )}^{4}(\eta +1)({\eta}^{3}+3{\eta}^{2}-3\eta +1)}{12\sqrt{{(1-{\eta}^{2})}^{3}}}{\epsilon}^{3}+O({\epsilon}^{5}).\end{array}$$

3.24

*Circular gait*:

$${T}_{\mathrm{p}}=E\left(\frac{{\eta}^{3}+3{\eta}^{2}-3\eta -1}{{\eta}^{2}(\eta +3)}\right)\sqrt{\frac{{c}_{\mathrm{t}}{l}^{3}}{3}{\eta}^{2}{(1-\eta )}^{3}(\eta +3)}\epsilon +O({\epsilon}^{3}).$$

3.25

The function *E*() in (3.25) denotes a complete elliptic integral of the second kind, whose definition is reviewed in the electronic supplementary material. The proof of this proposition is also given in the electronic supplementary material.

In order to validate the expansions of *T*_{p} in (3.24) and (3.25), figure 5 plots the *O*(*ε*) and *O*(*ε*^{3}) approximations of *T*_{p} as a function of stroke amplitude *ε*, compared with the exact value of *T*_{p} obtained by numerical integration of (2.28), for both gaits (figure 5*a*, square, and and55*b*, circular) under equal links length, $\eta =\frac{1}{3}$. It can be seen that for small amplitudes, the period time *T*_{p} indeed scales linearly with *ε* and that for intermediate values of *ε* the *O*(*ε*^{3}) correction term slightly improves the approximation accuracy (square gait only).

Using the expressions found in (3.24),(3.25) for the period times under constant power, along with previous results of net displacement (3.14), we now calculate several approximations of Lighthill's energetic efficiency *ξ* using (2.29). Then we obtain optimal values of length ratio *η* and stroke amplitude *ε* for maximal energetic efficiency, under both square and circular gaits. Using the expansions for *X* and *T*_{p}, the efficiency *ξ* can be written as

$$\xi =\frac{{X}^{2}}{{T}_{\mathrm{p}}^{2}}={\epsilon}^{2}\frac{{({X}^{(2)}+{\epsilon}^{2}{X}^{(4)}+\cdots )}^{2}}{{({T}_{\mathrm{p}}^{(1)}+{\epsilon}^{2}{T}_{\mathrm{p}}^{(3)}+\cdots )}^{2}}.$$

4.1

Calculating Taylor expansion of (4.1) in powers of *ε*, one obtains

$$\xi ={\epsilon}^{2}\frac{{({X}^{(2)})}^{2}}{{({T}_{\mathrm{p}}^{(1)})}^{2}}+2{\epsilon}^{4}{X}^{(2)}\frac{{T}_{\mathrm{p}}^{(1)}{X}^{(4)}-{X}^{(2)}{T}_{\mathrm{p}}^{(3)}}{{({T}_{\mathrm{p}}^{(1)})}^{3}}+O({\epsilon}^{6}).$$

4.2

Explicit approximations for the efficiency *ξ* under square and circular gaits are discussed below.

Substituting (3.14), (3.15) and (3.24) into (4.2), the first two elements in the expansion of Lighthill's energetic efficiency under the square gait are given by

$${\xi}^{(2)}=\frac{{({X}^{(2)})}^{2}}{{({T}_{\mathrm{p}}^{(1)})}^{2}}=\frac{3{\eta}^{2}{(1-\eta )}^{3}{(\eta +3)}^{2}}{32{\mathrm{c}}_{\mathrm{t}}l{(\eta +1)}^{3}}$$

4.3

and

$$\begin{array}{r}{\xi}^{(4)}=\frac{{\eta}^{2}{\epsilon}^{4}{(\eta -1)}^{3}(\eta +3)({\eta}^{9}+5{\eta}^{8}-3{\eta}^{7}-39{\eta}^{6}-37{\eta}^{5}+71{\eta}^{4}+263{\eta}^{3}+371{\eta}^{2}-48\eta +56)}{256{c}_{\mathrm{t}}l{(\eta +1)}^{5}}.\end{array}$$

4.4

Figure 6*a* plots the *O*(*ε*^{2}) and *O*(*ε*^{4}) approximations of *ξ* as a function of *ε* in dashed and dash-dotted curves, respectively, under the square gait with length ratio of $\eta =\frac{1}{3}$. For comparison, the exact value of *ξ* obtained from numerical integration is plotted in solid curve. While the *O*(*ε*^{2}) approximation grows monotonically with *ε*, adding the next-order correction of *O*(*ε*^{4}) also captures a local maximum in the efficiency, which is obtained at an intermediate stroke amplitude of *ε*≈1. Nevertheless, this approximation of *ξ* becomes negative for larger values of *ε*, which is non-physical according to the definition of *ξ* in (2.29). A more reasonable approximation is obtained by plugging the *O*(*ε*^{4}) and *O*(*ε*^{3}) approximations for *X* and *T*_{p}, respectively, into (2.29). With a slight abuse of notation, this approximation is denoted here as *O*(*ε*^{4/3}), and is given by

$${\xi}^{(4/3)}=\frac{{({\epsilon}^{2}{X}^{(2)}+{\epsilon}^{4}{X}^{(4)})}^{2}}{{(\epsilon {T}_{\mathrm{p}}^{(1)}+{\epsilon}^{3}{T}_{\mathrm{p}}^{(3)})}^{2}}.$$

4.5

Note that this approximation is different from the direct expansion of (4.1) into (4.2), as demonstrated below. A plot of this approximation is overlayed as a dotted line in figure 6*a*. It can be seen that the *O*(*ε*^{4/3}) approximation is always positive, and correctly captures the general trend of *ξ* as a function of *ε*. Finally, it can be seen from the plot that the exact value of *ξ* attains a higher maximum value at *ε*≈3, where the direction of swimming is reversed. This second maximum is not captured by any of the approximate expressions mentioned above. Nevertheless, as already mentioned in [13], this maximum at large strokes where the joint angles approach ±*π* is often impractical, due to possible collisions between the links. Moreover, resistive force theory which assumes negligible hydrodynamic interaction between the links is no longer valid in this range [30]. Thus, we focus here on approximations of the first maximum of *ξ* which is attained at moderate amplitudes, as discussed next.

Lighthill's energetic efficiency *ξ* under the square gait: (*a*) plot of *ξ* versus *ε* for $\eta =\frac{1}{3}$. (*b*) Plot of *ξ* versus *η* for *ε*=1. Solid curves, exact (numeric) computation. Dashed, dash-dotted and dotted **...**

We now use the approximations of Lighthill's energetic efficiency in order to obtain locally optimal values of both length ratio *η* and stroke amplitude *ε*. Using only the leading-order approximation of the efficiency in (4.3) gives an optimal length ratio of *η**=0.327. Taking also the next-order correction term into account and implementing the same calculation steps shown above in §3c, the optimal geometry and amplitude are obtained as

$${\eta}^{\ast}=0.2853,\phantom{\rule{1em}{0ex}}{\epsilon}^{\ast}=0.908\hspace{0.17em}\text{rad}={52}^{\circ}\phantom{\rule{1em}{0ex}}\mathrm{and}\phantom{\rule{1em}{0ex}}\xi =\mathrm{0.5838.}$$

4.6

Using *O*(*ε*^{(4/3)}) approximation, the optimal geometry and amplitude are obtained by solving a system of two polynomial equations as

$${\eta}^{\ast}=0.2754,\phantom{\rule{1em}{0ex}}{\epsilon}^{\ast}=1.1038\hspace{0.17em}\text{rad}={63.2}^{\circ}\phantom{\rule{1em}{0ex}}\mathrm{and}\phantom{\rule{1em}{0ex}}\xi =\mathrm{0.7294.}$$

4.7

Finally, numerical calculation of *ξ* under the square gait and conducting optimization using Matlab's function **fmincon** gives optimal values of

$${\eta}^{\ast}=0.2879,\phantom{\rule{1em}{0ex}}{\epsilon}^{\ast}=1.133\hspace{0.17em}\text{rad}={64.9}^{\circ}\phantom{\rule{1em}{0ex}}\mathrm{and}\phantom{\rule{1em}{0ex}}\xi =\mathrm{0.7709.}$$

4.8

Figure 6*b* shows a plot of the approximations of *ξ* as a function of *η* for a large amplitude of *ε*=1, compared with the exact computation of *ξ* obtained numerically. While there are large discrepancies in the value of the efficiency *ξ* (best captured by the *O*(*ε*^{4/3}) approximation), all approximations predict optimal values around *η*=0.28.

For the circular gait, the constant-power period time *T*_{p} has been approximated only to first order *O*(*ε*) in (3.25). Using (3.14), the expansion of Lighthill's efficiency in (4.2) can thus be obtained only to leading order, as

$$\xi ={\epsilon}^{2}\frac{3{(1-\eta )}^{3}(\eta +3)}{{c}_{t}lE\left(\frac{{\eta}^{3}+3{\eta}^{2}-3\eta -1}{{\eta}^{2}(\eta +3)}\right)}+O({\epsilon}^{3}).$$

4.9

Similar to (4.5), we define the *O*(*ε*^{4/1}) approximation as

$${\xi}^{(4/1)}=\frac{{({\epsilon}^{2}{X}^{(2)+}{\epsilon}^{4}{X}^{(4)})}^{2}}{{({T}_{\mathrm{p}}^{(1)})}^{2}}.$$

4.10

Substituting (3.14) and (3.25) into (4.10) gives the expression for the *O*(*ε*^{4/1}) approximation (not shown for brevity). The efficiency *ξ* and its approximations are plotted in figure 7*a* as a function of the stroke amplitude *ε* under equal links lengths $\eta =\frac{1}{3}$, and in figure 7*b* as a function of *η* under amplitude of *ε*=1. In both plots, *O*(*ε*^{2}) approximations appear in dashed curves, *O*(*ε*^{4/1}) approximations appear in dotted curves, and the exact values computed numerically are shown in solid curves. As before, the leading-order *O*(*ε*^{2}) approximation does not capture the optimum with respect to amplitude but both approximations capture optimum with respect to *η*.

Lighthill's energetic efficiency *ξ* under the circular gait: (*a*) plot of *ξ* versus *ε* for $\eta =\frac{1}{3}$. (*b*) Plot of *ξ* versus *η* for *ε*=1. Solid curves, exact (numeric) computation. Dashed and dotted curves, approximations **...**

Using the leading-order approximation for *ξ* in (4.9) and numerically searching for maximum with respect to *η*, optimal geometry is obtained as *η**=0.3139, with an efficiency of *ξ*^{(2)}=1.5565. Using the *O*(*ε*^{4/1}) approximation in (4.10), optimal values of both *ε* and *η* are obtained numerically as

$${\eta}^{\ast}=0.2379,\phantom{\rule{1em}{0ex}}{\epsilon}^{\ast}=1.3821\hspace{0.17em}\text{rad}\phantom{\rule{1em}{0ex}}\mathrm{and}\phantom{\rule{1em}{0ex}}\xi =1.22.$$

4.11

Finally, numerical calculation of *ξ* under the circular gait and conducting optimization using Matlab's **fmincon** function give optimal values of

$${\eta}^{\ast}=0.2684,\phantom{\rule{1em}{0ex}}{\epsilon}^{\ast}=1.3747\hspace{0.17em}\text{rad}\phantom{\rule{1em}{0ex}}\mathrm{and}\phantom{\rule{1em}{0ex}}\xi =1.2777.$$

4.12

Again, the *O*(*ε*^{4/1}) expression in (4.10) achieves a reasonable approximation of the optimum, as also seen from figure 7*a*,*b*.

An important observation is that even in leading-order *O*(*ε*^{2}) approximation of Lighthill's energetic efficiency *ξ*, the optimal geometry *η** for maximizing *ξ* actually depends on the shape of the chosen gait trajectory. This can clearly be seen in the different dependence of *ξ*^{(2)} on *η* in (4.3) and (4.9). This is a substantial difference from the case of maximizing displacement *X*, where dependence of the leading-order approximation in *η* is independent of the gait's shape, as manifested in theorem 3.2.

In this work, we have analysed the motion of the three-link ‘Purcell swimmer’. We provided a systematic method for deriving an expansion of the velocities of the swimmer and presented leading-order expressions and next-order corrections for the net displacement over one period in the cases of a square and circular gait. Examination of the correction terms confirms that there is a reversal in the direction of the net displacement at high amplitudes, a result which has previously only been shown numerically. The gait amplitude and swimmer geometry that optimize the displacement over one period were approximated using the first two terms in the expansion. Additionally, by writing asymptotic expressions for the period time under constant power expenditure, we were able to write, for the first time, leading-order term for the energetic efficiency of the square and circular gaits as well as a next-order correction term for the square gait. Once again, we used the obtained expressions in order to find the energetically optimal gait amplitude and swimmer geometry. The results demonstrate the utility of perturbation methods for obtaining approximate explicit expressions for nonlinear dynamics of locomotion systems, which enable analysis and optimization of their performance.

We now briefly discuss some limitations of this work and list possible directions for future extension of the research. First, the swimmer's dynamic equations have been formulated using the simplification of resistive force theory [26,27], in which hydrodynamic interaction between the links is neglected. It is well known (cf. [15,31]) that this assumption holds only for highly slender links, and for small stroke amplitudes where the gap between the links remains large even in the vicinity of the joints. Hydrodynamic interactions may be accounted for by using more refined models of slender body theory as in [30,18]. The decoupled relations in (2.5) should then be modified to account for inter-link resistance, and the analysis will probably become purely numerical due to the added complexity of the dynamics. Nevertheless, the structure and geometric symmetries of the dynamics in (2.11) will be maintained without qualitative changes. In fact, according to the numerical investigation in [30,18], no large quantitative changes in the results are expected even for moderate stroke amplitudes (*ε*≈2*rad*).

Second, the gait optimization conducted in this work was limited to varying the amplitude of predefined shape trajectories and did not address the possibility of other trajectories which may have better performance. Shape optimization of swimmers has been extensively studied, cf. [32–34]. One possible approach [18,35,36] is to represent a subset of all possible periodic shape changes by a finite set of variables (e.g. coefficients of a truncated Fourier series) and perform numerical optimization over this discrete set of variables. Another approach is optimal control theory [37], which is based on calculus of variations. This approach has been exploited by Alouges *et al.* [38] for obtaining energy-optimal gaits of unidirectional axisymmetric swimmers and by Giraldi *et al.* [21] for gait optimization of multi-link microswimmers under input constraints. In our recent work [39], we have used optimal control in order to obtain unconstrained optimal gaits for maximal displacement of Purcell's swimmer, which exactly reproduce the gaits obtained numerically in [18]. Furthermore, we are currently working on extending the results of [38] to planar motion of multi-link microswimmers.

Third, experiments conducted on a macro-scale swimmer prototype in a highly viscous fluid [14] have recently demonstrated the symmetry properties of Purcell's swimmer performing certain symmetric gaits. Similar experiments can be done to verify some results from this paper, e.g. the dependency of net displacement on swimmer geometry and gait amplitude. Finally, in many practical situations, the actuation of robotic swimmers will likely be by applying torques at the joints rather than prescribing the joint angles directly [40]. This calls for the formulation of an asymptotic expansion of the relation between the joint torques and the joint angles.

Both authors contributed equally to the writing of this paper and both authors gave final approval for publication.

We declare we have no competing interests.

This work is supported by the Israel Science Foundation (ISF) under grant no. 567/14.

1. Lauga E, Powers TR
2009.
The hydrodynamics of swimming microorganisms. *Rep. Prog. Phys.*
72, 096601 (doi:10.1088/0034-4885/72/9/096601)

2. Lauga E.2016. Bacterial hydrodynamics. *Annu. Rev. Fluid Mech.* (http://arxiv.org/abs/1509.02184. )

3. Abbott JJ, Lagomarsino MC, Zhang L, Dong L, Nelson BJ
2009.
How should microrobots swim?
*Int. J. Robot. Res.*
28, 1434–1447. (doi:10.1177/0278364909341658)

4. Kósa G, Jakab P, Hata N, Jólesz F, Neubach Z, Shoham M, Zaaroor M, Székely G
2008.
Flagellar swimming for medical micro robots: theory, experiments and application. In *2nd IEEE RAS & EMBS Int. Conf. on Biomedical Robotics and Biomechatronics (BioRob)*, pp. 258–263. New York, NY: IEEE.

5. Jang B.
*et al.*
2015.
Undulatory locomotion of magnetic multilink nanoswimmers. *Nano Lett.*
15, 4829–4833. (doi:10.1021/acs.nanolett.5b01981) [PubMed]

6. Gao W.
*et al.*
2012.
Cargo-towing fuel-free magnetic nanoswimmers for targeted drug delivery. *Small*
8, 460–467. (doi:10.1002/smll.201101909) [PubMed]

7. Peyer KE, Zhang L, Nelson BJ
2013.
Bio-inspired magnetic swimming microrobots for biomedical applications. *Nanoscale*
5, 1259–1272. (doi:10.1039/C2NR32554C) [PubMed]

8. Sitti M, Ceylan H, Hu W, Giltinan J, Turan M, Yim S, Diller E
2015.
Biomedical applications of untethered mobile milli/microrobots. *Proc. IEEE*
103, 205–224. (doi:10.1109/JPROC.2014.2385105)

9. Happel J, Brenner H
1965.
*Low Reynolds number hydrodynamics*. Englewood Cliffs, NJ: Prentice-Hall.

10. Purcell EM.
1977.
Life at low Reynolds number. *Am. J. Phys.*
45, 3–11. (doi:10.1119/1.10903)

11. Berg HC, Anderson RA
1973.
Bacteria swim by rotating their flagellar filaments. *Nature*
245, 380 (doi:10.1038/245380a0) [PubMed]

12. Taylor G.
1951.
Analysis of the swimming of microscopic organisms. *Proc. R. Soc. Lond. A*
209, 447–461. (doi:10.1098/rspa.1951.0218)

13. Becker LE, Koehler SA, Stone HA
2003.
On self-propulsion of micro-machines at low Reynolds number: Purcell's three-link swimmer. *J. Fluid Mech.*
490, 15–35. (doi:10.1017/S0022112003005184)

14. Gutman E, Or Y
2016.
Symmetries and gaits for Purcell's three-link microswimmer model. *IEEE Trans. Robot.*
32, 53–69. (doi:10.1109/TRO.2015.2500442)

15. Avron JE, Raz O
2008.
A geometric theory of swimming: Purcell's swimmer and its symmetrized cousin. *New J. Phys.*
10, 063016 (doi:10.1088/1367-2630/10/6/063016)

16. Hatton RL, Ding Y, Choset H, Goldman DI
2013.
Geometric visualization of self-propulsion in a complex medium. *Phys. Rev. Lett.*
110, 078101 (doi:10.1103/PhysRevLett.110.078101) [PubMed]

17. Hatton RL, Choset H
2013.
Geometric swimming at low and high Reynolds numbers. *IEEE Trans. Robot.*
29, 615–624. (doi:10.1109/TRO.2013.2251211)

18. Tam D, Hosoi AE
2007.
Optimal stroke patterns for Purcell's three-link swimmer. *Phys. Rev. Lett.*
98, 068105 (doi:10.1103/PhysRevLett.98.068105) [PubMed]

19. Shapere A, Wilczek F
1989.
Geometry of self-propulsion at low Reynolds number. *J. Fluid Mech.*
198, 557–585. (doi:10.1017/S002211208900025X)

20. Nayfeh AH.
2008.
*Perturbation methods*. New York, NY: John Wiley & Sons.

21. Giraldi L, Martinon P, Zoppello M
2015.
Optimal design of Purcell's three-link swimmer. *Phys. Rev. E*
91, 023012 (doi:10.1103/PhysRevE.91.023012) [PubMed]

22. Passov(Gutman) E, Or Y
2012.
Dynamics of Purcell's three-link microswimmer with a passive elastic tail. *Eur. Phys. J. E*
35, 1–9. (doi:10.1140/epje/i2012-12001-6) [PubMed]

23. Gutman E, Or Y
2014.
Simple model of a planar undulating magnetic microswimmer. *Phys. Rev. E*
90, 013012 (doi:10.1103/PhysRevE.90.013012) [PubMed]

24. Gutman E, Or Y
2016.
Optimizing an undulating magnetic microswimmer for cargo towing. *Phys. Rev. E*
93, 063105 (doi:10.1103/PhysRevE.93.063105) [PubMed]

25. Man Y, Lauga E
2013.
The wobbling-to-swimming transition of rotated helices. *Phys. Fluids*
25, 071904 (doi:10.1063/1.4812637)

26. Cox RG.
1970.
The motion of long slender bodies in a viscous fluid part 1. General theory. *J. Fluid Mech.*
44, 791–810. (doi:10.1017/S002211207000215X)

27. Gray J, Hancock G
1955.
The propulsion of sea-urchin spermatozoa. *J. Exp. Biol.*
32, 802–814.

28. Childress S.
2012.
A thermodynamic efficiency for Stokesian swimming. *J. Fluid Mech.*
705, 77–97. (doi:10.1017/jfm.2011.561)

29. Lighthill J.
1975.
*Mathematical biofluiddynamics*. Philadelphia, PA: SIAM.

30. Huber G, Koehler SA, Yang J
2011.
Micro-swimmers with hydrodynamic interactions. *Math. Comput. Model.*
53, 1518–1526. (doi:10.1016/j.mcm.2010.04.002)

31. Alouges F, DeSimone A, Heltai L, Lefebvre-Lepot A, Merlet B
2013.
Optimally swimming stokesian robots. *Discrete Continuous Dyn. Syst. B*
18, 1189–1215. (doi:10.3934/dcdsb.2013.18.1189)

32. Delfour MC, Zolésio J-P
2011.
*Shapes and geometries: analysis, differential calculus, and optimization*, 2nd edn
Advances in design and control, vol. 4
Philadelphia, PA: SIAM.

33. Moubachir M, Zolésio J-P
2006.
*Moving shape analysis and control: applications to fluid structure interactions*. Pure and applied mathematics, vol. 277
London, UK: Chapman and Hall.

34. Walker SW.
2015.
*The shapes of things: a practical guide to differential geometry and the shape derivative*. Advances in design and control, vol. 28
Philadelphia, PA: SIAM

35. Avron JE, Gat O, Kenneth O
2004.
Optimal swimming at low Reynolds numbers. *Phys. Rev. Lett.*
93, 186001 (doi:10.1103/PhysRevLett.93.186001) [PubMed]

36. Shapere A, Wilczek F
1989.
Efficiencies of self-propulsion at low Reynolds number. *J. Fluid Mech.*
198, 587–599. (doi:10.1017/S0022112089000261)

37. Bryson A, Ho Y-C
1975.
*Applied optimal control: optimization, estimation, and control*. Washington DC: Hemisphere Pub. Corp.

38. Alouges F, DeSimone A, Lefebvre A
2009.
Optimal strokes for axisymmetric microswimmers. *Eur. Phys. J. E*
28, 279–284. (doi:10.1140/epje/i2008-10406-4) [PubMed]

39. Wiezel O, Or Y
In preparation
Using optimal control to obtain maximum displacement gait for Purcell's three-link swimmer. In *Proc. 2016 IEEE Conf. on Decision and Control (CDC)*.

40. Or Y.
2012.
Asymmetry and stability of shape kinematics in microswimmers’ motion. *Phys. Rev. Lett.*
108, 258101 (doi:10.1103/PhysRevLett.108.258101) [PubMed]

Articles from Proceedings. Mathematical, Physical, and Engineering Sciences are provided here courtesy of **The Royal Society**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |