Optimal Strokes of Low Reynolds Number Linked-Sphere Swimmers

Optimal gait design is important for micro-organisms and micro-robots that propel themselves in a fluid environment in the absence of external force or torque. The simplest models of shape changes are those that comprise a series of linked-spheres that can change their separation and their sizes. We examine the dynamics of three existing linked-sphere types of modeling swimmers in low Reynolds number Newtonian fluids using asymptotic analysis, and obtain their optimal swimming strokes by solving the Euler-Lagrange equation using the shooting method. The numerical results reveal that (1) with the minimal 2 degrees of freedom in shape deformations, the model swimmer adopting the mixed shape deformation modes strategy is more efficient than those with a single-mode of shape deformation modes, and (2) the swimming efficiency mostly decreases as the number of spheres increases, indicating that more degrees of freedom in shape deformations might not be a good strategy in optimal gait design in low Reynolds number locomotion.


Introduction
Swimming by shape changes at low Reynolds number (LRN) is widely used in biology and micro-robotic design. In this flow regime, inertial effects are negligible, and the micro-organisms or micro-robots propel themselves by exploiting the viscous resistance of the fluid. For example, while a scallop that can only open or close its shell can swim in the ocean by accelerating the surrounding water, such a swimming strategy does not work at LRN, which is generalized by the principle: any time-reversible swimming stroke leads to no net translation at LRN Newtonian fluid, known as the scallop theorem [1].
It is important to understand how the performance of swimming depends on the geometric patterns of shape deformations for micro-swimmers. In microbiology, to fight the viscous resistance, different microorganisms adopt various propulsion mechanisms and directed locomotion strategies for searching for food and running from predators. For example, individual cells such as bacteria find food by a combination of taxis and kinesis using a flagellated or ciliated mode of swimming [5,4,2,3]. Recently, it was discovered that Dd cells can occasionally detach from the substrate and stay completely free in suspension for a few minutes before they slowly sink; during the free suspension stage, cells continue to form pseudopods that convert to rear-ward moving bumps, thereby propelling the cell through the surrounding fluid in a totally adhesion-free fashion [6]. Human neutrophils can swim to a chemoattractant fMLP (formyl-methionylleucyl-phenylalanine) source at a speed similar to that of cells migrating on a glass coverslip under similar conditions [7]. Most recently and equally striking, Drosophila fat body cells can actively swim to wounds in an adhesion-independent motility mode associated with actomyosin-driven, peristaltic cell shape deformations [8]. In micro bio-engineering, medical microrobots revolutionize many aspects of medicine in recent years, which make existing therapeutic and diagnostic procedures less invasive [9]. Different LRN swimming models and micro-robots have been designed since Purcell's two-hinge model was advanced [1]. In particular, various linked-sphere types of models have appeared, since their simple geometry permits both analytical and computational results [10,12,11,15,14,16,17,13]. These analytical and numerical results have greatly inspired the designs of micro robotic devices, for example, swimmers with 2 and 3 rotatory cylinders have been built to study the hydrodynamic interaction between a wall and an active swimmer [18,19]. Other micro-robots inspired by analytical/numerical works include the Quadroar swimmer, which consists of rotating disks and a linear actuator [20], as well as the Purcell's two-hinge model [21].
An important problem in LRN swimming is to find the optimized swimming stroke of the microswimmer, either (1) with respect to time, i.e., the stroke in one swimming cycle that moves the cell farthest, or (2) with respect to energy, i.e., among all strokes with designated starting and end points, find the one that consumes the least energy. These are usually called the time optimal control and the energy optimal control problems, respectively. Both optimal problems have attracted substantial interest in optimization as well as geometry [30,27,33,22,23,25,26,29,24,31,32,28]. Moreover, recently techniques from machine learning have been applied to linked-sphere types of model gait design, which allows incorporating environmental influences on the micro-swimmer's swimming behavior, including noise and a frictional medium [34].
Here by investigating the optimal strokes of a group of linked-sphere types of LRN modeling swimmers, we study the efficiencies of propelling mechanisms at LRN of different types of shape deformations. We start from three linked-sphere swimmers-Najafi-Golestanian (NG) 3-sphere accordion model ( Figure  1a) [10,12,11], pushmepullyou (PMPY) 2-sphere model (Figure 1b) [15], and the volume-exchange (VE) 3-sphere model (Figure 1c) [17]. All three models have only 2 degrees of freedom in their shape deformations, which is a minimal requirement that enables the swimmer to propel itself at LRN, according to the scallop theorem [1]. In particular, the shape deformations in the three models can be generalized as body elongation and/or mass transportation, and we investigate the efficiencies of these two shape deformation modes on the swimming performances. Then we generalize our results into a linear chain of spheres. The results are presented as follows: in Section 2 we present a brief introduction of the LRN swimming problem; in Section 3 we review the three existing linked-sphere types of swimmers: NG 3-sphere, PMPY 2-sphere and VE 3-sphere models and discuss their optimization problems in Section 4; finally in Section 5 we discuss the optimization problem of models consisting of a chain of spheres. Figure 1: Linked-sphere LRN swimmers. (a) a NG 3-sphere model [10,12,11]. (b) a PMPY 2sphere model [15]. (c) a VE 3-sphere model [17]. Figures show example strokes of each swimmer. Only shape deformations are presented, net translation of the swimmer is not shown.

