B Bonnard

O Cots

J.-C Faugère

A Jacquemard
email: alain.jacquemard@u-bourgogne.fr

J Rouot

M Safey

E L Din

T Verron

ALGEBRAIC-GEOMETRIC TECHNIQUES FOR THE FEEDBACK CLASSIFICATION AND ROBUSTNESS OF THE OPTIMAL CONTROL OF A PAIR OF BLOCH EQUATIONS WITH APPLICATION TO MAGNETIC RESONANCE IMAGING *

Keywords: 1991 Mathematics Subject Classification. 49K15, 14Q20, 81Q93. July 4, 2017 optimal control, singular trajectories, symbolic computation, Bloch equations, Magnetic Resonance Imaging

The aim of this article is to classify the singular trajectories associated with the optimal control problems of a pair of controlled Bloch equations. The motivation is to analyze the robustness of the optimal solutions to the contrast and the time-minimal saturation problem, in magnetic resonance imaging, with respect to the parameters and B 1 -inhomogeneity. For this purpose, we use various computer algebra algorithms and methods to study solutions of polynomial systems of equations and inequalities which are used for classification issues: Gröbner basis, cylindrical algebraic decomposition of semi-algebraic sets, Thom's isotopy lemma.

Résumé. L'objectif de cet article est de classifier les trajectoires singulières associées aux problèmes de contrôle optimaux d'une paire d'équations de Bloch contrôlées. La motivation est d'analyser la robustesse de la solution optimale du problème de contraste et de multisaturation en temps minimal, en imagerie par résonance magnétique nucléaire, par rapport aux paramètres et les inhomogénéités B 1 . On utilise le calcul symbolique pour étudier les solutions de systèmes polynomiaux d'équations et d'inéquations dans les problèmes de classification : base de Gröbner, décomposition algébrique cylindrique des ensembles semi-algébriques et le lemme d'isotopie de Thom.

INTRODUCTION

Optimal control algorithms were introduced in Nuclear Magnetic Resonance (NMR) to improve the control field very recently [START_REF] Conolly | Optimal control solutions to the magnetic resonance selective excitation problem[END_REF] and at the end of the nineties, new methods appeared in optimal control of NMR systems both from the analytical and numerical points of view [START_REF] Khaneja | Sub-Riemannian geometry and time optimal control of three spin systems: quantum gates and coherence transfer[END_REF][START_REF] Skinner | Application of optimal control theory to the design of broadband excitation pulses for high-resolution nmr[END_REF]. More recently, under the impulse of S. Glaser, the combination of geometrical optimal control based on the Maximum Principle [START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF] and related numerical algorithms (gradient methods [START_REF] Khaneja | Optimal control of coupled spin dynamics: design of NMR pulse sequences by gradient ascent algorithms[END_REF], shooting and continuation methods [START_REF] Cots | Contrôle optimal géométrique : méthodes homotopiques et applications[END_REF]) lead to sophisticated results starting from a complete solution to the time-minimal saturation of a single spin [START_REF] Lapert | Singular extremals for the time-optimal control of dissipative spin 1 2 particles[END_REF], and applications to the contrast problem in magnetic resonance imagery (MRI), see [START_REF] Bonnard | Singular trajectories and the contrast imaging problem in nuclear magnetic resonance[END_REF][START_REF] Lapert | Exploring the Physical Limits of Saturation Contrast in Magnetic Resonance Imaging[END_REF]. They are the basis to the numeric computations of robust optimal controls with respect to B 0 and B 1 inhomogeneities and were validated vert recently by in vitro (Figure 1) and in vivo experiments [START_REF] Bonnard | Geometric optimal control of the contrast imaging problem in nuclear magnetic resonance[END_REF][START_REF] Van Reeth | Optimal control design of preparation pulses for contrast optimization in {MRI}[END_REF][START_REF] Sugny | Optimal control theory for applications in magnetic resonance imaging[END_REF] (Figure 2).

The mathematical model which very accurately is suitable for analyzing such problems is to consider an ensemble of spins, each spin being described by a magnetization vector M = (M X , M Y , M Z ) in the laboratory frame 0XY Z whose evolution satisfies the so-called Bloch equation

dM dt = γM ∧ B + R(M), (1)
where γ is the gyromagnetic ratio, B(t) is the total magnetic field applied to the system which decomposes into

B(t) = B 0 + B 1 (t)
where B 0 is a strong constant magnetic field oriented along the Z axis, and B 1 (t) is the control RF-field in the transverse (X,Y )-plane. The R(M) term represents the dissipation, of the form

R(M) = M X T 2 , M Y T 2 , M Z -M 0 T 1
where M 0 is the equilibrium magnetization, which can be normalized to 1 using the rescaling M → M/M 0 , and T 1 , T 2 are the relaxation parameters which are the chemical signatures of the observed species. The control components are denoted u(t) = -γB Y , v(t) = -γB X and up to a time rescaling, one can impose the control bounds:

u 2 + v 2 = 1.
The Bloch equations can be written in a rotating frame 0xyz with M(t) = S(t)q(t), q = (x, y, z), S(t) = exp(ωtΩ z ), Ω z = 0 -1 0 1 0 0 0 0 0 as:

d dt   x y z   =   -1/T 2 -∆ω u 2 ∆ω -1/T 2 -u 1 u 2 u 1 -1/T 1     x y z   +   0 0 1/T 1   (2)
where ∆ω = ω 0 -ω is the resonance offset, ω 0 = -γB 0 is the resonance frequency and the RF-control field is represented as u 2 = u cos(ωt) -v sin(ωt) u 1 = u sin(ωt) + v cos(ωt) [START_REF] Bonnard | Feedback equivalence for nonlinear systems and the time optimal control problem[END_REF] which preserves the control bound u 2 1 + u 2 2 ≤ 1. Finally, in the moving frame the system takes the normalized form

           dx dt = -Γx -∆ωy + u 2 z dy dt = ∆ωx -Γy -u 1 z dz dt = γ(1 -z) + u 1 y -u 2 z. (4)
The resonant case is when ∆ω = 0 which will be assumed in the sequel of this article. The system is then written as

dq dt = F(q) + u 1 G 1 (q) + u 2 G 2 (q), u 2 1 + u 2 2 ≤ 1
Thanks to the symmetry of revolution along the Z axis related to the B 0 -orientation, in many cases it is physically relevant to consider the situation with u 2 = 0 and the spin system is a 2-dimensional system, with q = (y, z) and it is simply written as dq dt = F(q) + uG(q), q = (y, z), |u| ≤ 1.

