Integrability methods in the time minimal coherence transfer for Ising chains of three spins

The objective of this article is to analyze the integrability properties of extremal solutions of Pontryagin Maximum Principle in the time minimal control of a linear spin system with Ising coupling in relation with conjugate and cut loci computations. Restricting to the case of three spins, the problem is equivalent to analyze a family of almost-Riemannian metrics on the sphere S2, with Grushin equatorial singularity. The problem can be lifted into a SR-invariant problem on SO(3), this leads to a complete understanding of the geometry of the problem and to an explicit parametrization of the extremals using an appropriate chart as well as elliptic functions. This approach is compared with the direct analysis of the Liouville metrics on the sphere where the parametrization of the extremals is obtained by computing a Liouville normal form. Finally, an algebraic approach is presented in the framework of the application of differential Galois theory to integrability.


1.
Introduction. Over the past decade, the application of geometric optimal control techniques to the dynamics of spin systems with applications to Nuclear Magnetic Resonance (NMR) spectroscopy and quantum information processing [17] has been an intense research direction. In particular, a series of articles focus on the time optimal control of a linear chain of spins with Ising couplings [14,23]. Using Hamiltonian Formalism and an adapted rotating frame, the control system is defined by H = H d + H c , where H d is the internal Hamiltonian, H d = i,j J ij I iz I jz , with J ij representing the coupling between the spins, I kα = 1 ⊗ ... ⊗ I α ⊗ ... ⊗ 1, I α , α ∈ x, y, z being the Pauli matrix and the system is controlled by external radio-frequency pulse on resonance to each spin defining H c = i (u i1 I ix + u i2 I iy ).
Restricting to the case of three spins, the objective of this article is to provide the preliminary work to compute the optimal solutions parametrized by Pontryagin Maximum Principle. We here focus on the integrability aspects of the problem by using three different approaches.
A first point of view which already appears in the pioneering work [14], see also [9] in a different context, consists in lifting the problem on a sub-Riemannian invariant (SR-invariant) problem on SO(3) that depends on a parameter k representing the ratio of the coupling constants J 12 , J 23 between the spins. Such a metric is a limit case of invariant Riemannian metrics on SO(3), the so-called Euler-Poinsot rigid body problem in mechanics. Using the seminal work in [13], we define a chart that identifies locally SO(3) to S 2 × S 1 which enlightens the geometry of the problem and leads to an explicit computation of the extremals using elliptic functions.
Another approach consists in integrating the system directly on S 2 . In this context the problem is equivalent to analyze a family of 2D− Liouville metrics on S 2 with an equatorial singularity. The integrability properties is equivalent to the calculation of the Liouville normal form [4,5] using the additional first integral. In our case it corresponds to the Hamiltonian of the round metrics on S 2 , induced by the Casimir function on SO (3). This point of view is very important to analyze the optimality properties related to the conjugate and cut loci of the metrics, indeed it is related to similar calculations on Liouville surfaces that generalizes the case of ellipsoids [10,11].
Finally, the third approach consists in using our problem as a test-bed platform to apply the algebraic framework of Galois differential theory in integrability [18] to compute the solutions. First, the optimal control is calculated using the Jacobi elliptic functions and inserted in the equations. This reduces the computations to the integration of a time-depending linear equation whose coefficients are expressed in terms of the Jacobi elliptic functions. The Picard Vessiot extension and the associated Galois group are computed to parametrize the extremal solutions. Different computations were done independently but the parametrizations are compared in the conclusion. Also we briefly discuss the application to the computations of the conjugate and cut loci, following the method in [7] combining geometric analysis on Liouville surfaces and numerical computations. More details can be found in [12].
2. The mathematical model. In this paper, we detail the presentation of the problem in the case of three spins, but the problem can be easily generalized to a chain of n spins. We follow the presentation of [14,22,23].
We introduce the spin 1/2 matrices I α related to the Pauli matrices by a 1/2 factor. Such matrices satisfy: The Hilbert space L of the system is the space formed by the tensor product of the three two-dimensional spin 1/2 Hilbert space. Assuming a single input system, the Hamiltonian of the system decomposes into: We consider the time evolution of the vector X = ( denoting the expectation value. To compute the dynamics, we introduce a 8 × 8 matrix ρ ∈ L, called the density matrix, which satisfies: Using the definition of the expectation value of a given operator: O = Tr(Oρ), one gets: Hence, we deduce that: By rescaling the time by a factor J 12 , it becomes: Similar computations lead to the evolution of X given by: The optimal control problem is to transfer in minimum time (1, 0, 0, 0) to (0, 0, 0, 1) . It is an intermediate step to realize the transfer in minimum time from I 1x to I 3x . Indeed, it connects the first spin to the third one by controlling the second spin.
Introducing the following coordinates: the system becomes: where r = (r 1 , r 2 , r 3 ) ∈ S 2 and u 1 = −k sin(α), u 3 = − cos(α) are the components of the control. In those coordinates, the minimum time problem is equivalent to determine the fastest transfer on the sphere from (1, 0, 0) to (0, 0, 1).

Connection with invariant metrics on SO(3) and integration.
3.1. Lifting procedure. A first approach to analyze the optimal control problem and parametrize the extremals consists in lifting the problem onto SO(3). We introduce the matrix R(t) = (r ij (t)) on SO (3) where the third row is identified to the unit vector r(t) defined previously: r 31 = r 1 , r 32 = r 2 , r 33 = r 3 , and we consider the right-invariant control system: where the last column of R describes the evolution of the vector r. Our optimal control problem can then be stated as: min u(.) T 0 (I 1 u 2 1 (t) + I 3 u 2 3 (t))dt for the right-invariant control system with the following boundary conditions: which consist in steering the third axis of the frame R from e 1 to e 3 , where (e i ) is the canonical basis of R 3 . Similarly, it can be transformed into a left-invariant control problem to use the geometric framework and the computations in [13] : with the corresponding boundary conditions. This defines a left-invariant SR-problem on SO(3) depending upon the parameter k 2 = I 1 /I 3 . Upon an appropriate limit process I 2 → +∞, this is related to the Euler-Poinsot rigid body motion [2] : T 0 (I 1 u 2 1 (t) + I 2 u 2 2 (t) + I 3 u 2 3 (t))dt which is well-known model for left-invariant metrics on SO(3), depending on two parameters e.g. the ratios I 2 /I 1 , I 3 /I 1 . Two special cases are: • The bi-invariant case I 1 = I 2 = I 3 where the geodesics are the rotations of SO(3). • The case of revolution where I 1 = I 3 .
3.2. Extremal equations and Integration. The optimal solutions to our problem can be parametrized by the Pontryagin Maximum Principle [19], and thanks to the explicit formula given in [13], the solutions can be computed in both the Riemannian and the sub-Riemannian cases using elliptic functions.
We introduce: that satisfy the Lie brackets relations: Consider now the following optimal control problem on SO (3): where the I i 's are the principal momenta of inertia of the body. The extremal equation will be derived using appropriate coordinates. Let λ be an element of T * R SO(3) and denote H i = λ( A i (R)), i = 1, 2, 3 the symplectic lift on the vector fields A i . The pseudo-Hamiltonian associated to the problem takes the form: The extremal control is computed using the relation ∂H ∂u i = 0, and we obtain Plugging this expression for the u i into H, we get the Hamiltonian : The SR-case is obtained by setting u 2 = 0 which corresponds to take I 2 → +∞, and leads to the Hamiltonian : The evolution of the vector H = (H 1 , H 2 , H 3 ) is given by the Euler equation: where {,} denotes the Poisson bracket. Using the relation between Poisson and Lie brackets: , we obtain the Euler-equation: • SR-case (I 2 → +∞) The extremals equations are a classical example of (super) integrable system which is a consequence of the following proposition [13]. To provide details on the quadratures we introduce the following.
• r is the third row of the matrix R, • Φ 1 , Φ 2 , Φ 3 are the Euler-angles computed with the convention:

Useful formulas.
We recall the following: Hamiltonian using Euler-angles. Expressed in terms of the Euler angles the Hamiltonian for the Euler Poinsot rigid body motion takes the form: where p i is the canonical impulse associated to Φ i , i = 1, 2, 3. Observe that Φ 1 is a cyclic variable. As previously, the SR-case is obtained by taking I 2 → +∞ .
In both cases we have the following crucial proposition [13].
1. The angles Φ 2 and Φ 3 can be obtained from the relations: We have, with r = (r 1 , r 2 , r 3 ) ∈ S 2 (third row of R): While the Euler equation can be integrated using H n and the Casimir function , the angle Φ 1 can be computed by quadrature using [13]. Proposition 3. In the invariant case, Φ 1 is solution of the equation This leads to an uniform integration procedure in the invariant case using elliptic functions that we detail in the SR-case.

Integration of the Euler equation in the SR-case. We fix the level set for the
2 , and we introduce the angle β as follows: Using the Euler equation, we deduce that β is solution of the pendulum equation: We introduce ν = 2β, and we obtain the equation: We can assume I 3 > I 1 and I 3 = 1 and use [16]. We define and according to the constant C we have two types of generic solutions.
• Oscillating solutions where m is the modulus defined by C and we denote 4K(m) the period of sn, and the control is given by : We can deduce the following proposition.

Proposition 4.
Setting ω 2 = 1 k 2 − 1, the components of the control are given by: • In the oscillating case: • In the rotating case: To compute Φ 1 , we need the elliptic integral of the third kind See [16] for the result in the Riemannian case and [8] in the SR-case. 4. Direct integration on S 2 using Liouville theory.

4.1.
2D-Riemannian metrics whose geodesics flows are integrable by means of linear and quadratic integrals in momenta. We first recall some results from [4,5].
4.1.1. Linear Case. Let g be a real analytic Riemannian metric g(x, y) = a(x, y)dx 2 +2b(x, y)dxdy + c(x, y)dy 2 on a surface M and assume that the extremal flow possesses a non zero linear (in momenta) integral F . Then, there exists local coordinates u and v in which the metric has the polar form du 2 + m(u)dv 2 and p v is a first integral (Clairaut relation). This case is called the case of revolution.

4.1.2.
Liouville case. If the metric g admits an additional first integral F quadratic in momenta, the surface M is called a Liouville surface. This case is more intricate, and we present in details the algorithm to compute a normal form to integrate the extremal flow. First, we introduce isothermal coordinates (x, y) such that the metric takes the form: g = λ(x, y)(dx 2 + dy 2 ). If we denote p x , p y the adjoint variables, the first integral is given by: According to [4] this mapping R is holomorphic.
which preserves the isothermal form and the orientation satisfies the Cauchy-Riemann relations: Expressing F using the (u, v) coordinates, we obtain: ). An easy computation provides: where φ = (ϕ u + iψ u ). We choose the change of coordinates such that S = 1. Hence, we must solve the equation In the new coordinates, the metric takes the Liouville normal form To integrate, we use [5], Theorem 6.5. 1. The equations of the extremals for the Liouville metric can be written as: INTEGRABILITY FOR ISING 3-SPINS CHAINS 9 2. The extremals themselves are defined by the relation:

