Inverse Laplace Transformation Analysis of Stretched Exponential Relaxation

We investigate the effectiveness of the Inverse Laplace Transform (ILT) analysis method to extract the distribution of relaxation rates from nuclear magnetic resonance data with stretched exponential relaxation. Stretched-relaxation is a hallmark of a distribution of relaxation rates, and an analytical expression exists for this distribution for the case of a spin-1/2 nucleus. We compare this theoretical distribution with those extracted via the ILT method for several values of the stretching exponent and at different levels of experimental noise. The ILT accurately captures the distributions for $\beta \lesssim 0.7$, and for signal to noise ratios greater than $\sim 40$; however the ILT distributions tend to introduce artificial oscillatory components. We further use the ILT approach to analyze stretched relaxation for spin $I>1/2$ and find that the distributions are accurately captured by the theoretical expression for $I=1/2$. Our results provide a solid foundation to interpret distributions of relaxation rates for general spin $I$ in terms of stretched exponential fits.


Introduction
One of the most important quantities measured in magnetic resonance is the spin lattice relaxation rate, T −1 1 . This quantity probes the interaction between the nuclear spins and their environment, and reveals information about the local dynamics at the nuclear spin site. In conductors, the dominant contribution to T −1 1 arises from the hyperfine coupling to the electron spins, and in this case T −1 1 can be directly related to the dynamical spin susceptibility of the electrons [1]. This relationship has been extensively utilized to study a number of correlated electron systems, ranging from high temperature superconducting cuprates [2][3][4][5][6] , heavy fermion materials [7][8][9], iron-based superconductors [10][11][12] and other exotic materials [8,13]. A common issue in correlated electron systems is the presence of electronic inhomogeneity [14][15][16][17][18][19]. Even in a single crystal, inhomogeneity may arise intrinsically via frustration among competing orders, from the presence of impurities, or via inhomogeneous electronic responses such as in the mixed state of a type II superconductor. If the inhomogeneity is purely static, then the nature of the inhomogeneous distribution can often be studied via the effect on the NMR spectrum [20,21]. If the inhomogeneity fluctuates, then the spectrum may be motionally narrowed, precluding such investigations. However, in many cases the dynamics of the inhomogeneity may be reflected in T −1 1 , in which case there is a distribution of relaxation rates, rather than a single homogeneous T −1 1 . If each nucleus in a crystal relaxes with a different relaxation rate, W 1 (r), where r describes the position of the nucleus in the lattice, then the total magnetization measured experimentally exhibits a complicated relaxation curve. If the distribution of relaxation rates, P (W 1 ), is sufficiently narrow, then the NMR magnetization recovery data can be fit to a exponential form, exp[−t/T 1 ] (for a spin I = 1/2 nucleus), where t is the recovery time, and T −1 1 is the median of P (W 1 ). On the other hand, if P (W 1 ) is wide, then this exponential form will not fit the magnetization recovery data well, making it difficult to extract a meaningful value of T −1 1 . It is common practice to fit the recovery to a stretched exponential form, exp[−(t/T 1 ) β ], where β is the so-called stretching exponent that provides a rough measures of the width of P (W 1 ), and satisfies 0 < β ≤ 1 [22,23]. In many glassy systems, it has been observed that β decreases from unity as the temperature is reduced [24,25]. As β is reduced, the distribution P (W 1 ) grows by several decades, and β is related to the logarithmic width of the distribution. Analyzing magnetization recovery data with a stretched exponential is straightforward to implement and requires only one extra fitting parameter.
A significant disadvantage to this approach, however, is that it makes an implicit assumption about the shape of P (W 1 ) that may or may not accurately reflect the true distribution. A more direct method to extract the distribu-tion is desirable, but is in fact an ill-posed problem and is non-trivial to implement. The magnetization decay curve, M (t), is related to P (W 1 ) via a Fredholm integral of the first kind [26,27]. For a spin-1/2 nucleus, P (W 1 ) reduces to the inverse Laplace transform of M (t). Such problems are notoriously difficult because even small experimental errors in the values of M (t) give rise to large variations in P (W 1 ), and there is often no unique solution for a given data set. Different approaches have been developed to extract P (W 1 ), such as the maximum entropy method [28], and via linearization methods such as Tikhonov regularization [29]. The latter approach was adopted early on by researchers studying pore size distribution of rocks in the petrochemical industry [30][31][32][33], for the investigation of dielectric spectra in glasses [34], and recently has been used to analyze the glassy NMR behavior of high temperature superconductors [35,36]. This technique holds promise to shed light on many physical systems of interest, but several questions concerning the limits of validity of this approach remain outstanding. To better understand these limits, we have conducted numerical studies comparing the inverse Laplace transform (ILT) for stretched exponential decays for several different nuclear spins (I = 1/2, · · · , 9/2), different levels of signal to noise ratios, and different numbers of measured time points. We find that the ILT algorithm reproduces the theoretical distribution for a spin 1/2 nucleus for small stretching exponents β ≤ 0.8 when the distribution is not narrowly peaked. For higher spin nuclei, we find that P (W 1 ) for stretched relaxation is independent of I as long as the stretched relaxation curve is properly defined. These results provide important guidance for setting up experiments with sufficient signal to noise to properly extract the distribution of relaxation rates, and for interpreting the distribution when the relaxation can be described by stretched exponentials.

