Bernard Bonnard
email: bernard.bonnard@u-bourgogne.fr.

Olivier Cots
email: olivier.cots@irit.fr.

Boris Wembe
email: boris.wembe@irit.fr.

Keywords: 49K15, 53C60, 70H05 Helhmoltz-Kirchhoff N vortices model, Zermelo navigation problem, Geometric optimal control, Conjugate and cut loci, Clairaut-Randers metric with polar singularity

come

A Zermelo navigation problem with a vortex singularity

A ZERMELO NAVIGATION PROBLEM WITH A VORTEX SINGULARITY * Bernard Bonnard 1 , Olivier Cots 2 and Boris Wembe 3 Abstract. Helhmoltz-Kirchhoff equations of motions of vortices of an incompressible fluid in the plane define a dynamics with singularities and this leads to a Zermelo navigation problem describing the ship travel in such a field where the control is the heading angle. Considering one vortex, we define a time minimization problem which can be analyzed with the technics of geometric optimal control combined with numerical simulations, the geometric frame being the extension of Randers metrics in the punctured plane, with rotational symmetry. Candidates as minimizers are parameterized thanks to the Pontryagin Maximum Principle as extremal solutions of a Hamiltonian vector field. We analyze the time minimal solution to transfer the ship between two points where during the transfer the ship can be either in a strong current region in the vicinity of the vortex or in a weak current region. The analysis is based on a micro-local classification of the extremals using mainly the integrability properties of the dynamics due to the rotational symmetry. The discussion is complex and related to the existence of an isolated extremal (Reeb) circle due to the vortex singularity. The explicit computation of cut points where the extremal curves cease to be optimal is given and the spheres are described in the case where at the initial point the current is weak.

Résumé. Les équations d'Helhmoltz-Kirchhoff pour le mouvement tourbillonaire d'un fluide incompressible dans le plan définissent une dynamique hamiltonienne à singularités localisées aux tourbillons. Cela conduit à définir un problème de Zermelo décrivant le mouvement d'un navire où le contrôle est l'angle de cap et le critère à minimiser est le temps de transfert entre deux points du plan. Dans cet article, on se limite au cas d'un seul tourbillon localisé en zéro et le problème est analysé avec les techniques du contrôle optimal géométrique combinées à des simulations numériques, le contexte géométrique étant l'extension des métriques de Randers dans le plan possédant une symétrie de révolution. On doit en effet considérer le cas d'un courant faible, mais aussi d'un courant fort localisé au voisinage du tourbillon. Les trajectoires candidates à minimiser le temps sont paramétrées en utilisant le Principe du Maximum de Pontryaguine comme des extrémales solutions d'un système hamiltonien dont les projections sur l'espace d'état sont les géodésiques. L'analyse du problème optimal repose sur la classification microlocale des solutions extrémales utilisant l'intégrabilité de la dynamique. La discussion est complexe et repose sur l'existence d'un cercle géodésique dit de Reeb, conséquence de la singularité tourbillonaire. La discussion est complétée par l'évaluation des points de coupure, les points où les extrémales cessent d'être optimales. Les sphères sont décrites dans le cas d'un état initial à courant faible.

Introduction

Helhmoltz and Kirchhoff originated the model of the displacement of particles in a two dimensional fluid, see [START_REF] Helmholtz | On integrals of hydrodynamics equations, corresponding to vortex motions[END_REF][START_REF] Kirchhoff | Vorlesungen uber mathematische Physik[END_REF] for the original articles and [START_REF] Arnold | Topological Methods in Hydrodynamics[END_REF][START_REF] Newton | The N-Vortex Problem -Analytical Techniques[END_REF][START_REF] Saffman | Vortex dynamics[END_REF] for a modern presentation of Hamiltonian dynamics. In this model, the vorticity of the fluid is concentrated at points z i := (x i , y i ), i = 1, • • • , N , with circulation parameters k i and the configuration space is R 2N with coordinates (x 1 , y 1 , • • • , x N , y N ) endowed with the symplectic form ω := N i=1 k i dy i ∧ dx i . The dynamics is given by the Hamiltonian canonical equations

k i ẋi = ∂H ∂y i , k i ẏi = - ∂H ∂x i , (1)
1 ≤ i ≤ N , where the Hamiltonian function H is

H := - 1 π i<j k i k j ln z i -z j , (2)
where z i -z j is the Euclidean distance. In this article, we consider a motionless single vortex given by eq. ( 1), that we fix at the origin of the reference frame. This corresponds to set z 1 = (x 1 , y 1 ) at (0, 0). We consider a particle as a (point) vortex with zero circulation, setting k 2 = 0 with N = 2, under the influence of the current generated by the vortex and given by the vector field with a singularity at the origin. The current is defined by [START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF] and denoted, omitting indices, F 0 (x, y) := X 1 (x, y) ∂ ∂x + X 2 (x, y) ∂ ∂y , with z := (x, y). These classical notations being not adapted for later considerations, we will denote by x := (x 1 , x 2 ) the position of the particle (instead of z = (x, y)) and by k the circulation parameter of the vortex. In the following we will consider a (Zermelo) time minimization problem.

To define a Zermelo navigation problem, following [START_REF] Carathéodory | Calculus of Variations and Partial Differential Equations of the First Order, Part 1, Part 2[END_REF][START_REF] Zermelo | Über das Navigations problem bei ruhender oder veränderlicher wind-verteilung[END_REF] and see [START_REF] Bryson | Applied optimal control[END_REF] for the optimal control frame related to Zermelo's problems, we consider the particle as the ship of the navigation problem and the control is defined by α the heading angle of the ship axis. Hence the control field is given by u := u max (cos α, sin α) where u max is the maximal amplitude and this leads to a control system written as:

dx dt = F 0 (x) + u 1 F 1 (x) + u 2 F 2 (x), (3)
with F 1 := ∂/∂x 1 , F 2 := ∂/∂x 2 , x = (x 1 , x 2 ) and u := (u 1 , u 2 ) bounded by u ≤ u max . We then consider the associated time minimal control problem to transfer the ship from an initial configuration x 0 to a target x f , where x 0 and x f are two points of the punctured plane R 2 \ {0}. By a rescaling we can assume u max = 1 and denoting by g the Euclidean metric on the plane, u ≤ 1 bounds the control amplitude by 1 and we have two cases: the case F 0 g < 1 of weak current versus the case F 0 g > 1 of strong current. In the weak case, the time minimal problem defines a Randers metric in the plane, which is a specific Finsler metric, see [START_REF] Bao | Zermelo navigation on Riemannian manifolds[END_REF] for this geometric frame. In the neighborhood of the vortex we have F 0 g > 1, hence, due to the singularity we have a non trivial extension of the classical case.

