Central Limit Theorems via Analytic Combinatorics in Several Variables

The field of analytic combinatorics is dedicated to the creation of effective techniques to study the large-scale behaviour of combinatorial objects. Although classical results in analytic combinatorics are mainly concerned with univariate generating functions, over the last two decades a theory of analytic combinatorics in several variables (ACSV) has been developed to study the asymptotic behaviour of multivariate sequences. In this work we survey ACSV from a probabilistic perspective, illustrating how its most advanced methods provide efficient algorithms to derive limit theorems, and comparing the results to past work deriving combinatorial limit theorems. Using the results of ACSV, we provide a SageMath package that can automatically compute (and rigorously verify) limit theorems for a large variety of combinatorial generating functions. To illustrate the techniques involved, we also establish explicit local central limit theorems for a family of combinatorial classes whose generating functions are linear in the variables tracking each parameter. Applications covered by this result include the distribution of cycles in certain restricted permutations (proving a limit theorem stated as a conjecture in recent work of Chung et al.), integer compositions, and $n$-colour compositions with varying restrictions and values tracked. Key to establishing these explicit results in arbitrary dimension is an interesting symbolic determinant, which we compute by conjecturing and then proving an appropriate $LU$-factorization. It is our hope that this work provides readers a blueprint to apply the powerful tools of ACSV to prove central limit theorems in their own work, making them more accessible to combinatorialists, probabilists, and those in adjacent fields.

Let (X n ) be a sequence of random variables.A limit theorem (or limit law ) for X n is an approximation of the cumulative distribution functions P(X n ≤ k) as n → ∞.A local limit theorem is an approximation of the exact probabilities P(X n = k) as n → ∞.A central limit theorem (CLT) or local central limit theorem (LCLT) compares these probabilities to a normal density function The classical CLT states that if A 1 , A 2 , . . . is a sequence of independent identically distributed random variables with an expected value µ and a finite variance σ 2 > 0 then the sequence of random variables converges in distribution, after rescaling, to the standard normal distribution, so that for all real x.Note that we rescale the sequence X n so that its limit is a fixed distribution instead of one that varies with n.
The classical CLT has a long history (see Section 1 below for a brief recap), and has been generalized to weaken its assumptions, give explicit rates of convergence, and work in more abstract settings.In ddimensions, a multivariate (L)CLT compares probabilities to multivariate normal density functions The classical multivariate CLT states that if (A n ) is a sequence of d-dimensional independent identically distributed random variables with an expected value vector µ and a positive definite covariance matrix Σ then the sequence of d-dimensional random variables converges in distribution to the multivariate normal distribution with density function ϕ 0,Σ (x).
In where c i,n denotes the number of objects in C with size n and parameter value i. Well-established combinatorial theories exist to derive algebraic, differential, or functional equations satisfied by the generating functions of many types of combinatorial classes (see, for instance, [21,27,47,48]).Thus, one often has a multivariate generating function C(z, t) encoded in some way, and wants to determine the limiting behaviour of the d-dimensional array (c i,n ) i∈N d as a function of i when n → ∞.
Strikingly, the theory of analytic combinatorics shows how central limit theorems for many combinatorial parameters can be derived directly from a study of the singular behaviour of C(z, t).Flajolet and Sedgewick [21, Chapter IX] gives a detailed introduction, and Hwang [32] lists more than 25 limit laws proven by the late Philippe Flajolet using this framework, including patterns in words, cost analyses of various sorting algorithms, parameters in different kinds of trees, polynomials with restricted coefficients, and ball-urn models.
The classical theory of analytic combinatorics is chiefly concerned with the derivation of asymptotics of univariate generating functions.In typical cases, one applies transfer theorems to deduce the asymptotic behaviour of a generating function's coefficients from series expansions of the generating function near its singularities.For many classes of generating functions, it is essentially algorithmic (on examples occurring in practice) to compute the singularities determining asymptotics, apply transfer theorems, and find asymptotic behaviour.General multivariate results in analytic combinatorics can be traced back to work of Bender, Richmond, and collaborators in the 1980s (see Section 1 below for more details).Although such results apply to many applications in combinatorics and the analysis of algorithms, they capture only one type of multivariate singularity and require one to verify certain properties of multivariate generating functions (such as grouping singularities by moduli), which can be computationally expensive.
These issues led to the development of analytic combinatorics in several variables (ACSV) [41,44] starting in the early 2000s.Following the framework of univariate analytic combinatorics, the goal is to derive effective methods to take an encoding of a multivariate generating function and return asymptotics or limit theorems of its series coefficients.Most results in ACSV focus on sequences with multivariate rational generating functions: such sequences appear naturally in many applications, and many functions with more complicated singular behaviour can be encoded by rational functions in a higher number of variables (for instance, every algebraic function in d variables can be encoded as an explicit subseries of a rational function in 2d variables [12]).The crux of this paper is to illustrate both the power of ACSV for proving limit theorems and the (in many cases completely) automatic nature of the required computations.

Paper Structure
In addition to summarizing the existing literature on combinatorial limit theorems, and clarifying how the methods of ACSV build on and modify these results, we hope to illustrate how explicit and algorithmic modern methods for verifying central limit theorems are.To this end, we provide easy to use theorems2 , employing the ACSV template which can be applied by researchers in their own works.In addition to illustrating these methods on a collection of applications, including resolving a conjecture stated in Chung et al. [9] on the behaviour of cycles in a class of restricted permutations3 , we provide a SageMath package that can rigorously compute and prove central limit theorems for a large variety of combinatorial classes tracking parameters with multivariate rational generating functions, available at https://github.com/ACSVMath/Limit-Theorems-For-Combinatorial-ParametersSection 1 gives a brief history of central limit theorems, paying particular attention to the use of analytic methods.After surveying results in a probabilistic context, we describe the development of CLTs in combinatorics.Section 2 then gives an overview of the methods of ACSV, illustrating how multivariate LCLTs follow from the general theory.Section 3 illustrates how these results lead to a general limit theorem with explicit constants holding for a number of combinatorial objects, and Section 4 proves the explicit constants appearing in this limit theorem using a guess-and-check method to compute a symbolic determinant through LU -factorization.

