VORTEX PAIRS ON A TRIAXIAL ELLIPSOID AND KIMURA’S CONJECTURE

. We consider the problem of point vortices moving on the surface of a triaxial ellipsoid. Following Hally’s approach, we obtain the equations of motion by constructing a conformal map from the ellipsoid into the sphere and composing with stereographic projection. We focus on the case of a pair of opposite vortices. Our approach is validated by testing a prediction by Kimura that a (inﬁnitesimally close) vortex dipole follows the geodesic ﬂow. Poincar´e sections suggest that the global ﬂow is non-integrable.

1. Introduction. The equations describing the motion of N -point vortices on an ideal planar fluid were introduced in 1867 by Helmholtz [15] and described as an Hamiltonian system in 1876 by Kirchhoff [23]. Equations for point vortices on the two dimensional sphere were derived independently by I. Gromeka [13] and by E. Zermelo [38] and rediscovered in 1977 by Bogomolov [2]. In 1999 Kimura [22] studied all complete surfaces with constant curvature. Kimura conjectured that on any surface a pair of infinitesimally close opposite vortices would move along a geodesic. For the hyperbolic plane a recent study was carried out by Montaldi [29].
In 1980 D. Hally [14] wrote the equations for the point vortex dynamics on a simply connected compact surface (i.e, surfaces diffeomorphic to spheres) using isothermal coordinates. For such a surface with metric ds 2 = h 2 (z, z)|dz| 2 where z ∈ C ∪ ∞ represents stereographic coordinates on the sphere, Hally's equations arė where Γ k represents the k-th vortex intensity and z k its coordinates. The notation (z, z) should be familiar to the reader, meaning Rez = (z + z)/2, Imz = (z − z)/(2i).
The topological constraint of the surface being compact imposes Bounded domains inside curved surfaces, simply or multiply connected, can be studied on its planar image via Theorem 2 of C.C.Lin's classical paper [27]. Presently, there are powerful methods to produce conformal mappings to the unit circle or the unit circle with circular holes [9].
Recent work. The case of a compact Riemann surface S of any genus endowed with an arbitrary metric was addressed by Boatto and Koiller [1]. The constraint (2) can be relaxed. In fact, for compact surfaces, the Green function G(s 1 , s 2 ) of the Laplace-Beltrami operator governing the vortex-vortex interactions, also encodes a background counter-vorticity, uniformly distributed with respect to the metric. The Robin function (desingularization of G) accounts for the self interactions. For vortices in the round sphere there is a sizeable literature (for a fairly complete list see [1]). We now review some work on vortices moving on surfaces with nonconstant curvature. In 2008 Castilho and Machado [7] wrote Hally's equations as an Hamiltonian system, with Hamiltonian function and symplectic form Ω = N n=1 Γ n h 2 (z n , z n )dz n ∧ dz n .
Using perturbation theory they obtained first order approximations for Hally's equations for an ellipsoid of revolution x 2 R 2 + y 2 R 2 + z 2 R 2 (1+ ) = 1, for small values of . The ellipsoid's symmetry was used to reduce the dimension of the problem. In 2010 Kim [20] obtained the full equations for any ellipsoid of revolution. Several other surfaces of revolution were considered in [10]. Kimura's conjecture for vortex pairs (Γ 2 = −Γ 1 ) was first tested in [25]. As for numerical methods: San Miguel [28], used least-squares fitting to obtain discretized conformal mappings between ovaloids and the sphere. He integrated the vortex pair equations using the Gaussian collocation method. The first study on a genus 2 surface was done by C. Ragazzo [32]. Based on a relation between the Laplace-Beltrami Green function and the heat kernel, an algorithm is presented to determine the motion of a single vortex which is governed by Robin's function. The method is applied to compute the motion of a vortex on the Bolza surface, namely a constant curvature genus 2 surface whose fundamental domain is a regular octagon.
Summary of the paper. We study the motion of a vortex pair on the triaxial ellipsoid E 2 (a, b, c) : with 0 < a < b < c . In section 2 we review Jacobi's confocal conics coordinates λ 1 , λ 2 , with a ≤ λ 1 ≤ b ≤ λ 2 ≤ c that parametrize one octant of the ellipsoid. Jacobi also derived a system of conformal coordinates, but they become singular at the four umbilical points that belong to the ellipse with semiaxis a, c (y = 0), corresponding to λ 1 = λ 2 = b. In section 3 we recall the sphero-conical coordinates µ 1 , µ 2 for the unit sphere S 2 , that depends on three affine parameters I 1 < I 2 < I 3 , with I 1 < µ 1 < I 2 < µ 2 < I 3 .
In section 4 we construct a conformal map between the two surfaces, with the help of a simple, but useful lemma 4.1. The sphero-conical parameters I j are chosen in such a way that one octant of the ellipsoid is mapped exactly into one octant of the sphere under a common isothermal parametrization. There are two relations between the affine triples a, b, c and I 1 , I 2 , I 3 , given in terms of complete elliptic integrals. Moreover, each λ i is an elliptic function of their corresponding µ i , i = 1, 2. Dupin's lines of curvature of the ellipsoid are mapped into a topologically equivalent system of curves in the sphere. The umbilics The conformal factor is (λ 2 − λ 1 )/(µ 2 − µ 1 ), which is 0/0 at the umbilical points. We computed the limit in §5: it is equal to [ In section 6 we write the vortex pair equations on the ellipsoid and present our methodology to numerically integrate them. Composing the ellipsoid to sphere map with the stereographic projection from the sphere into the z-complex plane we get the conformal factor h(z,z) required for Hally's equation.
In section 7 our approach is validated by verifying Kimura's conjecture [22] about the relation between the dipole dynamics and the geodesic flow. We also compute exploratory Poincaré maps for the flow suggesting that it is chaotic. Directions for future research are presented in section 8. In Appendix A (following [1]) two proofs of Kimura's conjecture are outlined. Appendix B outlines Carlson's method for numerically computing elliptic integrals [6].
They satisfy Fig. 2. Clearly, the semiaxis extremes correspond to:

ADRIANO REGIS RODRIGUES, CÉSAR CASTILHO AND JAIR KOILLER
There are four umbilical points located in the middle ellipse (y = 0): corresponding to λ 1 = λ 2 = b. The following result is classical: The metric ds 2 induced by the embedding of the ellipsoid in R 3 is of Liouville type ( [4]) (note that the second term is positive).
In a famous paper Jacobi showed that the geodesics on the triaxial ellipsoid are integrable ( [17],1839). In §28 of his Vorlesungen ( [18], 1866), he presented a derivation using the (now called) Hamilton-Jacobi PDE, that separates using confocal quadrics coordinates 1 . In §28 Jacobi also constructed a local conformal map from the triaxial ellipsoid to the plane (pp. 215-217 of second edition). Jacobi's map was implemented recently in [31] and [19]. Isothermal coordinates (u, v) on an octant of E 2 can be constructed using elliptic integrals of the third kind Π (see [5]) We define the functions u = P (λ 1 ) (increasing) and v = Q(λ 2 ) (decreasing) by Proposition 2. The metric (11) in the ellipsoid (5) induced by its embedding in the euclidian space has, in the first octant, the isothermal coordinates are given by (13) and (15), The conformal map becomes singular (because the conformal factor vanishes) when λ 2 (v) = λ 1 (u) which occurs only for u = K 1 , v = K 2 , precisely the umbilical point.
The lines of curvature (of both families) in Fig.1 cannot be given a consistent direction along the middle ellipse 2 . In Fig.2, the top panel shows the transformation 1 The triaxial ellipsoid geodesics problem also separates in sphero-conical coordinates in §3 c, but we preferred to use Jacobi's confocal coordinates). 2 Source https://en.wikipedia.org/wiki/Umbilical_point from confocal coordinates (λ 1 , λ 2 ) to Jacobi's (u, v) on every octant. The following Proposition explains the middle and bottom panels. An underlined letter means that in the planar map the region appears flipped by 180 • but overlap perfectly in the ellipsoid. This observation certainly has not escaped to Jacobi, to whom it must have appeared so trivial that he did not even bothered to put in print.

