AN OPTIMAL CONTROL APPROACH TO CILIARY LOCOMOTION

. We consider a class of low Reynolds number swimmers, of prolate spheroidal shape, which can be seen as simpliﬁed models of ciliated microor- ganisms. Within this model, the form of the swimmer does not change, the propelling mechanism consisting in tangential displacements of the material points of swimmer’s boundary. Using explicit formulas for the solution of the Stokes equations at the exterior of a translating prolate spheroid the govern- ing equations reduce to a system of ODE’s with the control acting in some of its coeﬃcients (bilinear control system). The main theoretical result asserts the exact controllability of the prolate spheroidal swimmer. In the same geo-metrical situation, we consider the optimal control problem of maximizing the eﬃciency during a stroke and we prove the existence of a maximum. We also provide a method to compute an approximation of the eﬃciency by using ex- plicit formulas for the Stokes system at the exterior of a prolate spheroid, with some particular tangential velocities at the ﬂuid-solid interface. We analyze the sensitivity of this eﬃciency with respect to the eccentricity of the considered spheroid and show that for small positive eccentricity, the eﬃciency of a prolate spheroid is better than the eﬃciency of a sphere. Finally, we use numerical optimization tools to investigate the dependence of the eﬃciency on the number of inputs and on the eccentricity of the spheroid. The “best” numerical result obtained yields an eﬃciency of 30 . 66% with 13 scalar inputs. In the limiting case of a sphere our best numerically obtained eﬃciency is of 30 . 4%, whereas the best computed eﬃciency previously reported in the literature is of 22%.


