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.

1 Introduction

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):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\mathbf{U}& =e^{-i\mathit{\gamma }\left({\mathit{\tau }}_{\mathrm{p}}\right)I_z}{e}^{-i\mathit{\beta }\left({\mathit{\tau }}_{\mathrm{p}}\right)I_x}{e}^{-i\mathit{\alpha }\left({\mathit{\tau }}_{\mathrm{p}}\right)I_z}\\ & =\left[\begin{array}{rr}{e}^{-i{\mathit{\chi }}^{+}\left({\mathit{\tau }}_p\right)}\cos \left(\frac{\mathit{\beta }\left({\mathit{\tau }}_p\right)}{\mathrm{2}}\right)& -i{e}^{i{\mathit{\chi }}^{-}\left({\mathit{\tau }}_p\right)}\sin \left(\frac{\mathit{\beta }\left({\mathit{\tau }}_p\right)}{\mathrm{2}}\right)\\ -i{e}^{-i{\mathit{\chi }}^{-}\left({\mathit{\tau }}_p\right)}\sin \left(\frac{\mathit{\beta }\left({\mathit{\tau }}_p\right)}{\mathrm{2}}\right)& e^{i{\mathit{\chi }}^{+}\left({\mathit{\tau }}_p\right)}\cos \left(\frac{\mathit{\beta }\left({\mathit{\tau }}_p\right)}{\mathrm{2}}\right)\end{array}\right]\\ & =\prod _{k=\mathrm{1}}^K{\mathbf{U}}_k,\end{array}\end{array}$$

in which χ^±(τ_p) is given as follows:

$$\begin{array}{}\text{(2)}& {\mathit{\chi }}^{\pm }\left({\mathit{\tau }}_p\right)=\left[\mathit{\alpha }\right({\mathit{\tau }}_p)\pm \mathit{\gamma }({\mathit{\tau }}_p\left)\right]/\mathrm{2},\end{array}$$

I_k 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:

$$\begin{array}{}\text{(3)}& {\mathbf{U}}_k=\left[\begin{array}{rr}{\mathit{\kappa }}^{*}\left({\mathit{\tau }}_p\right)& -i{\mathit{\eta }}^{*}\left({\mathit{\tau }}_p\right)\\ -i\mathit{\eta }\left({\mathit{\tau }}_p\right)& \mathit{\kappa }\left({\mathit{\tau }}_p\right)\end{array}\right],\end{array}$$

in which κ(τ_p) and η(τ_p) are given by the following:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& \mathit{\kappa }\left({\mathit{\tau }}_p\right)=\cos ({\mathit{\omega }}_{ek}\mathrm{\Delta }{t}_k/\mathrm{2})+i\cos {\mathit{\theta }}_k\sin ({\mathit{\omega }}_{ek}\mathrm{\Delta }{t}_k/\mathrm{2})\\ & \mathit{\eta }\left({\mathit{\tau }}_p\right)=e^{i{\mathit{\varphi }}_k}\sin {\mathit{\theta }}_k\sin ({\mathit{\omega }}_{ek}\mathrm{\Delta }{t}_k/\mathrm{2}).\end{array}\end{array}$$

In Eq. (4), ω_1k, ϕ_k, and Δt_k 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); ${\mathit{\omega }}_{ek}=({\mathit{\omega }}_{\mathrm{1}k}^{\mathrm{2}}+{\mathrm{\Omega }}_k^{\mathrm{2}}{)}^{\mathrm{1}/\mathrm{2}}$ is the effective field; and ${\mathit{\theta }}_k={\tan }^{-\mathrm{1}}({\mathit{\omega }}_{\mathrm{1}k}/{\mathrm{\Omega }}_k)$ 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):

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\mathrm{d}f\left(t\right)}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\omega }}^{+}\left(t\right)f^{\mathrm{2}}\left(t\right)+i\mathrm{\Omega }\left(t\right)f\left(t\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\omega }}^{-}\left(t\right),\end{array}$$

in which the following applies:

$$\begin{array}{}\text{(6)}& f\left(t\right)=\tan \left({\displaystyle \frac{\mathit{\beta }\left(t\right)}{\mathrm{2}}}\right)e^{i\mathit{\gamma }\left(t\right)},\end{array}$$

where ${\mathit{\omega }}^{\pm }\left(t\right)={\mathit{\omega }}_x\left(t\right)\pm i{\mathit{\omega }}_y\left(t\right)$ 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:

$$\begin{array}{}\text{(7)}& \mathit{\alpha }\left({\mathit{\tau }}_{\mathrm{p}}\right)=\underset{\mathrm{0}}{\overset{{\mathit{\tau }}_{\mathrm{p}}}{\int }}\mathrm{d}t\mathit{\left\{}{\mathit{\omega }}_x\right(t\left)\sin \right[\mathit{\gamma }\left(t\right)]-{\mathit{\omega }}_y(t\left)\cos \right[\mathit{\gamma }\left(t\right)\left]\mathit{\right\}}/\sin \left[\mathit{\beta }\right(t\left)\right].\end{array}$$

The Riccati equation can be transformed into a second-order differential equation as follows:

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}y\left(t\right)}{\mathrm{d}{t}^{\mathrm{2}}}}-\left[{\displaystyle \frac{\mathrm{d}\ln \left[{\mathit{\omega }}^{-}\left(t\right)\right]}{\mathrm{d}t}}+i\mathrm{\Omega }\left(t\right)\right]{\displaystyle \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left|\mathit{\omega }\right(t\left){|}^{\mathrm{2}}y\right(t)=\mathrm{0},\end{array}$$

by using the following definition:

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\mathrm{d}\ln \left[y\left(t\right)\right]}{\mathrm{d}t}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\omega }}^{-}\left(t\right)f\left(t\right).\end{array}$$

A more compact form is obtained by defining the following:

$$\begin{array}{}\text{(10)}& {\widehat{\mathit{\omega }}}^{-}\left(t\right)=\exp \left[i\underset{\mathrm{0}}{\overset{t}{\int }}\mathrm{\Omega }\left(t^{\prime }\right)\mathrm{d}{t}^{\prime }\right]{\mathit{\omega }}^{-}\left(t\right),\end{array}$$