Analytic Methods for Central Limit Theorems
The application of analytic methods to prove limit theorems in probability and combinatorics has a long and wide-ranging history.It is impossible to give a full account of such a vast topic in this space, so our presentation is a broad overview tailored to the context of our work.

Probabilistic CLTs
The proto-history of the CLT goes back at least as far 4 as a 1733 offprint Approximatio ad summam terminorum binomii (a + b) n in seriem expansi of Abraham de Movire (printed in English in the 1738 edition of his seminal Doctrine of Chances [11]).Motivated by the computation of explicit bounds for the Law of Large Numbers, de Moivre used Stirling's approximation for n! (which de Moivre independently approximated around the same time as Stirling) to deduce5  Perhaps the first systematic uses of analytic methods to derive CLT-like results are due to Laplace.Laplace's approach to the CLT built off of his ground-breaking work [34] on the approximation of parameterized integrals in the 1770s.Still influential to this day, Laplace's method (in its classical form) is used to approximate integrals of the form b a A(x)e −nϕ(x) dx where A and ϕ are analytic functions with ϕ minimized over [a, b] at a unique point c ∈ (a, b) such that ϕ has a Taylor expansion at x = c that begins with a quadratic term.In an 1810 memoir, Laplace [35] introduced the concept of characteristic functions by making a change of variable t = e ix in integral representations for the sum of certain independent and identically distributed random variables.He then (formally) used his method to asymptotically approximate the probability that this sum lies between two factors at √ n-scale by integrating a normal density.Although not working with a modern standard of rigor, Laplace's techniques continue to be influential to this day.
Many famous mathematicians in the nineteenth century worked on topics related to central limit theorems, including Gauss [25] (who derived the now sometimes-eponymous Gaussian function in the context of probabilistic error analysis), Poisson [45] (who stated a CLT for a normalized sum of random variables, and gave explicit conditions on characteristic functions for a CLT to hold), Dirichlet (who gave proofs with more rigor, including correct truncation bounds to prove errors arising in Laplace's method go to zero), and Cauchy [8] (who gave an updated proof of the CLT using characteristic functions).As pointed out by Fischer [19], it is also instructive from a historical perspective to reflect that Cauchy and Dirichlet were working at the time period when mathematics and probability were starting to move away from a discipline chiefly concerned by modelling observations of the physical world to a more abstract logic-based subject.A Russian school, involving mathematicians such as Hausdorff, Chebyshev, and Markov, also developed CLT-like results in the late nineteenth and early twentieth century using techniques such as the method of moments and moment generating functions.The first "modern" treatment of the CLT (as a general mathematical result not dedicated to specific applications or to illustrate analytic methods) is often considered to be work of Lyapunov [39,40] around the turn of the twentieth century.The term central limit theorem was likely coined by Pólya [46] (who studied various aspects of CLTs, and when sequences of distribution functions converge) in 1920.Lévy [38] and Feller [18] gave necessary and sufficient conditions for the CLT to hold in 1935.
Modern probability theory texts typically prove the CLT using the continuity theorem developed by Lévy [36,37] in the 1920s, which states that if the characteristic functions φ n (t) = E e itXn of a sequence (X n ) of random variables converge pointwise to a continuous function φ(t) then φ is the characteristic function of some random variable X and the sequence (X n ) converges in distribution to X (see [7,Theorem 26.3] or [16,Theorem 3.3.6]for proofs).If A n is a sequence of independent random variables with a common characteristic function φ(t) then the characteristic function of Assuming that the A j have a mean 0 and a finite variance σ > 0, it is possible to get an expansion as t → 0, so the characteristic function of for each t as n → ∞.The CLT in the mean zero case then follows from the continuity theorem, and the result for a general mean follows from a shift of the X n .A complete proof of the general Lindeberg-Feller CLT can be found in [7,Theorem 27.2] and [16, Theorem 3.4.5].Under slightly stronger conditions (such as finiteness of the third moments of the A j ) the Berry-Esseen theorem [6,17] describes the rate of convergence between the cumulative distribution functions of the X n and the cumulative distribution function of the normal distribution (see [16,Theorem 3.4.9]for a modern presentation).Multivariate limit theorems, including the multivariate CLT, can be established from univariate results through the classical Cramér-Wold theorem [10]: a sequence (X n ) of d-dimensional random variables converges in distribution to a random variable Y if and only if the sequence of univariate random variables θ • X n converges to θ • Y for all θ ∈ R d (see [7,Theorem 29.4] for a proof).

Combinatorial CLTs
As described above, in combinatorial contexts one has a multivariate sequence defined implicitly by a multivariate generating function tracking some parameter χ, and wishes to extract a limit theorem for the random variables (X n ) taking the values i ∈ N d with probabilities In the 1-dimensional case, a direct computation verifies where C z denotes the partial derivative of C with respect to z and the operator [t n ] extracts the coefficient of t n from a power series expansion.Similarly, the kth moment of X n can be expressed in terms of the first k derivatives (with respect to z) of C(z, t), and limit theorems for X n can be established by extracting coefficients from these expressions.The theory of analytic combinatorics shows how, instead of explicitly computing coefficients, their behaviour can be approximated -to a degree often sufficient to establish limit theorems -implicitly by studying the analytic behaviour of C(z, t) in t when z ≈ 1 (for CLTs) or when |z| ≈ 1 (for LCLTs).
A similar approach can be used for higher-dimensional parameters, except that the translation of the analytic behaviour of C to the limit behaviour of c i,n is more delicate.Classical treatments typically proceed by requiring that the coefficient slices encoding the objects of size n are approximated by quasi-powers.We briefly describe some of these results, and illustrate the methods involved in their proof, to compare with the limit theorems derived using the framework of ACSV in Section 2.
Theorem 1 (Bender and Richmond [5,Theorem 1]).Suppose that ϕ n (z) ∼ f (n)g(z)λ(z) n uniformly in a neighbourhood of z = 1 where g(z) is uniformly continuous and λ(z) has a quadratic Taylor series expansion with error term O(( |x k − 1|) 3 ).If the Hessian matrix H of log λ(e s ) at s = 0 is non-singular then where p n (i) is the scaled coefficient defined in (2) above, m is the gradient of log λ(e s ) at s = 0, and the inequalities are taken coordinate-wise.
Theorem 1 and its corollaries below were first established in the 1-dimensional case by Bender [4].Generalizations of Theorem 1 in the 1-dimensional case were given by Hwang [29,30,31] and Heuberger and Kropf [28] (see also the treatment in Flajolet and Sedgewick [21,Chapter IX]).
Proof Sketch.The random variable X k with probability distribution function p n has characteristic function .
Shifting X k by its mean nm and scaling by 1/ √ n gives a new random variable with a characteristic function ] for all fixed s.The 1-dimensional random variable s • X n thus satisfies a CLT for any s by the continuity theorem, meaning X n satisfies a multivariate CLT by the Cramér-Wold theorem.
The key insight for combinatorial CLTs is that it is often possible to prove that ϕ n satisfies the assumptions of Theorem 1 directly from the analytic behaviour of C(z, t) near some of its singularities.
Corollary 2 (Bender and Richmond [5, Corollary 1]).Let C(e s , t) be the generating function described in (1).Suppose that there exists ϵ, δ > 0, a number q ∈ Q \ {−1, −2 . . .}, and functions A(s), B(s, t), and r(s) such that Proof.For any fixed s with ∥s∥ < ϵ, bounding the modulus of the Cauchy integral by the length of the curve of integration times the modulus of the integrand implies for some C > 0. Thus, taking the coefficient of t n in (3) gives, after making use of the general binomial expansion if ϵ is small enough such that |r(0) − r(s)| < δ/3 whenever ∥s∥ < ϵ.This final sum is a geometric series with ratio less than 1, so the asymptotic behaviour of binomial coefficients implies uniformly in a sufficiently small neighbourhood of z = 1, and the stated CLT follows from Theorem 1.
At s = 0, further algebraic manipulation shows Corollary 2 applies with q = 0 and A(s) = (1 − z/r)C, giving a central limit theorem with ◁ Further CLTs of this type using the theory of singularity analysis [20] were derived (in the 1-dimensional case) by Flajolet and Soria [22,23] and, in the multivariate setting, by Gao and Richmond [24].In particular, the work of Gao and Richmond allows C(z, t) to have algebraic and logarithmic-type singularities near z = 1.Drmota [13,14] used an analytic approach (based on techniques such as the saddle-point method) to derive central limit theorems from generating functions defined implicitly by (systems of) functional equations; see also the treatment in Drmota [15,Section 2.2].A detailed survey of analytic methods for the derivation of (largely 1-dimensional) CLTs is given in Flajolet and Sedgewick [21,Chapter IX].
Central limit theorems of the type described in Corollary 2 are derived from knowledge of the generating function C(z, t) in a neighbourhood of z = 1.When further information about the behaviour of C is known for all points with the same coordinate-wise modulus, LCLTs can also be produced.
Theorem 4 (Bender and Richmond [5,Corollary 2]).Let R be a compact subset of (−∞, ∞) d and suppose that there exists ϵ > 0, a number q ∈ Q \ {−1, −2 . . .}, and functions A(s), B(s, t), and r(s) such that where (i) A(s) is continuous and non-zero in an ϵ-neighbourhood of R, (ii) r(s) is non-zero and has continuous third-order partial derivatives in an ϵ-neighbourhood of R, (iii) B(s, t) is analytic and bounded for s in an ϵ-neighbourhood of R and |t| < |r(s)|(1 + ϵ), (iv) the Hessian matrix H(s) of λ(s) = 1/r(s) is nonsingular for all s in an ϵ-neighbourhood of R, and uniformly for all k such that k/n = −∇ log r(w) has a solution for w ∈ R. Furthermore, there is a local limit theorem uniformly, where u = (j − k)/ √ n.
The key to establishing the LCLT in Theorem 4 is Condition (v), which implies that the modulus of the singularity t = r(s) of C(z, t) is uniquely minimized among the points z = e s with fixed coordinate-wise modulus when z has positive real coordinates.For instance, in one variable, information about C(z, t) near the point z = 1 is not sufficient to establish local limit theorems, one must also verify that C(z, t) has no other singularities with |z| = 1 and the same value of |t|.The advanced methods of ACSV, however, show that the multivariate situation is more complicated: in at least two variables there can be singularities that do not form 'obstructions' to deforming domains of integration, making these singularities irrelevant to determining asymptotic behaviour.We return to this in the context of ACSV in Section 2 below.
Proof Sketch.The Cauchy integral formula implies uniformly for s in an ϵ-neighbourhood of R, and following the proof of Theorem 1 we shift the mean of X k by −n∇ log r(w) and scale by 1/ √ n to get a new random variable with characteristic function for certain values of w.Tracing through the definitions of the characteristic functions, and making a polar change of variables, to establish the LCLT from the Cauchy integral above it becomes sufficient to prove that for w given in the statement of the theorem.Near the origin (up to roughly √ n-scale) it can be shown that this difference is small using (5), and the second integral, which has an explicit integrand, is o(1) when bounded sufficiently away from the origin.Roughly speaking, Condition (v) implies that |ϕ n (z)| is small when z is away from the positive real axis, so the first integral is also o(1) when bounded sufficiently away from the origin and the claimed limit theorem holds.
An extension of Theorem 4 to functions C(z, t) with algebraic and logarithmic singularities was given by Gao and Richmond [24].

