Rigorous Computation of Diffusion Coefficients for Expanding Maps

For real analytic expanding interval maps, a novel method is given for rigorously approximating the diffusion coefficient of real analytic observables. As a theoretical algorithm, our approximation scheme is shown to give quadratic exponential convergence to the diffusion coefficient. The method for converting this rapid convergence into explicit high precision rigorous bounds is illustrated in the setting of Lanford’s map x↦2x+12x(1-x)(mod1)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x\mapsto 2x +\frac{1}{2}x(1-x) \pmod 1 $$\end{document}.


Introduction
For a real analytic 1 expanding interval map T : X → X with absolutely continuous invariant probability measure μ, and a real analytic function g : X → R, the corresponding diffusion coefficient (or variance) is the quantity σ 2 μ (g) defined by The quantity σ 2 μ (g) plays a role in the central limit theorem: it is well known (see e.g. [11]) that provided g is not equal to a coboundary plus a constant then 1 √ n n−1 i=0 g • T i − n g dμ converges in law to a normal distribution with mean zero and variance σ 2 μ (g) > 0. A difficult problem of practical interest is to calculate, or to approximate, the diffusion coefficient σ 2 μ (g), noting that (1) is only rarely amenable to direct evaluation. Bahsoun et al. [2] recently gave a method for the rigorous approximation of diffusion coefficients, including error bounds, based on Ulam's method. They illustrated this approach with the particular map introduced by Lanford [8], and the function g(x) = x 2 , showing that 0.3458 ≤ σ 2 μ (g) ≤ 0.4152.
In this paper we develop an alternative algorithm for approximating diffusion coefficients of expanding interval maps. In general the method uses the periodic points of T , and exploits the real analyticity of the map T and the function g. The method gives highly accurate approximations to the diffusion coefficient, both at the level of a theoretical algorithm converging with a given asymptotic speed (namely quadratic exponential convergence, as described in Theorem 1 below), and, most importantly, at the level of completely rigorous certified error bounds (see Theorems 2 and 3). The real analyticity assumptions will be crucial in establishing both the theoretical asymptotics and the concrete error bounds, since explicit use is made of the holomorphic extensions of the function g and (the inverse branches of) the map T to certain regions of the complex plane. The general asymptotic speed of our algorithm is as follows: Theorem 1 Let T : X → X be a real analytic expanding interval map with absolutely continuous invariant probability measure μ, and suppose g : X → R is real analytic. There exists a sequence σ 2 n → σ 2 μ (g), where each σ 2 n can be explicitly computed in terms of periodic points of period up to n. The rate of convergence is quadratic exponential, in the sense that there exist constants C > 0 and κ ∈ (0, 1) such that The constants C and κ of Theorem 1 can be rendered explicit, a procedure which involves consideration of holomorphic extensions to regions in the complex plane. A more challenging problem, in the context of a specific map T and function g, is to establish effective error bounds on |σ 2 μ (g) − σ 2 n |, preferably of very high accuracy; a key purpose of this article is to show that in such practical settings there is considerable scope for sharpening our optimal version of the simple asymptotic form (3) so as to obtain effective high quality bounds on the diffusion coefficient. As a model case we shall orient our discussion of this problem around the specific example considered in [2], namely Lanford's map T , and the function g(x) = x 2 , both of which are real analytic; henceforth we refer to this as the model problem. The problem of obtaining high accuracy rigorous estimates on σ 2 μ (g) involves both theoretical and computer programming elements, and any proof of such bounds will invariably be computer-assisted. As a starting point we note that, using only a modern desktop computer, it is possible to locate all the periodic points of the Lanford map T up to period P, for some 2 20 ≤ P ≤ 30. Choosing maximum period P = 25 yields the sequence of approximations to σ 2 μ (g) given in Table 2, which at the level of non-rigorous empirical observation suggests that σ 2 μ (g) = 0.36010948619916067289882418682857674924166999779722 ± 10 −50 , and indeed a more optimistic interpretation of Table 2 suggests the slightly more accurate σ 2 μ (g) = 0.3601094861991606728988241868285767492416699977972288644 ± 10 −55 .
The task is to now harness these computed approximate values σ 2 n (and in particular the last of these computed approximations, σ 2 P ) so as to produce a fully rigorous approximation to σ 2 μ (g), together with an error bound. Any naive expectation that the theoretical asymptotic (3), together with specific values for κ and C, would automatically yield an effective error bound on |σ 2 μ (g) − σ 2 n | is tempered by the realisation that, for the model problem, κ is reasonably close to the value 3 1, and C is extremely large. 4 Although, as noted above, the value P = 25 is deemed to be the maximum such that all (σ 2 n ) P n=1 can be explicitly evaluated, a finer analysis 5 of the estimates yielding the asymptotic (3) suggests that a good quality rigorous effective estimate on σ 2 μ (g) remains out of reach for P ≤ 30. In order to obtain high quality effective estimates on |σ 2 μ (g) − σ 2 n | we therefore develop a hybrid approach, consisting of three distinct types of computation, the first type being the exact evaluation of σ 2 n (see Sect. 3 for the formulae defining σ 2 n ) for all sufficiently 2 In general the specific value of P will depend on available hardware, on the computer programming implementation of our algorithm, and on the time available to make the computation. For the Lanford map T we found it possible to locate points up to period 20 in less than an hour, while locating points of period up to 25 took around a day (computations were performed in an arbitrary precision environment, giving several hundred correct decimal digits); note that since T is a 2-branch map, incrementing the maximum period by one entails an approximate doubling of the computer run time. 3 For any two branch expanding map, our techniques yield a value of κ lying in the range [2 −1/2 , 1), while for the Lanford map itself our optimal value is κ ≈ 0.927734 (this is the square root of the quantity θ defined in (57)). Note that although the term κ n 2 is approximately 4.3 × 10 −21 when n = 25, the value of C in (3) is too large for the asympotic estimate |σ 2 μ (g) − σ 2 n | ≤ Cκ n 2 to be effectively used until n is significantly larger (and, crucially, above the maximum value of n for which all 2 n period-n points can be located using the computational resources at our disposal). 4 The size of C will depend on κ, and C becomes larger the closer κ is chosen to the optimal value of approximately 0.927734 (see Footnote 3). As an indication of its order of magnitude, we use the fact (see Footnote 5) that |σ 2 μ (g) − σ 2 n | is related to (and in fact somewhat larger than) the quantity K n 1/20 E n (θ ) We are at liberty to work with any κ ∈ (θ 1/2 , 1) = (0.927734, 1), and for example with the concrete choice κ = 0.95 we can compute sup n∈N K n 1/20 E n (θ )/0.95 n 2 ≈ 4.440429 × 10 10 (the supremum is attained at n = 26), so that K n 1/20 E n (θ ) ≤ C κ n 2 for C = 4.5 × 10 10 . It follows, after some additional calculations (along the lines of those detailed in Sect. 8), that the value of C in (3) could be chosen to be of the order of 10 11 when κ = 0.95. 5 This finer analysis consists of using what we call Euler bounds, with the quality of the estimate on |σ 2 μ (g) − σ 2 n | closely related to the size of the quantities K n t E n (θ ) = K n t ( n i=1 (1−θ i )) −1 θ n(n+1)/2 given in Appendix Tables 5 and 6 (for t = 0 and t = 1/20 respectively), where θ = κ 2 ≈ 0.860691, K 0 ≈ 3.378, K 1/20 ≈ 3.631. We note that for sufficiently small values of n, the quadratic exponential decay of the term θ n(n+1)/2 is swamped by the exponential increase of the term K n t , and the strong increase of ( n i=1 (1 − θ i )) −1 (though this latter term is bounded, by ( ∞ i=1 (1 − θ i )) −1 ≈ 8876.45). In particular, for n = 20 both K n t E n (θ ) terms are greater than 1 (hence n = 20 represents a hopeless case for this naive method), while if n = 25 then K n 0 E n (θ ) ≈ 0.000084 and K n 1/20 E n (θ ) ≈ 0.00051, which in fact can be used (via arguments similar to those used in the proof of Theorem 3) to justify only a single decimal digit of σ 2 μ (g). small values of n (i.e. for all 1 ≤ n ≤ P, where e.g. P = 25 for the model problem), using exact locations of periodic points (i.e. evaluated to a given precision, typically several hundred decimal places). We next make the observation (see Corollary 1(b)) that σ 2 μ (g) can be expressed in terms of certain infinite series; it turns out that there are five such series, which for convenience we denote here 6 n of these series, and each of these is O(κ P 2 ) as P → ∞, a result which incidentally leads to the proof of Theorem 1. A consequence is that the task of obtaining a concrete bound on |σ 2 μ (g) − σ 2 P | reduces to bounding each tail ∞ n=P+1 s ( j) n , and here we note that the previously described difficulties in bounding |σ 2 μ (g) − σ 2 P | (for e.g. P = 25 in our model problem) stem from the natural upper bounds on the terms s ( j) n being insufficiently sharp for n ≈ P.
Our resolution of this problem of insufficiently sharp bounds consists of splitting the tails n we require some new techniques, whose justification (see Sect. 6) stems from the theory of eigenvalues and approximation numbers applied to a certain auxiliary (transfer) operator; these techniques require a non-trivial amount of computation, though a key point is that the computational effort is relatively light in comparison to that required for locating the 2 n period-n points for some high value of n (e.g. n ≈ P). The coefficients s ( j) n are related to the Taylor series for the determinant of the transfer operator, and can be bounded in terms of the approximation numbers of the operator. These approximation numbers can in turn be bounded by making a judicious choice of basis for an underlying Hilbert space whose inner product is defined by Lebesgue integration, and explicitly computing the norms of the images of (finitely many of) these basis elements under the transfer operator yields a bound on the approximation numbers which implies a bound on the s ( j) n for P + 1 ≤ n ≤ Q (see Sect. 7 for further details).
In Sect. 8 we combine all of these various ingredients, in the context of the model problem, to obtain the following rigorous bound on the diffusion coefficient, noting that it represents a significant improvement 8 on the estimate 0.3458 ≤ σ 2 μ (g) ≤ 0.4152 established in [2] for the same combination of function g(x) = x 2 and Lanford map T . 6 In terms of the later notation, these series correspond (see Corollary 1(b)) to ∞ n=1 nc n (0), ∞ n=1 n(n − 1)c n (0), ∞ n=1 c n (0), ∞ n=1 nc n (0), and ∞ n=1 c n (0), which themselves correspond to partial derivatives of the determinant of a (transfer) operator. 7 As will become clear, one virtue of this method is that it perfectly feasible, from a computational point of view, to choose Q rather large (e.g. some value well over 100), a choice which may be important for expanding maps T for which the expansion is rather mild, corresponding to significant inertia in the quadratic exponential decay of the terms s ( j) n , stemming from a value κ ∈ (0, 1) being close to 1. 8 While the rigorous estimate of [2] is less accurate than that of Theorem 2, the general strategy of [2] is based on Ulam's discretization method [16] and can be applied to a wider class of maps T and functions g for which there is no analyticity assumption (see [2] for details and references, and e.g. [10] for a further guide to the literature on numerical computations in the context of piecewise expanding maps).