(Communicated by Enrique Fernández-Cara)
Abstract. We consider a class of low Reynolds number swimmers, of prolate spheroidal shape, which can be seen as simplified models of ciliated microorganisms. Within this model, the form of the swimmer does not change, the propelling mechanism consisting in tangential displacements of the material points of swimmer's boundary. Using explicit formulas for the solution of the Stokes equations at the exterior of a translating prolate spheroid the governing equations reduce to a system of ODE's with the control acting in some of its coefficients (bilinear control system). The main theoretical result asserts the exact controllability of the prolate spheroidal swimmer. In the same geometrical situation, we consider the optimal control problem of maximizing the efficiency during a stroke and we prove the existence of a maximum. We also provide a method to compute an approximation of the efficiency by using explicit formulas for the Stokes system at the exterior of a prolate spheroid, with some particular tangential velocities at the fluid-solid interface. We analyze the sensitivity of this efficiency with respect to the eccentricity of the considered spheroid and show that for small positive eccentricity, the efficiency of a prolate spheroid is better than the efficiency of a sphere. Finally, we use numerical optimization tools to investigate the dependence of the efficiency on the number of inputs and on the eccentricity of the spheroid. The "best" numerical result obtained yields an efficiency of 30.66% with 13 scalar inputs. In the limiting case of a sphere our best numerically obtained efficiency is of 30.4%, whereas the best computed efficiency previously reported in the literature is of 22%.
1. Introduction and statement of the main results. Ciliates are swimming microorganisms which exploit the bending of a large number of small and densely packed organelles, termed cilia, in order to propel themselves in a viscous fluid. In this work we consider an envelope model for such ciliary locomotion, where the dynamics of the individual cilia are replaced by time periodic tangential displacements of the points on the boundary of the microorganism. We refer to the works of Taylor [29], Blake [6], Childress [9], Ishikawa, Simmonds and Pedley [15] or to the recent review paper of Lauga and Powers [18] for a detailed description of this model. The aim of this work is to study the locomotion mechanism of these organisms, combining classical tools of low Reynolds number flow theory with the modern techniques of optimal control. The main novelty we bring is the analysis of a non spherical case, which is closer to the elongated form of most ciliated microorganisms than the spherical one previously considered in the literature.
Our main theoretical result gives the controllability of the considered system (roughly speaking, this means that for every initial position P 1 and every final position P 2 there exists a time-periodic motion of the cilia steering the organism from P 1 to P 2 ). We also consider the optimal control problem, with state and control constraints. This problem is to maximize the efficiency of the swimmer. Tackling this optimal control problem with a reasonable computational cost requires semiexplicit formulas for the solutions of the Stokes system in the exterior of a prolate spheroid, with boundary conditions which correspond either to a translation of the spheroid or to an arbitrary tangential velocity field on the surface. The first type of boundary conditions have been tackled in classical reference such as Happel and Brenner [14] whereas tackling the tangential boundary conditions requires very technical calculations involving Gegenbauer functions. Moreover, the study of the problem for spheroids close to a sphere requires a careful asymptotic analysis of these special functions. Our results can be seen as extensions of San Martín, Takahashi and Tucsnak [25], Sigalotti and Vivalda [27]. In these two papers, the controllability is obtained by imposing the boundary velocity in (1.10), but this velocity is not related to a displacement of the points on the boundary, and in particular, (1.14) is not imposed. From an optimization perspective, our results are an extension of the study done in Lauga and Michelin [22] in the spherical case.
Related results, namely considering swimmers with shapes which change during the stroke, have been obtained in Shapere and Wilczek [26], Alouges, DeSimone and Lefebvre [3,4], Lohéac, Scheid and Tucsnak [21], Lohéac and Munnier [19]. Let us quote also the following results for the controllability of a three-spheres swimmers in the presence of boundaries in the fluid domain: Alouges and Giraldi [5], Gérard-Varet and Giraldi [13]. Concerning the optimal control problem associated with the swimming, Lohéac and Scheid [20] is interested in the time optimal control problem whereas Alouges, De Simone and Heltai [2] deals with a numerical method based on boundary integral formulation in order to optimize a stroke in the case of a three-spheres swimmer.
In the case of very high Reynolds number (perfect fluid), one can also consider the swimming problem and the corresponding control problem can be dealt in a similar way as in the case of very low Reynolds number: see Chambrion and Munnier [7,8], Kanso, Marsden, Rowley, and Melli-Huber [17], Munnier [23], Munnier and Pinçon [24].
To give the precise statement of our results we need some notation from spheroidal geometry. We first recall that a spheroid is a revolution ellipsoid. In Cartesian coordinates (x 1 , x 2 , x 3 ), the equation of a spheroid with Ox 3 as the symmetry axis is The spheroid is called prolate if a 1 < a 3 and oblate if a 1 > a 3 . To study the Stokes equations at the exterior of a prolate spheroid, it is convenient to use prolate spheroidal coordinates, denoted by (η, θ, ϕ). Following, for instance, Dassios et al. [11] or [14], these coordinates are related to the usual Cartesian ones by where the focal distance c is a fixed positive number and The corresponding orthonormal basis is given by e ϕ = −(sin ϕ) e 1 + (cos ϕ) e 2 . The level lines θ = θ 0 , η > 0, ϕ = 0 and θ ∈ [0, π], η = η 0 , ϕ = 0 are described in Figure 1, for various values of θ 0 , η 0 and for a fixed value of the focal distance c.   Figure 2 represents a family of spheroids for various values of c and of η 0 such that a 3 = c cosh η 0 = 1 and a 3 = c cosh η 0 → 1. In the limit case c → 0 (so that η 0 → ∞), we recover the unit sphere of R 3 .
Assume that at t = 0 the swimmer occupies the closed bounded set S 0 delimited by the spheroid, centered at the origin, of equation η = η 0 , for some η 0 > 0. For t > 0 we denote by S(t) the closed set occupied by the swimmer at instant t. Within Blake's envelope model, the propelling mechanism of ciliates consists of time periodic tangential displacements of the points on the swimmer's boundary.
In this article, we make the assumption that these displacements are azimuthally symmetric with respect to the axis Ox 3 of the unit vector e 3 . This implies, in particular, that at each instant t, the position of the center of mass of the swimmer is h(t)e 3 , where h a real valued function, and that the domain S(t) occupied by the swimmer at instant t is simply obtained by a translation of the vector h(t)e 3 of S 0 , i.e., we have Our approach cannot be directly adapted for displacements of the boundary points which are not azimuthally symmetric with respect to an axis. Indeed, we need closed form solutions of the Stokes system which are not available, in general, for velocity fields on the boundary without axial symmetry.
Let v (respectively p) be the Eulerian velocity (respectively the pressure) field in the fluid at instant t. These fields are defined, for each instant t, in the time dependent domain R 3 \ S(t). For the remaining part of this work, we choose to work, instead of v and p, with the fields defined by (1.4) This allows us to take the space variable x in the fixed domain R 3 \ S 0 .
The model that we consider for the displacements of the boundary points of S 0 is the following: for each t ≥ 0, the point x ∈ ∂S 0 , of prolate spheroidal coordinates (η 0 , ξ, ϕ), is displaced to a point of ∂S 0 whose prolate spheroidal coordinates are and where χ(·, t) is a diffeomorphism from [0, π] onto itself. These displacements at the boundary lead to a tangential velocity V χ e θ where (e η , e θ , e ϕ ) is the orthonormal basis associated with the prolate spheroidal coordinates (η, θ, ϕ): The full system describing the coupled motion of the swimmer and of the surrounding fluid writes In the above system, we have denoted by µ the viscosity coefficient of the fluid. We have denoted by σ(v, p) the field of Cauchy stress tensors in the fluid defined by Moreover, as in the remaining part of this work, we have denoted by n the unit normal vector field to ∂S 0 directed to the interior of S 0 .
Remark 1.1. Note that, in the system (1.7)-(1.11), we have used the low Reynolds number approximation, in which all the inertia terms are neglected. Moreover, the equilibrium condition for the angular momentum is automatically satisfied due to the azimuthal symmetry of the problem.
A sufficient condition for χ(·, t) to be a diffeomorphism from [0, π] onto itself is The controllability result corresponds to the following problem: assume h 0 , h 1 ∈ R, find T and χ such that χ(·, 0) = χ(·, T ), (1.14) and such that the solution (v, p, h) of the system (1.7)-(1.11) with h(0) = h 0 satisfies In order to solve the above control problem, we can suppose that the function χ in (1.5) can be written as 16) where N is a positive integer and the given functions (g i ) 1≤i≤N are supposed to be smooth, with g i (0) = g i (π) = 0 for every i ∈ {1, . . . , N }. We also add to (1.7)-(1.11) the equationα (1.17) We can also use the following initial conditions (1.18) Here and in what follows α stands for the vector valued function (α i ) 1≤i≤N . In that case, equation (1.14) can be written as This means that we want to control the final value of the functions (α i ), and it is thus more convenient to use (β i ) as inputs of our system. As we shall see in We are now in a position to state our main controllability result.
where the functions (g i ) have been introduced in (1.16).
Then for every h 1 ∈ R and ε > 0, there exist T > 0 and β ∈ C ∞ [0, T ], R N such that the solution (v, p, h, α) of (1. In the above theorem and in all what follows, we denote by | · | the Euclidean norm of R k , k ≥ 0. Note that, for ε small enough the third condition in the above theorem implies condition (1.13). An important ingredient in the proof of Theorem 1.1 is to write h as a function of χ. In order to do this, we use the particular shape of our swimmer and the fact that in that case we can have an explicit formula for the solution of the following Stokes problem by using the Gegenbauer functions defined in Appendix A.
Remark 1.2. In the above theorem, we use the symmetries to consider only the problem of a translation in one direction. The same question could be tackled for an arbitrary shape of the swimmer and for any translation and any rotation, but we would need explicit formulas for Stokes' systems such as (1.19)-(1.20) and the corresponding proof should be more involved (more Lie brackets should be computed, see the proof of Theorem 1.1).
The second part of this work is devoted to the maximization of the swimmer's efficiency, which is classically defined (see, for instance, [6]) as the ratio between the average power that an external force would spend to translate the system rigidly at the same average speed and the average power expended by the swimmer during a stroke starting and ending at the non-deformed shape. This study is more complicated than in the spherical case (see [6] and [22]) since there is no explicit formula of the solution of the Stokes system at the exterior of a prolate spheroid with a tangential boundary condition for the velocity. Therefore, an important part of our effort has been devoted to obtain explicit formulas for an approximation of this solution.
Let us consider (v, p, h, χ) the solution of (1.6), (1.7)-(1.11), (1.14), (1.15) and h(0) = 0. Then the efficiency eff of a motion driving the mass center from the origin to the final position h 1 e 3 , with h 1 = 0, in time T is given by (1.21) To derive the formula for the numerator in the right hand side of (1.21), we have used the fact that the mean velocity necessary to swim from the origin to h 1 e 3 in time T is h1 T e 3 , so that the work needed to translate the prolate spheroid from the origin to a position centered in h 1 e 3 in time T is given by Let us note that, this efficiency depends on the final position h 1 and on the final time T , but also on the geometric parameters c and η 0 (see (1.3)). An important problem is to maximize the efficiency eff with respect to the input χ which varies in a set of admissible controls. In the present case, given T > 0 and h 1 ∈ R, we choose the set M ε,r of admissible controls to be the subset of all the functions χ satisfying χ − Id L ∞ (0,T ;H 2 (0,π)) ≤ ε, ∂χ ∂t L 2 (0,T ;W 1,∞ (0,π)) ≤ r, for ε small enough and given r > 0. The result in Theorem 1.1 implies that for every h 1 there exists T > 0 such that M ε,r = ∅. We show the existence of controls maximizing the efficiency, in the sense below: Assume M ε,r is defined as above. Then there exists χ * ∈ M ε,r such that eff(χ * ) = sup We then consider a numerical method to compute the efficiency defined by (1.21). The idea is to decompose the solution as v =ḣ(t)v (0) + v (1) , p =ḣ(t)p (0) + p (1) . (1.24) where (v (1) , p (1) ) is the solution of the Stokes system The above problem can be approximated by using the Gegenbauer functions defined in Appendix A. More precisely, such an approximation can be obtained by inverting a tridiagonal finite-dimensional matrix L. Since this matrix may not be always invertible, we only consider the case of e := c 2 , small enough. We derive the formulas when a 3 = 1 or equivalently (according to (1.3)) In that case, the tridiagonal matrix L is invertible and we obtain an approximation for (v (1) , p (1) ) (see Proposition 5.1) and for the efficiency (see Section 6). We investigate in Section 7 the influence of the focal distance c on the swimmers efficiency. We show that the approximated efficiency of the swimmer increases when the sphere deforms to a prolate spheroid. Finally, in Section 8 we propose a methodology to solve the above mentioned optimal control problem numerically and we provide some of the results of our simulations. One of the main difficulty we have to tackle is the presence of constraints on the state and on the control. The maximum efficiency obtained numerically is of 30.66%, in the case in which 13 scalar input functions are used. In the limiting case of a sphere our best numerically obtained efficiency is of 30.4%, whereas the best computed efficiency previously reported in the literature, see [22], is of 22%.
2. Notation and preliminaries. In this paper, the usual notation is used for Lebesgue and Sobolev spaces on a domain Ω. In particular, the homogeneous Sobolev space of order (k, q), k ≥ 1, 1 < q < ∞, on Ω is denoted by D k,q (Ω) with associated seminorm |u| k,q,Ω = |i|=k D i u L q (Ω) . Whenever confusion will not arise, we shall omit Ω in the above norms and seminorms. Throughout the paper we shall use the same font style to denote scalar, vector and tensor-valued functions and corresponding function spaces. Consider the Stokes system, with Dirichlet boundary conditions, posed in Ω = R 3 \ S 0 , where S 0 can be a general compact connected set, with non empty interior. The results of this section will be applied in the case where S 0 is the closed bounded set delimited by a prolate spheroid. The following well-posedness result for the exterior Stokes system is well-known (see, for instance, Galdi

This solution satisfies
(1 + |x|)v ∈ L ∞ (Ω) Moreover, there exists c(Ω) > 0 such that: Using the above proposition, we deduce the following classical result. Its proof is standard and is thus skipped.
In particular, Remark 2.1. The above result implies in particular that one can write the efficiency defined by (1.21) as Consequently the inverse of the efficiency is a quadratic functional of the velocity of the fluid but its dependence on the control χ and on the rigid velocityḣ is more complicated.
In the case of an obstacle S 0 which is the closure of the interior of a prolate spheroid, more detailed information on the solution (v, p) of (2.1) can be obtained by using special functions.
To this aim, following Dassios et al. [11] and [14], we modify the prolate spheroidal coordinates which have been introduced in (1.2) and we set With this modification, the change of coordinates (1.2) becomes ) with x i given by (2.7), the natural basis corresponding to the above coordinate system is given by

JORGE SAN MARTÍN, TAKÉO TAKAHASHI AND MARIUS TUCSNAK
Moreover, the Lamé coefficients (also called scale factors) corresponding to the change of coordinates (2.7) are given by It can be easily checked that {g ζ , g τ , g ϕ } form an orthogonal basis in R 3 . Writing we obtain a right-handed orthonormal basis {e ζ , e τ , e ϕ } expressed by Note that the prolate spheroid equation (1.1) can be rewritten as where τ 0 > 1 is fixed. A simple calculation shows that the surface element of the prolate spheroid obtained by taking τ = τ 0 in (2.7) is given by In this work we are interested in solutions of (2.1) which are axisymmetric with respect to Ox 3 . This means, in particular that the velocity field v writes where e ζ and e τ have been introduced in (2.10) and (2.11) and the pressure satisfies p = p(ζ, τ ). It is not difficult to check that the (τ, ζ) and (ζ, τ ) components of the gradient of v are given by where the Lamé coefficients Ξ τ and Ξ ζ are given by (2.9). Inserting the expression of Ξ τ and Ξ ζ from (2.9) in the last formula, we obtain It follows that the component σ τ,ζ := (σe τ ) · e ζ of the corresponding stress tensor is given by The axisymmetric case is characterized (see [14]) by the existence of a stream function ψ. According to [11] this function is related to the velocity field v by Moreover, the stream function ψ lies in a kernel of a relatively simple differential operator. More precisely, consider the differential operator E 2 which associates to each smooth function f of variables ζ and τ the function E 2 f defined by Then (see, for instance, [14, p.104]) the fact that (v, p) satisfies the Stokes system is equivalent to the fact that ψ satisfies the equation where E 4 is the square of the operator E 2 introduced in (2.19). More precisely, v is given through formulas (2.17)-(2.18), (2.14), and p can be obtained from the integration of the relation ∇p = −µ∆v (see [14, 4-15, p.116]).
In the case of a prolate spheroid translating with a constant velocity equal to e 3 we have the following result.
is given, using prolate spheroidal coordinates, by where G 2 and H 2 are the second order Gegenbauer functions defined in Appendix A. Moreover, the corresponding stress tensor σ (0) satisfies, on ∂S 0 , (2.22) Finally, the force exerted by the solid on the fluid is

23)
Proof. We refer to [14,Section 4.30] for the formula of the stream function (2.21) and for the formula (2.23). More precisely, it is proved in [14,Section 4.30] that (2.24) satisfies (2.20) and that the velocity v (0) associated with By inserting (2.25) in the last term in the right hand side of (2.16) we obtain Concerning the first term in the right hand side of (2.16), we first tackle the part which can be evaluated directly using (2.25), i.e., we write Finally, we note that (2.17) implies that where ψ (0) is given by (2.21). Using next (2.21), together with the fact (following from (A.28)) that On the other hand, (2.17) and (2.25) imply that Using the last two formulas in (2.28), it follows that  .23), it follows that the work needed to translate the prolate spheroid from the origin to a position centered in (0, 0, h 1 ) in time T is given by 3. Proof of Theorem 1.1. The aim of this section is to prove Theorem 1.1. As a first step, we rewrite the system (1.6), (1.7)-(1.11), (1.16) and (1.17) and prove it well-posedness. First, we justify the formula (1.6) of V χ (see (1.10)). This formula uses the prolate spheroidal coordinates (1.2) and, by using the prolate spheroidal coordinates (2.7), it can be written as Denoting by x 1 (t), x 2 (t) and x 3 (t) the Cartesian coordinates corresponding to the points moving at the boundary of S 0 , we have, using (1.2), that Differentiating the above formula with respect to t, we obtain that the boundary velocity of these moving points is given by V χ (θ, ϕ, t)e θ , with the formula (1.6).
By using (2.23) and the symmetry, the above equation is equivalent tȯ We can then use Proposition 2.2 to write Since n = −e τ and e θ = −e ζ , we deduce from the above relation that We can combine (2.13) (2.22), (3.1) and the above equation to deduce The above equation and (3.5) yielḋ Relation (3.3) follows now by making the change of variables ξ = χ −1 (arccos ζ, t)) in the above formula.
In the case χ is given by (1.16), Proposition 3.1 can be used to write the system in a form which allows the application of some classical control theoretical tools. (3.14) We note that Relation (3.9) follows now by an integration by parts in (3.14).  (3.15) where (A n ) n≥1 is a sequence in L ∞ ([0, ∞); 2 ). Plugging (3.15) into (3.7), and using (A.19), we obtainḣ Since we have reduced the governing equations to a finite dimensional affine control system without drift, we are now in a position to use a well-known result due to Chow. We recall this result below and we refer, for instance, to Trélat [30, chap. 5, Proposition 5.14, p. 89], Jurdjevic [16] or Coron [10, Theorem 3.19, p.135]) for more information on this theorem. Theorem 3.1. Let m, n ∈ N * and let (f i ) i=1,...,n be C ∞ vector fields on R n . Consider the control system, of state trajectory X, In the above theorem, the Lie-bracket of two N -dimensional smooth vector fields f 1 , f 2 is a new vector field defined by We recall that a Lie-algebra is a space closed for the Lie-bracket [·, ·] and that Lie{f 1 , . . . , f m }, [·, ·] is the smallest Lie-algebra containing {f 1 , . . . , f m }. Moreover, Lie X {f 1 , . . . , f m } is the subspace of R n spanned by all the values in X of the vector fields in Lie{f 1 , . . . , f m }.
We are now in a position to prove Theorem 1.1.
Proof of Theorem 1.1. It clearly suffices to prove the result for N = 2. In that case, we only need to assume π 0 [g 2 (ξ)g 1 (ξ) − g 1 (ξ)g 2 (ξ)] sin 2 ξ dξ = 0. (3.19) We first remark that, for N = 2, equations (3.9)-(3.10) can be written in the form (3.17), provided that we set X = h α and where the functions F 1 and F 2 have been defined in (3.12) and (3.13). It is easily checked that The above formula, (3.18) and (3.20) yield that On the other hand, going back to (3.12) and (3.13) and using the fact that F (x) = − sin 2 x 2 , we obtain that The above formula and (3.21) imply that Using (3.19) and (3.20) and [f 1 , f 2 ](X) span R 3 . By continuity, there exists ε ∈ (0, ε] such that the same property holds for X = h α , with |α| < ε. We can thus conclude the proof by using Proof. Let us rewrite the efficiency defined by (1.21) by using the decomposition (3.4).
(4.6) Using Sobolev embeddings, we have that for ε small enough, if χ ∈ M ε,r then χ is a bijection of [0, π] onto itself and there exists a constant C = C(ε) such that χ −1 L ∞ (0,T ;W 1,∞ (0,π)) ≤ C. In what follows, we fix ε small enough so that the above result holds true. In particular, the velocity V χ given by (1.6) is in H 1 (∂S 0 ). From Theorem 1.1, we already know that for T large enough, and for ε, r > 0, the set M ε,r is non empty.
We can define h ∈ H 1 (0, T ) bẏ Then, from (4.10) and (4.9), we deduce Let us define (v, p) as the solution of (1.7)-(1.10) associated with χ and h defined above. Then we see that (v, p, h) is a solution of the control problem and in particular, χ ∈ M ε,r .

