A Zermelo navigation problem with a vortex singularity

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.


Introduction
Helhmoltz and Kirchhoff originated the model of the displacement of particles in a two dimensional fluid, see [23,26] for the original articles and [2,29,35] for a modern presentation of Hamiltonian dynamics. The vorticity 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 , and the dynamics is given by the Hamiltonian canonical equations

where the Hamiltonian function H is
where z i − z j is the Euclidian distance. In this article, we consider a single vortex, motionless by eq. (1), that we fix at the origin of the reference frame. This corresponds for instance to fix z 1 = (x 1 , y 1 ) to (0, 0). The idea is therefore to consider a particle as a (point) vortex with zero circulation, considering for instance 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, defined by (1) 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.
To define a Zermelo navigation problem, following [14,38] and see [11] 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 and thus 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: 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 flat Riemannian 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 [4] for this geometric frame. Note that when the ship visits 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 [14,38] and their study is an important inspiration for our work. Extension of classical calculus of variation is optimal control, with the Hamiltonian formulation coming from the Pontryagin Maximum Principle [31], and 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 [28]. From the dynamics of such systems there is an intense activity research initiated by Poincaré [30] to compute periodic trajectories avoiding collisions and such techniques lead to the concept of choregraphy developed by [16] for the N -body problem and [13] for the N -vortex system, showing the relations between both dynamics in the Hamiltonian frame [28]. From the control point of view, there is a lot of development related to space navigation for the N -body problem, see [8], valuable in our study for ship navigation in the N -vortex problem. Optimal control problem in this aera was developed for this purpose in relation with Finsler geometry, see [10] or [36] for a general setting in the planar case. More general results about vortex control may be found in [32,37] for instance.

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 bẏ x 2 (t) r(t) 2 + u 1 (t),ẋ 2 (t) = k 2π with r(t) 2 := x 1 (t) 2 + x 2 (t) 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: with F 0 , F 1 and F 2 three real analytic (i.e. C ω ) vector fields, where the current (or drift) is given by 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) .
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 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 driving (4) in minimum time from the initial condition x 0 to the target x f : V (x 0 , x f , µ, u max ) := inf {T | (T, u) ∈ D x0 and x u (T, 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) We first prove that there exists an admissible trajectory connecting any pair of points in M .
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 [15,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. 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: where s := sign(µ), x 0 := (R cos θ 0 , R sin θ 0 ) and where the control u is associated to v by the relation 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 : 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.
For the second lemma, we introduce the following definition.
Definition 2.6. An admissible trajectory is a trajectory x associated to a pair (T, u) ∈ D x0 , that is such that x = x u (·, x 0 ), satisfying in addition the constraints x(T ) = x f . Let x 1 and x 2 denote two admissible trajectories and let T 1 and T 2 denote, respectively, the first time such that x 1 (T 1 ) = x f = x 2 (T 2 ). Then, we say that x 1 is better than x 2 if T 1 ≤ T 2 , and strictly better if T 1 < T 2 .
Let us fix (x 0 , x f , µ) ∈ M 2 × R * and introduce r 0 := x 0 and r f := x f . Then: There exists ε > 0 such that for any admissible trajectory having a contact with B(0, ε), one can construct a strictly better admissible trajectory 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 having a contact with B(0, ε) and associated to a pair denoted (T, u). Then, there exists two times 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 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. 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.7 and R 2 := r 0 + T * , with r 0 := x 0 . By lemma 2.7, 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.
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.10. All the extremals are of order zero, that is there are no switching times.
We have the standard following result: Proposition 2.11. 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 .
#-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.12. Let (x(·), p(·), p 0 , u(·)) be an extremal. Denoting z := (x, p), then, we have over [0 , T ]: 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: 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 ].
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.14. We define the exponential mapping by 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 . Remark 2.16. The previous definition is related to the more classical Definition 2.22. 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 and the extremals are solution of the following Hamiltonian system: Introducing the Mathieu transformation then, in polar coordinates, the Hamiltonian is given by (we still denote by p the covector in polar coordinates) 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: 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 and since u = (cos α, sin α) = Φ/ Φ and Φ = (H 1 , H 2 ) never vanishes along any extremal, we have along any extremal. Hence, the strict Legendre-Clebsch condition is satisfied and we are in the regular case [6], but with a free final time T .
Definition 2.17. Let z(·) be a reference extremal curve solution ofż = #-H(z) given by (12). The variational equationi is called a Jacobi equation. A Jacobi field is a non trivial solution J of (16). It is said to be vertical

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 [6], we make the following generic assumptions on the reference extremal in order to derive second order optimality conditions: The reference extremal is normal and strict.
Definition 2.18. 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 [9].
Theorem 2. 19. 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: See [6,17] for more details about algorithms to compute conjugate times in a more general setting and [12] 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 [11]. This corresponds precisely to the so-called Goh transformation for the analysis of singular trajectories in optimal control, see for instance the reference [7]. This is presented next, in relation with the problem, to derive sufficient C 0 -optimality conditions under generic assumptions, see [9].

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 becomeṡ The associated pseudo-Hamiltonian reads 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: and denoting H X , H Y the Hamiltonian lifts: the Poisson bracket reads: Differentiating twice p(·) · G(x(·)) with respect to the time t, one gets: Lemma 2.20. Singular extremals (z(·), u(·)) are solution of the following equations: 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: 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:

The case of dimension 3 applied to the Zermelo problem
Consider the following affine control system: 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: where we denote: Moreover, let us introduce D := det(G, [G, F ], F ). In the vortex problem, one has: In our problem with the Goh extension, one orients p(·) using the convention of the maximum principle: and we introduce the following: Under assumptions (B1), (B2) and (B3), a geodesic is called: Note that the condition DD ≥ 0 amounts to the generalized Legendre-Clebsch condition ∂ ∂u and according to the higher-order maximum principle [27], this condition is a necessary (small) time minimization condition. See [9] 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.23. 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.18 for the projection of x(·) on the (x 1 , x 2 ) plane. In the exceptional case, the reference geodesic is C 0 -time minimizing.
2.5. Influence of the circulation 2.5.1. Influence of the circulation on the drift Then, we have (noticing that if u ∈ U, then −u ∈ U): This leads to introduce the following definition.
Remark 2.25. Note that if the drift could have been weak at any point x ∈ M , then we would have been in the Finslerian case [3] 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.26. 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.  Figure 2. 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.27. 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( and since p(t) = 0, the result follows.
Remark 2.28. 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: 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: 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: The normal extremals solution of the PMP are parameterized by the set 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 .

Resolution of the shooting equation
We introduce the shooting mapping where x f is the target and exp is the exponential mapping defined by (13). The shooting mapping is defined for 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 [12], 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 [12,17] 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.  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. (17). 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 α . Since one has similar numerical results in the strong drift case, we make the following conjecture:  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 curve is the graph of the smallest singular value σ min (t) for one α ∈ Θ(r 0 , µ) ∩ ∪ N i=1 {α i } and for t ∈ [0 , 60] on the left subgraph, and for one α ∈ Λ(r 0 , µ) ∩ ∪ N i=1 {α i } on the right subgraph. In this case, we stop the numerical integration when r(t) ≤ 10 −2 which explains why the minimal value of σ min is around 10 −5 .

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 (15), one obtains: where λ 1 := (2µc + p θ )p θ , λ 2 := 2(µp θ ) 2 , λ 3 := µc + p θ and λ 4 := µ 2 p θ . The aim is to get a polynomial one dimensional distribution to integrate, so we use the time reparameterization: Denoting by the derivative with respect to s, the system is written: We use Poincaré compactification where the sphere is identified to x = 1 and this leads to the following dynamics: which can be projected onto the 3-sphere r 2 + θ 2 + p 2 θ + x 2 = 1, up to a time reparameterization. Integrating the one-dimensional distribution from the system (20), 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: (20), one can obtain s = f (r), θ = g(r). This gives the integration and one gets a one-dimension foliation on the 3-sphere.

Using quadratures from
Let us end this section recalling the following. Letq = H(q), q = (x 1 , x 2 , x 3 , x 4 ), be an homogeneous polynomial vector field of degree k: H(λq) = λ k H(q). It can be projected onto the 3-sphere, the geometry projection being represented on the following figure: We set R 2 := x 2 1 + x 2 2 + x 2 3 + x 2 4 , v := q/|q| and up to reparameterization,q = H(q) decomposes into: It can be computed using projective charts. Consider: together with the vector field ∂/∂θ associated to the symmetry which respect to θ-rotation: We use the projective coordinates y 1 = r, y 2 = θ/r, y 3 = p r /r and y 4 = x/r so that the equatorial plane x = 0 is given by y 4 = 0. Reparameterizing, this leads to the dynamics: while θ = 1 projects ontoẏ 2 = 1,ẏ 3 =ẏ 4 = 0. Denote by G 1 and G 2 the corresponding vector fields, one has: which is colinear to G 2 . Again the two dimensional foliation on the the 3-sphere defined by ω(G 1 ) = ω(G 2 ) = 0 is given by the one form ω associated to the (r, p r ) dynamics represented as:ẏ 3 = λ 1 y 4 4 − λ 2 y 6 4 − y 2 3 ,ẏ 4 = −y 3 y 4 . To illustrate the integrability properties, we only give the formula in the abnormal case for the r component:

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 subsysteṁ 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 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), µ).
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) 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 α .
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 [22]. 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 [22] 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 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: Writing dt = 4µ 2 + r 2 4µ 2 − r 2 dr, then by integration we have with a constant c := ln K. Introducing 2µρ 2 := r, ρ ∈ [0 , 1], then one can define the function 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 [22,34]. 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 : 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 [22]. Since F (x, t + 1) = e F (x, t), then, R 2 defines a foliation on the annulus + 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: Example. A simple example which is detailed in [22] 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. On the left subgraph of Figure 12    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: 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 [22].
Let us exhibit now a function f from the Vortex application that satisfies the condition (7). We fix µ > 0. Coming back to eq. (24), one can write Integrating we have: with a constant c := ln K. Introducing for r ∈ [0 , 2µ]: 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 13. 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 14 shows the foliation on the punctured disk of radius 2µ in the (x 1 , x 2 )-plane of the Vortex application.   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(·),θ(·), pr(·), pθ) := (r(·), 2θ 0 − θ(·), −p r (·), p θ ). Then,z(·) is solution oḟ This means thatz(·) 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(·), θ(·)) andz 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: 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 15.
where V is the value function defined by (7). 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: 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: Let us fix ε > 0 and consider x ∈ B(x f , η ε ). By definition of the value function and by (26), we have 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 ceases 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: 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: Then, the cut locus from x 0 is simply defined by: 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: 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).
We have the following similar result: 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.
The splitting line is then given by solving F split = 0 since we have 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 16, we can see that the wavefront 3 W(x 0 , 3.5) 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 Figure 17 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. Let us describe the evolution of the wavefronts according to Figure 17. 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 18, 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 16.
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.29, 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:   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 19 represents four balls that we can group into three categories (see also Figure 17): • 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 19 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 19 is of type B. One can see on Figure 20, 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 19 represents a ball of type C. One can see on Figure 21 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 22.

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 Figure 20. Three balls with a zoom around the vortex (red dot) at times t ∈ {2.86, 2.88, 2.9}. Compare to Figure 18. For t = 2.9 the ball B(x 0 , t) is of type C and has a small hole with a fetus shape.  Figure 21. Spheres S(x 0 , 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.
of a strong current at the initial condition, in particular in relation with the abnormal extremals. A possible further extension concerns the case of several vortices by analogy with the N-body problem.
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. Figure 22. 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.