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

ρ2=x2+y2+z2,θ=arctan(y/x),ϕ=arccos(z/ρ). \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 q˙=X(q)+i=1muiYi(q)\dot{q} = X(q) + \sum_{i=1}^{m} u_i Y_i(q), where XX, Y1Y_1, ..., YmY_m are vector fields. Considering that we are working in local coordinates, we can define the matrix Y(q)=[Y1(q)Ym(q)]Rn×mY(q) = [Y_1(q) \cdots Y_m(q)] \in \R^{n\times m}, where nn is the dimension of the state space and mm the dimension of the control space. Let φ ⁣:RnRn\varphi\colon \R^n \to \R^n be a diffeomorphism and qq' the new coordinates given by q=φ(q)q' = \varphi(q). The dynamics in the new coordinates with the affine feedback u=A(q)v+b(q)u = A(q)\, v + b(q), A(q)Rm ×mA(q) \in \R^{m \times m} invertible, b(q) Rmb(q) \in \R^m, is given by the new vector fields

{X=φX+(φY)(bφ1),Y=(φY)(Aφ1), \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 φX\varphi_* X is the image (or push-forward) of XX by φ\varphi defined by φX(q)=dφ(φ1(q))X(φ1(q))\varphi_* X (q') = \mathrm{d} \varphi (\varphi^{-1}(q')) \cdot X(\varphi^{-1}(q')) and we define in a similar way φY(q)=[φY1(q)φYm(q)]Rn×m\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 ρ˙0\dot{\rho} \le 0 for every point on the sphere, i.e. such that ρ=1\rho = 1. On such points we have:

ρ˙=δcos2ϕ+γcosϕΓ=δX2+γXΓ=P(X),X=cosϕ. \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=1X=1, and its discriminant is (2Γγ)20(2\Gamma - \gamma)^2 \ge 0, so PP has at most another root X2=Γ/δX_2 = {\Gamma}/{\delta}. If δ<0\delta < 0, the polynomial function is negative between its roots, so PP is non-positive over [1,1][-1, 1] if and only if X21X_2 \le -1, that is: Γ/δ1γ0{\Gamma}/{\delta} \le -1 \Leftrightarrow \gamma \ge 0. If δ>0\delta > 0, we need to have X21X_2 \ge 1, that is Γ/δ1γ2Γ{\Gamma}/{\delta} \ge 1 \Leftrightarrow \gamma \le 2\Gamma. If δ=0\delta = 0 then PP is non-positive over [1,1][-1, 1] if and only if γ=Γ0\gamma = \Gamma \ge 0.

Hence, the Bloch ball is invariant under the Bloch equations if and only if 0γ2Γ0 \le \gamma \le 2 \Gamma but it is the smallest invariant ball if and only if we also have γ0\gamma \ne 0 since if γ=0\gamma = 0 then every ball centered at the origin is invariant under the Bloch equations and finally if 0<γ2Γ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)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 (ρ0,ϕ0,θ0)(\rho_0, \phi_0, \theta_0) in the Bloch ball. Let vˉ()=(vˉ1(),vˉ2())\bar{v}(\cdot) = (\bar{v}_1(\cdot), \bar{v}_2(\cdot)) denotes an optimal control which steers (ρ0,ϕ0,θ0)(\rho_0, \phi_0, \theta_0) to the target ρ=0\rho = 0 in time tˉf\bar{t}_f. The control law satisfies vˉ12(t)+vˉ22(t)1\bar{v}_1^2(t) + \bar{v}_2^2(t) \le 1 for all t[0,tˉf]t \in [0, \bar{t}_f]. Since θ\theta is a cyclic variable and since the component v1v_1 appears only in the θ\theta dynamics, then, one can define the control law v()=(0,vˉ2())v^*(\cdot) = (0, \bar{v}_2(\cdot)) which satisfies the constraint and also steers the system from the initial point (ρ0,ϕ0,θ0)(\rho_0, \phi_0, \theta_0) to the target ρ=0\rho = 0 in the same time tˉf\bar{t}_f, and thus, the control v()v^*(\cdot) is also optimal.

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

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

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