Fast logarithmic Fourier-Laplace transform of nonintegrable functions

TL;DR

This paper introduces a fast, flexible numerical Fourier-Laplace transform method that efficiently handles nonintegrable functions with power-law tails, demonstrating exponential convergence and broad applicability in physical models.

Contribution

It extends the logarithmic Fourier transform to nonintegrable functions, providing a method with proven exponential convergence and practical guidelines for applications involving power-law tails.

Findings

Achieves exponential convergence for a broad class of functions.

Effectively handles functions with power-law asymptotic behavior.

Demonstrates applicability in physical models near criticality.

Abstract

We present an efficient and very flexible numerical fast Fourier-Laplace transform, that extends the logarithmic Fourier transform (LFT) introduced by Haines and Jones [Geophys. J. Int. 92(1):171 (1988)] for functions varying over many scales to nonintegrable functions. In particular, these include cases of the asymptotic form $f(\nu\to0)\sim\nu^a$ and $f(|\nu|\to\infty)\sim\nu^b$ with arbitrary real $a>b$. Furthermore, we prove that the numerical transform converges exponentially fast in the number of data points, provided that the function is analytic in a cone $|\Im{\nu}|<\theta|\Re{\nu}|$ with a finite opening angle $\theta$ around the real axis and satisfies $|f(\nu)f(1/\nu)|<\nu^c$ as $\nu\to 0$ with a positive constant $c$, which is the case for the class of functions with power-law tails. Based on these properties we derive ideal transformation parameters and discuss how the logarithmic Fourier transform can be applied to convolutions. The ability of the logarithmic Fourier transform to perform these operations on multiscale (non-integrable) functions with power-law tails with exponentially small errors makes it the method of choice for many physical applications, which we demonstrate on typical examples. These include benchmarks against known analytical results inaccessible to other numerical methods, as well as physical models near criticality.

Equations

$$f(\nu )=\mathcal{F}(\hat{f})(\nu )={\int }_{-\infty }^{\infty }dt\hat{f}(t)e^{i\nu t},$$

$$\hat{f}(t)={\mathcal{F}}^{-1}(f)(t)={\int }_{-\infty }^{\infty }\frac{d\nu }{2\pi }f(\nu )e^{-i\nu t},$$

$$\nu =\sigma \bar{\nu }{e}^{\omega }\text{and}t=\eta \bar{t}{e}^{\tau },$$

$$\begin{array}{rl}& \hat{f}(\eta |t|)=e^{-k\tau }\\ & {\times \sum _{\sigma =\pm 1}\int \frac{d\omega }{2\pi }f(\sigma e^{\omega })e^{k(\omega +\tau )-i\sigma \eta \exp (\omega +\tau )}{e}^{(1-k)\omega }|}_{\tau =\ln |t|},\end{array}$$

$$\hat{f}(\eta |t|)=\frac{e^{-k\tau }}{2\pi }\sum _{\sigma =\pm 1}{\mathcal{F}}_{s\to \tau }\left[{\mathcal{F}}_{\omega \to s}\left(f(\sigma e^{\omega })e^{(1-k)\omega }\right)(s){\mathcal{F}}_{x\to s}^{-1}\left(e^{kx-i\sigma \eta \exp (x)}\right)(s)\right]\left(\tau =\ln |t|\right),$$

$$F_{\sigma }(\omega ):=f(\sigma e^{\omega })e^{(1-k)\omega }\in L^1,$$

$$g_{\sigma }(s):={\mathcal{F}}_{\omega \to s}\left(f(\sigma e^{\omega })e^{(1-k)\omega }\right)(s)$$

$${\int }_0^{\infty }d\nu \left|f(\sigma \nu )\right||\nu {|}^{-k}<\infty .$$

$$1+b<k<1+a.$$

$$\begin{array}{rl}{h}_{\sigma \eta }(s):=& {\mathcal{F}}_{x\to s}^{-1}\left(e^{kx-i\eta \sigma \exp (x)}\right)(s)\\ =& \frac{1}{2\pi }{\left(i\sigma \eta \right)}^{is-k}\Gamma (k-is),\end{array}$$

$$\int dx{e}^{kx-i\eta \sigma \exp \left(x\right)}{e}^{-isx}={\left(i\sigma \eta \right)}^{is-k}\Gamma (k-is)$$

