Large $z$ Asymptotics for Special Function Solutions of Painlevé II in the
Complex Plane

Large $\boldsymbol{z}$ Asymptotics for Special Function Solutions
of Painlevé II in the Complex Plane††This paper is a contribution to the
Special Issue on Painlevé Equations and Applications in Memory of Andrei
Kapaev. The full collection is available at
https://www.emis.de/journals/SIGMA/Kapaev.html

Alfredo DEAÑO

A. Deaño

School of Mathematics, Statistics and Actuarial Science, University of Kent,
UK<EMAIL_ADDRESS>
Received April 17, 2018, in final form September 22, 2018; Published online
October 03, 2018

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é 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.

Painlevé equations; asymptotic expansions; Airy functions

34M55; 34E05; 33C10; 30E10

_Dedicated to the memory of Andrei A. Kapaev_

##  1 Introduction and motivation

The six Painlevé differential equations, ${\rm P}_{\rm I}$–${\rm P}_{\rm 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, 17, 18, 21], as well as
the digital library of mathematical functions [30, Chapter 32], for more
details and the complete list of the Painlevé equations.

Generic solutions of ${\rm P}_{\rm I}$–${\rm P}_{\rm VI}$ are sometimes called
Painlevé transcendents, or nonlinear special functions, and they cannot be
expressed in terms of elementary or even classical special functions. However,
for specific values of the parameters in the differential equations, it is
known that families of rational and special function solutions exist for ${\rm
P}_{\rm II}$–${\rm P}_{\rm VI}$; these appear for instance in the theory of
semiclassical orthogonal polynomials, see for example [33] and the recent
monograph [34], and in random matrix theory, see [18, 19, 20].

These rational and special function solutions of Painlevé equations can be
constructed from a given _seed function_ $\varphi(z)$, with a suitable initial
value of the parameters; successive application of Bäcklund transformations
[12, Section 4], [30, Section 32.7] leads to a sequence of _tau functions_ ,
that as shown by Okamoto [28] (see also Forrester and Witte [19, 20] or
Clarkson [12]) have the form of $n\times n$ Wronskian determinants

$\displaystyle\tau_{n}(z)=\det\left(\frac{\textrm{D}^{j+k}}{\textrm{D}z^{j+k}}\varphi(z)\right)_{j,k=0,1,\ldots,n-1}$

with initial values $\tau_{0}(z)=1$ and $\tau_{1}(z)=\varphi(z)$. Here D is a
differential operator that depends on the particular Painlevé equation that we
are considering. The solution of the Painlevé equation (and of other
associated equations) can then be written directly in terms of these tau
functions.

In this paper we are interested in the special function solutions of the
second Painlevé equation, denoted ${\rm P}_{\rm II}$,

$\displaystyle q^{\prime\prime}=zq+2q^{3}+\alpha,\qquad\alpha\in\mathbb{C}.$
(1.1)

These special function solutions of ${\rm P}_{\rm II}$ can be written in terms
of standard Airy functions; more precisely, the general seed function is given
by

$\displaystyle\varphi(z)=C_{1}\operatorname{Ai}\big{(}{-}2^{-1/3}z\big{)}+C_{2}\operatorname{Bi}\big{(}{-}2^{-1/3}z\big{)},$
(1.2)

where $C_{1}$ and $C_{2}$ are constants, and the tau function is the $n\times
n$ Wronskian determinant

$\displaystyle\tau_{n}(z)=\det\left(\frac{{\rm d}^{j+k}}{{\rm
d}z^{j+k}}\varphi(z)\right)_{j,k=0,1,\ldots,n-1},\qquad n\geq 1,$ (1.3)

with $\tau_{0}(z):=1$. As shown in [13, Theorems 4 and 7], the functions

$\displaystyle p_{n}(z)=-2\frac{{\rm d}^{2}}{{\rm
d}z^{2}}\log\tau_{n}(z),\qquad q_{n}(z)=\frac{{\rm d}}{{\rm
d}z}\log\frac{\tau_{n-1}(z)}{\tau_{n}(z)},\qquad\sigma_{n}(z)=\frac{{\rm
d}}{{\rm d}z}\log\tau_{n}(z)$ (1.4)

are special function (Airy) solutions of the ${\rm P}_{\rm 34}$ equation

$\displaystyle p_{n}\frac{{\rm d}^{2}p_{n}}{{\rm
d}z^{2}}=\frac{1}{2}\left(\frac{{\rm d}p_{n}}{{\rm
d}z}\right)^{2}+2p_{n}^{3}-zp_{n}^{2}-\frac{\left(\alpha+\tfrac{1}{2}\right)^{2}}{2},$
(1.5)

the ${\rm P}_{\rm II}$ equation (1.1) and the symmetric
$\textrm{S}_{\textrm{II}}$ equation

$\displaystyle\left(\frac{\textrm{d}^{2}\sigma_{n}}{\textrm{d}z^{2}}\right)^{2}+4\left(\frac{\textrm{d}\sigma_{n}}{\textrm{d}z}\right)^{3}+2\frac{\textrm{d}\sigma_{n}}{\textrm{d}z}\left(z\frac{\textrm{d}\sigma_{n}}{\textrm{d}z}-\sigma_{n}\right)=\frac{1}{4}\left(\alpha+\tfrac{1}{2}\right)^{2},$
(1.6)

respectively, with $\alpha=n-\tfrac{1}{2}$.

The Airy solutions are used in the recent work of Clarkson, Loureiro and Van
Assche [15], in the asymptotic analysis of the partition function and free
energy in the cubic Hermitian random matrix model by Bleher, Deaño and
Yattselev [3, 4, 5], and in the study of multiple orthogonal polynomials with
a cubic potential in the complex plane by Van Assche, Filipuk and Zhang [35].
In random matrix theory, the pure $\operatorname{Ai}$ case of ${\rm P}_{\rm
II}$ special function solutions arises in the calculation of averages of
powers of the characteristic polynomial in the GUE (Gaussian unitary
ensemble), see [19, Proposition 28]. The asymptotic behavior of these Airy
solutions has 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 [17, Chapter 11], and Its and Kapaev in [23,
Proposition 4.3] characterize rational and Airy solutions of ${\rm P}_{\rm
II}$ as those that do not exhibit asymptotic elliptic behavior in the complex
plane.

The aim of this paper is to obtain large $z$ asymptotic approximations for
$\tau_{n}(z)$ corresponding to special function solutions of ${\rm P}_{\rm
II}$. As a direct consequence of (1.4), the analysis of $\tau_{n}(z)$ leads
directly to asymptotics for the Painlevé functions $p_{n}(z)$, $q_{n}(z)$ and
$\sigma_{n}(z)$. The methodology is related to ideas used for the large $n$
asymptotics for rational solutions (of ${\rm P}_{\rm II}$–${\rm P}_{\rm IV}$)
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]. Next we summarise the main steps:

  1. 1.

Relate the entries of the Wronskian determinant (1.3), which are Airy
functions, with the moments of a suitably chosen weight function $w(t,z)$,

$\displaystyle\mu_{m}(z)=\int_{\Gamma}t^{m}w(t,z){\rm d}t,\qquad
m=0,1,2,\ldots,$

where $\Gamma$ is a suitable contour on the real line or in the complex plane.
Namely, if

$\displaystyle\frac{{\rm d}^{m}}{{\rm
d}z^{m}}\varphi(z)=\beta^{m}\int_{\Gamma}t^{m}w(t,z){\rm
d}t=\beta^{m}\mu_{m}(z),$

for $m\geq 0$ and some constant $\beta$, then we can write

$\displaystyle\tau_{n}(z)=\det\big{(}\beta^{j+k}\mu_{j+k}(z)\big{)}_{j,k=0}^{n-1}=\beta^{n(n-1)}\det\left(\mu_{j+k}(z)\right)_{j,k=0}^{n-1}=\beta^{n(n-1)}D_{n}(z),$

for $n\geq 1$, 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. 2.

Apply the classical Heine’s formula [22, Corollary 2.1.3], [31, Section 2.2],
that gives a multiple integral representation for the Hankel determinant
obtained before:

$\displaystyle
D_{n}(z)=\frac{1}{n!}\int_{\Gamma^{n}}\Delta_{n}(t)^{2}\prod_{k=1}^{n}w(t_{k},z){\rm
d}t_{k},\qquad\Delta_{n}(t)=\prod_{1\leq j<k\leq n}(t_{k}-t_{j}).$

  3. 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 in one variable can be
found in many references, for instance [7, 27, 29, 32], and in the
multivariate context, we refer the reader to [6], [7, Chapter 8] or [16]. For
convenience, we detail the calculation instead of just writing the leading
term in the asymptotic expansion given in [16, equation (1.61)]; the main
technical details will depend on the different sectors where the variable $z$
grows large in $\mathbb{C}$, which will condition the deformation of $\Gamma$
that is needed, as well as the value of the constants in the seed function,
that lead to very different asymptotic behaviors.

######  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
$\tau_{n}(z)$. This is a general type of identity that relates consecutive tau
functions, and that for ${\rm P}_{\rm II}$ reads

$\displaystyle\frac{\tau_{n+1}(z)\tau_{n-1}(z)}{\tau_{n}(z)^{2}}=\frac{{\rm
d}^{2}}{{\rm d}z^{2}}\log\tau_{n}(z),\qquad n\geq 0,$ (1.7)

