A Complete Solution of the Bloch Equation

TL;DR

This paper presents a comprehensive, explicit analytical solution to the Bloch equation, addressing previous limitations and providing clear physical insights into two-level system dynamics.

Contribution

It offers a complete, general propagator solution for the Bloch equation, explicitly handling singular parameter cases and simplifying the roots calculation.

Findings

Explicit time dependence of the system is derived.

The solution accommodates division by zero in previous models.

Models of system dynamics as relaxation and rotation are developed.

Abstract

The Bloch equation is the fundamental dynamical model applicable to arbitrary two-level systems. Analytical solutions to date are incomplete for a number of reasons that motivate further investigation. The solution obtained here for the propagator, which generates the time evolution of the system and embodies all the system dynamics, is compact and completely general. The parameter space that results in division by zero in previous treatments is explicitly defined and accommodated in the solution. Polynomial roots required for the solution are expressed in terms of a single real root obtained using simple functional forms. A simple graphical rendition of this root is developed that clarifies and characterizes its dependence on the physical parameters of the problem. As a result, the explicit time dependence of the system as a function of its physical parameters is immediately evident. Several intuitive models of system dynamics are also developed. In particular, the Bloch equation is separable in the proper coordinate system, written as the sum of a relaxation operator and either a null operator or a commuting rotation. The propagator thus drives either pure relaxation or relaxation followed by a rotation. The paper provides a basis for increased physical insight into the Bloch equation and its widespread applications.

Equations

$$\Gamma=\left(\begin{array}[]{ccc}R_{1}&\mbox{$\omega_{3}$}&-\mbox{$\omega_{2}$}\\ -\mbox{$\omega_{3}$}&R_{2}&\mbox{$\omega_{1}$}\\ \mbox{$\omega_{2}$}&-\mbox{$\omega_{1}$}&R_{3}\end{array}\right)$$

$$\dot{\mathbf{M}}(t)+\Gamma \mathbf{M}(t)={\mathbf{M}}_0{R}_3.$$

$${\omega }_{12}^2$$

$${\omega }_e^2$$

$$\mathbf{M}(t)$$

$$=$$

$${\mathcal{L}}^{-1}[f(s)]$$

$$A(s)=s\mbox{$\openone$}+\Gamma,$$

$$\mathrm{adj}\,A(s)=A_{0}+A_{1}\,s+\mbox{$\openone$}\,s^{2},$$

$$(A_0{)}_{ij}$$

$$(A_1{)}_{ij}$$

$$p(s)=c_0+c_{1}s+c_2{s}^2+s^3$$

$$c_0$$

$$c_1$$

$$c_2$$

$$A^{-1}(s)=(s\mbox{$\openone$}+\Gamma)^{-1}=\dfrac{\mathrm{adj}\,A(s)}{p\,(s)}$$

$$e^{-\Gamma t}=\sum _{\mathrm{res}}\frac{\mathrm{adj}A(s)}{p(s)}{e}^{st}.$$

$$\mathrm{res}(s_0)=\frac{1}{(k-1)!}\underset{s\to s_0}{\lim }\frac{d^{k-1}}{d{s}^{k-1}}[(s-s_0{)}^{k}g(s)]$$

$${\mathbf{M}}_{\infty }$$

$$p(z-c_2/3)$$

$$a$$

$$b$$

$$z_{2,3}$$

$$\varpi^{2}=3\big{[}(z_{1}/2)^{2}+a/3\big{]},$$

$$\gamma =\frac{|b/2|}{|a/3{|}^{3/2}},$$

$$s_i=z_i-c_2/3,$$

$$\frac{c_2}{3}=\frac{1}{3}\sum _i{R}_i\equiv \bar{R}$$

$$\Gamma$$

$$R_{ip}$$

$$a$$

$$b$$

$$e^{-\Gamma t}=e^{-\bar{R}t}{e}^{-{\Gamma }_{\mathrm{p}}t}.$$

$$\mathrm{adj}{A}_{\mathrm{p}}(z)$$

Full text

broadband decoupling; heteronuclear decoupling; IS spin system; J coupling; optimal control theory

I Introduction

The Bloch equation needs little formal introduction. It was proposed originally as a classical, phenomenological model for the dissipative dynamics observed in magnetic resonance Bloch (1946). However, its impact has been more widespread as a result of its applicability to quantum two-level systems Feynman et al. (1957). The field of quantum optics is a particularly noteworthy example of its significance beyond the realm of magnetic resonance Fano (1983). More recently, the burgeoning field of quantum computing provides additional incentive for understanding Bloch equation dynamics due to the necessity of preserving coherence in the presence of relaxation Nielsen and Chuang (2000).

A thorough treatment of this fundamental dynamical model for arbitrary two-level quantum systems, including relaxation, is therefore relevant to a host of important physical systems. Yet, there is to date no complete and general solution for the Bloch equation. Previous solutions Torrey (1949); Madhu and Kumar (1995, 1997); Bain (2010) share some or all of the following limitations. They (i) divide by zero at particular, but unspecified, values of the field and relaxation rates; (ii) are not valid for arbitrary field orientation, which can be important in many applications; (iii) are cumbersome, linked to tables of multiply nested variables with obscure connection to the physical parameters of the problem; (iv) depend unavoidably on the roots of a cubic polynomial, roots that are only qualitatively specified or are expressed as complicated functions of the physical parameters; (v) provide only a small measure of the physical insight that might be expected from an analytical solution.

The present work derives a compact, complete solution to the Bloch equation. The solutions are valid for arbitrary (constant) fields and relaxation rates. The paper begins with a theoretical overview. The intent is to provide a fairly complete general understanding of the problem and the formal simplicity of the solution. The next section is devoted to the explicit form of the solutions. Their dependence on the physical parameters is fully characterized in a following section. The spin-spin (transverse) relaxation rate minus the spin-lattice (longitudinal) rate provides a convenient and particularly useful frequency scale for representing system parameters in the analysis of the Bloch equation. Conditions that result in division by zero in previous solutions are identified, yielding quantitative bounds for oscillatory (underdamped) and non-oscillatory (critically damped and underdamped) dynamics. The polynomial roots required in the solution are reduced to a single root with a fully characterized dependence on the physical parameters that admits a simple graphical representation. Next, a simple vector model is developed that reveals the underlying simplicity of the dynamics, with a modified system of relaxation rates playing a role analogous to longitudinal and transverse relaxation. The modified rates result from the interaction/coupling between the fields and relaxation processes. Additionally, and incidentally, a method for finding eigenvectors emerges that does not appear to have been considered previously. The Bloch equation is also shown to represent a system of three mutually coupled harmonic oscillators, providing additional perspective. Details of the calculations are provided in appendices. The concluding appendix applies the solutions to a representative set of cases yielding solutions that can be straightforwardly verified by other methods.

II Theoretical Overview

We first summarize the basic framework of the Bloch equation to recollect and also define the fundamental parameters of the problem. The equation describes the dynamics of a magnetization $\bm{M}$ subjected to a static (polarizing) magnetic field $\bm{H}_{0}=H_{0}\,\bm{\hat{z}}$ and a sinusoidally alternating field $2H_{a}\cos\omega_{a}t$ applied orthogonal to $\bm{H}_{0}$. For $H_{a}\ll H_{0}$, the equilibrium magnetization is not appreciably affected by the applied field and is therefore, to a good approximation, the time-independent value $\bm{M}_{0}=\chi H_{0}\,\bm{\hat{z}}$ produced by the polarizing field.

One then considers a reference frame rotating about $\bm{H}_{0}$ at an angular frequency $\omega_{a}$ equal to the frequency of the applied field. In this frame, the resulting effective field $\bm{H}_{e}$ is also time-independent. The evolution of the magnetization in this frame, neglecting dissipative effects, is simply a precession about the field at the Larmor frequency $\bm{\omega}_{e}=-\gamma\bm{H}_{e}$ due to the torque $\gamma\bm{M}\times\bm{H}_{e}$ on $\bm{M}$, with $\bm{H}_{e}=(H_{a}\cos\phi,H_{a}\sin\phi,H_{0}-\omega_{a}/\gamma)$. The phase $\phi$ relative to the $x$-axis in the rotating frame is arbitrary in the context of a single applied field and can be set equal to zero. However, the relative phase is required for problems involving sequentially applied fields.

Relaxation rates $R_{i}$ are then assigned to each component $M_{i}$ to include dissipative processes. The torque can be written as a matrix-vector product (Jaynes***), which, together with relaxation, gives the matrix

[TABLE]

comprised of the rates and the components of $\bm{\omega}_{e}$. Including the initial polarization $M_{0}$ then gives the Bloch equation in the rotating frame as

[TABLE]

In what follows, both $\bm{H}_{e}$ and $\bm{\omega}_{e}$ are referred to as fields, and we further define the transverse field $\bm{\omega}_{12}$ as a component of the total field $\bm{\omega}_{e}$, with respective magnitudes (squared)

[TABLE]

Framing the problem most generally to include arbitrary $\bm{\omega}_{e}$ and distinct relaxation rates might be expected to complicate the solution compared to previous treatments. However, the symmetry introduced provides concise expressions for \Gamma_{ij}=\varepsilon_{ijk}\,\mbox{\omega_{k}} in terms of the usual Levi-Civita tensor (summed over repeated indices) and $\Gamma_{ii}=R_{i}$, resulting in a compact and rather simple solution.

II.1 Bloch equation solution

Multiplying Eq. (2) by $e^{\Gamma t}$ and integrating the resulting exact differential over the interval $[0,t\,]$ gives the solution

[TABLE]

as a function of the steady-state $\bm{M}_{\infty}\equiv\Gamma^{-1}\bm{M}_{0}R_{3}$ and transient $\bm{M}(0)$ responses. The crux of the problem is a solution for the propagator $e^{-\Gamma t}$.

II.2 The propagator \bm{e^{-\Gamma\,t}}$

The Laplace transform $\mathcal{L}$ of $e^{-at}$ is equal to ($s+a)^{-1}$ for constant $a$. The matrix exponential $e^{-\Gamma t}$ for constant $\Gamma$ is then the inverse Laplace transform \mathcal{L}^{-1}\,[\,(s\mbox{\openone}+\Gamma)^{-1}\,], where $\openone$ is the identity element. The inverse Laplace transform of a function $f(s)$ can be written in terms of the Bromwich integral as (cf. Arfken (1970))

[TABLE]

where the real constant $\gamma$ is chosen such that Re$\,(s)<\gamma$ for all singularities of $f(s)$. Closing the contour by an infinite semicircle in the left half plane ensures convergence of the integral for $t>0$. The desired $F(t)$ is then the sum of the residues of the integrand.

For f(s)=(s\mbox{\openone}+\Gamma)^{-1}, recall the textbook theorem for the inverse of a matrix $A$, with terms defined as follows:

• (i)

$A(i|j)$ is the matrix obtained by deleting row $i$ and column $j$ of $A$.

• (ii)

The cofactor of $A_{ij}$ is $C_{ij}=(-1)^{i+j}\det A(i|j)$.

• (iii)