4.2.
Computations of the Liouville normal form. The Hamiltonian of the metric g can be written as: where H 0 is the Hamiltonian of the Grushin case (k = 1): and We can interpret S 2 as the homogeneous space SO(3)/SO (2) where SO (2) is the Lie subgroup leaving e 3 invariant. In this interpretation, the Casimir function corresponds to the bi-invariant case. On the homogeneous space, this defines the round sphere with constant curvature +1 whose metric in spherical coordinates takes the form: and corresponds to the Hamiltonian: A direct computation provides the following proposition.

4.2.1.
Integration in the Grushin case. This situation corresponds to a case of revolution and the integration is standard. The metric is already in the polar form and in our problem it is interesting to interpret the Grushin case as a deformation of the round case using the following homotopy: where m λ (ϕ) = sin 2 ϕ/(1 − λ sin 2 ϕ), λ ∈ [0, 1]. Indeed, it allows us to use the geometric framework developed in [6]. Except for the meridian solutions, the ϕ−variable is periodic and denoting ψ = π 2 − ϕ, the evolution of ψ is given by: We denote X = sin ψ and X + and X − the positive and negative roots of: . Introducing Y as X = X + Y , we have the following proposition.
1. The period of ψ is given by: 2. The ψ-variable in the normalized coordinates is given by: and we get
Moreover : Hence in the isothermal coordinates (x, y) : x − 2xyp x p y + (k 2 − y 2 )p 2 y and using the notation of section 4.1.2, one has : The solution of (4) is Case k<1 We take The dual variables are given by : Hence we have Therefore we have the following proposition : Proposition 7. In the case k < 1, the Liouville form is given by : Case k>1 Similarly using z = i √ k 2 − 1 sin w, we get : In both cases, the integration of the geodesics using Theorem 4.1 gives us elliptic integrals.
5. Algebraic techniques. In this section we present the techniques to integrate the equations using differential Galois techniques, see [3,18] for an introduction.

