<figure><img src="image_1.png" /></figure>

HAL
open science

# Energy-optimal strokes for multi-link microswimmers: Purcell's loops and Taylor's waves reconciled

Oren Wiezel, Laetitia Giraldi, Antonio Desimone, Yizhar Or, François Alouges

► To cite this version:

Oren Wiezel, Laetitia Giraldi, Antonio Desimone, Yizhar Or, François Alouges. Energy-optimal strokes for multi-link microswimmers: Purcell's loops and Taylor's waves reconciled. 2024. hal-01970740

HAL Id: hal-01970740

https://inria.hal.science/hal-01970740

Preprint submitted on 16 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
PAPER • OPEN ACCESS

# Energy-optimal strokes for multi-link microswimmers: Purcell's loops and Taylor's waves reconciled

To cite this article: François Alouges et al 2019 *New J. Phys.* 21 043050

View the [article online](#) for updates and enhancements.

You may also like

- Physics of microswimmers—single particle motion and collective behavior: a review
J Elgeti, R G Winkler and G Gompper

- Navigation control of flagellated magnetic microswimmer by parametric excitation
T Sonamani Singh and R D S Yadava

- Simulation of microswimmer hydrodynamics with multiparticle collision dynamics
Andreas Zöttl
PAPER • OPEN ACCESS

# Energy-optimal strokes for multi-link microswimmers: Purcell's loops and Taylor's waves reconciled

To cite this article: François Alouges et al 2019 *New J. Phys.* 21 043050

View the [article online](#) for updates and enhancements.

You may also like

- Physics of microswimmers—single particle motion and collective behavior: a review
J Elgeti, R G Winkler and G Gompper

- Navigation control of flagellated magnetic microswimmer by parametric excitation
T Sonamani Singh and R D S Yadava

- Simulation of microswimmer hydrodynamics with multiparticle collision dynamics
Andreas Zöttl
<figure><img src="image_1.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED

23 January 2019

ACCEPTED FOR PUBLICATION

20 March 2019

PUBLISHED

29 April 2019

Original content from this
work may be used under
the terms of the Creative
Commons Attribution 3.0
licence.

Any further distribution of
this work must maintain
attribution to the
author(s) and the title of
the work, journal citation
and DOI.

<figure><img src="image_2.png" /></figure>

PAPER

# Energy-optimal strokes for multi-link microswimmers: Purcell's loops and Taylor's waves reconciled

François Alouges¹, Antonio DeSimone² Ⓓ, Laetitia Giraldi³ Ⓓ, Yizhar Or⁴,⁵ Ⓓ and Oren Wiezel⁴ Ⓓ

¹ CMAP, Ecole Polytechnique, Cedex, France

² The BioRobotics Institute, Scuola Superiore Sant'Anna, Pisa, Italy

³ Université Côte d'Azur, Inria, CNRS, LJAD, France

⁴ Faculty of Mechanical EngineeringTechnion - Israel Institute of Technology, Israel

⁵ Author to whom any correspondence should be addressed

E-mail: izi@technion.ac.il

**Keywords:** microswimmer dynamics, optimal strokes, energy efficiency, asymptotic analysis

## Abstract

Micron-scale swimmers move in the realm of negligible inertia, dominated by viscous drag forces. In this paper, we formulate the leading-order dynamics of a slender multi-link (N-link) microswimmer assuming small-amplitude undulations about its straight configuration. The energy-optimal stroke to achieve a given prescribed displacement in a given time period is obtained as the largest eigenvalue solution of a constrained optimal control problem. Remarkably, the optimal stroke is an ellipse lying within a two-dimensional plane in the (N-1)-dimensional space of joint angles, where N can be arbitrarily large. For large N, the optimal stroke is a traveling wave of bending, modulo edge effects. If the number of shape variables is small, we can consider the same problem when the prescribed displacement in one time period is large, and not attainable with small variations of the joint angles. The fully nonlinear optimal control problem is solved numerically for the cases N = 3 (Purcell's three-link swimmer) and N = 5 showing that, as the prescribed displacement becomes small, the optimal solutions obtained using the small-amplitude assumption are recovered. We also show that, when the prescribed displacements become large, the picture is different. For N = 3 we recover the non-convex planar loops already known from previous studies. For N = 5 we obtain non-planar loops, raising the question of characterizing the geometry of complex high-dimensional loops.

## 1. Introduction

The analysis of biological and bio-inspired swimming at microscopic scales has attracted considerable attention in the recent literature, starting from the seminal work by Taylor [1], Lighthill [2], and Purcell [3]. One of the reasons is that swimming of unicellular organisms is at the root of many fundamental processes in biology: reproduction through the swimming of sperm cells is just one example [4–7]. Moreover, biology inspires the design of bio-mimetic artificial devices that may have important applications in medicine as drug delivery, diagnostic or therapeutic devices (for example: smart endoscopic capsules), see e.g. [8–10]. The large size of the recent literature makes it impossible to provide an exhaustive survey of the state of the art. Specialized reviews have appeared in recent years, such as [6, 11–13], and several monographs are available [2, 14, 15]. The interested reader may find in these works and in the references cited therein several hundred papers to explore the subject.

At the scale of a single cell, viscous forces dominate inertia in fluid flows, which are then governed by the (steady) Stokes equations [16]. These arise from the Navier–Stokes equation, in the (formal) limit of zero Reynolds and Womersley numbers. The linearity of the Stokes system has the consequence that net propulsion through periodic shape actuation is only possible through histories of shape changes that are not reciprocal in time, a fact popularized as ‘the scallop theorem’ by Purcell [3].

Non-reciprocal periodic shape changes can arise in a number of ways: thanks to a non-trivial topology of the space of shapes, through non-trivial loops in the space of shapes, thanks to the propagation of traveling waves of
<figure><img src="image_2.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED

23 January 2019

ACCEPTED FOR PUBLICATION

20 March 2019

PUBLISHED

29 April 2019

Original content from this
work may be used under
the terms of the Creative
Commons Attribution 3.0
licence.

Any further distribution of
this work must maintain
attribution to the
author(s) and the title of
the work, journal citation
and DOI.

<figure><img src="image_3.png" /></figure>

PAPER

# Energy-optimal strokes for multi-link microswimmers: Purcell's loops and Taylor's waves reconciled

François Alouges¹, Antonio DeSimone² Ⓓ, Laetitia Giraldi³ Ⓓ, Yizhar Or⁴,⁵ Ⓓ and Oren Wiezel⁴ Ⓓ

¹ CMAP, Ecole Polytechnique, Cedex, France

² The BioRobotics Institute, Scuola Superiore Sant'Anna, Pisa, Italy

³ Université Côte d'Azur, Inria, CNRS, LJAD, France

⁴ Faculty of Mechanical EngineeringTechnion - Israel Institute of Technology, Israel

⁵ Author to whom any correspondence should be addressed

E-mail: izi@technion.ac.il

**Keywords:** microswimmer dynamics, optimal strokes, energy efficiency, asymptotic analysis

## Abstract

Micron-scale swimmers move in the realm of negligible inertia, dominated by viscous drag forces. In this paper, we formulate the leading-order dynamics of a slender multi-link (N-link) microswimmer assuming small-amplitude undulations about its straight configuration. The energy-optimal stroke to achieve a given prescribed displacement in a given time period is obtained as the largest eigenvalue solution of a constrained optimal control problem. Remarkably, the optimal stroke is an ellipse lying within a two-dimensional plane in the (N-1)-dimensional space of joint angles, where N can be arbitrarily large. For large N, the optimal stroke is a traveling wave of bending, modulo edge effects. If the number of shape variables is small, we can consider the same problem when the prescribed displacement in one time period is large, and not attainable with small variations of the joint angles. The fully nonlinear optimal control problem is solved numerically for the cases N = 3 (Purcell's three-link swimmer) and N = 5 showing that, as the prescribed displacement becomes small, the optimal solutions obtained using the small-amplitude assumption are recovered. We also show that, when the prescribed displacements become large, the picture is different. For N = 3 we recover the non-convex planar loops already known from previous studies. For N = 5 we obtain non-planar loops, raising the question of characterizing the geometry of complex high-dimensional loops.

## 1. Introduction

The analysis of biological and bio-inspired swimming at microscopic scales has attracted considerable attention in the recent literature, starting from the seminal work by Taylor [1], Lighthill [2], and Purcell [3]. One of the reasons is that swimming of unicellular organisms is at the root of many fundamental processes in biology: reproduction through the swimming of sperm cells is just one example [4–7]. Moreover, biology inspires the design of bio-mimetic artificial devices that may have important applications in medicine as drug delivery, diagnostic or therapeutic devices (for example: smart endoscopic capsules), see e.g. [8–10]. The large size of the recent literature makes it impossible to provide an exhaustive survey of the state of the art. Specialized reviews have appeared in recent years, such as [6, 11–13], and several monographs are available [2, 14, 15]. The interested reader may find in these works and in the references cited therein several hundred papers to explore the subject.

At the scale of a single cell, viscous forces dominate inertia in fluid flows, which are then governed by the (steady) Stokes equations [16]. These arise from the Navier–Stokes equation, in the (formal) limit of zero Reynolds and Womersley numbers. The linearity of the Stokes system has the consequence that net propulsion through periodic shape actuation is only possible through histories of shape changes that are not reciprocal in time, a fact popularized as ‘the scallop theorem’ by Purcell [3].

Non-reciprocal periodic shape changes can arise in a number of ways: thanks to a non-trivial topology of the space of shapes, through non-trivial loops in the space of shapes, thanks to the propagation of traveling waves of
deformation. Corresponding biological examples are, respectively, the rotary motion of helical bacterial flagellar bundles, the different shapes of cilia in the power and in the recovery part of one stroke, the beating of a eukaryotic flagellum causing the propagation of bending waves along the length of the flagellum. In fact, bending waves propagating along cilia/flagella and shape modulation during one stroke are used not only for propulsion by micro-organisms, but also for transport inside organs in humans and other higher organisms [17–19]. The second example is the one with more connections with the study of minimal artificial swimmers (i.e. swimmers with only two internal degrees of freedom to control shape such as the three-link swimmer of Purcell [3, 20], the three-sphere swimmer of Najafi and Golestanian [21], and others). The third example is possibly the most thoroughly exploited paradigm in the fabrication of micro-swimmer prototypes, often through the action of an external magnetic field on a flexible filament [8, 10, 22, 23].

Patterns of optimal actuation have been investigated independently, for each of these three examples, with a variety of analytical and numerical methods. For the case of flagellar and ciliary propulsion, these include [2, 24–28] leading to recognizing, for example, helical shapes as the optimal ones for filaments in three dimensions, and smoothed saw-tooth traveling waves as the optimal wave forms for the planar beating of a one-dimensional flagellum or for a planar sheet. In the limit of small amplitudes, the latter reduces to Taylor’s traveling sine waves [1]. For the study of optimal strokes for minimal artificial swimmers the reader is referred, e.g. to [20, 29–34].

More recently, it has been recognized that, in the presence of external forces or torques, the scallop paradigm has to be reconsidered. This is a consequence of the fact that, in this case, the governing equations can no longer be cast as an affine control systems without drift. Similar remarks apply to the presence of parts of the swimmer body whose shape is not directly controlled and it instead emerges from the balance between elastic restoring forces and viscous resistances. Since here we will not pursue this issue any further, we refer the reader to some of the relevant literature [35–40].

In spite of all the recent progress, several questions remain open for investigation as testified by the growing pace at which research articles on low Reynolds number swimmers are being published. New aspects of the very same fundamental swimming problem are thus continuously emerging. Rather than adding one more, we are motivated here by the quest for unifying perspectives over this vast literature. For example: what is the connection, if any, between optimal actuation by traveling waves inspired by Taylor’s swimming sheet [1] and optimal actuation by closed loops in the space of shapes shown in Tam and Hosoi’s gaits [29] for Purcell’s swimmer and further discussed in [41]? The two main paradigms for producing non-reciprocal shape changes have remained mostly disconnected in the recent literature, confined to two independent streams of research.

Motivated by this question, we focus on one specific example, namely, a planar swimmer consisting of *N* equal links of fixed length *L* representing the *N*-dimensional generalization of Purcell’s famous three-link swimmer, and use geometric control theory to identify some unifying principles. In particular, we consider the problem of prescribing an admissible displacement in one shape cycle, and look for the gait that minimizes a suitable cost functional giving a measure of the expended energy during that cycle.

The problem of optimal control is nonlinear in the shape parameters, and determining (even numerically) the optimal gait explicitly becomes very hard as soon as the number of shape parameters becomes large. Under the assumption of small-amplitude joint angles, by considering the approximation at leading order in the shape parameters, we obtain an affine control system which can be analyzed in full detail. In previous studies, asymptotic analysis in this small-angle regime has been used for Purcell’s and other swimmers where the number *N* of shape degrees of freedom is small [42–44], but never for large *N*. We find that optimal gaits are always two-dimensional elliptical loops, independent of *N*. These gaits bridge Purcell’s loops for the two-dimensional shape space associated with *N* = 3, to gaits that, modulo edge effects, can be identified with Taylor’s traveling waves of bending for large *N*. Interestingly, this analysis can be done without any explicit reference to a concrete model for the interaction between the swimmer and the surrounding fluid, which can be modeled with the full detail of Stokes hydrodynamics, or with any of the simplified models to treat the viscous drag for a slender swimmer. The result only depends on structural properties and symmetries of the governing equations, which in turn reflect the geometric symmetry of the physical problems. In fact, in this regime of small-amplitude joint angles, the perturbations from the rectilinear geometry of the reference configuration are small, and a slender one-dimensional swimmer with homogenous geometric and mechanical properties that interacts with a homogeneous surrounding medium is a system which is essentially invariant under shifts along the body axis. This is exactly true for an infinite or a periodic system and approximately true, modulo edge effects, for a system of finite length. The relevance of traveling waves as optimal gaits is therefore naturally suggested by the geometric symmetry of the system and, in fact, it emerges naturally from the symmetries of the governing equations.

Our result can be seen as a unifying principle in the sense that many different concrete problems, where different geometries of micro-swimmers obtained as *N* is varied, and different models for the fluid-structure interactions are employed, lead to different governing equations, but their behavior can be interpreted on the basis of the same single principle, thanks to the fact that the governing equations all share the same abstract
deformation. Corresponding biological examples are, respectively, the rotary motion of helical bacterial flagellar bundles, the different shapes of cilia in the power and in the recovery part of one stroke, the beating of a eukaryotic flagellum causing the propagation of bending waves along the length of the flagellum. In fact, bending waves propagating along cilia/flagella and shape modulation during one stroke are used not only for propulsion by micro-organisms, but also for transport inside organs in humans and other higher organisms [17–19]. The second example is the one with more connections with the study of minimal artificial swimmers (i.e. swimmers with only two internal degrees of freedom to control shape such as the three-link swimmer of Purcell [3, 20], the three-sphere swimmer of Najafi and Golestanian [21], and others). The third example is possibly the most thoroughly exploited paradigm in the fabrication of micro-swimmer prototypes, often through the action of an external magnetic field on a flexible filament [8, 10, 22, 23].

Patterns of optimal actuation have been investigated independently, for each of these three examples, with a variety of analytical and numerical methods. For the case of flagellar and ciliary propulsion, these include [2, 24–28] leading to recognizing, for example, helical shapes as the optimal ones for filaments in three dimensions, and smoothed saw-tooth traveling waves as the optimal wave forms for the planar beating of a one-dimensional flagellum or for a planar sheet. In the limit of small amplitudes, the latter reduces to Taylor’s traveling sine waves [1]. For the study of optimal strokes for minimal artificial swimmers the reader is referred, e.g. to [20, 29–34].

More recently, it has been recognized that, in the presence of external forces or torques, the scallop paradigm has to be reconsidered. This is a consequence of the fact that, in this case, the governing equations can no longer be cast as an affine control systems without drift. Similar remarks apply to the presence of parts of the swimmer body whose shape is not directly controlled and it instead emerges from the balance between elastic restoring forces and viscous resistances. Since here we will not pursue this issue any further, we refer the reader to some of the relevant literature [35–40].

In spite of all the recent progress, several questions remain open for investigation as testified by the growing pace at which research articles on low Reynolds number swimmers are being published. New aspects of the very same fundamental swimming problem are thus continuously emerging. Rather than adding one more, we are motivated here by the quest for unifying perspectives over this vast literature. For example: what is the connection, if any, between optimal actuation by traveling waves inspired by Taylor’s swimming sheet [1] and optimal actuation by closed loops in the space of shapes shown in Tam and Hosoi’s gaits [29] for Purcell’s swimmer and further discussed in [41]? The two main paradigms for producing non-reciprocal shape changes have remained mostly disconnected in the recent literature, confined to two independent streams of research.

Motivated by this question, we focus on one specific example, namely, a planar swimmer consisting of *N* equal links of fixed length *L* representing the *N*-dimensional generalization of Purcell’s famous three-link swimmer, and use geometric control theory to identify some unifying principles. In particular, we consider the problem of prescribing an admissible displacement in one shape cycle, and look for the gait that minimizes a suitable cost functional giving a measure of the expended energy during that cycle.

The problem of optimal control is nonlinear in the shape parameters, and determining (even numerically) the optimal gait explicitly becomes very hard as soon as the number of shape parameters becomes large. Under the assumption of small-amplitude joint angles, by considering the approximation at leading order in the shape parameters, we obtain an affine control system which can be analyzed in full detail. In previous studies, asymptotic analysis in this small-angle regime has been used for Purcell’s and other swimmers where the number *N* of shape degrees of freedom is small [42–44], but never for large *N*. We find that optimal gaits are always two-dimensional elliptical loops, independent of *N*. These gaits bridge Purcell’s loops for the two-dimensional shape space associated with *N* = 3, to gaits that, modulo edge effects, can be identified with Taylor’s traveling waves of bending for large *N*. Interestingly, this analysis can be done without any explicit reference to a concrete model for the interaction between the swimmer and the surrounding fluid, which can be modeled with the full detail of Stokes hydrodynamics, or with any of the simplified models to treat the viscous drag for a slender swimmer. The result only depends on structural properties and symmetries of the governing equations, which in turn reflect the geometric symmetry of the physical problems. In fact, in this regime of small-amplitude joint angles, the perturbations from the rectilinear geometry of the reference configuration are small, and a slender one-dimensional swimmer with homogenous geometric and mechanical properties that interacts with a homogeneous surrounding medium is a system which is essentially invariant under shifts along the body axis. This is exactly true for an infinite or a periodic system and approximately true, modulo edge effects, for a system of finite length. The relevance of traveling waves as optimal gaits is therefore naturally suggested by the geometric symmetry of the system and, in fact, it emerges naturally from the symmetries of the governing equations.

Our result can be seen as a unifying principle in the sense that many different concrete problems, where different geometries of micro-swimmers obtained as *N* is varied, and different models for the fluid-structure interactions are employed, lead to different governing equations, but their behavior can be interpreted on the basis of the same single principle, thanks to the fact that the governing equations all share the same abstract
<figure><img src="image_3.png" /><figcaption>**Figure 1.** The N-link swimmer. Represented here with $N = 5$ links. The position $(x, y)$ corresponds to the midpoint of the central link, and $\theta$ is the angle that this link makes with the horizontal axis. The shape parameters $(\phi_1, \dots, \phi_{N-1})$ are the angles between two neighboring links at each joint.</figcaption></figure>

structure and symmetries. To reinforce this point even further, we note that a similar result has been obtained in the different but related context of crawling of one-dimensional objects on solid surfaces or within a solid matrix [45], where it is shown that peristaltic waves used by earthworms for locomotion can be interpreted as optimal gaits when only small deformations are allowed. This is true provided, of course, that both mechanical properties of the crawler and frictional interactions with the surroundings are invariant with respect to shifts along the body axis.

To remove the assumption of small-amplitude joint angles and analyze the case in which large displacements are prescribed (and large joint angle variations are allowed) we can only resort to numerics, and our current capabilities limit the size of the system we can handle (in this paper, $N = 5$). This is nevertheless sufficient to reveal that optimal gaits of large amplitude have interesting and nontrivial geometry. They may differ both qualitatively and quantitatively from the planar elliptical loops representing their small amplitude limits. In particular, we find non-convex planar loops for $N = 3$, as it was already well-known from the analysis of optimal gaits for Purcell's swimmer, but also non-planar closed space curves for $N = 5$. However, all these nontrivial gaits duly converge to planar elliptical loops when the size of the prescribed displacement in one cycle becomes small.

Seen from the perspective of our results on the small-amplitude angle regime, these results are not surprising. In particular, when the restricted setting of small perturbations from the rectilinear geometry is abandoned, and large shape changes are considered, then invariance under shifts along the body axis is lost and traveling waves are no longer a natural basis for the study of the properties of optimal gaits. In addition, the optimal loops are no longer planar: the study of their geometric properties emerges in this way as a new and completely open field for future investigations that will require the development of new theoretical and numerical tools.

The rest of the paper is organized as follows. We describe the mathematical setting of the optimal control problem for the N-link swimmer in section 2. In section 3 we consider the regime of small-amplitude joint angles, hence of small deviations from the rectilinear shape, derive the general result that optimal loops are planar ellipses, and we solve numerically some specific examples which, for large N, show the emergence of traveling waves as optimal gaits. In section 4 we move to the regime of large amplitude joint angles, hence of large deviations from the rectilinear shape, and discuss some numerical results obtained via direct numerical optimization using the software BOCOP [46]: we exhibit non-planar complex loops for the case $N = 5$ and show that they converge to planar ellipses as the size of the prescribed displacement in one cycle becomes small.

## 2. Dynamics of the N-link swimmer

We develop our analysis at leading-order of the dynamics of the N-link neutrally buoyant swimmer using only the assumptions following from the linearity of Stokes equations, their invariance properties, and the symmetries of the swimmer. These assumptions are valid for several classic models such as Purcell's swimmer [3], the N-link swimmer in [47], and Taylor's swimming sheet [1]. They also hold irrespectively of whether the forces exerted by the fluid surrounding the swimmer are calculated by solving explicitly the Stokes equations, or they are evaluated by using any of the several approximate methods that have been used for slender swimmers (resistive force theory [4], slender body theory [48] etc).

The swimmer is planar and composed of *N* identical segments connected by rotational joints, see figure 1. We assume in the following that *N* is odd, and denote by (x, y) the position of the midpoint of the central link and *θ* its orientation (the angle that the central link makes with the horizontal axis), see figure 1. The case of even *N* simply requires a different choice of coordinates. The swimmer is immersed in an unbounded domain occupied by a homogeneous viscous fluid governed by Stokes equations and is driven by shape changes, i.e. the joint angles $\phi = (\phi_1, \phi_2, \dots, \phi_{N-1})^T$ are given functions of time t. In particular, we will be concerned with *T*-periodic
<figure><img src="image_4.png" /><figcaption>**Figure 1.** The N-link swimmer. Represented here with $N = 5$ links. The position $(x, y)$ corresponds to the midpoint of the central link, and $\theta$ is the angle that this link makes with the horizontal axis. The shape parameters $(\phi_1, \dots, \phi_{N-1})$ are the angles between two neighboring links at each joint.</figcaption></figure>

structure and symmetries. To reinforce this point even further, we note that a similar result has been obtained in the different but related context of crawling of one-dimensional objects on solid surfaces or within a solid matrix [45], where it is shown that peristaltic waves used by earthworms for locomotion can be interpreted as optimal gaits when only small deformations are allowed. This is true provided, of course, that both mechanical properties of the crawler and frictional interactions with the surroundings are invariant with respect to shifts along the body axis.

To remove the assumption of small-amplitude joint angles and analyze the case in which large displacements are prescribed (and large joint angle variations are allowed) we can only resort to numerics, and our current capabilities limit the size of the system we can handle (in this paper, $N = 5$). This is nevertheless sufficient to reveal that optimal gaits of large amplitude have interesting and nontrivial geometry. They may differ both qualitatively and quantitatively from the planar elliptical loops representing their small amplitude limits. In particular, we find non-convex planar loops for $N = 3$, as it was already well-known from the analysis of optimal gaits for Purcell's swimmer, but also non-planar closed space curves for $N = 5$. However, all these nontrivial gaits duly converge to planar elliptical loops when the size of the prescribed displacement in one cycle becomes small.

Seen from the perspective of our results on the small-amplitude angle regime, these results are not surprising. In particular, when the restricted setting of small perturbations from the rectilinear geometry is abandoned, and large shape changes are considered, then invariance under shifts along the body axis is lost and traveling waves are no longer a natural basis for the study of the properties of optimal gaits. In addition, the optimal loops are no longer planar: the study of their geometric properties emerges in this way as a new and completely open field for future investigations that will require the development of new theoretical and numerical tools.

The rest of the paper is organized as follows. We describe the mathematical setting of the optimal control problem for the N-link swimmer in section 2. In section 3 we consider the regime of small-amplitude joint angles, hence of small deviations from the rectilinear shape, derive the general result that optimal loops are planar ellipses, and we solve numerically some specific examples which, for large N, show the emergence of traveling waves as optimal gaits. In section 4 we move to the regime of large amplitude joint angles, hence of large deviations from the rectilinear shape, and discuss some numerical results obtained via direct numerical optimization using the software BOCOP [46]: we exhibit non-planar complex loops for the case $N = 5$ and show that they converge to planar ellipses as the size of the prescribed displacement in one cycle becomes small.

## 2. Dynamics of the N-link swimmer

We develop our analysis at leading-order of the dynamics of the N-link neutrally buoyant swimmer using only the assumptions following from the linearity of Stokes equations, their invariance properties, and the symmetries of the swimmer. These assumptions are valid for several classic models such as Purcell's swimmer [3], the N-link swimmer in [47], and Taylor's swimming sheet [1]. They also hold irrespectively of whether the forces exerted by the fluid surrounding the swimmer are calculated by solving explicitly the Stokes equations, or they are evaluated by using any of the several approximate methods that have been used for slender swimmers (resistive force theory [4], slender body theory [48] etc).

The swimmer is planar and composed of *N* identical segments connected by rotational joints, see figure 1. We assume in the following that *N* is odd, and denote by (x, y) the position of the midpoint of the central link and *θ* its orientation (the angle that the central link makes with the horizontal axis), see figure 1. The case of even *N* simply requires a different choice of coordinates. The swimmer is immersed in an unbounded domain occupied by a homogeneous viscous fluid governed by Stokes equations and is driven by shape changes, i.e. the joint angles $\phi = (\phi_1, \phi_2, \dots, \phi_{N-1})^T$ are given functions of time t. In particular, we will be concerned with *T*-periodic
<figure><img src="image_4.png" /><figcaption>Figure 2. Symmetry with respect to the x-axis. The swimmer at position $(x, y, \theta)$ and with shape $(\phi_1, \dots, \phi_{N-1})$ is symmetrized into a swimmer at position $(x, -y, -\theta)$ with shape $(-\phi_1, \dots, -\phi_{N-1})$.</figcaption></figure>

shape changes $\phi(t)$, which we call strokes, satisfying $\phi(t) = \phi(T + t)$, $\forall t > 0$. Our analysis aims at resolving the leading-order terms of the swimmer's dynamics, and is restricted to up to second order $O(\epsilon^2)$ terms in the shape parameters $\phi_i(t) = O(\epsilon)$, which are assumed to be small, i.e. $\epsilon \ll 1$.

The dynamics of the swimmer's planar motion is governed by

$$ \begin{cases} \dot{x} = f(\theta, \phi) \cdot \dot{\phi}, \\ \dot{y} = g(\theta, \phi) \cdot \dot{\phi}, \\ \dot{\theta} = h(\phi) \cdot \dot{\phi}. \end{cases} \qquad (1) $$

These equations come from the balance of viscous force and torque, which are linear in $\dot{x}, \dot{y}, \dot{\theta}$ and $\dot{\phi}$. Details of the derivation in the case of resistive force theory [4] are given in [43, 47, 49]. The special structure of (1), namely the fact that functions $f, g$, and $h$ are independent of $x$ and $y$, and that the last one is also independent of $\theta$, are consequences of the translational and rotational invariance of the problem.

Due to the geometric structure of the swimmer, the functions $f, g$ and $h$ that define the dynamics satisfy further relations that are deduced from symmetries of the system. Among them, the symmetry with respect to the x-axis, depicted in figure 2, transforms the swimmer parametrized by $\phi$ at position $(x, y, \theta)$ to the one parametrized by $-\phi$, at position $(x, -y, -\theta)$. The invariance of the dynamics under such a transformation leads to

$$ \begin{cases} f(\theta, \phi) \cdot \dot{\phi} = f(-\theta, -\phi) \cdot (-\dot{\phi}), \\ g(\theta, \phi) \cdot \dot{\phi} = -g(-\theta, -\phi) \cdot (-\dot{\phi}), \\ h(\phi) \cdot \dot{\phi} = -h(-\phi) \cdot (-\dot{\phi}). \end{cases} \qquad (2) $$

## 2.1. Dynamics at leading order

For $\theta$ and $\phi$ of order $\epsilon$, we expand the dynamics (1) to second order in $\theta$, $\phi$ and $\dot{\phi}$ by expanding $f$ to first order as

$$ f(\theta, \phi) \sim F_0 + F_0\theta + F_0\phi + O(\epsilon^2), \qquad (3) $$

where $F_0 = f(0, 0) \in \mathbb{R}^{N-1}$, $F_0 = \frac{\partial f}{\partial \theta}(0, 0) \in \mathbb{R}^{N-1}$ and $F_0 = \frac{\partial f}{\partial \phi}(0, 0) \in \mathbb{R}^{N-1,N-1}$. We also expand similarly $g$ and $h$. Using the symmetry relations (2), we get for $f, g$ and $h$ respectively

$$ \begin{aligned} & (F_0 + F_0\theta + F_0\phi) \cdot \dot{\phi} = -(F_0 - F_0\theta - F_0\phi) \cdot \dot{\phi} \\ & (G_0 + G_0\theta + G_0\phi) \cdot \dot{\phi} = (G_0 - G_0\theta - G_0\phi) \cdot \dot{\phi} \\ & (H_0 + H_0\phi) \cdot \dot{\phi} = (H_0 - H_0\phi) \cdot \dot{\phi}, \end{aligned} $$

from which we deduce, $\phi$ and $\dot{\phi}$ being arbitrary, that $F_0 = G_0 = 0$ and $G_0 = H_0 = 0$. The original system (1) therefore reduces to
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Symmetry with respect to the x-axis. The swimmer at position $(x, y, \theta)$ and with shape $(\phi_1, \dots, \phi_{N-1})$ is symmetrized into a swimmer at position $(x, -y, -\theta)$ with shape $(-\phi_1, \dots, -\phi_{N-1})$.</figcaption></figure>

shape changes $\phi(t)$, which we call strokes, satisfying $\phi(t) = \phi(T + t)$, $\forall t > 0$. Our analysis aims at resolving the leading-order terms of the swimmer's dynamics, and is restricted to up to second order $O(\varepsilon^2)$ terms in the shape parameters $\phi_i(t) = O(\varepsilon)$, which are assumed to be small, i.e. $\varepsilon \ll 1$.

The dynamics of the swimmer's planar motion is governed by

$$ \begin{cases} \dot{x} = f(\theta, \phi) \cdot \dot{\phi}, \\ \dot{y} = g(\theta, \phi) \cdot \dot{\phi}, \\ \dot{\theta} = h(\phi) \cdot \dot{\phi}. \end{cases} \qquad (1) $$

These equations come from the balance of viscous force and torque, which are linear in $\dot{x}, \dot{y}, \dot{\theta}$ and $\dot{\phi}$. Details of the derivation in the case of resistive force theory [4] are given in [43, 47, 49]. The special structure of (1), namely the fact that functions $f, g$, and $h$ are independent of $x$ and $y$, and that the last one is also independent of $\theta$, are consequences of the translational and rotational invariance of the problem.

Due to the geometric structure of the swimmer, the functions $f, g$ and $h$ that define the dynamics satisfy further relations that are deduced from symmetries of the system. Among them, the symmetry with respect to the x-axis, depicted in figure 2, transforms the swimmer parametrized by $\phi$ at position $(x, y, \theta)$ to the one parametrized by $-\phi$, at position $(x, -y, -\theta)$. The invariance of the dynamics under such a transformation leads to

$$ \begin{cases} f(\theta, \phi) \cdot \dot{\phi} = f(-\theta, -\phi) \cdot (-\dot{\phi}), \\ g(\theta, \phi) \cdot \dot{\phi} = -g(-\theta, -\phi) \cdot (-\dot{\phi}), \\ h(\phi) \cdot \dot{\phi} = -h(-\phi) \cdot (-\dot{\phi}). \end{cases} \qquad (2) $$

## 2.1. Dynamics at leading order

For $\theta$ and $\phi$ of order $\varepsilon$, we expand the dynamics (1) to second order in $\theta$, $\phi$ and $\dot{\phi}$ by expanding $f$ to first order as

$$ f(\theta, \phi) \sim F_0 + F_0\theta + F_0\phi + O(\varepsilon^2), \qquad (3) $$

where $F_0 = f(0, 0) \in \mathbb{R}^{N-1}$, $F_0 = \frac{\partial f}{\partial \theta}(0, 0) \in \mathbb{R}^{N-1}$ and $F_0 = \frac{\partial f}{\partial \phi}(0, 0) \in \mathbb{R}^{N-1, N-1}$. We also expand similarly $g$ and $h$. Using the symmetry relations (2), we get for $f, g$ and $h$ respectively

$$ \begin{aligned} & (F_0 + F_0\theta + F_0\phi) \cdot \dot{\phi} = -(F_0 - F_0\theta - F_0\phi) \cdot \dot{\phi} \\ & (G_0 + G_0\theta + G_0\phi) \cdot \dot{\phi} = (G_0 - G_0\theta - G_0\phi) \cdot \dot{\phi} \\ & (H_0 + H_0\phi) \cdot \dot{\phi} = (H_0 - H_0\phi) \cdot \dot{\phi}, \end{aligned} $$

from which we deduce, $\phi$ and $\dot{\phi}$ being arbitrary, that $F_0 = G_0 = 0$ and $G_0 = H_0 = 0$. The original system (1) therefore reduces to
$$ \begin{cases} \dot{x} = (\mathbf{F}_{\theta}\theta + \mathbf{F}_{\phi}\phi) \cdot \dot{\phi} + O(\epsilon^3), \\ \dot{y} = \mathbf{G}_0 \cdot \dot{\phi} + O(\epsilon^3), \\ \dot{\theta} = \mathbf{H}_0 \cdot \dot{\phi} + O(\epsilon^3). \end{cases} \qquad (4) $$

Integrating in time we obtain

$$ \theta(t) = \int_{0}^{t} \dot{\theta} \, dt = \int_{0}^{t} \mathbf{H}_{0} \cdot \dot{\phi} \, dt = \mathbf{H}_{0} \cdot (\phi(t) - \phi(0)) + O(\epsilon^3) \qquad (5) $$

and

$$ y(t) = \int_{0}^{t} \dot{y} \, dt = \int_{0}^{t} \mathbf{G}_{0} \cdot \dot{\phi} \, dt = \mathbf{G}_{0} \cdot (\phi(t) - \phi(0)) + O(\epsilon^3). $$

We deduce that, up to second order, the swimmer experiences no global rotation or transverse translation after one complete stroke. Indeed, $\phi$ being T-periodic, we have

$$ \Delta\theta = \theta(T) - \theta(0) = O(\epsilon^3), \qquad \Delta y = y(T) - y(0) = O(\epsilon^3). \qquad (6) $$

Moreover, using the expression for $\theta$ in (5) in the first equation of (4) and integrating over one period leads to

$$ \Delta x = x(T) - x(0) = \int_{0}^{T} (\mathbf{H}_{0}^{T} \otimes \mathbf{F}_{\theta} + \mathbf{F}_{\phi}) \phi \cdot \dot{\phi} \, dt + O(\epsilon^3). $$

Calling $\mathbf{F}_{\phi}^{\theta}$ the matrix $\mathbf{H}_{0}^{T} \otimes \mathbf{F}_{\theta} + \mathbf{F}_{\phi}$, integrating by parts the last expression, and using again the T-periodicity of $\phi$ leads to

$$ \int_{0}^{T} \mathbf{F}_{\phi}^{\theta} \phi \cdot \dot{\phi} \, dt = - \int_{0}^{T} \mathbf{F}_{\phi}^{\theta} \dot{\phi} \cdot \phi \, dt = - \int_{0}^{T} \mathbf{F}_{\phi}^{\theta T} \phi \cdot \dot{\phi} \, dt. $$

Finally we get the net longitudinal displacement in one cycle as

$$ \Delta x = \int_{0}^{T} \frac{1}{2} (\mathbf{F}_{\phi}^{\theta} - \mathbf{F}_{\phi}^{\theta T}) \phi \cdot \dot{\phi} \, dt + O(\epsilon^3). \qquad (7) $$

**Remark 2.1.** Our analysis at order two in the joint-angle amplitudes proves that the net lateral displacement $\Delta y$ and the net rotation both vanish at order two, as given in equation (6), while $\Delta x$ may be non-zero at order two, as shown in equation (7). The fundamental physical explanation of this result lies in the fact that the swimmer has symmetry about its straightened configurations, and only time-periodic strokes consisting of small-amplitude deviations about this configuration are analyzed. These are zero-mean trajectories, which, due to the swimmer's symmetry, cancel out all displacements at leading order, except for the $x$ displacement which is along the swimmer's axis of symmetry. This fact has also been previously explained in some related works [44, 50, 51]. This is also related to the fact the Lie-bracket vector field for each pair of the joint angle inputs, evaluated at the zero (straight) configuration, gives only $x$-displacement while all other motions only appear in higher order Lie brackets, see [52].

**Remark 2.2.** Other symmetries may be also used to obtain further information about the remaining coefficients of the system. For instance, using the rotational invariance of the system one can show that

$$ \begin{cases} f(\theta, \phi) = \cos(\theta)f(0, \phi) - \sin(\theta)g(0, \phi), \\ g(\theta, \phi) = \sin(\theta)f(0, \phi) + \cos(\theta)g(0, \phi), \end{cases} $$

from which one easily obtains $\mathbf{F}_{\theta} = -\mathbf{G}_{0}$ (and $\mathbf{G}_{\theta} = \mathbf{F}_{0} = 0$). Similarly, using the symmetry of the system with respect to the $y$-axis, it can be shown that both $\mathbf{F}_{\theta}$ and $\mathbf{G}_{0}$ are 'even' vectors (in the sense that $\mathbf{F}_{\theta,i} = \mathbf{F}_{\theta,N-i}$ and similarly for $\mathbf{G}_{0}$) while $\mathbf{H}_{0}$ is an 'odd' vector (for which $\mathbf{H}_{0,i} = -\mathbf{H}_{0,N-i}$).

## 2.2. Power expended

We now derive a leading-order expression for the mechanical energy expended by the swimmer during one cycle. The instantaneous power (i.e. rate of mechanical work) is given by the scalar product between force and velocity densities integrated over its surface (see e.g. [2]). Moreover, due to the linearity of Stokes equations, both the forces and the velocities acting on the swimmer depend linearly on $\dot{\phi}$. Thus, the instantaneous power density expended by the swimmer is a quadratic form in $\dot{\phi}$, with coefficients depending on $\phi$ (see [43] for concrete expressions using resistive force theory). We may therefore write the total energy $E$ expended by the swimmer during the stroke as
$$ \begin{cases} \dot{x} = (\mathbf{F}_{\theta}\theta + \mathbf{F}_{\phi}\phi) \cdot \dot{\phi} + O(\epsilon^3), \\ \dot{y} = \mathbf{G}_0 \cdot \dot{\phi} + O(\epsilon^3), \\ \dot{\theta} = \mathbf{H}_0 \cdot \dot{\phi} + O(\epsilon^3). \end{cases} \qquad (4) $$

Integrating in time we obtain

$$ \theta(t) = \int_{0}^{t} \dot{\theta} dt = \int_{0}^{t} \mathbf{H}_{0} \cdot \dot{\phi} dt = \mathbf{H}_{0} \cdot (\phi(t) - \phi(0)) + O(\epsilon^3) \qquad (5) $$

and

$$ y(t) = \int_{0}^{t} \dot{y} dt = \int_{0}^{t} \mathbf{G}_{0} \cdot \dot{\phi} dt = \mathbf{G}_{0} \cdot (\phi(t) - \phi(0)) + O(\epsilon^3). $$

We deduce that, up to second order, the swimmer experiences no global rotation or transverse translation after one complete stroke. Indeed, $\phi$ being T-periodic, we have

$$ \Delta\theta = \theta(T) - \theta(0) = O(\epsilon^3), \qquad \Delta y = y(T) - y(0) = O(\epsilon^3). \qquad (6) $$

Moreover, using the expression for $\theta$ in (5) in the first equation of (4) and integrating over one period leads to

$$ \Delta x = x(T) - x(0) = \int_{0}^{T} (\mathbf{H}_{0}^{T} \otimes \mathbf{F}_{\theta} + \mathbf{F}_{\phi}) \phi \cdot \dot{\phi} dt + O(\epsilon^3). $$

Calling $\mathbf{F}_{\phi}^{\theta}$ the matrix $\mathbf{H}_{0}^{T} \otimes \mathbf{F}_{\theta} + \mathbf{F}_{\phi}$, integrating by parts the last expression, and using again the T-periodicity of $\phi$ leads to

$$ \int_{0}^{T} \mathbf{F}_{\phi}^{\theta} \phi \cdot \dot{\phi} dt = - \int_{0}^{T} \mathbf{F}_{\phi}^{\theta} \dot{\phi} \cdot \phi dt = - \int_{0}^{T} \mathbf{F}_{\phi}^{\theta T} \phi \cdot \dot{\phi} dt. $$

Finally we get the net longitudinal displacement in one cycle as

$$ \Delta x = \int_{0}^{T} \frac{1}{2} (\mathbf{F}_{\phi}^{\theta} - \mathbf{F}_{\phi}^{\theta T}) \phi \cdot \dot{\phi} dt + O(\epsilon^3). \qquad (7) $$

**Remark 2.1.** Our analysis at order two in the joint-angle amplitudes proves that the net lateral displacement $\Delta y$ and the net rotation both vanish at order two, as given in equation (6), while $\Delta x$ may be non-zero at order two, as shown in equation (7). The fundamental physical explanation of this result lies in the fact that the swimmer has symmetry about its straightened configurations, and only time-periodic strokes consisting of small-amplitude deviations about this configuration are analyzed. These are zero-mean trajectories, which, due to the swimmer's symmetry, cancel out all displacements at leading order, except for the $x$ displacement which is along the swimmer's axis of symmetry. This fact has also been previously explained in some related works [44, 50, 51]. This is also related to the fact the Lie-bracket vector field for each pair of the joint angle inputs, evaluated at the zero (straight) configuration, gives only $x$-displacement while all other motions only appear in higher order Lie brackets, see [52].

**Remark 2.2.** Other symmetries may be also used to obtain further information about the remaining coefficients of the system. For instance, using the rotational invariance of the system one can show that

$$ \begin{cases} f(\theta, \phi) = \cos(\theta)f(0, \phi) - \sin(\theta)g(0, \phi), \\ g(\theta, \phi) = \sin(\theta)f(0, \phi) + \cos(\theta)g(0, \phi), \end{cases} $$

from which one easily obtains $\mathbf{F}_{\theta} = -\mathbf{G}_{0}$ (and $\mathbf{G}_{\theta} = \mathbf{F}_{0} = 0$). Similarly, using the symmetry of the system with respect to the $y$-axis, it can be shown that both $\mathbf{F}_{\theta}$ and $\mathbf{G}_{0}$ are 'even' vectors (in the sense that $\mathbf{F}_{\theta,i} = \mathbf{F}_{\theta,N-i}$ and similarly for $\mathbf{G}_{0}$) while $\mathbf{H}_{0}$ is an 'odd' vector (for which $\mathbf{H}_{0,i} = -\mathbf{H}_{0,N-i}$).

## 2.2. Power expended

We now derive a leading-order expression for the mechanical energy expended by the swimmer during one cycle. The instantaneous power (i.e. rate of mechanical work) is given by the scalar product between force and velocity densities integrated over its surface (see e.g. [2]). Moreover, due to the linearity of Stokes equations, both the forces and the velocities acting on the swimmer depend linearly on $\dot{\phi}$. Thus, the instantaneous power density expended by the swimmer is a quadratic form in $\dot{\phi}$, with coefficients depending on $\phi$ (see [43] for concrete expressions using resistive force theory). We may therefore write the total energy $E$ expended by the swimmer during the stroke as
$$E(\phi) = \int_{0}^{T} \langle \mathbf{P}(\phi) \dot{\phi}, \dot{\phi} \rangle dt, \qquad (8)$$

where $\mathbf{P}(\phi) \in \mathcal{M}_{N,N}(\mathbb{R})$ is a symmetric and positive definite matrix, and the bracket $\langle \cdot, \cdot \rangle$ stands for the scalar product in $\mathbb{R}^N$. Expanding this expression gives the $O(\epsilon^2)$ leading-order term of the energy as

$$E_0(\phi) = \int_0^T \langle \mathbf{P}_0 \dot{\phi}, \dot{\phi} \rangle dt, \qquad (9)$$

where $\mathbf{P}_0 := \mathbf{P}(\phi = 0)$ is symmetric and positive definite.

### 3. Optimal strokes of small amplitude

In this section, we study energy-optimal strokes while focusing on the case in which the joint angles remain small (small amplitude approximation), so that only small perturbations of the rectilinear geometry are allowed.

#### 3.1. Analysis of the optimal control problem

We begin by reviewing the well-known criterion of optimal energy efficiency due to Lighthill, and discuss its relation with the energy-optimal strokes studied here. Lighthill's efficiency of swimming is defined as the ratio of $E_{\text{drag}}/E_{\text{swim}}$, where $E_{\text{drag}}$ is the energy needed by an external actor to pull the swimmer during a time $T$ at an average speed $\Delta x/T$, while $E_{\text{swim}}$ is the energy expended by the swimmer during a stroke to propel itself at the same average speed. As the drag force is proportional to the velocity, $E_{\text{drag}}$ is proportional to $(\Delta x)^2/T$ whereas $E_{\text{swim}}$ behaves like $\Delta x$ as (7) and (8) indicate $(\Delta x$ and $E$ being quadratic in $(\phi, \dot{\phi})$ are of the same order $O(\epsilon^2))$. Thus, Lighthill's efficiency is proportional to $\Delta x$, and maximizing it amounts to increase $\Delta x$, and thus the amplitude of the stroke as much as possible, beyond the range of small angles. This implies that optimal strokes that maximize Lighthill's efficiency typically involve large-angle trajectories [29, 43]. On the other hand, if the amplitude of the stroke is constrained, by e.g. the maximum allowed amplitude of the angles, the displacement will be constrained as well, and it is expected that this constraint will be saturated when attempting to maximize Lighthill's efficiency, making the resulting optimal strokes of limited general interest. We therefore turn our attention to the problem of maximizing the efficiency (i.e. minimizing the energy) for a displacement achieved during one stroke which is bounded, or more simply fixed. Optimal strokes are thus defined as $T$-periodic strokes $\phi(t)$ that expend the minimal energy $E$, among all strokes that achieve a given displacement $\Delta x$ in a given time period $T$.

We showed in the previous section that, at leading order, the $y$ and $\theta$ displacements are negligible with respect to the $x$ displacement. Optimal strokes that provide a (longitudinal) displacement $\Delta x$ are thus sought as solutions to the constrained minimization problem

$$\min_{\phi \in \mathcal{D}} E_0(\phi), \qquad (10)$$

where, using (7), we have

$$\mathcal{D} := \left\{ \phi \in (\mathcal{C}^1(0, T))^{N-1} | \phi(0) = \phi(T) \text{ and } \frac{1}{2} \int_0^T (\mathbf{F}_{\phi}^{\theta} - \mathbf{F}_{\phi}^{\theta T}) \phi \cdot \dot{\phi} dt = \Delta x \right\}$$

In the following, we show that the solution of this optimization problem describes a planar ellipse in shape space, which is $(N-1)$-dimensional. We also give a method for its computation. Note that discrepancies between the solutions of the optimal control problem defined above under the leading-order approximation and under the exact nonlinear dynamics in (1) and (8) become smaller and smaller by decreasing the displacement $\Delta x$ and making the stroke amplitude $\epsilon$ of order $O(\sqrt{\Delta x})$.

In order to find the Euler–Lagrange first-order necessary condition associated with the constrained optimization problem (10), we set

$$\tilde{\mathcal{E}}_0(\phi) = \mathcal{E}_0(\phi) + \lambda K_0(\phi),$$

where $K_0(\phi) = \frac{1}{2} \int_0^T (\mathbf{F}_{\phi}^{\theta} - \mathbf{F}_{\phi}^{\theta T}) \phi \cdot \dot{\phi} dt$ and $\lambda \in \mathbb{R}$ is the Lagrange multiplier associated with the constraint $K_0(\phi) = \Delta x$. Writing that the first order functional derivative of $\tilde{\mathcal{E}}_0$ vanishes at the optimal stroke $\phi^*$

$$\delta \tilde{\mathcal{E}}_0(\phi^*) = 0$$

amounts to writing

$$\frac{d}{dt} \frac{\partial L}{\partial \dot{\phi}}(\phi^*, \dot{\phi}^*) - \frac{\partial L}{\partial \phi}(\phi^*, \dot{\phi}^*) = 0 \qquad (11)$$
$$E(\phi) = \int_{0}^{T} \langle \mathbf{P}(\phi) \dot{\phi}, \dot{\phi} \rangle dt, \qquad (8)$$

where $\mathbf{P}(\phi) \in \mathcal{M}_{N,N}(\mathbb{R})$ is a symmetric and positive definite matrix, and the bracket $\langle \cdot, \cdot \rangle$ stands for the scalar product in $\mathbb{R}^N$. Expanding this expression gives the $O(\epsilon^2)$ leading-order term of the energy as

$$E_0(\phi) = \int_0^T \langle \mathbf{P}_0 \dot{\phi}, \dot{\phi} \rangle dt, \qquad (9)$$

where $\mathbf{P}_0 := \mathbf{P}(\phi = 0)$ is symmetric and positive definite.

### 3. Optimal strokes of small amplitude

In this section, we study energy-optimal strokes while focusing on the case in which the joint angles remain small (small amplitude approximation), so that only small perturbations of the rectilinear geometry are allowed.

#### 3.1. Analysis of the optimal control problem

We begin by reviewing the well-known criterion of optimal energy efficiency due to Lighthill, and discuss its relation with the energy-optimal strokes studied here. Lighthill's efficiency of swimming is defined as the ratio of $E_{\text{drag}}/E_{\text{swim}}$, where $E_{\text{drag}}$ is the energy needed by an external actor to pull the swimmer during a time $T$ at an average speed $\Delta x/T$, while $E_{\text{swim}}$ is the energy expended by the swimmer during a stroke to propel itself at the same average speed. As the drag force is proportional to the velocity, $E_{\text{drag}}$ is proportional to $(\Delta x)^2/T$ whereas $E_{\text{swim}}$ behaves like $\Delta x$ as (7) and (8) indicate $(\Delta x$ and $E$ being quadratic in $(\phi, \dot{\phi})$ are of the same order $O(\epsilon^2))$. Thus, Lighthill's efficiency is proportional to $\Delta x$, and maximizing it amounts to increase $\Delta x$, and thus the amplitude of the stroke as much as possible, beyond the range of small angles. This implies that optimal strokes that maximize Lighthill's efficiency typically involve large-angle trajectories [29, 43]. On the other hand, if the amplitude of the stroke is constrained, by e.g. the maximum allowed amplitude of the angles, the displacement will be constrained as well, and it is expected that this constraint will be saturated when attempting to maximize Lighthill's efficiency, making the resulting optimal strokes of limited general interest. We therefore turn our attention to the problem of maximizing the efficiency (i.e. minimizing the energy) for a displacement achieved during one stroke which is bounded, or more simply fixed. Optimal strokes are thus defined as $T$-periodic strokes $\phi(t)$ that expend the minimal energy $E$, among all strokes that achieve a given displacement $\Delta x$ in a given time period $T$.

We showed in the previous section that, at leading order, the $y$ and $\theta$ displacements are negligible with respect to the $x$ displacement. Optimal strokes that provide a (longitudinal) displacement $\Delta x$ are thus sought as solutions to the constrained minimization problem

$$\min_{\phi \in \mathcal{D}} E_0(\phi), \qquad (10)$$

where, using (7), we have

$$\mathcal{D} := \left\{ \phi \in (C^1(0, T))^{N-1} \middle| \phi(0) = \phi(T) \text{ and } \frac{1}{2} \int_0^T (\mathbf{F}_{\phi}^{\theta} - \mathbf{F}_{\phi}^{\theta T}) \phi \cdot \dot{\phi} dt = \Delta x \right\}$$

In the following, we show that the solution of this optimization problem describes a planar ellipse in shape space, which is $(N-1)$-dimensional. We also give a method for its computation. Note that discrepancies between the solutions of the optimal control problem defined above under the leading-order approximation and under the exact nonlinear dynamics in (1) and (8) become smaller and smaller by decreasing the displacement $\Delta x$ and making the stroke amplitude $\epsilon$ of order $O(\sqrt{\Delta x})$.

In order to find the Euler-Lagrange first-order necessary condition associated with the constrained optimization problem (10), we set

$$\tilde{\mathcal{E}}_0(\phi) = \mathcal{E}_0(\phi) + \lambda K_0(\phi),$$

where $K_0(\phi) = \frac{1}{2} \int_0^T (\mathbf{F}_{\phi}^{\theta} - \mathbf{F}_{\phi}^{\theta T}) \phi \cdot \dot{\phi} dt$ and $\lambda \in \mathbb{R}$ is the Lagrange multiplier associated with the constraint $K_0(\phi) = \Delta x$. Writing that the first order functional derivative of $\tilde{\mathcal{E}}_0$ vanishes at the optimal stroke $\phi^*$

$$\delta \tilde{\mathcal{E}}_0(\phi^*) = 0$$

amounts to writing

$$\frac{d}{dt} \frac{\partial L}{\partial \dot{\phi}}(\phi^*, \dot{\phi}^*) - \frac{\partial L}{\partial \phi}(\phi^*, \dot{\phi}^*) = 0 \qquad (11)$$
where the Lagrangian L is such that

$$ \tilde{\mathcal{E}}_0(\phi) = \int_0^T L(\phi, \dot{\phi}) dt. $$

Here, we have

$$ L(\phi, \dot{\phi}) = \langle \mathbf{P}_0 \dot{\phi}, \dot{\phi} \rangle + \frac{\lambda}{2} (\mathbf{F}_{\phi}^{\theta} - \mathbf{F}_{\phi}^{\theta T}) \phi \cdot \dot{\phi} $$

and therefore (11) becomes

$$ \mathbf{P}_0 \ddot{\phi}^* = \frac{\lambda}{2} (\mathbf{F}_{\phi}^{\theta T} - \mathbf{F}_{\phi}^{\theta}) \dot{\phi}^*. \quad (12) $$

We denote by M the skew symmetric matrix

$$ \mathbf{M} = \frac{1}{2} \mathbf{P}_0^{-\frac{1}{2}} (\mathbf{F}_{\phi}^{\theta T} - \mathbf{F}_{\phi}^{\theta}) \mathbf{P}_0^{-\frac{1}{2}}, \quad (13) $$

and decompose the equation (12) along the eigen-elements of M. Eigenvectors of skew symmetric matrices go by pairs, associated with conjugate and purely imaginary eigenvalues. We therefore set $(\mathbf{v}_j^\pm)_{1 \le j \le N'} (N' = \lfloor \frac{N}{2} \rfloor)$ the (complex and orthonormal) eigenvectors associated with the purely imaginary eigenvalue $\pm i\mu_j$ with $\mu_j \ge 0$:

$$ \mathbf{M} \mathbf{v}_j^\pm = \pm i \mu_j \mathbf{v}_j^\pm. $$

Projecting $\mathbf{P}_0^1 \phi^*$ on the $(\mathbf{v}_j^\pm)_{1 \le j \le N'}$ as

$$ \mathbf{P}_0^{\frac{1}{2}} \phi^*(t) = \sum_{j=1}^{N'} \psi_j^{\pm}(t) \mathbf{v}_j^{\pm}, $$

we deduce from (12) that

$$ \ddot{\psi}_j^{\pm} = \pm i \lambda \mu_j \dot{\psi}_j^{\pm}, $$

or $\psi_j^{\pm}(t) = \frac{\alpha_j^{\pm}}{\pm i\lambda\mu_j} \exp(\pm i\lambda\mu_j t) + C_j^{\pm}$ where $C_j^{\pm}$ is a constant that we may take equal to 0. The solution of (12) is thus expressed as

$$ \phi^*(t) = \sum_{j=1}^{N'} \frac{\alpha_j^{\pm} e^{\pm i \lambda \mu_j t}}{\pm i \lambda \mu_j} \mathbf{P}_0^{-1/2} \mathbf{v}_j^{\pm}, \quad \text{for } t \in [0, T]. \quad (14) $$

Since we focus on periodic strokes, $\phi(0) = \phi(T)$, we must have

$$ \lambda \mu_j T = 2\pi k_j \quad 1 \le j \le N', \quad k_j \in \mathbb{N}. \quad (15) $$

By plugging the solution (14) into (9), we find

$$ E_0(\phi^*) = T \sum_{j=1}^{N'} (|\alpha_j^+|^2 + |\alpha_j^-|^2), \quad (16) $$

while the x-displacement is given by

$$ \Delta x = T \sum_{j=1}^{N'} \frac{(|\alpha_j^+|^2 + |\alpha_j^-|^2)}{\lambda} = \frac{E_0(\phi^*)}{\lambda}. \quad (17) $$

Since $\Delta x$ is fixed in the optimization problem, minimizing the energy requires choosing $\lambda$ as small as possible. In view of (15), this is achieved if $\lambda = \pm \frac{2\pi}{\mu_M T}$, where $\mu_M = \max\{\mu_j, 1 \le j \le N'\}$, the direction of the translation depending on the sign of $\lambda$, and $\psi_j^\pm = 0$ if $\mu_j \ne \mu_M$. Assuming $\lambda > 0$ ($\lambda < 0$ is handled similarly), we deduce that $\lambda = 2\pi/T\mu_M$ and the solution has only two modes corresponding to $\mathbf{v}_M^\pm$

$$ \phi^*(t) = \frac{\alpha_M^+ e^{2i\pi t}}{2i\pi} \mathbf{P}_0^{-1/2} \mathbf{v}_M^+ - \frac{\alpha_M^- e^{-2i\pi t}}{2i\pi} \mathbf{P}_0^{-1/2} \mathbf{v}_M^-, \quad \text{for } t \in [0, T], \quad (18) $$

and is therefore an ellipse drawn in the plane $(\mathbf{P}_0^{-1/2} \mathbf{v}_M^+, \mathbf{P}_0^{-1/2} \mathbf{v}_M^-)$. Furthermore, noticing that $\phi^*$ is real and $\mathbf{v}_M^+ = \overline{\mathbf{v}_M}$, we must have $\alpha_M^+ = \overline{\alpha_M}$. The net-displacement achieved by the gait (18) will be $\Delta x = \frac{2T}{\lambda} |\alpha_M^+|^2$. Therefore, the optimal gait amplitude scales as $\sqrt{\Delta x}$, as expected from (7).

Let us point out the important observation that the optimal gait (18) is always an ellipse lying within a two-dimensional plane regardless of the number of links N, which makes the dimension of the shape space arbitrarily large. In addition, the optimal gait in (18) can be written equivalently as sinusoidal inputs:
where the Lagrangian L is such that

$$ \tilde{\mathcal{E}}_0(\phi) = \int_0^T L(\phi, \dot{\phi}) dt. $$

Here, we have

$$ L(\phi, \dot{\phi}) = \langle \mathbf{P}_0 \dot{\phi}, \dot{\phi} \rangle + \frac{\lambda}{2} (\mathbf{F}_{\phi}^{\theta} - \mathbf{F}_{\phi}^{\theta T}) \phi \cdot \dot{\phi} $$

and therefore (11) becomes

$$ \mathbf{P}_0 \ddot{\phi}^* = \frac{\lambda}{2} (\mathbf{F}_{\phi}^{\theta T} - \mathbf{F}_{\phi}^{\theta}) \dot{\phi}^*. \quad (12) $$

We denote by **M** the skew symmetric matrix

$$ \mathbf{M} = \frac{1}{2} \mathbf{P}_0^{-\frac{1}{2}} (\mathbf{F}_{\phi}^{\theta T} - \mathbf{F}_{\phi}^{\theta}) \mathbf{P}_0^{-\frac{1}{2}}, \quad (13) $$

and decompose the equation (12) along the eigen-elements of **M**. Eigenvectors of skew symmetric matrices go by pairs, associated with conjugate and purely imaginary eigenvalues. We therefore set $(\mathbf{v}_j^\pm)_{1 \le j \le N'} (N' = \lfloor \frac{N}{2} \rfloor)$ the (complex and orthonormal) eigenvectors associated with the purely imaginary eigenvalue $\pm i\mu_j$ with $\mu_j \ge 0$:

$$ \mathbf{M} \mathbf{v}_j^\pm = \pm i \mu_j \mathbf{v}_j^\pm. $$

Projecting $\mathbf{P}_0^1 \phi^*$ on the $(\mathbf{v}_j^\pm)_{1 \le j \le N'}$ as

$$ \mathbf{P}_0^{\frac{1}{2}} \phi^*(t) = \sum_{j=1}^{N'} \psi_j^\pm(t) \mathbf{v}_j^\pm, $$

we deduce from (12) that

$$ \ddot{\psi}_j^\pm = \pm i \lambda \mu_j \dot{\psi}_j^\pm, $$

or $\psi_j^\pm(t) = \frac{\alpha_j^\pm}{\pm i \lambda \mu_j} \exp(\pm i \lambda \mu_j t) + C_j^\pm$ where $C_j^\pm$ is a constant that we may take equal to 0. The solution of (12) is thus expressed as

$$ \phi^*(t) = \sum_{j=1}^{N'} \frac{\alpha_j^\pm e^{\pm i \lambda \mu_j t}}{\pm i \lambda \mu_j} \mathbf{P}_0^{-1/2} \mathbf{v}_j^\pm, \quad \text{for } t \in [0, T]. \quad (14) $$

Since we focus on periodic strokes, $\phi(0) = \phi(T)$, we must have

$$ \lambda \mu_j T = 2\pi k_j \quad 1 \le j \le N', \quad k_j \in \mathbb{N}. \quad (15) $$

By plugging the solution (14) into (9), we find

$$ E_0(\phi^*) = T \sum_{j=1}^{N'} (|\alpha_j^+|^2 + |\alpha_j^-|^2), \quad (16) $$

while the x-displacement is given by

$$ \Delta x = T \sum_{j=1}^{N'} \frac{(|\alpha_j^+|^2 + |\alpha_j^-|^2)}{\lambda} = \frac{E_0(\phi^*)}{\lambda}. \quad (17) $$

Since $\Delta x$ is fixed in the optimization problem, minimizing the energy requires choosing $\lambda$ as small as possible. In view of (15), this is achieved if $\lambda = \pm \frac{2\pi}{\mu_M T}$, where $\mu_M = \max\{\mu_j, 1 \le j \le N'\}$, the direction of the translation depending on the sign of $\lambda$, and $\psi_j^\pm = 0$ if $\mu_j \ne \mu_M$. Assuming $\lambda > 0$ ($\lambda < 0$ is handled similarly), we deduce that $\lambda = 2\pi/T\mu_M$ and the solution has only two modes corresponding to $\mathbf{v}_M^\pm$

$$ \phi^*(t) = \frac{\alpha_M^+ e^{2i\pi t}}{2i\pi} \mathbf{P}_0^{-1/2} \mathbf{v}_M^+ - \frac{\alpha_M^- e^{-2i\pi t}}{2i\pi} \mathbf{P}_0^{-1/2} \mathbf{v}_M^-, \quad \text{for } t \in [0, T], \quad (18) $$

and is therefore an ellipse drawn in the plane $(\mathbf{P}_0^{-1/2} \mathbf{v}_M^+, \mathbf{P}_0^{-1/2} \mathbf{v}_M^-)$. Furthermore, noticing that $\phi^*$ is real and $\mathbf{v}_M^+ = \overline{\mathbf{v}_M}$, we must have $\alpha_M^+ = \overline{\alpha_M}$. The net-displacement achieved by the gait (18) will be $\Delta x = \frac{2T}{\lambda} |\alpha_M^+|^2$. Therefore, the optimal gait amplitude scales as $\sqrt{\Delta x}$, as expected from (7).

Let us point out the important observation that the optimal gait (18) is always an ellipse lying within a two-dimensional plane regardless of the number of links N, which makes the dimension of the shape space arbitrarily large. In addition, the optimal gait in (18) can be written equivalently as sinusoidal inputs:
$$ \phi_k^*(t) = a_k \sin(\omega t + p_k), \quad k = 1, \dots, N-1, \qquad (19) $$

where $\omega = 2\pi/T$. The amplitudes $a_k$ in (19) are of order $O(\epsilon)$, and scale as $\sqrt{\Delta x}$. Finally, an important property of the optimal solution $\phi^*(t)$ is that it satisfies

$$ P(t) = \langle \mathbf{P}_0 \dot{\phi}(t), \dot{\phi}(t) \rangle = 2|\alpha_M^+|^2 = \text{const}, \qquad (20) $$

where $P(t)$ is the mechanical power generated by the swimmer. Equation (20) implies that the optimal gait generates a constant power over the entire cycle. A similar result for the three-sphere swimmer was proved in [31] and this agrees with a fundamental observation made in [20], which states that for any given trajectory in shape space, the time-parametrization associated with constant power is the one that minimizes the total energy expenditure.

## 3.2. Numerical results: from Purcell's loops to Taylor's waves

In this section, we present numerical results for the problem of computing optimal gaits for slender swimmers by using the derivation described above and leading to the Euler-Lagrange equations (11). In order to derive a simple and concrete formulation of the dynamics, we use here the local drag approximation of resistive force theory [4, 48] for slender links. Our computations are made using MATLAB. A similar approach based on the Euler-Lagrange equations (11) was used in [31] for the three-sphere swimmer where, however, nonlocal hydrodynamic forces were fully resolved by solving the Stokes system for the surrounding fluid.

Our numerical code is based on the computation of $\phi^*(t)$ from formula (18). The first steps do not rely on the small amplitude approximation, and they will be used in later sections as well.

First, we derive the dynamics of the swimmer in (1) using resistive force theory [4, 48]. This theory states that the viscous drag force $f_i$ and torque $m_i$ on the ith slender link of length $l$ under planar motion are proportional to its linear and angular velocities, respectively. Thus, one can write the expression for the drag force and torque exerted on each link:

$$ \begin{aligned} f_i &= -c_i l (\mathbf{v}_i \cdot \mathbf{t}_i) \mathbf{t}_i - c_n l (\mathbf{v}_i \cdot \mathbf{n}_i) \mathbf{n}_i \\ m_i &= -\frac{1}{12} c_n l^3 \omega_i, \end{aligned} \qquad (21) $$

where $\mathbf{v}_i$, $\omega_i$ are the linear and angular velocities of the ith link, and $\mathbf{t}_i$, $\mathbf{n}_i$ are unit vectors in its axial and normal directions. The resistance coefficients in (21) are $c_n = 2c_t = 4\pi\eta/\log(l/a)$ where $\eta$ is the dynamic viscosity of the fluid and $a \ll l$ is the radius of the slender links' cross-section. Using (21), the swimmer's dynamic equations can be derived from force and torque balance, see for instance [43] for Purcell's swimmer and [47] for the general N-link swimmer.

Next, we exploit the assumption of small amplitude approximation and we calculate the matrices $F_\phi^\theta$ and $P_0$ in (7), (9) associated with the leading-order expressions for the dynamics and the mechanical energy expenditure. Then it only remains to compute the eigenvalues and eigenvectors of the matrix $M$ in (13) in order to assemble the expression of the optimal gait $\phi^*(t)$ in (18). This computation is done numerically in MATLAB.

## 3.3. Purcell's swimmer and N-link swimmers for small N

We use the small amplitude approximation and solve equation (12) numerically to compute optimal strokes for Purcell's three link swimmer and for other cases with small N. Purcell's three-link swimmer is the minimal model of a linked microswimmer. We used the formula (18) in order to compute the energy-optimal gait, which is shown in figure 6(a) as a trajectory in the plane of joint angles ($\phi_1, \phi_2$) for $\Delta x = 0.1l$. Snapshots of the swimmer's configuration at quarter-period times are also illustrated on the plot, and we recover the key ingredient that enables the propulsion of this system, namely a phase shift between the two actuators.

The case of N = 5 link swimmer is also computed and solutions are shown in figure 7. Four snapshots of the swimmer obtained again at quarter-period times are presented corresponding to a sequence of shapes along the optimal stroke.

Both cases (N = 3 and N = 5) will be discussed further and compared with those obtained without assuming that the prescribed displacements are small and joint angles have small amplitude in section 4.

## 3.4. Optimal gaits for swimmers with many links

We now show results of optimal gaits for swimmers with $N = 11$ and $N = 101$, obtained by solving the eigenvalue problem (12). Figure 3(a) shows a single snapshot of both swimmers under the optimal gait for a prescribed displacement $\Delta x = 0.1l$. The figure indicates that the optimal gaits look like a traveling wave. In order to test this observation quantitatively, we rewrite the optimal gaits using the sinusoidal representation (19), and compute the joint amplitudes $a_k$ and phase difference between consecutive joints
$$ \phi_k^*(t) = a_k \sin(\omega t + p_k), \quad k = 1, \dots, N-1, \qquad (19) $$

where $\omega = 2\pi/T$. The amplitudes $a_k$ in (19) are of order $O(\epsilon)$, and scale as $\sqrt{\Delta x}$. Finally, an important property of the optimal solution $\phi^*(t)$ is that it satisfies

$$ P(t) = \langle \mathbf{P}_0 \dot{\phi}(t), \dot{\phi}(t) \rangle = 2|\alpha_M^+|^2 = \text{const}, \qquad (20) $$

where $P(t)$ is the mechanical power generated by the swimmer. Equation (20) implies that the optimal gait generates a constant power over the entire cycle. A similar result for the three-sphere swimmer was proved in [31] and this agrees with a fundamental observation made in [20], which states that for any given trajectory in shape space, the time-parametrization associated with constant power is the one that minimizes the total energy expenditure.

## 3.2. Numerical results: from Purcell's loops to Taylor's waves

In this section, we present numerical results for the problem of computing optimal gaits for slender swimmers by using the derivation described above and leading to the Euler-Lagrange equations (11). In order to derive a simple and concrete formulation of the dynamics, we use here the local drag approximation of resistive force theory [4, 48] for slender links. Our computations are made using MATLAB. A similar approach based on the Euler-Lagrange equations (11) was used in [31] for the three-sphere swimmer where, however, nonlocal hydrodynamic forces were fully resolved by solving the Stokes system for the surrounding fluid.

Our numerical code is based on the computation of $\phi^*(t)$ from formula (18). The first steps do not rely on the small amplitude approximation, and they will be used in later sections as well.

First, we derive the dynamics of the swimmer in (1) using resistive force theory [4, 48]. This theory states that the viscous drag force $f_i$ and torque $m_i$ on the ith slender link of length $l$ under planar motion are proportional to its linear and angular velocities, respectively. Thus, one can write the expression for the drag force and torque exerted on each link:

$$ \begin{aligned} f_i &= -c_i l (\mathbf{v}_i \cdot \mathbf{t}_i) \mathbf{t}_i - c_n l (\mathbf{v}_i \cdot \mathbf{n}_i) \mathbf{n}_i \\ m_i &= -\frac{1}{12} c_n l^3 \omega_i, \end{aligned} \qquad (21) $$

where $\mathbf{v}_i$, $\omega_i$ are the linear and angular velocities of the ith link, and $\mathbf{t}_i$, $\mathbf{n}_i$ are unit vectors in its axial and normal directions. The resistance coefficients in (21) are $c_n = 2c_t = 4\pi\eta/\log(l/a)$ where $\eta$ is the dynamic viscosity of the fluid and $a \ll l$ is the radius of the slender links' cross-section. Using (21), the swimmer's dynamic equations can be derived from force and torque balance, see for instance [43] for Purcell's swimmer and [47] for the general N-link swimmer.

Next, we exploit the assumption of small amplitude approximation and we calculate the matrices $F_\phi^\theta$ and $P_0$ in (7), (9) associated with the leading-order expressions for the dynamics and the mechanical energy expenditure. Then it only remains to compute the eigenvalues and eigenvectors of the matrix $M$ in (13) in order to assemble the expression of the optimal gait $\phi^*(t)$ in (18). This computation is done numerically in MATLAB.

## 3.3. Purcell's swimmer and N-link swimmers for small N

We use the small amplitude approximation and solve equation (12) numerically to compute optimal strokes for Purcell's three link swimmer and for other cases with small N. Purcell's three-link swimmer is the minimal model of a linked microswimmer. We used the formula (18) in order to compute the energy-optimal gait, which is shown in figure 6(a) as a trajectory in the plane of joint angles $(\phi_1, \phi_2)$ for $\Delta x = 0.1l$. Snapshots of the swimmer's configuration at quarter-period times are also illustrated on the plot, and we recover the key ingredient that enables the propulsion of this system, namely a phase shift between the two actuators.

The case of $N = 5$ link swimmer is also computed and solutions are shown in figure 7. Four snapshots of the swimmer obtained again at quarter-period times are presented corresponding to a sequence of shapes along the optimal stroke.

Both cases ($N = 3$ and $N = 5$) will be discussed further and compared with those obtained without assuming that the prescribed displacements are small and joint angles have small amplitude in section 4.

## 3.4. Optimal gaits for swimmers with many links

We now show results of optimal gaits for swimmers with $N = 11$ and $N = 101$, obtained by solving the eigenvalue problem (12). Figure 3(a) shows a single snapshot of both swimmers under the optimal gait for a prescribed displacement $\Delta x = 0.1l$. The figure indicates that the optimal gaits look like a traveling wave. In order to test this observation quantitatively, we rewrite the optimal gaits using the sinusoidal representation (19), and compute the joint amplitudes $a_k$ and phase difference between consecutive joints
<figure><img src="image_5.png" /><figcaption>Figure 3. Optimal strokes for swimmers with $N = 11$ (red) and $N = 101$ (blue and blow-up) for $\Delta x = 0.01l$, where $l$ is the length of a single link. (a) A single snapshot of both swimmers during the optimal stroke. (b) Joint angle's amplitudes $a_k$ versus (scaled) joint number $k$. (c) Phase differences $\Delta p_k$, defined by (22), versus (scaled) joint number $k$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4. Optimal strokes for multi-link swimmers with large $N$ for $\Delta x = 0.01l$, where $l$ is the length of a single link. (a) Joint angle's amplitude $a_k$ averaged across all joints, as a function of $N$. (b) Phase difference $\Delta p_k$ averaged across all joints, as a function of $N$.</figcaption></figure>

$$ \Delta p_k = p_{k+1} - p_k \mod[360^\circ], \qquad (22) $$

where $p_k$ are given by (19). The results of amplitudes and phase differences across the joints for both swimmers are plotted in figures 3(b) and (c), respectively. Remarkably, the amplitudes $a_k$ for both swimmers display a nearly identical and slightly non-uniform symmetric distribution along the swimmer's body, with a peak at the center (and up to 14% decrease towards the ends), see figure 3(b). The phase difference between joints is very close to a uniform value of 120°, indicating a traveling wave with wavelength of three links. This is further confirmed when computing the optimal gait for increasing numbers of $N$ links for $\Delta x = 0.01l$. Figures 4(a) and (b) plot the gait's amplitude and phase difference, respectively, averaged across all swimmer's joints. Note that for a swimmer with discrete links, a waveform with length of integer number of $m$ links corresponds to a phase difference of $\Delta\varphi = 2\pi/m$. Taking $m = 1$ or $m = 2$ results in time-reversible motion, hence the smallest possible integer is $m = 3$. That is, the energy-optimal gait has the shortest possible wavelength. This result is in agreement with Taylor's observation in [1].

In hindsight, the emergence of traveling waves as optimal solutions for the control problem of a slender homogeneous swimmer whose shape is assumed to remain close to the rectilinear one, and swimming in a homogeneous fluid is not surprising. Consider for a moment the case in which the swimmer is of infinite length, but possesses a periodic shape (with period $Nl$). In particular, we look for strokes that are periodic in both space (with period $Nl$) and time (with period $T$). This system is translationally invariant with respect to shifts along the body axis. As a consequence, the initial problem (1), the optimal control problem (10), the governing operators, and the resulting solutions possess the same symmetry. In particular, this entails that the matrices $P_0$ and $F_\phi^\theta - F_\phi^\theta T$ (and therefore $M$) are circulant matrices, i.e. matrices $X$ that are constant along their diagonals and whose entries $X_{ij}$ only depend on $i-j \mod[N-1]$. Such matrices share the same eigenvectors that are given by $v_k^\pm = \frac{1}{\sqrt{N-1}}(1, \omega_k^\pm, (\omega_k^\pm)^2, \dots, (\omega_k^\pm)^{N-2})^T$ with $\omega_k^\pm = e^{\pm 2ik\pi/(N-1)}$ for $k = 1, \dots, (N-1)/2$. Using (18) shows
<figure><img src="image_6.png" /><figcaption>Figure 3. Optimal strokes for swimmers with $N = 11$ (red) and $N = 101$ (blue and blow-up) for $\Delta x = 0.01l$, where $l$ is the length of a single link. (a) A single snapshot of both swimmers during the optimal stroke. (b) Joint angle's amplitudes $a_k$ versus (scaled) joint number $k$. (c) Phase differences $\Delta p_k$, defined by (22), versus (scaled) joint number $k$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 4. Optimal strokes for multi-link swimmers with large $N$ for $\Delta x = 0.01l$, where $l$ is the length of a single link. (a) Joint angle's amplitude $a_k$ averaged across all joints, as a function of $N$. (b) Phase difference $\Delta p_k$ averaged across all joints, as a function of $N$.</figcaption></figure>

$$ \Delta p_k = p_{k+1} - p_k \mod[360^\circ], \qquad (22) $$

where $p_k$ are given by (19). The results of amplitudes and phase differences across the joints for both swimmers are plotted in figures 3(b) and (c), respectively. Remarkably, the amplitudes $a_k$ for both swimmers display a nearly identical and slightly non-uniform symmetric distribution along the swimmer's body, with a peak at the center (and up to 14% decrease towards the ends), see figure 3(b). The phase difference between joints is very close to a uniform value of 120°, indicating a traveling wave with wavelength of three links. This is further confirmed when computing the optimal gait for increasing numbers of $N$ links for $\Delta x = 0.01l$. Figures 4(a) and (b) plot the gait's amplitude and phase difference, respectively, averaged across all swimmer's joints. Note that for a swimmer with discrete links, a waveform with length of integer number of $m$ links corresponds to a phase difference of $\Delta\varphi = 2\pi/m$. Taking $m = 1$ or $m = 2$ results in time-reversible motion, hence the smallest possible integer is $m = 3$. That is, the energy-optimal gait has the shortest possible wavelength. This result is in agreement with Taylor's observation in [1].

In hindsight, the emergence of traveling waves as optimal solutions for the control problem of a slender homogeneous swimmer whose shape is assumed to remain close to the rectilinear one, and swimming in a homogeneous fluid is not surprising. Consider for a moment the case in which the swimmer is of infinite length, but possesses a periodic shape (with period $Nl$). In particular, we look for strokes that are periodic in both space (with period $Nl$) and time (with period $T$). This system is translationally invariant with respect to shifts along the body axis. As a consequence, the initial problem (1), the optimal control problem (10), the governing operators, and the resulting solutions possess the same symmetry. In particular, this entails that the matrices $P_0$ and $F_\phi^\theta - F_\phi^\theta T$ (and therefore $M$) are circulant matrices, i.e. matrices $X$ that are constant along their diagonals and whose entries $X_{ij}$ only depend on $i-j \mod[N-1]$. Such matrices share the same eigenvectors that are given by $v_k^\pm = \frac{1}{\sqrt{N-1}}(1, \omega_k^\pm, (\omega_k^\pm)^2, \dots, (\omega_k^\pm)^{N-2})^T$ with $\omega_k^\pm = e^{\pm 2ik\pi/(N-1)}$ for $k = 1, \dots, (N-1)/2$. Using (18) shows
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Log-log plots of the optimal energy $E^*$ (solid line) and scaled energy $Q$ (dashed line) for $\Delta x = 0.01L$ as a function of $N$. Here $L = Nl$ is the total length of the swimmer and the energy is normalized by $c_t L^3/T$.</figcaption></figure>

that the solution of the optimal control problem is therefore a traveling sinusoidal wave. In the finite length case, the same is true modulo edge-effects, as it is apparent from figure 3. Notice that a similar result has been obtained in the different but related context of crawling of one-dimensional objects on solid surfaces or within a solid matrix [45], where the optimality of actuation strategies based on peristaltic waves is discussed.

Finally, we fix the total length $L = Nl$ of a swimmer, compute the energy-optimal gait for moving a given displacement of $\Delta x = 0.01L$, and obtain the energy $E^*$ along this gait. Figure 5 shows a log-log plot of $E^*$ as a function of the links number $N$. Remarkably, it can be seen that for large $N$, the energy $E^*$ decays to zero as $1/N$. Analyzing this seemingly counter-intuitive behavior for large $N$ more closely (see appendix) reveals that the optimal energy indeed scales as $E^* \sim \frac{\Delta x l^2}{N}$ for large $N$. This suggests that a more suitable performance measure would be the scaled optimal energy defined as $Q = \frac{NE^*}{\Delta x l^2}$. This quantity, (multiplied by $\Delta x$ for better graphical visibility), is also plotted as a function of $N$ as the dashed line in figure 5, which indicates that it converges to a finite non-zero value at the limit of large $N$.

# 4. Optimal strokes of large amplitude

In this section we remove the assumption of small angle amplitudes, and consider the problem of determining energetically-optimal strokes when the prescribed displacement in one cycle is no longer small. At the present state of our knowledge, this can only be done numerically.

We compute numerically the optimal gaits by utilizing the BOCOP toolbox of optimal control, which uses direct optimization methods by discretizing times and states [46]. BOCOP applies numerical integration of the full nonlinear dynamic equations (1) and the energy formula (8), without assuming small-amplitude strokes. We solve the problem for large prescribed displacements. In addition, we consider the behavior of the optimal solutions when the prescribed displacements become smaller and smaller, in order to obtain an independent check of our analysis of section 3, which is based on leading-order approximation when imposed displacements and angle amplitudes are small. The discretized nonlinear optimization problem is solved by the IPOPT solver [53] with MUMPS [54], while the derivatives are computed by sparse automatic differentiation with CPPAD [55]. In the numerical experiments, we used a midpoint (implicit 1st order) discretization with 100 time steps.

## 4.1. Optimal gaits for Purcell's three-link swimmer

The optimal gaits computed within the framework of the small angle amplitudes of section 3.3 are shown in figure 6(a) as a trajectory in the plane of joint angles ($\phi_1, \phi_2$) for $\Delta x = 0.1l$. Snapshots of the swimmer's configuration at quarter-period intervals are also illustrated on the plot.

For comparison, we computed energy-optimal gaits using BOCOP for different displacements $\Delta x$, and the resulting trajectories scaled by $\sqrt{\Delta x}$ are shown in figure 6(b). It can be seen from the plot that for small displacements $\Delta x$ the optimal gait obtained by BOCOP converges (after scaling) to the one obtained by the formula (18). Note that we did not constrain the initial conditions of the swimmer in this computation with BOCOP, nor required any symmetry relations of the gait. The only constraint is on zero net rotation and net displacement of magnitude $\Delta x$, without a specified direction. The agreement between the two methods of computation confirms the validity of our small-amplitude analysis. When the displacement $\Delta x$ is increased further, the energy-optimal gaits obtained with BOCOP begin to deviate significantly from the small-amplitude one. For $\Delta x = 0.26$, the optimal gait (dashed) coincides with the gait obtained by Tam and Hosoi [29] that
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Log-log plots of the optimal energy $E^*$ (solid line) and scaled energy $Q$ (dashed line) for $\Delta x = 0.01L$ as a function of $N$. Here $L = NI$ is the total length of the swimmer and the energy is normalized by $c_t L^3/T$.</figcaption></figure>

that the solution of the optimal control problem is therefore a traveling sinusoidal wave. In the finite length case, the same is true modulo edge-effects, as it is apparent from figure 3. Notice that a similar result has been obtained in the different but related context of crawling of one-dimensional objects on solid surfaces or within a solid matrix [45], where the optimality of actuation strategies based on peristaltic waves is discussed.

Finally, we fix the total length $L = NI$ of a swimmer, compute the energy-optimal gait for moving a given displacement of $\Delta x = 0.01L$, and obtain the energy $E^*$ along this gait. Figure 5 shows a log-log plot of $E^*$ as a function of the links number $N$. Remarkably, it can be seen that for large $N$, the energy $E^*$ decays to zero as $1/N$. Analyzing this seemingly counter-intuitive behavior for large $N$ more closely (see appendix) reveals that the optimal energy indeed scales as $E^* \sim \frac{\Delta x l^2}{N}$ for large $N$. This suggests that a more suitable performance measure would be the scaled optimal energy defined as $Q = \frac{NE^*}{\Delta x l^2}$. This quantity, (multiplied by $\Delta x$ for better graphical visibility), is also plotted as a function of $N$ as the dashed line in figure 5, which indicates that it converges to a finite non-zero value at the limit of large $N$.

# 4. Optimal strokes of large amplitude

In this section we remove the assumption of small angle amplitudes, and consider the problem of determining energetically-optimal strokes when the prescribed displacement in one cycle is no longer small. At the present state of our knowledge, this can only be done numerically.

We compute numerically the optimal gaits by utilizing the BOCOP toolbox of optimal control, which uses direct optimization methods by discretizing times and states [46]. BOCOP applies numerical integration of the full nonlinear dynamic equations (1) and the energy formula (8), without assuming small-amplitude strokes. We solve the problem for large prescribed displacements. In addition, we consider the behavior of the optimal solutions when the prescribed displacements become smaller and smaller, in order to obtain an independent check of our analysis of section 3, which is based on leading-order approximation when imposed displacements and angle amplitudes are small. The discretized nonlinear optimization problem is solved by the IPOPT solver [53] with MUMPS [54], while the derivatives are computed by sparse automatic differentiation with CPPAD [55]. In the numerical experiments, we used a midpoint (implicit 1st order) discretization with 100 time steps.

## 4.1. Optimal gaits for Purcell's three-link swimmer

The optimal gaits computed within the framework of the small angle amplitudes of section 3.3 are shown in figure 6(a) as a trajectory in the plane of joint angles $(\phi_1, \phi_2)$ for $\Delta x = 0.1l$. Snapshots of the swimmer's configuration at quarter-period intervals are also illustrated on the plot.

For comparison, we computed energy-optimal gaits using BOCOP for different displacements $\Delta x$, and the resulting trajectories scaled by $\sqrt{\Delta x}$ are shown in figure 6(b). It can be seen from the plot that for small displacements $\Delta x$ the optimal gait obtained by BOCOP converges (after scaling) to the one obtained by the formula (18). Note that we did not constrain the initial conditions of the swimmer in this computation with BOCOP, nor required any symmetry relations of the gait. The only constraint is on zero net rotation and net displacement of magnitude $\Delta x$, without a specified direction. The agreement between the two methods of computation confirms the validity of our small-amplitude analysis. When the displacement $\Delta x$ is increased further, the energy-optimal gaits obtained with BOCOP begin to deviate significantly from the small-amplitude one. For $\Delta x = 0.26$, the optimal gait (dashed) coincides with the gait obtained by Tam and Hosoi [29] that
<figure><img src="image_8.png" /><figcaption>Figure 6. Optimal strokes for Purcell’s swimmer, plotted as loops in the plane of joint angles ($\phi_1, \phi_2$). (a) Optimal stroke obtained by our small-amplitude analysis, with snapshots of the swimmer’s configurations. (b) Comparison with optimal strokes computed using BOCOP for different displacements $\Delta x$, angles $\phi_i$ are scaled by $\sqrt{\Delta x}$.</figcaption></figure>

maximizes Lighthill’s efficiency. This was calculated by considering the first two terms in a Fourier series expansion of the two joint angles. Notably, the long axis of this ellipse-like-shaped gait is ‘skewed’ with respect to the long axis of the exact ellipse representing our small-amplitude energy-optimal gait. The reason for this fundamental difference is the fact that optimization of Lighthill’s efficiency as in [29] does not involve a constraint on the traveled distance. Therefore, it produces large displacements through large amplitude strokes, as discussed at the beginning of section 3.

When $\Delta x$ is further increased to an upper limit of $\Delta x = 0.306$, the optimal gait (dashed-dotted) deforms into the famous peanut-shaped loop obtained in [29] as the maximal-displacement gait. Note that for large displacements, we had to add constraints on symmetries of the gait and initial conditions, otherwise BOCOP began to search for different gaits which are not ‘simple loops’.

We now further analyze the motion of Purcell’s swimmer using leading-order terms as studied in [43]. The sinusoidal gait in (19) for the two joint angles can be rewritten as

$$ \phi_1(t) = \varepsilon \sin(t + \varphi/2) \quad \phi_2(t) = \varepsilon \sin(t - \varphi/2). \qquad (23) $$

That is, it has a stroke amplitude of $\varepsilon$ and phase difference of $\varphi$. Using resistive force theory and the leading order expansion in $\varepsilon$, as explained in [43], the leading-order expressions for displacement $\Delta x$ and energy $E$ in a cycle under the gait (23) are obtained as:

$$ \Delta x = \frac{5\pi}{81}\varepsilon^2 l \sin(\varphi), \quad E = \frac{2\pi c_t}{81}\varepsilon^2(11 \cos(\varphi) + 16). \qquad (24) $$

Therefore, for sinusoidal gaits of the form (23), our constrained optimization of minimizing the energy to cover a given displacement reduces to minimizing $E/\Delta x$. This gives an optimization problem for a scalar function of $\varphi$. Using elementary calculus, the minimum of this function is obtained at the optimal phase difference of $\varphi^* = \cos^{-1}(-11/16) = 133.43^{\circ}$, and the resulting optimal gait (23) is exactly identical to the one obtained by using the eigenvalue formulation in (18), which is shown in figure 6(a).

## 4.2. Optimal gaits for 5-link swimmer

We now show a comparison between our analytical formulation of optimal gaits in (18) and numerical optimization using BOCOP for the five-link swimmer, whose space of joint angles is four-dimensional. Snapshots at quarter-period times of the swimmer’s configuration for the analytical optimal gait corresponding to $\Delta x = 0.1l$ are shown in figure 7 and compared to the numerical solutions obtained for $\Delta x = l$. According to our small-amplitude analysis, the optimal gait is planar (it lies in the two-dimensional linear subspace $S$ spanned by eigenvectors associated with the pair of imaginary eigenvalues of $M$ corresponding to the maximal magnitude $\mu_M$). In order to compare further our analytical optimal gaits with the results of BOCOP computations, we first plot in figure 8(a) the projections onto the planar subspace $S$ of the four-dimensional optimal gaits. These optimal gaits were numerically computed by BOCOP for different values of $\Delta x$ and then scaled by $\sqrt{\Delta x}$ for comparison with the analytical optimal gaits. In order to test the theoretical prediction that optimal gaits should lie within the two-dimensional subspace $S$, we computed the maximal Euclidean distance $d$ of each optimal gait in $\mathbb{R}^4$ from the plane $S$. Figure 8(b) shows a log-log plot of this distance $d$ as a function of $\Delta x$.
<figure><img src="image_9.png" /><figcaption>Figure 6. Optimal strokes for Purcell’s swimmer, plotted as loops in the plane of joint angles ($\phi_1, \phi_2$). (a) Optimal stroke obtained by our small-amplitude analysis, with snapshots of the swimmer’s configurations. (b) Comparison with optimal strokes computed using BOCOP for different displacements $\Delta x$, angles $\phi_i$ are scaled by $\sqrt{\Delta x}$.</figcaption></figure>

maximizes Lighthill’s efficiency. This was calculated by considering the first two terms in a Fourier series expansion of the two joint angles. Notably, the long axis of this ellipse-like-shaped gait is ‘skewed’ with respect to the long axis of the exact ellipse representing our small-amplitude energy-optimal gait. The reason for this fundamental difference is the fact that optimization of Lighthill’s efficiency as in [29] does not involve a constraint on the traveled distance. Therefore, it produces large displacements through large amplitude strokes, as discussed at the beginning of section 3.

When $\Delta x$ is further increased to an upper limit of $\Delta x = 0.306$, the optimal gait (dashed-dotted) deforms into the famous peanut-shaped loop obtained in [29] as the maximal-displacement gait. Note that for large displacements, we had to add constraints on symmetries of the gait and initial conditions, otherwise BOCOP began to search for different gaits which are not ‘simple loops’.

We now further analyze the motion of Purcell’s swimmer using leading-order terms as studied in [43]. The sinusoidal gait in (19) for the two joint angles can be rewritten as

$$ \phi_1(t) = \varepsilon \sin(t + \varphi/2) \quad \phi_2(t) = \varepsilon \sin(t - \varphi/2). \qquad (23) $$

That is, it has a stroke amplitude of $\varepsilon$ and phase difference of $\varphi$. Using resistive force theory and the leading order expansion in $\varepsilon$, as explained in [43], the leading-order expressions for displacement $\Delta x$ and energy $E$ in a cycle under the gait (23) are obtained as:

$$ \Delta x = \frac{5\pi}{81}\varepsilon^2 l \sin(\varphi), \quad E = \frac{2\pi c_t}{81}\varepsilon^2(11 \cos(\varphi) + 16). \qquad (24) $$

Therefore, for sinusoidal gaits of the form (23), our constrained optimization of minimizing the energy to cover a given displacement reduces to minimizing $E/\Delta x$. This gives an optimization problem for a scalar function of $\varphi$. Using elementary calculus, the minimum of this function is obtained at the optimal phase difference of $\varphi^* = \cos^{-1}(-11/16) = 133.43^{\circ}$, and the resulting optimal gait (23) is exactly identical to the one obtained by using the eigenvalue formulation in (18), which is shown in figure 6(a).

## 4.2. Optimal gaits for 5-link swimmer

We now show a comparison between our analytical formulation of optimal gaits in (18) and numerical optimization using BOCOP for the five-link swimmer, whose space of joint angles is four-dimensional. Snapshots at quarter-period times of the swimmer’s configuration for the analytical optimal gait corresponding to $\Delta x = 0.1l$ are shown in figure 7 and compared to the numerical solutions obtained for $\Delta x = l$. According to our small-amplitude analysis, the optimal gait is planar (it lies in the two-dimensional linear subspace $S$ spanned by eigenvectors associated with the pair of imaginary eigenvalues of $M$ corresponding to the maximal magnitude $\mu_M$). In order to compare further our analytical optimal gaits with the results of BOCOP computations, we first plot in figure 8(a) the projections onto the planar subspace $S$ of the four-dimensional optimal gaits. These optimal gaits were numerically computed by BOCOP for different values of $\Delta x$ and then scaled by $\sqrt{\Delta x}$ for comparison with the analytical optimal gaits. In order to test the theoretical prediction that optimal gaits should lie within the two-dimensional subspace $S$, we computed the maximal Euclidean distance $d$ of each optimal gait in $\mathbb{R}^4$ from the plane $S$. Figure 8(b) shows a log-log plot of this distance $d$ as a function of $\Delta x$.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Snapshots of five-link swimmer during the optimal stroke. (a)–(d): analytical solutions for $\Delta x = 0.1l$; (e)–(h): numerical solutions for $\Delta x = l$ where $l$ is the length of a single link.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Comparison of the small-amplitude optimal stroke of the five-link swimmer with BOCOP computations for different displacements $\Delta x$. (a) Projection of the strokes onto the plane $S$, angles $\phi_i$ are scaled by $\sqrt{\Delta x}$ ($V_1, V_2 \in \mathbb{R}^4$ are two orthonormal basis vectors spanning the plane $S$). (b) Log-log plot of the maximal distance $d$ of optimal strokes obtained with BOCOP from the plane $S$, as a function of displacement $\Delta x$.</figcaption></figure>

It can be seen that for small displacements this distance decays to zero as $(\Delta x)^{3/2}$, indicating that the optimal gaits obtained numerically using BOCOP are indeed converging to planar loops lying within $S$, in agreement with the prediction of our asymptotic analysis.

Moving to the discussion of the optimal gaits of large amplitude, the fact that we are considering a system whose shape space is four-dimensional enables us to address and explore the deviations from planarity of the optimal loops. The computed optimal strokes are shown in the bottom panel of figure 7 and in figure 9, where we plot several different projections, by plotting the joint evolution of several pairs of shape variables. While, in the regime of small displacements, we do recover the elliptical loops predicted by our theory, the shape of the strokes for large displacements shows large discrepancies, and the geometry of the strokes is difficult to understand because of their high-dimensionality. One notable feature of figure 7 is that, in the large angle amplitude regime, optimal strokes consist of undulating shapes (wave-forms) such that the distance between successive peaks is of the order of the whole swimmer length, see figure 7(h), as it is in fact observed in biological systems such as sperm cells [6, 7].

In order to give a concrete and visual representation of the emergence of non-planarity when angles and displacements become large, we derive from the computed four-dimensional loops a dimensionally-reduced
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Snapshots of five-link swimmer during the optimal stroke. (a)–(d): analytical solutions for $\Delta x = 0.1l$; (e)–(h): numerical solutions for $\Delta x = l$ where $l$ is the length of a single link.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Comparison of the small-amplitude optimal stroke of the five-link swimmer with BOCOP computations for different displacements $\Delta x$. (a) Projection of the strokes onto the plane $S$, angles $\phi_i$ are scaled by $\sqrt{\Delta x}$ ($V_1, V_2 \in \mathbb{R}^4$ are two orthonormal basis vectors spanning the plane $S$). (b) Log-log plot of the maximal distance $d$ of optimal strokes obtained with BOCOP from the plane $S$, as a function of displacement $\Delta x$.</figcaption></figure>

It can be seen that for small displacements this distance decays to zero as $(\Delta x)^{3/2}$, indicating that the optimal gaits obtained numerically using BOCOP are indeed converging to planar loops lying within $S$, in agreement with the prediction of our asymptotic analysis.

Moving to the discussion of the optimal gaits of large amplitude, the fact that we are considering a system whose shape space is four-dimensional enables us to address and explore the deviations from planarity of the optimal loops. The computed optimal strokes are shown in the bottom panel of figure 7 and in figure 9, where we plot several different projections, by plotting the joint evolution of several pairs of shape variables. While, in the regime of small displacements, we do recover the elliptical loops predicted by our theory, the shape of the strokes for large displacements shows large discrepancies, and the geometry of the strokes is difficult to understand because of their high-dimensionality. One notable feature of figure 7 is that, in the large angle amplitude regime, optimal strokes consist of undulating shapes (wave-forms) such that the distance between successive peaks is of the order of the whole swimmer length, see figure 7(h), as it is in fact observed in biological systems such as sperm cells [6, 7].

In order to give a concrete and visual representation of the emergence of non-planarity when angles and displacements become large, we derive from the computed four-dimensional loops a dimensionally-reduced
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Optimal strokes of the five-link swimmer obtained with BOCOP computations for different displacements $\Delta x$. (a) $\phi_1$ versus $\phi_2$. (b) $\phi_2$ versus $\phi_3$. (c) $\phi_1$ versus $\phi_4$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Three-dimensional representation of the optimal loops for the 5-link-swimmer as a function of the amplitude of the prescribed displacement $\Delta x$, in normalized coordinates $(\pi_1, \pi_2, \pi_3)$. The 3D picture is produced using the ISOMAP approach and the optimal strokes for different prescribed displacements are plotted in black. It is therefore a three-dimensional nonlinear projection of the original angles in 4D. As the displacement goes to 0, the strokes converge to a planar ellipse, as the one obtained in the small displacement limit. In turn, when the prescribed displacement increases, the strokes are non-planar and their geometry is complex.</figcaption></figure>

representation in three dimensions, using the ISOMAP approach [56, 57]. Figure 10 shows clearly that the loops are non-planar, and that interesting features of their two-dimensional projections such as cusps and crossings are really an outcome of the projection of higher dimensional, non self-intersecting curves. This interesting behavior is entirely new with respect to the results currently available in the literature for the $N = 3$ case of Purcell's swimmer whose shape space is two-dimensional, and whose optimal strokes are necessarily planar curves. In addition, we see clearly that, as the imposed displacement becomes small, the loops converge to the planar ellipses of the small amplitude approximation of section 3.

These results open the new and unexplored question of characterizing the geometry of optimal loops for swimmers with shape space of large dimensions, when the amplitude of the joint angles is allowed to become large. Conducting numerical computation with BOCOP for a number of links $N > 5$ is currently beyond its limitations of memory allocation. There is however no such limitations for the solution of the small amplitude version of the optimal control problem (18) presented here.

**5. Conclusion**

In this work, we have studied optimal periodic strokes of multi-link micro-swimmers that achieve a given prescribed displacement in a given time with minimum energy. Exploiting the linearity of Stokes flows and geometric symmetries of the swimmer, the optimization has been formulated as a constrained variational problem, where leading-order expansion of the dynamics leads to the optimal solution of an eigenvalue problem. Remarkably, it is proven that energy-optimal strokes for N-link swimmers reduce to ellipses lying in a two-dimensional subspace of the space of shapes.

For large $N$, the optimal strokes become traveling waves with the shortest possible wavelength of three links, in agreement with the observation made in Taylor's classic work [1]. A noticeable difference from [1] due to the
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Optimal strokes of the five-link swimmer obtained with BOCOP computations for different displacements $\Delta x$. (a) $\phi_1$ versus $\phi_2$. (b) $\phi_2$ versus $\phi_3$. (c) $\phi_1$ versus $\phi_4$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Three-dimensional representation of the optimal loops for the 5-link-swimmer as a function of the amplitude of the prescribed displacement $\Delta x$, in normalized coordinates $(\pi_1, \pi_2, \pi_3)$. The 3D picture is produced using the ISOMAP approach and the optimal strokes for different prescribed displacements are plotted in black. It is therefore a three-dimensional nonlinear projection of the original angles in 4D. As the displacement goes to 0, the strokes converge to a planar ellipse, as the one obtained in the small displacement limit. In turn, when the prescribed displacement increases, the strokes are non-planar and their geometry is complex.</figcaption></figure>

representation in three dimensions, using the ISOMAP approach [56, 57]. Figure 10 shows clearly that the loops are non-planar, and that interesting features of their two-dimensional projections such as cusps and crossings are really an outcome of the projection of higher dimensional, non self-intersecting curves. This interesting behavior is entirely new with respect to the results currently available in the literature for the $N = 3$ case of Purcell's swimmer whose shape space is two-dimensional, and whose optimal strokes are necessarily planar curves. In addition, we see clearly that, as the imposed displacement becomes small, the loops converge to the planar ellipses of the small amplitude approximation of section 3.

These results open the new and unexplored question of characterizing the geometry of optimal loops for swimmers with shape space of large dimensions, when the amplitude of the joint angles is allowed to become large. Conducting numerical computation with BOCOP for a number of links $N > 5$ is currently beyond its limitations of memory allocation. There is however no such limitations for the solution of the small amplitude version of the optimal control problem (18) presented here.

## 5. Conclusion

In this work, we have studied optimal periodic strokes of multi-link micro-swimmers that achieve a given prescribed displacement in a given time with minimum energy. Exploiting the linearity of Stokes flows and geometric symmetries of the swimmer, the optimization has been formulated as a constrained variational problem, where leading-order expansion of the dynamics leads to the optimal solution of an eigenvalue problem. Remarkably, it is proven that energy-optimal strokes for N-link swimmers reduce to ellipses lying in a two-dimensional subspace of the space of shapes.

For large $N$, the optimal strokes become traveling waves with the shortest possible wavelength of three links, in agreement with the observation made in Taylor's classic work [1]. A noticeable difference from [1] due to the
<figure><img src="image_13.png" /></figure>

finite length of the swimmer is the slightly non-uniform distribution of energy-optimal wave amplitudes, which decay symmetrically from swimmer's center towards its ends.

Numerical optimization for the fully nonlinear problem, in which the amplitude of the prescribed displacements and the excursions of the joint angles are allowed to be large, is obtained using the optimization toolbox BOCOP for multi-link swimmers with three and five links. When the prescribed displacements are small, the numerical results show excellent agreement with those obtained using the small amplitude approximation. When the imposed displacements are large, the picture changes significantly. For Purcell's 3-link swimmer, we obtain non-convex closed curves similar to the ones previously obtained by Tam and Hosoi for maximizing displacement or Lighthill's efficiency. For the 5-link swimmer, we obtain non-planar loops of complex and intriguing geometry. A more precise characterization of the properties of these highly-dimensional closed loops seems an interesting open problem. For both $N = 3$ and $N = 5$, all these complex shapes converge to the planar ellipses of the small amplitude approximation when the prescribed displacements become small.

Possible directions for future extension of this research are comparison with measured strokes of biological swimming microorganisms as in [25, 58], generalization of the model to include elasticity as in [36–39], and the study of optimal control of magnetically-actuated microswimmers [10].

## Acknowledgments

OW wishes to thank the Technion-France Association for their support in his research internship visit to Paris in 2015, in which this work has been initiated. We gratefully acknowledge support by the Israel Science Foundation (OW and YO, grant no. 567/14), by the European Research Council (ADS, grant no. 340685-MicroMotility) and by the Centre National de la Recherche Scientifique (LG, CNRS Défi INFITI, projet COMM). We all wish to warmly thank Dr Pierre Martinon for his help and technical support with the BOCOP software.

## Appendix—scaling laws for large N

In order to study scaling of the displacement and energy for large $N$, we now consider a swimmer with $N$ links of fixed lengths $l = 1$, under the gait of traveling wave:

$$ \phi_k(t) = \varepsilon \sin(\omega t + k\Delta\varphi), \quad \text{for } k = 1 \dots N, \text{ where } \omega = 2\pi/T. \qquad (25) $$

We consider gaits with fixed amplitude $\varepsilon$ and phase difference of $\Delta\varphi = 2\pi/3$. Using our formulation based on resistive-force theory (21), we numerically calculate the displacement $\Delta x$ and energy expenditure $E$ along a cycle, as a function of $N$. The results are shown on log-log scale plots in figures A1(a) and (b). It can be seen that for large $N$, the displacement $\Delta x$ converges to a constant. This is analogous to the 'infinite sheet' limit of Taylor's net swimming speed [1]. In addition, $\Delta x$ is linear in $l$ and scales quadratically with the stroke amplitude $\varepsilon$ as in [1]. The energy expenditure $E$, however, grows linearly with $N$. This is because for large $N$, the interaction between links becomes negligible and each link contributes an equal amount of viscous dissipation. In addition, computation for fixed $N$ and varying links' length $l$ reveals that the energy expenditure scales cubically, as $l^3$ (plot
<figure><img src="image_14.png" /></figure>

finite length of the swimmer is the slightly non-uniform distribution of energy-optimal wave amplitudes, which decay symmetrically from swimmer's center towards its ends.

Numerical optimization for the fully nonlinear problem, in which the amplitude of the prescribed displacements and the excursions of the joint angles are allowed to be large, is obtained using the optimization toolbox BOCOP for multi-link swimmers with three and five links. When the prescribed displacements are small, the numerical results show excellent agreement with those obtained using the small amplitude approximation. When the imposed displacements are large, the picture changes significantly. For Purcell's 3-link swimmer, we obtain non-convex closed curves similar to the ones previously obtained by Tam and Hosoi for maximizing displacement or Lighthill's efficiency. For the 5-link swimmer, we obtain non-planar loops of complex and intriguing geometry. A more precise characterization of the properties of these highly-dimensional closed loops seems an interesting open problem. For both $N = 3$ and $N = 5$, all these complex shapes converge to the planar ellipses of the small amplitude approximation when the prescribed displacements become small.

Possible directions for future extension of this research are comparison with measured strokes of biological swimming microorganisms as in [25, 58], generalization of the model to include elasticity as in [36–39], and the study of optimal control of magnetically-actuated microswimmers [10].

## Acknowledgments

OW wishes to thank the Technion-France Association for their support in his research internship visit to Paris in 2015, in which this work has been initiated. We gratefully acknowledge support by the Israel Science Foundation (OW and YO, grant no. 567/14), by the European Research Council (ADS, grant no. 340685-MicroMotility) and by the Centre National de la Recherche Scientifique (LG, CNRS Défi INFITI, projet COMM). We all wish to warmly thank Dr Pierre Martinon for his help and technical support with the BOCOP software.

## Appendix—scaling laws for large N

In order to study scaling of the displacement and energy for large $N$, we now consider a swimmer with $N$ links of fixed lengths $l = 1$, under the gait of traveling wave:

$$ \phi_k(t) = \varepsilon \sin(\omega t + k\Delta\varphi), \quad \text{for } k = 1 \dots N, \text{ where } \omega = 2\pi/T. \qquad (25) $$

We consider gaits with fixed amplitude $\varepsilon$ and phase difference of $\Delta\varphi = 2\pi/3$. Using our formulation based on resistive-force theory (21), we numerically calculate the displacement $\Delta x$ and energy expenditure $E$ along a cycle, as a function of $N$. The results are shown on log-log scale plots in figures A1(a) and (b). It can be seen that for large $N$, the displacement $\Delta x$ converges to a constant. This is analogous to the 'infinite sheet' limit of Taylor's net swimming speed [1]. In addition, $\Delta x$ is linear in $l$ and scales quadratically with the stroke amplitude $\varepsilon$ as in [1]. The energy expenditure $E$, however, grows linearly with $N$. This is because for large $N$, the interaction between links becomes negligible and each link contributes an equal amount of viscous dissipation. In addition, computation for fixed $N$ and varying links' length $l$ reveals that the energy expenditure scales cubically, as $l^3$ (plot
not shown). This is explained by the following observation. The mechanical power dissipation of each link scales as $P_i = f_i v_i$ where $f_i$ is the viscous drag force and $v_i$ is the link's linear velocity. The viscous force $f_i$ scales as $f_i = Rv_i$ where R is a viscous resistance coefficient (see [4, 16, 48]). The link's velocity scales as $v_i \sim \varepsilon l/T$, while the resistance coefficient R scales linearly with the link's length l. Since total energy dissipation scales as $NP_i T$, we deduce the following scaling laws for large N:

$$ \Delta x \sim \varepsilon^2 l, \quad E \sim \varepsilon^2 Nl^3/T. \qquad (26) $$

Consider now the case where the required displacement $\Delta x$ as well as the total length of the swimmer $L = Nl$, are held fixed while N is varied. That is, the links' length scales as $l = L/N$. The scaling laws in (26) will now change to

$$ \Delta x \sim \varepsilon^2 \frac{L}{N}, \quad E \sim \varepsilon^2 \frac{L^3}{N^2 T}. \qquad (27) $$

Fixing the displacement $\Delta x$, it is deduced from (26) that for large N, the amplitude $\varepsilon$ increases as $\sqrt{\frac{N\Delta x}{L}}$. Moreover, the energy scales as $E \sim \frac{\Delta x l^2}{N T}$. This scaling relation explains the decay rate of $E^*$ in figure 5.

## ORCID iDs

Antonio DeSimone https://orcid.org/0000-0002-2632-3057

Laetitia Giraldi https://orcid.org/0000-0003-2684-0203

Yizhar Or https://orcid.org/0000-0002-9091-9357

Oren Wiezel https://orcid.org/0000-0003-0097-959X

## References

[1] Taylor G I 1951 Analysis of the swimming of microscopic organisms Proc. R. Soc. A 209 447–61

[2] Lighthill J 1975 Mathematical biofluiddynamic Society for Industrial and Applied Mathematics () (Philadelphia, PA: Saunders)

[3] Purcell E M 1977 Life at low Reynolds number Am. J. Phys. 45 3–11

[4] Gray J and Hancock G J 1955 The propulsion of sea-urchin spermatozoa J. Exp. Biol. 32 802–14 http://jeb.biologists.org/content/32/4/802

[5] Dresdner R D and Katz D F 1981 Relationships of mammalian sperm motility and morphology to hydrodynamic aspects of cell function Biol. Reprod. 25 920–30

[6] Gaffney E, Gadélha H, Smith D, Blake J and Kirkman-Brown J 2011 Mammalian sperm motility: observation and theory Annu. Rev. Fluid Mech. 43 501–28

[7] Ishimoto K, Gadélha H, Gaffney E, Smith D and Kirkman-Brown J 2017 Coarse-graining the fluid flow around a human sperm Phys. Rev. Lett. 118 124501

[8] Dreyfus R, Baudry J, Roper M L, Fermigier M, Stone H A and Bibette J 2005 Microscopic artificial swimmers Nature 437 862–5

[9] Feng J and Cho S 2014 Mini and micro propulsion for medical swimmers Micromachines 5 97–113

[10] Alouges F, DeSimone A, Giraldi L and Zoppello M 2015 Can magnetic multilayers propel artificial microswimmers mimicking sperm cells? Soft Robot. 2 117–28

[11] Lauga E and Powers T R 2009 The hydrodynamics of swimming microorganisms Rep. Prog. Phys. 72 096601

[12] Guasto J, Rusconi R and Stoker R 2012 Fluid mechanics of planktonic microorganisms Annu. Rev. Fluid Mech. 44 373–400

[13] Goldstein R E 2015 Green algae as model organisms for biological fluid dynamics Annu. Rev. Fluid Mech. 47 343–75

[14] Childress S 1981 Mechanics of Swimming and Flying (Cambridge: Cambridge University Press)

[15] Bonnard B, Chyba M and Rouot J 2018 Geometric and Numerical Optimal Control: Application to Swimming at Low Reynolds Number and Magnetic Resonance Imaging (Berlin: Springer)

[16] Brenner H and Happel J 1965 Low Reynolds Number Hydrodynamics: With Special Applications to Particulate Media (Berlin: Springer)

[17] Blake J R 1971 A spherical envelope approach to ciliary propulsion J. Fluid Mech. 46 199–208

[18] Satir P and Sleigh M A 1990 The physiology of cilia and muco-ciliary interactions Annu. Rev. Physiol. 52 137–55

[19] Feriani J, Juenet M, Fowler C, Bruot N, Chioccioli M, Holland S, Bryant C and Cicuta P 2017 Assessing the collective dynamics of motile cilia in cultures of human airway cells Biophys. J. 113 109–19

[20] Becker L E, Koehler S A and Stone H A 2003 On self-propulsion of micro-machines at low Reynolds number: Purcell's three-link swimmer J. Fluid Mech. 490 15–35

[21] Najafi A and Golestanian R 2004 Simple swimmer at low Reynolds number: three linked spheres Phys. Rev. E 69 062901

[22] Ghosch H and Fischer P 2009 Controlled propulsion of artificial magnetic nanostructured propellers Nano Lett. 9 2243–5

[23] Hu W, Lum G, Mastrangeli M and Sitti M 2018 Small-scale soft-bodied robot with multimodal locomotion Nature 554 81–5

[24] Pironneau O and Katz D F 1974 Optimal swimming of flagellated micro-organisms J. Fluid Mech. 66 391–415

[25] Berman R S, Kenneth O, Sznitman J and Leshansky A M 2013 Undulatory locomotion of finite filaments: lessons from Caenorhabditis elegans New J. Phys. 15 075022

[26] Montenegro-Johnson F and Lauga E 2014 Optimal swimming of a sheet Phys. Rev. E 89 060701(R)

[27] Tam D and Hosoi A E 2011 Optimal feeding and swimming gaits of biflagellated organisms Proc. Natl Acad. Sci. USA 108 1001–6

[28] Osterman N and Vilfan A 2011 Finding the ciliary beating pattern with optimal efficiency Proc. Natl. Acad. Sci. USA 108 15727–32

[29] Tam D and Hosoi A E 2007 Optimal stroke patterns for Purcell's three-link swimmer Phys. Rev. Lett. 98 068105

[30] Alouges F, DeSimone A, Heltai L, Lefebvre A and Merlet B 2013 Optimally swimming Stokesian robots Discrete Continuous Dyn. Syst. B 18 1189–215

[31] Alouges F, DeSimone A and Lefebvre A 2008 Optimal strokes for low Reynolds number swimmers: an example J. Nonlinear Sci. 18 277–302
not shown). This is explained by the following observation. The mechanical power dissipation of each link scales as $P_i = f_i v_i$ where $f_i$ is the viscous drag force and $v_i$ is the link's linear velocity. The viscous force $f_i$ scales as $f_i = R v_i$ where R is a viscous resistance coefficient (see [4, 16, 48]). The link's velocity scales as $v_i \sim \varepsilon l/T$, while the resistance coefficient R scales linearly with the link's length l. Since total energy dissipation scales as $NP_i T$, we deduce the following scaling laws for large N:

$$ \Delta x \sim \varepsilon^2 l, \quad E \sim \varepsilon^2 Nl^3/T. \qquad (26) $$

Consider now the case where the required displacement $\Delta x$ as well as the total length of the swimmer $L = Nl$, are held fixed while N is varied. That is, the links' length scales as $l = L/N$. The scaling laws in (26) will now change to

$$ \Delta x \sim \varepsilon^2 \frac{L}{N}, \quad E \sim \varepsilon^2 \frac{L^3}{N^2 T}. \qquad (27) $$

Fixing the displacement $\Delta x$, it is deduced from (26) that for large N, the amplitude $\varepsilon$ increases as $\sqrt{\frac{N \Delta x}{L}}$. Moreover, the energy scales as $E \sim \frac{\Delta x l^2}{N T}$. This scaling relation explains the decay rate of $E^*$ in figure 5.

## ORCID iDs

Antonio DeSimone https://orcid.org/0000-0002-2632-3057

Laetitia Giraldi https://orcid.org/0000-0003-2684-0203

Yizhar Or https://orcid.org/0000-0002-9091-9357

Oren Wiezel https://orcid.org/0000-0003-0097-959X

## References

[1] Taylor G I 1951 Analysis of the swimming of microscopic organisms Proc. R. Soc. A 209 447–61

[2] Lighthill J 1975 Mathematical biofluiddynamic Society for Industrial and Applied Mathematics () (Philadelphia, PA: Saunders)

[3] Purcell E M 1977 Life at low Reynolds number Am. J. Phys. 45 3–11

[4] Gray J and Hancock G J 1955 The propulsion of sea-urchin spermatozoa J. Exp. Biol. 32 802–14 http://jeb.biologists.org/content/32/4/802

[5] Dresdner R D and Katz D F 1981 Relationships of mammalian sperm motility and morphology to hydrodynamic aspects of cell function Biol. Reprod. 25 920–30

[6] Gaffney E, Gadélha H, Smith D, Blake J and Kirkman-Brown J 2011 Mammalian sperm motility: observation and theory Annu. Rev. Fluid Mech. 43 501–28

[7] Ishimoto K, Gadélha H, Gaffney E, Smith D and Kirkman-Brown J 2017 Coarse-graining the fluid flow around a human sperm Phys. Rev. Lett. 118 124501

[8] Dreyfus R, Baudry J, Roper M L, Fermigier M, Stone H A and Bibette J 2005 Microscopic artificial swimmers Nature 437 862–5

[9] Feng J and Cho S 2014 Mini and micro propulsion for medical swimmers Micromachines 5 97–113

[10] Alouges F, DeSimone A, Giraldi L and Zoppello M 2015 Can magnetic multilayers propel artificial microswimmers mimicking sperm cells? Soft Robot. 2 117–28

[11] Lauga E and Powers T R 2009 The hydrodynamics of swimming microorganisms Rep. Prog. Phys. 72 096601

[12] Guasto J, Rusconi R and Stoker R 2012 Fluid mechanics of planktonic microorganisms Annu. Rev. Fluid Mech. 44 373–400

[13] Goldstein R E 2015 Green algae as model organisms for biological fluid dynamics Annu. Rev. Fluid Mech. 47 343–75

[14] Childress S 1981 Mechanics of Swimming and Flying (Cambridge: Cambridge University Press)

[15] Bonnard B, Chyba M and Rouot J 2018 Geometric and Numerical Optimal Control: Application to Swimming at Low Reynolds Number and Magnetic Resonance Imaging (Berlin: Springer)

[16] Brenner H and Happel J 1965 Low Reynolds Number Hydrodynamics: With Special Applications to Particulate Media (Berlin: Springer)

[17] Blake J R 1971 A spherical envelope approach to ciliary propulsion J. Fluid Mech. 46 199–208

[18] Satir P and Sleigh M A 1990 The physiology of cilia and muco-ciliary interactions Annu. Rev. Physiol. 52 137–55

[19] Feriani J, Juenet M, Fowler C, Bruot N, Chioccioli M, Holland S, Bryant C and Cicuta P 2017 Assessing the collective dynamics of motile cilia in cultures of human airway cells Biophys. J. 113 109–19

[20] Becker L E, Koehler S A and Stone H A 2003 On self-propulsion of micro-machines at low Reynolds number: Purcell's three-link swimmer J. Fluid Mech. 490 15–35

[21] Najafi A and Golestanian R 2004 Simple swimmer at low Reynolds number: three linked spheres Phys. Rev. E 69 062901

[22] Ghosch H and Fischer P 2009 Controlled propulsion of artificial magnetic nanostructured propellers Nano Lett. 9 2243–5

[23] Hu W, Lum G, Mastrangeli M and Sitti M 2018 Small-scale soft-bodied robot with multimodal locomotion Nature 554 81–5

[24] Pironneau O and Katz D F 1974 Optimal swimming of flagellated micro-organisms J. Fluid Mech. 66 391–415

[25] Berman R S, Kenneth O, Sznitman J and Leshansky A M 2013 Undulatory locomotion of finite filaments: lessons from Caenorhabditis elegans New J. Phys. 15 075022

[26] Montenegro-Johnson F and Lauga E 2014 Optimal swimming of a sheet Phys. Rev. E 89 060701(R)

[27] Tam D and Hosoi A E 2011 Optimal feeding and swimming gaits of biflagellated organisms Proc. Natl Acad. Sci. USA 108 1001–6

[28] Osterman N and Vilfan A 2011 Finding the ciliary beating pattern with optimal efficiency Proc. Natl. Acad. Sci. USA 108 15727–32

[29] Tam D and Hosoi A E 2007 Optimal stroke patterns for Purcell's three-link swimmer Phys. Rev. Lett. 98 068105

[30] Alouges F, DeSimone A, Heltai L, Lefebvre A and Merlet B 2013 Optimally swimming Stokesian robots Discrete Continuous Dyn. Syst. B 18 1189–215

[31] Alouges F, DeSimone A and Lefebvre A 2008 Optimal strokes for low Reynolds number swimmers: an example J. Nonlinear Sci. 18 277–302
[32] Alouges F, DeSimone A and Heltai L 2011 Numerical strategies for stroke optimisation of axisymmetric micro-swimmers *Math. Models Methods Appl. Sci.* **21** 361–87

[33] Wiezel O and Or Y 2016 Using optimal control to obtain maximum displacement gait for Purcell's three-link swimmer 2016 IEEE 55th Conf. on Decision and Control (CDC) pp 4463–8

[34] Giuliani N, Heltai N and DeSimone A 2018 Predicting and optimizing micro-swimmer performance from the hydrodynamics of its components: the relevance of interactions *Soft Robot.* **5** 410–24

[35] Burton L, Hatton R, Smith D, Choset H and Hosoi A 2010 Two-link swimming using buoyant orientation *Phys. Fluids* **22** 091703

[36] Spagnolie S E and Lauga E 2010 The optimal elastic flagellum *Phys. Fluids* **22** 455

[37] Passov (Gutman) E and Or Y 2012 Dynamics of Purcell's three-link microswimmer with a passive elastic tail *Eur. Phys. J. E* **35** 1–9

[38] Montino A and DeSimone A 2015 Three-sphere low-Reynolds-number swimmer with a passive elastic arm *Eur. Phys. J. E* **38** 1–10

[39] Cicconofri G and DeSimone A 2016 Motion planning and motility maps for flagellar microswimmers *Eur. Phys. J. E* **39** 72

[40] Ishimoto K and Gaffney E 2018 An elastohydrodynamical simulation study of filament and spermatozoan swimming driven by internal couples *IMA J. Appl. Math.* **83** 655–79

[41] Raz O and Avron J 2008 Comment on optimal stroke patterns for Purcell's three-link swimmer *Phys. Rev. Lett.* **100** 029801

[42] Giraldi L, Martinon P and Zoppello M 2015 Optimal design of Purcell's three-link swimmer *Phys. Rev. E* **91** 023012

[43] Wiezel O and Or Y 2016 Optimization and small-amplitude analysis of Purcell's three-link microswimmer model *Proc. R. Soc. A* **472** 20160425

[44] Alouges F and Di Fratta G 2018 Parking 3-sphere swimmer: I. Energy minimizing strokes *Discrete Cont. Dyn. Syst.* **B** **23** 1797–817

[45] Agostinelli D, Alouges F and DeSimone A 2018 Peristaltic waves as optimal gaits in metameric bio-inspired robots *Frontiers Robot. AI* **5** 99

[46] Bonnans F, Martinon P and Grélard V 2012 Bocop—a collection of examples *PhD Thesis* INRIA https://hal.inria.fr/hal-00726992

[47] Alouges F, DeSimone A, Giraldi L and Zoppello M 2013 Self-propulsion of slender micro-swimmers by curvature control: n-link swimmers *Int. J. Nonlinear Mech.* **56** 132–41

[48] Cox R G 1970 The motion of long slender bodies in a viscous fluid: I. General theory *J. Fluid Mech.* **44** 791–810

[49] Giraldi L, Martinon P and Zoppello M 2013 Controllability and optimal strokes for N-link microswimmer 52nd IEEE Conf. on Decision and Control (CDC) pp 3870–5

[50] Avron J E and Raz O 2008 A geometric theory of swimming: Purcell's swimmer and its symmetrized cousin *New J. Phys.* **10** 063016

[51] Gutman E and Or Y 2016 Symmetries and gaits for Purcell's three-link microswimmer model *IEEE Trans. Robot.* **32** 53–69

[52] Koiller J, Ehlers K and Montgomery R 1996 Problems and progress in microswimming *J. Nonlinear Sci.* **6** 507–41

[53] Wächter A and Biegler L T 2006 On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming *Math. Program.* **106** 25–57

[54] Amestoy P R, Duff I S, LExcellent J-Y and Koster J 2001 A fully asynchronous multifrontal solver using distributed dynamic scheduling *SIAM J. Matrix Anal. Appl.* **23** 15–41

[55] Bell B M 2012 CppAD: a package for C++ algorithmic differentiation *Computation. Infrastruct. Operations Res.* **57** 10

[56] Tenenbaum J, Silva V and Langford J 2000 A global geometric framework for non-linear dimensionality reduction *Science* **290** 2319–23

[57] Arroyo M, Milan D, Heltai L and DeSimone A 2012 Reverse engineering the euglenoid movement *Proc. Natl Acad. Sci. USA* **109** 17874–9

[58] Rossi M, Cicconofri G, Beran A, Noselli G and DeSimone A 2017 Kinematics of flagellar swimming in Euglena gracilis: Helical trajectories and flagellar shapes *Proc. Natl Acad. Sci. USA* **114** 13085–90
[32] Alouges F, DeSimone A and Heltai L 2011 Numerical strategies for stroke optimisation of axisymmetric micro-swimmers *Math. Models Methods Appl. Sci.* **21** 361–87

[33] Wiezel O and Or Y 2016 Using optimal control to obtain maximum displacement gait for Purcell's three-link swimmer 2016 IEEE 55th Conf. on Decision and Control (CDC) pp 4463–8

[34] Giuliani N, Heltai N and DeSimone A 2018 Predicting and optimizing micro-swimmer performance from the hydrodynamics of its components: the relevance of interactions *Soft Robot.* **5** 410–24

[35] Burton L, Hatton R, Smith D, Choset H and Hosoi A 2010 Two-link swimming using buoyant orientation *Phys. Fluids* **22** 091703

[36] Spagnolie S E and Lauga E 2010 The optimal elastic flagellum *Phys. Fluids* **22** 455

[37] Passov (Gutman) E and Or Y 2012 Dynamics of Purcell's three-link microswimmer with a passive elastic tail *Eur. Phys. J. E* **35** 1–9

[38] Montino A and DeSimone A 2015 Three-sphere low-Reynolds-number swimmer with a passive elastic arm *Eur. Phys. J. E* **38** 1–10

[39] Cicconofri G and DeSimone A 2016 Motion planning and motility maps for flagellar microswimmers *Eur. Phys. J. E* **39** 72

[40] Ishimoto K and Gaffney E 2018 An elastohydrodynamical simulation study of filament and spermatozoan swimming driven by internal couples *IMA J. Appl. Math.* **83** 655–79

[41] Raz O and Avron J 2008 Comment on optimal stroke patterns for Purcell's three-link swimmer *Phys. Rev. Lett.* **100** 029801

[42] Giraldi L, Martinon P and Zoppello M 2015 Optimal design of Purcell's three-link swimmer *Phys. Rev. E* **91** 023012

[43] Wiezel O and Or Y 2016 Optimization and small-amplitude analysis of Purcell's three-link microswimmer model *Proc. R. Soc. A* **472** 20160425

[44] Alouges F and Di Fratta G 2018 Parking 3-sphere swimmer: I. Energy minimizing strokes *Discrete Cont. Dyn. Syst.* **B** **23** 1797–817

[45] Agostinelli D, Alouges F and DeSimone A 2018 Peristaltic waves as optimal gaits in metameric bio-inspired robots *Frontiers Robot. AI* **5** 99

[46] Bonnans F, Martinon P and Grélard V 2012 Bocop—a collection of examples *PhD Thesis* INRIA https://hal.inria.fr/hal-00726992

[47] Alouges F, DeSimone A, Giraldi L and Zoppello M 2013 Self-propulsion of slender micro-swimmers by curvature control: n-link swimmers *Int. J. Nonlinear Mech.* **56** 132–41

[48] Cox R G 1970 The motion of long slender bodies in a viscous fluid: I. General theory *J. Fluid Mech.* **44** 791–810

[49] Giraldi L, Martinon P and Zoppello M 2013 Controllability and optimal strokes for N-link microswimmer 52nd IEEE Conf. on Decision and Control (CDC) pp 3870–5

[50] Avron J E and Raz O 2008 A geometric theory of swimming: Purcell's swimmer and its symmetrized cousin *New J. Phys.* **10** 063016

[51] Gutman E and Or Y 2016 Symmetries and gaits for Purcell's three-link microswimmer model *IEEE Trans. Robot.* **32** 53–69

[52] Koiller J, Ehlers K and Montgomery R 1996 Problems and progress in microswimming *J. Nonlinear Sci.* **6** 507–41

[53] Wächter A and Biegler L T 2006 On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming *Math. Program.* **106** 25–57

[54] Amestoy P R, Duff I S, LExcellent J-Y and Koster J 2001 A fully asynchronous multifrontal solver using distributed dynamic scheduling *SIAM J. Matrix Anal. Appl.* **23** 15–41

[55] Bell B M 2012 CppAD: a package for C++ algorithmic differentiation *Computation. Infrastruct. Operations Res.* **57** 10

[56] Tenenbaum J, Silva V and Langford J 2000 A global geometric framework for non-linear dimensionality reduction *Science* **290** 2319–23

[57] Arroyo M, Milan D, Heltai L and DeSimone A 2012 Reverse engineering the euglenoid movement *Proc. Natl Acad. Sci. USA* **109** 17874–9

[58] Rossi M, Cicconofri G, Beran A, Noselli G and DeSimone A 2017 Kinematics of flagellar swimming in Euglena gracilis: Helical trajectories and flagellar shapes *Proc. Natl Acad. Sci. USA* **114** 13085–90