The adjugate of $A$ is the matrix $(\mathrm{adj}\,A)_{ij}=C_{ji}$, i.e., the transpose of the cofactor matrix for $A$, which is the same as the cofactors of $A$ transpose.

• (iv)

Then $A^{-1}=\mathrm{adj}\,A/\det A$.

For

[TABLE]

the elements of $\mathrm{adj}\,A(s)$, are simple $(2\times 2)$ determinants, giving

[TABLE]

a polynomial in $s$ with coefficient matrices

[TABLE]

Thus, each element of $\mathrm{adj}\,A$ can be constructed from one of two relatively simple elementary forms, with diagonal or off-diagonal elements, respectively, related by cyclic permutation of indices. The matrices for the results above can be found in Appendix D.

The determinant of $A(s)$ is the characteristic polynomial of $(-\Gamma)$. The resulting third degree polynomial $p\,(s)$ is of the form

[TABLE]

with coefficients

[TABLE]

One therefore has, simply,

[TABLE]

The desired solution for $F(t)=e^{-\Gamma t}$ is then the sum of the residues of the integrand in Eq. (5),

[TABLE]

Recall for reference in what follows that for a function $g(s)$ with a pole of order $k$ at $s=s_{0}$, the coefficient of $(s-s_{0})^{-1}$ in the Laurent series expansion of $g(s)$ about $s=s_{0}$, i.e., the residue at $s_{0}$, is

[TABLE]

The poles clearly occur at the roots of $p(s)$, i.e., the eigenvalues of $-\Gamma$.

II.3 Steady state solution

The steady state response $\bm{M}_{\infty}$ defined in Eq. (4) is equal to $\Gamma^{-1}\bm{M}_{0}R_{3}$, with $\Gamma^{-1}$ obtained from Eq. (11) as $\textrm{adj}A(0)/p(0)$. The dependence on $\textrm{adj}A$ is only in the third column, since $\bm{M}_{0}$ is along $\bm{\hat{z}}$, with $p(0)$ given by $c_{0}$ in Eq. (10). Then

[TABLE]

Letting $R_{1}=R_{2}=1/T_{2}$ and $R_{3}=1/T_{1}$ gives

[TABLE]

which reduces to Bloch’s result Bloch (1946), obtained for $\omega_{2}=0$.

III Solutions for the Propagator

The results obtained so far provide the basis for a complete, compact, general solution of the Bloch equation. The solution is valid for all values of the system parameters. Degenerate roots of the characteristic polynomial, which give rise to division by zero in previous solutions, are fully addressed. In a subsequent section, the solution provided in the most general form for the case of three unequal relaxation rates is applied to the more typical case $R_{1}=R_{2}$ for isotropic media.

III.1 Roots of the characteristic polynomial

The solution for $e^{-\Gamma\,t}$ given in Eq. (12) requires the roots of the characteristic polynomial $p(s)$. As is well known, the substitution $s=z-c_{2}/3$ reduces Eq. (9) to the standard canonical form

[TABLE]

where

[TABLE]

Solutions for the roots $z_{i}$ are then available as functions of $a$ and $b$ from standard formulas. However, the roots are complicated functions of the polynomial coefficients (and hence, the physical parameters in the Bloch equation), which hinders physical insight. In Appendix A, simpler expressions are derived for the roots that reduce their complexity compared to previous treatments. The fundamental results are summarized below.

Any polynomial with real coefficients has at least one real root, assigned here to $z_{1}$. The solutions can then be consolidated in a convenient form that does not appear to have been employed before. The other two roots are written as a function of $z_{1}$,

[TABLE]

in terms of a discriminant

[TABLE]

which will be positive, negative, or zero depending on the value of $z_{1}$, the sign of $a$, and their relative magnitudes.

The roots are further characterized here in terms of the positive parameter

[TABLE]

leading to the following delineation of the roots:

• (i)

$a>0$ or $a<0$ and $\gamma>1$

• 3 distinct roots (1 real, 2 complex conjugate)
• (ii)

$a<0$ and $\gamma<1$

• 3 distinct real roots
• (iii)

$a<0$ and $\gamma=1$

• 2-fold degenerate roots $z_{+}=z_{-}=-\tfrac{1}{2}z_{1}$
• (iv)

$a=0=b$

• 3-fold degenerate roots $z_{i}=0$

The physical parameters that define these effective domains for the roots are derived for the standard case $R_{1}=R_{2}$ in a subsequent section.

In addition, the sign of $b$ determines the sign of $z_{1}$. Thus, in all cases, the set of three roots for a given $b<0$ is equal and opposite to the set obtained for parameters that flip the sign of $b$. The case $b=0$ (i.e., $\gamma=0$) reduces simply to $z_{1}\sim\mbox{sgn}(0)=0$. From Eqs. (26) and (27), there are then two additional real or imaginary roots depending on the sign of $\varpi^{2}$.

The roots of $p\,(s)$ needed for the solution of $e^{-\Gamma\,t}$ in Eq. (12) are then

[TABLE]

where, referring to Eq. (10),

[TABLE]

is the average of the relaxation rates.

III.2 A convenient matrix partitioning

We first seek to avoid transforming the characteristic polynomial to canonical form, solving for these roots, then transforming back to obtain the roots of the original polynomial. The result of this endeavor leads to additional simplifications in what follows.

Partition $\Gamma$ as the sum of commuting matrices

[TABLE]

where the diagonal elements of $\Gamma_{\mathrm{p}}$ are

[TABLE]

This partitioning gives $c_{2}=\sum_{i}R_{ip}=0$. Therefore, the characteristic polynomial for $-\Gamma_{\mathrm{p}}$ is in the standard canonical form $q(z)$ of Eq. (24) with coefficients

[TABLE]

from Eq. (10). We then have

[TABLE]

The focus henceforth will be the solution for $e^{-\Gamma_{\mathrm{p}}\,t}$, obtained from Eq. (12) with roots $s_{i}=z_{i}$ given in Eq. (129) and $\mathrm{adj}\,A\rightarrow\mathrm{adj\,}A_{\mathrm{p}}$ obtained from $\Gamma_{\mathrm{p}}$. Using Eqs.(7) and (LABEL:adjA_coeff) with relaxation rates $R_{i}\rightarrow R_{i\mathrm{p}}$ gives

[TABLE]

The partitioning also allows the coefficient matrices to be found rather easily in terms of $\Gamma_{p}$ as

[TABLE]

The result for $A_{1\mathrm{p}}$ is clear by inspection of the off-diagonal elements and confirmed using $\sum_{i}R_{ip}=0$ for the diagonal elements. One expands the $R_{i\mathrm{p}}$ according to Eq. (35) to obtain the expression for $A_{0\mathrm{p}}$. More generally, as can be verified by direct substitution,

[TABLE]

The derivation is fairly straightforward and readily generalized to higher dimensional matrices. However, these topics exceed the scope of the present work.

Either Eq. (39) or (40) can be used to obtain a concise implementation of the Cayley-Hamilton theorem, which states that every square matrix is a solution to its characteristic equation. Thus, $\Gamma_{\mathrm{p}}$ is a solution of Eq. (24). One can solve for $\Gamma_{\mathrm{p}}^{3}$, and subsequently for all higher powers of $\Gamma_{\mathrm{p}}$, in terms of the set \{\mbox{\openone},\Gamma_{\mathrm{p}},\Gamma_{\mathrm{p}}^{2}\}. The series expansion of $e^{-\Gamma_{\mathrm{p}}\,t}$ can then be expressed in terms of the same set. The net result is a relatively simple means for obtaining the scalar coefficients $a_{i}(t)$ in a solution of the form

[TABLE]

III.3 Simple pole solution

In the case that the roots $z_{i}$ of $q\,(z)$ are distinct, the residues are due to simple first-order poles. Factoring $q\,(z)$ as $\prod_{i}(z-z_{i})$ gives $(z-z_{n})/q\,(z)=\prod_{i\neq n}(z-z_{i})$, as needed to evaluate the residue of a first order pole at $z=z_{n}$. The derivative $q^{\prime}(z)=\sum_{j}\prod_{i\neq j}(z-z_{i})$ evaluated at $z_{n}$ is also equal to $\prod_{i\neq n}(z_{n}-z_{i})$, since the other terms in the sum vanish at $z=z_{n}$.

The matrix exponential $e^{-\Gamma_{\mathrm{p}}t}$ can then be written simply and succintly as

[TABLE]

There are thus independent contributions from each root, with shared dependencies on the fields $\omega_{i}$ and relaxation rates $R_{ip}$ which comprise $\mathrm{adj}\,A_{\mathrm{p}}$ according to Eq. (LABEL:adjA_coeff).

Performing the sum, with substitutions from Eqs. (LABEL:adjAp_Poly) and (39), and collecting terms according to Eq. (41) gives a result that can be written compactly as matrix multiplication in the form

[TABLE]

The derivative of the characteristic polynomial can be calculated from either the factored form involving the roots or the polynomial form in Eq. (9). Each provides information that might be useful for different applications.

For parameter values

• (i)

$a>0$ or $a<0$ and $\gamma>1$,

$\varpi$ is real from Eqs. (129) and (129), so two of the roots are complex conjugates. Although Eq. (42) or (54) is the most straightforward form of the solution and readily used in numerical calculations, the individual terms are complex. A more transparently real-valued expression is obtained by performing the sum in Eq. (42) after rationalizing complex denominators and writing the roots $z_{2,3}$ in terms of $z_{1}$ using Eqs. (26) and (27), as detailed in Appendix B. The result is Eq. (54) with

[TABLE]

The coefficient $a$ can be found in terms of the roots $z_{i}$ upon expanding the factored form for $q(z)$ to obtain $a=z_{1}z_{2}+z_{1}z_{3}+z_{2}z_{3}$. The solution for the matrix exponential is thus separable into a term the depends directly on the physical parameters of the problem through $\Gamma_{p}$, a term that depends on the roots $z_{i}$, and a term that gives the time dependence, which in turn is solely a function of the roots.

For the case

• (ii)

$a<0$ and $\gamma<1$,

$\varpi$ is imaginary, as given by Eq. (129), so there are three real roots. There is no oscillatory behavior in the straightforward result given in Eq. (54). The solution can written alternatively in terms of $\mu=|\varpi|$ using Eq. (62), with $\varpi=i\mu$ giving $\cos\varpi\,t\rightarrow\cosh\mu\,t$ and $\sin\varpi\,t/\varpi\rightarrow\sinh\mu\,t/\mu$.

III.4 Second-order pole solution

For the case

• (iii)

$a<0$ and $\gamma=1$,

two of the three real roots are equal, giving a doubly degenerate root $z_{2}=z_{3}=-z_{1}/2$, since $\varpi=0$ when $\gamma=1$ in either Eq. (129) or Eq. (129). The characteristic polynomial $q\,(z)\rightarrow(z-z_{1})(z-z_{2})^{2}$. The contribution from the first-order pole at $z_{1}$ is obtained as before from the $i=1$ term of Eq. (42). The residue at $z_{2}$ is calculated in Appendix B, leading to a solution that can be written in the form

[TABLE]

