P Bettiol

B Bonnard

J Rouot

Optimal strokes at low Reynolds number: a geometric and numerical study of Copepod and Purcell swimmers

In this article we provide a comparative geometric and numerical analysis of optimal strokes for two different rigid links swimmer models at low Reynolds number: the Copepod swimmer (a symmetric swimmer recently introduced by Takagi) and the long-standing three-link Purcell swimmer by Purcell. The design of strokes satisfying standard performance criteria leads to investigate optimal control problems which can be analyzed in the framework of sub-Riemannian geometry. In this context nilpotent approximations allow to compute strokes with small amplitudes, which in turn can be used numerically to obtain more general strokes. For the Copepod model a detailed analysis of both abnormal and normal strokes is also described. First and second order optimality conditions, combined with numerical analysis, allow us to detect optimal strokes for both the Copepod and the Purcell swimmers. C 1 -optimality is investigated using the concept of conjugate point. Direct and indirect numerical schemes are implemented in Bocop and HamPath software to perform numerical simulations, which are crucial to complete the theoretical study and evaluate the optimal solutions.

Introduction

Swimming models at low Reynolds number employing rigid links in the swimmer design have been introduced in the 1950s (cf. [START_REF] Purcell | Life at low Reynolds number[END_REF]), and have been widely studied since they are of great interest for applications (cf. [START_REF] Alouges | Optimal strokes for low Reynolds number swimmers: an example[END_REF][START_REF] Chambrion | Optimal Strokes for Driftless Swimmers: A General Geometric Approach[END_REF][START_REF] Takagi | Swimming with stiff legs at low Reynolds number[END_REF], and see also [START_REF] Lighthill | Note on the swimming of slender fish[END_REF] for earlier literature). The concept of optimal stroke in swimming modeling conveys the idea of periodic 'best deformation strategy' (minimizing the mechanical energy dissipated by the drag forces, for instance) performed by the swimmer body interacting with a fluid, and can be contextualized within the periodic optimal control framework (see for instance [START_REF] Alouges | Optimal strokes for low Reynolds number swimmers: an example[END_REF][START_REF] Chambrion | Optimal Strokes for Driftless Swimmers: A General Geometric Approach[END_REF] and the references therein).

The problem of finding periodic solutions minimizing a given functional is a classical one, and goes back to Hamiltonian mechanics. A number of significant approaches and different techniques are now available. For instance the standard problem of determining periodic solutions of Hamiltonian vector fields, introduced by Poincaré in relation with the N-body problem, has been studied using continuation and variational methods: the aim was to show the existence of one-parameter family of periodic trajectories emanating from an equilibrium point [23, Thm. 9.2.1], and to apply direct methods to compute periodic solutions, in relation with the class of homotopy associated with the topology induced by collisions. These are valuable contributions still widely employed in direct and indirect numerical schemes like in the Bocop and HamPath software [START_REF] Bonnans | Bocop -A collection of examples[END_REF][START_REF] Cots | Contrôle optimal géométrique : méthodes homotopiques et applications[END_REF], which have been used here for our numerical analysis.

In this paper we consider two different rigid-link swimming models: the Copepod swimmer (recently introduced in [START_REF] Takagi | Swimming with stiff legs at low Reynolds number[END_REF] to interpret the behaviour of particular zoo-plankton species) and the Purcell threelink swimmer [START_REF] Purcell | Life at low Reynolds number[END_REF]. Both dynamics are modeled using the Resistive force theory (see [START_REF] Gray | The propulsion of sea-urchin spermatozoa[END_REF], [START_REF] Childress | Mechanics of swimming and flying[END_REF]Chapter 5]). Assuming that the displacement is performed minimizing a quadratic cost (w.r.t. the controls), the optimal strokes can be determined employing techniques coming from sub-Riemannian (SR-) geometry. In this area a number of important concepts and useful results, which are relevant for our study, have been recently provided. For instance, strokes with small amplitude can be investigated invoking the concept of nilpotent approximation (see [START_REF] Bellaïche | The tangent space in sub-Riemannian geometry[END_REF] or [START_REF] Jean | Control of Nonholonomic Systems: from Sub-Riemannian Geometry to Motion Planning[END_REF] for an introduction on this construct, and cf. [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF][START_REF] Chakir | Small Subriemannian Balls on R 3[END_REF] for explicit computations of spheres with small radius). In the SR-geometry context important features are related to the role of normal and abnormal geodesics (cf. [START_REF] Montgomery | A tour of subriemannian geometries, their geodesics and applications[END_REF][START_REF] Bonnard | On the role of abnormal minimizers in sub-Riemannian geometry[END_REF]), the smoothness of the minimizers (cf. [START_REF] Hakavuori | Non-minimality of corners in sub-Riemannian geometry[END_REF]), and crucial issues arise in the analysis of conjugate points in relation with C 1 -optimality (for the fixed initial and final points problem) in both, normal and abnormal, cases (cf. [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF]Chapter 6], [START_REF] Agrachev | Strong minimality of abnormal geodesics for 2-distributions[END_REF]). First and second order optimality conditions combined with numerical simulations become essential tools in the study of optimal strokes, which, due to the particular periodic structure of shape variables, typically can be embedded in families of strokes having the same cost performance. The non-uniqueness of periodic minimizers in general prevents employing standard second-order sufficient conditions. However, refined second-order conditions covering general classes of problems with non-unique minimizers have been provided in literature (cf. [START_REF] Wang | Necessary and sufficient conditions for local optimality of a periodic process[END_REF] and [START_REF] Gavriel | Second order sufficient conditions for optimal control problems with nonunique minimizers: an abstract framework[END_REF]). The results obtained by [START_REF] Gavriel | Second order sufficient conditions for optimal control problems with nonunique minimizers: an abstract framework[END_REF] are particularly appropriate for the nature of swimming problems here considered, where only one part of the variables are periodic, whereas the 'displacement variables' usually do not enjoy this feature.

The starting point for our work on the Copepod swimmer is the paper [START_REF] Takagi | Swimming with stiff legs at low Reynolds number[END_REF] in which two types of strokes are exhibited invoking merely geometric arguments (which do not involve optimal control theory): sinusoidal oscillations with a phase delay between adjacent links producing a regular curve (stroke), and a sequence of two separate links movements followed by a simultaneous recovery links paddling in unison generating a 'triangular stroke'. The aim in [START_REF] Takagi | Swimming with stiff legs at low Reynolds number[END_REF] was to demonstrate the possibility of swimming for this model, and subsequently to compare the obtained strokes with respect to a (locomotion) efficiency, which corresponds to a physical criterion for swimming problems widely adopted in literature. Under this perspective the conclusion of [START_REF] Takagi | Swimming with stiff legs at low Reynolds number[END_REF] is that (regular) strokes provided by sinusoidal oscillations have a better performance than the triangular one. Introducing the concept of geometric efficiency (which is proportional to the notion of efficiency in [START_REF] Takagi | Swimming with stiff legs at low Reynolds number[END_REF]), our analysis clarifies this point in the context of optimal control theory. We show that the triangular stroke is the unique abnormal stroke and traces a set (a triangle on the shape plane) which is characterized by the so-called Martinet points in SR-geometry, exception made for the corners of the triangle. Though the abnormal (triangular) stroke provides a large displacement, its length is also very big due to the corners, and as a consequence its efficiency is low. We provide a classification of normal strokes, and numerical results show that only (regular) simple loops, laying in the interior of the domain circumscribed by the triangle, can be regarded as candidates for minimizers. Considering the problem in which we fix the strokes displacement, the necessary optimality conditions related to the notion of conjugate points are applicable in the normal case (see [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF]) and do select only simple loop (normal) strokes as possible optimizers. Particular symmetry properties of the problem allow us to complement these necessary conditions with sufficient conditions [START_REF] Gavriel | Second order sufficient conditions for optimal control problems with nonunique minimizers: an abstract framework[END_REF] taking into account the non-uniqueness of minimizers in our problem.

Swimming models (at low Reynolds number) with three rigid links are well-known and studied in literature and can be described in terms of equations written in explicit form [START_REF] Passov | Supplementary notes to: Dynamics of Purcell's three-link microswimmer with a passive elastic tail[END_REF]. The Purcell three-link swimmer model can be interpreted as a control system which involves functions with long and complicated expressions, and the nilpotent approximation in SR-geometry is particularly appropriate for investigating this model. The 2-distribution provided by the associated vector fields has a particular structure (with growth vector (2, 3, 5), cf. [START_REF] Bettiol | The three links Purcell swimmer and some geometric problems related to periodic optimal controls[END_REF]) called the flat Cartan case (see [START_REF] Cartan | Les systèmes de Pfaff a cinq variables et les équations aux derivées partielles du second ordre[END_REF][START_REF] Yu | Symmetries of flat rank two distributions and sub-Riemannian structures[END_REF] for details on this geometric structure). The nilpotent approximation analysis leads to consider three crucial cases: the degenerate, the rotating and the oscillating case. In a previous paper [START_REF] Bettiol | The three links Purcell swimmer and some geometric problems related to periodic optimal controls[END_REF] only the oscillating case has been explored in detail. Here we complete the picture showing the integrability of the associated normal extremals (for all three cases) in the class of elliptic functions. We obtain analytical expressions also for abnormal extremals, proving that they are straight lines in the shape variables plane. Numerical simulations using Bocop and HamPath software allow us to compute more general strokes and test their optimality owing to appropriate algorithms which determine possible conjugate points for both normal and abnormal cases. Combining nilpotent approximation in SR-geometry with numerical continuation methods proves to be an extremely useful approach to compute strokes of large amplitude for the Purcell swimmer starting from strokes with small amplitude. The typical outcome for the oscillating case is an 'eight' stroke for which a conjugate point does appear in the reference time period. On the other hand, we can find simple loop strokes, generated in the degenerate and the rotating case, without any conjugate points (in the reference time interval). This analysis supplemented again with symmetry results allows us to individuate a normal extremal which satisfies also the refined second order sufficient conditions test for optimality.

The paper is organized as follows. In Sec. 2, we briefly introduce the two swimming models, some basic concepts and results needed for our study. Sec. 3 presents the combination of our geometric and numeric analysis to determine optimal strokes of the Copepod swimmer. Sec. 4 is devoted to the investigation of optimal strokes for the three-link Purcell swimmer. (x, y) θ 2

θ 1 α x y l 0 l 1 l 2 Figure 1: (left) Copepod swimmer, (right) Purcell swimmer.
Copepod swimmer. It is formed joining two pairs of symmetric links together (see [START_REF] Takagi | Swimming with stiff legs at low Reynolds number[END_REF]). Each pair of symmetric links can be seen as a scallop. We shall assume that all links length is normalized to l = 1. This swimming model allows displacements along a straight line and the swimmer velocity at position x ∈ R is given by:

θ1 = u 1 , θ2 = u 2 , ẋ = θ1 sin(θ 1 ) + θ2 sin(θ 2 ) ∆(θ) , (1)
where θ = (θ 1 , θ 2 ) ∈ R 2 and ∆(θ) := 2 + sin 2 (θ 1 ) + sin 2 (θ 2 ). Write q = (θ 1 , θ 2 , x) ∈ R 3 , F 1 (q) := (1, 0, sin θ 1 /∆(θ)) and F 2 (q) := (0, 1, sin θ 2 /∆(θ)). Therefore, the control system associated with the Copepod model takes the form

q = u 1 F 1 (q) + u 2 F 2 (q) . (2)
The mechanical power for this model is given by a positive quadratic form, and taking this quantity as a performance criterion yields to minimize the quadratic cost:

T 0 L(q(t), u(t))dt (3)
in which L(q, u) := a(q)u 2 1 + 2b(q)u 1 u 2 + c(q)u 2 2 , with

a(q) := 1 3 - sin 2 θ 1 2∆(θ) , b(q) := -sin θ 1 sin θ 2 2∆(θ) , c(q) := 1 3 - sin 2 θ 2 2∆(θ) .
Purcell swimmer. We shall consider also the Purcell three-link swimmer represented in Fig. 1 (right): the point (x, y) ∈ R 2 provides the position of the center of the body (central stick) which has length l 0 and forms an angle α with the x-axis. The two external arms have length l 1 and l 2 , forming relative (w.r.t. the central stick) angles θ 1 and θ 2 respectively. This model is much more complicated than the