Theorem 2 For the Lanford map T , with absolutely continuous invariant probability mea-
The organisation of this article is as follows. Section 2 consists of preliminary material drawn from the ergodic theory of expanding maps, thermodynamic formalism, and Hilbert spaces of holomorphic functions. Our algorithm is described in Sect. 3, together with various reformulations of the diffusion coefficient. The rapid convergence of the algorithm is illustrated in Sect. 4 for certain cases where σ 2 μ (g) is known explicitly, and in Sect. 5 for the model problem, where σ 2 μ (g) does not have a (known) closed form. The key theoretical tools for deriving rigorous error estimates, based on the theory of eigenvalues and approximation numbers, are developed in Sects. 6 and 7. These tools are then applied in detail to the model problem in Sect. 8, proving a result (Theorem 3) that is slightly stronger than Theorem 2, and concluding with a proof of Theorem 1. Some of the numerical data used in the proof of Theorem 3 is collected as an Appendix.

Ergodic Theory of Expanding Interval Maps
Suppose the unit interval X = [0, 1] is partitioned as Given T : X → X , we shall always assume that T | X i is real analytic, for each 1 ≤ i ≤ l. We say T is piecewise expanding if there exists λ > 1 such that |T | | X i ≥ λ for all 1 ≤ i ≤ l. We say that T is Markov if there exists a d × d matrix A (the transition matrix for T ) with each entry either 0 or 1, such that is called the Markov partition for T . T is topologically mixing if some power of the transition matrix A is a strictly positive matrix.
It is well known (see [9]) that any topologically mixing piecewise C ω expanding Markov map admits a unique ergodic absolutely continuous invariant probability measure, and we shall denote this measure by μ. Our results are valid for all such maps, though to simplify the exposition we shall always assume that T is a full branch expanding map. In other words, each T | X i is assumed to be a surjection onto X , or, equivalently, every entry of the corresponding transition matrix A is a 1. For each 1 ≤ i ≤ l we write as the collection of inverse branches of T . Since T is expanding, each inverse branch is a contraction mapping on X ; indeed the real analyticity 9 of T ensures that the inverse branches have a holomorphic extension to some common complex neighbourhood of X on which they are all contraction mappings.
. . , T n−1 (x)) ∈ X n : T n (x) = x} denote the collection of periodic orbits of (not necessarily least) period n, considered as ordered ntuples. For x ∈ O n and g : X → R, define and for n ≥ 1, t ∈ C, define a n (t) := a g,n (t) = 1 n x∈O n For a continuous function f : X → R, its pressure P( f ) = P( f, T ) is defined (see e.g. [13]) by