5.1.
Step 1 : Computation of the control. To illustrate the techniques it is sufficient to consider one case. Hence for simplicity we assume the following normalizations : I 3 = 1 and I 1 > I 3 .
The components of the control are given in proposition 4.

5.2.
Step 2 : Computation of the position. One must integrate the time depending linear equation : There are two possible controls We first begin with the controls (6). This is a system of differential equations, and its coefficients belong to the function field K = C(sn(ωt, m), cn(ωt, m), dn(ωt, m)). This field is a differential field, i.e. for any f ∈ K, ∂ t f ∈ K. Let us consider the resolvent matrix of equation (5), and the differential field extension L = K (a 1,1 , . . . , a 3,3 ) generated by the entries of the resolvent matrix.  (5). We say that equation (5) is solvable if there exist a tower of field extensions K 0 = K ⊂ K 1 ⊂ · · · ⊂ K n = L such that for all i = 0 . . . n − 1 where ∂ta a ∈ K i . The Galois group of equation (5) is the group of differential automorphisms of L stabilizing K, i.e. automorphisms σ of L such that ∂ t σ = σ∂ t and σ | K = id, see ( [21],1.25,1.26,1.27) for details. (5) is isomorphic to a Lie subgroup of GL 3 (C). The equation (5) is solvable if and only if G is virtually solvable, i.e. its identity component is solvable.