to yield the following:

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}y\left(t\right)}{\mathrm{d}{t}^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{d}\ln \left[{\widehat{\mathit{\omega }}}^{-}\left(t\right)\right]}{\mathrm{d}t}}{\displaystyle \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left|\widehat{\mathit{\omega }}\right(t\left){|}^{\mathrm{2}}y\right(t)=\mathrm{0}.\end{array}$$

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.

2 Theory

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 $q\in [\mathrm{0},\mathrm{1}]$.

This relationship is established by constructing the homotopy as follows (Liao, 2012):

$$\begin{array}{}\text{(12)}& \begin{array}{rl}\mathcal{H}\left[\mathrm{\Phi }(t;q):q\right]=(\mathrm{1}-q)& \mathcal{L}\left[\mathrm{\Phi }(t;q)-y_{\mathrm{0}}\left(t\right)\right]\\ & -q{c}_{\mathrm{0}}H\left(t\right)\mathcal{N}\left[\mathrm{\Phi }(t;q)\right],\end{array}\end{array}$$

in which ℒ[ ] is a linear (differential) operator, and 𝒩[ ] is an (nonlinear differential) operator satisfying the following:

$$\begin{array}{}\text{(13)}& {\displaystyle }\mathcal{L}\left[\mathrm{0}\right]& {\displaystyle }=\mathrm{0}\text{(14)}& {\displaystyle }\mathcal{N}\left[y\left(t\right)\right]& {\displaystyle }=\mathrm{0},\end{array}$$

where y₀(t) is an initial approximation for the desired solution y(t), c₀≠0 is a convergence control parameter, and H(t)≠0 is an auxiliary function (vide infra). When q=0, the homotopy becomes as follows:

$$\begin{array}{}\text{(15)}& \mathcal{H}\left[\mathrm{\Phi }(t;\mathrm{0}):\mathrm{0}\right]=\mathcal{L}\left[\mathrm{\Phi }(t;\mathrm{0})-y_{\mathrm{0}}\left(t\right)\right].\end{array}$$

Therefore, when $\mathcal{H}\left[\mathrm{\Phi }(t;\mathrm{0}):\mathrm{0}\right]=\mathrm{0}$, Eq. (13) requires $\mathrm{\Phi }(t,\mathrm{0})=y_{\mathrm{0}}\left(t\right)$. Similarly, when q=1, the homotopy becomes as follows:

$$\begin{array}{}\text{(16)}& \mathcal{H}\left[\mathrm{\Phi }(t;\mathrm{1}):\mathrm{1}\right]=-c_{\mathrm{0}}H\left(t\right)\mathcal{N}\left[\mathrm{\Phi }(t;\mathrm{1})\right].\end{array}$$

Therefore, when $\mathcal{H}\left[\mathrm{\Phi }(t;\mathrm{1}):\mathrm{1}\right]=\mathrm{0}$, Eq. (14) requires $\mathrm{\Phi }(t,\mathrm{1})=y\left(t\right)$. Stated more succinctly, as q increases from 0→1, Φ(t;q) deforms from the initial approximation y₀(t) to the exact solution y(t). To proceed, the Maclaurin series for Φ(t;q) is assumed to exist as follows:

$$\begin{array}{}\text{(17)}& \mathrm{\Phi }(t;q)=\sum _{n=\mathrm{0}}^{\mathrm{\infty }}{y}_n\left(t\right)q^n,\end{array}$$

in which y_n(t) is given by the following:

$$\begin{array}{}\text{(18)}& y_n\left(t\right)={\displaystyle \frac{\mathrm{1}}{n\mathrm{!}}}{\displaystyle \frac{{\mathrm{d}}^n\mathrm{\Phi }(t;q)}{\mathrm{d}{q}^n}}{|}_{q=\mathrm{0}}.\end{array}$$

Conditions concerning convergence of the series are discussed by Liao (2012). Equation (17) has the desired property $\mathrm{\Phi }(t;\mathrm{0})=y_{\mathrm{0}}\left(t\right)$ and yields the following:

$$\begin{array}{}\text{(19)}& \mathrm{\Phi }(t;\mathrm{1})=y\left(t\right)=\sum _{n=\mathrm{0}}^{\mathrm{\infty }}{y}_n\left(t\right).\end{array}$$

HAM then consists of successively determining y_n(t), beginning with the initial approximation y₀(t), until y(t) is approximated to the desired accuracy. The choices of ℒ[ ], y₀(t), c₀, 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:

$$\begin{array}{}\text{(20)}& \mathcal{N}\left[g\left(t\right)\right]={\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}g\left(t\right)}{\mathrm{d}{t}^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{d}\ln \left[{\widehat{\mathit{\omega }}}^{-}\left(t\right)\right]}{\mathrm{d}t}}{\displaystyle \frac{\mathrm{d}g\left(t\right)}{\mathrm{d}t}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left|\widehat{\mathit{\omega }}\right(t){|}^{\mathrm{2}},\end{array}$$

in which g(t) is an arbitrary function. The linear operator is chosen to be the following:

$$\begin{array}{}\text{(21)}& \mathcal{L}\left[g\left(t\right)\right]={\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}g\left(t\right)}{\mathrm{d}{t}^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{d}\ln \left[{\widehat{\mathit{\omega }}}^{-}\left(t\right)\right]}{\mathrm{d}t}}{\displaystyle \frac{\mathrm{d}g\left(t\right)}{\mathrm{d}t}},\end{array}$$

and the initial approximation is y₀(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):

$$\begin{array}{}\text{(22)}& (\mathrm{1}-q)\mathcal{L}\left[\mathrm{\Phi }(t;q)-y_{\mathrm{0}}\left(t\right)\right]=q{c}_{\mathrm{0}}\mathcal{N}\left[\mathrm{\Phi }(t;q)\right].\end{array}$$

The derivative of Eq. (22), with respect to q, yields the first-order deformation equation as follows:

$$\begin{array}{}\text{(23)}& \begin{array}{rl}& -\mathcal{L}\left[\mathrm{\Phi }(t;q)-y_{\mathrm{0}}\left(t\right)\right]+(\mathrm{1}-q)\mathcal{L}\left[{\displaystyle \frac{\mathrm{d}\mathrm{\Phi }(t;q)}{\mathrm{d}q}}\right]\\ & =c_{\mathrm{0}}\mathcal{N}\left[\mathrm{\Phi }(t;q)\right]+q{c}_{\mathrm{0}}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}q}}\mathcal{N}\left[\mathrm{\Phi }(t;q)\right].\end{array}\end{array}$$

The limit q→0 gives the following:

$$\begin{array}{ll}{\displaystyle }-\mathcal{L}\left[\mathrm{\Phi }(t;\mathrm{0})-y_{\mathrm{0}}\left(t\right)\right]+\mathcal{L}\left[{\displaystyle \frac{\mathrm{d}\mathrm{\Phi }(t;q)}{\mathrm{d}q}}{|}_{q=\mathrm{0}}\right]& {\displaystyle }=c_{\mathrm{0}}\mathcal{N}\left[\mathrm{\Phi }(t;\mathrm{0})\right]\\ \text{(24)}& {\displaystyle }\mathcal{L}\left[y_{\mathrm{1}}\left(t\right)\right]& {\displaystyle }=c_{\mathrm{0}}\mathcal{N}\left[y_{\mathrm{0}}\left(t\right)\right],\end{array}$$

in which the second line is obtained using $\mathrm{\Phi }(t;\mathrm{0})=y_{\mathrm{0}}\left(t\right)$ and Eq. (18). Substituting for 𝒩[ ], ℒ[ ], and y₀(t) yields the following:

