# Inverse Laplace Transformation Analysis of Stretched Exponential Relaxation

H. Choi111Department of Materials Science and Engineering, Northwestern
University, Evanston, IL 60208, USA<EMAIL_ADDRESS>College of
Nanoscale Science and Engineering, State University of New York Polytechnic
Institute, New York 12203, USA  I. Vinograd  C. Chaffey  N. J. Curro
Department of Physics, University of California, Davis, California 95616, USA

###### Abstract

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.

###### keywords:

Spin lattice relaxation, inhomogeneity, distribution, inverse Laplace
transform

††journal: Journal of Magnetic Resonance

##  1 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}(\mathbf{r})$, where $\mathbf{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})^{\beta}]$, where
$\beta$ is the so-called stretching exponent that provides a rough measures of
the width of $P(W_{1})$, and satisfies $0<\beta\leq 1$ [22, 23]. In many
glassy systems, it has been observed that $\beta$ decreases from unity as the
temperature is reduced [24, 25]. As $\beta$ is reduced, the distribution
$P(W_{1})$ grows by several decades, and $\beta$ 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
distribution 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,\cdots,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 $\beta\leq
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.

##  2 Methods

When a spin $I=1/2$ nucleus at lattice position $\mathbf{r}$ is not in thermal
equilibrium, the magnetization component along the quantization axis
(typically the magnetic field direction) relaxes as:

$m_{z}(\mathbf{r},t)=m_{0}\left(1-\phi e^{-W_{1}(\mathbf{r})t}\right),$ (1)

where $m_{0}$ is the equilibrium magnetization, and $\phi$ is a parameter that
describes the initial condition and depends on the pulse sequence employed in
the measurement. We assume that $W_{1}(\mathbf{r})$ depends on position,
$\mathbf{r}$. All of the nuclei contribute to the measured signal, so that

$\displaystyle M(t)$ $\displaystyle=$
$\displaystyle\int_{V}m_{z}(\mathbf{r},t)d\mathbf{r}$ (2) $\displaystyle=$
$\displaystyle M_{0}\int_{0}^{\infty}K(W_{1},t)P(W_{1})dW_{1},$

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-\phi
e^{-W_{1}t}$. 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$,

$(M(t)-M_{0})/\phi=\int_{0}^{\infty}e^{-W_{1}t}P(W_{1})dW_{1}$ (3)

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:

$M_{i}=M(t_{i})=\sum_{j}K_{ij}P_{j}+e_{i}$ (4)