The Diffusion Coefficient
Suppose g : X → R is real analytic. Its diffusion coefficient (or variance) σ 2 μ (g) is defined by The diffusion coefficient can be expressed in terms of pressure as follows: Lemma 1 Let T : X → X be a real analytic expanding interval map with absolutely continuous invariant probability measure μ, and suppose g : X → R is real analytic. If p(t) := P(tg − log |T |), then the integral of g with respect to μ is given by and the diffusion coefficient is given by

Holomorphic Extensions
As noted in Sect. 2.1, the inverse branches of the real analytic expanding map T extend as contraction mappings to some common (simply connected) complex neighbourhood U of X . If g : X → R is real analytic then U may be chosen so that g is holomorphic on a neighbourhood of U . By the Riemann mapping theorem, no generality is lost by assuming that U can be chosen to be a disc D, and henceforth we make this assumption: an open disc D ⊂ C containing X will be called admissible (for the map T and function g) if g has a holomorphic extension to a neighbourhood of D, and each inverse branch τ i has a holomorphic extension to D such that ∪ l i=1 τ i (D) ⊂ D. This will allow consideration of transfer operators acting on certain Hilbert spaces of holomorphic functions.

Transfer Operators and Determinants
For a real analytic function g : X → R, an important ingredient in our method of approximating the diffusion coefficient σ 2 μ (g) is the function g : for sufficiently small values of z, and by analytic continuation to the whole of C 2 . It can be shown that (7) defines an entire function (see Corollary 4), with Taylor series expansion (where we write c n (t) for c g,n (t) whenever g is understood), from which we deduce the recurrence relation For T and g with holomorphic extensions to D as in Sect. 2.3, the corresponding transfer for z ∈ X , and by holomorphic continuation for z ∈ D, where f = tg − log |T |. The function g is the determinant det(I − zL g,t ) (see [12]), and its zeros are precisely the reciprocals of the eigenvalues of L g,t . The leading (i.e. largest in modulus) eigenvalue of L g,t is e p(t) = e P(tg−log |T |) .