5.
The Stokes problem at the exterior of a prolate spheroid with a tangential velocity on the boundary. In order to tackle the second term in the right hand side of the boundary condition (1.10) we investigate in this section the problem (2.1), with S 0 being the prolate spheroid of equations (2.7) and with the boundary velocity v * is azimuthally symmetric with respect to e 3 and tangential to ∂S 0 . Unlike in the case of a spherical swimmer, studied in [6], we have no longer a separation of variables (see [11]). Therefore there is no simple way of expressing the solution of (2.1) in function of the coefficients occurring in the decomposition of v * in an appropriate basis. This leads us to consider an approximation of the exact solution involving several tridiagonal matrices.
To describe this approximation we use the prolate spheroidal coordinates introduced in (2.6), (2.7), (2.8) to write the boundary condition (1.10) as In Proposition 2.3 we give a formula for the solution of the Stokes system (2.1) with v * = e 3 and the main aim of this section is to provide an approximate solution of (2.1) with v * = u(ζ)e ζ ζ ∈ [−1, 1], (5.1) where u is a smooth function.
We first give, following [11], an exact formula for a class of fields (v, p) satisfying the Stokes equation outside S 0 and with the normal component of v vanishing on ∂S 0 . We next show that for every small enough value of the the geometric parameter e = c 2 we can choose one of the fields (v, p) in the above mentioned family such that the tangential trace of v on ∂S 0 has, up to an error term of order O(e), is equal to u in (5.1).
Let us consider the functions ψ n defined by for n ≥ 4 and ). Using the above notation we provide below, by a slight variation of arguments in [11], an explicit formula for the solution of (2.1) with a particular boundary condition.
Proof. Let us first prove that ψ n defined by (5.2) satisfies (2.20). In order to do this, we use (A.22) and (A.28) to deduce that In particular, for n ≥ 2, Using (5.9), we obtain, for n ≥ 4, Then, combining the above equality with (5.10) and (5.11) we obtain The recurrence relations (A.38) and (A.39) transform the above relation into from which we deduce (2.20) (using again (5.9)). For n = 2, 3, similar calculations give the same result (we use (A.37) instead of (A.38)).
If the shape of the prolate is close to a ball, one can obtain asymptotic formulas of the above results. More precisely, in what follows, we set e := c 2 , and we assume cτ 0 = 1 so that We next study the above formulas and their consequences when e → 0. We first normalize the functions ψ n defined in ( (1 − eζ 2 )(1 − ζ 2 ) a n G n−2 (ζ) + b n G n (ζ) + c n G n+2 (ζ) e ζ , (5.14) with a 2 = a 3 = 0, a n = −eα n H n (1/ √ e)C n−2 (e) (n ≥ 4), , we obtain that .
In particular, we deduce from (5.12) that v n (ζ, τ 0 ) with (5.15)-(5.17). Now, we want to consider a linear combination of the v n in order to approximate the boundary condition (1.10). More precisely, we will choose (B n ) n∈{2,...,K} so that K n=2 B n v n is close to velocity produced by the cilia motions, that is In order to do this, we use that (G n ) n≥2 is an orthogonal basis of L 2 (−1, 1) (for the weight w(x) = 1/(1 − x 2 )) to write A n G n (ζ). (5.23) With this decomposition, the velocity produced by the cilia motions can be written as More precisely, one has the following result (whose proof is classical and omitted).
Lemma 5.3. Assume (5.23). Then the coefficients A n are given by B n ψ n is a stream function corresponding to a solution (v (1) , p (1) ) of where the remaining term R K (e) is such that Proof. Using (A.37) and (A.38), we deduce that for K ≥ 4, where α n = α n for n ≥ 4 and α 2 = α 3 = 0. The above relation can be transformed as A n−2 β n−2 G n (ζ).
6. Approximation of the swimming efficiency. In this section, we consider an approximation of efficiency for the prolate in the case e small. As in the previous section, we consider the case We also recall (see Section 1) that the motion of the points on the swimmers surface is defined, using prolate spheroidal coordinates, by and that the input functions of the considered system are the functions Given the function (β i ) 1≤i≤N and an eccentricity e, the corresponding efficiency eff(χ) is defined by (1.21) In all what follows, we write eff (β; e) since χ can be obtained from β and to emphasize the dependence on e.
For K ∈ N we approximate the efficiency defined in (1.21) by the function β e → eff K (β; e): (1) ) is the solution of (5.36)-(5.37) considered in Proposition 5.1 and that is obtained from (A n ) 2≤n≤K ; • v (1) A n G n (ζ)e ζ is a truncation of the trace of v (1) on ∂S 0 ; • h is defined by (3.16) and h(0) = 0.
The approximate character of eff K is due to the replacement of the last occurrence of v (1) by its truncation v (1) K in the definition (1.21) of eff (β; e). This allows us to express eff K using the coefficients (A n ) 2≤n≤K defined by (5.25). To accomplish this goal we first prove the following result.
(6.6) On the other hand, using (2.17) and (6.1), it follows that Using (2.13) and the fact that n = −e τ , it follows from the above relation that For the second integral in the right hand side of (6.8), we use (A.21) and (5.37), A 2 n n(n − 1)(2n − 1) .
Inserting the last formula and (6.9) in (6.8) we obtain the conclusion.
The result below provides an useful expression for the efficiency defined in (6.3).
Proposition 6.2. Let (β i ) 1≤i≤N be a family of functions in L ∞ (0, ∞). Then the swimmer's efficiency, defined by (6.3), is given by: Proof. We first note that, applying (2.23), the numerator of the right hand side of (1.21) writes On the other hand, from (3.16) it follows that

