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

Bernard Bonnard, Olivier Cots, Jean-Charles Faugère, Alain Jacquemard,
Jérémy Rouot, Mohab Safey El Din, Thibaut Verron

► To cite this version:

Bernard Bonnard, Olivier Cots, Jean-Charles Faugère, Alain Jacquemard, Jérémy Rouot, et al. 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. 2017. [hal-01556806](https://hal-inria HAL Id: hal-01556806)

HAL Id: hal-01556806

https://inria.hal.science/hal-01556806

Submitted on 5 Jul 2017

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# 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 *

B. BONNARD<fnref n="1" />,<fnref n="2" />, O. COTS<fnref n="3" />, J.-C. FAUGÈRE<fnref n="4" />, A. JACQUEMARD<fnref n="1" />, J. ROUOT<fnref n="5" />,
M. SAFEY EL DIN<fnref n="4" /> AND T. VERRON<fnref n="3" />

**Abstract.** 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.

**1991 Mathematics Subject Classification.** 49K15,14Q20,81Q93.

July 4, 2017.

## 1. INTRODUCTION

Optimal control algorithms were introduced in Nuclear Magnetic Resonance (NMR) to improve the control field very recently [20] 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 [33,48]. More recently, under the impulse of S. Glaser, the combination of geometrical optimal control based on the Maximum Principle [43] and related numerical algorithms (gradient methods [34],

<footnote>Keywords and phrases: optimal control,singular trajectories,symbolic computation,Bloch equations,Magnetic Resonance Imaging</footnote>

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

<footnote n="1">¹ Institut de Mathématiques de Bourgogne 9 avenue Savary, 21078 Dijon, France;
e-mail: bernard.bonnard@u-bourgogne.fr & alain.jacquemard@u-bourgogne.fr</footnote>

<footnote n="2">² Inria Sophia Antipolis</footnote>

<footnote n="3">³ Toulouse Univ., INP-ENSEEIHT-IRIT, UMR CNRS 5505, 2 rue Camichel, 31071 Toulouse, France;
e-mail: olivier.cots@enseeiht.fr & thibaut.verron@enseeiht.fr</footnote>

<footnote n="4">⁴ Sorbonne Universités, UPMC Univ Paris 06, CNRS, INRIA Paris Center, LIP6, Equipe PolSys, F-75005, Paris, France;
e-mail: jean-charles.faugere@inria.fr & mohab.safey@lip6.fr</footnote>

<footnote n="5">⁵ LAAS-CNRS, Univ. de Toulouse, LAAS, 31031 Toulouse cedex 4, France; e-mail: jeremy.rouot@laas.fr</footnote>
<figure><img src="image_2.png" /><figcaption>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.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIGURE 2. Contrast optimization in a *in vivo* setting. Species : brain – parietal muscle</figcaption></figure>

shooting and continuation methods [21]) lead to sophisticated results starting from a complete solution to the time-minimal saturation of a single spin [38], and applications to the contrast problem in magnetic resonance imagery (MRI), see [6,37].

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 [9, 45, 49] (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* OXYZ whose evolution satisfies the so-called *Bloch equation*

$$ \frac{dM}{dt} = \gamma M \wedge B + R(M), \qquad (1) $$

where $\gamma$ 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) = \left( \frac{M_X}{T_2}, \frac{M_Y}{T_2}, \frac{M_Z - M_0}{T_1} \right) $$

where $M_0$ is the equilibrium magnetization, which can be normalized to 1 using the rescaling $M \mapsto 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) = -\gamma B_Y, v(t) = -\gamma 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 Oxyz with $M(t) = S(t)q(t)$, $q = (x, y, z)$, $S(t) = \exp(\omega t \Omega_z)$, $\Omega_z = \begin{bmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}$ as:

$$ \frac{d}{dt} \begin{bmatrix} x \\ y \\ z \end{bmatrix} = \begin{bmatrix} -1/T_2 & -\Delta\omega & u_2 \\ \Delta\omega & -1/T_2 & -u_1 \\ u_2 & u_1 & -1/T_1 \end{bmatrix} \begin{bmatrix} x \\ y \\ z \end{bmatrix} + \begin{bmatrix} 0 \\ 0 \\ 1/T_1 \end{bmatrix} \qquad (2) $$
where $\Delta\omega = \omega_0 - \omega$ is the resonance offset, $\omega_0 = -\gamma B_0$ is the resonance frequency and the RF-control field is represented as

$$ \begin{cases} u_2 = u \cos(\omega t) - v \sin(\omega t) \\ u_1 = u \sin(\omega t) + v \cos(\omega t) \end{cases} \qquad (3) $$

which preserves the control bound $u_1^2 + u_2^2 \le 1$.

Finally, in the *moving frame* the system takes the normalized form

$$ \begin{cases} \frac{dx}{dt} = -\Gamma x - \Delta \omega y + u_2 z \\ \frac{dy}{dt} = \Delta \omega x - \Gamma y - u_1 z \\ \frac{dz}{dt} = \gamma (1-z) + u_1 y - u_2 z. \end{cases} \qquad (4) $$

The resonant case is when $\Delta\omega = 0$ which will be assumed in the sequel of this article. The system is then written as

$$ \frac{dq}{dt} = F(q) + u_1 G_1(q) + u_2 G_2(q), u_1^2 + u_2^2 \le 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

$$ \frac{dq}{dt} = F(q) + uG(q), q = (y, z), |u| \le 1. \qquad (5) $$

Also note that the Bloch ball $|q| \le 1$ is invariant for the dynamics provided that the parameters satisfy the physical constraint $2\Gamma \ge \gamma \ge 0$.

Our study will concern the following optimal control problems.

First of all, consider a single spin system described by (5). 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 [36].

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

$$ \dot{z} = \vec{H}_F(z) + u \vec{H}_G(z), z = (q, p) \qquad (6) $$

$H_F$ and $H_G$ being the Hamiltonians $p \cdot F(q)$ and $p \cdot G(q)$, and moreover the optimal control is given by the maximization condition

$$ H_F(z) + uH_G(z) = \max_{|v| \le 1} H_F(z) + vH_G(z) \qquad (7) $$

and it is a concatenation of *bang arcs* $\delta_+$, $\delta_-$ 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 = -\frac{\gamma}{2(\Gamma - \gamma)}$ and the interesting situation is when $2\Gamma > 3\gamma$ so that this horizontal line intersects the Bloch ball. More precisely, following [36] we recall the following:

**Proposition 1.1.** Provided that $2\Gamma > 3\gamma$, the time-minimal solution to the saturation problem in the experimental situation is of the form $\delta_+ \delta_{SH} \delta_+^B \delta_{SV}$ where $\delta_{SV}$ and $\delta_{SH}$ are vertical and horizontal singular arcs and the intermediate arc $\delta_-^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| \to \infty$ when $y \to 0$.

(Note: in experimental situation, the line ($\delta_{SH}$) is accessible from the North pole).

All the information about this study is reported on Fig. 4. Note that in the case $2\Gamma \le 3\gamma$, the standard inversion sequence $\delta_+ \delta_{SV}$ is time-minimal.
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Dynamics under the control $u \equiv 0$</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 4. Time-minimal saturation: left: $2\Gamma > 3\gamma$, right: $2\Gamma \le 3\gamma$</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 5. Algebraic sets involved in saturation singular arcs</figcaption></figure>

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 (5), with respective parameters $(\Gamma_1, \gamma_2)$, $(\Gamma_2, \gamma_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):

$$ \left\{ \begin{array}{l} \frac{d q_1}{d t} = F(q_1) + u G(q_1) \\ \\ \frac{d q_2}{d t} = F(q_2) + u (1 - \varepsilon) G(q_2) \end{array} \right. \qquad |u| \le 1 $$

where $\varepsilon$ 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)) = \alpha |q_1(t_f)|^2 - \beta |q_2(t_f)|^2 $$

where $\alpha, \beta$ 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(\bullet)} c(q(t_f))$ subject to $\dot{q} = F + uG$, $|u| \le 1$, with some prescribed boundary conditions $q(t_f) \in M$, with $M$ the terminal manifold defined by $f(q) = 0$. Again the maximum principle described by (6), (7) 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 \frac{\partial c}{\partial q}(q(t_f)) + \delta \frac{\partial f}{\partial q}(q(t_f)), p_0 \le 0, \delta: \text{constant vector} \quad (8)$$

Clearly, the optimal control problems boil down to analyze the so-called singular extremals given by

$$\dot{z} = H_F(z) + uH_G(z), H_G = 0 \quad (9)$$

where $u$ has to be admissible, that is $|u| \le 1$. The surface $\Sigma$ 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) = \frac{\partial F}{\partial q}(q)G(q) - \frac{\partial G}{\partial q}F(q)$$

and the Poisson bracket of $H_F, H_G$ is given by

$$\{H_F, H_G\}(z) = dH_F(\vec{H}_G(z)) = p \cdot [F, G](q).$$

Hence differentiating $t \mapsto H_G(z(t)) = 0$, one gets

$$z \in \Sigma' : H_F(z) = \{H_F, H_G\}(z) = 0 \quad (10)$$

and the singular control is defined by solving

$$\{\{H_G, H_F\}, H_F\}(z) + u\{\{H_G, H_F\}, H_G\} = 0. \quad (11)$$

If $F, G$ are real analytic, this defines a *meromorphic equation* of the form

$$\frac{dz}{dt} = \frac{X(z)}{g(z)}, z \in \Sigma'. \quad (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 $((\Gamma_1, \gamma_1), (\Gamma_2, \gamma_2))$ related to the classification of solutions of (12).

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 [7].

The organization of this article is the following. In Section 2, we present based on [3] 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).

<footnote n="4">⁴The necessary computations are detailed in Maple code files which are available on http://mercurey.gforge.inria.fr/</footnote>
## 2. 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.

### 2.1. The classification problem of singular extremals

In this section, we recall the results of [3] which will justify all our algebraic computations.

**Definition 2.1.** Let E and F be two $\mathbb{R}$-vector spaces, and let G be a group acting linearly on E and F. An homomorphism $\mathcal{X}: G \to \mathbb{R}^0$ is called a character. Let $\mathcal{X}$ be a character. A semi-invariant of weight $\mathcal{X}$ is a map $\lambda: E \to \mathbb{R}$ such that for all $g \in G$ and all $x \in E$, $\lambda(g.x) = \mathcal{X}(g)\lambda(x)$. It is an invariant if $\mathcal{X}=1$. A map $\lambda: E \to F$ is a semi-covariant of weight $\mathcal{X}$ for all $g \in G$ and for all $x \in E$, $\lambda(g.x) = \mathcal{X}(g)g.\lambda(x)$. It is called a covariant if $\mathcal{X}=1$.

Notations and definitions. Let $U \subset \mathbb{R}^n$ be an open subset and let $V^\omega(U)$ be the set of real analytic vector fields on $U$ identified to maps from $U$ to $U$. We denote by $C^\omega(U, \mathbb{R})$ the set of analytic maps from $U$ to $\mathbb{R}$. Let $G_d$ be the group of $C^\omega$ diffeomorphisms of $U$. The coordinates of $TU \simeq U \times \mathbb{R}^n$ are denoted $z = (q, p)$ and we endow $T^*U$ with its canonical symplectic structure defined by $\Omega = \sum_{i=1}^n dq_i \wedge dp_i$. Let $Z$ be a $C^\omega$ vector field on $U$, the Hamiltonian lift of $Z$ is defined by the Hamiltonian $H_Z(z) = p \cdot Z(q)$. If $\varphi \in G_d: q = \varphi(Q)$, one can lift $\varphi$ into a symplectic diffeomorphism $\vec{\varphi}$, called a *Mathieu transformation* defined by:

$$ \vec{\varphi}: q = \varphi(Q), p = P \frac{\partial \varphi^{-1}}{\partial Q} $$

