Large $z$ Asymptotics for Special Function Solutions of Painlev\'e II in the Complex Plane

In this paper we obtain large $z$ asymptotic expansions in the complex plane for the tau function corresponding to special function solutions of the Painlev\'e II differential equation. Using the fact that these tau functions can be written as $n\times n$ Wronskian determinants involving classical Airy functions, we use Heine's formula to rewrite them as $n$-fold integrals, which can be asymptotically approximated using the classical method of steepest descent in the complex plane.


Introduction and motivation
The six Painlevé differential equations, P I -P VI , have attracted a great deal of attention in the last decades. They feature in a large and increasing number of areas in mathematics, ranging from random matrix theory to integrable systems (continuous as well as discrete), orthogonal polynomials, partial differential equations and combinatorics. We refer the reader to standard references such as [12,18,20,23], as well as the Digital Library of Mathematical Functions, [17,Chapter 32], for more details.
The standard list of the Painlevé differential equations is the following: where α, β, γ, δ ∈ C and = d dz . We refer the reader to [23] for general properties of the equations and their solutions in the complex domain.
Generic solutions of P I -P VI are sometimes called Painlevé transcendents, or nonlinear special functions, and they cannot be expressed in terms of elementary or even classical special functions. As a consequence, asymptotic analysis of the solutions in the complex domain has been a key element in their study from the early days of the theory, and it Table 1: Classical special functions corresponding to special function solutions of the Painlevé equations P II -P VI .
has generated an enormous literature. In particular, the methodology of Riemann-Hilbert problems and the Deift-Zhou method of steepest descent have proved to be powerful approach in this direction, see for instance the monograph [18].
From the point of view of integrable systems, the Hamiltonian formulation gives more insight into the structure of the Painlevé functions than the standard form of the ODEs given above. It is known, see [12,21,22] and also [17, §32.6] and references therein, that each Painlevé equation can be written as a system: in terms of a suitable Hamiltonian H(q, p, z; v), where v is a vector of parameters present in the differential equation. Then the function q(z) satisfies the corresponding Painlevé equation, p(z) satisfies a related one (that sometimes, as in the case of P II , has independent interest), and the Hamiltonian itself, in some cases modified with elementary functions, satisfies the symmetric σ equation. Among all solutions of the Painlevé equations, families of rational and special function solutions are known to exist, for specific values of the parameters in the differential equations, for P II -P VI . These families of solutions are constructed in terms of elementary or classical special functions, and they appear in the theory of semiclassical orthogonal polynomials, see for example [34] and the recent monograph [35], and in random matrix theory, see [20,21,22].
The standard methodology to generate such special function solutions is to impose that, together with the corresponding Painlevé equation, the function u(z) satisfies a Riccati equation of the form This leads to conditions on the parameters of the Painlevé equations, and to specific forms for the coefficients p 0 (z), p 1 (z) and p 2 (z). Linearisation of this Riccati equation via u(z) = ±ϕ (z)/ϕ(z) produces a second order linear equation of classical type. The function ϕ(z), which is a polynomial or a classical special function, such as Airy, Bessel or hypergeometric, is usually called the seed function.
This procedure gives the first element of a family of rational or special function solutions. This family can be expanded applying Bäcklund transformations, which generate other solutions with modified parameters in a certain direction in the parameter space. Remarkably, in the case of special function solutions, the n-th element of the family can be written as an n × n Wronskian determinant involving the seed function ϕ(z). This fact was shown by Okamoto [30], and then exploited as a consequence of the Hamiltonian structure by Forrester and Witte [21,22], see also a summary in Clarkson [12]. Thus, we have a general formula: with initial values τ 0 (z) = 1 and τ 1 (z) = ϕ(z). Here D is a differential operator that depends on the particular Painlevé equation that we are considering. We observe that in general the function τ n (z) is analytic in the variable z (for P II and P IV ), or has fixed algebraic singularities (for P III , P V and P VI ), being a combination of classical functions and derivatives. Despite this explicit representation, the asymptotic behavior and the zero/pole patterns of these tau functions and related families of Painlevé special function solutions in C are still highly non-trivial. For rational solutions (of P II −P IV ) this has been studied by several authors, including Balogh, Bertola and Bothner [1], Bertola and Bothner [2], Bothner, Miller and Sheng [8], Buckingham [9], Buckingham and Miller [10,11]. In these references, the emphasis is in the large n behavior, with this parameter possibly coupled with others. Using the fact that these Wronskian determinants can be rewritten as Hankel/Toeplitz determinants constructed with the moments of a suitably chosen weight function, the machinery of Riemann-Hilbert problems and the Deift-Zhou method of steepest descent then gives the asymptotics of the underlying (complex) orthogonal polynomials and related quantities.
The aim of this paper is to illustrate how asymptotic approximations for τ n (z) corresponding to special function solutions can be obtained using similar techniques, but in a different asymptotic regime, namely with the size of the determinants n fixed and the variable z tending to infinity in the complex plane. The main steps of this approach are the following: 1. Relate the entries of the Wronskian determinant (1.4), which are classical special functions, with the moments of a suitably chosen weight function w(t, z), µ m (z) = Γ t m w(t, z)dt, m = 0, 1, 2, . . . , (1.5) which are defined on a certain contour Γ on the real line or in the complex plane.
for m ≥ 0 and some constant β, then we can write with D 0 (z) := 1.
It is worth mentioning that often there are several possible weight functions and contours, and some choices may be simpler and/or impose restrictions on parameters.
2. Apply the classical Heine's formula [ In the standard theory of orthogonal polynomials, the weight function is required to be positive along the integration contour. In the sequel this condition will not hold, as w(t, z) will be complex valued in general, but the integral formula remains valid without any technical difficulty, as long as the weight function decays suitably to make the integral convergent. We also point out that (1.8) is of great importance when connecting the theory of orthogonal polynomials with other areas such as random matrix theory, where the Hankel determinant can be connected with the partition function of a Hermitian random ensemble, we refer the reader to the work of Bleher [3] or Deift [16] for more details.
3. Apply the (classical) method of steepest descent to this n-fold integral, to obtain the leading asymptotic behavior in different sectors of the complex z plane. The details of this classical asymptotic method can be found in many references, for instance [7,29,31,33]; in the present context, the technical details will depend on several aspects: first, the different sectors where the variable z grows large in C, which will condition the deformation of Γ that is needed; secondly, the value of the constants in the seed functions, that may lead to very different asymptotic behaviors; lastly, the multidimensionality that makes it necessary to study different contributions to the total integral coming from the n variables. Also, rotational symmetries of the seed function, which are obtained from classical formulas for special functions, can be very useful in restricting the sectors in C where the asymptotic analysis is needed.
In this paper, we are interested in the special function (Airy) solutions of P II , as a relatively simple example of this methodology. These Airy solutions have been investigated recently by Clarkson in [13], but the asymptotic results are restricted to the real line, and only even values of n in the oscillatory regime are included. Asymptotic results for the seed case can also be found in [18,Chapter 11], and Its and Kapaev in [25,Proposition 4.3] characterize rational and Airy solutions of P II as those that do not exhibit asymptotic elliptic behavior in the complex plane. It is worth mentioning the numerical work done by Fornberg and Weideman in [19], which includes the special function solutions of P II and is extremely useful to understand the pole distribution of the Painlevé functions.
In the field of orthogonal polynomials, the connection of P II does not really arise from deformation of classical weights on the real line, as is the case for P III − P VI ; the fact that the exponential factor is cubic in the integral representation of the Airy functions makes it necessary to work with orthogonality that is defined on suitable contours in the complex plane. Nevertheless, these complex orthogonal polynomials have been recently studied by Clarkson, Loureiro and Van Assche in [15], and we note also the paper of Van Assche, Filipuk and Zhang [36], which uses a similar weight function, but in the context of multiple orthogonal polynomials.
In random matrix theory, the pure Ai case of P II solutions arises in the calculation of averages of powers of the characteristic polynomial in the GUE (Gaussian Unitary Ensemble), see [21,Proposition 28]. Also, the complex orthogonal polynomials mentioned before have been used in the asymptotic analysis (as n → ∞) of the partition function and free energy in the cubic Hermitian random matrix model, [4,5,6]. Remark 1.1. We observe that in principle the leading term in the asymptotic expansion could also be obtained using the Toda equation satisfied by the tau functions τ n (z): This kind of identity relates consecutive tau functions in a certain direction in the parameter space, and it is usually obtained as a consequence of the Hamiltonian formulation of the Painlevé equations, see for example Forrester and Witte [21,22]. Here D is a differential operator that can be just the derivative d/dz or a different one, C n is a constant (in z) and the original tau function may need modifications by elementary functions, see [20,Section 8.2.5].
Making a suitable ansatz of the leading term in the large z asymptotic expansion allows to construct a proof by induction, see [14,Proposition 5.2] for an example in P IV . However, we find that this methodology does not give a very precise estimate of the subleading terms or the order of the error terms, and also it poses problems in the oscillatory regime, where the leading term as z → ∞ is usually the result of combinations of different exponential contributions, that are difficult to keep track of in this Toda equation. For this reason, we prefer to calculate the expansions using steepest descent of the integral arising from Heine's formula. Once the general structure of the asymptotic expansion for τ n (z) is proved, (1.9) may be used to identify the coefficients therein.