Proposition 8. The Galois group G of equation
Until now, all these concepts are theoretical. But we can classify all Lie subgroup of GL 3 (C) and make some test to know if the Galois group is a virtually solvable one, this has been done in all generality for third order operators in [20]. The first thing to notice is that equation (5) conserves the Euclidean norm. Indeed, r 2 1 + r 2 2 + r 2 3 is constant because the matrix associated to the differential system (5) is in so (3). So the Galois group G is in fact a subgroup of the (complex) group SO 3 (C).
The Lie subgroups of SO 3 (C) are well known. All strict subgroups of SO 3 (C) are finite or stabilize one axis. So the Galois group is either finite, either a finite extension of SO 2 (C), or SO 3 (C). All these groups are virtually solvable, except the last one SO 3 (C) which is simple. So, if we want to try to solve this equation, the first thing to do is to know if G stabilize one axis. Proof. We begin by a variable change z = sn(wt, m). The system now becomes The coefficients are now in , which is an algebraic field and so will be easier to manipulate. We denote abusively again L the Picard Vessiot field of this system. We now use the cyclic vector method to build a third order differential equation for which r 1 is solution: for any i ∈ N, ∂ i t r 1 is a linear form in r 1 , r 2 , r 3 . Thus these linear forms for i = 0 . . . 3 are not independent. So there exist a linear combination of ∂ i t r 1 , i = 0 . . . 3 equal to zero. The relation between these linear forms gives a differential equation for r 1 (z) (taking into account the relation between parameters w 2 = 1/k 2 − 1) Let us remark the following • We can express r 2 , r 3 as linear combinations in K 0 of derivatives of r 1 . Thus L is generated by the solutions of (8) and their derivatives. So L is also the Picard Vessiot field of equation (8). • The equation (8) has now rational coefficients due to simplifications. The relation w 2 = 1/k 2 − 1 allowed to only use the parameters m, k.
Now the expsols routine of package DEtools of Maple find r 1 (z) = √ 1 − m 2 z 2 as solution (which is dn before the variable change). So there is a solution in K 0 , and this implies that the Galois group stabilize an axis (and is equal to id on this axis). So we already know that the Galois group is a subgroup of O 2 (C), and thus that our equation is solvable. We can now reduce the order of the equation, posing r 1 (z) = √ 1 − m 2 z 2 f (z)dz. We obtain a second order differential equation, on which we use the Kovacic algorithm [15] to obtain the solutions for r 1 : This implies that the integral belongs to L. This is an integral over an element of K 0 , and this integral is an elliptic integral, [16]. So it does not belong to K 0 , and thus the Galois group G is of dimension at least 1. Now the only two remaining possibilities for G are SO 2 (C) C * or O 2 (C). In the first case, the Galois group would be connected, and so no algebraic extensions would be necessary to express the solutions. This is indeed the case, as we can put the square root √ k 2 m 2 + m 2 z 2 − k 2 − m 2 inside the integral in the exponential. So expressing the solutions only uses an exponential integral of an element of K 0 . So G C * .
Remark that the Galois group we computed is over the base field K 0 . However, the third order differential equation for r 1 has rational coefficients. So it makes sense to compute the Galois group over C(z). Over C(z), (1 − z 2 )(1 − m 2 z 2 ) is an extension of degree 2, and so the Galois group is G D ∞ = O 2 (C).