Copepod model and allows displacements on the plane. Normalizing the links length (l 0 = 2, l 1 = l 2 = 1), the system dynamics takes the form:

  ẋ ẏ α  = 1 G   cos(α) -sin(α) 0 sin(α) cos(α) 0 0 0 1     g 11 g 12 g 21 g 22 g 31 g 32   θ1 θ2 , θ = u = S(θ)τ, (4)
where τ represents the physical control torque, and g ij , G and S are functions depending only on (θ 1 , θ 2 ) which have long and involved expressions (cf. [START_REF] Passov | Supplementary notes to: Dynamics of Purcell's three-link microswimmer with a passive elastic tail[END_REF] for details). For the Purcell three-link swimmer, the expanded mechanical power is given by

T 0 τ u dt = T 0 u S -1 (θ)u dt. ( 5
)
This quantity can be considered as a cost function to minimize and can be written in the form of a quadratic form like (3). Denoting q = (θ 1 , θ 2 , x, y, α) ∈ R 5 the state variable, we can write the swimming control system (4) as

q = u 1 F 1 (q) + u 2 F 2 (q) . (6)
(The expressions of the control vector fields F 1 , F 2 and of the mechanical energy can be found in [START_REF] Rouot | Méthodes géométriques et numériques en contrôle optimal et applications au transfert orbital poussée faible et à la nage à faible nombre de Reynolds[END_REF]).

Definition 2.1. A stroke of period T consists of a periodic motion (trajectory) in the shape variables.

For both models the two angular variables θ = (θ 1 , θ 2 ) are called the shape variables, and the period is taken T = 2π.

Remark 2.2. Standard existence theorems (cf. [START_REF] Vinter | Optimal control, Systems & Control: Foundations & Applications[END_REF]Sec. 2.8]) guarantee that for the optimal control problems here considered minimizers always exist.

Discussion on state constraints. Notice that a formulation of the control systems, corresponding to the swimmer models, should take account also of natural (physical) state constraints for the angular variables (θ 1 , θ 2 ):

• Copepod case. One has: θ 1 , θ 2 ∈ [0, π], θ 1 ≤ θ 2 .
• Purcell case. The state constraints here depend on the assumptions about the length l 0 of the body, and l 1 , l 2 which are respectively the lengths of the leg and the arms. We shall perform our computations assuming l 0 = 2 and l 1 = l 2 = 1. Hence we have the amplitude bounds:

θ 1 , θ 2 ∈ [-π, π].
We will find candidates to be optimal strokes which remain (strictly) inside the interior of the state constraints for the variables (θ 1 , θ 2 ), therefore these may be considered somehow non-active.

Elements of sub-Riemannian geometry

The reference optimal control problem for the swimmer models here considered is written in the SRgeometry context as

q = 2 i=1 u i F i (q), min u(.) J(u(.)) := T 0 (u R(q)u) dt, ( 7
)
where R is a positive definite quadratic form. The set of admissible controls U is the set of bounded measurable mappings. The length of an admissible trajectory q(.) of ( 7) associated with a control u is l(q(.))

= T 0 (u (t)R(q(t))u(t)) 1/2 dt .
Invoking standard concepts of SR-geometry (cf. [START_REF] Jean | Control of Nonholonomic Systems: from Sub-Riemannian Geometry to Motion Planning[END_REF]), we recall the following facts. Let D be the distribution span{F 1 , F 2 }. Using a feedback transformation u = β(q)v, we can find a local orthonormal frame (G 1 , G 2 ) such that the integral cost function reduces to T 0 v v dt. Near a point fixed q 0 , one can choose the so-called privileged coordinates so that the distribution D can be approximated by a nilpotent distribution denoted D = span{ F1 , F2 }. In our case, it is a crucial tool to initialize the numerical algorithms to compute strokes of the Purcell swimmer. Similarly, one can choose a nilpotent orthonormal frame denoted { Ĝ1 , Ĝ2 } to approximate the SR-problem and which plays a central role also to estimate the volume of small SR-balls [START_REF] Jean | Control of Nonholonomic Systems: from Sub-Riemannian Geometry to Motion Planning[END_REF]. We underline the fact that the nilpotent approximation is particularly relevant for the Purcell case where the complexity of the (non-approximated) functions involved leads to mathematical expressions which cannot be always explicitly written, nor numerically handled in a satisfactory way. In this paper we shall adopt the following convention for the Lie bracket of two vectors fields

F 1 , F 2 [F 1 , F 2 ](q) = ∂F 1 ∂q (q)F 2 (q) - ∂F 2 ∂q (q)F 1 (q).
We recall that if H 1 , H 2 are respectively the Hamiltonian lifts of F 1 and F 2 : p•F 1 , p•F 2 . Writing z = (q, p) (where q and p are respectively the state and adjoint variables) the symplectic coordinates (and -→ H the Hamiltonian vector field associated with the Hamiltonian H), one defines the Poisson bracket as:

{H 1 , H 2 }(z) = dH 1 ( - → H 2 )(z) = p • [F 1 , F 2 ](q).

Maximum Principle and computations of geodesic equations

The Maximum Principle is used to compute the geodesic equations. Assuming that the couple (F 1 , F 2 ) provides an orthonormal frame, the pseudo-Hamiltonian (also called unmaximized Hamiltonian) takes the form:

H(z, u) = 2 i=1 u i H i (z) + λ 0 2 i=1 u 2 i
where H i (.) is the Hamiltonian lift of F i , for i = 1, 2, and λ 0 is a constant which can be normalized to λ 0 = -1/2 (normal case) or λ 0 = 0 (abnormal case). We recall that both cases can be approached using standard techniques involving the maximality (also called Weierstrass) condition, obtaining important information about a minimizer candidate:

• Normal case: from ∂H ∂u = 0, we obtain u i = H i and, plugging these expressions of u i into H leads to the (true) Hamiltonian in the normal case:

H n = 1/2 (H 2 1 + H 2 2
). The corresponding solutions which satisfy the Maximum Principle are called normal extremals and their projections on the qspace are called normal geodesics. A normal geodesic is called strict if it is not the projection of an abnormal extremal.

• Abnormal case: we have the conditions H i (z) = 0, i = 1, 2. The corresponding solutions of the Maximum Principle are called abnormal extremals and their projections are called abnormal geodesics.

The concept provided by the following definition will play a crucial role in our context (cf Prop. 3.3 below for a link with optimality conditions).

Definition 2.3. Assuming the arc-length parameterization H n = 1/2, the exponential mapping is written:

exp q0 : (t, p(0)) → Π(exp(t -→ H n (z(0))))
, where z(0) = (q 0 , p(0)) and Π is the projection: z → q. A conjugate time (normal case) is a time t c such that the function exp q0 is not of full rank at t c and the corresponding point is called a conjugate point along the geodesic with initial condition z(0). We denote t 1c the first conjugate point.

We also recall that given the left endpoint condition q(0) = q 0 , the conjugate locus C(q 0 ) is the set of first conjugate points considering all normal geodesics emanating from q 0 . The sphere S(q 0 , r) is formed by the set of points at SR-distance r from q 0 .

Concepts of optimal control theory adapted to the swimmer problem

The study of minimizers for the Copepod and the Purcell models has to take into account two important features: the presence of state constraints and the periodicity of strokes. We shall overcome the first difficulty considering the family of all possible (not necessarily constrained) trajectories, and showing that the candidates to be minimizers stay in the interior of the state constraints, and therefore the latter can be in fact considered non-active. As a consequence we shall not enter state constraints in system [START_REF] Bonnans | Bocop -A collection of examples[END_REF]. On the other hand, the boundary conditions related to strokes provide periodic properties on controls and trajectories, which will play a crucial role in the selection of minimizers. We shall also consider an additional efficiency term (related to a standard physical criterion for swimming problems) in the cost which allows to establish that abnormal extremals do not have the best performance and, then, are not good candidates to be minimizers. This requires to consider optimal control problems involving a terminal cost function and specific boundary conditions. The proposition below summarizes standard first order optimality necessary conditions (see for instance [START_REF] Vinter | Optimal control, Systems & Control: Foundations & Applications[END_REF]Thm. 6.2.1]). We introduce the extra state variable q 0 which allows to rewrite the integral term in the cost functional as a final cost term, obtaining an extended state variable q = (q, q 0 ). This formulation is flexible enough to investigate the swimmers models considered here coupled with different cost functionals (such as mechanical energy, simplified SR energy, etc.).

Proposition 2.4. Write q = (q, q 0 ). Consider a minimizer (q * , u * ) for the optimal control problem min J(q, u) := g(q(T )) subject to bounded measurable functions u :

[0, T ] → R 2 and arcs q ∈ W 1,1 ([0, T ]; R n+1 ) satisfying q = 2 i=1 u i F i (q), q0 = L(q, u) (q(0), q(T )) ∈ C (8)
where T is fixed and

C is a closed subset of R n+1 × R n+1 . Write the Hamiltonian H(q, p, u) = p•F (q, u) = p• 2 i=1 u i F i (q)+p 0 L(q, u), p = (p, p 0 ).
Then, there exist an adjoint arc p * (•) and λ ≥ 0 such that the following necessary conditions are satisfied:

(p(0), λ) = 0, q * = ∂H ∂ p , ṗ * = - ∂H ∂ q a.e.
H(q * , p * , u * ) = max u∈U H(q * , p * , u) a.e .

together with the transversality condition:

(p * (0), -p * (T )) ∈ λ∇ q(T ) g(q * (T )) + N C (q * (0), q * (T ))

where N C is the limiting normal cone. Moreover there exists r ∈ R such that H(q * (t), p * (t), u * (t)) = r a.e..

Transversality condition ( 9) above provides useful information on periodic components of the optimal control problem, in particular when these do not intervene in the final cost g(•). It is expressed in terms of the limiting normal cone of the end-point constraint C, to summarize in a unified way all the different boundary conditions which we shall consider in our paper (for the notion on normal cones we refer the reader for instance to [START_REF] Vinter | Optimal control, Systems & Control: Foundations & Applications[END_REF]Chapter 4]).

Boundary conditions associated with periodicity. We split the state variable q into (q , q ) where q represents the periodic part of q. Let p = (p , p ) be the associated splitting of the adjoint vector. Imposing periodic conditions (as end-point conditions) on q , q (0) = q (T ), and assuming that g(.) does not depend on the q -component, the Maximum Principle leads to the transversality condition:

p (0) = p (T ). ( 10
)
For the swimmers models here considered the periodic component q typically includes the shape variables θ. However we shall consider a class of problems for the Purcell swimmer in which also α conditions has to satisfy the periodic boundary conditions.

Definition 2.5. A normal (resp. abnormal) stroke is a stroke corresponding to a normal (resp. abnormal) extremal. A piecewise smooth abnormal stroke is a piecewise smooth stroke such that each smooth sub-arc corresponds to an abnormal arc.

Shooting equation. To define the shooting equation, one restricts the flow to normal extremals, solution of -→ H n , with the following boundary conditions associated with the state variables splitting:

• q (0) = q 0 , q (T ) = q T , where q 0 , q T are fixed (and includes the non-periodic end point),

• q (0) = q (T ), p (0) = p (T ).

In the framework of SR-geometry and in relation with the underlying fixed end-points we have the following two properties [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF][START_REF] Chakir | Small Subriemannian Balls on R 3[END_REF].

Property 2.6. The shooting mapping fails to be locally injective if there exists conjugate points.

Property 2.7. The shooting mapping is defined on the cylinder and it fails to be proper if there exists abnormal extremals.

We also introduce the following notion of efficiency (cf. [START_REF] Alouges | Optimal strokes for low Reynolds number swimmers: an example[END_REF][START_REF] Chambrion | Optimal Strokes for Driftless Swimmers: A General Geometric Approach[END_REF][START_REF] Takagi | Swimming with stiff legs at low Reynolds number[END_REF]).

Definition 2.8. The geometric efficiency E of a stroke γ is defined as