see other examples in [19, 20].

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 ${\rm P}_{\rm 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\to\infty$ 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 $\tau_{n}(z)$ is proved,
(1.7) may be used to identify the coefficients therein.

##  2 Main results

In the study of the asymptotic behavior of the tau function (1.3) 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)|<\frac{\pi}{3}$, since we can use the following rotational
symmetries of the seed function (1.2):

######  Lemma 2.1.

The Airy seed function $\varphi(z)$ given by (1.2) satisfies the following
identities:

$\displaystyle\varphi\big{(}{\rm e}^{\pm\frac{2\pi{\rm
i}}{3}}z\big{)}=\widetilde{C}_{1\pm}\operatorname{Ai}\big{(}{-}2^{-1/3}z\big{)}+\widetilde{C}_{2\pm}\operatorname{Bi}\big{(}{-}2^{-1/3}z\big{)},$

with new constants

$\displaystyle\widetilde{C}_{1\pm}=\frac{C_{1}}{2}{\rm e}^{\pm\frac{\pi{\rm
i}}{3}}+\frac{3C_{2}}{2}{\rm e}^{\mp\frac{\pi{\rm
i}}{6}},\qquad\widetilde{C}_{2\pm}=\frac{C_{1}}{2}{\rm e}^{\mp\frac{\pi{\rm
i}}{6}}+\frac{C_{2}}{2}{\rm e}^{\pm\frac{\pi{\rm i}}{3}}.$

The tau function then satisfies

$\displaystyle\tau_{n}\big{(}{\rm e}^{\pm\frac{2\pi{\rm i}}{3}}z\big{)}={\rm
e}^{\pm\frac{2\pi{\rm i}}{3}n(n-1)}\tau_{n}(z).$

###### Proof.

The proof is a straightforward manipulation of standard formulas for Airy
functions, in particular

$\displaystyle\operatorname{Bi}(z)={\rm e}^{-\frac{\pi{\rm
i}}{6}}\operatorname{Ai}\big{(}{\rm e}^{-\frac{2\pi{\rm i}}{3}}z\big{)}+{\rm
e}^{\frac{\pi{\rm i}}{6}}\operatorname{Ai}\big{(}{\rm e}^{\frac{2\pi{\rm
i}}{3}}z\big{)},$ $\displaystyle\operatorname{Ai}\big{(}{\rm
e}^{\mp\frac{2\pi{\rm i}}{3}}z\big{)}=\tfrac{1}{2}{\rm e}^{\mp\frac{\pi{\rm
i}}{3}}\left(\operatorname{Ai}(z)\pm{\rm i}\operatorname{Bi}(z)\right),$

