## Abstract

We consider the random continued fractionwhere *s*_{n} are independent random variables with the same gamma distribution. Every realization of the sequence defines a Stieltjes function that can be expressed asfor some measure *σ* on the positive half-line. We study the convergence of the finite truncations of the continued fraction or, equivalently, of the diagonal Padé approximants of the function *S*. Using the Dyson–Schmidt method for an equivalent one-dimensional disordered system and the results of Marklof *et al.*, we obtain explicit formulae (in terms of modified Bessel functions) for the almost sure rate of convergence of these approximants, and for the almost sure distribution of their poles.

## 1. Introduction

Let be a sequence of positive real numbers and consider the analytic continued fraction(1.1)This continued fraction defines a *Stieltjes function*; it can be represented in integral form as(1.2)for some measure *σ* supported on the non-negative half-line such that the *moments*(1.3)exist. By an obvious use of the geometric series, every Stieltjes function can be expanded formally in powers of *t,*(1.4)Hence, *S* is the *moment-generating function* of the measure *σ*. It is a well-known fact of great practical importance that, given the first *n* of the moments, one may construct the rational function(1.5)whereThis truncation of the continued fraction (1.1) has a MacLaurin expansion whose *n*th partial sum agrees with that of the series (1.4). Hence, *S*_{n} is a diagonal (if *n* is odd) or near-diagonal (if *n* is even) *Padé approximant* of *S*.

Now suppose that *S*_{n} are independent positive random variables with the same distribution, say *μ*. We shall consider the following questions.

What are the almost sure analytic properties of these Stieltjes functions?

What is the almost sure leading asymptotic behaviour of the error as

*n*→∞?

These questions are of interest because Padé approximation is widely used in applied mathematics as a practical means of accelerating the convergence of the partial sums of series obtained by perturbation methods. As pointed out by Bender & Orszag (1978), the consideration of many particular cases where the *S*_{n} are *deterministic* reveals a wide range of large-*n* behaviours. Our motivation for studying the random case is to gain some insight into the asymptotic behaviour of Padé approximation in the ‘generic’ case.

Our study of diagonal Padé approximation reduces to aspects of the large-*n* behaviour of the denominators *Q*_{n}. For this reason, as in the deterministic case, the cornerstone of the analysis is the three-term recurrence relation(1.6)(The *P*_{n} satisfy the same recurrence relation, albeit with different initial conditions.) This recurrence relation makes a link between Padé approximation and a rich set of other mathematical entities, such as orthogonal polynomials, products of random matrices and discrete Schrödinger-like operators. By exploiting results that are well known in these related fields, one may obtain—for a very large class of distributions *μ* of the coefficients *S*_{n}—some partial answers to the questions stated earlier. Our contribution in the present paper is to elaborate the particular case where *μ* is the gamma distribution. More precisely, we obtain explicit formulae (in terms of Bessel functions) for the leading term in the asymptotic behaviour of the error of Padé approximation and for the asymptotic distribution of the poles, as well as the location of the essential spectrum of the measure *σ*.

In the remainder of this introductory section, we describe briefly the key ideas underlying the analysis. Then we summarize our main results in the form of a theorem.

### (a) The moment problem

The Stieltjes moment problem is, given a sequence , to determine whether or not there exists a measure *σ* such that equation (1.3) holds for every *n*. Historically, mathematical objects such as the analytic continued fraction (1.1), orthogonal polynomials and Padé approximants were introduced as tools in the study of this moment problem (Akhiezer 1961; Nikishin & Sorokin 1988; Simon 1998). Stieltjes (1894) showed that a necessary and sufficient condition for the *existence* of a measure *σ* with the prescribed moments is(1.7)He also showed that(1.8)is a necessary and sufficient condition for the *uniqueness* of the measure.

Let us assume that condition (1.7) holds and describe in very broad terms one way of ‘reconstructing’ *σ* from its moments (see Akhiezer (1961) and Nikishin & Sorokin (1988) for a detailed treatment).