There is thus a term linear in the time, $t$. Note that Eq. (72) is also the limit of Eq. (62) as $\varpi\rightarrow 0$, which requires $a\rightarrow-3(z_{1}/2)^{2}$ according to Eq. (27), providing an independent verification of the simple-pole result.

III.5 Third-order pole solution

The case

• (iv)

$a=0=b$

gives a triply degenerate, real root $z_{1}=0$ for $q\,(z)\rightarrow z^{3}$. The residue is one-half the second derivative of $\mathrm{adj}\,A_{\mathrm{p}}(z)\,e^{zt}$ with respect to $z$, evaluated at $z=0$, giving

[TABLE]

There is now a term that is quadratic in the time. The same result is obtained from Eq. (72) in the limit $z_{1}\rightarrow 0$ upon series expansion of the exponential terms.

In addition, the Cayley-Hamilton theorem is simple to apply directly in this case, since $q(\Gamma_{\mathrm{p}})=0=\Gamma_{\mathrm{p}}^{3}$. The series expansion of $e^{-\Gamma_{\mathrm{p}}\,t}$ is therefore truncated, giving the Eq. (73) result directly and verifying the self-consistency of the solutions.

As mentioned at the beginning of the section, the solutions can be further simplified when $R_{1}=R_{2}$ to provide increased insight into the nature of the solutions and the constraints that determine root multiplicities.

IV Characterization of the Solutions

Substituting $R_{1}=R_{2}$ gives rates $R_{ip}$ which can be written in the simple form

[TABLE]

where

[TABLE]

since the transverse relaxation rate $R_{2}$ is greater than or equal to the longitudinal rate $R_{3}$ in physical systems. The coefficients of the characteristic polynomial for $-\Gamma_{\mathrm{p}}$ then simplify to

[TABLE]

The rate $R_{\delta}$ provides a convenient and simplifying frequency scale for characterizing the solutions in the sections which follow.

IV.1 Criteria for the existence of degenerate roots

The resulting simpler form for the polynomial coefficients makes possible a straightforward analysis of the conditions for which there are degeneracies in the roots. As discussed in section III.1, there is a two-fold degeneracy in the roots for $a<0$ and $\gamma=1$. This is equivalent, using Eq. (28) for $\gamma$, to

[TABLE]

The trivial solution $a=0=b$ gives a three-fold degenerate root $z_{i}=0$.

Details are deferred to Appendix C, where the existence of degenerate roots is characterized in terms of

[TABLE]

For each $\omega_{3}$ defined by the range $0\leq\lambda_{3}\leq 1$, one finds two solutions for $\lambda_{12}$ that satisfy $D(a,b)=0$ and give real values for $\omega_{12}$. Thus, for each $\omega_{3}\in[\,0,R_{\delta}^{2}/3\,]$, there are two values of $\omega_{12}$ that produce degeneracies in the roots $z_{i}$. The two solutions for $\lambda_{12}$ can be expressed concisely in the form

[TABLE]

for $\lambda_{b}=\tfrac{3}{4}(\sqrt{3}-\tfrac{5}{3})$. The solutions converge at $\lambda_{3}=1$ to $\eta_{1}=\eta_{2}=27/4$, giving $\omega_{12}^{2}=8(R_{\delta}^{2}/3)$. Then $a=0=b$ from Eq. (36), giving a three-fold degenerate root $z_{i}=0$ of Eq. (24).

The following simple and explicit criteria define the characteristics of the roots:

• (i)

$\omega_{3}^{2}>R_{\delta}^{2}/3$

there is no real-valued solution for $\omega_{12}^{2}$ such that $\gamma=1$, i.e., $D(a,b)=0$, and hence no degenerate roots $z_{i}$. One then has the simple-pole solution of Eq. (62).

• (ii)

$\omega_{3}^{2}<R_{\delta}^{2}/3$

there are two different real-valued solutions for $\omega_{12}^{2}$ as a function of $\lambda_{3}$ that each give a two-fold degeneracy in the roots $z_{i}$, requiring the second-order pole solution of Eq. (72).

• (iii)

$\omega_{3}^{2}=R_{\delta}^{2}/3$

gives $\omega_{12}^{2}=8(R_{\delta}^{2}/3)$ for $\lambda_{3}=1$, resulting in a three-fold degenerate root $z_{i}=0$ which requires the third-order pole solution of Eq. (73).

IV.2 Characterization of the damping

Solutions for the roots $z_{i}$ are characterized according to whether the discriminant $\varpi^{2}$ of Eq. (27) is positive, negative, or zero, and can be described, respectively, as underdamped, overdamped, or critically damped, analogous to a damped harmonic oscillator.

The solution for the propagator in the case of degenerate roots ($a<0,\gamma=1$) has a term linear in time, characteristic of a critically damped harmonic oscillator. For a three-fold degeneracy in the roots, there is an additional term that is quadratic in the time. The allowed values of $\omega_{3}^{2}$, as discussed in the previous section, are restricted to the narrow range parameterized according to $0\leq\lambda_{3}\leq 1$. The two solutions $\omega_{12,1}^{2}$ and $\omega_{12,2}^{2}$ for each $\omega_{3}^{2}$, as determined from Eqs. (78) and (79), are the solid curves plotted in Fig. 1.

Using the same scaling of $\omega_{3}$ and $\omega_{12}$ as in Eq. (78), we also have

[TABLE]

Solutions in the range $\omega_{12,1}^{2}<\omega_{12}^{2}<\omega_{12,2}^{2}$ bounded by the critical damping parameters give $a<0$ and $\gamma<1$, resulting in three distinct real roots and overdamped evolution. The range of bounding values is fairly narrow, becoming increasingly so with increasing $\lambda_{3}$ and converging to a single value $\omega_{12}^{2}=8R_{\delta}^{2}/3$ as $\lambda_{3}\rightarrow 1$, as shown in the figure.

Underdamped, oscillatory solutions are obtained for all other field values, either $\omega_{3}^{2}>R_{\delta}^{2}/3$ (i.e., $\lambda_{3}>1$) or $\omega_{12}^{2}\geq\omega_{12,1}^{2}$ and $\omega_{12}^{2}\leq\omega_{12,2}^{2}$ for $\lambda_{3}\leq 1$.

IV.3 Characterization of the roots

The solution to the Bloch equation has a relatively simple form and can be expressed in terms of a single root, $z_{1}$, of the characteristic polynomial for $-\Gamma_{p}$. Although the solutions for $z_{1}$ have also been expressed in relatively simple functional form, these forms provide little physical insight. It remains to shed some light on the dependence of this root on the field $\bm{\omega}_{e}$ and the relaxation rates.

IV.3.1 Physical limits of the roots

The roots $z_{i}$, being functions of $a,b$ and $\gamma$, also scale as $R_{\delta}$. The associated decay rates are $\mathrm{Re}(s_{i})=\mathrm{Re}(z_{i})-\bar{R}$, from Eq. (29). Defining

[TABLE]

and using Eq. (75) for $R_{\delta}$ gives the decay rates

[TABLE]

The limiting rates are $R_{2}$ and $R_{3}$, which therefore constrains $\lambda_{z}$ to the range

[TABLE]

The damping has equal contributions from $R_{2}$ and $R_{3}$ for $\lambda_{z}=1/2$, with a larger contribution from either $R_{2}$ or $R_{3}$ if $\lambda_{z}$ is less than or greater than 1/2, respectively.

The dependence of $z_{1}$ on $\bm{\omega}_{e}$ and $R_{\delta}$, calculated according to Eqs. (129), is shown in Fig. 2, where contours of $\lambda_{z}$ are plotted as a function of $\lambda_{12}$ and $\lambda_{3}$. As discussed earlier, there is only one real root for $\lambda_{3}>1$. When $\lambda_{3}\leq 1$, there is also a single real root for values of $\lambda_{12}$ outside the narrow bounds that define critical damping. Within these bounds where the solutions represent overdamping, any of the three real roots can be designated as $z_{1}$, with $z_{\pm}$ from Eq. (129) giving the other two. For $\omega_{12}=0$, the relaxation rate is $R_{3}$ (i.e., $\lambda_{z}=2$), independent of the offset parameter $\lambda_{3}$, as is well-known. As $\omega_{12}$ increases for fixed $\omega_{3}$, the relaxation rate approaches $R_{2}$ ($\lambda_{z}=-1$), with the drop-off from $\lambda_{z}=2$ becoming increasingly steep at lower values of $\omega_{3}$. For the other roots in which $\mathrm{Re}(z_{\pm})=-1/2\,z_{1}$, the upper limit in Eq. (83) becomes 1/2.

IV.3.2 A linear relation for the roots

Equation (24) evaluated at the real root $z_{1}$ yields the linear relation

[TABLE]

for coefficient $b$ that will satisfy Eq. (24) as a function of a given coefficient $a$, with slope and intercept determined by $z_{1}$. Substituting the expressions for $a$ and $b$ in Eq. (80), rearranging and collecting terms after writing $9\lambda_{z}=6\lambda_{z}+3\lambda_{z}$ gives

[TABLE]

with slope $m_{\mathrm{s}}$ and intercept $\lambda_{12}^{\mathrm{int}}$

[TABLE]

There is thus a simple graphical representation for the value of the root $z_{1}$ as a function of the physical parameters $\,\omega_{12},\omega_{3},R_{\delta}$. There are a continuum of field values for a given $R_{\delta}$ that give the same $z_{1}$. Lines of constant $z_{1}$ as a function of $\lambda_{12}$ and $\lambda_{3}$ become hyperbolas when Eq. (85) is rewritten in terms of $\omega_{12}^{2},\omega_{3}^{2},R_{\delta}^{2}$ using Eq. (78).

V Intuitive Representations of System Dynamics

In most cases, the parameters of the Bloch equation yield three distinct roots for the characteristic polynomial $p(s)$ of Eq. (9), described as cases (i) and (ii) in section III.1. Exceptions were considered in more detail in section IV for the condition $R_{1}=R_{2}$. To provide additional physical insight, we develop a straightforward vector model of the time evolution for $M(t)$ given in Eq. (4). This requires the action of the propagator $e^{-\Gamma\,t}$ on an arbitrary vector. An alternative vector model is also considered, followed by a coupled oscillator model.

The eigensystem for $\Gamma$ is considered in sections that follow, but one can substitute notation for the partitioned matrix $\Gamma_{\mathrm{p}}$ in the expressions which are derived, since, as defined in Eq. (34), the matrices differ by a constant $\bar{R}$ times the identity matrix. The difference in the eigenvalues is also $\bar{R}$, from Eqs. (29) and (30). Thus $-\Gamma$ and $-\Gamma_{\mathrm{p}}$ have the same eigenvectors $\bm{s}_{i}\equiv\bm{z}_{i}$. Simple analytical expressions for the eigenvectors and other constituents of the model are derived in Appendix D. Each (unnormalized) eigenvector, which can assume different analytical forms depending on the scaling, is found to comprise the columns of $\mathrm{adj}\,A(s_{i})=\mathrm{adj}\,A_{\mathrm{p}}(z_{i})$, as discussed in Appendix E, providing an alternative method for calculating an eigenvector.

V.1 Existing models specific to simple limiting cases