see [30, formulas (9.2.10) and (9.2.11)]. The transformation for $\tau_{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 $\mathbb{C}$:

  * •

the non-oscillatory sector $|\arg(-z)|<\tfrac{\pi}{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 2.1.

Our main result about the asymptotic behavior of Airy-type solutions of ${\rm
P}_{\rm 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 $\operatorname{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.

###  2.1 Non-oscillatory regime

######  Theorem 2.2.

For $n\geq 1$, if we define

$\displaystyle
K_{n}=\frac{2^{-\frac{3n^{2}}{4}-\frac{n}{6}}}{\pi^{\frac{n}{2}}},$ (2.1)

then the function $\tau_{n}(z)$ has the following asymptotic behavior as
$|z|\to\infty$:

  1. $1.$

If $C_{2}\neq 0$, in the sector $|\arg(-z)|<\frac{\pi}{3}$,

$\displaystyle\tau_{n}(z)=K_{n}(-z)^{-\frac{n^{2}}{4}}\sum_{r=0}^{n}\mathbf{A}_{n,r}(z){\rm
e}^{\frac{\sqrt{2}}{3}(n-2r)(-z)^{\frac{3}{2}}},$ (2.2)

where

$\displaystyle\mathbf{A}_{n,r}(z)\sim(-z)^{\frac{3}{2}r(n-r)}\sum_{j=0}^{\infty}\frac{a_{n,r}^{(j)}}{(-z)^{3j/2}}$

and the leading coefficient is

$\displaystyle
a_{n,r}^{(0)}=(-1)^{\left\lfloor{r/2}\right\rfloor}2^{(n-r)(\frac{5r}{2}+1)}C_{1}^{r}C_{2}^{n-r}G(r+1)G(n-r+1),$

in terms of the Barnes $G$ function, $G(n)=\prod\limits_{k=0}^{n-2}k!$, see
[30, Section 5.17].

  2. $2.$

If $C_{2}=0$, in the sector $|\arg(-z)|<\pi$,

$\displaystyle\tau_{n}(z)=K_{n}(-z)^{-\frac{n^{2}}{4}}\mathbf{A}_{n,n}(z){\rm
e}^{-\frac{\sqrt{2}}{3}n(-z)^{\frac{3}{2}}}.$ (2.3)

From this theorem and the symmetry relations in Lemma 2.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 ${\rm
P}_{\rm II}$:

######  Corollary 2.3.

If $C_{2}\neq 0$, then the Airy solutions of ${\rm P}_{\rm II}$ are tronquée
solutions $($asymptotically free of poles$)$ in the sectors

$\displaystyle
S_{k}=\left\\{-\frac{\pi}{3}+\frac{2k\pi}{3}<\arg(-z)<\frac{\pi}{3}+\frac{2k\pi}{3},\,k\in\mathbb{Z}\right\\}.$

If $C_{2}=0$, then the Airy solutions of ${\rm P}_{\rm II}$ are tronquée
solutions in the sector

$\displaystyle S=\\{|\arg(-z)|<\pi\\}.$

This result proves a conjecture by Clarkson [13, p. 99]. Using these
asymptotic expansions, we can determine the asymptotic behavior of the
Painlevé functions in (1.4).

######  Corollary 2.4.

For $n\geq 1$, the functions $q_{n}(z)$, $p_{n}(z)$ and $\sigma_{n}(z)$ in
(1.4) admit asymptotic expansions of the following form:

  1. $1.$

If $C_{2}\neq 0$, as $|z|\to\infty$ with $|\arg(-z)|<\frac{\pi}{3}$,

$\displaystyle\sigma_{n}(z)=-\frac{n(-z)^{1/2}}{\sqrt{2}}-\frac{n^{2}}{4z}+\frac{\sqrt{2}n(4n^{2}+1)}{32(-z)^{5/2}}+\mathcal{O}\big{(}(-z)^{-4}\big{)},$
$\displaystyle
p_{n}(z)=-\frac{n}{\sqrt{2}(-z)^{1/2}}-\frac{n^{2}}{2z^{2}}-\frac{5n(4n^{2}+1)\sqrt{2}}{32(-z)^{7/2}}+\mathcal{O}\big{(}(-z)^{-5}\big{)},$
$\displaystyle
q_{n}(z)=\frac{(-z)^{1/2}}{\sqrt{2}}+\frac{2n-1}{4z}-\frac{\sqrt{2}(12n^{2}-12n+5)}{32(-z)^{5/2}}+\mathcal{O}\big{(}(-z)^{-4}\big{)}.$

  2. $2.$

If $C_{2}=0$, as $|z|\to\infty$ with $|\arg(-z)|<\pi$,

$\displaystyle\sigma_{n}(z)=\frac{n(-z)^{1/2}}{\sqrt{2}}+\frac{n^{2}}{4z}-\frac{\sqrt{2}n(4n^{2}+1)}{32(-z)^{5/2}}+\mathcal{O}\big{(}(-z)^{-4}\big{)},$
$\displaystyle
p_{n}(z)=\frac{n}{\sqrt{2}(-z)^{1/2}}-\frac{n^{2}}{2z^{2}}+\frac{5n(4n^{2}+1)\sqrt{2}}{32(-z)^{7/2}}+\mathcal{O}\big{(}(-z)^{-5}\big{)},$
$\displaystyle
q_{n}(z)=-\frac{(-z)^{1/2}}{\sqrt{2}}+\frac{2n-1}{4z}+\frac{\sqrt{2}(12n^{2}-12n+5)}{32(-z)^{5/2}}+\mathcal{O}\big{(}(-z)^{-4}\big{)}.$

###  2.2 Oscillatory regime

The behavior of $\tau_{n}(z)$ in the oscillatory regime is particularly
interesting. Figs. 1 and 2 show the $\tau_{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 2.5 makes
this idea more precise.

######  Theorem 2.5 (oscillatory regime).

For $n\geq 1$ and $z\in\mathbb{R}^{+}$, the function $\tau_{n}(z)$ has the
following asymptotic behavior as $z\to\infty$:

$\displaystyle\tau_{2s}(z)=K_{2s}z^{-s^{2}}\left[\mathbf{B}_{2s,s}(z)+2\sum_{r=0}^{s-1}\left(\mathbf{B}_{2s,r}(z)\cos(\psi_{2s,r}(z))+\mathbf{D}_{2s,r}(z)\sin(\psi_{2s,r}(z))\right)\right],$
$\displaystyle\tau_{2s-1}(z)=2K_{2s-1}z^{-\frac{(2s-1)^{2}}{4}}\sum_{r=0}^{s-1}\left(\mathbf{B}_{2s-1,r}(z)\cos(\psi_{2s-1,r}(z))+\mathbf{D}_{2s-1,r}(z)\sin(\psi_{2s-1,r}(z))\right),$

where $K_{n}$ is given by (2.1), the phase function is

$\displaystyle\psi_{n,r}(z)=(n-2r)\left(\frac{\sqrt{2}z^{3/2}}{3}+\frac{n\pi}{4}\right)$
(2.4)

and

$\displaystyle\mathbf{B}_{n,r}(z)\sim
M_{n,r}z^{\tfrac{3}{2}r(n-r)}\sum_{j=0}^{\infty}\frac{b_{n,r}^{(j)}}{z^{\tfrac{3j}{2}}},\qquad\mathbf{D}_{n,r}(z)\sim
M_{n,r}z^{\tfrac{3}{2}r(n-r)}\sum_{j=0}^{\infty}\frac{d_{n,r}^{(j)}}{z^{\tfrac{3j}{2}}},$

with

$\displaystyle M_{n,r}=(-1)^{r(n-r)}2^{\frac{5}{2}r(n-r)}G(r+1)G(n-r+1),$
(2.5)

and leading terms

$\displaystyle
b_{n,r}^{(0)}=\sum_{p=0}^{\left\lfloor{n/2}\right\rfloor}(-1)^{p}C_{1}^{2p}C_{2}^{n-2p}\sum_{q=\max(0,2p-r)}^{\min(2p,n-r)}(-1)^{q}{r\choose
2p-q}{n-r\choose q},$ $\displaystyle
d_{n,r}^{(0)}=\sum_{p=0}^{\lfloor{n/2}\rfloor-1}(-1)^{p-1}C_{1}^{2p+1}C_{2}^{n-2p-1}\sum_{q=\max(0,2p+1-r)}^{\min(2p+1,n-r)}(-1)^{q}{r\choose
2p+1-q}{n-r\choose q}.$ (2.6)

Figure 1: Plots of $\tau_{1}(z)$ (left), $\tau_{3}(z)$ (centre), $\tau_{5}(z)$
(right). In all cases $C_{2}=0$.

Figure 2: Plots of $\tau_{2}(z)$ (left), $\tau_{4}(z)$ (centre), $\tau_{6}(z)$
(right). In all cases $C_{2}=0$.

######  Remark 2.6.

The general formula for the coefficients $b_{n,r}^{(0)}$ and $d_{n,r}^{(0)}$
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 $\operatorname{Ai}$ function in
the seed): in this case, only the term $p=\lfloor{n/2}\rfloor$ survives,
$b_{2s-1,r}^{(0)}=d_{2s,r}^{(0)}=0$, and

$\displaystyle
b_{2s,r}^{(0)}=(-1)^{s}C_{1}^{2s}\sum_{q=\max(0,2s-r)}^{2s-r}(-1)^{q}{r\choose
2p-q}{n-r\choose q},$ $\displaystyle
d_{2s-1,r}^{(0)}=(-1)^{s}C_{1}^{2s-1}\sum_{q=\max(0,2s-r-1)}^{2s-r-1}(-1)^{q}{r\choose
2p+1-q}{n-r\choose q}.$

As before, from this result about the tau function, we can deduce the
asymptotic behavior of the Painlevé functions in the oscillatory regime.

######  Corollary 2.7.

For $n\geq 1$ and $s\geq 1$, the functions $q_{n}(z)$, $p_{n}(z)$ and
$\sigma_{n}(z)$ in (1.4) admit the following asymptotic expansions as
$z\to\infty$, excluding arbitrarily small but fixed neighborhoods of the poles
of the leading terms in the approximations:

  1. $1.$

If $C_{2}\neq 0$, then

$\displaystyle\sigma_{2s}(z)=\frac{s^{2}}{2z}+\frac{s}{2z}\frac{b_{2s,s-1}^{(0)}\sin(\psi_{2s,s-1}(z))-d_{2s,s-1}^{(0)}\cos(\psi_{2s,s-1}(z))}{b_{2s,s}^{(0)}}+\mathcal{O}\big{(}z^{-5/2}\big{)},$
$\displaystyle\sigma_{2s-1}(z)=\sqrt{\frac{z}{2}}\frac{-b_{2s-1,s-1}^{(0)}\sin(\psi_{2s-1,s-1}(z))+d_{2s-1,s-1}^{(0)}\cos(\psi_{2s-1,s-1}(z))}{b_{2s-1,s-1}^{(0)}\cos(\psi_{2s-1,s-1}(z))+d_{2s-1,s-1}^{(0)}\sin(\psi_{2s-1,s-1}(z))}+\mathcal{O}\big{(}z^{-1}\big{)},$

with coefficients given by (2.6) and phase function (2.4). Also,

$\displaystyle
p_{2s}(z)=-\frac{\sqrt{2}s}{\sqrt{z}}\frac{b_{2s,s-1}^{(0)}\cos(\psi_{2s,s-1}(z))+d_{2s,s-1}^{(0)}\sin(\psi_{2s,s-1}(z))}{b_{2s,s}^{(0)}}+\mathcal{O}\big{(}z^{-2}\big{)},$
$\displaystyle
p_{2s-1}(z)=z\left[1+\left(\frac{(b_{2s-1,s-1}^{(0)}\sin(\psi_{2s-1,s-1}(z))-d_{2s-1,s-1}^{(0)}\cos(\psi_{2s-1,s-1}(z))}{b_{2s-1,s-1}^{(0)}\cos(\psi_{2s-1,s-1}(z))+d_{2s-1,s-1}^{(0)}\sin(\psi_{2s-1,s-1}(z))}\right)^{2}\right]$
$\displaystyle\hphantom{p_{2s-1}(z)=}{}+\mathcal{O}\big{(}z^{-1/2}\big{)},$
(2.7)

and

$\displaystyle
q_{n}(z)=(-1)^{n}\sqrt{\frac{z}{2}}\frac{-b_{2s-1,s-1}^{(0)}\sin(\psi_{2s-1,s-1}(z))+d_{2s-1,s-1}^{(0)}\cos(\psi_{2s-1,s-1}(z))}{b_{2s-1,s-1}^{(0)}\cos(\psi_{2s-1,s-1}(z))+d_{2s-1,s-1}^{(0)}\sin(\psi_{2s-1,s-1}(z))}$
$\displaystyle\hphantom{q_{n}(z)=}{}+\mathcal{O}\big{(}z^{-1}\big{)},$ (2.8)

where $s=\left\lceil{\frac{n}{2}}\right\rceil$.

  2. $2.$

If $C_{2}=0$, then

$\displaystyle\sigma_{2s}(z)=\frac{s^{2}}{2z}-\frac{s}{2z}\sin(\psi_{2s,s-1}(z))+\mathcal{O}\big{(}z^{-5/2}\big{)},$
$\displaystyle\sigma_{2s-1}(z)=\sqrt{\frac{z}{2}}\cot(\psi_{2s-1,s-1}(z))+\mathcal{O}\big{(}z^{-1}\big{)}.$

Also,

$\displaystyle
p_{2s}(z)=\frac{\sqrt{2}s}{\sqrt{z}}\cos(\psi_{2s,s-1}(z))+\mathcal{O}\big{(}z^{-2}\big{)},$
$\displaystyle
p_{2s-1}(z)=\frac{z}{\sin^{2}(\psi_{2s-1,s-1}(z))}+\mathcal{O}\big{(}z^{-1/2}\big{)},$
(2.9)

and

$\displaystyle
q_{n}(z)=(-1)^{n}\sqrt{\frac{z}{2}}\cot(\psi_{2s-1,s-1}(z))+\mathcal{O}\big{(}z^{-1}\big{)},\qquad
s=\left\lceil{\frac{n}{2}}\right\rceil.$ (2.10)

We note that this is in accordance with the results in [13, Theorems 5, 8 and
9], 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 [26], Kuijlaars, Its and
Östensson consider the asymptotics (for real $z$) of a one parameter family of
solutions of ${\rm P}_{\rm 34}$, depending on a parameter that is related to
$\alpha$ in (1.5). This family is relevant in the analysis of critical edge
behavior in unitary random matrix ensembles [25], and it includes the Airy
solutions as a particular case. The results in (2.9) agree with Theorem 1.2 in
[26], for the asymptotic behavior of the tronquée solutions of the ${\rm
P}_{\rm 34}$ equation111In the notation of [26], we have $b=s_{2}=(-1)^{n}$,
in terms of the standard Stokes multipliers for ${\rm P}_{\rm II}$. Since the
restriction $b>0$ applies in the steepest descent analysis of the
Riemann–Hilbert problem, we only recover the case of even $n$ for Airy
solutions, with $b=1$ and $\beta_{0}=0$. We also note that the parameter
$\alpha$ in [26] corresponds to $-\alpha/2-1/4=-n/2$ in our notation, and a
change of variables is needed as well in the solutions of ${\rm P}_{\rm 34}$..

######  Remark 2.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 $\sigma_{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 $\operatorname{Ai}$ solution, we refer the reader to the discussion in
Appendix A.

##  3 Proof of Theorem 2.2

We recall the classical integral representations of the Airy functions

$\displaystyle\operatorname{Ai}(z)=\frac{1}{2\pi{\rm i}}\int_{\infty{\rm
e}^{-\pi{\rm i}/3}}^{\infty{\rm e}^{\pi{\rm
i}/3}}\exp\big{(}\tfrac{1}{3}t^{3}-zt\big{)}\mathrm{d}t,$
$\displaystyle\operatorname{Bi}(z)=\frac{1}{2\pi}\left(\int_{-\infty}^{\infty{\rm
e}^{-\pi{\rm i}/3}}+\int_{-\infty}^{\infty{\rm e}^{\pi{\rm
i}/3}}\right)\exp\big{(}\tfrac{1}{3}t^{3}-zt\big{)}\mathrm{d}t,$ (3.1)

for $z\in\mathbb{C}$, see [30, formulas (9.5.4) and (9.5.5)]. Bearing in mind
the form of the seed function (1.2), we take the weight function

$\displaystyle w(t,z)=\exp\big{(}\tfrac{1}{3}t^{3}+2^{-1/3}zt\big{)},$ (3.2)

and for $m\geq 0$ we define the moments

$\displaystyle\mu^{\operatorname{Ai}}_{m}(z)=\frac{C_{1}}{2\pi{\rm
i}}\int_{\infty{\rm e}^{-\pi{\rm i}/3}}^{\infty{\rm e}^{\pi{\rm
i}/3}}t^{m}w(t,z){\rm d}t,$
$\displaystyle\mu^{\operatorname{Bi}}_{m}(z)=\frac{C_{2}}{2\pi}\left(\int_{-\infty}^{\infty{\rm
e}^{-\pi{\rm i}/3}}+\int_{-\infty}^{\infty{\rm e}^{\pi{\rm
i}/3}}\right)t^{m}w(t,z){\rm d}t.$

Then, we have
$\varphi(z)=\mu^{\operatorname{Ai}}_{0}(z)+\mu^{\operatorname{Bi}}_{0}(z)$,
cf. (1.2), and as a direct consequence,

$\displaystyle\frac{{\rm d}^{k}}{{\rm
d}z^{k}}\varphi(z)=2^{-\frac{k}{3}}\big{[}\mu^{\operatorname{Ai}}_{k}(z)+\mu^{\operatorname{Bi}}_{k}(z)\big{]},\qquad
k\geq 0.$

The Wronskian determinant constructed from the seed function and the Hankel
determinant for the weight function (3.2) are therefore related as follows:

$\displaystyle\tau_{n}(z)=\det\left(\frac{{\rm d}^{j+k}}{{\rm
d}z^{j+k}}\varphi(z)\right)_{j,k=0}^{n-1}=2^{-\frac{n(n-1)}{3}}\det\big{(}\mu^{\operatorname{Ai}}_{j+k}(z)+\mu^{\operatorname{Bi}}_{j+k}(z)\big{)}_{j,k=0}^{n-1}.$

Consider now that we have a total of $r$ Airy $\operatorname{Ai}$ integrals in
the determinant, $0\leq r\leq n$, and consequently $n-r$ Airy
$\operatorname{Bi}$ integrals. Then

$\displaystyle\tau_{n}(z)=\frac{2^{-\frac{n(n-1)}{3}}}{(2\pi)^{n}}\sum_{r=0}^{n}\frac{C_{1}^{r}C_{2}^{n-r}}{{\rm
i}^{r}}D_{n,r}(z),$ (3.3)

where the Hankel determinants $D_{n,r}(z)$ can be written, following the
standard theory [22, Corollary 2.1.3], as $n$-fold integrals:

$\displaystyle
D_{n,r}(z)=\frac{1}{n!}\int\cdots\int\Delta_{n}(\mathbf{t})^{2}\prod_{k=1}^{n}w(t_{k},z){\rm
d}t_{k},\qquad\Delta_{n}(\mathbf{t})=\prod_{1\leq j<k\leq n}(t_{k}-t_{j}),$
(3.4)

where $r$ integrals are taken along the path corresponding to
$\operatorname{Ai}$ and $n-r$ along the one for $\operatorname{Bi}$, and we
denote $\mathbf{t}=(t_{1},t_{2},\ldots,t_{n})$. 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
$\operatorname{Ai}$ and the last $n-r$ integrals to $\operatorname{Bi}$, and
then sum over the possible permutations.

Using the scaled variables $\mathbf{t}=2^{-1/6}\sqrt{\rho}\mathbf{u}$, where

$\displaystyle z=-\rho{\rm e}^{{\rm i}\alpha},\qquad\rho\geq
0,\qquad|\alpha|<\frac{\pi}{3},$

cf. [7, Section 7.3], we obtain

$\displaystyle
D_{n,r}(z)=\frac{\left(2^{-1/6}\rho^{1/2}\right)^{n^{2}}}{n!}{n\choose
r}\int_{\Gamma_{\alpha}}\cdots\int_{\Gamma_{\alpha}}\Delta_{n}(\mathbf{u})^{2}\exp\left(\frac{\rho^{3/2}}{\sqrt{2}}\phi(\mathbf{u})\right)\prod_{k=1}^{n}{\rm
d}u_{k},$

where the phase function is

$\displaystyle\phi(\mathbf{u})=\sum_{k=1}^{n}\left[\frac{1}{3}u_{k}^{3}-{\rm
e}^{{\rm i}\alpha}u_{k}\right],$

and $\Gamma_{\alpha}$ is any smooth infinite path that joins the sectors
$\infty{\rm e}^{-\pi{\rm i}/3}$ and $\infty{\rm e}^{\pi{\rm i}/3}$ (for the
$\operatorname{Ai}$ case) and the sectors $-\infty$ and $\infty{\rm
e}^{\pm\pi{\rm i}/3}$ (for the $\operatorname{Bi}$ case).

Clearly, the gradient and Hessian of this function are

$\displaystyle\nabla\phi(\mathbf{u})=\big{(}u_{1}^{2}-{\rm e}^{{\rm
i}\alpha},\ldots,u_{n}^{2}-{\rm e}^{{\rm i}\alpha}\big{)},\qquad
H\phi(\mathbf{u})=2\operatorname{diag}(u_{k})_{k=1,\ldots,n}.$

The stationary points

$\displaystyle\mathbf{u}^{*}=(u_{\pm},\ldots,u_{\pm}),\qquad u_{\pm}=\pm{\rm
e}^{{\rm i}\alpha/2}.$

(with any combination of signs) are non-degenerate, since $\det
H\phi(\mathbf{u}^{*})\neq 0$.

For $|\alpha|<\frac{\pi}{3}$, the main contribution to each
$\operatorname{Bi}$ integral is given by the stationary point $u_{k-}$, since
$\textrm{Re}(\phi(u_{k-}))>0$ and $\textrm{Re}(\phi(u_{k+}))<0$, and it will
appear doubled because of the two paths joining $-\infty$ and $\infty{\rm
e}^{\pm\pi{\rm i}/3}$ in (3.1); for the $\operatorname{Ai}$ integrals, the
relevant stationary point is $u_{k+}$, since path deformation through $u_{k-}$
would change the asymptotic behavior (from exponentially decreasing to
exponentially increasing). Therefore, for the asymptotic analysis we need to
consider stationary points of the form

$\displaystyle\mathbf{u}^{(r)}=(\underbrace{u_{+},\ldots,u_{+}}_{r\,\text{times}},\underbrace{u_{-},\ldots,u_{-}}_{n-r\,\text{times}}).$
(3.5)

We add and subtract the value of the phase function at this point, to obtain

$\displaystyle
D_{n,r}(z)=\frac{\left(2^{-1/6}\rho^{1/2}\right)^{n^{2}}}{r!(n-r)!}\exp\left(\frac{\rho^{3/2}}{\sqrt{2}}\phi(\mathbf{u}^{(r)})\right)$
$\displaystyle\hphantom{D_{n,r}(z)=}{}\times\int_{\Gamma_{\alpha}}\cdots\int_{\Gamma_{\alpha}}\Delta_{n}(\mathbf{u})^{2}\exp\left(\frac{\rho^{3/2}}{\sqrt{2}}\big{[}\phi(\mathbf{u})-\phi\big{(}\mathbf{u}^{(r)}\big{)}\big{]}\right)\prod_{k=1}^{n}{\rm
d}u_{k}.$

By analyticity of the integrand in all the variables $u_{k}$, the precise
structure of these global paths is not relevant for the analysis, as long as
they connect the correct sectors in the complex plane, since we can deform the
contours sequentially in the different variables. We will use this freedom to
integrate along the paths of steepest descent through the points
$\mathbf{u}^{(r)}$. These paths, denoted here by $\Gamma_{\alpha}$, are
implicitly given by the equation

$\displaystyle\operatorname{Im}\phi(\mathbf{u})=\operatorname{Im}\phi\big{(}\mathbf{u}^{(r)}\big{)},$

and then it follows that the function
$\phi(\mathbf{u})-\phi(\mathbf{u}^{(r)})$ is real-valued on $\Gamma_{\alpha}$.
The paths of steepest descent are difficult to describe globally for general
values of $\alpha$, see the discussion in [7, Section 7.3], however, for
asymptotic approximations we only need their existence around the stationary
points.

We isolate the stationary points by fixing $\delta>0$ and two discs
$D\big{(}\mathbf{u}^{(r)},\delta\big{)}$ of radius $\delta$; then we define

$\displaystyle\Gamma_{\alpha,\delta}=\Gamma_{\alpha}\cap
D\big{(}\mathbf{u}^{(r)},\delta\big{)},$

that is, small portions of the steepest descent path around the stationary
points. Then as $z\to\infty$ in the sector that we are considering, we have

$\displaystyle
D_{n,r}(z)=\exp\left(\frac{\sqrt{2}}{3}(n-2r)(-z)^{3/2}\right)[I_{n,r}(z)+E_{n,r}(z)],$
(3.6)

where

$\displaystyle
I_{n,r}(z)=\frac{\left(2^{-1/6}\rho^{1/2}\right)^{n^{2}}}{{r!(n-r)!}}\int_{\Gamma_{\alpha,\delta}}\Delta_{n}(\mathbf{u})^{2}\exp\left(\frac{\rho^{3/2}}{\sqrt{2}}\big{[}\phi(\mathbf{u})-\phi\big{(}\mathbf{u}^{(r)}\big{)}\big{]}\right)\prod_{k=1}^{n}{\rm
d}u_{k},$

and $E_{n,r}(z)$ is the remainder.

It is important to note that we need to choose $\delta>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 (3.6), that is,
$E_{n,r}(z)=o\big{(}\exp\big{(}{-}\frac{\sqrt{2}n}{3}(-z)^{3/2}\big{)}\big{)}$
as $|z|\to\infty$. Computing such a $\delta$ 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.

We apply a final change of variables to transform the exponential terms in the
integral into Gaussians, which is a particular case of Morse lemma in the
literature [16, Chapter 1, Section 2]: for each $1\leq k\leq n$, we define

$\displaystyle\phi(\mathbf{u})-\phi\big{(}\mathbf{u}^{(r)}\big{)}=-\sum_{k=1}^{n}v^{2}_{k}.$

This change of variable can be written by components, and as $v_{k}\to 0$ we
have

$\displaystyle u_{k}=u_{+}+{\rm i}{\rm e}^{-\frac{\alpha{\rm
i}}{4}}v_{k}+\frac{{\rm e}^{-{\rm
i}\alpha}}{6}v_{k}^{2}+\mathcal{O}\big{(}v_{k}^{3}\big{)},\qquad
u_{k}=u_{-}+{\rm e}^{-\frac{\alpha{\rm i}}{4}}v_{k}+\frac{{\rm e}^{-{\rm
i}\alpha}}{6}v_{k}^{2}+\mathcal{O}\big{(}v_{k}^{3}\big{)}.$

This maps the contour $\Gamma_{\alpha,\delta}$ onto
$[-\varepsilon,\varepsilon]^{n}$ on the real axis, for some $\varepsilon>0$.
Then, with an exponentially small error again, we can extend the integrals to
the whole real axis, using a standard estimate: for $\varepsilon>1$,
$C=\frac{\rho^{3/2}}{\sqrt{2}}>0$, and $f$ analytic, even and with at most
polynomial growth, we have

$\displaystyle\int_{-\infty}^{\infty}f(v){\rm e}^{-Cv^{2}}{\rm
d}v-\int_{-\varepsilon}^{\varepsilon}f(v){\rm e}^{-Cv^{2}}{\rm
d}v=2\int_{\varepsilon}^{\infty}f(v){\rm e}^{-Cv^{2}}{\rm d}v\leq
2\int_{\varepsilon}^{\infty}f(v){\rm e}^{-Cv}{\rm d}v.$

The last integral can be written in terms of incomplete Gamma functions (see
[30, Sections 8.2 and 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 $\varepsilon$ large enough, we have a remainder
that is exponentially small with respect to all the terms in (3.6).

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:

$\displaystyle\Delta_{n}(\mathbf{u})^{2}=\prod_{1\leq j<k\leq
r}(u_{k}-u_{j})^{2}\prod_{r+1\leq j<k\leq n}(u_{k}-u_{j})^{2}\prod_{1\leq
j\leq r,r+1\leq k\leq n}(u_{k}-u_{j})^{2}.$

This term, together with the differentials, becomes

$\displaystyle\Delta_{n}(\mathbf{u})^{2}\prod_{k=1}^{n}{\rm
d}u_{k}=2^{2r(n-r)}{\rm e}^{\frac{\pi{\rm i}}{2}r^{2}-\frac{\alpha{\rm
i}}{4}(n^{2}+6r(n-r))}$
$\displaystyle\hphantom{\Delta_{n}(\mathbf{u})^{2}\prod_{k=1}^{n}{\rm
d}u_{k}=}{}\times\prod_{1\leq j<k\leq r}(v_{k}-v_{j})^{2}\prod_{r+1\leq
j<k\leq n}(v_{k}-v_{j})^{2}\Xi(\mathbf{v})\prod_{k=1}^{n}{\rm d}v_{k},$

where

$\displaystyle\Xi(\mathbf{v})=1+\sum_{k=1}^{n}a_{k}v_{k}+\sum_{j,k=1}^{n}b_{j,k}v_{j}v_{k}+\cdots,\qquad|v|\to
0,$

for some coefficients $a_{k}$ and $b_{j,k}$ whose exact form is not relevant
for the leading term in the asymptotic expansion. This leads to two decoupled
Selberg integrals, and writing everything together we obtain

$\displaystyle
I_{n,r}(z)=\big{(}2^{-1/6}\rho^{1/2}\big{)}^{n^{2}}\exp\left(\frac{\sqrt{2}(n-2r)}{3}(-z)^{3/2}-\left(\frac{n^{2}}{4}+\frac{3r(n-r)}{2}\right)\alpha{\rm
i}+\frac{r^{2}\pi{\rm i}}{2}\right)$
$\displaystyle\hphantom{I_{n,r}(z)=}{}\times
2^{(n-r)(2r+1)}\frac{S_{r}(z)S_{n-r}(z)}{r!(n-r)!}\big{(}1+\mathcal{O}\big{(}\rho^{-3/2}\big{)}\big{)},$
(3.7)

where

$\displaystyle S_{d}=\int_{\mathbb{R}^{d}}\Delta_{d}(v)^{2}\prod_{k=1}^{d}{\rm
e}^{-Cv_{k}^{2}}{\rm
d}v_{k}=\frac{(2\pi)^{d/2}}{(2C)^{d^{2}/2}}\prod_{k=1}^{d}k!=\frac{(2\pi)^{d/2}}{(2C)^{d^{2}/2}}G(d+2),$
(3.8)

for $d\geq 1$ and $\operatorname{Re}C>0$, in terms of the Barnes $G$ function,
see [30, Section 5.17]. Identifying $C=\rho^{3/2}/\sqrt{2}$ and $z=-\rho\,{\rm
e}^{{\rm i}\alpha}$, we have

$\displaystyle
I_{n,r}(z)=2^{-\frac{5n^{2}}{12}+\frac{n}{2}+(n-r)(\frac{5}{2}r+1)}\pi^{\frac{n}{2}}(-z)^{-\frac{n^{2}}{4}-\frac{3}{2}r(n-r)}G(r+1)G(n-r+1)$
$\displaystyle\hphantom{I_{n,r}(z)=}{}\times\exp\left(\frac{\sqrt{2}(n-2r)}{3}(-z)^{3/2}+\frac{r^{2}\pi{\rm
i}}{2}\right)\big{(}1+\mathcal{O}\big{(}(-z)^{-3/2}\big{)}\big{)}.$

Combining the powers of $2$ and $\pi$ with the prefactor in (3.3) and summing
over $r$, we arrive at the leading term in (2.2).

This calculation gives very relevant information about the remainder as well:
the order of the error in (3.7) 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

$\displaystyle\int_{\mathbb{R}}v_{k}^{2}{\rm e}^{-Cv_{k}^{2}}{\rm
d}v_{k}=\frac{1}{2C}\int_{\mathbb{R}}{\rm e}^{-Cv_{k}^{2}}{\rm
d}v_{k},\qquad\operatorname{Re}C>0,$

will contribute to an error of order $\mathcal{O}\big{(}\rho^{-3/2}\big{)}$ as
$\rho\to\infty$, since in our situation we have $C=\rho^{3/2}/\sqrt{2}$. This
is true for higher order terms as well, so each exponential level in (2.2)
contains in fact a full asymptotic expansion in inverse powers of
$(-z)^{3/2}$, which appears in the coefficients $\mathbf{A}_{n,r}(z)$ in
(2.2).

In the case $C_{2}=0$ we have no $\operatorname{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 (2.2), which gives (2.3); furthermore, following the standard asymptotic
theory of the Airy $\operatorname{Ai}$ function, see for example [27, Section
4.7], the expansion holds in the larger sector $|\arg(-z)|<\pi$.

##  4 Proof of Theorem 2.5

In the oscillatory regime, the asymptotic analysis is similar, but slightly
more complicated because both $\operatorname{Ai}$ and $\operatorname{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.

We start the calculation with (3.3) and (3.4) as before. Assuming that $z>0$,
we make the change of variables $\mathbf{t}=2^{-1/6}z^{1/2}\mathbf{u}$, and we
have

$\displaystyle
D_{n,r}(z)=\frac{\left(2^{-1/6}z^{1/2}\right)^{n^{2}}}{n!}{n\choose
r}\int_{\Gamma_{\alpha}}\cdots\int_{\Gamma_{\alpha}}\Delta_{n}(\mathbf{u})^{2}\exp\left(\frac{\rho^{3/2}}{\sqrt{2}}\phi(\mathbf{u})\right)\prod_{k=1}^{n}{\rm
d}u_{k},$ (4.1)

where the phase function is now

$\displaystyle\phi(\mathbf{u})=\sum_{k=1}^{n}\left[\frac{1}{3}u_{k}^{3}+u_{k}\right],$

with stationary points $\mathbf{u}^{*}=(\pm{\rm i},\ldots,\pm{\rm i})$. By
symmetry, we can consider without loss of generality the stationary point
$\mathbf{u}^{(r)}$ in (3.5) again, and sum in $r$ over the ${n\choose r}$
possible permutations. Adding and subtracting the phase function at (3.5), we
obtain

$\displaystyle
D_{n,r}(z)=\frac{\left(2^{-1/6}z^{1/2}\right)^{n^{2}}}{n!}{n\choose r}$
$\displaystyle\hphantom{D_{n,r}(z)=}{}\times\exp\left(\frac{\sqrt{2}{\rm
i}}{3}(n-2r)z^{3/2}\right)\int_{\Gamma^{n}}\Delta_{n}(\mathbf{u})^{2}\exp\left(\frac{z^{3/2}}{\sqrt{2}}\left[\phi(\mathbf{u})-\phi(\mathbf{u}^{(r)})\right]\right)\prod_{k=1}^{n}{\rm
d}u_{k}.$

As before, we isolate these points by fixing $\delta>0$ and two discs
$D\big{(}\mathbf{u}^{(r)},\delta\big{)}$ of radius $\delta$ around the
stationary points in each of the variables $u_{k}$. We take
$\Gamma_{\delta}=\Gamma\cap D\big{(}\mathbf{u}^{r},\delta\big{)}$, where
$\Gamma$ is the corresponding path of steepest descent. Then as $z\to\infty$
we have

$\displaystyle D_{n,r}(z)=\exp\left(\frac{\sqrt{2}{\rm
i}}{3}(n-2r)z^{3/2}\right)[I_{n,r}(z)+E_{n,r}(z)],$

where

$\displaystyle
I_{n,r}(z)=\frac{\left(2^{-1/6}z^{1/2}\right)^{n^{2}}}{r!(n-r)!}\int_{\Gamma_{\delta}}\Delta_{n}(\mathbf{u})^{2}\exp\left(\frac{z^{3/2}}{\sqrt{2}}\big{[}\phi(\mathbf{u})-\phi\big{(}\mathbf{u}^{(r)}\big{)}\big{]}\right)\prod_{k=1}^{n}{\rm
d}u_{k},$

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 $\operatorname{Ai}$ and $\operatorname{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 $\operatorname{Ai}$ and $u_{-}$ for $\operatorname{Bi}$). In the
oscillatory regime, however, we have (independently of the parameter $r$), $p$
integrals of type $\operatorname{Ai}$ and $n-p$ integrals of type
$\operatorname{Bi}$, with $0\leq p\leq n$, and any integral around the
stationary point $-{\rm i}$ has different orientation depending if we are
integrating along the $\operatorname{Ai}$ or the $\operatorname{Bi}$ contour:
the contour for an integral along an $\operatorname{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 $\operatorname{Ai}$ type in the last $n-r$
cases (where the point $-{\rm i}$ is taken into account) and therefore $p-q$
integrals of $\operatorname{Ai}$ type in the first $r$ cases, for any value of
$q$ from $0$ to $p$. We define

$\displaystyle H_{n,r,p}=\sum_{q=\max(0,p-r)}^{\min(p,n-r)}{r\choose
p-q}{n-r\choose q}(-1)^{q},$

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

$\displaystyle I_{n,r}(z)=(-4)^{r(n-r)}\exp\left[(n-2r){\rm
i}\left(\frac{\sqrt{2}z^{3/2}}{3}+\frac{n\pi}{4}\right)\right]\frac{S_{r}(z)S_{n-r}(z)}{r!(n-r)!}\sum_{p=0}^{n}H_{n,r,p}$
$\displaystyle\hphantom{I_{n,r}(z)=}{}\times\big{(}2^{-1/6}z^{1/2}\big{)}^{n^{2}}\big{(}1+\mathcal{O}\big{(}z^{-3/2}\big{)}\big{)},$

in terms of Selberg integrals (3.8) again. This leads to

$\displaystyle
I_{n,r}(z)=2^{-\frac{5n^{2}}{12}+\frac{n}{2}}\pi^{\frac{n}{2}}\,z^{-\frac{n^{2}}{4}-\frac{3}{2}r(n-r)}\exp\left[(n-2r){\rm
i}\left(\frac{\sqrt{2}z^{3/2}}{3}+\frac{n\pi}{4}\right)\right]M_{n,r}\sum_{p=0}^{n}H_{n,r,p}$
$\displaystyle\hphantom{I_{n,r}(z)=}{}\times\big{(}1+\mathcal{O}\big{(}z^{-3/2}\big{)}\big{)},$

where the coefficient $M_{n,r}$ is given by (2.5). The powers of $2$ and $\pi$
can then be combined and simplified using (3.3), which relates $\tau_{n}(z)$
and $D_{n}(z)$, and (4.1).

Finally, the asymptotic expansion can be written in terms of sines and
cosines, instead of complex exponentials, noting that
$H_{n,n-r,p}=(-1)^{p}H_{n,r,p}$, for $0\leq r\leq n$, and grouping terms
depending on the parity of $p$. The proof of this symmetry relation can be
obtained by writing the different cases: if $p\leq n-r$ and $p-r\leq 0$, then

$\displaystyle H_{n,r,p}=\sum_{q=0}^{p}{r\choose p-q}{n-r\choose q}(-1)^{q},$

and

$\displaystyle H_{n,n-r,p}=\sum_{q=0}^{p}{n-r\choose p-q}{r\choose
q}(-1)^{q}=\sum_{q=0}^{p}{n-r\choose q}{r\choose
p-q}(-1)^{p-q}=(-1)^{p}H_{n,r,p},$

where we have reversed the order inside the sum (i.e., $q\mapsto p-q$). If
$p\leq n-r$ and $p-r\geq 0$, then

$\displaystyle H_{n,r,p}=\sum_{q=p-r}^{p}{r\choose p-q}{n-r\choose
q}(-1)^{q},$

and

$\displaystyle H_{n,n-r,p}=\sum_{q=0}^{r}{n-r\choose p-q}{r\choose
q}(-1)^{q}=\sum_{q=0}^{r}{n-r\choose p-r+q}{r\choose r-q}(-1)^{r-q}$
$\displaystyle\hphantom{H_{n,n-r,p}}{}=\sum_{s=p-r}^{p}{n-r\choose s}{r\choose
p-s}(-1)^{p-s}=(-1)^{p}H_{n,r,p},$

reversing the sum in the first step and shifting $q=s-p+r$ in the second one.

The case $p>n-r$ can be proved in a similar way.

##  5 Proof of Corollaries 2.4 and 2.7

In the non-oscillatory regime, we can derive asymptotic expansions for the
Painlevé functions $\sigma_{n}(z)$, $p_{n}(z)$ and $q_{n}(z)$ in (1.4) quite
straightforwardly. We will use the asymptotic expansions for $\tau_{n}(z)$ and
(1.4), observing that differentiation is permitted since $\tau_{n}(z)$ is an
analytic function of $z$ inside the relevant sectors, see [29, Section 1.8].
In the oscillatory regime, one would need to open a sector around the positive
real axis; this calculation is similar to the oscillatory case, writing
$z=\rho{\rm e}^{{\rm i}\alpha}$, $|\alpha|<\frac{\pi}{3}$ and considering both
stationary points.

Suppose first that $C_{2}\neq 0$. Instead of working directly with the
asymptotic expansion (2.2), it is simpler to pick the leading term therein,
corresponding to $r=0$, i.e.

$\displaystyle\tau_{n}(z)=K_{n}\mathbf{A}_{n,0}(z){\rm
e}^{\frac{\sqrt{2}}{3}n(-z)^{\frac{3}{2}}}+\mathcal{O}\Big{(}{\rm
e}^{\frac{\sqrt{2}}{3}(n-2)(-z)^{\frac{3}{2}}}\Big{)}$

as $|z|\to\infty$ in the sector $|\arg(-z)|<\frac{\pi}{3}$. Using this result,
we can deduce the form of the asymptotic expansion for the function
$\sigma_{n}(z)$, which is

$\displaystyle\sigma_{n}(z)=\frac{\tau_{n}^{\prime}(z)}{\tau_{n}(z)}=\sum_{k=0}^{K-1}s_{n,k}(-z)^{\tfrac{1}{2}-\tfrac{3k}{2}}+\mathcal{O}\Big{(}(-z)^{\tfrac{1}{2}-\tfrac{3K}{2}}\Big{)},\qquad
K\geq 1.$ (5.1)

Then, substituting this expansion into the differential equation (1.6), we can
identify the coefficients $s_{n,k}$. Then, we can compute

$\displaystyle
p_{n}(z)=-2\sigma^{\prime}_{n}(z)=-2\sum_{k=0}^{K-1}p_{n,k}(-z)^{-\tfrac{1}{2}-\tfrac{3k}{2}}+\mathcal{O}\left((-z)^{\tfrac{1}{2}-\tfrac{3K}{2}}\right),$

where the coefficients $p_{n,k}$ follow easily from $s_{n,k}$. Finally, from
(1.4) again, we have

$\displaystyle
q_{n}(z)=\sigma_{n-1}(z)-\sigma_{n}(z)=\sum_{k=0}^{K-1}q_{n,k}(-z)^{-\tfrac{1}{2}-\tfrac{3k}{2}}+\mathcal{O}\Big{(}(-z)^{\tfrac{1}{2}-\tfrac{3K}{2}}\Big{)},$

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

$\displaystyle\tau_{n}(z)=K_{n}\mathbf{A}_{n,n}(z){\rm
e}^{-\frac{\sqrt{2}}{3}n(-z)^{\frac{3}{2}}}.$

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 (5.1). The reason for this is that in Theorem 2.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
symbolic software, but for brevity we will give only the leading terms.

The leading terms in Theorem 2.5 clearly correspond to $r=s$ if $n=2s$ is
even, and to $r=s-1$ if $n=2s-1$ is odd. For $s\geq 1$, we have

$\displaystyle\sigma_{2s}(z)=\frac{\tau_{2s}^{\prime}(z)}{\tau_{2s}(z)}=\frac{s^{2}}{2z}+\frac{s}{2z}\,\frac{b_{2s,s-1}^{(0)}\sin(\psi_{2s,s-1}(z))-d_{2s,s-1}^{(0)}\cos(\psi_{2s,s-1}(z))}{b_{2s,s}^{(0)}}+\mathcal{O}\big{(}z^{-5/2}\big{)},$

with coefficients given by (2.6) and phase function (2.4).

If $C_{2}=0$, this expression simplifies considerably, since only the term
$p=s$ survives, and then $b_{2s,s}^{(0)}=-b_{2s,s-1}^{(0)}=C_{1}^{2s}$ and
$d_{2s,s-1}^{(0)}=0$, so

$\displaystyle\sigma_{2s}(z)=\frac{s^{2}}{2z}-\frac{s}{2z}\sin(\psi_{2s,s-1}(z))+\mathcal{O}\big{(}z^{-5/2}\big{)}.$

This is in agreement with [13, Theorem 9], correcting the phase function. In
the odd case, we have

$\displaystyle\sigma_{2s-1}(z)=\sqrt{\frac{z}{2}}\frac{-b_{2s-1,s-1}^{(0)}\sin(\psi_{2s-1,s-1}(z))+d_{2s-1,s-1}^{(0)}\cos(\psi_{2s-1,s-1}(z))}{b_{2s-1,s-1}^{(0)}\cos(\psi_{2s-1,s-1}(z))+d_{2s-1,s-1}^{(0)}\sin(\psi_{2s-1,s-1}(z))}+\mathcal{O}\big{(}z^{-1}\big{)}.$

If $C_{2}=0$, this expression simplifies again

$\displaystyle\sigma_{2s-1}(z)=\sqrt{\frac{z}{2}}\cot(\psi_{2s-1,s-1}(z))+\mathcal{O}\big{(}z^{-1}\big{)}.$

In both cases, the asymptotic approximation is valid away from the zeros of
the denominators that appear in the leading terms, we refer the reader to [29,
Chapter 1, Section 8.1] for the general theory.

Similar calculations lead to the asymptotic expansions for $p_{n}(z)$,
although care is needed because the leading terms may come from subleading
ones before, as a result of differentiation of the trigonometric functions.
Straightforward manipulations, using the fact that
$p_{n}(z)=-2\sigma_{n}^{\prime}(z)$ and
$q_{n}(z)=\sigma_{n-1}(z)-\sigma_{n}(z)$, for $n\geq 1$, which follows from
(1.4), lead to (2.7), (2.8), (2.9), and (2.10).

Higher order terms in these asymptotic expansions can be computed by using the
corresponding differential equations. This is quite straightforward in the
non-oscillatory regime, but more involved in the oscillatory one, since
several trigonometric functions with different phase functions intervene.

##  Appendix A Stokes phenomenon for the $\boldsymbol{\operatorname{Ai}}$
solution

The case $C_{2}=0$, where the seed function only contains the Airy
$\operatorname{Ai}$ functions, is especially relevant both because of its
asymptotic behavior and in applications. It is interesting to observe that in
this case we have strong asymptotics for $\tau_{n}(z)$ in the whole cut plane
$|\arg(-z)|<\pi$ involving only one exponential factor, see (2.3); however, we
can also use the rotational symmetry given by Lemma 2.1, to study the
asymptotics in the regions $0<\arg z<\frac{2\pi}{3}$ and $-\frac{2\pi}{3}<\arg
z<0$. This calculation uses a different seed function in the original sector
$|\arg(-z)|<\frac{\pi}{3}$ and it gives subdominant exponential terms. These
do not affect the leading asymptotic behavior, but they give a non-linear
Stokes phenomenon for this family of solutions of ${\rm P}_{\rm II}$, in the
spirit of Its and Kapaev [24], see also [17, Chapter 11].

For instance, if we consider the sector $0<\arg z<\frac{2\pi}{3}$ and the seed
function with $\widetilde{C}_{1}=1$ and $\widetilde{C}_{2}=0$ (corresponding
to the pure Airy $\operatorname{Ai}$ function in this rotated sector), we have
from (2.3) the following asymptotic approximation

$\displaystyle\tau_{n}(z)=K_{n}(-z)^{-\frac{n^{2}}{4}}\mathbf{A}_{n,n}(z){\rm
e}^{-\frac{\sqrt{2}}{3}n(-z)^{\frac{3}{2}}}$
$\displaystyle\hphantom{\tau_{n}(z)}{}=(-1)^{\left\lfloor{n/2}\right\rfloor}K_{n}G(n+1)(-z)^{-\frac{n^{2}}{4}}\big{(}1+\mathcal{O}\big{(}(-z)^{-\frac{3}{2}}\big{)}\big{)}{\rm
e}^{-\frac{\sqrt{2}}{3}n(-z)^{\frac{3}{2}}}.$

However, applying Lemma 2.1, we can use $\tau_{n}(z)={\rm e}^{\frac{2\pi{\rm
i}}{3}n(n-1)}\tau_{n}\big{(}z{\rm e}^{\frac{2\pi{\rm i}}{3}}\big{)}$, and then
the seed function in the region $|\arg(-z)|<\frac{\pi}{3}$ is
$\varphi\big{(}z{\rm e}^{\frac{2\pi{\rm i}}{3}}\big{)}$, with constants

$\displaystyle C_{1}=\tfrac{1}{2}{\rm e}^{-\frac{\pi{\rm i}}{3}},\qquad
C_{2}=\tfrac{1}{2}{\rm e}^{\frac{\pi{\rm i}}{6}}.$

After simplification, the first two coefficients $\mathbf{A}_{n,0}\big{(}z{\rm
e}^{\frac{2\pi{\rm i}}{3}}\big{)}$ and $\mathbf{A}_{n,1}\big{(}z{\rm
e}^{\frac{2\pi{\rm i}}{3}}\big{)}$ give

$\displaystyle\tau_{n}(z)=(-1)^{\left\lfloor{n/2}\right\rfloor}K_{n}G(n+1)(-z)^{-\frac{n^{2}}{4}}{\rm
e}^{-\frac{\sqrt{2}}{3}n(-z)^{\frac{3}{2}}}$
$\displaystyle\hphantom{\tau_{n}(z)=}{}\times\left[1+(-1)^{n}\frac{2^{\frac{5n-7}{2}}{\rm
i}(-z)^{\frac{3}{2}(n-1)}}{\Gamma(n)}{\rm
e}^{\frac{2\sqrt{2}}{3}(-z)^{\frac{3}{2}}}\right]\big{(}1+\mathcal{O}\big{(}(-z)^{-\frac{3}{2}}\big{)}\big{)}.$

Other exponentially small contributions can be calculated in a similar way. As
a consequence, we have a similar Stokes phenomenon for the Painlevé functions:

$\displaystyle\sigma_{n}(z)=-\frac{n(-z)^{1/2}}{\sqrt{2}}\left[1+(-1)^{n+1}\frac{2^{\frac{5n-5}{2}}{\rm
i}}{\Gamma(n+1)}(-z)^{\frac{3}{2}(n-1)}{\rm
e}^{\frac{2\sqrt{2}}{3}(-z)^{\frac{3}{2}}}\right]\big{(}1+\mathcal{O}\big{(}(-z)^{-\frac{3}{2}}\big{)}\big{)},$
$\displaystyle
p_{n}(z)=-\frac{n}{\sqrt{2}(-z)^{1/2}}\left[1+(-1)^{n+1}\frac{2^{\frac{5n}{2}-2}{\rm
i}}{\Gamma(n)}(-z)^{\frac{3}{2}n}{\rm
e}^{\frac{2\sqrt{2}}{3}(-z)^{\frac{3}{2}}}\right]\big{(}1+\mathcal{O}\big{(}(-z)^{-\frac{3}{2}}\big{)}\big{)},$

and

$\displaystyle
q_{n}(z)=\frac{(-z)^{1/2}}{\sqrt{2}}\left[1+(-1)^{n+1}\frac{2^{\frac{5n-5}{2}}{\rm
i}}{\Gamma(n)}(-z)^{\frac{3}{2}(n-1)}{\rm
e}^{\frac{2\sqrt{2}}{3}(-z)^{\frac{3}{2}}}\right]\big{(}1+\mathcal{O}\big{(}(-z)^{-\frac{3}{2}}\big{)}\big{)},$

for $n\geq 1$.

Similar calculations can be carried out in the sector $\frac{4\pi}{3}<\arg
z<2\pi$, using the seed function $\varphi\big{(}z{\rm e}^{-\frac{2\pi{\rm
i}}{3}}\big{)}$ with constants

$\displaystyle C_{1}=\tfrac{1}{2}{\rm e}^{\frac{\pi{\rm i}}{3}},\qquad
C_{2}=\tfrac{1}{2}{\rm e}^{-\frac{\pi{\rm i}}{6}}.$

As mentioned before, these extra exponential terms are not strictly needed in
the non-oscillatory regime $|\arg(-z)|<\pi$, but they become relevant when one
examines the transition between exponential and trigonometric behavior.

An alternative way to calculate these subleading exponential terms would be to
extend the asymptotic expansion obtained in the oscillatory regime, with
trigonometric functions, to a wider sector in $\mathbb{C}$ around the positive
real axis. Then in the overlapping region between non-oscillatory and
oscillatory behavior the two asymptotic expansions should coincide, and it
would be possible to extract the subleading terms.

### Acknowledgements

The author acknowledges financial support from the EPSRC grant “Painlevé
equations: analytical properties and numerical computation”, reference
EP/P026532/1, and from the project MTM2015-65888-C4-2-P from the Spanish
Ministry of Economy and Competitivity. The author wishes to thank M.
Fasondini, D. Huybrechs, A. Iserles, A.R. Its, A.B.J. Kuijlaars, A.F.
Loureiro, C. Pech and W. Van Assche for stimulating discussions on the topic
and scope of this paper, as well as the organisers of the workshop “Painlevé
Equations and Applications” held at the University of Michigan, August 25–29,
2017, for their hospitality. The comments, remarks and corrections of the
anonymous referees have lead to an improved version of the paper, and they are
greatly appreciated.

## References

  * [1] Balogh F., Bertola M., Bothner T., Hankel determinant approach to generalized Vorob’ev–Yablonski polynomials and their roots, Constr. Approx. 44 (2016), 417–453, arXiv:1504.00440.
  * [2] Bertola M., Bothner T., Zeros of large degree Vorob’ev–Yablonski polynomials via a Hankel determinant identity, Int. Math. Res. Not. 2015 (2015), 9330–9399, arXiv:1401.1408.
  * [3] Bleher P., Deaño A., Topological expansion in the cubic random matrix model, Int. Math. Res. Not. 2013 (2013), 2699–2755, arXiv:1011.6338.
  * [4] Bleher P., Deaño A., Painlevé I double scaling limit in the cubic random matrix model, Random Matrices Theory Appl. 5 (2016), 1650004, 58 pages, arXiv:1310.3768.
  * [5] Bleher P., Deaño A., Yattselev M., Topological expansion in the complex cubic log-gas model: one-cut case, J. Stat. Phys. 166 (2017), 784–827, arXiv:1606.04303.
  * [6] Bleistein N., Saddle point contribution for an $n$-fold complex-valued integral, 2012.
  * [7] Bleistein N., Handelsman R.A., Asymptotic expansions of integrals, 2nd ed., Dover Publications, Inc., New York, 1986.
  * [8] Bothner T., Miller P.D., Sheng Y., Rational solutions of the Painlevé-III equation, arXiv:1801.04360.
  * [9] Buckingham R.J., Large-degree asymptotics of rational Painlevé-IV functions associated to generalized Hermite polynomials, Int. Math. Res. Not., to appear, arXiv:1706.09005.
  * [10] Buckingham R.J., Miller P.D., Large-degree asymptotics of rational Painlevé-II functions: noncritical behaviour, Nonlinearity 27 (2014), 2489–2578, arXiv:1310.2276.
  * [11] Buckingham R.J., Miller P.D., Large-degree asymptotics of rational Painlevé-II functions: critical behaviour, Nonlinearity 28 (2015), 1539–1596, arXiv:1406.0826.
  * [12] Clarkson P.A., Painlevé equations – nonlinear special functions, in Orthogonal polynomials and special functions, Lecture Notes in Math., Vol. 1883, Springer, Berlin, 2006, 331–411.
  * [13] Clarkson P.A., On Airy solutions of the second Painlevé equation, Stud. Appl. Math. 137 (2016), 93–109, arXiv:1510.08326.
  * [14] Clarkson P.A., Jordaan K., The relationship between semiclassical Laguerre polynomials and the fourth Painlevé equation, Constr. Approx. 39 (2014), 223–254, arXiv:1301.4134.
  * [15] Clarkson P.A., Loureiro A.F., Van Assche W., Unique positive solution for an alternative discrete Painlevé I equation, J. Difference Equ. Appl. 22 (2016), 656–675, arXiv:1508.04916.
  * [16] Fedoryuk M.V., Asymptotic methods in analysis, in Analysis I, Encyclopaedia of Mathematical Sciences, Vol. 13, Springer, Berlin, Heidelberg, 1989, 83–191.
  * [17] Fokas A.S., Its A.R., Kapaev A.A., Novokshenov V.Yu., Painlevé transcendents. The Riemann–Hilbert approach, Mathematical Surveys and Monographs, Vol. 128, Amer. Math. Soc., Providence, RI, 2006.
  * [18] Forrester P.J., Log-gases and random matrices, London Mathematical Society Monographs Series, Vol. 34, Princeton University Press, Princeton, NJ, 2010.
  * [19] Forrester P.J., Witte N.S., Application of the $\tau$-function theory of Painlevé equations to random matrices: PIV, PII and the GUE, Comm. Math. Phys. 219 (2001), 357–398, math-ph/0103025.
  * [20] Forrester P.J., Witte N.S., Application of the $\tau$-function theory of Painlevé equations to random matrices: $\rm P_{V}$, $\rm P_{III}$, the LUE, JUE, and CUE, Comm. Pure Appl. Math. 55 (2002), 679–727, math-ph/0201051.
  * [21] Gromak V.I., Laine I., Shimomura S., Painlevé differential equations in the complex plane, De Gruyter Studies in Mathematics, Vol. 28, Walter de Gruyter & Co., Berlin, 2002.
  * [22] Ismail M.E.H., Classical and quantum orthogonal polynomials in one variable, Encyclopedia of Mathematics and its Applications, Vol. 98, Cambridge University Press, Cambridge, 2009.
  * [23] Its A.R., Kapaev A.A., The irreducibility of the second Painlevé equation and the isomonodromy method, in Toward the Exact WKB Analysis of Differential Equations, Linear or Non-Linear (Kyoto, 1998), Kyoto University Press, Kyoto, 2000, 209–222.
  * [24] Its A.R., Kapaev A.A., Quasi-linear Stokes phenomenon for the second Painlevé transcendent, Nonlinearity 16 (2003), 363–386, nlin.SI/0108010.
  * [25] Its A.R., Kuijlaars A.B.J., Östensson J., Critical edge behavior in unitary random matrix ensembles and the thirty-fourth Painlevé transcendent, Int. Math. Res. Not. 2008 (2008), Art. ID rnn017, 67 pages, arXiv:0704.1972.
  * [26] Its A.R., Kuijlaars A.B.J., Östensson J., Asymptotics for a special solution of the thirty fourth Painlevé equation, Nonlinearity 22 (2009), 1523–1558, arXiv:0811.3847.
  * [27] Miller P.D., Applied asymptotic analysis, Graduate Studies in Mathematics, Vol. 75, Amer. Math. Soc., Providence, RI, 2006.
  * [28] Okamoto K., Studies on the Painlevé equations. III. Second and fourth Painlevé equations, $P_{{\rm II}}$ and $P_{{\rm IV}}$, Math. Ann. 275 (1986), 221–255.
  * [29] Olver F.W.J., Asymptotics and special functions, AKP Classics, A K Peters, Ltd., Wellesley, MA, 1997.
  * [30] Olver F.W.J., Olde Daalhuis A.B., Lozier D.W., Schneider B.I., Boisvert R.F., Clark C.W., Miller B.R., Saunders B.V. (Editors), NIST digital library of mathematical functions, Release 1.0.16 of 2017-09-18, available at http://dlmf.nist.gov/.
  * [31] Szegő G., Orthogonal polynomials, Colloquium Publications, Vol. 23, 4th ed., Amer. Math. Soc., Providence, R.I., 1975.
  * [32] Temme N.M., Asymptotic methods for integrals, Series in Analysis, Vol. 6, World Sci. Publ. Co. Pte. Ltd., Hackensack, NJ, 2015.
  * [33] Van Assche W., Discrete Painlevé equations for recurrence coefficients of orthogonal polynomials, in Difference Equations, Special Functions and Orthogonal Polynomials, World Sci. Publ., Hackensack, NJ, 2007, 687–725, math.CA/0512358.
  * [34] Van Assche W., Orthogonal polynomials and Painlevé equations, Australian Mathematical Society Lecture Series, Vol. 27, Cambridge University Press, Cambridge, 2018.
  * [35] Van Assche W., Filipuk G., Zhang L., Multiple orthogonal polynomials associated with an exponential cubic weight, J. Approx. Theory 190 (2015), 1–25, arXiv:1306.3835.