Full asymptotic expansion for orbit-summable quadrant walks and discrete polyharmonic functions

Enumeration of walks with small steps in the quadrant has been a topic of great interest in combinatorics over the last few years. In this article, it is shown how to compute exact asymptotics of the number of such walks with fixed start- and endpoints for orbit-summable models with finite group, up to arbitrary precision. The resulting representation greatly resembles one conjectured by Chapon, Fusy and Raschel for walks starting from the origin (AofA 2020), differing only in terms appearing due to the periodicity of the model. We will see that the dependency on start- and endpoint is given by discrete polyharmonic functions, which are solutions of $\triangle^n v=0$ for a discretisation $\triangle$ of a Laplace-Beltrami operator. They can be decomposed into a sum of products of lower order polyharmonic functions of either the start- or the endpoint only, which leads to a partial extension of a recent theorem by Denisov and Wachtel (Ann. Prob. 43.3).


Introduction and Motivation
Enumeration of lattice paths has by now become a standard problem in combinatorics.In particular the asymptotics of lattice path counting problem in the half-and quarter plane have been addressed in a variety of works.To do so, the by now standard approach as used in [2] is to consider the generating function Q(x, y; t) = n∈N k,l∈Z x k y l t n q(0, B; n), where q(0, B; n) denotes the (possibly weighted) number of lattice paths starting at (0, 0) and ending at a point B = (k, l).While for walks in the entire plane this contains negative powers of x and y and is thus not a generating function in these two variables, it is one in t.Utilizing a functional equation for Q(x, y; t), it is then often possible to extract information about the asymptotic behaviour of q(0, B; n), as was done for instance in [2], where the authors gain various asymptotic expressions for the case of a directed model (i.e. the first component of a step is always positive) closely tied to the zeros of the so-called kernel.In the undirected case, restricting ourselves to walks in the quarter plane, asymptotics of the coefficients of Q(1, 1; t) or Q(0, 0; t) (that is, of the walks starting at the origin and ending either anywhere or back at the origin) have been computed for some cases in [26], using a complex boundary value problem, and for a family of models related to the Gouyou-Beauchamps model in [20] via Analytic Combinatorics in Several Variables (ACSV, [38]).In [8], the authors give hypergeometric expressions for the generating function in the 19 cases where the group is finite (see [15]) and the model not algebraic, and show that it is transcendental.Additionally, they find explicit expressions for some specializations of the generating function for the 4 algebraic models.To do so, they use orbit summation paired with creative telescoping.This goes in a similar direction as [10], and in general the question about the algebraic properties of the generating function Q(x, y; t).This function is known to be D-finite in x, y, t if and only if the group is finite, which was first conjectured in [15] and subsequently proven by various authors, for instance [9,41,35,32].
There have also been extensions of the usual setting, where either steps are not required to be small or not to be homogeneous [27,34], or where the setting is not in the quarter plane but some other domain, e.g. in 3 or more dimensions [6,16,40,29,3].
Instead of counting all walks of a certain length, however, one could also sort them by their endpoint.For some standard models, like for example the Simple Walk (S = {←, ↑, →, ↓}), or the Tandem Walk (S = {←, ↑, ց}), explicit formulas for q(0, B; n) are given in [15], making use of coefficient extraction and orbit summation.One of the main results in the aforementioned article is that one can make use of orbit-summation in order to obtain an expression for the generating function Q(x, y; t) in order to show it is D-finite, but it turns out that this approach can be utilized to compute asymptotics of q(0, B; n) to arbitrary high order instead.For highly symmetric models, this is used in [39] in order to find diagonal representation of the generating functions, and then use ACSV to compute the asymptotics of the coefficients of Q(1, 1; t).In [40], this is then refined to all orbit-summable models with finite group.Since the principal idea behind the saddle point method and ACSV is in fact rather similar (even though the latter is more general and uses much heavier machinery), it seems likely that their approach could be generalized to compute exact asymptotics of q(k, l; n) as well.A general first order approximation with fixed start-and endpoint is given in [22], where using coupling with Brownian motion the authors show that, under some moment conditions, we have where q(A, B; n) is the number of paths with starting-and endpoints A = (u, v) and B = (k, l) respectively, and V (A), Ṽ (B) are discrete harmonic functions in u, v and k, l respectively.Note that these moment conditions are always satisfied for bounded step sets, which includes all cases considered in this article.
In [18], the authors used methods from [4] to show that for a Brownian motion, the asymptotic expansion of the heat kernel with respect to time contains (continuous) polyharmonic functions.They noticed as well that a discrete analogue of this exists for the Simple Walk, and conjectured that this might also be the case for other models.Using orbit summation together with a saddle point method, we will show this to be true in a more general context in Thm.3.1 and Thm.4.1.In particular, for all orbit-summable models (see Section 2), we obtain for any given m ∈ N asymptotics of the form where c ∈ N, γ ∈ R and the α i , β i , ζ i are roots of unity.The coefficients v p (x, y) are so-called discrete polyharmonic functions (for a formal definition see Section 2) of order p in B and, in reversed direction, in A. If the model has no drift then they are polynomials, otherwise they contain additional exponential factors.Compared to (1.3), the expansion appearing in [18] is different in the lack of the α i , β i , ζ i , which appear due to periodicity properties: at certain points we might have q(A, B; n) = 0, and at those points some kind of cancellation needs to occur.The structure of the polyharmonic functions v p (x, y) will be described more closely in Thm.4.2, where it will be shown that we can write where the index p is the same as in (1.3) and the h i (A) and g i (B) are (in the case of y adjoint) polyharmonic of degree at most p in A and B respectively (this will be detailed in Thm.4.2), and the number of summands k depends on p.For p = 1 this reduces precisely to (1.2), so in a certain sense one could view Thm.4.2 as an extension (albeit under much stronger conditions) of ( [22,Thm. 6]).
While the results of this article are stated for orbit-summable small step models in two dimensions, they do generalize to higher dimensions andwith a slight condition -to orbit-summable large step models as well, see the remarks after Thm.3.1 as well as Appendix A and Appendix B. This article is structured as follows: 1.In Section 2, a short overview over some definitions and the tools utilized will be given.