Special function solutions of P II
Instead of focusing directly on solutions of the differential equations in (1.1), we use the Hamiltonian formulation, which provides a broader perspective of the different functions involved. For P II , the Hamiltonian is given by where p = p(z) and q = q(z), see [12, §3], [23, §18] whereas q(z) satisfies the P II equation: and σ(z) = H II (q, p, z) satisfies the symmetric S II equation: It is well known, see for instance [12, 6) and the seed function is a combination of Airy functions: where C 1 and C 2 are constants. The tau function then has the following form: with τ 0 (z) := 1. These tau functions can be used as building blocks to construct special function solutions of P II as well as of the associated equations. Namely, we have the following result: 13], Theorems 4 and 7). Consider the n × n determinant (2.8) for n ≥ 1, with τ 0 (z) := 1 and seed function given by (2.7), then the functions are special function (Airy) solutions of the P 34 , P II and the σ II equations respectively, see (2.2), (2.3) and (2.4), with parameter α = n − 1 2 . In view of the previous theorem, in this paper we will investigate the asymptotic properties of the τ functions directly, and then derive corresponding results for the solutions of the Painlevé equations.

Main results
In the study of the asymptotic behavior of the tau function (2.8) we distinguish, much like in the case of classical Airy functions, two regimes: non-oscillatory (exponential) and oscillatory (trigonometric). Furthermore, in the first case it is enough to obtain asymptotics in the sector | arg(−z)| < π 3 , since we can use the following rotational symmetries of the seed function (2.7): Lemma 3.1. The Airy seed function ϕ(z) given by (2.7) satisfies the following identities: with new constants

