B Bonnard

J.-B Caillau

G Picot

Geometric and numerical techniques in optimal control of the two and three-body problems

The objective of this article is to present geometric and numerical techniques developed to study the orbit transfer between Keplerian elliptic orbits in the two-body problem or between quasi-Keplerian orbits in the Earth-Moon transfer when low propulsion is used. We concentrate our study on the energy minimization problem. From Pontryagin's maximum principle, the optimal solution can be found solving the shooting equation for smooth Hamiltonian dynamics. A first step in the analysis is to find in the Kepler case an analytical solution for the averaged Hamiltonian, which corresponds to a Riemannian metric. This will allow to compute the solution for the original Kepler problem, using a numerical continuation method where the smoothness of the path is related to the conjugate point condition. Similarly, the solution of the Earth-Moon transfer is computed using geometric and numerical continuation techniques.

Introduction

In this article we consider the orbit transfer in the two and three-body problem, using low propulsion. In the first case, the model is given by Kepler equation

q = - q |q| 3 + u m
where m represents the mass of the satellite, subject to ṁ = -δ|u| modelling fuel consumption. The control satisfies the constraint |u| ≤ ε where ε is a small parameter.

The physical optimal control has to maximize the final mass which leads to min u(.) t f 0 |u|dt * Mathematics Institute, Bourgogne University & CNRS, 9 avenue Savary, F-21078 Dijon (bernard.bonnard@u-bourgogne.fr).

† Same address (jean-baptiste.caillau@u-bourgogne.fr), supported by Conseil Régional de Bourgogne (contract no. 2009-160E-160-CE-160T).

‡ Same address (gautier.picot@u-bourgogne.fr), supported by CNRS (contract no. 37244) and Conseil Régional de Bourgogne (contract no. 079201PP02454515).

where t f is the fixed transfer time. From Pontryagin maximum principle [START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF], fixing the boundary conditions-e.g. a transfer from a low eccentric to a geostationary orbit-, an optimal solution can be numerically computed using a shooting algorithm. This leads to a complicated numerical problem. In [START_REF] Gergaud | Homotopy method for minimum consumption orbit transfer problem[END_REF], the following numerical scheme is proposed. One computes the optimal solution using the convex homotopy min u(.)

t f 0 [λ|u| 2 + (1 -λ)|u|]dt, λ ∈ [0, 1],
which amounts to regularizing the L 1 -minimization problem into an L 2 -problem. This was the starting point of the use of continuation methods in orbital transfer, when low propulsion is applied-see also [START_REF] Caillau | Contribution à l'étude du contrôle en temps minimal des transferts orbitaux[END_REF] for the use of the continuation method in the time minimal control problem, where the homotopy parameter is the bound of the maximal amplitude of the thrust. From the mathematical point of view, the original dynamics associated with the optimal flow is replaced by another Hamiltonian one, and a continuation is made to solve the shooting equation.

The first motivation of this article is to present a neat geometric result from [START_REF] Bonnard | Geodesic flow of the averaged controlled Kepler equation[END_REF]: Neglecting the mass variation, restricting to coplanar transfer (the inclination being considered as a homotopy parameter) and replacing the L 1problem by an averaged L 2 -problem, one can substitute the Hamiltonian vector field defined by the maximum principle with

H = 1 2n 5/3 9n 2 p 2 n + 5 2 (1 -e 2 )p 2 e + 5 -4e 2 2 p 2 θ e 2 ,
where n is the mean motion, e the eccentricity, and θ the angle of the pericenter (the singularity e = 0 corresponds to circular orbits). Coordinates (n, e, θ) are moreover orthogonal coordinates for the Riemannian metric associated with H, g = dn 2 9n 1 3

+ 2n

5 3

5(1 -e 2 ) de 2 + 2n

5 3 5 -4e 2 dθ 2 .
Such a metric is isometric to

g = dr 2 + r 2 c 2 (dϕ 2 + G(ϕ)dθ 2 )
where r = 2 5 n 5/6 , ϕ = arcsin e and c = 2/5, G(ϕ) = 5 sin 2 (ϕ) 1 + 4 cos 2 (ϕ)

•

The Hamiltonian flow -→ H is Liouville integrable and the metric in the above normal form captures the main properties of the averaged orbital transfer. Indeed, one can extract from g the following two-dimensional Riemannian metrics:

g 1 = dr 2 + r 2 dψ 2 which is associated with the orbital transfer where θ is kept fixed (this encompasses the case of circular targets). Such a metric is flat and geodesics are straight lines in suitable coordinates.

g 2 = dϕ 2 + G(ϕ)dθ 2 which represents the restriction of the metric to r 2 = c 2 and describes by homogeneity the orbit transfer in the general case.

A generalization of the results of [START_REF] Sinclair | The cut locus of a two-sphere of revolution and Toponogov's comparison theorem[END_REF] will allow to compute for the metric g 2 the conjugate and cut loci and to get a global optimality solution for the averaged optimal control problem. This is the starting point to analyze the original optimal control problem using a continuation method.

A second motivation of this article is to present some results from geometric control theory connected to our analysis with adapted numerical codes developed to compute the solutions. First of all, the maximum principle is only a necessary optimality condition. In order to get sufficient optimality conditions under generic assumptions one must define the concept of conjugate point, associated with the energy minimization problem. This concept was already introduced in the standard litterature of calculus of variations [5]. If the Hamiltonian optimal dynamics is described by a smooth Hamiltonian vector field -→ H , conjugate points are the image of the singularities of the exponential mapping: exp x(0) : p(0) -→ Π x exp t f -→ H (x(0), p(0)) where Π x : (x, p) -→ x is the canonical projection. Such points can be numerically computed using the code [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF]. An important remark, in view of the use of the (smooth) continuation method in optimal control is to observe that the shooting equation is precisely to find p(0) such that exp x(0) (p(0)) = x 1 where x 1 is the terminal condition and the derivative is generated using the variational equation of -→ H . This will lead to convergence results for the smooth continuation method in optimal control, related to estimates of conjugate points.

The last section is devoted to the Earth-Moon transfer, using low propulsion. The model is the standard circular restricted model [START_REF] Meyer | Introduction to Hamiltonian dynamical systems and the N -body problem[END_REF] where the two primaries are fixed in a rotating frame. Up to a normalization the system can be written in Hamiltonian form,

ẋ = -→ H 0 (x) + u 1 -→ H 1 (x) + u 2 -→ H 2 (x)
where x = (q, p) ∈ R 4 and the drift -→ H 0 is given by

H 0 (x) = 1 2 (p 2 1 + p 2 2 ) + p 1 q 2 -p 2 q 1 - 1 -µ ρ 1 - µ ρ 2 ,
q being the position of the spacecraft, ρ 1 representing the distance to the Earth with mass 1-µ located at (-µ, 0), and ρ 2 the distance to the Moon with mass µ located at (1 -µ, 0), µ ≃ 1.2153e -2 being a small parameter. The Hamiltonian fields associated with the control are given by

H i (x) = -q i , i = 1, 2,
and the control bound is |u| ≤ ε. The parameter µ is small and this remark was used by Poincaré to study the dynamics of the free motion described by -→ H 0 by making a deformation of the case µ = 0 which corresponds to Kepler equation in rotating coordinates [START_REF] Meyer | Introduction to Hamiltonian dynamical systems and the N -body problem[END_REF]. Inspired by this approach and using our preliminary geometric analysis, we propose a simple solution to the Earth-Moon transfer using low propulsion for the energy minimization problem.

2 Geometric and numerical methods

Maximum principle

We consider the energy minimization problem min u(.) t f 0 |u| 2 dt, for a smooth control system of the form

ẋ = F 0 (x(t)) + m i=1 u i (t)F i (x(t)) = F (x(t), u(t)), x ∈ X.
The set of admissible controls is the subset U of measurable bounded mappings u(.) with corresponding trajectory x(.) defined on the whole interval [0, t f ]. Pontryagin maximum principle [START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF] tells us that Proposition 2.1. If (x, u) is an optimal pair on [0, t f ], there exists a non trivial pair (p 0 , p), p 0 ≤ 0 and p an absolutely continuous adjoint vector valued in T * X, such that on [0, t f ] we have

ẋ = ∂H ∂p (x, p, u), ṗ = - ∂H ∂x (x, p, u), (1)
and

H(x, p, u) = max v∈R m H(x, p, v) (2)
where

H(x, p, u) = p 0 m i=1 u 2 i + < p, F (x, u) >.
Definition 2.1. The mapping H from T * X × R m to R is called the pseudo-Hamiltonian. A triple (x, p, u) solution of (1-2) is called an extremal trajectory.

Computation of extremals

From the maximization condition (4), one deduces that ∂H/∂v = 0, and there are two types of extremals:

-Abnormal extremals. They correspond to the situation p 0 = 0 and are implicitely defined by the relations H i = 0, i = 1, . . . , m, where H i =< p, F i (x) > are the Hamiltonian lifts.

-Normal extremals. If p 0 < 0, it can be normalized to -1/2 by homogeneity.

From ∂H/∂v = 0, one deduces u i = H i for i = 1, . . . , m, and plugging such H i into H defines a true smooth Hamiltonian

H n = H 0 + 1 2 m i=1
H 2 i whose solutions are the normal extremals.

The concept of conjugate point

Definition 2.2. Let z = (x, p) be a normal reference extremal defined on [0, t f ].

The variational equation

δz(t) = d -→ H n (z(t))δz(t)
is called the Jacobi equation. A Jacobi field is a non-trivial solution δz = (δx, δp). It is said to be vertical at time t if δx(t) = dΠ x (z(t))δz(t) = 0 where Π x is the projection (x, p) → x.

The following standard geometric result is crucial.

Proposition 2.2. Let L 0 be the fiber T * x0 X and L t = exp t ( -→ H n )(L 0 ) be its image by the one-parameter subgroup generated by -→ H n . Then L t is a Lagrangian submanifold whose tangent space at z(t) is generated by the Jacobi fields which are vertical at t = 0. Definition 2.3. We fix x 0 = x(0) and define for t

∈ [0, t f ] the exponential mapping exp x0,t (p 0 ) = Π x (z(t, z 0 ))
where z(t, z 0 ), with z 0 = (x 0 , p 0 ), denotes the normal extremal departing from z 0 when t = 0.

Definition 2.4. Let z = (x, p) be the reference normal extremal. The time t c ∈ [0, t f ] is called conjugate if the mapping exp x0,tc is not an immersion at p(0). The associated point x(t c ) is said to be conjugate to x 0 . We denote by t 1,c the first conjugate point and by C(x 0 ) the conjugate locus formed by the set of first conjugate points when we consider all normal extremals starting from x 0 .

The conjugate time notion admits the following generalization.

Definition 2.5. Let M 1 be a regular submanifold of M , and let us define

M ⊥ 1 = {(x, p) ∈ T * M | x ∈ M 1 , p ⊥ T x M 1 }. Then t foc ∈ [0, t f ]
is called a focal time if there exists a Jacobi field J = (δx, δp) such that δx(0) = 0 and J(t foc ) is tangent to M ⊥ 1 .

Remark 2.1. The concept of conjugate point is related to the necessary and sufficient optimality conditions, under generic assumptions, see for instance [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF].