AN OPTIMAL CONTROL APPROACH TO CILIARY LOCOMOTION 321
The last two formulas imply that To compute the denominator of the right hand side of (1.21) we first note that (1) ,ḣp (0) + p (1) )n dΓ = 0 (6.17) from (1.11).
On the other hand, by linearity, The last formula, (6.16) and (6.3) clearly imply that which yields the conclusion.
7. Sensitivity with respect to the focal distance. In this section, we study the behavior of the function eff K defined in (6.3) when e → 0. This means that we consider a family of prolate spheroids like those in Figure 2, i.e., with e = c 2 → 0, cosh η 0 = τ 0 = 1 c .
Theorem 7.1. The efficiency eff K defined by (6.3) satisfies with eff is the efficiency of the spherical swimmer. More precisely, if S 0 is the sphere centered at the origin and of radius 1, then one can consider the same problem as for the prolate: σ(v, p)n dΓ = 0 (t ≥ 0), (7.12) h(0) = 0, α(0) = 0, (7.13) , t ≥ 0). (7.14) Then, one can decompose A n (t)G n (cos θ) (7.15) and in that case one can obtain (see, for instance, [22]) that the efficiency defined by (1.21) satisfies Let us note that in [22], the tangential velocity is written by using the Legendre function of first kind θ → P 1 n (cos(θ)) (n ≥ 1), but that the obtained formulas can be seen as a particular case of the corresponding results for the prolate spheroidal case (see (A.5)). Indeed, it suffices to note that, using (A.20) and (A.22), we have P 1 n (cos(θ)) = n(n + 1) sin(θ) G n+1 (cos(θ)) (n ≥ 1).
Finally, it is easy to check that the efficiency in (7.16) is less than 1/2 (see [22]). Obtaining such a theoretical upper bound is an interesting open question in the prolate spheroidal case. From the numerical point of view, see Section 8, getting close to the theoretical upper bound seems difficult because of the singular behavior of the optimal stroke (we refer again to [22]).
8. Numerical optimization of the prolate spheroidal swimmer efficiency. We first recall the main steps to obtain the efficiency eff K defined in (6.3). We fix g i , i = 1, . . . , N such that the hypothesis of Theorem 1.1 holds. Then, for β ∈ L ∞ (0, T ; R N ), we compute Then we use formula (5.25) π 0 g i (ξ)G n (cos(χ(ξ, t))) ∂χ ∂ξ (ξ, t) dξ that gives in particular (see (3.16))  One of the difficulties in computing the efficiency of the prolate spheroidal swimmer is to evaluate formulas containing H n (1/ √ e) which, for small e, may lead to important numerical errors. Therefore, for small e, we use the asymptotic formula (derived in Theorem 7.1), eff K (β; e) ≈  To evaluate eff K , for given β and e, we discretize the time interval [0, T ] and the angle interval [0, π] by using uniform meshes of size N T and N ξ . With this full discretization the original optimal control problem, of input function β, reduces to an optimal control problem where the unknown is a matrix of size N × N T (we search for an approximation of β i (t j )). The constraints in this problem are the injectivity of χ and the periodicity of α (i.e. α(T ) = 0). This maximization is performed using the IPOPT (Interior Point Optimizer) package ( [31]).
We take g i (ξ) = sin(iξ), i ≥ 1. Table 1 contains the results of numerical computations of the optimal efficiency for various values of the number of scalar controls and with the discretization parameters K = 50, N T = 250, N ξ = 500. The above values of K, N T and N ξ are large enough (for the considered values of N ) to make the results stable (not more than 0.05% of modification of the efficiency) with respect to coherent augmentations of these parameters. In the above statement the term "coherent augmentations" signifies that increasing the number K of basis functions implies (due to high frequency oscillations) a significant augmentation of N ξ .
In Figure 3 we describe the evolution of h with respect to time and the evolution of the cilia (of χ) at several times: t = 0, t = T /16, t = T /8 et t = 3T /16.
We also represent in Figure 4 the evolution for N = 5 of h with respect to time and the evolution of the cilia (of χ) at several times: t = 0, t = T /4, t = T /2 et t = 3T /4. The same graphics are drawn in Figure 5 for N = 13.  Table 2. Maximum displacement for the optimal eff K for different values of N .  Our first conclusion is that, in the spherical case, we obtain numerically an efficiency of 30.4%, with 13 scalar input functions and a maximum displacement of 104 • . We remark in Figure 5 that the obtained displacement field has similar features with the "shock structure" obtained in [22] (where the maximal obtained efficiency is 22% with a maximal displacement 52.6 • ). We believe that obtaining numerically, in the spherical case, an efficiency closer to the theoretical bound of 50% requires to increase the number of modes N and, consequently, refining the discretization. Accomplishing this program by using the full-discretization and the standard optimization algorithms proposed in this paper seems a difficult computational issue.  For the prolate spheroidal ciliate we have shown, both analytically and numerically, that the optimal efficiency increases by small augmentations of the eccentricity e. It would be interesting to perform similar computations for larger values of e. From a theoretical view-point, the main difficulty in this case would be to establish that the matrix L in (5.31)-(5.34) is invertible, whereas computationally one should use "exact" formulas instead of the asymptotic ones in order to evaluate the Gegenbauer functions.
Finally, let us mention that, according to the left pictures in Figures 3, 4, 5 the position of the mass center of the ciliate goes backward at a certain stage of the stroke. We think that this feature should appear independently of the number of modes (although in a less obvious manner, as we can see in Figure 5) and that a mathematical proof of this assertion is an interesting challenge.
Appendix A. Some background on Legendre and Gegenbauer functions. We gather in this section some definitions and properties of functions of Legendre and Gegenbauer type. For details on this rich topic, including the proofs of the results we state here, we refer, for instance, to Abramowitz and Stegun [1] or Whittaker and Watson [32].

(A.2)
It is well-known that the Legendre polynomials satisfy the differential equation (1 − x 2 )P n (x) = −n(n + 1)P n (x) (A. 3) and that they satisfy the orthogonality conditions