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

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.


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,5,6,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,9,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,12,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 travelling 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,18,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,25,26,27] 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,28,29,30,31,32,33].
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 [34] and [35,36,37,38,39].
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 [28] for Purcell's swimmer and further discussed in [40]? 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 non-linear 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 smallangle regime has been used for Purcell's and other swimmers where the number N of shape degrees of freedom is small [41,42,43], 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 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 [44], 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 [45]: we exhibit nonplanar 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.

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 [46], 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 [47] etc.).
The swimmer is planar and composed of N identical segments connected by In particular, we will be concerned with T −periodic shape changes φ(t), which we call strokes, satisfying φ(t) = φ(T +t), ∀t > 0. Our analysis aims at resolving the leadingorder terms of the swimmer's dynamics, and is restricted to up to second order O(ε 2 ) terms in the shape parameters φ i (t) = O(ε), which are assumed to be small, i.e. ε 1. The dynamics of the swimmer's planar motion is governed by These equations come from the balance of viscous force and torque, which are linear inẋ,ẏ,θ andφ. Details of the derivation in the case of resistive force theory [4] are given in [42,46,48]. 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 θ, 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 Fig. 2, transforms the swimmer parametrized by φ at position (x, y, θ) to the one parametrized by −φ, at position (x, −y, −θ). The invariance of the dynamics under such a transformation leads to x-axis