• Copepod swimmer: E = x(T )/l(γ), • Purcell swimmer: E = x(T ) 2 + y(T ) 2 /l(γ)
that is the ratio between the euclidean displacement along (part of ) the state variable and the sub-Riemannian length of the stroke.

Remark 2.9. For the Copepod swimmer, in the concept of efficiency, E is proportional to

E := x 2 (T )/E(γ)
, where E(γ) is the energy of the stroke γ. E has been already used in the literature in the case of periodic strokes (see [START_REF] Alouges | Optimal strokes for low Reynolds number swimmers: an example[END_REF][START_REF] Takagi | Swimming with stiff legs at low Reynolds number[END_REF], and cf. [START_REF] Bettiol | Sub-Riemannian geometry and swimming at low Reynolds number: the Copepod case[END_REF] for further discussions).

Criteria imposing efficiency maximization can be viewed in terms of final cost minimization providing additional information on the necessary optimality conditions.

Optimality conditions and Geometric efficiency. Consider the optimal control problem of the form ( 8) where g = -E (for either the Copepod or the Purcell swimmer). If A(q(0), T ) is the accessibility set from q(0) = (q(0), 0) at time T and if (q * , p * , u * ) is an optimal solution for (8), then the geometric optimality conditions become q * (T ) ∈ ∂A(q(0), T ) .

Suppose also that we have a free right end-point condition for the non-periodic component, then the transversality condition (9) yields p * (T ) = λ∇ q g(q * ) and p (0) = p (T ) .

The Copepod swimmer

We start recalling two types of strokes which have been obtained in [START_REF] Takagi | Swimming with stiff legs at low Reynolds number[END_REF] employing just geometric arguments (not invoking optimal control techniques) showing that the model generates displacements; subsequently these are compared with respect to a locomotion efficiency term. These strokes constitute two reference cases for our analysis of the Copepod model, contributing as a motivation for our study.

First case

The two couples of legs are paddling in sequence followed by a recovery stroke performed in unison. In this case, the first step is to steer θ 2 from 0 to π, followed by θ 1 ; and the unison sequence corresponds to a displacement from π to 0 keeping θ 1 = θ 2 . Note that it corresponds to a triangle stroke on the (θ 1 , θ 2 )-plane. Moreover θ 1 and θ 2 stay on the boundary of the domain

{(θ 1 , θ 2 , x) : 0 ≤ θ 1 ≤ θ 2 ≤ π}.

Second case

The two legs are assumed to oscillate sinusoidally with period 2π according to

θ 1 (t) = Φ 1 + a cos(t), θ 2 (t) = Φ 2 + a cos(t + k) ( 11
)
where a, Φ 1 , Φ 2 and k are parameters.

Abnormal geodesics for the Copepod swimmer

We recall (cf. [START_REF] Zhitomirskiȋ | Typical singularities of differential 1-forms and Pfaffian equations[END_REF]) some basic facts about two-dimensional distributions in R 3 (D = span{F 1 , F 2 } associated with the reference control system in our case), which provide useful information on abnormal curves.

Local classification of two-dimensional distributions in dimension three and abnormal extremals. Denoting H i (z) = p • F i (q), for i = 1, 2, the Hamiltonian lifts of the vector field F i , the abnormal extremals are characterized by

H 1 (z) = H 2 (z) = 0.
Differentiating with respect to the time variable and using the condition,

ż = u 1 - → H 1 (z) + u 2 - → H 2 (z),
we obtain the relations

{H 1 , H 2 }(z) = 0 u 1 {{H 1 , H 2 }, H 1 }(z) + u 2 {{H 1 , H 2 }, H 2 }(z) = 0,
which allow to characterize the corresponding abnormal controls. Next, we introduce only the two stable models which are relevant for our study.

Contact case. A reference point q 0 ∈ R 3 is a contact point if span{F 1 , F 2 , [F 1 , F 2 ]
} is of dimension 3 at q 0 . At a contact point, identified with 0, there exists a system of local coordinates q = (x, y, z) such that D = ker(α) where α := dz + (xdy -ydx).

The corresponding nilpotent frame is

F1 = ∂ ∂x + y ∂ ∂z , F2 = ∂ ∂y -x ∂ ∂z .
Taking F1 , F2 as an orthonormal frame, we obtain the Heisenberg model in SR-geometry. Observe that dα = -2 dy ∧ dx is the Darboux form and that ∂ ∂z is the characteristic direction of dα.

The Martinet case. A point q 0 is a Martinet point if at q 0 we have the following property:

[F 1 , F 2 ] ∈ D = span{F 1 , F 2 } but { [[F 1 , F 2 ], F 1 ], [[F 1 , F 2 ], F 2 ] } D.
Then, there exist local coordinates q = (x, y, z) near q 0 (which can be identified as 0 in these new coordinates) such that D = ker ω, where ω := dz -y 2 / 2dx.

The corresponding nilpotent frame takes the form

F1 = ∂ ∂x + y 2 2 ∂ ∂z , F2 = ∂ ∂y .
Moreover, we obtain

F3 = [ F1 , F2 ] = y ∂ ∂z , [[ F1 , F2 ], F1 ] = 0, [[ F1 , F2 ], F2 ] = ∂ ∂z . (12)
Locally, at a Martinet point, the surface of equation y = 0, where F1 , F2 , [ F1 , F2 ] are linearly dependent, is called the Martinet surface and is foliated by abnormal curves, solutions of the vector field ∂ ∂x . In particular, through the origin it corresponds to the curve t → (t, 0, 0). Taking F1 , F2 as an orthonormal frame, it corresponds to the so-called flat Martinet case.

We apply now these constructs to the Copepod model. Recalling that

F i = sin(θ i ) ∆(θ) ∂ ∂x + ∂ ∂θ i , i = 1, 2 , and ∆(θ) := 2 + sin 2 (θ 1 ) + sin 2 (θ 2 ), (13)
we obtain

F 3 := [F 1 , F 2 ] = f (θ 1 , θ 2 ) ∂ ∂x and [[F 1 , F 2 ], F i ] = ∂f ∂θ i (θ 1 , θ 2 ) ∂ ∂x , for i = 1, 2 , (14)
where

f (θ 1 , θ 2 ) := 2 sin(θ 1 ) sin(θ 2 )(cos(θ 1 ) -cos(θ 2 )) / ∆ 2 (θ). ( 15
)
Consider the subset of R 3

Σ := {(θ 1 , θ 2 , x) | sin θ 1 sin θ 2 (cos θ 1 -cos θ 2 ) = 0 and 0 ≤ θ 1 ≤ θ 2 ≤ π}.
Observe that Σ represents the boundary of the physical domain (state constraints)

{(θ 1 , θ 2 , x) | 0 ≤ θ 1 ≤ θ 2
≤ π} which has on the (θ 1 , θ 2 )-plane the shape of a triangle having edges

Σ 1 := {(θ 1 , θ 2 , x) ∈ Σ | θ 1 = 0, 0 < θ 2 < π}, Σ 2 := {(θ 1 , θ 2 , x) ∈ Σ | θ 2 = π, 0 < θ 1 < π}, Σ 3 := {(θ 1 , θ 2 , x) ∈ Σ | 0 < θ 1 = θ 2 < π}.
From ( 13)-( 15) we easily deduce the following result.

Proposition 3.1.

i) The vector fields F 1 , F 2 and [F 1 , F 2 ] are linearly dependent on the set Σ which contains no contact points. Each point of

Σ 1 ∪ Σ 2 ∪ Σ 3 is a Martinet point.
ii) The triangle stroke is abnormal and the associated adjoint vector can be normalized to

p = (1, sin(θ 1 ) / ∆(θ), sin(θ 2 ) / ∆(θ)) .
Remark 3.2. The previous proposition provides the interpretation of the policy described in the 'First case' of Sec. 3, which corresponds exactly to the abnormal stroke. Notice that it provides in the (θ 1 , θ 2 )-plane the boundary of the physical domain for the Copepod model. A recent contribution [START_REF] Hakavuori | Non-minimality of corners in sub-Riemannian geometry[END_REF], applying SR-geometry arguments, proves that such an abnormal curve with corners cannot be optimal (in [START_REF] Hakavuori | Non-minimality of corners in sub-Riemannian geometry[END_REF] the authors consider problems in absence of state constraints). Introducing the efficiency term (see Def. 2.8) allow us to exclude the triangle stroke from the optimal trajectories (cf. the discussion at the end of Sec. 3).

The normal case

In the previous section we have discussed the abnormal case, which provide strokes having necessarily a 'triangular' shape. The 'Second case' (mentioned at the beginning of Sec. 3) suggests to investigate also strokes which are not abnormal and can be possibly described in terms of regular functions. This requires dealing with the class of normal extremals. All the numerical simulations are performed using the mechanical energy (in fact the results are not senstive to the choice of the cost, see [START_REF] Bettiol | Sub-Riemannian geometry and swimming at low Reynolds number: the Copepod case[END_REF]). Nevertheless, in order to interpret Takagi's result in the next subsection, we have to consider the simplified cost.

Geometric structure for the simplified cost and relation with Takagi's sinusoidal

paddling [START_REF] Bonnard | On the role of abnormal minimizers in sub-Riemannian geometry[END_REF] Note that if the cost is simplified to

T 0 (u 2 1 +u 2
2 )dt, some geometric computations can be made, in relation with the Heisenberg case (assuming F 1 , F 2 orthonormal) and can be used in the numerical implementation, in particular to determine strokes with small amplitudes. In this case, we have

H n = 1/2 H 2 1 + H 2 2 .
Using the Poincaré coordinates (q, H), H = (H 1 , H 2 , H 3 ) and H i = p•F i (q), for i = 1, 2, 3, straightforward computations (evaluated in the triangular domain) yield:

Ḣ1 = dH 1 ( H n ) = {H 1 , H 2 }H 2 = H 2 H 3 , Ḣ2 = dH 2 ( H n ) = {H 2 , H 1 }H 1 = -H 1 H 3 , Ḣ3 = dH 3 ( H n ) = {H 3 , H 1 }H 1 + {H 3 , H 2 }H 2 , with {H 3 , H 1 }(z) = p • [[F 1 , F 2 ], F 1 ](q) and {H 3 , H 2 }(z) = p • [[F 1 , F 2 ], F 2 ](q).
At a contact point, F 1 , F 2 , F 3 form a frame, therefore we obtain

[[F 1 , F 2 ] , F 1 ] (q) = 3 i=1 λ i (q)F i (q), [[F 1 , F 2 ] , F 2 ] (q) = 3 i=1 λ i (q)F i (q)
where

λ 1 = λ 2 = 0, ∂f ∂θ1 = λ 3 f, λ 1 = λ 2 = 0, ∂f ∂θ2 = λ 3 f . We conclude that Ḣ1 = H 2 H 3 , Ḣ2 = -H 1 H 3 , Ḣ3 = H 3 (λ 3 H 1 + λ 3 H 2 ) . (16)
The associated one dimensional distribution can be analyzed setting ds = H 3 dt and we obtain

dH 1 ds = H 2 , dH 2 ds = -H 1 , dH 3 ds = λ 3 H 1 + λ 3 H 2 . ( 17
)
In particular, differentiating one more time the first relation in ( 17) with respect to s and using the second relation in [START_REF] Gavriel | Second order sufficient conditions for optimal control problems with nonunique minimizers: an abstract framework[END_REF], we have the harmonic oscillator H 1 + H 1 = 0. Furthermore H 3 can be analyzed using the remaining equation ( 16). Observe that with the approximation λ 3 , λ 3 constant, the equation takes the form dH 3 ds = A cos(s + ρ).

with A, ρ constant. In these computations, we recognize the Heisenberg case, corresponding to λ 3 = λ 3 = 0. Observe that when q is not a contact point (that is

F 2 = [F 1 , F 2 ] ∈ span{F 1 , F 2 })
, in order to deal with the Martinet case, we can choose the frame F 1 , F 2 and F 3 , where

F 1 = F 1 , F 2 = F 2 and F 3 = ∂ ∂x .

A strokes classification and numerical results (mechanical cost)