Conjugate points and smooth continuation method

Smooth continuation is a general numerical method to solve a system of equations F (x) = 0 where F : R n -→ R n is a smooth mapping, see [START_REF] Allgower | Figure 22: Norm of extremal control, Earth-Moon transfer, µ = 1.2153e -2[END_REF]. The principle is to construct a homotopy path h(x, λ) such that h(x, 0) = G(x) and h(x, 1) = F (x) where G(x) is a map having known zeros, or where the zeros can be easily computed using a Newton type algorithm. The zeros along the path can be calculated by different methods, the simplest being a discretization

0 = λ 0 < λ 1 < • • • < λ n = 1
of the homotopy parameter where, at step i + 1, the zero computed at step i is used to initialize Newton algorithm. The approach has to be adapted to optimal control problems: The shooting equation comes from the projection of a symplectic mapping, the Jacobian can be computed using Jacobi fields and one must consider the central extremal fields associated with the problem (see [START_REF] Bonnard | The smooth continuation method in optimal control with an application to control systems[END_REF]). A short description of the method is given below in our case study.

Shooting equation

We consider a family H λ , λ ∈ [0, 1], of smooth Hamiltonians on T * X associated with normal extremals of an energy minimization problem. We fix the boundary conditions x 0 , x 1 and the transfer time t f . This leads to a family exp λ x0,t f (p 0 ) of exponential mappings. Using the notation E λ : p 0 -→ exp λ x0,t f (p 0 ), one must solve the shooting equation E λ (p 0 ) = x 1 .

Proposition 2.3. For each λ, the shooting equation is of maximal rank if and only if the point x 1 is not conjugate to x 0 for the corresponding λ. Moreover, in this case, the solutions of the shooting equation contain a smooth branch, which can be parameterized by λ and the derivative E ′λ can be generated integrating the Jacobi equation.

From the above proposition, to ensure convergence of the method one must control -the distance to the conjugate loci, -that the branch has is defined on the whole interval [0, 1].

The second point is related to two standard problems in optimal control: Existence of Lipschitzian minimizers-hence solutions of the maximum principle [START_REF] Guerra | Existence and Lipschitzian regularity for relaxed minimizers[END_REF]-, and compactness of the domain of the exponential mapping. Next we present a nice geometric situation for which convergence of the method is ensured.

Definition 2.6. Consider the normal extremal field -→ H n of an energy minimization problem with fixed final time t f . Given an initial condition x 0 , the separating locus L(x 0 ) is the set of points where two distinct normal extremal curves intersect with same cost. The cut point along a normal extremal is the first point where it ceases to be optimal. The cut locus Cut(x 0 ) is the set of such points when we consider all extremals initiating from x 0 and losing optimality exactly at time t f .

Convergence of the continuation method in the Riemannian case

We first recall that the Riemannian problem can be, at least locally, reset in the following framework. Let F 1 , . . . , F n be a set of smooth vector fields on a manifold X and assume that they are linearly independent. One can define a Riemannian metric on X by asserting that {F 1 , . . . , F n } form an orthonormal frame. Introducing the control system dx(t)/dt = i . Fixing the level set H n = 1/2 parameterizes trajectories by arc length. For the energy minimization problem, the transfer time can be arbitrarily prescribed.

n i=1 u i (t)F i (x(t)), the length of the curve x(•) is l(x) = T 0 n i=1 (u 2 i (t))
Theorem 2.1. Let g λ , λ ∈ [0, 1], be a smooth family of complete Riemannian metrics on X. Let us fix the initial point x 0 . Denote i λ (x 0 ) the distance from x 0 to the cut locus Cut λ (x 0 ) and by i λ = inf x0 i λ (x 0 ) the injectivity radius of the corresponding metric. Then, -for length shorter than inf λ i λ (x 0 ), the continuation method with initial condition of the shooting equation at x 0 converges, -for length shorter than inf λ i λ the continuation method converges for every initial condition of the shooting equation.

Remark 2.2. In the Riemannian case, the situation is neat. Completeness leads to existence of smooth normal minimizers, the domain of the exponential mapping is a sphere, and estimates of the injectivity radius are related to the curvature tensor. In general, such estimates are a difficult problem and a pragmatic point of view is to have numerical approximations [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF].

3 The energy minimization problem in orbital transfer with low thrust

Preliminaries

Neglecting the mass variation and restricting to the coplanar case, the system is represented in Cartesian coordinates by q = -q |q| 3 + u where q = (q 1 , q 2 ) is the position and x = (q, q) ∈ R 4 is the state. We denote by H 0 (q, q) = (1/2) q2 -1/|q| the Hamiltonian of the free motion. We have the following first integrals:

-C = q ∧ q (momentum), -L = -q/|q| + q ∧ C (Laplace integral).

Proposition 3.1. The domain Σ e = {(q, q) | H < 0, C = 0}, called the elliptic domain, is filled by elliptic orbits and to each orbit (C, L) corresponds a unique (oriented) ellipse.

To represent the space of ellipses, one introduces the following geometric coordinates:

-the semi-major axis of the ellipse a, related to the semi-latus rectum P by the relation

a = P/ √ 1 -e 2 ,
-the argument of the pericenter θ,

-the eccentricity e, the eccentricity vector being (e x , e y ) = (e cos θ, e sin θ).

To represent the position of the satellite we use the longitude l ∈ S 1 , while l ∈ R takes into account the rotation number and is called the cumulated longitude.

Observe that e = 0 corresponds to circular orbits. The control u can be decomposed into moving frames attached to the satellite, the two standard frames being -the radial-orthoradial frame {F r , F or } where F r = (q/|q|) ∂/∂ q,

-the tangential-normal frame {F t , F n } where F t = ( q/| q|) ∂/∂ q.

Averaging of periodic sub-Riemannian problems

Let X be an n-dimensional smooth manifold and let F i (l, x), i = 1, . . . , m be smooth vector fields parameterized by l ∈ S 1 that set up a constant rank m distribution on X (see also Remark 3.2). Given a positive pulsation ω on S 1 ×X relating the time t and the angle l according to

dl = ω(l, x)dt, (3)
one defines a periodic sub-Riemannian problem as follows: Given two points x 0 and x f on the manifold, minimize the L 2 dt control norm of trajectories connecting the two points and associated with the previous vector fields,

dx dl = m i=1 u i F i (l, x), u ∈ R m , min u(.) t f 0 |u| 2 dt = l f 0 |u| 2 dl ω(l, x) (here |u| 2 = m i=1 |u i | 2 ).
The total angular length, l f > 0 is fixed, implicitly defining t f through (3).

Pontryagin maximum principle asserts that minimizing trajectories are projection of Hamiltonian curves on the cotangent bundle (extremals), z = (x, p), such that

dx dl = ∂H ∂p (l, x, p, u), dp dl = - ∂H ∂x (l, x, p, u),
where

H(l, x, p, u) = p 0 |u| 2 + m i=1 u i H i (l, x, p), H i = p, F i (l, x) , i = 1, . . . , m,
is the Hamiltonian parameterized by l ∈ S 1 , u ∈ R m , and a non-positive constant p 0 . We restrict the analysis to the normal case, p 0 < 0, and pass to affine coordinates in (p 0 , p) setting p 0 = -1/2. Pontryagin maximization condition that, almost everywhere along an extremal,

H(l, z(l), u(l)) = max v∈R m H(l, z(l), v),
implies that the control is given by the dynamic feedback

u(l, z) = ω(l, x)(H 1 , . . . , H m )(l, z), (4)
and that z is an integral curve of the maximized Hamiltonian

H n (l, z) = ω(l, x) 2 m i=1 H 2 i (l, z).
We are interested in the behaviour of solutions for large angular length, so we set ε = 1/l f as the small parameter and renormalize the problem in the following manner-typical of systems with two time scales. In contrast to l which is the fast time, define the slow time s = εl in [0, 1], and renormalize variables on the cotangent bundle by x = x, p = p/ε. Lemma 3.1. In the renormalized variables,

d z ds = -→ H n (s/ε, z). (5)
Proof. Obvious since H n (l, x, •) is quadratic in the adjoint state p.

For fixed z, H n (•, z) is a smooth function on S1 and can be expanded into its Fourier series. In particular, one can define its first coefficient or average,

H(z) = 1 2π 2π 0 H n (l, z) dl.
Since taking symplectic vector field and averaging readily commute, it is well known [2, 20] that H trajectories are good approximations of those of H n as ε goes to zero.

Theorem 3.1. Given any initial condition, the solution z ε of (5) converges uniformly on [0, 1] towards the solution z of the averaged Hamiltonian.

Remark 3.1. The same rate of convergence holds for the cost which can be added as a new state in the augmented system. But as the integrand depends on the fast variable l through the pulsation ω, there is a priori no higher order approximation of the performance index [4,[START_REF] Chaplais | Averaging and deterministic optimal control[END_REF].

As for H n (s/ε, z) which converges to H in the large space of Schwartz distributions, 1 one can only expect weak convergence on the control u ε of the original problem on [0, 1/ε]. We address two questions: first, can we describe the fast oscillating control by means of Fourier series with slowly varying coefficients and, secondly, what is the asymptotic behaviour of u ε ∞ when ε → 0? The importance of such an estimate relies on the fact that one has in practice to estimate l f so as to meet a given requirement on the L ∞ -norm of the control.

We first note that the average system provides a sub-Riemannian approximation of the original one.

Proposition 3.2. The averaged Hamiltonian is a nonnegative quadratic form in the adjoint. If constant, its rank is not less than m, and H can locally be written as a sum of squares, thus defining a sub-Riemannian (Riemannian if k = n) problem.

Proof. Clearly, H(x, •) remains quadratic nonnegative in p by linearity and positivity of the integral and

ker H(x, •) = l∈S 1 ker H n (l, x, •)
so that, if constant, the rank k is at least equal to m. In this case, the quadratic form can always be decomposed into a sum of k squares in a chart by taking the square root of the associated symmetric nonnegative matrix.

Remark 3.2. (i) The increase in the rank is related to the generation of Lie brackets of the following system with drift: Set x = (l, x) and define

F0 (x) = ω(x) ∂ ∂l , Fi (x) = ω(x)F i (x), i = 1, . . . , m.
The fast oscillations (with respect to l) of the control in

dx dt = F0 (x) + m i=1 u i Fi (x)
generate new directions, namely

(ad j F0 ) Fi , j ≥ 0, i = 1, . . . , m,
a natural requirement being that the distribution { F0 , F1 , . . . , Fm } on the augmented space S 1 × X be bracket generating. This is equivalent to the bracket condition on S 1 × X for the distribution {∂/∂l, F 1 , . . . , F m }.

(ii) The assumption of rank constancy is sufficient to get a decomposition into a smooth sum of squares [3]. That the assumption is crucial is illustrated by the fact that it cannot be removed, even in the analytical category. Eigenvalues (and associated projectors) are indeed analytic functions on the set of matrices in the neighbourhoud of a simple (hence diagonalizable) endomorphism [START_REF] Kato | Perturbation theory for linear operators[END_REF]Theorem II.5.16]. Avoiding semi-simple eigenvalues is necessary as is clear considering x 1

x 2 x 2 -x 1 whose eigenvalues are not differentiable at (0, 0). But even in the simple symmetric nonnegative analytic case, existence of a differentiable square root matrix may fail as illustrated by

