Analysis of PFG Anomalous Diffusion via Real-Space and Phase-Space Approaches

Pulsed-field gradient (PFG) diffusion experiments can be used to measure anomalous diffusion in many polymer or biological systems. However, it is still complicated to analyze PFG anomalous diffusion, particularly the finite gradient pulse width (FGPW) effect. In practical applications, the FGPW effect may not be neglected, such as in clinical diffusion magnetic resonance imaging (MRI). Here, two significantly different methods are proposed to analyze PFG anomalous diffusion: the effective phase-shift diffusion equation (EPSDE) method and a method based on observing the signal intensity at the origin. The EPSDE method describes the phase evolution in virtual phase space, while the method to observe the signal intensity at the origin describes the magnetization evolution in real space. However, these two approaches give the same general PFG signal attenuation including the FGPW effect, which can be numerically evaluated by a direct integration method. The direct integration method is fast and without overflow. It is a convenient numerical evaluation method for Mittag-Leffler function-type PFG signal attenuation. The methods here provide a clear view of spin evolution under a field gradient, and their results will help the analysis of PFG anomalous diffusion.


Introduction
Anomalous dynamic behavior exists in many polymer or biological systems [1][2][3][4][5][6].These systems often consist of various molecules such as macromolecules and small penetrant molecules like water.These molecules often demonstrate anomalous dynamics behavior: their rotational motion time constants usually obey a stretched exponential distribution, namely, the Kohlrausch-Williams-Watts (KWW) function distribution [7,8], and their translational diffusion jump time and jump length could follow power law distributions [2].Nuclear magnetic resonance (NMR) is one of the important techniques used to detect these anomalous dynamic behaviors.For instance, anomalous diffusion can be detected by pulsed-field gradient (PFG) NMR experiments.
PFG diffusion NMR [9][10][11][12][13][14][15] has been a powerful tool for measuring normal diffusion.The history of using a field gradient to measure diffusion can be tracked back to Hahn, who first observed the influence of the molecule diffusion under a gradient magnetic field upon echo amplitudes in 1950 [9].PFG diffusion measurements have broad applications in NMR and magnetic resonance imaging (MRI) [13][14][15][16].For instance, the water diffusion difference in different parts of the brain can be used as an important contrast factor to build imaging of acute strokes [16].The PFG technique, being an ultra-valuable tool for normal disunion, could play an increasingly important role in monitoring anomalous diffusion occurring in many polymer and biological systems.
PFG anomalous diffusion [17][18][19][20][21][22][23] is different from PFG normal diffusion.Unlike normal diffusion, anomalous diffusion has a non-Gaussian probability distribution [2,24], and its mean square displacement is not linearly proportional to diffusion time [2].These non-Gaussian characteristics make it complicated to analyze PFG anomalous diffusion.Although PFG anomalous diffusion can be approximately analyzed by traditional methods such as the apparent diffusion coefficient method, PFG theories based on the fractional derivative could not only improve the analysis accuracy on the diffusion domain size, diffusion constant, and other variables [17,18,[25][26][27][28], but also yield additional information such as the time derivative order α and space derivative order β that are related with diffusion jump time and jump length distributions determined by material properties.
Much effort has been devoted to studying PFG anomalous diffusion based on fractional calculus, which includes the propagator method [29], the modified Bloch equation method [17,25,30,31], the effective phase-shift diffusion equation (EPSDE)method [18], the instantaneous signal attenuation method [26], the modified-Gaussian or non-Gaussian distribution method [27], etc.Additionally, PFG anomalous diffusions in restricted geometries such as plate, sphere, and cylinder have been investigated [28].These theoretical methods analyze PFG anomalous diffusion from different angles.Therefore, each of them can have its own advantages in handling certain types of PFG anomalous diffusion.To better apply the PFG technique to studying anomalous diffusion, it is still valuable to develop new theoretical treatments for PFG anomalous diffusion.
In this paper, two methods based on the fractional derivative [1,2,6,[32][33][34] are proposed to give general analytical PFG signal attenuation expressions for anomalous diffusion.The first method is the recently developed EPSDE method [18].This method describes the spin phase evolution by an effective phase diffusion process in virtual phase space [18].Solving the EPSDE gives valuable information about the phase evolution process such as the phase probability distribution function and the moment of mean phase displacement.Meanwhile, other conventional methods render it difficult to get this phase information, and usually assume an approximate phase distribution such as Gaussian phase distribution [13].In this paper, it will be shown that a solvable PFG signal attenuation equation can be derived by applying a Fourier transform to the effective phase-shift equation.The second method is to observe the signal intensity at the origin, and is an ultra-simple new method.For a homogeneous diffusion spin system, although the magnetization amplitude attenuates because of the gradient magnetic field effect, the phase of magnetization keeps constant at the origin of the gradient field.Such a specific phase property is employed to derive a PFG signal attenuation equation in this paper.The above two methods give the same signal attenuation equation, from which the general PFG signal attenuation expression can be derived by the Adomian decomposition method [35][36][37][38][39].Besides the Adomian decomposition method, a direct integration method was proposed for the numerical evaluation of the PFG signal attenuation, which is a fast and simple method.The results include the finite gradient pulse width (FGPW) effect [13][14][15], namely, the signal attenuation during each gradient pulse application period.Theoretically, during a short gradient pulse, the PFG signal attenuation can be neglected; nevertheless, the gradient pulse used in a clinical MRI is usually long [16,40].Additionally, a longer gradient pulse allows the measuring of slower diffusion under the same gradient maximum intensity, which matters in the study of polymer and biological systems where the molecule diffusion is often slow.Therefore, we need to consider the FGPW effect in many real applications.The two methods here give the same results in terms of the FGPW effect.These results agree with reported results from some other methods [18,26,27] and the continuous time random walk simulation [26,29].Furthermore, PFG anomalous diffusion of intramolecular multiple quantum coherence (MQC) is also discussed [27].The MQC effect has the benefit of enhancing the gradient effect on PFG signal attenuation [41].The two methods provide complementary views of PFG anomalous diffusion from both the real space and phase space.

The Phase-Space Method: The Effective Phase-Shift Diffusion Method
For the sake of simplicity, only one-dimensional anomalous diffusion is considered in this paper.In PFG experiments, the field gradient pulse results in a time-and space-dependent magnetic field B(z, t) = B 0 + g(t) • z, where B 0 is the exterior magnetic field, z is the position, and g(t) is the time-dependent gradient [13][14][15].The magnetic field exerts a torque on each spin moment.The torque changes the spin angular momentum direction, which leads the spin to precess about the magnetic field with Larmor frequency ω(z, t) = −γB(z, t).In a rotating frame with angular frequency ω 0 = −γB 0 , a diffusing spin at position z(t ) has a time-dependent angular frequency γg(t ) • z(t ), and its phase accumulated along the diffusion path is [13][14][15]19] where φ(t) is the net-accumulating phase.The range of φ(t) is −∞ < φ(t) < ∞ rather than −π ≤ φ(t) ≤ π, and cos(φ(t)) is the projection factor of the spin magnetization to the observing coordinate axis.The PFG signal comes from the ensemble contribution from all spins by averaging over all possible phases [13][14][15]: where S(0) is the signal intensity at the beginning of the first gradient pulse, S(t) is the signal intensity at time t, and P(φ, t) is the accumulating phase probability distribution function.In a symmetric diffusion system, Equation ( 2) can be further written as Because the NMR signal comes from numerous spins, the percentage of spins with potentially infinite φ(t) is very small, and can therefore be neglected.Additionally, −1 ≤ cos(φ(t)) ≤ 1; therefore, 0 ≤ |S(t)| ≤ S(0), and S(t) is a finite quantity that can be measured in NMR experiments.A possibly infinite "mean square phase displacement" |φ(t)| β is not necessarily an obstacle in PFG NMR measurement, although it may have a significant effect in other academic fields.The PFG signal attenuation can be obtained based on Equation (2) or (3), but the phase evolution process based on the phase path integral Equation (1) is complicated.Nevertheless, the recently developed EPSDE method provides a simple way to describe the phase evolution based on an effective phase diffusion process [18].In the following, the effective phase diffusion equation method [18] will be briefly reintroduced.Additionally, a general solution including the FGPW effect will be given that has not been reported in [18].
For simplicity, only diffusion with a symmetric probability distribution is studied here.The self-diffusion process can be described by a random walk, which consists of a sequence of independent random jumps with waiting times ∆t 1 , ∆t 2 , ∆t 3 , . . ., ∆t n , and corresponding displacement lengths ∆z 1 , ∆z 2 , ∆z 3 , . . ., ∆z n .Based on the random walk, Equation (1) can be rewritten as [18] where is the wavenumber [13][14][15].In Equation (4), is negligible, and the temporal and spatial summation orders are interchanged because ∆t i g(t i ) and ∆z m are often noncorrelated.In most PFG experiments, K(t tot ) = 0 where t tot is the time at the end of the rephasing gradient pulse, so Equation ( 4) can be further simplified to [18] The −∑ m K(t m ) • ∆z m term in Equation ( 6) is an effective phase random walk with a jump length −K(t m ) • ∆z m and a jump waiting time ∆t m The phase random walk can be treated as an effective phase diffusion, and can be obtained from the corresponding spin diffusion in the real space by scaling the jump length ∆z(t) by a factor −K(t).Because |∆z| β / ∆t α ∝ D f and |K(t)∆z| β / ∆t α ∝ D φe f f (t), [18], where D f and D φe f f (t) are the fractional diffusion coefficients for the real-space diffusion and the effective phase diffusion, respectively, and their units are m β /s α and rad β /s α .As both the effective phase diffusion and the real spin particle diffusion belong to the same type of diffusion, they should obey the same type of diffusion equation.The one-dimensional real-space diffusion equation based on the fractional derivative is [6,[31][32][33][34] where M xy (z, t) is the spin magnetization in PFG experiments, and t D α * is the Caputo fractional derivative defined as [6,[31][32][33][34] and ∂ β ∂|z| β is the space fractional derivative [6,[31][32][33][34] defined in Appendix A. By replacing coordinate z with φ, and D f with D φe f f (t), respectively, from real-space diffusion Equation ( 9), the effective phase diffusion equations can be obtained [18] as As F{ f (x, t)} = f (q, t) and F ∂ β ∂|x| β f (x, t); q = −q β f (q) [31][32][33][34] where q is the wavenumber, applying a Fourier transform F{ f ( which is a q-space phase diffusion equation.Here, q is the wavenumber dedicated to the above Fourier transform, which differs from the field gradient pulse-induced wavenumber K(t) defined by Equation (5).When q = 1, P(q, t) in Equation ( 10) equals the attenuated signal amplitude because, in PFG experiments, the signal attenuation can be calculated based on Equation (2) as [18] The normalized value of S(0) equaling 1 will be used and dropped out throughout this paper.Substituting Equation (11) and q = 1 into Equation ( 10), we get which is the PFG signal attenuation equation.The solution of Equation ( 12) will be given in Section 2.3 as the same equation will be obtained by observing the signal intensity at the origin in the subsequent section.

Observing the Signal Intensity at the Origin in Real Space
For a nondiffusing spin in a rotating frame with angular frequency ω 0 = −γB 0 , its precession phase can be described as [13][14][15] where ψ(z, t) is the precession phase, which is different from the accumulated phase φ(t) described by (1).The nondiffusing spin system has a time-and space-modulated magnetization, which can be described as where S(t) is the magnetization amplitude, which can be referred to as signal intensity, and exp(−iK(t) • z) is the phase term.For a pulsed gradient spin echo (PGSE) or pulsed gradient stimulated echo (PGSTE) experiment with a constant gradient, as shown in Figure 1, the wavenumber is where δ is the gradient pulse length, and ∆ is the diffusion delay starting from the beginning of the first dephasing gradient pulse to the beginning of the last rephrasing gradient pulse.Because of the two counteracting gradient pulses-a dephasing pulse and a rephasing pulse-K(t) returns to 0 at the end of the rephasing pulse, and the phase of nondiffusing spins will be refocused.Therefore, the signal of the nondiffusing spin system does not attenuate when the T 2 relaxation is neglected.
For the spin diffusion in a homogeneous sample, its magnetization still can be described by Equation ( 14), M xy (z, t) = S(t) exp(−iK(t) • z).At a random position z, the possibilities of spins diffusing to opposite directions z ± ∆z are equal [26].The opposite movements yield exp(−iK(t) , and have no effect on the phase, exp(−iK(t) • z), but do make the signal intensity decay by a factor of cos(K(t) • ∆z) for these spins.By substituting M xy (z, t) = S(t) exp(−iK(t) • z) into diffusion Equation (7), and applying At the origin, z = 0, iγgz = 0, exp(−iK(t) • z) = 1, and S(t) exp(−iK(t) • z) = S(t); we thus have Equation ( 17) is identical to Equation ( 12) obtained previously by the EPSDE method in Section 2.1.
where δ is the gradient pulse length, and Δ is the diffusion delay starting from the beginning of the first dephasing gradient pulse to the beginning of the last rephrasing gradient pulse.Because of the two counteracting gradient pulses-a dephasing pulse and a rephasing pulse-( ) K t returns to 0 at the end of the rephasing pulse, and the phase of nondiffusing spins will be refocused.Therefore, the signal of the nondiffusing spin system does not attenuate when the T2 relaxation is neglected.
(a) (b) This method of observing the signal intensity at the origin could be understood in another way: if only the signal intensity or amplitude is observed from a selected position at different diffusion times t i , i = 1, 2, 3 • • • in a PFG anomalous diffusion experiment, the observed signal intensity or amplitude shall still obey Equation ( 17) regardless of whether the phase is observed or not.It is reasonable that both ways can give the same signal attenuation expression, because the signal amplitude is homogeneous throughout the sample at each instant, and selecting the origin at a random position in the sample does not affect the result.

Analytical Solution by the Adomian Decomposition Method
The same equation as Equation ( 12) or (17) has been obtained by the fractional integral modified Bloch equation method.The three different methods get the same PFG signal attenuation equation.Reference [30] employed the Adomian decomposition method to give a general analytical solution to Equation (12).In the following, the general solution from the Adomian decomposition will be given in more detail for the cases of normal diffusion and general anomalous diffusion.

Normal Diffusion
First, the signal attenuation equation-Equation ( 12) or ( 17)-can be used to obtain the PFG signal attenuation expression for normal diffusion, which is a specific case of anomalous diffusion with α = 1 and β = 2.When α = 1 and β = 2, Equation (17) reduces to where D is the diffusion constant of normal diffusion.The solution of Equation ( 18) is which is the PFG signal attenuation expression for normal diffusion.The same result as expression (19) was obtained in [42].For PGSE or PGSTE experiments, the PFG signal attenuation calculated based on Equation ( 19) is exp −D f γ 2 g 2 δ 2 (∆ − δ/3) , which is a routinely used expression for PFG normal diffusion [13][14][15].

General Anomalous Diffusion
Second, the PFG signal attenuation for general anomalous diffusion can be obtained.The same equation as Equation ( 17) and its solution have been obtained by the modified Bloch equation method in [30].According to [30], Equations ( 12) or ( 17) can be written equivalently as where Based on the Adomian decomposition method [35][36][37][38][39], the solution of Equation ( 20) is [30,31] where and In PGSE or PGSTE experiments, the time t can be separated into three periods: 0 < t ≤ δ, δ < t ≤ ∆, and ∆ < t ≤ ∆ + δ.If the initial condition, S (1) (0 + ) = 0, is used [31-34] for free diffusion in a homogeneous sample, we can get the following: (a) PFG signal attenuation under short gradient pulse (SGP) approximation: δ is short enough and the diffusion inside each gradient pulse can be neglected.We get where S 0 (∆) = 1, and . Equation ( 23) replicates the SGP approximation result obtained in references [18,26,27].(b) Single pulse attenuation: This is an ideal situation-the first gradient pulse is regular, but the second gradient pulse is infinitely narrow and has the purpose of counteracting the first gradient pulse.We get where , and ) is a Mittag-Leffler-type function [43].Equation ( 24) is consistent with the results obtained by the modified Bloch equation proposed in [25].(c) General PFG signal attenuation: The PGSE or PGSTE experiment includes three periods: 0 < t 1 ≤ δ, δ < t 2 ≤ ∆, and ∆ < t 3 ≤ ∆ + δ.The integration in Equation (22c) during these three periods is tedious, but can be calculated with computer assistance.Nevertheless, we can get the first and second terms of Equation (22a) as the following: where B(x,y) and B(a;x,y) are the Beta function and incomplete Beta function, respectively.When t = ∆ + δ, Equation (25b) gives Equation ( 26) agrees with the signal attenuation expression S(t) = E α,1 − t 0 D f K β (t )dt α obtained by the instantaneous signal attenuation method [26].At small signal attenuation, be approximately expanded as Equation ( 27) is the same as that given by the combination of Equations (25a) and ( 26) (note that B(x, y) = B(y, x) and B(a; x, y) = B(1 − a; y, x)); however, at large signal attenuation, E α,1 − t 0 D f K β (t )dt α may deviate from the combination of Equations (25a) and ( 26) as the higher-order terms in the expansion cannot be omitted.This agreement can be explained by the following: in typical PGSE or PGSTE experiments, both the gradient intensity and pulse length of the dephasing pulse and the rephrasing pulse are identical, When 0 < α ≤ 2, β = 2, the general anomalous diffusion reduces to time-fractional diffusion and Equation (25b) can be further calculated as When α = 1, 0 < β ≤ 2, the general anomalous diffusion reduces to space-fractional diffusion.Equation ( 22) reduces to S(t) = exp − t 0 D f K β (t )dt , which is the same as that obtained by [18,26,30,31].

Numerical Evaluation of Mittag-Leffler Function-Based PFG Signal Attenuation by the Direct Integration Method
Although the Adomian decomposition method can be used to numerically evaluate the PFG signal attenuation, its calculation speed is slow and it could cause overflow at large signal attenuation.References [30,31,44] proposed an alternative method: a direct integration method that can give the same numerical results, but with a much faster speed and without overflow.From Equation ( 20), if we where , Equation ( 30) can be further written as [44] Equation ( 30) can be rewritten in a discrete form as [44] where t j = j ∑ k=1 ∆t k .Based on Equation (32), the PFG signal attenuations S(t 1 ), S(t 2 ), . . ., S(t n ) can be calculated step by step starting from S(t 1 ) through to S(t n ).Similarly, the Mittag-Leffler-type function and its derivative can be numerically evaluated by [44] This method is simple, and there is no overflow issue in the calculation.

Continuous-Time Random Walk (CTRW) in a Lattice Model
Random walk simulation is a powerful numerical method that employs a stochastic jump process to model normal and anomalous diffusion in physics, chemistry, biology, and many other disciplines.It can be used to simulate the PFG signal attenuation in mathematically tractable or intractable systems.In this paper, the continuous-time random walk (CTRW) simulation method developed in [26] is used to verify the theoretical results; it is based on two models: the CTRW model [45] and the Lattice model [46,47].In the simulation, a sequence of independent random waiting time and jump length combinations (∆t 1 , ∆ξ 1 ), (∆t 2 , ∆ξ 2 ), (∆t 3 , ∆ξ 3 ), . . ., (∆t n , ∆ξ n ) is produced by the computer program.The individual waiting time ∆t and jump length ∆ξ are given according to [26] by the following expressions: and where η t and η z are scale constants, Φ = π(V − 1/2), and U, V ∈ (0, 1) are two independent random numbers.The CTRW simulation based on the above waiting time and jump length satisfies the time-space fractional diffusion equation under the diffusive limit, providing a simple way to simulate anomalous diffusion in various research areas, such as physics and economics [45].
Although continuous waiting time and jump length are used in the simulation, the accumulating spin phase associated with the diffusion path is recorded in a discrete manner.Such a discrete phase recording manner is convenient and reasonable.The simulation can be viewed as a numerical PFG experiment, whose observables such as phase can be observed in a discrete time selected by an experiment observer; the observing or recording manner will not affect the fundamental numerical experiment process, namely the production of the continuous random walk sequence (∆t 1 , ∆ξ 1 ), (∆t 2 , ∆ξ 2 ), (∆t 3 , ∆ξ 3 ) , . . ., (∆t n , ∆ξ n ).The spin phase was recorded by the lattice model developed in [46,47], which has been applied to simulate PFG diffusion in polymer systems [48].The spin phase in the simulation is where ∆ξ j , and t j is the discrete record time, which takes place between the lth and (l + 1)th steps of the random walk.The second term γg(t j )ξ(t l+1 )(t j − t l ) on the right-hand side of Equation ( 37) corresponds to the partial phase evolution of the (l + 1)th jump step that needs to be recorded at the time t j .The PFG signal attenuation in the simulation can be obtained by averaging over all the walkers in the simulation [26,47,48]: where N walks is the total number of walks.The total number of walks used in each simulation is 1,000,000.

Results and Discussion
This paper uses two different methods to describe PFG anomalous diffusion: the EPSDE method [18] and the method of observing the signal intensity at the origin.The major results in this paper are summarized in Table 1.The EPSDE method [18] uses an effective phase diffusion process to describe the phase evolution, while observing the signal intensity at the origin only monitors the real-space signal amplitude change at the origin where the phase of the magnetization remains constant.Each of these two approaches has its unique advantages.The EPSDE method can directly obtain the phase distribution, which greatly simplifies PFG diffusion analysis [18]; in contrast, the traditional methods often need to get the real-space spin particle probability distribution first, which is then used to obtain the phase distribution.The ultra-simple method of observing the signal intensity at the origin obtains the signal attenuation equation by substituting the magnetization M xy (z, t) = S(t) exp(−iK(t) • z) directly into the diffusion equation and using the phase property exp(−iK(t) • z) = 1 and iγgz = 0 at the origin.These two approaches view the spin evolution process from two different spaces-the real space and the phase space-providing a clear picture of spin dynamics in PFG diffusion experiments.
Table 1.Comparison of pulsed-field gradient (PFG) anomalous diffusion results by the effective phase-shift diffusion method, the observing the signal intensity at the origin method and other methods.

Virtual Phase-Space Diffusion Real-Space Spin Diffusion
Effective phase-shift diffusion equation [18] Observing the signal intensity at the origin Diffusion equations: [32]  Solving methods: Fourier transform Substitute M xy (z, t) = S(t) exp(−iK(t) • z) into equation PFG signal attenuation equation*: attenuation expression by the Adomian decomposition method [30,31]*: Under short gradient pulse (SGP) approximation: At small attenuation: Other methods: (1) (3) [18,26,27] * The PFG signal attenuation equation and expression obtained in this paper are the same as that obtained by the integral modified Bloch equation method [30].
These two methods agree with each other.They get the same PFG signal attenuation Equation ( 11) or (16), which can be solved by the Adomian decomposition method.The solution gives PFG signal attenuation expressions (21a)-(21c), which include the FGPW effect.Understanding the FGPW effect is important as clinical MRI applications often use a long gradient pulse.Additionally, using long gradient pulses allows researchers to monitor slower diffusion in polymer or biological systems where the molecules or ions often diffuse slowly.The consistent results from two different methods help us better understand the FGPW effect.
Additionally, the two methods agree with other methods.First, the PFG signal attenuation Equation ( 17) obtained here is the same as that obtained by the integral modified Bloch equation method [30,31].Second, the results agree with the CTRW simulation as shown in Figure 2. In Figure 2a, the fractional diffusion constant D f = 6.6 × 10 −11 m β /s α was obtained from fitting the curve of mean square displacement versus t by equation |z| β = 2D f t α /Γ(1 + α). Figure 2b shows the PFG signal attenuation of anomalous diffusion obtained by Equations ( 22a)-(22c) based on the Adomian decomposition method, Equation (32) based on the direct integration method, and the CTRW simulation.Other parameters used in Figure 2 are α = 0.5, β = 2, g = 0.1 T/m, ∆ − δ equaling 0 ms, 20 ms.Interested readers can refer to [30] for additional comparisons between the theoretical prediction and the CTRW simulation.Third, the PFG signal attenuation expressions (22a)-(22c) are also consistent with those obtained by various approximation methods such as the non-Gaussian phase distribution (nGPD) method [27] and the instantaneous signal attenuation method [26].The same SGP approximation results, Equation (23), can be obtained by all these five different methods: the EPSDE method [18], the instantaneous signal attenuation method [26], the nGPD method [27], the method to observe the signal intensity at the origin, and the modified Bloch equation method [30].Moreover, Equation (27), S(t) ≈ E α,1 − t 0 D f K β (t )dt α , can be approximately obtained in this paper; this agrees with the result obtained from the instantaneous signal attenuation method [26].
The numerical evaluation of PFG signal attenuation by Equation ( 20) or ( 30) can be conveniently performed by the direct integration method.Figure 2b shows that the results of the direct integration method agree with those obtained from the Adomian decomposition method.Compared to the Adomian decomposition method, the computing speed of the direct integration method improves by orders of magnitudes.The speed improvement is because in the direct integration method, each S(t j ) in Equation (32) is only a single time integration, while in the Adomian decomposition method, S(t j ) = ∞ ∑ n=0 S n (t j ) is a superposition of many terms of time integration.Additionally, from Equations ( 33) and (34), the direct integration method can be used to calculate Mittag-Leffler-type functions and their derivatives.Figure 3 shows that the MLF calculated from the direct integration method agrees with that calculated by other methods in [49,50].The Pade approximation method in [50] is only for subdiffusion.Because the direct integration method does not cause overflow, it can be a useful method for calculating Mittag-Leffler-type functions.The Fortran code for MLF calculation by the direct integration method can be obtained from the following link: https://github.com/GLin2017/Mittag-Leffler-function-calculated-by-Direct-Integration.Matlab code for calculating the MLF function by other methods to a desired accuracy can be found at www.mathworks.com/matlabcentral/fileexchange/8738-mittag-leffer-function,2005.This direct integration method will be a great help to the application of current theoretical results to PFG anomalous diffusion in NMR and MRI.
In the current results, only PFG anomalous diffusion of single quantum coherence is considered.The results here can be easily extended to handle PFG anomalous diffusion of intramolecular MQC [27,40].As the gradient field-induced phase evolution of an n-order intramolecular MQC is n times faster than the corresponding single quantum coherence, the PFG signal attenuation of intramolecular MQC is the same as that of single quantum coherence with an effective gradient intensity ng [42].Therefore, for intramolecular MQC in PGSE or PGSTE experiments, the PFG signal attenuation equation is where From Equation (40), it is obvious that the a MQC (t) in the MQC equals n β a(t).
Mathematics 2018, 6, x 13 of 16 for calculating the MLF function by other methods to a desired accuracy can be found at www.mathworks.com/matlabcentral/fileexchange/8738-mittag-leffer-function,2005.This direct integration method will be a great help to the application of current theoretical results to PFG anomalous diffusion in NMR and MRI.33) from the direct integration method, the Pade approximation method [50], and the method used in [49].
In the current results, only PFG anomalous diffusion of single quantum coherence is considered.The results here can be easily extended to handle PFG anomalous diffusion of intramolecular MQC [27,40].As the gradient field-induced phase evolution of an n-order intramolecular MQC is n times faster than the corresponding single quantum coherence, the PFG signal attenuation of intramolecular MQC is the same as that of single quantum coherence with an effective gradient intensity ng [42].Therefore, for intramolecular MQC in PGSE or PGSTE experiments, the PFG signal attenuation equation is for calculating the MLF function by other methods to a desired accuracy can be found at www.mathworks.com/matlabcentral/fileexchange/8738-mittag-leffer-function,2005.This direct integration method will be a great help to the application of current theoretical results to PFG anomalous diffusion in NMR and MRI.33) from the direct integration method, the Pade approximation method [50], and the method used in [49].
In the current results, only PFG anomalous diffusion of single quantum coherence is considered.The results here can be easily extended to handle PFG anomalous diffusion of intramolecular MQC [27,40].As the gradient field-induced phase evolution of an n-order intramolecular MQC is n times faster than the corresponding single quantum coherence, the PFG signal attenuation of intramolecular MQC is the same as that of single quantum coherence with an effective gradient intensity ng [42].Therefore, for intramolecular MQC in PGSE or PGSTE experiments, the PFG signal attenuation equation is from the direct integration method, the Pade approximation method [50], and the method used in [49].
This paper only considers the spin self-diffusion that can be described by the time-space fractional diffusion equation based on the fractional derivative.In general, the two methods used in this paper can be applied to other types of anomalous diffusions such as that described by the time-space diffusion equation based on the fractal derivative [51], etc.Additionally, only the symmetric anomalous diffusion in homogeneous spin systems is studied here.In real applications, the anomalous diffusion can take place in complicated systems such as inhomogeneous systems, restricted geometries [28], anisotropic systems [44], nonsymmetric systems, etc.Additionally, the gradient field may be nonlinear [52].The current methods may face challenges in these complicated situations, which reminds us that much effort is required for the study of PFG anomalous diffusion.

Figure 3 .
Figure 3. Comparing the Mittag-Leffler function (MLF) calculated by different methods: Equation (33) from the direct integration method, the Pade approximation method[50], and the method used in[49].

Figure 3 .
Figure 3. Comparing the Mittag-Leffler function (MLF) calculated by different methods: Equation (33) from the direct integration method, the Pade approximation method[50], and the method used in[49].