Recall that *x*′ is a *point of increase* of the measure *σ* ifThe *spectrum* of *σ* is the set of its points of increase and will be denoted spec(*σ*). For *x*>0, we shall denote by *δ*_{x} the probability measure on whose only point of increase is *x*.

Given the *m*_{n}, we may define an inner product, say , on the space of polynomials as follows: if *p* and *q* are two polynomials with coefficients *p*_{i} and *q*_{i}, respectively, thenKnowing the moments, we may compute *s*_{n} and *S*_{n}. We remark that *P*_{2n} and *Q*_{2n} are polynomials of degree *n*−1 and *n,* respectively. Set(1.9)Then the fact that *S*_{2n} matches the moment-generating series to *O*(*t*^{2n}) implies that *ψ*_{n} is orthogonal to every polynomial of degree less than *n*, in the sense of the inner product . It follows (see Akhiezer 1961, ch. 1) that the roots of *ψ*_{n} are simple and lie in ; denote them byGaussian quadrature then defines a discrete measure(1.10)that converges weakly to a measure *σ* that solves the moment problem. In particular, the inner product coincides with the inner product in .

The moment problem can also be approached from the point of view of operator theory. The recurrence relation (1.6) for *Q*_{n} implies the following recurrence relation for *ψ*_{n}:whereThese numbers may be used to define a certain Jacobi operator, say , with a domain contained in the Hilbert space , as follows: first, we consider sequences with only finitely many terms and set(1.11)Given condition (1.8), it is then possible to extend this definition uniquely to obtain an essentially self-adjoint operator; we use the same symbol to refer to this extension. It may then be proved that the moments of the spectral measure of the operator are precisely the *m*_{n}, and so this spectral measure coincides with *σ* (see Nikishin & Sorokin 1988).

### (b) The density of states

Consider the finite-dimensional truncationof the operator . The spectrum of is the set of zeros of the polynomial *ψ*_{n} defined by equation (1.9). By the Chebyshev–Markov–Stieltjes theorem (see Nikishin & Sorokin 1988, §2.8), between any two zeros of *ψ*_{n}, there is a point of increase of *σ*; so we are led to studying the distribution of *λ*_{n,j}.

Define a measure *κ*_{n} on by(1.12)This measure is the normalized eigenvalue counting measure of the matrix . Indeed, we have

The normalized counting measure *κ*_{n} has a weak limit, say *κ*, as *n*→∞, and so there is a function *N*, called the *integrated density of states* of , defined byIf *κ* is absolutely continuous, one can also speak of the *density of states*, say *ϱ*, defined by(1.13)Although the measures *κ* and *σ* may be very different, their essential spectra are the same. In the context of Padé approximation, the integrated density of states describes the distribution of the poles of the approximants.

### (c) Krein's string

There is an interpretation, due to Krein, of the spectrum of the operator in terms of the characteristic frequencies of a vibrating string (Akhiezer 1961). Consider a weightless, infinite, perfectly elastic string, tied at one endpoint *x*=0, along which some beads are distributed. Let *s*_{2n} be the mass of the *n*th bead, and denote by (*x*_{n}, *y*_{n}) its position in the *xy*-plane. We assume that the *x*_{n} are fixed and given by the recurrence relationFor a string of uniform unit tension, the small vertical motion is then described by the discrete wave equation(1.14)To study the characteristic frequencies of the string, we setThen equation (1.14) reduces toandwhere *n*≥1. Comparing this with the definitions of and *ψ*_{n} given earlier, it is readily seen that

### (d) The complex Lyapunov exponent

Dyson (1953) developed a method for studying the characteristic frequencies of the one-dimensional disordered chain(1.15)Here *c*_{2n−1} and *c*_{2n−2} are the ratios of the elastic modulus of the *n*th spring and of the mass of the two particles attached to it. Disorder may be modelled in many ways; for instance, by assuming that the *c*_{n} are independent and identically distributed. The approach was later simplified by Schmidt (1957) and applied to the tight-binding Anderson model for a one-dimensional crystal with impurities.