The Diffusion Coefficient in Terms of Derivatives of the Determinant
The reformulation (6) of the diffusion coefficient σ 2 μ (g) in terms of pressure, together with the fact that e − p(t) is a zero of g (·, t), suggests the possibility of representing σ 2 μ (g) in terms of partial derivatives of g . In order to establish such a representation, as Proposition 1 below, we first adopt the following notational conventions:

Notation 2
We write first partial derivatives as and second partial derivatives as Proposition 1 If T : X → X is a real analytic expanding map with absolutely continuous invariant probability measure μ, the function g : X → R is real analytic, and the determinant g is defined by (7), then the diffusion coefficient σ 2 μ (g) can be expressed as , so that p(0) = 0 and therefore Differentiating gives so Lemma 1 gives Differentiating again gives so evaluating at t = 0 and using Lemma 1 gives or in other words The zeros of g (·, t) are the reciprocals of the eigenvalues of L g,t , and since e p(t) = z(t) −1 is the leading eigenvalue of L g,t then so differentiating (15) with respect to t gives and therefore Combining (12) and (17) gives Differentiating (16) with respect to t gives and evaluating this at t = 0 then using (11), (12) and (14), gives in other words which is the required expression (10).
Corollary 1 Under the same hypotheses as Proposition 1, .
can be expressed as can be expressed as .
(e) If g : X → R is real analytic such that g dμ = 0, andσ 2 N is defined bŷ Proof Part (a) follows from Proposition 1, by substituting (18) into (10). Since the Taylor series around 0 for g (·, t) is written (cf. (8)) as g (z, t) = 1 + ∞ n=1 c n (t)z n , termwise differentiation yields (b). Part (d) is a special case of formula (10) in Proposition 1, together with (b), while parts (c) and (e) follow directly from the definitions of σ 2 N andσ 2 N .
Remark 1 A consequence of Corollary 1 is that if g is known to have integral zero with respect to the absolutely continuous invariant probability measure μ, then there is a choice of sequence of approximants to the corresponding diffusion coefficient: both the sequence σ 2 N and the sequenceσ 2 N converge to σ 2 μ (g).

Periodic Orbit Formulae
The quantities σ 2 N approximating the diffusion coefficient σ 2 μ (g) are accessible to us in terms of those periodic points of T of period up to N . Recall from (4) that a n (t) = a g, We are interested in derivatives up to order 2, evaluated at t = 0, so for n ≥ 1 define α n := a n (0), β n := a n (0), γ n := a n (0), in other words