2.
In Section 3, it will be shown that models with finite group allowing for orbit summation in a manner similar as in [15] satisfy an asymptotic relation of the form (1.3), and the example of the Gouyou-Beauchamps model will be worked out in detail.In particular it will be explicitly shown how to compute the constants γ, c, α i , β i , ζ i and the functions v p (B) in (1.3).
3. In Section 4, we will consider the same problem where instead of (0, 0), we start our paths at an arbitrary point (u, v).It turns out, maybe not surprisingly, that due to the symmetry of the problem the resulting solution is similar to that obtained in Section 3, leading to Thm. 4.1.
Continuing from there, we can then find a decomposition of the v p as in (1.4) in Thm.4.2.The resulting decomposition of the asymptotic terms for the Simple Walk is given in Section 4.2, using an explicit basis of polyharmonic functions constructed in [42].

Walks in the Quarter Plane
In order to count lattice paths in the quarter plane, first of all we need a set S ⊆ Z 2 of permissible steps, together with a family of weights (ω s ) s∈S .A lattice path of length n is then a sequence of points (x 0 , . . ., x n ) ⊆ Q n , with Q := Z ≥0 × Z ≥0 being the quarter plane, such that x k − x k−1 ∈ S for all k ∈ {1, . . ., n}.We will count such a path by its weight, that is, the product n k=1 ω x k −x k−1 .We will in the following assume that: 1. our step set consists of small steps only, i.e. S ⊆ {−1, 0, 1} 2 \ {(0, 0)}, 2. our step set is non-degenerate, i.e. there is no (possibly rotated) halfplane containing all allowed steps.
In order to keep the notation short, we will in the following sometimes denote by S not only the set of allowed steps, but also the associated weights (ω s ) s∈S .Let q((u, v), (k, l); n) be the (weighted) number of paths of length n from (u, v) to the endpoint (k, l).It can then be shown (see e.g.[15, Lemma 4]) that the generating function satisfies the functional equation where where S(x, y) is the step-counting Laurent polynomial For ease of notation, if (u, v) = (0, 0), then we will just write q(k, l; n) and Q(x, y; t).
Given that we consider paths in the quarter plane only, we see that the series Q u,v (x, y; t) are indeed power series in x, y and t; in particular they are power series in t with polynomial coefficients in x, y.Given some power series A(x, y) = i,j∈N a i,j x i y j , we denote by [x n ] the linear operator extracting the n-th coefficient, i.e. [x n ]A(x, y) = j∈N a n,j y j , and [x n y m ]A(x, y) = a n,m .In the same manner, we define for Laurent series the operator [x > ] := ∞ n>0 x n [x n ] extracting all positive powers with their coefficients, and the operator [x ≥ ] extracting all nonnegative ones.As we allow small steps only and our step set is non-degenerate, we can write the kernel K(x, y) as with a(x), b(x), c(x), ã(y), b(y), c(y) all being non-zero.Consider now the two birational transformations Ψ : (x, y) → x −1 c(y) ã(y) , y . (2.8) These two transformations, which clearly depend on our step set S, generate the so-called group of a model, denoted by G.This group can be either finite or infinite.If G is finite, then as both Φ, Ψ are involutions, any element g ∈ G can be written as either (Φ • Ψ) k , or as Ψ • (Φ • Ψ) k for some k.We define sgn(g) = 1 in the first, and sgn(g) = −1 in the second case.The study of this group has been central to many results about lattice walks in the quarter plane, see e.g.[15,32,23,26].The main reason for this is that if we let then it follows that k(x, y) is invariant under G.This observation is the starting point for orbit summation methods as is done in [15]: we rewrite (2.2) (leaving out the dependency on t for easier readability) as Assuming now that the group is finite, and picking any g ∈ G, we hence obtain In the following, it will be convenient to consider models with zero drift, that is, where If this is not the case, then we can utilize the Cramer-transformation: we multiply each weight ω i,j by a factor of α i β j , where we choose α, β such that the drift will be 0. The existence of such an α, β is ensured for non-singular models, see e.g.[22, 1.5].The reason why this substitution is very convenient combinatorially is fairly simple; given the number q((i, j), (k, l); n) of paths from (i, j) to (k, l) with n steps weighted by old weights ω i,j , then for the equivalent q using the new weights we have q((i, j), (k, l); n) = α k−i β j−l q((i, j), (k, l); n). (2.14) Additionally, it turns out that the group of the model is, in a certain sense, invariant under this transformation, and in particular orbit-summability is preserved.
Lemma 2.2.Let S be a model and Ŝ be a Cramer-transform of S, with weights α, β.Let G, Ĝ be the respective groups with generators Φ, Ψ and Φ, Ψ. Lastly, let ι be the mapping (x, y) → (αx, βy).Then we have: Ŝ is orbit-summable if and only if S is orbit-summable.
Proof.We show the first part for Φ only, the statement for Ψ will then follow by symmetry.We can define ĉ(x) and â(x) as in (2.5).One finds that, by definition of Ŝ, we will have â(x) = βa(αx) and ĉ(x) = β −1 c(αx).From this it follows that The isomorphism of G and Ĝ follows immediately and is given by g → ι −1 • g • ι ∈ Ĝ.By this isomorphism, we also see immediately that positive and negative powers of x, y are preserved, which is in turn all that matters for orbit-summability, hence we are done.
Given a model with zero drift, one can consider the correlation coefficient θ (see [26,47]) given by θ = arccos It turns out that this correlation coefficient is closely linked to the properties of the group (its restriction to the surface C := {(x, y) : K(x, y) = 0} is finite if and only if π/θ ∈ Q [25]), and also has a geometric interpretation (see [47]).Additionally, in [46], [33], it is shown that π/θ is tied to the growth of the positive harmonic function.Given a model with non-zero drift, we can first apply a Cramer-transform to get rid of the drift, and then define θ as above.