Luck (1992) gave a very readable, well-motivated account of the Dyson–Schmidt approach; in brief, it builds on the intimate connection between second-order difference equations, continued fractions and Markov chains. For our purpose, it will be convenient to work with the random difference equation(1.16)where *t* is a parameter in , and is the branch of the square-root function defined on that returns a number with a non-negative real part. The following lemma is obvious.

*For every* ,*where u*_{n} *solves the difference equation* *(1.16)* *with u*_{−1}*=*0 *and u*_{0}=1.

The relevant continued fraction is(1.17)and we write(1.18)for its truncation. Let *u*_{−1} and *u*_{0} be complex random variables. Then equation (1.16) defines a sequence of general terms *u*_{n} by recurrence. The distribution *ν*_{μ} of the random variable *Z* is a stationary distribution for the Markov chain(1.19)(In the terminology of iterated random maps, the random variables *Z*_{n} and are, respectively, the *backward* and *forward iterates* associated with the continued fraction; when *u*_{−1}=0 and *u*_{0}=1, they have the same distribution, but their asymptotic behaviours are very different; see Diaconis & Freedman 1999.)

It follows that the growth of *u*_{n} may be quantified by means of the *complex Lyapunov exponent* defined by(1.20)where ln denotes the principal branch of the logarithm. Indeed, standard results from the theory of Markov chains imply that if has a unique stationary distribution, then(1.21)for almost every realization of ** s**, independently of the choice of (Breiman 1960; Furstenberg 1963; Meyn & Tweedie 1993). In particular, we have the formulaEquation (1.21) is also central to the study of the integrated density of states of the operator introduced earlier (Dyson 1953; Schmidt 1957; Luck 1992). We shall see that, under a very mild assumption on the distribution

*μ*of

*s*

_{n},

### (e) Furstenberg's theorem

In order to carry out this programme, we shall also make use of the connection between the Markov chain and the product of random matrices(1.22)where(1.23)The distribution *μ* from which the *s*_{n} are drawn induces, via equation (1.23), a distribution on the group of unimodular 2×2 matrices. The fundamental results of Furstenberg & Kesten (1960) and Furstenberg (1963), which are commonly referred to as ‘Furstenberg's theorem’, imply in particular that, under very mild assumptions on *μ*, there is a unique measure on the group of unimodular matrices that is invariant under the action of the random matrix (1.23). Furthermore, the number(1.24)which quantifies the growth of the product and is independent of the choice of matrix norm , may be shown to be *strictly positive*. These results are of great relevance to our problem for two reasons: firstly, any invariant measure yields a measure that is stationary for the Markov chain , and vice versa; secondly,Hence, we deduce at once the uniqueness of the measure , as well as the exponential growth of *u*_{n}.

*Let the s*_{n} *be the independent random variables in* *with a common distribution μ that has at least two points of increase*. *Let* *be the sequence defined by the recurrence* *(1.16)*. *Suppose that**for some ϵ*>0. *Then, for almost every realization of the sequence* ** s**,

*the following holds independently of the starting values u*

_{0}≠0

*and u*

_{−1}:

*for Lebesgue almost every*,

*Furthermore*,

*if we set*,

*then for Lebesgue almost every x*>0

*,*

The proof of this proposition is provided in the electronic supplementary material, appendix A.

### (f) The gamma distribution: statement of the main result

Using such machinery, we are able, for a very wide class of distributions *μ*, to deduce the almost sure exponential nature of the convergence of diagonal Padé approximation and also to deduce the almost sure singularity of the measure *σ*. A more quantitative study requires the calculation of the complex Lyapunov exponent, but there are very few known instances where it can be expressed in terms of familiar functions.