From the micro-local point of view, there exists a variety of different types of normal strokes for the Copepod (e.g. simple loops, 'eights', 'limaçons') in relation with the classification of planar periodic curves [START_REF] Berger | La taxonomie des courbes[END_REF]. Our aim in this section is first to establish numerically which kind of normal strokes can be realized by the Copepod swimmer. Subsequently we compute the conjugate points along a strict normal stroke to select minimizer candidates for the reference problem (with given end-point conditions). We underline the fact that at this stage we are not taking into account the state constraints, for we want to classify normal strokes among all possible normal geodesics, and provide (numerically) a selection of candidates to satisfy optimality criteria in terms of this classification. The exploration is purely numerical and does not pretend to be exhaustive (we refer to [START_REF] Bettiol | Sub-Riemannian geometry and swimming at low Reynolds number: the Copepod case[END_REF] for a complementary analysis, in particular, it can be proved that the only strokes with small amplitudes in the interior of the triangle are simple loops and limaçons).

Numerical methods. The period T is fixed at 2π in our simulations. We use the HamPath software [START_REF] Cots | Contrôle optimal géométrique : méthodes homotopiques et applications[END_REF] at two levels:

1. The shooting equations associated with the problem are

x(0) = 0, x(2π) = x f , θ j (0) = θ j (2π) j = 1, 2, p k (0) = p k (2π) k = 2, 3. (18)
2. We consider a normal stroke and we test its optimality by showing the non-existence of conjugate points using the variational equation to compute Jacobi fields. Recall (cf. Def. 2.3 and [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF]) that given a reference curve (q(t), p(t)) solution of the vector field -→ H n , a time t c ∈]0, 2π] is a conjugate time if there exists a Jacobi field δz = (δq, δp), that is a non-zero solution of the variational equation

δz(t) = ∂ H n ∂z (q(t), p(t))δz(t) (19)
such that δq(0) = δq(t c ) = 0. We denote δz i = (δq i , δp i ), i = 1, ..., 3, three linearly independent solutions of ( 19) with initial condition δq(0) = 0. At time t c we have the following rank condition

rank{δq 1 (t c ), δq 2 (t c ), δq 3 (t c )} < 3. ( 20
)
The following proposition is an immediate consequence of well-known results establishing a fundamental relationship between weak minimizers and existence of conjugate points (see [START_REF] Agrachev | Strong minimality of abnormal geodesics for 2-distributions[END_REF][START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF][START_REF] Wang | Necessary and sufficient conditions for local optimality of a periodic process[END_REF]).

Proposition 3.3. A necessary optimality condition for a strict normal stroke to provide a weak minimizer is the non-existence of a conjugate point on ]0, 2π[.

Complexity of optimal policies and selection of normal strokes. Fig. 2 illustrates the behaviour of four different types of strokes (we are not taking into account the state constraints at this stage) highlighting the complexity of the model. These can be regarded as examples covering the generic classification of periodic planar curves [START_REF] Berger | La taxonomie des courbes[END_REF]. Conjugate points are also computed to check the second order optimality conditions (cf. Prop. 3.3 above). There is no conjugate point on ]0, 2π] just in the case of the simple loop, but they do appear for the limaçon case, the eight case. Further simulations, taking into account more complicated shapes (combining two eights for instance), confirm the presence of conjugate points on ]0, 2π[. From these numerical simulations we identify only one class of curves in which the conjugate point is absent: the simple loops. We deduce that this class provides good candidates for optimality. Continuation method. We compute numerically simple loops from the stroke computed by Bocop of Fig. 3 and the latter is used to initialize the continuation method with HamPath to generate a oneparameter family of simple loops. We compare the efficiency of these strokes to deduce the most efficient stroke in this family (see Fig. 5). (We refer the reader to [START_REF] Bettiol | Sub-Riemannian geometry and swimming at low Reynolds number: the Copepod case[END_REF] for a microlocal and SR-normal forms study, which allows continuation method starting from a stroke with small amplitude).

Optimal curves confined in the triangular domain (state constraint)

In this section we enter the state constraints, in the sense that, owing to numerical analysis, we are able to select the curves (strokes) which are good candidates to be minimizers for our reference problem. More precisely, numerical simulations provide evidence that only simple loops confined in the interior of the state constraints (the domain having a triangular shape on the (θ 1 , θ 2 )-variables) can be considered as minimizers. In our study, we use a combination of the Bocop and HamPath software.

• Bocop : This software implements a direct method which can take into account state constraints and is used to initialize the shooting algorithm of the HamPath software (see [START_REF] Bonnans | Bocop -A collection of examples[END_REF]).

• HamPath : This software presented in [START_REF] Cots | Contrôle optimal géométrique : méthodes homotopiques et applications[END_REF] implements an indirect method based on the Pontryagin Maximum Principle. A shooting algorithm and continuation method are used to solve the boundary value problem. Note that it cannot be directly applied to compute the optimal solution with state constraints, due to the complexity of the implementation of the maximum principle with state constraints [START_REF] Bonnard | Mécanique céleste et contrôle des véhicules spatiaux. Mathématiques & Applications[END_REF]. Fixing the energy level H n = 1/2, the domain of the exponential map is not compact (it is a cylinder). It turns out that the shooting problem, consisting in finding an initial adjoint vector, is ill-conditioned since the exponential map is non proper (|p| → +∞) as illustrated in Fig. 4 when we compute normal extremals near the abnormal extremal (the triangle stroke).

Second order sufficient conditions for the Copepod strokes. The aim of this section is to check second order sufficient conditions for normal extremals represented by simple loop strokes of the Copepod swimmer. We employ here particular second order sufficient conditions which can be used in presence of non-uniqueness of minimizers [START_REF] Gavriel | Second order sufficient conditions for optimal control problems with nonunique minimizers: an abstract framework[END_REF]. We consider the optimal control problem in which we minimize the mechanical cost (3) over trajectories satisfying [START_REF] Alouges | Optimal strokes for low Reynolds number swimmers: an example[END_REF] and such that

x(0) = 0, x(T ) = x T , θ j (0) = θ j (T ) j = 1, 2. ( 21
)
Proposition 3.4. Take I = (-ε, ε) for some ε > 0 small enough and let (q(.), p(.), ū(.)) be a normal extremal (on [0, T ]) where q = ( θ1 , θ2 , x), p = (p 1 , p 2 , p 3 ) and ū = (ū 1 , ū2 ). Write θj (.), j = 1, 2, p(•), x(.) and ū(•) their corresponding periodic extensions. For all a ∈ I and t ∈ [0, T ], we define

q a (•) = (θ a 1 (•), θ a 2 (•), x a (•)), u a 1 (•), u a 2 (•) and p a (•) by θ a j (t) = θj (t + a), u a j (t) = ūj (t + a) for j = 1, 2, x a (t) = x(t + a) -x(a), p a (t) = (p 1 (t), p 2 (t + a), p 3 (t + a)).
Then, for ε > 0 small enough, the normal extremal (q(.), p(.), ū(.)) is continuously embedded in the family of extremals (q a (.), p a (.), u a (.)) a∈I .

In what follows, we denote by (q, ū) the simple loop of Fig. 3 (left) satisfying ( 21) with the associated adjoint vector p.

Numerical result. The simple loop (q, ū) is weak-locally optimal.

To confirm this claim, we shall invoke the second order test provided by [START_REF] Gavriel | Second order sufficient conditions for optimal control problems with nonunique minimizers: an abstract framework[END_REF]Thm. 4.2]. Full details are provided in [START_REF] Rouot | Méthodes géométriques et numériques en contrôle optimal et applications au transfert orbital poussée faible et à la nage à faible nombre de Reynolds[END_REF]. The data of the control systems ( 2)-(3) and Prop. 3.4 ensure that all the assumptions required to apply [17, Thm. 4.2] are satisfied. Numerical arguments allow to conclude that there is no conjugate point on [0, T ] for the normal extremal (q(.), p, ū(.)) (cf. Sec. 3.2.2 and Fig. 3 (left)). It implies that the Riccati equation associated with the accessory problem has a global symmetric solution (cf. [START_REF] Gavriel | Second order sufficient conditions for optimal control problems with nonunique minimizers: an abstract framework[END_REF]Equation (7)].

Using the Isoda integrator from the FORTRAN library odepack, we compute numerically the matrix W , defined in [START_REF] Gavriel | Second order sufficient conditions for optimal control problems with nonunique minimizers: an abstract framework[END_REF]Equation (8)]. Consider the linear subspace

L s := { (y 0 , y T ) ∈ R 3 × R 3 | ∇ q0,q T c(q 0 , q T ) (y 0 y T ) = 0 } ,
where c(q 0 , q T ) = (x(0), x(T ), θ 1 (0) -θ 1 (T ), θ 2 (0) -θ 2 (T )) . We introduce the matrix M s such that L s = Im(M s ). Standard second order sufficient conditions of [START_REF] Gavriel | Second order sufficient conditions for optimal control problems with nonunique minimizers: an abstract framework[END_REF]Thm. 3.4] lead to check that W is definite-positive on L s . This is equivalent to check that W s = M s (W +W )M s is positive-definite. Due to the non-uniqueness of the extremal (q, p, ū) (cf. the symmetry highlighted by Prop. 3.4), W s is not definite-positive (see Table 1). For the refined second order sufficient conditions of [START_REF] Gavriel | Second order sufficient conditions for optimal control problems with nonunique minimizers: an abstract framework[END_REF]Thm. 4.2], we consider the vector Γ = (0, ẋ2 (0), ẋ3 (0), 0, ẋ2 (0), ẋ3 (0)) the linear subspace

L r = L s ∩ {(y 0 , y T ) ∈ R 3 × R 3 | Γ (y 0 y T ) = 0}
and the matrix M r such that L r = Im(M r ), hence M r = 0 -ẋ3 (0) ẋ2 (0) 0 -ẋ3 (0) ẋ2 (0) . Numerical simulations confirm that W r = M r (W +W )M r is positive-definite and taking different relative tolerances for the integrator, W s has zero eigenvalue associated with the vector Γ (see Table 1). Therefore the refined second order conditions are applicable confirming the (local) optimality of the reference simple loop (q, ū). Table 1: The standard second order sufficient conditions fail: W s has zero eigenvalue and the refined second order sufficient conditions are satisfied: W r is positive-definite.

Comparison of strokes geometric efficiency. To compare normal and abnormal solutions corresponding to different displacements, in Fig. 5 we represent the efficiency ratio E = x(T )/l(γ) concerning solutions obtained for a given displacement x(T ) and l(γ) is the length of the stroke γ (this quantity does not depend upon the parameterization). For the triangle, a displacement along the vertical or horizontal edge gives x = 2

√

q • qdt and is given by 2π √ 2H n where H n is the energy level. The efficiency curve is displayed in Fig. 5 where the normal strokes corresponding to the maximal efficiency is also represented. Observe that from computations above obtained for the abnormal stroke, the efficiency turns out to be 5.56e-2. This result has been compared with the efficiency of normal strokes (see Fig. 5) establishing the optimality of normal strokes (in terms of efficiency).

Conclusions. From our analysis we deduce that the (triangle) abnormal stroke is not optimal. Indeed, one can choose a normal stroke (inside the triangle) such that the displacement is x/2 with x = x ab (T ) 2.742.10 -1 and the length is less than l ab /2 where l ab is the length of the triangle. Applying twice the normal stroke, we obtain the same displacement x than with the abnormal stroke but with a length strictly lower than l. The contribution [START_REF] Hakavuori | Non-minimality of corners in sub-Riemannian geometry[END_REF] proves that such abnormal trajectory with a corner is not C 0 -optimal in the unconstrained case. We were not able to prove the same result for our abnormal stroke taking into account the triangle constraint.

The three-Link Purcell swimmer

We restrict attention now to the optimal control problem associated with the Purcell swimmer model, which has been presented in Sec. 2.

A symmetry property

