Approximate representations of shaped pulses using the homotopy analysis method

Abstract 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.


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;Kupce et al., 1995;Cavanagh et al., 2007).The Euler angles for an arbitrary shaped pulse can be extracted from a numer-ical 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): U = e −iγ (τ p )I z e −iβ(τ p )I x e −iα(τ p )I z = e −iχ + (τ p ) cos( 2 ) e iχ + (τ p ) cos( in which χ ± (τ p ) is given as follows: I k are the Cartesian spin operators, and the product is timeordered from right to left.The propagator for the kth rectangular pulse segment is given as follows: Published by Copernicus Publications on behalf of the Groupement AMPERE.
T. Crawley and A. G. Palmer III: HAM for shaped pulses in which κ(τ p ) and η(τ p ) are given by the following: κ(τ p ) = cos(ω ek t k /2) + i cos θ k sin(ω ek t k /2) η(τ p ) = e iφ k sin θ k sin(ω ek t k /2). (4) 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); ω ek = (ω 2 1k + 2 k ) 1/2 is the effective field; and θ k = tan −1 (ω 1k / 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): in which the following applies: where ω ± (t) = ω x (t) ± iω y (t) 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.

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 ∈ [0, 1].
This relationship is established by constructing the homotopy as follows (Liao, 2012): in which L[ ] is a linear (differential) operator, and N [ ] is an (nonlinear differential) operator satisfying the following: where y 0 (t) is an initial approximation for the desired solution y(t), c 0 = 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 H [ (t; 0) : 0] = 0, Eq. ( 13) requires (t, 0) = y 0 (t).Similarly, when q = 1, the homotopy becomes as follows: Therefore, when H [ (t; 1) : 1] = 0, Eq. ( 14) requires (t, 1) = y(t).Stated more succinctly, as q increases from 0 → 1, (t; q) deforms from the initial approximation y 0 (t) to the exact solution y(t).To proceed, the Maclaurin series for (t; q) is assumed to exist as follows: in which y n (t) is given by the following: Conditions concerning convergence of the series are discussed by Liao (2012).Equation ( 17) has the desired property (t; 0) = y 0 (t) and yields the following: HAM then consists of successively determining y n (t), beginning with the initial approximation y 0 (t), until y(t) is approximated to the desired accuracy.The choices of L[ ], y 0 (t), c 0 , and H (t) provide considerable flexibility in finding approximate solutions to differential equations.For simplicity in the following, the auxiliary function is 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 y 0 (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 (t; 0) = y 0 (t) and Eq. ( 18).Substituting for N [ ], L[ ], and y 0 (t) yields the following: in which the final line is obtained using dy 0 (t)/dt = 0.This differential equation does not contain a term proportional to y 1 (t).Hence, the homogenous equation (setting the righthand 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 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): 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 Results obtained using y 0 (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 L[ ] and y 0 (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 y 0 (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 righthand side to 0) are y ± (t) = e ±iδ(t)/2 .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 y n (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:

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 0 (t) = 1 in method 1 leads to obtaining simple iterative solutions that can be calculated very efficiently.The form of y 0 (t) given in Eq. ( 36) could also be used in Eq. ( 24) to obtain an alternative expression for y 1 (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 0 (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 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 0 = −1 are nearly indistinguishable from the exact result of Eq. (3) (using τ p = τ k ) over the range of resonance offsets from 0 to /ω 1 = 15 1/2 .The first-order result provides a highly accurate estimate of γ (τ p ) but overestimates β(τ p ).The role of the convergence control parameter c 0 is illustrated in Fig. 1c, d.A value of c 0 = −0.925was chosen, using Eqs.( 46) and ( 30), to scale the first-order result for β(τ p ) to be equal to π/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 0 works well for a rectangular 180 • pulse (not shown) and 90 • EBURP-2, 90 • Q5, and 180 • RE-BURP 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 0 = −1, even to the third order in , is accurate for β(τ p ) only to slightly more than /ω 1 = 1.When c 0 = −0.925, the third-order power series has improved accuracy for resonance offsets up to nearly /ω 1 = 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 /ω 1 = 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.
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 (Kupce 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 0 = −1, the fifthorder approximation is extremely accurate compared with numerical calculations using Eqs.( 1)-(3) (Fig. 2a-c).With c 0 = −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 powerseries 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 c 0 = −1 and c 0 = −0.925.The choice of c 0 = −0.925has 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 (Kupce and Freeman, 1995) pulses.As for the EBURP-2 pulse, the fifth-order approximation for the REBURP pulse is highly accurate for both c 0 = −1 and c 0 = −0.925.The third-order approximation also is highly accurate when c 0 = −0.925.The WURST-20 pulse uses a linear frequency shift, generated by applying a quadratic phase shift during the pulse, and is an examhttps://doi.org/10.5194/mr-2-175-2021 Magn.Reson., 2, 175-186, 2021 ple of a phase-modulated or complex waveform.Again, the more complicated waveform requires higher order approximation, but 11th-order, with c 0 = −1, or ninth-order, with c 0 = −0.925,results are highly accurate.Method 2 yields a power series for y(t).If c 0 = −1, the resulting series is identical to the power series expansion obtained from perturbation theory (Li et al., 2014), while c 0 = −0.925provides 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 0 = −1, 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 τ α /τ p and τ γ /τ 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 directly to the original Riccati equation of Eq. (5).For example, by analogy to method 2, choosing the following: in which f 0 (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:  df 3 (t) dt 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 ra-diofrequency 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 M k (t) = e −R 2 t Mk (t) are the Cartesian components of the magnetization, and M 0 is the equilibrium magnetization.The use of transformed variables M(t) rather than M(t) simplifies the following discussion.The linear operator is chosen as L[g(t)] = dg(t)/dt, in which g(t) = [g x (t), g y (t), g z (t)] T is an arbitrary vector function.The nonlinear operator is as follows: The initial zeroth-order approximations for HAM are M0 (t) = M 0 (t) and are given by the solutions to the Bloch equations in the absence of exchange for initial equilibrium magnetization [0, 0, M 0 ] T .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 Mx1 (t) = My1 (t) = 0, and Mz1 (t) 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 R 2 , while relaxation of M z (t) depends on the average z magnetization during the pulse (calculated in the absence of relaxation).For macromolecules, R 2 R 1 typically, and the term proportional to R 1 /R 2 is small.The nth-order deformation equation leads to the following expression for n > 1: If c 0 = −1, 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 R 2 , except for cases in which < Mz0 (t) >= 0, such as the onresonance 360 • pulse.The above expressions display the fundamental dependence of relaxation during a pulse applied to equilibrium magnetization on the time-average z magnetization.

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.Special issue statement.This article is part of the special issue "Robert Kaptein Festschrift".It is not associated with a conference.
Acknowledgements.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.
Financial support.This research has been supported by the National Institutes of Health (NIH; grant nos.R35GM130398 and P41GM118302).
Review statement.This paper was edited by Malcolm Levitt and reviewed by Fabien Ferrage and one anonymous referee.