In his seminal paper on the disordered chain (1.15), Dyson studied in some detail the particular case where the *c*_{n} are independent and gamma distributed. Dyson found the invariant distribution of the continued fractionin the particular case where *t*>0; he then used analytic continuation to obtain an expansion for the complex Lyapunov exponent at *t*<0, and hence for the distribution of the characteristic frequencies. Dyson's continued fraction is not equivalent to ours (compare equations (1.15) and (1.14)), and so his analytical results do not transfer to our problem. However, the continued fraction (1.17) with independent gamma-distributed *s*_{n}, i.e. where(1.25)was examined also by Letac & Seshadri (1983); they obtained the probability distribution and found an explicit formula for the corresponding Lyapunov exponent for *t*>0. In a recent paper, we generalized this result by finding and the real part of the complex Lyapunov exponent for *every* complex *t* (see Marklof *et al.* in press); a straightforward extension of these calculations leads to the remarkably simple formula(1.26)In this expression, denotes differentiation with respect to *a*, and *K*_{a} is the modified Bessel function of the second kind (see electronic supplementary material, appendix B). The following summarizes the key results of this paper.

*Suppose that the s*_{n} *are independent draws from the gamma distribution with parameters a*>0 *and b*>0. *Then*, *for almost every realization of the sequence* ** s**,

*the following holds*:

*The density of states is given explicitly by the formula**and its absolutely continuous part is empty*.*For Lebesgue almost every*,

### (g) Relation to other work

As stated earlier, our focus in the present paper is on the performance of diagonal Padé approximation, viewed as a method of summing some random series. Related questions have been considered in the past, in different contexts. Foster & Pitcher (1974) studied the convergence of random *T*-fractions; these are continued fraction expansions which are in a one-to-one correspondence with the space of formal power series, but whose convergents are not Padé approximants. Foster & Pitcher (1974) showed that under very general conditions on the distribution of the coefficients, the difference between two successive convergents tends to zero exponentially fast, and that the exponent is twice the Lyapunov exponent associated with an infinite product of random matrices. Geronimo (1993) studied the random measure (on the unit circle) generated by a three-term recurrence relation with random identically distributed coefficients; he showed the positivity of the corresponding Lyapunov exponent and deduced that the random measure is singular with respect to the Lebesgue measure. This list is not exhaustive (see also Csordas *et al.* (1973) and Mannion (1993)).

The question of the nature of the measure *σ* has a counterpart in the theory of disordered systems which has been studied extensively in the context of Anderson localization. For example, the tight-binding Anderson model uses a discretized version of the Schrödinger equation with a potential that takes random identically distributed values at every point in a doubly infinite lattice. The resulting operator has a second-order finite-difference form like that of the operator —in which, more precisely, *h*_{n}=1 and *v*_{n} are independent and identically distributed—but it acts on sequences in . For a very wide choice of the distribution of the potential values *v*_{n}, the Lyapunov exponent of the discretized Schrödinger operator is strictly positive, so that, by Ishii's formula, the absolutely continuous spectrum is empty. A more refined study (see, for instance, Carmona & Lacroix (1990) and Pastur & Figotin (1992)) reveals that these operators have a *pure point* spectrum, i.e. the spectrum is the closure of the discrete spectrum. Such operators are said to exhibit the *localization property* because, for this type of spectrum, the generalized eigenfunctions decay exponentially fast as |*n*|→∞. The rigorous extension of such detailed results to the semi-infinite case would involve technicalities which are beyond the scope of the present paper.

Our analysis exploits a number of ideas and techniques found in these earlier studies. *We view our main contribution as that of exhibiting an interesting example of a class of random Stieltjes functions for which the leading behaviour of the error of diagonal Padé approximation and the density of states of the corresponding Jacobi operator are given explicitly in terms of special functions*.