Computer Implementation
Although in certain special cases (e.g. the doubling map of Sect. 4) the periodic points of T are rational and known explicitly, more generally a non-trivial aspect of our algorithm is to locate these periodic points (to within a specified precision). 10 For this, note that for 1 ≤ i ≤ l the inverse branch τ i : X → X i , defined as τ i = (T | X i ) −1 , is uniformly contracting. For each ξ ∈ {1, . . . , l} n the composition τ := τ ξ 1 • . . . • τ ξ n is also uniformly contracting, and the set of period-n points for T is precisely the set of fixed points of such compositions τ . The fixed point for the contraction mapping τ can be determined using standard techniques (e.g. choose x 0 ∈ X and evaluate x : Having located the period-n points of T , and formed the collection O n , for all 1 ≤ n ≤ N , the calculation of orbit sums a n (0) and their derivatives a n (0), a n (0) is then possible (using (4) and (22)) for n = 1, . . . , N . Differentiation of the recurrence relation (9) yields recurrence relations for the derivatives c n (0) and c n (0) which can then be computed for n = 1, . . . , N , and substitution into (20) gives the approximant σ 2 N .

Test Cases: Approximation of Known Diffusion Coefficients
For certain combinations of map T and function g, the diffusion coefficient is known exactly. While for these cases there is clearly no need for a numerical algorithm to approximate σ 2 μ (g), it is nonetheless instructive to consider them, by way of a warm-up exercise.

Perfect Approximation of the Diffusion Coefficient via Periodic Orbits
As a simple first example we describe here an expanding map T and function g whose diffusion coefficient σ 2 μ (g) is exceedingly well approximated by the sequence σ 2 n : in fact it turns out that each σ 2 n is equal to σ 2 μ (g). Let T : X → X be the doubling map, defined by T (x) = 2x (mod 1) on [0, 1) and T (1) = 1, its absolutely continuous invariant probability measure μ being Lebesgue measure itself. Consider the function g : X → R defined by g(x) = 2x − 1, which clearly satisfies g dμ = 0. In fact g is cohomologous to the function h defined by and it is easily seen that 1 n ( n−1 i=0 h • T i ) 2 dμ = 1 for all n ≥ 1, so the corresponding diffusion coefficient is given by the exact formula While the existence of an exact formula for σ 2 μ (g) means there is no need for numerical approximations, this example has the noteworthy feature that our approximations σ 2 n are perfect for each value of n: Proposition 2 For T : X → X the doubling map, and g(x) = 2x − 1, so in particular Proof If n ≥ 1 and x ∈ O n then m x = 2 n , and O n has cardinality 2 n , so α n = 2 n n(2 n −1) .
Now g (z, t) = exp − ∞ n=1 a g,n (t)z n for z of sufficiently small modulus, therefore and setting t = 0, so that a g,n (0) = β n = 0 for all n ≥ 1 by (26), gives for z of sufficiently small modulus.
Comparing (28) and (29), which are valid for z of sufficiently small modulus, gives which is in fact valid for all z ∈ C, by analytic continuation, since both sides of the equation are entire functions of z. Writing g (z, t) = 1 + ∞ n=1 c g,n (t)z n we deduce from (30) that ∞ n=1 c g,n (0)z n = ∞ n=1 nc g,n (0)z n , and the required equality (24), and hence (25), follows by comparing coefficients.

Remark 2
The setting of Proposition 2 allows an explicit illustration of the quadratic exponential decay of the coefficients c n (0) (and hence of the c n (0) = nc n (0)). Writing we see that , and therefore 11 so in particular

Rapid Approximation
Suppose, as in Sect. 4.1, that T : X → X is the doubling map, and now define g : X → R by g(x) = x 2 . Clearly the integral of g is known explicitly, namely g dμ = 1/3, and if 11 In fact a slight sharpening of (31) gives |c n (0) f • T n f dμ (see e.g. [3]) gives (33) More generally, for T the doubling map, we note in passing that there are exact formulae for the diffusion coefficient of monomials x k (e.g. for g(x) = x 3 it can can be shown that σ 2 μ (g) = 2783 11760 ), and indeed for general polynomials, which can be derived from the following result: Proposition 3 Let T : X → X be the doubling map, and μ Lebesgue measure. If B k denotes the k-th Bernoulli polynomial, then its diffusion coefficient is given by Proof The Bernoulli polynomial B k is an eigenvector of the Perron-Frobenius operator L, with corresponding eigenvalue 2 −k , since it is readily checked that the generating function LG(x, y) = G(x, y/2) (see [6]). Now σ 2 g. [1]), so the result follows.