Proposition 3. The real elliptic functions
obtained by inverting respectivelly (13) and (15), give rise to a double (branched) covering of the ellipsoid by a flat torus. The lattice has fundamental domains of sizes 4K 1 and 4K 2 . Each of the 16 rectangles of sizes K 1 × K 2 in the plane (u, v) corresponds to an octant: the ellipsoid is covered twice.
For the numerical work will not need to use their explicit formulas. Instead, we will be only computing elliptic integrals. Four rectangles in the plane surrounding a point marked U cover twice the sector formed by two octants with a common umbilical point. Rectangles with vertices corresponding to the umbilical points ±U 1 , ±U 2 with centers B (or −B) are mapped to sides y > 0 or y < 0 of the ellipsoid (5). In the next sections we will remedy the defect at the umbilical points.
Remark 1. In contradistinction with geodesics,which depends on the local metric, vortex motion depends on non local effects, which are encoded in the Green function of the Laplace-Beltrami operator of the metric [1]. An early attempt to use Jacobi's coordinates (u, v) in our numerical experiments for the vortex pair problem on the triaxial ellipsoid was not satisfactory. However, for initial conditions where the distance between the vortices was small, we show in §6 that the dynamics tends to move along a geodesic of the ellipsoid, as predicted by Kimura's conjecture. Therefore, in order to study point vortices on the full triaxial ellipsoid a global map from the ellipsoid to the sphere is needed.
Our construction of a conformal map from the ellipsoid to the sphere makes use of confocal coordinates in E and sphero-conical coordinates in S 2 , that has four "fake" singular points. The distribution of coordinate lines on the two surfaces correspond: in the ellipsoid they are the principal curvature lines, that have the umbilical points as singularities. The umbilics on the ellipsoid are mapped into the singular points of the sphero-conical coordinates. Hence the map will be global. 3. Sphero-conical coordinates (µ 1 , µ 2 ) on the sphere S 2 [3]. This coordinate system depends on three arbitrary parameters I 1 < I 2 < I 3 , that will be chosen so that one octant of the ellipsoid (5) gets mapped exactly into one octant of the sphere (6). This will permit to extend the conformal transformation between the whole surfaces, in such as way that the coordinate lines distributions correspond.
The parametrization of (γ 1 , γ 2 , γ 3 ) ∈ S 2 is . Taking µ 1 = µ 2 = I 2 one gets four distinguished points in the sphere that will be paired with the ellipsoid umbilics (10) in the conformal transformation between the surfaces. A quick derivation of this coordinate system on an octant of the sphere comes indirecly from the following problem (see [26]). Let A = diag(I 1 , I 2 , I 3 ). Diagonalize Ax, x , x ∈ R 3 , restricted to the subset x, γ = 0; that is, find extremals of Ax, x constrained to x 2 = 1 and to x, γ = 0. The extremals can be located by considering the function f (x) = 1 2 Ax, x , constrained to the subset of S 2 defined by ϕ −1 (0, 0) where Let µ and k be Lagrange multipliers. We look for solutions of That is Considering that x 2 j = 1 and For each choice of value µ, equation (22) represents an elliptical cone in the Euclidian space of the (γ 1 , γ 2 , γ 3 ). The intersection of these cones with the sphere γ, γ = 1, are curves that represent an orthogonal system of coordinates, since the extremal vectors (21), are extremals of the quadratic form Ax, x , and are also parallel to the gradient of (22). Equation (22) has two roots I 1 < µ 1 (γ, I 1 , I 2 , I 3 ) < I 2 , I 2 < µ 2 (γ, I 1 , I 2 , I 3 ) < I 3 that are explicitly obtained from a quadratic equation. Conversely, equations (19) for (γ 1 , γ 2 , γ 3 ) can be obtained by solving the linear system Proposition 4. The standard metric ds 2 of S 2 can be written in terms of the sphero-conical coordinates (19) with parameters (I 1 , I 2 , I 3 ) as For a derivation, see [3].

4.
Constructing the conformal map from E 2 (a, b, c) to the unit sphere. Two conformal maps from a closed simply connected surface to the sphere differ by Moebius transformations in S 2 . For the triaxial ellipsoid we found two references: Schering [35] in 1857 and Craig [8] in 1880. Both are quite intricate analytically, so we opted do to an ab initio construction, that may have its own interest due to its simplicity. Our map is equivalent to theirs. As this paper was being revised, we found a post in a cartography forum with a similar idea, combining Jacobi's projection of the y > 0 side of the ellipsoid, with a projection due to Goyou of half the sphere on a common rectangle in the plane 3 . Such a map will also makes explicit the (unique) complex structure in E 2 (a, b, c), via the global isothermic coordinates z,z obtained by stereographic projection of the sphere over the complex plane 4 .