ACSV and Multivariate CLTs
The multivariate results in Section 1 are derived from, in the words of Flajolet and Sedgewick [21, Page 768], "a perturbative theory of one-variable complex function theory."In contrast, the methods of ACSV in its modern form are based around the theory of complex analysis in several variables.Although the ACSV framework draws from a much larger (and more recently developed) collection of mathematical techniques 6 , it is possible to package the end results in ways that allow them to be used without understanding most of this background.
Perhaps the main advantages of ACSV are a unified framework to study the asymptotics of multivariate generating functions beyond those that exhibit the quasi-power behaviour described above (although we mainly stick to the 'smooth' case with this behaviour for our treatment here) and, as alluded to previously, a deeper understanding of which singularities contribute to coefficient asymptotics (which can greatly simplify computations).The results are also explicit to the point of being completely automated for large classes of combinatorial generating functions (such as in the software package developed for this article), allow for the computation of asymptotic expansions to arbitrary order, and work under many different sets of assumptions.
Remark 5.There are several natural ways to take the sequence indices of f i,k to infinity and search for limit theorems: for instance, one could examine coefficient slices with i For us it makes the most combinatorial sense to take the final index to infinity (as in the results above) and study coefficient behaviour in the first d indices, but the methods of ACSV adapt naturally to other situations.Work extending the types of limit theorems derived by our software package is ongoing.
In this section we give an overview of the methods of ACSV, and survey the explicit limit theorems they provide.Applications of these results are given in Section 3.