Swimming at LRN by Shape Changes-the Exterior Problem
The Navier-Stokes equations for an incompressible fluid of density ρ, viscosity µ, and velocity u are where f ext is the external force field. Herein we assume that the swimmer is self-propelled and does not rely on any exterior force, and therefore we require that f ext = 0. The Reynolds number based on a characteristic length scale L and speed scale U is Re = ρLU/µ, and when converted to dimensionless form and the symbols re-defined, the equations read Here Sl = ωL/U is the Strouhal number and ω is a characteristic frequency of the shape changes. When Re 1 the convective momentum term in Equation (1) can be neglected, but the time variation requires that ReSl = ωL 2 /ν 1. When both terms are neglected, which we assume throughout, the flow is governed by the Stokes equations: We also consider the propulsion problem in an infinite domain and impose the condition u| x→∞ = 0 on the velocity field.
In the LRN regime time does not appear explicitly, momentum is assumed to equilibrate instantaneously, and bodies move by exploiting the viscous resistance of the fluid. As a result, time-reversible deformations produce no motion, which is the content of the scallop theorem [1]. For a self-propelled swimmer, there is no net force or torque, and therefore movement is a purely geometric process: the net displacement of a swimmer during a stroke is independent of the rate at which the stroke is executed, as long as the Reynolds and Strouhal numbers remain small enough.
Suppose that a swimmer occupies the closed compact domain Ω(t) ⊂ R 3 , at time t, and let ∂Ω(t) denote its prescribed time-dependent boundary. A swimming stroke γ is specified by a time-dependent sequence of the boundary ∂Ω(t), and it is cyclic if the initial and final shapes are identical, i.e., ∂Ω(0) = ∂Ω(T ) where T is the period [35]. The swimmer's boundary velocity V relative to fixed coordinates can be written as a part v that defines the intrinsic shape deformations, and a rigid motion U. If u denotes the velocity field in the fluid exterior to Ω, then a standard LRN self-propulsion problem is: Given a cyclic shape deformation γ(t) that is specified by v, solve the Stokes equations subject to: Force-free condition: Torque-free condition: In order to treat general shape changes of a swimmer defined by Ω(t) ⊂ R 3 with boundary ∂Ω(t), one must solve the exterior Stokes Equations (2) for u, with a prescribed velocity v(t) on ∂Ω(t) and subject to the decay conditions u ∼ 1/r and p ∼ 1/r 2 as r → ∞. The solution has the representation: where G is the free-space Green's function, T is the associated third-rank stress tensor, n is the exterior normal and f is the force on the boundary [36]. The constraints that the total force and the total torque vanish determine the center-of-mass translational and angular velocities. When x ∈ ∂Ω(t), this is an integral equation for the force distribution on the boundary, the solution of which determines the forces needed to produce the prescribed shape changes.