x 2 1 + x 2 2 0 0 1 . (6)
To get a positive result with non constant rank, one must actually restrict to the symmetric nonnegative case with analytic dependence on one real variable only.

Eigenvalues and eigenvectors are then analytic on the real line [Ibid., Theorem II.6.1 and §II.6.2], and nonnegativeness ensures analyticity of square roots of the eigenvalues.

(iii) In the case of periodic sub-Riemannian systems, the loss of regularity may originate in averaging. The analytical distribution on

S 1 × R 2 F 1 (l, x) = √ 2(x 1 cos l + x 2 sin l) ∂ ∂x 1 , F 2 (l, x) = ∂ ∂x 2 has quadratic form 2(x 1 cos l + x 2 sin l) 2 0 0 1
whose averaged is (6) which does not define an even differentiable sub-Riemannian system on R 2 . (iv) A much stronger requirement is the existence of a change of coordinates on X (inducing a symplectic transformation on the cotangent) so that the averaged quadratic form be diagonal. For a given z ∈ T * X, let

u(l, z) = k∈Z c k (z)e k (l), e k (l) = e ikl ,
denote the Fourier series of the control (4). One has the following convergence result.

Proposition 3.3. For any positive ε, the normal optimal control is

u ε (l) = ε k∈Z c k • z ε (εl) e k (l)
and the series converges pointwisely. Moreover, for any k ∈ Z,

c k • z ε → c k • z uniformly on [0, 1] as ε → 0+
where the sub-Riemannian extremal z depends only on the boundary conditions x 0 , x f on X.

Proof. For l ∈ [0, l f ], u ε (l) = u(l, x ε (εl), ε p ε (εl)) = εu(l, x ε (εl), p ε (εl)) = ε k∈Z c k ( z ε (εl)) e k (l)
thanks to the pointwise convergence of the Fourier series at l for z = z ε (l). For any k ∈ Z,

c k (z) = 1 2π 2π 0 u(l, z)e k (l)dl
and the dependence on z is continuous since the integrand is bounded in a compact neighbourhood of the image of z.

Let us finally denote λ(l, x) the biggest (positive) eigenvalue of the nonnegative quadratic form in p

|u(l, x, p)| 2 = ω 2 (l, x) m i=1 H 2 i (l, z) = 2ω(l, x)H n (l, z),
and define λ(x) = max l∈S 1 λ(l, x). As a supremum of continuous functions, x → λ(x) is only lower semi-continous and we must assume the existence of a continous upper bound in the subsequent result.

Proposition 3.4. Let σ be a continuous function such that σ 2 dominates λ in a neighbourhood of x. Then

lim sup ε→0+ u ε ∞ ε ≤ σ • x • |p| ∞ . Proof. On [0, l f ], one has |u ε (l)| 2 = ε 2 |u(l, z ε (εl))| 2 ≤ ε 2 σ 2 • x ε (εl) • | p ε (εl)| 2 ,
so the result holds since the square root of the right hand side converges uniformly to σ • x • |p| as ε goes to zero by continuity of σ.

Computations in Kepler case

Using the neat geometric coordinates of §3.1 on the three dimensional space of ellipses, we get

H n (l, z) = ω(l, x) 2 (H 2 1 + H 2 2 )(l, z), u(l, z) = ω(l, x)(H 1 , H 2 )(l, z), |u(l, z)| 2 = ω 2 (l, x)(H 2 1 + H 2 2 )(l, z), (7)
with

F 1 (l, x) = P 2 W 2 sin l ∂ ∂e x -cos l ∂ ∂e y , F 2 (l, x) = P 2 W 2 2P W ∂ ∂P +(cos l + e x + cos l W ) ∂ ∂e x + (sin l + e y + sin l W ) ∂ ∂e y ,
and ω(l, x) = W 2 P 3/2 , W = 1 + e x cos l + e y sin l. Introducing mean motion, n = a -3/2 , and true anomaly, τ = l -θ, one gets

F 1 (l, x) = P 2 W 2 - 3en 1 -e 2 ∂ ∂n + sin τ ∂ ∂e -cos τ 1 e ∂ ∂θ , (8)
F 2 (l, x) = P 2 W 2 - 3W n 1 -e 2 ∂ ∂n + (cos τ + e + cos τ W ) ∂ ∂e + (sin τ + sin τ W ) 1 e ∂ ∂θ , (9)
with

W = 1 + e cos τ.
As a result, the computation of Fourier series of H n (or u) are performed with respect to τ rather than l so, as θ only appears through τ ,

H n (l, z) = k∈Z c k (n, e, p)e k (l -θ) = k∈Z c k (z)e k (l),
and coefficients verify c k (z) = c k (n, e, p)e k (θ).

Proposition 3.5. The adjoint p θ is a linear first integral of H.

Proof. According to the previous remark, θ is cyclic in the averaged Hamiltonian.

The remarkable feature of the set (n, e, θ) of coordinates is the following [START_REF] Edelbaum | Optimal low-thrust rendez-vous and station keeping[END_REF][START_REF] Geffroy | Généralisation des techniques de moyennation en contrôle optimal. Application aux problèmes de transfert et rendez-vous orbitaux à poussée faible[END_REF]7].

Proposition 3.6. The averaged Hamiltonian is Riemannian and orthogonal in (n, e, θ) coordinates,

H(z) = 1 2n 5/3 9n 2 p 2 n + 5 2 (1 -e 2 )p 2 e + 5 -4e 2 2 p 2 θ e 2 .
Fourier coefficients of the control are obviously obtained from (8-9), noting that Lemma 3.2. One has

e + cos τ W = -z - 2 √ 1 -e 2 e k≥1 z k cos kτ, sin τ W = - 2 e k≥1 z k sin kτ, where z = -e/(1 + √ 1 -e 2 ) is the only pole in the open unit disk of W = 1 + e cos τ .
We finally provide a continuous upper bound of the eigenvalues of the quadratic form associated with the control norm, allowing us to estimate precisely u ε ∞ as ε → 0+.

Proposition 3.7. Eigenvalues of the quadratic form (7) are uniformly dominated by

σ 2 (n, e) = 4(1 -e 2 ) n 2/3 (1 + e) 2 n 4/3 + 1 + e n 2/3 e + 1 -e 2 .
Proof. Since we have a rank two distribution (8-9) of vector fields (parameterized by l ∈ S 1 ) on the three-dimensional manifold X, a simple computation shows that the maximum eigenvalue is

λ(l, x) = ω 2 (l, x) 2 F 2 1 + F 2 2 + (F 2 1 -F 2 2 ) 2 + 4(F 1 |F 2 ) 2 (l, x),
hence the result.

The resulting estimate provided by Proposition 3.4 depends only on the geodesic connecting the two prescribed points on the manifold. Complete quadrature for these geodesics are computed in the next section.

Analysis of the averaged system

The main step in the analysis is to use further normalizations to obtain a geometric interpretation.

Proposition 3.8. In the elliptic domain, we set r = 2 5 n 5/6 , ϕ = arcsin e, and the metric is isometric to

g = dr 2 + r 2 c 2 (dϕ 2 + G(ϕ)dθ 2 ) with c = 2/5 and G(ϕ) = 5 sin 2 ϕ 1 + 4 cos 2 ϕ •
Geometric interpretation. This normal form captures the main properties of the averaged orbital transfer. Indeed, we extract from g two Riemannian metrics in dimension two g 1 = dr 2 + r 2 dψ 2 with ψ = ϕ/c which is associated with orbital transfer where θ is kept fixed, and also

g 2 = dϕ 2 + G(ϕ)dθ 2
which represents the restriction to r 2 = c 2 .

Analysis of g 1

When p θ vanishes, θ is constant. The corresponding extremals are geodesics of the Riemannian problem in dimension two defined by dθ = 0. We extend the elliptic domain to the meridian half-planes all isometric to

Σ 0 = {n > 0, e ∈] -1, +1[}.
In polar coordinates (r, ψ), Σ 0 is defined by {r > 0, ψ ∈]-π/(2c), π/(2c)[}. This extension allows to go through the singularity corresponding to circular orbits. Geometrically, this describes transfers where the angle of the pericenter is kept fixed and p θ = 0 corresponds to the transversality condition. Such a policy is clearly associated with steering the system towards circular orbits where the angle θ of the pericenter is not prescribed. An important physical subcase is the geostationary final orbit.

In the domain Σ 0 , the metric g 1 = dr 2 + r 2 dψ 2 is a polar metric isometric to the flat metric dx 2 + dz 2 if we set x = r sin ψ and z = r cos ψ. We deduce the following proposition.

Proposition 3.9. The extremals of the averaged coplanar transfer in Σ 0 are straight lines in suitable coordinates, namely x = 2 3/2 5 n 5/6 sin(c -1 arcsin e), z = 2 3/2 5 n 5/6 cos(c -1 arcsin e) with c = 2/5. Since c < 1, the domain is not convex and the metric g 1 is not complete.

Proof. The axis e x = 0 corresponds to circular orbits. Among the extremals, we have two types (see Fig. 1): complete curves of type 1, and non-complete curves of type 2 that meet the boundary of the domain. The domain is not geodesically convex and in subdomain II, the existence theorem fails. For each initial condition, there exists a separatrix S which corresponds to a segment line in the orbital coordinates which is meeting n = 0 in finite time. Its length gives the bound for a sphere to be compact.

In order to complete the analysis of g and to understand the role of g 2 , we present now the integration algorithm.

Integrability of the extremal flow

The integrability property is a consequence of the normal form only,

g = dr 2 + r 2 (dϕ 2 + G(ϕ)dθ 2 ),
and the associated Hamiltonian is decomposed into

H = 1 2 p 2 r + 1 r 2 H ′ , H ′ = 1 2 (p 2 ϕ + p 2 θ G(ϕ)
).

Lemma 3.3. The Hamiltonian vector field -→ H admits three independent first integrals in involution, H, H ′ , p θ , and is Liouville integrable.

To get a complete parameterization, we proceed as follows. We use the (n, e, θ) coordinates and write

H = 1 4n 5/3 [18n 2 p 2 n + H ′′ ]
with

H ′′ = 5(1 -e 2 )p 2 e + 5 -4e 2 e 2 p 2 θ .
Lemma 3.4. Let s = n 5/3 then s(t) is a polynomial of degree 2, s(t) = c 1 t 2 + ṡ(0)t + s(0) with s(0) = n 5/3 (0), ṡ(0) = 15n(0)p n (0) and c 1 = 25H/2.

Lemma 3.5. Let dT = dt/4n 5/3 . If H ′′ (0) = 0, then

T (t) = 1 2 |∆| [arctan L(s)] t 0 where L(t) = (2at + b)/ |∆|, a = c 1 , b = ṡ(0) and ∆ = -25H ′′ (0)/2 is the discriminant of s(t).
This allows to make the integration. Indeed if H ′′ = 0, p e = p θ = 0 and the trajectories are straight lines (the line S in Fig. 1). Otherwise, we observe that n 5/3 (t) is known and depends only upon n(0), p n (0) and H which can be fixed to 1/2 by parameterizing by arc length. Hence, it is sufficient to integrate the flow associated with H ′′ using the parameter dT = dt 4n 5/3 where T is given by the previous lemma. Let H ′′ = c 2 3 and p θ = c 2 . Using p e = ė/10(1 -e 2 ), we obtain