Background for ACSV
Let F (z, t) = G(z, t)/H(z, t) be a ratio of complex-valued functions G and H analytic in a domain D ⊂ C d+1 containing the origin, and suppose that F has a power series expansion valid in some neighbourhood of the origin (meaning, in particular, that H(0, 0) ̸ = 0).To derive limit theorems we study the asymptotic behaviour of coefficients Asymptotic arguments typically start with the Cauchy integral representation where T is any product of circles |z j | = |t| = ε sufficiently close to the origin.The methods of ACSV manipulate the domain of integration T to convert the Cauchy integral (6) into something that can be asymptotically approximated.As in the more classical univariate case, this process depends heavily on the singular set of the generating function F .Because F is a ratio, its singularities form a subset of the analytic variety V = {(z, t) ∈ C d+1 : H(z, t) = 0} defined by the vanishing of the denominator H, and includes all points where H vanishes and the numerator G does not.In many applications F is a rational function, in which case we may assume that G and H are coprime polynomials and the singular set of F equals the algebraic variety V defined by the vanishing of H (a similar characterization holds for general meromorphic functions, but one must introduce the notion of coprime germs of holomorphic functions).
Univariate meromorphic functions always have a finite set of dominant singularities (the singularities with minimal modulus) dictating the asymptotic behaviour, and explicit expressions for their coefficient asymptotics can be determined by adding up contributions given by each of these points.In contrast, if F is not a polynomial and the dimension d ≥ 2 then the set V is infinite and the geometry of V plays a large role in determining coefficient behaviour.In order to characterize the singularities determining asymptotics, we make the following definitions.
there is no element of V which is coordinate-wise closer to the origin, i.e., there does not exist (y, q) with |y j | < |w j | for all 1 ≤ j ≤ d and |q| < |s| such that H(y, q) = 0; • a strictly minimal point if it is minimal and no other point of V has the same coordinate-wise modulus; • a finitely minimal point if it is minimal and only a finite number of points in V have the same coordinatewise modulus; • a smooth critical point in the direction (r, 1) ∈ R d+1 >0 if it satisfies the system of equations and one of these partial derivatives does not vanish (which, in fact, implies that all of the derivatives do not vanish).
Remark 7. If H and all of its partial derivatives simultaneously vanish at a point w then either w is a zero of H with multiplicity greater than one or V is not a manifold near w.In the first case, H can be replaced by its square-free part near w to determine critical points (when H is a polynomial this means replacing it by the product of its distinct irreducible factors).In the second case, when V has non-smooth points, critical points can be defined by stratifying V into a finite collection of smooth manifolds that 'fit together nicely' and calculating critical points on each stratum.In general, if H is a polynomial then the critical points on each stratum are defined by a finite collection of polynomial equalities and inequalities that can be computed automatically from H (see [44,Section 8.2]).To simplify our presentation here we state our main results for the smooth case with zeroes of multiplicity one.