Linked-Sphere LRN Swimmers
To date, various simple linked-sphere models for which both analytical and computational results can be obtained have appeared. In particular, there are three linked-sphere models that are designed to investigate the effects of changing body length vs. mass transportation on LRN swimming behavior. The three models are: the NG 3-sphere accordion model ( Figure 1a) [10,12,11], the PMPY 2-sphere model ( Figure 1b) [15], and the VE 3-sphere model (Figure 1c) [17].

Fundamental Solutions for Translating and Radially Deforming Spheres
In either of the three models, shape deformations either come from linking-arm length change (the NG and the PMPY models) or radial deformations of the spheres (the PMPY and the VE models). When a sphere of radius a(t) is pulled through a quiescent fluid with a steady force f under no-slip conditions at the surfaces, the resulting flow field is given by [37]: where x 0 is the position of the center of the sphere and r = x − x 0 , and the resulting velocity of the sphere is given by Stokes' law A second fundamental solution is the velocity field u produced by a radially-expanding sphere, which can be generated by a point source at the center x 0 of the sphere. The corresponding velocity is [36]: where v = 4πa 3 /3 is the volume of the sphere andr = r/r. By combining the two fundamental solutions (3) and (5), we obtain the solution of the flow velocity generated by a sphere of radius a, centered at x 0 , subject to a pulling/pushing force f , and radially deforming at the rateȧ: Next, the fundamental solution involves the interaction between two spheres. Suppose that at time t the ith (i = 1, 2) sphere has radius a i (t), centered at x i (t), subjected to force f i (t) and with radially deforming rateȧ i (t). In the scenario that the two spheres are far apart, i.e., a i /l 1, where l = |x 1 −x 2 |, the translational velocity of the ith sphere can be obtained by the Reflection Method [37,41,39,42,38,40]: where δU t i is due to the translation of the other sphere, which arises from the flow given by Equation (3) and is where l = x i − x j . The velocity δU e i is resulted from a flow generated by the radial deformation of the other sphere, defined by Equation (5): Finally, the power consumption P of a sphere with radius a, dragged by a force f , translating at velocity U and radially expanding/contacting at a rate ofȧ comprises two parts: P t = f · U that results from the drag force on the sphere, and P e = 16πµaȧ 2 that results from the radial deforming [15,40]. Therefore the total power expended is:

Non-Dimensionalization of the System
Let A, L be the characteristic sphere radius and link-arm length. The reflection method holds in the regime ε := A/L 1. Shape deformations of the linked-sphere models (see Sections 3.3-3.5) come from the expanding/contracting of the link-arms (l), or the radial expanding/contracting of the spheres (ȧ).
For the shape deformation scales, we require thatȧ/l ∼ O(1), and |δa|, |δl| A, where δa and δl denote the shape deformations in a and l, respectively.
We non-dimensionalize the radii and rod lengths by A and L: We non-dimensionalize other length scales by A as well, and time by T -the swimming stroke period, thus we obtain: We assume that the deformations δa, δl in a and l are of the same order, both of which are scaled by and we have the approximation: The drag force f exerted on a sphere of radius a in a quiescent fluid is related to the sphere velocity U via Equation (4), which leads to the following scaling for forces: Finally, we non-dimensionalize the power P as so that while Power = Force · Velocity in dimensional form, after non-dimensionalization we still have P * = f * · U * . Thus the non-dimensional version of Equation (9) is:

NG 3-Sphere Swimmer
The NG swimmer ( Figure 1a) consists of three spheres with fixed radii a i (i = 1, 2, 3) and two linkingarms with adjustable length l i (t) (i = 1, 2) [10, 12,11]. While the spheres are rigid, the linking-arms can stretch or contract. The geometry of the model is shown in Figure 1a, where the three spheres align along the x-axis, numbered from sphere 1 to 3 from left to right. Let e be the unit vector pointing in the positive x direction, and we have e = (x i − x j )/ x i − x j for any pair of i, j with i > j. Thus we have f i = f i e and U i = U i e for i = 1, 2, 3. For simplicity, we consider the simple geometry a 1 = a 2 = a 3 = A. Assume initial shape l 1 (0) = l 2 (0) = L. The (non-dimensional) velocity of each sphere is given by: is the scaled distance between spheres i and j. The velocities are related via the following relations: The system is force-free, ad the velocity and power of the swimmer are: The PMPY swimmer (Figure 1b) consists of two spheres with radii a i (t) (i = 1, 2) and one linking-arm with length l(t) [15]. The spheres can expand or contract in the radial direction, and the linking-arm can stretch or contract. The (non-dimensional) velocity of each sphere is given by: The velocities are related by: and the force-free condition is: Assume initial conditions a 1 (0) = a 2 (0) = A and l(0) = L. There is no mass exchange between the swimmer and surrounding fluid, thus the total volume of the two spheres is conserved: The velocity and power of the whole swimmer are:

VE 3-Sphere Swimmer
The VE swimmer (Figure 1c) consists of three spheres with radii a i (t) (i = 1, 2, 3) and two rigid linkingarms with fixed length l i (i = 1, 2) [17]. The spheres can only expand or contract in the radial direction, and each can only exchange mass with its neighboring sphere(s), that is, between sphere 1 and 2, or between sphere 2 and 3, but not between sphere 1 and 3. The (non-dimensional) velocity of each sphere is given by: Assume initial conditions a 1 (0) = a 2 (0) = a 3 (0) = A and l 1 = l 2 = L. The system is force-free and total volume is conserved: The velocity and power of the whole swimmer are: Notice that the linking-arms are rigid.

Euler-Lagrange Equation for Optimal Strokes of LRN Swimmers
(Hereafter we will continue our discussion using the non-dimensional quantities but omitting the * notation for simplicity). As is introduced in Section 2, a swimming stroke γ(t) is a t-series of swimmer shapes.
In the case when the swimmer has only finitely many degrees of freedom in their shape deformations, γ can be considered as a function from [0, T ] to R m , where m is the number of degrees of freedom of the swimmer, and each component in γ(t) defines one degree of freedom. For example, all three linked-sphere swimmers in Sections 3.3-3.5 have only 2 degrees of freedom, thus we have γ NG (t) = (l 1 (t), l 2 (t)) T for a NG swimmer, γ PMPY (t) = (l(t), a 1 (t)) T for a PMPY swimmer, and γ VE (t) = (a 1 (t), a 3 (t)) T for a VE swimmer.
The linear and driftless properties of Stokes flows make the LRN swimming problem a classic driftless controllable system [35,43,31,32,28,25,44]. Consider a cyclic stroke γ(t) of a LRN swimmer, where t ∈ [0, 1], γ(0) = γ(1). The net translation X(γ) and energy dissipation E(γ) within a stroke can be represented as: where for NG, PMPY and VE equations, the operators F and G can be calculated from the equations from Sections 3.3-3.5, respectively. For a given initial shape which we denote by γ 0 and a given net translation X 0 , let Γ{γ 0 , X 0 } be the set of all cyclic strokes that starts and ends with shape γ 0 and results in a net translation of X 0 in a cycle: The optimal stroke problem can be formulated as follows: Given an initial shape γ 0 and a net translation X 0 , find the stroke γ in Γ{γ 0 , X 0 } that minimizes the energy dissipation E: The Euler-Lagrange equation for this optimization problem is [22,23]: where m is the number of degrees of freedom in the swimmer's shape deformation, and γ = (γ 1 , γ 2 , · · · , γ m ) T ∈ R m . Solutions of Equation (13) give the geodesics connecting (γ 0 , 0) and (γ 0 , X 0 ), which we refer to as the geodesic strokes. The optimal stroke is the geodesic stroke that consumes the least energy [22,23].
We use the Shooting Method (SM) to numerically solve the governing Equation (13) for the optimal stroke problem. The numerical methods are developed in [22,23] and have been applied for NG with equal sizes and PMPY swimmers; recent research works also report the efficiency of NG swimmers with unequal sized sphere [45]. In the following we use SM for all three (NG, PMPY and VE) linking-arm swimmers and compare their swimming performances. In Appendix A we present an outline of SM for the convenience of the readers of this paper.