4.1.
A simple lemma. The following is immediate: local coordinates on surfaces S andS, respectively. Assume that the respective metrics can be written as Then the correspondence defines a conformal map between the surfaces with conformal factor h given by In other words, we use the coordinate patch (ξ 1 , ξ 2 ) ∈ R = [0, r 1 ] × [0, r 2 ] as common isothermal parameters for the two surfaces, We apply this Lemma using S = S 2 with the sphero-conical coordinates (µ 1 , µ 2 ) andS = E 2 with Jacobi confocal coordinates (λ 1 , λ 2 ) with the corresponding metrics given by (24) and (11). It is VERY important that given a < b < c we chose I 1 < I 2 < I 3 such that condition (27) holds. On both coordinate systems and in each octant the coordinate curves meet the great circles perpendicularly, except at the point corresponding to the umbilic points at the ellipsoid.
Identifying the corresponding points of the two surfaces on the double coverings by the lattice in the complex plane w = u + iv, a global map from E 2 to S 2 results. w, w are common isothermal coordinates. The umbilical points in the ellipsoid (and the special points in the sphere) are ramification points of order 1/2.

Technical details.
Computing each function λ i (µ i ), i = 1, 2 requires one elliptic integral and one elliptic integrals inversion. Computing the parameters I 1 , I 2 , I 3 of the sphero-conical coordinates involves a nonlinear system involving two complete elliptic integrals of the first kind. Recall that the standard metric on S 2 written in the sphero-conical coordinates is given by (24), Define the functions with φ = arcsin Here F is the elliptic integral of the first kind [5] F Let P and Q as (13) and (15) respectively. The relations (28) become We must impose the condition (27) in Lemma 4.1 that assures that one octant of the sphere is mapped exactly over an octant of the ellipsoid. This amounts to Given a, b and c, finding I 1 , I 2 and I 3 satisfying (40). Note that both triples (a, b, c) and (I 1 , I 2 , I 3 ) can be considered as projective quantities. Let be the complete integral of the first type. From (34) and (36) it follows that Observing that k 1 and k 2 are complementary, we get After solving this system for k 1 , k 2 , the parameters I 1 , I 2 and I 3 are obtained from 5. Conformality at the umbilical points. The conformal factor is given by (32), with the coordinates implicitly related by (39). Umbilics correspond to λ 1 = λ 2 = b, and at the sphere we have µ 1 = µ 2 = I 2 in view of (40). Therefore we get a 0/0 indeterminacy. Let α, β > 0, with α + β = 1, and take Let's examine the limit Using l'Hôpital, we should investigate the two limits in the combination We will now show that both are equal to We do it for the first: where dP/dλ 1 is computed at λ 1 = P −1 S(µ 1 ). Therefore Let's now try to compute this limit as µ 1 → I 2 . Recall that P −1 S(I 2 ) = b. Pulling out (if we may) the factors that have a direct limit we get The limit in the right is, by a stroke of luck, √ L. The second limit is computed analogously and gives, seemingly by another stroke of luck, the same result, but his is indeed what we expect to happen by Riemann surfaces theory.
In passing, we have also shown that for µ 1 < I 2 < µ 2 , both close to I 2 , we have where a, b, c and I 1 , I 2 , I 3 are related by (44). We defined a global map from E 2 to S 2 by identifying the corresponding points on the two surfaces (ellipsoid and the sphere) on their double coverings by the same lattice in the complex w-plane (w = u + iv). Summarizing: Proposition 5. Global isothermal coordinates (z, z) on E 2 (except for the point corresponding to z = ∞) are obtained by stereographic projection from the γ-sphere to the complex z-plane, namely: Moreover, z = z(w) is a (complex) elliptic function with quarter periods K 1 , iK 2 .
6. The vortex pair on the triaxial ellipsoid E 2 . The case of a vortex pair is given when N = 2 and Γ 1 = −Γ 2 = Γ in Hally's equations (1). The symplectic form and Hamiltonian function are respectively Hamilton's equations are Let us now discuss how to do the actual computations. Let (z, z) denote the stereographic coordinates through the south pole S on S 2 \ {S}.
The ellipsoid metric in stereographic coordinates is given by (59) To simplify notation we define ) and r(z, z) ≡ 1 + zz.
(60) We also write λ i (z, z) = λ i (µ i (z, z)), i = 1, 2. Therefore, we get Proposition 6. The conformal factor of the ellipsoid over the z-complex plane is given by These determine equations (53). L z and M z , are given by that are obtained through implicit diferentiation of (39).
On the other hand, the geodesic equations in local coordinates (u, v) arë The first fundamental form coefficients are and we obtain the Christoffel symbols