The remainder of the paper is devoted to a detailed proof of theorem 1.3: the first statement follows immediately from Dyson's formula for the density of states, which is derived in §2. In §3, we deduce the singularity of the measure *σ* from the positivity of the real part of the Lyapunov exponent. In §4 we show that the error of diagonal Padé approximation is inversely proportional to the square of *u*_{n}; this yields the third statement in the theorem. Finally, in §5, we provide a numerical illustration of our results.

## 2. The formula for the density of states

So that we can use proposition 1.2, we shall henceforth suppose that the *s*_{n} are independent draws from a distribution *μ* on such that

*μ*has at least two points of increase.There exists

*ϵ*>0 such that

To avoid needless repetitions, we shall not mention these particular assumptions explicitly again in the statement of the intermediate results leading to theorem 2.4. We begin by relating the growth of to the complex Lyapunov exponent.

*For almost every realization of the sequence* ** s**,

*we have*,

*for Lebesgue almost every*,

*and*,

*for Lebesgue almost every*,

Let . By definition,Hence, by lemma 1.1,where *u*_{n} solves the difference equation (1.16) with . This yieldsThe first statement in the proposition then follows from the electronic supplementary material, appendix A, corollary 5.3. The proof of the second statement is identical. ▪

Next, we examine the implications of the lemma for the distribution of . By virtue of the recurrence relation satisfied by , we can writeLetand let . Then(2.1)

*Suppose that**Then for almost every realization of the sequence* ** s**,

*for Lebesgue almost every*,

Consider the real part in equation (2.1). ThenBy lemma 2.1, the first term on the right tends tothe second term tends to zero and, by the ergodic theorem, the third term tends to ▪

*Under the same assumption, for almost every realization of the sequence* ** s**,

*the sequence*

*has a weak limit*,

*say κ*,

*which is a probability measure on*.

*In particular*,

The proof is a specialization of that given by Goldsheid & Khoruzhenko (2005) in the more general case of a non-Hermitian Jacobi matrix. See electronic supplementary material, appendix C for details. ▪

*Let the s*_{n} *be independent random variables in* *with a common distribution μ that has at least two points of increase*. *Suppose also that**for some ϵ*>0. *Then, for almost every realization of the sequence* ** s**,

*for Lebesgue almost every*,

Let . ThenHence, we have the identity(2.2)Consider the imaginary part; the result then follows from lemma 2.1 and corollary 2.3. ▪

*Suppose that the s*_{n} *are independent and gamma distributed with parameters a and b*. *Then*, *for λ*>0*, we have the following formula for the density of states*:

Let *λ*>0 and setFor , we haveHence,and soWe deduce the formulae(2.3)and(2.4)The proposition, together with the last equation, then yieldsDifferentiating both sides with respect to *λ*, we findWe obtain the desired result by changing the order of differentiation on the right-hand side, and making use of the identity ▪

*Under the same assumption, for almost every realization of the sequence* ** s**,

By differentiating the identity (see Watson 1966, §13.73)with respect to *a*, we deduce that *ϱ* is strictly positive. ▪

## 3. Singularity of the spectrum

Corollary 2.6 implies in particular that the radius of convergence of the generating series of the moments is zero almost surely; in other words, the random Stieltjes functions that we have constructed are not analytic at the origin.

The problem of determining the nature of the spectrum is more delicate. Every measure may be decomposed into three disjoint parts: its absolutely continuous, singular continuous and discrete parts, denoted by , and , respectively. Ishii (1973) and Yoshioka (1973) showed that the spectrum of the absolutely continuous part is given by the formula(3.1)This result may be established by examining the resolvent of the Jacobi operator . The following result is essentially equivalent:

*Let the s*_{n} *be independent random variables in* *with a common distribution μ that has at least two points of increase*. *Assume that**for some ϵ*>0. *Then, for almost every realization of the sequence* ** s**,

By lemma 2.1 and electronic supplementary material, appendix A, corollary 5.5, for almost every realization of ** s**, for Lebesgue almost every , we have(3.2)Now, let