Discrete Polyharmonic Functions
Given a step set S and a discrete function f defined on the quarter plane Q, we can define the Markov operator If we take the discrete random walk (X n ) with the transition probabilities given by the reverse of S, and the induced Markov chain M n := f (X n ), then we can interpret the operator P as the expectation One can then proceed to look at the expected change during a time step, weighted by a parameter t, which is given by The operator △ can be viewed as the discrete equivalent of a Laplace-Beltrami operator.We call a function If t = 1, then we simply speak of harmonic and polyharmonic functions respectively.Note that strictly speaking, the polyharmonicity is given by the first conditions, while the second ones are Dirichlet boundary conditions.Due to the underlying combinatorics, however, one immediately sees that we are only interested in polyharmonic functions satisfying the Dirichlet problem as above.In particular (1-)harmonic and biharmonic functions have a variety of applications in the theory of stochastic processes and physics, see e.g.[36,1,19].The occurrence of t-polyharmonic functions in the asymptotics of path counting problems, however, was noted only fairly recently in [18].In [42], it is shown that in the zero drift case, the space of 1-polyharmonic functions of order n is isomorphic to C[[z]] n , and a basis consisting of polyharmonic functions with rational generating function is given.Given a model S, we can also define a discrete Laplacian △ for the model with directions reversed S. Since △ is the adjoint operator to △ on the space L 2 (Z 2 ), we will call it the adjoint Laplacian.In Section 4.1, we will encounter functions of the form f (x, y), for x, y ∈ Z 2 , which are polyharmonic in x and adjoint polyharmonic in y.We will call such a function f (x, y) multivariate polyharmonic of order p if for all 0 ≤ k ≤ p.Note that one can verify by an easy computation using linearity that the ordering of the Laplacians does not matter as they commute, i.e. we have Lastly, discrete polyharmonic functions behave well with respect to Cramer transformations: via a short computation, one can show Lemma 2.3.Let Ŝ be a Cramer-transform of S, with ωi,j = α i β j ω i,j .Let △, △ be the associated Laplacians.We then have This directly implies that we have a bijection between the polyharmonic functions w.r.t.△ and those w.r.t.△, given by adding a factor of α −k β −l .

Saddle Points
In the following we will use the saddle point method as described for instance in [30,VIII].In particular, we will be interested in saddle points of S(x, y), which is a Laurent polynomial with only positive coefficients.We call a point (x 0 , y 0 ) a dominant saddle point if (x 0 , y 0 ) is a minimizer of S(x, y) in R + × R + .Note that this minimization property directly implies the usual defining property of a saddle point, namely that the orthogonal directional derivatives vanish.Additionally, note that S(x 0 , y 0 ) > 0, due to the positivity of coefficients.
Lemma 2.4.For any non-degenerate model, a dominant saddle point exists.
Proof.Since the model is non-degenerate, we know that S(x, y) is coercive on R + × R + (it goes to infinity wherever we approach the boundary; see also [22, 1.5]).Therefore, S(x, y) must attain its minimum at some point (x 0 , y 0 ), which is then by definition a dominant saddle point.

By a short computation, one can check that
Lemma 2.5.The dominant saddle point s 0 = (x 0 , y 0 ) is equal to (1, 1) if and only if the drift of the model is 0.
When using a Cauchy type integral to compute asymptotics later on, we will integrate over the torus {(x, y) : |x|= x 0 , |y|= y 0 }.We will call a point (x i , y i ) on this torus a saddle point associated with (x 0 , y 0 ), if |S(x i , y i )|= S(x 0 , y 0 ).Note that by the positivity of the coefficients of S(x, y), we know that S(x, y) takes the maximum (absolute) value on the aforementioned torus at (x 0 , y 0 ).Indeed, for any other point (x i , y i ) we can have S(x i , y i ) = S(x 0 , y 0 ) only if there is some From here, it is not difficult to see that α i , β i and ζ i must be roots of unity.It is also clear that there can only be finitely many such ζ i , in a one-to-one correspondence with finitely many pairs (x i , y i ), 0 ≤ i ≤ l (the maximum l appearing for the 19 unweighted orbit-summable models is 3, see [43, App.D]).One can check directly that the (x i , y i ) are indeed saddle points of S(x, y), moreover, as we will see in Lemma 3.1, the local behaviour of S(x 0 , y 0 ) and S(x i , y i ) is the same up to the factor ζ i .By the same reasoning as in [30,VIII], it turns out that when computing the asymptotics of the coefficients via the Cauchy formula, the main contributions to the contour integral come from the points (x i , y i ), as the modulus of S(x, y) will be smaller elsewhere, leading to exponentially smaller terms.Note that while this article considers mainly the 2-dimensional case, all the definitions above extend to more dimensions in a natural manner; only the structure of the group becomes more complicated as we will have more than two transformations, see e.g.[6].