A neat treatment of the historical Zermelo navigation was made by [START_REF] Carathéodory | Calculus of Variations and Partial Differential Equations of the First Order, Part 1, Part 2[END_REF][START_REF] Zermelo | Über das Navigations problem bei ruhender oder veränderlicher wind-verteilung[END_REF] and their study is an important inspiration for our work. Optimal control, with the Hamiltonian formulation coming from the Pontryagin Maximum Principle [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes, Translated from the Russian[END_REF], forms the frame that we shall use in our analysis, combined with recent development concerning Hamiltonian dynamics to deal with N vortices or N bodies dynamics, see [START_REF] Meyer | Introduction to Hamiltonian Dynamical Systems and the N-Body Problem[END_REF]. An intense research activity was lead by H. Poincaré on the dynamics of such systems [31] to compute periodic trajectories avoiding collisions and such techniques lead to the concept of choregraphy developed by [START_REF] Chenciner | Simple Choreographic Motions of N Bodies: A Preliminary Study[END_REF] for the N -body problem and [START_REF] Calleja | Choreographies in the n-vortex problem[END_REF] for the N -vortex system, showing the relations between both dynamics in the Hamiltonian frame [START_REF] Meyer | Introduction to Hamiltonian Dynamical Systems and the N-Body Problem[END_REF]. From the control point of view, there is a lot of development related to space navigation for the N -body problem, see [START_REF] Bonnard | Mécanique céleste et contrôle des véhicules spatiaux[END_REF], valuable in our study for ship navigation in the N -vortex problem. Optimal control problem in this area was developed for this purpose in relation with Finsler geometry, see [START_REF] Bonnard | Optimal Control with Applications in Space and Quantum Dynamics[END_REF] or [START_REF] Serres | Géométrie et classification par feedback des systèmes de contrôle non linéaires de basse dimension[END_REF] for a general setting in the planar case. More general results about vortex control may be found in [START_REF] Protas | Vortex dynamics models in flow control problems[END_REF][START_REF] Vainchtein | Optimal control of a co-rotating vortex pair: averaging and impulsive control[END_REF] for instance.

From the time minimal point of view, using the Maximum Principle, we lift the control dynamics (3) by defining the pseudo-Hamiltonian H(x, p, u) := H 0 (x, p)

+ 2 i=1 u i H i (x, p)
where H i (x, p) := p • F i (x), i = 1, 2, 3, are the Hamiltonian lifts of F 0 , F 1 and F 2 and the Maximum Principle leads to analyze the extremal curves solution of the Hamiltonian vector field

#- H := ∂H ∂p ∂ ∂x - ∂H ∂x ∂ ∂p ,
defined by the maximized Hamiltonian H(x, p) := max u ≤1 H(x, p, u). In the context of geometric optimal control, see [START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF][START_REF] Jurdjevic | Geometric Control Theory[END_REF] for a general reference, time minimal solutions are obtained by a micro-local analysis of #-H combined with the computation of the cut point along extremal curves, that is the first point where an extremal curve ceases to be globally optimal. Fixing the initial point x(0) = x 0 of the navigation problem, the set of such points is called the cut locus. Even in the Riemannian geometry, the determination of the cut locus is a very complicated problem and there are only a few results. A major recent contribution concerns the case of ellipsoids solving the Jacobi conjecture [START_REF] Itoh | The cut loci and the conjugate loci on ellipsoids[END_REF]. Parallel development were obtained in the frame of space navigation where geometric analysis is combined with numerical methods, see again [START_REF] Bonnard | Optimal Control with Applications in Space and Quantum Dynamics[END_REF] for a general reference for such contributions. In the frame of Zermelo navigation problem with a small current, called Randers metrics, some results were obtained recently for sphere of revolutions [START_REF] Hama | The cut locus of a Randers rotational 2-sphere of revolution[END_REF].

Our aim is to extend those results for the navigation problem, with a single vortex, which combines different new phenomena in particular the existence of a singularity localized at the vortex position which leads to the strong current case and the need of extending the Finsler case [START_REF] Bao | An Introduction to Riemann-Finsler Geometry[END_REF][START_REF] Bao | Zermelo navigation on Riemannian manifolds[END_REF]. The ultimate goal is to analyze the regularity of the value function: x f → V (x 0 , x f , µ) where V (x 0 , x f , µ) is the minimum time from x 0 to reach the point x f of the punctured plane, in presence of a vortex with a circulation k := 2πµ. The parameter µ is introduced later for practical convenience. The regularity of x f → V (x 0 , x f , µ) is analyzed, in particular, in relation with Legendrian and Lagrangian singularities [START_REF] Goryunov | Lagrangian and Legendrian Singularities[END_REF] associated to the Hamiltonian dynamics #-H. Regularity of the value function in relation with the Hamilton-Jacobi-Bellman equation leads to sufficient global optimality conditions, see the seminal reference [START_REF] Boltyanskii | Sufficient conditions for optimality and the justification of the dynamic programming method[END_REF].

The organization of this article is the following. In Section 2, we present the existence theorem to transfer in minimum time any two points of the punctured plane. We state the Maximum Principle to parameterize the minimizers as extremal curves of a smooth Hamiltonian vector field. This leads to define a shooting method used to compute candidates as minimizers. Extremal curves are classified using generic assumptions into hyperbolic and abnormal curves, candidate as time minimizing curves and elliptic curves candidate as time maximizers. Conjugate points where an extremal curve ceases to be optimal for the C 1 -topology are calculated numerically and leads to the conclusion of the absence of such points, hence, the optimality problem boils down to compute cut points in the case of an empty conjugate locus. The final Section 3 is the main contribution of this article with the existence result. Thanks to the integrability of the extremal flow due to the rotational symmetry, the micro-local classification of extremals is presented. The two important phenomena is the existence of abnormal minimizers and of a single extremal circle trajectory called a Reeb circle. Using this classification, the cut points can be computed along any extremal to determine the time minimal value function, combining geometric analysis and numerical simulations using the HamPath code. We present in details the case where at the initial point the current is weak. This gives a nontrivial extension of the Finsler situation.

Existence results and Pontryagin Maximum Principle

Existence of time minimal solutions

We consider a single vortex centered in the reference frame and thus the control system of our Zermelo navigation problem is given by 2 the square distance of the ship, that is the particle, to the origin and where k is the circulation of the vortex. This control system may be written in the following form:

ẋ1 (t) = - k 2π x 2 (t) r(t) 2 + u 1 (t), ẋ2 (t) = k 2π x 1 (t) r(t) 2 + u 2 (t), with r(t) 2 := x 1 (t) 2 + x 2 (t)
ẋ(t) = F 0 (x(t)) + 2 i=1 u i (t)F i (x(t)), (4)
with F 0 , F 1 and F 2 three real analytic (i.e. C ω ) vector fields, where the current (or drift) is given by

F 0 (x) := µ x 2 1 + x 2 2 Å -x 2 ∂ ∂x 1 + x 1 ∂ ∂x 2 ã , (5)
with µ := k/2π, and where the control fields are F 1 := ∂/∂x 1 and F 2 := ∂/∂x 2 . Considering the polar coordinates (x 1 , x 2 ) =: (r cos θ, r sin θ) and an adapted rotating frame for the control, v := u e -iθ , the control system (4) becomes

ṙ(t) = v 1 (t), θ(t) = µ r(t) 2 + v 2 (t) r(t) . (6)
We give hereinafter some classical definitions and refer to [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF] for more details. We consider admissible control laws in the set U := {u : [0 , +∞) → U | u measurable} , where the control domain U := B(0, u max ) ⊂ R 2 denotes the Euclidean closed ball of radius u max > 0 centered at the origin. Since the drift introduces a singularity at the origin, we define by M := R 2 \ {0} the state space, and for any u ∈ U and x 0 ∈ M , we denote by x u (•, x 0 ) the unique solution of (4) associated to the control u such that x u (0, x 0 ) = x 0 . We introduce for a time T > 0 and an initial condition x 0 ∈ M , the set U T,x0 ⊂ U of control laws u ∈ U such that the associated trajectory x u (•, x 0 ) is well defined over [0 , T ], and we denote by A T,x0 := Im E T,x0 the atteignable set (or reachable set) from x 0 in time T , where we have introduced the endpoint mapping

E T,x0 : U T,x0 -→ M u -→ x u (T, x 0 ).
Then, we denote by A x0 := ∪ T ≥0 A T,x0 the atteignable set from x 0 . We recall that the control system is said to be controllable from x 0 if A x0 = M and controllable if A x0 = M for any x 0 ∈ M . Now, for a given pair (x 0 , x f ) ∈ M 2 and some parameters u max ∈ R * + and µ ∈ R, we define the problem of steering (4) in minimum time from the initial condition x 0 to the target x f :

(P ) V (x 0 , x f , µ, u max ) := inf {T | (T, u) ∈ D x0 and x u (T, x 0 ) = x f } ,
where

D x0 := {(T, u) ∈ [0 , +∞) × U | u ∈ U T,x0 }.
We emphasize the fact that the value function V depends on the initial condition x 0 , the target x f and the parameters u max and µ. The first main result is the following:

Theorem 2.1. For any (x 0 , x f , µ, u max ) ∈ M 2 × R * × R * + , the problem (P ) admits a solution. Remark 2.2. Note that when µ = 0, the result is clearly true in R 2 but false in M = R 2 \ {0}.

Up to a time reparameterization τ := t u max and a rescaling of µ, one can fix u max = 1 and we have the relation V (x 0 , x f , µ, u max ) = V (x 0 , x f , µ/u max , 1)/u max . We thus fix from now u max = 1 and write the value function (with a slight abuse of notation)

V (x 0 , x f , µ) := V (x 0 , x f , µ, 1). ( 7
)
We first prove that there exists an admissible trajectory connecting any pair of points in M .

Lemma 2.3. The system (4) is controllable.

Proof. Let consider a pair (x 0 , x f ) ∈ M 2 . We introduce r 0 := x 0 and r f := x f . From x 0 , we can apply a constant control v(t) = (±1, 0) (depending on whether r f is smaller or greater than r 0 ) until the distance r f is reached and then apply a constant control v(t) = (0, sign(µ)) until the target x f is reached.

Remark 2.4. The controllability gives us that the value function is finite while the existence of solutions implies that the value function is lower semi-continuous.

The existence of time-optimal solutions relies on the classical Filippov's theorem [16, Theorem 9.2.i] and the main idea is to prove that the problem (P ) is equivalent to the same problem with the restriction that the trajectories remain in a compact set. To prove this, we need a couple of lemmas. Let us introduce some notations for the first lemma: for a trajectory-control pair (x, u), we associate the pair (q, v) with q := (r, θ) the polar coordinates and v = u e -iθ . We denote by q v (•, q 0 ) the solution of ( 6) with control v such that q v (0, q 0 ) = q 0 . We define for (ε, R, µ) ∈ R + × R * + × R * and θ 0 ∈ R, two optimization problems: (a) The minimum time to make a complete round at a distance R to the vortex:

T θ (R, θ 0 , µ) := inf {T | (T, u) ∈ D x0 and q v (T, (R, θ 0 )) = (R, θ 0 + s 2π)} ,
where s := sign(µ), x 0 := (R cos θ 0 , R sin θ 0 ) and where the control u is related to v by u = ve iθ . (b) The minimum time to reach the circle of radius ε from a distance R to the vortex:

T r (ε, R, θ 0 , µ) := inf {T | (T, u) ∈ D x0 and r v (T, (R, θ 0 )) = ε} .
Since it is clear, due to the rotational symmetry of the problem, that T θ (R, •, µ) and T r (ε, R, •, µ) are invariant with respect to θ 0 , one can fix θ 0 = 0 and set T θ (R, µ) := T θ (R, 0, µ) and T r (ε, R, µ) := T r (ε, R, 0, µ). Besides, from the proof of lemma 2.5, one can notice that T r does not depend on µ, hence, one can define T r (ε, R) := T r (ε, R, 0). Under these considerations, we have the following comparison between T θ and T r :

Lemma 2.5. For any (ε, R, µ) s.t. µ = 0, 0 < R < R µ and 0 ≤ ε < ε µ,R , with R µ := |µ| 2π -1 and ε µ,R := R Å 1 - 2πR |µ| + R ã , we have 0 ≤ ε < ε µ,R < R < R µ and T θ (R, µ) < T r (ε, R),
that is the minimum time to make a complete round at a distance R to the vortex is strictly smaller than the minimum time to reach the circle of radius ε < R.

Proof. It is clear from (6) that T θ (R, µ) is given by the control v(t) = (0, s). This gives by a simple calculation

T θ (R, µ) = 2πR 2 /(|µ| + R). It is also clear that T r (ε, R, µ) is given by v(t) = (-1, 0), whence T r (ε, R, µ) = R -ε = T r (ε, R
) and indeed T r does not depend on µ. Fixing ε = 0, we have T θ (R, µ) = T r (0, R) ⇔ R = |µ|/(2π -1) =: R µ . Besides, we have

T θ (R, µ) < T r (ε, R) ⇔ ε < R Å 1 - 2πR |µ| + R ã =: ε µ,R
but also we have 0 < ε µ,R ⇔ R < R µ , whence the conclusion.

Next an admissible trajectory x associated to a pair (T, u) ∈ D x0 is such that x(T ) = x u (T, x 0 ) = x f . Let us fix (x 0 , x f , µ) ∈ M 2 × R * and introduce r 0 := x 0 and r f := x f . Then: Lemma 2.6. There exists ε > 0 such that any optimal trajectory is contained in M \ B(0, ε).

Proof. Let consider (ε, R) s.t. 0 < R < min{R µ , r 0 , r f } and 0 < ε < ε µ,R . Let us recall that ε < ε µ,R < R since R < R µ and consider an admissible trajectory x intersecting B(0, ε) and associated to a pair denoted (T, u). Then, there exists two times 0

< t 1 ≤ t 2 < T s.t. x([t 1 , t 2 ]) ⊂ B(0, ε). Since 0 < ε < R < min{r 0 , r f }, there exists also 0 < t in < t 1 ≤ t 2 < t out < T s.t.
x([t in , t out ]) ⊂ B(0, R) and s.t. x(t in ) and x(t out ) belong to ∂B(0, R) = S(0, R), the sphere of radius R centered at the origin. In all generality, one can assume that ∀ t ∈ [0 , t in ) ∪ (t out , T ], x(t) / ∈ B(0, R). Let consider the circular arc from x(t in ) to x(t out ) obtained with a control v = (0, s), s := sign(µ), realized in a time denoted τ > 0. It is clear that τ ≤ T θ (R, µ) since T θ (R, µ) is the time to make a circular arc of angle 2π. It is also clear from lemma 2.5 and from the definition of T r that τ ≤ T θ (R, µ) < T r (ε, R) ≤ t out -t in . Let us replace the part x([t in , t out ]) by the circular arc. Then, the new trajectory associated to the pair denoted (T , u ) is still admissible and is strictly better than x since T = T -(t out -t in ) + τ < T . Moreover, this new trajectory is by construction contained in M \ B(0, ε). Whence the conclusion.

® ® ® ε R x f x 0 ® ® ® ε R x f x 0 Figure 1.
Illustration of the construction of a strictly better admissible trajectory. The vortex is represented by a red ball, while the trajectories are the solid black lines. One can see on the left, a trajectory crossing the ball of radius ε. This trajectory is replaced on the right subgraph by a strictly better admissible trajectory.

We are now in position to prove Theorem 2.1.

Proof of Theorem 2.1. By lemma 2.3, there exists an admissible trajectory x. Let T * denote the first time s.t.

x(T * ) = x f . Let us introduce R 1 := ε from lemma 2.6 and R 2 := r 0 + T * , with r 0 := x 0 . By lemma 2.6, the problem (P ) is equivalent to the same problem with the additional constraint R 1 ≤ r(t). Since ṙ(t) = v 1 (t) and v 1 (t) ≤ 1, then for any t ∈ [0 , T * ] we have r(t) ≤ r 0 + T * . The problem (P ) is thus equivalent to the same problem with the additional constraints: R 1 ≤ r(t) ≤ R 2 . The trajectories of the equivalent problem are contained in the compact set B(0, R 2 ) \ B(0, R 1 ). The result follows from the Filippov's existence theorem.

Classification of the extremal curves

In this section, we recall concepts and results from [START_REF] Bonnard | Optimal Control with Applications in Space and Quantum Dynamics[END_REF]. Let x 0 ∈ M and (T, u) ∈ D x0 be an optimal solution of problem (P ) with x := x u (•, x 0 ) the associated optimal trajectory. According to the Pontryagin Maximum Principle [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes, Translated from the Russian[END_REF], there exists an absolutely continuous function p : [0 , T ] → R 2 satisfying the adjoint equation almost everywhere over [0 , T ]:

ṗ(t) = -∇ x H(x(t), p(t), u(t)), (8)
where1 H(x, p, u)

:= p • (F 0 (x) + u 1 F 1 (x) + u 2 F 2 (x)
) is the pseudo-Hamiltonian associated to the problem (P ). Besides, there exists p 0 ≤ 0 such that: the pair (p(•), p 0 ) never vanishes [START_REF] Bonnard | Mécanique céleste et contrôle des véhicules spatiaux[END_REF] and such that the optimal control satisfies the maximization condition almost everywhere over [0 , T ]:

H(x(t), p(t), u(t)) = max w∈U H(x(t), p(t), w) = -p 0 . ( 10
)
Definition 2.7. An extremal is a 4-uplet (x(•), p(•), p 0 , u(•)) satisfying ( 4) and ( 8)- [START_REF] Bonnard | Théorie des singularités de l'application entrée/sortie et optimalité des trajectoires singulières dans le problème du temps minimal[END_REF]. It is said abnormal whenever p 0 = 0 and normal whenever p 0 = 0. It is called strict if p(•) is unique up to a factor. An extremal

(x(•), p(•), p 0 , u(•)) is called a BC-extremal if x(0) = x 0 and if there is a time T ≥ 0 s.t. x(T ) = x f .
Let us introduce the Hamiltonian lifts H i (x, p) := p • F i (x), i = 0, 1, 2, the function Φ := (H 1 , H 2 ) and the switching function ϕ defined for t ∈ [0 , T ] by ϕ(t) := Φ(z(t)) = p(t), z(•) := (x(•), p(•)). The maximization condition [START_REF] Bonnard | Théorie des singularités de l'application entrée/sortie et optimalité des trajectoires singulières dans le problème du temps minimal[END_REF] implies for a.e. t ∈ [0 , T ]:

u(t) = ϕ(t) ϕ(t) = p(t) p(t) , whenever ϕ(t) = 0. Introducing the switching surface Σ := z ∈ M × R 2 Φ(z) = 0 = M × {0} and denoting z := (x, p) ∈ M × R 2
, one can define outside Σ the Hamiltonian:

H(z) := H(z, Φ(z)/ Φ(z) ) = H 0 (z) + Φ(z) = H 0 (z) + p . (11)
Definition 2.8. An extremal (x(•), p(•), p 0 , u(•)) contained outside the switching surface Σ is called of order zero.

Let us recall that a switching time 0 < t < T is a time s.t. ϕ(t) = 0 and s.t. for any ε > 0 (small enough) there exists a time τ ∈ (t -ε, t + ε) ⊂ [0 , T ] s.t. ϕ(τ ) = 0. We can show that the extremals are only of order zero, and thus are smooth: Proposition 2.9. All the extremals are of order zero, that is there are no switching times. Proof. Let (x(•), p(•), p 0 , u(•)) be an extremal. If there exists a time t s.t. ϕ(t) = 0, then p(t) = 0 and we have H(x(t), p(t), u(t)) = 0 = -p 0 which is impossible by [START_REF] Bonnard | Mécanique céleste et contrôle des véhicules spatiaux[END_REF].

We have the standard following result: Proposition 2.10. The extremals of order zero are smooth responses to smooth controls on the boundary of u ≤ 1. They are singularities of the endpoint mapping E T,x0 for the L ∞ -topology when u is restricted to the unit sphere S 1 .

For any Hamiltonian H(z), resp. pseudo-Hamiltonian H(z, u), we denote by #-

H(z) := (∇ p H(z), -∇ x H(z)), resp.
#-H(z, u) := (∇ p H(z, u), -∇ x H(z, u)), its associated Hamiltonian vector field, resp. pseudo-Hamiltonian vector field. With these notations, we have the following classical but still remarkable fact: Proposition 2.11. Let (x(•), p(•), p 0 , u(•)) be an extremal. Denoting z := (x, p), then, we have over [0 , T ]:

ż(t) = #- H(z(t), u(t)) = #- H(z(t)) = # - H 0 (z(t)) + Å p(t) p(t) , 0 ã , (12)
that is the extremals are given by the flow of the Hamiltonian vector field associated to the maximized Hamiltonian H.

Proof. Since the extremal is of order zero, the control t → u(t) is smooth and the adjoint equation ( 8) is satisfied all over [0 , T ]. Besides, denoting (with a slight abuse of notation) u(z) := Φ(z)/ Φ(z) , we have:

H (z) = ∂H ∂z (z, u(z)) + ∂H ∂u (z, u(z)) • u (z) = ∂H ∂z (z, u(z)) + Φ(z) T Å I 2 Φ(z) - Φ(z)Φ(z) T Φ(z) 3 ã =0 • Φ (z) = ∂H ∂z (z, u(z)) = H 0 (z) + Å 0, p p ã .
This proposition shows the importance of the true Hamiltonian H which encodes all the information we need and gives a more geometrical point of view: we will thus refer to trajectories as geodesics. Besides, from the maximum principle, optimal extremals are contained in the level set {H ≥ 0}. Let z(•, x 0 , p 0 ) := (x(•, x 0 , p 0 ), p(•, x 0 , p 0 )) be a reference extremal curve solution of ż = #-H(z) with initial condition z(0, x 0 , p 0 ) = (x 0 , p 0 ) and defined over the time interval [0 , T ]. Lemma 2.12. One has x(t, x 0 , λp 0 ) = x(t, x 0 , p 0 ) and p(t, x 0 , λp 0 ) = λp(t, x 0 , p 0 ). Thanks to this lemma, and since p never vanishes, we can fix by homogeneity p 0 = 1. We thus introduce the following definition that gives us a way to parameterize the extremals of order zero. Definition 2.13. We define the exponential mapping by

exp x0 (t, p 0 ) := Π • e t #- H (x 0 , p 0 ), (13)
where

e t #- H (x 0 , p 0 ) is the solution at time t of ż(s) = #- H(z(s)), z(0) = (x 0 , p 0 ).
It is defined for small enough nonnegative time t and we can assume that p 0 belongs to S 1 . Definition 2.14. Let z(•) := (x(•), p(•)) be a reference extremal of order zero, defined on [0 , T ]. Let H be the Hamiltonian defined by [START_REF] Bonnard | Optimal Control with Applications in Space and Quantum Dynamics[END_REF]. The associated geodesic x(•) is called exceptional if H = 0, hyperbolic if H > 0 and elliptic if H < 0, along the reference extremal z(•).

Remark 2.15. The previous definition is related to the more classical definition 2.21. Even if the elliptic geodesics are not optimal according to the PMP, they still play a role in the analysis of the optimal synthesis, in particular in the computation of the cut locus when the current (or drift) is strong, see Section 3.

In cartesian coordinates, the Hamiltonian writes

H(x 1 , x 2 , p 1 , p 2 ) = µ x 2 1 + x 2 2 (-p 1 x 2 + p 2 x 1 ) + » p 2 1 + p 2 2 ,
and the extremals are solution of the following Hamiltonian system:

ẋ1 = -µ x 2 r 2 + p 1 p , ẋ2 = µ x 1 r 2 + p 2 p , ṗ1 = - µ r 4 2 x 1 x 2 p 1 -(x 2 1 -x 2 2 ) p 2 , ṗ2 = µ r 4 (x 2 1 -x 2 2 ) p 1 -2 x 1 x 2 p 2 .
Introducing the Mathieu transformation

Å p r p θ ã = Å cos θ sin θ -r sin θ r cos θ ã Å p 1 p 2 ã (14)
then, in polar coordinates, the Hamiltonian is given by (we still denote by p the covector in polar coordinates)

H(r, θ, p r , p θ ) = p θ µ r 2 + p r ,
where p r := » p 2 r + p 2 θ /r 2 . It is clear from H that θ is a cyclic variable and thus the problem has a symmetry of revolution and by Noether theorem, the adjoint variable p θ is a first integral. This relation p θ = constant corresponds to the Clairaut relation on surfaces of revolution. Hence, we can fix θ(0) = 0 and consider p θ has a parameter of the associated Hamiltonian system in polar coordinates:

ṙ = p r p r , θ = 1 r 2 Å µ + p θ p r ã , ṗr = p θ r 3 Å 2µ + p θ p r ã , ṗθ = 0. ( 15
)
2.3. C 1 -second order necessary conditions in the regular case

Since the extremals are of order zero, one can restrict u(t) to the 1-sphere S 1 . Writing u(t) =: (cos α(t), sin α(t)), we have with some abuse of notations

H = H 0 + u 1 H 1 + u 2 H 2 = H 0 + cos α H 1 + sin α H 2 ,
with α the new control. Differentiating twice with respect to α, we have

∂H ∂α = -sin α H 1 + cos α H 2 , ∂ 2 H ∂α 2 = -(cos α H 1 + sin α H 2 ),
and since u = (cos α, sin α) = Φ/ Φ and Φ = (H 1 , H 2 ) never vanishes along any extremal, we have

∂ 2 H ∂α 2 = - » H 2 1 + H 2 2 < 0
along any extremal. Hence, the strict Legendre-Clebsch condition is satisfied and we are in the regular case [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF], but with a free final time T .

Definition 2.16. Let z(•) be a reference extremal curve solution of ż = #-H(z) given by [START_REF] Bryson | Applied optimal control[END_REF]. The variational equation i

δz(t) = #- H (z(t)) • δz(t), (16
) is called a Jacobi equation. A Jacobi field is a non trivial solution J of [START_REF] Cesari | Optimization-theory and applications: problems with ordinary differential equations[END_REF]. It is said to be vertical at time t if δx(t) := Π (z(t)) • J(t) = 0, where Π : (x, p) → x is the standard projection.

Let z(•, x 0 , p 0 ) := (x(•, x 0 , p 0 ), p(•, x 0 , p 0 )) with p 0 ∈ S 1 , be a reference extremal curve solution of ż = #-H(z) with initial condition z(0, x 0 , p 0 ) = (x 0 , p 0 ) and defined over the time interval [0 , T ]. Following [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF], we make the following generic assumptions on the reference extremal in order to derive second order optimality conditions:

(A1) The trajectory x(•, x 0 , p 0 ) is a one-to-one immersion on [0 , T ].

(A2) The reference extremal is normal and strict.

Definition 2.17. Let z = (x, p) be the reference extremal defined hereinabove. Under assumptions (A1) and (A2), the time 0 < t c ≤ T is called conjugate if the exponential mapping is not an immersion at (t c , p 0 ). The associated point exp x0 (t c , p 0 ) = x(t c , x 0 , p 0 ) is said to be conjugate to x 0 . We denote by t 1c the first conjugate time.

The following result is fundamental, see [START_REF] Bonnard | Théorie des singularités de l'application entrée/sortie et optimalité des trajectoires singulières dans le problème du temps minimal[END_REF].

Theorem 2.18. Under assumptions (A1) and (A2), the extremities being fixed, the reference geodesic x(•) is locally time minimizing (resp. maximizing) for the L ∞ -topology on the set of controls up to the first conjugate time in the hyperbolic (resp. elliptic) case.

Algorithm to compute conjugate times. Writing the reference trajectory x(t) := x(t, x 0 , p 0 ) and considering a Jacobi field J(•) := (δx(•), δp(•)) along the reference extremal, which is vertical at the initial time (i.e. δx(0) = 0) and normalized by p 0 • δp(0) = 0 (since p 0 is restricted to S 1 ), then t c is a conjugate time if and only if t c is a solution of the following equation:

t → det(δx(t), ẋ(t)) = 0. ( 17
)
See [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF][START_REF] Cots | Contrôle optimal géométrique : méthodes homotopiques et applications[END_REF] for more details about algorithms to compute conjugate times in a more general setting and [START_REF] Caillau | Differential continuation for regular optimal control problems[END_REF] for details about the numerical implementation of these algorithms into the HamPath software.

The Zermelo-Carathéodory-Goh point of view

From the historical point of view in the Zermelo navigation problem, Zermelo and Carathéodory use for the parameterization of the geodesics the derivative of the heading angle α instead of the angle α itself, see [START_REF] Bryson | Applied optimal control[END_REF]. This corresponds precisely to the so-called Goh transformation for the analysis of singular trajectories in optimal control, see for instance the reference [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF]. This is presented next, in relation with the problem, to derive sufficient C 0 -optimality conditions under generic assumptions, see [START_REF] Bonnard | Théorie des singularités de l'application entrée/sortie et optimalité des trajectoires singulières dans le problème du temps minimal[END_REF].

Goh transformation

Retricting to extremals of order zero, the Goh transformation amounts to set (we use the same notation u for the new control but no confusion is possible): α = u, that is to take α as the control of the ship. Note that such a transformation transforms L ∞ -optimality conditions on the set of controls into C 0 -optimality conditions on the set of trajectories. For the geodesics computations, this amounts only to a reparameterization of extremal curves of order zero. Considering the vortex problem in cartesian coordinates (x 1 , x 2 ), we introduce x := (x 1 , x 2 , x 3 ), x 3 := α, and the control system becomes

ẋ = F (x) + u G(x), with F (x) := Å F 0 (x 1 , x 2 ) + cos x 3 F 1 (x 1 , x 2 ) + sin x 3 F 2 (x 1 , x 2 ) 0 ã and G = ∂ ∂x 3 .
The associated pseudo-Hamiltonian reads

H(x, p, u) := p • (F (x) + u G(x)), p := (p 1 , p 2 , p 3 ),
and relaxing the bound on the new control, the maximization condition implies p • G = 0 along any extremal. These extremals are called singular and the associated control is also called singular. Let us recall how we compute the singular extremals. First we need to introduce the concepts of Lie and Poisson brackets. The Lie bracket of two C ω vector fields X, Y on an open subset V ⊂ R n is computed with the convention:

[X, Y ](x) := ∂X ∂x (x) Y (x) - ∂Y ∂x (x) X(x),
and denoting H X , H Y the Hamiltonian lifts:

H X (z) := p • X(x), H Y (z) := p • Y (x), with z := (x, p) ∈ V × R n , the Poisson bracket reads: {H X , H Y } := H X • #- H Y = p • [X, Y ](x), where #- H X = ∂H ∂p ∂ ∂x - ∂H ∂x ∂ ∂p .
Differentiating twice p(•) • G(x(•)) with respect to the time t, one gets:

Lemma 2.19. Singular extremals (z(•), u(•)) are solution of the following equations:

H G (z(t)) = {H G , H F }(z(t)) = 0, {{H G , H F }, H F }(z(t)) + u(t) {{H G , H F }, H G }(z(t)) = 0.
If {{H G , H F }, H G } = 0 along the extremal, then the singular control is called of minimal order and it is given by the dynamic feedback:

u s (z(t)) := - {{H G , H F }, H F }(z(t)) {{H G , H F }, H G }(z(t)) .
Plugging the control u s in feedback form into the pseudo-Hamiltonian leads to define a true Hamiltonian denoted H s (z) := H(z, u s (z)), and one has: Lemma 2.20. Singular extremals of minimal order are the solutions of

ż(t) = #- H s (z(t)), with the constraints H G (z(t)) = {H G , H F }(z(t)) = 0.

The case of dimension 3 applied to the Zermelo problem

Consider the following affine control system: ẋ = F (x) + u G(x), with u ∈ R and x ∈ R 3 , where F , G, are C ω vector fields. Let z(•) := (x(•), p(•)) be a reference singular extremal curve on [0 , T ]. We assume the following:

(B1) The reference geodesic t → x(t) is a one-to-one immersion on [0 , T ].

(B2) F and G are linearly independent along x(•).

(B3) G, [G, F ], [[G, F ], G] are linearly independent along x(•).
From (B3), p is unique up to a factor and the geodesic is strict and moreover u s can be computed as a true feedback:

u s (x) = - D (x) D(x) ,
where we denote:

D := det(G, [G, F ], [[G, F ], G]), D := det(G, [G, F ], [[G, F ], F ]).
Moreover, let us introduce D := det(G, [G, F ], F ). In the vortex problem, one has:

α = - D (x) D(x) .
In our problem with the Goh extension, one orients p(•) using the convention of the maximum principle: p(t) • F (x(t)) ≥ 0 on [0 , T ] and we introduce the following definition consistent with definition 2.14:

Definition 2.21. Under assumptions (B1), (B2) and (B3), a geodesic is called:

• hyperbolic if DD > 0, • elliptic if DD < 0, • abnormal (or exceptional ) if D = 0.
Note that the condition DD ≥ 0 amounts to the generalized Legendre-Clebsch condition ∂ ∂u

d 2 dt 2 ∂H ∂u (z(t)) ≥ 0
and according to the higher-order maximum principle [START_REF] Krener | The high order maximal principle and its application to singular extremals[END_REF], this condition is a necessary (small) time minimization condition. See [START_REF] Bonnard | Théorie des singularités de l'application entrée/sortie et optimalité des trajectoires singulières dans le problème du temps minimal[END_REF] for the general frame relating the optimal control problems using the Goh transformation and applicable to our study and for the following result.

Theorem 2.22. Under assumptions (B1), (B2) and (B3), a reference hyperbolic (resp. elliptic) geodesic x(•) defined on [0 , T ] is time minimizing (resp. maximizing) on [0 , T ] with respect to all trajectories contained in a C 0 -neighborhood of x(•) if T < t 1c where t 1c is the first conjugate time along x(•) as defined by 2.17 for the projection of x(•) on the (x 1 , x 2 ) plane. In the exceptional case, the reference geodesic is C 0 -time minimizing and maximizing.

F(x, µ) := F 0 (x, µ) + 2 i=1 u i F i (x) u := (u 1 , u 2 ) ∈ U . (18)
Then, we have (noticing that if u ∈ U, then -u ∈ U):

0 ∈ F(x, µ) ⇔ ∃ u := (u 1 , u 2 ) ∈ U s.t. F 0 (x, µ) = 2 i=1 u i F i (x) ⇔ F 0 (x, µ) g = F 0 (x, µ) ≤ 1 ⇔ |µ| ≤ x = r.
This leads to introduce the following definition.

Definition 2.23. The drift F 0 (x, µ) is said to be weak at the point

x if F 0 (x, µ) < 1, strong at x if F 0 (x, µ) > 1 and moderate at x if F 0 (x, µ) = 1.
Remark 2.24. Note that if the drift could have been weak at any point x ∈ M , then we would have been in the Finslerian case [START_REF] Bao | An Introduction to Riemann-Finsler Geometry[END_REF] with a metric of Randers type. However, this is not possible for µ = 0 (the case µ = 0 is trivially Euclidean on R 2 ), since in this case, we have for any x ∈ M ∩ B(0, |µ|) = ∅ that the drift F 0 (x, µ) is not weak.

Remark 2.25. One can also notice that at the initial time, the strength of the drift depends on the ratio |µ|/r 0 . See Figure 2 for an illustration of the different possible strengths of the drift. The vortex is placed at the origin and marked by a black dot, as the initial point x 0 := (2, 0) in cartesian coordinates. The black circle corresponds to the set of initial directions x 0 + F(x 0 , µ) and the thick blue vector is the initial direction associated to the initial control u(0) := (cos α, sin α), with α = 7π/8. It is decomposed as the sum of the drift (oriented vertically) and the control field. On the left, we have a weak drift (µ = 0.5 r 0 ) at x 0 , in the middle we have a moderate drift (µ = r 0 ) and on the right a strong drift (µ = 2 r 0 ).

Influence of the circulation on the abnormal extremals

Proposition 2.26. Let (x(•), p(•), p 0 , u(•)) be an abnormal extremal, that is p 0 = 0. Then, the drift is strong or moderate all along the geodesic.

Proof. Since the abnormal extremal is of order 0, then all along the extremal we have

H(x(t), p(t)) = p(t) • F 0 (x(t), µ) + p(t) = 0. So, the Cauchy-Schwarz inequality gives p(t) = |p(t) • F 0 (x(t), µ)| ≤ p(t) F 0 (x, µ)
and since p(t) = 0, the result follows.

Remark 2.27. According to the previous proposition, the abnormal geodesics (that is the projection of the abnormal extremals on the state manifold) are contained in the ball B(0, |µ|).

According to the PMP, p 0 = 0 for the abnormal extremals while p 0 < 0 for the normal extremals. In the normal case, by homogeneity, one can fix p 0 = -1 and the initial adjoint vector p 0 := p(0) of normal extremals lives in the one dimensional space p ∈ R 2 H(x 0 , p) = 1 . This parameterization is very classical. Another possibility is to set (p 0 , p 0 ) = 1 since by the PMP, the pair (p(•), p 0 ) does not vanish. Finally, in our case, since all the extremals are of order zero, that is since p does not vanish, we can also fix p 0 = 1. We consider this third possibility but in polar coordinates, that is, denoting p 0 := (p r (0), p θ ) (recalling that p θ is constant) we parameterize the initial adjoint vector by:

p 0 ∈ p ∈ R 2 p r0 = 1 .
We thus introduce α ∈ [0 , 2π) such that p r (0) = cos α, p θ = r 0 sin α, which gives the initial control v(0) = (p r (0), p θ /r 0 ) = (cos α, sin α). According to the Mathieu transformation ( 14), one has in cartesian coordinates that p x (0) = cos θ 0 cos α -sin θ 0 sin α and p y (0) = sin θ 0 cos α + cos θ 0 sin α, so in the particular case θ 0 = 0, we have u(0) = (p x (0), p y (0)) = (cos α, sin α). This parameterization has the advantage to cover the normal and the abnormal extremals. According to the PMP, we have at the initial time:

H(q 0 , p 0 ) = p θ µ r 2 0 + p 0 r0 = p θ µ r 2 0 + 1 = -p 0 ≥ 0,
with q 0 = (r 0 , θ 0 ). Introducing (with a slight abuse of notation) H(α) := µ sin α/r 0 + 1, then, the abnormal extremals are characterized by H(α) = 0 ⇔ sin α = -r 0 /µ. We have three cases:

• If the drift is weak at the initial point, then this equation has no solution which explains why there is no abnormal extremals. In this case, H(α) > 0 for any α and thus there are only hyperbolic geodesics. • If the drift is moderate at the initial point, that is if |µ| = r 0 , then the single abnormal extremal is parameterized by α = π/2 if µ < 0 and by α = 3π/2 if µ > 0, • In the last case when the drift is strong, then for a given µ, the equation H(α) = 0 has two distinct solutions α a 1 < α a 2 in [0 , 2π). If µ < 0, then α a 1 and α a 2 are contained in (0 , π) while if µ > 0, then α a 1 and α a 2 are contained in (π , 2π). We have in addition the following symmetry:

α a 2 = π -α a 1 if µ < 0 and α a 2 = 3π -α a 1 if µ > 0.
The normal extremals solution of the PMP are parameterized by the set {α ∈

[0 , 2π) | H(α) > 0} = [0 , α a 1 ) ∪ (α a 2 , 2π) while for α ∈ (α a 1 , α a 2 )
we have H(α) < 0. One can see on Figure 3, the two abnormal directions with two hyperbolic and elliptic directions (that is resp. associated to hyperbolic and elliptic geodesics). The two abnormal directions define the boundary of the cone of admissible directions and reveal a lack of accessibility in the neighborhood of x 0 .

Numerical results

Resolution of the shooting equation

We introduce the shooting mapping

S(T, p 0 ) := exp x0 (T, p 0 ) -x f , (19)
where x f is the target and exp is the exponential mapping defined by [START_REF] Caillau | Differential continuation for regular optimal control problems[END_REF]. The shooting mapping is defined for (T, p 0 ) ∈ (T, p 0 ) ∈ R + × S 1 T < t p0 where t p0 ∈ R * + ∪ {+∞} is the maximal time such that exp x0 (•, p 0 ) is well defined over [0 , t p0 ). Let (T, p 0 ) be a solution of S = 0 such that the associated extremal is normal. The shooting mapping is differentiable at (T, p 0 ) and if T is not a conjugate time, then its Jacobian is of full rank at (T, p 0 ), which is a necessary condition to compute numerically the BC-extremals by means of Newton-like algorithms. We present in the following some examples of hyperbolic geodesics fixing the initial condition to x 0 := (2, 0) and solving the shooting equations S = 0 thanks to the HamPath code [START_REF] Caillau | Differential continuation for regular optimal control problems[END_REF], for different final conditions and for different strengths of the drift.

HamPath code. A Newton-like algorithm is used to solve the shooting equation S(T, p 0 ) = 0. Providing H and S to HamPath, the code generates automatically the Jacobian of the shooting function. To make the implementation of S easier, HamPath supplies the exponential mapping. Automatic Differentiation is used to produce #-H and is combined with Runge-Kutta integrators to assemble the exponential mapping and the variational equations ( 16) used to compute conjugate times. See [START_REF] Caillau | Differential continuation for regular optimal control problems[END_REF][START_REF] Cots | Contrôle optimal géométrique : méthodes homotopiques et applications[END_REF] for more details about the code.

Example 1. For this first example we want to steer the particle from x 0 to x f := (-2, 0) with µ := 2 x 0 (strong drift). In this case, we obtain a final time T ≈ 1.641 and the shooting equation S = 0, is solved with a very good accuracy of order 1e -12 (it is the same for the others examples but it won't be mentioned anymore). The associated hyperbolic geodesic is portrayed on Figure 4. The point vortex is represented by a black dot as the initial condition. The initial velocity ẋ(0) is given with the boundary (the black circle) of x 0 + F(x 0 , µ), cf. eq. ( 18). One can see that the drift is strong since x 0 ∈ x 0 + F(x 0 , µ). Example 2. To emphazise the influence of the final condition, let us take again µ := 2 x 0 and set x f := (2.5, 0). We can note from Figure 4 that the solution turns around the point vortex and profits from the circulation. In this case we obtain a final time T ≈ 2.821. Examples 3-4. Here we want to observe what happens for a weak drift. We set µ := 0.5 x 0 and present two cases with x f := (-2, 0) (cf. left subgraph of Figure 5) and x f := (2.5, 0) (cf. right subgraph of Figure 5). When x f = (-2, 0), the final condition is the same as in the example 1 but since the drift is weaker, the final time is longer. This is because the particle takes advantage of the vortex circulation. On the other hand, for x f = (2.5, 0) (same final condition as example 2) and considering a weak drift, then the particle does not turn around the vortex.

Numerical results on the absence of conjugate points

According to Section 3.2, there are two types of geodesics reaching the boundary of the domain: the bounded geodesics that reach the vortex and the unbounded geodesics. There exists also a unique geodesic separating these two cases, which is called the separatrix. Fixing θ 0 = 0 and using the parameterization u(0) = (cos α, sin α), α ∈ [0 , 2π), from Section 2.5.2, we define for a pair (r 0 , µ) ∈ R * + × R * , the set Λ(r 0 , µ) of parameters α such that the associated geodesics converge to the vortex and Θ(r 0 , µ) the set of parameters such that the associated geodesics go to infinity (in norm). The sets Λ(r 0 , µ) and Θ(r 0 , µ) depend on the strength of the drift and the sign of µ, and they are given in Section 3.2. One can find in Figure 6 the smallest singular value, denoted σ min (•), of det( ẋ(•), δx(•)) over the time, see eq. [START_REF] Chenciner | Simple Choreographic Motions of N Bodies: A Preliminary Study[END_REF]. For a fixed α ∈ [0 , 2π), we denote by t α ∈ R * + ∪ {+∞} the maximal time such that the associated geodesic is well defined over [0 , t α ). If there exists a time t c ∈ (0 , t α ) such that σ min (t c ) = 0, then the time t c is a conjugate time. If not, the geodesic has no conjugate time. One can see from the left subgraph of Figure 6, that for any weak drift and for any α ∈ Θ(r 0 , µ), that there are no conjugate times and σ min (t) → 1 when t → t α = +∞. 2 From the right subgraph of Figure 6, it is clear that for any weak drift and for any α ∈ Λ(r 0 , µ), that there are no conjugate times and that σ min (t) → 0 when t → t α . On the two figures, the red curve corresponds to the separatrix. It is clear that for the separatrix σ min (t) → 0.5 when t → t α = +∞. Since one has similar numerical results in the strong drift case, we make the following conjecture:

Conjecture 2.28. For any (x 0 , µ) ∈ M × R * , the conjugate locus is empty. i=1 {α i } on the right subgraph. In this case, we stop the numerical integration when r(t) ≤ 10 -3 which explains why the minimal value of σ min is between 10 -5 and 10 -6 . The red curve on each plot corresponds to the separatrix.

Micro-local analysis and properties of the value function

Poincaré compactification on S 3 of the extremal dynamics and integrability results

The objective of this section is to provide the geometric frame to analyze the extremals of order zero. Reparameterizing, the flow defines a polynomial vector field wich can be compactified using Poincaré method to analyze the behaviors of extremal curves converging either to the origin or to the infinity, see Section 2.2. Thanks to the rotational symmetry, this foliation can be integrated which is crucial to define in the next section the concept of Reeb circle.

Introducing c := -p 0 , then from the condition H = p θ µ/r 2 + p r = -p 0 , one gets r 2 p r = cr 2 -µp θ . Plugging this into [START_REF] Carathéodory | Calculus of Variations and Partial Differential Equations of the First Order, Part 1, Part 2[END_REF], one obtains:

r 3 (cr 2 -µ p θ ) ṙ = r 5 p r , r 3 (cr 2 -µ p θ ) θ = λ 3 r 3 -λ 4 r, r 3 (cr 2 -µ p θ ) ṗr = λ 1 r 2 -λ 2 ,
where λ 1 := (2µc + p θ )p θ , λ 2 := 2(µp θ ) 2 , λ 3 := µc + p θ and λ 4 := µ 2 p θ . The aim is to get a polynomial system to integrate, so we use the time reparameterization: dt r 3 (cr 2 -µp θ ) = ds.

Denoting by the derivative with respect to s, the system is written:

r = r 5 p r , θ = λ 3 r 3 -λ 4 r, p r = λ 1 r 2 -λ 2 .
We use Poincaré compactification where the sphere is identified to x = 1 and this leads to the following dynamics:

r = r 5 p r , θ = λ 3 r 3 x 3 -λ 4 rx 5 , p r = λ 1 r 2 x 4 -λ 2 x 6 , x = 0, (20)
which can be projected onto the 3-sphere r 2 + θ 2 + p 2 θ + x 2 = 1, up to a time reparameterization. Integrating the system [START_REF] Godbillon | Feuilletages, études géométriques[END_REF], one has x = c 1 and we get:

dr dp r = r 5 p r λ 1 r 2 x 4 -λ 2 x 6
which can integrated by separating the variables: dr r 5 λ 1 r 2 x 4 -λ 2 x 6 = p r dp r .

Hence, we obtain:

λ 2 c 6 1 4r 4 - λ 1 c 4 1 2r 2 = 1 2 p 2 r + K 1 .
Using quadratures from [START_REF] Godbillon | Feuilletages, études géométriques[END_REF], one can obtain s = f (r), θ = g(r). To illustrate the integrability properties, we only give the formula in the abnormal case for the r component:

Proposition 3.1. For α ∈ {α a
1 , α a 2 } (see Section 2.5.2 for the definition of α a 1 and α a 2 ), one has:

r(t) = Å s(t) 2 - 2 cos α r 0 s(t) + 1 r 2 0 ã -1 2
, where s(t) := sin α r 0 tan Å sin α r 0 t + arctan -cos α sin α ã + cos α r 0 .

Micro-local analysis of the extremal solutions

Since θ is a cyclic variable, and so p θ is a parameter, then we can focus the analysis on the subsystem

ṙ = p r p r , ṗr = p θ r 3 Å 2µ + p θ p r ã , (21)
to determine the different types of extremals. First of all, it is clear that the unique equilibrium point satisfying (p r , p θ ) = (0, 0) and r ≥ 0 is given by (r e , p r,e ) := (2|µ|, 0), and p θ must satisfy sign(p θ ) = -sign(µ). Now, introducing

a(p θ , µ) := µ 2 p 2 θ and b(p θ , r 0 , µ) := -p θ Å p θ + 2 µ + 2 p θ µ 2 r 2 0 ã , one can define ϕ(r, p θ , r 0 , µ) := a(p θ , µ) Å 1 r 4 - 1 r 4 0 ã + b(p θ , r 0 , µ) Å 1 r 2 - 1 r 2 0 ã + 1 - p 2 θ r 2 0
and one has the following result: Lemma 3.2. Along any extremal parameterized by (p r (0), p θ ) r(0) = 1 holds p r (t) 2 = ϕ(r(t), p θ , r(0), µ).

Proof. From ( 21), one has along the extremal (omitting the time variable):

p r ṗr = p θ r 3 (p θ + 2µ p r ) ṙ. ( 22
)
Using the parameterization p 0 r0 = 1 (with p 0 := (p r (0), p θ ) and r 0 := r(0)), one has

p(t) r(t) = - Å p 0 + p θ µ r(t) 2 ã = p θ µ Å 1 r 2 0 - 1 r(t) 2 ã + 1, (23)
since along the extremal the Hamiltonian is constant and equal to -p 0 . Putting ( 23) in ( 22), one has

p r ṗr = - Å 2a r 5 + b r 3 ã ṙ.
Integrating this equation, we have

1 2 p 2 r = 1 2 Å a r 4 + b r 2 + c ã , c := p r (0) 2 - a r 4 0 - b r 2 0 ,
and since p r (0) 2 = 1 -p 2 θ /r 2 0 , we get the conclusion. Let us fix p θ , r 0 and µ and introduce the polynomial function of degree two P (X) := aX 2 + bX + c, with

∆ := b 2 -4ac = p 3 θ Å p θ Å 1 + 4 µ 2 r 2 0 ã + 4µ
ã , its discrimant. Then, we have ϕ(r, p θ , r 0 , µ) = P (1/r 2 ) and the sign of the discrimant plays a crucial role in the analysis.

Remark 3.3. If p θ = 0 then P (X) = 1 and if µ = 0 then P (X) = -p 2 θ X + c. The polynomial P is of degree two if and only if p θ = 0 and µ = 0. The case µ = 0 is the simple Euclidean case that we do not develop. If p θ = 0 then ṗr = 0, ṙ = ±1 and this case is clear: these are the fastest geodesics converging either to the vortex ( ṙ = -1) or going to infinity ( ṙ = +1).

Let us introduce now

p * θ := -

4 µ 1 + 4 µ 2 r 2 0 ,
such that, assuming p θ = 0, µ > 0 and r 0 > 0, we have:

∆(p θ , r 0 , µ)      < 0 if p θ ∈ (p * θ , 0), = 0 if p θ = p * θ , > 0 if p θ ∈ R \ [p * θ , 0]. Remark 3.4. If µ < 0, then p * θ > 0 and we have ∆(p θ , r 0 , µ) < 0 if p θ ∈ (0 , p * θ ), ∆(p * θ , r 0 , µ) = 0 and ∆(p θ , r 0 , µ) > 0 if p θ ∈ R \ [0 , p * θ ]
. Remark 3.5. Note that with our parameterization, we are only interested in parameters p θ satisfying |p θ | ≤ r 0 . We can check that for any r 0 > 0 and µ = 0, we have |p * θ | ≤ r 0 . The particular case |p * θ | = r 0 is given by r 0 = 2|µ|. Thus, setting p θ = r 0 sin α, we can find two parameters 0

< α * 1 ≤ α * 2 < 2π such that p * θ =: r 0 sin α * 1 = r 0 sin α * 2 . This gives us sin α * 1 = sin α * 2 = - 4 λ 1 + 4λ 2 , λ := µ r 0 ,
and we can see again that the ratio µ/r 0 is of particular interest. Note that α * 1 = α * 2 if r 0 = 2|µ|. Let us analyze the three different cases depending on the sign of the discrimant, considering that p θ = 0 and fixing the parameters µ > 0 and r 0 > 0.

• Case ∆ < 0. In this case, the polynomial P has no roots and so p r (t) 2 = P (1/r(t) 2 ) never vanishes. See Figure 7 for an illustrative example. Thus, p r is of constant sign such as ṙ and so r is monotone, even strictly monotone since p r is nonzero. This case is given by p θ ∈ (p * θ , 0), thus, in this case |p θ | ≤ |p * θ | ≤ r 0 and setting (p r (0), p θ ) = (cos α, r 0 sin α), then r(t) → 0 (i.e. to the vortex) when t → t α if cos α < 0 and r(t) → +∞ when t → t α if cos α > 0, where t α ∈ R * + ∪ {+∞} is the maximal positive time such that the associated geodesic is defined over [0 , t α ). Denoting

α * := arcsin Å p * θ r 0 ã ∈ [- π 2 , 0)
and introducing (cf. remark 3.5)

α * 1 := π -α * and α * 2 = 2π + α * = 3π -α * 1 , then, for α ∈ (π , α * 1 )
, we have r(t) → 0, while for α ∈ (α * 2 , 2π), we have r(t) → +∞, when t → t α . • Case ∆ = 0. Since we assume p θ = 0 then ∆(p θ , r 0 , µ) = 0 if and only if p θ = p * θ . In this case P (1/r 2 ) = 0 if and only if r = r * := 2µ = r e . See Figure 8 for an illustrative example. We have four possibilities:

-If r 0 > r * and p r (0) > 0, that is p r (0) = cos α * 2 , then r(t) → +∞ when t → t α * 2 .

-If r 0 < r * and p r (0) < 0, that is p r (0) = cos α * 1 , then r(t) → 0 when t → t α * 1 .

-If r 0 > r * and p r (0) = cos α * 1 < 0, or if r 0 < r * and p r (0) = cos α * 2 > 0, then r(t) → r * . -If r 0 = r * and p r (0) = cos α * 1 = cos α * 2 = 0, then r(t) = r * for any time t and thus, in this case, the geodesic describes a circle. Remark 3.6. We can notice that (r e , p r,e ) = (2µ, 0) is an unstable hyperbolic fixed point and its associated linearized system has ±1/2µ as eigenvalues with associated eigenvectors (±µ(4µ 2 +r 2 0 )/r 2 0 , 1). Remark 3.7. Let us consider the θ dynamics:

θ = 1 r 2 Å µ + p θ p r ã .
Fixing p θ = -4µr 2 /(r 2 + 4µ 2 ) and using the parameterization p r = 1 (this amounts to consider r as the initial distance to the vortex), then we have θ = 0 ⇔ µ + p θ = 0 ⇔ r = 2|µ|/ √ 3. See also eq. ( 24).

• Case ∆ > 0. In this case, P has two roots X -< X + given by:

X ± := p 2 θ r 2 0 + 2 µ 2 p 2 θ + 2 µ p θ r 2 0 ± r 0 » p 3 θ (4 p θ µ 2 + 4 µ r 2 0 + p θ r 2 0 ) 2 µ 2 p 2 θ r 2 0 . Figure 8. Case ∆ = 0. Graph of r → P (1/r 2 ). Setting: µ = 2, r 0 = 3µ, p θ = p * θ .
By a tedious calculation, one can prove that X -≥ 0 and X -= 0 if and only if

p θ = -r 2 0 /µ < p * θ < 0. Let us introduce r 1 := 1 X + and r 2 := 1 X - ,
with the convention 1/0 = +∞. Then, for any r ∈ (r 1 , r 2 ), we have P (1/r 2 ) < 0 and thus it clear that r 0 ∈ (r 1 , r 2 ), since by definition P (1/r 2 0 ) = p r (0) 2 ≥ 0. (1) Subcase p θ > 0. Since p θ is positive we have 0 < X -< X + and thus 0 < r 1 < r 2 < +∞. We can easily prove that p θ > 0 ⇒ r 1 < r 0 since:

r 1 < r 0 ⇔ 1 < r 2 0 X + ⇔ 0 < p 2 θ r 2 0 + 2 µ p θ r 2 0 + r 0 » p 3 θ (4 p θ µ 2 + 4 µ r 2 0 + p θ r 2 0 ).
This can be seen on Figure 9. Since r 0 ∈ (r 1 , r 2 ), then p θ > 0 ⇒ r 2 ≤ r 0 and thus r(t) will go to +∞ in any case. More precisely, setting (p r (0), p θ ) = (cos α, r 0 sin α), then we have:

-If α ∈ (0 , π/2), then ṙ(0) > 0 and thus r 2 < r 0 . In this case we have ṙ(t) > 0 for any t ≥ 0 and thus r is strictly increasing and r(t) → +∞ when t → t α .

-If α = π/2, then ṙ(0) = 0, that is r 2 = r 0 , but r(t) is stricty increasing for t > 0 and we still have r(t) → +∞ when t → t α .

-If α ∈ (π/2 , π), then ṙ(0) < 0 and again r 2 < r 0 . In this case, r(t) is decreasing over [0 , tα ], where tα > 0 is the time such that P (1/r( tα ) 2 ) = 0, that is such that r( tα ) = r 2 . Since (r 2 , 0) is not an equilibrium point (a necessary condition is sign(p θ ) = -sign(µ)), and since along the geodesic holds p r (t) 2 = P (1/r(t) 2 ), then necessarily, r(t) is increasing over [ tα , t α ) (the sign of p r changes at t = tα and no more changing may occur after), even strictly increasing for t > tα and we still have r(t) → +∞ when t → t α .

(2) Subcase p θ < p * θ . Here we have -r 0 ≤ p θ < p * θ < 0. Setting (p r (0), p θ ) = (cos α, r 0 sin α), this case corresponds to α ∈ (α * 1 , α * 2 ). Let us recall that if r 0 = r * = 2µ, then p * θ = -r 0 (i.e. α * 1 = α * 2 = 3π/2), and so if r 0 = r * , then this last case is empty. We thus have either r 0 < r * or r 0 > r * . Let us prove now that r 1 < r * . To do this, we need the first following result:

p θ < p * θ ⇔ p θ r 2 0 + 4 µ 2 p θ + 4 µ r 2 0 < 0.
Besides,

r 1 < r * = 2 µ ⇔ 1 < 4 µ 2 X + ⇔ 0 < p θ p θ r 2 0 + 4 µ 2 p θ + 4 µ r 2 0 + 2 √ ∆, Figure 9. Case ∆ > 0 and p θ > 0. Graph of r → P (1/r 2 ). Setting: µ = 2, r 0 = 3µ, p θ = r 0 /2.
which is true since p θ < 0 and p θ < p * θ . We are now in position to conclude: -If r 0 > r * then necessarily r 0 ≥ r 2 (see Figure 10) since r 0 ∈ (r 1 , r 2 ) and thus r(t) → +∞ when t → t α for any α -If r 0 < r * , using the fact that |p θ | ≤ r 0 < r * = 2µ, then one can prove by a tedious calculation that r 0 ≤ r 1 (see Figure 11). Then, we have r(t) → 0 when t → t α for any α

∈ (α * 1 , α * 2 ). More precisely, r(t) is strictly increasing if α ∈ (3π/2 , α * 2 ), it is increasing if α = 3π/2 ( ṙ(0) = 0) and has one oscillation if α ∈ (α * 1 , 3π/2).
∈ (α * 1 , α * 2 ). More precisely, r(t) is strictly decreasing if α ∈ (α * 1 , 3π/2), it is decreasing if α = 3π/2 ( ṙ(0) = 0) and has one oscillation if α ∈ (3π/2 , α *
2 ). Remark 3.8. According to Section 2.5.2, the abnormal extremals are given by p a θ := -r 2 0 /µ, that is by α = α a 1 and α = α a 2 (defined in Section 2.5.2). Since (recalling that µ > 0)

p a θ -p * θ = - r 4 0 µ(r 2 0 + 4µ 2 ) < 0, it is clear that p a θ < p * θ , which gives α * 1 < α a 1 ≤ α a 2 < α * 2
, and thus the abnormal case is contained in this last case of ∆ > 0. Besides, the abnormal extremals exist only if the drift is strong or moderate, that is if r 0 ≤ µ < r * . To end this remark, let us mention that if the drift is strong, then for any α ∈ (α a 1 , α a 2 ) = ∅ the associated extremal is elliptic and from the previous analysis, one can say that all the elliptic geodesics converge to the vortex, that is r(t) → 0 when t → t α . We want now to group together the geodesics converging to the vortex and separate them from the ones going to infinity. We thus introduce for r 0 > 0 and µ = 0, the following sets contained in [0 , 2π):

Λ(r 0 , µ) :=          [π , α * 1 ) if r 0 ≥ 2µ and µ > 0 [π , α * 2 ) if r 0 < 2µ and µ > 0 (α * 2 , π] if r 0 ≥ 2|µ| and µ < 0 (α * 1 , π] if r 0 < 2|µ| and µ < 0 Θ(r 0 , µ) :=          [0 , π) ∪ (α * 1 , 2π) if r 0 > 2µ and µ > 0 [0 , π) ∪ (α * 2 , 2π) if r 0 ≤ 2µ and µ > 0 [0 , α * 2 ) ∪ (π , 2π) if r 0 > 2|µ| and µ < 0 [0 , α * 1 ) ∪ (π , 2π) if r 0 ≤ 2|µ| and µ < 0 Ψ(r 0 , µ) := ® {α * 1 } if (r 0 > 2µ and µ > 0) or if (r 0 ≤ 2|µ| and µ < 0) {α *
2 } if (r 0 ≤ 2µ and µ > 0) or if (r 0 > 2|µ| and µ < 0) Definition 3.9. The geodesics parameterized by α ∈ Ψ(r 0 , µ) are called separating geodesics.

We can now summarize the analysis:

Theorem 3.10. Let us fix x 0 ∈ M and µ = 0. Then, to any α ∈ [0 , 2π) is associated a unique zero order extremal z(•) := (x(•), p(•)) solution of the Hamiltonian system ż = #-H(z) given by eq. ( 12), satisfying x(0) = x 0 , and parameterized in polar coordinates by p(0) r0 = 1 (with r 0 := x 0 ). Besides, we have:

∀ α ∈ Λ(r 0 , µ) : lim t→tα r(t) = 0, ∀ α ∈ Θ(r 0 , µ) : lim t→tα r(t) = +∞ with t α = +∞, ∀ α ∈ Ψ(r 0 , µ) : lim t→tα r(t) = r * = 2|µ| with t α = +∞.

Reeb foliations

3-dimensional compact Reeb component

A concrete example of foliation of codimension 1 on the 3-sphere S 3 is given by the Reeb foliation [START_REF] Godbillon | Feuilletages, études géométriques[END_REF]. To construct such a foliation, the idea is to glue together two 3-dimensional foliations on two solid tori along their boundary. Let us recall only the construction of such a foliation on the solid torus D 2 × S 1 , where D 2 := (x, y) ∈ R 2 x 2 + y 2 ≤ 1 = B(0, 1), following the presentation of [START_REF] Hector | Les feuilletages de Reeb[END_REF]. Let f : D 2 → R be a smooth mapping such that the following conditions are satisfied:

(1) f is of the form: f (x, y) =: ϕ(x 2 + y 2 ), that is in polar coordinates, f (r, θ) = ϕ(r 2 ), (2) f (∂D 2 ) = {0} and f (x, y) > 0 for (x, y) ∈ ∂D 2 = S 2 , (3) f has no critical points on ∂D 2 .

Let consider the smooth mapping F : D 2 × R → R defined by F (x, y, t) := f (x, y) e t , where (x, y) ∈ D 2 and t ∈ R. F is a submersion from D 2 × R to R and defines a foliation of D 2 × R whose leaves are the level sets of F : F -1 (a), a ≥ 0. This foliation is called the 3-dimensional non compact Reeb component and it is denoted by R 3 in [START_REF] Hector | Les feuilletages de Reeb[END_REF]. Since F (x, y, t + 1) = e F (x, y, t), then, R 3 defines a foliation on the quotient manifold D 2 × S 1 = D 2 × R / (x, y, t) ∼ (x, y, t + 1). This foliation being called the 3-dimensional compact Reeb component and denoted R 3 c . This foliation admits a unique compact leaf: the boundary of D 2 × R diffeomorphic to the torus T 2 = S 1 × S 1 . All the others leaves are diffeomorphic to R 2 and are also given by the quotient of the graphs of functions of the form (in polar coordinates): (r, θ) ∈ B(0, 1) → -ln(ϕ(r 2 )) + b, b ∈ R. Remark 3.11. A simple example which is detailed in [START_REF] Hector | Les feuilletages de Reeb[END_REF] is given by the mapping f (r, θ) = 1 -r 2 . Note that it is possible to construct a Reeb component on the solid torus with a function F which is a submersion only on B(0, 1) × R. Such an example is given by the function f (r, θ) = exp(-exp(1/(1 -r 2 ))). This possibility will be useful in the following section to construct a 2-dimensional Reeb component from the Vortex application.

Let us exhibit now a function f in the Vortex application that can be used to construct a 3-dimensional compact Reeb component and thus a Reeb foliation on S 3 . To this end, we focus on the separating geodesics, that is the geodesics associated to parameters α ∈ Ψ(r 0 , µ). These geodesics are parameterized by p θ = p * θ = -4µr 2 0 /(r 2 0 + 4µ 2 ) and p(0) r0 = 1, where p(0) = (p r (0), p θ ) is the initial covector. Let us fix the circulation µ > 0 and the initial distance r 0 > 0. Let α ∈ Ψ(r 0 , µ) denotes the single element contained in this set. Let us denote by (r(•), θ(•)) the associated separating geodesic and by γ its orbit, that is γ := {(r(t), θ(t)) | t ∈ (t α , t α )}, where t α < 0 is the minimal negative time such that the geodesic is well defined by backward integration. The time t α > 0 being the maximal positive time such that the geodesic is well defined by forward integration. Let us consider now an arbitrary point (r 1 , θ 1 ) ∈ γ and define σ as the orbit associated to the unique separating geodesic when we consider (r 1 , θ 1 ) as the new initial state. Then, it is clear that γ = σ since again the set Ψ(r 1 , µ) has only one element. Thus, considering that at any time t, r(t) is the initial distance to the vortex, leads to reparameterize the separating geodesics by setting p θ = -4µr(t) 2 /(r(t) 2 + 4µ 2 ) and p(t) r(t) = 1, which gives

p r (t) 2 = 1 - p 2 θ r(t) 2 = (r(t) 2 -4µ 2 ) 2 (r(t) 2 + 4µ 2 ) 2 .
In the case 0 < r 0 < 2µ, the parameter α ∈ Ψ(r 0 , µ) is given by α * 2 and thus p r (t) > 0 for any time t (see Section 3.2). In this case, the dynamics of the separating geodesics reduces to the 2-dimensional differential equations:

ṙ = p r p r = p r = 4µ 2 -r 2 4µ 2 + r 2 , θ = 1 r 2 Å µ + p θ p r ã = 1 r 2 (µ + p θ ) = µ r 2 Å 4µ 2 -3r 2 4µ 2 + r 2 ã . (24)
Writing

dt = 4µ 2 + r 2 4µ 2 -r 2 dr,
then by integration we have

t(r) = 4µ atanh Å r 2µ ã -r + c = 2µ ln Å 2µ + r 2µ -r ã -r + c = ln Ç Å 2µ + r 2µ -r ã 2µ e -r K å ,
with a constant c := ln K. Introducing 2µρ 2 := r, ρ ∈ [0 , 1], then one can define the function

f (ρ, θ) := ϕ(ρ 2 ) := Å 1 -ρ 2 1 + ρ 2 ã 2µ e 2µρ 2
which satisfies the conditions (1)-( 2)-( 3) to construct a 3-dimensional compact Reeb component from F (ρ, θ, t) := f (ρ, θ) e t .

2-dimensional compact Reeb component

In the Vortex application, one can find a foliation in the 2-dimensional state space given by the separating geodesics. To present this foliation, we need to introduce a generalized 2-dimensional Reeb component following the presentation of [START_REF] Hector | Les feuilletages de Reeb[END_REF]. We thus first recall what is a 2-dimensional Reeb component and then we introduce the generalization. Let f : [-1 , 1] → R be a smooth mapping such that the following conditions are satisfied:

(4) f is of the form: f (x) =: ϕ(x 2 ), ( 5) f ({-1, 1}) = {0} and f ((-1 , 1)) > 0, (6) f has no critical points on {-1, 1}.

Let consider the smooth mapping

F : [-1 , 1] × R → R defined by F (x, t) := f (x) e t , where x ∈ [-1 , 1] and t ∈ R. F is a submersion from [-1 , 1] × R to R and defines a foliation of [-1 , 1] × R whose leaves are the level sets of F : F -1 (a), a ≥ 0. This foliation is called the 2-dimensional non compact Reeb component and it is denoted by R 2 in [23]. Since F (x, t + 1) = e F (x, t), then, R 2 defines a foliation on the annulus [-1 , 1] × S 1 = [-1 , 1] × R / (x, t) ∼ (x, t + 1)
. This foliation being called the 2-dimensional compact Reeb component and denoted R 2 c . This foliation admits two compact leaves: the two boundaries of the annulus, diffeomorphic to S 1 . All the others leaves are diffeomorphic to R and each of their two extremities wind up around one of the two compact leaves. These non compact leaves are also given by the quotient of the graphs of functions of the form: x ∈ (-1 , 1) → -ln(ϕ(x 2 )) + b, b ∈ R. Remark 3.12. A simple example which is detailed in [START_REF] Hector | Les feuilletages de Reeb[END_REF] is given by the function f (x) = 1-x 2 . The maximum of this function is given by solving f (x) = 2x = 0, that is for x = 0. Note that it is possible to construct a Reeb component on the annulus with a function F which is a submersion only on (-1 , 1) × R considering for instance f (x) = exp(-exp(1/(1 -x 2 ))). Remark 3.13. Note that we can construct a 3-dimensional Reeb component from the 2-dimensional Reeb component if f satisfies the condition (4)-( 5)- [START_REF] Bolsinov | Integrable Hamiltonian Systems, Geometry, Topology, Classification[END_REF]. See [START_REF] Hector | Les feuilletages de Reeb[END_REF] for details.

Let us slightly generalize the construction of a Reeb component on the annulus. To do so, we consider a function f : [a , b] → R smooth on (a , b), a < b in R, and such that: c . This foliation, like R 2 c , has two compact leaves diffeomorphic to S 1 and all the others leaves are diffeomorphic to R and given by the quotient of the graphs of functions of the form:

(7) f ({a, b}) = {0}, f ((a , b)) > 0. Then, F (x, t) := f (x) e t is a submersion from (a , b) × R to R
x ∈ (a , b) → -ln(f (x)) + c, c ∈ R.
Remark 3.14. The most important difference with the classical case is that we do not impose that f (x) is of the form ϕ(x 2 ). Even in the case a = -b, if f (x) is not of the form ϕ(x 2 ), then we do not have necessarily f (-x) = f (x) and if we do not have this symmetry, then it is not possible to construct a 3-dimensional Reeb component from a 2-dimensional one using the construction detailed in [START_REF] Hector | Les feuilletages de Reeb[END_REF].

Let us exhibit now a function f from the Vortex application that satisfies the condition [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF]. We fix µ > 0. Coming back to eq. ( 24), one can write

dθ = µ r 2 Å 4µ 2 -3r 2 4µ 2 -r 2 ã dr.
Integrating we have:

θ(r) = atanh Å r 2µ ã + µ r + c = 1 2 ln Å 2µ + r 2µ -r ã + µ r + c = ln Ç Å 2µ + r 2µ -r ã 1 2 e µ r K å ,
with a constant c := ln K. Introducing for r ∈ [0 , 2µ]:

f (r) := Å 2µ -r 2µ + r ã 1 2 e -µ r , (25)
with f (0) = 0, then the separating geodesics inside the punctured disk of radius 2µ are exactly given by the level sets of the function F (r, θ) := f (r) e θ . Besides, the function f satisfies the condition ( 7) and thus, it defines a generalized 2-dimensional compact Reeb component on the annulus [0 , 2µ] × S 1 , see Figure 12. The separating geodesics inside the punctured disk of radius 2µ correspond to the non compact leaves of the foliation while the unique circle geodesic, that is the separating geodesic with r 0 = 2µ, corresponds to one of the two compact leaves. The vortex corresponds to the second compact leaf. Figure 13 shows the foliation on the punctured disk of radius 2µ in the (x 1 , x 2 )-plane of the Vortex application. Definition 3.15. In the Vortex application, the circle of radius 2|µ| (assuming µ = 0) is called the Reeb circle.

Symmetries of the extremal curves

There exists three continuous symmetries and one discrete symmetry of main interest. The first continuous symmetry was considered at the beginning of the paper when we have fixed u max = 1. The second symmetry is the rotational symmetry that was mentioned at the end of Section 2.2. As a consequence of this symmetry of revolution, one can easily prove that:

∀∆θ ∈ R : V (r 0 , θ 0 , r f , θ f , µ) = V (r 0 , θ 0 + ∆θ, r f , θ f + ∆θ, µ),
where V is the value function [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF] expressed in polar coordinates. Hence, one can fix for instance θ 0 = 0 and reduce the set of variables of the value function. Note that the Reeb circle is the unique geodesic invariant by the rotational symmetry. The third continuous symmetry comes from the invariance of the set of geodesics by homothety combined with time/circulation dilatation. Let us detail this third symmetry. We start by introducing the change of variables (r, θ) := (λr, θ) for a given λ ∈ R * + . The associated Mathieu transformation leads to p := (p r , p θ ) = p r λ , p θ and the Hamiltonian system [START_REF] Carathéodory | Calculus of Variations and Partial Differential Equations of the First Order, Part 1, Part 2[END_REF] in polar coordinates becomes

ṙ = λ ṙ = λ p r p r = λ p r λ p r = λ p r p r , θ = θ = 1 r 2 Å µ + p θ p r ã = λ 2 r2 Å µ + p θ λ p r ã = λ r2 Å λµ + p θ p r ã
for the state and

ṗr = ṗr λ = p θ λr 3 Å 2µ + p θ p r ã = λ 2 p θ r3 Å 2µ + p θ λ p r ã = λp θ r3 Å 2λµ + p θ p r ã , ṗθ = 0,
for the covector, since

p 2 r = p r λ , p θ 2 r = p 2 r λ 2 + p 2 θ λ 2 r 2 = 1 λ 2 p 2 r .
Introducing now μ := λµ, the time reparameterization t = ψ(τ ) =: τ /λ, and the notation (r, θ, pr , pθ )(τ ) := (r, θ, p r , p θ )(ψ(τ )), then the system in the time τ writes

d dτ r(τ ) = d dτ r(ψ(τ )) = ṙ(ψ(τ ))ψ (τ ) = ṙ(t) λ = pr (τ ) p(τ ) r(τ )
,

d dτ θ = 1 r2 Å μ + pθ p r ã , d dτ pr = pθ r3 Å 2μ + pθ p r ã , d dτ pθ = 0,
which is equivalent to the original system [START_REF] Carathéodory | Calculus of Variations and Partial Differential Equations of the First Order, Part 1, Part 2[END_REF]. As a conclusion, the set of solutions of the problem (P ) is invariant by the dilatation (r(•), θ(•), t, µ) → σ λ (r(•), θ(•), t, µ) := (λr(•), θ(•), λt, λµ), for any λ ∈ R * + , and we have 1 λ V (λr 0 , θ 0 , λr f , θ f , λµ) = V (r 0 , θ 0 , r f , θ f , µ).

It is thus possible to fix for instance r 0 and again to reduce the set of variables of the value function.

Remark 3.16. Setting λ := 1/r 0 and considering also the invariance by rotation, we can write:

V (r 0 , θ 0 , r f , θ f , µ) = r 0 V (1, 0, ∆r, ∆θ, µ r 0 ),
with ∆r := r f /r 0 and ∆θ := θ f -θ 0 . This shows that the value function depends in fact only on the three variables ∆r, ∆θ and µ/r 0 .

We end this section presenting the discrete symmetry on the set of extremals. Let (r(•), θ(•), p r (•), p θ ) be a reference extremal on [0 , T ] and introduce θ 0 := θ(0). Let consider the following discrete symmetry: z(•) := (r(•), θ(•), p r (•), p θ ) := (r(•), 2θ 0 -θ(•), -p r (•), p θ ). Then, z(•) is solution of ż = -#-H(z), z(0) = (r(0), θ 0 , -p r (0), p θ ).

This means that z(•) is obtained by backward integration from the same initial condition in the state space than z(•) but with the initial covector (-p r (0), p θ ). Let us extend the reference extremal to obtain a maximal solution of ż = #-H(z), z(0) = (r(0), θ 0 , p r (0), p θ ). In this case, its orbit in the state space is given by γ := Im(r(•), θ(•)) and z is also a maximal solution of its associated differential equation, and we denote by γ := Im(r(•), θ(•)) its orbit in the state space. Then, we have:

Å 1 0 0 -1 ã γ + Å 0 2θ 0 ã = γ,
that is γ is obtained by an affine reflection of axis θ = θ 0 . In cartesian coordinates, one has a linear reflection of axis Rx 0 , see Figure 14.

3.5. Properties of the value function and its level sets in the weak case Let us introduce the mapping V µ (x 0 , x f ) := V (x 0 , x f , µ), where V is the value function defined by [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF]. We are first interested in the continuity of x f → V µ (x 0 , x f ) at points x f where the drift is weak, that is for target x f such that x f > |µ|. This situation is more general than the Finslerian case with Randers metric since the drift is not necessarily weak all along the geodesics. We have: Figure 14. Geodesics in the state space in the (x 1 , x 2 )-plane. The vortex is placed at the origin and represented by a black dot. The initial point, denoted x 0 , is the black dot at an intersection of both geodesics. The geodesic in green is obtained by a reflection of axis Rx 0 from the blue geodesic and vice-versa. Assuming that the initial point is given at time t = 0, then the plain lines are obtained for times t > 0 while the dashed parts are for negative times. On the left subgraph is represented a normal geodesic with both extremities going to infinity while on the right subgraph is represented the two abnormals (in the case of a strong drift) with a cusp at r = |µ|. Both extremities of the abnormal geodesics go the vortex. Proposition 3.17. Let µ = 0 be the vortex circulation and x 0 ∈ M the initial condition. Then, V µ (x 0 , •) is continuous at any point x f such that x f > |µ|.

Proof. Since we consider that the drift is weak at the target x f then local controllability around x f holds, that is we have the following:

∀ ε > 0, ∃ η ε > 0 s.t. ∀ (x 1 , x 2 ) ∈ B(x f , η ε ) × B(x f , η ε ) : V µ (x 1 , x 2 ) ≤ ε. (26)
Let us fix ε > 0 and consider x ∈ B(x f , η ε ). By definition of the value function and by [START_REF] Jurdjevic | Geometric Control Theory[END_REF], we have

V µ (x 0 , x) -V µ (x 0 , x f ) ≤ V µ (x f , x) ≤ ε and V µ (x 0 , x f ) -V µ (x 0 , x) ≤ V µ (x, x f ) ≤ ε,
which prove the result.

We are now interested in the construction of the optimal synthesis associated to problem (P ) for given values of µ and x 0 , that is we want to get for any x f ∈ M the value of V µ (x 0 , x f ) together with the optimal geodesic joining the points x 0 and x f . To construct such an optimal synthesis, a classical approach is to compute the level sets of V µ (x 0 , •), which corresponds in Riemmanian or Finslerian geometry to compute the spheres centered in x 0 . To compute the spheres, one preliminary work is to compute the cut locus, that is the set of cut points where the geodesics cease to be optimal. Let us introduce some notations before characterizing the cut points. Following Section 2.5.2 and fixing x 0 ∈ M , we introduce the notation in polar coordinates p 0 (α) := (cos α, r 0 sin α), with α ∈ [0 , 2π) and where r 0 := x 0 . Let t ∈ R + , we introduce the wavefront from x 0 at time t by

W(x 0 , t) := exp x0 (t, p 0 (α)) α ∈ [0 , 2π) s.t. t < t α ,
where the exponential mapping is given by ( 13) and depends on the circulation µ supposed to be given, and where we recall that t α ∈ R * + ∪ {+∞} is the maximal positive time such that the associated geodesic is defined over [0 , t α ). Then, we define what corresponds to the balls and spheres in Riemannian or Finslerian geometry:

B(x 0 , t) := {x ∈ M | V µ (x 0 , x) < t} , B(x 0 , t) := {x ∈ M | V µ (x 0 , x) ≤ t} , S(x 0 , t) := {x ∈ M | V µ (x 0 , x) = t} .
We use the terminology ball and sphere by analogy with Geometry. Note that the spheres S(x 0 , •) are the level sets of V µ (x 0 , •). To construct the optimal synthesis, one of the main important loci to compute is the cut locus.

We define first the cut times as:

t cut (x 0 , α) := sup t ∈ R * + exp x0 (•, p 0 (α)) is optimal over [0 , t] .
Then, the cut locus from x 0 is simply defined by:

Cut(x 0 ) := x ∈ M ∃ (t, α) ∈ R * + × [0 , 2π) s.t. x = exp x0 (t, p 0 (α)) and t = t cut (x 0 , α) .
Any point of the cut locus is called a cut point. From the cut times, we define also the injectivity radius from x 0 as:

t inj (x 0 ) := inf α∈[0,2π) t cut (x 0 , α).
Remark 3.18. This time is of particular interest since if the drift is weak at the initial point x 0 , then exp x0 is a diffeomorphism from (0 , t inj (x 0 )) × S 1 to B(x 0 , t inj (x 0 )) \ {x 0 }. This is similar to the Riemannian and Finslerian cases and for times 0 ≤ t ≤ t inj (x 0 ), we have S(x 0 , t) = W(x 0 , t).

In Finslerian geometry [START_REF] Bao | An Introduction to Riemann-Finsler Geometry[END_REF] (as in Riemannian geometry [19, Proposition 2.2, Chapter 13]), the cut locus is a part of the union of the conjugate locus and the splitting locus, which is defined by:

Split(x 0 ) := ¶ x ∈ M ∃(t, α 1 , α 2 ) ∈ R * + × [0 , 2π) 2 s.t. α 1 = α 2 and exp x0 (t, p 0 (α 1 )) = exp x0 (t, p 0 (α 2 )) = x © .
We have the following similar result:

Proposition 3.19. Let x(•) be a geodesic. If x c := x(t c ) is a cut point along x(•) such that V µ (x 0 , •) is continuous at x c , then: (a) either x c is a conjugate point along x(•), (b) or x c is a splitting point along x(•).
Conversely if a) or b) holds, then there exists a time 0 ≤ t 1 ≤ t c such that x(t 1 ) is a cut point along x(•).