The tau function then satisfies
Proof. The proof is a straightforward manipulation of standard formulas for Airy functions, in particular see [17, 9.2.10-11]. The transformation for τ n (z) follows directly from the properties of the seed function.
Having this result, it is enough to study the tau functions in two different sectors of C: • The non-oscillatory sector | arg(−z)| < π 3 .
• The real axis, where oscillatory behavior occurs.
Once the asymptotic behavior is determined in these sectors, the results in the rotated ones follows directly by changing the constants suitably, according to Lemma 3.1.
Our main result about the asymptotic behavior of Airy-type solutions of P II is given in the two theorems below. We present the general result together with a number of consequences, and we highlight the case C 2 = 0 (when the seed function contains only Ai functions), which is particularly important in applications, because it has a distinguished asymptotic behavior, an extended non-oscillatory sector and particular importance in applications coming from orthogonal polynomials and random matrix theory.
(ii) If C 2 = 0, in the sector | arg(−z)| < π, From this theorem and the symmetry relations in Lemma 3.1, we can draw a number of consequences. Firstly, we can determine asymptotically free of poles regions in the complex plane for the special function solutions of P II : Corollary 3.3. If C 2 = 0, then the Airy solutions of P II are tronquée solutions (asymptotically free of poles) in the sectors Using these asymptotic expansions, we can determine the asymptotic behavior of the Painlevé functions in (2.9).