First of all due to the structure of equations ( 4) we have the following result. Lemma 4.1. Let q(t) = ( θ1 (t), θ2 (t), x(t), ȳ(t), ᾱ(t)) and q(t) = (θ 1 (t), θ 2 (t), x(t), y(t), α(t)) be the solutions associated with u(.) with respective initial conditions (θ 10 , θ 20 , 0, 0, 0) and (θ 10 , θ 20 , 0, 0, α 0 ). Then θ j (t) = θj (t) j = 1, 2, α(t) = ᾱ(t) + α 0 , x(t) = cos(α 0 )x(t) -sin(α 0 )ȳ(t), y(t) = sin(α 0 )x(t) + cos(α 0 )ȳ(t).

(

) 22
Straightforward computations yield the following results.

Proposition 4.2. Let H n be the extremal normal Hamiltonian associated with any quadratic cost 2π 0 (a(q)u 2 1 + 2b(q)u 1 u 2 + c(q)u 2 2 )dt. Denoting (p θ , p x , p y , p α ) the adjoint components, we have the following first integrals

I 1 = p x , I 2 = p y , I 3 = H n , I 4 = (p x y -p y x) -p α Corollary 4.3. The shooting conditions • θ, α 2π-periodic, • x(0) = y(0) = 0 and (x 2 + y 2 )(2π) = r 2 , can be completed considering • p θ , p α 2π-periodic, • (p x y -xp y )(2π) = 0 i.e. (p x , p y )(2π) is normal to S(r) : x 2 + y 2 = r 2 ).
Observe that, since I 4 = (p x y -xp y )-p α is a first integral, and at t = 0 and t = 2π we have p x y -p y x = 0, we deduce that p α (0) = p α (2π). Hence the assertion above 'p α is 2π-periodic' is equivalent to p is normal to S(r) and one of the conditions can be relaxed and be replaced by α(0) = 0 to determine the solution.

Nilpotent approximation

Due to the mathematical complexity of the expressions of the vector fields F 1 and F 2 employed in this model, the nilpotent approximation will play a crucial role in our analysis. First of all, owing to the integrability of the associated normal extremals in the class of elliptic functions, the nilpotent approximation will allow to make a micro-local analysis of the different kinds of strokes and to establish the existence of conjugate points using a suitable time rescaling. Second, the abnormal extremals forming piecewise smooth strokes can be easily computed in this approximation and, then, the optimality of these strokes can be studied using the concept of the (corresponding) conjugate point.

The flat nilpotent model

The Purcell system (4) can be written as a control system of the form q = F u = 2 i=1 u i F i , where q = (θ 1 , θ 2 , x, y, α) ∈ R 5 . Though the vectors fields F 1 , F 2 have a complicated expression, they provide a 2-distribution with growth (2, 3, 5) (cf. [START_REF] Bettiol | The three links Purcell swimmer and some geometric problems related to periodic optimal controls[END_REF]). There exists a unique nilpotent model associated with a 2-dimensional distribution in dimension 5 with growth vector [START_REF] Alouges | Optimal strokes for low Reynolds number swimmers: an example[END_REF][START_REF] Bellaïche | The tangent space in sub-Riemannian geometry[END_REF][START_REF] Bettiol | The three links Purcell swimmer and some geometric problems related to periodic optimal controls[END_REF], see [START_REF] Cartan | Les systèmes de Pfaff a cinq variables et les équations aux derivées partielles du second ordre[END_REF][START_REF] Yu | Symmetries of flat rank two distributions and sub-Riemannian structures[END_REF]. Definition 4.4. We call the flat Cartan model the 2-dimensional distribution in dimension five defined by the two vector fields:

F1 (q) = ∂ ∂ q1 , F2 (q) = ∂ ∂ q2 + q1 ∂ ∂ q3 + q3 ∂ ∂ q4 + q2 1 ∂ ∂ q5 . ( 23
)
where q = (q 1 , q2 , q3 , q4 , q5 ) are the privileged coordinates with the following weights: 1 for q1 and q2 , 2 for q3 , and 3 for q4 and q5 .

Computation of the nilpotent approximation

The nilpotent approximation of the Purcell model is computed at q = 0. It provides a nilpotent approximation for the SR-problem with the simplified cost 2π 0 (u 2 1 + u 2 2 )dt. The two-jets of F 1 and F 2 at q = 0 are

F 1 (q) = ∂ ∂q 1 + - 1 6 q 5 - 4 27 q 1 - 2 27 q 2 ∂ ∂q 3 + 1 6 - 1 12 q 5 2 - 2 27 q 5 q 2 - 4 27 q 5 q 1 - 1 27 q 1 2 - 1 27 q 1 q 2 - 1 36 q 2 2 ∂ ∂q 4 + - 7 27 + 2 81 q 1 2 - 2 81 q 1 q 2 - 5 162 q 2 2 ∂ ∂q 5 + O(|q| 3 ). F 2 (q) = ∂ ∂q 2 + 1 6 q 5 + 4 27 q 2 + 2 27 q 1 ∂ ∂q 3 + - 1 6 + 1 12 q 5 2 + 4 27 q 5 q 2 + 2 27 q 5 q 1 + 1 36 q 1 2 + 1 27 q 1 q 2 + 1 27 q 2 2 ∂ ∂q 4 + - 7 27 - 5 162 q 1 2 - 2 81 q 1 q 2 + 2 81 q 2 2 ∂ ∂q 5 + O(|q| 3 ) .
The local diffeomorphism ϕ, which transforms F 1 , F 2 into the nilpotent approximation F1 , F2 , can be explicitly written using a sequence ϕ = ϕ N o ... o ϕ 1 : R 5 → R 5 , where N = 13 (cf. [START_REF] Bettiol | The three links Purcell swimmer and some geometric problems related to periodic optimal controls[END_REF]). This leads to a complicated transformation whose role is to relate the privileged coordinates to the physical coordinates (θ 1 , θ 2 , x, y, α) in particular we have a 'stability' property for the shape variables:

Proposition 4.5. The shape variables θ = (θ 1 , θ 2 ) corresponds to the (q 1 , q2 ) coordinates.

Integration of normal extremal trajectories

Computing with ( 23), one has

F1 (q) = ∂ ∂ q1 , (q) = ∂ ∂ q2 + q1 ∂ ∂ q3 + q3 ∂ ∂ q4 + q2 1 ∂ ∂ q5 , [ F1 , F2 ](q) = - ∂ ∂ q3 -2q 1 ∂ ∂ q5 , [[ F1 , F2 ], F1 ](q) = -2 ∂ ∂ q5 , [[ F1 , F2 ], F2 ](q) = ∂ ∂ q4 .
All brackets of length greater than 3 are zero. Write ẑ = (q, p). Employing the corresponding Hamiltonian lifts, one has:

H 1 (ẑ) = p • F1 (q) = p1 , H 2 (ẑ) = p • F2 (q) = p2 + p3 q1 + p4 q3 + p5 q2 1 , H 3 (ẑ) = p • [ F1 , F2 ](q) = -p 3 -2q 1 p5 , H 4 (ẑ) = p • [[ F1 , F2 ], F1 ](q) = -2p 5 , H 5 (ẑ) = p • [[ F1 , F2 ], F2 ](q) = p4 . The SR-Cartan flat case is q = 2 i=1 u i Fi , min u 2π 0 (u 2 1 + u 2 2 )dt.
and the normal Hamiltonian takes the form

H n = 1/2 (H 2 1 + H 2 2 ). (24)
More precisely, using the Poincaré coordinates, the control system is written

q1 = H 1 , q2 = H 2 , q3 = H 2 q1 , q4 = H 2 q3 , q5 = H 2 q2 1 . (25)
Deriving with respect to the time variable, we have

Ḣ1 = dH 1 ( H n ) = {H 1 , H 2 }H 2 = p • [ F1 , F2 ](q)H 2 = H 2 H 3 , Ḣ2 = -H 3 H 1 , Ḣ3 = H 1 H 4 + H 2 H 5 , Ḣ4 = 0 hence H 4 = c 4 , Ḣ5 = 0 hence H 5 = c 5 .
Fixing the level energy,

H 2 1 + H 2 2 = 1 we set H 1 = cos(ϑ) and H 2 = sin(ϑ). Ḣ1 = -sin(ϑ) θ = H 2 H 3 = sin(ϑ)H 3 . Hence θ = -H 3 and θ = -(H 1 c 4 + H 2 c 5 ) = -c 4 cos(ϑ) -c 5 sin(ϑ) = -ω 2 sin(ϑ + φ) (26)
where ω and φ are constant. More precisely, we have

ω = p4 (0) 2 + 4p 5 (0) 2 1/4 , φ = arctan (-2 p5 (0)/p 4 (0)) .
• First, we consider the degenerate case which corresponds to θ = 0. Therefore, ϑ(t) = p30 t + ϑ 0 where ϑ 0 , p30 are constant and for p30 = 0, the solutions of (25), q1 , q2 , are expressed as

q1 (t) = q10 + 1/p 3 (0) sin(p 3 (0) t + ϑ 0 ), q2 (t) = q20 -1/p 3 (0) cos(p 3 (0) t + ϑ 0 ). ( 27
)
where q10 , q20 are constant.

• Second, the case where θ = 0 leads to a pendulum equation. Indeed, taking ψ = ϑ+φ, (26) becomes

1/2 ψ2 -ω 2 cos(ψ) = B, (28)
where B is the constant

B = 1/2 (p 3 (0) + 2 q1 (0)p 5 (0)) 2 -p1 (0) p4 (0) -2 p5 (0) p2 (0) -2 p5 (0) p4 (0) x3 (0).
We have the following two possible cases.

-Oscillating case. We introduce

k 2 = 1/2 + B/(2 ω 2 ) with 0 < k < 1 so that (28) becomes ψ2 = 4ω 2 k 2 -sin 2 (ψ/2)
and, using standard relations on elliptic functions (cf. [21, Chapter 5]), we obtain sin(ψ/2) = k sn(u, k), cos(ψ/2) = dn(u, k) where u = ωt + ϕ 0 . H 1 and H 2 are elliptic functions of the first kind and the solutions of (25), q1 , q2 , are expressed as

ω q1 (u) = ω q10 + -2 k sin (φ) cn (u ) + (-u + 2 E (u )) cos (φ) ω q2 (u) = ω q20 + -2 k cos (φ) cn (u ) + (u -2 E (u )) sin (φ) (29)
where q10 and q20 are constant, cn, sn, dn are elliptic functions of the first kind and E(.) is the elliptic integral of the second kind.

-Rotating case. We introduce k 2 = 2 ω 2 /(B + ω 2 ) with 0 < k < 1 so that (28) becomes

ψ2 = 4ω 2 /k 2 1 -k 2 sin 2 (ψ/2) .
q(.), there must hold H 2 (q(.), p(.)) = 0 and, differentiating with respect to t, {H 1 , H 2 }(q(.), p(.)) = 0 must hold too. Deriving one more time, the extremals associated with the controls

u a (q, p) = {H 1 , {H 2 , H 1 }}(q, p) / {H 2 , {H 1 , H 2 }}(q, p) = 2 p5 / p4
satisfy the relation H 2 = {H 1 , H 2 } = 0 along (q(.), p(.)) and are solutions of q = ∂Ha ∂ p , ṗ = -∂Ha ∂ q where H a is the true Hamiltonian H a (q, p) = H 1 (q, p) + u a H 2 (q, p) = p1 + 2 p5 p2 + p3 q1 + p4 q3 + p5 q2 1 /p 4 .

From the Pontryagin maximum principle, we also have H 1 (q(.), p(.)) = 0. The extremal system subject to the constraints

H 1 = H 2 = {H 1 , H 2 } = 0 is
integrable, in particular the abnormals can be written as polynomial functions (see [START_REF] Rouot | Méthodes géométriques et numériques en contrôle optimal et applications au transfert orbital poussée faible et à la nage à faible nombre de Reynolds[END_REF]).

Lemma 4.6. The θ-projection of abnormals are straight lines and form triangular strokes:

q1 (t) = t + q1 (0), q2 (t) = 2 p5 (0)/p 4 (0)t + q2 (0)
where (q(0), p(0)) are constant satisfying the constraints

H 1 = H 2 = {H 1 , H 2 } = 0.

Numerical results

Computations of conjugate points