Numerical Results
The numerical results of optimal strokes of a NG 3-sphere, PMPY 2-sphere and VE 3-sphere swimmers are given in Table 1 and Figure 2, where we take ε = A/L = 0.2 in all simulations. For different values of the net translation X, the optimal strokes γ of the three swimmers are computed by the SM method described in Section 4.1 and Appendix A, the energy dissipation of the optimal strokes E(γ) is recorded in Table 1, and the γ-paths (solid lines) and the X(γ(t)), t ∈ [0, 1] trajectories (dashed lines) are given in Figure 2. The simulation results show that PMPY is the most efficient swimmer among the three, next is NG, the last is VE: where we define the efficiency of an LRN swimmer as the net translation to energy dissipation ratio of the optimal stroke of the swimmer: where γ is the optimal stroke.
The dimensional form of Eff has the unit of force −1 . In addition, from Table 1 we find that for single loop optimal strokes, the efficiency for each swimmer is about the same for different values of X(γ), with a slight increase in the efficiency for NG and PMPY swimmers as X(γ) increases: As X increases, geodesic strokes of multiple loops show up for NG and PMPY swimmers. However, we do not catch multi-loop solutions in our numerical simulations for VE swimmers, as large shape deformations of VE models easily lead to negative volume of the spheres. Comparing to single loop geodesic strokes, multi-loop geodesic strokes are less efficient (see Appendix B).
We summarize our main conclusions from the above numerical results as follows: 1. A PMPY model that adopts a mixed-mode of shape deformations is the most efficient among the three; next is NG and the last is VE, both of which adopt a single-mode of shape deformations (see Remark below).
2. Single-loop geodesic strokes are more efficient than multi-loop geodesic strokes.
3. The efficiency of a given LRN swimmer is almost the same for different X, as long as the optimal stroke is single-looped.
Finally, we numerically compare the optimal strokes of the three model swimmers to stepwise square strokes. The results are presented in Appendix C, which clearly shows the improvement of swimming efficiency in the optimal strokes comparing to the stepwise square strokes.
Remark. At LRN, different microorganisms adopt various propulsion mechanisms and directed locomotion strategies for searching for food and running from predators. While many microorganisms including bacteria use a flagellated or ciliated mode of swimming [53,47,48,51,49,46,5,4,2,3,50,52], some cells use an amoeboid swimming strategy, which relies on the generation, protrusion and even travel of pseudopodia or blebs [7,54,6]. Such shape deformations propagating over the whole cell body can usually be generalized as a combination of two modes: stretching of the cell body, and mass transportation along the cell body. Comparing with the linked-sphere types of models, cell body stretching can be considered as model while mass transportation as modeȧ. According to the scallop theorem [1], a minimal LRN swimmer has at least two degrees of freedom in its shape deformations, therefor we have three possible combinations: two single-mode of shape changes: both in body stretching (both arel-type of shape deformations), both in mass transportation (both areȧ-type of shape deformations), and a mixed-modes of shape change: body stretching combined with mass transportation (l,ȧ). It is clearly seen that the three combinations correspond to the three linked-sphere swimmers (NG, VE and PMPY) we discussed here, and our numerical results indicate that:  (Figure 2a Comparing to single-mode shape deformations, mixed-modes of shape deformations, i.e., body stretching combined with mass transportation is more efficient when swimming at LRN. An asymptotic analysis study on the three swimmers that derives the same conclusion can be found in [55], which illustrates that a PMPY swimmer has a net translation X(γ) ∼ O(1), comparing to X(γ) ∼ O(ε 2 ) for a NG or a VE swimmer. Therefore the mixed control strategy as is adopted by a PMPY swimmer is superior to uni-mode of control strategies as are adopted by a NG or a VE swimmer.

Remark.
A more popular definition of efficiency of NG swimmer in its dimensional form is given by [30]: which measures the ratio of the average energy dissipation in the stroke γ to that applied to a rigid sphere of radius A traveling the same distance X. However, this definition does not apply well to our discussion in that (1) we are considering non-rigid spheres where the spheres are allowed to deform radially, and (2) we are comparing swimmers of different numbers of spheres. Therefore we adopt the definition given by Equation (14) instead.

Optimal Strokes of (m + 1)-Linked-Sphere LRN Swimmers
In this section we will discuss (m + 1)-linked-sphere LRN swimmers, starting from the NG-type of swimmers (Section 5.1), where rigid spheres are considered and shape deformations are only allowed to be in linking-arms length changes. Next we consider the PMPY-type of swimmers (Section 5.2), where both spheres and linking-arms are allowed to deform. Finally we discuss the effect of sphere separation distance on the swimmer's swimming efficiency (Section 5.3).
According to the scallop theorem [1], a minimal LRN swimmer should have at least 2 degrees of freedom. When it comes to LRN swimmers with more than 2 degrees of freedom, the question is: is it better to use all degrees of freedom or to take a subset of them so to reach more efficient swimming strategy? In our simulations, we find that the optimal strokes for different numbers of linking-arms have all made use of all degrees of freedom, for example, Figure 3d,e show the γ-paths of the optimal strokes for 5 and 10 arms, respectively, which show that all linking-arms performed length deformations in the optimal stroke, none is disabled.
The energy dissipation of the optimal stroke reaches a minimum min m E(γ) = 1.70 × 10 −2 at m = 1, i.e., a 2-sphere PMPY swimmer, which is the one we discussed in Section 3.4. E(γ) increases as m increases (Figure 3b, red dots). Correspondingly, the most efficient swimmer is the 2-sphere NG swimmer with Eff = 29.5 × 10 −3 , and the efficiency decreases as m increases (Figure 3c, red dots). Comparing to NG type of swimmers with the same structure (i.e., same numbers of spheres and linking-arms), PMPY swimmers consumes less energy and is thus more efficient (Figure 3b,c).
Similar to NG swimmers, for PMPY swimmers, the optimal strokes for different numbers of linkingarms have also made use of all degrees of freedom. Figure 3fg show the γ-paths of the optimal strokes for 5 and 10 arms, respectively, which show that all linking-arms performed length deformations and all spheres performed radial deformations in the optimal stroke. In addition, we find that for large m, the amplitude of sphere radial deformations a i (t) are smaller than that of the linking-arms length changes l i (t) (Figure 3f,g); on the other hand, in the most efficient PMPY type of swimmers, i.e., the 2-sphere PMPY swimmer, the amplitude of sphere radial deformation a 1 (t) is larger than that of the linking-arm length change l(t) (Figure 2a,c).

(m+1)-Linked-Sphere NG and PMPY Swimmers With Widely Separated Spheres
In this part we consider swimmers with widely separated spheres, where ε = A/L = 0.01, comparing to ε = 0.2 in Sections 5.1 and 5.2. The simulation results are given in Figure 4 for NG and PMPY types of swimmers with different numbers of linking-arms m, where X = 5 × 10 −4 as before.
First, the most efficient NG type of swimmer is again the 4-sphere NG swimmer, with a minimum energy dissipation min m E(γ) = 5745.09 × 10 −2 ; the most efficient PMPY type of swimmer is again the 2-sphere PMPY swimmer, with a minimum energy dissipation min m E(γ) = 2.03 × 10 −2 . Comparing to ε = 0.2 swimmers, we find that the optimized energy dissipation of NG swimmer with widely separated sphere increases greatly while that of PMPY swimmer is almost the same (Table 2).  Next, similar to ε = 0.2 swimmers, when m increases, the energy dissipation E(γ) increases for both NG and PMPY types of swimmers (Figure 4a), resulting in fast decreasing efficiency (Figure 4b). However, we should point out that while for a 2-sphere PMPY swimmer, E(γ) is almost the same despite the separation distance of the spheres (Table 2), for large m, the energy dissipation of PMPY type of swimmers increases much faster with widely separated spheres (ε = 0.01) comparing to those with closer spheres (ε = 0.2) (Figures 3b and 4a, Table 3). Remark. In [55], using asymptotic analysis, it is shown that for NG 3-sphere and VE 3-sphere swimmers, the net translation X(γ) ∼ O(ε 2 ); while for PMPY 2-sphere swimmers, X(γ) ∼ O(1). Therefore as the separation distance L increases, the efficiency of NG and VE 3-sphere swimmers quickly decreases, while that of PMPY 2-sphere does not change much, which makes PMPY swimmers superior to NG and VE swimmers from an efficiency point of view, considering that all have a minimal 2-degree of freedom in shape deformations. The asymptotic analysis results supports our numerical results for optimal strokes (Table 2). However, when there are more degrees of freedom in a swimmer's shape deformations, the efficiency behavior of PMPY types of swimmers clearly depends on the separation distance L (Table 3). Further asymptotic analysis in this direction might be worthwhile.

Discussion
A successful and efficient locomotory gait design is important for micro-organisms and micro-robot who live or present in a viscous fluid environment. While bacteria often adopt a flagellated or ciliated mode of swimming strategy, amoeboid cells deform their cell bodies to propel themselves, resisting the viscous resistance from surrounding fluid. Such shape deformations can be generalized into two modes: stretch or elongation of a part of or the whole cell body, and mass transportation along the the cell body which does not greatly change the cell length. Interestingly, in bio-engineering, there are three linked-sphere LRN swimming models (NG 3-sphere, PMPY 2-sphere and VE 3-sphere models) that adopt a uni-or mixed modes of the aforementioned two shape deformation changes. By analyzing the optimal strokes of the three swimmers, we show that PMPY which adopts the mixed control is the most efficient among the three. We also consider models consisting a chain of spheres, and again we find that the PMPY-type of swimmers that use mixed controls are more efficient than NG-type of swimmers which use uni-controls of length change. We also find that generally speaking, the swimming efficiency decreases as the number of spheres increases, implying that more degrees of freedom in shape deformations is not a good strategy in optimal gait design. When the sphere separation distance increases, the efficiencies of NG type of swimmers greatly decrease, while the efficiencies of PMPY type of swimmers decrease in swimmers with many spheres, but are not affected much for swimmers with less spheres. The findings in our paper can be potentially applied in the design of micro-robots with more complex structures. In addition, the numerical scheme presented here can be applied to more advanced LRN swimming systems. For example, it can be applied for general 2D and 3D swimmers, when the swimmer shapes can be represented by conformal mappings or spherical harmonics [35,58,57,56].
Moreover, for micro-organisms, their swimming behaviors are also significantly affected by environmental factors in biological media, therefore complex rheology of the surrounding fluid should be considered. A starting point might be bringing the linked-sphere swimmers into viscoelastic medium. We would like to point out that asymptotic analysis results for NG 3-sphere swimmer in linearized viscoelastic fluid have been obtained [13].

A Shooting Method
We use the SM to solve the optimal stroke problem Equation (12), where the Euler-Lagrange Equation (13) is solved using the 2nd order Runge-Kutta method [22,23]. For a LRN swimmer with finitely many degrees of freedom in its shape deformations, given an initial shape γ 0 and an initial shape deformation denoted by v 0 , there exists a unique solution γ(t) (t ∈ [0, 1]) to the Euler-Larange Equation (13) and satisfies the initial condition: γ(0) = γ 0 ,γ(0) = v 0 . For the given initial shape γ 0 , define a function Φ γ0 as: where λ is the Lagrange multiplier in Equation (13), γ is the solution to the Euler-Larange Equation (13), and X(γ) is the net translation resulted from the stroke γ (see Equation (10)). Then the optimal stroke problem Equation (12) becomes: Given initial shape γ 0 and a net translation X 0 , find The corresponding stroke is the optimal stroke γ.
Remark. From the geometric point of view, let Q be the set of all possible configurations of the LRN swimmer. If the swimmer has only finitely many degrees of freedom in this shape deformations, then Q is a finite dimensional smooth manifold and has the geometric structure Q = S × G, where S is the set of all unlocated shapes of the swimmer and G is the group of rigid motions. All the linked-sphere swimmers we considered are rotation-free and can only perform translation along one direction, therefore G ∼ = R. A swimming stroke γ(t) is a path in S, andγ(t) ∈ T γ(t) S, that is, a shape deformationγ(t) is a tangent vector in the tangent space T γ(t) S on the shape γ(t). Therefore Equation (15) defines a function Φ γ0 : The geometric structure Q = S×G makes Q the trivial principal fiber bundle over S. For more discussions regarding the geometric properties of a LRN swimming system, please refer to [35,43,31,32,28,25,44].

B Geodesic Strokes
In [22], Alouges et al identified three types of geodesic strokes by their shapes, referred to as "drop, bean and pretzel". While "drops" and "beans" are single loop strokes, pretzel can be considered as the general 2-loop strokes. In our simulations of NG 3-sphere and PMPY 2-sphere models, we also capture the three types of geodesic strokes ( Figure 5, Tables 4 and 5), where X = 0.045, ε = 0.2 for the NG swimmer and X = 0.1, ε = 0.2 for the PMPY swimmer. In the NG swimmer, the "drop" stroke gives the optimal stroke among the three (Figure 5a,b, Tables 4), and we note that in [22] it is also reported that the "drop" stroke is the most efficient one. On the other hand, in the PMPY swimmer, the "drop" and "bean" strokes consume similar energy (Figure 5c,d, Tables 5). In either NG or PMPY model, single loop geodesic strokes are clearly more efficient then the 2-loop geodesic strokes, or say, the "pretzel" strokes (Tables 4 and 5). Next, with large X, we catch multi loop geodesic strokes with more than 2 loops in the PMPY swimmer ( Figure 6, Table 6). Finally, We do not catch multi-loop geodesic strokes for VE swimmers, possibly because that large shape deformations of VE models easily lead to negative volume of the spheres.

C Comparison between Optimal and Square Strokes
Here we numerically compare the optimal stokes γ for NG 3-sphere, PMPY 2-sphere and VE 3-sphere models as obtained in Section 4.2 with stepwise square strokes γ sq . The energy dissipation and efficiency are shown in Table 7, the γ-paths and the X(γ(t)), t ∈ [0, 1] trajectories are given in Figure 7. To make a fair comparison, we take ε = 0.2 and X(γ) = 0.001 in all simulations. Data in Table 7 shows the improvement of swimming efficiency in the optimal strokes in all three model swimmers comparing to the stepwise square strokes. Table 7: Energy dissipation E(γ) and efficiency Eff of the optimal strokes γ and square strokes γ sq . All simulations are taken with X = 0.001 and ε = 0.2. The γ-path of the optimal stroke γ of a NG swimmer with 5 and 10 linked arms, respectively, where γ = (l 1 , l 2 , · · · , l m ) T . (f,g) The γ-path of the optimal stroke γ of a PMPY swimmer with 5 and 10 linked arms, respectively, where γ = (l 1 , l 2 , · · · , l m , a 1 , a 2 , · · · , a m ) T , where the solid lines give the l i -trajectories and dashed lines the a i -trajectories. ε = 0.2 and X = 5 × 10 −4 in all simulations.     . (a,c,e) give the γ-paths of the optimal strokes (γ) and the square strokes (γ sq ) of the NG, PMPY and VE swimmers, all with X(γ) = 0.001 and ε = 0. 2. (b,d,f ) give the stroke trajectories X(γ(t)), t ∈ [0, 1] (dashed lines) corresponding to the γ-paths (solid lines).