$$\begin{array}{}\text{(25)}& \begin{array}{rl}& {\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}{y}_{\mathrm{1}}\left(t\right)}{\mathrm{d}{t}^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{d}\ln \left[{\widehat{\mathit{\omega }}}^{-}\left(t\right)\right]}{\mathrm{d}t}}{\displaystyle \frac{\mathrm{d}{y}_{\mathrm{1}}\left(t\right)}{\mathrm{d}t}}\\ & =c_{\mathrm{0}}\left({\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}{y}_{\mathrm{0}}\left(t\right)}{\mathrm{d}{t}^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{d}\ln \left[{\widehat{\mathit{\omega }}}^{-}\left(t\right)\right]}{\mathrm{d}t}}{\displaystyle \frac{\mathrm{d}{y}_{\mathrm{0}}\left(t\right)}{\mathrm{d}t}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left|\widehat{\mathit{\omega }}\right(t\left){|}^{\mathrm{2}}{y}_{\mathrm{0}}\right(t)\right)\\ & ={\displaystyle \frac{c_{\mathrm{0}}}{\mathrm{4}}}\left|\widehat{\mathit{\omega }}\right(t){|}^{\mathrm{2}},\end{array}\end{array}$$

in which the final line is obtained using $\mathrm{d}{y}_{\mathrm{0}}\left(t\right)/\mathrm{d}t=\mathrm{0}$. This differential equation does not contain a term proportional to y₁(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:

$$\begin{array}{}\text{(26)}& y_{\mathrm{1}}\left(t\right)={\displaystyle \frac{c_{\mathrm{0}}}{\mathrm{4}}}\underset{\mathrm{0}}{\overset{t}{\int }}{\widehat{\mathit{\omega }}}^{-}\left(t^{\prime }\right)\underset{\mathrm{0}}{\overset{t^{\prime }}{\int }}{\widehat{\mathit{\omega }}}^{+}\left(t^{\prime \prime }\right)\mathrm{d}{t}^{\prime \prime }\mathrm{d}{t}^{\prime }.\end{array}$$

The higher-order approximations y_n(t) are obtained in a similar fashion. The nth derivative, with respect to q of Eq. (22), yields the following (for n>1):

$$\begin{array}{}\text{(27)}& \begin{array}{rl}& -n\mathcal{L}\left[{\displaystyle \frac{{\mathrm{d}}^{n-\mathrm{1}}\mathrm{\Phi }(t;q)}{\mathrm{d}{q}^{n-\mathrm{1}}}}\right]+(\mathrm{1}-q)\mathcal{L}\left[{\displaystyle \frac{{\mathrm{d}}^n\mathrm{\Phi }(t;q)}{\mathrm{d}{q}^n}}\right]\\ & =n{c}_{\mathrm{0}}{\displaystyle \frac{{\mathrm{d}}^{n-\mathrm{1}}}{\mathrm{d}{q}^{n-\mathrm{1}}}}\mathcal{N}\left[\mathrm{\Phi }(t;q)\right]+q{c}_{\mathrm{0}}{\displaystyle \frac{{\mathrm{d}}^n}{\mathrm{d}{q}^n}}\mathcal{N}\left[\mathrm{\Phi }(t;q)\right].\end{array}\end{array}$$

Executing the derivatives, taking the limit q→0, and dividing both sides of the equation by n! gives the following:

$$\begin{array}{}\text{(28)}& \begin{array}{rl}{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}{y}_n\left(t\right)}{\mathrm{d}{t}^{\mathrm{2}}}}& -{\displaystyle \frac{\mathrm{d}\ln \left[{\widehat{\mathit{\omega }}}^{-}\left(t\right)\right]}{\mathrm{d}t}}{\displaystyle \frac{\mathrm{d}{y}_n\left(t\right)}{\mathrm{d}t}}\\ & =(c_{\mathrm{0}}+\mathrm{1})\mathit{\{}{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}{y}_{n-\mathrm{1}}\left(t\right)}{\mathrm{d}{t}^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{d}\ln \left[{\widehat{\mathit{\omega }}}^{-}\left(t\right)\right]}{\mathrm{d}t}}{\displaystyle \frac{\mathrm{d}{y}_{n-\mathrm{1}}\left(t\right)}{\mathrm{d}t}}\mathit{\}}\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{c}_{\mathrm{0}}\left|\widehat{\mathit{\omega }}\right(t\left){|}^{\mathrm{2}}{y}_{n-\mathrm{1}}\right(t),\end{array}\end{array}$$

with the solution obtained by the same approach as for Eq. (26), as follows:

$$\begin{array}{}\text{(29)}& \begin{array}{rl}{y}_n\left(t\right)=(c_{\mathrm{0}}+& \mathrm{1}\left)y_{n-\mathrm{1}}\right(t)\\ & +{\displaystyle \frac{c_{\mathrm{0}}}{\mathrm{4}}}\underset{\mathrm{0}}{\overset{t}{\int }}{\widehat{\mathit{\omega }}}^{-}\left(t^{\prime }\right)\underset{\mathrm{0}}{\overset{t^{\prime }}{\int }}{\widehat{\mathit{\omega }}}^{+}\left(t^{\prime \prime }\right)y_{n-\mathrm{1}}\left(t^{\prime \prime }\right)\mathrm{d}{t}^{\prime \prime }\mathrm{d}{t}^{\prime }.\end{array}\end{array}$$

Successive use of Eqs. (26) and (29) allows y(t) and, hence, f(t) to be determined to arbitrary accuracy, as follows:

$$\begin{array}{}\text{(30)}& f\left(t\right)=\left({\displaystyle \frac{-\mathrm{2}}{{\mathit{\omega }}^{-}\left(t\right)}}\right){\displaystyle \frac{\mathrm{d}\ln \left[y\left(t\right)\right]}{\mathrm{d}t}}=\left({\displaystyle \frac{-\mathrm{2}}{{\mathit{\omega }}^{-}\left(t\right)}}\right){\displaystyle \frac{{\sum }_{m=\mathrm{0}}^N\frac{\mathrm{d}{y}_m\left(t\right)}{\mathrm{d}t}}{{\sum }_{n=\mathrm{0}}^N{y}_n\left(t\right)}},\end{array}$$

in which N is the order of approximation. For completeness, the derivatives of Eqs. (26) and (29) are, respectively, as follows:

$$\begin{array}{}\text{(31)}& {\displaystyle \frac{\mathrm{d}{y}_{\mathrm{1}}\left(t\right)}{\mathrm{d}t}}={\displaystyle \frac{c_{\mathrm{0}}}{\mathrm{4}}}{\widehat{\mathit{\omega }}}^{-}\left(t\right)\underset{\mathrm{0}}{\overset{t}{\int }}{\widehat{\mathit{\omega }}}^{+}\left(t^{\prime }\right)\mathrm{d}{t}^{\prime }\text{(32)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{y}_n\left(t\right)}{\mathrm{d}t}}=(c_{\mathrm{0}}& +\mathrm{1}){\displaystyle \frac{\mathrm{d}{y}_{n-\mathrm{1}}\left(t\right)}{\mathrm{d}t}}\\ & +{\displaystyle \frac{c_{\mathrm{0}}}{\mathrm{4}}}{\widehat{\mathit{\omega }}}^{-}\left(t\right)\underset{\mathrm{0}}{\overset{t}{\int }}{\widehat{\mathit{\omega }}}^{+}\left(t^{\prime }\right)y_{n-\mathrm{1}}\left(t^{\prime }\right)\mathrm{d}{t}^{\prime }.\end{array}\end{array}$$