Remark 3.20. The proposition above gives a characterization of the cut points where the mapping V µ (x 0 , •) is continuous. In Riemannian and Finslerian geometry, the mapping V µ (x 0 , •) is replaced by the distance function which is continuous with respect to its second argument and so all the cut points are either conjugate or splitting points. By Proposition 3.17, if the drift is weak at the cut points, then V µ (x 0 , •) is continuous at this point and so the characterization holds. In the strong case, we can have abnormal minimizers and so we can have cut points which are neither conjugate nor splitting points. However, if the drift is weak at the initial point x 0 , then there are no abnormals and we can expect that the mapping V µ (x 0 , •) is continuous at any point, even where the drift is strong. This result would be useful to conclude in the following part since we are interested in the construction of the optimal synthesis for a given initial condition x 0 where the drift is weak. Note that the numerical experiments of the following section suggest that V µ (x 0 , •) is continuous at any point of M .

Level sets of the value function and optimal synthesis

In this section, we fix for the numerical experiments x 0 := (3, 0) and µ := 0.6 x 0 . The drift is thus weak at the initial position x 0 . We decompose the construction of the optimal synthesis in three steps. In the first step, we compute the splitting locus. In a second time we give the cut locus and we finish by the construction of the spheres and balls.

Step 1: Computation of the splitting locus. Let us introduce the mapping