Normal case. In the normal case, we consider the extremal system given by the normal Hamiltonian [START_REF] Montgomery | A tour of subriemannian geometries, their geodesics and applications[END_REF]. In Sec. 4.2.3, we described three types of extremals. For each case, we have computed solutions using HamPath , representing the control, the state and the adjoint variables as functions of time (see Fig. 6, Fig. 7, Fig. 8). We also illustrate the conjugate points evaluated according to the algorithm [START_REF] Bonnard | Mécanique céleste et contrôle des véhicules spatiaux. Mathématiques & Applications[END_REF], as well as the smallest singular value for the rank test.

Property on the first conjugate point. For the normal extremals in the oscillating case and the rotating case presented in Sec. 4.2.3, we take a large number of random initial adjoint vectors p(0) such that H 1 (q(0), p(0)) 2 + H 2 (q(0), p(0)) 2 = 1 and such that 0 < k(p(0)) < 1 where k is given by [START_REF] Wang | Necessary and sufficient conditions for local optimality of a periodic process[END_REF] for the oscillating case and by [START_REF] Zhitomirskiȋ | Typical singularities of differential 1-forms and Pfaffian equations[END_REF] for the rotating case. Then we numerically integrate the extremal system. We compute the first conjugate time t 1c , the pulsation ω = (p 4 (0) 2 + 4 p5 (0) 2 ) 1/4 , and the complete elliptic integral K(k), where k is the modulus given by [START_REF] Wang | Necessary and sufficient conditions for local optimality of a periodic process[END_REF] in the oscillating case or by [START_REF] Zhitomirskiȋ | Typical singularities of differential 1-forms and Pfaffian equations[END_REF] in the rotating case. Let γ(.) be a normal extremal starting at t = 0 from the origin and defined on [0, +∞[. As illustrated on Fig. 9, there exists a first conjugate point along γ corresponding to a conjugate time t 1c satisfying the inequalities: 0.34 ω t 1c -0.4 < K(k) < 0.53 ω t 1c -0.8 for the oscillating case, 0.33 ω t 1c + 0.16 < K(k) < 0.55 ω t 1c -1.27 for the rotating case.

Abnormal case. Fig. 10 illustrates the time evolution of the state variables. We check numerically the second order optimality conditions [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF]. Both the determinant test and the smallest singular value for the rank condition indicate that there is no conjugate time for abnormal extremals (Fig. 11). • The analytical expressions of θ 1 (t), θ 2 (t), given for the degenerate case and for the oscillating case in section 4.2.3, allow us to compute the strokes with small amplitudes of the nilpotent model. The results employing the Singular Decomposition Value (SVD) test (cf. [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF]) for conjugate points are also illustrated (see Fig. 6 and Fig. 7) showing that the simple loop have no conjugate points on [0, T ] while the eight stroke have a first conjugate point on [0, T ].

• The previous solutions are used to compute strokes for the Purcell swimmer with the T 0 (u 2 1 + u 2 2 )dt cost. More precisely, the initial adjoint vector p(0) of the nilpotent model gives a good initialization of the shooting algorithm used by HamPath to solve the following boundary value problem.

           q = ∂ Hn ∂p , ṗ = -∂ Hn ∂q , θ j (T ) = θ j (0) j = 1, 2, x(0) = y(0) = α(0) = 0, x(T ) 2 + y(T ) 2 = c 1 , α(T ) = c 2 , p θj (T ) = p θj (0) j = 1, 2, p α (0) = p α (T ) (33)
where T, c 1 , c 2 are fixed constants and Hn is the normal Hamiltonian associated with the

T 0 (u 2 1 + u 2
2 )dt cost. Then, with T fixed to 2π and c 2 to 0, we perform a discrete homotopy on the radius c 1 to obtain stroke with larger amplitudes (see Fig. 12). • We consider now the following optimal control problem

           q = ∂Hn ∂p , ṗ = -∂Hn ∂q , θ j (T ) = θ j (0) j = 1, 2, x(0) = y(0) = α(0) = 0, x(T ) 2 + y(T ) 2 = c 1 , α(T ) = c 2 , p θj (T ) = p θj (0), j = 1, 2, p α (0) = p α (T ) (34)
where H n = 1/2 a(q)u 2 1 + 2b(q)u 1 u 2 + c(q)u 2 2 is the true Hamiltonian associated with the mechanical cost, and u 1 , u 2 are the optimal controls. We take an extremal of (33) to initialize a discrete homotopy with parameter λ ∈ [0, 1], of the following optimal control problem

           q = ∂H λ ∂p , ṗ = -∂H λ ∂q , θ j (T ) = θ j (0) j = 1, 2, x(0) = y(0) = α(0) = 0, x(T ) 2 + y(T ) 2 = c 1 , α(T ) = c 2 , p θj (T ) = p θj (0), j = 1, 2, p α (0) = p α (T ) (35)
where H λ = λH n +(1-λ) Hn . When λ reaches the value 1, we obtain an extremal of (34). Since the latter homotopy is discrete, we may not follow a unique branch and obtain many kinds of strokes: Fig. 14 and Fig. 15 are two different strokes solutions of (34) and the SVD rank condition show that the only candidates for optimality are the simple loops. Then we perform a second homotopy on the radius c 1 to have a one-parameter family of strokes. Fig. 16 represents two one-parameter families of solutions of (34) corresponding respectively to the strokes of Fig. 14 and Fig. 15. To compare these two families of strokes, we compute in Fig. 17 their geometric efficiencies and we conclude that for a given radius r = c 1 , the corresponding stroke of the family for λ = 1 is more efficient. Conclusions. Two one-parameter families of simple loop for the mechanical cost appear to be good candidates for minimizers, and their respective efficiency is compared in Fig. 17. Note that the efficiency increases with the radius of the circle c 1 .

4.4 Sufficient second order conditions for the Purcell strokes Lemma 4.1 shows one symmetry for the Purcell swimmer. This model has actually a further symmetry: any time translation of the shape variables (θ 1 , θ 2 ) and the orientation variable α provides a stroke, which has the same cost. The presence of (two) symmetries requires (as in the case of the Copepod model) applying particular second order sufficient conditions [START_REF] Gavriel | Second order sufficient conditions for optimal control problems with nonunique minimizers: an abstract framework[END_REF]. For this purpose, we provide numerical results on second order sufficient conditions for normal extremals of the Purcell swimmer. We consider the optimal control problem in which we minimize the mechanical cost (5) over trajectories satisfying [START_REF] Bettiol | Sub-Riemannian geometry and swimming at low Reynolds number: the Copepod case[END_REF] and such that θ j (0) = θ j (T ) j = 1, 2, α(0) = α(T ),

x(0) = 0, y(0) = 0, x(T ) 2 + y(T ) 2 = r (r is fixed) .

The symmetry in the Purcell case is shown in the following proposition.

Proposition 4.7. Take I = (-ε, ε) for ε > 0 small enough. Consider a = (φ, σ) ∈ I × I. Let (q, ū, p) where q = ( θ1 , θ2 , x, ȳ, ᾱ) and ū = (ū 1 , ū2 ) be a normal extremal. Write θj (.), j = 1, 2, p(.), x(.), ȳ(.) and ū(.) their corresponding periodic extensions. For all a ∈ I and t ∈ [0, T ], we define q a (.) = (θ a 1 (.), θ a 2 (.), x a (.), y a (.), α a (.)) where u a j (t) = ūj (t + σ), θ a j (t) = θj (t + σ) j = 1, 2, x a (t) = cos(φ) (x(t + σ) -x(σ)) -sin(φ) (ȳ(t + σ) -ȳ(σ)), y a (t) = sin(φ) (x(t + σ) -x(σ)) + cos(φ) (ȳ(t + σ) -ȳ(σ)), α a (t) = ᾱ(t + σ) + φ.

(37)

Then, the normal extremal (q(.), ū(.), p(.)) is continuously embedded in the family of extremals (q a (.), u a (.), p a (.)) a∈I where p a (.) is the adjoint vector associated with (q a (.), u a (.)).

Consider the simple loop of Fig. 15 satisfying (36), written below (q(.), ū(.)), with the associated adjoint vector p(.). We have the following result.

Numerical result. The simple loop (q, ū) is a local weak minimizer. We prove this claim using again results by [START_REF] Gavriel | Second order sufficient conditions for optimal control problems with nonunique minimizers: an abstract framework[END_REF]. Straightforward computations and Prop. 4.7 above show that the Purcell model satisfies all the assumptions required by [START_REF] Gavriel | Second order sufficient conditions for optimal control problems with nonunique minimizers: an abstract framework[END_REF]Thm. 4.2]. Owing to the fact that there is no conjugate points for the normal extremal (q(.), ū(.), p) (cf. Sec. 4.3.2), the Riccati equation has a global symmetric solution. As for the Copepod, we compute numerically the matrix W , defined in [START_REF] Gavriel | Second order sufficient conditions for optimal control problems with nonunique minimizers: an abstract framework[END_REF]Equation (8)]. We set L s = { (y 0 , y T ) ∈ R 5 × R 5 | ∇ q0,q T c(q 0 , q T ) (y 0 y T ) = 0 } where c(q 0 , q T ) = (θ 1 (0) -θ 1 (T ), θ 2 (0) -θ 2 (T ), x(0), y(0), α(0) -α(T ), x(T ) 2 + y(T ) 2 -r). We take the matrix N s such that L s = Im(N s ). Now from the symmetry described in Prop. 4.7, we define Γ φ := ∇ φ q a (0) ∇ φ q a (T ) φ=0 , Γ σ := ∇ σ q a (0) ∇ σ q a (T ) σ=0 and Γ r := (Γ φ Γ σ ) .

We consider the linear subspaces We take two different relative tolerances for the integrator used to compute numerically the matrix W . Table 2 shows that the matrices W s = N s (W + W ) N s , Wφ = N φ (W + W ) N φ and Wσ = N σ (W + W ) N σ have zero eigenvalues (whose eigenvectors are Γ φ and Γ σ ). In particular, W s is not definite-positive hence the standard sufficient second order conditions [START_REF] Gavriel | Second order sufficient conditions for optimal control problems with nonunique minimizers: an abstract framework[END_REF]Thm. 3.4] are not applicable.

On the other hand, the refined second order sufficient conditions of [START_REF] Gavriel | Second order sufficient conditions for optimal control problems with nonunique minimizers: an abstract framework[END_REF]Thm. 4.2] are satisfied since the eigenvalues of W r = N r (W + W ) N r are positive. This confirms the claim stating that the reference extremal (q, ū) is indeed a local weak minimizer. Observe also that when we apply the refined second order conditions we have to take account of all the possible symmetries which generate the embedding in a family of weak normal extremals. Indeed, if we consider just one symmetry the refined second order sufficient conditions are not applicable, as the analysis of the matrices Wσ and Wφ above clearly show. Table 2: The standard second order sufficient conditions fail: W s has zero eigenvalues and the refined second order sufficient conditions are satisfied: W r is positive-definite.

Figure 2 :

 2 Figure 2: Normal strokes: simple loop (left), limaçon with inner loop (right) and eight case, and a two self-intersecting case (bottom). First conjugate points on [0, 2π] are computed with a svd test and they appear with a cross.

Figure 3 :

 3 Figure 3: Normal stroke where the constraints are satisfied: simple loop with no conjugate point on [0, T ].

Figure 4 :

 4 Figure 4: Norm of the initial adjoint vector p = (p 1 , p 2 , p 3 ) and value of p 1 (0) 2 for normal strokes such that H n = 1/2 and having different displacements, illustrating the non properness of the exponential mapping.

  along the hypotenuse x = -√ 2arctanh √ 2/2 and the total displacement is x ab (T ) 2.742.10 -1 . The length of a normal stroke γ is l(γ) = 2π 0

Figure 5 :

 5 Figure5: Efficiency curve (left) and the corresponding minimizing curve with the best performance (right). Note that the efficiency of the abnormal curve is 5.56e -2 vs of order 8.89e -2 for normal strokes.

Figure 10 :

 10 Figure10: Abnormal case: state variables for q(0) = (1, 0, 1, 0, 0), p(0) = (0, 0, -2, 1, 1).