Results obtained using y₀(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 y_n(t) have similarities to the Fourier integrals obtained from average Hamiltonian theory by Warren (1984).

The above choices of ℒ[ ] and y₀(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:

$$\begin{array}{}\text{(33)}& {\displaystyle }\begin{array}{rl}\mathcal{N}\left[g\left(t\right)\right]={\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}g\left(t\right)}{\mathrm{d}{t}^{\mathrm{2}}}}& -\left[{\displaystyle \frac{\mathrm{d}\ln \left[\mathit{\omega }\left(t\right)\right]}{\mathrm{d}t}}+i\mathrm{\Omega }\right]{\displaystyle \frac{\mathrm{d}g\left(t\right)}{\mathrm{d}t}}\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\mathit{\omega }}^{\mathrm{2}}\left(t\right)g\left(t\right)\end{array}\text{(34)}& {\displaystyle }\mathcal{L}\left[g\left(t\right)\right]={\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}g\left(t\right)}{\mathrm{d}{t}^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{d}\ln \left[\mathit{\omega }\left(t\right)\right]}{\mathrm{d}t}}{\displaystyle \frac{\mathrm{d}g\left(t\right)}{\mathrm{d}t}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\mathit{\omega }}^{\mathrm{2}}\left(t\right)g\left(t\right)\text{(35)}& {\displaystyle }{y}_{\mathrm{0}}\left(t\right)=\cos \left[{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\delta }\left(t\right)\right],\text{(36)}& {\displaystyle }\mathit{\delta }\left(t\right)=\underset{\mathrm{0}}{\overset{t}{\int }}\mathit{\omega }\left(t^{\prime }\right)\mathrm{d}{t}^{\prime }.\end{array}$$

This choice of y₀(t) satisfies the following:

$$\begin{array}{}\text{(37)}& {\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}{y}_{\mathrm{0}}\left(t\right)}{\mathrm{d}{t}^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{d}\ln \left[\mathit{\omega }\left(t\right)\right]}{\mathrm{d}t}}{\displaystyle \frac{\mathrm{d}{y}_{\mathrm{0}}\left(t\right)}{\mathrm{d}t}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\mathit{\omega }}^{\mathrm{2}}\left(t\right)y_{\mathrm{0}}\left(t\right)=\mathrm{0},\end{array}$$

and is the exact on-resonance solution for y(t). Consequently, the first-order deformation equation leads to the following:

$$\begin{array}{}\text{(38)}& \begin{array}{rl}& {\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}{y}_{\mathrm{1}}\left(t\right)}{\mathrm{d}{t}^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{d}\ln \left[\mathit{\omega }\left(t\right)\right]}{\mathrm{d}t}}{\displaystyle \frac{\mathrm{d}{y}_{\mathrm{1}}\left(t\right)}{\mathrm{d}t}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\mathit{\omega }}^{\mathrm{2}}\left(t\right)y_{\mathrm{1}}\left(t\right)\\ & =-i{c}_{\mathrm{0}}\mathrm{\Omega }{\displaystyle \frac{\mathrm{d}{y}_{\mathrm{0}}\left(t\right)}{\mathrm{d}t}}.\end{array}\end{array}$$

The solutions to the homogeneous equation (setting the right-hand side to 0) are $y^{\pm }\left(t\right)=e^{\pm i\mathit{\delta }\left(t\right)/\mathrm{2}}$. The method of variation of parameters then gives the inhomogeneous solution as follows:

$$\begin{array}{}\text{(39)}& y_{\mathrm{1}}\left(t\right)=-i{c}_{\mathrm{0}}\mathrm{\Omega }\underset{\mathrm{0}}{\overset{t}{\int }}\sin \left[{\displaystyle \frac{\mathit{\delta }\left(t\right)}{\mathrm{2}}}-{\displaystyle \frac{\mathit{\delta }\left(t^{\prime }\right)}{\mathrm{2}}}\right]\sin \left[{\displaystyle \frac{\mathit{\delta }\left(t^{\prime }\right)}{\mathrm{2}}}\right]\mathrm{d}{t}^{\prime }.\end{array}$$

The nth-order deformation equation for n>1 is as follows:

$$\begin{array}{}\text{(40)}& \begin{array}{rl}{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}{y}_n\left(t\right)}{\mathrm{d}{t}^{\mathrm{2}}}}& -\left[{\displaystyle \frac{\mathrm{d}\ln \left[\mathit{\omega }\left(t\right)\right]}{\mathrm{d}t}}+i\mathrm{\Omega }\right]{\displaystyle \frac{\mathrm{d}{y}_n\left(t\right)}{\mathrm{d}t}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\mathit{\omega }}^{\mathrm{2}}\left(t\right)y_n\left(t\right)\\ & =(\mathrm{1}+c_{\mathrm{0}})\mathit{\{}{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}{y}_{n-\mathrm{1}}\left(t\right)}{\mathrm{d}{t}^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{d}\ln \left[\mathit{\omega }\left(t\right)\right]}{\mathrm{d}t}}{\displaystyle \frac{\mathrm{d}{y}_{n-\mathrm{1}}\left(t\right)}{\mathrm{d}t}}\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\mathit{\omega }}^{\mathrm{2}}\left(t\right)y_{n-\mathrm{1}}\left(t\right)\mathit{\}}-i{c}_{\mathrm{0}}\mathrm{\Omega }{\displaystyle \frac{\mathrm{d}{y}_{n-\mathrm{1}}\left(t\right)}{\mathrm{d}t}},\end{array}\end{array}$$

with the following solution:

$$\begin{array}{}\text{(41)}& \begin{array}{rl}{y}_n\left(t\right)=& (\mathrm{1}+c_{\mathrm{0}})y_{n-\mathrm{1}}\left(t\right)\\ & -i{c}_{\mathrm{0}}\mathrm{\Omega }\underset{\mathrm{0}}{\overset{t}{\int }}{\displaystyle \frac{\mathrm{2}}{\mathit{\omega }\left(t^{\prime }\right)}}\sin \left[{\displaystyle \frac{\mathit{\delta }\left(t\right)}{\mathrm{2}}}-{\displaystyle \frac{\mathit{\delta }\left(t^{\prime }\right)}{\mathrm{2}}}\right]{\displaystyle \frac{\mathrm{d}{y}_{n-\mathrm{1}}\left(t^{\prime }\right)}{\mathrm{d}{t}^{\prime }}}\mathrm{d}{t}^{\prime }.\end{array}\end{array}$$

Each y_n(t) is proportional to Ωⁿ, and these results yield the following power series in Ω for y(t):

$$\begin{array}{}\text{(42)}& y\left(t\right)=y_{\mathrm{0}}\left(t\right)+\sum _{n=\mathrm{1}}^N(\mathrm{2}+c_{\mathrm{0}})y_n\left(t\right),\end{array}$$

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:

$$\begin{array}{}\text{(43)}& {\displaystyle \frac{\mathrm{d}{y}_{\mathrm{1}}\left(t\right)}{\mathrm{d}t}}=-i{c}_{\mathrm{0}}\mathrm{\Omega }{\displaystyle \frac{\mathit{\omega }\left(t\right)}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{t}{\int }}\cos \left[{\displaystyle \frac{\mathit{\delta }\left(t\right)}{\mathrm{2}}}-{\displaystyle \frac{\mathit{\delta }\left(t^{\prime }\right)}{\mathrm{2}}}\right]\sin \left[{\displaystyle \frac{\mathit{\delta }\left(t^{\prime }\right)}{\mathrm{2}}}\right]\mathrm{d}{t}^{\prime }\text{(44)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{y}_n\left(t\right)}{\mathrm{d}t}}=(\mathrm{1}+& c_{\mathrm{0}}){\displaystyle \frac{\mathrm{d}{y}_{n-\mathrm{1}}\left(t\right)}{\mathrm{d}t}}-i{c}_{\mathrm{0}}\mathrm{\Omega }{\displaystyle \frac{\mathit{\omega }\left(t\right)}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{t}{\int }}{\displaystyle \frac{\mathrm{2}}{\mathit{\omega }\left(t^{\prime }\right)}}\\ & \times \cos \left[{\displaystyle \frac{\mathit{\delta }\left(t\right)}{\mathrm{2}}}-{\displaystyle \frac{\mathit{\delta }\left(t^{\prime }\right)}{\mathrm{2}}}\right]{\displaystyle \frac{\mathrm{d}{y}_{n-\mathrm{1}}\left(t^{\prime }\right)}{\mathrm{d}{t}^{\prime }}}\mathrm{d}{t}^{\prime }.\end{array}\end{array}$$

2.1 Methods

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 ${\mathit{\omega }}_{\mathrm{1}}/\left(\mathrm{2}\mathit{\pi }\right)=\mathrm{25}\mathrm{000}$ Hz and an on-resonance 90^∘ pulse length of 10.0 µs or ${\mathit{\omega }}_{\mathrm{1}}/\left(\mathrm{2}\mathit{\pi }\right)=\mathrm{250}$ 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 ${\mathit{\omega }}_{\mathrm{1}}/\left(\mathrm{2}\mathit{\pi }\right)=\mathrm{9000}$ Hz and 90^∘ pulse lengths of 455.2 and 504.9 µs, respectively. REBURP (Geen and Freeman, 1991) pulses were simulated using a maximum ${\mathit{\omega }}_{\mathrm{1}}/\left(\mathrm{2}\mathit{\pi }\right)=\mathrm{10}\mathrm{000}$ Hz and a 180^∘ pulse length of 626.5 µs. WURST-20 (Kupc̆e and Freeman, 1995) pulses were simulated using a maximum ${\mathit{\omega }}_{\mathrm{1}}/\left(\mathrm{2}\mathit{\pi }\right)=\mathrm{9512}$ Hz, a frequency sweep of 50 000 Hz, and a pulse length of 440.0 µs.

Equation (7) can be recast as follows:

$$\begin{array}{}\text{(45)}& \begin{array}{rl}\mathit{\alpha }\left({\mathit{\tau }}_{\mathrm{p}}\right)={\displaystyle \frac{i}{\mathrm{4}}}\underset{\mathrm{0}}{\overset{{\mathit{\tau }}_{\mathrm{p}}}{\int }}& \mathrm{d}t\mathit{\left\{}{\mathit{\omega }}^{+}\right(t\left)f^{*}\right(t)-{\mathit{\omega }}^{-}(t\left)f\right(t\left)\mathit{\right\}}\\ & \times \mathit{\{}\mathrm{1}+|f\left(t\right){|}^{\mathrm{2}}\mathit{\}}/|f\left(t\right){|}^{\mathrm{2}},\end{array}\end{array}$$

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 y₀(t)=1 in method 1 leads to obtaining simple iterative solutions that can be calculated very efficiently. The form of y₀(t) given in Eq. (36) could also be used in Eq. (24) to obtain an alternative expression for y₁(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 y₀(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 ω₁. For example, using Eq. (26) yields the following:

$$\begin{array}{}\text{(46)}& y_{\mathrm{1}}\left(t\right)={\displaystyle \frac{c_{\mathrm{0}}{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{4}{\mathrm{\Omega }}^{\mathrm{2}}}}\left(\mathrm{1}-e^{i\mathrm{\Omega }t}\right)+i{\displaystyle \frac{c_{\mathrm{0}}{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}t}{\mathrm{4}\mathrm{\Omega }}};\end{array}$$

however, analytic calculations of higher order y_n(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 $c_{\mathrm{0}}=-\mathrm{1}$ are nearly indistinguishable from the exact result of Eq. (3) (using τ_p=Δτ_k) over the range of resonance offsets from 0 to $\mathrm{\Omega }/{\mathit{\omega }}_{\mathrm{1}}={\mathrm{15}}^{\mathrm{1}/\mathrm{2}}$. The first-order result provides a highly accurate estimate of γ(τ_p) but overestimates β(τ_p). The role of the convergence control parameter c₀ is illustrated in Fig. 1c, d. A value of $c_{\mathrm{0}}=-\mathrm{0.925}$ was chosen, using Eqs. (46) and (30), to scale the first-order result for β(τ_p) to be equal to $\mathit{\pi }/\mathrm{2}$ 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 c₀ 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 $c_{\mathrm{0}}=-\mathrm{1}$, even to the third order in Ω, is accurate for β(τ_p) only to slightly more than $\mathrm{\Omega }/{\mathit{\omega }}_{\mathrm{1}}=\mathrm{1}$. When $c_{\mathrm{0}}=-\mathrm{0.925}$, the third-order power series has improved accuracy for resonance offsets up to nearly $\mathrm{\Omega }/{\mathit{\omega }}_{\mathrm{1}}=\mathrm{2}$. 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 $\mathrm{\Omega }/{\mathit{\omega }}_{\mathrm{1}}=\mathrm{3.5}$ 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.

Figure 1 HAM approximations for on-resonance 90^∘ rectangular pulse with ${\mathit{\omega }}_{\mathrm{1}}/\left(\mathrm{2}\mathit{\pi }\right)=\mathrm{25}\mathrm{000}$ Hz. Exact calculation of Euler angles β(τ_p) and γ(τ_p) (black line). For a rectangular pulse, α(τ_p)=γ(τ_p). First-order (blue dotted line), second-order (reddish-purple dashed line), and third-order (orange dash–dot–dotted line) HAM results using method 1. The third-order result is found using the power series of method 2 (green dash-dotted line). Results are shown for (a, b) $c_{\mathrm{0}}=-\mathrm{1}$ and (c, d) $c_{\mathrm{0}}=-\mathrm{0.925}$. The exact, second-order HAM and third-order HAM curves for method 1 are virtually indistinguishable.

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 $c_{\mathrm{0}}=-\mathrm{1}$, the fifth-order approximation is extremely accurate compared with numerical calculations using Eqs. (1)–(3) (Fig. 2a–c). With $c_{\mathrm{0}}=-\mathrm{0.925}$ (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.

Figure 2 HAM approximations for 90^∘ EBURP-2 pulse. Numerical calculation of Euler angles α(τ_p), β(τ_p), and γ(τ_p) using Eqs. (1)–(3) (black line). First-order (blue dotted line), second-order (reddish-purple dashed line), third-order (green dash-dotted line), and fifth-order (orange dash–dot–dotted line) HAM results, using method 1. Results are shown for (a, b, c) $c_{\mathrm{0}}=-\mathrm{1}$ and (d, e, f) $c_{\mathrm{0}}=-\mathrm{0.925}$. The numerical calculation and fifth-order HAM curves are nearly indistinguishable.

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 $c_{\mathrm{0}}=-\mathrm{1}$ and $c_{\mathrm{0}}=-\mathrm{0.925}$. The choice of $c_{\mathrm{0}}=-\mathrm{0.925}$ has a remarkable effect in terms of increasing the accuracy the fifth-order approximation to nearly that of the seventh-order result.

Figure 3 HAM approximations for 90^∘ Q5 pulse. Numerical calculation of Euler angles α(τ_p), β(τ_p), and γ(τ_p) using Eqs. (1)–(3) (black line). Fifth-order (blue dotted line) and seventh-order (orange dash–dot–dotted line) HAM results, using method 1. Results are shown for (a, b, c) $c_{\mathrm{0}}=-\mathrm{1}$ and (d, e, f) $c_{\mathrm{0}}=-\mathrm{0.925}$. The numerical calculation and seventh-order HAM curves are nearly indistinguishable.

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 $c_{\mathrm{0}}=-\mathrm{1}$ and $c_{\mathrm{0}}=-\mathrm{0.925}$. The third-order approximation also is highly accurate when $c_{\mathrm{0}}=-\mathrm{0.925}$. 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 $c_{\mathrm{0}}=-\mathrm{1}$, or ninth-order, with $c_{\mathrm{0}}=-\mathrm{0.925}$, results are highly accurate.

Figure 4 HAM approximations for (a, b) REBURP and (c, d) WURST-20 inversion pulses. Numerical calculation of Euler angle β(τ_p) using Eqs. (1)–(3) (black line). (a, b) First-order (blue dotted line), second-order (reddish-purple dashed line), third-order (green dash-dotted line), and fifth-order (orange dash–dot–dotted line) HAM results, using method 1. (c, d) Fifth-order (blue dotted line), seventh-order (reddish-purple dashed line), ninth-order (green dash-dotted line), and 11th-order (orange dash–dot–dotted line) HAM results, using method 1. Results are shown for (a, c) $c_{\mathrm{0}}=-\mathrm{1}$ and (b, d) $c_{\mathrm{0}}=-\mathrm{0.925}$. The numerical calculation and (a, b) fifth-order and (c, d) 11th-order HAM curves are nearly indistinguishable.

Method 2 yields a power series for y(t). If $c_{\mathrm{0}}=-\mathrm{1}$, the resulting series is identical to the power series expansion obtained from perturbation theory (Li et al., 2014), while $c_{\mathrm{0}}=-\mathrm{0.925}$ 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 $c_{\mathrm{0}}=-\mathrm{1}$, yields the following:

$$\begin{array}{ll}{\displaystyle }f\left(t\right)=& {\displaystyle \frac{\sin \left[\frac{\mathit{\delta }\left(t\right)}{\mathrm{2}}\right]+i\mathrm{\Omega }\underset{\mathrm{0}}{\overset{t}{\int }}\cos \left[\frac{\mathit{\delta }\left(t\right)}{\mathrm{2}}-\frac{\mathit{\delta }\left(t^{\prime }\right)}{\mathrm{2}}\right]\sin \left[\frac{\mathit{\delta }\left(t^{\prime }\right)}{\mathrm{2}}\right]\mathrm{d}{t}^{\prime }}{\cos \left[\frac{\mathit{\delta }\left(t\right)}{\mathrm{2}}\right]-i\mathrm{\Omega }\underset{\mathrm{0}}{\overset{t}{\int }}\sin \left[\frac{\mathit{\delta }\left(t^{\prime }\right)}{\mathrm{2}}-\frac{\mathit{\delta }\left(t^{\prime }\right)}{\mathrm{2}}\right]\sin \left[\frac{\mathit{\delta }\left(t^{\prime }\right)}{\mathrm{2}}\right]\mathrm{d}{t}^{\prime }}}\\ \text{(47)}& {\displaystyle }& {\displaystyle }\approx \tan \left[{\displaystyle \frac{\mathit{\delta }\left(t\right)}{\mathrm{2}}}\right]\left(\mathrm{1}+i{\displaystyle \frac{\mathrm{\Omega }}{\sin \left[\mathit{\delta }\left(t\right)\right]}}\underset{\mathrm{0}}{\overset{t}{\int }}\sin \left[\mathit{\delta }\left(t^{\prime }\right)\right]\mathrm{d}{t}^{\prime }\right),\end{array}$$

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:

$$\begin{array}{}\text{(48)}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\mathit{\alpha }}={\displaystyle \frac{\mathrm{1}}{\sin \left[\mathit{\delta }\left({\mathit{\tau }}_{\mathrm{p}}\right)\right]}}\underset{\mathrm{0}}{\overset{{\mathit{\tau }}_{\mathrm{p}}}{\int }}\sin \left[\mathit{\delta }({\mathit{\tau }}_{\mathrm{p}}-t^{\prime })\right]\mathrm{d}{t}^{\prime }\text{(49)}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\mathit{\gamma }}={\displaystyle \frac{\mathrm{1}}{\sin \left[\mathit{\delta }\left({\mathit{\tau }}_{\mathrm{p}}\right)\right]}}\underset{\mathrm{0}}{\overset{{\mathit{\tau }}_{\mathrm{p}}}{\int }}\sin \left[\mathit{\delta }\left(t^{\prime }\right)\right]\mathrm{d}{t}^{\prime }.\end{array}$$

For a 90^∘ pulse, the above equations can be written compactly as follows:

$$\begin{array}{}\text{(50)}& {\mathit{\tau }}_{\mathit{\alpha }}+i{\mathit{\tau }}_{\mathit{\gamma }}=\underset{\mathrm{0}}{\overset{{\mathit{\tau }}_{\mathrm{p}}}{\int }}{e}^{i\mathit{\delta }\left(t^{\prime }\right)}\mathrm{d}{t}^{\prime }.\end{array}$$

The ratios ${\mathit{\tau }}_{\mathit{\alpha }}/{\mathit{\tau }}_{\mathrm{p}}$ and ${\mathit{\tau }}_{\mathit{\gamma }}/{\mathit{\tau }}_{\mathrm{p}}$ 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:

$$\begin{array}{}\text{(51)}& {\displaystyle }& {\displaystyle }\mathcal{N}\left[g\right(t\left)\right]={\displaystyle \frac{\mathrm{d}g\left(t\right)}{\mathrm{d}t}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\omega }\left(t\right)g^{\mathrm{2}}\left(t\right)-i\mathrm{\Omega }g\left(t\right)-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\omega }\left(t\right)\text{(52)}& {\displaystyle }& {\displaystyle }\mathcal{L}\left[g\right(t\left)\right]={\displaystyle \frac{\mathrm{d}g\left(t\right)}{\mathrm{d}t}}-i\mathrm{\Omega }g\left(t\right)\text{(53)}& {\displaystyle }& {\displaystyle }{f}_{\mathrm{0}}\left(t\right)=\tan \left[{\displaystyle \frac{\mathit{\delta }\left(t\right)}{\mathrm{2}}}\right],\end{array}$$

in which f₀(t) is the exact solution for Ω=0, yields the following series solution:

$$\begin{array}{}\text{(54)}& f\left(t\right)=\tan \left[{\displaystyle \frac{\mathit{\delta }\left(t\right)}{\mathrm{2}}}\right]+\sum _{n=\mathrm{1}}^N{f}_n\left(t\right).\end{array}$$

The result obtained from the first-order deformation equation is as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\mathrm{d}{f}_{\mathrm{1}}\left(t\right)}{\mathrm{d}t}}-i\mathrm{\Omega }{f}_{\mathrm{1}}\left(t\right)=-i{c}_{\mathrm{0}}\mathrm{\Omega }{f}_{\mathrm{0}}\left(t\right)\\ \text{(55)}& {\displaystyle }& {\displaystyle }{f}_{\mathrm{1}}\left(t\right)=-i{c}_{\mathrm{0}}\mathrm{\Omega }{e}^{i\mathrm{\Omega }t}\underset{\mathrm{0}}{\overset{t}{\int }}{e}^{-i\mathrm{\Omega }t}\tan \left[{\displaystyle \frac{\mathit{\delta }\left(t^{\prime }\right)}{\mathrm{2}}}\right]\mathrm{d}{t}^{\prime }.\end{array}$$

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 Φ²(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:

$$\begin{array}{}\text{(56)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}{f}_{\mathrm{2}}\left(t\right)}{\mathrm{d}t}}& -i\mathrm{\Omega }{f}_{\mathrm{2}}\left(t\right)\\ & =c_{\mathrm{0}}\mathit{\{}{\displaystyle \frac{\mathrm{d}{f}_{\mathrm{1}}\left(t\right)}{\mathrm{d}t}}-i\mathrm{\Omega }{f}_{\mathrm{1}}(t)-\mathit{\omega }(t\left)f_{\mathrm{0}}\right(t\left)f_{\mathrm{1}}\right(t\left)\mathit{\right\}}.\end{array}\end{array}$$