F split (t, α 1 , x, α 2 ) := (x -exp x0 (t, p 0 (α 1 )), x -exp x0 (t, p 0 (α 2 ))) ∈ R 4 .
The splitting line is then given by solving F split = 0 since we have

Split(x 0 ) = {x ∈ M | ∃ (t, α 1 , α 2 ) s.t. α 1 = α 2 and F split (t, α 1 , x, α 2 ) = 0} .
Introducing y := (t, α 1 , x) ∈ R 4 and λ := α 2 ∈ R, we have to solve F split (y, λ) = 0 ∈ R 4 . Numerically, we compute the splitting line with the differential path following method (or homotopy method) of the HamPath software. Under some regularity assumptions, the set F -1 split (0) is a disjoint union of differential curves, each curve being called a path of zeros. To compute a path of zeros, we look for a first point on the curve by fixing the homotopic parameter λ to a certain value λ and calling a Newton method to solve F split (•, λ) = 0. Then, we use a Prediction-Correction (PC) method to compute the differential curve. The HamPath code implements a PC method with a high-order Runge-Kutta scheme with variable step-size for the prediction, hence the discretization grid of the homotopic parameter is computed by the numerical integrator. Besides, the Jacobian of F split is computed by automatic differentiation combined with variational equations of the exponential mapping. Remark 3.21. If there exists several paths of zeros, then each path has to be found manually and compared.