Notation 3 Writing the Lanford map T as
we see that its real analytic inverse branches τ i : X → X i are given by 13 Table 2 Quadratic exponential convergence of approximations σ 2 n (formed using periodic points of period up to n) to the diffusion coefficient σ 2 μ (g) for the Lanford map T (x) = 2x + 1 2 x(1 − x) (mod 1) with absolutely continuous invariant probability measure μ and function g(x) = x 2 . Convergence is O(κ n 2 ) as n → ∞ for some κ < 1, and it appears that κ may be chosen to be approximately equal to location of the period-n points for T , since they are fixed points of suitable compositions τ ξ 1 •· · ·•τ ξ n . The computational procedure for locating the collection of period-n points is very swift for smaller values of n; high level software packages such as Mathematica or Matlab may be used for this purpose, though the exponential growth in the number of period-n points makes it advantageous to use an imperative programming language for larger n. For this paper the location of periodic points for the Lanford map was carried out on a personal computer with Fortran F07 compiler using MPFUN MPFR packages by D. Bailey (allowing for thread-safe arbitrary precision computations), for all periods n up to n = 25, following the algorithm as described in Sect. 3. (b) It is noteworthy that in all cases studied (those of Sect. 4 as well as this section), the approximations to σ 2 μ (g) are rather inaccurate (e.g. not correct to 2 decimal places) until points of period at least 5 are incorporated into the approximation. (c) In [2, §4.6], a non-rigorous experiment is performed, which seems to suggest that σ 2 μ (g) lies in [0.361, 0.363]. This contrasts with our sequence of approximations in Table 2, and in particular with our best approximation σ 2 25 . It follows from Theorem 3 that the approximation error of the experiment in [2] is at least 10 −3 . (d) Note that min x∈X T (x) = T (1) = 3/2, corresponding to the fact that 2/3 is the largest value attained on X by the derivatives of the inverse branches τ i . The value 2/3 appears to be significant concerning the rate at which the approximants σ 2 n approach the diffusion coefficient σ 2 μ (g). Assuming σ 2 25 to be approximately equal to σ 2 μ (g), so that δ n := |σ 2 n − σ 2 25 | ≈ |σ 2 n − σ 2 μ (g)|, we note that the values ε n := exp(n −2 log δ n ) are close to √ 2/3 (e.g. ε 22 ≈ √ 0.668617, ε 23 ≈ √ 0.667478, ε 24 ≈ √ 0.666508), and we therefore tabulate l n := [log √ 2/3 |σ 2 n − σ 2 25 |] in Table 2 to illustrate the quadratic exponential convergence. In view of this, it is unsurprising that the value √ 2/3 (or some value rather close to it) also appears to dictate the quadratic exponential decay of the coefficients c n (t) of the determinants g (z, t) = 1 + ∞ n=1 c n (t)z n : for example the c n (0) in Appendix Table 7 are such that the terms κ n := exp(n −2 log |c n (0)|) appear to be converging to a value close to, or equal to, √ 2/3 (e.g. for 20 ≤ n ≤ 25 the κ n are approximately √ 0.674246, √ 0.673346, √ 0.672570, √ 0.671899, √ 0.671313, √ 0.670801 respectively). On the basis of the observed behaviour for the Lanford map, and for the doubling map in Sect. 4, one might speculate that for more general real analytic maps T : X → X and functions g : as n → ∞ for all κ > κ T , and that lim n→∞ exp(n −2 log |c g,n (t)|) = κ T for all t ∈ C. This would constitute a strengthening of our results that if θ ∈ (κ 2 T , 1) is the contraction ratio for an admissible disc, then |σ 2 n − σ 2 μ (g)| = O(κ n 2 ) as n → ∞ for all κ > θ 1/2 , and lim sup n→∞ exp(n −2 log |c g,n (t)|) ≤ θ 1/2 for all t ∈ C.

Eigenvalues and Approximation Numbers
In this section we recall the definition of approximation numbers s n (L g,t ) for the transfer operator L g,t , and introduce a sequence α n (t) of upper bounds for s n (L g,t ), which we call approximation bounds. By then defining the associated contraction ratio θ ∈ (0, 1) we are able to establish (see Corollary 2) the exponential bound s n (L g,t ) ≤ α n (t) ≤ K t θ n , for a certain explicit constant K t > 0, which in particular will facilitate (see Corollary 3 in Sect. 7) a proof of the quadratic exponential decay of the Taylor coefficients for the associated determinant.
Let D ⊂ C be an open disc of radius centred at c, and let {λ n (t)} ∞ n=1 denote the eigenvalue sequence for the operator L g,t : H 2 (D) → H 2 (D), with the convention that eigenvalues are ordered by decreasing modulus and repeated according to their algebraic multiplicities. The Taylor coefficients c n (t) of g (·, t) then satisfy (see e.g. [15,Lem. 3.3]) the identity c n (t) = If, for k ≥ 0, we define m k : D → C by then {m k } ∞ k=0 constitutes an orthonormal basis for H 2 (D). For n ≥ 1 we can define the corresponding nth approximation bound α n (t) by and these values yield a simple upper bound on the approximation numbers of the transfer operator: Lemma 2 For n ≥ 1, the nth approximation number of L g,t : g,t := L g,t P n , where P n : and the Cauchy-Schwarz inequality then implies g,t has rank n − 1, so the required inequality (39) follows.

