Asymptotics of Bivariate Analytic Functions with Algebraic Singularities

In this paper, we use the multivariate analytic techniques of Pemantle and Wilson to derive asymptotic formulae for the coefficients of a broad class of multivariate generating functions with algebraic singularities. Flajolet and Odlyzko (1990) analyzed the coefficients of a class of univariate generating functions with algebraic singularities. These results have been extended to classes of multivariate generating functions by Gao and Richmond (1992) and Hwang (1996, 1998), in both cases by immediately reducing the multivariate case to the univariate case. Pemantle and Wilson (2013) outlined new multivariate analytic techniques and used them to analyze the coefficients of rational generating functions. These same multivariate techniques can be used to analyze functions with algebraic singularities.


Introduction
For several decades, singularity analysis has been used to derive asymptotic formulae for coefficients of univariate generating functions.In 1990, for example, Flajolet and Odlyzko found asymptotics for a large class of univariate functions with algebraic singularities in [FO90].Examining the coefficients of multivariate generating functions is notoriously more difficult and technical.Pemantle and Wilson developed techniques to tackle multivariate rational generating functions in [PW13] and previous work, where they rely on the multivariate Cauchy integral, identifying and analyzing critical regions in the domain of integration that contribute to the integral's asymptotics through Morse theory.In this paper, we will look at the coefficients of H(x, y) −β , where H is an analytic function and β ∈ Z ≤0 is a real number.Under some assumptions about the zero set of H, we will find an asymptotic approximation for the coefficients [x r y s ]H(x, y) −β as r and s approach infinity with r s in a nearly-fixed ratio, as described in Theorem 1.
Flajolet and Odlyzko's 1990 results relied on using the Cauchy integral formula and explicit contour manipulations.Later in the 1990s, Gao, Richmond, Bender, and Hwang extended these results to classes of bivariate functions by temporarily fixing a variable and applying univariate results, which required special restrictions on the bivariate functions.(See Section 2.2 below for more details.)In this paper, we instead rely on the multivariate techniques that Pemantle and Wilson developed, manipulating the multivariate Cauchy integral formula directly.More details of these techniques are in Section 2.1 below.By using a combination of the Pemantle and Wilson techniques and the contour manipulations of the original Flajolet and Odlyzko work, we avoid using Morse theory.The algebraic singularities lead to manipulations of the torus on a Riemann surface instead of multidimensional complex space.However, this does not change the main methods of the asymptotic analysis, except requiring careful tracking of the argument of some expressions.
In Section 3, we state our main result (Theorem 1), which we prove in subsequent sections.Then, in Section 7, we look at examples of our results, including an application of Theorem 1 to a generating function that encodes properties of the stationary distributions of random colorings on the complete graphs, as found in [BCCG15].
An extended abstract of this paper, [Gre16], appeared in the proceedings of the 28th International Conference on Formal Power Series and Algebraic Combinatorics.

Historical Background
In this section, we provide some information about previous results in singularity analysis on which this work relies.

Multivariate Analytic Combinatorics of Rational Functions
In [PW13], Pemantle and Wilson outline a program which greatly extends the results of previous work on multivariate generating function analysis.Although many of the technical details of the program are not needed to prove the results in this paper, Pemantle and Wilson's work still lays the foundation for our approach.In the simplest case, Pemantle and Wilson begin with a rational function, F (z) = G(z)/H(z), where G and H are polynomials with real coefficients in the variables z 1 , . . ., z d , and where F (z) is analytic near the origin.We write z = (z 1 , . . ., z d ) and z r = z r 1 1 • • • z r d d .Then, F (z) has the series representation: The multivariate Cauchy integral formula tells us for r ∈ Z d ≥0 : Here, the torus is small enough that it does not enclose any singularities of F (z).The goal is to approximate z nr F (z) for some fixed unit vector r ∈ R d ≥0 as n approaches infinity, or if r has irrational components, to approximate coefficients [z nsn ]F (z) for large n with s n tending towards r.
To analyze the Cauchy integral, the torus T can be expanded into a cycle C which gets stuck on some chosen subset of the singularities of F (z) (which are the zeroes of H(z)), and expands beyond them elsewhere.Due to the z −r term in the integrand, we expect that as r → ∞, the integrand will decay exponentially faster in the regions of C away from the singularities of F , since the magnitude of z is larger in these regions.In this case, we can approximate the integral by analyzing the integrand near the singularities, since the rest of the integral decays too quickly to contribute to the asymptotics.However, the method of expanding T needs to be chosen carefully in order to ensure this works.
To expand T successfully, we need to minimize the maximum modulus of z −r along our contour C. The reason for this is as follows: we want to find a contour where the integrand attains its maximum modulus over some small interval, and then decays rapidly away from this interval.At a point where the maximum modulus is not minimized, the argument of the z −r term will oscillate rapidly as r tends to infinity, which leads to cancellation near the singularity.However, when the maximum modulus is minimized, we can approximate the integral in this region by using saddle point methods.
To minimize the maximum modulus, we consider the height function, h(z) := −r • Re log z.Although this excludes the contribution from F (z) in the integrand, F (z) is bounded on compact sets, so h still approximates the log modulus of the integrand as r approaches infinity in the direction of r.With the goal of expanding the torus T until it hits a singularity of F , we consider the values of h on V := {z : H(z) = 0}.On a cycle where the maximum of h is minimized, the points where the maximum of h is attained are saddle points of h.Thus, the critical points of h restricted to V will be candidates for the singularities that will contribute to the asymptotics.
To find the critical points of h, we consider a stratification of the space V, restricting our attention to critical points within a certain stratum S. When V is a smooth manifold near a critical point, the critical point is called smooth.In this case, the stratum S is of dimension d − 1, and we can find d − 1 equations (in addition to H = 0) that characterize the location of the smooth critical points: When H is a polynomial, the above critical point equations form a system of polynomial equations.In this case, Gröbner bases can help compute the critical points.In general, it is not necessarily true that all critical points will contribute to the leading term of the Cauchy integral.In this paper, we require that the critical points be minimal (described in Section 3 below), which guarantees that they do.To see examples of identifying minimal critical points, see Section 7 below.
After determining which critical points are candidates for contributing to the asymptotics, we still must expand the torus T into a cycle C which hugs V near these points.Goresky and MacPherson show in [GM88] how Morse theory can lead to an explicit description of the domain of integration near a critical point.In the case of generating functions without algebraic singularities, this machinery can be used to evaluate the residues of the integrals near each critical point, quickly leading to asymptotic expansions for the coefficients.However, in the case where H −β has algebraic singularities, we rely on specific homotopies of the contour, and hence we do not need to use Morse theory to determine the domain of integration.The portion of the contour near a particular critical point is called a quasi-local cycle.The asymptotics of the coefficients are thus given by a sum of integrals over these quasi-local cycles.Below, we find the leading-term asymptotics for the coefficients.However, these integral analyses can be used to find complete asymptotic expansions of coefficients, as in Raichev and Wilson's work in [RW08].