ADRIANO REGIS RODRIGUES, CÉSAR CASTILHO AND JAIR KOILLER
Hence, the geodesic equations arė The numerical integration was performed as follows. Let v(0) and q(0) be initial conditions for some geodesic. If γ(t) is its projection over the ellipsoid, integrating numerically determine p 1 = γ( /2) and p 2 = γ(− /2). q 1 and q 2 are the initial conditions for the vortex dipole. The geodesics initial conditions are q(0) andṽ(0) whereṽ(0) is a π/2 positive rotation of v(0) suitably normalized. In the figures the dashed lines represent the geodesics and solid lines represent the vortices. In the top figures u is the horizontal axis and v is the vertical axis. See Figs. 3 and 4.
Exploratory Poincaré Sections were produced using Henon's method [16]. We depicted the stroboscopic positions of one of the vortices. See Figs. 5 to 7.
8. Final comments. In this paper we make a first study about point vortices moving on the surface of a triaxial ellipsoid. We focused on the case of a pair of opposite vortices. Our methodology was validated by testing Kimura's conjecture on close by pairs. The numerical experiments we presented on global behavior are exploratory, and we are preparing a more thorough study for a sequel paper.
1. Domains (simply or multiply connected) in the ellipsoid can be mapped to topologically equivalent domains in the plane. Then Theorem 2 of [27] as geometrized in [1] allow to study confined vortex motion in the planar image. 2. Equilibria and their stability. It is geometrically evident that vortex pairs placed at the ends of the principal axis are equilibria. Which of the three configurations are stable? 3. One of the referees suggested superposing Poincaré sections of the ellipsoid geodesics system with the sections of the vortex pair. Can the vortex pair system be regarded as a KAM perturbation of the integrable geodesic system on the triaxial ellipsoid when the vortices are sufficiently close? For a surface S, a numerical construction for the symplectomorphism S × S − diagonal ≡ T * S − zero section is needed (see [1]) but this is not a trivial task. 4. What would be the behavior of a closeby vortex pair such that its midpoint is initially on a geodesic passing through an umbilical point? It is known that the geodesic path will pass successively through the opposite umbilical point, approaching the middle ellipse as the time t → ±∞. That is, for Jacobi's integrable geodesic problem the middle ellipse is a periodic orbit with coinciding stable and unstable manifolds. Can transversality be shown for the vortex pair system regarded as a perturbation of the geodesic system? How about vortex pairs placed at opposite umbilical points? We conjecture that they will traverse a periodic orbit passing through the other pair. 5. We plan to pursue a more thorough investigation of Poincaré sections in the future. One of the referees suggests to make parameters move away from circular symmetry (two equal axis). Chaotic regions should became more and more visible (e.g. take a=1, b=1, 1.1, 1.2,1.3,..., c=9). 6. The single vortex problem on a compact surface S (moving on an uniform countervorticity field) has the Robin function R as Hamiltonian   result (see [36]) holds only for genus zero: where K is the Gaussian curvature function of the surface and ∆ the Laplace-Beltrami operator of the metric. For the triaxial ellipsoid it should be possible to solve this Poisson equation using the confocal coordinates. 7. Three point vortices on the sphere are integrable due to the SO(3) symmetry and were studied in [21] and [34]. How is the motion of three vortices affected on a nearly spherical oblate or prolate ellipsoid with two equal axis? 8. We used a Runge-Kutta method in the simulations. A symplectic integrator would be more adequate for very long times. Indeed, one could apply directly a symplectic integrator on S 2 ×S 2 , instead of the plane, just using the conformal map from the ellipsoid to the sphere. A new numerical method for point vortices in the sphere was proposed in [37]. Numerical methods for point vortices on surfaces could be produced using techniques of discrete differential geometry, implementing numerically the conformal maps and Green functions of a discretized Laplacian [12].