Approximate representations of shaped pulses using the homotopy analysis method
The evolution of nuclear spin magnetization during a radiofrequency pulse in the absence of relaxation or coupling interactions can be described by three Euler angles. The Euler angles, in turn, can be obtained from the solution of a Riccati differential equation; however, analytic solutions exist only for rectangular and hyperbolic-secant pulses. The homotopy analysis method is used to obtain new approximate solutions to the Riccati equation for shaped radiofrequency pulses in nuclear magnetic resonance (NMR) spectroscopy. The results of even relatively low orders of approximation are highly accurate and can be calculated very efficiently. The results are extended in a second application of the homotopy analysis method to represent relaxation as a perturbation of the magnetization trajectory calculated in the absence of relaxation. The homotopy analysis method is powerful and flexible and is likely to have other applications in magnetic resonance.
Numerous aspects of nuclear magnetic resonance (NMR) spectroscopy are formulated in terms of differential equations, few of which have closed-form analytical solutions. In an era characterized by ever-increasing computational capabilities, numerical solutions to such differential equations are always possible and are frequently the preferred approach for applications such as data analysis. However, approximate solutions can provide useful formulas and insights that are difficult to discern from purely numerical results.
As one example, the net evolution of the magnetization of an isolated spin during a radiofrequency pulse, i.e. in the absence of relaxation and scalar or other coupling interactions, can be described by three rotations with Euler angles α(τp), β(τp), and γ(τp), in which τp is the pulse length (Zhou et al., 1994; Siminovitch, 1997a, b). Shaped pulses, in which the amplitude (Rabi frequency), phase, or radiofrequency are time dependent, are widely applied in modern NMR spectroscopy and other magnetic resonance techniques (Geen and Freeman, 1991; Emsley and Bodenhausen, 1992; Kupc̆e et al., 1995; Cavanagh et al., 2007). The Euler angles for an arbitrary shaped pulse can be extracted from a numerical calculation in which the shaped pulse is represented by a series of K short rectangular pulses with appropriate amplitudes and phases. Thus, the propagator for a shaped pulse expressed in the Cartesian basis is given by the following (Siminovitch, 1995):
in which χ±(τp) is given as follows:
Ik are the Cartesian spin operators, and the product is time-ordered from right to left. The propagator for the kth rectangular pulse segment is given as follows:
in which κ(τp) and η(τp) are given by the following:
In Eq. (4), ω1k, ϕk, and Δtk are the radiofrequency field strength, phase angle, and duration of the kth pulse segment; Ωk is the resonance offset in the rotating frame of reference during the kth pulse segment (and is constant if the offset is fixed); is the effective field; and is the tilt angle. Values of α(τp), β(τp), and γ(τp) are then obtained from the matrix elements of U.
Alternatively, the Euler angles can be determined from the solution of the following Riccati equation (Zhou et al., 1994):
in which the following applies:
where and ωx(t) and ωy(t) are the Cartesian amplitude components of the radiofrequency field in the rotating frame of reference. After solution of the Riccati equation, β(τp) and γ(τp) are obtained from the magnitude and argument of f(τp), and the value of α(τp) is obtained by integration as follows:
The Riccati equation can be transformed into a second-order differential equation as follows:
by using the following definition:
A more compact form is obtained by defining the following:
to yield the following:
The Riccati differential equation only can be solved analytically for a single rectangular or hyperbolic-secant pulse (Silver et al., 1985; Zhou et al., 1994; Rourke, 2002). Approximate solutions for arbitrary shaped pulses have been derived by the perturbation theory for the limits of small, using Eq. (11), and large, using Eq. (5), resonance offsets (Li et al., 2014); however, perturbation theory is unwieldy to apply to a high order and, obviously, depends on the perturbation parameters being small in some respects.
The homotopy analysis method (HAM) is a fairly recent development, first reported in 1992 (Liao, 1992), for approximating solutions to differential equations, particularly nonlinear ones. HAM does not depend on small parameters, unlike perturbation theory, and has proven powerful in a number of applications outside of NMR spectroscopy (Liao, 2012). The present paper illustrates HAM by application to the solutions of Eqs. (5) and (11) and, subsequently, by extension to the Bloch equations, including relaxation.
In topology, a pair of functions defining different topological spaces are said to be homotopic if the shape defined by one function can be continuously transformed (deformed in the lexicon of topology) into the shape defined by the other. Analogously, the essence of HAM is to map a function of interest, here y(t), to a second function, Φ(t;q), which has a known solution and is a function of both t and an embedding parameter .
This relationship is established by constructing the homotopy as follows (Liao, 2012):
in which ℒ[ ] is a linear (differential) operator, and 𝒩[ ] is an (nonlinear differential) operator satisfying the following:
where y0(t) is an initial approximation for the desired solution y(t), c0≠0 is a convergence control parameter, and H(t)≠0 is an auxiliary function (vide infra). When q=0, the homotopy becomes as follows:
Therefore, when , Eq. (13) requires . Similarly, when q=1, the homotopy becomes as follows:
Therefore, when , Eq. (14) requires . Stated more succinctly, as q increases from 0→1, Φ(t;q) deforms from the initial approximation y0(t) to the exact solution y(t). To proceed, the Maclaurin series for Φ(t;q) is assumed to exist as follows:
in which yn(t) is given by the following:
Conditions concerning convergence of the series are discussed by Liao (2012). Equation (17) has the desired property and yields the following:
HAM then consists of successively determining yn(t), beginning with the initial approximation y0(t), until y(t) is approximated to the desired accuracy. The choices of ℒ[ ], y0(t), c0, and H(t) provide considerable flexibility in finding approximate solutions to differential equations. For simplicity in the following, the auxiliary function is H(t)=1.
The iterative algorithm in HAM is illustrated by application to the second-order differential form of the Riccati equation. In the first example, the nonlinear operator is obtained from Eq. (11) as follows:
in which g(t) is an arbitrary function. The linear operator is chosen to be the following:
and the initial approximation is y0(t)=1.
From the relationships of Eqs. (13) and (14) embedded in the initial homotopy, Eq. (12), the zeroth-order deformation equation is defined as follows (Liao, 2012):
The derivative of Eq. (22), with respect to q, yields the first-order deformation equation as follows:
The limit q→0 gives the following:
in which the second line is obtained using and Eq. (18). Substituting for 𝒩[ ], ℒ[ ], and y0(t) yields the following:
in which the final line is obtained using . This differential equation does not contain a term proportional to y1(t). Hence, the homogenous equation (setting the right-hand side to 0) can be solved by two successive integrations, and the inhomogeneous solution is obtained by the technique of variation of parameters (Arfken et al., 2013). The solution is as follows:
The higher-order approximations yn(t) are obtained in a similar fashion. The nth derivative, with respect to q of Eq. (22), yields the following (for n>1):
Executing the derivatives, taking the limit q→0, and dividing both sides of the equation by n! gives the following:
with the solution obtained by the same approach as for Eq. (26), as follows:
Successive use of Eqs. (26) and (29) allows y(t) and, hence, f(t) to be determined to arbitrary accuracy, as follows:
in which N is the order of approximation. For completeness, the derivatives of Eqs. (26) and (29) are, respectively, as follows:
Results obtained using y0(t)=1, together with Eqs. (26), (29), and (30), will be called method 1 in the following discussion. The iterated form of the above expressions for yn(t) have similarities to the Fourier integrals obtained from average Hamiltonian theory by Warren (1984).
The above choices of ℒ[ ] and y0(t) are not unique. Different choices lead to different series approximations and, hence, to different qualitative and quantitative results. As a second example, Ω(t)=Ω is assumed to be fixed and only amplitude-modulated pulses ω(t) with x phase are considered (these assumptions can be relaxed as needed). Returning to Eq. (8), the homotopy is defined as follows:
This choice of y0(t) satisfies the following:
and is the exact on-resonance solution for y(t). Consequently, the first-order deformation equation leads to the following:
The solutions to the homogeneous equation (setting the right-hand side to 0) are . The method of variation of parameters then gives the inhomogeneous solution as follows:
The nth-order deformation equation for n>1 is as follows:
with the following solution:
Each yn(t) is proportional to Ωn, and these results yield the following power series in Ω for y(t):
which is substituted into Eq. (30) to obtain f(t). Results using Eqs. (39), (41), and (42) will be called method 2 in the following discussion. For completeness, the derivatives of Eqs. (39) and (41) are as follows:
Numerical integration was performed using the trapezoid method implemented in Python 3.6. Pulse shapes were discretized in 1000 increments. Rectangular pulses were simulated using Hz and an on-resonance 90∘ pulse length of 10.0 µs or Hz and an on-resonance 90∘ pulse length of 1 ms. EBURP-2 (Geen and Freeman, 1991) and Q5 (Emsley and Bodenhausen, 1992) pulses were simulated using a maximum Hz and 90∘ pulse lengths of 455.2 and 504.9 µs, respectively. REBURP (Geen and Freeman, 1991) pulses were simulated using a maximum Hz and a 180∘ pulse length of 626.5 µs. WURST-20 (Kupc̆e and Freeman, 1995) pulses were simulated using a maximum Hz, a frequency sweep of 50 000 Hz, and a pulse length of 440.0 µs.
Equation (7) can be recast as follows:
for numerical calculations. Alternatively, α(τp) can be obtained from the argument of f(τp) calculated for the time-reversed pulse (Li et al., 2014). The latter is more computationally demanding, but more numerically stable, and was used for the results presented herein.
2.2 Results and discussion
In the present applications, HAM converts the second-order Riccati differential equation, Eq. (8), which cannot be solved directly, into a series of second-order differential equations that have convenient solutions. The choice of y0(t)=1 in method 1 leads to obtaining simple iterative solutions that can be calculated very efficiently. The form of y0(t) given in Eq. (36) could also be used in Eq. (24) to obtain an alternative expression for y1(t) to be substituted into Eqs. (29) and (30). The resulting first-order expressions for y(t) are usually more accurate than the first-order results obtained using y0(t)=1, but this advantage becomes less pronounced at higher orders of approximation and comes at increased computational cost. Thus, Eqs. (26), (29), and (30) are most suitable in practice.
A first example of the results of the above analysis is given for a rectangular 90∘ pulse in Fig. 1. The integrals in Eqs. (26) and (29) can be performed analytically for a rectangular pulse with amplitude ω1. For example, using Eq. (26) yields the following:
however, analytic calculations of higher order yn(t) do not have advantages over numerical integration. As shown in Fig. 1a, b, the second- and third-order results obtained with method 1 and are nearly indistinguishable from the exact result of Eq. (3) (using τp=Δτk) over the range of resonance offsets from 0 to . The first-order result provides a highly accurate estimate of γ(τp) but overestimates β(τp). The role of the convergence control parameter c0 is illustrated in Fig. 1c, d. A value of was chosen, using Eqs. (46) and (30), to scale the first-order result for β(τp) to be equal to at Ω=0. As shown, the resulting first-order result, using method 1, is now nearly exact at all resonance offsets. In the present application, adjusting the convergence control parameter provides accuracy equivalent to 1 or 2 additional higher orders of approximation. Remarkably, this same value of c0 works well for a rectangular 180∘ pulse (not shown) and 90∘ EBURP-2, 90∘ Q5, and 180∘ REBURP and WURST inversion pulses (vide infra).
In contrast to the results of method 1, the power series for y(t) obtained, using method 2 with , even to the third order in Ω, is accurate for β(τp) only to slightly more than . When , the third-order power series has improved accuracy for resonance offsets up to nearly . However, further increases in accuracy at larger resonance offsets require very large orders of approximation N in Eq. (42). For example, extending the accuracy of the power series for β(τp) to offsets requires N=50. The differences between the results of method 1 and method 2 reflect the inevitable shortcomings of power series and perturbation approaches when the expansion parameter is not small.
A more challenging example is given by the 90∘ EBURP-2 pulse (Geen and Freeman, 1991). In principle, the integrals in Eqs. (26) and (29) can be performed analytically because the pulse shape is expressed as a Fourier series (as are other pulses in the BURP (Geen and Freeman, 1991) and SNOB (Kupc̆e et al., 1995) families). In practice, the number of terms that must be calculated becomes very large, and numerical integration is much more efficient. Calculations using method 1 are shown in Fig. 2. With , the fifth-order approximation is extremely accurate compared with numerical calculations using Eqs. (1)–(3) (Fig. 2a–c). With (Fig. 2d–f), even the small deviations observed for the fifth-order HAM approximation are eliminated, and the third-order result is accurate, except at the edge of the excitation band. In contrast, perturbation theory or power-series expansions (method 2) are extremely poor at reproducing β(τp), essentially failing as soon as Ω is non-zero (not shown). The accuracy of the method 1 approximations over the full range of resonance offsets shows that HAM, with an appropriate choice of linear operator and starting functions, can provide approximate solutions valid far beyond the range of perturbation theory.
The Gaussian Q5 90∘ pulse (Emsley and Bodenhausen, 1992) has a more complicated amplitude modulation profile than the EBURP-2 pulse and requires higher orders of approximation to obtain accurate results. Results obtained for method 1 with fifth- and seventh-order approximations are shown in Fig. 3. The seventh-order results are highly accurate for both and . The choice of has a remarkable effect in terms of increasing the accuracy the fifth-order approximation to nearly that of the seventh-order result.
The application of HAM is not limited to 90∘ pulses or to amplitude-modulated pulses. Figure 4 shows the performance of method 1 for the 180∘ REBURP (Geen and Freeman, 1991) and WURST-20 inversion (Kupc̆e and Freeman, 1995) pulses. As for the EBURP-2 pulse, the fifth-order approximation for the REBURP pulse is highly accurate for both and . The third-order approximation also is highly accurate when . The WURST-20 pulse uses a linear frequency shift, generated by applying a quadratic phase shift during the pulse, and is an example of a phase-modulated or complex waveform. Again, the more complicated waveform requires higher order approximation, but 11th-order, with , or ninth-order, with , results are highly accurate.
Method 2 yields a power series for y(t). If , the resulting series is identical to the power series expansion obtained from perturbation theory (Li et al., 2014), while provides additional accuracy. However, as noted above, the power series requires very high orders N to obtain accuracy comparable to results from modest orders using method 1. Thus, method 1 is much more powerful for general calculations; however, the power series leads to a convenient expression for the near-resonance phase shift γ(τp). The first-order power series for y(t), assuming , yields the following:
in which the second equality is the expansion to first order in Ω, and the resulting trigonometric functions have been simplified. This result is identical to the previously reported result from the first-order perturbation theory (Li et al., 2014). The argument of the first-order approximation of f(t) is a good estimate of the phase γ(τp) of the transverse magnetization following the pulse. As noted above, the phase α(t) is obtained by repeating the calculation with the time-reversed pulse. Therefore, as concluded from the perturbation theory, an amplitude-modulated shaped pulse acts as an ideal rotation of the angle β(τp) preceded and followed by time delays τα and τγ over the frequency range for which the first-order approximation holds (Lescop et al., 2010; Li et al., 2014), as in the following:
For a 90∘ pulse, the above equations can be written compactly as follows:
The ratios and are the average projections of a unit vector onto the z axis and −y axis, respectively, over the duration of the pulse (for a vector is oriented along the z axis at time 0).
The above explications have focused on solutions to the transformed Riccati equation in Eq. (8). However, HAM also could be applied directly to the original Riccati equation of Eq. (5). For example, by analogy to method 2, choosing the following:
in which f0(t) is the exact solution for Ω=0, yields the following series solution:
The result obtained from the first-order deformation equation is as follows:
However, additional terms in the series lack the simple iterative structure shown in Eqs. (29) and (41) because of the increasing complexity of the higher derivatives of Φ2(t;q) that must be calculated for the nth order deformation equation. For example, the differential equations for the next two terms in the series for f(t) become the following:
In addition, results obtained using Eq. (30) to obtain f(t) from y(t) generally are more accurate than results obtained by direct calculation of f(t) at the same order of approximation. Thus, in this particular application, use of HAM with the transformed Riccati equation, Eq. (8), yields more convenient expressions. Nonetheless, this example demonstrates the particular power of HAM in directly converting the solution of a nonlinear differential equation into a series of linear first-order differential equations, which always can be solved by integration (Liao, 2012).
For many applications, the Euler angles for a shaped pulse are easily obtained from Eqs. (1)–(3). However, calculations for method 1 using Eqs. (26), (29), and (30) are extremely efficient. In Python 3.6, the seventh-order HAM approximation for the Q5 pulse is approximately 20 times faster than direct calculations using Eqs. (1)–(3). Thus, these approximations may be particularly useful for the computational design of radiofrequency pulses in which many iterations of a search or optimization routine are necessary (Gershenzon et al., 2008; Li et al., 2011; Nimbalkar et al., 2013; Asami et al., 2018).
The Euler angle representation is particularly convenient because, once calculated, the Euler angles can be used to determine the outcome of a shaped pulse applied to arbitrary initial magnetization. The Riccati equation can be extended to incorporate radiation damping, but not relaxation, as discussed by Rourke (2002). However, the Euler angles can serve to generate the initial approximations for a second application of HAM to obtain approximate solutions to the Bloch equations for particular initial conditions, including relaxation. In the following, Ω(t)=Ω is assumed to be fixed, and only amplitude-modulated pulses ω(t) with x phase are considered (these assumptions can be relaxed as needed). The Bloch equations for a pulse applied with x phase can be written in the following form:
in which are the Cartesian components of the magnetization, and M0 is the equilibrium magnetization. The use of transformed variables rather than M(t) simplifies the following discussion. The linear operator is chosen as , in which is an arbitrary vector function. The nonlinear operator is as follows:
The initial zeroth-order approximations for HAM are and are given by the solutions to the Bloch equations in the absence of exchange for initial equilibrium magnetization . The initial approximations are calculated by using the Euler angles determined from method 1 described above. Thus, in the following:
and the system of first-order deformation equations yield the following:
The solution to this equation yields , and is given by the following:
The first-order approximation of magnetization during the pulse is given by the following:
At this level of approximation, relaxation of transverse magnetization depends simply on R2, while relaxation of Mz(t) depends on the average z magnetization during the pulse (calculated in the absence of relaxation). For macromolecules, R2≫R1 typically, and the term proportional to is small.
The nth-order deformation equation leads to the following expression for n>1:
If , the above recursive expressions can be written compactly as follows:
For a rectangular pulse applied to equilibrium magnetization (with magnitude set to unity for convenience), the initial approximations are as follows:
In this case, Γ(t)=Γ and the series of approximations given in Eq. (67) can be summed to give the following:
and this yields identical results as a direct integration of the Bloch equations. Equations (65), (67), and (70) explicitly show the effect of relaxation as a perturbation of the evolution of magnetization in the absence of relaxation.
Figure 5 shows the magnetization components for rectangular 90∘, 180∘, 270∘, and 360∘ nominal on-resonance pulses in the absence and presence of relaxation. Calculations were performed in the absence of relaxation using Eq. (68) and in the presence of relaxation using the HAM approximations, i.e., Eqs. (69) and (70). The first-order HAM approximation is surprisingly accurate for moderate values of R2, except for cases in which , such as the on-resonance 360∘ pulse. The above expressions display the fundamental dependence of relaxation during a pulse applied to equilibrium magnetization on the time-average z magnetization.
Fast, accurate methods for solving differential equations have widespread application in NMR spectroscopy. The present work has illustrated the homotopy analysis method (Liao, 2012) for approximating solutions for differential equations by application to the Riccati differential equation for the Euler angle representation of radiofrequency pulse shapes and to solutions of the Bloch equations incorporating relaxation. The freedom to select the linear operator, lowest-order approximate solution, convergence control parameter, and auxiliary function is powerful in obtaining series solutions that are highly accurate for low orders of approximation and efficient to calculate or that provide qualitatively convenient series allowing physical insight. It can be expected that homotopy analysis method will find other applications in NMR spectroscopy.
RMarkdown and bibtex files are provided in the Supplement. The RMarkdown file contains the Python code for all calculations described in the paper.
The supplement related to this article is available online at: https://doi.org/10.5194/mr-2-175-2021-supplement.
AGP conceived the project. The theoretical derivations, numerical calculations, and writing of the paper were performed by TC and AGP.
The authors declare that they have no conflict of interest.
This article is part of the special issue “Robert Kaptein Festschrift”. It is not associated with a conference.
Some of the work presented here was conducted at the Center on Macromolecular Dynamics by NMR Spectroscopy located at the New York Structural Biology Center. Arthur G. Palmer is a member of the New York Structural Biology Center. This paper is dedicated to Robert Kaptein on the occasion of his 80th birthday.
This research has been supported by the National Institutes of Health (NIH; grant nos. R35GM130398 and P41GM118302).
This paper was edited by Malcolm Levitt and reviewed by Fabien Ferrage and one anonymous referee.
Arfken, G. B., Weber, H. J., and Harris, F. E.: Mathematical Methods for Physicists, Academic Press, Boston, MA, 7th edn., 2013. a
Asami, S., Kallies, W., Günther, J. C., Stavropoulou, M., Glaser, S. J., and Sattler, M.: Ultrashort broadband cooperative pulses for multidimensional biomolecular NMR experiments, Angew. Chem. Int. Ed., 57, 14498–14502, 2018. a
Cavanagh, J., Fairbrother, W. J., Palmer, A. G., Rance, M., and Skelton, N. J.: Protein NMR Spectroscopy: Principles and Practice, Academic Press, San Diego, CA, 2nd edn., 2007. a
Emsley, L. and Bodenhausen, G.: Optimization of shaped selective pulses for NMR using a quaternion description of their overall propagators, J. Magn. Reson., 97, 135–148, 1992. a, b, c
Geen, H. and Freeman, R.: Band-selective radiofrequency pulses, J. Magn. Reson., 93, 93–141, 1991. a, b, c, d, e, f
Gershenzon, N. I., Skinner, T. E., Brutscher, B., Khaneja, N., Nimbalkar, M., Luy, B., and Glaser, S. J.: Linear phase slope in pulse design: Application to coherence transfer, J. Magn. Reson., 192, 235–243, 2008. a
Kupc̆e, E. and Freeman, R.: Adiabatic pulses for wideband inversion and broadband decoupling, J. Magn. Reson., Ser. A, 115, 273–276, 1995. a, b
Kupc̆e, E., Boyd, J., and Campbell, I.: Short selective pulses for biochemical applications, J. Magn. Reson., Ser. B, 106, 300–303, 1995. a, b
Lescop, E., Kern, T., and Brutscher, B.: Guidelines for the use of band-selective radiofrequency pulses in hetero-nuclear NMR: Example of longitudinal-relaxation-enhanced BEST-type 1H-15N correlation experiments, J. Magn. Reson., 203, 190–198, 2010. a
Li, J. S., Ruths, J., Yu, T.-Y., Arthanari, H., and Wagner, G.: Optimal pulse design in quantum control: A unified computational method, P. Natl. Acad. Sci. USA, 108, 1879–1884, 2011. a
Li, Y., Rance, M., and Palmer, A. G.: Rotation operator propagators for time-varying radiofrequency pulses in NMR spectroscopy: applications to shaped pulses and pulse trains, J. Magn. Reson., 248, 105–114, 2014. a, b, c, d, e
Liao, S.: The proposed homotopy analysis technique for the solution of nonlinear problems, PhD thesis, Shanghai Jiao Tong University, Shanghai, China, 1992. a
Liao, S.: Homotopy Analysis Method in Nonlinear Differential Equations, Springer, Berlin, Heidelberg, 2012. a, b, c, d, e, f
Nimbalkar, M., Luy, B., Skinner, T. E., Neves, J. L., Gershenzon, N. I., Kobzar, K., Bermel, W., and Glaser, S. J.: The Fantastic Four: A plug 'n' play set of optimal control pulses for enhancing NMR spectroscopy, J. Magn. Reson., 228, 16–31, 2013. a
Rourke, D. E.: Solutions and linearization of the nonlinear dynamics of radiation damping, Concepts Magn. Reson., 14, 112–129, 2002. a, b
Silver, M. S., Joseph, R. I., and Hoult, D. I.: Selective spin inversion in nuclear magnetic resonance and coherent optics through an exact solution of the Bloch-Riccati equation, Phys. Rev. A, 31, 2753–2755, 1985. a
Siminovitch, D. J.: Comment on “Rotation operator approach to spin dynamics and the Euler geometric equations” [J. Chem. Phys. 101, 6424 (1994)], J. Chem. Phys., 103, 2766–2768, 1995. a
Siminovitch, D. J.: Rotations in NMR: Part I. Euler-Rodrigues parameters and quaternions, Concepts Magn. Reson., 9, 149–171, 1997a. a
Siminovitch, D. J.: Rotations in NMR: Part II. Applications of the Euler-Rodrigues parameters, Concepts Magn. Reson., 9, 211–225, 1997b. a
Warren, W. S.: Effects of arbitrary laser or NMR pulse shapes on population-inversion and coherence, J. Chem. Phys., 81, 5437–5448, 1984. a
Zhou, J., Ye, C., and Sanctuary, B. C.: Rotation operator approach to spin dynamics and the Euler geometric equations, J. Chem. Phys., 101, 6424–6429, 1994. a, b, c