To compute a splitting curve (that is a path of zeros), we need to get a first point on the curve. This step is easy since the splitting points are self-intersections of the wavefronts. On the right subgraph of Figure 15, we can see that the wavefront3 W(x 0 , 3.5) has at least three self-intersections labelled 1, 2 and 3. Figure 15. Four wavefronts with different scaling at times t ∈ {0.5, 2, 2.8, 3.5}. The black dot is x 0 . For t = 3.5, the wavefront has at least three self-intersections labelled 1, 2 and 3.

Each self-intersection of the wavefront W(x 0 , 3.5) is a point of a splitting curve. We denote by Σ 1 , Σ 2 and Σ 3 the three splitting curves associated respectively to the self-intersections 1, 2 and 3, and we have Σ 1 ∪ Σ 2 ∪ Σ 3 ⊂ Split(x 0 ). On Figure 16 is represented the graph of the function α 2 → t(α 2 ) along the three splitting curves for times t ≤ 3.2, the homotopic parameter α 2 being strictly monotone along the splitting curves. Since Σ 1 is below Σ 2 and Σ 3 on the figure, 4 it is clear that only Σ 1 has a chance to be part of the cut locus, compared to Σ 2 and Σ 3 . Note that there may exist others splitting curves but the numerical experiments suggest that Σ 1 is below any of them. Remark 3.22. Note that the extremities of a splitting curve are the vortex and infinity according to the numerical experiments. Hence, we cannot compute all the curve. Numerically, we use an option of the HamPath code to stop the homotopy if the curve goes out the annulus of smaller circle of radius 0.005 and of larger one of radius 100. Let us describe the evolution of the wavefronts according to Figure 16. For t = 0, the wavefront is reduced to the singleton {x 0 }. For times t ∈ (0 , t inj (x 0 )) with t inj (x 0 ) = inf t(•) along Σ 1 , the wavefronts have no selfintersections. For t = t inj (x 0 ), the wavefront has one self-intersection, while for t ∈ (t inj (x 0 ) , t vor (x 0 )), where t vor (x 0 ) := x 0 is the minimal time to reach the vortex, the wavefronts have two self-intersections contained in Σ 1 , see Figure 17, and up to six self-intersections in Σ 1 ∪ Σ 2 ∪ Σ 3 . Finally, for times t ≥ t vor (x 0 ), the wavefronts have one self-intersection in Σ 1 and three in Σ 1 ∪ Σ 2 ∪ Σ 3 , see the right subgraph of Figure 15.