Oscillatory regime
The behavior in the oscillatory regime is particularly interesting. Figures 1 and 2 show the τ n (z) functions for different values of n, on the positive real axis, with C 2 = 0. It is apparent that in the even case there is a leading algebraic term, with oscillations of small amplitude superimposed, whereas in the odd case the leading term is itself oscillatory with increasing amplitude (except when n = 1). Theorem 3.5 makes this idea more precise.  Theorem 3.5. (Oscillatory regime). For n ≥ 1 and z ∈ R + , the function τ n (z) has the following asymptotic behavior as z → ∞: (B 2s,r (z) cos(ψ 2s,r (z)) + D 2s,r (z) sin(ψ 2s,r (z))) , (B 2s−1,r (z) cos(ψ 2s,r (z)) + D 2s−1,r (z) sin(ψ 2s,r (z))) , (3.11) where K n is given by (3.2), the phase function is (3.14) Remark 3.6. The general formula for the coefficients b n,r is cumbersome, but it can be easily evaluated with symbolic software, and in many cases the sums reduce to just a few terms. Furthermore, several important simplifications can be made if C 2 = 0 (pure Ai function in the seed): in this case, only the term p = n/2 survives, b As before, from this result about the tau function, we can deduce the asymptotic behavior of the Painlevé functions in the oscillatory regime.
Corollary 3.7. For n ≥ 1 and s ≥ 1, the functions q n (z), p n (z) and σ n (z) in (2.9) admit the following asymptotic expansions of the following form as z → ∞: (3.15) with coefficients given by (3.14) and phase function (3.12). Also,  and q n (z) = (−1) n z 2 cot(ψ 2s−1,s−1 (z)) + O(z −1 ), (3.20) We note that this is in accordance with the results in [13, Theorems 5 and 8], but it also extends the asymptotic results to the complex plane, and it includes the case of n odd in the oscillatory regime. Furthermore, it provides more precise estimates of the remainder terms.
It is worth mentioning that in the reference [28], Kuijlaars, Its andÖstensson consider the asymptotics (for real z) of a one parameter family of solutions of P 34 , depending on a parameter that is related to α in (2.2). These family is relevant in the analysis of critical edge behavior in unitary random matrix ensembles [27], and it includes the Airy solutions as a particular case. The results in (3.19) agree with Theorem 1.2 in [28], for the asymptotic behavior of the tronquée solutions of the P 34 equation 1 .
Remark 3.8. In order to give the leading asymptotic behavior it is enough to keep a few terms in the previous expansion, namely r = 0 in the non-oscillatory regime, and r = s (if n = 2s is even), and r = s − 1 (if n = 2s − 1 is odd) in the oscillatory regime. We have opted to present the full expansion because in order to examine the asymptotic behavior of solutions of the Painlevé equations σ n (z), p n (z) and q n (z), in particular in the oscillatory regime, higher order terms are needed. Also, these extra terms give exponential contributions that are important in the Stokes phenomenon for the Ai solution, we refer the reader to the discussion in Appendix A.