ė2 = 20(1 -e 2 ) e 2 [c 3 e 2 -(5 -4e 2 )c 2 2 ].
To integrate, we set w = 1 -e 2 for e ∈]0, 1[, so the equation takes the form

dw dT = Q(w)
where

Q(w) = 80w[(c 2 3 -c 2 2 ) -(c 2 3 + 4c 2 2 )w] with positive discriminant. Hence the solution is w = 1 2 c 2 3 -c 2 2 c 2 3 + 4c 2 2 1 + sin( 4 √ 5 c 2 3 + 4c 2 2 )T + K ,
K being a constant. We deduce that

θ(T ) = θ(0) + 2c 2 T 0 1 + 4w(s) 1 -w(s) ds
where θ(0) can be set to 0 by symmetry. To conclude, we must integrate (1 + 4w(s))/(1 -w(s)) with w = K 1 (1 + sin x) and x = (4s/ √ 5) c 2 3 + 4c 2 2 + K. Therefore, we must evaluate an integral of the form

A + B sin x C + D sin x dx•
More precisely, the formula is

A + B sin x C + D sin x dx = B D x + AD -BC dx C + D sin x with dx C + D sin x = 2 √ C 2 -D 2 arctan C tan(x/2) + D √ C 2 -D 2 and C 2 -D 2 > 0.
The previous lemmas and computations give Proposition 3.10. For H ′′ = 0, the solutions of -→ H can be computed using elementary functions and

n(t) = 25 2
Ht 2 + 15n(0)p n (0)t + n 5/3 (0)

3/5 , e(t) = 1 -K 1 (1 + sin K 2 (t)), θ(t) = θ(0) + p θ 2|p θ | K 3 -4x + 10 K 3 arctan (1 -K 1 ) tan(x/2) -K 1 K 3 K2(t) K , with K = arcsin 1 -e(0) 2 K 1 -1 , K 1 = 1 2 H ′′ (0) -p 2 θ H ′′ (0) + 4p 2 θ , K 2 (t) = 4 √ 5 T (t) H ′′ (0) + 4p 2 θ + K , K 3 = 5p 2 θ H ′′ (0) + 4p 2 θ •
For H ′′ = 0, they are straight lines.

Remark 3.3. The above formulas give the complete solution of the associated Hamilton-Jacobi Equation.

Geometric properties of g 2

The previous integration algorithm shows that the extremals of this metric describe the evolution of the angular variables θ and ϕ parameterized by dT = dt/r(t) 2 where r(t) 2 is a second order polynomial whose coefficients depend only upon the energy level H fixed to 1/2, r(0) and p r (0). We now give some basic properties of g 2 .

Lemma 3.6. The metric g 2 can be extended to an analytic metric on the whole sphere S 2 , where θ and ϕ are spherical coordinates with two polar singularities at ϕ = 0 or π corresponding to e = 0, whereas the equator corresponds to e = 1; θ is an angle of revolution. The meridians are projections on S 2 of the extremals of g 1 .

Lemma 3.7. The two transformations (ϕ, θ) → (ϕ, -θ) and (ϕ, θ) → (π -ϕ, θ) are isometry of g 2 . This induces the following symmetries for the extremal flow:

-if p θ → -p θ then we have two extremals of same length symmetric with respect to the meridian θ = 0, -if p ϕ → -p ϕ then we have two extremals of same length intersecting on the antipodal parallel, ϕ = π -ϕ(0).

Such properties (illustrated on Fig. 3.4.3) are shared by the following oneparameter family of metrics.

Metrics induced by the flat metric on oblate ellipsoid of revolution.

We consider the flat metric of R 3 , g = dx 2 +dy 2 +dz 2 , restricted to the ellipsoid defined by

x = sin ϕ cos θ, y = sin ϕ sin θ, z = µ cos ϕ where µ ∈]0, 1[. A simple computation leads to

E µ (ϕ)dϕ 2 + sin 2 ϕ dθ 2
for the restricted metric, where E µ (ϕ) = µ 2 + (1 -µ 2 ) cos 2 ϕ, and we can write

g 2 = 1 E µ (ϕ) (E µ (ϕ)dϕ 2 + sin 2 ϕ dθ 2 )
where µ = 1/ √ 5. We deduce the following lemma.

Lemma 3.8. The metric g 2 is conformal to the flat metric restricted to an oblate ellipsoid of revolution with parameter µ = 1/ √ 5.

A global optimality result with application to orbital transfer

In this section, we consider an analytic metric on

R + × S 2 g = dr 2 + (dϕ 2 + G(ϕ)dθ 2 )
and let H be the associated Hamiltonian. We fix the parameterization to arc length by restricting to the level set H = 1/2. Let x 1 , x 2 be two extremal curves starting from the same initial point x 0 and intersecting at some positive t. We get the relations

r 1 (t) = r 2 (t), ϕ 1 (t) = ϕ 2 (t), θ 1 (t) = θ 2 (t),
and from Lemma 3.4, we deduce that Lemma 3.9. Both extremals x 1 and x 2 share the same p r (0) and for each t, r 1 (t) = r 2 (t).

If we consider now the integral curves of H ′ where H = (1/2)p 2 r + H ′ /r 2 on the fixed induced level and parameterize these curves using dT = dt/r 2 , we deduce the following characterization.

Proposition 3.11. The following conditions are necessary and sufficient to characterize extremals of H ′ = 0 intersecting with same length ϕ 1 (T ) = ϕ 2 (T ) and θ 1 (T ) = θ 2 (T ) together with the compatibility condition

T = t 0 dt r 2 (t) = 2 √ ∆ arctan L(t) t t=0 .
Theorem 3.2. A necessary global optimality condition for an analytic metric on R + × S 1 normalized to

g = dr 2 + r 2 (dϕ 2 + G(ϕ)dθ 2 )
is that the injectivity radius be greater than or equal to π on the sphere r = 1, the bound being reached by the flat metric in spherical coordinates.

Proof. We observe that in the flat case, the compatibility condition cannot be satisfied. Moreover, the injectivity radius on S 2 is π corresponding to the halflength of a great circle. For the analytic metric on S 2 under consideration, the injectivity radius is the length of the conjugate point at minimum distance or the half-length of a closed geodesic [START_REF] Do Carmo | Riemannian geometry[END_REF]. The conjugate point is, in addition, a limit point of the separating line. Hence, if the injectivity radius is smaller than π, we have two minimizers for the restriction of the metric on S 2 which intersects with a length smaller than π. We shall show that it corresponds to a projection of two extremals x 1 and x 2 which intersect with same length.

For such extremals r(0) = 1, we set p r (0) = ε, H = 1/2 and we get

2H ′ = p 2 ϕ (0) + p 2 θ (0) G(ϕ(0)) = λ 2 (ε), λ(ε) = 1 -ε 2 .
If t 1 is the injectivity radius on the level set H ′ = 1/2, for H ′ = λ 2 (ε)/2 and p r (0) = ε, it is rescaled to T 1 = t 1 /λ(ε). The compatibility relation for T = T 1 then gives

T 1 = arctan t + ε λ(ε) -arctan ε λ(ε) •
Clearly, the maximum of the right member is π, taking ε < 0, |ε| → 1. Hence, it can be satisfied since t 1 < π. The flat case shows that it is the sharpest bound.

By homogeneity, we deduce the following corollary.

Corollary 3.1. If the metric is normalized to dr 2 + (r 2 /c 2 )(dϕ 2 + G(ϕ)dθ 2 ), then the bound for the injectivity radius on r 2 = c 2 is cπ.

Riemannian curvature and injectivity radius in orbital transfer

Using standard formulae from Riemannian geometry [START_REF] Do Carmo | Riemannian geometry[END_REF], we have the following proposition.

Proposition 3.12. Let g be a smooth metric of the form dr 2 + r 2 (dϕ 2 + G(ϕ)dθ 2 ). The only non-zero component of the Riemann tensor is

R 2323 = r 2 - G ′′ (ϕ) 2 -G(ϕ) + G ′ (ϕ) 2 4G(ϕ)
which takes the form R 2323 = -r 2 F (F ′′ + F ) if we set G(ϕ) = F 2 (ϕ). We have therefore R 2323 = 0 if and only if F (ϕ) = A sin(ϕ + ϕ 0 ) which is induced by the flat case in spherical coordinates.

Hence, the main non-zero sectional curvature of the metric is

K = R 2323 | ∂ ∂θ ∧ ∂ ∂ϕ | 2
and computing this term in the case of orbital transfer, we get:

Lemma 3.10. The sectional curvature in the plane (ϕ, θ) is given by

K = (1 -24 cos 2 ϕ -16 cos 4 ϕ) r 2 (1 + 4 cos 2 ϕ) 2
and K → 0 as r → +∞.

Proposition 3.13. The Gauss curvature of the metric of

g 2 = dϕ 2 + G(ϕ)dθ 2 with G(ϕ) = (sin 2 ϕ)/(1 + 4 cos 2 ϕ) is K = 5(1 -8 cos 2 ϕ) (1 + 4 cos 2 ϕ) 2 •
Theorem 3.3. The Gauss curvature of g 2 is negative near the poles and maximum at the equator. The injectivity radius is π/ √ 5 and is reached by the shortest conjugate point along the equator.

Proof. Clearly K is maximum and equal to five along the equator which is an extremal solution. Hence a direct computation gives that the shortest conjugate point is on the equator with length π/ √ 5. It corresponds to the injectivity radius if the half-length of a shortest periodic extremal is greater than π/ √ 5. Simple closed extremals are computed in [START_REF] Bonnard | Geodesic flow of the averaged controlled Kepler equation[END_REF] using the integrability property and a simple reasoning gives that the shortest corresponds to meridians whose length is 2π. Hence the result is proved.

Corollary 3.2. Since π/ √ 5 < π 2/5, the necessary optimality condition of Theorem 3.3 is not satisfied in orbital transfer for the extension of the metric to R + × S2 .

Cut locus on S 2 and global optimality results in orbital transfer

From the previous section, the computation of the injectivity radius for the metric on S 2 is not sufficient to conclude about global optimality. A more complete analysis is necessary to evaluate the cut locus. This analysis requires numerical simulations. The main results are [START_REF] Bonnard | Geodesic flow of the averaged controlled Kepler equation[END_REF][START_REF] Bonnard | Conjugate and cut loci of a two-sphere of revolution with application to optimal control[END_REF]: Proposition 3.14. For the metric g 2 on S 2 , they are exactly five simple closed geodesics modulo rotations around the poles, the shortest being meridians with length 2π, the longest the equator with length 2π √ 5.

Theorem 3.4. Except for poles, the conjugate locus is a deformation of a standard astroid with axial symmetry and two cusps located on the antipodal parallel.

With the exception of poles, the cut locus is a simple segment, located on the antipodal parallel, with axial symmetry, and whose extremities are cusps points of the conjugate locus. For a pole, the cut locus is reduced to the antipodal pole.

Proof. The proof is made by direct analysis of the extremal curves. The main problem is to prove that the separating line is given by points on the antipodal parallel where, because of the isometry ϕ → π -ϕ, two extremals curves with same length intersect. This property cannot occur before. The results are represented Fig. 3.

Geometric interpretation and comments. The metric is conformal to the restriction of the flat metric to an oblate ellipsoid of revolution. For such a metric, the cut locus is known since Jacobi and is similar to the one represented on Fig. 3. It is a remarkable property that there is no bifurcation of the cut locus when the metric is deformed by the factor E µ (ϕ) In orbital transfer for instance, the Gauss curvature is not positive. On S 2 , relations between the conjugate and cut loci allow to deduce the cut locus from the conjugate locus. 2

The conjugate locus can also easily be computed using the code [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF]. It can also be deduced by inspecting the extremal flow only, the conjugate locus being an envelope. Finally, we observe that in order to have intersecting minimizers, we must cross the equator ϕ = π, that is e = 1. The same is true for conjugate points. Hence we deduce Theorem 3.5. Conjugate loci and separating lines of the averaged Kepler metric are always empty in the spaces of ellipses where e ∈ [0, 1[.

The averaged system in the tangential case

An interesting question is to analyze if the averaged system in the tangential case where the control is oriented along F t retains similar properties [START_REF] Bonnard | Energy minimization of single-input orbit transfer by averaging and continuation[END_REF]. The first step is to compute the corresponding averaged system. Proposition 3.15. If the control is oriented along F t only, the averaged Hamiltonian associated with energy minimization is

H t = 1 2n 5/3 9n 2 p 2 n + 4(1 -e 2 ) 3/2 1 + √ 1 -e 2 p 2 e + 4(1 -e 2 ) 1 + √ 1 -e 2 p 2 θ e 2
and corresponds to the Riemannian metric

g t = dn 2 9n 1/3 + n 5/3 4 1 + √ 1 -e 2 (1 -e 2 ) 3/2 de 2 + 1 + √ 1 -e 2 (1 -e 2 ) e 2 dθ 2
where (n, e, θ) remain orthogonal coordinates.

Construction of the normal form

We proceed as in Section 3.4 and set r = 2 5 n 5/6 , e = sin ϕ 1 + cos 2 ϕ .

The metric becomes

g t = dr 2 + r 2 c 2 (dϕ 2 + G(ϕ)dθ 2 ), c = 2 5 < 1, and
G(ϕ) = sin 2 ϕ 1 -(1/2) sin 2 ϕ 1 -sin 2 ϕ 2 .
Hence the normal form is similar to the full control case. As before, we introduce the metrics

g 1 = dr 2 + r 2 dψ 2 , ψ = ϕ/c, and
g 2 = dϕ 2 + G(ϕ)dθ 2 .
The main difference with the full control case will concern the singularities of the function G.

Metrics g 1 and g 2

The metric g 1 corresponds again to transfer to circular orbits and is the polar form of the flat metric dx 2 + dz 2 , if x = r sin ψ and z = r cos ψ.

The normal form reveals the same homogeneity property between the full control and the tangential case, so the metric g 2 can be used to make a similar optimality analysis, evaluating the conjugate and cut locus. But the metric g 2 cannot be interpreted as a smooth metric on S 2 . This can be seen by computing the Gauss curvature.

Proposition 3.16. The Gauss curvature of g 2 is given by

K = (3 + cos 2 ϕ)(cos 2 ϕ -2) (1 + cos 2 ϕ) cos 2 ϕ •
In particular K → -∞ when ϕ → π/2, and since K < 0, the conjugate locus of any point is empty.

Nevertheless, the extremals can be smoothly extended through the singular boundary of the domain, the equator ϕ = π/2.

Integration of the extremal flow

The algorithm based on the normal form is similar to the bi-input case, but we compare the respective transcendence. The Hamiltonian is written

H = 1 4n 5/3 [18n 2 p 2 n + H ′′ ]
where H ′′ now takes the form

H ′′ = 8(1 -e 2 ) 3/2 1 + √ 1 -e 2 p 2 e + 8(1 -e 2 ) 1 + √ 1 -e 2 p 2 θ e 2 •
We set H ′′ = c 2 3 , p θ = c 2 , and from

p e = 4n 5/3 (1 + √ 1 -e 2 )e 16(1 -e 2 ) 3/2 we obtain dw dT 2 = Q(w) (1 + w) 2
where w = √ 1 -e 2 , T is as in the bi-input case, and Q is the fourth-order polynomial Hence, the integration requires the computation of the elliptic integral

Q(w) = 32w[c 2 3 (1 -w 2 )(1 + w) -8c 2 2 w 2 ].
dw(1 + w) Q(w)
which has an additional complexity. It is related to the pole of order 2 of the metric at the equator. See [START_REF] Bonnard | Singular metrics on the two-sphere in space mechanics[END_REF] for both aspects.

Conclusion in both cases

The previous analysis shows that the full control case and the tangential one admit a uniform representation in coordinates (θ, ϕ). In particular, it allows to make a continuation between the respective Hamiltonians, i.e., between the respective functions G. A correction has to be made between orbit elements e which are respectively defined by e = sin ϕ versus e = sin ϕ 1 + cos 2 ϕ.

The flows in the two cases are presented on Fig. 4 and reveal the similar structure. The conjugate locus is reached after having crossed the equator. On Fig. 5 and 6 we present a first continuation result, in the tangential case, showing the convergence of the method from the averaged to the non averaged trajectory for specific transfer orbits.

Conjugate and cut loci on a two-sphere of revolution

The problem of computing the conjugate and cut loci in orbital transfer is connected to a very old geometric problem which goes back to Jacobi and is briefly introduced next.

Definition 3.1. The two-sphere S 2 endowed with a smooth metric of the form dϕ 2 + G(ϕ)dθ 2 in spherical coordinates is called two-sphere of revolution.

Many of them can be realized as Riemannian surfaces of revolution embedded in R 3 by rotating on a smooth curve homeomorphic to a half-circle. The classical examples are the following: The averaged trajectories are clearly nice approximations of the optimal ones of the original system. Hence, convergence of the underlying shooting method to compute the non-averaged minimizing trajectory is easily obtained. -Round sphere S 2 . It is constructed restricting the Euclidian metric to S 2 and given by dϕ 2 + sin 2 (ϕ)dθ 2 .

-Oblate ellipsoid of revolution O(µ). If we restrict the Euclidian metric to the surface x = sin ϕ cos θ, y = sin ϕ sin θ, z = µ cos ϕ, with µ < 1, the metric is (1 -(1 -µ 2 ) sin 2 ϕ)dϕ 2 + sin 2 (ϕ)dθ 2 which can be set in the form dϕ 2 + G(ϕ)dθ 2 using a quadrature.

We recall some basic properties on the ellipsoid of revolution.

Proposition 3.17. On an oblate ellipsoid of revolution, -the Gauss curvature is monotone increasing from the North pole to the equator, -the cut locus of a point which is not a pole is a subarc of the antipodal parallel,

-the conjugate locus of a point which is not a pole has a standard astroid shape with four cusps.

The simple structure of the cut locus is a consequence of [START_REF] Sinclair | The cut locus of a two-sphere of revolution and Toponogov's comparison theorem[END_REF].

Theorem 3.6. Let dϕ 2 + G(ϕ)dθ 2 be a metric on a two-sphere of revolution.

We assume:

-The transformation ϕ -→ π -ϕ is an isometry i.e. G(π -ϕ) = G(ϕ),

-The Gauss curvature K is monotone non decreasing along a meridian from the North pole to the equator.

Then, the cut locus of a point not a pole is a simple branch located on the antipodal parallel.

Application. Let g λ be the family of analytic metrics on S 2 defined by

g λ = dϕ 2 + G λ (ϕ)dθ 2 , G λ (ϕ) = (1 + λ) 2 sin 2 ϕ (1 + λ cos 2 ϕ) , λ ≥ 0.
The Gauss curvature is given by

K λ = - 1 √ G ∂ 2 √ G ∂ 2 ϕ = (1 + λ)(1 -2λ cos 2 ϕ) (1 + λ cos 2 ϕ) 2 •
Hence if 0 < λ ≤ 2, then K λ is monotone non decreasing from the North pole to the equator and the previous theorem asserts that the one parameter family has a cut locus reduced to a simple branch for λ ∈]0, 2].

If λ > 2 the Gaussian curvature K λ is not monotone and the result cannot be applied. In particular the orbit transfer with full control corresponds to λ = 4, while at the limit case λ = +∞ a singularity appears. The Riemannian metric has a pole at the equator, the situation being similar to the one occuring when the thrust is only tangential. Hence Theorem 3.6 has to be refined to deal with such situations. The final result is coming from [START_REF] Bonnard | Singular metrics on the two-sphere in space mechanics[END_REF].

We consider a metric of the form g = ϕ 2 + G(ϕ)dθ 2 where G ′ is non zero on ]0, π 2 [ and G(π -ϕ) = G(ϕ). Definition 3.2. The first return mapping to the equator is the map

∆θ : p θ ∈]0, √ G(π/2)[ → ∆θ(p θ )
that is the θ-variation of the extremal parameterized by arc length and associated with the adjoint vector component p θ . The extremal flow is called tame if the first return mapping is monotone non-increasing for p θ ∈]0, √ G(π/2)[. Theorem 3.7. In the tame case, the cut locus of a point different from a pole is a subset of the antipodal parallel. If moreover ∆θ ′ (p θ ) < 0 < ∆θ ′′ (p θ ) on ]0, √ G(π/2)[, then the conjugate locus of such a point has exactly four cusps.

Remark 3.4. This result can be extended to the singular case where the metric has poles at the equator.

4 Energy minimization in the Earth-Moon space mission with low thrust

The N-body problem

In this section, we follow mainly [START_REF] Meyer | Introduction to Hamiltonian dynamical systems and the N -body problem[END_REF] (see also [START_REF] Pollard | Mathematical introduction to celestial mechanics[END_REF] and [START_REF] Szebehely | Theory of orbits: The restricted problem of three bodies[END_REF]). Consider N point masses m 1 , . . . , m N moving in a Galilean reference system R 3 where the only forces acting are their mutual attractions. If q = (q 1 , . . . , q N ) ∈ R 3N is the state and p = (p 1 , . . . , p N ) is the momentum vector, the equations of motion are

q = ∂H ∂p , ṗ = - ∂H ∂q
where the Hamiltonian is

H = N i=1 |p i | 2 2m i -U, U (q) = N 1≤i<j≤N Gm i m j |q i -q j | •
A subcase is the coplanar situation where the N masses are in a plane R 2 . In this case the Galilean reference frame can be replaced by a rotating frame defined by

K = 0 1 -1 0 , exp(ωtK) = cos ωt sin ωt -sin ωt cos ωt ,
and introducing a set of coordinates which uniformly rotates with frequency ω, one defines the symplectic transformation

u i = exp(ωtK)q i , v i = exp(ωtK).
A standard computation gives the Hamiltonian of the N -body problem in rotating coordinates,

H = N i=1 |v| 2 2m i - N i=1 w t u i Kv i - N 1≤i<j≤N Gm i m j |q i -q j | •
In particular, the Kepler problem in rotating coordinates up to a normalization has the following Hamiltonian

H = |p| 2 2 -t qKp - 1 |q| •
-Euler points. They are the collinear points denoted L 1 , L 2 and L 3 located on the line y = 0 defined by the primaries. For the Earth-Moon problem they are given by

x 1 ≃ -1.0051, x 2 ≃ 8.369e -1, x 3 ≃ 1.1557.
-Lagrange points. The two points L 4 and L 5 form with the two primaries an equilateral triangle.

Some important information about stability of the equilibrium points are provided by the eigenvalues of the linearized system. The linearized matrix evaluated at points L 1 , L 2 or L 3 admits two real eigenvalues, one being strictly positive, and two imaginary ones. The collinear points are consequently not stable. In particular, the eigenvalues of the linearized matrix evaluated at L 2 with µ = 1.2153e -2 are approximately ±2.931837 and ±2.334248i. When it is evaluated at L 4 or L 5 , the linearized matrix has two imaginary eigenvalues since µ < µ 1 = (1/2)(1 -√ 69/9) in the Earth-Moon system. Points L 4 and L 5 are thus stable according to Arnold stability theorem [START_REF] Meyer | Introduction to Hamiltonian dynamical systems and the N -body problem[END_REF]. See also [START_REF] Marsden | New methods in celestial mechanics and mission design[END_REF] for a review of mission design techniques using the equilibrium points.

The continuation method

The mathematical continuation method in the restricted circular problem is omnipresent in Poincaré's work, in particular for the continuation of circular orbits [START_REF] Poincaré | [END_REF]. Geometrically, it is simply a continuation of trajectories of Kepler problem into trajectories of the 3-body problem. It amounts to considering µ as a small parameter-the limit case µ = 0 being Kepler problem in the rotating frame-writing

H 0 = |p| 2 2 -t qKp - 1 |q| + o(µ).
The approximation for µ is valid, a neighbourhood of the primaries being excluded. In the Earth-Moon problem, since µ is very small, the Kepler problem is a good approximation of the motion in a large neighbourhood of the Earth. This point of view is important in our analysis, as indicated by the status report of the SMART-1 mission since most of the mission time is spent under the influence of the Earth attraction only, see [START_REF] Racca | SMART-1: The first time of Europe to the Moon[END_REF][START_REF] Racca | SMART-1 mission description and development status[END_REF].

The control problem

The control system in the rotating frame is deduced from the previous model and can be written in Hamiltonian form

dx dt = -→ H 0 (x) + u 1 -→ H 1 (x) + u 2 -→ H 2 (x)
where x = (q, p), -→ H 0 is the free motion and -→ H 1 , -→ H 2 are given by -→ H i = -q i , i =1,2. As for the Kepler problem, the mass variation of the satellite can be introduced in the model dividing u i by m(t) and adding the equation ṁ = -δ|u|. Again, it will be not taken into account here. Moreover we still restrict our analysis to the energy minimization problem with fixed final time t f and control valued in R 2 . The physical problem which is to maximize the final mass can be analyzed using a continuation method.

A lunar mission using low propulsion called SMART-1 was realized by ESA and the practical details of the mission-in particular the description of the trajectory-are reported in [START_REF] Racca | SMART-1 mission description and development status[END_REF]. We present a trajectory analysis based on our geometric and numerical techniques. For simplicity, we have fixed the boundary conditions to circular orbits, the one around the Earth corresponding to the geostationary one. But everything can be applied to other boundary conditions, for instance those described in the report status of the SMART-1 mission. A trajectory comparison is discussed in the final section (see also [6]).

Our analysis is based on a numerical continuation, including second-order optimality check, where µ is the parameter of the continuation. The averaged system is finally applied to get an approximate energy minimizing trajectory for the phase of the mission that starts from the circular Earth orbit and aims at a quasi-elliptic orbit where the apogee is about 338000 Kilometers. In the restricted problem approximation, the effect of the inclination (around 30 degrees in phase two of the mission) is not modelled here.

Numerical continuation for the Earth-L 2 transfer

As a first approach we choose to simulate the Earth-L 2 transfer in the restricted 3-body problem. Indeed, in the limit case µ = 0, the Moon and the point L 2 are identical. Moreover, in the Earth-Moon system, the point L 2 and the Moon are very close. As a result, the first phase of an Earth-Moon transfer is comparable to an Earth-L 2 transfer. Solving the shooting function associated with the Earth-L 2 transfer is consequently useful to provide a good approximation of the solution of the Earth-Moon transfer shooting function.

Using a circular orbit around the Earth for the geostationary one, we set as initial condition x 0 = (1 -µ + 1.099e -1, 0, 0, 2.8792). The point L 2 is located on (x L2 , 0) with x L2 solution of the equation

x - (1 -µ)(x + µ) ρ 3 1 - µ(x -1 + µ) ρ 3 2 = 0.
Since we want to reach L 2 with a zero speed, we fix the target x f = (x L2 , 0, 0, 0). By making the parameter µ vary from zero, one builds up a family (S µ ) µ of shooting functions which connects the Kepler and the 3-body problem. The numerical continuation method can be applied to deduce low thrust extremal trajectories of the Earth-L 2 transfer from the Kepler ones.

In accordance with the report status of ESA, we fix the transfer time to 121 time units of the restricted 3-body problem, which approximately corresponds to the transfer time from the Earth to the point L 2 during the SMART-1 mission (about 17 months). In addition, we consider a constant spacecraft mass of 350 Kilograms, see [START_REF] Racca | SMART-1: The first time of Europe to the Moon[END_REF][START_REF] Racca | SMART-1 mission description and development status[END_REF]. Setting µ = 0, we compute an extremal using the shooting method, then increase µ up to 1.2153e-2 with a discrete continuation. At each step, the first conjugate time along the extremal is computed to ensure convergence of the continuation method. The Euclidian norm of the extremal control is plotted Fig. 14 to draw a comparison between the control bound and the maximum thrust allowed by electro-ionic engines. Figs. 7 to 14 present the computed spacecraft trajectories in both rotating and fixed frames, as well as the first conjugate time and the norm of the control along trajectories in Kepler and Earth-Moon systems. The numerical continuation method, considering µ as a small parameter, leads to deduce an extremal trajectory of the energy minimization Earth-L 2 transfer problem from one corresponding to the Kepler case, in accordance with Poincaré's work. The second order optimality condition check ensure that the computed extremals are locally energy minimizing in L ∞ ([0, t f ]). We also note that in both cases µ = 0 and µ = 1.2153e -2, the maximum value reached by the norm of the extremal control is inferior to the bound |u| ≤ 8e -2 of the SMART-1 electro-ionic engine (7.3e -2 Newtons), while transfer times are comparable.

Numerical continuation method for the Earth-Moon transfer

The second part of our trajectory analysis is devoted to the Earth-Moon transfer. We use the same dynamics and initial conditions as previously. In this case, the circular orbit around the Moon, denoted O L and defined by (x 1 -1+µ) 2 +x 2 2 = 1.7e-3, x 2 3 +x 2 4 = 2.946e-1, (x 1 -1+µ, x 2 ) ⊥ (x 3 , x 4 ), is chosen as the target orbit. In addition to the former necessary conditions, the maximum principle provides the transversality condition p(t f ) ⊥ T x(t f ) O L . The shooting equation is modified accordingly, and one has to check local optimality computing focal instead of conjugate points. The transfer time is fixed to 124 time units of the restricted 3-body problem. The extremal trajectory is obtained using the shooting method and initializing p 0 with the initial costate vector of the Earth-L 2 transfer. The Earth-Moon trajectories in both rotating and fixed frames, the first focal time and the norm of the extremal control are presented from Fig. 15 to Fig. 22 for µ varying between 0 and 1.2153e -2. In the two cases µ = 0 and µ = 1.2153e -2, the first focal time along extremals t 1,foc is higher than 3t f /2, ensuring local optimality. The maximal bound of the norm of the extremal control is about 4.5e -2 which approximately corresponds to the half of the maximal thrust allowed during SMART-1 mission.

It is interesting to notice that the Earth-L 2 Keplerian trajectory differs from the Earth-Moon Keplerian trajectory. The first phase of the Earth-Moon transfer matches the Earth-L 2 transfer. It underlines the crucial role of the neighbourhood of the point L 2 where the attractions of the two primaries compensate each other. By treating the first phase as a Keplerian transfer from the geostationary orbit (GEO) to a geostationary transfer orbit (GTO), we are able to use the averaging results of section 3.2 in order to estimate the maximal bound of the control associated with energy minimizing trajectories.

The first phase of the Earth-Moon trajectory with µ = 1.2153e -2 can be approximated by a Keplerian transfer from the GEO to a GTO orbit with semimajor axis a ≃ 5.84e -1, eccentricity e ≃ 3.96e -1, and argument of pericenter θ = 8π/7. Using the system of coordinates introduced in Proposition 3.9, one can compute the appropriate geodesic and estimate the final longitude required to achieve a prescribed maximal bound on the norm of the control thanks to Proposition 3.7. [2] Arnold, V. I. Mathematical methods of classical mechanics. Springer, New-York, 1989.

[3] Bellaïche, A. The tangent space in sub-Riemannian geometry. Sub-Riemannian geometry, 1-78, Progr. Math. 144, Birkhäuser, 1996.

Figure 1 :

 1 Figure 1: Geodesics of the metric g 1 in (n, e x ) and flat coordinates.

Figure 2 :

 2 Figure 2: Action of the symmetry group on the extremals.

Figure 3 :

 3 Figure 3: Conjugate and cut loci in averaged orbital transfer.

Figure 4 :

 4 Figure 4: Extremal flow of g 2 in the full control and tangential cases, in the (ϕ, θ) coordinates, starting from ϕ = π/6.

Figure 5 :

 5 Figure 5: Computation by continuation of the non-averaged solution.The averaged trajectories are clearly nice approximations of the optimal ones of the original system. Hence, convergence of the underlying shooting method to compute the non-averaged minimizing trajectory is easily obtained.

Figure 6 :

 6 Figure6: Trajectory of a non-averaged solution for ε = 1e-2 with (e(0), n(0)) = (7.5e-1, 5e-1) and (e(t f ), n(t f )) = (5e-2, 3e-1). Dashed ellipses are averaged ellipses and provide a good approximation of the motion.

Figure 7 :

 7 Figure 7: Earth-L 2 trajectory in the rotating frame, µ = 0.

Figure 8 :

 8 Figure 8: Earth-L 2 trajectory in the fixed frame, µ = 0.

Figure 9 :

 9 Figure 9: First conjugate time, Earth-L 2 transfer, µ = 0.

Figure 10 :

 10 Figure 10: Norm of the extremal control, Earth-L 2 transfer, µ = 0.

Figure 11 :

 11 Figure 11: Earth-L 2 trajectory in the rotating frame, µ = 1.2153 -2.

Figure 12 :

 12 Figure 12: Earth-L 2 trajectory in the fixed frame, µ = 1.2153e -2.

Figure 13 :

 13 Figure 13: First conjugate time, Earth-L 2 transfer, µ = 1.2153e -2.

Figure 14 :

 14 Figure 14: Norm of the extremal control, Earth-L 2 transfer, µ = 1.2153e -2.

Figure 15 :Figure 16 :Figure 17 :

 151617 Figure 15: Earth-Moon trajectory in the rotating frame, µ = 0.

Figure 18 :Figure 19 :

 1819 Figure 18: Norm of extremal control, Earth-Moon transfer, µ = 0.

[ 4 ]

 4 Bensoussan, A. Regular perturbations in optimal control. Singular perturbations in systems and control (Udine, 1982), 169-183, CISM Courses and Lectures, 280, Springer, Vienna, 1983. [5] Bliss, G. A. Lectures on the Calculus of Variations. Chicago, 1946.

Figure 20 :Figure 21 :

 2021 Figure 20: Earth-Moon trajectory in the fixed frame, µ = 1.2153e -2.

[ 6 ]

 6 Bombrun, A.; Chetboun, J.; Pomet, J.-B. Transfert Terre-Lune en poussée faible par contrôle feedback -La mission SMART-1. INRIA Research report (2006), no. 5955, 1-27.

[ 7 ]

 7 Bonnard, B.; Caillau, J.-B. Riemannian metric of the averaged energy minimization problem in orbital transfer with low thrust. Ann. Inst. H. Poincaré

Topological dual of the space of smooth compactly supported functions on the real line.

For instance, a domain bounded by two intersecting minimizing curves must contain a conjugate point.

The circular restricted 3-body problem in Jacobi coordinates

The following representation of the Earth-Moon problem fits in the so-called circular restricted 3-body problem. In the rotating frame, the Earth which is the biggest primary planet with mass 1 -µ is located at (-µ, 0) while the Moon with mass µ is located at (1 -µ, 0) (the small parameter being µ ≃ 1.2153e -2).

We note z = x + iy the position of the spacecraft, and ρ 1 , ρ 2 the distances to the primaries,

The equation of motion is

where -V is the potential of the system defined by

The system can be written using Hamiltonian formalism setting q 1 = x, q 2 = y, p 1 = ẋ -y, q 2 = ẏ + x, and the Hamiltonian describing the motion writes

Jacobi integral, Hill regions and equilibrium points

The Jacobi integral using Hamiltonian formalism is simply the Hamiltonian H 0 which gives

where

The Hill domain for the value h is the region where the motion can occur, that is {(x, y) ∈ R 2 | Ω(x, y) + h ≥ 0}.

The equilibrium points of the problem are well known. They split into two different types:


M Lesur

E Gravier

K Lim

C Djerroud

M Idouakass

IMPURITY TRANSPORT IN COLLISIONLESS TRAPPED-PARTICLE-DRIVEN TURBULENCE

The paper reports on properties of the radial turbulent transport of light to heavy impurities due to trapped-particledriven turbulence in the core of tokamaks, based on a gyrokinetic bounce-averaged model. The focus is on the effects of impurity concentration, charge, mass, gradientson impurity transport. The validity of the usual passive treatment of impurities is found to be limited to concentrations of W 40+ tungsten (respectively, C 6+ carbon) below 2x10 -4 (1.4%). The diffusive impurity transport decreases slightly with increasing mass number. It depends strongly on the charge number Z, but the effect depends on the nature of turbulence: it increases with Z for turbulence dominated by Trapped Electron Modes (TEM) but decreases with Z for turbulence dominated by Trapped Ion Modes (TIM). The direction of thermo-diffusion also depends on the nature of the turbulence, and on the impurity temperature profile. Finally, the direction of the curvature pinch depends on magnetic shear: it is inward for a positive magnetic shear, while a negative magnetic shear prevents impurity core accumulation, in agreement with quasilinear theory.

INTRODUCTION

In the context of magnetic confinement fusion plasmas, ion species other than hydrogen isotopes are referred to as impurities. Tokamak plasmas often contain impurities such as helium, tungsten, beryllium, carbon, argon, or neon. Accumulation of impurities in the core threatens the viability of fusion. One issue is the dilution of fuel, which degrades the efficiency of fusion reactors. Another urgent issue is radiative loss, which can be particularly strong for heavy impurities. All in all, the efficiency of fusion is very sensitive to core impurity concentration [START_REF] Pütterich | Calculation and experimental test of the cooling factor of tungsten[END_REF]. Improving theories of impurity transport is crucial to address these urgent issues.

To account for both kinetic and global effects, to which turbulent transport is sensitive, we adopt a global gyrokinetic approach. However, since we base the present exploration on many parameter scans, which would require hundreds of millions of hours of parallel computing if performed with a global 5D gyrokinetic code, we adopt a reduced approach. Namely, we start from a 3D bounce-averaged gyrokinetic model, focusing on particles trapped deeply by the magnetic field inhomogeneity in toroidal geometry. The orbits of these particles live in a 2D phase-space parametrized by an energy invariant. This approach is meant to explore the basic mechanisms of turbulent transport in collisionless trapped-particledriven turbulence, rather than to provide quantitative predictions since this in some sense we isolate a part of the complex system of core tokamak turbulence. Here we investigate the impacts of impurity concentration, charge, mass and gradients.

MODEL

In conventional 5D gyrokinetics, particle motion is described in a 4D phase-space (the 3 coordinates of the gyrocenter + parallel velocity or momentum) parameterized by the gyrokinetic counterpart of the magnetic moment. In contrast, here, we limit the kinetic description to the bounce-averaged dynamics of deeply trapped IAEA-CN-286/1160 particles. Such a particle is parametrized by its precession angle α, its average poloidal magnetic flux 𝜓 (a radial coordinate), and an invariant 𝐸 akin to a kinetic energy. The distribution function f s (α, 𝜓, 𝐸, 𝑡) of these particles, for each species 𝑠, is evolved according to a Vlasov equation,

∂f s ∂t + [J 0 ϕ, f s ] α,ψ + ω d,s ∂f s ∂α = 0
where [J 0 ϕ, f s ] α,ψ are Poisson brackets, J 0 is a gyrokinetic-gyrobounce operator, ω d,s is the precession frequency, which depends on both 𝐸 and the charge number Z.

The other types of particles in the plasma are described via their density response to electric potential fluctuations ϕ, which forms the left-hand-side of the quasineutrality equation,

∑ eZ s 2 n eq,s T eq,s s [ 1 -f t f t (ϕ -ϵ ϕ,s ⟨ϕ⟩) -Δ ̅ s ϕ] = 4π√2 ∑ Z s m s 3/2 s ∫ J 0,s f s E 1/2 dE ∞ 0
Here f t is the fraction of trapped particles, ϵ ϕ,s ∈ [0,1] controls the response of passing particles. Finally, Δ ̅ s ϕ is a polarization term, where the Laplacian operator depends on the Larmor radius ρ 𝑐 and the banana width δ b .

The semi-Lagrangian simulation code TERESA has been developed based on this model. The role of kinetic effects was discussed in Refs. [START_REF] Lesur | A simple model for electron dissipation in trapped ion turbulence[END_REF][3] [START_REF] Medina | Test particle dynamics in low-frequency tokamak turbulence[END_REF]. TERESA was applied in 2018 to a deuterium plasma with non-trace concentrations of tungsten, showing that peaked or hollow impurity density profile can change the nature of turbulence from TEM to TIM [START_REF] Idouakass | Impurity density gradient influence on trapped particle modes[END_REF]. This paper reports on more recent works based on TERESA simulations with a deuterium plasma and various impurities.

LIMIT OF VALIDITY OF THE PASSIVE TREATMENT

In gyrokinetic simulations, trace impurity species can be either treated self-consistently (what we refer to as active treatment), or, on the contrary, without including them in Maxwell equations (what we refer to as passive treatment). The latter approach is consistent with the assumption that impurities with low enough concentrations are impacted by turbulence generated by electrons and main ions, but do not impact it significantly in return. Indeed, the relative contribution of impurity density in quasi-neutrality is of the order of C.Z 2 , where C is the impurity concentration, and Z its charge number. In the "trace" limit (C.Z 2 << 1), the impurity behaves as a passive species, or test-particle, which means that its presence does not affect the turbulent state.

In this section, which is based on Ref. [START_REF] Lesur | Validity limits of the passive treatment of impurities in gyrokinetic tokamak simulations[END_REF], we relax this assumption, and investigate the impacts of impurity, on impurity transport, as a function of its concentration, in the presence of trapped-particle-driven turbulence. We run a series of simulations with self-consistent (active) treatment of impurities for a wide range of concentrations and compare the results with a series of simulations with passive treatment of impurities.

Impact of tungsten concentration

For each species, the equilibrium distribution function is chosen as 2D Maxwellian with radial density profile such that κ n = |𝑑𝑛/𝑑𝜓| =const. and similarly for the temperature profile with gradient 𝜅 𝑇 . The impacts of tungsten concentration depend on the relationship between equilibrium density gradient and temperature gradient [START_REF] Lesur | Validity limits of the passive treatment of impurities in gyrokinetic tokamak simulations[END_REF]. Here we focus on a case where these gradients are equal. The gradients and the simulation parameters are given in Table 1. Thermal baths are imposed at both radial boundaries. 1 summarizes the impact of W 40+ concentration. It includes the relative variation, due to the presence of impurity, of linear growth rate, turbulence intensity, and impurity density flux. For example, the relative variation of linear growth rate is (𝛾 -𝛾 0 )/𝛾 0 , where 𝛾 0 is the growth rate when C=0. A similar definition is used for turbulence intensity. As for the impurity density flux, it is normalised to C, and the variation is with respect to the impurity flux in the trace limit. We note 𝜔 0 the precession frequency of particles with velocity equal to the Deuterium thermal speed. The measured impact of impurity concentration on linear growth rate 𝛾 for C<10 -3 agrees with the analytic expression 𝛾 = 𝛾 0 -𝐶 𝑑𝛾 with 𝑑𝛾 = 2.2 × 10 -4 𝜔 0 , which was obtained at first order in C in the appendix 1 of Ref. [START_REF] Lesur | Validity limits of the passive treatment of impurities in gyrokinetic tokamak simulations[END_REF]. The stabilizing effect of impurity on TEM is consistent with Refs. [START_REF] Du | Coupling of ion temperature gradient and trapped electron modes in the presence of impurities in tokamak plasmas[END_REF][8], although our interpretation differs: we have shown that this is mainly due to a modification of the response of passing particles when 𝐶. 𝑍 2 ∼ 1. The measured impact of impurity concentration on turbulence intensity is in qualitative agreement with a mixing length estimate, 𝑒ϕ/𝑇 0 = γ/(𝑘 𝑟 ρ 𝑐,𝑖 𝑘 θ 𝑐 𝑠 ). However, the rms of the potential is twice more sensitive to C than predicted from the latter estimate.

Furthermore, the sensitivity of impurity density flux to C is even more threshold-like than expected from the sensitivity of turbulence intensity. As observed on FIG. 1, there is a transition for all linear growth rate, turbulence intensity, and impurity density and heat fluxes, centred around 𝐶. 𝑍 2 = 1.1. However, the transition is significantly steeper for the fluxes. This transition may be viewed as a threshold, which gives the range of validity of the passive treatment, C < 2.10 -4 . When the temperature gradient is finite and the density profile is flat, we obtain qualitatively similar results, except for the density flux which vanishes. For the heat flux, the range of validity of the passive treatment is slightly smaller than for the case with equal gradients, C < 10 -4 .

FIG. 1. Relative variation, due to the presence of impurity, of linear growth-rate, turbulence intensity (rms of the potential), and density flux, all versus tungsten concentration C.

Here the impurity flux is normalised to the impurity density, and the variation is relative to the trace limit.

Impact of carbon concentration

We have performed a similar analysis for C 6+ carbon. There are two qualitative differences compared with the tungsten case:

-The transition for the flux is centred around C.Z 2 = 0.5; -The transition for the flux is at lower concentration than the transition for turbulence intensity and growthrate, which is opposite to the tungsten case.

However, our main conclusion stands: the transition is steeper for the flux than for turbulence, and steeper for turbulence than for the growth-rate.

Synchronisation between impurity density and electric potential fluctuations

We have investigated the physical mechanism responsible for this difference in behaviour between turbulence intensity and transport. We observe that phase-synchronization between impurity density fluctuations and electric potential fluctuations occurs for high enough impurity concentrations, which quenches impurity transport.

IAEA-CN-286/1160

DIFFUSIVE TRANSPORT, THERMO-DIFFUSIVE PINCH, AND CURVATURE-DRIVEN PINCH

Hereafter, the impurity species are treated self-consistently but in the trace limit, so that impurity concentrations do not affect the nature of the turbulence.

Impact of charge and mass numbers on diffusion

To evaluate impurity transport, we plot the cumulative flux along the 𝜓-direction (meaning along the radial direction) given by

Γ cumul = ∫ Γ(t * )dt * t 0 with Γ(𝑡) = - 2 √π ∫ 𝑑α ∫ 𝑓 ∂(𝒥ϕ) ∂α 𝐸 1/2 ∞ 0 dE 2π 0
The impurity flux due to Trapped Electron Mode (TEM) turbulence increases with (A,Z) (FIG. 2). In contrast, the impurity flux is found to decrease with (A,Z) when the Trapped Ion Mode (TIM) dominated (FIG. 3). In this case the higher charge impurity has smaller impurity flux. In fact, we have shown [START_REF] Gravier | Diffusive impurity transport driven by trapped particle turbulence in tokamak plasmas[END_REF] that the dependency on mass number A is much weaker than the dependency on the charge number Z.

FIG. 2. Impurity cumulative flux plotted against time in the case of a TEM turbulence.

Helium diffuses from the core to the edge while other impurities diffuse from the edge to the core, according to the initial density gradients.

FIG. 3. Impurity cumulative flux plotted against time in the case of a TIM turbulence. Helium diffuses from the core to the edge while other impurities diffuse from the edge to the core, according to the initial density gradients.

Let us compare the theoretically predicted impurity flux and the nonlinear impurity transport simulations. Quasilinear theory yields:

∂n 0 ∂t + ∂Γ ψ ∂ψ = 0
with Γ 𝜓 the impurity transport in the radial direction:

Γ ψ = -κ n ∫ ∑|n𝒥ϕ n | 2 γ n (ω r,n -n Ω d Z E) 2 + γ n 2 F eq ψ=0 E 1/2 dE n E
Both 𝛾 𝑛 and 𝜔 𝑟,𝑛 have very weak dependency on both 𝑍 and 𝐴, since only trace impurities are considered, and therefore turbulence is governed by electrons and main ions. Therefore, from this last equation we can see that the radial impurity flux depends on the charge number 𝑍 as expected but depends weakly on the mass number. Indeed, the impurity mass only impacts the gyro-bounce-averaging operator, where the Larmor radius and the banana width are proportional to √A/Z. But this 𝐴-dependence of the gyro-bounce-averaging operator is weak and does not depend on the sign of 𝜔 𝑟 and the nature of the instability. Nevertheless, regarding the 𝐴dependence, the impurity flux is expected to decrease as 𝐴 increases, in both TEM and TIM cases. These results agree with the numerical simulations. This expression can be further evaluated by considering two different cases: negative and positive frequencies. The real frequencies 𝜔 𝑟,𝑛 and the linear growth rates 𝛾 𝑛 are needed and determined by solving the linear dispersion relation. The saturation level of the plasma potential is given by nonlinear simulations.

Here, we choose five different impurities with the same mass number (A = 20) but with 5 different charge numbers Z = 2, 4, 8, 10 and 12. We describe how impurity diffusive transport depends on the charge number, depending on the nature of the dominant instabilities: the diffusion coefficient increases with Z in the case of TEM turbulence [START_REF] Gravier | Diffusive impurity transport driven by trapped particle turbulence in tokamak plasmas[END_REF] (FIG. 4), while it decreases with Z in TIM-dominated turbulence, in good agreement with quasilinear theory.

FIG. 4. Impurity diffusion coefficient DTEM plotted against the impurity charge number Z in the case of TEM turbulence.

IAEA-CN-286/1160

Thermo-diffusive pinch in TIM and TEM turbulences

The total particle fluxconsisting of the diffusive and convective partcan be decomposed into

Γ 𝑧 𝑁 𝑧 = -𝐷 𝑧 ∇𝑁 𝑧 𝑁 𝑧 + 𝑉 𝑧 = +𝐷 𝑧 𝜅 𝑛𝑧 ⏟ 𝐷𝑖𝑓𝑓𝑢𝑠𝑖𝑜𝑛 + 𝑉 𝑧 ⏟ 𝐶𝑜𝑛𝑣𝑒𝑐𝑡𝑖𝑜𝑛 with 𝑉 𝑧 ∝ 𝐶 𝑇 𝜅 𝑇𝑧 + 𝐶 𝑃
The first term in the pinch velocity 𝑉 𝑧 corresponds to the thermo-diffusive pinch which stems from the inherent temperature gradient while the second term is the curvature pinch, proportional to magnetic field curvature.

Here, we use five different values of impurity temperature profiles to investigate the thermo-diffusive pinch with no magnetic curvature (𝐶 𝑃 = 0). FIG. 5 shows that, as expected, the pinch velocity 𝑉 𝑧 disappears when the temperature profile is flat. 𝑉 𝑧 changes its sign when the temperature gradient does. We observe that the sign of the coefficient 𝐶 𝑇 changes depending on the TIM or TEM nature of the turbulence at a given temperature gradient [START_REF] Lim | Impurity pinch generated by trapped particle driven turbulence[END_REF]. Here, 𝜅 𝑇𝑧 > 0 means that the core is hotter than the edge, therefore the right part of the figure is more relevant. In this case, curvature pinch is inward except for reversed magnetic shear, consistently with [START_REF] Futatani | Reversal of impurity pinch velocity in tokamaks plasma with a reversed magnetic shear configuration[END_REF].

FIG. 5. Thermodiffusion velocity VZ plotted against the impurity temperature gradient, for both TEM and TIM turbulences.

Impact of magnetic shear on curvature pinch

In a similar way that thermo-diffusive pinch was studied by imposing 𝐶 𝑃 = 0, the direction of magnetic curvature pinch can be investigated by imposing a flat temperature gradient for the impurity (𝜅 𝑇𝑧 = 0). Given that the Maxwellian distribution function 𝐹 𝑒𝑞 contains the effect of magnetic curvature through the Hamiltonian 𝐻 𝑒𝑞 = In the bounce-averaged gyrokinetic code, the Hamiltonian can be expressed in terms of toroidal precession frequency Ω 𝑑 which is a characteristic frequency of trapped particles.

𝐻 𝑒𝑞 = 𝐸(1 + 𝑒αΩ 𝑑 𝜓)

where α is an ad-hoc coefficient which determines the positive/negative magnetic shear according to its sign.

A positive value of α represents a positive magnetic shear while a negative value of α mimics a negative magnetic shear. The direction of the curvature pinch VZ with different values of α is represented in FIG. 6. Our simulations show that the direction of the impurity curvature pinch is inward in the case of a positive magnetic shear, while a negative magnetic shear switches the sign of VZ and therefore prevents impurity core accumulation. This is in qualitative agreement with quasilinear theory.

FIG. 6. Magnetic curvature pinch VZ plotted with different values of α (with flat temperature gradient of impurity).

The sign of α is representative of the sign of magnetic shear.

SUMMARY, CAVEATS AND PERSPECTIVES

To summarize, we first analysed the impact of impurity concentration on mode growth, saturation amplitude, and transport in TEM/TIM turbulence. There is a transition for all growth rate, turbulence intensity, and impurity density and heat fluxes, centered around C.Z 2 = 1.1 for W 40+ and 0.5 for C 6+ . However, the transition is significantly steeper for the fluxes due to phase-synchronization between electric potential fluctuations and impurity density fluctuations for high enough concentrations. The range of validity of the passive treatment of tungsten is limited to concentrations smaller than 10 -4 (in order of magnitude).

Secondly, the nonlinear simulations demonstrate that, in qualitative agreement with quasi-linear theory :

-Impurity diffusive transport decreases (increases) with Z in TIM (TEM) turbulence. The dependency on mass number is much weaker; -For standard impurity temperature gradient, thermo-diffusion brings impurities outwards (inwards) in TIM (TEM) turbulence. -Curvature pinch is inward except for reversed magnetic shear.

The model used in the above analyses is collisionless and limited to trapped-particle-driven turbulence. We applied it as a prototype of turbulent transport in toroidal fusion plasmas to draw qualitative trends. Although this is outside the scope of this paper, we have recently made progress using more advanced models. Firstly, ITG modes must be included for the quasi-steady-state turbulent spectrum to be more accurate. Secondly, in general, both collisional (neoclassical) and turbulent processes contribute to the radial transport of impurities, in some conditions in synergy [START_REF] Esteve | Self-consistent gyrokinetic modeling of neoclassical and turbulent impurity transport[END_REF][13] [START_REF] Idomura | Dynamics of enhanced neoclassical particle transport of tracer impurity ions in ion temperature gradient driven turbulence[END_REF]. It turns out that an accurate collision operator is crucial. Thirdly, sources and boundary conditions can have significant impacts. The new version of the GYSELA-X gyrokinetic code is well-adapted in these respects. Recently we have been applying it to investigate impurity transport [START_REF] Lim | Gyrokinetic modelling of light to heavy impurity transport in tokamaks[END_REF].

TABLE 1

 1

	. SIMULATION PARAMETERS
		𝑎/𝑅 0		f t	κ n = 𝜅 𝑇
		0.1		2/3	0.25
	Species	Z	ρ 𝑐 /𝑎	δ b /a	ϵ ϕ,s	Concentration	𝑇 𝑒 /𝑇 𝑖
	Electrons	-1	0.01	0.01	1	1	1
	Deuterium	1	0.03	0.1	0	1 -40 C	1
	Tungsten	40	0.0072	0.024	0	C	1

ACKNOWLEDGEMENTS

This work was carried out within the framework of the French Federation for Magnetic Fusion Studies (FR-FCM) and of the EUROfusion Consortium. It received funding from the Agence Nationale de la Recherche for the project GRANUL (ANR-19-CE30-0005), and from the Euratom research and training programmes 2014-2018 and 2019-2020 under Grant Agreement No. 633053 for the project WP17-ENRCEA-02. The views and opinions expressed herein do not necessarily reflect those of the European Commission. This work was granted access to the HPC resources of IDRIS under allocations 2019 and 2020 made by GENCI, of EXPLOR under Project No. 2017M4XXX0251, and of CINECA Marconi under projects GSNTIT and GSNTITS.