$$\begin{array}{}\text{(57)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}{f}_{\mathrm{3}}\left(t\right)}{\mathrm{d}t}}& -i\mathrm{\Omega }{f}_{\mathrm{3}}\left(t\right)\\ & =c_{\mathrm{0}}\mathit{\{}{\displaystyle \frac{\mathrm{d}{f}_{\mathrm{2}}\left(t\right)}{\mathrm{d}t}}-i\mathrm{\Omega }{f}_{\mathrm{2}}(t)-\mathrm{2}\mathit{\omega }(t\left)f_{\mathrm{0}}\right(t\left)f_{\mathrm{2}}\right(t)-\mathit{\omega }(t\left)f_{\mathrm{1}}^{\mathrm{2}}\right(t\left)\mathit{\right\}}.\end{array}\end{array}$$

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:

$$\begin{array}{}\text{(58)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\widehat{\mathit{M}}\left(t\right)=-\mathbf{\Gamma }\left(t\right)\widehat{\mathit{M}}\left(t\right)+\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ e^{R_{\mathrm{2}}t}{R}_{\mathrm{1}}{M}_{\mathrm{0}}\end{array}\right],\end{array}$$

$$\begin{array}{}\text{(59)}& \widehat{\mathit{M}}\left(t\right)=\left[\begin{array}{c}{\widehat{M}}_x\left(t\right)\\ {\widehat{M}}_y\left(t\right)\\ {\widehat{M}}_z\left(t\right)\end{array}\right],\end{array}$$

$$\begin{array}{}\text{(60)}& \mathbf{\Gamma }\left(t\right)=\left[\begin{array}{ccc}\mathrm{0}& \mathrm{\Omega }& \mathrm{0}\\ -\mathrm{\Omega }& \mathrm{0}& {\mathit{\omega }}_x\left(t\right)\\ \mathrm{0}& -{\mathit{\omega }}_x\left(t\right)& -(R_{\mathrm{2}}-R_{\mathrm{1}})\end{array}\right],\end{array}$$

in which $M_k\left(t\right)=e^{-R_{\mathrm{2}}t}{\widehat{M}}_k\left(t\right)$ are the Cartesian components of the magnetization, and M₀ is the equilibrium magnetization. The use of transformed variables $\widehat{\mathit{M}}\left(t\right)$ rather than M(t) simplifies the following discussion. The linear operator is chosen as $\mathcal{L}\left[\mathit{g}\right(t\left)\right]=\mathrm{d}\mathit{g}\left(t\right)/\mathrm{d}t$, in which $\mathit{g}\left(t\right)=\left[g_x\right(t),g_y(t),g_z(t){]}^T$ is an arbitrary vector function. The nonlinear operator is as follows:

$$\begin{array}{}\text{(61)}& \mathcal{N}\left[\mathit{g}\right(t\left)\right]={\displaystyle \frac{\mathrm{d}\mathit{g}\left(t\right)}{\mathrm{d}t}}+\mathbf{\Gamma }\left(t\right)\mathit{g}\left(t\right)-\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ e^{R_{\mathrm{2}}t}{R}_{\mathrm{1}}{M}_{\mathrm{0}}\end{array}\right].\end{array}$$

The initial zeroth-order approximations for HAM are ${\widehat{\mathit{M}}}_{\mathrm{0}}\left(t\right)={\mathit{M}}_{\mathrm{0}}\left(t\right)$ and are given by the solutions to the Bloch equations in the absence of exchange for initial equilibrium magnetization $[\mathrm{0},\mathrm{0},M_{\mathrm{0}}{]}^T$. The initial approximations are calculated by using the Euler angles determined from method 1 described above. Thus, in the following:

$$\begin{array}{}\text{(62)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\widehat{\mathit{M}}}_{\mathrm{0}}\left(t\right)+\mathbf{\Gamma }\left(t\right){\widehat{\mathit{M}}}_{\mathrm{0}}=\mathrm{0},\end{array}$$

and the system of first-order deformation equations yield the following:

$$\begin{array}{}\text{(63)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\widehat{\mathit{M}}}_{\mathrm{1}}\left(t\right)=-c_{\mathrm{0}}\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ (R_{\mathrm{2}}-R_{\mathrm{1}}){\widehat{M}}_{z\mathrm{0}}\left(t\right)+e^{R_{\mathrm{2}}t}{R}_{\mathrm{1}}{M}_{\mathrm{0}}\end{array}\right].\end{array}$$

The solution to this equation yields ${\widehat{M}}_{x\mathrm{1}}\left(t\right)={\widehat{M}}_{y\mathrm{1}}\left(t\right)=\mathrm{0}$, and ${\widehat{M}}_{z\mathrm{1}}\left(t\right)$ is given by the following:

$$\begin{array}{}\text{(64)}& \begin{array}{rl}& {\widehat{M}}_{z\mathrm{1}}\left(t\right)\\ & =-c_{\mathrm{0}}\left[(R_{\mathrm{2}}-R_{\mathrm{1}})\underset{\mathrm{0}}{\overset{t}{\int }}{\widehat{M}}_{z\mathrm{0}}\left(t^{\prime }\right)\mathrm{d}{t}^{\prime }+{\displaystyle \frac{R_{\mathrm{1}}}{R_{\mathrm{2}}}}(e^{R_{\mathrm{2}}t}-\mathrm{1})M_{\mathrm{0}}\right]\\ & =-c_{\mathrm{0}}\left[(R_{\mathrm{2}}-R_{\mathrm{1}})t<{\widehat{M}}_{z\mathrm{0}}\left(t\right)>+{\displaystyle \frac{R_{\mathrm{1}}}{R_{\mathrm{2}}}}(e^{R_{\mathrm{2}}t}-\mathrm{1})M_{\mathrm{0}}\right].\end{array}\end{array}$$

The first-order approximation of magnetization during the pulse is given by the following:

$$\begin{array}{}\text{(65)}& \mathit{M}\left(t\right)=e^{-R_{\mathrm{2}}t}\left[\begin{array}{c}{M}_{x\mathrm{0}}\left(t\right)\\ M_{y\mathrm{0}}\left(t\right)\\ M_{z\mathrm{0}}\left(t\right)+{\widehat{M}}_{z\mathrm{1}}\left(t\right)\end{array}\right].\end{array}$$

At this level of approximation, relaxation of transverse magnetization depends simply on R₂, while relaxation of M_z(t) depends on the average z magnetization during the pulse (calculated in the absence of relaxation). For macromolecules, R₂≫R₁ typically, and the term proportional to $R_{\mathrm{1}}/R_{\mathrm{2}}$ is small.

The nth-order deformation equation leads to the following expression for n>1:

$$\begin{array}{}\text{(66)}& {\widehat{\mathit{M}}}_n\left(t\right)=(\mathrm{1}+c_{\mathrm{0}}){\widehat{\mathit{M}}}_{(n-\mathrm{1})}\left(t\right)-c_{\mathrm{0}}\underset{\mathrm{0}}{\overset{t}{\int }}\mathbf{\Gamma }\left(t^{\prime }\right){\widehat{\mathit{M}}}_{(n-\mathrm{1})}\left(t^{\prime }\right)\mathrm{d}{t}^{\prime }.\end{array}$$

If $c_{\mathrm{0}}=-\mathrm{1}$, the above recursive expressions can be written compactly as follows:

$$\begin{array}{}\text{(67)}& \begin{array}{rl}{\widehat{\mathit{M}}}_n\left(t\right)=& (-\mathrm{1}{)}^{n-\mathrm{1}}\underset{\mathrm{0}}{\overset{t}{\int }}\mathbf{\Gamma }(t_{n-\mathrm{1}}\left)\mathrm{d}{t}_{n-\mathrm{1}}\underset{\mathrm{0}}{\overset{t_{n-\mathrm{1}}}{\int }}\mathbf{\Gamma }\right(t_{n-\mathrm{2}})\mathrm{d}{t}_{n-\mathrm{2}}\mathrm{\dots }\\ & \times \underset{\mathrm{0}}{\overset{t_{\mathrm{2}}}{\int }}\mathbf{\Gamma }\left(t_{\mathrm{1}}\right){\widehat{\mathit{M}}}_{\mathrm{1}}\left(t_{\mathrm{1}}\right)\mathrm{d}{t}_{\mathrm{1}}.\end{array}\end{array}$$

For a rectangular pulse applied to equilibrium magnetization (with magnitude set to unity for convenience), the initial approximations are as follows:

$$\begin{array}{}\text{(68)}& {\mathit{M}}_{\mathrm{0}}\left(t\right)=\left[\begin{array}{c}[\mathrm{1}-\cos ({\mathit{\omega }}_{\mathrm{e}}t\left)\right]\cos \mathit{\theta }\sin \mathit{\theta }\\ -\sin \left({\mathit{\omega }}_{\mathrm{e}}t\right)\sin \mathit{\theta }\\ \cos \left({\mathit{\omega }}_{\mathrm{e}}t\right){\sin }^{\mathrm{2}}\mathit{\theta }+{\cos }^{\mathrm{2}}\mathit{\theta }\end{array}\right],\end{array}$$

$$\begin{array}{}\text{(69)}& \begin{array}{rl}{\widehat{M}}_{z\mathrm{1}}\left(t\right)=-c_{\mathrm{0}}& (R_{\mathrm{2}}-R_{\mathrm{1}})\left[{\displaystyle \frac{\mathrm{1}}{{\mathit{\omega }}_{\mathrm{e}}}}\sin \left({\mathit{\omega }}_{\mathrm{e}}t\right){\sin }^{\mathrm{2}}\mathit{\theta }+t{\cos }^{\mathrm{2}}\mathit{\theta }\right]\\ & -c_{\mathrm{0}}{\displaystyle \frac{R_{\mathrm{1}}}{R_{\mathrm{2}}}}(e^{R_{\mathrm{2}}t}-\mathrm{1}).\end{array}\end{array}$$

In this case, Γ(t)=Γ and the series of approximations given in Eq. (67) can be summed to give the following:

$$\begin{array}{}\text{(70)}& \begin{array}{rl}\mathit{M}\left(t\right)=e^{-R_{\mathrm{2}}t}& \mathit{\left\{}{\mathit{M}}_{\mathrm{0}}\right(t)+\underset{\mathrm{0}}{\overset{t}{\int }}{e}^{\mathbf{\Gamma }(t-t^{\prime })}\\ & \times \left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ (R_{\mathrm{2}}-R_{\mathrm{1}})M_{z\mathrm{0}}\left(t^{\prime }\right)+e^{R_{\mathrm{2}}{t}^{\prime }}{R}_{\mathrm{1}}{M}_{\mathrm{0}}\end{array}\right]\mathrm{d}{t}^{\prime }\mathit{\}},\end{array}\end{array}$$

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 R₂, except for cases in which $<{\widehat{M}}_{z\mathrm{0}}\left(t\right)>=\mathrm{0}$, 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.

Figure 5 HAM approximations for rectangular (a, b, c) 90^∘, (d, e, f) 180^∘, (g, h, i) 270^∘, and (j, k, l) 360^∘ pulses applied to initial z magnetization. Values of (a, d, g, j) M_x(Ω), (b, e, h, k) M_y(Ω), and (c, f, i, l) M_z(Ω) are shown as functions of resonance offset Ω. Magnetization components in the absence of exchange, using Eq. (68) (black dotted line), first-order HAM approximation of the Bloch equations using Eq. (69) (reddish-purple dashed line), and exact HAM solution of the Bloch equations using Eq. (70) (blue solid line). Calculations used ${\mathit{\omega }}_{\mathrm{1}}/\left(\mathrm{2}\mathit{\pi }\right)=\mathrm{250}$ Hz, $R_{\mathrm{1}}=\mathrm{2}{\mathrm{s}}^{-\mathrm{1}}$, $R_{\mathrm{2}}=\mathrm{100}{\mathrm{s}}^{-\mathrm{1}}$, and $c_{\mathrm{0}}=-\mathrm{1}$.

3 Conclusion

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.