As a point of departure, consider first the simple limiting cases for which the dynamics is already well known and readily visualized. In the absence of relaxation, i.e., all $R_{i}=0$, any magnetization vector $\mathcal{M}$ rotates about the total effective field $\bm{\omega}_{e}$ at constant angular frequency $\omega_{e}$. The time evolution of a vector under the action of the propagator has a simple solution in a coordinate system rotated to align one of the axes with the effective field. The component of $\mathcal{M}$ along $\bm{\omega}_{e}$ is constant, and the components in the plane perpendicular to $\bm{\omega}_{e}$ rotate at angular frequency $\omega_{e}$ in the plane. By constrast, the solution for each component $\mathcal{M}_{i}(t)$ in the standard $(x_{1},x_{2},x_{3})$-coordinate system is more complicated, and it is not immediately apparent by inspection that the solution is a rotation.

If the relaxation is switched on with equal rates $R_{i}=R$ on the diagonal, the relaxation matrix R\mbox{\openone} commutes with the remaining rotation matrix, and the solution is a dynamic scaling $e^{-Rt}$ of the rotating vector $\mathcal{M}$. In addition, for $\omega_{12}=0$ and $R_{1}=R_{2}\neq R_{3}$, the relaxation matrix still commutes with the rotation about nonzero $\omega_{3}$. The evolution is then a scaling $e^{-R_{2}t}$ of the transverse component $\mathcal{M}_{12}$, which rotates at angular frequency $\omega_{3}$ in the plane perpendicular to $\omega_{3}$, along with exponential decay $e^{-R_{3}t}$ of component $\mathcal{M}_{3}$, as illustrated in Fig. 3a. In the case of pure relaxation, with all the field components $\omega_{i}=0$, the solution is a non-oscillatory exponential decay $e^{-R_{i}t}$ for each component $\mathcal{M}_{i}$ along coordinate axis $x_{i}$.

V.2 A more general model

With the exception of the above simple cases, there has been no analogous picture of system dynamics when the rotation and relaxation do not commute. The combined, noncommutative action of arbitrary fields and dissipation rates appears to require something more complex. Yet, the simple visual model shown in Fig. 3a, which is comprised of independent relaxation and rotation elements, is readily extended to the general case when viewed in an appropriate coordinate system.

V.2.1 One real, two complex conjugate roots

The solution for each component $\mathcal{M}_{i}$ is known to be a combination of oscillation and bi-exponential decay Torrey (1949), as is also evident from the propagator derived in Eq. (12). The underlying simplicity of the system dynamics can be demonstrated starting with the eigensystem for $\Gamma$ (or, alternatively, $\Gamma_{\mathrm{p}}$, as noted above).

The real eigenvalue $s_{1}$ of $-\Gamma$ has a real eigenvector $\bm{s}_{1}$ which can be used as one axis of a physical coordinate system, but the complex roots $s_{+}$ and $s_{-}=s_{+}^{*}$ have associated complex eigenvectors $\bm{s}_{+}$ and $\bm{s}_{-}=\bm{s}_{+}^{*}$. The eigenvectors are most generally not orthogonal, but they are linearly independent, given the distinct eigenvalues.

Define the real vectors

[TABLE]

The set $\{\bm{\tilde{s}}_{1},\bm{\tilde{s}}_{2},\bm{\tilde{s}}_{3}\}$ can then be used as an alternative basis for describing the system evolution. System states and operators are transformed between bases in the usual fashion by a matrix $P$ comprised of the $\{\bm{\tilde{s}}_{i}\}$, entered as column vectors. Vector $\tilde{\mathcal{M}}$ and matrix $\tilde{\Gamma}$ in the new basis are given by

[TABLE]

with $P$ invertible since the $\bm{\tilde{s}}_{i}$ are linearly independent.

The potentially tedious process of calculating $e^{-\tilde{\Gamma}\,t}$ from Eq. (88) can be bypassed, with $e^{-\tilde{\Gamma}\,t}$ deduced from the action of $\Gamma$ on its eigenvectors (see Appendix D). In terms of constants

[TABLE]

the solution $\tilde{\mathcal{M}}(t)=e^{-\tilde{\Gamma}t}\tilde{\mathcal{M}}(0)$ for the time dependence of state vector $\tilde{\mathcal{M}}$ in the new basis is found to be

[TABLE]

Viewed in the $\{\bm{\tilde{s}}_{i}\}$ coordinate system, the component of $\mathcal{M}$ along $\bm{\tilde{s}}_{1}$ (i.e., $\tilde{\mathcal{M}}_{1}$) decays at the rate $\bar{R}-z_{1}$, while components in the $(\bm{\tilde{s}}_{2},\bm{\tilde{s}}_{3})$-plane rotate in the plane and decay at the rate $\bar{R}+z_{1}/2$. Thus, even in the most general case of three unequal rates $R_{1},R_{2},R_{3}$, there emerges a single “planar” relaxation rate $R_{2s}$ and a new “longitudinal” relaxation rate $R_{1s}$ defined as

[TABLE]

Defining $\tilde{\mathcal{M}}(t)$ as the state $\bm{M}(t)-\bm{M}_{\infty}$ expressed in the $\{\bm{\tilde{s}}_{i}\}$ coordinates and working backwards from Eq. (97) gives the Bloch equation in this basis as

[TABLE]

The diagonal matrix consisting of the relaxation rates $R_{is}$ commutes with the matrix of off-diagonal elements, which generates a rotation about $\bm{\tilde{s}}_{1}$, and one immediately obtains the solution given in Eq. (97).

One therefore has considerable latitude in the choice of $\bm{\tilde{s}}_{2}$ and $\bm{\tilde{s}}_{3}$, since all components in the plane they define decay at the same rate. Rotating these coordinate axes in the plane by any angle results in an equally valid set of axes for representing the dynamics. The vectors $\bm{\tilde{s}}_{2}$ and $\bm{\tilde{s}}_{3}$ constructed from a particular column in the coefficient matrices of Eq. (154) are related by such a rotation to the axes constructed from one of the other columns (excepting when one of the columns returns the irrelevant zero vector). By contrast, $\bm{\tilde{s}}_{1}$ defines the unique axis for longitudinal decay, so the $\bm{\tilde{s}}_{1}$ chosen from different columns must be related by a scale factor.

Note also that the rotation in the plane is not at a constant angular frequency $\varpi$ unless $\bm{\tilde{s}}_{2}$ and $\bm{\tilde{s}}_{3}$ are orthogonal. A component aligned with $\bm{\tilde{s}}_{2}$ rotates to $\bm{\tilde{s}}_{3}$ during a time defined by the condition $\varpi t=\pi/2$, then rotates from there to $-\bm{\tilde{s}}_{2}$ in the same time. In an oblique coordinate system, the rotations are through different angles in the same time, so clearly the angular frequency of the rotation is not constant.

V.2.2 Three real roots

In this case, all the eigenvectors are real and the new basis is simply the eigenbasis $\{\bm{s}_{1},\bm{s}_{2},\bm{s}_{3}\}$ obtained from the roots

[TABLE]

defined in Eq. (29). The real roots $z_{i}$ are obtained for $\varpi^{2}<0$ in Eq. (27). Substituting $\varpi\rightarrow i\mu$ in Eq. (26) gives $z_{2,3}=-1/2\,z_{1}\mp\mu$.

The matrix $\Gamma$ is obviously diagonal in its eigenbasis, and, by extension, so is the propagator in this basis. Thus

[TABLE]

Each component of $\mathcal{M}$ along $\bm{\tilde{s}}_{i}$ decays at the rate determined by $s_{i}$. In contradistinction to the rates that emerge from the oscillatory solutions, here, even in the typical case of equal transverse rates $R_{1}=R_{2}$ and longitudinal rate $R_{3}$, we find three distinct rates

[TABLE]

due to the coupling of the field with the relaxation processes.

Given $e^{-\tilde{\Gamma}t}$ as obtained in Eq. (97) or (107), the propagator in the standard coordinate basis is $e^{-\Gamma t}=Pe^{-\tilde{\Gamma}t}P^{-1}$ from Eq. (88). One obtains a simple, factored solution for the propagator and a correspondingly simple physical interpretation of the dynamics, with oscillation frequencies and decay rates hinging upon the primary real root $z_{1}$. The dependence of this root on the fields and relaxation rates has been shown previously in Fig. 2.

V.2.3 Degenerate roots

The vector model approach to obtaining the propagator is only applicable to the case of distinct eigenvalues. Degenerate eigenvalues do not give the linearly independent eigenvectors necessary to define a new coordinate system. However, the degeneracies are a relatively trivial component of the parameter space, at least for $R_{1}=R_{2}$, as shown in Fig. 1. Moreover, the solution has to be continuous as the degeneracies are approached, with a smooth transition from oscillatory, decaying solutions to pure decay as one crosses the parameter-space boundary identifying the degenerate solutions.

V.3 Discussion and representative examples

The solutions of section III are represented in the standard coordinate system, expressed in general form for the case of three unequal relaxation rates. Here, they are applied to specific physical examples, with $R_{1}=R_{2}$. The trajectories of initial states under the action of the propagator are plotted to illustrate the underlying simplicity of the dynamics and corroborate the alternative coordinate system that defines the vector model. Parameters for the examples are chosen to demonstrate the damping and rotation that are characteristic of the dynamics for all but a small region of the parameter space. A purely damped solution and model dynamics given by Eq. (107) is rather featureless, by comparison. Unless stated otherwise, the first column of $\mathrm{adj}\,A_{\mathrm{p}}$ is chosen to calculate the coordinate basis $\{\bm{\tilde{s}}_{i}\}$.

V.3.1 Free precession, $\bm{\omega}_{e}=(0,0,\omega_{3})$

When the only field in the rotating frame is the offset from resonance, $\omega_{3}$, the matrix $\Gamma_{\mathrm{p}}$ is the sum of a diagonal relaxation matrix and the matrix which generates a rotation about $\omega_{3}$. Since they commute, the propagator factors into the product of exponential decay and a rotation, leading to the standard interpretation of the dynamics discussed previously. This example also provides a simple context for applying the more general vector model. The eigenvalues are easily obtained as $z_{1}=2R_{\delta}$ and $z_{\pm}=-R_{\delta}\pm i\omega_{3}$. Then Eq. (146) gives, upon identifying $\varpi\equiv\omega_{3}$ and eliminating common factors in individual columns,

[TABLE]

As noted earlier, there is always only one unique nonzero result for $\bm{\tilde{s}}_{1}$, with any apparent differences between columns simply a matter of scale. The nonzero columns for $\bm{\tilde{s}}_{2}$ are orthogonal, as are those of $\bm{\tilde{s}}_{3}$. The columns thus differ, as expected, by a rotation in the $(\bm{\tilde{s}}_{2},\bm{\tilde{s}}_{3})$-plane, in this case by $90^{\circ}$. Choosing the second column and a left-handed rotation by $\phi=\tan^{-1}(3R_{\delta}/\omega_{3})$ or the first column and a right-handed rotation by $90-\phi$ gives the more typical result $\bm{\tilde{s}}_{2}=(0,1,0)$ and $\bm{\tilde{s}}_{3}=(1,0,0)$ depicted in Fig. 3a. The model dynamics for an initial state $\mathcal{M}_{0}$ is a spiral about $\bm{\omega}_{e}$, which is aligned along the $z$-axis, with rotation at constant angular frequency $\omega_{e}$ in the $(x,y)$-plane, as required. The relaxation rate obtained from Eq. (82) or Eq. (89) for $z_{1}=2R_{\delta}$, with $\lambda_{z}=2$, is $R_{1s}=R_{3}$, while the roots $z_{\pm}$ with $\lambda_{z}=-1$ give $R_{2s}=R_{2}$.