Figure 17. Three wavefronts with the same scaling and a zoom around the vortex (red dot) for t ∈ {2.86, 2.88, 2.9}. For t = 2.9 > t inj (x 0 ) ≈ 2.889, the wavefront has two self-intersections.

Step 2: Computation of the cut locus. According to Proposition 3.17, at a final point x f where the drift is weak the mapping V µ (x 0 , •) is continuous. In this case, from Proposition 3.19, if x f is a cut point, then either it is a conjugate point or a splitting point. Since, there are no conjugate points from Conjecture 2.28, x f is a splitting point. Finally, from step 1, we can conclude that if x f is a cut point then it belongs to Σ 1 . In this section, we consider moreover that the drift is weak at the initial condition x 0 . In this case, we conjecture the following: Conjecture 3.23. If the drift is weak at x 0 , then V µ (x 0 , •) is continuous at any point of M .

This conjecture means that the characterization of cut points from Proposition 3.19 holds in all the state space M . Hence, we have Cut(x 0 ) ⊂ Σ 1 . Now, by the converse part of Proposition 3.19 and according to the previous conjecture, we can conclude that:

Cut(x 0 ) = Σ 1 .

Step 3. Computations of the spheres and balls. Figure 18 represents four balls that we can group into three categories (see also Figure 16):