Definition 2
Let D be the smallest disc, concentric with D, such that ∪ l i=1 τ i (D) ⊂ D , and , the respective radii of D, D . The corresponding contraction ratio θ = θ D is defined to be θ = θ D := / . Lemma 3 Let D be an admissible disc, with contraction ratio θ = θ D . If g : (0, 1) → R has a holomorphic continuation to a bounded function on D, and each |τ i (·)| has a holomorphic continuation to a bounded function on D, then for all k ≥ 0, Now each τ i (D) is contained in the disc D , with the same centre c as D, and of radius = θ , therefore |C i (m k )(z)| = −k |τ i (z) − c| k < −k ( ) k = θ k for all z ∈ D. It follows that C i (m k ) ≤ θ k and combining this with (41) gives the required inequality (40).

Corollary 2 Under the hypotheses of Lemma
Proof Combining (38) and Lemma 3 gives while the inequality s n (L g,t ) ≤ α n (t) is the content of Lemma 2.

Euler Bounds and Computed Bounds
In this section we introduce two different kinds of bound on the Taylor series coefficients of the determinant g (·, t). The first of these, the Euler bound, has a simple closed form and is readily seen to converge to zero at a quadratic exponential rate. This implies the quadratic exponential decay of the Taylor coefficients (see Corollary 3), and hence that the determinant is an entire function (see Corollary 4); importantly, the inequality proved in Corollary 3 is subsequently used in Sect. 8 to rigorously bound one part of the error term in our diffusion coefficient approximation. The second kind of bound on the Taylor coefficients of the determinant is based on the approximation bounds α n (t) introduced in Sect.
In view of the following bound (44), and the fact that the identity in (43) was first given by Euler (cf. [5,Ch. 16]), we shall refer to the quantity K n t E n (θ ) as the Euler bound on the n th Taylor coefficient of the determinant g (·, t).
Proof The asymptotic (45) is immediate from (44), and this in particular implies that the Taylor coefficients of g (·, t) tend to zero faster than any exponential, hence the function is entire.
In order to exploit Lemma 2, which asserts that we require a practical means of computing the approximation bound α n (t). This will consist of bounding ∞ k=n−1 L g,t (m k ) 2 by the sum of an exactly computed long finite sum N k=n−1 L g,t (m k ) 2 (the H 2 (D) norms of the summands can be evaluated using numerical integration, 14 since each L g,t m k is known in closed form) and a rigorous upper bound on ∞ k=N +1 L g,t (m k ) 2 using (40). With this in mind, for n, N ∈ N with n ≤ N , we define the lower computed approximation bound and the upper computed approximation bound α g,n,N , Lemma 4 For t ∈ C, and n, N ∈ N with n ≤ N , Proof The inequality α n,N , 1−θ 2 , and the inequality follows. To prove that α n,N ,+ (t) ≤ K t (1 + θ 2(N +2−n) ) 1/2 θ n , note that combining (42) with α n,N ,− (t) ≤ α n (t) gives α n,N ,− (t) ≤ K t θ n , so (47) gives where the sum is over those i = (i 1 , . . . , i n ) ∈ N n which satisfy i 1 < i 2 < . . . < i n , and the sequence (α M n,N ,+ (t)) ∞ n=1 is defined by: Note that from (42), (48) and (50) we have s n (L g,t ) ≤ α M n,N ,+ (t), which combined with (36) establishes that the Taylor bounds β M n,N ,+ (t) are indeed bounds on the modulus of the n th Taylor coefficient of g (·, t): As computable approximations to β M n,N ,+ (t) we then define the lower computed Taylor bound by β M,− g,n,N , In practice the sum on the righthand side of (53) will be extremely small, though is sufficient for the upper computed Taylor bound to be an upper bound on |c n (t)|:

Validated Numerics: The Lanford Map
With the theory of Sects. 6 and 7 in hand, we are finally in a position to rigorously justify the quality of our computed approximation (see Sect. 5) to the diffusion coefficient σ 2 μ (g) in the case of our model problem, namely the case of T the Lanford map, μ its absolutely continuous invariant probability measure, and g : X → R the function g(x) = x 2 . In Sect. 8.1 we choose a suitable disc D, compute the associated contraction ratio θ and constants K t , and make choices of the natural numbers M, N , Q which arise in connection with the computed Taylor bounds of Sect. 7. In Sect. 8.2 we establish (see Proposition 5) rigorous bounds on the tails of five series which arise in the formula for σ 2 μ (g) derived in Corollary 1; each series represents a certain derivative of the determinant, and the bounds are established via our Euler bounds and computed Taylor bounds on its Taylor coefficients c n (t). In Sect. 8.3 these tail estimates are combined with the exact evaluations of the corresponding truncated series obtained via periodic point calculations (as described in Sect. 5) to prove a rigorous bound on σ 2 μ (g) (see Theorem 3). In Sect. 8.4 we prove Theorems 1 and 2, which were stated in Sect. 1; Theorem 2 is seen to be a minor variant of Theorem 3, while the more abstract Theorem 1 is established by combining the techniques used to prove Theorem 3 with the Taylor series asymptotic (45) from Corollary 4.