Asymptotic Results
The earliest techniques of ACSV were derived using an explicit surgery method that, in the smooth setting, yields limit theorems similar to the results of Bender and Richmond described above.ACSV determines asymptotic behaviour by manipulating the domain of integration T in (6), splitting it into some regions where the Cauchy integral is negligible and other regions where the integral can be approximated with analytic techniques.Roughly speaking, in the smooth setting one can push out the domain of integration T in (6) to approach the set of singularities of F , take a residue in one variable to reduce to a d − 1 dimensional integral lying 'on the singular set,' and then (hopefully) determine asymptotics of this lower dimensional integral using the saddle-point method.Minimal points are those to which T can be easily deformed, while critical points are those where a saddle-point analysis can be performed locally after computing a residue.
The surgery approach to ACSV (introduced for the smooth case in [43]) applies in the presence of finitely minimal critical points.The assumption of finite minimality allows one to make explicit residue computations by fixing the moduli of the z variables and varying only the modulus of the t variable.This makes the surgery method a fairly straightforward analogue of the techniques from univariate analytic combinatorics, however finite minimality is difficult to verify computationally and is a stronger condition than necessary.The surgery method is covered in detail in Melczer [41,Chapter 5], yielding asymptotic results relying on one further quantity.
Definition 8. Let (w, s) ∈ C d+1 * be a smooth critical point in the direction (r, 1) and suppose H t (w, s) ̸ = 0.The phase Hessian H(w, s) of H at (z, t) = (w, s) is the d × d matrix H with entries where U i,j = wiwj Hz i z j (w,s) tHt(w,s) and U d+1,d+1 = tHtt(w,s) Ht(w,s) .The point (w, s) is called nondegenerate if the phase Hessian matrix H(w, s) has non-zero determinant.
Theorem 9. Suppose that the rational function F (z, t) = G(z, t)/H(z, t) admits a nondegenerate strictly minimal smooth critical point (w, s) ∈ C d+1 * in the direction r ∈ R d+1 * , such that H t (w, s) ̸ = 0. Then for any nonnegative integer M there exist computable constants C 0 , . . ., C M such that where H = H(w, s) is the phase Hessian matrix and The asymptotic expansion (9) holds uniformly in neighbourhoods R ⊂ R d+1 * of r where there is a smoothly varying nondegenerate strictly minimal critical point such that H t does not vanish.
The fact that the asymptotic behaviour described in Theorem 9 varies smoothly with the direction r under small perturbations allows (with a small amount of extra analysis) for the derivation of LCLTs.
Proposition 10 (Melczer [41, Proposition 5.10]).Suppose F (z, t) = G(z, t)/H(z, t) has a power series expansion F (z, t) = (i,k)∈N d+1 f i,k z i t k at the origin such that f i,k is non-negative for all but a finite number of terms.Suppose further that, in some direction (m, 1), there is a strictly minimal critical point of the form (1, ρ) for some ρ > 0. If H t (w, ρ) and G(w, ρ) are non-zero, and the phase Hessian The requirement of strict minimality in Proposition 10 is analogous to Condition (v) in Theorem 4 above.The matrix in Definition 8 equals the H in Theorems 1 and 4 up to sign, with the entries now determined explicitly from evaluations of partial derivatives of H using (8).
Remark 11.Proposition 5.10 in Melczer [41] requires G and H to be polynomials to simplify its presentation, however the result continues to hold for analytic functions [41,Remark 5.14].Other generalizations that can be handled with the ACSV framework include (non-power series) Laurent expansions, non-smooth geometries, and non-minimal points when additional assumptions are verified.
The most computationally expensive hypothesis to verify in Proposition 10 is strict minimality (see Melczer and Salvy [42] for a complexity analysis of smooth ACSV methods).Generically7 V is smooth and the set of smooth critical point equations (7) admits a finite set of solutions, however verifying strict minimality of a point w requires examining whether V intersects the (always infinite) set of points with the same coordinate-wise moduli as w.The perturbative approach in Section 1.2 and the surgery ACSV method both require strict minimality (or finite minimality with some extra bounding), however more advanced multivariate methods show that this can be weakened.
Indeed, in the univariate setting it is the singularities of minimal modulus that contribute to dominant asymptotic behaviour, so it is tempting to assume that minimal singularities are the ones determining dominant multivariate asymptotics.In fact, as the theory of ACSV matured its methods were re-examined through more advanced mathematical frameworks, illustrating how critical points are the singularities dictating asymptotic behaviour.The most explicit results still hold for minimal critical points, however one only needs to verify that the (generically finite) set of critical points has no other elements with the same coordinate-wise modulus as a candidate minimal critical point.Without this strengthening, algorithms to compute asymptotics (or rigorously establish limit theorems) would not terminate in reasonable time even for examples with relatively low dimension and degree.
Theorem 12 (Pemantle, Wilson, and Melczer [44, Theorem 9.12]).Suppose that the rational function F (z, t) = G(z, t)/H(z, t) admits a nondegenerate minimal smooth critical point (w, s) ∈ C d+1 * in the direction r ∈ R d+1 * , such that H t (w, s) ̸ = 0 and no other critical point has the same coordinate-wise modulus as (w, s).Then the conclusions of Theorem 9 hold, where the expansion (9) now holds uniformly over neighbourhoods where there is a smoothly varying nondegenerate minimal critical point such that H t does not vanish and no other critical point has the same coordinate-wise modulus.
Theorem 12 was first discussed in Pemantle and Baryshnikov [3] using cones of hyperbolicity.To briefly summarize, if (z, t) is a minimal point then the tangent plane to the modified function H(e x1 , . . ., e x d , e q ) at (log(z), log(t)) is defined by a normal vector that is a multiple of a real vector v.If (z, t) is not critical then v is not parallel to the direction vector r and this can be used to locally deform a domain of integration near (z, t) into a region of complex space where it can be bounded and shown to be negligible.The framework of hyperbolic cones shows that these local deformations can be done in a consistent manner (away from critical points) and also generalizes to non-smooth cases.
Example 13.The rational function admits (x, t) = (1, 1) as a minimal critical point in the direction r = (1, 1), however this point is not finitely minimal as (−1, t) is a singularity for any t ∈ C, and none of Theorem 4, Theorem 9, or Proposition 10 directly apply.Since there are no other critical points with the same coordinate-wise modulus as (1, 1), Theorem 12 does apply and we can compute In fact, the local central limit theorem When dealing with combinatorial generating functions, the importance of critical points can be seen directly.Indeed, the idea underpinning cones of hyperbolicity (that the tangent space near any smooth minimal point has a normal vector that is a multiple of a real vector) implies that every minimal point is critical in some direction.When the series under consideration has nonnegative coefficients, this implies that any minimal point with the same coordinate-wise modulus as a critical point with positive coordinates is also critical [41, Corollary 5.5], giving the following result on which we base our algorithm for automatically finding and verifying LCLTs.Proposition 14. Suppose F (z, t) = G(z, t)/H(z, t) has a power series expansion F (z, t) = (i,k)∈N d+1 f i,k z i t k at the origin such that f i,k is non-negative for all but a finite number of terms.Suppose further that, in some direction (m, 1), there is a minimal critical point of the form (1, ρ) for some ρ > 0 and no other critical point has the same coordinate-wise modulus.If H t (w, ρ) and G(w, ρ) are non-zero, and the phase Hessian H of H at (1, ρ) is nonsingular, then the LCLT (10) holds.
Going even further, Baryshnikov, Pemantle, and Melczer [2] show how asymptotic behaviour can be characterized in the absence of minimal critical points using techniques from stratified Morse theory.This requires introducing the notion of critical points at infinity and gives asymptotics as a linear combination of asymptotic expansions with (generally unknown) integer coefficients, so here we stick to the explicit case of (not necessarily strictly) minimal critical points covered by Theorem 12.
Finally, although it is beyond the scope of the present article, we give a non-smooth example from forthcoming work that illustrates the differences from the smooth cases discussed above.
Example 15.Consider the generating function .
The methods of ACSV for transverse multiple points (see [41,Chapter 9] or [44,Chapter 10]) imply that the asymptotic behaviour of f λn,n varies with λ ∈ (0, 1) in a computable manner.Indeed, if 0 < λ < 1/3 or 1/2 < λ < 1 then the dominant asymptotic behaviour of f λn,n is dictated by a smooth minimal critical point vanishing on only one of the denominator factors of F , and f λn,n → 0 exponentially quickly.In contrast, when λ ∈ (1/3, 1/2) then the dominant asymptotic behaviour is determined by the non-smooth minimal critical point (1, 1) where both denominator factors vanish, and the methods of ACSV imply f λn,n → 6.Thus, unlike the smooth case where a central limit theorem is obtained, here we have a range of indices where the same minimal critical point determines asymptotics, giving a limiting distribution with a flat plateau (see the middle plot in Figure 1).A minor modification of the arguments in [?, Section 6.2]developed there only for generating functions whose denominator factors are linear -shows how to capture the transition between the different limiting regimes, which occurs on a square-root scale.Indeed, if Characterize the smallest positive solution t = ρ of H(1, t) = 0 and define m = ρHt(1,ρ) , . . ., 2. Show that H(s, . . . ,s, sρ) ̸ = 0 for all 0 < s < 1 3. Show that no other critical point has the same coordinate-wise modulus as (1, ρ)