Asymptotics Involving Algebraic Singularities
In their 1990 paper [FO90], Flajolet and Odlyzko described how to compute the asymptotics of a class of univariate generating functions with algebraic singularities.They considered functions of the form, where α, γ, δ, and K are arbitrary real numbers, along with other related classes of functions.Their results differed from previous results both in the class of generating functions covered, and in their method of proof.Because we will use similar techniques in our proofs later, we take a moment to summarize their proof here.Flajolet and Odlyzko relied on the univariate Cauchy integral formula: Here, [z n ] g(z) represents the coefficient of z n in the power series expansion of g, and C is any positively-oriented contour around the origin which does not enclose any singularities of g(z).circle around the origin, the authors expanded C in hopes of finding a contour which is easier to analyze.In order to expand C, Flajolet and Odlyzko also require the extra assumption that f is analytic within the expanded contour C * .As C expands, it must avoid not only the singularity at 1, but also the branch cut emanating from this point.They expand the contour so it looks like the contour C * , as shown in Figure 1.Like a Hankel contour, this contour wraps around the branch cut of g and extends beyond the singularity at 1, although C * does not extend to infinity.From here, the contour is broken up into segments, γ 1 , γ 2 , γ 3 , and γ 4 .

Starting with any function
As n approaches infinity, f (z)/z n+1 , the integrand in the Cauchy integral formula, decays exponentially faster on γ 4 than it does on γ 1 .For this reason, the integral over γ 4 is negligible in the asymptotic expansion of [z n ] f (z).Likewise, the contribution along most of γ 2 and γ 3 is negligible, meaning that the asymptotics of [z n ] f (z) are controlled by the integrand near z = 1.However, near z = 1, f (z) = O(|z − 1| α ), which means that f is bounded along the contours near the critical point, leading to the bound, [z n ]f (z) = O(n −α−1 ).Flajolet and Odlyzko then extended their results to functions g(z) with the form in Equation (2).Before moving on to further developments with algebraic singularities, we highlight the connection between the proof outlined above and the results in this paper.To analyze bivariate generating functions, we deform a torus in two complex dimensions.Let (p, q) be a point contributing to the asymptotics of such a bivariate generating function.Then, in the proof below, one of the circles in the torus is expanded to a circle of radius |q|, while the other circle of the torus is expanded until it wraps around the singularity at p, similarly to how the Flajolet-Odlyzko contour wraps around the singularity in the univariate case.
After Flajolet and Odlyzko published their results in 1990, other researchers extended these results to classes of multivariate generating functions.Bender and Richmond, [BR83], had already considered the asymptotics of multivariate generating functions with poles in 1983.In 1992, Gao and Richmond, [GR92], considered classes of bivariate generating functions F (z, x) which are of a form they called algebraico-logarithmic, which includes some generating functions with algebraic singularities.These algebraico-logarithmic functions could be reduced to univariate generating functions where the results of Flajolet and Odlyzko can be applied.
Then, in his 1996 and 1998 papers, [Hwa96] and [Hwa98], Hwang expanded upon the multivariate results, using a probability framework and deriving large deviation theorems.In 1996, Hwang considered sequences of random variables {X n }.Assuming that the moment generating functions of the X n were of a particular form, Hwang proved a central limit theorem for {X n }.Then, he considered a class of bivariate generating functions P (w, z) such that after approximating [z n ] P (w, z) with Flajolet and Odlyzko's univariate results, [z n ] P (w, z) satisfied the same conditions he required previously of the moment generating functions of X n .Applying his central limit theorem gave asymptotic results for a new class of bivariate generating functions.In 1998, Hwang extended his results by using univariate saddle point methods to approximate integrals.