V.3.2 On resonance, $\bm{\omega}_{e}=(\omega_{1},\omega_{2},0)$

On resonance, the root $z_{1}=-R_{\delta}$, and $\varpi^{2}=\omega_{e}^{2}-(3/2R_{\delta})^{2}$ from Eq. (189). The associated eigenvector $\bm{\tilde{s}}_{1}$ is obtained by inspection from Eq. (145), with $\bm{\tilde{s}}_{2}$ and $\bm{\tilde{s}}_{3}$ obtained from Eqs. (146) and (154), giving

[TABLE]

Thus, on resonance, the propagator still generates a spiral about the effective field $\bm{\omega}_{e}=\bm{\tilde{s}}_{1}$ with precession in the $(\bm{\tilde{s}}_{2},\bm{\tilde{s}}_{3})$-plane orthogonal to $\bm{\tilde{s}}_{1}$. However, as considered in section V.2.1, the rotation frequency driven by $\varpi$ is not constant, since $\bm{\tilde{s}}_{2}$ is not perpendicular to $\bm{\tilde{s}}_{3}$. The deviation from orthogonality, determined by the third component of $\bm{\tilde{s}}_{2}$, is small for fields that are large compared to $R_{\delta}$. The respective decay rates $R_{1s}$ and $R_{2s}$ are $R_{2}$ and $1/2(R_{2}+R_{3})$, using $\lambda_{z}=-1$ and $\lambda_{z}=1/2$ as determined from $z_{1}$ and $-z_{1}/2$. Components along $\bm{\tilde{s}}_{1}$, i.e., in the $(x,y)$-plane, decay at the usual spin-spin relaxation rate, as would be expected. Components rotating in the plane orthogonal to $\bm{\tilde{s}}_{1}$ experience equal influence, on average, from their projection onto the longitudinal $z$-axis defining $\omega_{3}$ and their projection into the $(x,y)$-plane, so one might predict from the model that they decay at the average of the usual spin-spin and longitudinal relaxation rates. These values for the decay rates have been obtained previously as elements of the solution in the standard coordinate system Torrey (1949) without the physical interpretation presented here.

The trajectory for an initial state ${\cal M}_{0}$ due to the action of propagator $e^{-\Gamma t}$ with $\bm{\omega}_{e}=(\omega_{1},0,0)$ and nonzero relaxation is shown in Fig. 3b. Values of the parameters are given in the caption. For nonzero $\omega_{2}$, the figure is simply rotated about the $z$-axis by angle $\phi=\tan^{-1}(\omega_{2}/\omega_{1})$. The state ${\cal M}_{0}$ has been chosen with equal components parallel and orthogonal to $\bm{\omega}_{e}$ to most clearly illustrate the dynamics predicted by the vector model. The slight misalignment between $\bm{\tilde{s}}_{2}$ and the $y$-axis is evident in the figure and becomes more prominent as the magnitude of the field, $\omega_{12}$, is reduced relative to $R_{\delta}$.

V.3.3 Off resonance, general $\bm{\omega}_{e}$

Most generally, $\bm{\tilde{s}}_{1}$ is not aligned with $\bm{\omega}_{e}$. Dividing column $j$ of the matrix in Eq. (145) by (nonzero) $\omega_{j}$ quantifies the degree to which $\bm{\tilde{s}}_{1}$ deviates from $\bm{\omega}_{e}$, due to the coupling between the fields and the relaxation. The result is an expression of the form $\bm{s}_{1}=\bm{\omega}_{e}+\delta\bm{v}$, where vector $\delta\bm{v}$ is comprised of the second term in each row of the $j^{\mathrm{th}}$ column divided by $\omega_{j}$.

In addition, $\bm{\tilde{s}}_{1}$ is typically not orthogonal to the $(\bm{\tilde{s}}_{2},\bm{\tilde{s}}_{3})$-plane. One then has to further modify intuitions developed from orthogonal coordinate systems. For example, in Fig. 3c, ${\cal M}_{0}$ is aligned with the normal to the $(\bm{\tilde{s}}_{2},\bm{\tilde{s}}_{3})$-plane. It therefore has no orthogonal projection in the plane and might naively be expected to have no evolution in the plane. However, $\bm{\tilde{s}}_{1}$ is distinctly different than the normal, and ${\cal M}_{0}$ is the vector sum of a component along $\bm{\tilde{s}}_{1}$ and a component parallel to the plane, which are the quantities relevant for the vector model. As shown in the figure, the parallel component rotates and decays in the plane while the component along $\bm{\tilde{s}}_{1}$ strictly decays. Similarly, $\mathcal{M}_{0}$ orthogonal to $\bm{\tilde{s}}_{1}$ as in Fig. 3d nonetheless has a component along $\bm{\tilde{s}}_{1}$ in the oblique coordinates that decays to generate the spiral shown in the figure.

By constrast, the dynamics viewed in standard coordinates is oscillation of each component ${\cal M}_{i}(t)$ combined with relaxation at two separate rates. As in simpler examples, it can be decoupled into two independent dynamical systems, one of which rotates in a plane and decays at one rate and another which decays along a fixed axis, albeit in an oblique coordinate system.

The deviation of $\bm{\tilde{s}}_{1}$ from the normal to the plane is quantified in Appendix D off resonance for $\bm{\omega}_{12}$ of either $x$- or $y$-phase and the case $\omega_{1}=\omega_{2}=\omega_{3}$.

V.4 Alternative vector model

The Bloch equation is typically represented in vector form but can be conveniently packaged in matrix form, which is the approach taken here. The physics of its solution—the torque on a magnetic moment in a magnetic field subject to relaxation of the magnetization—can be made more explicit by returning to the original vector operations, motivated by the treatment in Jaynes (1955) for the rotation of a vector about the field.

Partition $\Gamma_{p}$ into its diagonal elements $R_{ip}$ and off-diagonal $\omega_{i}$, writing $\Gamma_{p}=\mathcal{R}_{\mathrm{p}}+\Omega$. The diagonal matrix $\mathcal{R}_{\mathrm{p}}$ scales each component $\mathcal{M}_{i}$ of a vector $\mathcal{M}$ by $R_{ip}$, and $\Omega$ implements the cross product $(-\bm{\omega}_{e}\times\ )$. According to Eq. (41), the propagator acting on $\mathcal{M}$ generates three separate vectors $\bm{v}_{n}=\Gamma_{p}^{n}\mathcal{M},(n=0,1,2)$, which can be represented starting with $\bm{v}_{0}=\mathcal{M}$ as

[TABLE]

Each succeeding $\bm{v}_{n}$ is a nonuniform scaling of the previous $\bm{v}_{n-1}$ added to a vector ($\bm{v}_{n-1}\times\bm{\omega}_{3}$) that is orthogonal to $\bm{v}_{n-1}$. The time dependence of $\bm{v}_{n}$ is given by the associated term $a_{n}(t)e^{-\bar{R}t}$ found in Eqs. (62–73). The $a_{n}(t)$ are factored as the product of a matrix $W(z_{1})$ and vector $\bm{u}(t)$. Each $a_{n}(t)$ is merely a different linear combination of the same three simple functions $u_{i}(t)$ that comprise the components of $\bm{u}$, weighted according to the corresponding elements from row $n$ of the matrix $W$. A given $\bm{v}_{n}(t)$ thus maintains a fixed orientation, changing length with a time dependence consisting of the different weightings of the $u_{i}(t)$ for different $\bm{v}_{n}$. The trajectory $\mathcal{M}(t)=\sum_{n}\bm{v}_{n}(t)$ can thus be represented in terms of the decaying oscillations of three vectors fixed in place.

Alternatively, expand (\mbox{\openone},\Gamma_{p},\Gamma_{p}^{2})\,W(z_{1})\bm{u}(t) and group terms of the same time dependence $u_{i}(t)$, as, for example, in Eq. (134). The propagator applied to $\mathcal{M}$ gives three different linear combinations of the $\bm{v}_{n}$, with a time dependence $u_{i}(t)$ for the $i^{\mathrm{th}}$ combination. The resulting interpretation of $\mathcal{M}(t)$ is similar to the previous paragraph, but the functional form of the decaying oscillations is simpler using this different set of vectors.

V.5 The Bloch equation as a system of coupled oscillators

Any quantum N-level system can be represented as a system of coupled harmonic oscillators Skinner (2013), albeit requiring negative or even antisymmetric couplings. The Bloch equation is perhaps particularly interesting, since it incorporates dissipation for the most elementary case, i.e., 2-level systems.

Expressing Eq. (2) in terms of $\mathcal{M}(t)=\bm{M}(t)-\bm{M}_{\infty}$ yields a homogeneous first-order differential equation and an alternative route to the solution, Eq. (4b). Differentiating again with respect to time and substituting $\mathcal{\dot{M}}=-\Gamma\mathcal{M}$ gives

[TABLE]

with

[TABLE]

As considered previously Skinner (2013), damping is provided by the antisymmetric part of $\Gamma^{2}$ in addition to the terms $R_{i}^{2}$ on the diagonal. For the system of three coupled oscillators illustrated in Fig. 4, the displacement $r_{i}$ of mass $m_{i}$ from equilibrium is equal to $\mathcal{M}_{i}$. We can write the coupling constants $k_{ij}=\kappa_{ij}+\sigma_{ij}$ in terms of symmetric $\kappa_{ij}$ and anti-symmetric $\sigma_{ij}$ connected in parallel. Then, by inspection, $\kappa_{ij}=\omega_{i}\omega_{j}$ and $\sigma_{ij}=\varepsilon_{ijk}(R_{i}+R_{j})\omega_{k}$, assuming unit masses, and $k_{ii}=-(\Gamma^{2})_{ii}-\sum_{j}k_{ij}$. For a given positive $\sigma_{ij}$, a positive displacement of mass $m_{j}$ results in a positive force on $m_{i}$. The resulting positive displacement of $m_{i}$ provides a negative force on $m_{j}$ due to $\sigma_{ji}<0$ which opposes the original displacement of $m_{j}$ and damps the motion. Stated differently, energy transferred from $m_{j}$ to $m_{i}$ is not reciprocally transferred back from $m_{i}$ to $m_{j}$, and the motion is quenched. An antisymmetric coupling acts as a negative feedback mechanism that curbs system oscillations.

The usual representation of damped oscillators employs a velocity-dependent friction force. The above implementation is frictionless. It provides an alternative model for investigating dissipative processes with the potential for new insights within the well understood context of coupled harmonic oscillations.

VI Conclusion