Methods
When a spin I = 1/2 nucleus at lattice position r is not in thermal equilibrium, the magnetization component along the quantization axis (typically the magnetic field direction) relaxes as: where m 0 is the equilibrium magnetization, and φ is a parameter that describes the initial condition and depends on the pulse sequence employed in the measurement. We assume that W 1 (r) depends on position, r. All of the nuclei contribute to the measured signal, so that where V is the volume of the sample, and M 0 = N 0 m 0 , where N 0 is the number of nuclei in the crystal. In the second line, rather than integrating over real space we express the integral as a distribution over a normalized distribution of W 1 and kernel function K(W 1 , t) = 1 − φe −W1t . This kernel function changes for higher spins, I > 1/2, as described below, however the general approach to solving for P (W 1 ) remains the same. For I = 1/2, is equivalent to the Laplace transform of P (W 1 ). Thus in principle, the distribution can be obtained by simply taking the inverse Laplace transform of the measured data.
The ILT approach offers a powerful method to determine P (W 1 ), however it requires a number of assumptions. To determine the distribution, the problem is first linearized: where i ∈ {1, · · · , N } are the measured time points, K ij = K(W 1,j , t i ), e i are experimental errors, P j = P (W 1j ), and j ∈ {1, · · · , L} with L > N are the points in the distribution. Since L > N , there are in fact more points in the distribution than experimentally measured, and the vector P is underdetermined. Tikhonov regularization [29] offers a method to obtain a solution by minimizing the functional: subject to the condition that every element P j ≥ 0. Here α is the Tikhonov regularization parameter that enforces P to have a stable solution. This procedure ensures that the distribution is positive definite, hence physically realistic, but has the effect of broadening and smoothing the distribution, depending on the choice of α [27,37]. The solution of (5) is P =H ·K † · c, where † means transpose, and the matrixH has elements where I is the identity matrix. This equation can be solved iteratively [36]. The sum of the residuals is: The distribution clearly depends on the choice of α, and becomes broader and smoother as α increases. The optimal value of α is usually determined by the so-called self-consistency method [27,36], in which α is chosen as the minimum of either α 1 or α 2 , where α 1 = | e|, the sum of the experimental errors of the measurements of M , and α 2 satisfies: We use the IGOR Pro software environment to solve for c numerically using a set of M data, for various numbers of N data points, and with M = 128 logarithmically-spaced values of W 1 . By computing χ(α) for a broad range of α, we find the optimal regularization parameter and use this to determine the distribution P for a given data set,

Results
To determine the effectiveness of the ILT method, it is valuable to test the algorithm to extract known distributions from test data sets. It is also instructive to determine the optimal experimental conditions to get the most accurate measurement of P (W 1 ). For example, the number of data points in a typical experiment lies between N ∼ 5 − 20. However, the major constraint is the total experimental time, t expt = N t 1 , where t 1 is the measurement time for a single point. For Gaussian noise, the measurement error e i ∼ t −1/2 1 ∼ N 1/2 , for a fixed t expt , therefore fewer points would result in lower measurement noise. An interesting question is whether it is better to have more points, N , with higher noise, or fewer points with lower noise, in order to determine P (W 1 ) with the best fidelity.