$$\left|\Gamma (k-is)(i\sigma \eta {)}^{is-k}\right|\propto \{\begin{array}{ll}\sqrt{2\pi }|s{|}^{k-1/2}{e}^{-\pi |s|}& \sigma \eta s\to \infty \\ \sqrt{2\pi }|s{|}^{k-1/2}& \sigma \eta s\to -\infty \end{array},$$

$$n:=\max (0,\lceil k+1/2\rceil )$$

$$\hat{f}(\eta |t|)=\frac{e^{-k\tau }}{(2\pi {)}^2}\sum _{\sigma =\pm 1}{\mathcal{F}}_{s\to \tau }\left[{\mathcal{F}}_{\omega \to s}\left(f(\sigma e^{\omega })e^{(1-k)\omega }\right)(i\sigma \eta {)}^{is-k}\Gamma (k-is)\right]\left(\tau =\ln |t|\right),$$

$$\mathcal{F}\mathcal{L}(f)(t)={\int }_{-\infty }^{\infty }\frac{d\nu }{2\pi }f(\nu )\exp (e^{i\varphi }\nu t),$$

$$\begin{array}{rl}{\nu }_n& =e^{{\omega }_n}\text{and}{\nu }_{-n}=-e^{{\omega }_n}\text{with}{\omega }_n=\Delta \omega \left(n+{\omega }_s\right)\\ t_n& =e^{{\tau }_n}\text{and}{t}_{-n}=-e^{{\tau }_n}\text{with}{\tau }_n=\Delta \tau \left(n+{\tau }_s\right)\end{array}$$

$$\begin{array}{rl}{\hat{f}}_N(t_{\eta n})=& e^{-k{\tau }_n}\sum _{\sigma =\pm 1}\sum _{l=1}^N\frac{\Delta s}{2\pi }{e}^{i{s}_l{\tau }_n}(i\sigma \eta {)}^{i{s}_l-k}\Gamma (k-i{s}_l)\\ & \cdot \sum _{m=1}^N\frac{\Delta \omega }{2\pi }f(\sigma e^{{\omega }_m})e^{(1-k){\omega }_m}{e}^{i{\omega }_m{s}_l}\end{array}$$

$$I_1^{\sigma }(s)=\int \frac{d\omega }{2\pi }f(\sigma e^{\omega })e^{(1-k)\omega }{e}^{is\omega },$$

$$\begin{array}{rl}& S_1^{\sigma }(s)=\sum _{m\in \mathbb{Z}}\frac{\Delta \omega }{2\pi }f(\sigma e^{{\omega }_n})e^{(1-k){\omega }_n}{e}^{is{\omega }_n}\\ & =\{\begin{array}{ll}\oint \frac{dz}{2\pi i}f(\sigma e^{-iz})(1+n_B(z))e^{-(1-k+is)iz}& \Re (s)<0\\ \oint \frac{dz}{2\pi i}f(\sigma e^{-iz})(n_B(z))e^{-(1-k+is)iz}& \Re (s)>0\end{array}.\end{array}$$

$$\begin{array}{rl}& E_1^{\sigma }(s)=S_1^{\sigma }(s)-I_1^{\sigma }(s)=\\ & ={\int }_{R^{(1)}-i\infty }^{R^{(1)}-+i\infty }\frac{dz}{2\pi i}f(\sigma e^{-iz})(n_B(z))e^{-(1-k+is)iz}\\ & +{\int }_{-R^{(1)}+i\infty }^{-R^{(1)}-i\infty }\frac{dz}{2\pi i}f(\sigma e^{-iz})(1+n_B(z))e^{-(1-k+is)iz},\end{array}$$

$$E_1^{\sigma }(s)=e^{R^{(1)}(-{\beta }_{\omega }+|s|)}{F}_1^{\sigma }(s,R^{(1)}).$$

$$I_2^{\sigma }(\tau )=\int \frac{ds}{2\pi }\Gamma (k-is)(i\sigma \eta {)}^{is-k}{I}_1^{\sigma }(s)e^{is\tau }.$$

$$\begin{array}{rl}{E}_2^{\sigma }(\tau )=& I_2^{\sigma }(\tau )\\ & -\sum _{|l|\le N}\frac{\Delta s}{2\pi }\Gamma (k-i{s}_l)e^{i{s}_l\tau }(-i\sigma \eta {)}^{i{s}_l-k}{S}_1^{\sigma }(s_l).\end{array}$$

$$E_2^{\sigma }(\tau )=e^{R^{(2)}(-{\beta }_s+|\tau |)}{F}_2^{\sigma }(\tau ,R^{(2)})+E_{2,\Gamma }^{\sigma }(\tau ).$$

$$\begin{array}{rl}& E_{2\Gamma }^{\sigma }(\tau )=-\sum _{\begin{array}{c}m=\lceil -k\rceil \\ \wedge m\ge 0\end{array}}^{\lfloor R^{(2)}-k\rfloor }\frac{(-1{)}^m}{m!}\frac{e^{(k+m)\tau }(i\sigma \eta {)}^m{I}_1^{\sigma }(-i(m+k))}{e^{{\beta }_s(k+m)}-1}\\ & -\sum _{\begin{array}{c}m=\lceil -R^{(2)}-k\rceil \\ \wedge m\ge 0\end{array}}^{\lfloor -k\rfloor }\frac{(-1{)}^m}{m!}\frac{e^{(k+m)\tau }(i\sigma \eta {)}^m{I}_1^{\sigma }(-i(k+m))}{1-e^{-{\beta }_s(k+m)}}.\end{array}$$

$$1\le \frac{2\Delta \omega }{{\pi }^2}\ln \left(\sqrt{2\pi }\frac{ϵ}{\delta }\right).$$

$$\left|f(\sigma |\nu |)f\left(\frac{\sigma }{|\nu |}\right)\right|\underset{\nu \to 0}{<}|\nu {|}^{c_{\sigma }}.$$

$$f(\sigma e^{\omega })=\sum _{n=0}^{\infty }{\left(\frac{d^n}{d{\omega }^n}f(\sigma e^{\omega })\right)}_{{\omega }_0}\frac{(\omega -{\omega }_0{)}^n}{n!}$$

$$R^{(1)}=\underset{\begin{array}{c}{\omega }_0\in \bar{\mathbb{R}}\\ \sigma =\pm 1\end{array}}{\inf }{R}_{{\omega }_0}^{\sigma }.$$

$$\theta (R^{(1)})=\min (\pi /2,R^{(1)})>0.$$

Full text

I Introduction

In physics, one is often confronted with the need to Fourier transform or convolve functions that are either only numerically available or whose exact transformation is not known. Since the reinvention of the fast Fourier transform (FFT) by Cooley and Tukey FFT1965 , which reduces the numerical cost for both of these operations from $\mathcal{O}(N^{2})$ to $\mathcal{O}(N\log_{2}{N})$, where $N$ denotes the number of grid points, the FFT has been established as the standard method for most situations. However, it necessarily requires an equidistant grid, which is quite inconvenient for many applications in theoretical physics. There, one frequently has to deal with slowly (i.e. algebraically) decaying functions, while the opposite limit of small arguments contains a lot of physical information. An example is provided by Green’s functions in many-body problems with short-range interactions zwer14varenna . To implement an FFT under such circumstances, it is necessary to use a fine grid for small arguments that extends to very high frequencies, which is of course not very practicable due to the huge number of required data points. Consequently, a number of alternative methods have been introduced in the literature: Sometimes, sufficient knowledge about the asymptotic behavior at large arguments can be gained, subtracted and treated separately, such that the remainder of the function under consideration decays fast enough to be amenable to the application of an FFT BDMC2011 ; BDMC2013a ; BDMC2013b . More often, however, it is necessary to waive the advantages of the FFT in favor of a more flexible sampling, specifically adapted to the problem. This, however requires to apply a discrete Fourier transform (DFT) with $\mathcal{O}(N^{2})$ numerical complexity Num_recipes .

A combination of the best of both worlds, i.e. an $N\log_{2}{N}$ scaling on a logarithmic grid, which is able to cover all physically relevant orders of magnitude, has first been proposed by Haines and Jones in form of the logarithmic Fourier transform (LFT), which they have applied in a geophysical context LFT1988 . In its original form however, the LFT is only applicable under very restrictive assumptions on the properties of the function $f(\nu)$ under consideration (e.g. $f(0)=0$) and on the allowed range of the trade-off parameter, which is necessary to adjust the LFT according to the asymptotics of $f(\nu)$.

The aim of this work is to present a generalized version of the logarithmic Fourier-Laplace transformation that in particular applies to functions with nonintegrable power-law tails. We give the corresponding definition in section II and show how the original restrictions can be lifted to extend to generalized functions GelfandBook . Moreover, in section III we give a proof that the LFT converges exponentially fast in the number of grid points used for the numerical evaluation, provided the function satisfies certain analyticity conditions. Furthermore, we discuss how the theorem can be applied for practical purposes and in particular show that functions with algebraic tails are perfectly amenable to the LFT. In section IV, we find an ideal set of the trade-off parameters, based on the asymptotic behavior of the input data and extend the excellent performance of the LFT to convolutions in section V. In section VI we provide several classes of mathematical examples highlighting the advantages of LFTs over FFTs and discuss possible optimizations. Finally, we show in section VII how the LFT can be applied to typical multiscale problems in physics on the example of a density-density correlation function and a simple variant of mode-coupling theory. We conclude in section VIII.

II Definition

II.1 Mathematical Formulation

Following the standard convention in the physics literature, we define the Fourier transform of a function $\hat{f}(t)$ in the time domain as

[TABLE]

while the inverse transform to frequency $\nu$ is given by

[TABLE]

for both $f,\hat{f}\in L^{1}[\mathbb{R},\mathbb{C}]$. In the following, we utilize the LFT to extend the set of argument functions to include certain distributions, the precise properties of which we state below. We introduce the logarithmic frequency and time coordinates $\omega$ and $\tau$ via

[TABLE]

where $\sigma=\pm 1=\eta$ are necessary to distinguish between the positive and negative real axis, while the prefactors $\bar{\nu}$ and $\bar{t}$ are required for dimensional purposes and will be set to unity in the remainder of this paper. With these definitions, the inverse Fourier transform (2) can be written as a convolution for every $t\in\mathbb{R}$:

[TABLE]

where $k\in\mathbb{R}$ denotes the trade-off parameter LFT1988 . By the help of the convolution theorem of Fourier analysis (see also Eq. (39) below), the integral in (4) can be reformulated in terms of the product of two Fourier transforms

[TABLE]

provided that $k$ is chosen such that each of the three Fourier integrals converges, the conditions for which we will detail now.

Since we ultimately aim for a numerical implementation, the LFT can in general only be applied if

[TABLE]

such that the Fourier transformation

[TABLE]

exists in the integral sense of Eq. (1). Regarding the original function $f(\nu)$ this statement is equivalent to

[TABLE]

In the particular case of a power-law behavior, i.e. $f(\nu)\to\nu^{a}$ for $|\nu|\to 0$ and $f(\nu)\to\nu^{b}$ for $|\nu|\to\infty$, the trade-off parameter has to be chosen according to

[TABLE]

As a result, for theses functions the LFT even admits a pole of $f$ located at the origin or a branch cut beginning just there, as well as nonintegrable, algebraically growing asymptotics, provided that they can be controlled by an appropriate value of $k$.

Applying the definition of the $\Gamma$ function the $f$-independent inverse Fourier transform in Eq. (5) can be formally rewritten as

[TABLE]

for $k\in\mathbb{R}\setminus\mathbb{Z}_{0}^{-}$, where the exclusion of nonpositive integers is due to the poles of the Gamma function $\Gamma(k-is)$. We point out that this result has to be considered as the analytic continuation of the integral representation

[TABLE]

that, indeed, only holds if $0<k<1$, as emphasized by Haines and Jones LFT1988 .

Finally, we have to consider the transformation $\mathcal{F}_{s\to\tau}(g_{\sigma}(s)h_{\sigma\eta}(s))$ from the auxiliary variable $s$ to $\tau$ in Eq. (5). Since in any practical implementation the factor $g_{\sigma}(s)$ will only be known in an approximate, discretized form, no analytic continuation can be applied and we have to demand that $g_{\sigma}\cdot h_{\sigma\eta}\in L^{1}$. Given the asymptotics of the product frei05book

[TABLE]

we conclude that $g_{\sigma}\cdot h_{\sigma\eta}\in L^{1}$ requires $g_{\sigma}$ to satisfy $\lim_{|s|\to\infty}|s|^{k+1/2}g_{\sigma}(s)=0$. According to the lemma of Riemann-Lebesgue for differentiable functions koer89book , the latter condition is fulfilled if $F_{\sigma}(\omega)$ is at least

[TABLE]

times differentiable with the derivatives F^{({\color[rgb]{0,0,1}l})}_{\sigma}(\omega)\in L_{1}, for $0\leq l\leq n$. With respect to the original function $f(\nu)$ this implies that $f^{(n)}(\nu)$ exists, while the integrability condition on $F^{(n)}(\omega)$ reduces to Eq. (8), as can be shown by partial integration.

All in all, the logarithmic Fourier transform reads

[TABLE]

which can be applied with a given value of the trade-off parameter $k\in\mathbb{R}\setminus\mathbb{Z}_{0}^{-}$ to all functions $f(\nu)$, that satisfy the summability criterion (8) and are $n$ times differentiable, with $n$ set by Eq. (13).

The computation of the Fourier transformation $f(\sigma|\nu|)=\mathcal{F}(\hat{f}(t))$ via the LFT follows analogously, with the same conditions on $\hat{f}(t)$. It yields the same result as in Eq. (14), yet, with the replacements $f\to\hat{f}$, $\sigma\leftrightarrow\eta$, $\tau\leftrightarrow\omega$ and $(i\sigma\eta)^{is-k}\to(-i\sigma\eta)^{is-k}$, as well as an additional factor of $2\pi$ on the right-hand side of Eq. (14).

We remark that the LFT can readily be generalized to Fourier-Laplace transforms of the form

[TABLE]

with $\phi\in[0,2\pi[$. The result in (14) still holds, merely with the factor $(i\sigma\eta)^{is-k}$ substituted by $(e^{i\phi}\sigma\eta)^{is-k}$. Half-sided transforms, which correspond to the standard definition in case of the Laplace transformation, are simply obtained by restricting the outermost sum of Eq. (14) to $\sigma=1$. Clearly, in most applications Laplace transforms, that is $\phi=\pi$, involve quickly converging integrals. Therefore we will focus on the most critical case of Fourier transforms ($\phi=\pi/2$) and inverse Fourier transforms ($\phi=3\pi/2$), where the transformation kernel entails no exponential suppression of large arguments. Nonetheless, we stress, that even for exponentially decaying integrals the logarithmic transforms are orders of magnitude faster than equidistant grids as is highlighted by the trivial example $f(\nu)=e^{-|\nu|}$ in section VI.

After having discussed the mathematical framework of the LFT let us briefly comment on the role of $k$ (see also Ref. LFT1988 ). The term trade-off parameter refers to the fact that $k>0$ ($k<0$) suppresses both the integrand of the Fourier transformation $\mathcal{F}_{\omega\to s}$ in the definition (14) (cf. also Eq. (8)) for large $\omega\to\infty$ (small $\omega\to-\infty$) and the result in the $\tau\to\infty$ ($\tau\to-\infty$) limit, which corresponds to $t\to\pm\infty$ ($t\to 0$), due to the overall prefactor. Simultaneously, the convergence in the opposite limits is diminished. This dependence on $k$ can be utilized to tune the properties of the LFT to suit the asymptotics of $f$.

A different perspective on the LFT is opened by the interpretation of the trade-off parameter as a shift of the final integration over $s$ to a contour in the complex plane. The discussion of the LFT in terms of contour integrals, which are sensitive to the analytic structure of $f$, is crucial to understand the convergence properties of the LFT detailed in section III

II.2 Numerical implementation

So far, we have only utilized exact analytical reformulations of the problem. However, upon introducing exponential grids with index $\;n\in\{1,2,...,N\}$ in both frequency and time

[TABLE]

and discretizing the auxiliary space via $s_{n}=\Delta s(n+s_{s})$, the usefulness of the form (14) is immediately revealed: The equidistant grids in $\tau$ and $\omega$ make it possible to take advantage of the efficient FFT algorithm – even for Laplace transforms, where fast algorithms otherwise require more elaborate methods from approximation theory Rokhlin1988 ; Strain1992 – while covering low frequencies and short times with a high density of points, as opposed to a reduced sampling density at large arguments. Since in many physical applications the high-energy or frequency range shows an algebraic behavior, this covering of the frequency (momentum) and time (position) domain will be very favorable under many circumstances. Important physical examples include generic correlation functions in frequency and momentum space, while in a critical theory algebraic tails appear in the position and time argument. For instance the momentum distribution $n(k)$ of ultracold Fermions in the vicinity of an open-channel dominated Feshbach resonance chin10Feshbach obeys the Tan energy theorem Tan08energy : For momenta $k$ that exceed any intrinsic inverse length scale $n(k)$ decays like $\mathcal{C}/k^{4}$, where $\mathcal{C}$ is the observable Tan contact density kuhn11contact ; hoin13contact . On the other hand, the phase transition to the superfluid is signaled by an instability of the pair propagator in the low-momentum limit. In this system, the LFT has been applied to study the phase diagram in the presence of a finite spin imbalance Frank2018 . Furthermore, similar challenges arise in the efficient simulation of analog low/high pass filters Christensen1990 , in the context of signal processing as well as in the numerical solution of differential equations Boyd2001 .

In addition to the convenient distribution of points, the grid (16) acquires a high degree of flexibility as the step sizes $\Delta\omega,\Delta\tau$ and $\Delta s$ and the linear shifts $\omega_{s}$ and $\tau_{s}$, that play the role of the prefactors $\bar{t}$ and $\bar{\nu}$, together with $s_{s}$ can be chosen at will. This allows to adjust the method to the asymptotics of various functions as is shown in section VI. The standard choice for all the shifts is $-N/2$ in order to cover positive and negative exponents equally. We will return to the question of how to determine the ideal transformation parameters in section IV.

Before continuing, we remark however, that functions with important features on intermediate scales which cannot be considered as part of the asymptotics of small or large arguments (not even by using the entire set of parameters available in Def. (16)), will yield no advantage over an ordinary FFT. Such cases appear for double-peak structures whose centers are too far apart to be scaled to the high grid-density at $\omega\to 0$ without including an impractically large $N$. Similarly, functions that oscillate uniformly on all scales with a fixed frequency $\bar{\omega}$ will be inevitably undersampled by the given grid at frequencies $\omega\gtrsim\ln(\bar{\omega}/\Delta\omega)$.

III Convergence properties

III.1 Theoretical perspective

Now we address the issue of how efficiently a function $f(\nu)$ that obeys the properties stated below Eq. (14) can be sampled and then Fourier transformed on the exponential grid. This requires to answer the question of how quickly the sum

[TABLE]

representing the numerical, discrete approximation on the grids defined in Eq. (16) converges towards the exact integral (14) as $N\to\infty$. First of all, we note that Eq. (17) indeed approaches the LFT from Eq. (14). To see this one has to consider the limits of the largest values $|\omega_{\pm N}|,|s_{\pm N}|\to\infty$ at vanishing stepsizes $\Delta\omega,\Delta s\to 0.$ Taking the latter limit yields well-defined integrals on finite intervals, since all terms represent measurable functions. In particular, the sum in the second line can be interpreted as Fourier coefficient $F_{l}$ of the periodic function $F_{\sigma}(\omega)$ with period $2N\Delta\omega$. The differentiability of $F_{\sigma}(\omega)$ then implies the asymptotic behavior $F_{l}\lesssim Cs^{(n+1)}$ with a positive constant $C$ koer89book , such that the limit $|\omega_{\pm N}|,|s_{\pm N}|\to\infty$ exists and by its uniqueness we recover the definition of the LFT.

Beyond the mere existence, we now show that under conditions satisfied in many relevant application, the LFT converges exponentially fast in the number of grid points.

Theorem: Let $f(\sigma e^{\omega})$ be a function that is analytic in a closed strip of width $R^{(1)}>0$ around $\bar{\mathbb{R}}$, i.e. the affinely extended real axis of the logarithmic argument $\omega$, and whose asymptotic behavior can be controlled by a suitable choice of the trade-off parameter $k$, such that $F_{\sigma}(\omega)\in\mathcal{S}(\mathbb{R})$, where $\mathcal{S}$ denotes the space of Schwartz functions GelfandBook . Then the deviation of the approximation Eq. (17) from the exact expression Eq. (14) vanishes exponentially in the number of grid points $N$.

Proof: The rate of this convergence will not depend on the exact values of the centers of the grids $\omega_{s}$, $\tau_{s}$ and $s_{s}$. To keep the notation simple, we will in the following assume them to be given by integers.

Obviously, $\mathcal{S}\subset L^{1}$ and the Schwartz functions satisfy the differentiability condition (13) by definition. Furthermore, the integral

[TABLE]

is finite and itself a Schwartz function since the integrand is an element of $\mathcal{S}(\mathbb{R})$. By virtue of the Paley-Wiener theorem titc75book $I^{\sigma}_{1}(s)$ is analytic in a strip around the real $s$-axis, whose width $R^{(2)}>0$ is determined by the asymptotic decrease of $F_{\sigma}(\omega)$, which at least is exponential. Furthermore, the truncation error due to the finite summation interval, scales like $|F_{\sigma}(\omega_{\pm N})|$ and thus in any case merely gives rise to exponential corrections. Therefore, we consider right away the infinite sum. The latter can be replaced by a contour integral around the imaginary axis in the mathematically positive direction, which reads 111This method of rewriting sums in terms of contour integrals is also used to evaluate Matsubara sums in the context of finite temperature quantum field theory abri75 ; altl2010book ; fett71

[TABLE]

Here $n_{B}(z)=1/(\exp(2\pi z/\Delta\omega)-1)$ is the Bose-Einstein distribution with ”inverse temperature” $\beta_{\omega}=2\pi/\Delta\omega$, whose simple poles at $\beta_{\omega}\omega_{n}$ make sure that one recovers the original series with the help of the residue theorem. Subtracting the exact integral $I_{1}^{\sigma}(s)$, which is also taken along the imaginary axis, one obtains for the difference

[TABLE]

where we have made use of Cauchy’s theorem to deform the integration contour such that it remains within boundary of the analytic domain of $f(\sigma\exp(\omega))$ (cf. Fig. 1).

Extracting the dominant exponential behavior we can write for the error

[TABLE]

The function $F_{1}$ arises from the remaining integrals and is the Fourier transformation of an integrable function, which in particular implies that the exponential prefactor indeed yields the leading asymptotic behavior for $|s|\to\infty$ due to the lemma of Riemann and Lebesgue. Moreover, by increasing $\beta_{\omega}\sim N$ the error becomes exponentially small in the number of grid points, provided $|s|<\beta_{\omega}$. From a theoretical perspective, we can take the limit $N\to\infty$ and achieve exponential convergence uniformly in $s$.

Next, one has to investigate the convergence properties of the remaining sum over the auxiliary variable $s$ in (17). First, we focus on the properties of the exact integral

[TABLE]

The asymptotics of the product (12) gives at most rise to an algebraic growth, while for $|s|\to\infty$ the integral $I_{1}^{\sigma}(s)\in\mathcal{S}$ decays exponentially fast, thus rendering $I_{2}^{\sigma}(\tau)$ well-defined. To estimate the error arising from replacing the analytic expression in Eq. (14) by a discretized numerical approximation we consider the difference

[TABLE]

Since the difference $E_{1}^{\sigma}(s)=S_{1}^{\sigma}(s)-I_{1}^{\sigma}(s)$ becomes exponentially small with decreasing $\Delta\omega$, we can replace the sum by the exact integral $I_{1}^{\sigma}(s)$. Furthermore, the sum can be extended to $l\in\mathbb{Z}$ because the truncation to $|l|\leq N$ neglects only terms that are exponentially suppressed due to the asymptotics of $I_{1}^{\sigma}(s)$ that, provided $s_{\pm N}$ are large enough, determines the exponential tails of $S_{1}^{\sigma}(s)$. Then, the error $E_{2}^{\sigma}(\tau)$ can be treated analogously to Eq. (20) in terms of complex contour integrals by introducing the Bose distribution $n_{B}(z)$, which now involves the inverse temperature $\beta_{s}=2\pi/\Delta s$. In this step we have to analytically continue the integral $I_{1}^{\sigma}(s\to z)$, which is in general not known in closed form for genuine complex arguments $z$. Yet, performing this continuation numerically is not required as we need the expression only on a formal level to determine the error. Like above, we shift the contour integrals into the complex plane to the fixed finite real parts $\pm R^{(2)}$ within the width of the analytic domain of $I_{1}^{\sigma}(z)$, but have to take into account the simple poles of the $\Gamma$ function located at the nonpositive integers enclosed by the modified contour. Altogether, we can estimate the error by

[TABLE]

Once again $F_{2}^{\sigma}(\tau,R^{(2)})$ is the Fourier transform of an integrable function and does not overcome the leading exponential such that the first term vanishes exponentially for all $|s|$ in the limit $\beta_{s}\sim 1/\Delta s\sim N\to\infty$. The second term summarizes the contributions from the residues $\text{Res}(\Gamma,-m)=(-1)^{m}/m!$ for $m\in\mathbb{N}_{0}$ and reads

[TABLE]

Note that the Bose functions control the exponential function $\exp((k+m)\tau)$ via $\beta_{s}$ in analogy to the first term in Eq. (24).

In total, we observe that both $E_{1}^{\sigma}$ and $E_{2}^{\sigma}$ decrease exponentially in the limit $\Delta_{\omega},\Delta s\to 0$ for all $s$, while all intermediate steps do not violate this scaling. $\square$

Note that for $\phi\in]\pi/2,3\pi/2[$, which includes the half-sided Laplace transform, convergence is even faster. This, however, takes its toll, when considering the inverse transform, where the exponential growth of $|\Gamma(k-is)|^{-1}\sim e^{\pi|s|/2}$ in the case of the Laplace transform severely magnifies errors and limits the useful interval in $s$ and thereby the attainable precision. For the special case of $k=1/2$, this has been analyzed in detail by Epstein and Schotland Epstein2008 , who, given a noise $\delta$ on the input, also derive a bound on the maximum resolution $\Delta\omega$ and precision $\epsilon$ of the numeric inverse Laplace transform:

[TABLE]

As we will see below, this bound can be significantly improved by the use of the theorem in Sec. III.1 in combination with the knowledge of $R^{(1)}$.

III.2 Practical aspects

Regarding the implementation of the LFT for practical purposes, it is helpful to relate the analytic properties of the function $f$ in the logarithmic argument $\omega$ to the original argument $\nu$, since $f(\nu)$ corresponds to the form that is given in most applications. From this we find bounds for $R^{(1)}$ for some relevant classes of functions. First of all, the analyticity of $f(\sigma e^{\omega})$ implies that $f(\nu)$ is also analytic for $\nu\in\mathbb{R}\setminus\{0\}$. Moreover, the exponential decay of $F_{\sigma}(\omega\to\pm\infty)$ requires that a constant $c_{\sigma}>0$ exists, such that

[TABLE]

To connect the analyticity properties of $f$ with respect to the variables $\nu$ and $\omega$, we first note that the analytic strip of width $R^{(1)}$ in $\omega$ is equivalent to the statement that for each $\omega_{0}\in\bar{\mathbb{R}}$ there is an $R^{\sigma}_{\omega_{0}}>0$, such that the Taylor series of $f(\sigma e^{\omega})$

[TABLE]

converges for all $|\omega-\omega_{0}|\leq R^{\sigma}_{\omega_{0}}$. The width of the strip is given by

[TABLE]

To estimate the minimal size of the analytic domain of $f(\nu)$, we first determine the image of the line $\omega_{0}+i\lambda R^{(1)}$, with $\lambda\in[-1,1]$ centered around the real $\omega_{0}$ under the mapping (3). This yields circular sectors of radius $\nu_{0}=\exp{\omega_{0}}$ that are symmetric around the real $\sigma\nu$ half-axis and centered at $\nu=0$. The half opening angle is given by $\min(\pi/2,R^{(1)})$. The restriction of the angle arises from the separation of the variable $\nu$ with respect to $\sigma=\text{sign}(\Re(\nu))$. Taking the union over all $\omega_{0}\in\mathbb{R}$, which yields the domain of analyticity of $f(\nu)$ gives rise to an infinite cone $|\Im(\nu)|<\tan{(\theta)}|\Re(\nu)|$ with opening angle

[TABLE]

Note that for asymptotically large linear-frequency arguments the required width in $\nu$ space grows linearly, but admits nonanalyticities close to the origin. This is to be expected since nonsmooth variations of $f(\nu)$ that are related to nonanalyticities are very well captured by the exponentially dense grid in the vicinity of the origin. On the other hand, the same behavior at large frequencies will cause deteriorated numerical results due to severe undersampling of sharp features.

In view of these arguments it becomes apparent that algebraic functions with the asymptotic behavior $f(\nu\to 0)\sim\nu^{a}$ and $f(\nu\to\infty)\sim\nu^{b}$ with $a>b$ are ideal candidates for the application of the LFT, since they trivially satisfy condition (27) and the transformation $\nu\to\sigma\exp(\omega)$ removes any nonanalyticity located at the origin, which arises from any $a\in\mathbb{R}\setminus\mathbb{N}_{0}$. We recall that the above considerations do not make any reference to the integrability properties of $f(\nu)$, which in case of algebraic functions can be controlled by the trade-off parameter.

Another class of asymptotic behavior is given by exponential functions, which we discuss with the help of the simple example of a single, dominant exponent $f(\nu)\sim\exp(\alpha(\sigma\nu)^{c})$ for $\sigma\nu\to 0,\infty$, with $c\in\mathbb{R}\setminus\{0\}$ and $\alpha\in\mathbb{C}\setminus\{0\}$ to include oscillating functions. The prefactor is allowed to contain an arbitrary algebraic function to which the problem would be reduced for $\alpha=0$ or $c=0$. After the mapping to $\omega=\Re\omega+i\Im\omega$ we have

[TABLE]

For real $\omega$ the limit $c\,\Re\omega\to\infty$ requires $\Re\alpha<0$, since otherwise the inner Fourier transform $\mathcal{F}_{\omega\to s}$ of the LFT in Eq. (14) is not defined. Note that the trade-off parameter cannot be used to remedy the super-exponential growth for $\Re\alpha>0$, which is also beyond the scope of the notion of generalized Fourier transformations GelfandBook . Furthermore, we observe that even if $\Re\alpha<0$, the boundary of the analytic strip cannot overcome the constraint

[TABLE]

because $f(\sigma e^{\omega})$ diverges exponentially for $\Im\omega$ beyond that value if $c\,\Re\omega\to\infty$. Considering an oscillatory example $f(\nu)\sim e^{i\nu}$, we obtain $R^{(1)}=0$, irrespective of an integrable algebraic prefactor. Therefore, the convergence of the LFT is degraded to an algebraic one on a fundamental level as mentioned at the end of Sec. II.

More generally the width of the analytic strip $R^{(1)}$ can be determined from the positions of singularities of $f(\nu)$ in the complex plane. Modeling $f(\nu)$ in the vicinity of a nonanalyticity at $\hat{\nu}=\nu_{0}\pm iR^{\sigma}_{\nu_{0}}$, where $\nu_{0}\in\mathbb{R}\setminus\{0\}$, $\sigma=\text{sgn}(\nu_{0})$ and $R^{\sigma}_{\nu_{0}}>0$, in the form

[TABLE]

with $\lambda\in\mathbb{C}$ and $\alpha\in\mathbb{R}\setminus\mathbb{N}_{0}$, we find for the $k^{\text{th}}$ derivative

[TABLE]

This form can be used to extract $R_{\nu_{0}}^{\sigma}$ and therefore $R^{(1)}$ and is accessible even for numerical data via finite difference approximations.

Regarding the application of the LFT in practice, we note that it can be implemented with any existing library of FFTs. However, one additionally has to evaluate the Gamma function for complex arguments. Fortunately, these do not depend on $f$ and can be tabulated if several transformations have to be performed. In any case, all modern programming languages include fast algorithms (e.g. Spouge’s approximation spou94 ) to compute $\Gamma$ and the LFT will, therefore, not be severely slowed down compared to an FFT with the same number of data points.

IV Optimal parameter choices

So far, we have shown that the LFT converges exponentially fast towards the exact Fourier transform if the analytic structure of $f$ satisfies the conditions of the above theorem. However, we have not yet made any rigorous statements regarding how many points have to be used to reach the desired precision or how the trade-off parameters should be adjusted to obtain an optimal convergence. As we have seen in section III.1, the differences between the exact Fourier transform and the LFT approximant are well known, such that generic statements about the ideal parameter choices, as well as significant improvements to the results, can be made.

In the following, we will focus mainly on algebraic functions, whose convergence properties can be influenced by the trade-off parameter in contrast to exponential functions. First, we recall that the only restriction on $k$ is given by the convergence of the inner integral in the LFT (14), thus the trade-off parameter has to be chosen according to $1+b<k<1+a$ and $\notin\mathbb{Z}_{0}^{-}$. From a numerical perspective one has to keep in mind that $k$ has to avoid nonpositive integers by a finite margin to bypass the poles of the $\Gamma$ function. In practice, a deviation of 0.01 turns out to be sufficient. Apart from this constraint it is desirable to choose $k$ close to $k_{\text{opt}}=1+(a+b)/2$ that symmetrizes the asymptotic behavior of the integrand in $I^{\sigma}_{1}(s)$ (see Eq. (18)) on both ends of the $\omega$-interval, such that the smallest truncation errors are achieved for the standard value for the centers of the grids (see Eq. (16)). With this trade-off parameter, using $N$ data points, a truncation error no larger than $\epsilon$ requires

[TABLE]

As argued below equation (12) the decay of the integrand in (22) is dominated by $S_{1}^{\sigma}(s)$, the asymptotic behavior of which can be estimated from the contour integral (19), which gives rise to the asymptotic behavior $S_{1}^{\sigma}(s)\sim e^{-R^{(1)}|s|}$. However, approximately at $s_{\text{max}}=\pi/\Delta\omega$ this function drops below the error $E_{1}^{\sigma}(s)$ from Eq. (21). To compute the remaining Fourier transform with truncation errors that are consistent with the previous steps one consequently demands, that $S_{1}^{\sigma}(s_{\text{max}})<\epsilon$ from which one concludes

[TABLE]

Together with (35) this fixes the lower bound of points necessary for an absolute accuracy of roughly $\epsilon e^{-k_{\text{opt}}\tau}=\epsilon e^{-(1+(a+b)/2)\tau}$ to

[TABLE]

Remarkably this scales only logarithmically with the desired precision. We emphasize that this statement holds in general, even if the optimal choice of the trade-off parameter is prohibited. In this case only the prefactor increases. However, one might suspect that $N$ will be drastically increased by the requirement $\Delta s\ll 1$ in order to control the error $E_{2}^{\sigma}(\tau)$, according to equation (24). Since the closest non-analyticity of $S_{1}^{\sigma}(z)$ to the real axis appears at a distance $R^{(2)}\simeq\min\{k-1-b,1+a-k\}$, we infer from the asymptotics $I_{2}^{\sigma}(\tau)\sim\exp(-R^{(2)}\tau)$ and $E_{2}^{\sigma}(\tau)\sim\exp(-(\beta_{s}+|\tau|)R^{(2)})$ that it makes only sense to include values $|\tau_{n}|\leq\tau_{\max}=\pi/\Delta s=\beta_{s}/2$.Therefore, $\Delta s\approx(b-a)\pi/\ln{(\epsilon)}$ is possible.

In addition, the error $E_{2\Gamma}^{\sigma}(\tau)$ due to the proximity between the integration contour and poles of the $\Gamma$ function does not influence $\Delta s$ and apart from a prefactor their $\tau$-dependence is exactly known, as can be seen from Eq. (25). Thus, if the numerical error in the time domain is dominated by these contributions one can fit the residues in $E^{\sigma}_{2\Gamma}$ to the limits $\tau\to\tau_{\pm N}$. In these regimes only numerical noise remains, because by virtue of the lemma by Riemann and Lebesgue the exact function $I_{2}^{\sigma}(\tau)$ has decreased below the desired precision threshold. This procedure works particularly well, since the exponential terms are known exactly (see also the examples in section VI). Moreover, depending on $k$, only a few dominant terms have to be subtracted, while the remaining terms are negligible due to their strictly monotonically decreasing exponents.

Note that this discussion does not include round-off errors. These will give rise to an additional limitation of the attainable precision, since the finite accuracy of the internal numerical operations sets a bound to the possible precision of the LFT and in particular determines how well the decay of $I_{2}^{\sigma}(\tau)$ for $|\tau|\to\infty$ can be resolved. Furthermore, the final multiplication with $e^{-k\tau}$ also affects the error estimate. For negative values of $k$ exceptional precision can be achieved in the regime of small $\tau$. For strongly negative values of $k$, however, these come at the price of enhanced errors at large arguments. As long as round-off errors are ignored, these are a minor problem, since the final Fourier transform in (3) will typically decay to zero at large arguments, which allows to remove the previously discussed systematic errors from $E_{2\Gamma}^{\sigma}$ (cf. section VI). In practice, round-off errors, unfortunately, dominate $I^{\sigma}_{2}(\tau)$ at large arguments, which sets a lower boundary to the useful interval of the trade-off parameter. Positive values of $k$ on the other hand, which are necessary to treat non-integrable functions, that is those with $b>-1$, will result in undesirably enhanced errors near the origin in the image space (here $t$). Removal of these errors will typically involve fitting the asymptotic behavior to the numerical data within the range of $t$-values, where $\hat{f}(t)>\epsilon e^{-k\tau}$ can be satisfied and extrapolating it towards $t\to 0$. The dynamical compression, that is the diminishing length of this $t$-interval as $k$ increases, is the price to pay for numerically Fourier transforming non-integrable functions.

If more than just the minimal number of data points necessary for a given precision are available, one can use them to compensate the enhanced truncation errors and perform several transformations with different trade-off parameters. Larger values of $k$ increase precision for $\tau>0$, while smaller values enhance the accuracy at negative $\tau$. If $k$ can be varied over a wide interval without creating too large truncation errors, this procedure can, for example, be used to mitigate the impact of dynamical compression (see section VI). To improve the results further one can use $\omega_{s}$ to shift the list of $\omega_{n}$ points to optimally sample the asymptotics of the function. For an arbitrary value of $k$ that is compatible with the convergence requirements the condition $F_{\sigma}(\omega_{\pm N})<\epsilon$ translates to

[TABLE]

which reduces to equation (37) and $\omega_{s}=-N/2$ if $k=k_{\text{opt}}$. Similarly, $s_{s}$ can be used to find the best distribution of the auxiliary space points $s_{n}$.

Finally, if a precision close to the round-off limit is required, $e^{-k\tau}$ cannot be orders of magnitude larger than $\hat{f}(\eta e^{\tau})$ but instead should stay as close as possible to $-\ln(\hat{f}(\eta e^{\tau}))/\tau$ in the range of $\tau$ arguments of interest.

In case of the inverse Laplace transform $\hat{f}(t)\to f(\nu)$ we can use the same procedure for the optimization of the transformation parameters on an input with multiplicative noise of amplitude $\delta$. Given a good estimate of $R^{(1)}$, which equals the width of the analytic strip of $\hat{f}(e^{\tau})$ minus $\pi/2$, this allows us to enhance the relative precision of the result near $\nu=1$ to $\delta^{R^{(1)}/(\pi/2+R^{(1)})}$, where, as opposed to Eq. (26), no constraint on $\Delta\omega$ is required. Note that only for $R^{(1)}\ll 1$ one is struck by the dreaded exponential enhancement of errors.

V LFT Convolutions

One of the most important applications of Fourier transforms in theoretical physics relies on the efficient calculation of convolutions. Due to the convolution theorem

[TABLE]

the computational complexity $\propto N^{2}$ of the direct discretized evaluation of the integral can be reduced to $\propto N\log_{2}(N)$ when using the FFT algorithm. There are, however, many situations, where the convolution theorem may not be utilized. For example, if one of the factors $f(\nu)$ or $g(\nu)$ decays too slowly to be integrable (that is no faster than $1/\nu$), its Fourier transform can no longer be understood as an integral and the identity (39) cannot be used to improve performance, since the FFT will fail to correctly determine either $\hat{f}(t)$ or $\hat{g}(t)$. Despite these complications, the convolution may still be defined as an ordinary integral (at least as long as the product $f(\nu)g(\nu)$ decays faster than $1/\nu$) and numerical evaluation is cumbersome but straightforward.

Here the LFT really excels. On the one hand, it can be evaluated much faster than any direct evaluation of the convolution (even if performed on an optimized grid) due to its superior scaling in the number of data points. On the other, it is able to cope with non-integrable functions as long as a suitable trade-off parameter exists. In other words, the LFT is much less plagued by convergence problems than the FFT. Slowly decaying functions, as we have already pointed out in the last section, can be numerically transformed with the LFT at the price of dynamical compression that inevitably reduces the signal-to-noise contrast at small arguments in the image space. However, for convolutions this problem is slightly less pronounced: If possible, setting the trade-off parameter of the final $\tau\to\omega$ back-transform in (39) to $k_{\text{back}}=1-k_{1}-k_{2}$, where $k_{1}$ and $k_{2}$ are the optimized parameters for the transforms of $f(\eta e^{\omega})$ and $g(\eta e^{\omega})$, respectively, renders the result unaffected by any dynamical compression in $\tau$ because the noise level in the time domain stays constant at $\mathcal{O}{(\epsilon)}$, as the problematic exponential prefactors cancel. In the exotic case of strongly divergent integrals, this value of $k_{\text{back}}$ may not be useful if $e^{k_{\text{back}}\omega_{\text{max}}}$ becomes much larger than the expected result, in which case more data points and the less aggressive, symmetric choice $k_{1}=k_{2}=k_{\text{back}}$ are typically better suited (see last example in VI).

VI Examples and optimizations

To benchmark the LFT and to illustrate the role of the transformation parameters, in particular of $k$, we compute the Fourier transforms for several examples and compare the results to the exact solutions. Without loss of generality, we focus on functions that are centered around the origin and that vary on a characteristic scale of unity. Deviations from that behavior can be remedied by preprocessing the function with a variable transformation which combines a shift of the original argument followed by rescaling it with a proper $\bar{\nu}$.

Let us begin with a benign example, a Lorentzian curve

[TABLE]

which could also be transformed with an ordinary fast Fourier transform. However, the slow convergence of the integral implies that reaching a global precision of $\epsilon=10^{-12}$ with the FFT requires roughly $10^{13}$ data points, which exceeds numerical feasibility by several orders of magnitude. In contrast, to achieve the same accuracy with the LFT only a little more than 300 points suffice, as can be deduced from Eq. (37) with $R^{(1)}=\pi/2$. Indeed, Fig. 2 has been obtained with 360 points and the (quasi-)optimal parameter $k=-1/100$, since $k_{\text{opt}}=0$ is prohibited by the $\Gamma$ function. Setting up the LFT in this way, the precision is no longer limited by the finite resolution, but by double-precision floating point arithmetic and error-propagation therein. In addition, the interval in $t$ can be chosen arbitrarily by adjusting $\Delta\tau$ and $\tau_{s}$ with no influence on the error level. As discussed in Sec. IV, due to the proximity to the pole of the $\Gamma$ function the last point $\hat{f}(t_{\text{max}})$ has to be subtracted which corresponds to the constant that arises from the leading $m=0$ contribution to $\exp(-k\tau)E_{2\Gamma}^{\sigma}(\tau)$, see Eq. (25). The same procedure, which amounts to nothing else than a trivial subtraction of a one-parameter fit has been used for all other plots (except Fig. 6) as well.

A significantly more demanding example (on the branch where $\sqrt{-1}=i$) is given by the function

[TABLE]

whose Fourier transform has to be understood in the sense of tempered distributions. For positive arguments $t>0$ it reads

[TABLE]

According to section IV, the optimized trade-off parameter is close to $k=1$, which removes the divergent behavior from the numerical integrals at the price of reduced precision at very small values of $t$. To demonstrate how the ideal choice of $k$ might depend on the data range of interest in the image space, Fig. 3 depicts the result for the ideal trade-off parameter $k=1.01$ and the suboptimal $k=0.71$. In order to achieve errors of $10^{-12}$ at $\tau=0$ in both cases, which requires roughly $N=600$ in the optimized setting, the grid size has been increased to $N=1000$. In agreement with the general discussion, values of $k$ larger than $1$ reduce errors at large arguments, while those smaller than unity increase precision close to $t=0$.

The next example shows

[TABLE]

which transforms into

[TABLE]

and is again correctly described by the LFT on only 560 points (see Fig. 4). However, the divergence of $f(\nu)$ as $\nu\to\infty$ results in an even stronger dynamical suppression than before.

The well-behaved function $f(\nu)=e^{-|\nu|}$, which in fact corresponds to the inverse transformation of the very first example (40), could also be treated by an FFT. However, the same accuracy as demonstrated in Fig. 5 with $N=480$ on the logarithmic grid would require more than a million data points on a linear grid, illustrating that even for the most benign functions the LFT can outperform the direct application of an FFT.

The convolution of

[TABLE]

with itself, that appears frequently in the evaluation of Feynman diagrams with non-relativistic propagators abri75 ; fett71 , cannot be treated by FFTs, as the integral over $f(\nu)$ does not exist. Using the LFT remedies this issue by the help of the trade-off parameter. For instance, using the optimized value of $k$ given in section IV, a constant error of roughly $10^{-12}$, which is limited only by the internal floating point precision, is obtained with only $N=560$ points. In Fig. 6 the two leading contributions to $E^{\sigma}_{2\Gamma}$ from Eq. (25) for the transformation from $t$ to $\nu$ with $m=0,1$ were subtracted by fitting the two corresponding parameters $I_{1}^{\sigma}(-i(k+m))$ to the high-frequency range near $\nu=10^{28}$.

Encouraged by these results, one can try to convolve some more exotic functions, for example $f(\nu)=\ln(\nu^{2}+1)$ with itself. In this case neither the convolution as an integral, nor the product of the distributions in the time domain is in general well-defined GelfandBook . Nevertheless, employing a cutoff $e^{-\delta|\nu|}$, with $\delta>0$, in the logarithmic frequency space and sending $\delta$ to zero at the end one finds analytically

[TABLE]

As Fig. 7 demonstrates, this result is again very accurately recovered by means of an LFT with $N=560$. Here, due to the large frequency interval used, round-off errors multiplied by $e^{-k_{\text{back}}\omega}$ are the limiting factor at large $\nu$. This affects the choice of trade-off parameters for the forward and backward LFTs, which is expected on general grounds, as remarked at the end of section V on convolutions. Furthermore, the first sub-leading divergence $\propto e^{2\omega}$ due to the next pole $(m=2)$ of the $\Gamma$ function beyond the constant has been fitted with a single parameter against the raw result at $\omega_{N}$ and subtracted.

VII Application to physical examples

Following the purely mathematical discussion, we now provide physical examples to highlight the real-world advantages of the LFT. These demonstrations are deliberately chosen to be simple, yet of relevance to current research and with apparent generalizations to more challenging problems.

VII.1 Polarization function

The polarization function of the one-dimensional Bose gas is given by Mahan_book ; Kamenev_book

[TABLE]

with $n_{B}(k)=1/(e^{\beta\xi_{k}}-1)$ the Bose-Einstein distribution for the inverse temperature $\beta=1/(k_{B}T)=1$ and the free dispersion $\xi_{k}=k^{2}/(2m)-\mu$ with momentum $k$, mass $m=1/2$ and chemical potential $\mu$. $\Pi(\omega,q)$ describes density-density correlations and in general has no known closed expression. One therefore has to rely on a numerical evaluation of the integral. In case of the fugacity $z=\exp{(\beta\mu)}$ approaching unity from below, the density fluctuations in the regime of long wave lengths proliferate, which is reflected in the singular behavior $n_{B}(0)\sim-1/(\beta\mu)$. As a consequence, the direct evaluation of the polarization function becomes numerically expensive. However, precisely these low temperature correlation functions are a common ingredient in quantum many-body theories, in particular: quantum critical transport sach11 , response near phase transitions in ultracold atoms endr12 and Bose gases in optical cavities rits13 .

In the following, we will show that an efficient evaluation of $\Pi(\omega,q)$ with unrivaled precision is possible by means of the LFT. The integral in Eq. (47) can be rewritten as a convolution, which allows for an efficient treatment that requires only two one-dimensional (half-sided) Fourier transforms:

[TABLE]

The computation of the polarization function on a two-dimensional $(\omega,q)$-grid of size $N\times N$ therefore requires only $\mathcal{O}(N^{2})$ operations for the evaluation of $\sigma(y(\omega,q))$ in contrast to $\mathcal{O}(N^{3})$ for the direct approach. In fact, the actual fast Fourier transform in Eq. (48b) results only in subleading corrections to the overall complexity. The main advantage of the LFT over other Fourier transforms, however, lies in its accuracy. We highlight this in Fig. 8 by comparing the absolute error obtained for an ordinary FFT, the LFT and a discrete Fourier transform (DFT) in combination with a cubic spline interpolation. All algorithms use the same number of data points as well as optimized transformation parameters, with the DFT and LFT operating on a common grid. For all relevant values the LFT outperforms the other methods by several orders of magnitude, furthermore, in contrast to the FFT a much larger interval can be sampled.

VII.2 Glass transition

The glassy, mechanically rigid state of amorphous materials and its realization by supercooling liquids has been studied for a very long time debe01 . Nevertheless, a theoretical description of this state is difficult since one has to deal with density fluctuations on various length scales and very slow relaxation processes, as observed in experiments stef94 . One theoretical approach to this problem is mode coupling theory reic05 ; goet08 which is based on an effective equation of motion for the dynamical structure factor $S(\mathbf{k},t)$. In the following, we illustrate how the LFT, which by construction is capable of dealing with multi-scale problems, can be applied to these kinds of models. Here we focus on a simplified variant of mode-coupling theory due to Leutheusser leut84 and Bengtzelius et al. beng84 and remark on the advantageous properties of the LFT for more generic problems of this kind.

To study the relaxation of density distortions one introduces an effective temporal order-parameter $\Phi(t)$, whose phenomenological time evolution is given by leut84

[TABLE]

The left side of the equation is a simple harmonic oscillator with damping rate $\gamma$ and frequency $\Omega_{0}$. The correlations responsible for the glass transition are incorporated in the memory integral on the right-hand side and weighted by the dimensionless, positive coupling constant $\lambda$. The initial conditions $\Phi(0)=1$ and $\dot{\Phi}(0)=0$ model the original deviation from the equilibrium state $\Phi\equiv 0$. Despite its simplicity, the above equation takes the basic properties of the glass transition into account which we briefly review before presenting the solution based on the LFT.

The physical order parameter distinguishes between two phases via the long-time limit $\Phi(t\to\infty)=c$: The ergodic phase is characterized by perfect relaxation corresponding to $c=0$, in contrast to the glass phase where the initial distortion never disappears completely and thus $c>0$ . To gain further insight into the phase diagram we apply the half-sided Fourier transformation

[TABLE]

Furthermore, we reparametrize the function $\Phi(t)=\delta\Phi(t)+c$ to obtain the asymptotics

[TABLE]

irrespective of the phase. The constant value $c$ gives rise to a term proportional to $\delta(\nu)$ in the Fourier transform of Eq. (49) that has to be canceled to satisfy $\delta\Phi(t\to\infty)=0$. This is the case if

[TABLE]

which not only determines the asymptotic value of the order parameter but also identifies the critical coupling for the glass transition $\lambda_{c}=1$. As has been shown in Refs. leut84 ; beng84 by analytic means, the approach to the phase boundary is characterized by a divergent low-frequency limit of $\delta\hat{\Phi}(\nu)$ that follows the power law

[TABLE]

In the ergodic phase the half-sided Fourier transform of Eq. (49) yields the self-consistent relation

[TABLE]

Similarly, in the glass phase one obtains a quadratic equation with the solution

[TABLE]

where in the first line one has to choose the branch that yields a positive $\Re\hat{\Phi}(\nu)$, since this function represents a retarded, bosonic correlation function fett71 . The equations (54) or (55) can be solved in an iterative manner, which requires repeated Fourier transformations between the time and frequency spaces. To reliably compute the critical behavior in the vicinity of the glass transition, however, one has to include very long times, especially to capture the extremely slow relaxation when approaching the instability from the ergodic phase. Yet, this is exactly the scenario the LFT has been devised for.

Note that the formulation of Eq. (49) in terms of $\delta\Phi$ is not only helpful for analyzing the problem in further detail, but in order to apply the LFT it is also mandatory because $\delta\Phi$ satisfies the condition (9). Figure 9 shows $\Phi(t)$ at the values $\lambda=1\pm 2^{-21}\approx 1\pm 4.7\cdot 10^{-7}$ in the immediate vicinity of the glass transition. In the ergodic phase the plateau, which characterizes the so-called regime of $\beta$-relaxation, reaches times of $\mathcal{O}(10^{12})$ before the final $\alpha$-relaxation to $\Phi=0$ sets in. In addition to the global features of the dynamics, the LFT reproduces the critical exponents from Eq. (53) with a numerical error on the order of $10^{-4}$ (see inset of Fig. 9). To achieve such small errors one can profoundly benefit from the flexibility of the LFT: While the LFTs have to operate on identical $\omega$ and $\tau$ grids, irrespective of the direction of the transformation, the auxiliary $s$ space can be sampled for the two directions $\omega\to\tau$ and $\tau\to\omega$ independently. Moreover, one can introduce two LFTs for $\Re\delta\Phi(\nu)$ and $\Im\delta\Phi(\nu)$ separately without altering the overall computational cost. In total, the numerical effort to obtain $\Phi(t)$ at all times for a given $\lambda$ scales like $N_{\text{it}}\cdot N\log N$ where $N_{\text{it}}$ denotes the number of iterations needed to reach convergence. Using optimized trade-off parameters and a grid of length $N=2400$ suffices to produce the results shown in Fig. 9. The small number of data points used in the LFT allows to find the converged order parameter in a couple of seconds, even close to the phase transition, where $N_{\text{it}}\sim 10^{4}$.

A direct numerical solution of the integro-differential equation (49) on a discretized time axis requires step sizes $\Delta t\lesssim 1$ independent of the magnitude of $t$ in order to compute the cancellations between the various terms with sufficient precision. These are responsible for the slow evolution and the formation of the plateau over six orders of magnitude in time. Larger $\Delta t$ would lead to an instability of the numerical solution that diverges away from the physical $\Phi(t)$. Due to the scaling $\mathcal{O}(N_{\text{dir}}^{2})$ of a direct approach and limited step size it is numerically completely unfeasible to reach times of order $10^{14}$. This, however, is necessary to reliably determine $\Phi(\omega=0)$ and consequently the critical exponents 222It is worth to mention, that in combination with a partially analytical solution, step sizes $\Delta t\gg 1$ are possible. The implementation of such a method is, however, far more complicated and time-consuming than the simple, direct approach presented here Goetze1988 . .

Regarding more complicated versions of mode-coupling theory that resolve the dependence on the length scales, thereby considering the structure factor $S(\mathbf{k},t)$ instead of $\Phi(t)$, the grid size of $N\sim 10^{3}$ used here is still small enough to incorporate a second argument without running into memory limitations. If the coupling between different wave vectors can be written in terms of convolutions, which usually is the case reic05 ; goet08 , the LFT can also be applied to simplify the spatial dependence. Very similar problems appear in the context of approximate equations of motions of correlation functions in quantum field theory, which typically include algebraic decays in frequency and momentum space. As mentioned earlier, an example of the application of the LFT in the context of ultracold Fermi gases can be found in Ref. Frank2018 .

VIII Conlcusion

We have shown rigorously that the LFT can be used to numerically transform nonintegrable functions, as long as their asymptotics can be controlled by the trade-off parameter. Furthermore, we have proven that one can achieve exponential convergence in the number of data points if the function is analytic in a cone with finite opening angle around the real axis in the original argument. Finally, we have given several examples that benchmark the superior convergence of the LFT compared to the FFT including functions that have to be considered within the concept of generalized Fourier transformations.

Acknowledgments The authors thank Wilhelm Zwerger for fruitful discussions and comments on the manuscript and Christian Johansen for pointing out erroneous factors of $2\pi$. This work has been supported by the Nanosystems Initiative Munich (NIM).