4.
Prove that the determinant of the phase Hessian H is nonzero at (1, ρ) in the direction r = (m, 1)

5.
Prove that H t (1, ρ) and G(1, ρ) are nonzero and read off the local central limit theorem from (10) Figure 2: A schema to prove LCLTs using Proposition 14.
is the Gaussian error function then when t grows sufficiently slower than √ n (for instance, when t bounded).See Figure 1 for an illustration.◁

Verifying Minimality
Because verifying criticality is easier than verifying minimality, most ACSV algorithms for asymptotics fix a direction, compute critical points in this direction, and then study the critical points to determine which are minimal points.In contrast, because Proposition 10 requires a minimal critical point of the form (1, t), to prove an LCLT it is often easiest to use (7) to discover a direction r = (m, 1) with critical points of this form and then verify the required conditions.Determining minimality is easier for points with positive coefficients when F (z) has only a finite number of non-negative coefficients, as it does under our assumptions.The following result should be seen as a multivariate generalization of the well-known Vivanti-Prinsheim theorem in the univariate case, and the approach to proving an LCLT that it suggests when combined with Proposition 14 is summarized in Figure 2.
Although proving strict minimality is difficult in general, there is one case arising in some combinatorial applications where it is automatic.Definition 17.A power series S(z) = n∈N d p n z n is called aperiodic if every element of Z d can be written as an integer linear combination of the exponents {n ∈ N d : p n ̸ = 0} appearing in S.
Proposition 18 (Melczer [41,Proposition 5.5]).Suppose F (z) = G(z)/H(z) is a ratio of functions G and H.If H(z) = 1 − S(z) for some aperiodic power series S with non-negative coefficients then every minimal point that is within the domain of convergence of the power series expansion of F (z) is strictly minimal and has positive real coordinates.
Proof.Suppose that w is a minimal point and for each 1 ≤ j ≤ d write w j = x j e iθj with x j > 0 and θ j ∈ R. Let s n denote the coefficient of z n in S(z).Then

A Family of Limit Theorems
In this section, and the next one, we illustrate how to apply the approach of Figure 2 to prove LCLTs.We begin with the example that originally motivated our work.Definition 19.For any d ∈ N, let F d (n) be the set of permutations σ on {1, . . ., n} such that i − d ≤ σ(i) ≤ i + 1 for all i.Let F d be the union of sets F d (n) for all n ∈ N.
Note that every element of F d , when written in disjoint cycle notation, has cycles of length at most d + 1.

Proposition 20 (Chung et al. [9, Theorem 1]). The number of permutations in
t n ]F (z, t) in the power series expansion of the rational function Chung et al. [9] prove Proposition 2 by considering the set of perfect matchings of the graph associated with F d and using this to find an explicit recurrence that can be manipulated.One of the motivations for their study of this family is a relationship to the determination of sample sizes required for sequential importance sampling of certain random perfect matchings in classes of bipartite graphs.Experimentally checking the properties of Proposition 14 in low dimension using a computer algebra system, we were initially surprised to find that some of the properties did not hold!To better understand the behaviour of the coefficients we thus plotted the coefficients of [t 150 ]F (z, t) for d = 1, shown on the left of Figure 3.
Figure 3 shows the problem establishing a limit theorem on the coefficients of F (z, t): the coefficients approach a normal distribution, but the distribution is supported on a d-dimensional slice of R d+1 .Indeed, if the size n of one of our restricted permutations is fixed, and the number i k of k cycles it contains is specified for all k ≥ 2, then its number of one cycles (i.e., fixed points) is uniquely determined as We thus prove Conjecture 21 true by setting z 1 = 1 and applying the techniques of ACSV to the coefficients of F (1, z 2 , . . ., z d , t). where • and H is the non-singular matrix with entries Theorem 22 follows directly from Proposition 14 after some direct calculations, the most difficult of which is computing a closed form for the determinant of the Hessian matrix H.We complete these calculations in Section 4, using a guess-and-check method for symbolic determinants which is a useful tool for establishing limit theorems such as these with parameterized dimension.In fact, we will prove the following more general result.
Theorem 23 (Main Theorem).Let F (z, t) = G(z,t) H(z,t) be a ratio of functions where , and let ρ be the smallest positive root of P (t).Suppose that • each of the q i (t) is a non-zero polynomial vanishing at the origin and q(t) is a complex-valued analytic function for |t| ≤ ρ that vanishes at the origin, • the power series expansion of S(z, t) = 1 − H(z, t) = q(t) + k q k (t)z k at the origin has non-negative coefficients, • the exponents appearing in the power series q(t) have greatest common divisor 1, and As n → ∞, the maximum coefficient of [t n ]F (z, t) as a polynomial in z 1 , . . ., z d approaches where H is the non-singular d × d matrix whose determinant is non-zero.Furthermore, Remark 24.The Sage package accompanying this article automatically proves central limit theorems for all explicit rational functions satisfying the conditions of Proposition 14, which includes all those satisfying the conditions of Theorem 23.
Although the form of H in Theorem 23 might seem restrictive, generating functions of this form appear quite frequently when tracking parameters in combinatorial classes using, for instance, the symbolic method framework described in Flajolet and Sedgewick [21].We end this section by describing some other combinatorial limit theorems it captures.

Strings with Tracked Letters
Let A = {a 1 , . . ., a ℓ } be an alphabet with ℓ letters and let Ω = {ω 1 , . . ., ω d } ⊂ A be a subset of d letters we wish to track.The multivariate generating function enumerating such strings is which trivially satisfies the hypotheses of Theorem 23 when ℓ > d (if ℓ = d then all coefficients live in a ddimensional hyperplane of R d+1 , because adding the number of occurrences of each letter gives the length of the string).Thus, as n → ∞, the maximum coefficient of [t n ]F (z, t) as a polynomial in z 1 , . . ., z d approaches where m = 1 ℓ , . . ., 1 ℓ and H is the non-singular d × d matrix with off-diagonal entries − 1 ℓ 2 and diagonal entries ℓ−1 ℓ 2 .Note that tracking the number of 1s in binary strings (where ℓ = 2 and d = 1) recovers the classical central limit theorem,

