Stieltjes moment sequences for pattern-avoiding permutations

A small subset of combinatorial sequences have coefficients that can be represented as moments of a nonnegative measure on $[0, \infty)$. Such sequences are known as Stieltjes moment sequences. This article focuses on some classical sequences in enumerative combinatorics, denoted $Av(\mathcal{P})$, and counting permutations of $\{1, 2, \ldots, n \}$ that avoid some given pattern $\mathcal{P}$. For increasing patterns $\mathcal{P}=(12\ldots k)$, we recall that the corresponding sequences, $Av(123\ldots k)$, are Stieltjes moment sequences, and we explicitly find the underlying density function, either exactly or numerically, by using the Stieltjes inversion formula as a fundamental tool. We show that the generating functions of the sequences $\, Av(1234)$ and $\, Av(12345)$ correspond, up to simple rational functions, to an order-one linear differential operator acting on a classical modular form given as a pullback of a Gaussian $\, _2F_1$ hypergeometric function, respectively to an order-two linear differential operator acting on the square of a classical modular form given as a pullback of a $\, _2F_1$ hypergeometric function. We demonstrate that the density function for the Stieltjes moment sequence $Av(123\ldots k)$ is closely, but non-trivially, related to the density attached to the distance traveled by a walk in the plane with $k-1$ unit steps in random directions. Finally, we study the challenging case of the $Av(1324)$ sequence and give compelling numerical evidence that this too is a Stieltjes moment sequence. Accepting this, we show how rigorous lower bounds on the growth constant of this sequence can be constructed, which are stronger than existing bounds. A further unproven assumption leads to even better bounds, which can be extrapolated to give an estimate of the (unknown) growth constant.


Context and motivation
Many important distributions in probability theory have moments that happen to coincide with some classical counting sequences in combinatorics. For instance, the moments of the standard exponential distribution E(1) with rate parameter 1, are equal to a n = n!, and thus count the number of permutations of S n = {1, . . . , n}.
Similarly, the moments of the standard Gaussian (normal) distribution N (0, 1), satisfy, for all n 0, a 2n = 1 · 3 · 5 · · · (2n − 1), a 2n+1 = 0, so that a 2n counts the number of pairings of 2n elements, or, equivalently, the number of perfect matchings of the complete graph K 2n . Such connections between probability and combinatorics turn out to be useful in a number of contexts. For instance, Billingsley's proof of the Central Limit Theorem by the "method of moments" [17, p. 408-410] arrives at a limiting Gaussian distribution (2), for sums of independent random variables, precisely thanks to the property of the corresponding moments to enumerate pairings (3).
In random graph theory, the Poisson distribution plays an important role. When the mean is 1, its moments a n := E(X n ), where Pr(X = k) = e −1 k! , count the number of set partitions of a collection of n elements, giving rise to the sequence of Bell numbers [50, p. 109-110]: a n = n! × coefficient of x n in e e x −1 .
It could be argued that the Poisson distribution arises, in many cases, from the emergence of set partitions, in accordance with Equations (4) and (5). Random matrix theory also abounds with examples, the most famous case being Wigner's semi-circle law. Here the Catalan numbers, (C n ) n 0 = (1, 1, 2, 5, 14, 42, . . .), known to enumerate a variety of trees and lattice paths and given by the formula C n = 1 n + 1 2n n (6) appear, in the asymptotic limit, as (renormalized) expectations of powers of traces of random matrices in various ensembles. Their moment representation x n x(4 − x) dx then implies convergence of the spectrum of the random matrices under consideration to the semicircle law with density 1 2π x(4 − x).

The main problem
There are thus good reasons to develop methods aimed at explicitly solving the so-called moment problem [112,3], which considers a sequence a = (a n ) n 0 of real nonnegative numbers, and searches to express its general term a n as the integral a n = Γ x n dρ(x) for all n 0, for some support Γ ⊆ R, and for some probability measure ρ. In most cases one can write dρ(x) = µ(x) dx, where µ(x) is a nonnegative function known as the probability density function of ρ (in short, density), in which case the above equation becomes a n = Γ x n µ(x) dx.
When in addition the support Γ is a subset of the half-line [0, ∞), the problem (7) is classically called the Stieltjes moment problem, and the sequence a ≡ (a n ) n a Stieltjes moment sequence 1 . The same problem is called the Hausdorff moment problem when Γ = [0, 1] and the Hamburger moment problem when Γ = R, see [68].
In this article, we will not consider the Hamburger moment problem, but only the Stieltjes and the Hausdorff problems. Note that a probability measure is entirely determined by its moments (i) in the Stieltjes case if |a n | n! for all n [17, Thm. 30.1], (ii) in the Hausdorff case, and more generally in the Stieltjes case with compact support [13], so that the problem (7) has a unique solution in those cases.
There are several necessary and sufficient conditions that the sequence a must satisfy in order to be a Stieltjes moment sequence, or equivalently, for a density function µ(x) to solve (7) for some support Γ ⊆ [0, ∞]. To express these conditions, one classically attaches to the sequence a the (infinite) Hankel matrices H ∞ n (a) =      a n a n+1 a n+2 . . . a n+1 a n+2 a n+3 . . . a n+2 a n+3 a n+4 . . . . . . . . . . . . . . .

Examples
Geometric sequences (τ n ) n 0 with growth rate τ > 0 are Stieltjes moment sequences, as τ n is equal to the n-th moment of the distribution with support Γ = {τ } which is τ with probability 1; equivalently, one can take α 0 = 1, α 1 = τ and α k = 0 for k 2 in (d). However, other very simple sequences with rational generating function are not Stieltjes; for instance, the Fibonacci sequence a = (1, 1, 2, 3, 5, 8, . . .) with generating function 1/(1 − x − x 2 ), is such that the second leading principal minor ∆ 2 1 (a) of H ∞ 1 (a) equals 1 × 3 − 2 × 2 = −1 < 0, and thus by (b) it cannot be Stieltjes. It is trivially a moment sequence for a linear combination of Dirac measures, hence it is a Hamburger moment sequence (with support Γ included in R but not in [0, ∞)). More generally, Stieltjes moment sequences with rational generating functions are well-understood: for instance, if (a n ) n 0 has a rational generating function with simple poles only, then a n is a finite sum of the form j c j p n j for some complex c j , p j , the corresponding measure ρ(x) is a weighted sum of Dirac measures j c j δ(x − p j ), and one has that (a n ) n 0 is Stieltjes if and only if all c j and p j are nonnegative real numbers.
A useful refinement of Theorem 1 (see e.g., the references cited before its statement), which essentially excludes Stieltjes moment sequences whose generating functions are rational, is that: (a ) a ≡ (a n ) n 0 is a Stieltjes moment sequence with a representing measure ρ having infinite support Γ (containing a dense subset of R) ⇐⇒ Perhaps the simplest example of a Stieltjes moment sequence with irrational generating function is a n = τ n+1 /(n + 1) with τ > 0, which corresponds to the uniform distribution U(0, τ ), with Γ = [0, τ ] and µ(x) ≡ 1. Observe that the generating function ∞ n=0 a n x n = −ln (1 − τ x)/x is not only irrational but already transcendental (i.e., nonalgebraic). However, it is D-finite, i.e., it satisfies a linear differential equation with polynomial coefficients. Another basic example of Stieltjes sequence is n!, corresponding to the standard exponential distribution E(1) for which, according to (1), one can take Γ = [0, ∞) and µ(x) = exp(−x). This can also be seen by using (d) and Euler's 1746 continued fraction expansion [47, §21], with α 0 = 1 and α 2k−1 = α 2k = k for k 1. Once again, the generating function n 0 n!x n is D-finite, but not algebraic. More generally, Euler [47, §26] showed that for the sequence with general term c (c + 1) · · · (c + n − 1) one can take α 2n−1 = c+n−1 and α 2n = n, so that the sequence is again Stieltjes for any real c > 0, and its generating function is transcendental D-finite and satisfies the linear differential equa- Here, and in all that follows, D x denotes the usual derivation operator ∂ ∂x . A slightly more involved example is 4 n (n + 1)! n!/(2n + 3)! which is Stieltjes using the support Γ = [0, 1] and the algebraic density µ(x) = √ 1 − x/4 and whose generating function is D-finite transcendental, equal to the electronic journal of combinatorics 27(4) (2020), #P4.20 Here, and in all that follows 2 F 1 (a, b; c; x) denotes the classical Gauss hypergeometric function ∞ n=0 (a)n(b)n (c)n x n n! defined for any a, b, c ∈ Q, −c / ∈ N, where (x) n denotes the Pochhammer symbol (x) n = x (x + 1) · · · (x + n − 1) for n ∈ N.
One of the simplest Stieltjes moment sequences with an algebraic but irrational generating function is the aforementioned Catalan sequence (6), which has a plethora of combinatorial interpretations [119]. As we will see in Section 2.5, in order to show that the Catalan sequence is Stieltjes, one can take in (a) the support Γ = [0, 4] and the density µ(x) = 4−x x . Equivalently, one can show that in (b), all the leading principal minors ∆ n 0 (a) and ∆ n 1 (a) of H ∞ 0 (a) and H ∞ 1 (a) are equal to 1, a property that uniquely characterizes the Catalan sequence [119, Pb. A35(b)]. Equivalently, all the α i 's in the continued fraction expansion (d) are equal to 1.