where $i\in\\{1,\cdots,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\in\\{1,\cdots,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 $\vec{P}$ is underdetermined.
_Tikhonov regularization_ [29] offers a method to obtain a solution by
minimizing the functional:

$\Phi(\vec{P}\,)=\frac{1}{2}|\tilde{\mathbf{K}}\cdot\vec{P}-\vec{M}|^{2}+\frac{1}{2}\alpha\left|\vec{P}\right|^{2}$
(5)

subject to the condition that every element $P_{j}\geq 0$. Here $\alpha$ is
the Tikhonov regularization parameter that enforces $\vec{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 $\alpha$ [27, 37]. The solution of
(5) is
$\vec{P}=\tilde{\mathbf{H}}\cdot\tilde{\mathbf{K}}^{\dagger}\cdot\vec{c}$,
where $\dagger$ means transpose, and the matrix $\tilde{\mathbf{H}}$ has
elements $H_{ij}=H((\tilde{\mathbf{K}}^{\dagger}\cdot\vec{c}\,)_{j})$ if $i=j$
and 0 otherwise, and $H(x)$ is the Heaviside function: $H(x)=1$ if $x>0$ and
$H(x)=0$ otherwise. The vector $\vec{c}$ satisfies:

$\vec{c}(\alpha)=\left(\tilde{\mathbf{K}}\cdot\tilde{\mathbf{H}}(\vec{c}(\alpha))\cdot\tilde{\mathbf{K}}^{\dagger}+\alpha\tilde{\mathbf{I}}\right)^{-1}\cdot\vec{M}$
(6)

where $\mathbf{I}$ is the identity matrix. This equation can be solved
iteratively [36]. The sum of the residuals is:

$\chi(\alpha)=|\vec{M}-\mathbf{K}\cdot\vec{P}|=\alpha|\vec{c}(\alpha)|.$ (7)

The distribution clearly depends on the choice of $\alpha$, and becomes
broader and smoother as $\alpha$ increases. The optimal value of $\alpha$ is
usually determined by the so-called self-consistency method [36, 27], in which
$\alpha$ is chosen as the minimum of either $\alpha_{1}$ or $\alpha_{2}$,
where $\alpha_{1}=|\vec{e}|$, the sum of the experimental errors of the
measurements of $\vec{M}$, and $\alpha_{2}$ satisfies:

$\left.\frac{d\ln\chi(\alpha)}{d\alpha}\right|_{\alpha_{2}}=0.1.$ (8)

We use the IGOR Pro software environment to solve for $\vec{c}$ numerically
using a set of $\vec{M}$ data, for various numbers of $N$ data points, and
with $M=128$ logarithmically-spaced values of $W_{1}$. By computing
$\chi(\alpha)$ for a broad range of $\alpha$, we find the optimal
regularization parameter and use this to determine the distribution $\vec{P}$
for a given data set, $\\{t_{i},M_{i}\\}$ with measurements errors $e_{i}$.

##  3 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\sim 5-20$. However, the major constraint
is the total experimental time, $t_{expt}=Nt_{1}$, where $t_{1}$ is the
measurement time for a single point. For Gaussian noise, the measurement error
$e_{i}\sim t_{1}^{-1/2}\sim 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.

###  3.1 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:

$K(t,W_{1})=1-2e^{-(W_{1}t)},$ (9)

with the corresponding magnetization recovery vector:

$M(t_{i})=1-2e^{-(W^{*}_{1}t_{i})^{\beta}}+e_{i},$ (10)

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_{\beta}(W_{1})$ can be expressed analytically as an infinite series:

$P_{\beta}(W_{1})=\frac{1}{\pi}\sum_{n=1}^{\infty}\frac{(-1)^{n+1}\Gamma(n\beta+1)}{n!(W_{1}/W_{1}^{*})^{n\beta+1}}\sin(n\pi\beta)$
(11)

where $\Gamma(x)$ is the Gamma function [23]. Figure 1 shows the distributions
corresponding to stretched relaxation for several values of $\beta$. Note that
because the $W_{1}$ values are distributed on a logarithmic scale, it is
necessary to multiply $P_{\beta}(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 $\beta\rightarrow 1$. As $\beta$ is reduced,
the distribution broadens considerably, and is several decades in width once
$\beta\leq 0.5$.

Figure 1:  (a) Linear-log plot of spin I = 1/2 M(t) with stretched exponential
form given by Eq. 10 with $N=15$ time points and $e_{i}=0$ for $\beta$ ranging
from 0.3 to 1 in 0.1 increments. (b) Linear-log plot of the theoretical
distribution $P_{\beta}(W_{1})$ given by Eq. 11 with corresponding $\beta$
values.

Fig. 2 shows the distributions extracted using the ILT method for several
values of $\beta$ with zero noise ($e_{i}=0$) and with $N=15$ time points.
There is relatively good agreement for intermediate values of $\beta$, but
once $\beta\gtrsim 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:

$S_{2}=\sum_{j=1}^{M}\left(P_{j}-P_{\beta}(W_{1j})\right)^{2}$ (12)

where $\vec{P}$ is determined by ILT. As shown in Fig. 3, $S_{2}$ generally
increases as $\beta$ approaches unity. However, $S_{2}$ has inflated values at
$\beta<0.5$ due to domain constraints. For lower $\beta$ 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 $\beta$ recoveries, resulting in larger $S_{2}$ values.

Figure 2:  Linear-log plot of $P(W_{1})$ ILT estimations of M(t) given by Eq.
10 with $N=15$ time points and $e_{i}=0$ and the theoretical distribution
$P_{\beta}(W_{1})$ given by Eq. 11 for $\beta$ values ranging from 0.3 to 1.0
in 0.1 increments. Figure 3:  Log-linear plot of $S_{2}$, given by Eq. 12, of
P(W1) ILT estimations of M(t) given by Eq. 10 with $N=15$ time points, and
$e_{i}=0$ for $\beta$ values ranging from 0.3 to 1.0 in 0.1 increments versus
$\beta$.

###  3.2 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 $\beta=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\gtrsim 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 $\beta$ 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.

Figure 4:  Linear-log plot of P(W1) ILT estimations of M(t) given by Eq. 10
with $N=6,15,30$ time points and $e_{i}=0$ and the theoretical distribution
$P_{\beta}(W_{1})$ given by Eq. 11 for = 0.8. Figure 5: Plot of $S_{2}$, given
by Eq. 12, of P(W1) ILT estimations of M(t) with $N$ ranging from 6 to 30 time
points and $e_{i}=0$ versus the number of M(t) recovery points for = 0.7, 0.8,
0.9.

###  3.3 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 $\sigma_{n}^{2}$,
such that the signal to noise ratio $SNR=\sigma_{n}^{-1}$. A sample set of
magnetization recovery points $(N=15)$ with $\beta=0.8$, and the corresponding
ILT distributions are shown in Fig. 6.

Figure 6: (a) Linear-log plot of M(t) given by Eq. 10 with $N=15$ time points,
SNR = 10, 20, 100, and = 0.8. (b) Linear-log plot of comparison between P(W1)
ILT estimations of M(t) in (a) and the theoretical distribution
$P_{\beta}(W_{1})$ given by Eq. 11 for = 0.8.

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}\sim SNR^{-s}$, where $s$ increases
as $\beta$ decreases. Extra structures, such as spurious peaks and shoulders
are apparent in the ILT distributions. These artifacts persist even up to SNR
levels of 100, although the overall shape of the distribution is qualitatively
correct.

Figure 7:  Log-log plot of $S_{2}$, given by Eq. 12 of P(W1) ILT estimations
of M(t) given by Eq. 10 with $N=15$ time points and SNR ranging from 10 to 100
versus SNR for = 0.7, 0.8, 0.9.

###  3.4 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\leftrightarrow-1/2$)
transition in the presence of magnetic fluctuations, $K(W_{1},t)=1-\phi
f(W_{1}t)$, where $f(x)$ is given by:

$f(x)=\sum_{j}c_{j}\exp({-\alpha_{j}x})$ (13)

and the coefficients $c_{j}$ are given in Table 1, and the $\alpha_{i}$ are
$\alpha_{1}=1$, $\alpha_{2}=6$, $\alpha_{3}=15$, $\alpha_{4}=28$ and
$\alpha_{5}=45$.

Table 1:  Coefficients in $f(x)$ for higher spins. | $I={3}/{2}$ | $I={5}/{2}$ | $I={7}/{2}$ | $I={9}/{2}$
---|---|---|---|---
$c_{1}$ | ${1}/{10}$ | ${1}/{35}$ | 1/84 | 1/165
$c_{2}$ | 9/10 | 8/45 | 3/44 | 24/715
$c_{3}$ | 0 | 50/63 | 75/364 | 6/65
$c_{4}$ | 0 | 0 | 1224/1716 | 1568/7293
$c_{5}$ | 0 | 0 | 0 | 7938/12155

Figure 8: (a) Comparison between the computed magnetization recovery $M(t)$
(Eq. 3) with $P_{\beta}(W_{1})$ given by Eq. 11 for various values of the
nuclear spin $I$ with N = 15 time points. (b) Corresponding residuals squared
between $M(t)$ and the stretched exponential recovery function (Eq. 14).

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:

$f_{\beta}(x)=\sum_{j}c_{j}\exp({-(\alpha_{j}x)^{\beta}}),$ (14)

or simply use the spin-1/2 expression of Eq. 9 [39, 24, 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_{\beta}(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_{\beta}(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_{\beta}(W_{1})$.

Table 2:  Coefficients in $f(x)$ for different satellite transitions of a spin 7/2 nucleus for magnetic fluctuations. | $|{1}/{2}|\leftrightarrow|3/2|$ | $|{3}/{2}|\leftrightarrow|5/2|$ | $|{5}/{2}|\leftrightarrow|7/2|$
---|---|---|---
$c_{1}$ | 196/429 | 49/429 | 4/429
$c_{2}$ | 49/132 | 49/132 | 3/44
$c_{3}$ | 1/1092 | 100/273 | 75/364
$c_{4}$ | 9/77 | 25/308 | 25/77
$c_{5}$ | 1/33 | 1/132 | 3/11
$c_{6}$ | 1/84 | 1/21 | 3/28
$c_{7}$ | 1/84 | 1/84 | 1/84

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}|\leftrightarrow|I_{z}|-1$, with
$1/2<|I_{z}|\leq 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
$\alpha_{i}$ are given by $\alpha_{1}=1$, $\alpha_{2}=3$, $\alpha_{3}=6$,
$\alpha_{4}=10$, $\alpha_{5}=15$, $\alpha_{6}=21$, and $\alpha_{7}=28$. Figure
10 compares the extracted distributions with $P_{\beta}(W_{1})$, and the
magnetization recovery using the convoluted $P_{\beta}(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$.

Figure 9: Distributions of $P(W_{1})$ extracted from stretched relaxation
curves using $f_{\beta}(x)$ (Eq. 14) for $I=1/2,3/2,5/2,7/2$ and 9/2, for
$\beta=0.8$. The dashed line is the theoretical distribution
$P_{\beta}(W_{1})$(Eq. 11).

##  4 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 $\beta\geq
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
$\beta\geq 0.8$).

Figure 10:  (Upper panel) Distributions extracted from stretched relaxation
curves using $f_{\beta}(x)$ (Eq. 14) for the central
($1/2\leftrightarrow-1/2$), first ($|1/2|\leftrightarrow|3/2|$), second
($|3/2|\leftrightarrow|5/2|$) and third satellites
($|5/2|\leftrightarrow|7/2|$) for $I=7/2$. (Lower panel) Comparison between
the computed magnetization recovery and the stretched expression (Eq. 14) for
the same transitions for $\beta=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 $\sim 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 $\vec{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 $\beta$. 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_{\beta}(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.

##  5 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.

Acknowledgment. We acknowledge stimulating discussions with W. Polonik, K.
Dahmen and E. Carlson. Work at UC Davis was supported by the NSF under Grants
No. DMR-1807889 and PHY-1852581.

## References

  * [1] T. Moriya, The effect of electron-electron interaction on the nuclear spin relaxation in metals, J. Phys. Soc. Jpn. 18 (4) (1963) 516–520.
URL http://dx.doi.org/10.1143/JPSJ.18.516

  * [2] P. C. Hammel, M. Takigawa, R. H. Heffner, Z. Fisk, K. C. Ott, Spin dynamics at oxygen sites in YBa2Cu3O7, Phys. Rev. Lett. 63 (18) (1989) 1992–1995.
URL http://dx.doi.org/10.1103/PhysRevLett.63.1992

  * [3] M. Bankay, M. Mali, J. Roos, I. Mangelschots, D. Brinkmann, Anisotropy and magnetic field dependence of the planar copper NMR spin-lattice relaxation rate in YBa2Cu4O8, Phys. Rev. B 46 (17) (1992) 11228–11231.
URL http://dx.doi.org/10.1103/PhysRevB.46.11228

  * [4] T. Imai, C. P. Slichter, K. Yoshimura, K. Kosuge, Low frequency spin dynamics in undoped and Sr-doped La2CuO4, Phys. Rev. Lett. 70 (7) (1993) 1002–1005.
URL https://doi.org/10.1103/PhysRevLett.70.1002

  * [5] S. Ohsugi, Y. Kitaoka, K. Ishida, K. Asayama, Cu NQR study of the spin dynamics in high-$T_{c}$ superconductor La2-xSrxCuO4, J. Phys. Soc. Jpn. 60 (7) (1991) 2351–2360.
URL journals.jps.jp/doi/abs/10.1143/JPSJ.60.2351

  * [6] M.-H. Julien, M. Horvatić, P. Carretta, C. Berthier, Y. Berthier, P. Ségransan, S. Loureiro, J.-J. Capponi, 63Cu and 199Hg NMR in overdoped HgBa2CaCu2O6+δ, Phys. C 268 (3-4) (1996) 197–204.
URL http://dx.doi.org/10.1016/0921-4534(96)00419-4

  * [7] R. R. Urbano, B.-L. Young, N. J. Curro, J. D. Thompson, L. D. Pham, Z. Fisk, Interacting antiferromagnetic droplets in quantum critical CeCoIn5, Phys. Rev. Lett. 99 (14) (2007) 146402.  doi:10.1103/PhysRevLett.99.146402.
URL http:/dx.doi.org/10.1103/PhysRevLett.99.146402

  * [8] Y. Kohori, T. Kohara, H. Shibai, Y. Oda, Y. Kitaoka, K. Asayama, Nuclear magnetic relaxation in the heavy-fermion superconductor UPt3, J. Phys. Soc. Jpn. 57 (2) (1988) 395–397.
URL http://dx.doi.org/10.1143/JPSJ.57.395

  * [9] G. Q. Zheng, N. Yamaguchi, H. Kan, Y. Kitaoka, J. L. Sarrao, P. G. Pagliuso, N. O. Moreno, J. D. Thompson, Coexistence of antiferromagnetic order and unconventional superconductivity in heavy-fermion CeRh1-xIrxIn5 compounds: Nuclear quadrupole resonance studies, Phys. Rev. B 70 (1) (2004) 014511.
URL http:/dx.doi.org/10.1103/PhysRevB.70.014511

  * [10] K. Kitagawa, N. Katayama, K. Ohgushi, M. Yoshida, M. Takigawa, Commensurate itinerant antiferromagnetism in BaFe2As2: 75As-NMR studies on a self-flux grown single crystal, J. Phys. Soc. Jpn. 77 (11) (2008) 114709.
URL http://jpsj.ipap.jp/link?JPSJ/77/114709/

  * [11] K. Kitagawa, N. Katayama, H. Gotou, T. Yagi, K. Ohgushi, T. Matsumoto, Y. Uwatoko, M. Takigawa, Spontaneous formation of a superconducting and antiferromagnetic hybrid state in SrF$e_{2}$As2 under high pressure, Phys. Rev. Lett. 103 (25) (2009) 257002\.
URL http://dx.doi.org/10.1103/PhysRevLett.103.257002

  * [12] T. Kissikov, R. Sarkar, M. Lawson, B. T. Bush, E. I. Timmons, M. A. Tanatar, R. Prozorov, S. L. Bud’ko, P. C. Canfield, R. M. Fernandes, N. J. Curro, Uniaxial strain control of spin-polarization in multicomponent nematic order of BaFe2As2, Nat. Commun. 9 (1) (2018) 1058.
URL http://dx.doi.org/10.1038/s41467-018-03377-8

  * [13] S. H. Baek, M. J. Graf, A. V. Balatsky, E. D. Bauer, C. Cooley, J. L. Smith, N. J. Curro, Antiferromagnetic patches and hidden order in URu2Si2 by impurity doping, Phys. Rev. B 81 (13) (2010) 132404.
URL http://dx.doi.org/10.1103/PhysRevB.81.132404

  * [14] J. M. Tranquada, B. J. Sternlieb, J. D. Axe, Y. Nakamura, S. Uchida, Evidence for stripe correlations of spins and holes in copper oxide, Nature 375 (6532) (1995) 561–563.
URL http://dx.doi.org/10.1038/375561a0

  * [15] S. A. Kivelson, E. Fradkin, V. J. Emery, Electronic liquid-crystal phases of a doped Mott insulator, Nature 393 (6685) (1998) 550–553.
URL http://dx.doi.org/10.1038/31177

  * [16] H. Alloul, J. Bobroff, M. Gabay, P. J. Hirschfeld, Defects in correlated metals and superconductors, Rev. Mod. Phys. 81 (1) (2009) 45.
URL http://link.aps.org/abstract/RMP/v81/p45

  * [17] J. Schmalian, P. G. Wolynes, Stripe glasses: Self-generated randomness in a uniformly frustrated system, Phys. Rev. Lett. 85 (4) (2000) 836–839.
URL http://dx.doi.org/10.1103/PhysRevLett.85.836

  * [18] S. Seo, X. Lu, J.-X. Zhu, R. R. Urbano, N. Curro, E. D. Bauer, V. A. Sidorov, L. D. Pham, T. Park, Z. Fisk, J. D. Thompson, Disorder in quantum critical superconductors, Nat. Phys. 10 (2014) 120–125.
URL http://dx.doi.org/10.1038/nphys2820

  * [19] A. P. Dioguardi, M. M. Lawson, B. T. Bush, J. Crocker, K. R. Shirer, D. M. Nisson, T. Kissikov, S. Ran, S. L. Bud’ko, P. C. Canfield, S. Yuan, P. L. Kuhns, A. P. Reyes, H.-J. Grafe, N. J. Curro, NMR evidence for inhomogeneous glassy behavior driven by nematic fluctuations in iron arsenide superconductors, Phys. Rev. B 92 (2015) 165116.
URL http://link.aps.org/doi/10.1103/PhysRevB.92.165116

  * [20] M. Horvatić, Y. Fagot-Revurat, C. Berthier, G. Dhalenne, A. Revcolevschi, NMR imaging of the soliton lattice profile in the spin-Peierls compound CuGeO3, Phys. Rev. Lett. 83 (2) (1999) 420–423.
URL http://dx.doi.org/10.1103/PhysRevLett.83.420

  * [21] Y. Nakai, Y. Hayashi, K. Kitagawa, K. Ishida, H. Sugawara, D. Kikuchi, H. Sato, Evidence of the bound states of the vortex state in an s-wave superconductor proved by NMR measurements, J. Phys. Soc. Japan 77 (Suppl.A) (2008) 333–335.
URL http://dx.doi.org/10.1143/jpsjs.77sa.333

  * [22] J. C. Phillips, Stretched exponential relaxation in molecular and electronic glasses, Rep. Progr. Phys. 59 (9) (1996) 1133–1207.
URL http://dx.doi.org/10.1088/0034-4885/59/9/003

  * [23] D. C. Johnston, Stretched exponential relaxation arising from a continuous sum of exponential decays, Phys. Rev. B 74 (2006) 184430.
URL http://link.aps.org/doi/10.1103/PhysRevB.74.184430

  * [24] V. F. Mitrović, M.-H. Julien, C. de Vaulx, M. Horvatić, C. Berthier, T. Suzuki, K. Yamada, Similar glassy features in the 139La NMR response of pure and disordered La1.88Sr0.12CuO4, Phys. Rev. B 78 (2008) 014504.
URL http://link.aps.org/doi/10.1103/PhysRevB.78.014504

  * [25] A. P. Dioguardi, T. Kissikov, C. H. Lin, K. R. Shirer, M. M. Lawson, H.-J. Grafe, J.-H. Chu, I. R. Fisher, R. M. Fernandes, N. J. Curro, NMR evidence for inhomogeneous nematic fluctuations in BaFe2(As1-xPx)2, Phys. Rev. Lett. 116 (2016) 107202.
URL http://link.aps.org/doi/10.1103/PhysRevLett.116.107202

  * [26] J. Butler, J. Reeds, S. Dawson, Estimating solutions of first kind integral equations with nonnegative constraints and optimal smoothing, SIAM J. Numer. Anal. 18 (3) (1981) 381–397.
URL https://doi.org/10.1137/0718025

  * [27] J. Honerkamp, J. Weese, Tikhonovs regularization method for ill-posed problems, Contin. Mech. Thermodyn. 2 (1) (1990) 17–30.
URL https://doi.org/10.1007/BF01170953

  * [28] A. Livesey, J. Brochon, Analyzing the distribution of decay constants in pulse-fluorimetry using the maximum entropy method, Biophysical journal 52 (5) (1987) 693–706.
URL https://doi.org/10.1016/S0006-3495(87)83264-2

  * [29] A. N. Tikhonov, Regularization of incorrectly posed problems, Soviet Mathematics Doklady, 1963.
  * [30] E. Fordham, A. Sezginer, L. Hall, Imaging multiexponential relaxation in the ($y$, Log${}_{e}T_{1}$) plane, with application to clay filtration in rock cores, Journal of Magnetic Resonance, Series A 113 (2) (1995) 139–150.
URL https://doi.org/10.1006/jmra.1995.1073

  * [31] L. Venkataramanan, Y.-Q. Song, M. D. Hurlimann, Solving Fredholm integrals of the first kind with tensor product structure in 2 and 2.5 dimensions, IEEE Transactions on signal processing 50 (5) (2002) 1017–1026.
URL https://doi.org/10.1109/78.995059

  * [32] Y.-Q. Song, L. Venkataramanan, M. Hürlimann, M. Flaum, P. Frulla, C. Straley, T1–T2 correlation spectra obtained using a fast two-dimensional Laplace inversion, J. Magn. Reson. 154 (2) (2002) 261–268.
URL http://dx.doi.org/10.1006/jmre.2001.2474

  * [33] J. Mitchell, T. Chandrasekera, L. Gladden, Numerical estimation of relaxation and diffusion distributions in two dimensions, Prog Nucl Mag Res Sp 62 (2012) 34–50.
URL http://dx.doi.org/10.1016/j.pnmrs.2011.07.002

  * [34] H. Schäfer, E. Sternin, R. Stannarius, M. Arndt, F. Kremer, Novel approach to the analysis of broadband dielectric spectra, Phys. Rev. Lett. 76 (1996) 2177–2180.
URL https://link.aps.org/doi/10.1103/PhysRevLett.76.2177

  * [35] A. Arsenault, T. Imai, P. Singer, K. Suzuki, M. Fujita, Magnetic inhomogeneity in charge-ordered La1.885Sr0.115CuO4 studied by NMR, Phys. Rev. B 101 (18) (2020) 184505.
URL https://doi.org/10.1103/PhysRevB.101.184505

  * [36] P. Singer, A. Arsenault, T. Imai, M. Fujita, La-139 NMR investigation of the interplay between lattice, charge, and spin dynamics in the charge-ordered high-Tc cuprate La1.875Ba0.125CuO4, Phys. Rev. B 101 (17) (2020) 174508.
URL https://doi.org/10.1103/PhysRevB.101.174508

  * [37] C. L. Lawson, R. J. Hanson, Solving least squares problems, SIAM, 1995.
  * [38] C. P. Slichter, Principles of Nuclear Magnetic Resonance, 3rd Edition, Springer-Verlag, 1992.
  * [39] N. Curro, P. Hammel, B. Suh, M. Hücker, B. Büchner, U. Ammerahl, A. Revcolevschi, Inhomogeneous low frequency spin dynamics in La1.65Eu0.2Sr0.15CuO4, Phys. Rev. Lett. 85 (3) (2000) 642–645.
URL http://dx.doi.org/10.1103/PhysRevLett.85.642

  * [40] A. P. Dioguardi, J. Crocker, A. C. Shockley, C. H. Lin, K. R. Shirer, D. M. Nisson, M. M. Lawson, N. apRoberts Warren, P. C. Canfield, S. L. Bud’ko, S. Ran, N. J. Curro, Coexistence of cluster spin glass and superconductivity in Ba(Fe1-xCox)2As2 for $0.060\leq x\leq 0.071$, Phys. Rev. Lett. 111 (2013) 207201.
URL http://link.aps.org/doi/10.1103/PhysRevLett.111.207201

  * [41] P. Sollich, F. Lequeux, P. Hébraud, M. E. Cates, Rheology of soft glassy materials, Physical review letters 78 (10) (1997) 2020.
URL https://doi.org/10.1103/PhysRevLett.78.2020

  * [42] H. Fujisaka, K. Egami, T. Yamada, Glassy dynamics in a spatially distributed dynamical system, Phys. Lett. A 174 (1-2) (1993) 103–110.
URL https://doi.org/10.1016/0375-9601(93)90550-J