• Type A. The ball is simply connected in R 2 with a smooth boundary;

• Type B. The ball is not simply connected in R 2 ;

• Type C. The ball is simply connected in R 2 with a singularity on its boundary. For times t ∈ (0 , t inj (x 0 )], we have W(x 0 , t) = S(x 0 , t) and the ball B(x 0 , t) is of type A. The two balls on the top of Figure 18 are of type A. For times t ∈ (t inj (x 0 ) , t vor (x 0 )), the balls have a hole around the vortex and so they are of type B. The ball at the bottom-left subgraph of Figure 18 is of type B. One can see on Figure 19, the evolution of the balls around the vortex for times close to t inj (x 0 ). The hole appears when the ball self-intersects, like a snake biting its tail. For times t ≥ t vor (x 0 ), the balls are simply connected in R 2 since the hole has vanished (the vortex is reached at t = t vor (x 0 )) and the spheres have a singularity at the cut point. In this case the balls have a shape of an apple and are of type C. The bottom-right subgraph of Figure 18 represents a ball of type C.

One can see on Figure 20 the evolution of the spheres S(x 0 , t) at times t ∈ {0.5, 1.5, 2.5, 3.5} together with the cut locus Σ 1 . One can notice that the singularity on S(x 0 , 3.5) in contained in the cut locus. Putting all together, we have constructed the optimal synthesis which is given on Figure 21.