Compositions with Tracked Summands
Fix a positive integer d and recall that an (integer) composition of size n is an ordered tuple of positive integers which sum to n.If C is the class of compositions, enumerated by size and the number of times each element of {1, . . ., d} occurs, then C has the multivariate generating function where where m = 1 4 , 1 8 , . . ., 1  2 d+1 and H is the d × d matrix with off-diagonal entries H i,j = −2 −i−j−2 (i + j − 3) and diagonal entries H j,j = 2 −2(j+1) (2 j+1 − 2j + 3).
The flexibility of Theorem 23 allows us to modify the restrictions on the elements appearing or tracked among the compositions under consideration.For instance, for any finite set Ω ⊂ Z >0 the multivariate generating function enumerating compositions by size and number of times each element of Ω = {ω 1 , . . ., ω d } occurs is where z k tracks the number of occurrences of ω k .The hypotheses of Theorem 23 still hold, so a local central limit theorem applies (although, of course, the parameters of the limiting distribution will depend on which summands are tracked).
Even more generally, we can restrict our compositions to summands in some set Λ ⊂ Z >0 and track occurrences of elements in a finite subset Ω = {ω 1 , . . ., ω d } ⊂ Λ.The relevant multivariate generating function becomes When the elements of Ω are coprime, Theorem 23 applies, so a local central limit theorem holds (when the elements of Ω have greatest common divisor larger than 1 then some intervention is needed to determine which singularities determine asymptotic behaviour).

n-Colour Compositions with Tracked Summands
Finally, we consider the class of n-colour compositions, introduced by Agarwal [1] and studied in Gibson et al. [26], where each integer i is coloured by one of i available colours (each summand must be coloured and different colourings give distinct n-colour compositions).
Example 27.The twenty-one n-colour compositions of 4 are where the subscripts denote the colours assigned to each summand.
Once again, we fix a positive integer d and track the number of times each element of {1, . . ., d} occurs.If S n is the class of positive integers where each integer i is coloured with one of i possible colours then is the multivariate generating function of positive integers where z k tracks the number of occurrences of k, so the corresponding multivariate generating function for n-colour compositions is As expected, the hypotheses of Theorem 23 hold, and we get a (messier) LCLT whose parameters are given explicitly in the Sage notebooks corresponding to this paper.Similar to our last example, we may further restrict which elements (or colours!) are allowed or tracked, and immediately get LCLTs.The proofs for such extensions follow similarly to those for integer compositions.

Proof of Theorem 23
We now prove Theorem 23 by applying the outline in Figure 2 to verify the conditions of Proposition 14.To that end, assume the hypotheses of Theorem 23 hold.
Step 1: Finding the correct direction.
Following the outline of Figure 2, we substitute w = (1, ρ) into the smooth critical point equations ( 7) to find that w is a critical point in the direction (m, 1) where Steps 2 and 3: Establishing minimality.
Our assumptions imply that the polynomial d k=1 q k (t) and series expansion of q(t) have non-negative coefficients, and 0 < q(ρ) < 1, which implies that |s| increases as 0 < t < ρ decreases.Thus H(s, . . ., s, sρ) ̸ = 0 for all 0 < s < 1, and w is minimal by applying Lemma 16 to the function 1/H(z, t) = 1/(1 − S(z, t)) whose series expansion at the origin has non-negative coefficients.Our assumptions imply that S(z, t) is an aperiodic power series with non-negative coefficients, so Proposition 18 implies that no other singularities (including no other critical points) have the same coordinate-wise modulus as w.
We now prove that the phase Hessian matrix defined by (8), which simplifies to (12) in our case, has nonzero determinant.To compute this symbolic determinant, we originally used the Sage computer algebra system to compute and factor the Hessian determinant for the permutation generating function described by Proposition 20 in small dimensions.Observing a pattern in the factors, we were able to conjecture, and then a posteriori prove, an LU -factorization for the Hessian matrix in this case, and then extend to the general case.This approach immediately gives the Hessian determinant: if HU = L for a lower triangular matrix L and an upper triangular matrix U with diagonal elements equal to 1 then det H is simply the product of the diagonal elements of L.
Remark 28.When the dimension d is fixed, the matrix equation HU = L defines an explicit system of equations in terms of the entries of U and L, so it is often possible to computationally determine suitable U and L in low dimension, conjecture their general form, then prove this inference.This approach to symbolic determinants is described in the well-known treatise of Krattenthaler [33], which attributes the popularization of such a 'guess-and-check' LU -factorization to George Andrews after he used it to great effect in a variety of papers starting in the 1970s.A companion Sage notebook to this paper gives a procedure to calculate U and L for any fixed d.Studying the numerator and denominator of the rational function entries we are able to deduce certain patterns, such as the denominators being constant down columns, which leads us to conjecture that in general dimension HU = L where Although these formulas are quite involved, we note that by keeping P ′ (ρ) and P ′′ (ρ) as symbolic parameters the entries of U and L are independent of d, which allows us to algorithmically verify that HU = L with the aid of a computer algebra system.We break this verification into three cases depending on the behaviour of the entries of U and L: see Figure 5 for an illustration of the different cases on a 4 × 4 matrix.
= Case 1: j > i Case 1 (j > i): As U is upper-triangular, we have where we split the sum in such a way that the summands in the indefinite series have a uniform definition.Using Sage for algebraic manipulations, we see that β ij = P ′ (ρ) q ′ j (ρ)ρ − q j (ρ) − P ′′ (ρ)q j (ρ)ρ + q ′ j (ρ)ρB j − q j (ρ)ρD j (15) Since α ij , β ij , γ ij , and the denominator in our above expression for H ia U aj are not dependent on a, Substitution of our above formulas then gives, after some algebraic simplification by Sage, that 1≤a<j a̸ =i in this case, as required.H ja U aj + H jj = L jj using the definitions above.Analogously to Case 1, + β jj q a (ρ) + γ jj q ′ a (ρ) ρ 2 P ′ (ρ) 3 r j where α jj , β jj and γ jj are defined in ( 14), ( 15) and ( 16) respectively.Since α jj , β jj , γ jj , and the denominator in this expression for H ja U aj are not dependent on a, we can algebraically simplify to get j−1 a=1 + β jj q a (ρ) + γ jj q ′ a (ρ) ρ 2 P ′ (ρ) 3 r j   = α jj D j + β jj A j + γ jj B j ρ 2 P ′ (ρ) 3 r j .
Substitution of our above formulas then gives, after some algebraic simplification by Sage, that j−1 a=1 H ja U aj + H jj − L jj = 0 in this case, as desired.
Case 3 (j < i): At this point, we have proven that HU is lower-triangular and described its diagonal entries, which gives the desired determinant.Although not needed to establish Theorem 23, for completeness we prove the claimed expression above for the below diagonal entries in our companion Sage notebook.
Step 3b: Computing the Hessian determinant Our LU -factorization expresses the Hessian determinant as the product of the diagonal entries of the lowertriangular matrix L. In fact, the form of these diagonal entries causes cancellation of many terms, leading to a relatively compact final answer, q k (ρ) = 0 and P ′ (ρ) = −q ′ (ρ) − d k=1 q ′ k (ρ) so that A d+1 = d k=1 q k (ρ) = 1 − q(ρ) and B d+1 = d k=1 q ′ k (ρ) = −q ′ (ρ) − P ′ (ρ).Making these substitutions into our definition of r d+1 gives, after some algebraic simplification, that the Hessian determinant has the form (13).
Remark 29.When each of the q k are monomials this inequality becomes an equality. Step