Full asymptotic expansion
The goal of this section is to compute the asymptotics of orbit-summable lattice walks from the origin to an arbitrary but fixed point in the quarter plane, and in particular to show the following: Theorem 3.1.Let S be a step set satisfying the general assumptions stated in Section 2 and such that S is orbit-summable, i.e.
Suppose that s 0 = (x 0 , y 0 ) is a dominant saddle point, with associated other saddle points s i , 1 ≤ i ≤ r (meaning that we consider r + 1 saddle points in total).Furthermore, let α i , β i , ζ i be the roots of unity as constructed in Section 2.Then, there is a constant c ∈ N, a constant γ > 0, and γpolyharmonic functions v p of degree p such that for any m ∈ N we have The polyharmonic functions v p (k, l) are polynomials precisely if the drift of the model is zero, else they contain an additional factor of x −k 0 y −l 0 , with (x 0 , y 0 ) the dominant saddle point.They can be computed explicitly via a Cauchy-type integral.Lastly, the constant c can be expressed using the correlation coefficient θ defined in (2.17) via c = π/θ, and for the exponential growth γ we have γ = min x,y∈R + S(x, y).
In order to keep the proof of Thm.3.1 reasonably concise, we will first start with two somewhat technical lemmas.In the first lemma, we will establish some periodicity properties of S(x, y) and the group, in the case where we have multiple saddle points.Lemma 3.1.Let (x 0 , y 0 ) be a dominant saddle point, and (x i , y i ) = (α i x 0 , β i x 0 ) be the associated ones, with S(x i , y i ) = ζ i S(x 0 , y 0 ).Define furthermore φ(x, y), ψ(x, y) such that Φ(x, y) = (x, φ(x, y)), Ψ(x, y) = (ψ(x, y), y) are the generators of the group as in (2.7),(2.8).We then have, for all x, y ∈ C: (3.4) Remark: Lemma 3.1 still holds true in more than two dimensions, with a completely analogous proof.Also note in particular that (3.2)-(3.4)hold true for all x, y ∈ C, and not only saddle points.
Proof.As (x 0 , y 0 ) is a dominant saddle point and (x i , y i ) := (α i x 0 , β i y 0 ) associated to it, we have |α i |= |β i |= |ζ i |= 1, and know that for each monomial x k y l appearing in S(x, y) we have (α i x) k (β i y) l = ζ i x k y l .Consequently, for all such k, l we have α k i β l i = ζ i , and thus (3.2) holds.Since Ψ by construction changes (x, y) to another pair (ψ(x, y), y) such that both S(x, y) = S(ψ(x, y), y), we know (remember that S(x, y) is quadratic in x, y) that if then we have ψ(x 1 , y) = x 2 .Now suppose we have a pair (x, y) such that ψ(x, y) = x.Then, using (3.2), we can rewrite: In this case, from the above we can conclude that (3.3) holds whenever ψ(x, y) = y.But ψ(x, y) = x can be true only when ∂S ∂x (•, y) = 0, i.e.where for a given y we have no solution other than x to S(•, y) = S(x, y).This means that we have ψ(x, y) = y almost everywhere.Since the functions on the leftand right-hand side of (3.3) are rational and agree almost everywhere, we know that they are indeed the same.By a symmetric argument, we show (3.4).
Lastly, we will show that given an asymptotic representation as in Thm.3.1 below, the functions appearing therein are indeed polyharmonic.Lemma 3.2.Suppose q(B; n) is a (combinatorial) quantity satisfying for all B ∈ Z ≥0 × Z ≥0 , n ≥ 0, and that at the same time it is of the form for all k ≥ 0, with the ζ i pairwise different roots of unity.Then, the v p,i are γ-polyharmonic of degree p. If, additionally, for a fixed point B we know that q(B; n) = 0 ∀n ∈ N, then v p,i (B) = 0 for all p, i.
Proof.Substituting (3.9) into (3.8)gives us Extracting the terms for p = 1 and noticing that the others are smaller by a factor of at least 1 n , we obtain Letting n go to ∞, we have All we need to do now is show that each of the v 1,i (x) is γ-harmonic by itself, i.e. that (3.11) holds for each summation index separately.To do so, let us forget for a moment the part O 1 n and solve the exact analogue of (3.11).Since by assumption the ζ i are all different, we know that the vectors are linearly independent (written as a matrix, they give a Vandermonde matrix with determinant i<j (ζ j − ζ i ) = 0).Therefore, the system has no nontrivial solutions, and hence All that remains to do now is to see that the error term O 1 n in (3.11) does not change this.To do so, suppose now that (3.14) is not satisfied for some i.Then we know that there are arbitrarily large n such that (3.13) does not hold, i.e. its right-hand side takes a value ε n = 0.As the ζ i are roots of unity, there are only finitely many values which ε n can take for different n, so we cannot have convergence of ε n to 0. But then, choosing n large enough, (3.10) cannot hold either; a contradiction.Thus, (3.14) must hold, and we know that the v 1,i (B) are harmonic.By induction, applying the discrete Laplacian △ to both sides of (3.10), we argue in the same fashion that the v k,i (B) must be polyharmonic of degree k.
To show that q(B; n) = 0 for all n implies that v p,i (B) = 0 for all n, let us assume the opposite.So suppose that q(B; n) = 0 for all n, but that at the same time we have p, i such that v p,i (B) = 0. Assume our p to be minimal with this property.Then, we know that because otherwise this would be a contradiction to (3.9) for large n.But then we can utilize independence of the vectors (ζ k 1 , ζ k 2 , . . ., ζ k l ) for k = 0, . . ., r − 1 as before and arrive at a contradiction.Finally, to show that the v p,i are polyharmonic of order p, we can apply the operator △ p−1 to both sides of (3.9), notice that the first p − 1 terms vanish by assumption and then repeat the above argument.