Conclusion

In this article, we have analyzed the Zermelo navigation problem with a vortex singularity. We have shown the existence of an optimal solution. Thanks to the integrability properties we made a micro-local classification of the extremal solutions, showing remarkably the existence of a Reeb foliation. The spheres and balls are described in the case where at the initial condition the current is weak. This gives us the optimal synthesis for any initial condition such that x 0 > |µ|.

Note that the current is not necessarily weak all along the geodesics and so the situation we have analyzed is more general than the Randers case from Finsler geometry. Still, this analysis has to be completed to the case of a strong current at the initial condition, in particular in relation with the abnormal extremals. Note that the historical work [START_REF] Carathéodory | Calculus of Variations and Partial Differential Equations of the First Order, Part 1, Part 2[END_REF] is a case study of the strong current situation which can be applied to our study. A possible further extension concerns the case of several vortices by analogy with the N-body problem.

Finally in a forthcoming article, the 2D-Zermelo navigation problem will be generalized to the case of arbitrary current and metric, with the symmetry of revolution. It will cover the historical problem and the vortex case. The analysis will be based on tools from Hamiltonian dynamics from [START_REF] Bolsinov | Integrable Hamiltonian Systems, Geometry, Topology, Classification[END_REF] and in particular the concept of Reeb graph.

We are grateful to Carlos Balsa for helpful discussions about vortex theory and for drawing authors' attention to this topic, and to Ulysse Serres for helpful discussions about Zermelo-like navigation problems.

2. 5 .

 5 Influence of the circulation 2.5.1. Influence of the circulation on the drift Denoting the drift (5) F 0 (x, µ) to emphasize the role of µ, one introduces for (x, µ) ∈ M × R * the set

Figure 2 .

 2 Figure2. The vortex is placed at the origin and marked by a black dot, as the initial point x 0 := (2, 0) in cartesian coordinates. The black circle corresponds to the set of initial directions x 0 + F(x 0 , µ) and the thick blue vector is the initial direction associated to the initial control u(0) := (cos α, sin α), with α = 7π/8. It is decomposed as the sum of the drift (oriented vertically) and the control field. On the left, we have a weak drift (µ = 0.5 r 0 ) at x 0 , in the middle we have a moderate drift (µ = r 0 ) and on the right a strong drift (µ = 2 r 0 ).

Figure 3 .

 3 Figure3. The vortex is placed at the origin and marked by a black dot, as the initial point x 0 := (2, 0) in cartesian coordinates. The drift is strong at the initial point since µ = 2 r 0 . The black circle represents the set of initial directions x 0 + F(x 0 , µ) and the blue vector is the initial hyperbolic direction associated to the initial control u(0) := (cos α, sin α), with α = 7π/8. The direction in green is elliptic while the two red directions are the abnormal directions located at the boundary of the cone of admissible directions.

2 Figure 4 . 2 Figure 5 .

 2425 Figure 4. Example 1 and 2. Geodesic with a strong drift at the initial point: µ = 2 x 0 . (Left: example 1) x 0 = (2, 0), x f = (-2, 0) and the final time is T ≈ 1.641. (Right: example 2) x 0 = (2, 0), x f = (2.5, 0) and the final time is T ≈ 2.821.

Figure 6 .

 6 Figure 6. Smallest singular value with respect to time. Setting: θ 0 = 0, µ = 2.0, r 0 = 4µ/3 (weak drift at the initial point). The set of parameters [0 , 2π) is uniformly discretized in N := 1000 sub-intervals and we define 0 =: α 1 < • • • < α N < 2π the associated parameters. Each blue curve is the graph of the smallest singular value σ min (t) for one α ∈ Θ(r 0 , µ) ∩ ∪ N i=1 {α i } and for t ∈ [0 , 50] on the left subgraph, and for one α ∈ Λ(r 0 , µ) ∩ ∪ Ni=1 {α i } on the right subgraph. In this case, we stop the numerical integration when r(t) ≤ 10 -3 which explains why the minimal value of σ min is between 10 -5 and 10 -6 . The red curve on each plot corresponds to the separatrix.

Figure 7 .

 7 Figure 7. Case ∆ < 0. Graph of r → P (1/r 2 ). Setting: µ = 2, r 0 = 3µ, p θ = p * θ /2 ∈ (p * θ , 0).

Figure 10 .

 10 Figure 10. Case ∆ > 0 and p θ < p * θ : r → P (1/r 2 ). Setting: µ = 2, r 0 = 3µ, p θ = (p * θ -r 0 )/2.

Figure 11 .

 11 Figure 11. Case ∆ > 0 and p θ < p * θ : r → P (1/r 2 ). Setting: µ = 2, r 0 = 4µ/3, p θ = (p * θ -r 0 )/2.

  and defines a foliation on [a , b] × R which is called the generalized 2-dimensional non compact Reeb component and denoted R 2,G . As in the two previous cases, R 2,G defines a foliation on the annulus [a , b] × S 1 which is called the generalized 2-dimensional compact Reeb component and denoted R 2,G

Figure 12 .

 12 Figure 12. (Left) Graphs of r ∈ (0 , 2µ) → -ln(f (r)) + c, c ∈ R, with f defined by eq. (25). The minimum is attained at r = 2µ/ √ 3 where f (r) = 0. (Right) Generalized 2-dimensional compact Reeb component of the annulus given by the function F (r, θ) = f (r) e θ . The intermediate dashed circle of radius 2µ/ √ 3 is orthogonal to all the non compact leaves in blue. The two compact leaves are the two red circles.

Figure 13 .

 13 Figure 13. Separating geodesics in the punctured disk of radius 2µ in the (x 1 , x 2 )-plane. The vortex is placed at the origin and represented by a red dot. The red circle is the Reeb circle of radius 2µ while the black and dashed circle is the circle of radius 2µ/ √ 3 where θ = 0 along the separating geodesics.

3. 5 . 1 .

 51 Continuity of the value function and characterization of the cut points

Figure 16 .

 16 Figure 16. Graph of the function t(α 2 ) along Σ 1 (blue), Σ 2 (red) and Σ 3 (green) for times t ≤ 3.2. The balls are of three types: A, B and C. See the description detailed in the step 3.

2 CFigure 18 .

 218 Figure 18. Four balls in the (x 1 , x 2 )-plane at times t ∈ {1.5, 2.8, 2.9, 3.5} with the initial condition x 0 represented by a black dot and the vortex by a red dot. The balls B(x 0 , 1.5) and B(x 0 , 2.8) on top are of type A. The ball B(x 0 , 2.9) (Bottom-Left) is of type B while B(x 0 , 3.5) (Bottom-Right) is of type C.

Figure 19 .

 19 Figure[START_REF] Do | Birkhäuser, Mathematics: Theory & applications[END_REF]. Three balls with a zoom around the vortex (red dot) at times t ∈ {2.86, 2.88, 2.9}. Compare to Figure17. For t = 2.9 the ball B(x 0 , t) is of type C and has a small hole with a fetus shape.

2 Figure 20 .

 220 Figure 20. Spheres S(x , t) at times t ∈ {0.5, 1.5, 2.5, 3.5} in black together with the cut locus Σ 1 in blue. The cut locus should go to the vortex but to gain in clarity, it is represented up to a distance of 0.12 to the vortex. The vortex being represented by a red dot while the initial condition is represented by a black dot.

Figure 21 .

 21 Figure 21. Optimal synthesis. The black curves represent the spheres of types A and C. The two blue curves is a single sphere of type B. Since the balls of type B are not simply connected and have one hole, the corresponding spheres have two connected components.

The standard inner product is written a • b or a, b .

It is clear that tα = +∞ for any α ∈ Θ(r 0 , µ), since | ṙ| ≤ 1.

The wavefronts are computed by homotopy. We compute only rough approximations of the extremities of the wavefronts which are not closed curves, like W(x 0 , 3.5). Hence, there may exists others self-intersections but which are not relevant for rest of the analysis.

We mean that the graph of α 2 → t(α 2 ) associated to Σ 1 is below the ones associated to Σ 2 and Σ 3 .

This research is supported by the French Ministry for Education, Higher Education and Research.