Figure 11 :

 11 Figure 11: Abnormal case: the second order sufficient condition indicates there is no conjugate point.

Figure 6 :

 6 Figure 6: (left) Control and state physical variables in the degenerate case of the nilpotent approximation with an simple loop. (right) SVD test of conjugate points (no conjugate point on [0, 2π]).

Figure 7 :Figure 8 :

 78 Figure 7: Control and state physical variables in the oscillating case of the nilpotent approximation with an eight shape (the cross stands for the first conjugate point).

Figure 9 :

 9 Figure 9: Conjugate points of normal extremals with constant energy H 2 1 + H 2 2 = 1 in the oscillating case (left) and in the rotating case (right).

4. 3 . 2

 32 Optimal strokes using a discrete numerical homotopy Method description.

Figure 12 :

 12 Figure 12: One parameter family of simple loop strokes of the Purcell swimmer with the T 0 (u 2 1 + u 2 2 )dt cost. The continuation is performed on the constant c 1 where we fixed T = 2π and c 2 = 0.

Figure 13 :

 13 Figure 13: (left) Simple loop stroke for the Purcell swimmer minimizing the cost T 0 (u 2 1 + u 2 2 )dt, taking T = 2π, c 1 = 0.068, c 2 = 0 and imposing the periodicity on α. (right) Test of conjugate points (no conjugate point on [0, 2π]).

Figure 14 :

 14 Figure 14: (left) State variables for the Purcell swimmer minimizing the mechanical cost, taking T = 2π, c 1 = 0.058 and c 2 = 0. (right) Test of conjugate points (no conjugate point on [0, 2π]).

Figure 15 :

 15 Figure 15: (left) State variables for the Purcell swimmer minimizing the mechanical cost, taking T = 2π, c 1 = 0.065 and c 2 = 0. (right) Test of conjugate points (no conjugate point on [0, 2π]).

2 Figure 16 :Figure 17 :

 21617 Figure 16: Two families of normal strokes for the Purcell swimmer minimizing the mechanical cost. We fixed T = 2π and c 2 = 0 and the family of strokes is obtained by a continuation on c 1 .

L

  φ := L s ∩ {(y 0 , y T ) ∈ R 5 × R 5 | Γ φ (y 0 y T ) = 0}, L σ := L s ∩ {(y 0 , y T ) ∈ R 5 × R 5 | Γ σ (y 0 y T ) = 0}, L r := L s ∩ {(y 0 , y T ) ∈ R 5 × R 5 | Γ r (y 0 y T ) = 0}and the matrices N φ , N σ and N r such thatL φ = Im(N φ ), L σ = Im(N σ ), L r = Im(N r ).

Invoking again elliptic functions properties ( [START_REF] Lawden | Elliptic functions and applications[END_REF]) we have sin(ψ/2) = sn(u/k, k), cos(ψ/2) = cn(u/k, k) where u = ωt + ϕ 0 . Still H 1 and H 2 are elliptic functions of the first kind. The solutions of [START_REF] Passov | Supplementary notes to: Dynamics of Purcell's three-link microswimmer with a passive elastic tail[END_REF], q1 , q2 , satisfy the relations

where q10 and q20 are constant, K(k), E(k) are respectively the complete elliptic integrals of the first and second kind, Z(.) is the Jacobi's Zeta function.

Computations of strokes with small amplitudes using the nilpotent approximation

We recall that the physical variables q are related to q using the transformation ϕ. The adjoint variables p are obtained by a Mathieu transformation associated with ϕ. Strokes with small amplitudes such that q(0) = 0 are computed from the nilpotent approximation in the following way:

• Degenerate case: The corresponding solutions qi (t), i = 1, 2 of ( 27) yield the periodic variables θ i (t) = qi (t), i = 1, 2 of period 2π/p 3 (0). Moreover, the constants q10 , q20 , ϑ 0 may be chosen so that q(0) = (θ 1 (0), θ 2 (0), x(0)) = 0.

• Oscillating case:

The modulus k can be expressed as

and, computing k(p(0)) such that the linear terms of θ 1 (t) = q1 (ωt + ϕ 0 ), θ 2 (t) = q2 (ωt + ϕ 0 ) of (29) vanish, leads to periodic strokes with eight shapes of period

.

The constants q10 , q20 are chosen such that ϑ(0) = 0. The initial adjoint vector p(0) has to verify the conditions H 1 (q(0), p(0)) 2 + H 2 (q(0), p(0)) 2 = 1, k(p(0)) ∈ (0, 1) and p4 (0) 2 + 4 p5 (0) 2 = 0. We integrate numerically the stroke in the physical variables starting from (q(0) = 0, p(0)) and compute the first conjugate points on [0, T ] (see Fig. 7).

• Rotating case: The modulus k can be expressed as

We have θ 1 (t) = q1 (ωt+ϕ 0 ), θ 2 (t) = q2 (ωt+ϕ 0 ) where q1 , q2 are explicitly written in [START_REF] Vinter | Optimal control, Systems & Control: Foundations & Applications[END_REF]. We choose p(0) so that H 1 (q(0), p(0)) 2 + H 2 (q(0), p(0)) 2 = 1, k(p(0)) ∈ (0, 1) and such that the denominator of k(p(0)) is nonzero. As k(p(0)) tends to 0, the linear terms of q1 (u), q2 (u) of ( 30) tend to 0. This is the case when p4 (0) → 0 and p5 (0) → 0, and at the limit, equation ( 26) reduces to the equation of the degenerate case: θ = 0.

Abnormal case. We can reduce the problem considering the minimal time problem for the single-input affine system (cf. [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF]) q(t) = F1 (q(t)) + u(t) F2 (q(t))

where u is now a scalar control. Denoting q(.) a reference minimum time trajectory, since we consider abnormal extremals, it follows from the Pontryagin maximum principle that along the extremal lift of


Laure Buhry
email: laure.buhry@inria.fr

Clément Langlet

Francesco Giovannini

Modelling the effects of propofol on neuronal synchronization in network of interneurons

Propofol is a chemical agent commonly used as an intravenous general anesthetic. At the cellular level, this short-acting ananesthetic positively modulates GABAergic inhibitory activity by targeting GABA-A receptors [START_REF] Vanlersberghe | Propofol[END_REF]. This type of receptors are widespread in the brain and can be present both within synaptic clefts, as well as on extrasynaptic locations along the dendrites and neuron membrane where they are responsible for tonic inhibition. At the macroscopic level of SEEG (deep Stereographic-Electro Encephalogram) or EEG (Electro Encephalogram) recordings, one observes, with certain doses of propofol, a paradoxical excitation phenomenon [START_REF] Brown | General Anesthesia and Altered States of Arousal: A Systems Neuroscience Analysis[END_REF] the generation mechanisms of which are not clearly understood. In this study, we suggest a potential mechanism for the appearance of paradoxical excitation occurring under propofol-induced general anaesthesia.

We show, with a model network of Hodgkin-Huxley neurons, that tonic inhibition -induced by the binding of propofol to extra-synaptic receptors -together with an increase of the synaptic time constant within an certain range [START_REF] Mcdougall | Propofol enhances both tonic and phasic inhibitory currents in second-order neurons of the solitary tract nucleus (NTS)[END_REF] can account for the phenomenon of paradoxical excitation. However, changes in the gain (or conductance) of the synaptic inhibition do not correspond to a sudden increase in neuronal population firing rate nor synchrony as described in the experiments [START_REF] Mcdougall | Propofol enhances both tonic and phasic inhibitory currents in second-order neurons of the solitary tract nucleus (NTS)[END_REF]. The action of propofol on extrasynaptic GABAergic receptors was modelled by varying the conductance g ton of a tonic current in the form I ton = g ton (V -E ton ) as described in [START_REF] Hutt | Study of GABAergic extra-synaptic tonic inhibition in single neurons and neural populations by traversing neural scales: application to propofol-induced anaesthesia[END_REF]. Figure 1 shows the evolution of the neuronal population firing rate and the coherence (or synchrony) of the network activity as the tonic inhibition and the synaptic conductance vary. The plots are given for different values of the synaptic time constants. The increase of these three variables, synaptic time constant and conductance, and tonic conductance reflect an increase in propofol doses.

  Network frequency for τi = 10 ms.

  Network synchronisation for τi = 10 ms.

  Network frequency for τi = 14 ms.

  Network synchronisation for τi = 14 ms.

  Network synchronisation for τi = 30 ms.

FrequencyFigure 1 :

 1 Figure 1: Propofol-enhanced tonic inhibition allows for tighter network synchronisation, regardless of the presence of stronger inhibitory synapses. In all the plots, the x axis represents the tonic conductance (gton) and the y axis represents the inhibitory synaptic weight (w i ). (A) Given τ i = 10 ms, the network frequency decelerates as tonic inhibition strengthens until a critical value at which it accelerates. (B) This acceleration is due to an abrupt increase in network synchronisation at gton ≥ 15 nS for all values of w i . (C) A longer synaptic time constant (τ i = 14 ms) shifts the network frequency bump towards lower values of gton. (D) Similarly, the network synchronisation bump shifts towards lower values of gton. (E) Extending the synaptic time constant (τ i = 30 ms) causes the bump-like pattern of the network frequency to disappear in favour of a linearly decelerating trend. (F) Similarly, the bump-like pattern of the network synchronisation disappears in favour of a linearly decelerating trend. The network was stimulated with a constant current I stim = 0.4 nA.

Acknowledgements

Laure Buhry thanks the CNRS for the PEPS JCJC 2016.


Laurent Cognet
email: lcognet@u-bordeaux1.fr

Dmitri A Tsyboulski

John-David R Rocha

Condell D Doyle

James M Tour

R Bruce

Stepwise Quenching of Exciton Fluorescence in Carbon Nanotubes by Single Molecule Reactions

The chemical reactions of single molecules with individual single-walled carbon nanotubes are observed, and luminescence quenching analysis reveals the diffusional range of excitons in semiconducting nanotubes.

Optical excitation of semiconducting single-walled carbon nanotubes (SWNTs) generates relatively strongly bound excitons that have spatial dimensions predicted to be a few nanometers (1)(2)(3). Experimental evidence of efficient exciton-exciton annihilation in nanotubes indicates that SWNT excitons have significant mobility along the tube axis (4)(5)(6). However their mobility is still experimentally and theoretically uncertain. A notable related effect is the strong suppression of photoluminescence (PL) when SWNT sidewalls are perturbed by chemical reactions (7)(8)(9). This quenching phenomenon has hampered the use of covalently derivatized SWNTs as near-infrared (near-IR) fluorophores.

We report here the use of single-nanotube microscopy to detect stepwise changes in SWNT PL intensity within segments of individual nanotubes while they are exposed to chemical reactants. These stepwise changes in PL intensity are caused by reactions of single molecules with one nanotube. Since the pioneering low-temperature experiments of Orrit and Moerner, single-molecule spectroscopy has proven to be a powerful tool that bypasses ensemble averaging in the study of static and dynamic nano-objects in various environments (10,11). Single-molecule approaches are especially appealing for SWNT fundamental studies and applications (12)(13)(14) because the bulk samples are highly heterogeneous. In the present work, the magnitudes of PL intensity steps caused by single molecule reactions reveals that the exciton excursion range in highly luminescent SWNTs is ~ 90 nm and essentially independent of nanotube structure. This room-temperature excitonic motion is deduced to be diffusional. Because each nanotube exciton visits a very large number of atomic sites during its lifetime, PL quenching provides an ultrasensitive method for sensing and studying certain types of chemical reactions with nanotube sidewalls at the single-molecule level.

