Fourier Transform of the Stretched Exponential Function: Analytic Error Bounds, Double Exponential Transform, and Open-Source Implementation libkww

The C library \texttt{libkww} provides functions to compute the Kohlrausch-Williams-Watts function, i.e.\ the Laplace-Fourier transform of the stretched (or compressed) exponential function $\exp(-t^\beta)$ for exponents $\beta$ between 0.1 and 1.9 with sixteen-digits accuracy. Analytic error bounds are derived for the low and high frequency series expansions. For intermediate frequencies the numeric integration is enormously accelerated by using the Ooura-Mori double exponential transformation. The source code is available from the project home page \url{http://apps.jcns.fz-juelich.de/doku/sc/kww}.


A. Purpose
Dynamic phenomena can be studied as function of time or frequency. Time-dependent relaxation functions and frequency-dependent spectra are related to each other by Fourier transforms. These transforms are only in special cases given by simple analytic expressions, like the Cauchy spectrum ("Lorentzian") for exponential relaxation. In general, the transforms must be computed numerically. Not rarely, the actual or perceived difficulty of this computation is a hurdle that prevents adequate data analysis.
Specifically, relaxation in disordered systems is often described by the stretched exponential function exp(−(t/τ ) β ). However, when the same dynamics is measured in the frequency domain, experimentalists often resort to other fit functions (like a sum of two Lorentzians or a Havrialiak-Negami function) that should be shaved away by Occam's razor since they introduce additional parameters without providing better fits, let alone deeper insights.
This paper describes the mathematical foundations of a new software that allows for speedy computation of the Fourier transform of the stretched exponential function. For low and high frequencies, well-known series expansions are used. For intermediate frequencies, the Fourier integral is calculated explicitely. This is the bottleneck that determines computing times. In the present implementation, the explicit integration is strongly accelerated thanks to a recent mathematical innovation, the double exponential transformation of Mori et al. [1,2]. * Electronic address: j.wuttke@fz-juelich.de

B. Uses of the stretched exponential
The stretched exponential function arises in different mathematical contexts, for instance as Lévy symmetric alpha-stable distribution, or as the complement of the cumulative Weibull distribution. In recent years, it has been found to provide a good fit to various socioeconomic statistics, like urban agglomeration sizes, currency exchange rate variations, or the 'success' of scientists, musicians, and Hollywood blockbusters [3,4,5].
In physics, the most important use of the stretched exponential function is the approximative description of relaxation in glass-forming liquids and amorphous polymers. In 1993, Böhmer et al. [6] listed stretching exponents for over 70 materials, obtained by viscoelastic, calorimetric, dielectric, optical, and other linear response measurements. As of 2008, that review has been cited over 900 times, indicating a huge increase in the use of the stretched exponential function for describing relaxation phenomena. An impressive body of literature on stretched exponential relaxation in molecular and electronic glasses has been reviewed in 1996 by Phillips [7]; I disagree, however, with his theoretical views, starting on the second page where he claims that β tends towards 1 for high temperatures, which is not confirmed by recent measurements in normal and slightly supercooled liquids (e. g. [8,9]).
Other physical applications of the stretched exponential function are the time dependence of luminiscence or fluorescence decays [10], and the concentration dependence of diffusion coefficients and viscosities [11]. In most applications, the exponent is restricted to values β ≤ 1. However, in recent years some uses of the "compressed" or "squeezed" exponential function with 1 < β < 2 have been proposed, mostly in protein kinetics [12,13,14], but also in magnetism [15].
The use of the Fourier transform to describe dynamic susceptibilities and scattering experiments has its foun-dations in linear response theory. The relations between response functions, relaxation functions, susceptibilities, correlation functions, and scattering laws are briefly summarized in Appendix A.

C. Historic Notes
In physics, the earliest known use of the stretched exponential function is by Rudolf Kohlrausch in 1854 to describe charge relaxation in a Leiden jar. He was followed by his son Friedrich Kohlrausch in 1863 who used the stretched exponential to describe torsional relaxation in glass fibers, thereby improving previous studies by Wilhelm Weber (1841) and Rudolf Kohlrausch (1847). In the modern literature, these early accomplishments are often confounded, and a majority of references to Poggendorff's Annalen der Physik und Chemie is incorrect [16].
In 1970, Williams and Watts introduced the Fourier sine and cosine transform of the stretched exponential function to describe dielectric response as function of frequency [17]. Their intuition is remarkable, since they were neither aware of earlier uses of the stretched exponential in the time domain, nor had they the technical means of actually computing the Fourier transform: based on analytic expressions for β = 1 and β = 0.5, they courageously extrapolated to β = 0.38. In a subsequent paper [18] it is quite obvious that the curves, perhaps drawn with a curving tool, are not really fits to the data.
It was noticed soon that series expansions can be used to calculate the Fourier transform of the stretched exponential in the limit of low or high frequencies [19,20,21]. Based on this, computer routines were implemented that complemented these series expansions by explicit integration for intermediate frequencies [22,23]. In actual fit routines, it was found more convenient to interpolate between tabulated values than to calculate the Fourier transform explicitely [24]. Other experimentalists fit their data with the Havriliak-Negami function, and use some approximations [25,26] to express their results as Kohlrausch-Williams-Watts parameters.

D. Claims
The present implementation improves upon previous work [22,23] in the following respects: (1) The regions where series expansions are applicable are redetermined, and error estimates are corrected. (2) A wider β range is covered. (3) The integration at intermediate frequency, the bottleneck of previous implementations, is improved in speed and accuracy by using the Oouri-Mori double exponential transform. (4) The implementation has a well defined accuracy of 10 −7 . (5) The implementation is made available as a C library libkww that is easy to retrieve, to install, and to use; the source code is released under the GNU General Public License so that long-term availability is ensured.
Claim (2-4) require some explanation. Dishon et al. published tables for values of β between 0.01 and 2 [22]. However, for β 0.3, these tables are useless because they do not cover the ω range where the Fourier transform is nontrivial. This will become clear in Sect. II D.
With respect to speed of calculation, one might oppose that given todays computing power this is no longer a serious concern. However, the Fourier transform of the stretched exponential is often embedded in a convolution of a theoretical model with an instrumental resolution function, which in turn is embedded in a nonlinear curve fitting routine. In such a situation, accelerating the innermost loop is still advantageous.
With respect to accuracy, one might argue that a numeric precision of 10 −3 or 10 −4 is largely sufficient for fitting spectroscopic data. However, violations of monotonicity at a level δ can trap a fit algorithm in a haphazard local minimum unless the minimum search step of the algorithm is correspondingly set to O(δ). By guaranteeing full single-precision accuracy δ = 1.2 · 10 −7 , the integration of libkww into existing fit routines is made easier and safer.

A. Notation
We write the stretched exponential function in dimensionless form as Motivated by the relations between relaxation, linear response, and dynamic susceptibility (Appendix A), we define the one-sided Fourier transform of f β as Strictly speaking, this is a Laplace transform. It can be read as two-sided Fourier transform if the relaxation function f β is multiplied by a Heavyside step function. The frequency ω can be construed as having a small positive imaginary part, ω + i0 + .
In most applications, one is interested in either the cosine or the sine transform, In physical applications, the stretched exponential function is usually written as Its Fourier transform F β,τ can be expressed quite simply by the dimensionless (2):

B. Small ω Expansion
For small and for large values of ω, F β (ω) can be determined from series expansions [19,20,21,22,23]. For small values of ω, one may expand the exp(iωt) term in (2). Substituting x = t β , and using the defining equation of the gamma function, one obtains a Taylor series in powers of ω (in Ref. [21] traced back to Cauchy 1853): with amplitudes For β ≥ 1, these series converge for all values of ω. In practice, however, the use of (7) must be restricted to small |ω|; for large |ω|, the alternating terms become prohibitively large. For β < 1, the series (7) constitute asymptotic expansions [27,28]. In Appendix B it is shown that the modulus of the remainder is not larger than that of the first neglected term. This improves upon a weaker and unproven estimate in Ref. [23]. To minimize the truncation error, the summation must be terminated just before the smallest term.

C. Large ω Expansion
A complementary series expansion for large ω can be derived by expanding the exp(−t β ) term in (2).
(where ± is given by sign a) one obtains a series in powers of ω −β (in Ref. [21] attributed to Wintner 1941 [29]): with amplitudes The convergence of these series is kind of symmetric to the one found in Sect. II B: For β ≤ 1, the series converge for all ω = 0; for β > 1, they (are asymptotic expansions. For small |ω| the are useless because the oscillating terms grow too much before they finally start to converge; for large |ω| they converge rapidly. In the following, we will explicitely consider only the case ω > 0.
To avoid numeric inaccuracies in evaluating the complex exponential or trigonometric factors for β → 2, these terms may be using the complementary exponent β := 2 − β: and similarly for F β and V β . At this occasion we note that there is a qualitative difference between Q β and V β in the special case β = 2 where the large ω expansion is not applicable for Q β because all terms are zero.
Concerning the computation of (10) for β > 1, two erroneous staments are made in Ref. [23], namely: (i) the most accurate results are obtained truncating the summation before the smallest term; and (ii) the truncation error is less than twice the first neglected term.
To see that statement (i) is wrong, let β be a "round" rational number like 3/2, and choose ω small so that B k decreases at least for the first few values of k. Then, for one of these k, the trigonometric factor in Q β (ω) or V β (ω) is zero. If statement (i) were true, (10) could be truncated with zero error to a finite sum. This is obviously not true. A correct truncation criterion must be based not on full terms, but only on amplitudes B k ; it must disregard oscillating prefactors like sin(kβπ/2).
To make sense of statement (ii), it must be reinterpreted as refering not to the first neglected term in its entirety, but only to its amplitude B n−1 . But even then, the statement is unfounded. In Ref. [23] it is underlaid by a reference to a specific page in a book on numerical analysis [30]. However, that page only says "the error committed is usually less than twice the first neglected term" (the emphasis is mine), followed by a reference to a specific page in a 1907 book on Celestial Mechanics [31]. Going back to that source, we find a rigorous theorem that holds only under very restrictive conditions that are not fulfilled here.
A valid estimate of the truncation error has been derived by Wintner [29] for the case β ≤ 1. In Appendix C, a simplified version of his argument is presented, and generalized to the case β > 1. The result is the following: When the sums in (10) are truncated after a term k = n − 1, then the modulus of the truncation error is not larger than For β ≤ 1, sin φ = 1. Therefore, the truncation error is not larger than the first neglected term, disregarding its oscillating prefactor. For 1 < β ≤ 2, the error bound increases with increasing β. In the worst case β = 2, one has a prefactor (sin φ) −1−nβ = 2 1/2+n . The truncation criterion is: stop summation with k = n − 1 if The leading-order terms in (7) and (10) are power-laws in ω. If we plot ln Q β or ln V β versus ln ω, these power-law asymptotes are straight lines. The intersection of a low-ω and a high-ω asymptote defines a cross-over frequency: Fig. 1 shows these cross-over frequencies as function of β. For β 0.5, the curves are very steep; for β → 0, both cross-over frequencies go rapidly to zero, the leading order being Because of this divergence, the limiting case of very small β values has no practical importance, at least in physical applications.
The divergence of the cross-over frequencies (17) also explains why previously published tables [22] of Q β (ω) and V β (ω) are useless for small exponents β 0.3: As these tables employ a fixed, linear ω grid, for β → 0 the nontrivial cross-over region is no longer covered; the tables actually contain no more than the asymptotic large ω power law.

E. Numeric Integration
Popular approaches to calculate numeric Fourier transforms include straightforward fast Fourier transform, and Tuck's simple "Filon-trapezoidal" rule [32]. Both methods evaluate f β (t) on an equidistant grid t k = k∆t. The Filon rule optimizes the weight of the grid points. For small β, the decay of f β extends over several decades. In that case one may either resort to brute force, increasing the number of grid points to the limits of the computers storage capacity, or use a decimation algorithm. The double-exponential transformation provides a simpler and more efficient alternative. Instead of optimizing weights, one optimizes the grid points t k , taking advantage of the zeros of the Fourier integrand.
The double-exponential transformation was first proposed by Takahasi and Mori in 1974 for the efficient evaluation of integrals with end-point singularities [1]. Afterwards, it was improved for the evaluation of oscillatory functions by Ooura and Mori [2]. Their algorithm is particularly efficient for calculating Fourier transforms at intermediate frequencies where series expansions become slow or fail altogether.
The double-exponential transform is a substitution of the integration variable t by a new variable k, using a monotonous function φ(k) that satisfies (for the case ω > 0) With the Fourier integral (2) becomes Applying the trapezoidal rule with stepwidth 1, we obtain an approximation with scale factors and complex phase factors These coefficients can be pre-computed independently of f and ω.
For large negative k, the conditions (18) and (19) ensure that c k goes rapidly towards 1 and b k towards 0. For large positive k, condition (20) implies that the argument of the complex exponential function in (26) tends towards iπk. At this point it is convenient to choose the summation limits according to In consequence, when calculating the cosine transform Q β , all k are half-integers, and therefore Re c k goes to 0. Similarly, for the sine transform V β , all k are integers, and Im c k → 0. So, in any case, the relevant part of c k tends rapidly towards 0 as k → ±∞. This allows us to approximate (22) by (23) with surprisingly low values of |N ± |. All transformations considered by Ooura and Mori [2] have the form .
The function χ(k) must fulfill the conditions Condition (30) guarantees that numerator and denominator of (28) have a zero at the same location k = 0 so that the singularity is removable, The slope at this point is These values are needed to compute the coefficients a 0 , b 0 for use in the sine transform. We note in passing that (26) becomes problematic for large positive k. Following Ref. [2], we use (28) to rewrite c k as For the calculation of Q β , one has to take the real part of (34), with k being half-integer values; for V β , the imaginary part, with integer k. The resulting expressions can both be written in the same form, using the floor function ⌊k⌋: To proceed, the function χ(k) must be specified. Originally, Ooura and Mori [33] had proposed where h is a parameter that controls the accuracy of the trapezoidal approximation. Later, after studying the influence of singularities near the real axis, they suggested [2] Since the stretched exponential function has no such singularities, there is no reason to prefer (37) over (36). Numeric experimentation shows that the integration is rather insensitive to the choice of χ.

A. Accuracy
Let as write S for either Q β (ω) or V β (ω). The numeric computation of S shall aim at a relative accuracy δ. We approximate S by the sums (7), (10), and (23), which we denote for short as In the cases of the small and large ω expansions (7) and (10), we have derived in the appendices B and C upper bounds for the truncation error: The sum (39) is computed incrementally until either the relative error e n /S n is smaller than δ or the computation fails for one of the reasons to be discussed in a moment. The accuracy of the trapezoidal approximation (23) depends on the step width, expressed by the parameter h in (36) or (37). For given h, N − and N + (and thereby n = N + − N − + 1) are chosen such that |b N± | < δ. The trapezoidal sum S n,h is calculated for decreasing step widths h j until either the relative change |(S n,hj − S n,hj−1 )/S n,hj | is smaller than δ or the computation fails.
Computations can fail at any k < n for one of the following reasons: (i) s k has grown beyond the largest floating-point number; (ii) s k has decreased below the smallest normalized floating-point number; (iii) in the case of an asymptotic series, e k+1 > e k : the terms start to grow before the desired acccuracy e n /S n < δ has been reached.
Furthermore, computations can be declared ex post to have failed: (iv) While incrementing (39), we also update a variable that contains the largest amplitude contributing to S n , z n = max k<n |s k |. If z n ≫ S n , then the accuracy of S n suffers from the cancellation of huge alternating terms. Let ǫ denote the accuracy of the floating-point data type used for computing the s k and S n . Then the accuracy of S n is at best ǫz n . We only accept S n if ǫz n < δS n .
In our implementation, we aim at single-precision accuracy δ = 2 −23 ≃ 1.2 · 10 −7 . Internally, s k and S n are double-precision variables with ǫ = 2 −52 ≃ 2.2 · 10 −16 . FIG. 2: Numeric Fourier transforms Q β (ω), V β (ω) over wide logarithmic ranges. Discrete symbols are obtained from the expansions (7) for small ω (blue) and (10) for large ω (red); solid lines are results of numeric quadrature using the double exponential transform. Fig. 2 gives a broad overview over numeric results; Fig. 3 concentrates on the crossover region between the two asymptotic power laws. Blue symbols have been obtained by the small ω expansion (7), red symbols by the large ω expansion (10). For some combinations of ω and β, both expansions fail for one of the four reasons enumerated in Sect. III A. In these cases, one must resort to numeric quadrature based on the double exponential transform (solid lines).   4 shows the large ω limit of the small ω expansion and the small ω limit of the large ω expansion for the case of the cosine transform as function of β. Results for the sine transform are very similar. These limits have been calculated by a simple interval division alogrithm, using a Ruby script kwwlimits.rb that is part of the source distribution. For intermediate values of β, there is only a tiny gap between the two expansions; for the sine transform near β ≃ 1.4, they even overlap. On the other hand, for β 0.3 the gap between the two series expansions extends over more than a decade in ω.

B. Results and Ranges
The data points used in Fig. 4 are copied from the output of kwwlimits.rb into the source of the high-level routines that implement Q β (ω) and V β (ω) for all values of ω and for 0.1 ≤ β ≤ 2. For given β, ω a simple interpolation is used to decide whether a series expansion is likely to work or not. In rare cases (for ω very close to a border line), the decision will be wrong. In case of a false negative, some computation time is wasted doing the numeric quadrature though the series expansion would have worked. In case of a false positive, the subroutine that implements the series expansion will return an error code; then the high-level routine will resort to calling the subroutine that implements the numeric quadrature.
The cosine transform requires special care in the case β → 2. In the limit β = 2, the transform is just a Gaussian, For β 2, Q β is very close to Q 2 for low frequencies, whereas for high frequencies it has a power-law tail In the cross-over range between these two regimes, the numeric quadrature fails to reach the required accuracy because of cancellation. This problem can be solved by quadrating not f β (t) but the difference f β (t)−f 2 (t). The analytic transform Q 2 (ω) is then simply added to the numeric Re FT[f β − f 2 ].

D. Public Interface
Routines for the computation of Q β (ω) and V β (ω) have been implemented in form of a small library libkww. In order to ensure maximum portability, the programming language C has been chosen. The source code is published under the terms of the GNU General Public License (GPL) [34]. The project website is http://www.messen-und-deuten.de/kww. The entry page summarizes the scope of libkww in the style of a Unix manual page. A link leads to the source code directory.
To prepare the source code for distribution, standard procedures of the free software project GNU have been used. After downloading a compressed source package, the standard command sequence $ tar zxvf kww-<version>.tgz $ cd kww-<version> $ ./configure $ make $ sudo make install is sufficient to install libkww on a Unix system.
Following Unix manual page conventions, the application programming interface (API) can be documented by the "synopsis" #include <kww.h> float kwwcf(float omega, float beta); float kwwsf(float omega, float beta); The letters c and s stand for cosine and sine transform, respectively; the routines return Q β (ω) and V β (ω). Following a convention of the C standard library libc, the letter f stands for the data type float of the return value. By choosing this single-precision data type instead of the more common double-precision type double, it is made clear to the user that no attempt has been made to achieve more than single-precision accuracy.
The use of the Fourier transform to describe dynamic susceptibilities and scattering experiments has its foundations in linear response theory. In this appendix, the relations between response functions, relaxation functions, susceptibilities, correlation functions, and scattering laws shall be briefly summarized.
The linear response B(t) to a perturbation A(t) can be written as Consider first the momentary perturbation A(t) = δ(t).
The response is B(t) = R(t). Therefore, the memory kernel R is identified as the response function. Consider next a perturbation A(t) = e ηt Θ(−t) that is slowly switched on and suddenly switched off (Θ is the Heavyside step function, η is sent to 0 + at the end of the calculation). For t > 0, one obtains B(t) = Φ(t) where Φ is the negative primitive of the response function, R(t) = −∂ t Φ(t). Since Φ describes the time evolution after an external perturbation has been switched off, it is called the relaxation function. Kohlrausch's stretched exponential function is a frequently used approximation for Φ(t).
Consider finally a periodic perturbation that is switched on adiabatically, A(t) = exp(−iωt + ηt), implying again the limit η → 0 + . The response can be written B(t) = χ(ω)A(t), introducing a dynamic susceptibility This motivates our definition (2) of the one-sided Fourier transform F (ω) of the relaxation function Φ(t). Since there is a simple relation between χ and F : In consequence, the imaginary part of the susceptibility, which typically describes the loss peak in a spectroscopic experiment, is given by the real part of the Fourier transform of the relaxation function, Im χ = ωRe F β (ω). Up to this point, the only physical input has been Eq. (A1). To make a connection with correlation functions, more substantial input is needed. Using the full apparatus of statistical mechanics (Poisson brackets, Liouville equation, Boltzmann distribution, Yvon's theorem), it is found [35] that for classical systems Pair correlation functions are typically measured in scattering experiments. For instance, inelastic neutron scattering at wavenumber q measures the scattering law S(q, ω), which is the Fourier transform of the density correlation function, In contrast to (2) and (A2), this is a normal, twosided Fourier transform. In consequence, if we let ρ(q, t) * ρ(q, 0) = Φ(t), then the scattering law S(q, ω) is proportional to the real part Re F (ω) of the one-sided Fourier transform of Φ(t).

APPENDIX B: TRUNCATION ERROR IN SMALL ω EXPANSION
In this appendix, we derive an upper bound for the error made by truncating the small ω expansion (7). We consider the cosine transform, and write the Taylor expansion with Lagrange remainder as with 0 ≤ ξ ≤ ω. From (7), we know that Therefore, the truncation error is not larger than the first neglected term. The same holds for the sine transform.

APPENDIX C: TRUNCATION ERROR IN LARGE ω EXPANSION
In this appendix, we derive an upper bound for the error made by truncating the high-ω expansion (10), simplifying an argument of Wintner [29], and generalizing it to cover not only the convergent case β ≤ 1 but also the asymptotic expansion for β > 1. Substituting s = ωt, the Fourier integral (2) takes the form For brevity, we discuss only the cosine transform Q β (ω) = Re F β (ω), which we rewrite as Q β (ω) = G(ω −β )/ω, introducing the functions G(x) = Re Now, we choose an integration path C in the complex plane, consisting of two line segments, s and se iφ , and two arcs, re iϕ and Re iϕ , with 0 < r ≤ s ≤ R < ∞ and 0 ≤ ϕ ≤ φ. The integral of γ along this path is zero: The contributions of the two arcs tend to 0 as r → 0 and R → ∞. Hence the contributions of the two line segments have equal modulus. This allows us to obtain the following bounds: that is independent of ξ. Evaluating the well-known integral (6) we obtain Only trivial changes are needed to adapt this argument to the sin trafo V β (ω).