Main Result: Bivariate Analytic Functions with Algebraic Singularities
In this paper, our goal is to find the asymptotics of the coefficients of H(x, y) −β , where H is an analytic function with real coefficients and β ∈ R is not a negative integer.Let us summarize notation in a bivariate setting.Let V be the zero set of the analytic function, H(x, y), where H(0, 0) = 0. We will approximate the coefficients [x r y s ] H(x, y) −β for a fixed β ∈ R as r and s approach infinity with their ratio approaching a constant, λ.Critical points in the direction of λ = r+O(1) s (as r and s approach infinity) are defined by: The critical points are smooth if the gradient of H does not vanish on V at the critical points.Let D be the domain of convergence of the power series of H −β that converges around the origin, (0, 0).Then, a critical point (p, q) is called minimal if (p, q) ∈ ∂D.A collection of critical points is called strictly minimal if there are no other zeroes of H on ∂D.For notational convenience, we will represent partial derivatives with subscripts, so that for instance, H x = ∂H ∂x .We will apply heuristics from Section 2.1 to prove the following: Theorem 1.Let H be an analytic function with exactly n strictly minimal critical points {(p i , q i )} n i=1 , all of which are smooth and lie on the same torus T * .(Hence, , and let λ = r+O(1) s as r, s → ∞ with r and s integers.Define χ 1 , χ 2 , and M i as follows (where χ 1 and χ 2 depend on i): For all i, assume p i , q i , H x (p i , q i ), and M i are nonzero, and assume that the real part of −q 2 i M i is strictly positive.Define x −β P as the value of x −β defined by using a ray from the origin of C as the branch cut of the logarithm.In this definition, choose any ray such that H(x, y) −β P = H(x, y) −β in a neighborhood of the origin in C 2 (as defined by the power series of H −β ), and such that this ray does not pass through −p i H x (p i , q i ) for any i.Let ω i be the signed number of times the curve H(tp i , tq i ) crosses this branch cut in a counterclockwise direction as t increases, 0 ≤ t < 1.Then, the following expression holds as r, s → ∞: Here, the square root in the denominator is taken to be the principal root.
Unfortunately, for general H, the formula in Theorem 1 is messy, as we must find how many times the image of H wraps around the origin along the path connecting (0, 0) to each critical point (p i , q i ), and it is difficult to determine the sign of the square root.Luckily, in the case where H has only real coefficients and there is a single smooth strictly minimal critical point, we can simplify the formula.
Corollary 2. Let H be an analytic function with a single smooth strictly minimal critical point (p, q), where p and q are real and positive.Let H have only real coefficients in its power series expansion about the origin.Assume H(0, 0) > 0, and consider H −β for β ∈ R with β ∈ Z ≤0 .Also, define H −β here with the standard branch chosen along the negative real axis, so that H(0, 0) −β > 0. Let λ = r+O(1) s as r, s → ∞ with r and s integers.Define the following quantities: Assume that H x (p, q) and M are nonzero.Then, the following expression holds as r, s → ∞: In the above expression, −H x (p, q)p will be a positive real number, and (−H x (p, q)p) −β will also be a positive real number.Additionally, −2πq 2 M is positive, so the positive square root is taken.
Unfortunately, even in the simplified setting of Corollary 2, it is challenging to verify that a given critical point is strictly minimal.Section 7 below gives a couple examples where the Corollary can be applied.The RAGlib Maple package, [Saf15], can help verify these conditions computationally.
Generating functions of the form in the Theorem and the Corollary are expected to appear in several contexts.For example, there are many ways of extending the Catalan numbers to multidimensional arrays, like the Fuss-Catalan numbers.The generating function for the Catalan numbers has a square root, and in multivariate extensions, the generating functions are still algebraic.Another example is in counting RNA secondary structures with various structural features, called motifs.RNA secondary structures can be analyzed using stochastic context free grammars, as in [PH14], where multiple variables can be used to track more than one type of motif in a secondary structure simultaneously.In such a context, the Theorem above would give asymptotics on the number of secondary structures with motifs in a fixed ratio, as the number of nucleotides in the sequence approaches infinity.

Proof Set-Up
To prove Theorem 1, we analyze the multivariate Cauchy integral formula, Equation (1).When reduced to two dimensions, the formula becomes the following: We can immediately reduce to the case where there is only one strictly minimal critical point, which we will label (p 1 , q 1 ) = (p, q).(Correspondingly, the i in the subscripts of χ 1,i , χ 2,i , ω i and M i will be dropped.)This reduction is possible because the critical points are discrete, so that the contribution from each critical point (p i , q i ) to the asymptotics is given by a quasi-local cycle disjoint from the other quasi-local cycles, and thus the contributions can be analyzed independently and then summed.An outline of the analysis is as follows: 1.In Section 4.1, we find a change of variables into (u, v) coordinates so that the analytic function H(x, y) essentially behaves as a linear function in u, with some minor error terms in v.This change of variables also allows us to choose a simpler quasi-local cycle near the critical point, (p, q), where the u and v components of the contour are independent of each other.
2. In Section 4.2, we describe an appropriate expansion of the torus T in the Cauchy integral formula.Inspired by the contour from the univariate Flajolet-Odlyzko results described above, we choose a contour which similarly wraps around p in the u-coordinate, while passing directly through q in the v-coordinate.The description of the contour is technical to ensure that the contour does not cross over the singularities of H −β .Next, in Section 4.3, we verify that the region of the contour away from (p, q) does not contribute to the asymptotics.
3. With the set-up complete, we are ready to analyze the Cauchy integral.In Section 5, we apply the change of variables and contour deformations from the previous step, and then justify that the integrand of the Cauchy integral is approximately the product of a function in u and a function in v.This step is by far the most tedious, taking many lemmas to justify, and requiring analyses along each part of the quasi-local cycle.
4. Finally, the Cauchy integral is broken up into the product of two univariate integrals, which are analyzed in Section 6.The u integral is approximately a univariate Cauchy integral, and can be related to binomial coefficients.The v integral is a standard Fourier-Laplace type integral.Multiplying the approximations of these integrals gives the final result.

A Convenient Change of Variables
In order to approximate H(x, y) as a univariate linear function near the critical point (p, q), it turns out it is sufficient that the power series expansion of H has no constant term, linear term, nor quadratic term in one of its two input variables.To transform H into this form, we define the following change of variables: Here, χ 1 and χ 2 are as defined in Theorem 1 above.Write H as a power series in u and v: Since H(p, q) = 0, we have that d 00 = 0. Notice that when (x, y) = (p, q), we also have that (u, v) = (p, q).We can easily verify that d 01 = d 02 = 0 by checking some derivatives of H.

Determining the Quasi-Local Cycle
For now, assume that there is a unique critical point, (p, q).Recall that the original domain of integration in Equation (3) is a torus T around the origin which encloses no singularities of H −β (x, y).To decrease the magnitude of the integrand exponentially as r and s approach infinity, we expand the torus T towards the minimal critical point, (p, q).Because (p, q) is a strictly minimal critical point, there cannot be any zeroes between the origin and (p, q) that would otherwise obstruct the deformation.Hence, we can expand the domain of integration through a homotopy until it is near the critical point.Before expanding T , it is the product of a small x circle and a small y circle.Begin the deformation by expanding the y component to the circle, |y| = |q|.The y portion of the quasi-local cycle, C y , will be the part of this circle where y = qe iθ for |θ| ≤ θ y , where θ y > 0 is a small constant.Note that q is not necessarily real.This contour is pictured on the left in Figure 2. Now, for each y ∈ C y , we will expand the x circle until it approaches the zero set of H near p.When y is close to q, we will wrap the x contour around the zero set of H.However, when y is further away from q, we will expand the x contour less, so that it does not come into contact with the zero set of H.
More explicitly, since H x (p, q) = 0 and H is analytic, the implicit function theorem guarantees that we can parameterize the variety V = {(x, y)|H(x, y) = 0} by a smooth function G(y), so that H(p + G(y), y) = 0 for all y ∈ C y with θ y sufficiently small.So, for y = qe iθ with |θ| ≤ θy 2 , we choose the x contour appearing on the right in Figure 2. It is not necessarily true that this contour avoids the branch cut of H −β : to account for this, we can view all of our order-of-magnitude computations as if they are on the Riemann surface of H −β .Then, we can readjust our arguments accordingly when analyzing the final form of our Cauchy integral.The equations for the pieces of the contour are as follows: Here, x > 0 is a small positive constant.Note that later on, the change of variables into (u, v) coordinates will allow us to drop the corresponding G(y) term in the contour when sufficiently close to the critical point (p, q).Now, as |θ| increases with |θ| ≥ θy 2 , we would like to find an interpolation of the x quasi-local contour, shrinking it until it no longer wraps around the zero set of H.To do this, notice that when y = qe iθyt for t ∈ −1, − 1 2 ∪ 1 2 , 1 , |p + G(y)| > |p| uniformly, since (p, q) is a strictly minimal critical point of H. Therefore, we can find a δ > 0 so that |p + G(y)| > |p| + δ for every t ∈ −1, − 1 2 ∪ 1 2 , 1 .For y = qe iθ with |θ| > θy 2 , we linearly interpolate the radius |x − G(y)| in γ 4 and γ 5 from |p| + x to |p| + δ as |θ| increases from θy 2 to θ y , while correspondingly adjusting each other part of the contour, γ i , to form a closed curve as necessary.This gradually shrinks the quasi-local contour until it no longer wraps around the zero set V. We will show that the integrand is small along all parts of this contour, so the details of how the γ i intersect are not important.
This completes the description of a possible quasi-local contour near (p, q), but we will morph it slightly so that it is more convenient.Consider applying the change of variables given in Section 4.1.Since v = y, the v portion of the contour is identical to the y portion of the contour.Then, since u = x + χ 1 (v − q) + χ 2 (v − q) 2 , each contour γ i (y) is translated by χ 1 (v − q) + χ 2 (v − q) 2 , so that it retains its overall shape but is centered at a new location.Additionally, the parameterization of the zero set G(y) changes to a new parameterization κ(v) (discussed below) such that H(p + κ(v), v) = 0. κ(v) oscillates slower than G(y) near the critical point.(Here, H(u, v) = H(x, y), as defined in Equation (4).)This will allow us to approximate our quasi-local cycle as a product contour near the critical point (p, q), because we will show that we can drop κ(v) from the contour for y close enough to q.
In summary, the final quasi-local cycle C(p, q) (in (u, v)-coordinates) near the critical point (p, q) has three regimes.The contour is an arc in v, and wraps around the zero set of H in u.Let v = qe iθ .When θ ≤ r − 2 5 , the u-contour wraps exactly around the point, p, and this portion of the contour is a product contour.When r − 2 5 ≤ θ ≤ θy 2 , the contour instead wraps around the point p + κ(v) for a suitably chosen interpolation κ.Finally, if θ ≥ θy 2 , then the u-contour gradually shrinks as θ increases, until it no longer intersects the zero set of H at all.
Using the chain rule, we have Hu (p, q) = H x (p, q).Since H x (p, q) = 0 by assumption, and since H is analytic, the implicit function theorem guarantees that there exists a smooth parameterization κ(v) of the zero set of H, so that H(p + κ(v), v) = 0 for v sufficiently close to q. Investigating κ(v) a little further, we use the power series expansion of H about (p, q), given in Equation (4), along with the facts that d 00 = d 01 = d 02 = 0, to obtain the following: Thus, κ(v) = O(v − q) 3 .However, by comparing κ(v) to G(y), we find that G(y) = −χ 1 (y − q) − χ 2 (y − q) 2 + O(y − q) 3 .This provides more insight into this change of variables: in addition to allowing us to write H as a nice power series with some vanishing coefficients, the change of variables also describes V near (p, q).By converting the contour into (u, v)coordinates, we are able to stabilize the u contours, slowing down the movement of the zero set of H when it is parameterized by v.We take advantage of this slow-down by morphing our contour slightly, as described in the following paragraph.
In order to break the 2-dimensional Cauchy integral into two one-dimensional integrals, we need the quasi-local contour to be a product contour near the critical point, (p, q).To achieve this goal, we will need to break into two cases: when |θ| ≤ r − 2 5 and when |θ| > r − 2 5 , for v = qe iθ .Let us first analyze |v − q| in these cases: In the second line, we use the power series expansion for e iθ , which holds uniformly as < 1 r for r sufficiently large.Hence, when |θ| ≤ r − 2 5 , |κ(v)| = O r − 6 5 .Therefore, for r sufficiently large, the point p + κ(v) is always within the circle of radius 1 r about the point p, and we can morph our u-contour so that it is centered exactly around the point p instead of the point p + κ(v).Thus, we will drop κ(v) from the definitions of all the γ i when θ ≤ r − 2 5 , which means that the u contour no longer depends on v when θ ≤ r − 2 5 .(Note that this corresponds to a similar shift in the original (x, y)-coordinates, which can be computed explicitly to justify that the original torus T can be morphed locally to this new contour.)The portion of the contour where θ ≤ r − 2 5 will yield the dominating contribution to the integral asymptotically.In the other regime, when θ ≥ r − 2 5 , we cannot simply eliminate κ(v).Instead, let κ(v) be 0 when θ ≤ r − 2 5 , let it be κ(v) when θ ≥ r − 7 20 , and let it linearly interpolate between 0 and κ(v) when r − 2 5 ≤ θ ≤ r − 7 20 .We replace κ(v) with κ(v) in the definition of the quasi-local cycle.Note that κ(v) = O(v − q) 3 as v tends to q.We will use this condition much later in the proof.This completes the description of the quasi-local cycle.

Away from the Quasi-Local Cycle
Let us justify that the integral over the quasi-local cycle provides the main contribution to the asymptotics of the coefficients, and that the remainder of the domain of integration contributes negligibly in comparison.In order to do so, we will find a way to expand the torus T away from the quasi-local cycle so that the integrand decays exponentially faster here when compared to the quasi-local cycle.To begin, consider the case where there is only one strictly minimal critical point, (p, q).Formally, by strictly minimal, we mean the following: Here again, V = {(x, y)|H(x, y) = 0}.Consider the torus, T (p,q) := {x : |x| = |p|} × {y : |y| = |q|}.From this torus, remove an open neighborhood N of the point (p, q), where N is so small that the angular sectors of the torus that it covers in x and y are smaller than the angular sectors of the torus that C covers in x and y.That is, the y component of N should only consist of y values whose arguments | arg(y) − q| < c < θ y for some constant c > 0. Similarly, for each y ∈ C, the arguments of the x values in N should not vary from p + G(y) more than the arguments of the x values in C. Now, T (p,q) \N is a closed set which does not intersect the closed set, V. Thus, there are open sets dividing these two sets.This implies that there is a neighborhood of T (p,q) \N which does not intersect V.There is some δ * > 0 such that the x arc of T (p,q) \N can be expanded by δ * without hitting V.
Then, at every point of this new cycle away from the critical point (p, q), we have that |x| ≥ |p| + δ * .This forces the Cauchy integral to decay exponentially faster away from the critical point than it does near (p, q), proving that the asymptotic contribution to the integral cannot come from T (p+δ * ,q) .After expanding T (p,q) to T (p+δ * ,q) := {x : |x| = |p| + δ * } × {y : |y| = |q|}, notice that T (p+δ * ,q) \N can be connected to the quasi-local cycle C by adding two short lines at the ends of the x contours connecting the circle of radius δ * to the ends of γ 4 and γ 5 .Because these lines are contained entirely within the region near (p, q) where the implicit function theorem holds for G(y), the lines cannot hit any zeroes of H. Also, the magnitude of xy 1 λ along these lines is always greater than the magnitude of pq 1 λ because |x| > |p|, which means that these lines also do not contribute to the asymptotics of the integral, and may be ignored.

Approximating with a Product Integral
With our quasi-local cycle and change of variables defined, we are ready to begin analyzing the Cauchy integral formula, Equation (3).Writing the quasi-local cycle near (p, q) as C(p, q), we apply the change of variables to the Cauchy integral formula restricted to C(p, q) to obtain the following integral: Here, we used the fact that the Jacobian of the transformation is 1.Our goal now is to show that this integral is essentially a product integral.The following lemma describes this precisely.
Lemma 3. The integral in Equation (6) is asymptotically equivalent to the following: The above estimate holds as r, s → ∞ with λ = r+O(1) s .Here, C (p, q) is the portion of C(p, q) where |θ| ≤ r − 2 5 .Hence, C (p, q) is a product contour.
The proof of this lemma involves two types of statements: near the critical point, where |u − p| and |v − q| are both sufficiently small, we will argue that the integrands are asymptotically the same.Away from the critical point, where at least one of |u − p| or |v − q| is sufficiently large, we will show that both integrands are small, and hence do not contribute asymptotically to either integral.(In the second integral, we need only show that the integrand is small when |u − p| is large, since |v − q| is always small in C (p, q).)

When |u − p| and |v − q| are Small
In order to match the two integrands when |v − q| and |u − p| are small, we rewrite the original: Here, K and L are correction factors with the following definitions: Thus, our goal is to show that K and L are asymptotically equivalent to 1.We can analyze K and L along each part of the contour where (u − p) and (v − q) are small, and verify that in every case, they are asymptotically equal to 1 + o(1) as r, s → ∞ with λ = r+O(1) s .Explicitly, we will show this for u in γ 1 , and for the parts of γ 2 and γ 3 sufficiently close to the critical point.Recall that C y is the v-portion of the quasi-local contour, shown on the left in Figure 2, where v = qe iθ for |θ| ≤ θ y with θ y > 0 a small constant.
Lemma 4. Assume v ∈ C y with |θ| ≤ r − 2 5 .Also, assume that either u ∈ γ 1 , or that u ∈ γ 2 ∪ γ 3 with u = p + ωt r and t ≤ r 3 10 .Then, the following holds uniformly as r, s → ∞ with λ = r+O(1) For u ∈ γ 1 , |u − p| = 1 r .Thus, we have 1 − u p = O (r −1 ), and 1 − u p < 1 for r sufficiently large.Hence, we can expand as a uniformly convergent geometric series for all u ∈ γ 1 .This yields the following: Now, we can replace the numerator in the base of K by the expression in Equation ( 7) to obtain the following: Equation ( 8) holds uniformly for |θ| ≤ r − 2 5 and u in the region of the lemma, as r → ∞. Between γ 1 and the regions of γ 2 and γ 3 described in the lemma, 1 − u p = O r − 7 10 .Also, from Equation (5), |v − q| = O r − 2 5 .Plugging these facts into Equation (8) yields the following: 1 We replace the base of K with this new expression and use the Taylor series for the natural logarithm to obtain: Next, we prove the corresponding statement for L(u, v) on γ 1 and the parts of γ 2 and γ 3 sufficiently close to p. Lemma 5. Assume v ∈ C y with |θ| ≤ r − 2 5 .Also, assume either that u ∈ γ 1 , or that u ∈ γ 2 ∪ γ 3 with u = p + ωt r and t ≤ r Proof.Recall that H(u, v) has a particularly nice power series, given in Equation ( 4): In this series, we have the restrictions, d 00 = d 01 = d 02 = 0. Hence, we can express H in the following manner: Here, we choose f, g, and h to be any functions satisfying Equation ( 9) such that ) , and h(u, v) = O(v − q) 3 , each uniformly as (u, v) approaches (p, q).Also, we recall that in the power series expansion of H, d 10 = H x (p, q).We now plug Equation ( 9) into the definition of L: In the region described in this Lemma, we have the restrictions, 1 r ≤ |u − p| ≤ r − 7 10 , and |v − q| = O r − 2 5 .Thus, we obtain the following expressions: Each of these statements holds uniformly over the region in the lemma as r → ∞.Plugging these into Equation (10) above yields the desired result: This completes the proof that our integrand is essentially a product integrand near the critical point.It remains to show that the contributions away from the critical point are negligible.