5.3.
Step 3 Simplifications of the solution. This expression is Liouvillian (finitely many integrals, exponentials and algebraic extensions), but is not the optimal one. Indeed, the Galois group has dimension one. This suggests that only one integral symbol should be necessary to write the solution. So a simple expression with no iterated integrals exists. Let us find it. As the Galois group over C(z) is D ∞ , there exists a basis of solutions of the form with F /F, G ∈ C(z). We first compute the symmetric power 2 of equation (8), i.e. the linear differential equation whose solutions are the products of two solutions of equation (8). We remark that this symmetric power should have F 2 in its space of solutions. We only have to compute the vector space of hyperexponential solutions to find "candidates" for F . The solutions are α + βz 2 . Thus F should be of the form F = α + βz 2 .
We can now inject F e √ Gdz in equation (8). This produces a large non linear differential equation in G(z). But G(z) also appears in the equation, and we can act the multivaluation of the square root G(z) −→ − G(z). This gives us another equation, and after simplification, we manage to obtain a linear equation Solving this equation with dsolve, we obtain only one solution in C(z) Injecting this in the original non linear equation gives Thus αk 2 + βk 2 − α = 0. Thus with α = k 2 , β = 1 − k 2 the equation is satisfied. So the solutions for r 1 are To conclude, let us do the same for the second equation. Proof. The proof is similar to the previous one for controls (6). We produce a third order differential equation for r 1 . This time, the expsols routine finds √ 1 − z 2 . Again, we search solutions under the form c 1 Gdz , and we find

5.4.
Step 4 Special values. For some special values of m, k, the Galois group simplifies. These exceptional parameter values can be seen as singularities or confluences on the formula of r 1 . Let us make an analysis for controls (6) (the other one is similar).
The main case appear when the integral in the exponential identically vanishes, i.e. k 2 m 2 − k 2 − m 2 = 0. There is one singularity for k = 0, but which is present in the original system and so is excluded. The other problematic cases are confluences: The elliptic integral simplifies, and it is no longer guaranteed that the exponential of it is transcendental. These are the cases k 2 = 1 and m 2 = 1.
The Galois group of (5) over K 0 is Proof. For k 2 m 2 − k 2 − m 2 = 0, the solutions can be obtained using a limiting process noting k 2 m 2 − k 2 − m 2 = and making a series expansion in The Galois group is then additive instead of multiplicative, G C + . If k 2 = 1, the equation becomes a first order equation, whose solutions are c 1 √ 1 − m 2 z 2 , so in K 0 , thus G = id. If m 2 = 1, the dsolve command returns an expression in which the only possibly transcendental term is ((z − 1)/(z + 1)) k 2 √ k 2 −1 . This term is algebraic if and only if 1/k ∈ E (and so G Z/nZ). If m = 0, the integral is no longer elliptic, and is given by The exponential of this expression is algebraic if and only if k ∈ E (and so G Z/nZ). The solution for c 1 = 0, c 2 = 1, c 3 = 0 is given by where Π denotes the elliptic integral of the third kind. The formula in the case c 1 = 0, c 2 = 0, c 3 = 1 is similar up to a sign. In the special case k 2 m 2 − k 2 − m 2 = 0, only an elliptic integral of the second kind is necessary, and in the even more special case of finite Galois groups, all the solutions are algebraic in z, and thus are algebraic in cos(ωt), sin(ωt) or exp(ωt). Note also that with the modulus m ∈]0, 1], k 2 m 2 − k 2 − m 2 = 0. 6. Conclusion.