A complete solution of the Bloch equation has been presented together with intuitive visual models of its dynamics. The solution is valid for arbitrary system parameters, yet is simpler than previous solutions. It can be expressed as the product of three separate terms: one which depends directly on the physical parameters of the problem through the matrix $\Gamma_{p}$, a term that depends on the roots of a cubic characteristic polynomial for the problem, and a term that gives the time dependence, which in turn is solely a function of the roots. Moreover, the time evolution of the system as a function of the physical parameters has been made more explicit and apparent.

The solutions depend critically on the three polynomial roots. Quantitative relations have been derived for the physical parameters that define the possible system dynamics: (i) oscillatory, underdamped evolution for one real and two complex-conjugate roots, (ii) non-oscillatory, overdamped evolution for three real roots, and (iii) non-oscillatory, critically damped evolution for doubly or triply degenerate (real) roots. The damping rates and the frequency driving the oscillatory behavior have been reduced to simple functions of a single root which is obtained as a straightforward function of the system parameters. In addition, a linear relation has been derived for the system parameters as a function of this real root, which provides a straightforward graphical realization of the damping rates and frequency for a given physical configuration.

An intuitive dynamical model developed here transforms the Bloch equation to a frame in which damping commutes with a rotation, providing a propagator for the time evolution of the system that is the product of a rotation times a decay, in either order. The decay rates in this frame result from interaction/coupling of the fields with the spin-lattice and spin-spin relaxation processes. The model was motivated by well-known visual models for simple cases such as equal relaxation rates or free precession (no fields transverse to the longitudinal, $z$-axis). The system state in such cases rotates about the effective field, with concurrent exponential decay of the longitudinal and transverse components. The extended model retains the same essential features: rotation, exponential decay of the invariant component in the rotation analogous to longitudinal relaxation, and a separate decay of the rotating components analogous to transverse relaxation. An alternative vector model has also been provided, as well as a representation of the Bloch equation as a system of coupled harmonic oscillators. The net result of the solutions and models is more direct physical insight into the dynamics of the Bloch equation.

Acknowledgements.

The author gratefully acknowledges support from the National Science Foundation under Grant CHE-1214006.

Appendix A Cubic Polynomials with Real Coefficients

The standard solutions for the three roots of Eq. (24), cast here in terms of

[TABLE]

are

[TABLE]

These solutions can be consolidated in a convenient form that does not appear to have been employed heretofore. Substituting (\Lambda_{+}-\Lambda_{-})=\big{[}(\Lambda_{+}+\Lambda_{-})^{2}-4\Lambda_{+}\Lambda_{-}\big{]}^{1/2} and noting $\Lambda_{+}\Lambda_{-}=-a/3$ gives

[TABLE]

in terms of a discriminant

[TABLE]

Any polynomial with real coefficients has at least one real root. Therefore $\varpi^{2}>0$ gives one real and two complex conjugate roots, with three real roots resulting from $\varpi^{2}\leq 0$.

One can then employ simple forms for $z_{1}$ Miura (1980); McKelvey (1984). The number of conditional dependencies relating the cited expressions for $z_{1}$ to the signs and relative magnitudes of $a$ and $b$ can be further simplified in terms of

[TABLE]

Then the roots can be calculated according to their domain of applicability as

[TABLE]

For $a>0$ or $a<0$ and $\gamma>1$, there is one real root and complex conjugate roots $z_{\pm}$. For ${a<0,\,\gamma<1}$, there are three real roots, with $\gamma=1$ in Eq. (129) or (129) giving $\varphi=0=\varpi$ and two degenerate roots $z_{+}=z_{-}$, while (129) reorders the roots relative to (129), so that the nondegenerate root for the case $\gamma=1$ is one of the $z_{\pm}$. Results for $a=0$ are straightforwardly obtained from Eqs. (125) and (26), or using the expressions in (129) and (129), with $\sinh^{-1}\gamma\rightarrow\cosh^{-1}\gamma\rightarrow\ln(2\gamma)$ in the limit $\gamma\rightarrow\infty$. Terms then result that are multiplied by $\sqrt{\alpha}$, cancelling the singularity at $a=0$. For the case $a=0=b$, there are three equal roots $z_{i}=0$.

Appendix B Calculation of $\bm{e^{-\Gamma_{\mathrm{p}}t}}$

B.1 First-order pole

Consider the case of one real root $z_{1}$ and two complex conjugate roots $z_{2,3}=-1/2z_{1}\pm i\,\varpi$, as given by Eq. (26), with $\varpi^{2}=3(z_{1}/2)^{2}+a>0$. Using Eq. (7) for $\mathrm{adj}A_{\mathrm{p}}$ in Eq. (42) gives

[TABLE]

with $q^{\prime}(z_{i})=\prod_{j\neq i}(z_{i}-z_{j})$, as discussed in section III.3, and A_{2\mathrm{p}}\equiv\mbox{\openone}.

Evaluating the $q^{\,\prime}(z_{i})$ and using Eq. (27) for $\varpi^{2}$ gives

[TABLE]

Since $z_{3}=z_{2}^{*}$, each coefficient $(A_{n})_{\mathrm{p}}$ in Eq. (42) multiplies a sum $S_{n}$

[TABLE]

Substituting $z_{2}=-z_{1}/2+i\,\varpi$ gives

[TABLE]

We then have, using Eq. (39) for the $(A_{n})_{\mathrm{p}}$,

[TABLE]

Arranging coefficients of (\mbox{\openone},-\Gamma_{p},\Gamma_{p}^{2}) in a matrix for each time-dependent term in the solution gives the result in Eq. (62). All three roots are real when $\varpi^{2}<0$, which is the case for $a<0$ and $\gamma<1$. Then $\varpi\rightarrow i\mu$ in Eq. (62), with $\mu^{2}=|3(z_{1}^{2}/2)+a|$ and $a=-|a|$.

B.2 Second-order pole

The case $\varpi=0$ gives doubly-degenerate real roots $z_{2}=z_{3}=-z_{1}/2$, and the characteristic polynomial $q\,(z)\rightarrow(z-z_{1})(z-z_{2})^{2}$. The residue at $s=z_{2}$ in Eq. (12) requires the derivative of $e^{st}\mathrm{adj}\,A_{\mathrm{p}}(s)/(s-z_{1})$ with respect to $s$, evaluated at $s=z_{2}$. Expanding $\mathrm{adj}\,A_{\mathrm{p}}$ using Eqs. (7) and (39) as above, utilizing a common denominator $(z_{2}-z_{1})^{2}$, and substituting $z_{2}=-z_{1}/2$ gives

[TABLE]

The contribution from the first-order pole at $z_{1}$ is obtained as before from the $i=1$ term of Eq. (62) with $a=-3(z_{1}/2)^{2}$, since $\varpi=0$, to yield the result of Eq. (72).

Appendix C Existence of Degenerate Roots

The characteristic polynomial for the case $R_{1}=R_{2}$ has degenerate roots for $D(a,b)=0$ (cf. Eq. (77)), which requires $a<0$. The special case $a=0=b$ discussed in section IV.1 gives $\omega_{3}^{2}=1$ and $\omega_{12}^{2}=8$, normalized to $R_{\delta}^{2}/3$. More generally, scale $\omega_{3}^{2}$ and $\omega_{12}^{2}$ in terms of the same normalization as

[TABLE]

where $\lambda_{3}\geq 0$, and

[TABLE]

Then $D(a,b)=0$ gives

[TABLE]

with

[TABLE]

The roots $\eta_{\,1}(\lambda_{3})$ and $\eta_{\pm}(\lambda_{3})$ of Eq. (138) can then be obtained using Eqs. (129) with the appropriate substitution of variables. Only those solutions such that $\omega_{12}^{2}\geq 0$ (i.e., is real) are of interest. The results, outlined in detail below, are that (i) there are no degenerate roots if $\omega_{3}^{2}>R_{\delta}^{2}/3$; and (ii) for each $\omega_{3}$ satisfying $0\leq\omega_{3}^{2}\leq R_{\delta}^{2}/3$, there are two values of $\omega_{12}^{2}$ that give degenerate roots.

• Note for use in what follows that
• $\cdot$

$a_{\eta}<0$ for all $\lambda_{3}\geq 0$

• $\therefore$ no Eq. (129) solutions for $\eta$
• $\cdot$

$\sqrt{\alpha_{\eta}}=\sqrt{|a_{\eta}/3|}=\frac{9}{4}\sqrt{8\lambda_{3}+1}$

• $\cdot$

$b_{\eta}=0$ for $\lambda_{3}\,=\,\frac{3}{4}(\sqrt{3}-\frac{5}{3})\equiv\lambda_{b}\approx 0.05$

• $\cdot$

$D(a_{\eta},b_{\eta})=\dfrac{3^{12}}{2^{\,6}}\,\lambda_{3}\,(\lambda_{3}-1)^{3}$

• $\cdot$

$\gamma_{\eta}(\lambda_{3})=\dfrac{|8\lambda_{3}^{2}+20\lambda_{3}-1|}{(8\lambda_{3}+1)^{3/2}}$   (cf. Eq. (128))

• $\gamma_{\eta}(0)=1,\quad\gamma_{\eta}(\lambda_{b})=0,\quad\gamma_{\eta}(1)=1$

If $\lambda_{3}>1$, then

• $\cdot$

$D(a_{\eta},b_{\eta})>0$, equivalent to $\gamma_{\eta}>1$

• $\cdot$

there is one real solution $\eta_{\,1}$ from Eq. (129)

• $\cdot$

Define $\varphi_{\eta}=\tfrac{1}{3}\cosh^{-1}\gamma_{\eta}$

• $\cdot$

$b_{\eta}>0$

• $\cdot$

$\eta_{\,1}=-2\sqrt{\alpha_{\eta}}\cosh\varphi_{\eta}$

• $\cosh\varphi_{\eta}\geq 1$ for all $\varphi_{\eta}$,
• $2\sqrt{\alpha_{\eta}}>\frac{9}{2}\,(3)$
• $\therefore\eta_{\,1}<-\frac{27}{2}$
• $\implies\omega_{12}^{2}\sim(\eta_{\,1}+\frac{9}{4}-\lambda_{3})<-\frac{45}{4}-\lambda_{3}<0$
• $\bullet$

No real $\omega_{12}$ such that Eq. (24) has degenerate roots for $\omega_{3}^{2}=\lambda_{3}\;R_{\delta}^{2}/3>R_{\delta}^{2}/3$

If $\lambda_{3}\leq 1$, then

• $\cdot$

$\omega_{12}^{2}\sim(\eta+\frac{9}{4}-\lambda_{3})\geq 0$ for $\eta\geq 0$

• $\cdot$

$D(a_{\eta},b_{\eta})\leq 0$, equivalent to $\gamma_{\eta}\leq 1$

• $\cdot$

there are three real solutions $\eta_{\,1},\eta_{\pm}$ from Eq. (129)

• $\cdot$

Define $\vartheta=\frac{1}{3}\,\sin^{-1}(\gamma_{\,\eta})$

• $0\leq\sin\vartheta\leq\frac{1}{2}$

• (a)

If $\lambda_{\,b}\leq\lambda_{3}\leq 1$, then