(

Also note that the Bloch ball |q| ≤ 1 is invariant for the dynamics provided that the parameters satisfy the physical constraint 2Γ ≥ γ ≥ 0.

Our study will concern the following optimal control problems. First of all, consider a single spin system described by [START_REF] Bonnard | Algebraic geometric classification of the singular flow in the contrast imaging problem in nuclear magnetic resonance[END_REF]. The saturation problem is to drive the magnetization vector from the north pole of the Bloch ball N = (0, 1) to the center O = (0, 0). A major result in NMR was to compute the time-minimal solution, using geometric optimal control based on the maximum principle [START_REF] Lapert | Développement de nouvelles techniques de contrôle optimal en dynamique quantique : de la Résonance Magnétique Nucléaire à la physique moléculaire[END_REF].

According to this principle, a time-minimal solution is given by the equations:

ż = H F (z) + u H G (z), z = (q, p) (6)
H F and H G being the Hamiltonians p • F(q) and p • G(q), and moreover the optimal control is given by the maximization condition

H F (z) + uH G (z) = max |v|≤1 H F (z) + vH G (z) (7)
and it is a concatenation of bang arcs δ + , δ -where u = +1 or u = -1 and the so-called singular arcs ; they are obtained by solving H G (z(t)) = 0. A straightforward computation shows that singular arcs are located on the set where G and [G, F] are colinear and are formed by the two lines y = 0 (axis of revolution) and z 0 = -γ 2(Γ-γ) and the interesting situation is when 2Γ > 3γ so that this horizontal line intersects the Bloch ball. More precisely, following [START_REF] Lapert | Développement de nouvelles techniques de contrôle optimal en dynamique quantique : de la Résonance Magnétique Nucléaire à la physique moléculaire[END_REF] we recall the following:. Proposition 1.1. Provided that 2Γ > 3γ, the time-minimal solution to the saturation problem in the experimental situation is of the form δ + δ SH δ B + δ SV where δ SV and δ SH are vertical and horizontal singular arcs and the intermediate arc δ B + is a bang arc connecting the two lines, called a bridge and is related to the property that the singular control as calculated along the singular line is such that |u S | → ∞ when y → 0.

(Note: in experimental situation, the line (δ SH ) is accessible from the North pole).

All the information about this study is reported on Fig. [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF]. Note that in the case 2Γ ≤ 3γ, the standard inversion sequence δ + δ SV is time-minimal. In MRI, we attribute to |q| a level of gray such that |q| = 1 corresponds to white while and saturation |q| = 0 to black. This led to consider the following optimal problems, taking a pair of spins, each governed by [START_REF] Bonnard | Algebraic geometric classification of the singular flow in the contrast imaging problem in nuclear magnetic resonance[END_REF], with respective parameters (Γ 1 , γ 2 ), (Γ 2 , γ 2 ) and controlled by the same RF-field u (coupling the dynamics).

Problem 1. Saturation of both spins in minimum time. Note that another problem is to saturate a couple of spin systems with the same parameters but taking into account B 1 -inhomogeneity and this leads to consider a pair of systems (5):

       dq 1 dt = F(q 1 ) + uG(q 1 ) dq 2 dt = F(q 2 ) + u(1 -ε)G(q 2 ) |u| ≤ 1
where ε is a small rescaling parameter related to the control field inhomogeneity.

Problem 2. The second problem which plays a central role in MRI is the contrast problem. One must separate in a fixed transfer time t f the contrast between the two species. This amounts to introduce a cost function

c(q(t f )) = α|q 1 (t f )| 2 -β |q 2 (t f )| 2
where α, β are weight parameters. The contrast by saturation is to saturate one spin, e.g.spin 1, imposing q

1 (t f ) = 0 while maximizing |q 2 (t f )| 2 .
Both problems correspond to Mayer problems in optimal control, that is an optimal problem of the form min u(•) c(q(t f )) subject to q = F +uG, |u| ≤ 1, with some prescribed boundary conditions q(t f ) ∈ M, with M the terminal manifold defined by f (q) = 0. Again the maximum principle described by ( 6), [START_REF] Bonnard | Geometric and numerical methods in the contrast imaging problem in nuclear magnetic resonance[END_REF] applies and is complemented by transversality conditions depending upon the cost function and boundary conditions and written as:

f (q(t f )) = 0, p(t f ) = p 0 ∂ c ∂ q (q(t f )) + δ ∂ f ∂ q (q(t f )), p 0 ≤ 0, δ : constant vector (8)
Clearly, the optimal control problems boil down to analyze the so-called singular extremals given by

ż = H F (z) + uH G (z), H G = 0 ( 9
)
where u has to be admissible, that is |u| ≤ 1. o The surface Σ defined by H G (z) = 0 is called the switching surface.

Computation of the singular extremals can be made explicit. If F, G are two vector fields, the Lie bracket is defined with the convention:

[F, G](q) = ∂ F ∂ q (q)G(q) - ∂ G ∂ q F(q)
and the Poisson bracket of H F , H G is given by

{H F , H G }(z) = dH F ( H G (z)) = p • [F, G](q). Hence differentiating t → H G (z(t)) = 0, one gets z ∈ Σ : H F (z) = {H F , H G }(z) = 0 (10)
and the singular control is defined by solving

{{H G , H F }, H F }(z) + u{{H G , H F }, H G } = 0. ( 11
)
If F, G are real analytic, this defines a meromorphic equation of the form

dz dt = X(z) g(z) , z ∈ Σ . (12)
Since F and G are affine vector fields and depend linearly on the parameters, this leads to algebraic computations which are presented in this article and performed using computer algebra algorithms. More precisely, they are related to the following problem.

Find rational invariants with respect to the parameters ((Γ 1 , γ 1 ), (Γ 2 , γ 2 )) related to the classification of solutions of [START_REF] Bruns | Determinantal rings[END_REF].

Clearly this algebraic study is the first step of a computational process to classify the optimal strategies and analyze their robustness with respect to the physical parameters [START_REF] Bonnard | Geometric and numerical methods in the contrast imaging problem in nuclear magnetic resonance[END_REF].

The organization of this article is the following. In Section 2, we present based on [START_REF] Bonnard | Feedback equivalence for nonlinear systems and the time optimal control problem[END_REF] the geometric frame relating the classification of singular extremals with the optimal control problems and the action of the feedback group on the system. They are used to generate rational invariants. The main contribution of this article is presented in Section 3 which contains the algebraic computation. The algebraic techniques are sketched : Gröbner basis, a cellular decomposition of semi-algebraic set called cylindrical algebraic decomposition, a fast decision procedure for real root solving polynomial systems based on the critical point method and real roots classification based on Thom's isotopy lemma. They are the computational tools necessary to the algebraic computation in the contrast and multisaturation problems. They are given in details for one specific problem: the classification of the set of singularities of the determinantal surface D = 0 4 . This shows the complexity of the study (Figs.8 to Figs.13).

GEOMETRIC FRAME

The aim of this section is twofold. First of all, recall the geometric framework relating the classification of singular extremals and the optimal problems. Second, a neat approach is proposed to generate all Lie brackets, in relation with semi-direct Lie algebras.

The classification problem of singular extremals

In this section, we recall the results of [START_REF] Bonnard | Feedback equivalence for nonlinear systems and the time optimal control problem[END_REF] which will justify all our algebraic computations. Definition 2.1. Let E and F be two R-vector spaces, and let G be a group acting linearly on E and F. An homomorphism

X : G → R 0 is called a character. Let X be a character. A semi-invariant of weight X is a map λ : E → R such that for all g ∈ G and all x ∈ E, λ (g.x) = X (g)λ (x). It is an invariant if X = 1. A map λ : E → F is a semi-covariant of weight X is for all g ∈ G and for all x ∈ E, λ (g.x) = X (g)g.λ (x). It is called a covariant if X = 1.
Notations and definitions. Let U ⊂ R n be an open subset and let V ω (U) be the set of real analytic vector fields on U identified to maps from U to U. We denote by C ω (U, R) the set of analytic maps from U to R. Let G d be the group of C ω diffeomorphisms of U. The coordinates of TU U × R n are denoted z = (q, p) and we endow T * U with its canonical symplectic structure defined by Ω = ∑ n i=1 dq i ∧ dp i . Let Z be a C ω vector field on U, the Hamiltonian lift of Z is defined by the Hamiltonian H Z (z) = p • Z(q). If ϕ ∈ G d : q = ϕ(Q), one can lift ϕ into a symplectic diffeomorphism ϕ, called a Mathieu transformation defined by:

ϕ : q = ϕ(Q), p = P ∂ ϕ -1 ∂ Q
writing p and P as row vectors. A system q = X + uY , X,Y ∈ V ω (U), is written as (X,Y ) and let S = {(X,Y )}. Take (X,Y ), (X ,Y ) in S, they are called feedback equivalent if there exists ϕ ∈ G d and a feedback u = α(q) + β (q)v, α, β ∈ C ω (U, R), β (q) non-zero for each q, such that

X = ϕ * X + ϕ * Y.α Y = ϕ * Y.β
where ϕ * Z is the image of Z defined by ϕ * Z = ∂ ϕ -1 ∂ q (Z • ϕ). This action defines a group structure on the set of triplets (ϕ, α, β ) and this group is the feedback group and denoted by G f . Singular extremals of order 2. Let (F, G) ∈ S, the singular extremals of order 2 are the solutions of

ż = H F (z) + u s H G (z), z ∈ Σ : H G = 0 such that {{H G , H F }, H G }(z)
is never vanishing and u s is given by [START_REF] Brown | Qepcad b: a program for computing with semi-algebraic sets using cads[END_REF] as

u s = - {{H G , H F }, H F }(z) {{H G , H F }, H G }(z) (13)
and they are contained in the surface Σ :

H G (z) = {H G , H F }(z) = 0. Let S : {{H G , H F }, H G } = 0. Plugging such u s in H(z, u) = H F (z) + uH G (z)
defines a true Hamiltonian denoted Ĥ(z). One has the following, see [START_REF] Bonnard | Feedback equivalence for nonlinear systems and the time optimal control problem[END_REF].

Proposition 2.2. The singular extremals of order 2 are the solutions of the Hamiltonian equation

q = ∂ Ĥ ∂ p , ṗ = - ∂ Ĥ ∂ q , (q, p) ∈ Σ \ S (14)
Moreover, an explicit representation is as follows:

Lemma 2.3. Assume that G and [G, F] are never colinear on U. Then (Σ \ S, Ω ) is a symplectic manifold, Ω denoting the restriction of the standard symplectic form and Ĥ on Σ \ S is a C ω Hamiltonian whose solutions are singular extremals of order 2.

The action on the feedback group G f on the set of singular extremals of order 2 is defined as follows: if (ϕ, α, β ) ∈ G f , then the feedback (α, β ) acts trivially and ϕ acts by the change of symplectic coordinates ϕ.

The key result in our analysis is the following [START_REF] Bonnard | Feedback equivalence for nonlinear systems and the time optimal control problem[END_REF].

Theorem 2.4. Let λ be the map (F, G) → ( Ĥ, Σ , S) (differential equation defined by [START_REF] Canny | The complexity of robot motion planning[END_REF]). Then for the respective action of the feedback group G F , λ is a covariant.

Moreover, in relation with our geometric study, it is worth extracting two more covariants, where G F acts on functions by ϕ action only. This is called a feedback invariant.

Proposition 2.5. The following sets are feedback invariants:

(1) The set C where F and G are colinear ;

(2) The set C where G and [F, G] are colinear ;

(3) On Σ , the set S defined by {{H G , H F }, H G } = 0.

Reparameterization. In order to study the meromorphic equation ( 14), one makes a change of parameterization:

ds = dt/{{H G , H F }, H G } (15)
and observes that the map:

(F, G) → {{H G , H F }, H G } is a semi-covariant if z ∈ Σ . Theorem 2.6. The map λ ; (F, G) → ( Ĥ{{H G , H F }, H G }, Σ ) is a semi-covariant.
Note that the associated vector field is C ω but not Hamiltonian. It amounts to relate the classification problem to the classification of smooth vector fields and is the basic tool to generate invariants. Exceptional and generic cases. The set of singular extremals of order 2 can be split according to [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF] into two subsets, recalling that the level sets Ĥ = c are invariants. Exceptional case. It corresponds to Ĥ = 0 and the corresponding singular extremals of order 2 are called exceptional. Their role in Mayer optimal problems is analyzed in [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF].

In our situation, in dimension 4, the constraint Ĥ = 0 leads to the additional relation H F = p • F = 0. Hence the singular control is defined by

H F = H G = {H G , H F } = 0 {{H G , H F }, H F } + u s {{H G , H F }, H G } = 0 Since p ∈ R 4 \ {0} this leads to the relation D (q) + u s D(q) = 0
where D, D are the associated determinants

D = det(F, G, [G, F], [[G, F], G]) D = det(F, G, [G, F], [[G, F], F])
and the corresponding control is given by a feedback u s = -D (q) D(q) and the vector field is denoted X e (q) = F(q)-D (q) D(q) G(q), which again can be analyzed using the reparameterization ds = dt/D(q(t)) and this gives the smooth vector field

X e r = DF -D G. (16)
In this case the ϕ-action is reduced to the standard action of diffeomorphisms ϕ on U acting on tensors, that is, on vector fields V by image ϕ * V and on mappings f : U → R by composition. Generic case. It corresponds to Ĥ = c = 0 and by homogeneity one can take c = 1.

Let

D = {{H G , H F }, H G }, D = {{H G , H F }, H F }, the singular control u S is defined by H G = {H G , H F } = 0 D + u S D = 0.

Lie brackets computations

Each spin system can be lifted on the semi-direct product S = GL(3, R) R 3 acting on the q-space using the action (A, a).q = Aq + a. The Lie bracket computation rule is [(A, a), (B, b)] = ([A, B], Ab -Ba). Our system is written (F, G) with Γ, -Γ, -γ) and a = (0, 0, γ) • G = Cq where C is the antisymmetric matrix C = E 32 -E 23 with E i j = (δ i j ). According to the Lie bracket computation rule on the semidirect Lie product, we can reduce the computations to matrix Lie brackets. Moreover, to make such computations we use the following standard results, see [START_REF] Helgason | Differential geometry, Lie groups, and symmetric spaces[END_REF] for the details.