Stretched Relaxation of a Spin 1/2
We first consider the case of stretched exponential relaxation of a spin I = 1/2, with the kernel function: with the corresponding magnetization recovery vector: where W * 1 is a characteristic rate scale. This function is shown in Fig. 1(a) on a linear-log scale with W * 1 = 1. The distribution P β (W 1 ) can be expressed analytically as an infinite series: where Γ(x) is the Gamma function [23]. Figure 1 shows the distributions corresponding to stretched relaxation for several values of β. Note that because the W 1 values are distributed on a logarithmic scale, it is necessary to multiply P β (W 1 ) by W 1 to properly normalize the distribution. The distributions are centered close to W * 1 , which is approximately equal to the median of the distribution. The distribution approaches a delta function as β → 1. As β is reduced, the distribution broadens considerably, and is several decades in width once β ≤ 0.5. Fig. 2 shows the distributions extracted using the ILT method for several values of β with zero noise (e i = 0) and with N = 15 time points. There is relatively good agreement for intermediate values of β, but once β 0.8, the distribution narrows and the ILT method fails to capture the narrow width. To measure the effectiveness of the approach, we compute the sum of the squares of the residuals, S 2 , defined as: where P is determined by ILT. As shown in Fig. 3, S 2 generally increases as β approaches unity. However, S 2 has inflated values at β < 0.5 due to domain constraints. For lower β recoveries, the complete magnetization recovery is not captured within the given time domain. As a result, artificial higher relaxation rates are produced in the ILT distributions of lower β recoveries, resulting in larger S 2 values.

Optimal Number of Measurements
To understand how well the algorithm behaves with different numbers of measured time points, we compare the extracted distribution for different values of N . As shown in Fig. 4 for β = 0.8, including a greater number of measured recovery points improves the quality of the 'fit' such that the ILT distribution more accurately reproduces the exact solution. In each case, the mean of the distribution is correct, but the width is too wide for N = 6 points. For 15 points, the agreement is better, but there is an oscillation present in the upper tail of the ILT distribution that is not present in the exact solution. The behavior of these oscillations depends on the N when N 12−15, but there are no obvious trends. In fact, for N = 30 the oscillations appear somewhat larger than for N = 15. To quantify the difference between the ILT and the exact distributions, we compute S 2 for various values of β and N as shown in Fig.  5. This quantity appears to reach an asymptotic value by approximately N = 12 to 15. For smaller values of N , S 2 oscillates between larger and smaller values for N odd or even values, respectively. The origin of this behavior is not understood.

Sensitivity to Noise
In order to understand the effect of experimental noise, we added random values e i to each M i value, where each e i is sampled from a Gaussian distribution centered at zero with second moment σ 2 n , such that the signal to noise ratio SN R = σ −1 n . A sample set of magnetization recovery points (N = 15) with β = 0.8, and the corresponding ILT distributions are shown in Fig. 6. Introducing noise clearly affects the ability of the ILT algorithm to accurately reproduce the known distribution. As shown in Fig. 7, there is an approximate power law relationship: S 2 ∼ SN R −s , where s increases as β decreases. Extra structures, such as spurious peaks and shoulders are apparent in the ILT distributions. These artifacts persist even up to SNR lev-

Relaxation of Higher Spins
For spins greater than I = 1/2, the kernel function consists of multiple exponential decays reflecting the normal modes relaxation of the spin system. If the nuclei experience only a Zeeman interaction, the energy splittings between the states are all equal. However, it is common that a quadrupolar interaction will further split these states such that the spectrum consists of 2I resonances [38]. For the central (I z = +1/2 ↔ −1/2) transition in the presence of magnetic fluctuations, and the coefficients c j are given in Table 1, and the α i are α 1 = 1, α 2 = 6, α 3 = 15, α 4 = 28 and α 5 = 45. For stretched relaxation, however, it is unclear how these functions should be modified. Past researchers have tended to either modify each exponential term with the same stretching exponent:  or simply use the spin-1/2 expression of Eq. 9 [24,39,40]. The problem with these ad hoc approaches is that P (W 1 ) should be independent of the nuclear spin so that it reflects the intrinsic dynamics of the environment, but it is unclear what fitting function should be used. In order to better understand the distributions for higher spin nuclei, we convoluted P β (W 1 ) in Eq. 11 for spin 1/2 with the various Kernel functions in Eq. 13 for different I, and compared with Eq. 14, as shown in Fig. 8. Surprisingly, there is near perfect agreement between the two curves, indicated by the low values of residuals squared across all time points. The larger values of residuals squared at the beginning and end of the time domain are most likely due to the limits of the computed magnetization recovery using Eq. 10 as the exponent term prohibits the recovery to completely reach values of -1 and 1. We further analyzed the various decay curves given by f β (x) (Eq. 14) with the ILT algorithm using the appropriate kernels (given by Eq. 13) to extract P (W 1 ) distributions for each value of I, as shown in Fig. 9. Although there are oscillations introduced by the ILT algorithm, the general shape of the distributions for all of the spins are similar to one another and well-described by P β (W 1 ). Although the relaxation function in Eq. 13 for higher spins is multiexponential, for the central transition the coefficients c j are such that the relaxation is dominated by one exponential. On the other hand, for the satellite transitions (|I z | ↔ |I z | − 1, with 1/2 < |I z | ≤ I), the relative weights of the different exponentials are more evenly distribution. The coefficients c j are given for the different satellite transitions in Table 2 for the case of I = 7/2. In this case the α i are given by α 1 = 1, α 2 = 3, α 3 = 6, α 4 = 10, α 5 = 15, α 6 = 21, and α 7 = 28. Figure 10 compares the extracted distributions with P β (W 1 ), and the magnetization recovery using the convoluted P β (W 1 ) with the stretched expression, Eq. 14, for the central and three satellite transitions for I = 7/2. Once again, there is good agreement. These studies indicate that Eq. 14 is the proper form for stretched exponential relaxation so that the distribution is independent of nuclear spin, I.