Remarks:
• In the proof of Thm.3.1 we will see that in our case, for different i the v p,i (k, l) differ only by a factor of α −k i β −l i .Here, given a dominant saddle point (x 0 , y 0 ), and an associated one (x i , y i ), then α i , β i are the numbers such that (x i , y i ) = (α i x 0 , β i y 0 ).This allows us to essentially talk about only a single polyharmonic function v p for any given p; namely the one defined by the dominant saddle point.
• By studying the exact shape of the functions v p more closely (in particular their degree as polynomials), we will see in Section 4.1 that in the context of this article, we always have △v p+1 = v p + r p−1 , where r p−1 is some polyharmonic function of degree at most p − 1.This hints at the fact that the polyharmonic functions appearing in the asymptotic expansions are not arbitrary, but do have some form of recursive structure.
We now have all ingredients ready for the proof of Thm.3.1.
(3.18)By Cauchy's formula, we have with Γ 1,2 being closed curves around the origin.To evaluate the asymptotics of this integral, we utilize the saddle point method, as described for instance in [30,Chapter VIII].
The main idea is to conveniently choose our contours Γ 1 , Γ 2 such that they make the integral as easy to compute as possible.
With this in mind, suppose that (x 0 , y 0 ) is a dominant saddle point, and pick Γ 1 = {|x|= |x 0 |}, Γ 2 = {|y|= |y 0 |}.We know that the modulus of S(x, y) on Γ 1 ∩ Γ 2 is maximal; and the only other points where it attains the same value are the associated saddle points (x i , y i ).At any other point, |S(x, y)| will be strictly smaller -hence, when n goes to infinity, it suffices to compute the integral locally around our saddle points, since the rest will grow exponentially slower.We could hypothetically run into issues if N(x, y) were to be infinite, but we will see that this is not the case at our saddle points.For any other points, one can easily verify using the definition of the group that we can have infinite values of N(x, y) only at lines of the form {x = ζ i } or {y = χ i } for finitely many ζ i , χ i .It follows from the residue theorem, applied for the two coordinates iteratively, that these singularities can safely be neglected.It is then easy to check that, given an ε > 0, the set {(x, y) : |S(x 0 , y 0 ) − S(x, y)|< ε} is contained in a domain of the form {(x, y) As previously mentioned, in order to find the asymptotics, the rest of the integral can be neglected, as it will be exponentially smaller.Changing our coordinates to x = x 0 e is/ √ n , y = y 0 e it/ √ n , this corresponds to a region of the form s To find the asymptotics, it therefore suffices to compute the integrals where F j (s, t, k, l, n) is the expression obtained by substituting x = x j e is/ √ n , y = y j e it/ √ n , for (x j , y j ) the relevant saddle points (i.e.(x 0 , y 0 ) and the ones associated to it, as outlined in Section 2).We will see that, given any fixed m ∈ N, each such integral can be written in the form where Q(s, t) is some (positive definite) quadratic form and the p j (s, t) are polynomials.Assume w.l.o.g. that δ 1 ≥ δ 2 .For s, t > δ 2 √ n, we can see that the integral over the remaining part of R is exponentially small in n, and therefore we can consider the integral Consequently, all we need to do is to consider the latter integral for all saddle points with maximum absolute value (of which there are finitely many), and by computing all the expressions up to a fixed p j (s, t), we will then have obtained the asymptotics of (3.19) and at the same time shown (3.1).
In the following, we will proceed in two steps: first, we pick a dominant saddle point and show that, locally around this point, everything works out smoothly.Then, we pick any associated saddle point and show that, up to powers of roots of unity, nothing changes from the first step.
1. Suppose (x 0 , y 0 ) is a dominant saddle point, let γ := S(x 0 , y 0 ) and fix an m ∈ N. First, we show that 0 = |N(x 0 , y 0 )|< ∞.This is due to (x 0 , y 0 ) being a saddle point: we have ∂S ∂y (x 0 , y 0 ) = 0 and thus y 0 is the unique solution to S(x 0 , •) = S(x 0 , y 0 ).Therefore, φ(y 0 ) = y 0 , and in the same manner we can see ψ(x 0 ) = x 0 .It follows immediately that the alternating orbit sum of xy evaluated at (x 0 , y 0 ) is 0; in particular it is finite.Our next step is to show that we can rewrite the integrand (that is, the one in (3.19)) as in (3.22).To do so, we substitute x → x 0 exp is √ n =: e s , y → y 0 exp it √ n =: e t , and then separate the integral into three parts: S(e s , e t ) n , N(e s , e t ) and the denominator 1/e k+1 s e l+1 t (note that a power in the denominator vanishes due to the substitution rule).
(a) Part 1: S(e s , e t ) n As (x 0 , y 0 ) is a dominant saddle point of S(x, y) and therefore ∂S ∂x (x 0 , y 0 ) = ∂S ∂y (x 0 , y 0 ) = 0, we have a Taylor expansion of S(x, y) around (x 0 , y 0 ) of the form with A(s, t, n) = n −3/2 2m−1 j=0 a j (s,t) , and the a j homogeneous polynomials.We know that Q(s, t) := u 4 s 2 + vst + w 4 t 2 is a positive definite quadratic form (because our saddle point is a local minimizer of S(x, y) in R + × R + ).Consequently, we can write log S(e s , e t ) = log γ − Q(s, t) − A(s, t, n) with once again B(s, t, n) = n −3/2 2m−1 j=0 b j (s,t) Note in particular that this last factor is the only part which depends on the endpoint (k, l).
Multiplying the power series together and sorting them by powers of n, we obtain as a result the contribution of this saddle point to (3.22) of the form where the q ′ p (s, t, k, l) are polynomials in s, t, k, l.One can easily see that q ′ p (s, t, k, l) is homogeneous of degree p; for odd p the double integral therefore vanishes by symmetry.Thus we can rewrite (3.30) as with q p := q ′ 2p .The factor 1 n c in (3.30) stems from the fact that we obtain a factor of 1/n by the substitution rule, and that the integral might vanish for small values of p.In particular, from this we can conclude that the c appearing in (3.1) is integer.
2. Suppose now that we have another saddle point (x j , y j ) = (α j x 0 , β j y 0 ) associated to (x 0 , y 0 ), and pick ζ j such that S(x j , y j ) = ζ j S(x 0 , y 0 ) (notice that we then have |ζ j |= 1, as discussed in Section 2).Our goal is now to describe the series expansion of the numerator around (x j , y j ) using the one around (x 0 , y 0 ).We substitute as before x → e ′ s := x j exp is √ n , y → e ′ t := y j exp it √ n .Due to Lemma 3.1, we can now conclude that the series representation (w.r.t.n at ∞) around S(e ′ s , e ′ t ) is the same as the one of ζ j S(e s , e t ), and the representation of N(e ′ s , e ′ t ) is the same as α j β j N(e s , e t ).Lastly, the expansion of 1/e ′(k+1) s e ′(l+1) t clearly changes only by adding a factor of α j k+1 β j l+1 as well.Therefore, we can conclude that the contribution of the saddle point (x j , y j ) is the same as the one of (x 0 , y 0 ) up to a factor of ζ n j α −k j β −l j .By Lemma 3.2, we deduce that the v p (k, l) are indeed γ-polyharmonic of degree p (note that for our dominant saddle point we have α = β = ζ = 1).Because the only point in the construction where k, l appear is in (3.29), one finds that the v p (k, l) are, up to a factor of x −k 0 y −l 0 , bivariate polynomials in k, l.Due to Lemma 2.5, they are therefore polynomials if and only if the model has zero drift.The fact that c = π/θ (with θ the correlation coefficient as defined in (2.17)) follows from [22,Thm. 6].To see that γ = min x,y∈R + S(x, y), it suffices to remember that at the beginning of the proof we had γ = S(x 0 , y 0 ) for a dominant saddle point (x 0 , y 0 ), which is a minimizer of S by its definition in Section 2.
The construction used in the proof allows us to give some further properties of the polyharmonic functions appearing in the asymptotic expansion.
Corollary 3.1.The degree of the polynomial part of v p (k, l) (that is, without the factor of x −k 0 y −l 0 ) is c + 2p − 1. Proof.By looking once again at the proof of Thm.3.1 and in particular (3.29).From there, the statement follows immediately.Corollary 3.2.For any orbit-summable model, we have π/θ ∈ N.
Proof.From Thm. 3.1, we know that c = π/θ.From Cor. 3.1, we know that c + 2p − 1 is the degree of the polynomial part of v p (k, l), and therefore integer.Since 2p − 1 is integer, the statement follows. Remarks: • Thm.3.1 holds true in higher dimensions as well, and indeed the proof translates directly.The one difference lies in the powers of n which appear: by the substitution rule dx i ds i = c i e is i / √ n , one obtains additional factors of n −1/2 .Thus, for even dimensions the constant c in (3.1) will be integer, whereas for odd dimensions it will be in 1  2 + N.An example case for three dimensions is treated in Appendix B.
• When looking at models with large steps, the one thing that could go wrong is that the numerator might have a singularity at a saddle point.Usually, this seems not to be the case, and for a given model this condition is very easy to check.An example is treated in Appendix A.
• A table of the first three asymptotic terms for the 19 unweighted orbitsummable models is given in [43, App.D].

Example: the Gouyou-Beauchamps model
In this section, we will illustrate the result of Thm.3.1 by computing the asymptotics for the Gouyou-Beauchamps model, and in doing so find an explicit formula for the polyharmonic functions appearing therein.The Gouyou-Beauchamps Walk is defined by the step set {տ, ց, ←, →}.Its step polynomial is S(x, y) = y x + x y + 1 x +x, and solving ∂S ∂x = ∂S ∂y = 0 yields the four solutions (x, y) = (±1, ±1).We find that S(1, 1) = 4, S(−1, 1) = −4, S(1, −1) = S(−1, −1) = 0. Verifying that the second derivatives do not vanish, we therefore have the dominant saddle point (1, 1) and one other associated to it, namely (−1, 1), to consider.Note that the appearance of two saddle points is not at all surprising here, due to parity (or, in more general terms, periodicity) considerations: if the first coordinate of a given point is even, then we can only hit it after an even or odd number of steps.Therefore we can already expect at this point the asymptotics resulting from the saddle points to be precisely the same up to a factor of (−1) k+n , which coincides with the statement of Thm.3.1.Hence, we will only consider the dominant saddle point here.
The alternating orbit sum of xy can be checked to be Out integrand is therefore of the form where after letting x → e s := e is/ √ n , y → e t := e it/ √ n , by the substitution rule we will end up with The first factor of −1/n is entirely harmless; we will therefore proceed to compute each of the factors in the second fraction separately, precisely as in the proof of Thm.3.1.
1. Series representation of S(e s , e t ) n : First, we can compute (we let m := √ n for readability) with From here we obtain where with the i j positive integers.Finally, we can compute S(e s , e t ) n = exp (n log [S(e s , e t )]) where C(s, t, m) = k≥1 c j j k with with the i j again positive integers.
2. Series representation of N(e s , e t ): Using (3.32), we find that with 3. Series representation of 1/e k+1 s e l+1 t : Lastly, we have Overall, we obtain as product of the three factors computed above for j 1 , j 2 , j 3 nonnegative integers.In particular, we notice that q p is homogeneous of degree p in s, t, and of degree p in k, l.In order to compute the contribution up to order O 1 n j of this saddle point to the asymptotics, all we need to do now is compute This gives us an explicit formula for the asymptotics of this model.In particular, we can check directly that all coefficients of 1 m k for odd k vanish.By computing the integrals, we see that The harmonic function v 1 (k, l) was already computed some time ago in [5,45,20].

Periodicity
Considering the combinatorial context, it is clear that in any asymptotic expansion as in (3.1), the discrete harmonic function v 1 (k, l) will always be positive (if it does not vanish).When looking at the computations in [43, App.D], it appears as if there is an even stronger pattern; namely that v p (k, l) is positive for even p, and negative for odd p.It turns out, however, that this is not generally true; a counterexample is given for instance by computing enough terms in the expansion of the Simple Walk [17].
It is also a direct consequence of Thm.3.1 that the number of saddle points is closely tied to the periodicity of the model.If we have a single saddle point, then clearly our model is aperiodic; but the number of saddle points also corresponds directly to the periodicity.
Lemma 3.3.Suppose that our model is irreducible (that is, the semi-group generated by step set S is all of Z 2 ), and that it is m-periodic.Then we have exactly m − 1 saddle points s 1 , . . ., s m−1 associated to our dominant saddle point s 0 = (x 0 , y 0 ).The ζ i corresponding to these saddle points (see Section 2) are -in some order -the m-th roots of unity.
natural because the underlying combinatorial problem is highly symmetrical: any path from A to B corresponds to a path from B to A with reversed steps.This, and a similar statement for the higher order terms, can be formalized in Thm.4.1 below and the following Thm.4.2.
Theorem 4.1.Suppose that S is a step set satisfying the general assumptions stated in Section 2 and that S is orbit-summable.Then, with (x 0 , y 0 ) being a dominant saddle point, then there are constants c, r ∈ N, γ ∈ R + , and α i , β i , ζ i ∈ C as well as functions v p (k, l, u, v) such that for any m ∈ N we have The v p (k, l, u, v) are polynomials precisely if the drift is zero (else they contain exponential factors).In this case they are of bidegree c + 2p − 1 in both (k, l) and (u, v), and of total degree 2c + 4p − 2. Each v p (k, l, u, v) is multivariate polyharmonic of degree p.
Remark: the constants α i , β i , ζ i stem from the saddle points associated to the dominant one via s i = (x i , y i ) = (α i x 0 , β i y 0 ).Also note that despite the appearance of complex numbers, all sums end up being real.
Before proving Thm.4.1, we first show the following lemma, which is a natural extension of Lemma 3.2.
Lemma 4.1.Suppose that q(A, B; n) is a (combinatorial) quantity satisfying s∈S ω s q(A − s, B; n − 1) = q(A, B; n), s∈S ω s q(A, B + s; n − 1) = q(A, B; n), for all A, B ∈ Z ≥ × Z ≥0 , n ≥ 0. Assume furthermore that q(A, B; n) has an asymptotic representation of the form for all m, with the ζ i pairwise different and of modulus 1.Then, each v p,i is multivariate polyharmonic of order p.
Proof.The proof works in the very same manner as the proof of Lemma 3.2; at each step we can choose whether to apply the identity (4.2) or (4.3), leading to an additional instance of △ or △ respectively.
Proof (of Thm.4.1).Analogous to the proof of Thm.3.1.We note that the contribution of N(e ′ s , e ′ y ) for an associated saddle point changes by a factor of α u+1 i β v+1 i instead of a factor of α i β i , and only need to keep track of the coefficients u, v throughout, which however behave exactly in the same manner as k, l.Finally, the polyharmonicity properties are a direct of Lemma 4.1.

Remark:
• While the starting point (u, v) and the end point (k, l) of the walk end up playing a very similar role (see also Thm. 4.2), which is not at all surprising from a combinatorial point of view, this is not at all obvious from the proof: the role of (k, l) is very easily summarized as these coefficients only appear in the integrand as a factor of x −k−1 y −l−1 , the starting point does not appear directly as factor x u+1 y v+1 , but instead as its orbit sum.A priori, without the combinatorial interpretation, it does not seem to be obvious that both of these occurrences lead to a symmetrical role in the result.
In the following, we will want to describe the coefficients v p (k, l, u, v) appearing in Thm.4.1 more precisely.We will consider the drift zero case, which is however not a real restriction as we can transform any other model to a zero drift one using the Cramer transformation discussed in Section 2. The goal will be to prove the following theorem: where the a i,j are constants, the h j i are polyharmonic of degree (at most) i, and the g j i are adjoint polyharmonic of degree (at most) p + 1 − i.
In order to prove Thm. 4.2, it turns out to be very useful to have a polynomial basis of the space of polyharmonic functions for any given model.As the group is finite by assumption and π/θ ∈ Z by Cor.3.2, we can use the basis given in [42], which consists of sequences h m n of n-polyharmonic functions satisfying: 2. the h m n (k, l) are bivariate polynomials of increasing degree in both n and m: the degree will increase by 2 for each step in n, whereas it will increase by at least 2 and at most c + 1 for each step in m.
Proof (of Thm.4.2).Taking u, v as parameters, we can for each (u, v) write v p (k, l, u, v) as a sum of the basis functions h m n (k, l), and obtain Since we know that v p (k, l, u, v) is a bivariate polynomial of bidegree c+2p−1 in both (k, l) and (u, v) (which is also where the conditions n, m ≤ p and n + m ≤ p + 1 come from), the only thing we need to show is that g m n is adjoint polyharmonic of degree p + 1 − n for any m.To do so, we utilize Lemma 4.1.First, consider n = p.We then have therefore v 1 p (u, v) is adjoint harmonic.As the discrete Laplacians are linear, we can now proceed by induction, in each step applying the same argument as above to all terms which are polyharmonic (in (k, l)) of order n, i.e. the multiples of h 1 n , . . ., h p+1−n n .The statement follows. Remarks: • Thm.4.2 tells us in particular that we can write each coefficient v p as a sum of products of polyharmonic and adjoint polyharmonic functions, so that the degree of the former and latter adds up to at most p + 1.This can be viewed as an extension of [22,Thm. 6], where it is shown that v 1 is the product of a harmonic and an adjoint harmonic function (albeit in a much more general setting).
• By simple degree considerations, the only base function h j p appearing in v p (k, l, u, v) can be h 1 p .From this it follows that △v p+1 (k, l, u, v) = v p (k, l, u, v) + r p−1 (k, l, u, v), with r p−1 polyharmonic of degree at most p − 1.
• If the step set is symmetric in the sense that ω s = ω −s , then one can easily see that a i,j = a j,i due to the symmetry of the underlying combinatorial problem (q(A, B; n) = q(B, A; n), where q denotes the paths with reversed steps).This holds true for an appropriate choice of basis in some other cases as well; for examples of this, see [43, App.C].
• A priori it is not at all clear which elements of the basis h j i constructed in [42] actually appear in connection with some combinatorial problems, i.e. if this basis is combinatorially reasonably chosen.As we will see in Section 4.2, this seems to be the case.
• We can in fact give an upper bound on the number of summands appearing in the decomposition (4.5).Denote by (h m n ) the basis of harmonic functions from [42] as discussed above.First of all, we can see by degree considerations that for any given p, all (n, m) such that h m n can appear in the decomposition are contained in the subset {(n, m) : n + m ≤ p + 1}.If one then writes these h m n as a table, this gives us a triangular shape of size increasing with p.Along the lines n + m = k, we will have functions of degree at least 2(k − 1).We can now count the number of possible products of a given function h j i with a base function hn m of the adjoint Laplacian.In order not to exceed the degree of v p , we can multiply h 1 1 with the entire triangle of adjoint base functions.For h 2 1 and h 1 2 , we cannot multiply them with any hn m with m + n = p + 1, and so on.All in all, this gives us a maximum of p(p+1)(p+2)(p+3) 24 summands.This maximum is achieved e.g. for the Simple Walk for v 1,2,3 .Generally, the larger the value of π/θ, the less summands we will have (since the degrees of the h j i increase more quickly).
In the following, we will see what this decomposition looks like in case of the Simple Walk.The Gouyou-Beauchamps model and the Tandem Walk are treated in [43, App.C].

Example: the Simple Walk
In the following, to keep the expressions a bit shorter, we will give the expressions as after the substitution k → k − 1, l → l − 1 etc. (i.e.we have kl instead of (k + 1)(l + 1)).This corresponds to a shift of the quarter plane, where instead of Z ≥0 × Z ≥0 we now consider Z >0 × Z >0 .For the Simple Walk, with S = {→, ↓, ←, ↑}, we then have (after rescaling by multiplicative constants) By symmetry, we can pick the base functions hj i (u, v) = h j i (u, v) for the adjoint Laplacian.For the first three asymptotic terms with arbitrary starting and ending points, we obtain (again up to multiplicative constants) One can check that Cor.4.2 takes the form As the degree of h j i is truly increasing by only 2 whenever we increase either i or j by one, it turns out that we have indeed 1, 5 and 15 different summands respectively.Due to the symmetry of this model, the second and third equations can be simplified a bit: for v 2 , letting g 2 := 4h 1 2 + 2h 2 1 (clearly, g 2 is then biharmonic) gives us For v 3 , letting in the same manner g 3 := 192 5 h 1 3 + 64 5 h 2 2 + 4h 3 1 , we have While one can view the definition of g 2 and g 3 as purely a crutch to make the resulting expressions shorter, they do in fact give in a sense a natural decomposition of the v i : g 2 consists of the highest order terms in v 2 , while g 3 consists of those in v 3 .

Remarks:
• By the above, when taking a scaling limit, then we have where α i is an appropriate scaling constant, which means that the g i already give us all the terms which will not vanish in this kind of limit.
• By [18, Thm.Noticing that this representation looks almost the same as the one for the discrete case in Thm.4.1, it is natural to compare the functions v p to their continuous counterparts f p .For v 1 , for instance, by [22,Lemma 13] we know that we will have (after appropriate scaling by a constant) v 1 → f 1 .However, this is not at all clear for p > 1.For the Simple Walk, one can check that v 2 → f 2 , but this fails for p = 3: we have f 3 (k, l, u, v) = kluv(3k 4 + 6k 2 l 2 + 3l 4 + 22k 2 u 2 +6l 2 u 2 + 3u 4 + 6k 2 v 2 + 22l 2 v 2 + 6u 2 v 2 + 3v 4 ), whereas the scaling limit of v 3 (k, l, u, v) turns out to be kluv(3k 4 + 6k 2 l 2 + 3l 4 + 10k 2 u 2 +6l 2 u 2 + 3u 4 + 6k 2 v 2 + 10l 2 v 2 + 6u 2 v 2 + 3v 4 ), where the coefficients of k 2 u 2 and l 2 v 2 do not match.While it might seem a bit disheartening that the asymptotics of the discrete case do not converge to those of the continuous case 'arbitrarily well', one could just as well argue that this is in fact not something that can be expected: [22] only gives us first-order convergence, so in a sense anything going beyond the first-order terms is already more than we could've hoped for.

Outlook
• We know by [22,Thm. 6] that, for any (not necessary orbit-summable) model with a finite number of steps, the first order term of the asymptotics of q(k, l; n) behaves as in Thm.3.1.However, not much is known about the higher order terms.In this context, one could view Thm.3.1 as a partial generalization for the case of orbit-summable models with finite group, but it would be interesting to see whether something similar holds true for a more general class of models, which could for example also include models with large steps, as in e.g.[28].While it appears at first that this might fail due to the more complicated structure of the group, in [7] the notion of an orbit is extended to a more general context, which might make an approach using a saddle point method feasible.
• While in the zero-drift case there is a unique combinatorially relevant harmonic function for a given model (namely the positive one), there is no equivalent in the polyharmonic case, nor a reasonable combinatorial interpretation.It would be very interesting to know if there is such an interpretation, or at least, in some sense, a 'canonical' p-polyharmonic function.Intuitively, one may think such a canonical function should coincide in its scaling limit with the (continuous) polyharmonic function appearing in the corresponding continuous heat kernel, but as the higher order asymptotics of the discrete and continuous cases do not coincide (see Section 4.2), this might not be an ideal choice.Maybe the functions g 2 , g 3 , . . .as defined at the end of Section 4.2 could serve as candidates for such canonical representatives instead.
• While the saddle point method as used in this article works well to compute precise asymptotics for any given starting point given that the model has finite group and is orbit-summable, it seems that it would be hard to apply when one of these conditions does not hold.If the path counting function is algebraic (which, at least in the unweighted case, implies a finite group), sometimes it is possible to obtain an explicit expression which can then be utilized to extract asymptotics [11,12,31].In the infinite group case, however, things seem to be more complicated.In the one-dimensional setting, a probabilistic approach seems to work in order to show an asymptotic expansion similar to (1.3), using only moment conditions [21].
In two dimensions, one can tackle this problem using an explicit parametrization of the zero-set of the kernel via elliptic functions (in particular Jacobi theta-functions as in [44]).This method can then be utilized to show for some cases with infinite group (for instance the model with steps {←, ↑, →, ↓, ր}) that the asymptotics of the number of walks returning to the origin contain logarithmic terms [24], which makes an asymptotic expansion as in (3.1) impossible.However, it turns out that the coefficients still have a similar structure in the sense that they consist of polyharmonic functions.
• Finally, one could pose similar questions for different domains, be it a higher dimension or a different cone, for instance the three quarter plane.See for example [13,49,48,7,14,50,6].
We can find the dominant saddle point of S(x, y, z) to be at s 0 = 2 3/4 , 2 1/2 , 1 , with 7 others associated to it.Using the notation as in (3.1), we can check that and that we have A and Appendix B, examples of the method applied to a model with large steps and a three-dimensional model are given.Additional examples of the decomposition of the polyharmonic coefficients as in Section 4.2 for the Gouyou-Beauchamps model and the Tandem Walk as well as the first three terms of the asymptotics for all the 19 unweighted, orbit-summable models are given in the preprint of this article [43, App.C, D].

Theorem 4 . 2 .
The polynomials v p (k, l, u, v) in the zero drift case of Thm.4.1 each have a representation of the form