Main objectives
In this article we are mainly interested in identifying combinatorial sequences that are Stieltjes moment sequences, and in finding an explicit description of the corresponding densities: either a closed formula, or a defining algebraic or differential (potentially nonlinear) equation, or a numerical approximation if none of those are available. We will mainly focus on sequences that occur in connection with pattern-avoiding permutations. The corresponding sequences, denoted Av(P), count permutations of S n = {1, 2, . . . , n} that avoid some given pattern P. In that context, the Catalan sequence emerges as the simplest case, corresponding to the enumeration of permutations avoiding any fixed pattern of length 3, such as (213) or (312) [113,73].
One reason for attempting to identify combinatorial sequences as Stieltjes (or Hausdorff) moment sequences is that such sequences are log-convex 2 . To see that a Stieltjes moment sequence is log-convex, it suffices to observe that for each n 1, the expression a n+1 a n−1 − a 2 n 0 is a minor of H ∞ 0 (a), so it is nonnegative by condition (c). Alternatively, one can deduce log-convexity directly from condition (a) as follows: observe that for any real numbers β and γ the expression β 2 a n−1 + 2βγ a n + γ 2 a n+1 = x n−1 (β + γx) 2 dρ(x) is nonnegative. Then setting β = a n and γ = −a n−1 yields a n+1 a n−1 − a 2 n 0. If the sequence a is positive, then its log-convexity implies that the ratios an a n−1 are lower bounds on the growth rate τ of the sequence 3 .
In this article, we will only deal with sequences having at most exponential growth. An additional characterization of the Stieltjes property [80,42] is available to study problems on such sequences, and is given in Theorem 2 below (taken from [42,Theorem 6.6]; the proof is based on [80,Corollary 1]). To state it, let us recall that a Nevanlinna function, (sometimes called a Pick function, an R function or a Herglotz function) is a complex analytic function on the open upper half-plane C + and has nonnegative imaginary part. A Nevanlinna function maps C + into itself, but is not necessarily injective or surjective.
Theorem 2. For a sequence of real numbers a = (a n ) n 0 , and a real number τ > 0, the following assertions are equivalent: (a) a is a Stieltjes moment sequence with exponential growth rate at most τ , (b) There exists a positive measure ρ on [0, τ ] such that, for each n 0, a n = τ 0 x n dρ(x), (c) The analytic continuation of the generating function f (z) = ∞ n=0 a n z n of a is a Nevanlinna function which is analytic and nonnegative on (−∞, 1/τ ).

Contributions and structure of the article
The paper is organized as follows. In Section 2 we begin by recalling a basic but extremely useful tool, the Stieltjes inversion formula. We review in §2.1 a first proof which is essentially the one given by Stieltjes himself in [120], and then give in §2.2 a second instructive proof based on a complex analytic approach, using the Hankel contour technique.
In §2.4 we discuss some structural consequences of the inversion formula, ensuring that algebraic properties of the generating function of moments are nicely (and algorithmically!) preserved by the density function. We illustrate this in several ways in §2.5 on the toy example of the Catalan sequence, and in §2.6 on a first example coming from the world of pattern-avoiding permutations of length greater than 3.
In Sections §3.5 and §3.7 we show that the density function for the Stieltjes moment sequence counting permutations that avoid the pattern (12 . . . k) is closely, but non-trivially, related to the density attached to the distance traveled by a walk in the plane with k − 1 unit steps in random directions.
In Section 4 we consider the much more challenging (and still unsolved!) case of permutations avoiding the pattern (1324). We provide compelling numerical evidence (but not a proof) that the corresponding counting sequence is a Stieltjes sequence. This is based on extensive enumerations given in [38], in which the first 50 terms of the generating function are found. Assuming that the sequence can indeed be expressed as a Stieltjes moment sequence, we numerically construct the density function, and obtain lower bounds to the growth constant that are better than existing rigorous bounds. By extrapolation, we make a rather precise conjecture as to the actual growth rate.