When (u − p) or (v − q) is Big
In order to finish justifying that the integral in Equation ( 6) is approximately a product integral, we must show both that the integral in Equation ( 6) and the product integral in Lemma 3 are small away from the critical point, (p, q).We show this by approximating the magnitude of the integrands on the portion of the contours away from (p, q), and seeing that the contributions to the integral away from the critical point are negligible in comparison to the final approximation for the coefficients of H −β .We look at the integrals separately.Lemma 6.Let C(p, q) represent the portion of C(p, q) where at least one of the following conditions holds: |θ| > r − 2 5 or |u − q| ≥ r − 7 10 .Then, the following holds uniformly as r, s → ∞ with λ = r+O(1) s for some constant d > 0: .
Proof.We bound the terms of the integrand separately.First, recall the nice power series, is the parameterization of the zero set of H(u, v) near (p, q), and that κ(v) is the interpolation of this parameterization that is zero sufficiently close to (p, q) and κ(v) further from (p, q), as described in Section 4.2.Define ū by ū = u − p − κ(v) and v by v = v − q.H(p + κ(v), v) can be represented as follows: With this in mind, we plug the point (p + κ(v) + ū, v) into the power series of H and extract the portion corresponding to H(p + κ(v), v): for all |θ| ≤ θ y .Additionally, |ū| ≥ 1 r on all parts of C. Therefore, for x , θ y , and θ x sufficiently small, all terms in the expansion of H are negligible except d 10 ū.Since |ū| ≥ 1 r and since H is bounded on compact sets, we have the following bound uniformly on C: Now, we turn to the remaining part of the integrand.Using the relation, s = r λ + O(1) as r, s → ∞, we have the following: Here, we also took a factor of p −r out of (u − χ q (v − q) − χ 2 (v − q) 2 ) −r and a factor of We can expand ϕ as a bivariate power series: This equation holds uniformly as (u, v) approaches (p, q).M is a constant in terms of the derivatives of H, as defined in the statement of Theorem 1. Rewriting ϕ(u, v) in terms of κ, ū, and v = qiθ + O(θ) 2 gives the following as ū, v → 0: From here, our goal is to bound e −rϕ in magnitude.To do so, we will investigate the real part of ϕ.Let d = Re − q 2 M 2 , which is a strictly positive number by assumption.We break into cases now.
Case 1: Consider the case where u is close to the critical point p, in the sense that either u ∈ γ 1 and |θ| ≥ r − 2 5 , or u ∈ γ 2 or γ 3 and |ū| ≤ r − 7 10 but |θ| ≥ r − 2 5 .Either |ū| = 1 r for u ∈ γ 1 , which is much smaller than θ 2 , or 1 p ū is a strictly positive real number if u ∈ γ 2 ∪ γ 3 .In either case, u at worst does not increase the real part of ϕ, and we obtain the following: The above inequality holds for r sufficiently large and for x and θ y small enough.Thus, we have for r sufficiently large: Figure 3: α must be small when θ x is small.
Case 2: Consider the case where u ∈ γ 2 or γ 3 and |ū| ≥ r − 3 10 .(This case is only relevant when |θ| is small enough for γ 2 and γ 3 to be part of the contour.)For sufficiently small x and θ y , the O(ūθ) term is dominated by the ū term.The remaining θ terms are dominated by the θ 2 term, so these θ terms can only increase the real part of ϕ.Thus, the real part of ϕ is at least half the 1 p ū term, and we have the following for r sufficiently large: Plugging this into the exponential yields the following: Case 3: Now, consider the case where u ∈ γ 4 or γ Also, for θ x sufficiently small (depending on x and |p|), the following holds: This statement should be clear graphically: p , and consider Figure 3. Clearly, as θ x tends to zero, α approaches zero as well, verifying Equation (15).Combining Equation (14) and Equation (15), we have the following: Just like in Case 2, in the expansion of ϕ, the O(ūθ) term is dominated by the ū term, and the remaining θ terms are dominated by the θ 2 term, which only adds to the real part of ϕ.Hence, for x , θ x , and θ y sufficiently small we have: This yields: This decay is much greater than in the other cases: this is because here, (u, v) is bounded away from (p, q) by a constant amount.
In every case, we have the following bound for x , θ x , y , and δ sufficiently small and r sufficiently large: Finally, notice that for x , θ y , θ x , and δ sufficiently small, Plugging Equation ( 16) and Equation ( 17) back into Equation (12) gives the following: Recognizing that the entire domain of integration has size bounded by a constant, we combine Equation (11) and Equation (18) to get the desired result: . Now, we examine the corresponding statement for the product integral.
Proof.This statement can be proved by using the same method as in Lemma 6, by rewriting the integrand and examining the analogue of the function ϕ(u, v) in Lemma 6.The corresponding function has exactly the same form as ϕ(u, v) in Equation (13).
We have nearly completed the proof of Lemma 3.However, it is not yet clear that the bounds we have found away from the critical point are small compared to the value of the whole integral.It turns out that the exponential term in these bounds, e − d 2 r 1 5 , will ensure that these bounds are small compared to the integral overall.To show this, it remains to evaluate the asymptotic contribution of the product integral, which will simultaneously show that the contributions to the integral away from the critical point are negligible.