Proof of Theorem 3.2
We recall the classical integral representations of the Airy functions: and for m ≥ 0 we define the moments Then, we have ϕ(z) = µ Ai 0 (z) + µ Bi 0 (z), cf. (2.7), and as a direct consequence, The Wronskian determinant constructed from the seed function and the Hankel determinant for the weight function (4.2) are therefore related as follows: Consider now that we have a total of r Airy Ai integrals in the determinant, 0 ≤ r ≤ n, and consequently n − r Airy Bi integrals. Then where the Hankel determinants D n,r (z) can be written, following the standard theory [24, Corollary 2.1.3], as n-fold integrals: where r integrals are taken along the path corresponding to Ai and n − r along the one for Bi. We note that the integrand is symmetric in the n variables (the Vandermonde determinant may change sign when permuting variables, but it appears squared), therefore we can suppose, without loss of generality, that the first r integrals correspond to Ai and the last n − r integrals to Bi, and then sum over the possible permutations. We write and we make the change of variables (cf. [7, Section 7.3]) t k = 2 −1/6 √ ρ u k , k = 1, . . . , n. (4.8) Then we obtain D n,r (z) = 2 −1/6 ρ 1/2 n 2 n! n r where Γ α is any smooth infinite path that joins the sectors ∞e −πi/3 and ∞e πi/3 (for the Ai case) and the sectors −∞ and ∞e ±πi/3 (for the Bi case), and the phase function is We note that for any u k , the stationary points are given by the solutions of φ (u k ) = u 2 k − e iα = 0. This gives For |α| < π 3 , the main contribution to the integrals is given by the stationary point u − , since Re(φ(u k− )) > 0 and Re(φ(u k+ )) < 0. However, deformation through u k− whilst respecting the asymptotic behavior is only possible for the integrals corresponding to Bi, since Ai is exponentially recessive in this regime of α.
We use the freedom in choosing the path of integration to integrate along the paths of steepest descent through the points u k± . These paths, denoted here by Γ α , are implicitly given by the equation Im φ(u) = Im φ(u ± ), which in this case are difficult to describe exactly for general values of α, see the discussion in [7,Section 7.3]. However, we only need to know their existence around the stationary points. Note that the function φ(u) − φ(u ± ) is real-valued on Γ α , by construction. We add and subtract the value of the phase function at the following stationary point: and we obtain (4.13) We isolate the stationary points by fixing δ > 0 and two discs D(u k± , δ) of radius δ; then we define Γ ± α,δ = Γ α ∩ D(u k± , δ), that is, small portions of the steepest descent path around the stationary points. Then as z → ∞ in the sector that we are considering, we have where I n,r (z) = 2 −1/6 ρ 1/2 n 2 r!(n − r)!
with Γ = (Γ + α,δ ) r ∪ (Γ − α,δ ) n−r , and E n,r (z) is the remainder. It is important to note that we need to choose δ > 0 in such a way that the remainder E n,r (z) is exponentially small with respect to all the exponential terms that are present in (4.14), that is, Computing such a δ explicitly is complicated in general, but it is clear that for large enough |z|, such a choice is always possible, given that the phase function is real (and decaying) along the path of steepest descent.
In order to obtain large z asymptotics of I n,r (z), we need to add r contributions from the stationary point u + (which also correspond to the Ai integrals), and n−r contributions from the stationary point u − (which also correspond to the Bi integrals). We also note that the contribution from u − is doubled, because of the path used for the Bi function, which gives an overall 2 n−r factor.
We apply a final change of variables to transform the exponential terms in the integral into Gaussians: for each 1 ≤ k ≤ n, we define This change of variable is again implicit, but by construction it maps u k± to v k = 0. By perturbation, we can write as v k → 0, (4.17) These changes of variable take each contour Γ ± α,δ onto a segment [−ε, ε] on the real axis, for some ε > 0. Then, with an exponentially small error again, we can extend the integrals to the whole real axis, using a standard estimate: for ε > 1, C = ρ 3/2 √ 2 > 0, and f analytic, even and with at most polynomial growth, we have The last integral can be written in terms of incomplete Gamma functions (see [17, §8.2, §8.11] for definitions and asymptotics) and a similar argument can be used in n variables, taking the Vandermonde determinant as the function f . If we make ε large enough, we have a remainder that is exponentially small with respect to all the terms in (4.14).
Since we have r integrals around the stationary point u + and n − r around u − , then the differentials become where the error term has the form: for some coefficients a k and b j,k whose exact form is not relevant for the leading term in the asymptotic expansion. In order to study the Vandermonde determinant, we split it in three parts, separating those terms that combine two u + or two u − values and then a final one that mixes positive and negative stationary points: Taking only the leading terms, we have 1≤j≤r,r+1≤k≤n with error terms Ξ 2 (v), Ξ 3 (v), Ξ 4 (v) similar in structure to (4.18).
The first two products give two decoupled Selberg integrals, and the third one contributes a constant to the leading term. Writing everything together, we obtain: for d ≥ 1 and Re C > 0, in terms of the Barnes G function, see [17, §5.17]. Identifying C = ρ 3/2 / √ 2 and z = −ρ e iα , we have Combining the powers of 2 and π with the prefactor in (4.6) and summing over r, we arrive at the leading term in (3.3).
This calculation gives very relevant information about the remainder as well: the order of the error term in (4.21) comes from the fact that any term beyond the leading one in the differentials or in the Vandermonde will produce linear terms in the components v k (which integrate to 0 against the Gaussian because of symmetry) and then quadratic terms, which, using the formula will contribute to an error of order O(ρ −3/2 ) as ρ → ∞, since in our situation we have C = ρ 3/2 / √ 2. This is true for higher order terms as well, so each exponential level in (3.3) contains in fact a full asymptotic expansion in inverse powers of (−z) 3/2 , which appears in the coefficients A n,r (z) in (3.3).
In the case C 2 = 0 we have no Bi integrals, so we take r = n, and all integrals will involve the stationary point u + only. The calculation is analogous and the leading term follows from this substitution into (3.3), which gives (3.6); furthermore, following the standard asymptotic theory of the Airy Ai function, see for example [29,Section 4.7], the expansion holds in the larger sector | arg(−z)| < π.