Dynamics at leading order
For θ and φ of order ε, we expand the dynamics (1) to second order in θ, φ andφ by expanding f to first order as where We also expand similarly g and h. Using the symmetry relations (2), we get for f , g and h respectively from which we deduce, φ andφ being arbitrary, that F 0 = G θ = 0 and G φ = H φ = 0. The original system (1) therefore reduces to Integrating in time we obtain and We deduce that, up to second order, the swimmer experiences no global rotation or transverse translation after one complete stroke. Indeed, φ being T -periodic, we have Moreover, using the expression for θ in (5) in the first equation of (4) and integrating over one period leads to .
integrating by parts the last expression, and using again the T -periodicity of φ leads to Finally we get the net longitudinal displacement in one cycle as Remark 2.1 Our analysis at order two in the joint-angle amplitudes proves that the net lateral displacement ∆y and the net rotation both vanish at order two, as given in Eq. (6), while ∆x may be non-zero at order two, as shown in Eq. (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 [49,50,43]. 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 [51].

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 from which one easily obtains F θ = −G 0 (and G θ = F 0 = 0). Similarly, using the symmetry of the system with respect to the y-axis, it can be shown that both F θ and G 0 are 'even' vectors (in the sense that F θ,i = F θ,N −i and similarly for

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φ. Thus, the instantaneous power density expended by the swimmer is a quadratic form inφ, with coefficients depending on φ (see [42] for concrete expressions using resistive force theory). We may therefore write the total energy E expended by the swimmer during the stroke as where P(φ) ∈ M N,N (R) is a symmetric and positive definite matrix, and the bracket < ·, · > stands for the scalar product in R N . Expanding this expression gives the O(ε 2 ) leading-order term of the energy as where P 0 := P(φ = 0) is symmetric and positive definite.

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.

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 drag /E swim , where E drag is the energy needed by an external actor to pull the swimmer during a time T at an average speed ∆x/T , while E 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 drag is proportional to (∆x) 2 /T whereas E swim behaves like ∆x as (7) and (8) indicate (∆x and E being quadratic in (φ,φ) are of the same order O(ε 2 )). Thus, Lighthill's efficiency is proportional to ∆x, and maximizing it amounts to increase ∆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 [28,42]. 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 φ(t) that expend the minimal energy E, among all strokes that achieve a given displacement ∆x in a given time period T . We showed in the previous section that, at leading order, the y and θ displacements are negligible with respect to the x displacement. Optimal strokes that provide a (longitudinal) displacement ∆x are thus sought as solutions to the constrained minimization problem where, using (7), we have 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)  In order to find the Euler-Lagrange first-order necessary condition associated with the constrained optimization problem (10), we set T )φ ·φ dt and λ ∈ R is the Lagrange multiplier associated with the constraint K 0 (φ) = ∆x. Writing that the first order functional derivative ofẼ 0 vanishes at the optimal stroke φ * δẼ 0 (φ * ) = 0 amounts to writing where the Lagrangian L is such that Here, we have and therefore (11) becomes We denote by M the skew symmetric matrix 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 (v ± j ) 1≤j≤N (N = N 2 ) the (complex and orthonormal) eigenvectors associated with the purely imaginary eigenvalue ±iµ j with µ j ≥ 0: Projecting P we deduce from (12) j is a constant that we may take equal to 0. The solution of (12) is thus expressed as Since we focus on periodic strokes, φ(0) = φ(T ), we must have By plugging the solution (14) into (9), we find while the x-displacement is given by Since ∆x is fixed in the optimization problem, minimizing the energy requires choosing λ as small as possible. In view of (15), this is achieved if λ = ± 2π µ M T , where µ M = max {µ j , 1 ≤ j ≤ N }, the direction of the translation depending on the sign of λ, and ψ ± j = 0 if µ j = µ M . Assuming λ > 0 (λ < 0 is handled similarly), we deduce that λ = 2π/T µ M and the solution has only two modes and is therefore an ellipse drawn in the plane (P The net-displacement achieved by the gait (18) will be ∆x = 2T λ |α + M | 2 . Therefore, the optimal gait amplitude scales as √ ∆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: φ * k (t) = a k sin(ωt + p k ), k = 1, . . . , N −1 where ω = 2π/T . The amplitudes a k in (19) are of order O(ε), and scale as √ ∆x. Finally, an important property of the optimal solution φ * (t) is that it satisfies 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 [30] 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.

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,47] for slender links. Our computations are made using Matlab. A similar approach based on the the Euler-Lagrange equations (11) was used in [30] 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 φ * (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 [47,4]. This theory states that the viscous drag force f i and torque m i on the i th 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: where v i , ω i are the linear and angular velocities of the i th link, and t i , n i are unit vectors in its axial and normal directions. The resistance coefficients in (21) are c n = 2c t = 4πη/ log(l/a) where η is the dynamic viscosity of the fluid and a 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 [42] for Purcell's swimmer and [46] for the general N -link swimmer.
Next, we exploit the assumption of small amplitude approximation and we calculate the matrices F θ φ 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 φ * (t) in (18). This computation is done numerically in Matlab.

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 energyoptimal gait, which is shown in Fig. 6a as a trajectory in the plane of joint angles (φ 1 , φ 2 ) for ∆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 Fig. 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.

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 3a shows a single snapshot of both swimmers under the optimal gait for a prescribed displacement ∆x = 0.1l. The figure indicates that the optimal gaits look like a travelling 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 where p k are given by (19). The results of amplitudes and phase differences across the joints for both swimmers are plotted in Figures 3b and 3c, 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 Fig  3b. The phase difference between joints is very close to a uniform value of 120 • , indicating a travelling wave with wavelength of three links. This is further confirmed when computing the optimal gait for increasing numbers of N links for ∆x = 0.01l. Figures 4a and 4b 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 ∆ϕ = 2π/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 N l). In particular, we look for strokes that are periodic in both space (with period N l) 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 θ φ − F θ φ 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 (18) shows 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 [44], where the optimality of actuation strategies based on peristaltic waves is discussed.  Finally, we fix the total length L = N l of a swimmer, compute the energyoptimal gait for moving a given displacement of ∆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 * ∼ ∆xL 2 N for large N . This suggests that the a more suitable performance measure would be the scaled optimal energy defined as Q = N E * ∆xL 2 . This quantity, (multiplied by ∆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 nonzero value at the limit of large N .

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 [45]. Bocop applies numerical integration of the full nonlinear dynamic equations (1) and the energy formula (8), without assuming smallamplitude 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 [52] with Mumps [53], while the derivatives are computed by sparse automatic differentiation with CppAD [54]. In the numerical experiments, we used a midpoint (implicit 1st order) discretization with 100 time steps.

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 Fig. 6a as a trajectory in the plane of joint angles (φ 1 , φ 2 ) for ∆x = 0.1l. Snapshots of the swimmer's configuration at quarterperiod intervals are also illustrated on the plot.
For comparison, we computed energy-optimal gaits using Bocop for different displacements ∆x, and the resulting trajectories scaled by √ ∆x are shown in Fig. 6b. It can be seen from the plot that for small displacements ∆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 ∆x, without a specified direction. The agreement between the two methods of computation confirms the validity of our small-amplitude analysis. When the displacement ∆x is increased further, the energy-optimal gaits obtained with Bocop begin to deviate significantly from the small-amplitude one. For ∆x = 0.26, the optimal gait (dashed) coincides with the gait obtained by Tam and Hosoi [28] that 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 [28] does not involve a constraint on the travelled distance. Therefore, it produces large displacements through large amplitude strokes, as discussed at the beginning of Section 3.
When ∆x is further increased to an upper limit of ∆x = 0.306, the optimal gait (dash-dotted) deforms into the famous peanut-shaped loop obtained in [28] 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 [42]. The sinusoidal gait in (19) for the two joint angles can be rewritten as That is, it has a stroke amplitude of ε and phase difference of ϕ. Using resistive force theory and the leading order expansion in ε, as explained in [42], the Therefore, for sinusoidal gaits of the form (23), our constrained optimization of minimizing the energy to cover a given displacement reduces to minimizing E/∆x. This gives an optimization problem for a scalar function of ϕ. Using elementary calculus, the minimum of this function is obtained at the optimal phase difference of ϕ * = cos −1 (−11/16) = 133.43 • , and the resulting optimal gait (23) is exactly identical to the one obtained by using the eigenvalue formulation in (18), which is shown in Fig. 6a.

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 ∆x = 0.1l are shown in Figure 7 and compared to the numerical solutions obtained for ∆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 µ M .) In order to compare further our analytical optimal gaits with the results of Bocop computations, we first plot in Figure 8a 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 ∆x and then scaled by √ ∆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 R 4 from the plane S. Figure 8b shows a log-log plot of this distance d as a function of ∆x. It can be seen that for small displacements this distance decays to zero as (∆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 Fig. 7 and in Fig. 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 Fig. 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 Fig. 7h, 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 nonplanarity when angles and displacements become large, we derive from the computed four-dimensional loops a dimensionally-reduced representation in three dimensions, using the Isomap approach [55,56]. 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 ∆x=0.1l ∆x=0.316l ∆x=0.548l ∆x=0.707l ∆x=0.894l ∆x=1l Analytic  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.

Conclusion
In this work, we have studied optimal periodic strokes of multi-link microswimmers 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 twodimensional 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 Figure 10: Three-dimensional representation of the optimal loops for the 5-linkswimmer as a function of the amplitude of the prescribed displacement ∆x, in normalized coordinates (π 1 , π 2 , π 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 non-linear projection of the origianl 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.
Taylor's classic work [1]. A noticeable difference from [1] due to the 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 multilink 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 nonconvex 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,57], generalization of the model to include elasticity as in [37,35,36,38], and the study of optimal control of magnetically-actuated microswimmers [10]. Fixing the displacement ∆x, it is deduced from (26) that for large N , the amplitude ε increases as N ∆x L . Moreover, the energy scales as E ∼ ∆xL 2 N T . This scaling relation explains the decay rate of E * in Figure 5.