Analyzing the Product Integral
Lemma 3 has reduced our work to computing the following: We break it up into two univariate integrals, to be analyzed separately: Above, U is the u-projection of the contour, C , which resembles the x contour in Figure 2, but with G(y) = 0. V is likewise the v-projection, which is the set, v : v = qe iθ , |θ| ≤ r − 2

5
. We analyze each integral in lemmas below.
Lemma 8.The following holds uniformly as r, s → ∞ with λ = r+O(1) Here, ω is defined to be the signed number of times the curve H(tp, tq) crosses the branch cut in the definition of the function x −β P , as described in the statement of the Theorem.The main idea behind the proof of this lemma is that the remaining integral is almost a univariate Cauchy integral, but with a shrunken domain of integration.Thus, the integral is approximately equal to the coefficients of [H x (p, q) • (u − p)] −β , which can be estimated using the binomial theorem and Stirling's approximation.To account for the branch cut in the original function H(x, y) −β , we add a term e −β(2πiω) , where ω counts the number of times the image of H wraps around the origin as the input increases from (0, 0) to (p, q).An example is illustrated in Figure 4.
Proof.The contour U is comprised of the segments γ i for 1 ≤ i ≤ 5 in the case where |v − q| ≤ r − 2 5 .The endpoints of the contour, at the beginning of γ 4 and end of γ 5 , both have magnitude |u| = |p| + x .We can attach these endpoints to a portion of the circle {u : |u| = |p| + x } to form a closed cycle C u that wraps around the origin and contains no singularities of [H x (p, q) • (u − p)] −β .Because u −r−1 is exponentially smaller on the circle {u : |u| = |p| + x } than it is near the critical point p, we have: Now, we can use the Cauchy integral formula to evaluate this integral.However, we finally must worry about how the analytic continuation of H −β is defined.H(0, 0) is nonzero by assumption, and the values of H −β are defined near the origin of C 2 by the generating function itself.Separately from the analytic continuation of H −β that we have used up to this point, we choose a branch of the logarithm with the following properties: the branch must agree with H −β on some small neighborhood of the origin, and its branch cut must be a line from the origin that is not the line (t) = −tH x (p, q)p for t ≥ 0, for any of the critical points (p, q).Define x −β P as the value of x −β obtained by using this branch of the logarithm.
Consider the curve H(tp, tq) in C, with t ∈ [0, 1).This curve may wrap around the origin several times, and in particular, may cross the branch cut described above.Recall the bivariate power series for H(x, y): Plugging in our parameterization yields: The above equations are true as t → 1. Recall the following conditions: H x (p, q) = 0, and H y (p, q) = p λq H x (p, q).Plugging this into our computations above yields: Thus, as t tends to 1, the curve H(tq, tq) is essentially linear, with quadratic error.As long as the branch cut chosen above is not the line (t) mentioned above, the curve will only cross the branch cut finitely many times.Let ω be the signed number of times the curve H(tp, tq) crosses the branch cut in the counter-clockwise direction for t ∈ [0, 1).That is, every time the curve crosses the branch cut in the counter-clockwise direction, add 1 to ω, and every time it crosses in the clockwise direction, subtract 1 from ω.If the curve only touches the branch cut without crossing it, leave ω unchanged.As t approaches 1, we have shown that H behaves essentially like H x (p, q)(u − p), and we have traced how the argument changes as we expand the two-dimensional torus towards the critical point.Now, in order to revert the integral over C u back to the appropriate coefficient of H x (p, q)(u − p) by using the Cauchy integral formula, we must follow the image of H x (p, q)(u − p) from u = p back to the origin u = 0.As u follows the line from p to 0, the H x (p, q)(u − p) will follow the line in C from 0 to −pH x (p, q), the point whose power we are trying to determine.Because this straight line is (t), it will not cross the branch cut we chose above.Thus, ω already accounts for the total number of times the branch cut is crossed.Figure 4 shows an example of this setup.In this example, ω = 1, because H(tp, tq) crosses the branch cut once in the counter-clockwise direction.
In conclusion, we have the following: In the last line, we used Stirling's approximation to complete the proof.
We turn our attention to the other integral, and find its asymptotic contribution.
Lemma 9.The following holds uniformly as r, s → ∞ with λ = r+O(1) s : Here, the square root is taken to be the principal root.
This integral is nearly a Fourier-Laplace integral, again with a shrunken domain of integration, and thus can be estimated using standard Fourier-Laplace approximations.
Above, we define ψ as ψ(v) := log 1 − χ 1 (v−q)+χ 2 (v−q) 2 p + 1 λ log v q .Next, we expand ψ(v) as a Taylor series about v = q: Also, since v = qe iθ and |θ| ≤ r − 2 5 , we have that v q O(1) = 1 + o(1).Plugging these expressions into the integral and rewriting it in terms of θ gives us the following: The last line is true because O(θ) 3 = O r − 6 5 implies that e −rO(θ) 3 = 1+o(1).The remaining integral is nearly a Fourier-Laplace integral, but it has a shrinking domain of integration.We justify that this can be replaced by a domain of integration of constant size.Specifically, we aim to show that for some > 0 small enough, the following holds: To see this, notice that if |θ| ≥ r − 2 5 , then we have: Above, as before, we define d = Re − q 2 M 2 .Then, we have the following:

This justifies Equation (21).
To analyze the remaining integral, we use the standard saddle point approximation, which is proved in Theorem 4.1.1 in [PW13].The amplitude A(θ) = e iθ and the phase φ(θ) = − q 2 M 2 θ 2 are both analytic functions near θ = 0, and Re(φ) ≥ 0 on the interval [− , ], with equality only at θ = 0. Thus, we have: In the above expression, the square root is the principal root.Plugging Equation ( 22) into the remaining integral in Equation (20) finishes the proof.
Plugging the results of Lemma 8 and Lemma 9 into Equation ( 19) gives us the final result: Unfortunately, for general H, the formula in the Theorem becomes quite messy, as we must find how many times the image of H wraps around the origin along the path connecting (0, 0) to each critical point (p, q).Additionally, the sign of the square root in the formula can cause headaches.Luckily, in the case where H has only real coefficients and there is a single smooth strictly minimal critical point, we can simplify the formula, as seen in Corollary 2 above.The proof of the Corollary is simple from here.
Proof.Since H has real coefficients and p and q are positive real numbers, we must have that −H x (p, q)p is real.The line H(tp, tq) for 0 ≤ t ≤ 1 is real and can't pass through the origin since (p, q) is minimal.Also, the line from 0 to −H x (p, q)p is real, and it approximates H(tp, tq) for t near 1, as described in Section 6 and Figure 4 above.This would mean that −H x (p, q)p is in fact positive.Additionally, the line H(tp, tq) cannot wrap around the origin, which forces ω = 0 in the statement of the original Theorem.As a result, {(−H x (p, q)p) −β } P is positive.With this term positive, the only other term with an unknown sign is −2πq 2 M .However, knowing that −2πq 2 M is real, in order for the coefficients of H −β to be real at all, −2πq 2 M must be negative.Thus, −2πq 2 M is always positive (since the principal root is taken), which forces the whole formula to be positive always.