Our studies required highly luminescent and relatively long nanotubes that were immobilized yet accessible to added reactant solutions. We therefore used very brief tip ultrasonication to disperse raw HiPco SWNTs in aqueous sodium dodecylbenzenesulfonate (SDBS) surfactant before mixing the dilute suspension with a melted agarose gel preparation [START_REF]Materials and methods are available as supporting material on Science Online[END_REF]. Aqueous gels are commonly used in molecular biology to provide an inert immobilizing environment and they have previously been applied in single-molecule optical measurements [START_REF] Dickson | [END_REF]. We captured near-IR fluorescence images of single nanotubes with a wide-field inverted microscope modified to include laser excitation and a low-noise InGaAs camera (14). Micrometerslong individual SWNTs were easily identified in these images and some of them displayed uniform and high PL intensity (Fig. 1A). The (n,m) structural assignment of each studied tube was deduced from the peak wavelength of its narrow Lorentzian emission spectrum (17), as measured with a dedicated spectrograph and multichannel near-IR detector (Fig. 1B). In this study, we selected the brightest nanotubes present in the weakly sonicated samples. We believe that these selected nanotubes have low defect densities and nonradiative relaxation rates that are dominated by intrinsic processes.

Most of these nanotubes displayed PL that varied linearly at the excitation intensities used here and remained stable for tens of minutes under continuous illumination (Fig. 1B andC). We used low laser intensities (~100 W/cm 2 ) to avoid exciton-exciton annihilation effects that would shorten the exciton lifetime and reduce the exciton's excursion range (4)(5)(6)(7). Furthermore, we verified that the nanotubes studied here showed no detectable spectral shifts, spectral broadening, or change in emission intensity during the typical duration of our experiments (Fig. 1B). Finally, we confirmed that the emission spectra of SWNTs embedded in agarose gels matched those of SWNTs in fluid surfactant suspension.

To initiate chemical reactions, we placed a 10-µL drop of reactant solution on the edge of the 20 mm by 20 mm sample gel while recording fluorescence images of a specific nanotube with stable and spatially uniform emission (Fig. S1). Two reactants were used. One was sulfuric acid (0.05 M), which causes pH-dependent, reversible quenching of SWNT fluorescence (7,8,18). The other was 4-chlorobenzenediazonium tetrafluoroborate (1 mg/mL in water), which is known to give irreversible covalent derivatization of carbon atoms on the nanotube sidewall (19). We recorded fluorescence image sequences at an 18 Hz frame rate and plotted the total intensity within selected 2 × 2 pixel regions to obtain PL intensity traces. The chosen spatial binning region (4 pixels) corresponded to 670 nm by 670 nm, which is essentially the diffraction-limit at the near-IR wavelengths.

The PL intensity trace from a 670 nm segment of an individual (8,6) nanotube after exposure to sulfuric acid is shown on Fig. 2A. That nanotube's luminescence spectrum is plotted in Fig. 2B. Coarsely, the PL time trace displays an exponential decay form. However, closer examination reveals a series of strikingly distinct steps between well-defined intensity levels, rather than a continuous decrease (see the inset of Fig. 2A). Moreover, the luminescence shows both upward as well as downward steps.

In order to analyze the distribution of step amplitudes, we compiled a histogram for the first 108 s period showing signal differences (δI) between successive image frames (Fig. 2C). Apart from the large main peak near zero, which arises mainly from signal noise, there are four distinct side bands distributed symmetrically around zero. These side bands are the signature of specific step amplitudes. The inner two, (labeled +1 and -1) correspond to positive and negative steps of the same unit intensity, whereas the outer ones (+2 and -2) correspond to double-amplitude steps. The presence of discrete (quantized) intensity changes is confirmed by the precisely linear relation linking the amplitudes of these intensity steps (Fig. 2D).

We attribute the observed discrete intensity steps to individual protonation reactions at the nanotube surface. The -1 and -2 steps would arise, respectively, from one or two independent sidewall reactions occurring within the 54 ms sampling interval, whereas the +1 and +2 steps arise from the reverse chemical processes. The mechanism through which acids quench nanotube luminescence has been described as a protonation at the sidewall of the nanotubes (7,8,18). In this scheme, a hole is injected into the nanotube π-system near the protonation site. If an exciton encounters such a chemically induced hole before it radiatively recombines, the exciton's luminescence will be efficiently quenched through nonradiative Auger processes (6).

Because the excitons have a limited excursion range, only a fraction of those generated within an observed 670 nm-long segment will be quenched by a single protonation-induced hole. If the density of such holes is low, then each one causes the same amount of quenching (corresponding to the observed -1 step amplitude), and each de-protonation restores the same amount of luminescence (the observed +1 step amplitude). However, when the typical distance between quenching sites becomes shorter or comparable to the exciton excursion range, the quenching impact of additional sites is reduced, giving smaller and eventually unresolved steps, as seen at later times in Fig. 2A.

To confirm that the observed stochastic stepwise quenching arises from reversible protonation, we have added ~10 μL of 1 M NaOH to nanotube samples that had previously been quenched by acid addition. The down triangle in Fig. 2E marks the point of acid addition, and the up triangle marks that of base addition. The period of luminescence restoration is expanded in Fig. 2F, which shows the expected stochastic stepwise intensity changes. Statistical analysis of upward and downward steps in a medium of known pH can in principle provide the pKa of the protonated nanotube. However, we could not perform this analysis on our data because of time-dependent pH gradients in the sample. Large intensity fluctuations are observed after luminescence has been restored by base addition. We attribute this behavior to local pH fluctuations near the nanotube within the agarose pores after the successive addition of acid and base without active mixing. Finally, the higher intensities seen in Fig. 2E after restoration likely reflect a higher final pH level at the nanotube (18).

To compare these exciton quenching effects with those caused by an irreversible chemical reaction, we performed similar measurements on samples exposed to diazonium salts (Fig. 3) (20). These traces also display distinct steps in luminescence intensity, but the great majority of the steps are negative (>95%), indicating a mainly irreversible process (9). Interestingly, adjacent diffraction-limited segments within a single nanotube exhibit very different dynamics, with little or no correlation between the time positions of the steps (compare each set of red and black traces in Fig. 3,A to C). Similar spatially uncorrelated behavior in acid quenching is vividly illustrated by the video in Supporting Online Materials [START_REF]Materials and methods are available as supporting material on Science Online[END_REF]. This observation demonstrates that the mean excursion range of excitons must be smaller than our optical resolution of 670 nm. The same conclusion can be drawn from the localized PL in bent nanotubes excited by linearly polarized light (( 14) and Fig. S2).

Comparison of the PL time traces within each frame of Fig. 3, A through C reveals that for different segments of one nanotube, the first few step heights are integer multiples of the same unit intensity change Δ (regions marked by a single star). This stepwise quenching pattern leads to emission intensity plateaus that match between distinct segments of the same tube. As was observed for acid quenching, the steps have identical amplitudes at low densities of quenching sites, and gradually smaller amplitudes at longer times when the typical distance between quenching sites becomes shorter than or comparable to the exciton excursion range. However, in contrast to the reversible protonation reaction of acid quenching, the single-molecule diazonium salt reactions are essentially irreversible. Here, the chlorobenzenediazonium ion (or its diazotate dimeric product) is thought to first adsorb onto the surface of a nanotube, forming a charge transfer complex. The complex then irreversibly decomposes to form a covalent bond with the nanotube surface (9,21).

We presume that either of these nanotube perturbations can cause efficient exciton nonradiative decay. A mechanism for luminescence quenching by sidewall derivatization is suggested by scanning tunneling spectroscopy results on nanotubes covalently functionalized with amido groups, obtained with 10 nm resolution. The semiconducting energy gap was locally filled at the derivatization site (22). Excitons encountering such electronic perturbations should undergo efficient nonradiative recombination.

For each nanotube studied, repeatable values of Δ are obtained at early stages of the reaction when only a few independent sites have been derivatized within the observed diffraction-limited segment. Significantly, we find that Δ values match for separate segments in a single nanotube (Fig. 3,A to C). As for the minor variations in intensity step heights found at early stages of the reaction, we attribute these to the finite length of the diffraction-limited segment, since quenching events near the edge of the region produce smaller steps. Interestingly, several simultaneous steps sometimes occur within a single 54 ms frame, (Fig. 3 A andB). This behavior probably results from nonuniform dilution of reactants in the sub-micrometer pores of the agarose gels, exposing the nanotube segments to large local concentration fluctuations.

Normalizing Δ by the initial luminescence intensity, I, gives a ratio representing the probability that a mobile exciton within the nanotube segment encounters a chemical reaction site during its lifetime. This probability also represents the fraction of the nanotube length explored by an exciton. The average exciton excursion range, Λ, is therefore simply given by ( )

I L Δ ×
, where L is the observed segment length. For L = 670 nm, we have measured Λ values for 19 different highly luminescent SWNTs. The resulting distribution is strikingly narrow, with a mean value of 91 nm (Fig. 4A). The small dispersion found for Λ values suggests that Λ does not depend strongly on nanotube structure. We plotted Λ values measured for 11 different (n,m) species (23) as a function of SWNT diameter (Fig. 4B) and found no systematic variation of exciton excursion range with diameter (for the 0.7 to 1.1 nm range studied) or with chiral angle.

As mentioned above, our determination of Λ assumes that the local perturbation induced by a reaction sufficiently modifies the electronic structure of these bright tubes (22) to completely quench excitons at that site. If this quenching probability, P, is in fact less than 1, then Λ would be given by 90 nm / P (but could not exceed 670 nm).

However, several observations support the estimate of P = 1. The results are independent of nanotube structure; the normalized initial step amplitudes are consistent for acid quenching and diazonium salt quenching; and Auger quenching of excitons is known to be efficient, Is the exciton motion is ballistic or diffusional? For the bright emissive SWNTs studied here, we can estimate the exciton recombination lifetime, τ, as 100 ps, which is among the larger values reported from time-resolved spectroscopy (24,25), If the motion is ballistic, then the exciton's kinetic energy would be given by If efficient exciton quenching occurs not only at chemical derivatization sites, but also at the ends of cut nanotubes, our value of the exciton diffusion range, Λ, naturally explains the strong reduction in PL reported for SWNTs shorter than ca. 100 nm (31). More generally, the measured diffusion range implies that one exciton will visit ca. 10 4 carbon atoms during its lifetime. Thus PL becomes remarkably sensitive to certain sidewall electronic perturbations, which may include chemical derivatizations and defects introduced during nanotube processing. PL may then prove useful for detecting local pH gradients in restricted environments, such as microfluidic channels or organelles inside biological cells.

τ

  effective mass. Using an estimated * exciton m of ~0.1 m e (1,26), we obtain ballistic E ~ 10 -7 eV. This value is orders of magnitude lower than the 0.012 eV (½ k B T) expected from rapid thermalization induced by the exciton-phonon coupling evident in SWNT spectroscopy (27). We conclude that the excitonic motion is instead diffusional, in agreement with a recent investigation on bulk samples (28). The corresponding diffusion coefficient for the exciton 1more than two orders of magnitude below a prior tentative estimate that was indirectly deduced from transient absorption experiments on ensembles of short SWNTs in annealed dried films (29). The exciton diffusion coefficient measured here for highly luminescent individual SWNTs may be used to interpret the ~150 cm -1 Lorentzian linewidths observed in emission spectroscopy of room-temperature SWNTs (see Figs. 1Band 2B). These widths arise from dephasing that is far more rapid (70 fs) than the ~100 ps exciton population lifetime τ.If dephasing instead occurs by diffusional hopping of the exciton along the tube axis, then the hopping step length d hop can be estimated from the relation represents the interval between hops, identified as the 70 fs dephasing time. Interestingly, this approach gives a hopping step length of ~2 nm, which closely matches theoretical estimates of the exciton size(1)(2)(3). We propose that the emission linewidth may originate from thermally assisted, short-range hopping of excitons along the nanotube axis. At cryogenic temperatures, suppression of this hopping would then remove the major source of line broadening and account for dramatically narrower spectral features (30).

Fig. 1 .

 1 Fig. 1. Fluorescence measurements on an individual SWNT immobilized in agarose gel.

Fig. 2 .

 2 Fig. 2. Reversible stepwise quenching of individual SWNTs by acid. (A) Luminescence

Fig. 3 .

 3 Fig. 3. Irreversible localized quenching of individual SWNTs with 4-

Fig. 4 .

 4 Fig. 4. Exciton excursion range (Λ) measured for various SWNTs. (A) The distribution (n = 19) of Λ values deduced from elementary step amplitudes (see text). (B) Plot of Λ