*η*>0 and consider the setBy equation (3.2), this set has Lebesgue measure zero almost surely. On the other hand, it follows from the Men'shov–Rademacher theorem (see Nikishin & Sorokin 1988, proposition 8.3) that for

*σ*almost every ,So, almost surely, for every

*σ measurable set A,*▪

## 4. The rate of convergence

Let be the sequence defined by the recurrence relation (1.16) with the starting values and . Also, denote by the shift operator on the space of complex sequences, i.e.In order to emphasize the dependence of the continued fraction (1.17) on the sequence ** s**, we shall sometimes write

*Z*(

**) and**

*s**Z*

_{n}(

**) instead of**

*s**Z*and

*Z*

_{n}. We have the following convenient representation of the error

Using the recurrence relations (1.6) and the identityit is straightforward to show (by induction on n) thatThen

Lemma 1.1 then yields the desired result. ▪

*Let the s*_{n} *be positive independent random variables with a common distribution μ that has at least two points of increase*. *Suppose also that there exists ϵ*>0 *such that**Then*, *for almost every realization of the sequence* ** s**,

*for Lebesgue almost every*

*,*

Let be fixed. We haveHence, by lemma 4.1,(4.1)Consider the almost sure limit of each of the three terms on the right-hand side of this equality as *n*→∞. The first term tends to zero. By proposition 1.2, the second term tends toFinally, it follows easily from the ergodic theorem that the third term tends to the same limit. The result follows from a standard argument involving the use of Fubini's theorem. ▪

## 5. A numerical illustration

Following the example of Dyson (1953), it is instructive to begin with an examination of the (deterministic) case whereSetThenwhereNote that *μ*^{∞} has only one point of increase, and so the hypothesis of proposition 3.1 is not satisfied. Indeed, for this choice of *μ*, the spectrum of the measure *σ* is absolutely continuous.

In this case, the continued fraction (1.17) reduces toand so the complex Lyapunov exponent isIn particular, an elementary calculation shows thatNext, let and denote by the gamma distribution with . As Dyson remarked, this distribution has mean 1 and variance 1/*a* and so we may, for large *a*, view it as a perturbation of . Indeed, using our explicit formula for together with the large-order expansions in Abramowitz & Stegun (1964), §9.7, we findLikewise settingAnd using the large-order expansions in Abramowitz & Stegun (1964), §9.3, we obtain, for This expansion breaks down at *λ*=4; as *a* increases, diverges to infinity there but tends to zero exponentially fast for *λ*>4.

The following computations were performed in multiple-precision floating-point arithmetic with the Maple software package; the eigenvalues and eigenvectors of the matrix _{n} were calculated using the `Eigenvals` function, which implements the QR algorithm. Figure 1 corresponds to the case where with *a*=8; it illustrates the convergence of the counting measure to the integrated density of states *N* as *n* increases—thus confirming the validity of our formula for *N*.

While the Lyapunov exponent and the density of states are non-random, the measure *σ* is random. We note that is absolutely continuous whereas, for every , almost every realization of *σ* is singular. Figure 2 shows the approximationof the integrated measure for particular realizations corresponding to *a*=8 and 64. Here is the discrete measure defined by the quadrature formula (1.10), and *n*=256.

## Acknowledgments

The authors gratefully acknowledge the support of the Engineering and Physical Sciences Research Council (United Kingdom) under grant GR/S87461/01 and of the Leverhulme Trust under a Philip Leverhulme Prize (J.M.). The second author would also like to thank Professor Alain Comtet for his hospitality during a visit to the Laboratoire de Physique Théorique et Modèles Statistiques, Université de Paris-Sud (Orsay), and for many useful discussions while some of this work was carried out.

## Footnotes

Electronic supplementary material is available at http://dx.doi.org/10.1098/rspa.2007.0014 or via http://www.journals.royalsoc.ac.uk.

- Received May 2, 2007.
- Accepted July 5, 2007.

- © 2007 The Royal Society