Discussion
The ILT algorithm appears to be most effective at accurately capturing the true distribution of relaxation rates when the distribution is sufficiently broad to begin with. For stretched exponential relaxation, when β ≥ 0.8, or when the width of the distribution is less than about one decade, the ILT algorithm overestimates the width. This observation reflects that fact that the Tikhonov regularization acts to smooth the distribution. Efforts to invert noisy data of ill-posed problems typically result in large fluctuations of the distribution function that are not physical, hence the effort to 'regularize' the solution [26,37].
Smoothing of a distribution is a necessary side-effect of the ILT algorithm, and will lead to overestimates of the distribution width when the distribution is intrinsically narrow (such that β ≥ 0.8). Choosing the optimal number of measured recovery points, N , is important to accurately capture the distribution, and our simulations indicate that N should be at least 12-15. Choosing a greater number of points improves the accuracy, but may lead to a reduction in signal to noise if the total experimental time is constrained. Signal to noise ratios above ∼ 40 are necessary to capture the salient features of a distribution, but we find that unphysical artifacts in the distribution persist even up to higher SNR values. This fact should be taken into account when interpreting distributions obtained from experimental systems.
The ILT algorithm artificially introduces discontinuities in dP (W 1 )/dW 1 when P (W 1 ) approaches zero because the Heaviside function forces P to vanish if it becomes negative. This behavior is somewhat problematic because it would be more physically-realistic if the tails of P (W 1 ) asymptotically approached zero smoothly. Such artificial cut-offs may not accurately capture the physics of glassy systems where the fluctuations are expected to exhibit power law distributions with long tails [41,42]. However, our analysis clearly demonstrates that fitting the magnetization recovery directly with the stretched exponential expression (Eq. 14) and inferring the distribution using the theoretical expression (Eq. 11) provides a straightforward description of the distribution of relaxation rates. This approach is valid for any nuclear spin, is easier to implement, and does not suffer from the introduction of artifacts. In such cases, all of the relevant physical information about the distribution is captured by the parameters W * 1 and β. Note that the ILT approach may still be necessary for cases in which the distribution is not well-described by a stretched exponential, for example a bimodal distribution of relaxation rates, or when the distribution is not expected to be well-approximated by P β (W 1 ). The latter distribution is biased towards high W 1 values, and thus stretched exponential fits would not be appropriate when the distribution is expected to be more symmetric or biased towards low W 1 values.

Conclusion
The ILT algorithm is a powerful method to extract distributions from time-series data sets, which has grown in popularity in recent years. A priori, this method makes no assumptions about the nature of the distribution, and is thus useful to study materials with complex inhomogeneous behavior. However, the algorithm does 'filter out' sharp features of a distribution, leading to artificial broadening and oscillatory components. These features are especially pronounced when the time-series data has significant levels of noise. On the other hand, many researchers have traditionally fit the time-series data directly with stretched exponentials of various forms, which are easier to implement and more direct. A drawback of this method has been poor understanding of the nature of the distribution, particularly for the case of I > 1/2. Our study indicates that the stretched exponential form described by Eq. 14 accurately captures the distribution independent of the spin I. This result implies that if the time-series data of any nucleus can be fit by this form, the full distribution P (W 1 ) can be inferred without the need to invert the data with a complicated algorithm such as the ILT.