Computed Approximation Bounds and Computed Taylor Bounds
We shall be particularly interested in the choices t = 0 and 16 t = 1/20, and in these cases the supremum norm on D for both functions w 0,t and w 1,t is attained by evaluating at z = c + , g(τ 0 (c + )) τ 0 (c + ) 15 We make this choice so as to minimise the error estimates arising from the computed Taylor bounds. 16 The choice t = 1/20 is close to optimal for the purpose of estimating c n (0) and c n (0) via Cauchy's integral formula in the proof of Proposition 5. This involves, respectively, the integration of c n (ζ )ζ −2 and c n (ζ )ζ −3 over a circular contour centred at 0, and for both integrands there is a tension between the bound on |c n (ζ )|,which increases with |ζ |, and the bound on |ζ −k | (for k = 2, 3), which decreases with |ζ |. = 0.575158859780423676330133482123073962520 . . . , We can then compute K 0 and K 1/20 to be We know that |c n (t)| ≤ K n t E n (θ ) for all n ≥ 1, and will be interested in those n which are large enough for this bound to be effective, for the cases t = 0 and t = 1/20. It is easily computed that in both of these cases, the Euler bound K n t E n (θ ) does not even become smaller than 1 until n > 20, and when n = 26 (the smallest value of n for which we do not have access to the period-n points) it is of the order of 10 −5 for small t, so the Euler bound by itself would only permit a bound on σ 2 μ (g) which is accurate to around 1 decimal place. It is therefore crucial that we use the computed Taylor bounds in order to yield the high accuracy bound on σ 2 μ (g) given in Theorem 3, and in the proof of that result we use the Euler bounds only for n > 40.
Henceforth let Q = 40, M = 300, N = 400 (so that in particular Q ≤ M ≤ N , as was assumed throughout Sect. 7), and consider the two cases t = 0 and t = 1/20. We first evaluate the H 2 (D) norms of the monomial images L g,t (m k ) for 0 ≤ k ≤ N = 400. These norms are decreasing in k, and Appendix Table 3 contains the first few evaluations, for 0 ≤ k ≤ 10. Using these norms L g,t (m k ) we then evaluate, for 1 ≤ n ≤ M = 300, the upper computed approximation bounds α n,N ,+ (t) = α n,400,+ (t) defined (cf. (47)) by 17 α n,N , These bounds are decreasing in n; Appendix Table 4

A Tale of Five Tails: The Ingredients for Validating the Diffusion Coefficient
The following Proposition 5 gives rigorous bounds on the tails of five series appearing in the formula for σ 2 μ (g) derived in Corollary 1.
Next we require an estimate on the terms c n (0). From Cauchy's integral formula where p is the positively oriented circle of radius p centred at 0, we see that |c n (0)| ≤ 1 p max t∈ p |c n (t)|, and making the choice p = 1/20 gives |c n (0)| ≤ 20 max t∈ 1/20 |c n (t)|. Therefore and using the values in Appendix Table 6 we can evaluate the finite sum and the values in Appendix Table 6 give The values in Appendix Table 6  which is the required bound (62).

The Rigorous Bound on the Diffusion Coefficient
In the proof of Theorem 3 we shall make repeated use of the following simple lemma, in settings where A and B are quantities which cannot be computed precisely, but where a and b are computable approximations, and errors α and β can be derived. A, B, a, b ∈ R and α, β > 0 satisfy |A − a| ≤ α and |B − b| ≤ β, then

Lemma 5 If
and We can now justify the quality of our computed approximation to σ 2 μ (g) as follows: Theorem 3 If T is the Lanford map, μ is its absolutely continuous invariant probability measure, and g(x) = x 2 , then the diffusion coefficient σ 2 μ (g) can be approximated by σ 2 25 , which is derived using T -periodic points of period up to 25, so that Proof For economy of notation, let us write where for our purposes N will equal either 25 or ∞, so in particular Corollary 1(b) gives Our periodic orbit calculations (as described in Sect. 5) yield the following:  (91) Using (77) again we see that (92) Writing we use (91) and (92) to see that and hence so from (90) and (94) we deduce |σ 2 μ (g) − σ 2 25 | ≤ η 3 + η 7 ≤ 1.48 × 10 −18 , and the desired bound (79) follows.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix: Numerical Data for the Model Problem
Here we include data (presented in truncated form) for various quantities used in the computation of the diffusion coefficient σ 2 μ (g) for T the Lanford map and g(x) = x 2 .   Table 5 Upper computed Taylor bounds β M,+ n,N ,+ (t) for Lanford map transfer operator L g,t with t = 0, g(x) = x 2 , M = 300, N = 400, and disc D centred at c = 0.664, of radius = 0.87 Table 5 continued n β M,+ n,N ,+ (0)