writing $p$ and $P$ as row vectors. A system $\dot{q} = X + uY$, $X, Y \in V^\omega(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 $\varphi \in G_d$ and a feedback $u = \alpha(q) + \beta(q)v$, $\alpha, \beta \in C^\omega(U, \mathbb{R})$, $\beta(q)$ non-zero for each $q$, such that

$$ \begin{cases} X' = \varphi * X + \varphi * Y.\alpha \\ Y' = \varphi * Y.\beta \end{cases} $$

where $\varphi * Z$ is the image of $Z$ defined by $\varphi * Z = \frac{\partial \varphi^{-1}}{\partial q}(Z \circ \varphi)$. This action defines a group structure on the set of triplets $(\varphi, \alpha, \beta)$ and this group is the *feedback group* and denoted by $G_f$.

Singular extremals of order 2. Let $(F, G) \in S$, the singular extremals of order 2 are the solutions of

$$ \dot{z} = \vec{H}_F(z) + u_s \vec{H}_G(z), z \in \Sigma: H_G = 0 $$

such that $\{\{H_G, H_F\}, H_G\}(z)$ is never vanishing and $u_s$ is given by (11) as

$$ u_s = - \frac{\{\{H_G, H_F\}, H_F\}(z)}{\{\{H_G, H_F\}, H_G\}(z)} \quad (13) $$

and they are contained in the surface $\Sigma': 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 $\hat{H}(z)$. One has the following, see [3].

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

$$ \dot{q} = \frac{\partial \hat{H}}{\partial p}, \dot{p} = -\frac{\partial \hat{H}}{\partial q}, (q, p) \in \Sigma' \setminus S \quad (14) $$

Moreover, an explicit representation is as follows:

**Lemma 2.3.** Assume that $G$ and $[G, F]$ are never colinear on $U$. Then $(\Sigma' \setminus S, \Omega')$ is a symplectic manifold, $\Omega'$ denoting the restriction of the standard symplectic form and $\hat{H}$ on $\Sigma' \setminus S$ is a $C^\omega$ 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 $(\varphi, \alpha, \beta) \in G_f$, then the feedback $(\alpha, \beta)$ acts trivially and $\varphi$ acts by the change of symplectic coordinates $\vec{\varphi}$.

The key result in our analysis is the following [3].

**Theorem 2.4.** Let $\lambda$ be the map $(F, G) \to (\hat{H}, \Sigma', S)$ (differential equation defined by (14)). Then for the respective action of the feedback group $G_F$, $\lambda$ is a covariant.

Moreover, in relation with our geometric study, it is worth extracting two more covariants, where $G_F$ acts on functions by $\vec{\varphi}$ 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 $\Sigma'$, 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\} \qquad (15)$$

and observes that the map: $(F, G) \mapsto \{\{H_G, H_F\}, H_G\}$ is a semi-covariant if $z \in \Sigma'$.

**Theorem 2.6.** The map $\lambda'$; $(F, G) \mapsto (\vec{H}\{\{H_G, H_F\}, H_G\}, \Sigma')$ is a semi-covariant.

Note that the associated vector field is $C^\omega$ 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 [4] into two subsets, recalling that the level sets $\hat{H} = c$ are invariants.

Exceptional case. It corresponds to $\hat{H} = 0$ and the corresponding singular extremals of order 2 are called exceptional. Their role in Mayer optimal problems is analyzed in [4].

In our situation, in dimension 4, the constraint $\hat{H} = 0$ leads to the additional relation $H_F = p \cdot F = 0$. Hence the singular control is defined by

$$\begin{cases} H_F = H_G = \{H_G, H_F\} = 0 \\ \{\{H_G, H_F\}, H_F\} + u_s\{\{H_G, H_F\}, H_G\} = 0 \end{cases}$$

Since $p \in \mathbb{R}^4 \setminus \{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 = -\frac{D'(q)}{D(q)}$ and the vector field is denoted $X^e(q) = F(q) - \frac{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_r^e = DF - D'G. \qquad (16)$$

In this case the $\vec{\varphi}$-action is reduced to the standard action of diffeomorphisms $\varphi$ on $U$ acting on tensors, that is, on vector fields $V$ by image $\varphi * V$ and on mappings $f: U \to \mathbb{R}$ by composition.
Generic case. It corresponds to $\hat{H} = c \neq 0$ and by homogeneity one can take $c = 1$.

Let $\mathcal{D}' = \{\{H_G, H_F\}, H_G\}, \mathcal{D} = \{\{H_G, H_F\}, H_F\}$, the singular control $u_S$ is defined by

$$ \begin{cases} H_G = \{H_G, H_F\} = 0 \\ \mathcal{D} + u_S \mathcal{D}' = 0. \end{cases} $$

## 2.2. Lie brackets computations

Each spin system can be lifted on the semi-direct product $S = GL(3, \mathbb{R}) \times \mathbb{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

* $F = (A, a), A = \text{diag}(-\Gamma, -\Gamma, -\gamma)$ and $a = (0, 0, \gamma)$

* $G = Cq$ where $C$ is the antisymmetric matrix $C = E_{32} - E_{23}$ with $E_{ij} = (\delta_{ij})$.

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 [32] for the details.

The Lie algebra $gl(n, \mathbb{R})$ of $n \times n$ matrices decomposes into $cI_n \oplus sl(n, \mathbb{R})$ where $sl(n, \mathbb{R})$ is the Lie algebra of matrices with zero trace. This algebra can be written as the direct sum $a \oplus b$ where $a$ is the Cartan subalgebra of diagonal matrices and $b = \bigoplus_{i \neq j} \mathbb{R}E_{ij}$. Moreover we have the Cartan decomposition $sl(n, \mathbb{R}) = so(n) \oplus h$ where $h$ is the subspace of symmetric matrices. If $A = \text{diag}(\lambda_1, \dots, \lambda_n) \in a$ then $[A, E_{ij}] = (\lambda_i - \lambda_j)E_{ij}$ and $\lambda_i - \lambda_j$ will form the nontrivial spectrum of $adA$ (where $adA$ is defined as $adA(B) = [A, B]$) with corresponding eigenvectors $E_{ij}$.

Application. Restricting to the $q = (y, z)$ space we have $A = \text{diag}(-\Gamma, -\gamma)$ and $a = (0, \gamma)$ and $A = cI_2 + \text{diag}(\lambda, -\lambda)$ where $c = -\frac{\Gamma+\gamma}{2}$ which is zero if and only if $\gamma = \Gamma = 0$ and $\lambda = \delta/2$ where $\delta = \gamma - \Gamma$. If $\delta \neq 0$, the nontrivial spectrum is $(\delta, -\delta)$.

In NMR, we have the following

**Lemma 2.7.** The case $\gamma = \Gamma$ is the case of water species.

Otherwise, an easy computation gives that the Lie algebra generated by $((A, a), (C, 0))$ is $gl(2, \mathbb{R}) \oplus \mathbb{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 = (-\Gamma y, \gamma(1 - z))$

* $G = (-z, y)$

Length 2.

* $[F, G] = (\gamma - \delta z, -\delta y)$

Length 3.

* $[[F, G], F] = (-\gamma(\gamma - 2\Gamma) + \delta^2 z, -\delta^2 y)$

* $[[F, G], G] = (-2\delta y, -\gamma + 2\delta z)$

Length 4.

* $[[[F, G], F], F] = (-\delta^3 z + \gamma(\gamma^2 - 3\gamma\Gamma + 3\Gamma^2), -\delta^3 y)$

* $[[[F, G], F], G] = [[[F, G], G], F] = (0, \gamma(\gamma - 2\Gamma))$

* $[[[F, G], G], G] = (-\gamma + 4\delta z, 4\delta y)$

Couples of spins.

* Contrast, multisaturation: pairs with $(\Gamma_1, \gamma_1), (\Gamma_2, \gamma_2)$

* Multisaturation with $B_1$ inhomogeneity: in this case $\Gamma_1 = \Gamma_2 = \Gamma$, $\gamma_1 = \gamma_2$ and for the second spin the vector field $G$ is rescaled into $(1 - \varepsilon)G$.
# 3. ALGEBRAIC COMPUTATIONS

## 3.1. Preliminaries

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

* In vitro four test cases [5] 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 [45]:

- $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 $\Pi_1, \Pi_2$ the respective projections on the first and second spind: $q \mapsto q_1$, $q \mapsto 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

$$ \Gamma_i y_i^2 = \gamma_i(1 - z_i)z_i, 0 \le z_i \le 1 $$

intersected with either of the sets

$$ \Gamma_1 y_1 z_2 = \Gamma_2 y_2 z_1 \\ \gamma_1 y_1 (z_2 - 1) = \gamma_2 y_2 (z_1 - 1) $$

for the contrast setting, or with either of the sets

$$ \Gamma(1 - \varepsilon)y_1z_2 = \Gamma y_2z_1 \\ \gamma y_1(z_2 - 1) = \gamma(1 - \varepsilon)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\delta z = \gamma (\text{if } \delta \neq 0, \text{ 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$, $\gamma_1 z_2 - \gamma_2 z_1 = (\delta_1 - \delta_2) z_1 z_2$;

- the 2 lines $y_i = 0$, $z_1 = \gamma_i/(\delta_1 + \delta_i)$, $z_2 = \gamma_i/(\delta_2 + \delta_i)$ ($i = 1, 2$);

- the plane $z_1 = \gamma_1/2\delta$, $z_2 = \gamma_2/2\delta$ if $\delta_1 = \delta_2 = \delta$.

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 = \gamma/2\delta$.

## 3.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 [22, 23, 29, 50] for more information about those tools and others.

### 3.2.1. Operations on polynomial ideals

Throughout the rest of this section, we let $n$ be a positive integer, and we consider polynomials in $A = \mathbb{C}[X_1, \dots, 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 $\langle F \rangle$; the set of common complex zeroes of $F$ is denoted by $V(F)$ or $V(\langle F \rangle)$, 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 \cap \mathbb{C}[X_{k+1}, \dots, 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 $\mathbb{C}^{n-k}$ with coordinates $X_{k+1}, \dots, 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^{\infty}) := \{g \in A \mid \exists m \in \mathbb{N}, gf^n \in I\}. $$

Geometrically, the zeroes of this ideal is the smallest algebraic set containing $V(I) \setminus V(f)$.

This can be computed by adding the polynomial $u \cdot 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 [1, 17], resultants [14–16, 40] or Gröbner bases [13, 25, 26]. 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) \cap V(G)$ is the algebraic set $V(F \cup G)$.

The radical of an ideal $I$ is the set

$$ \sqrt{I} := \{f \in A \mid \exists m \in \mathbb{N}, f^m \in 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 \supset I$ such that $\sqrt{I} = \sqrt{J}$. We refer to algorithms in [35] for computing Gröbner bases of the radical of an ideal given by a finite set of its generators.

### 3.2.2. 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^{\alpha_1} \cdots X_n^{\alpha_n}$ compatible with multiplication. In the sequel, we will use mostly 3 orderings:

* the degree reverse-lexicographic ordering (DRL), denoted $X_1 > X_2 > \cdots > X_n$;

* the $k$'th elimination ordering (elim-$k$), denoted $X_1 > \cdots > X_k \gg X_{k+1} > \cdots > X_n$, and defined as a block-DRL ordering;

* the lexicographical ordering lex, denoted $X_1 \gg X_2 \gg \cdots \gg 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 \cap \mathbb{C}[X_{k+1}, \dots, X_n]$ is a (DRL) Gröbner basis of the ideal $I \cap \mathbb{C}[X_{k+1}, \dots, 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 [28].

### 3.2.3. 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 \subset \mathbb{R}^n$ defined as a finite union $S = \cup_{i=1}^r S_i$ of sets defined by polynomial equations $P_i$ and inequalities $Q_i$ with real coefficients:

$$S_i = \{ \mathbf{x} \in \mathbb{R}^n \mid \begin{array}{l} \forall P \in \mathcal{P}_i, P(\mathbf{x}) = 0 \\ \forall Q \in \mathcal{Q}_i, Q(\mathbf{x}) > 0 \end{array} \}$$

and maps between those sets.

Given two semi-algebraic sets A and B, a map $A \to 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) [18,19,42] of $\mathbb{R}^n$ is a sequence $S_1, \dots, S_n$ such that for all $i \in \{1, \dots, n\}$, $S_i$ is a decomposition of $\mathbb{R}^i$ into connected semi-algebraic subsets, called cells, and such that:

* Cells of $S_1$ are all points or intervals;

* For any $i \in \{1, \dots, n\}$, and for any cell $S_i$ of $S_i$, there exists a finite number of continuous semi-algebraic functions

$$\xi_{S_1, \dots, S_{l_S}} : S \to \mathbb{R}$$

such that $S \times \mathbb{R} \subset \mathbb{R}^{i+1}$ is a disjoint union of cells of $S_{i+1}$, which are

- either the graph $\Gamma_{S,j}$ of a function $\xi_{S,j}$

- or a band $B_{S,j}$ bounded by $\xi_{S,j}$ and $\xi_{S,j+1}$ for $j \in \{0, \dots, l_S\}$, with the convention that $\xi_{S,0} = -\infty$ and $\xi_{S,l_S+1} = +\infty$:

$$B_{S,j} = \{(x', x_{i+1}) \in S \times \mathbb{R} \mid \xi_{S,j}(x') < x_{i+1} < \xi_{S,j+1}(x')\}.$$

Given a set of polynomials $\mathcal{P} \subset \mathbb{R}[X_1, \dots, X_n]$, a CAD $S = (S_1, \dots, S_n)$ is $\mathcal{P}$-invariant if, for any cell $S \in S_n$ and any polynomial $P \in \mathcal{P}$, P has constant sign on S.

In particular, given a semi-algebraic set defined by polynomial equations $\mathcal{P}$ and polynomial inequalities $\mathcal{Q}$, a $\mathcal{P} \cup \mathcal{Q}$-invariant CAD gives an extensive topological description of the semi-algebraic set by means of a partition into subsets homeomorphic to $]0, 1]^l$ for some $i \in \{0, \dots, n\}$.

Unfortunately, computing a $\mathcal{P}$-invariant cylindrical algebraic decomposition, for $\mathcal{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 \le 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 \ne 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 [11], the command CylindricalDecomposition in the computer algebra system Mathematica and the subpackage SemialgebraicsetTools of RegularChains in the computer algebra system Maple.

### 3.2.4. 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 \le i \le 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 [31] 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 [24, 30, 47]. These are implemented and available in the RAGlib (Real Algebraic Geometry library) Maple package [46]. In the sequel, we roughly describe these algorithms.

The input consists of a polynomial system of equations and inequalities

$$f_1 = \dots = f_p = 0, \quad h_1 > 0, \dots, h_r > 0.$$

defining a semi-algebraic set $S \subset \mathbb{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 = \dots = f_p = 0$, $h_1 \ge e, \dots, h_r \ge e$. Finding a suitable $e_0$ can be done using notions of critical values and asymptotic critical values (see e.g. [30] 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 \subset \mathbb{R}^n$ be a semi-algebraic set defined and $e_0$ be as above $C \subset S$ be a connected component of $S$. Then, for $0 < e < e_0$, there exists $\{i_1, \dots, i_\ell\} \subset \{1, \dots, r\}$ and a connected component $C'$ of the real algebraic set defined by

$$f_1 = \dots = f_p = 0, \quad h_{i_1} = \dots = h_{i_\ell} = e$$

such that $C' \subset 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 = \dots = 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 [24] is the following. Let $a = (a_1, \dots, a_n) \subset \mathbb{Q}^n$ be generically chosen and $\pi_a : (x_1, \dots, x_n) \to a_1x_1 + \dots + a_nx_n$. Then for any connected component $C'$ of the real algebraic set under study, either $C'$ contains a critical point of the map $\pi_a$ restricted to the set under study (then such points can be computed using e.g. Gröbner bases, see [27]), or for an arbitrary $\rho \in \mathbb{R}$, $\pi_a^{-1}(\rho)$ meets $C'$ (hence it is sufficient to perform a recursive call to the algorithm by adding the constraint $a_1X_1 + \dots + a_nX_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.
### 3.2.5. 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 $\mathbb{C}^{n+t}$, with coordinates $(\mathbf{X}, \mathbf{G}) = (X_1, \dots, X_n, G_1, \dots, G_t)$ and the projection onto the parameter space

$$ \pi : \begin{array}{ccc} \mathbb{C}^{n+t} & \rightarrow & \mathbb{C}^t \\ (\mathbf{x}, \mathbf{g}) & \mapsto & \mathbf{g} \end{array} $$

Assume that $V$ has dimension $t$. Let $B \in \mathbb{R}^{n+t}$ be a semi-algebraic set defined only by a set of polynomial inequalities $H \subset \mathbb{R}[\mathbf{X}, \mathbf{G}]$:

$$ B = \{ (\mathbf{x}, \mathbf{g}) \in \mathbb{R}^{n+t} \mid \forall h \in H, h(\mathbf{x}, \mathbf{g}) > 0 \}, $$

we consider the semi-algebraic set $V \cap B \subset \mathbb{R}^{n+t}$.

The classification problem is to identify a dense covering of the real parameter space $\mathbb{R}^t$ into open subsets $U_1, \dots, U_N$ such that on each $U_i$, the number of points in the fibers of $\pi$ restricted to $V \cap B$ is constant:

$$ \exists c_i, \forall \mathbf{g} \in U_i, \#\pi^{-1}(\mathbf{g}) \cap V \cap B = c_i. $$

This problem has been thoroughly studied, see [39, 51] 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 [10], is adapted from the more general strategy of described in [39, 51] 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, \dots, p_M\} \subset \mathbb{R}[\mathbf{G}]$ such that the dense covering $(U_i)$ given by the connected components of $\mathbb{R}^t \setminus (V(p_1) \cup V(p_2) \cup \dots \cup 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 $\pi$ restricted to $V$;

* the projection of points where $V$ meets the topological border $\partial B$ of the semi-algebraic set $B$;

* points $\mathbf{g}$ at which the fiber $\pi^{-1}(\mathbf{g}) \cap V$ has positive real dimension;

* points at which the projection $\pi$ restricted to $V \cap 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 := \bigcup_{h \in H} V(h). $$

Note in particular that $\partial B \subset B_0$. We denote by $\mathrm{sing}(V)$ the singular locus of $V$, $\mathrm{crit}(\pi, V)$ the critical locus of $\pi$ restricted to $V$, and $K(\pi, V) = \mathrm{sing}(V) \cup \mathrm{crit}(\pi, V)$.

The hypotheses are:

* $\mathcal{H}1$ There exists a nonempty Zariski-open subset $O_1$ of $\mathbb{C}^t$ such that for all $\mathbf{g} \in O_1$, the fiber $V \cap \pi^{-1}(\mathbf{g})$ is a nonempty finite subset of $\mathbb{C}^{n+t}$

* $\mathcal{H}2$ The restriction of the projection $\pi$ to $B$ is proper;

* $\mathcal{H}3$ The intersection $V \cap B_0$ has dimension at most $t-1$ in $\mathbb{C}^{nt}$;

* $\mathcal{H}4$ The algebraic set $V$ is equidimensional with dimension $t$.

**Lemma 3.2.** Let $V$ and $B$ satisfying hypotheses $\mathcal{H}1$, $\mathcal{H}2$, $\mathcal{H}3$ and $\mathcal{H}4$. Let $U$ be a nonempty connected open subset of $\mathbb{R}^t$ which does not meet $\pi(V \cap B_0) \cup \pi(K(\pi, V))$, and let $\mathbf{g} \in U$. Then $V \cap \pi^{-1}(\mathbf{g})$ is finite, and for any $\mathbf{g}' \in U$, $\#V \cap \pi^{-1}(\mathbf{g}') = \#V \cap \pi^{-1}(\mathbf{g})$.
Hence, under these hypotheses, it suffices to compute a nonzero polynomial whose zeroes cover $\pi(V \cap B_0)$ and $\pi(K(\pi, 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 \times k$ matrix with coefficients in $\mathbb{R}[\mathbf{X}, \mathbf{G}]$, $r \in \{0, \dots, k-1\}$, the $r$'th determinantal variety associated with $M$ is

$$V_r(M) = \{(\mathbf{x}, \mathbf{g}) \in \mathbb{C}^{n+t} \mid \mathrm{rank}(M(\mathbf{x}, \mathbf{g})) \le r\}.$$

These varieties can be defined algebraically in two ways. First, one can use the minors modelization:

$$V_r(M) = V((r+1)\text{-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(\mathbf{x}, \mathbf{g}) \cdot Y(\mathbf{y}) = 0$$

where $Y(\mathbf{y})$ is a full-rank $k \times k - r$ matrix whose entries are new variables $\mathbf{y} = y_{1,1}, \dots, 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(\pi, V_{r_0}(M))$ for $r_0 \in \{0, \dots, k-1\}$. Assume that the variety $V$ satisfies the following three properties:

H5 There exists a nonempty Zariski-open subset $O_2 \subset \mathbb{C}^t$ such that

$$V \cap \pi^{-1}(O_2) = V_{r_0-1} \cap \pi^{-1}(O_2);$$

H6 For any $r_1 < r_2 \in \{0, \dots, 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 \in \{0, \dots, k-1\}$, the variety $V_r$ is equidimensional with dimension $n+t-(k-r+1)^2$.

These properties are generic [12], [44]. 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 \setminus V_{r_0-1}$, which has dimension less than $t$;

(2) the projection of $K(\pi, V_{r_0-1})$;

(3) the projection of $B_0 \cap V_{r_0-1}$.

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

$$P^t \cdot M(\mathbf{x}, \mathbf{g}) = 0$$

where $P$ is a vector of dimension 4 and $M$ is a matrix of size $4 \times 3$. This incidence variety encodes an undefined linear system, so its projection onto the affine space $\mathbb{C}^{n+t}$ with coordinates $(\mathbf{X}, \mathbf{G})$ is the affine space itself. In this situation, generically, $\mathrm{sing}(V)$ is exactly $V_2(M)$.

### 3.3. Algebraic computations for multisaturation with $B_1$-inhomogeneity

#### 3.3.1. Conventions and notations

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

$$z_1 \leftarrow z_1 + 1, z_2 \leftarrow 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 \begin{bmatrix} -\Gamma y & -z-1 & \delta z - \Gamma & 2\delta y \\ -\gamma z & y & \delta y & -2\delta z + \Gamma - \delta \\ -\Gamma y & (1-\varepsilon)(-z-1) & (1-\varepsilon)(\delta z - \Gamma) & (1-\varepsilon)^2(2\delta y) \\ -\gamma z & (1-\varepsilon)y & (1-\varepsilon)\delta y & (1-\varepsilon)^2(-2\delta z + \Gamma - \delta) \end{bmatrix}$$

and

$$D' = \det(F, G, [G, F], [[G, F], F]) = \det \begin{bmatrix} -\Gamma y & -z-1 & \delta z - \Gamma & \gamma(\gamma - 2\Gamma) + \delta^2(z+1) \\ -\gamma z & y & \delta y & \delta^2 y \\ -\Gamma y & (1-\varepsilon)(-z-1) & (1-\varepsilon)(\delta z - \Gamma) & (1-\varepsilon)^2(\gamma(\gamma - 2\Gamma) + \delta^2(z+1)) \\ -\gamma z & (1-\varepsilon)y & (1-\varepsilon)\delta y & (1-\varepsilon)^2\delta^2 y \end{bmatrix}$$

The following polynomials will appear frequently in the remainder of the section:

* $P_{y_1} := y_1 - (1 - \varepsilon)y_2$

* $P_{y_2} := y_2 - (1 - \varepsilon)y_1$

* $P_{z_1} := 2(\Gamma - \gamma)z_1 + 2\Gamma - \gamma$

* $P_{z_2} := 2(\Gamma - \gamma)z_2 + 2\Gamma - \gamma$

The root of the univariate polynomials $P_{z_1}$ and $P_{z_2}$ is

$$z_s = \frac{\gamma - 2\Gamma}{2\Gamma - 2\gamma}.$$

### 3.3.2. Transfer time not fixed ($H_F = 0$)

Singularities of {$D = 0$}.

**Proposition 3.4.** The set of points satisfying $D = \frac{\partial D}{\partial y_1} = \frac{\partial D}{\partial z_1} = \frac{\partial D}{\partial y_2} = \frac{\partial D}{\partial z_2} = 0$ is given, generically on authorized values of $\gamma, \Gamma$, 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

$$\begin{cases} y_1 = (1 - \varepsilon)y_2 \\ z_1 = z_2 = z_s = \frac{\gamma - 2\Gamma}{2\Gamma - 2\gamma}. \end{cases}$$

If $\gamma = \Gamma$ (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 - \varepsilon)\tilde{D}$. The singularities of $D$ and those of $\tilde{D}$ are the same, so for the study, we consider the ideal

$$I := \left\langle \tilde{D}, \frac{\partial \tilde{D}}{\partial y_1}, \frac{\partial \tilde{D}}{\partial y_2}, \frac{\partial \tilde{D}}{\partial z_1}, \frac{\partial \tilde{D}}{\partial z_2} \right\rangle.$$

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 \gg z_2 > \varepsilon > \Gamma > \gamma$. This computation yields that

$$I \cap \mathbb{Q}[z_2, \varepsilon, \Gamma, \gamma] = \langle \varepsilon^2(\varepsilon - 2)^2(2\Gamma - \gamma)(\Gamma - \gamma)z_2^3 P_{z_2}^3 \rangle,$$
<figure><img src="image_7.png" /><figcaption>FIGURE 6. Structure of the study of the singularities of {$D = 0$}</figcaption></figure>

so singular points necessarily satisfy

$$ egin{cases} z_2 = 0 \\ \text{or} \\ P_{z_2} = 0 \end{cases} \iff z_2 = \frac{\gamma - 2\Gamma}{2\Gamma - 2\gamma} $$

If $\Gamma = \gamma$ (that is, if the matter is water), the second of these solutions does not exist. If $\gamma = 2\Gamma$ (which means that the matter is on the limit of the domain of validity $2\Gamma \ge \gamma$), 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 := \langle \mathrm{sqfr}(G), z_2 \rangle \\ I_2 := \langle \mathrm{sqfr}(G), P_{z_2} \rangle $$

where for any polynomial $f$, $\mathrm{sqfr}(f)$ is the square-free part of $f$ and $\mathrm{sqfr}(G)$ means that we apply $\mathrm{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 \gg z_1 > z_2 > \varepsilon > \Gamma > \gamma$, and we find that this ideal contains

$$ \gamma z_1^2 (\varepsilon - 1)^2 (2\Gamma - \gamma), $$

so $z_1 = 0$.

We then compute a Gröbner basis of $\langle \mathrm{sqfr}(G_1), z_1 \rangle$ with respect to the order $y_1 \gg y_2 > z_1 > z_2 > \varepsilon > \Gamma > \gamma$, and we find that this ideal contains

$$ \Gamma \gamma \varepsilon y_2^2 (\Gamma - \gamma) (2\Gamma - \gamma)^2 (\varepsilon - 2), $$

so $y_2 = 1$.

Finally, adding $y_2$ to the ideal yields that

$$ 0 = \Gamma y_1 (\varepsilon - 1) (2\Gamma - \gamma), $$
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 = (\gamma - 2\Gamma)/(2\Gamma - 2\gamma)$. We compute a Gröbner basis $G_2$ of $I_2$ with respect to the order $y_1 > y_2 \gg z_1 > z_2 > \varepsilon > \Gamma > \gamma$, and we find that the ideal contains

$$ z_2\gamma(z_1 - z_2)^2. $$

Since this case was already studied, we may assume that $z_2 \neq 0$, so

$$ z_1 = z_2 = \frac{\gamma - 2\Gamma}{2\Gamma - 2\gamma}. $$

Adding $z_1 - z_2$ to $\text{sqfr}(G_2)$ and computing a Gröbner basis for the order $y_1 \gg y_2 > z_1 > z_2 > \varepsilon > \Gamma > \gamma$, we find that the ideal contains

$$ \gamma^2 y_2(\varepsilon - 1)P_{y_1}(2\Gamma - \gamma), $$

so we have 2 new branches to consider.

If $y_2 \neq 0$, $y_1 = (1 - \varepsilon)y_2$. Otherwise, by adding $y_2 = 0$ to the system of equations, we find that the ideal contains

$$ \gamma y_1^2 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 = \frac{\gamma - 2\Gamma}{2\Gamma - 2\gamma} \qquad (17) $$

(2) the line

$$ \left\{ \begin{array}{l} y_1 = y_2 = 0 \\ z_1 = z_2 \end{array} \right. \qquad (18) $$

(3) the surface (parameterized by $y_1, y_2$)

$$ z_1 = z_2 = \frac{\Gamma P_{y_2}^2 (\gamma - 2\Gamma)}{2(\Gamma - \gamma)a_3} \qquad (19) $$

with

$$ a_3 = (\Gamma + \gamma)P_{y_1}^2 + \varepsilon(\varepsilon - 2)\Gamma(y_1 - y_2)(y_1 + y_2) $$

(4) the surface (parameterized by $y_1, z_2$)

$$ y_2 = \frac{y_1 z_2}{(1 - \varepsilon)z_1}, \quad z_1 = \frac{(2\Gamma - \gamma)z_2}{a_4} \qquad (20) $$

with

$$ a_4 = 2(\varepsilon - 2)(\Gamma - \gamma)\varepsilon z_2 + (2\Gamma - \gamma)(\varepsilon - 1)^2 $$
<figure><img src="image_8.png" /><figcaption>FIGURE 7. Structure of the study of {$D = D' = 0$}</figcaption></figure>

(5) the surface (parameterized by $y_2, z_2$)

$$z_1 = \frac{z_2 y_1}{(1 - \varepsilon)y_2}, \quad y_1 = \frac{(1 - \varepsilon)y_2 ((2\Gamma - \gamma)\Gamma y_2^2 + \gamma^2 z_2^2)}{a_5} \quad (21)$$

with

$$a_5 = \Gamma(2\varepsilon(\varepsilon - 2)(\Gamma - \gamma)z_2 + (\varepsilon - 1)^2(2\Gamma - \gamma))y_2^2 + \gamma^2 z_2^2$$

*Proof.* The structure of this proof is summarized in Fig. 7. The determinant $D'$ factors as

$$D' = 2\gamma^2(2\Gamma - \gamma)(\Gamma - \gamma)(z_1 - z_2)(\varepsilon - 1)((\varepsilon - 1)y_2z_1 + y_1z_2),$$

so we form the two ideals

$$\begin{aligned} I_1 &= \langle \tilde{D}, z_1 - z_2 \rangle \\ I_2 &= \langle \tilde{D}, (\varepsilon - 1)y_2 z_1 + y_1 z_2 \rangle. \end{aligned}$$

If $z_1 = z_2$, after substitution, $\tilde{D}$ has two factors depending on $y_1, y_2, z_2$: $P_{z_2}$ and

$$p_2 := 2(\Gamma - \gamma)((\Gamma + \gamma)P_{y_1}^2 + \varepsilon(\varepsilon - 2)\Gamma(y_1 - y_2)(y_1 + y_2))z_2 + \Gamma P_{y_2}^2(2\Gamma - \gamma)$$

The polynomial $P_{z_2}$ gives solution 1.

Let

$$a_3(y_1, y_2) = (\Gamma + \gamma)P_{y_1}^2 + \varepsilon(\varepsilon - 2)\Gamma(y_1 - y_2)(y_1 + y_2)$$

so that the coefficient of $z_2$ in $p_2$ is $2(\Gamma - \gamma)a_3$, it is homogeneous in $y_1, y_2$ with degree 2. Its discriminant in $y_2$ is

$$-4(\varepsilon - 2)^2 \varepsilon^2 y_1^2 \gamma \Gamma.$$

Since the parameters $\gamma, \Gamma$ 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 \neq 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 = \frac{\Gamma P_{y_2}^2 (\gamma - 2\Gamma)}{2(\Gamma - \gamma)a_3(y_1, y_2)}. $$

We now turn to the other branch, defined by $(\varepsilon - 1)y_2z_1 + y_1z_2 = 0$.

If $y_1 = z_1 = 0$, there are 2 curves of singular points defined (in $y_2, z_2$) by

$$ \Gamma(2\Gamma - \gamma)y_2^2 + \gamma^2 z_2^2 = 0. $$

Since $2\Gamma \ge \gamma$, the only solution is $z_2 = 0$ with either $2\Gamma = \gamma$ or $y_2 = 0$.

If $y_1 = 0$ and $z_1 \neq 0$, then (since $\varepsilon \neq 1$) we must have $y_2 = 0$. Furthermore, we may assume that $z_1 \neq z_2$ since this case was already studied. The remaining solutions form a curve defined by

$$ 0 = p_3 := (2(\varepsilon - 2)(\Gamma - \gamma)\varepsilon z_2 + (2\Gamma - \gamma)(\varepsilon - 1)^2)z_1 - (2\Gamma - \gamma)z_2. $$

Let $a_4(z_1, z_2)$ be the coefficient of $z_1$ in $p_3$, the solutions are given by either

$$ \begin{cases} a_4(z_2) \neq 0 \\ z_1 = \frac{2\Gamma - \gamma}{c_3(z_2)} z_2 \end{cases} $$

or (since by assumption $z_1 \neq 0$)

$$ z_2 = 2\Gamma - \gamma = 0. $$

So we may assume that $y_1 \neq 0$. We compute a Gröbner basis of $I_2 + \langle uy_1 - 1 \rangle$ for the order $u \gg z_1 > y_1 \gg z_2 > y_2 > \varepsilon > \gamma > \Gamma$. This basis contains a polynomial which factors as the product of

$$ p_4 = (2(\varepsilon - 2)(\Gamma - \gamma)\varepsilon z_2 + (\varepsilon - 1)^2(2\Gamma - \gamma))y_1 + (\varepsilon - 1)(2\Gamma - \gamma)y_2 $$

and

$$ p_5 = (\Gamma(2\varepsilon(\varepsilon - 2)(\Gamma - \gamma)z_2 + (\varepsilon - 1)^2(2\Gamma - \gamma))y_2^2 + \gamma^2 z_2^2)y_1 + (\varepsilon - 1)y_2((2\Gamma - \gamma)\Gamma y_2^2 + \gamma^2 z_2^2). $$

First, assume that $p_4 = 0$. We compute a Gröbner basis of $I_2 + \langle uy_1 - 1, p_4 \rangle$ for the order $u \gg z_1 > y_1 \gg z_2 > y_2 > \varepsilon > \gamma > \Gamma$, and we find that the last polynomial defining the ideal is $p_3$, whose solutions we already studied.

Finally, assume that $p_4 \neq 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)\Gamma\gamma^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

$$ \begin{cases} a_5(y_2, z_2) \neq 0 \\ y_1 = \frac{(1-\varepsilon)y_2((2\Gamma-\gamma)\Gamma y_2^2 + \gamma^2 z_2^2)}{\Gamma(2\varepsilon(\varepsilon-2)(\Gamma-\gamma)z_2 + (\varepsilon-1)^2(2\Gamma-\gamma))y_2^2 + \gamma^2 z_2^2)} \end{cases} $$

and

$$ \begin{cases} a_5(y_2, z_2) = 0 \\ y_2((2\Gamma - \gamma)\Gamma y_2^2 + \gamma^2 z_2^2) \end{cases} $$

which, as in the case $y_1 = z_1 = 0$, is only $y_2 = z_2 = 0$ if $2\Gamma > \gamma$. This partial solution completes into $y_1 = y_2 = z_1 = z_2 = 0$, which was already known. $\square$
Equilibrium positions.

**Lemma 3.6.** The equilibrium points of $\dot{X} = DF - D'G$ are all contained in $\{D = D' = 0\}$.

*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. $\square$

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

$$ \frac{d}{dt}(q + \delta q) = (DF - D'G)(q) + A(q) \cdot \delta q + R(\delta q). $$

where $A = \text{Jac}_q(DF - D'G)$, so that

$$ \frac{d}{dt}(\delta q) = A(q) \cdot \delta q + R(q)(\delta q). $$

We can compute $A(q)$ explicitly: Indeed, let $f = DF - D'G$. Its first derivative is

$$ df(q)(u) = dD(q)(u)F(q) + D(q)dF(q)(u) - D'(q)(u)G(q) - D'(q)dG(q)(u), \quad (22) $$

so

$$ A(q) = \nabla D(q) \cdot F(q) + D(q) \cdot \text{Jac}_q(F)(q) - \nabla D'(q) \cdot G(q) - D'(q) \cdot \text{Jac}_q(G)(q). $$

We examine the eigenvalue decomposition of $A(q)$.

Solution 1 (17). If $z_1 = z_2 = \frac{\gamma - 2\Gamma}{2\Gamma - 2\gamma}$, the characteristic polynomial of $A$ factors as

$$ T^2 \left( T - \gamma^2 (2\Gamma - \gamma)^2 (\varepsilon - 1) P_{y_1}^2 \right)^2 $$

The matrix $A(q)$ is diagonalizable.

Solution 2 (18). 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} \begin{bmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix} P $$

with the transition matrix

$$ P = \begin{bmatrix} 0 & 1 & -1 & 0 \\ \varepsilon\gamma^3(\varepsilon - 1)(\varepsilon - 2)z_1^2 P_{z_1} & 1 & 0 & 0 \\ 0 & 1 & 0 & 1 \\ \varepsilon\gamma^3(\varepsilon - 1)(\varepsilon - 2)z_1^2 P_{z_1} & 0 & 0 & 0 \end{bmatrix}. $$

Solution 3 (19). If $z_1 = z_2 = \Gamma P_{y_2}^2 (\gamma - 2\Gamma)/2(\Gamma - \gamma)a_3$, the characteristic polynomial of $A(q)$ factors as

$$ T^2 \left( T + \frac{b_3}{a_3} \right) \left( T - \frac{b_3}{a_3} \right) $$
with

$$b_3 = \Gamma\gamma^2(\varepsilon - 1)P_{y_1}P_{y_2}(2\Gamma - \gamma)^2.$$

The matrix $A(q)$ is diagonalizable.

Solution 4 (20). If $y_2 = \frac{y_1z_2}{(1-\varepsilon)z_1}$ and $z_1 = \frac{(2\Gamma-\gamma)z_2}{a_4}$, the characteristic polynomial of $A(q)$ factors as

$$T^2 \left( T - \frac{b_4}{a_4} \right) \left( T + \frac{b_4}{a_4} \right)$$

with

$$b_4 = 2\varepsilon^2\gamma^3 z_2^3 (\varepsilon - 1)(\varepsilon - 2)^2 (2\Gamma - \gamma)(\Gamma - \gamma)P_{z_2}.$$

The matrix $A(q)$ is diagonalizable.

Solution 5 (21). If $z_1 = \frac{z_2y_1}{(1-\varepsilon)y_2}$ and $y_1 = \frac{(1-\varepsilon)y_2((2\Gamma-\gamma)\Gamma_{y_2}^2+\gamma^2z_2^2)}{a_5}$, the characteristic polynomial of $A(q)$ factors as

$$T^2 \left( T - \frac{b_4(\Gamma y_2^2 + \gamma(z_2^2 + z_2))\Gamma y_2^2}{a_5} \right)^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$, $\nabla D = \nabla 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:

$$\frac{d}{dt}(N + \delta q) = R(N)(\delta q) = O(\|\delta 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)$:

$$\frac{d}{dt}(q + \delta q) = (DF - D'G)(q) + A(q)(\delta q) + Q(q)(\delta q) + O(\|\delta q\|^3),$$

with $\frac{d}{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 [5, Sec. 3.4]. Differentiating (22) along $q$ again, the second derivative of $f$ is

$$\begin{aligned} \frac{d^2 f(q)(u, v)}{dt^2} = & \frac{d^2 D(q)(u, v)F(q)}{dt} + \frac{d^2 D(q)(u)}{dt} \frac{dF(q)}{dt} + \frac{d^2 D(q)(v)}{dt} \frac{dF(q)}{dt} \\ & - \frac{d^2 D'(q)(u, v)G(q)}{dt} - \frac{d^2 D'(q)(u)}{dt} \frac{dG(q)}{dt} - \frac{d^2 D'(q)(v)}{dt} \frac{dG(q)}{dt} \end{aligned} \quad (23)$$

Note that second derivatives of $F$ and $G$ are 0, since their coordinates are affine in $q$.

We wish to compute $H_2(\delta q) = Q(S)(\delta q, \delta q) = \frac{1}{2}d^2 f(q)(\delta q, \delta q)$. Since $dD(S) = dD'(S) = 0$, we find in the end that

$$H_2(\delta q) = h_2(\delta q)F(S) - h'_2(\delta q)G(S),$$
with

$$
\begin{align*}
F(S) &= \left( 0, \frac{\gamma(2\Gamma - \gamma)}{2(\Gamma - \gamma)}, 0, \frac{\gamma(2\Gamma - \gamma)}{2(\Gamma - \gamma)} \right)^t \\
G(S) &= \left( \frac{\gamma}{2(\Gamma - \gamma)}, 0, \frac{(1 - \varepsilon)\gamma}{2(\Gamma - \gamma)}, 0 \right)^t \\
h_2(\delta q) &= \frac{1}{2} d^2 D(S)(\delta q, \delta q) = (1 - \varepsilon)(\delta z_1 - \delta z_2)(\delta z_1 - (1 - \varepsilon)^2 \delta z_2)(2\Gamma - \gamma)\gamma^2 \\
h'_2(\delta q) &= \frac{1}{2} d^2 D'(S)(\delta q, \delta q) = (1 - \varepsilon)(\delta z_1 - \delta z_2)(\delta y_2(\varepsilon - 1) + \delta y_1)(2\Gamma - \gamma)^2 \gamma^2.
\end{align*}
$$

Following [5] and [41], we study the projection of the differential equation $\dot{v} = H_2(v)$ on the sphere $S^3$. Let $w = v/\|v\|$ be this projection, it satisfies the differential equation

$$
\begin{align*}
\dot{w} &= \frac{1}{\|v\|^2} \left( \dot{v} \|v\| - v \frac{\langle v, \dot{v} \rangle}{\|v\|} \right) \\
&= \frac{H_2(v)}{\|v\|} - \frac{\langle v, H_2(v) \rangle}{\|v\|^3} v \\
&= \|v\| (H_2(w) - \langle w, H_2(w) \rangle w)
\end{align*}
$$

so we are to study the following differential equation on the sphere $S^3$:

$$
\dot{v} = H_2(v) - \langle v, H_2(v) \rangle v =: H_2^\pi(v).
$$

Invariants are related to the eigenvalues of the linearization of $H_2^\pi$ at points where $H_2^\pi(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 $\xi$ such that there exists $\lambda \in \mathbb{R} \setminus \{0\}$, $H_2(\xi) = \lambda \xi$.

We study the linearization of $H_2^\pi$ 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 $\delta z_1 = \delta z_2$;

(2) the projective line $\delta y_2 = (1 - \varepsilon) \delta y_1$, $\delta z_1 = (1 - \varepsilon)^2 \delta z_2$.

In the first case, the Jacobian of the system is nilpotent. In the second case, it is diagonalizable with non-zero eigenvalues:

$$
\frac{1}{2} \left( \delta \bar{y}_2 + 1 \pm \sqrt{ \delta \bar{y}_2^2 + (2\varepsilon - 1)^2 (\delta \bar{y}_2 + 1) - 4(\varepsilon - 1)^4 - 2\delta \bar{y}_2 + 1 } \right).
$$

*Proof.* First we study ray solutions. Let $\xi$ be a vector on a ray solution, such that

$$
H_2(\xi) = \lambda \xi.
$$

Let $\alpha\xi$ be another vector on the same line ($\alpha \in \mathbb{R}$), since $H_2$ is homogeneous with degree 2, one has

$$
H_2(\alpha\xi) = \alpha^2 H_2(\xi) = \alpha^2 \lambda \xi = \alpha \lambda(\alpha\xi).
$$

So each line or ray solutions contains a unique $\xi_0$ such that $H_2(\xi_0) = \xi_0$.

A Gröbner basis of the system $\langle H_2(\delta q) - \delta q \rangle$ is given by $\{\delta y_1, \delta z_1, \delta y_2, \delta z_2\}$, so there is no non-trivial ray.

This can also be seen in the following way: let $\delta q$ be a vector such that $H_2(\delta q) = \delta q$. By the structure of the vector $F(S)$, $\delta q$ satisfies $\delta z_1 = \delta z_2$, and so $h_2(\delta q) = h'_2(\delta q) = 0$, so $H_2(\delta q) = 0$, and which, by hypothesis, implies that $\delta 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^\pi$ in some neighborhood of these solutions in $S^3$. First we consider vectors $\delta q$ such that $\delta z_1 = \delta z_2$, we may perform the computations in the affine chart given by $\delta z_1 \neq 0$, with coordinates $\delta \bar{y}_1 = \delta y_1/\delta z_1$, $\delta \bar{y}_2 = \delta y_2/\delta z_1$, $\delta \bar{z}_2 = \delta z_2/\delta z_1$. The differential equation becomes

$$ \frac{d}{dt} \begin{pmatrix} \delta \bar{y}_1 \\ \delta \bar{y}_2 \\ \delta \bar{z}_1 \end{pmatrix} = \delta z_2 \bar{C} (\delta \bar{y}_1, \delta \bar{y}_2, \delta \bar{z}_1) $$

with $\bar{C}$ a polynomial vector field of degree 3.

At $\delta \bar{z}_1 = 1$, its Jacobian is nilpotent:

$$ \begin{bmatrix} 0 & 0 & (\varepsilon - 1)(\delta \bar{y}_1(\varepsilon(\varepsilon - 1) - 1) + \delta \bar{y}_2(1 - \varepsilon))(2\Gamma - \gamma)^2 \gamma^3 \\ 0 & 0 & (\varepsilon - 1)(\delta \bar{y}_1(2\varepsilon(\varepsilon - 2) + 1) + \delta \bar{y}_2(\varepsilon - 1))(2\Gamma - \gamma)^2 \gamma^3 \\ 0 & 0 & 0 \end{bmatrix} $$

Then we consider vectors $\delta q$ such that $\delta y_2 = (1 - \varepsilon)\delta y_1$, $\delta z_1 = (1 - \varepsilon)^2\delta z_2$. This time we use the chart $\delta z_2 \neq 0$ with coordinates $\delta \bar{y}_1 = \delta y_1/\delta z_2$, $\delta \bar{y}_2 = \delta y_2/\delta z_2$, $\delta \bar{z}_1 = \delta z_1/\delta z_2$. As above, we compute the differential equation in this chart, linearize the resulting vector field, and evaluate this Jacobian at $\delta \bar{y}_2 = (1 - \varepsilon)\delta \bar{y}_1$ and $\delta \bar{z}_1 = (1 - \varepsilon)^2$. This matrix has rank 2 and is diagonalizable with non-zero eigenvalues:

$$ \frac{1}{2} \left( \delta \bar{y}_2 + 1 \pm \sqrt{\delta \bar{y}_2^2 + (2\varepsilon - 1)^2 (\delta \bar{y}_2 + 1) - 4(\varepsilon - 1)^4 - 2\delta \bar{y}_2 + 1} \right). $$

□

Cubic approximation at $N$. We perform the same study at the North pole $N$. With expression (23), we can verify that the quadratic component of $R(N)$ is 0. Indeed, $F(N) = 0$ and $d^2D'(N) = 0$.

Further differentiating along $q$, we obtain

$$ H_3(\delta q) := \frac{1}{6} d^3 f(N)(\delta q, \delta q, \delta q) = \frac{1}{6} (3d^2D(N)(\delta q, \delta q)F(\delta q) - d^3D'(N)(\delta q, \delta 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 $\dot{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

$$ \left\{ \begin{array}{l} \delta z_1 = \delta z_2 = 0 \\ (\varepsilon - 1)\delta y_1 + \delta y_2 = \frac{1}{\Gamma(2\Gamma - \gamma)\sqrt{1 - \varepsilon}}; \end{array} \right. \qquad (24) $$

(2) the quadric

$$ \left\{ \begin{array}{l} \delta y_1 = \delta y_2 = 0 \\ ((\varepsilon - 1)\delta z_1 - \delta z_2)^2 - (\varepsilon - 2)\delta z_1 \delta z_2 = \frac{1}{\gamma^3(2\Gamma - \gamma)(1 - \varepsilon)}. \end{array} \right. $$

and three sets of real non-isolated singularities:
$$
(1) \textit{the plane} \qquad \left\{ \begin{array}{l} \delta z_1 = \delta z_2 \\ \delta y_1(1-\epsilon) = \delta y_2 \end{array} \right. \tag{25}
$$

$$
(2) \textit{the plane} \qquad \left\{ \begin{array}{l} \delta y_2 = (1-\epsilon)\delta y_1 \\ \delta z_2 = (1-\epsilon)^2 \delta z_1 \end{array} \right. \tag{26}
$$

(3) the surface defined by

$$
\left\{
\begin{array}{l}
0 = \Gamma(\epsilon - 1)(2\Gamma - \gamma)\delta y_1^2 + (2\Gamma - \gamma)\Gamma\delta y_1\delta y_2 + \gamma^2(\epsilon - 1)\delta z_1^2 - \gamma^2(\epsilon - 1)\delta z_1\delta z_2 \\
0 = \Gamma(\epsilon - 1)(2\Gamma - \gamma)\delta y_1\delta y_2 - \gamma^2\delta z_1\delta z_2 + (2\Gamma - \gamma)\Gamma\delta y_2^2 + \gamma^2\delta z_2^2 \\
0 = \delta y_1\delta z_2 + (\epsilon - 1)\delta y_2\delta z_1.
\end{array}
\right.
\tag{27}
$$

For points on the line (24), the linearization of $H_3^\pi$ 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:

$$
\begin{bmatrix}
0 & 0 & 0 \\
0 & * & 1 \\
0 & 0 & *
\end{bmatrix}
$$

For isolated singularities satisfying (26), the matrix is diagonalizable with 3 non-zero eigenvalues.

*Proof.* First we study ray solutions. Let $\xi$ be a vector on a ray solution, such that

$$
H_3(\xi) = \lambda \xi.
$$

Let $\alpha\xi$ be another vector on the same line ($\alpha \in \mathbb{R}$), one has

$$
H_3(\alpha\xi) = \alpha^3 H_3(\xi) = \alpha^3 \lambda\xi = \alpha^2 \lambda(\alpha\xi).
$$

So unlike in the quadratic case, a line of ray solutions contains 2 vectors $\xi_1, \xi_2$ such that either $H_3(\xi_1) = H_3(\xi_2) = 1$ or
$H_3(\xi_1) = H_3(\xi_2) = -1$.

In order to study ray solutions, we compute a Gröbner basis of the system $H_3(\delta q) - \iota \delta q = 0, \iota^2 = 1$, for an order eliminating $\alpha$. We find that the basis contains

$$
\{ \delta y_1 \delta z_1, \delta y_1 \delta z_2, \delta y_2 \delta z_1, \delta y_2 \delta z_2 \}
$$

so either $\delta y_1 = \delta y_2 = 0$ or $\delta z_1 = \delta z_2 = 0$.

If $\delta z_1 = \delta z_2 = 0$, computing a new Gröbner basis of the system, saturating with $\Gamma - \gamma$ and $\delta y_1$ (see Sec. 3.2.1), shows
that $\delta y_1$ and $\delta y_2$ must satisfy

$$
\left( \Gamma^2 (\epsilon - 1) (2\Gamma - \gamma)^2 \left( (\epsilon - 1) \delta y_1 + \delta y_2 \right)^2 \right)^2 - 1 = 0.
$$

Since $\epsilon - 1 < 0$, $2\Gamma - \gamma > 0$ and $\Gamma > 0$, this defines 2 lines of real solutions given by

$$
(\epsilon - 1)\delta y_1 + \delta y_2 = \pm \frac{1}{\Gamma(2\Gamma - \gamma)\sqrt{1 - \epsilon}}
$$

Those lines are equivalent in the projective space: each line of ray-solutions contains a vector in both lines.
If $\delta y_1 = \delta y_2 = 0$, the same technique shows that $\delta z_1$ and $\delta z_2$ must satisfy

$$ (\gamma^3(2\Gamma - \gamma)(\varepsilon - 1)((\varepsilon - 1)\delta z_1 - \delta z_2)^2 - (\varepsilon - 2)\delta z_1 \delta z_2))^2 - 1 = 0, $$

which defines 1 quadric

$$ ((\varepsilon - 1)\delta z_1 - \delta z_2)^2 - (\varepsilon - 2)\delta z_1 \delta z_2 = \frac{1}{\gamma^3(2\Gamma - \gamma)(1 - \varepsilon)}. $$

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(\delta q) = 0$, saturating by $\gamma - \Gamma$, $\gamma$, $\Gamma$, $\varepsilon - 1$ and $2\Gamma - \gamma$. Factoring the results, it appears that the solutions split into 5 cases:

(1) $\delta z_1 = \delta z_2$

(2) $\delta z_1 = -\delta z_2$

(3) $\delta z_1 = 0$

(4) $\delta z_2 = (1 - \varepsilon)^2 \delta z_1$

(5) otherwise.

In the case 1, a Gröbner basis is given by

$$ \left\{ \begin{array}{l} \delta y_1 (\delta y_1 (\varepsilon - 1) + \delta y_2)^2 \\ \delta y_2 (\delta y_1 (\varepsilon - 1) + \delta y_2)^2 \\ \delta z_2 (\delta y_1 (\varepsilon - 1) + \delta y_2)^2 \\ \delta z_1 - \delta z_2 \end{array} \right\} $$

and the solutions form the plane (25).

In the case 2, saturating by $\delta z_1 - \delta z_2$, a Gröbner basis is given by

$$ \left\{ \begin{array}{l} \delta y_1 + (1 - \varepsilon)\delta y_2 \\ \delta z_1 + \delta z_2 \\ \Gamma(\varepsilon^2 - 2\varepsilon + 2)(2\gamma - \gamma)\delta y_2^2 + 2\gamma^2 \delta z_2^2 \end{array} \right\} $$

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

In the case 3, saturating by $\delta z_1 - \delta z_2$ and $\delta z_1 + \delta z_2$, a Gröbner basis is given by

$$ \{ \delta y_1, \delta z_1, (2\Gamma - \gamma)\Gamma\delta y_2^2 + \gamma^2 \delta z_2^2 \} $$

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

In the case 4, saturating by $\delta z_1 \pm \delta z_2$ and $\delta z_1$, a Gröbner basis is given by

$$ \{ \delta y_1 \delta z_2 + (\varepsilon - 1) \delta y_2 \delta z_1, \delta y_2 + (\varepsilon - 1) \delta y_1, \delta z_2 - (\varepsilon - 1)^2 \delta z_1 \} $$

and the solutions form the plane (26).

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

$$ \left\{ \begin{array}{l} \Gamma(\varepsilon - 1)(2\Gamma - \gamma)\delta y_1^2 + (2\Gamma - \gamma)\Gamma\delta y_1\delta y_2 + \gamma^2(\varepsilon - 1)\delta z_1^2 - \gamma^2(\varepsilon - 1)\delta z_1\delta z_2 \\ \Gamma(\varepsilon - 1)(2\Gamma - \gamma)\delta y_1\delta y_2 - \gamma^2\delta z_1\delta z_2 + (2\Gamma - \gamma)\Gamma\delta y_2^2 + \gamma^2\delta z_2^2 \\ \qquad \delta y_1\delta z_2 + (\varepsilon - 1)\delta y_2\delta z_1 \\ \qquad (\Gamma(2\Gamma - \gamma)(\varepsilon - 1)^2\delta y_2^2 + \delta z_2^2\gamma^2)\delta z_1 - (2\Gamma - \gamma)\Gamma\delta z_2\delta y_2 - \gamma^2\delta z_2^3 \end{array} \right\} $$

The fourth polynomial is a combination of the other 3, and the solutions form the surface (27).
For the second part of the proposition, as in the quadratic case, we study the linearization of $H_3^\pi : \dot{v} = H_3(v) - \langle v, H_3(v) \rangle v$. In the affine chart with $\delta y_1 \neq 0$, with coordinates $\bar{\delta z_1} = \delta z_1/\delta y_1$, $\bar{\delta y_2} = \delta y_2/\delta y_1$ and $\bar{\delta z_2} = \delta z_2/\delta y_1$, the differential equation $\dot{v} = H_3(v)$ becomes

$$ \frac{d}{dt} \begin{pmatrix} \bar{\delta z_1} \\ \bar{\delta y_2} \\ \bar{\delta z_2} \end{pmatrix} = \delta y_1^2 \bar{Q} (\bar{\delta z_1}, \bar{\delta y_2}, \bar{\delta z_2}) $$

with $\bar{Q}$ a polynomial vector field of degree 4. In this chart, $H_3^\pi$ becomes

$$ \dot{v} = \bar{Q}(v). $$

We conclude by evaluating the Jacobian of $\bar{Q}$ at the relevant points.

<figure><img src="image_9.png" /></figure>

### 3.3.3. General case

Singularities of $\{\mathcal{D} = H_G = \{H_G, H_F\} = 0\}$.

**Proposition 3.9.** The set of singularities of $\{\mathcal{D} = H_G = \{H_G, H_F\} = 0\}$ is given, generically on authorized values of $\gamma, \Gamma$, 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 $\gamma = \Gamma$, solution 2 becomes a surface defined by $z_1 = z_2$, $y_2 = (1 - \varepsilon)y_1$.

*Proof.* By definition of $\mathcal{D} = \{\{H_G, H_F\}, H_G\}$, we want to study the zeroes of

$$ \begin{cases} 0 = p \cdot G \\ 0 = p \cdot [G, F] \\ 0 = p \cdot [[G, F], G] \end{cases} \qquad (D) $$

The singularities of this variety is the set of points at which the matrix

$$ \begin{bmatrix} G & [G,F] & [[G,F],G] \end{bmatrix} $$

has rank at most 2 (as per Remark 3.3).

We encode that with the incidence variety

$$ \begin{bmatrix} G & [G,F] & [[G,F],G] \\ v_1 & v_2 & v_3 \end{bmatrix} \cdot \begin{bmatrix} L_1 \\ L_2 \\ L_3 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 1 \end{bmatrix} $$

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, \Gamma, \gamma, \varepsilon$. We eliminate $L_1, L_2, L_3$ from the ideal in order to recover the projection, and we saturate by $1 - \varepsilon$.

We then compute a Gröbner basis for the elimination order $y_1 > y_2 > z_1 > z_2 \gg \gamma > \Gamma > \varepsilon$, this basis contains 10 polynomials, some of which have factors with multiplicity greater than 1 or are divisible by $1 - \varepsilon$ or $\gamma$. We take the square-free form of this basis, and we saturate by $1 - \varepsilon$ and $\gamma$ 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) ((\varepsilon - 1)P_{z_2}y_1 + P_{z_1}y_2). $$
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 \gg \gamma > \Gamma > \varepsilon$, take its square-free form, and recompute a Gröbner basis. The result contains the polynomial

$$P_{z_2}P_{y_2}.$$

The solutions decompose into 2 algebraic sets, defined by

$$z_1 = z_2 = z_S = \frac{2\Gamma - \gamma}{2\Gamma - 2\gamma}$$

and (adding $(\varepsilon - 1)y_1 + y_2$ to the ideal and saturating by $P_{z_2}$, $\varepsilon$ and $\varepsilon - 2$)

$$\begin{cases} 0 = y_1(\Gamma - \gamma) \\ 0 = y_2(\Gamma - \gamma) \\ 0 = P_{y_2} = (\varepsilon - 1)y_1 + y_2 \\ 0 = z_1 - z_2. \end{cases}$$

The latter is generically (if $\Gamma \neq \gamma$) defined by

$$\begin{cases} y_1 = y_2 = 0 \\ z_1 = z_2 \end{cases}$$

or if $\Gamma = \gamma$ (that is if the spin we consider is water), by

$$\begin{cases} y_1 = \frac{y_2}{\varepsilon - 1} \\ z_1 = z_2. \end{cases}$$

Then, starting again with the whole ideal, we add $(\varepsilon - 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 $(\Gamma, \gamma, \varepsilon)$, is an irreducible surface. $\square$

Locus of $\{\mathcal{D} = \mathcal{D}' = H_G = \{H_G, H_F\} = 0\}$.

**Proposition 3.10.** *The solutions form the union of the hyperplane defined by*

$$z_1 = z_2 \qquad (28)$$

and the hypersurface

$$y_1 = - \frac{y_2 P_{z_1}}{(\varepsilon - 1) P_{z_2}}. \qquad (29)$$

*Proof.* Points such that $\{\mathcal{D} = \mathcal{D}' = H_G = \{H_G, H_F\} = 0\}$ satisfy

$$\begin{cases} 0 = p \cdot G \\ 0 = p \cdot [G, F] \\ 0 = p \cdot [[G, F], G] \\ 0 = p \cdot [[G, F], F] \end{cases}$$

The projection of these points onto the space $(y_1, z_1, y_2, z_2)$ is given by the vanishing of the determinant $\Delta'$, defined as

$$\begin{align*} \Delta' &= \det \begin{bmatrix} G & [G, F] & [[G, F], G] & [[G, F], F] \\ (1-\varepsilon)G & [(1-\varepsilon)G, F] & [[[(1-\varepsilon)G, F], (1-\varepsilon)G] & [[[(1-\varepsilon)G, F], F]] \end{bmatrix} \\ &= (1-\varepsilon)^2 \det \begin{bmatrix} G & [G, F] & [[G, F], G] & [[G, F], F] \\ G & [G, F] & (1-\varepsilon)[[G, F], G] & [[G, F], F] \end{bmatrix} \end{align*}$$
This determinant factors as

$$-2 (\varepsilon - 1)^2 (2\Gamma - 1) (\Gamma - 1) (z_1 - z_2) [(\varepsilon - 1) P_{z_2} y_1 + P_{z_1} y_2]$$

□

Equilibrium positions.

**Lemma 3.11.** All equilibrium points of $\dot{Z} = \mathcal{D}\vec{H}_F - \mathcal{D}'\vec{H}_G$ satisfying $H_G = \{H_G, H_F\} = 0$ are contained in $\{\mathcal{D} = \mathcal{D}' = 0\}$.

*Proof.* Recall that $\vec{H}_F$ is defined as

$$\vec{H}_F = \left[ \frac{\partial H_F}{\partial p} \right] = \left[ \frac{F}{-p \cdot \frac{\partial H_F}{\partial q}} \right]$$

and $\vec{H}_G$ is defined in the same way. Let $z = (q, p)$ be a point such that

$$\mathcal{D}(z)\vec{H}_F(z) - \mathcal{D}'(z)\vec{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

$$\bullet X_D F(q) - X_{D'} G(q) = 0$$

$$\bullet H_G(z) = \{H_G, H_F\}(z) = 0$$

$$\bullet \mathcal{D}(z)\vec{H}_F(z) - \mathcal{D}'(z)\vec{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 $\langle 1 \rangle$, and so the associated system has no solution.

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

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

$$\mathcal{A} = \mathrm{Jac}(\mathcal{D}\vec{H}_F - \mathcal{D}'\vec{H}_G)$$

on equilibrium point, given as the union of points satisfying Eq. (28) and (29).

Solutions of Eq. (28). If $z_1 = z_2$, the matrix $\mathcal{A}$ has rank 2, and its characteristic polynomial is

$$T^6 \left( T^2 - \frac{\varepsilon(\varepsilon - 1)(\varepsilon - 2)(8(\Gamma - \gamma)^2 P_{y_1} y_1 y_2 + P_{z_2}(2P_{z_2} - \gamma)P_{y_2})}{2P_{y_1}(\Gamma - \gamma)y_1} T + \left( \frac{(\varepsilon - 1)\gamma P_{y_1}(2\Gamma - \gamma)}{y_1} \right)^2 \right)$$

Solutions of Eq. (29). If $y_2 = -\frac{y_2 p_{z_1}}{(\varepsilon - 1) p_{z_2}}$, the matrix $\mathcal{A}$ has rank 2, and its characteristic polynomial is

$$T^6 \left( T^2 - 4\varepsilon y_2 (\varepsilon - 1)(\varepsilon - 2)(\Gamma - \gamma)T + \left( \frac{2\gamma(z_1 - z_2)(\varepsilon - 1)^2(2\Gamma - \gamma)(\Gamma - \gamma)}{p_{z_1}} \right)^2 \right).$$

The discriminant of the degree 2 factor factors as

$$\frac{16(\Gamma - \gamma)^2(\varepsilon - 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_{z_1}^2}$$

with

$$a_6(z_1, y_2) = \varepsilon(\varepsilon - 2)P_{z_1}y_2$$

$$b_6(z_1, z_2) = (2\Gamma - \gamma)(\varepsilon - 1)\gamma(z_1 - z_2)$$
which induces the following classification of the eigenvalues of $\mathcal{A}$:

* if $|a(z_1, y_2)| > |b(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)

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 [10] 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\Gamma_2 \ge \gamma_2 > 0$, and that $(\gamma_1, \Gamma_1) \ne (\gamma_2, \Gamma_2)$. We normalize by means of $\Gamma_1 = \gamma_1 = 1$ (water), leaving free the two parameters $\Gamma_2, \gamma_2$ corresponding to the second matter. Under this normalization, the hypothesis becomes $(\gamma_2, \Gamma_2) \ne (1, 1)$. Let $\Pi$ be the plane $\{y_1 = y_2 = 0\}$.

The following theorems were proved in [10] using a real roots classification algorithm.

**Theorem 3.12.** Consider the 9 polynomials:

$$ 
\begin{aligned}
f_1 &= \Gamma_2 - 1 \\
f_2 &= 3\Gamma_2 - 2\gamma_2 - 1 \\
f_3 &= 3\Gamma_2^2 - 5\Gamma_2\gamma_2 + \gamma_2^2 + 2\Gamma_2 - 2\gamma_2 + 1 \\
f_4 &= 2\Gamma_2^2 - 5\Gamma_2\gamma_2 + 2\gamma_2^2 - 2\Gamma_2 + 3\gamma_2 \\
f_5 &= 2\gamma_2^3 - (3\Gamma_2 + 11)\gamma_2^2 + (9\Gamma_2 + 6 - 3\Gamma_2^2)\gamma_2 + 2\Gamma_2(\Gamma_2 + 2)(\Gamma_2 - 1) \\
f_6 &= \Gamma_2 - 2\gamma_2 + 1 \\
f_7 &= 2\Gamma_2 - \gamma_2 - 1 \\
f_8 &= \gamma_2 - 2 + \Gamma_2 \\
f_9 &= 2\Gamma_2^2 - 5\Gamma_2\gamma_2 + 2\gamma_2^2 + 1.
\end{aligned}
$$

The zeroes of their product divide the subset of $\mathbb{R}^2$ defined by $2\Gamma_2 > \gamma_2 > 0$ into connected components where the cardinality of $V_{\mathbb{R}} \cap \pi^{-1}(\gamma_2, \Gamma_2)$ is constant.

Using a Cylindrical Algebraic Decomposition, one can make this statement more precise by expliciting this cardinality. Let $\psi: (y_1, z_1, y_2, z_2) \mapsto (-y_1, z_1, -y_2, z_2)$ be the symmetry fixing $\Pi = \{y_1 = y_2 = 0\}$, and let us consider the semi-algebraic sets (see Fig. 8):

$$ 
\begin{aligned}
\mathcal{G}_1^- &= \{0 < \gamma_2 < 2\Gamma_2, \Gamma_2 < 1, f_2 > 0\}, \\
\mathcal{G}_1^+ &= \{0 < \gamma_2 < 2\Gamma_2, \Gamma_2 > 1, f_2 < 0, f_4 > 0\}, \\
\mathcal{G}_2^- &= \{0 < \Gamma_2 < 1, f_6 > 0, f_3 < 0\}, \\
\mathcal{G}_2^+ &= \{\Gamma_2 > 1, f_6 < 0, f_5 > 0\}, \\
\mathcal{G}_1 &= \mathcal{G}_1^+ \cup \mathcal{G}_1^-, \\
\mathcal{G}_2 &= \mathcal{G}_2^+ \cup \mathcal{G}_2^-, \\
\mathcal{G} &= \mathcal{G}_1 \cup \mathcal{G}_2.
\end{aligned}
$$
<figure><img src="image_10.png" /><figcaption>FIGURE 8. Curves involved in the definition of the semi-algebraic set $\mathcal{G}$. The blue (resp. green) sample points correspond to points in $\mathcal{G}_1^- \cup \mathcal{G}_1^+$ (resp. $\mathcal{G}_2^- \cup \mathcal{G}_2^+$). The circled numbers in each area correspond to the number of singularities in $\mathcal{B}$ for parameters in the area. Parameters in the red area are physically irrelevant.</figcaption></figure>

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

(1) if $(\gamma_2, \Gamma_2) \in \mathcal{G}_1$ there is one other singularity lying on $\Pi \cap \mathcal{B}$;

(2) if $(\gamma_2, \Gamma_2) \in \mathcal{G}_2$ there are two other singularities in $\mathcal{B}$, $\psi$-symmetric, outside $\Pi$.

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 and 9. Observe that the number of singularities inside $\mathcal{B}$ is an invariant of the contrast problem. Two of the pairs of biological matters studied in [8], water-cerebrospinal fluid (normalized parameters $[\gamma_2 = \frac{5}{4}, \Gamma_2 = \frac{25}{3}]$) and water-fat (normalized parameters $[\gamma_2 = \frac{25}{2}, \Gamma_2 = 25]$) correspond to points outside $\mathcal{G}$, and their invariant is 1 in both cases (see Fig. 9).

We can give a more precise description of those points. Assuming that $(\gamma_2 - \Gamma_2)(2\Gamma_2 - \gamma_2 - 2) \neq 0$, we define a point (with the coordinates $(y_1, z_1, y_2, z_2)$):

$$ \Omega: \left(0, \frac{\Gamma_2 - 1}{\Gamma_2 - \gamma_2}, 0, \frac{(1 - \Gamma_2)(2\Gamma_2 - \gamma_2)}{(2\Gamma_2 - \gamma_2 - 2)(\Gamma_2 - \gamma_2)}\right) $$

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.
<figure><img src="image_11.png" /><figcaption>FIGURE 9. Positions of the parameters corresponding to the pairs water-cerebrospinal fluid (red circle) and water-fat (red square) and the set $\mathcal{G}$ (with the same conventions as in Fig. 8). For both these pairs, there is only 1 singularity in $\mathcal{B}$.</figcaption></figure>

Let

$$ \delta(\gamma_2, \Gamma_2) = (\Gamma_2 - 1)(2\Gamma_2 - \gamma_2)(\Gamma_2 - 2\gamma_2 + 1)(\Gamma_2 - 2\gamma_2 + 2) \quad (30) $$

and

$$ \eta(\gamma_2, \Gamma_2) = (5\Gamma_2 - 4\gamma_2 + 1)(\Gamma_2 - \gamma_2)(\Gamma_2 + \gamma_2 - 1)(\Gamma_2 - 2\gamma_2 + 2). $$

If $\eta(\gamma_2, \Gamma_2) \neq 0$ we denote by $\Omega_a^+$ (resp. $\Omega_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 = \frac{\sqrt{3}\sqrt{\delta(\gamma_2, \Gamma_2)}}{(5\Gamma_2 - 4\gamma_2 + 1)(\Gamma_2 - \gamma_2)} $$

$$ z_1 = \frac{(\gamma_2 + \Gamma_2 - 1)(\Gamma_2 - 1)}{(5\Gamma_2 - 4\gamma_2 + 1)(\gamma_2 - \Gamma_2)} $$

$$ y_2 = \frac{\gamma_2 \sqrt{3}\sqrt{\delta(\gamma_2, \Gamma_2)}}{2(\gamma_2 + \Gamma_2 - 1)(\Gamma_2 - 2\gamma_2 + 2)(\Gamma_2 - \gamma_2)} $$

$$ z_2 = \frac{(2\Gamma_2 - \gamma_2)(\Gamma_2 - 1)}{2(\gamma_2 + \Gamma_2 - 1)(\gamma_2 - \Gamma_2)} $$

**Theorem 3.14.** For values of the parameters $(\gamma_2, \Gamma_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);

* $\Omega$ if $(\gamma_2 - \gamma_2)(2\Gamma_2 - \gamma_2 - 2) \neq 0$;
<figure><img src="image_12.png" /><figcaption>FIGURE 10. The curves involved in the decomposition of the parameter space (with the same conventions as in Fig. 8).</figcaption></figure>

• $\Omega_a^+$ and $\Omega_a^-$ if $\eta(\gamma_2, \Gamma_2) \neq 0$ and which are real if $\delta(\gamma_2, \Gamma_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 $(\gamma_2, \Gamma_2)$ outside of a finite union of curves, the singularities of $\{D=0\}$ inside the Bloch ball in the case of water are:

• $\emptyset$;

• $\Omega$ if and only if $(\gamma_2, \Gamma_2) \in \mathcal{G}_1$;

• $\Omega_a^+$ and $\Omega_a^-$ are real and inside the Bloch ball if and only if $(\gamma_2, \Gamma_2) \in \mathcal{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 $\mathcal{G}_1$ and $\mathcal{G}_2$ are disjoint.

*Proof.* The point $\Omega$ 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

$$ \begin{cases} \Gamma_2 - 1 > 0 \\ \Gamma_2 < \gamma_2 \\ 3\Gamma_2 - 2\gamma_2 - 1 = f_2 < 0 \end{cases} \quad \text{or} \quad \begin{cases} \Gamma_2 - 1 < 0 \\ \Gamma_2 > \gamma_2 \\ f_2 > 0 \end{cases} $$
<figure><img src="image_13.png" /><figcaption>FIGURE 11. Magnification of Fig. 10 near (1, 1).</figcaption></figure>

Note that in each case, the first and third condition imply the second one.

In the first case, since $\Gamma_2 > \gamma_2$, the inequalities $0 > z_2 > -2$ are equivalent to

$$ \begin{cases} 2\Gamma_2 - \gamma_2 - 2 < 0 \\ f_4 < 0 \end{cases} $$

Both inequalities are satisfied by assumption, recovering the defining inequalities of $G_1^-$. In the second case ($\Gamma_2 > \gamma_2$), the inequalities $0 > z_2 > -2$ are equivalent to

$$ \begin{cases} 2\Gamma_2 - \gamma_2 - 2 > 0 \\ f_4 > 0 \end{cases} $$

The first inequality is satisfied whenever the second is, so have recovered 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 $\Pi$ are:

* if $(2\Gamma_2 - \gamma_2)(\Gamma_2 - \gamma_2)(2\Gamma_2 - \gamma_2 - 2) \neq 0$: the center O of B, and the point $\Omega$;
* if $2\Gamma_2 = \gamma_2: O, \Omega$ and $N$;

* if $\Gamma_2 \neq 1$, and $(\Gamma_2 - \gamma_2)(2\Gamma_2 - \gamma_2 - 2) = 0: O$.

* if $\Gamma_2 = 1$ and $\gamma_2 \neq 1: O$ and $\Omega = N$.

*Proof.* We split the proof in seven cases, finding the solutions such that $z_2 = -1$, $2\Gamma_2 = \gamma_2$, $\Gamma_2 = 1$, $z_2 = 0$, $2\Gamma_2 = \gamma_2$, $\Gamma_2 = \gamma_2$, $2\Gamma_2 - \gamma_2 - 2 = 0$, and the remaining general case.

(1) Case $z_2 = -1$. We compute a basis of $S \cup \{y_1, y_2, z_2 + 1\}$ with the lexicographical order $y_1 \gg y_2 \gg z_1 \gg z_2 \gg \gamma_2 \gg \Gamma_2$. The polynomial $\gamma_2^2(z_1 + 1)$ is present, and since $\gamma_2 \neq 0$, we get $z_2 = -1$ and derive directly that the unique solution is $O$.

(2) Case $2\Gamma_2 = \gamma_2$. We compute a basis of $S \cup \{y_1, y_2, 2\Gamma_2 - \gamma_2\}$ with the same order $y_1 \gg y_2 \gg z_1 \gg z_2 \gg \gamma_2 \gg \Gamma_2$, and the result follows immediately.

(3) Case $\Gamma_2 = 1$. We found in the basis of $S \cup \{y_1, y_2, \Gamma_2 - 1\}$ with order $y_1 \gg y_2 \gg z_1 \gg z_2 \gg \gamma_2 \gg \Gamma_2$, the polynomial $z_2\gamma_2^2(\gamma_2 - 1)^2(z_2 + 1)$. The case $z_2 + 1$ is already treated, while the case $\gamma_2 = 0$ is excluded. The case $\gamma_2 = 1$ is also excluded, since that would mean $(\gamma_2, \Gamma_2) = (1, 1)$. In the case $z_2 = 0$, the polynomial $2\gamma_2^3 z_1(\gamma_2 - 2)$ appears in the basis of $S \cup \{y_1, y_2, \Gamma_2 - 1, z_2\}$. Since the case $\gamma_2 = 2 = 2\Gamma_2$ has already been studied, this implies $z_1 = 0$ hence the north pole $N$ is the unique solution. Note that for $\Gamma_2 = 1, N$ and $\Omega$ coincide.

(4) Case $z_2 = 0$. We compute a basis with the lexicographic $y_1 \gg y_2 \gg z_1 \gg z_2 \gg \gamma_2 \gg \Gamma_2$. The first polynomial is $(\Gamma_2 - 1)^2(2\Gamma_2 - \gamma_2)$. Since $2\Gamma_2 - \gamma_2 = 0$ and $\Gamma_2 = 1$ are already known, there is no extra solution.

(5) Case $\Gamma_2 = \gamma_2$. The basis contains the polynomial $\Gamma_2(\Gamma_2 - 1)^2(z_2 + 1)$, which leads to $z_2 = -1$ (already discussed).

(6) Case $2\Gamma_2 - \gamma_2 - 2 = 0$. The basis contains the polynomial $(\Gamma_2 - 1)^2(z_2 + 1)$, as in the previous case.

(7) Let us finally search the solutions such that $z_2(z_2 + 1)(\gamma_2 - \Gamma_2)(2\Gamma_2 - \gamma_2 - 2) \neq 0$. We consider the system

$$S \cup \{y_1, y_2, z_2 \xi_1 - 1, (z_2 + 1) \xi_2 - 1, (\Gamma_2 - \gamma_2) \lambda_1 - 1, (2\Gamma_2 - \gamma_2 - 2) \lambda_2 - 1\},$$

and compute a basis with order $y_1 \gg y_2 \gg z_1 \gg z_2 \gg \xi_1 \gg \xi_2 \gg \gamma_2 \gg \Gamma_2 \gg \lambda_1 \gg \lambda_2$. The last four polynomials of the basis are:

$$\Gamma_2 \lambda_1 + 2 \lambda_2 \lambda_1 + \lambda_1 - 2 \lambda_2 + z_2, -\Gamma_2 \lambda_1 + \lambda_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 $\Omega$ is a solution satisfying these conditions, hence it is the unique solution. $\square$

Solutions outside the symmetry plane $\{y_1 = y_2 = 0\}$.

**Proposition 3.18.** For values of the parameters $(\gamma_2, \Gamma_2)$ outside of a finite union of curves, the solutions of $S$ outside $\Pi$ are given by means of a two parameters family of $\Pi$-symmetric points $(\Omega_a^+(\gamma_2, \Gamma_2), \Omega_a^-(\gamma_2, \Gamma_2))$, which can be either real or complex depending on the sign of $\delta(\gamma_2, \Gamma_2)$.

*Proof of Proposition 3.18.*

*Proof.* The case $y_2 = 0, y_1 \neq 0$ is treated in Lemma 3.19. We assume now that $y_2 \neq 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 $\gamma_2 = h \Gamma_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 $\tilde{S}$ this new set of polynomials in the variables $(h, \Gamma_2, z_1, z_2, t, Y_2)$.
We first compute a Gröbner basis, considering the polynomials in $Q(\gamma_2, \Gamma_2)[y_1, z_1, y_2, z_2]$ with a lexicographical order, and find directly the rational expressions

$$ 
\begin{align*}
z_1 &= \frac{(1 - (1+h)\Gamma_2)(\Gamma_2 - 1)}{(4\Gamma_2h - 5\Gamma_2 - 1)\Gamma_2(h-1)} \\
z_2 &= \frac{1}{2} \frac{(\Gamma_2 - 1)(h-2)}{(1 - (1+h)\Gamma_2)(h-1)} \\
Y_2 &= \frac{3}{4} \frac{h^2\Gamma_2(\Gamma_2 - 1)(h-2)(2\Gamma_2h - \Gamma_2 - 1)}{(2\Gamma_2h - \Gamma_2 - 2)(\Gamma_2h + \Gamma_2 - 1)^2(h-1)^2} \\
t &= 2 \frac{(\Gamma_2h + \Gamma_2 - 1)(2\Gamma_2h - \Gamma_2 - 2)}{(4\Gamma_2h - 5\Gamma_2 - 1)h\Gamma_2}
\end{align*}
\tag{SoluGen} $$

which correspond to the points $\Omega_a^+$, $\Omega_a^-$. These points are symmetric with respect to $\Pi$ and real if and only if

$$ (\Gamma_2 - 1)(2\Gamma_2 - h\Gamma_2)(\Gamma_2 - 2h\Gamma_2 + 1)(\Gamma_2 - 2h\Gamma_2 + 2) \ge 0 $$

Recall that $h\Gamma_2 = \gamma_2$, so we have recovered the definition of $\delta$ (30). We define

$$ \bar{\eta}(\Gamma_2, h) = (4\Gamma_2h - 5\Gamma_2 - 1)(h-1)(2\Gamma_2h - \Gamma_2 - 2)(1 - (1+h)\Gamma_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 $\{(\gamma_2, \Gamma_2) \text{ s.t. } \bar{\eta}(\Gamma_2, h) = 0\}$. $\square$

**Lemma 3.19.** If $\gamma_2 \neq 2\Gamma_2$ and $\gamma_2 \neq 0$, the only solutions of the system $S \cup \{y_2\}$ are the solutions on $\Pi$ 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 \cup \{y_2\}$ with respect to the order $z_1 \gg z_2 \gg y_1 \gg y_2 \gg \gamma_2 \gg \Gamma_2$ is $\gamma_2^2 y_1 (2\Gamma_2 - \gamma_2)$. $\square$

### 3.4.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 $(\gamma_1, \Gamma_1, \gamma_2, \Gamma_2)$, which we may reduce to 3, normalizing by $\gamma_1 = 1$. The constraints are $2\Gamma_1 \ge 1$, $2\Gamma_2 \ge \gamma_2 > 0$, and we exclude the case of a single spin $(\Gamma_2, \gamma_2) = (\Gamma_1, 1)$.

In this case, real roots classification algorithms [10] allow to compute a generic classification of the parameter space according to the number of singularities.

**Theorem 3.20** ([10]). Splitting the subset of $\mathbb{R}^3$ defined by $2\Gamma_2 \ge \gamma_2 > 0$ and $2\Gamma_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.

**Theorem 3.21.** For values of the parameters $(\gamma_1, \Gamma_1, \gamma_2, \Gamma_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><img src="image_14.png" /><figcaption>FIGURE 12. All polynomials involved in the classification in the *in vivo* case</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIGURE 13. Classification regarding the second parameter in the *in vivo* case</figcaption></figure>

As an example, we show slices of this classification for $\Gamma_1 = 531/26 \approx 20.4$, which corresponds to the *in vivo* experimental setting. All polynomials involved are represented on Figure 12. In this case, there can be 1, 2, 3 or 4 singularities, and the corresponding sample points are represented respectively with gray dots, blue diamonds, green circles and orange boxes.

In the experimental setting, the parameters for the second matter were $\gamma_2 \approx 0.89$ and $\Gamma_2 \approx 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 $(\gamma_2, \Gamma_2) = (\gamma_1, \Gamma_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.
REFERENCES

[1] P. Aubry, D. Lazard, and M. Moreno Maza. On the theories of triangular sets. *Journal of Symbolic Computation*, 28(1-2):105-124, July 1999.

[2] S. Basu, R. Pollack, and M.-F. Roy. Algorithms in real algebraic geometry, volume 10 of algorithms and computation in mathematics, 2006.

[3] B. Bonnard. Feedback equivalence for nonlinear systems and the time optimal control problem. *SIAM J. Control Optim.*, 29(6):1300–1321, 1991.

[4] B. Bonnard and M. Chyba. *Singular trajectories and their role in control theory*, volume 40 of *Mathématiques & Applications (Berlin) [Mathematics & Applications]*. Springer-Verlag, Berlin, 2003.

[5] B. Bonnard, M. Chyba, A. Jacquemard, and J. Marriott. Algebraic geometric classification of the singular flow in the contrast imaging problem in nuclear magnetic resonance. *Math. Control Relat. Fields*, 3(4):397–432, 2013.

[6] B. Bonnard, M. Chyba, and J. Marriott. Singular trajectories and the contrast imaging problem in nuclear magnetic resonance. *SIAM J. Control Optim.*, 51(2):1325–1349, 2013.

[7] B. Bonnard, M. Claeys, O. Cots, and P. Martinon. Geometric and numerical methods in the contrast imaging problem in nuclear magnetic resonance. *Acta Appl. Math.*, 135:5–45, 2015.

[8] B. Bonnard and O. Cots. Geometric numerical methods and results in the contrast imaging problem in nuclear magnetic resonance. *Math. Models Methods Appl. Sci.*, 24(1):187–212, 2014.

[9] B. Bonnard, O. Cots, S. J. Glaser, M. Lapert, D. Sugny, and Y. Zhang. Geometric optimal control of the contrast imaging problem in nuclear magnetic resonance. *IEEE Trans. Automat. Control*, 57(8):1957–1969, 2012.

[10] B. Bonnard, J.-C. Faugère, A. Jacquemard, M. Safey El Din, and T. Verron. Determinantal sets, singularities and application to optimal control in medical imagery. In *Proceedings of the 2016 International Symposium on Symbolic and Algebraic Computation*, ISSAC '16, 2016.

[11] C. W Brown. Qepcad b: a program for computing with semi-algebraic sets using cads. *ACM SIGSAM Bulletin*, 37(4):97–108, 2003.

[12] W. Bruns and U. Vetter. *Determinantal rings*, volume 1327 of *Lecture Notes in Mathematics*. Springer-Verlag, Berlin, 1988.

[13] B. Buchberger. A theoretical basis for the reduction of polynomials to canonical forms. *ACM SIGSAM Bulletin*, 10(3):19–29, 1976.

[14] J. Canny. *The complexity of robot motion planning*, volume 1987 of *ACM Doctoral Dissertation Awards*. MIT Press, Cambridge, MA, 1988.

[15] J. Canny and I. Emiris. An efficient algorithm for the sparse mixed resultant. In *Applied algebra, algebraic algorithms and error-correcting codes* (*San Juan, PR, 1993*), volume 673 of *Lecture Notes in Comput. Sci.*, pages 89–104. Springer, Berlin, 1993.

[16] E. Cattani and A. Dickenstein. Introduction to residues and resultants. In *Solving polynomial equations*, pages 1–61. Springer Berlin Heidelberg, 2005.

[17] C. Chen and M. Moreno Maza. Algorithms for computing triangular decomposition of polynomial systems. *J. Symbolic Comput.*, 47(6):610–642, 2012.

[18] C. Chen and M. Moreno Maza. Quantifier Elimination by Cylindrical Algebraic Decomposition Based on Regular Chains. *J. Symbolic Comput.*, 75:74–93, 2016.

[19] G. E. Collins. Quantifier elimination for real closed fields by cylindrical algebraic decomposition. In *Automata theory and formal languages (Second GI Conf., Kaiserslautern, 1975)*, volume 33 of *Lecture Notes in Computer Science*, pages 134–183. Springer, Berlin, 1975.

[20] S. Conolly, D. Nishimura, and A. Macovski. Optimal control solutions to the magnetic resonance selective excitation problem. *IEEE transactions on medical imaging*, 5(2):106–115, 1986.

[21] O. Cots. *Contrôle optimal géométrique : méthodes homotopiques et applications*. PhD thesis, Université de Bourgogne, Dijon, France, 2012.

[22] D. A. Cox, J. Little, and D. O’Shea. *Using algebraic geometry*, volume 185 of *Graduate Texts in Mathematics*. Springer, New York, second edition, 2005.

[23] A. Dickenstein and I. Z. Emiris, editors. *Solving polynomial equations*, volume 14 of *Algorithms and Computation in Mathematics*. Springer-Verlag, Berlin, 2005. Foundations, algorithms, and applications.

[24] M. Safey El Din and É. Schost. Polar varieties and computation of one point in each connected component of a smooth real algebraic set. In J. Rafael Sendra, editor, *Symbolic and Algebraic Computation, International Symposium ISSAC 2003, Drexel University, Philadelphia, Pennsylvania, USA, August 3-6, 2003, Proceedings*, pages 224–231. ACM, 2003.

[25] J.-C. Faugère. A new efficient algorithm for computing Gröbner bases ($F_4$). *J. Pure Appl. Algebra*, 139(1-3):61–88, 1999. Effective methods in algebraic geometry (Saint-Malo, 1998).

[26] J.-C. Faugère. A new efficient algorithm for computing Gröbner bases without reduction to zero ($F_5$). In *Proceedings of the 2002 International Symposium on Symbolic and Algebraic Computation*, pages 75–83 (electronic), New York, 2002. ACM.

[27] J.-C. Faugère, M. Safey El Din, and P.-J. Spaenlehauer. Critical points and gröbner bases: the unmixed case. In J. van der Hoeven and M. van Hoeij, editors, *International Symposium on Symbolic and Algebraic Computation, ISSAC'12, Grenoble, France - July 22 - 25, 2012*, pages 162–169. ACM, 2012.

[28] J.-C. Faugère. FGb: A Library for Computing Gröbner Bases. In K. Fukuda, J. Hoeven, M. Joswig, and N. Takayama, editors, *Mathematical Software - ICMS 2010, volume 6327 of Lecture Notes in Computer Science*, pages 84–87, Berlin, Heidelberg, September 2010. Springer Berlin / Heidelberg.

[29] G.-M. Greuel and G. Pfister. *A Singular introduction to commutative algebra*. Springer, Berlin, extended edition, 2008. With contributions by Olaf Bachmann, Christoph Lossen and Hans Schönemann, With 1 CD-ROM (Windows, Macintosh and UNIX).

[30] A. Greuet and M. Safey El Din. Probabilistic algorithm for polynomial optimization over a real algebraic set. *SIAM Journal on Optimization*, 24(3):1313–1343, 2014.

[31] D. Yu. Grigor’ev and N.N. Vorobjov. Solving systems of polynomial inequalities in subexponential time. *Journal of Symbolic Computation*, 5(1):37–64, 1988.
[32] S. Helgason. *Differential geometry, Lie groups, and symmetric spaces*, volume 80 of *Pure and Applied Mathematics*. Academic Press, Inc. [Harcourt Brace Jovanovich, Publishers], New York-London, 1978.

[33] N. Khaneja, S. J. Glaser, and R. Brockett. Sub-Riemannian geometry and time optimal control of three spin systems: quantum gates and coherence transfer. *Phys. Rev. A* (3), 65(3, part A):032301, 11, 2002.

[34] N. Khaneja, T. Reiss, C. Kehlet, T. Schulte-Herbrüggen, and S. J. Glaser. Optimal control of coupled spin dynamics: design of NMR pulse sequences by gradient ascent algorithms. *J. of Magnetic Resonance*, 172(2):296–305, 2005.

[35] T. Krick and A. Logar. An algorithm for the computation of the radical of an ideal in the ring of polynomials. In *Applied algebra, algebraic algorithms and error-correcting codes (New Orleans, LA, 1991)*, volume 539 of *Lecture Notes in Comput. Sci.*, pages 195–205. Springer, Berlin, 1991.

[36] M. 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*. PhD thesis, Université de Bourgogne, Dijon, France, 2011.

[37] M. Lapert, Y. Zhang, M. A. Janich, S. J. Glaser, and D. Sugny. Exploring the Physical Limits of Saturation Contrast in Magnetic Resonance Imaging. *Scientific Reports*, 2:589, 2012.

[38] M. Lapert, Y. Zhang, M. Braun, S. J. Glaser, and D. Sugny. Singular extremals for the time-optimal control of dissipative spin $\frac{1}{2}$ particles. *Phys. Rev. Lett.*, 104(8, 083001), 2010.

[39] D. Lazard and F. Rouillier. Solving parametric polynomial systems. *Journal of Symbolic Computation*, 42(6):636–667, 2007.

[40] F. S. MacAulay. Some Formulae in Elimination. *Proc. London Math. Soc.*, S1-35(1):3.

[41] L. Markus. Quadratic differential equations and non-associative algebras. In *Contributions to the theory of nonlinear oscillations, Vol. V*, pages 185–213. Princeton Univ. Press, Princeton, N.J., 1960.

[42] S. McCallum and H. Hong. On using lizard's projection in CAD construction. *J. Symb. Comput.*, 72:65–81, 2016.

[43] L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze, and E. F. Mishchenko. *The mathematical theory of optimal processes*. Translated from the Russian by K. N. Trirogoff; edited by L. W. Neustadt. Interscience Publishers John Wiley & Sons, Inc. New York-London, 1962.

[44] K. Ranestad. Algebraic Degree in Semidefinite and Polynomial Optimization. In Miguel F. Anjos and Jean-Bernard Lasserre, editors, *Handbook on Semidefinite, Conic and Polynomial Optimization*, chapter 3, pages 61–75. Springer, 2012.

[45] E. Van Reeth, H. Ratiney, M. Tesch, D. Grenier, O. Beuf, S. J. Glaser, and D. Sugny. Optimal control design of preparation pulses for contrast optimization in {MRI}. *Journal of Magnetic Resonance*, 279:39–50, 2017.

[46] M. Safey El Din. Raglib (real algebraic geometry library). 2007. Maple package. http://www-polsys.lip6.fr/safey/RAGLib/.

[47] M. Safey El Din. Testing sign conditions on a multivariate polynomial and applications. *Mathematics in Computer Science*, 1(1):177–207, 2007.

[48] T. E. Skinner, T. O. Reiss, B. Luy, N. Khaneja, and S. J. Glaser. Application of optimal control theory to the design of broadband excitation pulses for high-resolution nmr. *Journal of magnetic resonance (San Diego, Calif. : 1997)*, 163(1):8–15, July 2003.

[49] D. Sugny, E. Van Reeth, H. Ratiney, S.J. Glaser, and M. Lapert. Optimal control theory for applications in magnetic resonance imaging. *Pacific Journal of Mathematics for Industry*, Submitted 2017.

[50] D. Wang. *Elimination methods*. Texts and Monographs in Symbolic Computation. Springer-Verlag, Vienna, 2001.

[51] L. Yang, X. Hou, and B. Xia. A complete algorithm for automated discovering of a class of inequality-type theorems. *Sci. China Ser. F*, 44(1):33–49, 2001.