Integrability methods
Euler angles and proposition 2 determine the third row of the matrix R(t) ∈ SO(3) using Jacobi elliptic functions (solutions of the pendulum equation) and a further quadrature gives the elliptic integral Π to parametrize the angle Φ 1 (the formulas in 3.2.1 gives exponential). The computation of all the rows of the matrix R(t) is associated to the computation of the Picard Vessiot extension. The method of section 5 consists in reducing the computation to a third order linear differential equation which is reparametrized using Jacobi elliptic functions to produce a specific algebraic solution corresponding to the third row of the matrix. Using this specific solution the third order equation is reduced to a second order equation which is solved using Kovacic algorithm. Concerning Liouville approach and theorem 4.1 the explicit computations of the integral (f (u) + a) −1/2 du and (g(v) − a) −1/2 dv superposes the Jacobi elliptic functions and the elliptic integral of the third kind.

Geometric applications
The geometric applications on the computations are the following. The Liouville approach is used to determine the conjugate and cut loci which is achieved in the next section. The Picard Vessiot extension is an important step to compute the conjugate and cut loci for the SR-invariant metrics on SO(3).
6.2. Optimality problem. Our work was a preliminary study to a complete analysis of the optimality problem which can be handled using the technical framework introduced in [7] combining geometric analysis and numerical simulations.
We use the following concepts. On the almost-Riemannian surface (M, g), the cut point along a geometric curve γ, projection of an extremal solution of the Maximum Principle, emanating from q 0 ∈ M is the first point where it ceases to be minimizing and we denote C cut (q 0 ) the set of such points forming the cut locus. The first conjugate point is the point where it ceases to be minimizing among the set of geodesics C 1 -close to γ emanating from q 0 and we denote C(q 0 ) the set of such points, forming the conjugate locus.
The computations of the conjugate and cut loci on surfaces is a very difficult problem and only recently [10] was proved the Jacobi conjecture : the conjugate locus of a non-umbilical point on ellipsoids has exactly four cusps. A consequence being that the cut locus is a segment. Also in [11] the result was extended to a wider class of Liouville surfaces which possess such simple cut and conjugate loci.
To generalize such computations in our case we can use two different techniques. First of all, we can make an explicit computations of the conjugate loci using our extremal parametrization or numerical simulations [7]. Secondly we can try to relate the simple structure of the conjugate loci to some monotonicity properties of Gaussian curvature or parametrized integrals related to extremals. Also in our case, one must take into account the existence of equatorial Grushin singularities of the metric.
We briefly present the method to determine the conjugate and cut loci of an one parameter family of Liouville metrics.
First, we have the following lemma.
The next step is to use the Grushin singularity resolution described in [6]. (This result is valid for every k.) Proposition 10. Near the equator point q 0 identified to 0, the conjugate and cut loci for the metric restricted to a neighborhood of 0 can be computed using the local model : dx 2 + dy 2 x 2 where x = 0 is the equator. The cut locus is a segment [− , ] minus 0 while the conjugate locus is formed by four symmetric curves of the form x = cy 2 , minus 0, and tangential to the meridian identified to θ 0 = 0. Note in particular that this computation allows to determine by continuation the conjugate and cut loci at an equatorial point. In particular we have, see also [9] for a similar result.
Proposition 11. For every k, the cut locus of an equatorial point is the equator minus this point.
Proof. A simple computation shows that the Gaussian curvature in each hemisphere is strictly negative. Hence there is no conjugate point for a geodesic starting from the equatorial point before returning to the equator. Due to the reflectional symmetry with respect to the equator two geodesics starting from an equatorial point intersect with same length when returning to the equator. This proves the result.
Finally, the simple structure of the conjugate and cut loci in the Grushin case is well known. Roughly spoken it can be simply deduced from the convexity of the period mapping p θ → T (p θ ) = 2π √ 1+p 2 θ given in proposition 6, see [6].
Proposition 12. In the Grushin case we have : 1. The cut and conjugate loci of a pole is the antipodal point.
2. The cut locus of an equatorial point q 0 is the whole equator minus this point, while the conjugate locus has a double heart shape, with four meridional singularities, two at q 0 and two cusps on the opposite meridian. 3. The conjugate locus of a point not a pole nor on the equator has only four cusps and the cut locus is a simple segment on the antipodal parallel.
The problem is to generalize this simple description to the case k = 1. This is the analog on the ellipsoid to generalize the case of revolution to a general ellipsoid.
From the theoretical point of view the main steps are the following : Step 1 : Liouville metrics on the sphere are classified using different notions of equivalence (Liouville equivalence, equivalence of geodesic flow or isometry of the metrics) see [5] and the relation between the case of ellipsoids and Euler-Poinsot rigid body motion is well understood.
Step 2 : The characterization of the metrics on the sphere satisfying the Jacobi conjecture on the ellipsoids : "The conjugate locus of a generic point has only four cusps and the cut locus is a segment" is the object of intense research activities [10,11] but the verification of the conditions is always a complicated task even in the case of revolution.
Due to this difficulty we adopt here a different point of view coming from [7] applied to the case of spins and to be compared to a similar analysis on the ellipsoid.