4: Checking final hypotheses
To conclude Theorem 23 from Proposition 10 it remains only to note that G(1, ρ) ̸ = 0, which is one of our assumptions, and that H t (w) = P ′ (ρ) ̸ = 0 because none of the q k (t) and q(t) have constant terms and we assume the denominator is not constant.

Conclusion
The methods of analytic combinatorics in several variables, while perhaps daunting to some outside users due to its reliance on a wide breadth of mathematical techniques, provide some of the most powerful tools for the study of multivariate generating functions.Although Bender and Richmond [5] already provided techniques to prove local central limit theorems for a variety of combinatorial generating functions, verifying required conditions on analytic regions for the generating functions is too expensive to implement in a general, practical, algorithm.Using results of ACSV, a better understanding of the singular sets of multivariate generating functions yields such an algorithm, in addition to providing a framework for further generalizations.It is our hope that by putting the results of ACSV into context with past probabilistic work, giving a simple outline of how to apply the results, implementing the results in a computer algebra package, and using the results to prove a family of limit theorems, that readers are inspired to look further into this growing area of combinatorics.
combinatorial contexts, one often has a combinatorial class (C, | • |) defined by a set of objects C and a size function |•| : C → N such that there are only a finite number of objects of any given size.A d-dimensional parameter on such a class is any map χ : C → N d , and the multivariate generating function with respect to the triple (C, | • |, χ) is the (d + 1)-dimensional power series 1 C(z, t) = n≥0   i∈N d c i,n z i   t n

where•
A(s) is continuous and non-zero, • r(s) is positive and has continuous third-order partial derivatives, and • B(s, t) is analytic and bounded for ∥s∥ < ϵ and |t| < |r(0)| + δ.If the Hessian matrix of λ(s) = 1/r(s) is nonsingular at s = 0 then the CLT in Theorem 1 holds.

Figure 1 :
Figure 1: The coefficients c(k) = f k,N (in blue) compared to their limiting behaviour (in red and orange) when N = 1000.The transitions in behaviour near N/3 and N/2 are shown on intervals of length 2 √ N in the left and right plots.
since w is within the domain of convergence of G and H and hence |S(w)| ≤ S(|w 1 |, . . ., |w d |).

Conjecture 21 (
Chung et al. [9, Page 45]).For fixed d the joint limiting distribution of the number of k-cycles approaches a multivariate normal distribution as n → ∞.

Figure 3 :
Figure 3: Left: The coefficients of [t 150 ]F (z, t) when d = 1, divided by the maximum coefficient and shifted to put the maximum at the origin.Right: The coefficients of [t 150 ]F (1, z 2 , t) compared to their limiting normal distribution and shifted to put the maximum at the origin.

Remark 26 .
Restricting compositions to use only the numbers 1, . . ., d + 1 and tracking occurrences of those numbers gives the same generating function as the family of restricted permutations discussed above, whose behaviour is described in Theorem 22.

Figure 4 :
Figure 4: Left: The coefficients of [t 100 ] 1 − z 1 t − t 2 /(1 − t) −1 divided by the projected maximum coefficient compared to their limiting normal distribution and shifted to put the maximum at the origin.Right:The coefficients of [t 100 ] 1 − z 1 t − z 2 t 2 − t 3 /(1 − t)−1 divided by the projected maximum coefficient and shifted to put the maximum at the origin.

Figure 5 :
Figure 5: The three cases, distinguished by colour, for d = 4.

Case 2 ( 1 a=1H
i = j): Since U is upper-triangular we again have(HU ) jj = d a=1 H ja U aj = j a=1 H ja U aj = j−ja U aj + H jj , = Case 2: i = jwhere we separate the case a = j since the entries of U have a different definition on the diagonal and above the diagonal.It is therefore sufficient to prove j−1 a=1 k (ρ) r d+1 P ′ (ρ) d ρ d r 1 = (−1) d d k=1 q k (ρ) r d+1 P ′ (ρ) d+2 ρ d+1 .Furthermore,P (ρ) = 1 − q(ρ) − d k=1 t is the multivariate generating function of positive integers where z k tracks the number of occurrences of k.Once again, it is easy to verify the hypotheses of Theorem 23 hold with ρ = 1/2, the smallest positive root ofP (t) = H(1, t) = 1 − t 1−t .Figure 4 illustrates the corresponding local central limit theorem when d = 1 and d = 2. Theorem 25.As n → ∞ the maximum coefficient of [t n ]C(z, t) as a polynomial in z 1 , . . ., z d approaches A n = 2 n n −d/2 2