$0\leq\gamma_{\eta}\leq 1$,

$0\leq\vartheta\leq\pi/6$,

$b_{\eta}\geq 0$

• $\cdot$

$\eta_{\,1}=2\sqrt{\alpha_{\eta}}\sin\vartheta$

• $\ \ \therefore\phantom{\eta_{1}}\eta_{\,1}\geq 0$
• $\implies\omega_{12}^{2}\,>0$
• $\cdot$

$\eta_{+}=-\sqrt{\alpha_{\eta}}\,\sin\vartheta\,+\,\sqrt{3}\,(\alpha_{\eta}-\alpha_{\eta}\sin^{2}\vartheta)^{1/2}$

• $\phantom{\eta_{+}}=2\sqrt{\alpha_{\eta}}\,\sin(\pi/3-\vartheta)$
• $\ \ \therefore\phantom{\eta_{1}}\eta_{+}\geq 0$
• $\implies\omega_{12}^{2}\,>0$
• (b)

If $0\leq\lambda_{3}\leq\lambda_{\,b}$, then

$1\geq\gamma_{\eta}\geq 0$,

$\pi/6\geq\vartheta\geq 0$,

$b_{\eta}\leq 0$

• $\cdot$

$\eta_{\,1}=-2\sqrt{\alpha_{\eta}}\sin\vartheta$

• $\ \ \therefore\ -\frac{9}{4}\,\leq\,\eta_{\,1}\leq\,0$
• $\implies\omega_{12}^{2}\sim\eta_{\,1}+\frac{9}{4}-\lambda_{3}\geq 0$,
• since $\eta_{\,1}\in[-\frac{9}{4},0\,]$ as $\lambda_{3}\in[0,\lambda_{b}]$
• $\cdot$

$\eta_{+}=\sqrt{\alpha_{\eta}}\,\sin\vartheta\,+\,\sqrt{3}\,(\alpha_{\eta}-\alpha_{\eta}\sin^{2}\vartheta)^{1/2}$

• $\phantom{\eta_{+}}=2\sqrt{\alpha_{\eta}}\,\sin(\pi/3+\vartheta)$
• $\therefore\eta_{+}\geq 0$
• $\implies\omega_{12}^{2}\,>0$
• $\bullet$

2 real $\omega_{12}^{2}$ such that Eq. (24) has degenerate roots for $0\,\leq\,\omega_{3}^{2}\,\leq\,R_{\delta}^{2}/3$

The solutions for $\omega_{12}^{2}$ become equal at $\omega_{3}^{2}=R_{\delta}^{2}/3$, as shown in Fig. 1, corresponding to the case $a=0=b$. There is then a three-fold degenerate root $z=0$ of Eq. (24). Recall that a solution to $D(a,b)=0$ for real $a,b$ requires $a=\omega_{12}^{2}+\omega_{3}^{2}-3R_{\delta}^{2}\leq 0$, which is readily verified for the solutions obtained above. Scaling $a$ according to Eqs. (136) and (137), dividing by $R_{\delta}^{2}/3$, and using the maximum value $\eta_{\mathrm{max}}=\sqrt{\alpha_{\eta}}=27/4$ at $\lambda_{3}=1$ gives

[TABLE]

Appendix D Vector Model

There is a simple physical interpretation for the action of the propagator $e^{-\Gamma\,t}$ when, as is most common, the matrix $\Gamma$ has three distinct eigenvalues. Supplementary details of the model introduced in section V.2 are presented here. Consider the case of one real eigenvalue and two complex conjugate eigenvalues. Results for the other possibility, that of three real eigenvalues, are obtained directly from Eq. (145) in what follows.

The eigenvalues of $-\Gamma$ are the roots $s_{1}=z_{1}-\bar{R}$ and $s_{2,3}\equiv s_{\pm}=-z_{1}/2\pm i\,\varpi-\bar{R}$, obtained from Eq. (29), with real $z_{1}$ given in Eqs. (129). The associated eigenvectors are $\bm{s}_{1}$ and the complex conjugate pair $\bm{s}_{\pm}$. The relation between $\bm{s}_{\pm}$ and the real vectors $\bm{\tilde{s}}_{2}$ and $\bm{\tilde{s}}_{3}$ defined in Eq. (LABEL:Real_s_2,3) is

[TABLE]

Defining $\bm{\tilde{s}}_{1}\equiv\bm{s}_{1}$ gives a set $\bm{\tilde{s}}_{i}$ of three linearly independent vectors that can be used as an alternative basis for representing arbitrary system states. We then have

[TABLE]

Similarly,

[TABLE]

These relations, together with $e^{-\Gamma\,t}\bm{\tilde{s}}_{1}=e^{s_{1}}\bm{\tilde{s}}_{1}$, yield the propagator $e^{-\tilde{\Gamma}\,t}$ for the evolution of states $\tilde{\mathcal{M}}=\sum_{i}\tilde{\mathcal{M}}_{i}\bm{\tilde{s}}_{i}$ expressed in the $\{\bm{\tilde{s}}_{i}\}$ basis, as given in Eq. (97).

As noted in Eq. (88), matrix $P$ generated from the $\{\bm{\tilde{s}}_{i}\}$ entered as column vectors transforms from the $\{\bm{\tilde{s}}_{i}\}$ basis to the standard basis, with $P^{-1}=\mathrm{adj}\,P/\det P$ giving the desired $\tilde{\mathcal{M}}$ starting with $\mathcal{M}$ in the standard basis. One easily shows that $\det P=\bm{\tilde{s}}_{1}\cdot(\bm{\tilde{s}}_{2}\times\bm{\tilde{s}}_{3})$, and row $i$, column $l$ of $\mathrm{adj}\,P$ is $(\bm{\tilde{s}}_{j}\times\bm{\tilde{s}}_{k})_{l}$ for cyclic permutation of $i=1,j=2,\mathrm{and}\ k=3$ to obtain

[TABLE]

The eigenvectors needed to construct the real basis are obtained in the usual fashion as solutions to (s\mbox{\openone}+\Gamma)\bm{s}=0 for each eigenvalue $s_{i}$. The solution for the three components of each eigenvector is overdetermined, by construction, so any one of the three equations is a linear combination of the other two and is redundant. We are free to assign any (nonzero) value to one of the components, leaving two equations and two unknowns. There are three different but equivalent forms for the eigenvector solution depending on which two equations are chosen. Setting the third component equal to one for simplicity gives an expression for the other two components involving a common denominator. Scaling the result by this factor gives the following result for eigenvector $\bm{s}_{i}$, with the left arrow signifying that the columns of the matrix map to $\bm{s}_{i}$:

[TABLE]

The different columns give equivalent results, as discussed in section V. In the absence of relaxation, the real root of Eq. (9) is $s_{1}=0$ with eigenvector $\bm{s}_{1}=(\omega_{1},\omega_{2},\omega_{3})$, which is the rotation axis for the resulting time evolution. In the case $\bm{\omega}_{e}=0$, in which $\Gamma$ is already diagonal, the coordinates reduce to the standard coordinate system as required.

One might recognize the righthand side of Eq. (145) as $\mathrm{adj\,}A(s_{i})$ from Eqs. (7–LABEL:adjA_coeff), with $\mathrm{adj}\,A(s_{i})=\mathrm{adj}\,A_{\mathrm{p}}(z_{i})$, since $s_{i}=z_{i}-\bar{R}$ and $R_{i}-\bar{R}=R_{i\mathrm{p}}$. We thus have the real basis vectors $\bm{\tilde{s}}_{2,3}\equiv\bm{\tilde{z}}_{2,3}$ equal to the respective real, imaginary parts of $\bm{z}_{+}=\mathrm{adj}\,A_{\mathrm{p}}(z_{+})$ according to Eq. (LABEL:Real_s_2,3), with $z_{+}=-z_{1}/2+i\,\varpi$. Then, using Eq. (LABEL:adjAp_Poly) for $\mathrm{adj}\,A_{\mathrm{p}}(z_{i})$ in polynomial form and eliminating common scale factors, the real basis vectors defining the oblique coordinate system can be written concisely as

[TABLE]

The result for $\bm{\tilde{z}}_{1}$ can be obtained directly from Eq. (145) with the substitutions $s_{i}\rightarrow z_{i}$ and $R_{i}\rightarrow R_{ip}$ for the corresponding parameters associated with $\Gamma_{\mathrm{p}}$. One can readily deduce the coefficient matrices $A_{0\mathrm{p}}$ and $A_{1\mathrm{p}}$ from Eq. (145) and the expression for $\bm{\tilde{s}}_{1}$ in Eq. (146) without recourse to the definitions for each element given in Eq. (LABEL:adjAp_Poly). The matrices are also given as simple functions of $\Gamma_{\mathrm{p}}$ in Eq. (39). For convenient reference, each coefficient matrix is written below.

[TABLE]

with $-R_{1\mathrm{p}}=R_{2\mathrm{p}}+R_{3\mathrm{p}}$ and cyclic permutations, since $\sum_{i}R_{i\mathrm{p}}=0$ by construction in the original matrix partitioning.

D.1 Measures of obliquity

Bloch equation dynamics are simple in the oblique coordinates of the model, consisting of independent rotation relaxation elements. This section provides examples that quantify the degree to which the plane of rotation is oblique to the axis $\bm{\tilde{z}}_{1}$ representing simple exponential decay. In what follows, the first column of $\mathrm{adj}\,A_{\mathrm{p}}$ is arbitrarily chosen to calculate the coordinate basis $\{\bm{\tilde{z}}_{i}\}$ in the case $R_{1}=R_{2}$. Similar results are obtained using any of the other columns.

D.1.1 Off resonance, $\bm{\omega}_{e}=(0,\omega_{2},\omega_{3})$

Off resonance, in contrast to the on-resonance example of section V.3.2, $\bm{\tilde{z}}_{1}$ is neither aligned with $\bm{\omega}_{e}$, nor is it orthogonal to the $(\bm{\tilde{z}}_{2},\bm{\tilde{z}}_{3})$-plane. Calculating the $\bm{\tilde{z}}_{i}$ as above provides the normal to the plane, $\bm{\tilde{n}}_{23}=\bm{\tilde{z}}_{2}\times\bm{\tilde{z}}_{3}$. We then have

[TABLE]

and

[TABLE]

which bears little resemblance to $\bm{\tilde{z}}_{1}$. Yet, scaling $\bm{\tilde{z}}_{1}$ by $f_{s}=-(\bm{\tilde{n}}_{23})_{1}/(\bm{\tilde{z}}_{1})_{1}$ from the first components gives, for components two and three, $f_{s}\bm{\tilde{z}}_{1}-\bm{\tilde{n}}_{23}\propto q(z_{1})$, the characteristic polynomial for $-\Gamma_{\mathrm{p}}$, which is zero when evaluated at its root $z_{1}$. Thus, within a scale factor or, equivalently, when both both vectors are normalized, we can write simply

[TABLE]

D.1.2 Off resonance, $\bm{\omega}_{e}=(\omega_{1},0,\omega_{3})$

Similarly, for $\omega_{2}=0$,

[TABLE]