Examples 7.1 A Generating Function with Multinomial Coefficients
In order to see how Theorem 1 can be applied, we look at an artificial example where the computations are particularly simple, and where it is easy to verify that there is a strictly minimal critical point.Consider the coefficients x r y r of the generating function, F (x, y) = 1 (1 − x − y) 1 2 .
Using the Multinomial Theorem, we have that the coefficient [x r y s ]F (x, y) = − 1 2 r − 1 2 −r s .Thus, using Theorem 1, we will find an approximation for − 1 2 r − 1 2 −r r for large r.Here, the direction λ = 1.So, the critical point equations become the following: The unique solution to these equations is (p, q) = (1/2, 1/2).Because the x partial derivative of 1 − x − y never vanishes, this point is a smooth critical point.Additionally, it is easy to verify that (1/2, 1/2) must be a strictly minimal critical point: on the zero set of 1 − x − y, we have that y = 1 − x, and it is clear that if |x| ≤ 1 2 , then |y| ≤ 1 2 , with equality only when x = y = 1 2 .Thus, we can apply Corollary 2 to approximate [x r y r ]F (x, y).With H(x, y) = 1 − x − y and λ = 1, we compute the following: Plugging these into Corollary 1 with β = 1/2 yields: [x r y r ]F (x, y) ∼ 2 2r−1/2 πr .