Related work
Stieltjes moment sequences. A basic tool is the so-called Stieltjes inversion formula, invented by Stieltjes in his pioneering work [120], see also the references cited in §2. It has been notably used to study the moment problem in many algebraic cases, that is for sequences whose generating function is an algebraic function, such as the Catalan sequence (6), and many variants, see the references cited at the end of §2.4. Various generalizations of the Catalan numbers are proved to be Stieltjes sequences e.g., in [96,79,86,85,58,76,66,77,78] to name just a few, using the inverse Mellin transform (a tool similar to the Stieltjes inversion) and Meijer G-functions.
At a higher level of complexity lie sequences whose generating function is D-finite but transcendental (that is, not algebraic). For instance, the sequence n! 3 is proved to be Stieltjes in [69,Eq. (64)], with an explicit density function involving Bessel functions [69,Eq. (65)]. Various interesting combinatorial sequences arise as binomial sums, and a natural question is to study if they are Stieltjes or not. For instance, the sequence g n = n k=0 2 k k n k 2 is proved to be Stieltjes in [136, p. 7]. The same paper proves that the Franel sequence f n = n k=0 n k 3 is Hamburger and the paper [77] conjectures that it is Stieltjes. In the same vein, the famous Apéry sequence A n = n k=0 n k 2 n+k k 2 has been conjectured by Sokal to be Stieltjes, and apparently proved to be so in an unpublished work by Edgar [116]. The sequence g n counts the moments of the distance from the origin of a 3-step random walk in the plane. More generally, Borwein et al. [20,21,22] studied the densities of uniform n-step random walks with unit steps in the plane, corresponding to the moment sequence W n (k) in (21). Among the aforementioned references, these papers by Borwein are the closest in spirit to the present article: they discover features of the corresponding densities, like modularity and representations via hypergeometric functions (for 3 and 4 steps), and perform a precise numerical study of their properties.
The tools used in [20,21,22] are again Meijer G-functions and inverse Mellin transforms. Note that Bessel functions emerge again, see also [7], a fact that we will also pop up in our article. Not all combinatorial sequences have D-finite generating functions; at an even higher level of complexity lie those whose generating function is D-algebraic (i.e., the solution of a nonlinear differential equation) or whose exponential generating function is D-algebraic, but non-D-finite.
For instance, the Bell sequence (B n ) n 0 defined as in (5) by n 0 B n x n /n! = exp(e x − 1) belongs to this latter class, as do the Euler and the Springer sequences (E n ) n 0 and (S n ) n 0 defined by n 0 E n x n /n! = sec x + tan x and n 0 S n x n /n! = 1/(cos x − sin x). The five integer sequences (B n ) n 0 , (E 2n ) n 0 , (E 2n+1 ) n 0 , (S 2n ) n 0 and (S 2n+1 ) n 0 have been proved to be Stieltjes moment sequences [98,76,116] and exponential versions of their generating functions are known to be Dalgebraic 4 . However, the full sequences (E n ) n 0 and (S n ) n 0 are not Stieltjes moment sequences [116] 5 . A complete picture, classifying combinatorial Stieltjes moment sequences 4 To be more precise, the five power series B n x n /n!, E 2n x 2n /(2n)!, E 2n+1 x 2n+1 /(2n + 1)!, S 2n x 2n /(2n)!, S 2n+1 x 2n+1 /(2n + 1)! are D-algebraic. 5 Interestingly, it is noted in [116] that (S n ) n 0 is a Hamburger moment sequence, but not (E n ) n 0 .
according to the (differential) transcendental nature of their generating functions is however still missing.
Longest increasing subsequences. The increasing subsequence problem covers a vast area, difficult to summarize in just one paragraph. It goes back at least to the 1930s, see for example [45], and was popularized in a wide-ranging article by Hammersley [60], in response to a query of Ulam [128]. Such subsequences are studied in many disciplines, including computer science [70,124,103], random matrix theory [39,104,6], representation theory of S n [110], combinatorics of Young tableaux [70], and physics [64,92,127] and [52,Chap. 10]. A recent book [107] is entirely devoted to the surprising and beautiful mathematics of the longest increasing subsequences. See also Stanley's comprehensive survey [118] of the vast literature on increasing subsequences in permutations.
Counting permutations whose longest increasing subsequences have length less that k is in close connection with counting permutations avoiding the pattern ( [107] mentioned above. In order to avoid possible confusions, we stress that counting consecutive patterns in permutations is a related but different area, in which an occurrence of a consecutive pattern in a permutation corresponds to a contiguous factor of the permutation. The corresponding generating functions are different from the ones for permutations avoiding the pattern (12 . . . k), and are quite well understood [67, §5] and [44].
General pattern avoidance. As mentioned before, avoiding the pattern (12 . . . k) corresponds to the longest increasing subsequences problem; avoiding more general patterns, or sets of patterns, is an even vaster topic. Whole books are dedicated to it [67,19] and many problems remain widely open, including in cases when the pattern is simple in appearance. This is already the case for the innocently looking pattern (1324), for which even the nature of the corresponding generating function remains unknown, as does the growth of the sequence [65,37,16,38]. Until very recently, it was believed by some experts that for any pattern, the corresponding generating function for avoiding permutations would be D-finite. This is known under the name of the Noonan-Zeilberger conjecture (1996). The conjecture was disproved by Garrabrant and Pak [54,93] using tools from complexity/computability theory; their striking result is that there exists a family of patterns which each have length 80, such that the corresponding generating function is not D-finite 6 . Note that this existential result is not yet complemented by any concrete example, ideally of a single avoided permutation, such as (1324). We might say that the current understanding is at a similar level as in (good old) times when the existence of transcendental numbers was proved, without being able to exhibit a single concrete transcendental number. Coming back to permutations avoiding the particular pattern (1324), there are compelling (but still empirical) arguments that the generating function is not D-finite and that its coefficients grow as (11.60 ± 0.01) n [38], the smallest rigorously proved interval containing the (exponential) growth constant being (10.27, 13.5) [16].

Transformation on sequences
Operation on densities Figure 1: The correspondence between basic transformations on sequences and operations on densities: here a n , u n , b n are respectively associated with densities µ(x), λ(x), ν(x).

Closure properties
The class of Stieltjes moment sequences S(Γ) with support Γ is a ring and is closed under forward and backward shift, under dilation and under division by the sequence (n + 1) n 0 . There are a few basic transformations on sequences that correspond to simple operations on densities. For instance, if b n = a n+1 is a sequence obtained by forward shift, then the density associated with b n is x · µ(x); similarly, the backward shift b n = a n−1 is associated with density function (µ(x)−µ(0))/x, but now suitable integrability conditions of µ at x ∼ 0 are needed. We list a few such correspondences in Fig. 1. We only stress the formal aspects, give the basic versions, and refrain from stating detailed validity conditions, as these are obvious consequences of change-of-variables formulae, partial integration, and similar elementary techniques. Many closure properties are proved by Bennett in [12, §2].

The moment problem and the Stieltjes inversion formula
In what follows, we restrict our attention to combinatorial sequences, that is, we assume that the coefficients a n of the sequence a are nonnegative integers, counting objects of size n in some fixed combinatorial family. We assume that the generating function f (z) = ∞ n=0 a n z n has radius of convergence z c ∈ (0, ∞), with z c = 1/τ , that f (z) is analytically continuable to the complex plane slit along [z c , ∞), and that f (z) = O(z −1 ) as z → ∞. Much of this section applies more generally, but to simplify matters, it is more convenient to consider this more restricted situation.
Assuming that a n is equal to the n-th moment of a probability measure with density µ, that is, assuming equation (7) holds, our aim is to solve the Stieltjes moment problem, that is, to obtain the density µ in terms of the generating function f of a. In other words, we need to solve the equation that expresses f in terms of µ, namely The basic tool used to solve the moment problem is the Stieltjes inversion formula, invented by Stieltjes [125,Chap. 2], and also Masson's article [82] and Simon's survey [114].

Via the Poisson kernel
We follow here the presentation in [90, Lecture 2]. If ρ is a probability measure with density µ with support Γ ⊆ R, then its Stieltjes transform (or, Cauchy transform) is the map G µ : where C + and C − denote respectively the upper and the lower half-planes. It is easy to show that G µ is analytic on C + , and moreover if the measure ρ is compactly supported (as assumed from the very beginning in our setting), with a support Γ (in our case, Γ = [0, τ ]) contained in the interval [−R, R], then where a n = Γ x n dρ(x) = Γ x n µ(x) dx is the n-th moment of ρ. (This follows, after multiplication by µ(x) followed by integration, from the uniformly convergent expansion 1/(z − x) = n 0 x n /z n+1 which holds for all |z| > R |x|.) Therefore, Equation (8) simply states: Starting from f , find µ such that The Stieltjes inversion formula is an effective way of solving (9), i.e. of recovering the density µ of the probability measure ρ from its Stieltjes transform G µ . Denoting this formula reads: the electronic journal of combinatorics 27(4) (2020), #P4.20 (Here, and below, " " stands for the operation of taking the imaginary part of a complex number.) The original proof of Stieltjes is based on the so-called Poisson kernel on the upper half plane defined by P (t) := 1 π t 2 + 2 . First, by definition, hence h (t) can be expressed as the convolution integral and the properties of the Poisson kernel permit one to conclude the proof of (11).
In the important case when G µ admits a continuous extension to C + ∪ I, for some interval I ⊆ R, the Stieltjes inversion formula simply reads: for all x ∈ I.
Putting things together, we conclude that, given a sequence a = (a n ) n 0 and a probability density µ(x), the following assertions are equivalent: • a n is the n-th moment of µ, i.e., a n = Γ x n µ(x) dx, for all n 0, • the generating function f (z) = n 0 a n z n is equal to in which case, moreover, the following integral representation holds:

Via the Hankel contour technique
The previous inversion process (13) can also be performed using Cauchy's coefficient formula in conjunction with special contours of integration known as Hankel contours. These contours are classical in complex analysis (e.g., to express the inverse of the Gamma function in the whole complex plane), in the inversion theory of integral transforms [40], and more recently in singularity analysis, see [49] and [50, Chap. VI]. Hankel contours come very close to the singularities then steer away: by design, they capture essential asymptotic information contained in the functions' singularities. The proof which follows is taken from [27], and does not appear to be published elsewhere in the rich literature on the moment problem, which is why we reproduce it here. To start with, Cauchy's coefficient formula provides an integral representation where the integration path C should encircle the origin and stay within the domain of analyticity of f (z). (For instance, any circle of radius < z c is suitable.) Due to the assumptions made on f (z), we can extend the contour C to be part of a large circle centered at the origin and with radius R, complemented by a Hankel-like contour that starts from R − i to z c − i , then winds around z c to the left, then continues from z c + i to R + i . As R → ∞, only the contribution of the Hankel part H of the contour survives, so that Since f (z) has real Taylor coefficients, it satisfies f (z) =f (z). This was originally true for |z| < z c and it survives in the split plane, by analytic continuation. Thus, if we let tend to 0, we obtain in the limit where f − (z), f + (z) are the lower and upper limits of f (z) as z → x ∈ R + from below and above, respectively. This is a typical process of contour integration, used here in a semi-classical way. By conjugacy, we also have since the contributions to the integral from the real part of f cancel out. By denoting we obtain the real integral representation a Then, under the change of variables t = 1/x, one again obtains the Stieltjes inversion formula (12), under the equivalent form

Stieltjes inversion formula for non-convergent series
In the case where the series f (z) = ∞ n=0 a n z n has radius of convergence 0, we cannot apply the Stieltjes inversion formula as described, as f (z) only makes sense as a formal power series, as the sum does not converge except at z = 0. Nonetheless, we mention a generalization of this method due to Hardy [61], which applies as long as there exists a constant C > 0 satisfying a n C n · (2n)! for all n > 0.
The idea is to introduce a new function which does converge for small s, and which extends to an analytic function on C \ R + . One can then prove that Finally the Stieltjes inversion formula can be applied to solve for µ.

Algebraic and algorithmic consequences of the Stieltjes inversion
An important consequence of the Stieltjes inversion formula is that the properties of (the generating function of) a sequence of moments are perfectly mirrored by the properties of the corresponding density function. More precisely, the following holds: Theorem 3. Assume as before that a = (a n ) n 0 is the moment sequence of the density function µ on a bounded domain Γ, i.e., a n = Γ x n µ(x) dx for all n 0. If the generating function f (z) = n 0 a n z n of a belongs (piecewise) to one of the following classes (i) algebraic, i.e., root of a polynomial equation P (z, f (z)) = 0, with P ∈ C[x, y], (ii) D-finite, i.e., solution of a linear ODE with polynomial coefficients, (iii) D-algebraic, i.e., solution of a nonlinear ODE with polynomial coefficients, then the same is true for the density function µ.
Moreover, one can effectively compute an (algebraic, resp. differential) equation satisfied by the density function µ starting from an equation for f . Instead of proving the general statement, we illustrate the proof on an example in the next subsection.
We expect that this theorem generalizes to some cases where the support is not bounded, perhaps to all sequences satisfying Hardy's condition (14). Note that the converse of (ii) holds true [10,31]. However, this is not the case for (i) and we do not believe it is the case for (iii).
For instance, concerning the converse of (i), if Γ = [0, 1] and µ( then the generating function F (z) of the sequence of moments a n = Γ x n µ(x)dx, is a transcendental function. See [94,108] for a study of the algebraicity of integrals of the form As for the converse of (iii), we do not have a counterexample for which the support Γ is bounded, but for unbounded support, we have the following counterexample: take the log-normal distribution µ(x) = exp(− ln(x) 2 /4)/x, which is clearly D-algebraic. Its moments over Γ = [0, ∞) are equal to a n = exp(n 2 ). The generating function f (z) of (a n ) n is not D-algebraic, since by a result of Maillet and Mahler [109,Eq. (2)], if f (z) were D-algebraic, then there would exist two positive constants K and C such that a n < K n! C for all n ∈ N, which is clearly not the case.
Another interesting example is the D-algebraic density µ( x , whose moment sequence e n = ∞ 0 x n µ(x)dx is the sequence (1, 1, 5, 61, 1385, 50521, . . .) of Euler's secant numbers (these are precisely the numbers E 2n mentioned on page 8). Although the exponential generating function n 0 e n z 2n (2n)! is sec(z) = 1 cos(z) , hence it is D-algebraic, the ordinary generating function f (z) = n 0 e n z n is not D-algebraic [24].
Finally, we give a potential counterexample to (iii) which has bounded support. The density µ(x) = ln(2 sin(x/2)) on the interval [0, 2π] is clearly D-algebraic. However, the corresponding moments a n = 2π 0 x n µ(x)dx are polynomials in π whose coefficients are linear combinations of odd zeta values ζ(2k + 1): If the generating series of (a n ) n were D-algebraic, this would contradict a commonly accepted number theoretic conjecture saying that the numbers π and ζ(2k + 1) are algebraically independent.

A basic, yet important example: The Catalan case
Here, we illustrate the above results and procedures on the case of Av (123), the permutations avoiding the pattern (123), which are counted by Catalan numbers C n = 1 n+1 2n n [113,73]. Let f (z) be the generating function of the Catalan sequence, f (z) = 1 + z + 2 z 2 + 5 z 3 + 14 z 4 + 42 z 5 + 132 z 6 + · · · 2.5.1 The first way: exploiting algebraicity It is well-known that f (z) is algebraic, being the root of the polynomial P (x, y) = xy 2 − y + 1. It follows that g(z) = f (1/z)/z is also algebraic, to be precise, it is a root of K(z, y) = 1 − zy + zy 2 .
Assume that C n is the n-th moment of a probability measure with density µ(x), so that As the exponential growth rate of C n is τ = 4, one may further assume the support of µ is Γ = [0, 4]. By the Stieltjes inversion formula, one has In this particular case, one can solve P by radicals, and obtain directly that We explain now a method that works more generally in the case when P cannot be solved by radicals. It yields, by a resultant computation, a polynomial equation satisfied by the density µ.
Writing g(x+iy) as A(x, y)+iB(x, y), where A and B are real functions, and denoting by ϕ(x) and ψ(x) the limits in y = 0 + of A and of B, we deduce 0 = K(z, g(z)) = K(x + iy, A(x, y) + iB(x, y)) and taking the limit y → 0 + yields K(x, ϕ(x) + iψ(x)) = 0. We thus have and by identification of real and imaginary parts, ϕ and φ satisfy In particular, ϕ(x) and of ψ(x) are algebraic; moreover, ψ(x) is a root of the resultant As µ(x) is nonnegative, the conclusion is that In other words, the previous procedure finds that Remark. Once deduced, this kind of equality can be proved algorithmically by the method of creative telescoping, see e.g., [5,71,25]; in our case, this method finds and proves that the integrand U (n, x) = x n 4−x x is a solution of the telescopic equation Integrating the last equality w.r.t. x between 0 and 4 implies a linear recurrence satisfied by f n = 4 0 U (n, x) dx: , the sequence f n /(2π) satisfies the same recurrence, and the same initial conditions, as the sequence (C n ). Therefore, they coincide, and this proves (16). Remark. Note that the Stieltjes inversion formula has allowed us to recover and to prove algorithmically the Marchenko-Pastur (or, free Poisson) density, arising in connection with the famous Wigner semicircle law in free probability and in random matrix theory [129,130]! The same approach allows one to prove that, for any s 2, the Fuss-Catalan numbers where the density function is a positive algebraic function of degree at most s(s−1)/2 [27], see also [97,84,9,79,86,87,88,89] for related results and generalizations; these can be obtained using the Stieltjes inversion formula. In all these references, the generating functions and the corresponding densities are algebraic functions. Remark. Catalan numbers count Dyck paths (among many other combinatorial objects). More generally, the previous method works mutatis mutandis for sequences that count lattice paths in the quarter plane whose stepset is contained in a half-plane. The generating function of such walks is known to be algebraic [8], hence the associated measure is algebraic as well. However, the support is generally not contained in [0, ∞), therefore the counting sequences are only Hamburger, and not Stieltjes, moment sequences. The simplest example is that of Motzkin paths, whose generating function 3]. The principal minor ∆ n 0 (a) is equal to 1 for all n but ∆ n 0 (a) is 1, 0, −1, −1, 0, 1 for n = 1, . . . , 6, repeating modulo 6 thereafter [2, Prop. 2], thus confirming that the Motzkin sequence is not Stieltjes, but only Hamburger. An alternative way of seeing this is via Flajolet's combinatorial continued fractions [48]. Indeed, the generating series for Motzkin paths has a Jacobi continued fraction [48,Prop. 5] instead of a Stieltjes continued fraction as required for a Stieltjes moment sequence by part (d) of Theorem 1. Some more examples are given in Appendix E.

The second way: exploiting D-finiteness
We can use the D-finiteness of the generating function f (z) = n 0 C n z n rather than its algebraicity. We give here the argument, since it will be used in the subsequent sections for other D-finite (but transcendental) generating functions. Recall that from the recurrence relation (n + 2) C n+1 − (4n + 2) C n = 0, for all n 0, the generating function f (z) is D-finite and it satisfies the linear differential equation It follows that g(z) = f (1/z)/z is also D-finite, and satisfies the differential equation: Recall that by the Stieltjes inversion formula, one has It follows, by taking the limit at y → 0 + , by taking the imaginary part, and by linearity, that ψ(x), and thus also the density µ(x), satisfy the homogeneous part of the previous differential equation, that is: Therefore, µ(x) is equal, up to a multiplicative constant λ, to 4−x x . The constant λ can be determined using The third way. There is an alternative way to find µ(x), without using the Stieltjes inversion formula. The method is also algorithmic, close to the one presented in [31]. The starting point is, once again, that the generating function f (t) = n 0 C n t n of the Catalan numbers satisfies a linear differential equation, Then, summation and differentiation imply that which is, after dividing through by t and reducing the pole order by Hermite reduction, Next, we integrate by parts, to reduce the derivative (1/(1 − xt)) to 1/(1 − xt): This provides an alternative proof for the differential equation (17), and the rest of the argument is similar. In the next subsection we give a more complex example, corresponding to pattern-avoiding permutations, which form the main objects studied in this article.

Av(1342)
The generating function of the sequence counting permutations of length n that avoid the pattern 1342 (https://oeis.org/A022558), starting is algebraic and is equal [18] to: Making the substitution x → 1 x in f (x) and multiplying by − 1 πx yields Since the imaginary part of the rational function is zero, the imaginary part of the first term is the density function. That is to say, the corresponding density is We plot µ(x) in Fig. 2, as well as the numerically constructed density function. To construct the density function numerically, we make a polynomial approximation µ(x) ≈ P (x) over the known range x ∈ [0, 8] with the property that the polynomial reproduces the initial moments of the density function -that is, just the coefficients of the sequence Av(1342), -constrained by P (8) = P (8) = 0, see §4.1 for details. Graphically the two curves are indistinguishable. In Fig. 3, we display the ratio between these two curves, which shows that the approximation is very good in the bulk of the distribution, but that it becomes less accurate as a ratio as the density approaches 0.
3 Increasing Av(123 · · · (k + 1)) sequences We apply the approach from the previous section, notably the Stieltjes inversion formula, to the explicit study of densities emerging from the classical "increasing subsequence problem", or equivalently, from the study of permutations that avoid an increasing pattern of the form (123 · · · (k + 1)). Later, in Section 4, we address one of the most interesting unsolved problems in the area of pattern-avoiding permutations, that of enumerating permutations avoiding the pattern (1324), by using a purely numerical approach to the Stieltjes inversion.
In random matrix theory, a classical result of Diaconis and Shashahani [39] states that, if U (k) represents the group of complex unitary matrices of size k, then This result was extended to the case n > k by Rains [104], who proved that From Rains's result, we can deduce that for any k, the counting sequence (f nk ) n∈N of Av(123 . . . k + 1) is a Stieltjes moment sequence. To see this, let X be a random variable with the same distribution as |Tr(M ) 2 | M ∈U (k) . Clearly X is supported on the subset [0, k 2 ] of R 0 ; indeed, for a k × k unitary matrix M , 0 |Tr(M )| 2 k 2 , as each entry of the matrix has modulus 1. Rains's result implies that f nk is equal to the n-th moment E(X n ) of the distribution of X. Hence (f nk ) n∈N is a Stieltjes moment sequence for any k. Moreover, the corresponding density function, that we will denote by µ k (x), is precisely the density function of |Tr(M ) 2 | M ∈U (k) 8 .
Equality (19) yields the following multiple integral representation for f nk : In particular, the generating function for Av(123 . . . (k + 1)) is D-finite, a result proved on the level of exponential generating functions and without using (20) by Gessel [55, p. 280] 9 . Moreover, Bergeron and Gascon [15] explicitly computed the differential equations corresponding to the exponential generating functions for k 11 (solving earlier conjectures made in [14]). Therefore, the results from the previous section imply that, for any k, the corresponding density µ k (x) is also a D-finite function.
In subsequent subsections we explicitly find this density. In the simple case k = 2, this corresponds to the Catalan numbers, already treated in §2.5. For k = 3, 4 we calculate the density function µ k (x) in terms of (pullbacks of) a rather special subset of 2 F 1 hypergeometric functions corresponding to classical modular forms.
Note the formal resemblance of (20) with the multiple integral which occurs in the theory of uniform random walk integrals in the plane, where at each step a unit step is taken in a random direction [20,21,22]. The integral (21) expresses the n-th moment of the distance to the origin after n steps. Of particular interest is the evaluation W n (1), which is equal to the expected distance after n steps. We will actually see that there is a nontrivial link between the generating functions of f nk and of W k (n), mirrored by an equally nontrivial link between the corresponding densities µ k (x) and p k (x). The relationship between classical modular forms and some particular 2 F 1 hypergeometric functions is discussed in a number of places in the literature (see for example [121], [81], [135]). One of the simplest illustrations of this intriguing relation can be seen in the identity [135,Eq. (74)] where E 4 (q) = 1 + 240 q + 2160 q 2 + · · · denotes the weight-four classical Eisenstein series, and where j(q) = q −1 + 744 + 196884 q + 21493760 q 2 + · · · is the j-invariant (of an elliptic curve). In our examples, the Hauptmodul 1728/j will be a rational function H(x) of a variable x. Modular forms (such as E 4 ) satisfy infinite order symmetries corresponding to the modular equations (see §3.2, or Appendix B and Appendix D). The pulled-back hypergeometric functions representing the classical modular forms 10 (like 2 F 1 1 12 , 5 12 ; 1; H(x) ) 9 Note that an integral identity due to Heine allows to prove, via the exponential generating function are annihilated by order-two linear differential operators. The corresponding nome q(x) (explicitly, the exponential of the ratio of two solutions of the order-two operator) is a globally bounded series, i.e., it can be recast into a series with integer coefficients after a rescaling x → N x.
The emergence of classical modular forms (closely associated with elliptic curves) in the density functions of Av(1234) and Av(12345) discussed below ( §3.1 and  §3.3), is a hint that the general Av(123 . . . (k + 1)) series could be related to integrable theory, an observation which is in agreement with the relation with random matrix problems, discussed above.

Av(1234)
The sequence Av(1234) admits nice closed forms ([55, p. 281], [118, p. 12]), Its generating function F 3 (x) starts 11 and is a solution of an order-three operator L 3 = L 1 ⊕ L 2 which is the direct sum 12 of an order-one operator L 1 with a simple rational solution R(x) and an order-two linear differential operator L 2 (recall that D x denotes the derivation d/dx): which can be shown to admit, as in [23], a pulled-back 2 F 1 hypergeometric function H(x) as solution. Consequently, in a neighborhood of x = 0, the series (23) can be written as the sum F 3 (x) = f 1 := R(x) + H(x), where: Note that the pulled-back 2 F 1 hypergeometric H(x) is not a Taylor series but a Laurent series. The Taylor series with positive integer coefficients (23) amounts to getting rid of the poles of the Laurent series for H(x). Note that the pulled-back 2 F 1 hypergeometric H(x) can also be rewritten using the identity as an order-one linear differential operator acting on a classical modular form since the 2 F 1 pulled-back hypergeometric function 2 F 1 3 4 , 3 4 ; 1; x is actually a classical modular form, see (B.1) in [1,Appendix B] and Theorem 3 in [106], see also §3.2. Equation (70) in Appendix B also underlies a x ↔ 1 9x involutive symmetry for the pulled-back hypergeometric function H(x) and thus also for the Av(1234) series (23).
The general theory previously displayed tells us that the density can be obtained from formula (10), no matter whether the series we consider is an algebraic series, a D-finite series (as in this case Av(1234)), a differentially algebraic series, or a not-evendifferentially algebraic series. Let us rewrite the Stieltjes inversion formula (10) giving the density as:

Naive approach
In the case of the Av(1234) series (23), the function θ(y) is the sum of 1/y · R(1/y) and 1/y · H(1/y), which reads: Since the imaginary part of a polynomial (such as y+5 6 ) evaluated at real values is zero, it is clear that finding the imaginary part of (29) gives just the imaginary part of the pulledback 2 F 1 hypergeometric function 1/y · H(1/y). However evaluating the imaginary part of a pulled-back 2 F 1 hypergeometric is not straightforward. This compels us to examine the function as a piecewise function.
The density µ 3 (y) is defined on the interval [0, 9]. We first note that the pullback in the pulled-back 2 F 1 hypergeometric function in (29), is small when y is small. Since we are dealing with 2 F 1 hypergeometric functions we need to know, for y ∈ [0, 9] when the pullback H = − The detailed piecewise analysis of this function and its imaginary part is given in Appendix A. The conclusion of this naive approach is that we find an expression which is the correct density, but divided by a factor of 3.

The magic of modularity
Actually, the origin of the factor 3 comes from the modular equation (73). More precisely, the following identity H(x) = H 2 (x) holds: where H 2 is the right-hand side of (30). This expresses the fact that the principal part of the generating function of Av(1234) corresponds to an order-one differential operator acting on a classical modular form.
Let us introduce G(y) := 1 y H(1/y) and G 2 (y) := 1 y H 2 (1/y). Naively, one could imagine that G(y) and G 2 (y) are equal for small y and possibly in some interval like [0,9]. In fact, this is not true: their series expansions around y = 0 satisfy 3 G(y) = G 2 (y). In the Stieltjes inversion formula, the density candidate gets multiplied by a factor of 3. Thus, the corresponding density becomes the correct one.

Less naive approach
Again, we need to perform a piecewise analysis of the slightly more complicated expression H 2 (y), as we did with H(y) in Appendix A. Not surprisingly, after performing this analysis, one gets the correct density. In Appendix C we sketch some arguments showing that this alternative expression H 2 (y) behaves, as expected, better than the simpler one using H(y). In Fig. 4 we plot this correct density function, as well as plotting the density function, numerically constructed as described in Sec. 4.1. Again, they are graphically indistinguishable.

Av(12345)
The positive integer coefficient series F 4 (x) for Av(12345) reads 13 : 1 + x + 2 x 2 + 6 x 3 + 24 x 4 + 119 x 5 + 694 x 6 + 4582 x 7 + 33324 x 8 + · · · (31) It is a solution of an order-four linear differential operator L 4 = L 1 ⊕ L 3 which is the direct sum of an order-one linear differential operator L 1 with a simple rational solution R(x) = 18 x 2 +10 x+1 12 x 3 and an order-three linear differential operator L 3 which can be seen, using van Hoeij's algorithm in [62] 14 , to be homomorphic (with ordertwo intertwiners) to the symmetric square of an order-two linear differential operator: which has a classical modular form solution that can be written as a pulled-back hypergeometric function. In Appendix D we show that this pulled-back hypergeometric function is a classical modular form. Appendix D also underlines another involutive symmetry x ↔ 1 64x on this new classical modular form and thus on the Av(12345) series (31).
Introducing F(x) = H(x) 2 , the square of this classical modular form one finds that the series F 4 (x) for Av(12345) can be written as the sum of the rational function R(x) (the solution of L 1 ) and an order-two linear differential operator acting on F(x) given by (34). Thus, the Av(12345) series (31) is actually the expansion of Alternatively, L 3 can also be seen as being homomorphic (with order-two intertwiners) to the symmetric square of another order two linear differential operator, with a pulledback hypergeometric solution: which emerges in the study of Domb numbers (counting 2n-step polygons on the diamond lattice, and also the moments of a 4-step random walk in two dimensions, cf. https: //oeis.org/A002895). In other words, the following identity holds: As we will see in the next subsection, the relation between the generating function F 4 (x) of Av(12345) and the Domb generating function is not fortuitous: it mirrors what happens on the level of the density functions, namely that the density function µ 4 (x) corresponding to Av(12345) is intimately related to the density function p 4 (x) of the distance travelled in 4 steps by a uniform random walk in the plane [21].

Density of Av(12345)
The Av(12345) density function µ 4 (x) is also the probability density function of the squared norm of the trace of a random 4 by 4 unitary matrix. We recall that the support of the measure is [0, 16]. The sequence (31) is a Stieltjes sequence, the density being given by the general formula (28). Again one has to change x → 1/x in (35) and take the imaginary part of the slightly involved hypergeometric result.
As far as finding the density µ 4 (x) is concerned (before undertaking the delicate piecewise analysis of taking the imaginary part) one finds, using (35) and Stieltjes inversion formula, that where: In this case the support of the measure µ 4 (x) is [0, 16] and the piecewise analysis will need to consider the intervals [0, 4] and [4,16]. The value 4 corresponds to the pullback 108 y (16 − y) (2 + y) 2 being equal to 1. We will not pursue these calculations in detail here, as they are similar to those in the previous section for Av(1234). We simply show in Fig. 5 a plot of the density function µ 4 (x) for Av(12345) obtained numerically, using the approach described in §4.1.
The linear differential operator annihilating the density µ 4 (x) is easily obtained from the pullback of the order-three operator L 3 in (32). It reads 15 : Similarly, µ 4 (x 2 ) is a solution of the linear differential operator

Relation between densities of Av(12345) and of 4-steps uniform random walks
An unexpectedly simple relation connects the density µ 4 (x) for Av(12345) and the density p 4 (x) of short uniform random walks studied in [21]. Recall that p 4 (x) denotes the density function of the distance travelled in 4 steps by a uniform random walk in Z 2 that starts at (0, 0) and consists of unit steps, each step being taken into a uniformly random direction. In [21,Thm. 4.7] it was shown that for x ∈ (2, 4), which satisfies the order-3 operator [21, Eq. (2.7)] Interestingly, the operator (38) satisfied by µ(x 2 ) is gauge equivalent (homomorphic) with A 4 , with a nontrivial order-2 intertwiner This shows that the density p 4 (x) and our density µ 4 (x 2 ), both with support [0, 4], are actually differentially related in a non-trivial way. More precisely, the operator (39) sends p 4 (x) to a solution of the operator (38) satisfied by µ 4 (x 2 ). Conversely, µ 4 (x 2 ) is mapped by a similar order-2 operator to a solution of the operator satisfied by p 4 (x). We can push a bit further the analysis made in [21]. The operator A 4 annihilating p 4 (x) is the symmetric square of which has a basis of solutions {y 1 Similarly, an explicit differential expression of order 2 in µ 4 (x 2 ) could also be written as a linear combination with constant coefficients of y 2 1 , y 2 2 and y 1 y 2 .

Av(123456): order-four operator homomorphic to its adjoint
The positive integer coefficient series F 5 (x) for Av(123456) reads 16 : 1 + x + 2 x 2 + 6 x 3 + 24 x 4 + 120 x 5 + 719 x 6 + 5003 x 7 + 39429 It is a solution of an order-five linear differential operator L 5 = L 1 ⊕ L 4 which is the direct sum of an order-one operator L 1 with a simple rational solution R(x) and an order-four linear differential operator L 4 . This is The linear differential operator annihilating the corresponding density µ 5 (x) can easily be obtained from the pullback of the order-four operator (42). Since we do not have any closed exact formula for its solutions 17 we cannot perform the kind of piecewise analysis we performed previously. Mezzarobba's algorithms [83] could in principle obtain (numerically) this density, by effective (numerical) analytic continuation. The order-four operator L 4 is MUM ("maximal unipotent monodromy") and is nontrivially homomorphic to its adjoint with an order-two intertwiner. Its exterior square has a rational solution of the form P (x)/(301 x 11 ), where P (x) is a polynomial with integer coefficients. The differential Galois group of the order-four operator L 4 is thus Sp(4, C), or a special subgroup of it [28] 18 .
Remark. A priori, we cannot exclude the fact that L 4 could be homomorphic to the symmetric cube of a second-order linear differential operator, or to a symmetric product of two second-order operators. Furthermore, it could also be, in principle, that these second-order operators admit classical modular forms as solutions (pullbacks of special 2 F 1 hypergeometric functions). However, these options can both be excluded by using some results from differential Galois theory [115], specifically from [100,Prop. 7,p. 50] for the symmetric cube case, and from [100,Prop. 10,p. 69] for the symmetric product case, see also [62, §3]. Indeed, if L 4 were either a symmetric cube or a symmetric product of order-two operators, then its symmetric square would contain a (direct) factor of order 3 or 1. This is ruled out by a factorization procedure which shows that the symmetric square of L 4 is (LCLM-)irreducible.
Still, we cannot exclude the fact that the solutions of L 4 (and in particular the generating function for Av(123456), minus the rational part (41)) could be written as an algebraic pullback of a 4 F 3 hypergeometric function.
We show in Fig. 6 a plot of the density function µ 5 (x) for Av(123456) obtained numerically, using the approaches described in §4.1. For k > 6 we find that our approximations are visually indistinguishable from the graph of e −x , so we do not display these plots separately.

Relation between densities of Av(123456) and of 5-steps uniform random walks
An unexpectedly simple relation connects the density µ 5 (x) for Av(123456) and the density p 5 (x) of short uniform random walks studied in [21]. Recall that p 5 (x) denotes the density function of the distance travelled in 5 steps by a uniform random walk in Z 2 that starts at (0, 0) and consists of unit steps, each step being taken into a uniformly random direction. In [21,Thm. 5.2] it was shown that for x ∈ (0, 5) \ {1, 3}, the density p 5 (x) satisfies the order-4 operator [21, Eq. (5.4)] Interestingly, the operator satisfied by our density µ 5 (x 2 ), is gauge equivalent (homomorphic) with A 5 , with a nontrivial order-3 intertwiner. This shows that the density p 5 (x) and our density µ 5 (x 2 ), both with support [0, 5], are actually differentially related in a non-trivial way. More precisely, the intertwiner sends p 5 (x) to a solution of the operator (43) satisfied by µ 5 (x 2 ). Conversely, µ 5 (x 2 ) is mapped by a similar order-2 operator to a solution of the operator satisfied by p 5 (x).

Av(1234567): order-five operator homomorphic to its adjoint
The positive integer coefficient series F 6 (x) for Av(1234567) reads 19 : 1 + x + 2 x 2 + 6 x 3 + 24 x 4 + 120 x 5 + 720 x 6 + 5039 x 7 + 40270 It is a solution of an order-six linear differential operator L 6 = L 1 ⊕ L 5 which is the direct sum of an order-one operator L 1 with a simple rational solution R(x) The linear differential operator L 5 is MUM and non-trivially homomorphic to its adjoint. The intertwiners have order 4. The symmetric square of L 5 has a rational solution of the form P (x)/x N , where P (x) is a polynomial. Its differential Galois group is thus a special subgroup of the orthogonal group SO(5, C): Once again, the operator satisfied by µ 6 (x 2 ) is gauge equivalent (homomorphic) with the operator A 6 satisfied by the density p 6 (x) in [21]. Therefore, p 6 (x) and our density µ 6 (x 2 ), both with support [0, 6], are actually differentially related in a non-trivial way.

Av(12345678): order-six operator homomorphic to its adjoint
The positive integer coefficient series F 7 (x) for Av(12345678) reads 20 : 1 + x + 2 x 2 + 6 x 3 + 24 x 4 + 120 x 5 + 720 x 6 + 5040 x 7 + 40319 x 8 + 362815 It is a solution of an order-seven linear differential operator L 7 = L 1 ⊕ L 6 which is the direct sum of an order-one operator L 1 with a simple rational solution R(x) and an order-six operator L 6 , The linear differential operator L 6 is non-trivially homomorphic to its adjoint. The intertwiners have order 4. Its differential Galois group is thus a special subgroup of the symplectic group Sp(6, C). As before, the operator satisfied by µ 7 (x 2 ) is gauge equivalent (homomorphic) with the operator A 7 satisfied by the density p 7 (x) in [21]. Therefore, p 7 (x) and our density µ 7 (x 2 ), both with support [0, 7], are actually differentially related in a non-trivial way.
From all the examples of the previous subsections, it is legitimate to expect that for any k, the generating function F k (x) of Av(1234. . .(k + 1)) is a series solution of a linear differential operator L k (x, D x ) of order k. The upper bound k can be proven in the spirit of [15,Prop. 1]. It also appears that L k is moreover the direct sum L k = L with a rational solution R k−2 (x) with a unique pole at x = 0, where the polynomial P k−2 (x) has degree k − 2 and is of the form Their coefficients are given in Table 1 for k 10. Note that they seem to have a very nice combinatorial interpretation: for fixed r, the sequence of coefficients ([x r ] P (x)) r+1 coincides with the sequence (Q r ( )) r+1 for some polynomial Q r (x) of degree 2r in x. For instance, when r = 1, this sequence is ( 2 + 1) 2 = (5, 10, 17, . . .), and when r = 2, this is 1 2 4 − 3 + 1 2 2 − + 3 3 = (18, 71, 198, . . .). Even more remarkably, the evaluation Q r (−d) of Q r (x) at negative integers x = −d appears to count permutations of length r+d that do contain an increasing subsequence of length d. In other words, the polynomials Q r (−d) are exactly the polynomials B d (d+r) of [41]. This surprising coincidence deserves a combinatorial explanation 21 .
The order-(k − 1) linear differential operator L (k) k−1 appears to be Fuchsian 22 , MUM and homomorphic to its adjoint. As a consequence of this, it has a differential Galois group which is a subgroup of the orthogonal group SO(k − 1, C) when k is even, and of the symplectic group Sp(k − 1, C) when k is odd.
The linear differential operator L (k) k−1 has a Laurent series solution k−1 (x) whose sum with the Laurent series expansion of the rational function R k−2 (x) is exactly the generating function F k (x) of Av(1234. . .(k + 1)): 21 Another intriguing property is that Q r (r) = r + 1 and Q r (r − 1) = r. 22 Its finite singularities are at 0, and at 1, 1 9 , . . . , 1 k 2 if k is odd and at 1 4 , . . . , 1 k 2 if k is even.
Therefore, the explicit calculation of the density function µ k (x) for the Stieltjes sequence Av(1234. . .(k + 1)) amounts to performing a piecewise analysis of the sum Because of the special form of R k−2 (x), the function 1/y · R k−2 (1/y) is a polynomial. The imaginary part of the evaluation of this polynomial for real values of y is zero. The density function for the Stieltjes sequence Av(1234. . .(k + 1)) is thus a piecewise D-finite function, and a solution of the order-(k − 1) linear differential operator Of course, this linear differential operator is again (non-trivially) homomorphic to its adjoint, and it has the same differential Galois group as L k−1 could be homomorphic to the (k−2)-th power of a second-order linear differential operator. However, we cannot exclude the fact that the solutions of L (k) k−1 , and in particular the generating function F k (x) of Av(12 . . . (k+1)), minus the rational part R k−2 (x) in (48), could be written as an algebraic pullback of a k−1 F k−2 hypergeometric function. If we were forced to make a bet, our guess would be that the solutions might belong to the world of (specializations of) multivariate hypergeometric functions.
In 1990 Gessel [55] had shown that the double Borel transform 24 of the generating function of Av(123 . . . (k + 1)) could be expressed as the determinant of a k × k matrix whose coefficients are Bessel functions. More precisely, recalling that the generating function of Av(123 . . . (k + 1)) is n f nk x n , with f nk as in (19) and (20), then the double Borel transform x n+k/2 n!(n + k)! .
Clearly, for fixed k, the power series I k (x) is D-finite, so an immediate consequence of Gessel's result is that Y k (x) is also D-finite. Bergeron and Gascon [15] showed that, moreover, Y k (x) can be expressed as homogeneous polynomials in the D-finite functions I 0 (x) and I 1 (x) 25 satisfying the order-two linear differential equations For instance, The underlying differential operators L 2 (x), L 3 (x), . . . , L 7 (x) for Y 2 (x), Y 3 (x), . . . , Y 7 (x) were derived Bergeron and Gascon in [15]. For instance, The solutions of the differential equations L k (y(x)) = 0 are not given by regular hypergeometric functions or solutions of Fuchsian linear differential operators, but rather by solutions of linear differential operators with an irregular singularity at ∞. The differential operators L k are again homomorphic to their adjoint. Moreover, since the operators (52) are homomorphic and the Y k are homogeneous polynomials in I 0 and I 1 , one deduces immediately that the L k 's are homomorphic to the k-th symmetric power of the order two operator L 1 = xD 2 x + D x − 1. It is known [32] that symmetric powers of a fixed second-order operator are related by a recurrence of order two with simple order-1 operators as coefficients: if M = D 2 x + 25 In Maple jargon, these are I 0 = BesselI(0, 2 √ x) and I 1 = BesselI(1, 2 √ x).
a(x)D x +b(x), and then its n-th symmetric power is given by the n-th term of the sequence of operators M 0 = 1, L 1 = D x , and for all k 1, Our operators L k are not symmetric powers, but are homomorphic to such symmetric powers. Consequently, we do not expect a recurrence as simple as the one on the M k 's. However, we remark that the L k 's still possess some remarkable features. For instance, By analogy with the emergence of these symmetric powers for L k , one could imagine that the L k 's associated with Av(12 . . . (k + 2)) might also be homomorphic to the symmetric (k − 1)-th power of an order 2 differential operator, which would hopefully correspond to a classical modular form.

A difficult case: Av(1324)
In this section, we investigate the density function for the most difficult Wilf class of length-4 pattern-avoiding permutations, namely Av(1324), whose generating function is unknown. The study of this (conjectured) density function is entirely numerical, based on the extensive exact enumerations of Conway, Guttmann and Zinn-Justin [38]. The series with positive integer coefficients for Av(1324) is 1 + x + 2 x 2 + 6 x 3 + 23 x 4 + 103 x 5 + 513 x 6 + 2762 x 7 + 15793 x 8 + . . . , and is known [38] up to order x 50 . The generating function is not known, and there are compelling arguments [38] that it is not D-finite. While it is known that the coefficients grow exponentially as λ n , the value of λ is not known. The best estimate is in [38] and is λ = 11.60 ± 0.01. The best rigorous bounds [16] are 10.27 < λ < 13.5, while the paper [36] gives the improved upper bound e π √ 2/3 ≈ 13.001954, but only under a certain conjecture about the number of inversions in a 1324-avoiding permutation.
The Hankel determinants are all positive and are monotonically increasing, which provides strong evidence for the conjecture that the series is a Stieltjes moment sequence.
Recall that if the sequence a is positive, then its log-convexity implies that the coefficient ratios an a n−1 are lower bounds on the growth rate µ of the sequence. In this way we obtain λ > 9.03.
In the case that a is a Stieltjes moment sequence, stronger lower bounds for µ can be calculated using a method first given by Haagerup, Haagerup and Ramirez-Solano [59], that we now summarize.
Using the coefficients a 0 , a 1 , . . . a n , one calculates the terms α 0 , α 1 , . . . α n in the continued fraction representation given in Theorem 1. It is easy to see that the coefficients of A(x) are nondecreasing in each α j . Hence A(x) is (coefficient-wise) bounded below by the generating function A n (x), defined by setting α n , α n+1 , α n+2 . . . , to 0. Therefore, the growth rate µ n of A n (x) is not greater than the growth rate µ of A(x). The growth rates µ 1 , µ 2 , . . . clearly form a non-decreasing sequence, and, since the coefficients of A n (x) are log-convex, µ n a n /a n−1 . It follows that this sequence µ 1 , µ 2 , . . . of lower bounds converges to the exponential growth rate µ of a.
If we assume further that the sequences α 0 , α 2 , α 4 . . . and α 1 , α 3 , α 5 . . . are nondecreasing, as we find empirically in many of the cases we consider, we can get stronger lower bounds for the growth rate by setting α n+1 , α n+3 . . . to α n−1 and α n+2 , α n+4 . . . to α n . For this sequence the exponential growth rate of the corresponding sequence a is ( √ α n + √ α n−1 ) 2 . By the method with which we constructed this bound, it is clear that Hence, the lower bounds ( √ α n + √ α n−1 ) 2 converge to the growth rate µ.
In particular, if α n α n+2 for each n and the limit lim n→∞ α n exists, then it is equal to µ/4. Extensive use of this result in studying the cogrowth sequences of various groups has been made in [59] and [43]. A plot of the known values of α n for Av(1324) is shown in Fig. 7  It appears that the coefficients α 2k and the coefficients α 2k+1 each form an increasing sequence which is roughly linear when plotted against n −2/3 . If this trend continues, the terms α n will certainly all be positive. If it could be proved that this is a Stieltjes moment sequence, these α values would imply a lower bound of λ > 10.302, which is an improvement on the best lower bound cited above. If we assume further that the sequences (α 2k ) k 0 and (α 2k+1 ) k 0 are both increasing, we get an improved lower bound λ > 10.607.
Extrapolating the two subsequences (α 2k ) k 0 and (α 2k+1 ) k 0 , it seems plausible that these both converge to the same constant c ≈ 2.9, which is consistent with the growth rate λ ≈ 11.60 predicted by Conway, Guttmann and Zinn-Justin [38].
From the numerical data, a histogram of the density function can be constructed, and this is shown in Fig. 8. We have no explanation for the little wiggle near the origin.

Numerics
We have used two different types of numerical construction of the density functions.
The histogram for Av(12 . . . k) (D-finite, discussed in §3.1) is constructed as follows: • take a random (k − 1) × (k − 1) matrix with entries which independently have standard complex normal distributions • orthogonalize it to create a random unitary matrix • compute the norm of the trace squared (by the Rains result [104] the n-th moment of this distribution is exactly the number of 12 . . . k-avoiding permutations of size n) • repeat the above steps 500 000 times and draw a histogram of the results, where each bar has width 0.02 Note that this procedure does not use the exact coefficients.
The polynomial approximations for the densities are constructed by computing the moments a 0 , a 1 , . . . , a n (that is, just the initial terms of the sequence), then approximating the density function by the unique polynomial P (t) of degree n + 2 over the known range [0, µ) with the following properties: • the initial moments of the distribution with density P (t) are a 0 , a 1 , . . . , a n • P (µ) = P (µ) = 0 For Av(1342), we have µ = 8. For Av(1324) we chose µ = 12, though for any value of µ between 11 and 13 the figure looks almost identical.

Final remarks
Remark 1. For each permutation π of length at most 4, we have proved or given strong evidence that the sequence Av n (π) n 0 , whose general term is the number of permutations of {1, . . . , n} that avoid the pattern π, is a Stieltjes moment sequence. Indeed, there is only one Wilf class of length 3, namely (123), considered in §2.5 and there are 3 Wilf classes of length 4, namely (1342), (1234) and (1324). We proved that the first three are Stieltjes moment sequences (in §2.5, §2.6 and §3.1) and presented numerical evidence (in §4.1) for the last one. Moreover, the same equally holds for any π of the form (12 . . . k), as showed in §3. This naturally suggests the further question: Open question: Is it the case that, for every permutation π, the sequence Av n (π) n 0 is a Stieltjes moment sequence?
We tested all known terms of these sequences and found they are consistent with being Stieltjes moment sequences. However, to our knowledge, substantial computations of the initial terms of these sequences have only been done in the cases discussed in our paper.
Remark 2. One of the referees suggested that the classes of triples and pairs of length-4 pattern-avoiding permutations enumerated by Albert et al. [4], and conjectured to be non-D-finite, should be similarly studied. These sequences are Av(4123, 4231, 4312), Av(4123, 4231), Av(4123, 4312) and Av(4231, 4321), and their first terms are given in the OEIS as A257562, A165542, A165545, A053617, respectively. Unfortunately, they all have Hankel determinants that become negative after a certain order, and so cannot be described as Stieltjes moment sequences.

Conclusion
We have shown how considering coefficients of combinatorial sequences as moments of a density function can be useful. Quite often in the literature the analysis of Stieltjes moment sequences with explicit densities is performed only for algebraic series which are slightly over-simplified degenerate cases.
We have studied here in some detail examples that are not algebraic series. For instance we have shown that Av(1234) is a Stieltjes sequence whose generating function corresponds, up to a simple rational function, to an order-one linear differential operator acting on a classical modular form represented as a pulled-back 2 F 1 hypergeometric function, and Av(12345) is a Stieltjes sequence whose generating function, up to a simple rational function, corresponds to an order-two linear differential operator acting on the square of a classical modular form represented as a pulled-back 2 F 1 hypergeometric function. The corresponding densities are of the same type, and do not involve the rational function. This scheme generalizes to all the Av(12345 . . . k) which are all Stieltjes moment sequences.
The linear differential operators annihilating such series are direct sums of an order-one operator with a simple rational solution where the denominator is just a power of x, and a linear differential operator homomorphic to its adjoint (thus corresponding to selected differential Galois groups). This last operator has a Laurent series solution, the generating function of the Stieltjes sequence corresponding to getting rid of the finite number of poles. The density is the solution of the x → 1/x pullback of this last operator, but finding the actual density requires delicate piecewise analysis which is difficult to perform when one does not have exact expressions for the generating function of the Stieltjes moment sequence.
We have shown that the density function for the Stieltjes moment sequence Av(12 . . . k) is closely, but non-trivially, related to the density attached to the distance traveled by a walk in the plane with k − 1 unit steps in random directions.
Finally, we have considered the challenging (and still unsolved!) case of Av(1324), and provided compelling numerical evidence that the corresponding counting sequence is a Stieltjes sequence. Assuming this, we obtained lower bounds to the growth constant that are better than existing rigorous bounds.
wishes to thank the ARC Centre of Excellence for Mathematical and Statistical Frontiers (ACEMS) for support.

A Density of Av(1234): Piecewise analysis
Discussion of the evaluation of 2 F 1 in Maple As far as the evaluation of a Gauss hypergeometric series like 2 F 1 (−1/4, 3/4; 1; H) goes, the hypergeometric series converges for |H| < 1, and 2 F 1 (a, b; c, H) is then defined for |H| 1 by analytic continuation. The point z = 1 is a branch point, and the interval [1, ∞] is the branch cut, with the evaluation on this cut taking the limit from below the cut. We need to look at [63,91,11]. The hypergeometric function 2 F 1 (a, b; c; H) can be rewritten 26 For H > 1, that the first term on the RHS of (60) corresponds to positive real values while the second term on the RHS of (60) corresponds to purely imaginary 27 values: Since the imaginary part of a polynomial (such as y+5 6 ) evaluated at real values is zero, the imaginary part of (29) is the imaginary part of the pulled-back hypergeometric term in (29) which reads: • For y ∈ [0, 1] (i.e. H ∈ [−∞, 0]), and 2 F 1 (−1/4, 3/4; 1; H) is given by (56) and is real and positive), the imaginary part of (29) (y − 1) 1/4 (y − 9) 3/4 can be rewritten, using (56), as the imaginary part of which reads, since y 2 + 18 y − 27 < 0: • For y ∈ [1,9] the imaginary part of (29) is thus the imaginary part of (y − 1) 1/4 (y − 9) 3/4 This can be rewritten, using (60), as the imaginary part of Note that the first term in (66) is just a complex number because of the factor (y−9) 3/4 (all the other factors are real and positive), while the second term in (66) is a complex number because of the factor (y − 9) 3/4 combined with the pure imaginary factor (1 − H) 1/2 (all the other factors are real and positive). After simplification, this can be written as This should give an imaginary part which reads, for y ∈ [1,9], B Av(1234) as a derivative of a classical modular form First note that The fact that 2 F 1 3 4 , 3 4 ; 1; is a classical modular form gives rise to this simple but infinite-order transformation.
For the non-classical modular form emerging naturally for Av(1234) a more complicated, but analogous, transformation exists. It is: where: The relation between these two Hauptmoduls is the modular equation: which is a representation of τ → 3 τ, where τ is the ratio of the two periods of the underlying elliptic function.

C Cut analysis for H 2 (y)
We have the formula f (x) = f 1 (x) := R(x) + H(x) near x = 0. We cannot directly apply the Stieltjes inversion formula to this, because is not analytic on C \ R >0 . Indeed H(x) has a cut wherever . These cuts are shown in Fig. 9. We can therefore only deduce that f 1 (x) coincides with the desired analytic extension f (x) inside the ring in Fig. 9. Alternatively we can use the following formula which can be derived from (71): This function is also not analytic on C\R >0 , as it has a cut wherever . These cuts are shown in Fig. 10. We can therefore only deduce that f 2 (x) coincides with the desired analytic extension f (x) inside the ring in Fig. 10. Nonetheless, this is an improvement on f 1 (x), as the ring in Fig. 10 contains the cut [1/9, 1]. As a consequence, applying the Stieltjes inversion formula directly to f (x) = R(x) + H 2 (x) yields the correct expression for the density on the interval [1,9].

D The classical modular form emerging in Av(12345)
Introducing U (p) 2 the linear operator U 2 given by (33) pulled-back by 64/x: One has the relations in (74) and in G(x) = F (x) √ x given by (34): Consequently the classical modular form (the square of which is F (x)) is which can also be written: The relation between the two (Hauptmodul) pullbacks is the (fundamental) modular equation which is a representation of τ → 2 τ . Note that one could have seen directly that L 3 , the order-three linear differential operator (32), is non-trivially homomorphic to its pullback by x → 1 64x .
E Densities for (Hamburger) moment sequences of walks Fig. 11 displays Hamburger moment sequences (not Stieltjes!), coming from the enumeration of walks restricted to the quarter plane. The context is the following: given a certain model (i.e., set of allowed steps), one forms the sequence (a n ) n 0 whose n-th term counts the walks confined to the quarter plane N 2 , starting at the origin, and consisting of exactly n steps. In general, the generating function F (z) = n a n z n is not algebraic, and not even D-finite. However, in Fig. 11 we consider only models leading to algebraic generating functions. These are basically of two types: either the support of the model is contained in a half-plane, e.g. for the first model (A001405), in which case the generating function is known to be algebraic [8]; or the support has a more general shape, and then algebraicity is not obvious, such as for 11th model (A001006) [30,Prop. 9], or for the last model (A151323) [30,Prop. 15], see also [26].
In each case we have computed the density function using the Stieltjes inversion formula. We describe below this process for the sixth model d s d E , as in this case the density involves a Dirac delta function, which has not appeared in our previous examples. Starting with the generating function for the model, we write G µ (z) = 1 z F ( 1 z ), as then the density µ(x) should be given by For x = 3, this yields the correct formula, however, G µ (x) has a pole at x = 3, as the limit does not exist at this point. There are several ways to resolve this problem (all of which, of course, yield the same density), the simplest being to remove the pole from F (z) before applying the Stieltjes inversion formula. To do this we write as thenF (z) has no poles. Applying the Stieljes inversion yields We then simply observe that where δ is the Dirac delta function. It follows that the density function µ(x) is Actually, much more can be said about the 13 sequences displayed in Fig. 11. For any sequence a among them, the leading principal minors ∆ n 0 (a) and ∆ n 1 (a) have very nice expressions: • in all cases except the last one, ∆ n 0 (a) is equal to q ( n 2 ) , where q is 1, 2, 3, or 4; • in all cases except the last one, the quotient sequence u n := ∆ n 1 (a)/∆ n 0 (a) is a linearly recurrent sequence with constant coefficients. which no longer satisfies a recurrence with constant coefficients, but still (conjecturally) satisfies a linear recurrence with polynomial coefficients u n+2 = 2 u n+1 − (4n + 9)(4n + 7) (2n + 5)(2n + 3) u n .
In cases 1-12, all these assertions can be proved using the explicit form of the corresponding (Stieltjes, or Jacobi) continued fractions, which can be determined since the generating function of a is an algebraic function of degree 2. This approach is very classical, and is for instance explained by Krattenthaler in [72, §2.7] and [74, §5.4]. It was applied in similar walk enumeration contexts by Tamm [123], by Brualdi and Kirkland [33], and by Chang, Hu, Lei and Yeh [34], to name just a few. Going beyond algebraicity degree 2 appears to be quite a challenging task. For degree 3, namely for some sequences in the orbit of the Catalan-Fuss sequence 1 3n+1 3n+1 n , explicit evaluations of Hankel determinants have been possible due to an approach introduced by Tamm [123] and extended by Gessel and Xin [57], see also [46] and Theorem 31 in [74]. Some of these evaluations are already connected to the enumeration of alternating sign matrices presenting (vertical) symmetries.
For the evaluations (79) and (80) of the Hankel determinants occurring in case 13, we can again use the J-fraction with γ 0 = 3, γ j = 2 for j > 0 and β j = (4j−1)(4j+1) (2j−1)(2j+1) . Similarly to cases 1-12, the Hankel determinants can easily be determined from these continued fraction coefficients. To prove that this is indeed the correct continued fraction, we adapt a method applied implicitly by Euler [47, §21] to S-fractions, which is described in detail in [101, §2.1] and [117]. The first step in our case is to set A 1 (x) := A(x)/β 0 , and define A j+1 (x) recursively by Then it suffices to show that each A j (x) ∈ R[x], as this immediately yields the continued fraction form (81) for A(x) = β 0 A 1 (x). The next step is to define another sequence {B n (x)} n∈N by B 0 (x) = 1 and B j (x) = A j (x)B j−1 (x). Writing A j (x) = B j (x)/B j−1 (x), the recursion (82) simplifies to Finally, we guess and prove the exact form of each series B j (x). First, the guess is that B j (x) is the unique series in R[x] with constant term 1 satisfying − 2j(2j + 1 + 6(j + 1)x) x(1 − 6x)(1 + 2x) B j (x) + 2j − (8j + 6)x − 24(j + 1)x 2 x(1 − 6x)(1 + 2x) B j (x) + B j (x) = 0. (84) To see that this equation uniquely defines B j (x), we observe that it is equivalent to the coefficients b j,0 , b j,1 , . . . of B j (x) being determined by the recurrence b j,n+1 = 2(2n + 2j + 1)(n + j)b j,n + 12(n + j)(n + j − 1)b j,n−1 (2j + n)(n + 1) , with the initial conditions b j,0 = 1 and b j,−1 = 0. Finally, we prove by induction that (84) holds for all j. The base cases are trivial to prove using the exact, algebraic forms of B 0 (x) and B 1 (x), so we proceed to the inductive step. We assume that (84) holds for j and j − 1, and prove that it holds for j + 1. Using (83), this reduces to proving which is true as the right hand side satisfies the differential equation (84) for B j−1 (x) and has constant term 1.
Model, sequence and OEIS tag Generating function F (z) Density µ(x) Support Γ