and

[TABLE]

Scaling $\bm{\tilde{z}}_{1}$ by $f_{s}=-(\bm{\tilde{n}}_{23})_{2}/(\bm{\tilde{z}}_{1})_{2}$ gives $f_{s}\bm{\tilde{z}}_{1}-\bm{\tilde{n}}_{23}\propto q(z_{1})$ for components one and three, so that

[TABLE]

within a scale factor.

D.1.3 $\bm{\omega_{1}=\omega_{2}=\omega_{3}\equiv\omega}$

In this case,

[TABLE]

and

[TABLE]

Scaling $\bm{\tilde{z}}_{1}$ by $f_{s}=(\bm{\tilde{n}}_{23})_{1}/(\bm{\tilde{z}}_{1})_{2}$ gives both $f_{s}(\bm{\tilde{z}}_{1})_{1}-(\bm{\tilde{n}}_{23})_{2}$ and $f_{s}(\bm{\tilde{z}}_{1})_{3}-(\bm{\tilde{n}}_{12})_{3}$ proportional to $q(z_{1})$, so that the vectors can be scaled to satisfy

[TABLE]

Appendix E An Alternative Method for Calculating

an Eigenvector

Equation (145) is simply \mathrm{adj}\,(s_{i}\mbox{\openone}+\Gamma) from Eqs. (6–LABEL:adjA_coeff). One therefore happens upon the modest result, apparently unrecognized, that an eigenvector $\bm{\upsilon}$ corresponding to a distinct eigenvalue $\upsilon$ of operator $\Upsilon$ can be obtained as

[TABLE]

seen as follows. Recall, the characteristic polynomial p(s)=\det(s\mbox{\openone}-\Upsilon) equals zero for eigenvalue $s=\upsilon$, and (s\mbox{\openone}-\Upsilon)^{-1}=\mathrm{adj}\,(s\mbox{\openone}-\Upsilon)/p(s) from Eq. (11). Then

[TABLE]

Only a single column of the adjugate matrix is required, so the method is fairly efficient. However, the trivial zero eigenvector solution can be one of the columns, requiring further completion of the adjugate to obtain the desired eigenvector.

For the case of degenerate eigenvalues, the method is incomplete. When the nullity (dimension of the null space) of (\upsilon\mbox{\openone}-\Upsilon) equals the order of the degeneracy, $k$ (i.e, the rank equals the dimension of the operator, $n$, minus $k$), there are $k$ distinct eigenvectors, but the method fails, returning only the zero eigenvector. If there is not a complete set of eigenvectors (the degenerate eigenvalue is defective in that the nullity is less than $k$), and the rank is greater than $n-k$), the method appears to return the eigenvectors that exist, but one rarely needs these, since the matrix $\Upsilon$ is not diagonalizable in this case.

Appendix F Limiting Cases

The solutions are evaluated and confirmed for $R_{1}=R_{2}$ and a representative set of limiting cases that can be readily solved by other methods.

F.1 Three distinct roots

Three examples are presented representing the separate cases $a=0$ and $b=0$.

• (i)

$b=0$, $a\neq 0$

According to the defining relations for $a$ and $b$ in Eq. (76), the condition $b=0$ implies $\omega_{12}^{2}=2R_{\delta}^{2}(1+\frac{1}{3}\lambda_{3})$, using Eq. (3) for $\omega_{e}^{2}$ and Eq. (78) for $\omega_{3}$. Then

[TABLE]

The roots of Eq. (24) are easily obtained, giving

[TABLE]

There are two cases, depending on the sign of $a$.

• • (1) $a>0$

Then Eq. (62) gives

[TABLE]

There is no exponential decay contribution due to this term, with the overall factor $e^{-\bar{R}t}$ in the final expression for $e^{-\Gamma t}$ providing a single system decay rate $\bar{R}$. Example (1)

• (***)

Choose $R_{\delta}=0$ to obtain

• (***)

$b=0$,   $a=\omega_{e}^{2}$,   $\varpi=\omega_{e}$

In this case, Eq. (168) represents a rotation about the field $\bm{\omega}_{e}$.

The propagator $U_{R}$ for a rotation about $\bm{\omega}_{e}$ is readily obtained by transforming to a coordinate system with new $z$-axis aligned with $\bm{\omega}_{e}$, rotating by angle $-\omega_{e}t$ about this axis, then transforming back to the original coordinates. Specifying the orientation of $\bm{\omega}_{e}$ in terms of polar angle $\theta$ and azimuthal angle $\phi$ relative to the $z$- and $x$-axes, respectively, one has $U_{R}=U_{z}(-\phi)U_{y}(-\theta)U_{z}(-\omega_{e}\,t)U_{y}(\theta)U_{z}(\phi)$ in terms of the elementary operators $U_{y}$ and $U_{z}$ for rotations about the $y$- and $z$- axes, respectively. Then $U_{R}$ provides a verification of the Eq. (168) result upon substituting $\cos\phi=\omega_{1}/\omega_{12}$, $\sin\phi=\omega_{2}/\omega_{12}$, $\cos\theta=\omega_{3}/\omega_{e}$, $\sin\theta=\omega_{12}/\omega_{e}$.

• • (2) $a<0$

for $\lambda_{3}<1$ gives $\varpi\rightarrow i\,\mu=i\,\sqrt{|a|}$ and

[TABLE]

Example (2)

• (***)

Choose $\omega_{1}^{2}=2R_{\delta}^{2}$, $\omega_{2}=0$, $\lambda_{3}=0$ to obtain

• (***)

$b=0$,   $a=-R_{\delta}^{2}$,   $\mu=R_{\delta}$

Equation (169) then gives

[TABLE]

For an independent calculation, the matrix $-\Gamma_{\mathrm{p}}$ can be diagonalized, with eigenvalues given by the $z_{i}$ and associated real-valued eigenvectors. The simple exponential of the diagonalized matrix is then transformed back to the original basis in the standard fashion using the matrix of eigenvectors and its inverse to obtain $e^{-\Gamma_{\mathrm{p}}\,t}$ as given above.

• (ii)

$a=0$. $b\neq 0$

The condition $a=0$ implies $\omega_{e}^{2}=3R_{\delta}^{2}$, leading to

[TABLE]

and root $z_{1}=-\mathrm{sgn}(b)|b|^{1/3}$ from Eq. (129). For $\mathrm{sgn}(b)=\pm 1$ and the definition $\tilde{\lambda}_{3}=|1-\lambda_{3}|^{1/3}$, we have accordingly

[TABLE]

Although the form of Eq. (62) does not simplify in this case as appreciably as for $b=0$, both the root $z_{1}$, which determines the decay rate, and the oscillatory frequency $\varpi$ are simple multiples of $R_{\delta}$. Example (3)

• (***)

Choose $\omega_{e}^{2}\rightarrow\omega_{1}^{2}=3R_{\delta}^{2}$,   $\omega_{2}=0=\omega_{3}$

Most off-diagonal elements of $\Gamma_{\mathrm{p}}$ are equal to zero in this case, and $\tilde{\lambda}_{3}=1$ for the Eq. (175) input parameters to Eq. (62). Defining $\kappa=(\sqrt{3}/2)R_{\delta}$ and combining the sums of trigonometric functions that appear on the diagonal gives the succinct form

[TABLE]

Again, the matrix $-\Gamma_{\mathrm{p}}$ is diagonalizable, providing a simple result for the matrix exponential in the eigenbasis and a straightforward means for calculating $e^{-\Gamma_{p}t}$ as obtained above. The associated eigenvectors are complex-valued in this case, making the algebra slightly more tedious. Alternatively, one can readily verify that $d/dt\,e^{-\Gamma_{\mathrm{p}}\,t}=-\Gamma_{\mathrm{p}}\,e^{-\Gamma_{\mathrm{p}}\,t}$.

F.2 Two equal roots

Degenerate roots require $\gamma=1$. For a given $\omega_{3}^{2}=\lambda_{3}R_{\delta}^{2}/3$, with $0\leq\lambda_{3}\leq 1$, there are two values $\omega_{12}^{2}$ that satisfy $\gamma=1$. Consider $\lambda_{3}=0$, in which case Eqs. (78) and (79) give

[TABLE]

• (i)

$\omega_{12}=0$

This is the case of pure relaxation, with $\Gamma_{\mathrm{p}}$ reduced to the diagonal elements $\{R_{\delta},R_{\delta},-2R_{\delta}\}$. We have $a=-3R_{\delta}^{2}$, $b=-2R_{\delta}^{3}<0$, and

[TABLE]

from Eq. (129). Thus, Eq. (72) gives the expected result

[TABLE]

• (ii)

$\omega_{12}^{2}=\frac{9}{4}R_{\delta}^{2}\rightarrow\omega_{1}^{2}$

Then $a=-3R_{\delta}^{2}/4<0$, $b=R_{\delta}^{3}/4>0$, and

[TABLE]

resulting in

[TABLE]

Verifying that $d/dt\,e^{-\Gamma_{\mathrm{p}}\,t}=-\Gamma_{\mathrm{p}}\,e^{-\Gamma_{\mathrm{p}}\,t}$ is fairly straightforward and represents the simplest test of the solution, since $\Gamma_{\mathrm{p}}$ is not diagonalizable.

F.3 Three equal roots

There is a three-fold degenerate root $z_{i}=0$ in the case $a=0=b$, since $q(z)\rightarrow z^{3}$. This requires $\omega_{e}^{2}=3R_{\delta}^{2}$ from Eq. (76), which then forces $\omega_{3}^{2}=R_{\delta}^{2}/3$ in the expression for $b$. As noted previously, the Cayley-Hamilton theorem is simple to apply directly in this case, since $q(\Gamma_{\mathrm{p}})=\Gamma_{\mathrm{p}}^{3}=0$. The series expansion of $e^{-\Gamma_{\mathrm{p}}\,t}$ is therefore truncated, giving the Eq. (73) result.

F.4 **On resonance **

When $\omega_{3}=0$, $b$ can be written in the form $R_{\delta}(a+R_{\delta}^{2})$ from Eq. (76), with $a\rightarrow\omega_{12}^{2}-3R_{\delta}^{2}$. The characteristic polynomial then becomes $z^{3}+R_{\delta}^{3}+a(z+R_{\delta})$, so that, by inspection,

[TABLE]

The solution for $e^{-\Gamma_{p}\,t}$ using Eq. (62) with the above parameters yields the solution for $e^{-\Gamma\,t}$ obtained originally in Torrey (1949) for the case $\varpi\neq 0$. As discussed above, if $\omega_{12}=3R_{\delta}/2$, there is a two-fold degeneracy in the roots, giving the solution in Eq. (188) for $e^{-\Gamma_{p}\,t}$.

For $\omega_{12}<3R_{\delta}/2$, the sinusoidal terms become the corresonding hyperbolic functions, as noted earlier, with $\cos\varpi\,t\rightarrow\cosh\mu\,t$ and $\sin\varpi\,t/\varpi\rightarrow\sinh\mu\,t/\mu$, where now $\mu=\sqrt{(\tfrac{3}{2}R_{\delta})^{2}-\omega_{12}^{2}}$.