Double exponential quadrature for fractional diffusion

We introduce a novel discretization technique for both elliptic and parabolic fractional diffusion problems based on double exponential quadrature formulas and the Riesz–Dunford functional calculus. Compared to related schemes, the new method provides faster convergence with fewer parameters that need to be adjusted to the problem. The scheme takes advantage of any additional smoothness in the problem without requiring a-priori knowledge to tune parameters appropriately. We prove rigorous convergence results for both, the case of finite regularity data as well as for data in certain Gevrey-type classes. We confirm our findings with numerical tests.


Introduction
The study of processes governed by fractional linear operators has gathered significant interest over the last few years [8,22,35] with applications ranging from physics [1] to image processing [1,15,16], inverse problems [19] and more. See [33] for an overview of applications in different fields. The goal is to solve problems of the form with parameters β,α ∈ (0, 1]. There are multiple (non-equivalent) ways of defining fractional powers of operators. We mention the integral fractional Laplacian and the spectral definition [22]. In this paper, we focus on the spectral definition which is equivalent to the functional calculus definition.
For discretization of such problems, both stationary and time dependent, multiple approaches have been presented. A summary of the most common can be found in [2,22]. They can be broadly distinguished into three categories. The first class of methods uses the Caffarelli-Silvestre extension to reformulate the problem as a PDE posed in one additional spatial dimension. This problem is then treated by standard finite element techniques [6,24,25,[27][28][29]. The second big class of discretization schemes, and the one our new scheme is part of, was first introduced in [7] and later extended to more general operators [5] and time dependent problems [3,4,26]. They are based on the Riesz-Dunford calculus (sometimes also referred to as Dunford-Taylor or Riesz-Taylor) and employ a sinc quadrature scheme to discretize the appearing contour integral. sinc quadrature, and overall sinc-based numerical methods are less well known than their polynomial based counterparts, but provide rapidly converging schemes [21,32] with very easy implementation. The quadrature relies on appropriate coordinate transforms in order to yield analytic, rapidly decaying integrands over the real line and then discretization using the trapezoidal quadrature rule. In [34] it was realized that by adding an additional sinh-transformation, it is possible to get an even faster convergence for certain integrals. Namely, writing N q for the number of quadrature points, instead of convergence of the form e − √ N q , it is possible to get rid of the square root and obtain rates of the form e − Nq ln Nq . Further developments in this direction are summarized in [23]. Such schemes are commonly referred to as double exponential quadrature or sinh-tanh quadrature. Thirdly there is the large class of methods based on rational approximation of the functions z −β and the Mittag Leffler-Function e γ,μ (z) (see (3.18) for the precise definition). As shown in [17], this class also encompasses the previous two approaches while also allowing some other methods, based on general rational approximation algorithms like Best-Uniform-Rational approximation (BURA) or the "Adaptive Antoulas-Anderson"-algorithm (AAA) from [30]. Finally, there exist some further methods based on reduced basis and rational Krylov methods [9,10,12,13] which are strongly related to rational approximation.
In this paper we investigate whether the discretization of the Riesz-Dunford integral can benefit from using a double exponential quadrature scheme instead of the more established sinc-quadrature. We present a scheme that retains all the advantages of [3][4][5] while delivering improved convergence rates. Namely, the scheme is very easy to implement if a solver for elliptic finite element problems is available. It is almost trivially parallelizable, as the main cost consists of solving a series of independent elliptic problems. In addition, it provides (compared to sinc-methods) superior accuracy over a wide range of applications and does not require subtle tweaking of parameters in order to get good performance. Instead it will automatically pick up any additional smoothness of the underlying problem to give improved convergence. Since for each quadrature point an elliptic FEM problem needs to be solved, reducing the number of quadrature points greatly increases performance of the overall method.
Compared to the BURA and AAA rational approximation methods, the sinc-and double-exponential quadrature based algorithms have several advantages. Firstly, the implementation is very simple with quadrature nodes that are known explicitly. The quadrature points are also independent of the spectrum of the operator L and no explicit bound on the largest eigenvalue is required. This makes them better suited for highly accurate but highly ill-conditioned discretizations like the hp-FEM scheme in [26].
Secondly, the quadrature points are independent of the function that is to be approximated. Most notably, when considering the time-dependent problem with inhomogeneus right-hand side in Sect. 3.3, all the linear systems that need to be solved are independent of the time t or the integration variable τ . This makes the full timedependent problem of the same cost (with respect to the number of systems that need solving) as the simple stationary problem. Thirdly, the quadrature based methods allow for very detailed analysis, as showcased in this article. In addition to the quadrature analysis, they also allow for detailed analysis of the error brought in by a discretization in space [26]. In practice, we also observed better numerical stability in the presence of rounding errors, as showcased in Fig. 4.
The paper is structured as follows. After fixing the model problem and notation in Sects. 1.1, 3 introduces the double exponential formulas in an abstract way and we collect some known properties. In addition, we provide one small convergence result which, to our knowledge, has not yet appeared in the literature; we show that the double exponential formulas at least provide comparable convergence of order e − √ N q even without requiring additional analyticity compared to standard sinc methods. The paper is structured as follows. In Sect. 1, we introduce the general setting and the functional calculus. Sect. 2 introduces the quadrature scheme as well as the model problems we are interested in. We also state the main convergence results. Sect. 3 is devoted to proving these results. Sect. 3.1 presents the abstract analysis for sinc methods and collects some known properties. In addition, we provide one small convergence result which, to our knowledge, has not yet appeared in the literature; we show that the double exponential formulas at least provide comparable convergence of order e − √ N q / ln(N q ) even without requiring additional analyticity compared to standard sinc methods. In Sect. 3.2, we look at the case of a purely elliptic problem without time dependence. It will showcase the techniques used and provide the building block for the more involved problems later on. In Sect. 3.3, we then consider what happens if we move into the time-dependent regime. Section 4 provides extensive numerical evidence supporting the theory. We also compare our new method to the standard sinc-based methods. Finally, Appendix A collects some properties of the coordinate transform involved. The proofs and calculations are elementary but somewhat lengthy and thus have been relegated to the appendix in order to not impact readability of the article. Throughout this work we will encounter two types of error terms. For those of the form e − γ k we will be content with not working out the constants γ explicitly. For the more important terms of the form e − γ √ k we will derive explicit constants γ which prove sharp in several examples of Sect. 4.
We close with a remark on notation. Throughout this text, we write A B to mean that there exists a constant C > 0, which is independent of the main quantities of interest like number of quadrature points N q or step size k such that A ≤ C B. The detailed dependencies of C are specified in the context. We write A ∼ B to mean A B and B A.

General setting and notation
In this paper, we consider problems of applying holomorphic functions f to selfadjoint operators, for example the Laplacian. The two large classes of problems treated in this paper stem from the study of fractional diffusion problems, both in the stationary as well as in the transient version. Since it does not incur additional difficulty compared to the explicit setting of Remark 1.2, we will work in the following abstract setting: Assumption 1.1 Let X be a Hilbert space and L be a positive definite, self adjoint operator on X such that there exists a sequence of eigenvalues λ j > 0 with associated eigenfunctions ϕ j ∈ X , j ∈ N 0 , such that (ϕ j ) ∞ j=0 is an orthonormal basis of X .
Given the eigenvalues and eigenfunctions of L, we define the spaces for β ≥ 0

Remark 1.2
The problem we have in mind for our applications is the following: given a bounded Lipschitz domain , we consider the space X := L 2 ( ) and the self adjoint operator where A ∈ L ∞ ( ; R d×d ) is uniformly symmetric and positive definite and c ∈ L ∞ ( ) satisfies c ≥ 0 almost everywhere. The domain of dom(L) is always taken to include homogeneous Dirichlet boundary conditions. In this case, the spaces H β (X ) correspond to the standard (fractional) Sobolev spaces often denoted by H β ( ) or H β ( ) in the literature. [5] considers an even more general class of operators, namely the class of "regular accretive operators". We expect some of the results of this article to carry over also to such a class, but since many of our proofs rely on the decomposition using real eigenvalues, such generalizations would be non-trivial.

Remark 1.3
The spaces H β are the natural setting for our regularity assumptions on the data. If we are interested in convergence beyond root-exponential rates, we need the following class of functions of Gevrey-type.
Compared to the standard Gevrey-class of functions, these spaces also include boundary conditions for the functions L n f for all n ∈ N. If the boundary conditions are met, we can then estimate Examples for such functions are those only containing a finite number of frequencies when decomposed into the eigenbasis of L, but also more complex functions such as smooth bump functions with compact support are admissible (see [31,Section 1.4]).
One natural way of defining a functional calculus for the operator L is based on the spectral decomposition.
An alternative definition for holomorphic functions, which will prove more useful for approximation is given in the following Definition. For simplicity, we restrict our considerations to decaying functions g. In this case, it can be shown (see also [3,Section 2]) that the operators resulting from Definitions 1.4 and 1.5 coincide.
where the integral is taken in the sense of Riemann, and C is the smooth path The parameter κ > 0 is taken sufficiently small such that κ < λ 0 , where λ 0 is the smallest eigenvalue of L. The parameters σ and θ can be used to tweak the discretization. We have observed the best behavior for σ := 1/2 and θ := 4; cf. Sect. 4.

Remark 1.6
The choice of path in Definition 1.5 is somewhat arbitrary. It is only required to encircle the spectrum of L with winding number 1. Throughout this paper, we will only ever use the same path and thus make it part of our definition.

Remark 1.7
One could also think to allow σ ∈ (0, 1). For the practical application of the scheme this does not make a big difference, but the analysis for σ = 1 in this paper makes heavy use of the half-angle formula. Therefore we restrict our view to the cases σ = 1 or σ = 1/2. In numerical experiments, methods with σ = 1/2 work, but we decided that the small difference in performance does not warrant the much greater complexity of analysis.

Model problems, discretization and results
In this section, we introduce the discretization methods and, in order to ease the reading of the article, we present the most important of the convergence results. All of the sometimes very technical proofs are relegated to Sect. 3. The main role in our discretization schemes will be played by the following coordinate transform which parametrizes the contour in Definition 1.5: We will focus on the cases σ ∈ 1 2 , 1 and θ ≥ 1. κ is again taken sufficiently small as in Definition 1.5.
Using this transformation, we can introduce the double exponential quadrature approximation of the Riesz-Dunford calculus in Definition 1.5. Because the discretization by quadrature will appear repeatedly for different functions and operators, we introduce the following notation: and Q L (g) := Q L (g, ∞) for the case where no cutoff is performed. The quadrature error will be denoted by where g(L) is given via the Riesz-Dunford integral 1.5. Again, we write E L (g) := E L (g, ∞).

Remark 2.2
In Definition 2.1, we will often work with the special case L = λ. This is taken to mean the scalar multiplication operator u → λu on the vector space X .
We apply the function to the following problems: For both model problems, we prove two convergence results, depending on the regularity of the data. In the case of "finite regularity", the data ( f or u 0 ) are assumed to be in a space H 2ρ for some ρ > 0. This results in bounds of root-exponential order N q / ln(N q ).
The second case is the one were the data are in the Gevrey-type classes G L introduced in (1.2). For such functions, the double-exponential discretization leads to an improved convergence of the form O(e − γ Nq ln(Nq ) ).

The elliptic problem
As our first model problem, we consider the following elliptic fractional diffusion problem: Using the Riesz-Dunford formula, this is equivalent to computing In order to get a discrete scheme, we replace the integral with the quadrature formula. Given N q ∈ N and k > 0, the approximation to (2.4) is then given by

Remark 2.3
Since in practice, the solution operator (L − z) −1 is not computable, one would in addition replace (L − z) −1 by a Galerkin solver in order to obtain a fully computable scheme. In the Setting of Remark 1.2, this means the following: given a closed subspace Given discretization parameters V h ⊆ H 1 ( ), N q ∈ N and k > 0, the fully discrete approximation to (2.4) is then given by In order to keep presentation to a reasonable length, we focus on the spatially continuous setting. We only remark that discretization in space can be easily incorporated into the analysis. For low order finite elements one can follow [3]; for an exponentially convergent hp-FEM scheme we refer to [26].
Remark 2. 4 We should point out that for the elliptic problem, there exist methods based on the Balakrishnan formula (see also Sect. 4) which do not require complex arithmetic. On the other hand, since we are only approximating real valued functions, we can exploit the symmetry of (2.2) to only solve for j ≥ 0, thus halving the number of linear systems. This results in (roughly) comparable computational effort for both the Balakrishnan and the double exponential schemes. Due to their better convergence the DE-schemes might therefore still be advantageous.
The convergence of the new method can be summarized in the following two theorems.
Theorem 2.5 Let u be the exact solution to (2.4) and assume f ∈ H 2ρ ( ) for some ρ ≥ 0. Let β ≥ β with β ∈ (0, 1] and u k := Q L (z −β , N q ) f denote the approximation computed using stepsize k > 0 and N q ∈ N quadrature points. Then, the following estimate holds for all ε ≥ 0 and r ∈ [0, β/2]: where the rate p(σ, θ ) is given by For ε > 0, the implied constant and γ may depend on ε, r , the smallest eigenvalue λ 0 of L, β, κ, θ and σ . But they are independent of ρ, β, k, and f . If ε = 0, the constants may in addition depend on ρ and β.

Remark 2.6
When comparing Theorem 2.5 to the estimates of the standard sincquadrature one might think that the double exponential method is inferior due to the √ k vs k behavior. This misconception can be cleared up by considering the better decay properties of the double-exponential formula. It allows to choose k ∼ ln(N q )/N q compared to the standard sinc-quadrature choice of k ∼ N −1/2 q without the cutoff error becoming dominant. Using this choice, the exponential term scales like N q / ln(N q ) for double exponential and N q for standard sinc respectively. As is shown in Sect. 4, the better constants in the exponential still often outweigh the presence of the ln-term for the double-exponential quadrature.

Remark 2.7
For most of the computation, the convergence rate is determined by the factor p(σ, θ ) in Corollary 3.8. We observe that for θ = 1, picking σ = 1/2 roughly doubles the convergence rate. Similarly, it often appears beneficial to pick larger values of θ . Especially for σ = 1, we get an asymptotic rate for θ → ∞, which is the same as in the case of σ = 1/2. But we need to point out that increasing θ means that we have to decrease the value d(θ ), which determines the rate in the higher orders terms of the form e −γ /k , thus leading to those terms dominating in a larger and larger preasymptotic regime. Overall, the method using σ = 1/2 and setting θ moderately large is expected to give the best convergence rates; cf. Sect. 4.
The previous theorem shows that in general, the convergence behaves like O(e − γ √ k ). It also shows that, if the function f in the right-hand side has some additional smoothness, the method automatically detects this and delivers an improved convergence rate. If the additional smoothness is in the right Gevrey type classes, we can establish convergence which is beyond the root exponential behavior. The details can be found in the following theorem: Theorem 2.8 Let u be the exact solution to (2.4) and assume that there exist constants Assume that β > β with β ∈ (0, 1]. Let u k := Q L (z −β , N q ) f denote the approximation computed using stepsize k ∈ (0, 1/2) and N q ∈ N quadrature points. Then, the following estimate holds: The implied constant and γ may depend on ω, the smallest eigenvalue λ 0 of L, κ, θ , σ , R f , β, and ω. If ω = 0, the logarithmic term may be removed.

The parabolic problem
The second model problem we consider is a time-dependent fractional diffusion problem of parabolic type. We fix α, β ∈ (0, 1] and a final time T > 0. Given an initial condition u 0 ∈ X and right-hand side f ∈ C([0, T ], X ) we seek u : [0, T ] → dom(L β ) satisfying where ∂ α t denotes the Caputo fractional derivative. Following [4], the solution u can be written using the Mittag-Leffler function e α,μ (see (3.18)) as (2.8) Here we again use either the spectral or, equivalently, the Riesz-Dunford calculus to define the operators. We discretize this problem by using our double exponential formula. Namely for k > 0 and using N q ∈ N quadrature points, (2.9)

Remark 2.9
In practice, in order to get a fully computable discrete scheme, one would again replace the resolvent by a Galerkin solver and the convolution in time by an appropriate numerical quadrature. For example, [4] presents a low order approximation scheme. In order to retain exponential convergence, [26] uses a scheme based on hp-FEM and hp-quadrature. We summarize the construction briefly. For a given degree p ∈ N 0 , and interval I , we denote the Gauss quadrature points and weights on (−1, 1) by (x [11,Section 2.7] for details. We then consider a geometric mesh on (0, 1) with grading factor σ ∈ (0, 1) and parameter L ∈ N, L ≤ p given by On each one of these elements, we apply a Gauss quadrature, reducing the order as we approach the singularity, i.e., we get the nodes and weights as The convolution in (2.7) is then replaced by In order to get a fully discrete scheme, this function is then discretized using the double exponential quadrature scheme: In order to not overwhelm the presentation of the paper, we do not consider these types of discretization errors. The analysis of such errors could be taken almost verbatim from the references [3,26].
The analysis of the method again comes in the form of two theorems, one for the case of finite regularity and one for regularity in the Gevrey-type classes G L (C f , R f , ω).

Error analysis
In this section, we analyze the quadrature error when applying a double exponential formula for discretizing certain integrals. For θ ≥ 1, δ > 0 we define the sets where for each θ , d(θ ) is a constant which is assumed sufficiently small in order for Lemmas A.3, A.4, and A.8 to hold. Since all the proofs analyzing the properties of ψ σ,θ are elementary but somewhat lengthy and cumbersome, they have been relegated to Appendix A. The most important properties are, that y → ψ σ,θ (y) for y ∈ R traces the contour in the definition of the Riesz Dunford calculus (see Definition 1.5), and that it is analytic in D d (θ) . The other important results concern the points where ψ σ,θ crosses the real axis, as these points correspond to (possible) poles in the integrand of Definition 1.5. The location of these points, as well as other important estimates are collected in Lemma A.8. Roughly summarizing, the finitely many points y satisfying ψ σ,θ (y) = λ have distance 1/ ln(λ) from the real axis. Away from such points ψ σ,θ (y) − λ λ holds and for y → ±∞ the function ψ σ,θ behaves doubly-exponential (Lemma A.4).

Abstract analysis of sinc-quadrature
In this section, we collect some results on sinc-quadrature formulas.

Remark 3.1
As is common in the literature, we define the sinc function as The following result is the main work-horse when analyzing sinc-quadrature schemes. In order to reduce the required notation, we use a simplified version of [32, Problem 3.2.6].

is a meromorphic function on the infinite strip D d(θ) . It is also continuous on
Denote by res(g; p ) the residue of g at p , and define γ (k; p ) := 1 sin(π p /k) . Then for all k > 0, using s := sign(Im( p )): (3.4) Proposition 3.2 requires certain decay properties for the integrand in a complex strip, and thus is not always applicable. As is shown in Appendix A, the transformation ψ σ,θ maps partly into the left-half plane. One can even show that the real part changes sign infinitely many times when evaluating along a line of fixed imaginary part. If we therefore consider the case when f (z) := e −z is the exponential function, this means that f • ψ is exponentially increasing in such regions. This puts showing estimates of the form required in Proposition 3.2 (iii) out of reach.
On the other hand, Lemma A.5 shows that for σ = 1, restricted to the domain D exp δ , the map ψ σ,θ stays in the right half-plane. Here the exponential function is decreasing. Similarly, the Mittag-Leffler function e α,μ is decreasing on slightly larger sectors, allowing for the choice of σ = 1/2 if α < 1. This motivates the following modification of Proposition 3.2.

Lemma 3.3 Assume that g : D
exp δ → C is holomorphic and is doubly-exponentially decreasing, i.e., there exist constants C g > 0, μ g > 0, such that g satisfies Then, for all 0 < ε < μ g/2, there exists a constant C > 0 which is independent of k, μ and g such that the following error estimate holds: Proof We closely follow the proof of [21, Theorem 2.13], but picking a different contour and later exploiting the strong decay properties of g.
For fixed t ∈ R, we fix N large enough such that t ∈ R N . By applying the residue theorem to the function one can show the equality Since asymptotically g(t) decreases doubly exponentially, while 1/ sin(π y/k) only grows exponentially along the path {(ξ, δ e −ξ ), ξ ∈ R}, we can pass to the limit N → ∞ to get the representation Integrating (3.7) over R and exchanging the order of integration gives: i sign(Im(y))π y k dy, (3.8) where in the last step we invoked [21,Lemma 2.19] to explicitly evaluate the integral. What remains to be done is bound the integral on the right-hand side. For simplicity, we focus on the upper-right half-plane. The other cases follow analogously. There, we can parameterize ∂ D exp δ as y = ξ + iδ e −ξ . We estimate For ε > 0, we can absorb the linear ξ -term into the first exponential, and estimate: where the second term will be used to regain integrability, whereas the first one will give us approximation quality. For ξ = 0 and ξ → ∞, we get sufficient bounds to prove (3.6). We thus have to look for maxima of the function with respect to ξ in between (0, ∞). Due to monotonicity of the exponential, we focus on the argument and set τ := e ξ . By setting its derivative to zero we get that the map .

Remark 3.4
It is also possible to admit meromorphic functions with finitely many poles into Lemma 3.3, as long as additional error terms analogous to (3.4) are introduced. Since we will not need this generalization we stay in the analytic setting.
While Lemma 3.3 provides a reduced rate of convergence compared to the more-standard sinc-quadrature of Proposition 3.2 (k −1/2 vs k −1 ), thus removing the advantage we want to achieve by using the double exponential transformation, we will later consider a class of functions which decay fast enough to allow us to tune the parameter μ ∼ k −1 to regain almost full speed of convergence.
Finally, we show how the transformation ψ σ,θ and the operator L enter the estimates. The next corollary also showcases how the cutoff error is controlled.

Corollary 3.5
Let O ⊆ C contain the right half-plane, and if σ = 1/2 also a sector Assume that g : O → C is analytic and satisfies the polynomial bound Then, for all ε > 0, s, r ∈ R such that μ − r + s − 2ε > 0, the quadrature errors can be bounded by: The constant C is independent of g, k,r ,s and β, but may depend on ε, σ , θ . The ratê γ depends on θ and ω. γ depends on σ .
Proof Let (λ j , v j ) ∞ j=0 denote the eigenvalues and eigenfunctions of the self-adjoint operator L. Following [3], plugging the eigen-decomposition of a function u into the Riesz-Dunford calculus, we can write the exact function g(L)u as and analogously for the discrete approximation Q L (g, N q )u. For the norm, as defined in (1.1), this means: We have thus reduced the problem to one of scalar quadrature, for which we aim to apply Lemma 3.3. We fix λ > λ 0 > κ. ψ σ,θ maps D exp δ analytically to O via Lemma A.5 (δ depends on θ and ω). What remains to be shown is a pointwise bound for the function By distinguishing the cases ψ σ,θ (y) < λ/2 and ψ σ,θ (y) ≥ λ/2 we get using either (A6) or Lemma A.5 We conclude using Lemma A.5: The double exponential growth of ψ σ,θ (see Lemma A.4) then gives after absorbing the cosh term by slightly adjusting ε: then gives, after readjusting ε: The cutoff error is handled easily, also using the estimate (3.10). We calculate where the last step follows by estimating the sum by the integral and elementary estimates.

The elliptic problem
In this section, we analyze the error when discretizing the elliptic fractional diffusion problem from Sect. 2.1. In order to analyze the quadrature error, we need to understand a specific scalar function. This is done in the next Lemma. Lemma 3.6 Fix λ > λ 0 > κ and β > 0. For y ∈ R, define the function Then the following statements hold: (i) g β λ can be extended to a meromorphic function on D d (θ) . It has finitely many poles. All poles p satisfy ψ σ,θ ( p) = λ and are all simple. For any ν ≥ 0, the number of poles within the strip can be bounded independently of ν, β and λ. The imaginary part of p can be bounded away from zero and for large λ, the following asymptotics hold: where the implied constants depend on θ , κ, and λ 0 . (ii) There exist constants C > 0, γ > 0, independent of λ and β and a value d λ ∈ The constant C may depend on β but can be chosen independently of λ and β.
Proof Proof of (i): We note that by Lemma A.3, ψ σ,θ is non-vanishing in D d(θ) . Since D d(θ) is simply connected, we may define It is easy to check that on R we have h(y) = ln(ψ σ,θ (y)) since the derivative as well as the value at y = 0 coincide. Thus, defining provides a valid meromorphic extension. The only poles are located where ψ σ,θ (z) = λ. By Lemma A.8 (i), the number of such poles within strips of width ln(λ) −1 is uniformly bounded. By Lemma A.3, ψ σ,θ has no zeros in the domain D d (θ) , which means all the poles are simple. The bound on the imaginary part follows from Lemma A.8 (ii).
Proof of (ii): We first note for y = a ± id λ , if λ < ψ σ,θ (y) /2, the trivial estimate Overall, we can estimate using Lemma A.4 where in the last step, we used that ψ σ,θ has the same asymptotic behavior as ψ σ,θ up to single exponential terms, which we absorb into the double exponential by slightly reducing γ .
Looking at |λ β 2 g β λ (y)|, one can calculate using two different ways to estimate ψ σ,θ (y) − λ: The integral bound then follows easily from the pointwise ones.
where the rate is given by Thus for k ∼ ln(N q )/N q we get (almost) exponential convergence: (3.14) The implied constants and γ may depend on λ 0 , β, σ , θ and κ.
Proof To cut down on notation, we only consider the case ln(λ/κ) ≥ c 1 > 1 so that the first term in the minimum of (3.11) dominates. If λ is small, the error can be absorbed into the e −γ /k term. The error E λ (z −β , N q ) corresponds to approximating g β λ by sinc quadrature. We split the error into two parts, the quadrature error and the cutoff error.
The term E c can be handled by the same argument as in Corollary 3.5. We therefore focus on the quadrature error E λ (z −β ) and apply Proposition 3.2. By Lemma 3.6(iii) it holds that N g β λ , D d(θ) < ∞. To satisfy assumption (ii), it suffices that (for sufficiently large y) the vertical strips do not contain any poles and we can use the asymptotics of Lemma 3.6(ii).
By Lemma 3.6, there are at most finitely many simple poles. The residue of the function at these poles can be easily calculated using the well-known rule provided that f is analytic and g (z 0 ) = 0. In our case this means, if ψ σ,θ (y λ ) = λ: where ζ ∈ N 0 denotes the branch of the complex logarithm picked by h. Thus, for a single pole y λ with s y λ := sign(Im(y λ )), recalling the definition of γ (k; y λ ) = 1 /sin(π y λ /k, we can estimate such that the number of elements in each bucket B is uniformly bounded (independently of λ, β and ). This allows us to calculate for the pole contribution in Proposition 3.2: where we used the elementary estimate 1 − e −2x min(x, 1) for x ≥ 0. Applying Proposition 3.2 and inserting this estimate for the pole-contributions gives: The bound from Lemma 3.6(iii) then completes the proof.
Proof We first show the estimate for ε > 0. We note that for ln(λ/κ) ≥ k −1 , we can bound the error in Theorem 3.7 by exp(−γ /k) (for an appropriate choice of constant γ ) due to the smallness of the term λ −β . Thus it remains to consider the case ln(λ/κ) < k −1 . Similarly, if ln(λ) ≤ max( c 2 ε , − ln(κ) p(σ,θ)−2ε ε , 1) =: μ 0 , the leading error term behaves like exp(−γ μ 0 k ). We are left to consider the remaining case. Writing μ := ln(λ), the error term can be estimated: We look for the minimum of the exponent. Setting the derivative of the map to zero, we get that the minimum satisfies Inserting this value into (3.16) gives the stated result (after slightly changing ε to get to the stated form). To see the case for ε = 0, we note that if ln(λ/κ) ≤ √ β+ρ−r , we can estimate for the leading term in Theorem 3.7: In the remaining case, we can estimate the higher order term in the ln(λ/κ)-asymptotics as We can also write λ −β+r −ρ = κ −β+r −ρ λ κ −β+r −ρ and continue as in the proof for δ > 0 but using μ := ln(λ/κ). This time we no longer have to compensate for the factors involving c 2 /μ and − ln(κ) by slightly reducing the rate. The price we pay is that the constant may blow up for ρ → ∞.
We can now leverage our knowledge about the function g β λ to gain insight into the discretization error for (2.5). This allows us to prove the two main theorems of this section. First we deal with the finite regularity case.
Proof of Theorem 2.5 Let (λ j , v j ) ∞ j=0 denote the eigenvalues and eigenfunctions of the self-adjoint operator L. Just as we did in the proof of Corollary 3.5, we plug the eigen-decomposition into the Riesz-Dunford calculus and Definition 2.1 to get for the discretization error: Applying Corollary 3.8 then gives for ρ ≥ 0 Next we prove the improved estimates for the case of G L (C f , R f , ω)-regularity.
Proof of Theorem 2.8 For simplicity of notation, we ignore the cutoff error, i.e., for now consider N q = ∞. The cutoff error can either be easily tracked throughout the proof or added at the end, analogously to Corollary 3.5.
We first note, that by Stirling's formula, we can estimate the derivatives of f by By assumption, we can apply Theorem 2.5 for any ρ ≥ 0. Picking ρ = δ k ln(k) 2 for δ sufficiently small and ε := p(σ, θ )/2 (because we need ρ-robust error estimates) gives: We need to show that the bracket in the exponential is positive. In order to do this, we expand the logarithmic term as ln 2δ This first term is negative, and for the others we note that 2ω √ δ |ln(k)| − ln(k) − 2 ln(| ln(k)|) + c 2 is uniformly bounded as | ln(| ln(k)|)| grows slower than | ln(k)| as k → 0. Due to the leading √ δ term, we can make δ small enough (independently of k) to ensure that the second term in the exponent of (3.17) is smaller than γ and the statement follows. If ω = 0, we don't have to compensate the factor e ωρ ln(ρ) , therefore picking ρ ∼ k −1 is sufficient and the improved statement follows.

The parabolic problem
Now that the stationary problem is well understood, we can move on to analyzing the discretization of the time dependent problem introduced in Sect. 2.2.

The Mittag Leffler function
The representation (2.8) hints that it is crucial to understand the Mittag-Leffler function if one wants to analyze the time dependent problem (2.7). We follow [20, Section 1.8].
For parameters α > 0, μ ∈ R, the Mittag-Leffler function is an analytic function on C and given by the power series . (3.18) We collect some important properties we will need later on. We start with the following decomposition result, also giving us asymptotic estimates.

Proposition 3.9
For 0 < α < 2, μ ∈ R and απ 2 < ζ < απ, we can decompose the Mittag-Leffler function as where R N α,μ is analytic away from zero and satisfies for a constant C > 0 depending only on z 0 and ζ .
where C can be taken as two rays {r ζ 0 : r ≥ 1}, {r ζ 0 : r ≥ 1} and a small circular arc connecting the two without crossing the negative real axis. ζ 0 is taken in the left halfplane such that the opening angle of C is sufficiently large in order to avoid possible poles of the integrand and ensure that the term (1 − t α /z) −1 is uniformly bounded. The stated result then follows easily by comparing the integral under consideration to the definition of the Gamma function.
Setting N = 1 in Proposition 3.9 and simple calculation yields the following estimates: For α = μ = 1, the Mittag-Leffler function e 1,1 is the usual exponential function. For the decomposition result, we can skip the terms involving powers z −n in this case as e z already decays faster than any polynomial.
Finally, we need a way of computing antiderivatives of the convolution kernel in (2.8).

Double exponential quadrature for the parabolic problem
The case of finite regularity In this section, we investigate the convergence of our method in the case that u 0 and f have finite H 2ρ -regularity for some ρ ≥ 0. It will showcase most of the new ingredients needed to go from the elliptic case to the time dependent one while keeping the technicalities to a minimum. The step towards Gevrey-regularity will then mainly consist of carefully retracing the argument and fine-tuning parameters. We start with the case if f = 0.

Proof
We start with N q = ∞ and split the Mittag-Leffler function according to (3.19).
We write For the first terms, we apply Theorem 2.5, and for the final term we use the decay estimate (3.20) and Corollary 3.5. Note that this is where we have to exclude the case α = β = 1 and σ = 1/2. If α < 1 the Mittag-Leffler function is contractive on a large enough sector. If β < 1, the map z → z β maps the required sector into the right half plane. Otherwise, the exponential function only decays in the right half-plane, not any slightly bigger sector. Thus, if σ = 1/2, Corollary 3.5 does not apply. Overall, we get the estimate: To simplify the calculations, we make use of the fact that β − r − 2ε ≥ β/2 − 2ε > 0 and ρ > 0. That way, the last term can be simplified to If η is an integer, we can pick N = η to get the statement for N q = ∞. For general η ≥ 1, we can interpolate between η and η + 1. The treatment of the cutoff error follows as in Corollary 3.5, exploiting that e α,μ (z) decays like (3.21) with s := β/2.
Picking η large enough, Lemma 3.11 shows that for fixed times t > 0 we get the same convergence rate as for the elliptic problem, though the approximation deteriorates as t gets small. Now that we understand the homogeneous problem, we can look at the case of allowing inhomogeneous right-hand sides f by using the representation formula (2.8), and finally prove the main result Theorem 2.10. We point out that naive application of Corollary (3.16) also inside the time-convolution integral would fail to give good rates, as the error may blow up faster than τ −α for small times, leading to a non-integrable function. Instead, the following proof relies on integration by parts and (3.22) to split the convolution into point evaluations similar to Lemma 3.11 and an integrable remainder term.

Proof of Theorem 2.10
As we have already estimated the error of the homogeneous part, we only consider the part corresponding to the inhomogenity, i.e., for now let u 0 = 0. We integrate by parts m times, using (3.22): Transferring this identity to the operator-valued setting, this means that we can analyze the quadrature error for these terms separately.
All the terms appearing are of the structure in Lemma 3.11. Most notably, the first m terms are evaluated at a fixed t > 0 thus we don't have to analyze them further and can just accept some t-dependence.
Investigating the remaining integral, we get by using η := m/α +q in Lemma 3.11: For q < 1, this is an integrable function (with respect to τ ) and the integral grows like t α(1−q) . We now focus on extracting the correct t dependencies. For small times, the dominating t-dependence in the estimates above can be found in the first term of (3.24), which behaves like t −mα(1−q) . If we put back the homogeneous contribution from Lemma 3.11, this term will dominate for small times like t −m−qα . For larger times, the initial error term in (3.24) is dominant, giving behavior T α . The cutoff error is treated like before, making use of the decay of e α,α . We just point out that the homogeneous cutoff error behaves like t −α/2 and the inhomogeneous part t α/2 . We crudely estimated both by max(t −m−qα , T α ) to simplify the statement of the theorem).

Remark 3.12
Corollary 2.10 shows that, as long as we assume that f is smooth enough in time we recover the same convergence rate p(σ, θ ) √ β + ρ − r as in the homogeneous and elliptic case.

The case of Gevrey-type regularity
If the data not only satisfies some finite regularity estimates but instead is even in some Gevrey-type class of functions, we can again improve the convergence rate, and almost get rid of the square root in the exponent. We go back to the homogeneous problem and assume that k < 1/2 so that the logarithmic terms can be written down succinctly.
Proof We go back to (3.23), but apply Theorem 2.8 to each of the first N terms, getting: We estimate the first N terms by For n ≤ δ k|ln(t )| 2 |ln(k)| 2 , we can estimate the exponent by (2) For δ small enough, depending on c 1 , α and γ , the term in brackets is uniformly positive (i.e., independently of t and k), we can thus estimate for some γ 1 > 0: The remainder term behaves like By picking N = δ k|ln(t )| 2 |ln(k)| 2 , the exponent be bounded up to a constant by By taking the factor δ sufficiently small, we get that the term in brackets stays uniformly positive, which shows The cutoff error can easily be dealt with as in the previous results, as the Mittag-Leffler function satisfies the decay bound (3.21) for s = 1/2.
Finally, we are in a position to also include the inhomogenity f into our treatment. This means we can prove the main result Theorem 2.11. Just as in Lemma 2.10, we use integration by parts to decompose the error into parts for positive times and a remainder integral with "nice enough" behavior with respect to τ .
Proof of Theorem 2. 11 We again work under the assumption u 0 = 0 and focus on the error when dealing with the inhomogenity f alone and also start with N q = ∞. We also for now take t ≤ 1.
Going back to (3.24) we get for N ∈ N 0 to be fixed later (3.26) For the first terms, we apply Lemma 3.13 to get exponential convergence, as long as f ( j) is in the right Gevrey-type class. Namely, we note that we can estimate  (3.27) Again restricting ω to absorb the factor N due to the summation.
For the remainder in (3.26), we look at the pointwise error at fixed 0 < τ < t, shortening f (N ) := f (N ) (t − τ ). Going back to (3.25), we can use the additional powers of t to get rid of the ln(t) term in the exponential: We then proceed as in the proof of Lemma 3.13, noting that since the τ -dependent terms can be bounded independently of N we can get by without the ln(t )-term in the exponent. Overall, we get by tuning N ∼ δ/(|ln(k)| 2 k) (also in (3.27)) appropriately: which easily gives the stated result. If t > 1, we can skip the integration by parts step for the integration over (1, t) and directly apply Lemma 3.13. The cutoff error is treated as always.

Numerical examples
In this section, we investigate, whether the theoretical results obtained in Sects. 3.2 and 3.3 can also be observed in practice. We compare the following quadrature schemes: (i) DE1: double exponential quadrature using σ = 1/2 and θ = 4, (ii) DE2: double exponential quadrature using σ = 1 and θ = 4, (iii) DE3: double exponential quadrature using σ = 1 and θ = 1, (iv) sinc: standard sinc quadrature (v) Balakrishnan: a quadrature scheme based on the Balakrishnan formula (vi) BURA: best uniform rational approximation For the double exponential quadrature schemes, we used k = 0.9 ln(r N q )/N q with r := 1 for β > 0.4 and r := 5 for β < 0.4. This makes the cutoff error decay like e −βr N 0.9 q , which is sufficiently fast to not impact the overall convergence rate. The factor 0.9 was observed to have some slightly improved stability compared to 1. The damping constant r was introduced to get good behavior for small β; see Sect. 4.3.
For the standard sinc-quadrature, the proper tuning of k and N q is more involved.
Following [4], we picked k = 2π d βN q with d = π/5. The Balakrishnan formula is only possible for the elliptic problem. It is described in detail in [5]. Following [5, Remark 3.1] we used where M is the number of negative quadrature points. This corresponds (in their notation) to taking s + := β/10, which was taken because it yielded good results (Fig. 4).

The pure quadrature problem
In this section, we focus on a scalar quadrature problem. Namely, we investigate how well our quadrature scheme can approximately evaluate the following functions using the Riesz-Dunford calculus (a) z −β and (b) e α,1 (−t α z β ) at different values  ∈ (4, ∞). This is equivalent to solving the elliptic and parabolic problem with data consisting of a single eigenfunction corresponding to the eigenvalue λ. Throughout, we used κ := 3. Theoretical investigations revealed, that the quadrature error is largest at ln(λ) ∼ k −1/2 (see the proof of Corollary 3.8). Therefore, we make sure that for each k under consideration, such a value of e 1 √ k is among the λ-values sampled. More precisely, the sample points consist of N max with k(N q ) = 0.9 ln(N q )/N q , and we consider the maximum error over all these samples. We used t := 1 for all experiments. We observe that for the most part, choosing σ = 1/2 and θ moderately large gives the best result. This agrees with our theoretical findings. This method fails to converge though if α = β = 1 is chosen as the parameters for the Mittag-Leffler function. This also agrees with the theory, because in this case, ψ σ,θ fails to map into the domain where e α,μ is decaying (see (3.21)). This shows that the restriction on σ in the theorems of Sect. 3.3 is necessary. If we only consider the elliptic problem, no such restriction is necessary, as the decay property is valid on all of the complex plane. All the other methods perform well in all of the cases. The straight-forward double exponential formula, i.e., σ = θ = 1, is often outperformed by the simple sinc quadrature scheme, (except in the α = β = 1 case of the exponential). For comparison, we've included the (rounded) predicted rate for the DE1 scheme in the plots. We observe that for several applications our estimates appear sharp. For f (z) = z −1 the scheme outperforms the prediction, but this might be due to a large preasymptotic regime. We note that for e −z β , we expect better estimates than the ones presented in this article to be possible due to the exponential decay. This is also true for the standard sinc methods, see [3].
Second, we look at the case of a single frequency λ and see how the convergence rate decays as λ → ∞. In order to better see the λ-dependence of the quadrature error, we consider the relative error of the quadrature, i.e., we look at E λ (z −β )/λ −β for β = 0.5. The theory from Theorem 3.7 predicts behavior of the form e − γ ln(λ)k , i.e., the rate drops like ln(λ). In Fig. 2, we can see this behavior quite well. In comparison, using standard sinc quadrature gives a λ-robust asymptotic rate, but only of order N q .

A 2d example
In order to confirm our theoretical findings in a more complex setting, we now look at a 2d model problem with more realistic data than a single eigenfunction. Namely, we work in the PDE-setting of Remark 1.2 using the unit square = (0, 1) 2 and the standard Laplacian with Dirichlet boundary conditions. We focus on two cases: first we look at what happens if the initial condition does not satisfy any compatibility condition, i.e., u 0 / ∈ H 2ρ for ρ ≥ 1/4. The second example is then taken such that the data is (almost) in the Gevrey-type class as required by Theorem 2.8 and Theorem 2.11. The inhomogenity in time is taken as f (t) := sin(t)u 0 , thus possessing analogous regularity properties. We computed the function at t = 0.1.
For the discretization in space and of the convolution in time of (2.8), we consider the scheme presented in [26]. It is based on hp-finite elements for the Galerkin solver and a hp-quadrature on a geometric grid in time for the convolution. As it is shown there, such a scheme delivers exponential convergence with respect to the polynomial degree and the number of quadrature points. Since we are not interested in these kinds of discretization errors, we fixed these discretization parameters in order to give good accuracy and only focus on the error due to discretizing the functional calculus. Namely, we used 5 layers of geometric refinement towards the boundary and vertices and a polynomial degree of p = 12.
Since the exact solution is not available, we computed a reference solution with high accuracy and compared our other approximations to it. The reference solution is computed by the DE1 scheme (as it outperformed the others) by using 8 additional quadrature points to the finest approximation present in the graph. As the DE1 scheme has finished convergence at this point, we can expect this to be a good approximation.
We start with the parabolic problem. The initial condition is given by For ω := 1, this function does not vanish near the boundary of and therefore only satisfies u 0 ∈ H 1/2−ε . We are in the setting of Lemma 2.10. By inserting ρ = 1/4 (up to ε) and r = 0, the predicted rates for DE1 and DE2 are roughly e − 6.13 √ k and e − 5.62 k respectively. Figure 3a contains our findings. We observe that all methods converge with exponential rate proportional to N q . The double exponential formulas outperforming the standard sinc quadrature. We also observe that picking σ = 1 and θ = 1 can greatly improve the convergence. The best results being delivered by DE1, i.e. σ = 1/2 and θ = 4. For DE1 and DE2, we observe that for a large part of the computation, the scheme outperforms the predicted asymptotic rate, but for DE2, the rate appears sharp for large values of N q .
As a second example, we used ω = 0.05. This function is almost equal to 0 in a vicinity of the boundary of . Thus we may hope to achieve the improved convergence rate of Theorem 2.11. Figure 3b shows that it is plausible that the exponential rate of order N q is achieved, and all the double exponential schemes greatly outperform the standard sinc quadrature. The best results are again achieved by DE1 and DE2, which also greatly outperform the predicted rate for the non-smooth case.

Elliptic problem and behavior for smallŤ
hus far, all our estimates worked under the assumption of β > β > 0. In order to shed some light on the behavior, and in addition gain insight into the behavior for the As geometry we again used the unit square. We chose f = 1 as the constant function. In this class, we also included the method based on the Balakrishnan formula as well as a rational approximation method, namely the one based on computing the best uniform rational approximation as described in [17]. Where we approximate z 1−β on [0, 1] using a rational function and then divide by z −1 and scale back to the interval [λ min , λ max ]. For computing the approximation we used the brasil algorithm described in [18], the implementation of which can be found in the baryrat python package [18]. To determine λ max , we used a simple power iteration with 10 iterations. This gave the estimate λ max ≈ 6 · 10 15 . For λ min we used the constant κ := 3 also used in the other algorithms.
For small β, preliminary experiments suggest a severe degrading of performance if the choice k := 0.9 ln(N q )/N q is made. Therefore it was necessary to introduce the constant r in our considerations. We point out that setting r := 1 for β > 0.4 is not fully necessary and only gives small improvements for larger values of β. Thus, if multiple values of β are of interest, in order to be able to reuse the approximate inverses (L − ψ σ,θ ( jk)) −1 , the choice of this damping factor should be according to the smallest value of β one is interested in.
In Fig. 4, we again observe that with θ = 4 and σ ∈ {0.5, 1}, the double exponential formulas significantly outperform the standard sinc based strategies, where σ = 0.5 again delivers the best performance. For comparison, we included the predicted rates for the DE1 and DE2 schemes into the graphics. We observe that asymptotically our estimates appear sharp, but with a large range of values, for which the scheme outperforms the predictions. The rational approximation method performs very well for small numbers of systems, but the performance degrades severely when higher accuracy is required. This instability with respect to numerical errors is most likely due to the requirement of rewriting the rational function in the partial fraction form to apply it to a matrix as described in [17] -even though a multiprecision library is Fig. 4 Comparison of approximation schemes for 2d elliptic problems with different parameter β used for computing the poles and residuals of the rational function. We also tried the method based on the AAA-algorithm [30], but there the numerical instability was even more problematic. If we talk about generic complex numbers without relation to any of the specific planes, we use the letter ζ instead.
We start out with some basic properties of sinh.

Lemma A.3 ψ σ,θ is analytic on the infinite strip D d(θ) . For d(θ ) sufficiently small, both ψ σ,θ and ψ σ,θ are non-vanishing on D d(θ) .
Proof The analyticity of ψ σ,θ is clear. In order to analyze the roots, we first rewrite for w = a + ib, separating the real and imaginary parts: We first focus on the case σ = 1. In this case, (A2) shows that any root y of ψ σ,θ must satisfy for w := π 2 sinh(y) =: a + bi: Since cosh has no roots, we get cos(b) = θ sin(b). As cos(b) = θ sin(b) and θ cos(b) = −sin(b) is impossible at the same time, we get that a = 0 and b = tan −1 (1/θ ) + π for some ∈ Z.
It remains to show that π 2 sinh(y) does not map to these points. Looking at the real part of sinh(y) we immediately deduce that if Im(y) ∈ (−π/2, π/2), in order to produce a purely imaginary result, it must hold that Re(y) = 0. For the imaginary part, we then get the equation: sin(Im(y)) = 2 + 2 tan −1 (1/θ ) π for some ∈ Z which is not possible for |Im(y)| ≤ d(θ ) < sin −1 ( 2 tan −1 (1/θ) π ). Next, we show that ψ σ,θ also does not vanish. A simple calculation shows Since the restriction θ ≥ 1 was not crucial for the proof, ψ 1, 1 θ and cosh have no roots in the symmetric (w.r.t. sign flip) domain D d (θ) . This shows that ψ 1,θ also is non-vanishing.
Just like we did when showing ψ 1 2 ,θ = 0 we can argue that a = 0. We get sin(b/2) = Im(t). Since t only depends on θ , we get that |b| > b 0 > 0 with a constant only depending on θ . We proceed as when showing ψ 1 2 ,θ = 0 to conclude that ψ 1 2 ,θ has no root in D d (θ) for d sufficiently small (depending on θ ).
We now look at how to adapt the proof to the case σ = 1/2. If |Re(w)| ≥ 2 ln(1 + √ 5), we get where in the last step we used the monotonicity of the expression and the fact that e |Re(w)| 2 − e |Re(w)/2| = 2 for Re(w) = 2 ln(1 + √ 5). The argument for the y − wtransformation stays the same. The upper bound also follows easily from the triangle inequality and the growth of sinh and cosh.
While on the full strip D d(θ) , the image of the transformation is difficult to study, the restriction to a certain subdomain is much better behaved.
Proof By Lemma A.2(iii) it is sufficient to consider the mapping of ϕ σ,θ restricted to small strips in the w-plane around the real axis. We start with the simpler case σ = 1. Going back to (A2) and writing w := π 2 sinh(y) =: a + ib, we note that if |b| is sufficiently small, we can guarantee that cos(b) − θ sin(b) > c > 0 for some constant c > 0 depending on θ .
Using (A6) then concludes the proof.
In order to apply the double exponential formulas for the Riesz-Dunford calculus, it is important to understand where ψ σ,θ (z) hits the real line. We start with the w-domain. Lemma A.6 Fix λ ≥ λ 0 > 1. Then the following holds for every w ∈ C with Re(w) = 0 and cosh(σ w) + iθ sinh(w) = λ : (A10) (i) There exist constants c 1 , c 2 , c 3 > 0 such that w satisfies log(λ) − c 1 ≤ |Re(w)| ≤ log(λ) + c 1 and where c 1 depends on λ 0 and θ , c 2 depends on λ 0 , and c 3 depends on θ . (ii) Given 0 < r < R, the number N w (λ, r , R) of points w satisfying (A10) with r ≤ |Im(w)| ≤ R is bounded uniformly in λ by The constant C depends only on θ . (iii) There exist at most four values p 1 , . . . , p 4 depending on λ, θ , and σ such that all points satisfying (A10) can be written as If w solves (A10) then −w does as well.
Proof We start with the simpler case σ = 1. By separating the real and imaginary parts as in (A2), we can observe that the critical points w = a + ib with a = 0 are located at This implies that |a| ∼ ln(λ), and we also see that for each , there are at most two such points, one in each half-plane. All the statements follow easily. Note that in (iii) only two families are needed. For the remainder of the proof we therefore focus on the case σ = 1/2. Proof of (i): We start with the bound on the real part and write w = a + ib. We note that if |a| > max 1, 2 ln 8 θ one can estimate using elementary considerations that e |a| /4 ≤ |sinh(w)| and e |a|/2 /θ ≤ e |a| /8. We then calculate: From this, the statement readily follows. The other direction is shown similarly: For |a| ≤ max 1, 2 ln 8 θ , we use the bound ϕ σ,θ (w) e |a| to see that λ = ϕ σ,θ (w) e |a| , giving that ln(λ) must be uniformly bounded. By taking c 1 large enough we can make the ln(λ) − c 1 negative, thus making the first estimate in (i) trivial. Since ln(λ) ≥ ln(λ 0 ) > 0, we can also immediately see | Re(w)| ≤ max 1, 2 ln( 8 /θ) + ln(λ).
The final bound on the real part of w then follows for c 1 := max ln(8/θ ), ln(9θ/8), c , where c is used to compensate for the case of small a.
Writing cosh 2 (w/2) = 1 + sinh 2 (w/2), we get that t := sinh(w/2) solves the quartic equation This means there can be at most 4 such values t 1 , . . . t 4 for any λ and it must hold that Here w j for j = 1, . . . , 4 is the solution to sinh(w j /2) = t j with Re(w j ) > 0 and minimal value of Im(w j ) . To see (ii), we note that for each t j at most ceil( R−r /4π) values lie in the sought after strip. Therefore we can estimate If λ ∈ [0, ] and |w| > C w := max(log(2 /c 1 ), 4 ln (2)), (where c 1 is the constant in (A7) or (A8)) we get: We therefore may from now on assume that λ is sufficiently large as we see fit. In preparation for the rest of the proof, we note that for ζ, μ ∈ R, w.l.o.g., |ζ | ≤ |μ|: Because it is much simpler, we start with the case σ = 1. We note that in this case M λ consists of the points mapped to ±λ. We distinguish three cases, depending on whether Re(w) is small and if Im(w) is close to a pole or not.
Proof of (iii): For d λ = d(θ ), we can not guarantee that ψ(ξ + i d λ ) does not hit the value λ. In this case, we have to modify d λ slightly to get robust estimates. For d ∈ R, consider the hyperbolas In the light of Lemma A.7 we need to ensure that dist(γ d λ , w p ) 1 for all w p ∈ M λ . We will be looking for d λ in a small strip around d(θ ). To simplify notation we define the length To make things symmetric with respect to the real axis, we consider M λ := M λ − M λ . It will therefore be sufficient to focus on the upper right quadrant of the complex plane. All other cases follow by symmetry. We write M y λ := sinh −1 ( 2 π M λ ) for the corresponding points in the y-domain. We start by noting that we can easily stay away from the problematic parts of the imaginary axis by making d(θ ) sufficiently small, as if |Re(sinh(y))| < ε we have |Im(sinh(y))| < (1 + ε) sin(Im(y)); thus for small real parts we can ensure to fit between (−b 0 , b 0 ) on the imaginary axis. This also means that we only consider points w λ ∈ M λ with |Re(w λ )| > ε > 0 since our path will have already positive distance to other possible poles.
By (i), the number of points y λ in M y λ in the strip d(θ ) − ω ≤ Im(y λ ) ≤ d(θ ) can be bounded by a constant N , independent of λ. In order to also avoid points in M y λ which are close but outside the critical strip we also avoid the boundary points d(θ ) − ω and d(θ ). Since N + 2 strips of width ω 2N +4 can not cover a strip of width ω, there exists a value d λ such that d(θ ) − ω ≤ d λ ≤ d(θ ) and |Im(y λ ) − d λ | ≥ ω 2(N + 2) ∀y λ ∈ M y λ .