Bloch equation

The mathematical model which 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 whose evolution satisfies the so-called Bloch equation

\[ \dot{M}(\tau) = \bar{\gamma}\, M(\tau) \wedge B(\tau) + R(M(\tau)),\]

where $\tau$ is the time, $\bar{\gamma}$ is the gyromagnetic ratio of the considered nucleus, $B(\tau) = (B_X(\tau), B_Y(\tau), B_Z(\tau))$ is the total magnetic field applied to the system which decomposes into

\[ B(\tau) = B_0 + B_1(\tau),\]

where $B_0 = B_Z\, e_Z$ is a constant and strong polarizing field in the $Z$-direction, while the control RF-field $B_1(\tau) = B_X(\tau)\, e_X + B_Y(\tau)\, e_Y$ is in the transverse plane $(X,Y)$. The vectors $e_X$, $e_Y$ and $e_Z$ form the standard basis of $\R^3$. The $R(M)$ term represents the dissipation and is 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, and $T_1$, $T_2$ are the relaxation times which are the chemical signatures of the observed species and satisfy the physical constraints $0 < T_2 \le 2\, T_1$. The control is denoted $\omega(\tau) = (\omega_X(\tau),\omega_Y(\tau)) \coloneqq (-\bar{\gamma} B_X(\tau), -\bar{\gamma} B_Y(\tau))$ and is bounded by $\omega_\mathrm{max}$, i.e. $\Vert{\omega(\tau)}\Vert\le \omega_\mathrm{max}$, where $\omega_\mathrm{max}$ is the maximal experimental intensity of the experiments and $\Vert{\cdot}\Vert$ is the Euclidean norm.

The Bloch equation can be written in a rotating frame where the equilibrium is normalized introducing the state $q=(x,y,z)$, the matrices

\[ \Omega_z = \begin{bmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}\]

and $S(\tau) = \exp(\tau \, \bar{\omega}\, \Omega_z)$, and writting $M(\tau) = M_0 S(\tau) q(\tau)$. In the rotating frame the Bloch equation becomes:

\[ \dot{q}(\tau) = \begin{bmatrix} -{1}/{T_2} & -\Delta\omega & \omega_y \\ \Delta\omega & -{1}/{T_2} & -\omega_x \\ -\omega_y & \omega_x & -{1}/{T_1} \\ \end{bmatrix}\, q(\tau) + \begin{bmatrix} 0 \\ 0 \\ {1}/{T_1} \end{bmatrix},\]

where $\Delta \omega = \omega_0 - \bar{\omega}$ is the resonance offset, $\omega_0 = -\bar{\gamma} B_Z$ is the resonance frequency and the new control is

\[ \begin{aligned} \omega_x(\tau) \coloneqq \omega_X(\tau) \cos(\bar{\omega} \tau) + \omega_Y(\tau) \sin(\bar{\omega} \tau), \\ \omega_y(\tau) \coloneqq \omega_Y(\tau) \cos(\bar{\omega} \tau) - \omega_X(\tau) \sin(\bar{\omega} \tau), \end{aligned}\]

which preserves the control bound $\omega_x^2(\tau) + \omega_y^2(\tau) \le \omega_\mathrm{max}$. Finally, we introduce the normalized control and time:

\[ u \coloneqq \frac{u_\mathrm{max}}{\omega_\mathrm{max}}\, \omega, \quad t \coloneqq \frac{\omega_\mathrm{max}}{u_\mathrm{max}}\, \tau,\]

such that the normalized control satisfies $\Vert{u(t)}\Vert\le u_\mathrm{max}$. In the sequel, we fix $\bar{\omega} = \omega_0$ which gives $\Delta \omega = 0$ (it is called the resonant case) and which leads to the final normalized Bloch equation:

\[ \left\{ \begin{aligned} \dot{x}(t) &= \displaystyle -\Gamma\, x(t) + u_2(t)\, z(t), \\[0.2em] \dot{y}(t) &= \displaystyle -\Gamma\, y(t) - u_1(t)\, z(t), \\[0.2em] \dot{z}(t) &= \displaystyle \gamma\,(1-z(t)) + u_1(t)\, y(t) - u_2(t)\, x(t), \end{aligned} \right.\]

where

\[ \gamma \coloneqq \frac{u_\mathrm{max}}{w_\mathrm{max}\, T_1}, \quad \Gamma \coloneqq \frac{u_\mathrm{max}}{w_\mathrm{max}\, T_2}, \quad 0 < \gamma \le 2\Gamma.\]

Time parameterization

In this setting, the normalized Bloch equation has 3 positive parameters: $\gamma$, $\Gamma$ and $u_\mathrm{max}$. However, one may choose one parameter to fix since for any $\lambda > 0$ and any triplet $(\gamma, \Gamma, u_\mathrm{max})$, both systems, in coordinates $(x, y, z)$, with parameters $(\gamma, \Gamma, u_\mathrm{max})$ or $(\lambda\gamma, \lambda\Gamma, \lambda u_\mathrm{max})$ are equivalent, up to a time reparameterization.

Note that for the constant control $u(\cdot) =(u_\mathrm{max},0)$, starting from $x=0$ and denoting $q = (y,z)$, the Bloch equation becomes

\[ \dot{x}(t) = -\Gamma\, x(t), \quad \dot{q}(t) = u_\mathrm{max} \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}\, q(t),\]

and the solution is given by $x(t) = 0$, $q(t) = R(u_\mathrm{max}\, t) \, q(0)$ where $R(\theta)$ is the rotation matrix of angle $\theta$. Thus, this trajectory is periodic of period $T$ given by:

\[ T = \frac{2 \pi}{u_\mathrm{max}}.\]

Classically, one may fix the period to $T=1$ and so $u_\mathrm{max} = 2 \pi$, or on the contrary one may fix $u_\mathrm{max} = 1$ and choose $T = 2 \pi$. From now, we choose to fix $u_\mathrm{max} = 1$.

Spherical coordinates

The Bloch equation in the spherical coordinates

\[ \begin{aligned} x = \rho\, \sin\phi \cos\theta, \quad %\\[0.1em] y = \rho\, \sin\phi \sin\theta, \quad %\\[0.1em] z = \rho\, \cos\phi, \end{aligned}\]

with the feedback $u = R(\theta)^{-1} v$ becomes:

\[ \left\{ \begin{aligned} \dot{\rho}(t) &= \gamma \cos \phi(t)\, (1 - \rho(t) \cos \phi(t)) - \Gamma \rho(t) \sin^2 \phi(t), \\[1.0em] \dot{\phi}(t) &= \delta \sin \phi(t) \cos \phi(t) - \frac{\gamma}{\rho(t)} \sin \phi(t) + v_2(t), \\[0.0em] \dot{\theta}(t) &= - \cot \phi(t)\, v_1(t), \end{aligned} \right.\]

where $\delta \coloneqq \gamma - \Gamma$ and with the control constraint $v_1^2 + v_2^2 = u_1^2 + u_2^2 \le 1$.

Details.

Inverting the spherical change of coordinates, we get

$\rho^2 = x^2 + y^2 + z^2, \quad %\\[0.1em] \theta = \arctan(y/x), \quad %\\[0.1em] \phi = \arccos(z/\rho).$

To get the dynamics with the new coordinates, we use the following procedure: let us assume that we have a multi-inputs affine control system of the form $\dot{q} = X(q) + \sum_{i=1}^{m} u_i Y_i(q)$ , where $X$ , $Y_1$ , ..., $Y_m$ are vector fields. Considering that we are working in local coordinates, we can define the matrix $Y(q) = [Y_1(q) \cdots Y_m(q)] \in \R^{n\times m}$ , where $n$ is the dimension of the state space and $m$ the dimension of the control space. Let $\varphi\colon \R^n \to \R^n$ be a diffeomorphism and $q'$ the new coordinates given by $q' = \varphi(q)$ . The dynamics in the new coordinates with the affine feedback $u = A(q)\, v + b(q)$ , $A(q) \in \R^{m \times m}$ invertible, $b(q) \in \R^m$ , is given by the new vector fields

$\left\{ \begin{aligned} X' &= \varphi_* X + (\varphi_* Y) \cdot (b \circ \varphi^{-1}), \\[0.2em] Y' &= (\varphi_* Y)\cdot (A\circ \varphi^{-1}), \end{aligned} \right.$

where $\varphi_* X$ is the image (or push-forward) of $X$ by $\varphi$ defined by $\varphi_* X (q') = \mathrm{d} \varphi (\varphi^{-1}(q')) \cdot X(\varphi^{-1}(q'))$ and we define in a similar way $\varphi_* Y (q') = [\varphi_* Y_1(q') \cdots \varphi_* Y_m(q') ] \in \R^{n\times m}$ .

The Bloch ball

Definition. The closed unit ball is called the Bloch Ball.

Proposition. The Bloch ball is invariant under the Bloch equations if and only if $0 \le \gamma \le 2 \Gamma$ and it is the smallest invariant closed ball centered at the origin if and only if $0 < \gamma \le 2\, \Gamma$.

Proof.

Consider the Bloch equations in spherical coordinates. We want to prove first that $\dot{\rho} \le 0$ for every point on the sphere, i.e. such that $\rho = 1$ . On such points we have:

$\dot{\rho} = -\delta \cos^2 \phi + \gamma \cos \phi - \Gamma = -\delta X^2 + \gamma X - \Gamma = P(X), \quad X = \cos \phi.$

This is a degree 2 polynomial vanishing at $X=1$ , and its discriminant is $(2\Gamma - \gamma)^2 \ge 0$ , so $P$ has at most another root $X_2 = {\Gamma}/{\delta}$ . If $\delta < 0$ , the polynomial function is negative between its roots, so $P$ is non-positive over $[-1, 1]$ if and only if $X_2 \le -1$ , that is: ${\Gamma}/{\delta} \le -1 \Leftrightarrow \gamma \ge 0$ . If $\delta > 0$ , we need to have $X_2 \ge 1$ , that is ${\Gamma}/{\delta} \ge 1 \Leftrightarrow \gamma \le 2\Gamma$ . If $\delta = 0$ then $P$ is non-positive over $[-1, 1]$ if and only if $\gamma = \Gamma \ge 0$ .

Hence, the Bloch ball is invariant under the Bloch equations if and only if $0 \le \gamma \le 2 \Gamma$ but it is the smallest invariant ball if and only if we also have $\gamma \ne 0$ since if $\gamma = 0$ then every ball centered at the origin is invariant under the Bloch equations and finally if $0 < \gamma \le 2 \Gamma$ , then the dynamics given by the normalized Bloch equation with a zero control (i.e. the free dynamics) is linear, its equilibrium point is the North pole $N=(0,0,1)$ and this equilibrium is globally asymptotically stable so any trajectory solution of the free dynamics would leave any ball centered at the origin strictly smaller than the unit ball. The result is proved.

Symmetry of revolution

Considering the Bloch equation in spherical coordinates, one can see that $\theta$ does not appear in the dynamics (it is a cyclic variable) so there exists a one dimensional symmetry group of translations on $\theta$. Hence, one may fix the initial value $\theta(0)$ to $\pi/2$ for instance (it implies $x(0) = 0$). Finally, imposing $\theta(0) = \pi/2$ and $u_2 = 0$, then, one can reduce (since $x(t)=0$ for any time $t$) the Bloch equation to the following two dimensional single-input control system:

\[ \left\{ \begin{aligned} \dot{y}(t) &= \displaystyle -\Gamma\, y(t) - u(t)\, z(t), \\[0.2em] \dot{z}(t) &= \displaystyle \gamma\,(1-z(t)) + u(t)\, y(t), \end{aligned} \right.\]

where $u \coloneqq u_1$ satisfies $|{u(t)}| \le 1$.

This control system is affine in the control and may be written in the form

\[ \dot{q}(t) = F_0(q(t)) + u(t)\, F_1(q(t)),\]

where $q \coloneqq (y,z)$ is the state and where the vector fields $F_0$ and $F_1$ are given by:

\[ F_0(q) = -\Gamma\, y \frac{\partial}{\partial y} + \gamma\, (1 - z) \frac{\partial}{\partial z}, \quad F_1(q) = -z \frac{\partial}{\partial y} + y \frac{\partial}{\partial z}.\]

The admissible control set $\mathcal{U}$ is then defined as

\[ \mathcal{U} \coloneqq \left\{ u(\cdot) \colon [0, \infty) \to [-1, 1] ~|~ u(\cdot) \text{ measurable} \right\}.\]

Restriction to the 2D case for the time-optimal saturation problem

Let us consider the Bloch equations in spherical coordinates and take an initial point $(\rho_0, \phi_0, \theta_0)$ in the Bloch ball. Let $\bar{v}(\cdot) = (\bar{v}_1(\cdot), \bar{v}_2(\cdot))$ denotes an optimal control which steers $(\rho_0, \phi_0, \theta_0)$ to the target $\rho = 0$ in time $\bar{t}_f$ . The control law satisfies $\bar{v}_1^2(t) + \bar{v}_2^2(t) \le 1$ for all $t \in [0, \bar{t}_f]$ . Since $\theta$ is a cyclic variable and since the component $v_1$ appears only in the $\theta$ dynamics, then, one can define the control law $v^*(\cdot) = (0, \bar{v}_2(\cdot))$ which satisfies the constraint and also steers the system from the initial point $(\rho_0, \phi_0, \theta_0)$ to the target $\rho = 0$ in the same time $\bar{t}_f$ , and thus, the control $v^*(\cdot)$ is also optimal.

On the contrary, let us consider an optimal control of the form $\bar{v}(\cdot) = (0, \bar{v}_2(\cdot))$ with $\bar{t}_f$ the optimal time. If there exists times $0 \le t_1 < t_2 \le \bar{t}_f$ such that $\bar{v}_2(t)^2 < 1$ on $[t_1, t_2]$ , then clearly, one can define a new control $v^*(\cdot) = (v_1^*(\cdot), \bar{v}_2(\cdot))$ such that $v_1^*(\cdot) \ne 0$ on $[t_1, t_2]$ , which satisfies the control constraint over $[0, \bar{t}_f]$ , and again this control $v^*(\cdot)$ steers the system to $\rho=0$ in time $\bar{t}_f$ and thus is optimal.

As a conclusion, if there exists an optimal control $\bar{v}(\cdot) = (\bar{v}_1(\cdot), \bar{v}_2(\cdot))$ , then there exists one of the form $v^*(\cdot) = (0, \bar{v}_2(\cdot))$ . Besides, if there exists an optimal control of the form $\bar{v}(\cdot) = (0, \bar{v}_2(\cdot))$ such that $\bar{v}_2(\cdot)^2 \ne 1$ , then there exist optimal controls of the form $v^*(\cdot) = (v_1^*(\cdot), \bar{v}_2(\cdot))$ with $v_1^*(\cdot) \ne 0$ . In any case, it is reasonable to consider only the optimal controls of the form $\bar{v}(\cdot) = (0, \bar{v}_2(\cdot))$ .

Besides, since $\theta$ is a cyclic variable, one can fix $\theta_0 = \pi/2$ and in the case $v_1(\cdot) = 0$ and $\theta(0) = \pi/2$ , we have $\theta(\cdot) = 0$ and $v_1(\cdot) = u_2(\cdot) = 0$ .