Proof of Theorem 3.5
In the oscillatory regime, the asymptotic analysis is similar, but slightly more complicated because both Ai and Bi integrals need to be evaluated both at u + and at u − , which for z > 0 have exponential contributions with the same real part. We will highlight the main differences with respect to the non-oscillatory case.
Assuming that z > 0, we make the change of variables and the phase function becomes cf. (4.10), with stationary points By symmetry, we can consider without loss of generality the stationary point (4.12) again, and sum in r over the n r possible permutations. Then, we have where D n,r (z) is given, analogous to (4.9) and adding and subtracting the phase function at (4.12), by (5.5) As before, we isolate these points by fixing δ > 0 and two discs D(u k± , δ) of radius δ around the stationary points in each of the variables u k : where Γ is the corresponding path of steepest descent. Then as z → ∞ we have D n,r (z) = [I n,r (z) + E n,r (z)] , where I n,r (z) = 2 −1/6 z 1/2 n 2 r!(n − r)!
7) and E n,r (z) is the remainder.
In the analysis of the Vandermonde determinant, it is convenient to split the different cases depending on which stationary point is considered, using the index r, and not in terms of Ai and Bi functions. Note that in the non-oscillatory case, both ideas are equivalent, since each Airy function requires only one of the stationary points (u + for Ai and u − for Bi). In the oscillatory regime, however, we have (independently of the parameter r), p integrals of type Ai and n − p integrals of type Bi, with 0 ≤ p ≤ n, and any integral around the stationary point −i has different orientation depending if we are integrating along the Ai or the Bi contour: the contour for an integral along an Ai contour is oriented from right to left, and each one of them adds a −1 factor. In order to quantify this, for any given p, we need to count all possible configurations where we have q integrals of Ai type in the last n − r cases (where the point −i is taken into account) and therefore p − q integrals of Ai type in the first r cases, for any value of q from 0 to p. We define H n,r,p = min(p,n−r) q=max(0,p−r) r p − q n − r q (−1) q , (5.8) where the limits of summation are set so that all binomial numbers are well defined. Writing together all the contributions and summing over p, we have 6 Proof of Corollary 3.4 In the non-oscillatory regime, we can derive asymptotic expansions for the Painlevé functions σ n (z), p n (z) and q n (z) in (2.9) quite straightforwardly. Suppose first that C 2 = 0. Instead of working directly with the asymptotic expansion (3.3), it is simpler to pick the leading term therein, corresponding to r = 0, i.e. τ n (z) = K n A n,0 (z)e as |z| → ∞ in the sector | arg(−z)| < π 3 . Using this result, we can deduce the form of the asymptotic expansion for the function σ n (z), which is for K ≥ 1 with coefficients that follow from s n,k once again.
In the case C 2 = 0, we can use a similar argument, but with the form τ n (z) = K n A n,n (z)e − The negative exponential term naturally leads to changes in the coefficients, but the form of the asymptotic expansion is the same.
In the oscillatory regime, calculations are more delicate, and the main difficulty is to establish a general pattern for the asymptotic expansion, as we did before in (6.1). The reason for this is that in Theorem 3.5 there is a clear leading term, but as we expand further, different trigonometric functions will be involved. These higher order terms can be computed using a symbolic program, but for brevity we will give only the leading terms.
If C 2 = 0, this expression simplifies considerably, since only the term p = s survives, and then b