• F = (A, a), A = diag(-
The Lie algebra gl(n, R) of n × n matrices decomposes into cI n ⊕ sl(n, R) where sl(n, R) is the Lie algebra of matrices with zero trace. This algebra can be written as the direct sum a ⊕ b where a is the Cartan subalgebra of diagonal matrices and b = i = j RE i j . Moreover we have the Cartan decomposition sl(n, R) = so(n) ⊕ h where h is the subspace of symmetric matrices. If A = diag(λ 1 , . . . , λ n ) ∈ a then [A, E i j ] = (λ i -λ j )E i j and λ i -λ j will form the nontrivial spectrum of adA (where adA is defined as adA(B) = [A, B]) with corresponding eigenvectors E i j . Application. Restricting to the q = (y, z) space we have A = diag(-Γ, -γ) and a = (0, γ) and A = cI 2 + diag(λ , -λ ) where c = -Γ+γ 2 which is zero if and only if γ = Γ = 0 and λ = δ /2 where δ = γ -Γ. If δ = 0, the nontrivial spectrum is (δ , -δ ).

In NMR, we have the following Lemma 2.7. The case γ = Γ is the case of water species.

Otherwise, an easy computation gives that the Lie algebra generated by ((A, a), (C, 0)) is gl(2, R) ⊕ R 2 . Moreover, all Lie brackets can be easily computed. They are listed next, up to order 4, needed in our computation. Length 1.

• F = (-Γy, γ(1 -z)) • G = (-z, y) Length 2. • [F, G] = (γ -δ z, -δ y) Length 3. • [[F, G], F] = (-γ(γ -2Γ) + δ 2 z, -δ 2 y) • [[F, G], G] = (-2δ y, -γ + 2δ z) Length 4. • [[[F, G], F], F] = -δ 3 z + γ γ 2 -3γΓ + 3Γ 2 , -δ 3 y • [[[F, G], F], G] = [[[F, G], G], F] = (0, γ(γ -2Γ)) • [[[F, G], G], G] = (-γ + 4δ z, 4δ y) Couples of spins.
• Contrast, multisaturation: pairs with (Γ 1 , γ 1 ), (Γ 2 , γ 2 ) • Multisaturation with B 1 inhomogeneity: in this case Γ 1 = Γ 2 = Γ, γ 1 = γ 2 and for the second spin the vector field G is rescaled into (1 -ε)G.

ALGEBRAIC COMPUTATIONS

Preliminaries

Four test cases. In the physical experiments different cases will be considered. For the contrast problem:

• In vitro four test cases [START_REF] Bonnard | Algebraic geometric classification of the singular flow in the contrast imaging problem in nuclear magnetic resonance[END_REF] with relaxation times in seconds:

-P 1 : water (T 1 = T 2 = 2.5 s) / cerebrospinal fluid (T 1 = 2 s, T 2 = 0.3 s) -P 2 : deoxygenated blood (T 1 = 1.35 s, T 2 = 0.05 s) / oxygenated blood (T 1 = 1.35 s, T 2 = 0.2 s) -P 3 : gray cerebral matter (T 1 = 0.92 s, T 2 = 0.1 s) / white cerebral matter (T 1 = 0.780 s, T 2 = 0.09 s) -P 4 : water (T 1 = T 2 = 2.5 s) / fat (T 1 = 0.2 s, T 2 = 0.1 s) • In vivo one test case [START_REF] Van Reeth | Optimal control design of preparation pulses for contrast optimization in {MRI}[END_REF]:

-Q 1 : brain (T 1 = 1.062 s, T 2 = 0.052 s) / parietal muscle (T 1 = 1.200 s, T 2 = 0.029 s) Colinear sets (feedback invariants). We denote by Π 1 , Π 2 the respective projections on the first and second spind: q → q 1 , q → q 2 . Each colinear set is a curve described by the projections corresponding to the colinear set of each spin and an additional relation.

• C 1 : F, G linearly dependent: the projections are the ovals defined by

Γ i y 2 i = γ i (1 -z i )z i , 0 ≤ z 1 ≤ 1
intersected with either of the sets

Γ 1 y 1 z 2 = Γ 2 y 2 z 1 γ 1 y 1 (z 2 -1) = γ 2 y 2 (z 1 -1)
for the contrast setting, or with either of the sets

Γ(1 -ε)y 1 z 2 = Γy 2 z 1 γy 1 (z 2 -1) = γ(1 -ε)y 2 (z 1 -1)
for the B 1 -inhomogeneity setting. • C 2 : G, [G, F] linearly dependent: for a single spin, the projections are the 2 lines defined by

y = 0 2δ z = γ (if δ = 0, that is if the matter is not water)
For two spins in the contrast setting, the projections are:

-the curve y 1 = y 2 = 0, γ 1 z 2 -γ 2 z 1 = (δ 1 -δ 2 )z 1 z 2 ; -the 2 lines y i = 0, z 1 = γ i /(δ 1 + δ i ), z 2 = γ i /(δ 2 + δ i ) (i = 1, 2); -the plane z 1 = γ 1 /2δ , z 2 = γ 2 /2δ if δ 1 = δ 2 = δ .
For two spins in the B 1 -inhomogeneity setting, the projections are:

-the line y 1 = y 2 = 0, z 1 = z 2 ; -the plane z 1 = z 2 = γ/2δ .

Algebraic techniques

In this section, we give a brief overview of algebraic tools and methods to solve the algebraic problems posed by the study. We refer the reader to various textbooks [START_REF] Cox | Using algebraic geometry[END_REF][START_REF]Solving polynomial equations[END_REF][START_REF] Greuel | A Singular introduction to commutative algebra[END_REF][START_REF] Wang | Elimination methods. Texts and Monographs in Symbolic Computation[END_REF] for more information about those tools and others.

Operations on polynomial ideals

Throughout the rest of this section, we let n be a positive integer, and we consider polynomials in

A = C[X 1 , . . . , X n ].
Our main object of study shall be the polynomial ideal generated by a finite family F of polynomials in A, i.e. the subset of A of algebraic combinations of elements in F. Their importance in algebraic geometry stems from the fact that this is the set of polynomials which vanish at all common solutions of the finite subset F of A under consideration. Such an ideal is denoted by F ; the set of common complex zeroes of F is denoted by V (F) or V ( F ), and called an algebraic set.

Our work will consist in successive transformations of finite sets of generators ideals of polynomial equations, in order to enable information extraction on the algebraic set associated to those ideals. These transformations are mostly of two forms: elimination of variables, and saturation of polynomials.

Eliminating the first k variables from the ideal I means computing a set of generators of

I k := I ∩ C[X k+1 , . . . , X n ].
Geometrically, the set of zeroes of the ideal I k is the smallest algebraic set containing the projection of the set of zeroes V (I) onto the affine subspace C n-k with coordinates X k+1 , . . . , X n .

Given an ideal I and a polynomial f , saturating I by f means computing a set of generators of the ideal

(I : f ∞ ) := {g ∈ A | ∃ m ∈ N, g f n ∈ I}.
Geometrically, the zeroes of this ideal is the smallest algebraic set containing V (I) \V ( f ). This can be computed by adding the polynomial u • f -1 (with a new variable u) to the generators of I, and eliminating u from that ideal.

This effectively reduces our two major computational tools to the polynomial elimination problem. Because this problem is of such importance in computational algebraic geometry, many tools have been developed for solving it, for example triangular sets and regular chains [START_REF] Aubry | On the theories of triangular sets[END_REF][START_REF] Chen | Algorithms for computing triangular decomposition of polynomial systems[END_REF], resultants [START_REF] Canny | The complexity of robot motion planning[END_REF][START_REF] Canny | An efficient algorithm for the sparse mixed resultant[END_REF][START_REF] Cattani | Introduction to residues and resultants[END_REF][START_REF] Macaulay | Some Formulae in Elimination[END_REF] or Gröbner bases [START_REF] Buchberger | A theoretical basis for the reduction of polynomials to canonical forms[END_REF][START_REF] Faugère | A new efficient algorithm for computing Gröbner bases (F 4 )[END_REF][START_REF] Faugère | A new efficient algorithm for computing Gröbner bases without reduction to zero (F 5 )[END_REF]. In the following computations, we shall mainly use Gröbner bases, because they will allow us to perform both the elimination steps and routine simplification steps.

In the next subsection we give more details about how Gröbner bases can be used in our setting. Before that, we conclude this subsection with a few algebraic and algorithmic techniques which we shall use in the following computations.

Given two systems of polynomial equations F and G, the set

V (F) ∩V (G) is the algebraic set V (F ∪ G).
The radical of an ideal I is the set

√ I := { f ∈ A | ∃ m ∈ N, f m ∈ I}.
and it has the same set of zeroes as I. Computing the radical of an ideal is a difficult problem, which we will not need to tackle here. However, we shall sometimes be interested in, given a set of generators of I, computing a set of generators of another ideal J ⊃ I such that √ I = √ J. We refer to algorithms in [START_REF] Krick | An algorithm for the computation of the radical of an ideal in the ring of polynomials[END_REF] for computing Gröbner bases of the radical of an ideal given by a finite set of its generators.

Gröbner bases

A Gröbner basis of an ideal I is a particular set of generators of I, with additional properties. We only detail the information needed to understand subsequent computations, and we refer the interested reader to [22, Chap. 2-3] for a comprehensive overview of Gröbner basis techniques and elimination theory. In some sense, Gröbner bases are to polynomial systems what triangular linear systems are to linear system solving. As triangular systems depend on an implicit order on the variables used by Gaussian elimination in the solving process, Gröbner bases depend on an ordering on the monomials (non-linearity implies that ordering the variables is not sufficient). Hence, Gröbner bases, as objects and algorithms for computing them depend on a monomial ordering, that is a total ordering on the monomials

X α 1 1 • • • X α n n
compatible with multiplication. In the sequel, we will use mostly 3 orderings:

• the degree reverse-lexicographic ordering (DRL), denoted

X 1 > X 2 > • • • > X n ; • the k'th elimination ordering (elim-k), denoted X 1 > • • • > X k X k+1 > • • • > X n
, and defined as a block-DRL ordering;

• the lexicographical ordering lex, denoted X 1 X 2 • • • X n .
Of course, the choice of monomial orderings impacts on the efficiency of the computations. Usually, the DRL ordering yields faster computations, but is nonetheless useful for simple cases, when one merely wants to simplify the generators of an ideal.

The elim-k ordering allows us to eliminate the first k variables: if G is a Gröbner basis of I with respect to the elim-k ordering, then

G ∩ C[X k+1 , . . . , X n ] is a (DRL) Gröbner basis of the ideal I ∩ C[X k+1 , . . . , X n ].
Finally, the lexicographical order has the elimination property for all indices k: from a single Gröbner basis, one can recover all the intermediate elimination bases. This is usually enough to obtain a full description of the solutions, but these bases are also the hardest to compute directly.

Most computer algebra systems include at least one toolbox for Gröbner basis computations. For example in Maple, one can use the built-in package Groebner as well or the third-party library FGb [START_REF] Faugère | FGb: A Library for Computing Gröbner Bases[END_REF].

Cylindrical algebraic decomposition

The last subsections are focused on tools for real algebraic geometry. The basic items of study of real algebraic geometry is semi-algebraic sets, that is a sets S ⊂ R n defined as a finite union S = r i=1 S i of sets defined by polynomial equations P i and inequalities Q i with real coefficients:

S i = x ∈ R n ∀ P ∈ P i , P(x) = 0 ∀ Q ∈ Q i , Q(x) > 0
and maps between those sets. Given two semi-algebraic sets A and B, a map A → B is called a semi-algebraic function if its graph is a semi-algebraic set. For example, polynomial or rational functions are semi-algebraic.

A Cylindrical Algebraic Decomposition (CAD) [START_REF] Chen | Quantifier Elimination by Cylindrical Algebraic Decomposition Based on Regular Chains[END_REF][START_REF] Collins | Quantifier elimination for real closed fields by cylindrical algebraic decomposition[END_REF][START_REF] Mccallum | On using lazard's projection in CAD construction[END_REF] of R n is a sequence S 1 , . . . , S n such that for all i ∈ {1, . . . , n}, S i is a decomposition of R i into connected semi-algebraic subsets, called cells, and such that:

• Cells of S 1 are all points or intervals;

• For any i ∈ {1, . . . , n}, and for any cell S of S i , there exists a finite number of continuous semi-algebraic functions ξ S,1 , . . . , ξ S,l S : S → R such that S × R ⊂ R i+1 is a disjoint union of cells of S i+1 , which are either the graph Γ S, j of a function ξ S, j or a band B S, j bounded by ξ S, j and ξ S, j+1 for j ∈ {0, . . . , l S }, with the convention that ξ S,0 = -∞ and ξ S,l S +1 = +∞:

B S, j = {(x , x i+1 ) ∈ S × R | ξ S, j (x ) < x i+1 < ξ S, j+1 (x )}.
Given a set of polynomials P ⊂ R[X 1 , . . . , X n ], a CAD S = (S 1 , . . . , S n ) is P-invariant if, for any cell S ∈ S n and any polynomial P ∈ P, P has constant sign on S.

In particular, given a semi-algebraic set defined by polynomial equations P and polynomial inequalities Q, a P ∪ Qinvariant CAD gives an extensive topological description of the semi-algebraic set by means of a partition into subsets homeomorphic to ]0, 1[ i for some i ∈ {0, . . . , n}.

Unfortunately, computing a P-invariant cylindrical algebraic decomposition, for P, is intrinsically difficult: the complexity of the algorithm and the size of the output, for most entries, is doubly exponential in n.

Nonetheless, this tool will prove useful for completing computations when n ≤ 2 or 3 needed for e.g. real root classification algorithms described hereafter. In this context, the main purpose of the CAD will be to compute one point per connected component of an open semi-algebraic set, defined by inequations P i = 0 and strict inequalities P j > 0. In the next paragraph, we describe an efficient alternative to cylindrical algebraic decision which we used for solving systems which were unreachable to Cylindrical Algebraic Decomposition because the number of variables was greater than 2.

Algorithms for computing Cylindrical Algebraic Decompositions are implemented in the C library QEPCAD [START_REF] Brown | Qepcad b: a program for computing with semi-algebraic sets using cads[END_REF], the command CylindricalDecomposition in the computer algebra system Mathematica and the subpackage SemialgebraicsetTools of RegularChains in the computer algebra system Maple.

Real root solving for polynomial systems

As we already said, the complexity (in terms of output size and run time) of Cylindrical Algebraic Decomposition is doubly exponential in the number n of variables. Recall also that a Cylindrical Algebraic Decomposition provides a partition of the semi-algebraic set under consideration into semi-algebraic pieces which are homeomorphic to ]0, 1[ i for 0 ≤ i ≤ n.

As we will see in the next paragraph, such a partition is not needed to solve classification problems over semi-algebraic sets. We only need to compute sample points in each connected component of a given semi-algebraic set. Thanks to the so-called critical point method initiated in [START_REF] Yu | Solving systems of polynomial inequalities in subexponential time[END_REF] and complexity improvements [2, Chapter 13], it has been proved that this can be done in time (sD) O(n) where s is the number of polynomials defining the semi-algebraic set under consideration and D is the maximum degree of those polynomials. The core geometric idea, on which this method is based, is to reduce the problem to a polynomial optimization one.

It took a decade to obtain algorithms based on this method which have a practical impact [START_REF] El Din | Polar varieties and computation of one point in each connected component of a smooth real algebraic set[END_REF][START_REF] Greuet | Probabilistic algorithm for polynomial optimization over a real algebraic set[END_REF][START_REF] Din | Testing sign conditions on a multivariate polynomial and applications[END_REF]. These are implemented and available in the RAGlib (Real Algebraic Geometry library) Maple package [START_REF] Din | Raglib (real algebraic geometry library[END_REF]. In the sequel, we roughly describe these algorithms.

The input consists of a polynomial system of equations and inequalities

f 1 = • • • = f p = 0, h 1 > 0, . . . , h r > 0.
defining a semi-algebraic set S ⊂ R n . We aim at computing sample points in each connected component of S.

The first step is a reduction to computing sample points in each connected component of a semi-algebraic set defined with non-strict inequalities. Indeed, Observe that there exists e 0 > 0 such that for any connected component C of S and all 0 < e < e 0 , there exists a connected component C e of the semi-algebraic set S e defined by

f 1 = • • • = f p = 0, h 1 ≥ e, . . . , h r ≥ e.
Finding a suitable e 0 can be done using notions of critical values and asymptotic critical values (see e.g. [START_REF] Greuet | Probabilistic algorithm for polynomial optimization over a real algebraic set[END_REF] and references therein).

This reduction to the study of closed semi-algebraic sets opens some perspectives. Indeed, the next step is to reduce this semi-algebraic problem to algebraic ones. The cornerstone of this reduction comes with the result below. Proposition 3.1. [2, Proposition 13.2] Let S ⊂ R n be a semi-algebraic set defined and e 0 be as above C ⊂ S be a connected component of S. Then, for 0 < e < e 0 , there exists {i 1 , . . . , i } ⊂ {1, . . . , r} and a connected component C of the real algebraic set defined by

f 1 = • • • = f p = 0, h i 1 = • • • = h i = e such that C ⊂ C.
Observe that using the above result we can reduce the problem of computing sample points in each connected component of S to the problem of computing sample points in several real algebraic sets and filter those points at which the inequalities used to define S are positive.

Let us mention that using standard computer algebra techniques, the number of real algebraic sets to be considered can be well-controlled and lies in s O(n) with s = p + r.

It remains to explain how we compute sample points in real algebraic sets. Hence, let us assume that a real algebraic set is given as the real solution set to the system

f 1 = • • • = f k = 0.
In order to keep this description as simple as it can be, we assume that any solution to this system, its jacobian matrix has full rank ; hence the real algebraic set is smooth.

The core idea developed in [START_REF] El Din | Polar varieties and computation of one point in each connected component of a smooth real algebraic set[END_REF] is the following. Let a = (a 1 , . . . , a n ) ⊂ Q n be generically chosen and π a :

(x 1 , . . . , x n ) → a 1 x 1 + • • • + a n x n .
Then for any connected component C of the real algebraic set under study, either C contains a critical point of the map π a restricted to the set under study (then such points can be computed using e.g. Gröbner bases, see [START_REF] Faugère | Critical points and gröbner bases: the unmixed case[END_REF]), or for an arbitrary ρ ∈ R, π -1 a (ρ) meets C (hence it is sufficient to perform a recursive call to the algorithm by adding the constraint a 1 X 1 + • • • + a n X n to the input). The overall cost of those computations lie in D O(n) .

Implementations of such algorithms are available in the RAGlib Maple package based on the C library FGb for computing Gröbner bases.

Real roots classification

Our last algebraic tool is a computational technique for solving the problem of classifying the real fibers of projection restricted to an algebraic set.

Consider an algebraic set V in C n+t , with coordinates (X, G) = (X 1 , . . . , X n , G 1 , . . . , G t ) and the projection onto the parameter space π : C n+t → C t (x, g) → g Assume that V has dimension t. Let B ∈ R n+t be a semi-algebraic set defined only by a set of polynomial inequalities

H ⊂ R[X, G]: B = (x, g) ∈ R n+t | ∀ h ∈ H, h(x, g) > 0 , we consider the semi-algebraic set V ∩ B ⊂ R n+t .
The classification problem is to identify a dense covering of the real parameter space R t into open subsets U 1 , . . . ,U N such that on each U i , the number of points in the fibers of π restricted to V ∩ B is constant:

∃ c i , ∀ g ∈ U i , #π -1 (g) ∩V ∩ B = c i .
This problem has been thoroughly studied, see [START_REF] Lazard | Solving parametric polynomial systems[END_REF][START_REF] Yang | A complete algorithm for automated discovering of a class of inequality-type theorems[END_REF] among other references. In our situation, the variety V shall be defined as the zero set of D and its partial derivatives along the X i 's, where D is the determinant of a matrix. The strategy that we use, as described in [START_REF] Bonnard | Determinantal sets, singularities and application to optimal control in medical imagery[END_REF], is adapted from the more general strategy of described in [START_REF] Lazard | Solving parametric polynomial systems[END_REF][START_REF] Yang | A complete algorithm for automated discovering of a class of inequality-type theorems[END_REF] but dedicated to the structure of the polynomial systems which appear in our context. This enables our variant algorithm to tackle those classification issues which were out of reach to the implementations of the general strategy.

We summarize first the key points of this general strategy. One wishes to compute a set of polynomials {p 1 , . . . , p M } ⊂ R[G] such that the dense covering (U i ) given by the connected components of R t \ (V (p 1 ) ∪V (p 2 ) ∪ • • • ∪V (p M )). These polynomials correspond to equations satisfied by points at which the cardinality of the fibers change. These points may be:

• the projection of singular points of V ;

• the critical values of π restricted to V ;

• the projection of points where V meets the topological border ∂ B of the semi-algebraic set B;

• points g at which the fiber π -1 (g) ∩V has positive real dimension;

• points at which the projection π restricted to V ∩ B is not proper. This is formalized in the following lemma, whose proof relies on Thom's isotopy lemma. We first set some notations and hypotheses. If B is a semi-algebraic set defined by polynomial inequalities H, we define

B 0 := h∈H V (h).
Note in particular that ∂ B ⊂ B 0 . We denote by sing(V ) the singular locus of V , crit(π,V ) the critical locus of π restricted to V , and

K(π,V ) = sing(V ) ∪ crit(π,V ).
The hypotheses are:

H1 There exists a nonempty Zariski-open subset O 1 of C t such that for all g ∈ O 1 , the fiber V ∩ π -1 (g) is a nonempty finite subset of C n+t H2 The restriction of the projection π to B is proper; H3 The intersection V ∩ B 0 has dimension at most t -1 in C nt ; H4 The algebraic set V is equidimensional with dimension t. Lemma 3.2. Let V and B satisfying hypotheses H1, H2, H3 and H4. Let U be a nonempty connected open subset of R t which does not meet π(V ∩ B 0 ) ∪ π(K(π,V )), and let g ∈ U. Then V ∩ π -1 (g) is finite, and for any g ∈ U,

#V ∩ π -1 (g ) = #V ∩ π -1 (g).
Hence, under these hypotheses, it suffices to compute a nonzero polynomial whose zeroes cover π(V ∩ B 0 ) and π(K(π,V )) to obtain the dense covering that we need.

The strategy that we use refines these computations by taking advantage of the fact that our variety V comes from a determinantal variety, defined as follows. Let M be a k × k matrix with coefficients in R[X, G], r ∈ {0, . . . , k -1}, the r'th determinantal variety associated with M is

V r (M) = (x, g) ∈ C n+t | rank(M(x, g)) ≤ r .
These varieties can be defined algebraically in two ways. First, one can use the minors modelization:

V r (M) = V ((r + 1)-minors of M).
Another way to define V r is as the projection of the incidence variety defined as the set of solutions of

M(x, g) •Y (y) = 0
where Y (y) is a full-rank k × kr matrix whose entries are new variables y = y 1,1 , . . . , y k,k-r .

We now describe more precisely how this property allows us to compute the classification for the MRI problem. The variety V that we will consider is defined as K(π,V r 0 (M)) for r 0 ∈ {0, . . . , k -1}. Assume that the variety V satisfies the following three properties:

H5 There exists a nonempty Zariski-open subset O 2 ⊂ C t such that V ∩ π -1 (O 2 ) = V r 0 -1 ∩ π -1 (O 2 );
H6 For any r 1 < r 2 ∈ {0, . . . , k -1}, for any matrix A of size r 1 in M, the ideal generated by the r 2 -minors of M containing A is radical; H7 For any r ∈ {0, . . . , k -1}, the variety V r is equidimensional with dimension n + t -(k -r + 1) 2 . These properties are generic [START_REF] Bruns | Determinantal rings[END_REF], [START_REF] Ranestad | Algebraic Degree in Semidefinite and Polynomial Optimization[END_REF]. If we assume that n = (k -r 0 + 1) 2 , then hypotheses H7 and H5 imply hypotheses H4 and H1 for V .

Under these hypotheses, the strategy consists in computing:

(1) the projection of V \V r 0 -1 , which has dimension less than t;

(2) the projection of K(π,V r 0 -1 );

(3) the projection of B 0 ∩V r 0 -1 .

Remark 3.3. In section 3.3.3, we will also be describing K(π,V ) where V is an incidence variety defined by the equations

P t • M(x, g) = 0
where P is a vector of dimension 4 and M is a matrix of size 4 × 3. This incidence variety encodes an underdefined linear system, so its projection onto the affine space C n+t with coordinates (X, G) is the affine space itself. In this situation, generically, sing(V ) is exactly V 2 (M).

3.3.

Algebraic computations for multisaturation with B 1 -inhomogeneity

Conventions and notations

In this section and the next one, we use the change of coordinates

z 1 ← z 1 + 1, z 2 ← z 2 + 1,
which places the center of the coordinates at the North pole of the Bloch ball. In this new system of coordinates, the center of the Bloch ball has coordinates (0, -1, 0, -1). With these coordinates, the determinants D and D defined in Sec. 2.1 are written

D = det(F, G, [G, F], [[G, F], G]) = det     -Γy -z -1 δ z -Γ 2δ y -γz y δ y -2δ z + Γ -δ -Γy (1 -ε)(-z -1) (1 -ε)(δ z -Γ) (1 -ε) 2 (2δ y) -γz (1 -ε)y (1 -ε)δ y (1 -ε) 2 (-2δ z + Γ -δ )     and D = det(F, G, [G, F], [[G, F], F]) =     -Γy -z -1 δ z -Γ γ(γ -2Γ) + δ 2 (z + 1) -γz y δ y δ 2 y -Γy (1 -ε)(-z -1) (1 -ε)(δ z -Γ) (1 -ε) 2 (γ(γ -2Γ) + δ 2 (z + 1)) -γz (1 -ε)y (1 -ε)δ y (1 -ε) 2 δ 2 y    
The following polynomials will appear frequently in the remainder of the section:

• P y 1 := y 1 -(1 -ε)y 2 • P y 2 := y 2 -(1 -ε)y 1 • P z 1 := 2(Γ -γ)z 1 + 2Γ -γ • P z 2 := 2(Γ -γ)z 2 + 2Γ -γ
The root of the univariate polynomials P z 1 and P z 2 is

z S = γ -2Γ 2Γ -2γ .
= ∂ D ∂ z 1 = ∂ D ∂ y 2 = ∂ D ∂ z 2
= 0 is given, generically on authorized values of γ, Γ, by

(1) the point y 1 = y 2 = z 1 = z 2 = 0, and (2) the curve defined by P y 1 = P z 1 = P z 2 = 0, which is parameterized by y 2 as

y 1 = (1 -ε)y 2 z 1 = z 2 = z S = γ-2Γ 2Γ-2γ . If γ = Γ (for
example if the matter is water), only the former solution exists.

Proof. The structure of this proof is summarized in Fig. 6. The determinant D can be factored as (1 -ε) D. The singularities of D and those of D are the same, so for the study, we consider the ideal

I := D, ∂ D ∂ y 1 , ∂ D ∂ y 2 , ∂ D ∂ z 1 , ∂ D ∂ z 2 .
In order to eliminate y 1 , y 2 and z 1 from the ideal I, we compute a Gröbner basis G of I with respect to the elimination ordering y 1 > y 2 > z 1 z 2 > ε > Γ > γ.This computation yields that

I ∩ Q[z 2 , ε, Γ, γ] = ε 2 (ε -2) 2 (2 Γ -γ)(Γ -γ)z 3 2 P 3 z 2 , D, ∂ D ∂ y 1 , ∂ D ∂ y 2 , ∂ D ∂ z 1 , ∂ D ∂ z 2 D, ∂ D ∂ y 1 , ∂ D ∂ y 2 , ∂ D ∂ z 1 , ∂ D ∂ z 2 z 2 = 0 z 1 = 0 y 2 = 0 y 1 = 0 P z 2 = 0 z 2 = 0 z 1 = z 2 y 1 = 0 P y 1 = 0 FIGURE 6.
Structure of the study of the singularities of {D = 0} so singular points necessarily satisfy

     z 2 = 0 or P z 2 = 0 ⇐⇒ z 2 = γ-2 Γ 2 Γ-2 γ . If Γ = γ (
that is, if the matter is water), the second of these solutions does not exist. If γ = 2 Γ (which means that the matter is on the limit of the domain of validity 2 Γ ≥ γ), both solutions coincide.

In all other cases, there are 2 distinct possible values for z 2 , and we consider both cases: we consider the two ideals

I 1 := sqfr(G), z 2 I 2 := sqfr(G), P z 2
where for any polynomial f , sqfr( f ) is the square-free part of f and sqfr(G) means that we apply sqfr to each element of G.

In order to lift the partial solution z 2 = 0, we compute a Gröbner basis G 1 of I 1 with respect to the ordering y 1 > y 2 z 1 > z 2 > ε > Γ > γ, and we find that this ideal contains

γz 2 1 (ε -1) 2 (2Γ -γ), so z 1 = 0.
We then compute a Gröbner basis of sqfr(G 1 ), z 1 with respect to the order y 1 y 2 > z 1 > z 2 > ε > Γ > γ, and we find that this ideal contains

Γγεy 2 2 (Γ -γ)(2Γ -γ) 2 (ε -2), so y 2 = 1.
Finally, adding y 2 to the ideal yields that 0 = Γy 1 (ε -1)(2Γ -γ), so the final solution is (y 1 , y 2 , z 1 , z 2 ) = (0, 0, 0, 0).

We now consider the partial solution z 2 = (γ -2Γ)/(2Γ -2γ). We compute a Gröbner basis G 2 of I 2 with respect to the order y 1 > y 2 z 1 > z 2 > ε > Γ > γ, and we find that the ideal contains

z 2 γ(z 1 -z 2 ) 2 .
Since this case was already studied, we may assume that z 2 = 0, so

z 1 = z 2 = γ -2Γ 2Γ -2γ .
Adding z 1 -z 2 to sqfr(G 2 ) and computing a Gröbner basis for the order y 1 y 2 > z 1 > z 2 > ε > Γ > γ, we find that the ideal contains γ 2 y 2 (ε -1)P y 1 (2Γ -γ), so we have 2 new branches to consider.

If

y 2 = 0, y 1 = (1 -ε)y 2 .
Otherwise, by adding y 2 = 0 to the system of equations, we find that the ideal contains

γy 2 1 z 2 ,
so y 1 = 0, and in particular, this point is on {P y 1 = 0}.

Locus of {D = D = 0}.

Proposition 3.5. The points of {D = D = 0} are given by:

(1) the plane

z 1 = z 2 = z S = γ -2Γ 2Γ -2γ (17)
(2) the line

y 1 = y 2 = 0 z 1 = z 2 (18)
(3) the surface (parameterized by y 1 , y 2 )

z 1 = z 2 = ΓP 2 y 2 (γ -2Γ) 2(Γ -γ)a 3 ( 19
)
with

a 3 = (Γ + γ)P 2 y 1 + ε(ε -2)Γ(y 1 -y 2 )(y 1 + y 2 )
(4) the surface (parameterized by y 1 , z 2 )

y 2 = y 1 z 2 (1 -ε)z 1 , z 1 = (2Γ -γ)z 2 a 4 ( 20
)
with solution 4)

a 4 = 2(ε -2)(Γ -γ)εz 2 + (2Γ -γ)(ε -1) 2 D = D = 0 D = 0 z 1 -z 2 = 0 P z 2 = 0 (solution 1) p 2 = 0 y 1 = y 2 = 0 (solution 2) a 3 = 0 (solution 3) D = 0 (ε -1)y 2 z 1 + y 1 z 2 = 0 y 1 = z 1 = 0 y 2 = z 2 = 0 (⊂ solution 2)
y 1 = 0 z 1 = 0 p 3 = 0 (
y 1 = 0 p 4 = 0 p 5 = 0 FIGURE 7.
Structure of the study of {D = D = 0}

(5) the surface (parameterized by y 2 , z 2 )

z 1 = z 2 y 1 (1 -ε)y 2 , y 1 = (1 -ε)y 2 (2Γ -γ)Γy 2 2 + γ 2 z 2 2 a 5 (21)
with

a 5 = Γ 2ε(ε -2)(Γ -γ)z 2 + (ε -1) 2 (2Γ -γ) y 2 2 + γ 2 z 2 2
Proof. The structure of this proof is summarized in Fig. 7. The determinant D factors as

D = 2γ 2 (2Γ -γ)(Γ -γ)(z 1 -z 2 )(ε -1)((ε -1)y 2 z 1 + y 1 z 2 ),
so we form the two ideals

I 1 = D, z 1 -z 2 I 2 = D, (ε -1)y 2 z 1 + y 1 z 2 .
If z 1 = z 2 , after substitution, D has two factors depending on y 1 , y 2 , z 2 : P z 2 and

p 2 := 2(Γ -γ) (Γ + γ)P 2 y 1 + ε(ε -2)Γ(y 1 -y 2 )(y 1 + y 2 ) z 2 + ΓP 2 y 2 (2Γ -γ)
The polynomial P z 2 gives solution 1. Let a 3 (y 1 , y 2 ) = (Γ + γ)P 2

y 1 + ε(ε -2)Γ(y 1 -y 2 )(y 1 + y 2 )
so that the coefficient of z 2 in p 2 is 2(Γγ)a 3 , it is homogeneous in y 1 , y 2 with degree 2. Its discriminant in y 2 is -4(ε -2) 2 ε 2 y 2 1 γΓ.

Since the parameters γ, Γ are necessarily positive, this discriminant is negative, and thus the only real root of a 3 (y 1 , y 2 ) is y 1 = y 2 = 0. If y 1 = y 2 = 0, p 2 vanishes regardless of z 2 .

If y 1 = 0, a 3 (y 1 , y 2 ) does not have any real root in y 2 , and z 2 is given by

(z 1 =)z 2 = ΓP 2 y 2 (γ -2Γ) 2(Γ -γ)a 3 (y 1 , y 2 )
.

We now turn to the other branch, defined by (ε -1)y 2 z 1 + y 1 z 2 = 0.

If y 1 = z 1 = 0, there are 2 curves of singular points defined (in y 2 , z 2 ) by

Γ(2Γ -γ)y 2 2 + γ 2 z 2 2 = 0.
Since 2Γ ≥ γ, the only solution is z 2 = 0 with either 2Γ = γ or y 2 = 0.

If y 1 = 0 and z 1 = 0, then (since ε = 1) we must have y 2 = 0. Furthermore, we may assume that z 1 = z 2 since this case was already studied. The remaining solutions form a curve defined by

0 = p 3 := 2(ε -2)(Γ -γ)εz 2 + (2Γ -γ)(ε -1) 2 z 1 -(2Γ -γ)z 2 .
Let a 4 (z 1 , z 2 ) be the coefficient of z 1 in p 3 , the solutions are given by either

     a 4 (z 2 ) = 0 z 1 = 2Γ -γ c 3 (z 2 ) z 2
or (since by assumption

z 1 = 0) z 2 = 2Γ -γ = 0.
So we may assume that y 1 = 0. We compute a Gröbner basis of I 2 + uy 1 -1 for the order u z 1 > y 1 z 2 > y 2 > ε > γ > Γ. This basis contains a polynomial which factors as the product of

p 4 = 2(ε -2)(Γ -γ)εz 2 + (ε -1) 2 (2Γ -γ) y 1 + (ε -1)(2Γ -γ)y 2 and p 5 = Γ 2ε(ε -2)(Γ -γ)z 2 + (ε -1) 2 (2Γ -γ) y 2 2 + γ 2 z 2 2 y 1 + (ε -1)y 2 (2Γ -γ)Γy 2 2 + γ 2 z 2 2 .
First, assume that p 4 = 0. We compute a Gröbner basis of I 2 + uy 1 -1, p 4 for the order u z 1 > y 1 z 2 > y 2 > ε > γ > Γ, and we find that the last polynomial defining the ideal is p 3 , whose solutions we already studied.

Finally, assume that p 4 = 0 and p 5 = 0. The discriminant in y 2 of the coefficient a 5 (y 2 , z 2 ) of y 1 in p 5 is

-4a 4 (z 2 )Γγ 2 z 2 2
where a 4 (z 2 ) is as above the coefficient of z 1 in p 3 . The last components of the solutions are given by Proof. Assume that at some point, either of the determinants D and D is non-zero, this implies that F and G are colinear. Since F and G form the first two columns of the matrices whose D and D are the respective determinants, D = D = 0 at that point. Linearization of the system at equilibrium points. For each of the components of the set of equilibrium points {D = D = 0} found in the previous paragraph, we inspect the behavior of the system in a neighborhood. Namely, for each equilibrium point q, we write d dt

   a 5 (y 2 , z 2 ) = 0 y 1 = (1-ε)y 2( (2Γ-γ)Γy 2 2 +γ 2 z 2 2 ) Γ(2ε(ε-2)(Γ-γ)z 2 +(ε-1) 2 (2Γ-γ))y 2 2 +γ 2 z 2
(q + δ q) = (DF -D G)(q) + A(q) • δ q + R(δ q).
where

A = Jac q (DF -D G), so that d dt (δ q) = A(q) • δ q + R(q)(δ q).
We can compute A(q) explicitely: Indeed, let

f = DF -D G. Its first derivative is d f (q)(u) = dD(q)(u)F(q) + D(q)dF(q)(u) -dD (q)(u)G(q) -D (q)dG(q)(u), (22)
so A(q) = ∇D(q).F(q) + D(q)Jac q (F)(q) -∇D (q).G(q) -D (q)Jac q (G)(q).

We examine the eigenvalue decomposition of A(q). Solution 1 [START_REF] Chen | Algorithms for computing triangular decomposition of polynomial systems[END_REF]. If

z 1 = z 2 = γ-2Γ
2Γ-2γ , the characteristic polynomial of A factors as

T 2 T -γ 2 (2Γ -γ) 2 (ε -1)P 2 y 1 2
The matrix A(q) is diagonalizable. Solution 2 [START_REF] Chen | Quantifier Elimination by Cylindrical Algebraic Decomposition Based on Regular Chains[END_REF]. If y 1 = y 2 = 0 and z 1 = z 2 , the characteristic polynomial of A(q) is T 4

The Jacobian matrix A(q) can be trigonalized as

A(q) = P -1     0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0     P
with the transition matrix

P =     0 1 -1 0 εγ 3 (ε -1)(ε -2)z 2 1 P z 1 1 0 0 0 1 0 1 εγ 3 (ε -1)(ε -2)z 2 1 P z 1 0 0 0     .
Solution 3 [START_REF] Collins | Quantifier elimination for real closed fields by cylindrical algebraic decomposition[END_REF]. If z 1 = z 2 = ΓP 2 y 2 (γ -2Γ)/2(Γγ)a 3 , the characteristic polynomial of A(q) factors as

T 2 T + b 3 a 3 T - b 3 a 3 with b 3 = Γγ 2 (ε -1)P y 1 P y 2 (2Γ -γ) 2 .
The matrix A(q) is diagonalizable. Solution 4 [START_REF] Conolly | Optimal control solutions to the magnetic resonance selective excitation problem[END_REF]. If

y 2 = y 1 z 2 (1-ε)z 1 and z 1 = (2Γ-γ)z 2 a 4
, the characteristic polynomial of A(q) factors as

T 2 T - b 4 a 4 T + b 4 a 4 with b 4 = 2ε 2 γ 3 z 3 2 (ε -1)(ε -2) 2 (2Γ -γ)(Γ -γ)P z 2 .
The matrix A(q) is diagonalizable.

Solution 5 [START_REF] Cots | Contrôle optimal géométrique : méthodes homotopiques et applications[END_REF]. If

z 1 = z 2 y 1 (1-ε)y 2 and y 1 = (1-ε)y 2( (2Γ-γ)Γy 2 2 +γ 2 z 2 2 ) a 5
, the characteristic polynomial of A(q) factors as

T 2 T - b 4 (Γy 2 2 + γ(z 2 2 + z 2 ))Γy 2 2 a 5 2 .
The matrix A(q) is diagonalizable. Special points. There are two points at which A vanishes: the North pole N = (0, 0, 0, 0) and S = (0, z S , 0, z S ). Both points are such that D = D = 0, ∇D = ∇D = 0, and additionally, at the North pole, F(N) = 0. The North pole is on solutions 2, 3, 4 and 5. The remainder at N is cubic:

d dt (N + δ q) = R(N)(δ q) = O( δ q 3 ).
The point S is the intersection of solutions 1 and 2. The remainder at S is quadratic. Higher order studies for the special points. Quadratic approximation at S. We now study the quadratic component H 2 = Q(S) of the remainder R(S):

d dt (q + δ q) = (DF -D G)(q) + A(q)(δ q) + Q(q)(δ q) + O( δ q 3 ),
with dq dt (S) = (DF -D G)(S) = 0 and A(S) = 0. We can compute Q by differentiating f = DF -D G again, as was done in [START_REF] Bonnard | Algebraic geometric classification of the singular flow in the contrast imaging problem in nuclear magnetic resonance[END_REF]Sec. 3.4]. Differentiating [START_REF] Cox | Using algebraic geometry[END_REF] along q again, the second derivative of f is

d 2 f (q)(u, v) = d 2 D(q)(u, v)F(q) + dD(q)(u)dF(q)(v) + dD(q)(v)dF(q)(u) -d 2 D (q)(u, v)G(q) -dD (q)(u)dG(q)(v) -dD (q)(v)dG(q)(u) (23)
Note that second derivatives of F and G are 0, since their coordinates are affine in q.

We wish to compute H 2 (δ q) = Q(S)(δ q, δ q) = 1 2 d 2 f (q)(δ q, δ q). Since dD(S) = dD (S) = 0, we find in the end that

H 2 (δ q) = h 2 (δ q)F(S) -h 2 (δ q)G(S), with
F(S) = 0, γ(2Γ -γ) 2(Γ -γ) , 0, γ(2Γ -γ) 2(Γ -γ) t G(S) = γ 2(Γ -γ) , 0, (1 -ε)γ 2(Γ -γ) , 0 t h 2 (δ q) = 1 2 d 2 D(S)(δ q, δ q) = (1 -ε)(δ z 1 -δ z 2 )(δ z 1 -(1 -ε) 2 δ z 2 )(2Γ -γ)γ 2 h 2 (δ q) = 1 2 d 2 D (S)(δ q, δ q) = (1 -ε)(δ z 1 -δ z 2 )(δ y 2 (ε -1) + δ y 1 )(2Γ -γ) 2 γ 2 .
Following [START_REF] Bonnard | Algebraic geometric classification of the singular flow in the contrast imaging problem in nuclear magnetic resonance[END_REF] and [START_REF] Markus | Quadratic differential equations and non-associative algebras[END_REF], we study the projection of the differential equation v = H 2 (v) on the sphere S 3 . Let w = v/ v be this projection, it satisfies the differential equation

ẇ = 1 v 2 v v -v v, v v = H 2 (v) v - v, H 2 (v) v 3 v = v H 2 (w) -w, H 2 (w) w
so we are to study the following differential equation on the sphere

S 3 : v = H 2 (v) -v, H 2 (v) v =: H π 2 (v).
Invariants are related to the eigenvalues of the linearization of H π 2 at points where H π 2 (v) = 0. Those points are: • lines of non-isolated singular points of H 2 , that is vectors v such that H 2 (v) = 0 • ray solutions, that is vectors ξ such that there exists λ ∈ R \ {0}, H 2 (ξ ) = λ ξ . We study the linearization of H π 2 in some neighborhood of these solutions in S 3 .

Proposition 3.7. The blow-up at point S has no ray solution, and two sets of non-isolated singularities:

(1) the projective plane δ z 1 = δ z 2 ;

(2) the projective line δ

y 2 = (1 -ε)δ y 1 , δ z 1 = (1 -ε) 2 δ z 2 .
In the first case, the Jacobian of the system is nilpotent. In the second case, it is diagonalizable with non-zero eigenvalues:

1 2 δ y 2 + 1 ± δ y 2 2 + (2ε -1) 2 ( δ y 2 + 1) -4(ε -1) 4 -2 δ y 2 + 1 .
Proof. First we study ray solutions. Let ξ be a vector on a ray solution, such that

H 2 (ξ ) = λ ξ .
Let αξ be another vector on the same line (α ∈ R), since H 2 is homogeneous with degree 2, one has

H 2 (αξ ) = α 2 H 2 (ξ ) = α 2 λ ξ = αλ (αξ ).
So each line or ray solutions contains a unique ξ 0 such that H 2 (ξ 0 ) = ξ 0 . A Gröbner basis of the system H 2 (δ q) -δ q is given by {δ y 1 , δ z 1 , δ y 2 , δ z 2 }, so there is no non-trivial ray. This can also be seen in the following way: let δ q be a vector such that H 2 (δ q) = δ q. By the structure of the vector F(S), δ q satisfies δ z 1 = δ z 2 , and so h 2 (δ q) = h 2 (δ q) = 0, so H 2 (δ q) = 0, and which, by hypothesis, implies that δ q = 0.

We now consider non-isolated singular points of H 2 , that is the zeroes of H 2 . Since F(S) and G(S) are linearly independent, those points are exactly the zeroes of h 2 and h 2 , as described in the statement of the proposition.

Then we study the linearization of H π 2 in some neighborhood of these solutions in S 3 . First we consider vectors δ q such that δ z 1 = δ z 2 , we may perform the computations in the affine chart given by δ z 1 = 0, with coordinates δ

y 1 = δ y 1 /δ z 1 , δ y 2 = δ y 2 /δ z 1 , δ z 2 = δ z 2 /δ z 1 . The differential equation becomes d dt   δ y 1 δ y 2 δ z 1   = δ z 2 C( δ y 1 , δ y 2 , δ z 1 )
with C a polynomial vector field of degree 3. At δ z 1 = 1, its Jacobian is nilpotent:

  0 0 (ε -1)( δ y 1 (ε(ε -1) -1) + δ y 2 (1 -ε))(2Γ -γ) 2 γ 3 0 0 (ε -1)( δ y 1 (2ε(ε -2) + 1) + δ y 1 (ε -1))(2Γ -γ) 2 γ 3 0 0 0   Then we consider vectors δ q such that δ y 2 = (1 -ε)δ y 1 , δ z 1 = (1 -ε) 2 δ z 2 .
This time we use the chart δ z 2 = 0 with coordinates δ

y 1 = δ y 1 /δ z 2 , δ y 2 = δ y 2 /δ z 2 , δ z 1 = δ z 1 /δ z 2 .
As above, we compute the differential equation in this chart, linearize the resulting vector field, and evaluate this Jacobian at δ y 2 = (1 -ε) δ y 1 and δ z 1 = (1 -ε) 2 . This matrix has rank 2 and is diagonalizable with non-zero eigenvalues:

1 2 δ y 2 + 1 ± δ y 2 2 + (2ε -1) 2 ( δ y 2 + 1) -4(ε -1) 4 -2 δ y 2 + 1 .
Cubic approximation at N. We perform the same study at the North pole N. With expression [START_REF]Solving polynomial equations[END_REF], we can verify that the quadratic component of R(N) is 0. Indeed, F(N) = 0 and d 2 D (N) = 0. Further differentiating along q, we obtain

H 3 (δ q) := 1 6 d 3 f (N)(δ q, δ q, δ q) = 1 6 3d 2 D(N)(δ q, δ q)F(δ q) -d 3 D (N)(δ q, δ q, δ q)G(N)
Note that since we centered the coordinates at the North pole, F is linear in q, so dF(q) = F, and G is affine in q, so dG(q) is constant.

As in the previous subsection, we study the projection of the differential equation v = H 3 (v) on the sphere S 3 , and its equilibrium points, which form lines of non-isolated singular points and ray solutions.

Proposition 3.8. The cubic blow-up at the North pole N, for admissible values of the parameters, has two sets of ray solutions:

(1) the projective line

     δ z 1 = δ z 2 = 0 (ε -1)δ y 1 + δ y 2 = 1 Γ(2Γ -γ) √ 1 -ε ; (24)
(2) the quadric

     δ y 1 = δ y 2 = 0 ((ε -1)δ z 1 -δ z 2 ) 2 -(ε -2)δ z 1 δ z 2 = 1 γ 3 (2Γ -γ)(1 -ε)
.

and three sets of real non-isolated singularities:

(1) the plane

δ z 1 = δ z 2 δ y 1 (1 -ε) = δ y 2 (25)
(2) the plane

δ y 2 = (1 -ε)δ y 1 δ z 2 = (1 -ε) 2 δ z 1 (26)
(3) the surface defined by

     0 = Γ(ε -1)(2Γ -γ)δ y 2 1 + (2Γ -γ)Γδ y 1 δ y 2 + γ 2 (ε -1)δ z 2 1 -γ 2 (ε -1)δ z 1 δ z 2 0 = Γ(ε -1)(2Γ -γ)δ y 1 δ y 2 -γ 2 δ z 1 δ z 2 + (2Γ -γ)Γδ y 2 2 + γ 2 δ z 2 2 0 = δ y 1 δ z 2 + (ε -1)δ y 2 δ z 1 . (27)
For points on the line (24), the linearization of H π 3 is diagonal: the vectors (1, 0, 0) and (0, 0, 1) are eigenvectors, with the same eigenvalue, and the vector (0, 0, 1) is in the kernel.

For isolated singularities satisfying (25), the matrix is not diagonalizable, its Jordan form has the following structure:

  0 0 0 0 * 1 0 0 *  
For isolated singularities satisfying (26), the matrix is diagonalizable with 3 non-zero eigenvalues.

Proof. First we study ray solutions. Let ξ be a vector on a ray solution, such that

H 3 (ξ ) = λ ξ .
Let αξ be another vector on the same line (α ∈ R), one has

H 3 (αξ ) = α 3 H 3 (ξ ) = α 3 λ ξ = α 2 λ (αξ ).
So unlike in the quadratic case, a line of ray solutions contains 2 vectors ξ 1 , ξ 2 such that either

H 3 (ξ 1 ) = H 3 (ξ 2 ) = 1 or H 3 (ξ 1 ) = H 3 (ξ 2 ) = -1.
In order to study ray solutions, we compute a Gröbner basis of the system H 3 (δ q) -ιδ q = 0, ι 2 = 1, for an order eliminating α. We find that the basis contains

{δ y 1 δ z 1 , δ y 1 δ z 2 , δ y 2 δ z 1 , δ y 2 δ z 2 } so either δ y 1 = δ y 2 = 0 or δ z 1 = δ z 2 = 0.
If δ z 1 = δ z 2 = 0, computing a new Gröbner basis of the system, saturating with Γγ and δ y 1 (see Sec. 3.2.1), shows that δ y 1 and δ y 2 must satisfy

Γ 2 (ε -1)(2Γ -γ) 2 ((ε -1)δ y 1 + δ y 2 ) 2 2 -1 = 0.
Since ε -1 < 0, 2Γ -γ > 0 and Γ > 0, this defines 2 lines of real solutions given by

(ε -1)δ y 1 + δ y 2 = ± 1 Γ(2Γ -γ) √ 1 -ε
Those lines are equivalent in the projective space: each line of ray-solutions contains a vector in both lines.

For the second part of the proposition, as in the quadratic case, we study the linearization of

H π 3 : v = H 3 (v) - v, H 3 (v) v.
In the affine chart with δ y 1 = 0, with coordinates δ

z 1 = δ z 1 /δ y 1 , δ y 2 = δ y 2 /δ y 1 and δ z 2 = δ z 2 /δ y 1 , the differential equation v = H 3 (v) becomes d dt   δ z 1 δ y 2 δ z 2   = δ y 2 1 Q( δ z 1 , δ y 2 , δ z 2 )
with Q a polynomial vector field of degree 4. In this chart,

H π 3 becomes v = Q(v).
We conclude by evaluating the Jacobian of Q at the relevant points.

General case

Singularities of {D = H G = {H G , H F } = 0}.
Proposition 3.9. The set of singularities of

{D = H G = {H G , H F } = 0} is given, generically on authorized values of γ, Γ, by (1)
the plane z 1 = z 2 = z S ;

(2) the line z 1 = z 2 , y 1 = y 2 = 0;

(3) an irreducible variety of dimension 5.

If γ = Γ, solution 2 becomes a surface defined by z 1 = z 2 , y 2 = (1 -ε)y 1 .
Proof. By definition of D = {{H G , H F }, H G }, we want to study the zeroes of

     0 = p • G 0 = p • [G, F] 0 = p • [[G, F], G] (D)
The singularities of this variety is the set of points at which the matrix

  G [G, F] [[G, F], G]  
has rank at most 2 (as per Remark 3.3).

We encode that with the incidence variety

    G [G, F] [[G, F], G] v 1 v 2 v 3     •   L 1 L 2 L 3   =     0 0 0 1    
with new variables L = L 1 , L 2 , L 3 and random numbers v 1 , v 2 , v 3 . This gives us a system of 4 polynomial equations in the 10 unknowns y 1 , z 1 , y 2 , z 2 , L 1 , L 2 , L 3 , Γ, γ, ε. We eliminate L 1 , L 2 , L 3 from the ideal in order to recover the projection, and we saturate by 1 -ε.

We then compute a Gröbner basis for the elimination order y 1 > y 2 > z 1 > z 2 γ > Γ > ε, this basis contains 10 polynomials, some of which have factors with multiplicity greater than 1 or are divisible by 1 -ε or γ. We take the square-free form of this basis, and we saturate by 1 -ε and γ before computing a new Gröbner basis for the same order. The result is a set of 11 polynomials which includes

(z 1 -z 2 ) ((ε -1)P z 2 y 1 + P z 1 y 2 ) .
This determinant factors as

-2 (ε -1) 2 (2 Γ -1) (Γ -1) (z 1 -z 2 ) (ε -1) P z 2 y 1 + P z 1 y 2 Equilibrium positions. Lemma 3.11. All equilibrium points of Ż = D H F -D H G satisfying H G = {H G , H F } = 0 are contained in {D = D = 0}. Proof.
Recall that H F is defined as

H F = ∂ H F ∂ p -∂ F ∂ q = F -p • ∂ H F ∂ q
and H G is defined in the same way. Let z = (q, p) be a point such that

D(z) H F (z) -D (z) H G (z) = 0,
by looking at the first 4 components of this system, we see that the vectors F(q) and G(q) are colinear.

Introduce new variables X D and X D and consider the ideal generated by

• X D F(q) -X D G(q) = 0 • H G (z) = {H G , H F }(z) = 0 • D(z) H F (z) -D (z) H G (z) = 0 saturated by X D and X D .
Computing a Gröbner basis of this ideal (for any order) yields that this ideal is actually 1 , and so the associated system has no solution.

Hence, at an equilibrium point, either D(z) or D (z) has to be 0. Since both O (the center of the Bloch ball) and N (the north pole) are on the hyperplane z 1 = z 2 = 0, which is contained in {D = D = 0}, we may assume that the point z is neither O nor N. Hence, F(z) and G(z) are non-zero, and so D(z) = D (z) = 0.

Eigenvalues of the linearization. We consider the eigenvalue decomposition of the matrix

A = Jac(D H F -D H G )
on equilibrium point, given as the union of points satisfying Eq. ( 28) and [START_REF] Greuel | A Singular introduction to commutative algebra[END_REF]. Solutions of Eq. [START_REF] Faugère | FGb: A Library for Computing Gröbner Bases[END_REF]. If z 1 = z 2 , the matrix A has rank 2, and its characteristic polynomial is

T 6 T 2 - ε(ε -1)(ε -2) 8(Γ -γ) 2 P y 1 y 1 y 2 + P z 2 (2 P z 2 -γ)P y 2 2P y 1 (Γ -γ)y 1 T + (ε -1)γP y 1 (2Γ -γ) y 1 2
Solutions of Eq. ( 29). If y 2 = -

y 2 p z 1 (ε-1)p z 2
, the matrix A has rank 2, and its characteristic polynomial is

T 6 T 2 -4εy 2 (ε -1)(ε -2)(Γ -γ)T + 2γ(z 1 -z 2 )(ε -1) 2 (2Γ -γ)(Γ -γ) p z 1 2
.

The discriminant of the degree 2 factor factors as

16(Γ -γ) 2 (ε -1) 2 (a 6 (z 1 , y 2 ) -b 6 (z 1 , z 2 )) (a 6 (z 1 , y 2 ) + b 6 (z 1 , z 2 )) p 2 z 1 with a 6 (z 1 , y 2 ) = ε(ε -2)P z 1 y 2 b 6 (z 1 , z 2 ) = (2Γ -γ)(ε -1)γ(z 1 -z 2 )
which induces the following classification of the eigenvalues of A: We consider the set of singularities of the hypersurface {D = 0} lying in the Bloch ball. We show in the following that these singularities can be computed by means of explicit rational expressions. These explicit solutions show how complex is the structure of the algebraic variety corresponding to the singularities of {D = 0}. The number of real singularities can be computed directly without knowing the exact expressions by the algorithm in [START_REF] Bonnard | Determinantal sets, singularities and application to optimal control in medical imagery[END_REF] using an optimized variant of the Cylindric Algebraic decomposition algorithm. Statement of the main results. We keep the same letters for the sake of simplification. The new coordinates of the center O of the Bloch ball are O : (0, -1, 0, -1). We consider the set of polynomials S = {D, D y 1 , D y 2 , D z 1 , D z 2 }. We recall that we assume that 2 Γ 2 ≥ γ 2 > 0, and that (γ 1 , Γ 1 ) = (γ 2 , Γ 2 ). We normalize by means of Γ 1 = γ 1 = 1 (water), leaving free the two parameters Γ 2 , γ 2 corresponding to the second matter. Under this normalization, the hypothesis becomes (γ 2 , Γ 2 ) = (1, 1). Let Π be the plane {y 1 = y 2 = 0}.

• if |a(z 1 , y 2 )| > |b(
The following theorems were proved in [START_REF] Bonnard | Determinantal sets, singularities and application to optimal control in medical imagery[END_REF] using a real roots classification algorithm.

Theorem 3.12. Consider the 9 polynomials:

f 1 = Γ 2 -1 f 2 = 3 Γ 2 -2 γ 2 -1 f 3 = 3 Γ 2 2 -5 Γ 2 γ 2 + γ 2 2 + 2 Γ 2 -2 γ 2 + 1 f 4 = 2 Γ 2 2 -5 Γ 2 γ 2 + 2 γ 2 2 -2 Γ 2 + 3 γ 2 f 5 = 2γ 3 2 -(3Γ 2 + 11) γ 2 2 + 9Γ 2 + 6 -3Γ 2 2 γ 2 + 2 Γ 2 (Γ 2 + 2) (Γ 2 -1) f 6 = Γ 2 -2 γ 2 + 1 f 7 = 2 Γ 2 -γ 2 -1 f 8 = γ 2 -2 + Γ 2 f 9 = 2 Γ 2 2 -5 Γ 2 γ 2 + 2 γ 2 2 + 1.
The zeroes of their product divide the subset of R 2 defined by 2 Γ 2 > γ 2 > 0 into connected components where the cardi-

nality of V R ∩ π -1 (γ 2 , Γ 2 ) is constant.
Using a Cylindrical Algebraic Decomposition, one can make this statement more precise by expliciting this cardinality. Let ψ : (y 1 , z 1 , y 2 , z 2 ) → (-y 1 , z 1 , -y 2 , z 2 ) be the symmetry fixing Π = {y 1 = y 2 = 0}, and let us consider the semialgebraic sets (see Fig. 8):

G - 1 = {0 < γ 2 < 2 Γ 2 , Γ 2 < 1, f 2 > 0}, G + 1 = {0 < γ 2 < 2 Γ 2 , Γ 2 > 1, f 2 < 0, f 4 > 0}, G - 2 = {0 < Γ 2 < 1, f 6 > 0, f 3 < 0}, G + 2 = {Γ 2 > 1, f 6 < 0, f 5 > 0}, G 1 = G + 1 ∪ G - 1 , G 2 = G + 2 ∪ G - 2 , G = G 1 ∪ G 2 . f 1 f 2 f 3 f 3 f 4 f 4 f 5 f 6 1 1 1 1 1
∪ G + 1 (resp. G - 2 ∪ G + 2 )
. The circled numbers in each area correspond to the number of singularities in B for parameters in the area. Parameters in the red area are physically irrelevant.

Theorem 3.13. For all (γ 2 , Γ 2 ) such that 2 Γ 2 > γ 2 > 0, the center O of the Bloch ball B is a singularity of {D = 0}. And if (γ 2 , Γ 2 ) ∈ G, there exist at most two other singularities in the interior of the Bloch ball:

(1) if (γ 2 , Γ 2 ) ∈ G 1 there is one other singularity lying on Π ∩ B;

(2) if (γ 2 , Γ 2 ) ∈ G 2 there are two other singularities in B, ψ-symmetric, outside Π.

Computing the polynomials listed in Th. 3.12 takes about 10 s using the strategy described in Section 3.2.5, using Gröbner basis computations with the Maple FGb library for elimination steps. Computing points in each cell of the parameter space defined by these polynomials, in order to prove Th. 3.13, takes 50 s using the CAD implementation from the RegularChains package in Maple.

The configuration is illustrated in Figs. 8 and9. Observe that the number of singularities inside B is an invariant of the contrast problem. Two of the pairs of biological matters studied in [START_REF] Bonnard | Geometric numerical methods and results in the contrast imaging problem in nuclear magnetic resonance[END_REF], water-cerebrospinal fluid (normalized parameters [γ 2 = 5 4 , Γ 2 = 25 3 ]) and water-fat (normalized parameters [γ 2 = 25 2 , Γ 2 = 25]) correspond to points outside G, and their invariant is 1 in both cases (see Fig. 9). . We can give a more precise description of those points. Assuming that (γ 2 -Γ 2 )(2 Γ 2 -γ 2 -2) = 0, we define a point (with the coordinates (y 1 , z 1 , y 2 , z 2 )):

Ω : 0, Γ 2 -1 Γ 2 -γ 2 , 0, (1 -Γ 2 ) (2 Γ 2 -γ 2 ) (2 Γ 2 -γ 2 -2) (Γ 2 -γ 2 )
In the same coordinates system, O : (0, -1, 0, -1) is the center of the Bloch ball, while N : (0, 0, 0, 0) is the North pole. and water-fat (red square) and the set G (with the same conventions as in Fig. 8). For both these pairs, there is only 1 singularity in B.

Let

δ (γ 2 , Γ 2 ) = (Γ 2 -1) (2 Γ 2 -γ 2 ) (Γ 2 -2 γ 2 + 1) (Γ 2 -2 γ 2 + 2) (30) and η(γ 2 , Γ 2 ) = (5 Γ 2 -4 γ 2 + 1) (Γ 2 -γ 2 ) (Γ 2 + γ 2 -1) (Γ 2 -2 γ 2 + 2) . If η(γ 2 , Γ 2 ) = 0 we denote by Ω + a (resp. Ω - a ) the points of coordinates (y 1 , z 1 , y 2 , z 2 ) (resp. (-y 1 , z 1 , -y 2 , z 2
)), with:

y 1 = √ 3 δ (γ 2 , Γ 2 ) (5 Γ 2 -4 γ 2 + 1) (Γ 2 -γ 2 ) z 1 = (γ 2 + Γ 2 -1) (Γ 2 -1) (5 Γ 2 -4 γ 2 + 1) (γ 2 -Γ 2 ) y 2 = γ 2 √ 3 δ (γ 2 , Γ 2 ) 2 (γ 2 + Γ 2 -1) (Γ 2 -2γ 2 + 2) (Γ 2 -γ 2 ) z 2 = (2 Γ 2 -γ 2 ) (Γ 2 -1) 2 (γ 2 + Γ 2 -1) (γ 2 -Γ 2 )
.

Theorem 3.14. For values of the parameters (γ 2 , Γ 2 ) outside of a finite union of curves, the following points are singular points of {D = 0} in the subcase of water:

• O (always in the interior of the Bloch ball); . The curves involved in the decomposition of the parameter space (with the same conventions as in Fig. 8).

• Ω if (γ 2 -γ 2 )(2Γ 2 -γ 2 -2) = 0;
• Ω + a and Ω - a if η(γ 2 , Γ 2 ) = 0 and which are real if δ (γ 2 , Γ 2 ) > 0. The remainder of this section will be devoted to the proof of Theorem 3.14. Beforehand, however, let us state this last theorem, which makes the connection between the results of Theorems 3.13 and 3.14, and in particular implies that generically, the points listed in Theorem 3.14 are the only singularities of {D = 0} in the Bloch ball. Theorem 3.15. For values of the parameters (γ 2 , Γ 2 ) outside of a finite union of curves, the singularities of {D = 0} inside the Bloch ball in the case of water are:

• O; • Ω if and only if (γ 2 , Γ 2 ) ∈ G 1 ;
• Ω + a and Ω - a are real and inside the Bloch ball if and only if (γ 2 , Γ 2 ) ∈ G 2 .

Remark 3.16. Note that generically the number of singularities of {D = 0} in the Bloch ball is at most 3, because the domains G 1 and G 2 are disjoint.

Proof. The point Ω is within the ball if and only if its coordinates satisfy

-2 < z 1 , z 2 < 0.
The inequalities in z 1 are equivalent to the disjonction of these two sets of inequalities Note that in each case, the first and third condition imply the second one.

     Γ 2 -1 > 0 Γ 2 < γ 2 3Γ 2 -2γ 2 -1 = f 2 < 0 or      Γ 2 -1 < 0 Γ 2 > γ 2 f 2 > 0 0.6 0.8 1
In the first case, since Γ 2 > γ 2 , the inequalities 0 > z 2 > -2 are equivalent to

2Γ 2 -γ 2 -2 < 0 f 4 < 0
Both inequalities are satisfied by assumption, recovering the defining inequalities of G - 1 . In the second case (Γ 2 > γ 2 ), the inequalities 0

> z 2 > -2 are equivalent to 2Γ 2 -γ 2 -2 > 0 f 4 > 0
The first inequality is satisfied whenever the second is, so have recoved the defining inequalities of G + 1 , proving the first point.

The second case is studied in the same way.

Proof of Theorem 3.14.

Proof. We consider separately those solutions which are on the symmetry plane {y 1 = y 2 = 0} (Proposition 3.17) and those which are not (Proposition 3.18). The result is an immediate consequence of those propositions.

Solutions on the symmetry plane {y 1 = y 2 = 0}.

Proposition 3.17. The singularities of D which are on the symmetry hyperplane Π are:

• if (2 Γ 2 -γ 2 )(Γ 2 -γ 2 )(2 Γ 2 -γ 2 -2) = 0:
the center O of B, and the point Ω;

• if 2 Γ 2 = γ 2 : O, Ω and N; • if Γ 2 = 1, and (Γ 2 -γ 2 )(2 Γ 2 -γ 2 -2) = 0: O. • if Γ 2 = 1 and γ 2 = 1: O and Ω = N.
Proof. We split the proof in seven cases, finding the solutions such that z

2 = -1, 2Γ 2 = γ 2 , Γ 2 = 1, z 2 = 0, 2 Γ 2 = γ 2 , Γ 2 = γ 2 , 2 Γ 2 -γ 2 -2 = 0
, and the remaining general case.

(1) Case z 2 = -1. We compute a basis of S ∪ {y 1 , y 2 , z 2 + 1} with the lexicographical order y

1 y 2 z 1 z 2 γ 2 Γ 2 .
The polynomial γ 2 2 (z 1 + 1) is present, and since γ 2 = 0, we get z 2 = -1 and derive directly that the unique solution is O.

(2) Case 2 Γ 2 = γ 2 . We compute a basis of S ∪ {y 1 , y 2 , 2 Γ 2 -γ 2 } with the same order y 1 y 2 z 1 z 2 γ 2 Γ 2 , and the result follows immediately. (3) Case Γ 2 = 1. We found in the basis of S ∪ {y 1 , y 2 , Γ 2 -1} with order y

1 y 2 z 1 z 2 γ 2 Γ 2 , the polynomial z 2 γ 2 2 (γ 2 -1) 2 (z 2 + 1)
. The case z 2 + 1 is already treated, while the case γ 2 = 0 is excluded. The case γ 2 = 1 is also excluded, since that would mean (γ 2 , Γ 2 ) = (1, 1). In the case z 2 = 0, the polynomial 2γ 3 2 z 1 (γ 2 -2) appears in the basis of S ∪ {y 1 , y 2 , Γ 2 -1, z 2 }. Since the case γ 2 = 2 = 2Γ 2 has already been studied, this implies z 1 = 0 hence the north pole N is the unique solution. Note that for Γ 2 = 1, N and Ω coincide. (4) Case z 2 = 0. We compute a basis with the lexicographic y 1 y 2 z 1 z 2 γ 2 Γ 2 . The first polynomial is (Γ 2 -1) 2 (2 Γ 2 -γ 2 ). Since 2 Γ 2 -γ 2 = 0 and Γ 2 = 1 are already known, there is no extra solution.

(5) Case Γ 2 = γ 2 . The basis contains the polynomial Γ 2 (Γ 2 -1) 2 (z 2 + 1), which leads to z 2 = -1 (already discussed). ( 6) Case 2 Γ 2 -γ 2 -2 = 0. The basis contains the polynomial (Γ 2 -1) 2 (z 2 + 1), as in the previous case. [START_REF] Bonnard | Geometric and numerical methods in the contrast imaging problem in nuclear magnetic resonance[END_REF] Let us finally search the solutions such that z 2 (z 2 + 1)(γ 2 -Γ 2 )(2 Γ 2 -γ 2 -2) = 0. We consider the system S ∪ {y 1 , y 2 , z 2 ξ 1 -1, (z 2 + 1)ξ 2 -1, (Γ 2 -γ 2 )λ 1 -1, (2 Γ 2 -γ 2 -2)λ 2 -1}, and compute a basis with order y 1 y 2 z 1 z 2 ξ 1 ξ 2 γ 2 Γ 2 λ 1 λ 2 . The last four polynomials of the basis are:

Γ 2 λ 1 + 2 λ 2 λ 1 + λ 1 -2 λ 2 + z 2 , -Γ 2 λ 1 + λ 1 + z 1 , y 2 , y 1 .
This implies that there is at most one solution under the above conditions. By plugging into the system, one verifies that Ω is a solution satisfying these conditions, hence it is the unique solution.

Solutions outside the symmetry plane {y 1 = y 2 = 0}. Proposition 3.18. For values of the parameters (γ 2 , Γ 2 ) outside of a finite union of curves, the solutions of S outside Π are given by means of a two parameters family of Π-symmetric points (Ω + a (γ 2 , Γ 2 ), Ω - a (γ 2 , Γ 2 )), which can be either real or complex depending on the sign of δ (γ 2 , Γ 2 ).

Proof of Proposition 3.18. . Proof. The case y 2 = 0, y 1 = 0 is treated in Lemma 3.19. We assume now that y 2 = 0. We transform the system S = {D, D y 1 , D y 2 , D z 1 , D z 2 } by means of the introduction of the new variables (t, h) such that y 1 = t y 2 , and γ 2 = h Γ 2 . After simplification by y 2 , one observes that the powers of y 2 in the simplified system are all even. So we perform the change of variables Y 2 = y 2 2 , and denote by S this new set of polynomials in the variables (h, Γ 2 , z 1 , z 2 ,t,Y 2 ).

We first compute a Gröbner basis, considering the polynomials in Q(γ 2 , Γ 2 )[y 1 , z 1 , y 2 , z 2 ] with a lexicographical order, and find directly the rational expressions

z 1 = (1 -(1 + h) Γ 2 ) (Γ 2 -1) (4 Γ 2 h -5 Γ 2 -1) Γ 2 (h -1) z 2 = 1 2 (Γ 2 -1) (h -2) (1 -(1 + h) Γ 2 ) (h -1) Y 2 = 3 4 h 2 Γ 2 (Γ 2 -1) (h -2) (2 Γ 2 h -Γ 2 -1) (2 Γ 2 h -Γ 2 -2) (Γ 2 h + Γ 2 -1) 2 (h -1) 2 t = 2 (Γ 2 h + Γ 2 -1) (2 Γ 2 h -Γ 2 -2) (4 Γ 2 h -5 Γ 2 -1) h Γ 2 (SoluGen)
which correspond to the points Ω + a , Ω - a . These points are symmetric with respect to Π and real if and only if

(Γ 2 -1) (2 Γ 2 -hΓ 2 ) (Γ 2 -2 hΓ 2 + 1) (Γ 2 -2 hΓ 2 + 2) ≥ 0
Recall that hΓ 2 = γ 2 , so we have recovered the definition of δ (30). We define

η(Γ 2 , h) = (4 Γ 2 h -5 Γ 2 -1) (h -1) (2 Γ 2 h -Γ 2 -2) (1 -(1 + h) Γ 2 ) ,
the lowest common multiple of the denominators of the solutions, so that those paired solutions are valid on a Zariski open set (including at least the complementary of {(γ 2 , Γ 2 ) s.t. η(Γ 2 , h) = 0}.

Lemma 3.19. If γ 2 = 2Γ 2 and γ 2 = 0, the only solutions of the system S ∪ {y 2 } are the solutions on Π given in Proposition 3.17. In particular, all these solutions are such that y 1 = 0.

Proof. The second polynomial of a basis of S∪{y 2 } with respect to the order z 1 z 2 y 1 y 2 γ 2 Γ 2 is γ 2 2 y 1 (2 Γ 2γ 2 ).

General case (3 projective parameters)

We conclude by giving classification results for the classification of the real singularities of {D = 0} in the Bloch ball, in the general case. We now have 4 parameters (γ 1 , Γ 1 , γ 2 , Γ 2 ), which we may reduce to 3, normalizing by γ 1 = 1. The constraints are 2 Γ 1 ≥ 1, 2 Γ 2 ≥ γ 2 > 0, and we exclude the case of a single spin (Γ 2 , γ 2 ) = (Γ 1 , 1).

In this case, real roots classification algorithms [START_REF] Bonnard | Determinantal sets, singularities and application to optimal control in medical imagery[END_REF] allow to compute a generic classification of the parameter space according to the number of singularities. Theorem 3.20 ( [10]). Splitting the subset of R 3 defined by 2Γ 2 ≥ γ 2 > 0 and 2Γ 1 > 1 into open subsets where the number of real singularities of {D = 0} in the Bloch ball in the fibers is constant, can be done by cutting out 12 irreducible surfaces, consisting of 5 planes, 3 quadrics, 2 surfaces of degree 9, and one of degree 14.

Note that this computation is more expensive than in the case of water (2.5 h), and the result is also larger: the degree 14 polynomial has 451 monomials.

Using critical points computations on the hypersurface defined by the product of these polynomials, we are able to compute at least one point per connected component of the complementary of this hypersurface. This second step of the computations is also more computationally expensive than in the case of water, it takes 48 h and results in 10 109 points.

We can then use these points to count the number of singularities in each cell. However, unlike the case of water, it is not possible to characterize the cells using only the signs of these polynomials. In the experimental setting, the parameters for the second matter were γ 2 0.89 and Γ 2 36.6. Scaling the graphic down to fit this point (Figure 13) shows that some of the polynomials can be eliminated for values in this order of magnitude. Several of the classification polynomials vanish at (γ 2 , Γ 2 ) = (γ 1 , Γ 1 ), which is a natural singularity. One can also note that unlike in the case of water, the areas where there is only 1 singularity are quite far from the actual values of the parameters.

FIGURE 1 .

 1 FIGURE 1. Experimental results: The inner circle shape sample mimics the deoxygenated blood, where T 1 = 1.3s and T 2 = 50ms; the outside moon shape sample corresponds to the oxygenated blood, where T 1 = 1.3s and T 2 = 200ms. (left) Without control, (right) Optimized contrast.

FIGURE 3 .FIGURE 4 .F

 34 FIGURE 3. Dynamics under the control u ≡ 0

FIGURE 5 .

 5 FIGURE 5. Algebraic sets involved in saturation singular arcs

3. 3 . 2 .Proposition 3 . 4 .

 3234 Transfer time not fixed (H F = 0) Singularities of {D = 0}. The set of points satisfying D = ∂ D ∂ y 1

2 and a 5 Lemma 3 . 6 .

 2536 (y 2 , z 2 ) = 0 y 2 (2Γ -γ)Γy 2 2 + γ 2 z 2 2which, as in the casey 1 = z 1 = 0, is only y 2 = z 2 = 0 if 2Γ > γ.This partial solution completes into y 1 = y 2 = z 1 = z 2 = 0, which was already known.Equilibrium positions. The equilibrium points of Ẋ = DF -D G are all contained in {D = D = 0}.

2 FIGURE 8 .

 28 FIGURE 8. Curves involved in the definition of the semi-algebraic set G. The blue (resp. green) sample points correspond to points in G -1 ∪ G + 1 (resp. G - 2 ∪ G + 2 ). The circled numbers in each area correspond to the number of singularities in B for parameters in the area. Parameters in the red area are physically irrelevant.

2 FIGURE 9 .

 29 Fig.6

2 FIGURE 10

 210 Fig.9

Theorem 3 . 21 .

 321 For values of the parameters (γ 1 , Γ 1 , γ 2 , Γ 2 ) outside of a finite union of hypersurfaces, the number of singularities of {D = 0} in the Bloch ball can be 1, 2, 3, 4 or 5.

FIGURE 12 .

 12 FIGURE 12. All polynomials involved in the classification in the in vivo case

FIGURE 13 .

 13 FIGURE 13. Classification regarding the second parameter in the in vivo case

  z 1 , z 2 )|: 2 single real eigenvalues;• if |a(z 1 , y 2 )| = |b(z 1 , z 2 )|: 1 double real eigenvalue; • if |a(z 1 , y 2 )| < |b(z 1 , z 2 )|: 2 single complex eigenvalues.

	3.4. Algebraic computations for the contrast problem
	3.4.1. Subcase of water (2 projective parameters)

The necessary computations are detailed in Maple code files which are available on http://mercurey.gforge.inria.fr/

* Work supported by ANR Project DFG Explosys (ANR-14-CE35-0013-01; GL203/9-1) and ERC-Advanced Grant for the TAMING project

If δ y 1 = δ y 2 = 0, the same technique shows that δ z 1 and δ z 2 must satisfy

which defines 1 quadric

.

We now consider non-isolated singularities, that is zeroes of H 3 . To this end, we compute a Gröbner basis of the system H 3 (δ q) = 0, saturating by γ -Γ, γ, Γ, ε -1 and 2Γ -γ. Factoring the results, it appears that the solutions split into 5 cases:

(1)

(5) otherwise. In the case 1, a Gröbner basis is given by

and the solutions form the plane [START_REF] Faugère | A new efficient algorithm for computing Gröbner bases (F 4 )[END_REF].

In the case 2, saturating by δ z 1 -δ z 2 , a Gröbner basis is given by

  which has no real non-zero solution for admissible values of the parameters.

In the case 3, saturating by δ z 1 -δ z 2 and δ z 1 + δ z 2 , a Gröbner basis is given by

which has no real non-zero solution for admissible values of the parameters.

In the case 4, saturating by δ z 1 ± δ z 2 and δ z 1 , a Gröbner basis is given by

and the solutions form the plane [START_REF] Faugère | A new efficient algorithm for computing Gröbner bases without reduction to zero (F 5 )[END_REF].

Finally, for the case 5, we compute a Gröbner basis, saturating by all the previous conditions. This basis is

The fourth polynomial is a combination of the other 3, and the solutions form the surface [START_REF] Faugère | Critical points and gröbner bases: the unmixed case[END_REF].

First, we add z 1 -z 2 to the ideal. Once again, we compute a Gröbner basis for the elimination order y 1 > y 2 > z 1 > z 2 γ > Γ > ε, take its square-free form, and recompute a Gröbner basis. The result contains the polynomial

The solutions decompose into 2 algebraic sets, defined by

and (adding (ε -1)y 1 + y 2 to the ideal and saturating by P z 2 , ε and ε -2)

The latter is generically (if Γ = γ) defined by

or if Γ = γ (that is if the spin we consider is water), by

Then, starting again with the whole ideal, we add (ε -1)P z 2 y 1 + P z 1 y 2 to the ideal and saturate by z 1 -z 2 . The result, generically on (Γ, γ, ε), is an irreducible surface.

Locus of {D

Proposition 3.10. The solutions form the union of the hyperplane defined by z 1 = z 2 [START_REF] Faugère | FGb: A Library for Computing Gröbner Bases[END_REF] and the hypersurface

The projection of these points onto the space (y 1 , z 1 , y 2 , z 2 ) is given by the vanishing of the determinant ∆ , defined as