Geometric framework Ellipsoids
The general case is obtained by gluing the oblate and prolate cases and the conjugate and cut loci are numerically determined using this analysis. For a non umbilical point the conjugate locus has four cusps and the cut locus is a segment. For an umbilical point they shrink to a single point.
Spin case The sphere is interpreted by gluing the two hemispheres at the equator associated to the Grushin singularity of the metric.
The Grushin case k = 1 is described in Proposition.12 and corresponds to the oblate case for an ellipsoid (outside the equator).
For the case k = 1 the first step is to compute the set η of points where the first integral F is proportional to the Hamiltonian H. For a general ellipsoid we We represent on fig.1 such points in each hemisphere in relation with the level sets of the Liouville coordinates.
Numerically we evaluate the conjugate and cut loci of a point in η which shrinks into a single point, see fig.2.
We compute numerically the conjugate and cut loci of a point not on the equator nor in η. It is represented on fig.3 confirming the simple structure of the conjugate and cut loci. 6.3. Link with conjugate and cut loci computation for an invariant Riemannian on sub-Riemannian metrics on SO(3). Note also that a side effect of our computations is a first step towards the computation of conjugate and cut loci for the metrics on SO(3), which is an important challenge. Also it can be generalized to the case of simple groups of dimension 3.
Controlling the two spins 2 and 3, one gets the system : where k 3 = J23 J12 , k 5 = J34 J12 . Using the notations : x 2 = r 2 cos α 2 , x 3 = r 2 sin α 2 , x 4 = r 3 cos α 3 , x 5 = r 3 sin α 3 , x 6 = r 4 .  From the control parametrization, the control set is the surface M = {u 2 1 + u 2 2 + u 2 3 − u 2 1 u 2 3 = 1}. The control constraint can be written : v r, dr dt = g 2 r, dr dt + g 4 r, dr dt = 1 with g 2 , g 4 respectively quadratic and quartic with respect to the speed. The optimal control problem is a time minimal transfer. The candidate as minimizers can be computed using the maximum principle. One can observe that the optimal problem is not a SR-problem unless we consider small controls and the control constraints can be approximated with g 4 = 0. In this case the problem is associated to an invariant SR-problem on SO(4). An interesting question in this context is to generalize the integrability issues.