Complete Graph Coloring
Using Corollary 2 above, we will find an asymptotic formula for the coefficients of the following bivariate generating function: This generating function describes the stationary distributions of a red/blue color-swapping algorithm on the complete graphs K r , where the coefficient of x r y s is proportional to the probability that s of the vertices of K r are blue in the stationary distribution, rescaled by a factor proportional to 2r−2 r .For more details, see [BCCG15].Because each y term is attached to an x of equal or greater power, the power series expansion of F will have no terms where the power of y is larger than the power of x.Thus, we will look at the asymptotics only the case where µ := λ −1 = s r ∈ (0, 1).(We switch to µ so that the range of possible directions is bounded.) The following paragraphs briefly describe how to check the conditions of Corollary 2 computationally.A Maple worksheet providing the code for these computations is available online:

Future Research
Possible future research directions include finding more complete asymptotic expansions for the coefficients of H −β , and also extending to more variables.Finding a more complete asymptotic expansion should be possible by analyzing the Cauchy integral using the same contour, but with more precise error handling.Extending the formula to more variables is challenging because the change of variables needed to approximate a multivariate function by a univariate linear function quickly becomes complicated.
Another research direction would be to look at other types of algebraic singularities.For example, one could study the coefficients of a function F (x, y) which is known to satisfy some polynomial equation but may not have the form H −β .One challenge in this case is to identify which singularities closest to the origin contribute to the asymptotics of the coefficients.If it is possible to determine the singularities, then the equation F satisfies can be used to estimate the behavior of F near the singularities.Assuming F is reasonably well-behaved near its singularities, the results in this paper could be applied like a transfer theorem to compute the asymptotics for the coefficients of F .In some cases, it is possible to find accurate computational estimates for the coefficients of F , which may be enough to determine where the contributing critical points for F are located.
Combining these results with other asymptotic techniques may yield stronger results and more complete asymptotic expansions, too.For example, creative telescoping methods take the generating function in question and find a partial differential equation that the function satisfies.By finding a basis of solutions to this differential equation, one can find complete asymptotic expansions to the coefficients of the generating function.Unfortunately, it is often difficult to find the correct coefficients of the solution to the PDE -this is referred to as the connection problem.However, if the leading-term asymptotics of the solution are known, the connection problem can often be solved.Thus, combining these creative telescoping methods with the first-order asymptotics results in this paper, one may be able to analyze generating functions without too many technical computations.

Figure 2 :
Figure 2: On the left, the y portion of the quasi-local contour.On the right, a close-up of the x portion of the quasi-local contour.