* [203] Zumbach, G. Volatility conditional on price trends.  Quantitative finance 10, 4 (2010), 431–442.
  * [204] Zumbach, G. Discrete time series, processes, and applications in finance, 2013 ed.  Springer, Berlin, Germany, 2012.


Painlevé IV Critical Asymptotics for Orthogonal Polynomials in the Complex
Plane

Painlevé IV Critical Asymptotics
for Orthogonal Polynomials 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

Marco BERTOLA †‡, José Gustavo ELIAS REBELO † and Tamara GRAVA †§

M. Bertola, J.G. Elias Rebelo and T. Grava

† Area of Mathematics, SISSA, via Bonomea 265 - 34136, Trieste, Italy
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>‡ Department of
Mathematics and Statistics, Concordia University,
‡ 1455 de Maisonneuve W., Montréal, Québec, Canada H3G 1M8 § School of
Mathematics, University of Bristol, UK

Received February 06, 2018, in final form August 14, 2018; Published online
August 30, 2018

We study the asymptotic behaviour of orthogonal polynomials in the complex
plane that are associated to a certain normal matrix model. The model depends
on a parameter and the asymptotic distribution of the eigenvalues undergoes a
transition for a special value of the parameter, where it develops a corner-
type singularity. In the double scaling limit near the transition we determine
the asymptotic behaviour of the orthogonal polynomials in terms of a solution
of the Painlevé IV equation. We determine the Fredholm determinant associated
to such solution and we compute it numerically on the real line, showing also
that the corresponding Painlevé transcendent is pole-free on a semiaxis.

orthogonal polynomials on the complex plane; Riemann–Hilbert problem; Painlevé
equations Fredholm determinant

34M55; 34M56; 33C15

To the memory of Andrei Kapaev, a master of Painlevé equation, a very generous
colleague and a special friend.

##  1 Introduction

In this work we consider the critical asymptotic behaviour for the orthogonal
polynomials that are associated to the matrix model described here. Over the
set of normal $n\times n$ matrices
$\mathcal{N}_{n}:=\\{M\colon[M,M^{\star}]=0\\}\subset\operatorname{Mat}_{n\times
n}(\mathbb{C})$ we consider the measure of the form

$\displaystyle M\to\frac{1}{\mathcal{Z}_{n,N}}{\rm
e}^{-N\operatorname{Tr}(W(M))}\mathrm{d}M,\qquad\mathcal{Z}_{n,N}=\int_{\mathcal{N}_{n}}{\rm
e}^{-N\operatorname{Tr}(W(M))}\mathrm{d}M,$ (1.1)

where ${\rm d}M$ stands for the induced volume form on $\mathcal{N}_{n}$,
invariant under conjugation by unitary matrices.

Since normal matrices are diagonalizable by unitary transformations, the
probability density (1.1) can be reduced to the form [35]

$\displaystyle\dfrac{1}{Z_{n,N}}\prod_{i<j}|\lambda_{i}-\lambda_{j}|^{2}{\rm
e}^{-N\sum\limits_{j=1}^{n}{W(\lambda_{i})}}\mathrm{d}A(\lambda_{1})\cdots\mathrm{d}A(\lambda_{n}),$
(1.2)

where $\lambda_{j}$ are the complex eigenvalues of the normal matrix $M$ [11]
and ${\rm d}A(\lambda)={\rm d}\operatorname{Re}(\lambda){\rm
d}\operatorname{Im}(\lambda)$ is the area element in the complex plane.

Associated to the above data we consider the sequence of (monic) orthogonal
polynomials $p_{n}(z)=z^{n}+\cdots$,

$\displaystyle\int_{\mathbb{C}}p_{n}(\lambda)\overline{p_{m}(\lambda)}{\rm
e}^{-NW(\lambda)}\mathrm{d}A(\lambda)=h_{n,N}\delta_{n,m},\qquad
h_{n,N}>0,\qquad n,m=0,1,\dots,$ (1.3)

where $W\colon\mathbb{C}\rightarrow\mathbb{R}$ is the external potential.

Likewise in the standard GUE, the orthogonal polynomial $p_{n}(z)$ is also the
expectation value of the characteristic polynomial of the random normal matrix
$M$: $p_{n}(z)=\langle\det(z{\mathbf{1}}-M)\rangle$, where the expectation is
with respect to the measure (1.1) (see, e.g., [18]). The statistical
quantities related to the eigenvalues of matrix models can be expressed in
terms of the associated orthogonal polynomials, $p_{n}(\lambda)$. In
particular, the probability measure (1.2) can be written as a determinantal
point field (sometimes also called “process”) with Christoffel–Darboux kernel

$\displaystyle K(\lambda,\eta)={\rm
e}^{-N(W(\lambda)+W(\eta))/2}\sum_{j=0}^{n-1}\frac{1}{h_{j,N}}p_{j}(\lambda)\overline{p_{j}(\eta)}.$

The average density of eigenvalues is

$\displaystyle\rho_{n,N}(\lambda)=\frac{1}{n}K(\lambda,\lambda)$

and in the limit

$\displaystyle n\to\infty,\qquad N\to\infty,\qquad\frac{N}{n}\to\frac{1}{T},$

it converges to the unique probability measure, $\mu^{*}$, in the plane which
minimizes the functional [19, 26]

$\displaystyle
I(\mu)=\iint\log|\lambda-\eta|^{-1}\mathrm{d}\mu(\lambda)\mathrm{d}\mu(\eta)+\frac{1}{T}\int
W(\lambda)\mathrm{d}\mu(\lambda).$

The functional $I(\mu)$ is the Coulomb energy functional in two dimensions and
the existence of a unique minimizer is a well-established fact under mild
assumptions on the potential $W(\lambda)$ [37]. If $W(\lambda)$ is twice
continuously differentiable then the equilibrium measure is given by

$\displaystyle\mathrm{d}\mu^{*}(\lambda)=\frac{1}{4\pi T}\Delta
W(\lambda)\chi_{D}(\lambda)\mathrm{d}A(\lambda),$

where $\chi_{D}$ is the characteristic function of the compact support set
$D=\operatorname{supp}(\mu^{*})$ and
$\Delta=4\partial_{\lambda}\bar{\partial}_{\lambda}$. For example, if
$W(\lambda)=|\lambda|^{2}$, one has the Ginibre ensemble [23] and the measure
$\mathrm{d}\mu^{*}(\lambda)$ is the uniform measure on the disk of radius
$\sqrt{T}$. Following [4], we consider the external potential to be of the
form

$\displaystyle
W(\lambda)=|\lambda|^{2d}-t\lambda^{d}-\bar{t}\bar{\lambda}^{d},\qquad\lambda\in\mathbb{C},\qquad
d\in{\mathbb{N}},\qquad t\in{\mathbb{C}}^{*}.$ (1.4)

The case $d=1$ is reducible to Hermite polynomials and we will consider only
$d\geq 2$. Up to a rotation $\lambda\mapsto{\rm e}^{{\rm i}\theta}\lambda$ we
can assume without loss of generality that $t\in{\mathbb{R}}_{+}$, leading to
the external potential and orthogonality measure

$\displaystyle
W(\lambda)=|\lambda|^{2d}-t\big{(}\lambda^{d}+\bar{\lambda}^{d}\big{)},\qquad\lambda\in\mathbb{C},\qquad{\rm
e}^{-N(|\lambda|^{2d}-t(\lambda^{d}+\bar{\lambda}^{d}))}\mathrm{d}A(\lambda).$
(1.5)

There remains a residual discrete rotational $\mathbb{Z}_{d}$-symmetry which
allows us to further reduce the problem as observed in [5] and [20]. The
potential $W$ in (1.5) can be written as
$W(\lambda)=\frac{1}{d}Q\big{(}\lambda^{d}\big{)}$, and hence the equilibrium
measure for $W$ can be obtained from the equilibrium measure associated to $Q$
by an unfolding procedure. Considering the particular case of (1.4), the
potential $Q$ can be rewritten as

$\displaystyle Q(u)=d|u|^{2}-dt(u+\bar{u})=d|u-t|^{2}-t^{2}d,$

which corresponds to the Ginibre ensemble whose equilibrium measure is the
normalized area measure of the disk centered at $u=t$ and of radius $t_{c}$

$\displaystyle|u-t|\leq t_{c},\qquad t_{c}=\sqrt{\frac{T}{d}}.$

Pulling back to the $\lambda$-plane, we obtain the equilibrium measure of $W$

$\displaystyle\mathrm{d}\mu_{W}=\frac{d}{\pi
t^{2}_{c}}|\lambda|^{2(d-1)}\chi_{D}\mathrm{d}A(\lambda),$

where $\mathrm{d}A$ is the area measure and $\chi_{D}$ is the characteristic
function of the domain $D$ whose boundary is a lemniscate

$\displaystyle
D:=\big{\\{}\lambda\in\mathbb{C},\,\big{|}\lambda^{d}-t\big{|}\leq
t_{c}\big{\\}}.$ (1.6)

Figure 1: Domain $D$ for different values of $t$: on the left $t<t_{c}$, in
the center $t=t_{c}$ and on the right $t>t_{c}$.

A plot of the domain $D$ for different values of $t$ is shown in Fig. 1.
Depending on the radius $t_{c}$, the domain $D$ is either connected for $t\leq
t_{c}$ or has $d$ connected components for $t>t_{c}$. The cases $t>t_{c}$ and
$t<t_{c}$ were analysed in [4]. In this manuscript we study the critical
regime $t=t_{c}$. Boundary points where $\Delta W>0$ have been classified by
Sakai [38] and fall in the category of regular points, cusps or double points.
Details about the application of Sakai’s theory to the setting when $\Delta
W>0$ near a singular point can be found in several papers, e.g., [33]. The
classification of boundary points when $\Delta W=0$ is still an open problem
to the best of our knowledge. In this manuscript we take a first step in
considering a specific example when $\Delta W=0$. In particular we will study
the asymptotic behaviour of the orthogonal polynomials associated to this
model.

The asymptotic analysis of orthogonal polynomials when the measure is of
planar type, is much harder than the orthogonality on contours. When the
orthogonality planar measure is not varying there are several studies like in
[36] or [24]. In a very recent work [27] the asymptotic of orthogonal
polynomials for a varying planar measure was obtained in the regular case.
Critical regimes have been considered in [34], in [9, 32] studying a normal
matrix model with cut-off, and in [1] where new types of determinantal point
fields, have emerged. The equilibrium problem for a class of potentials to
which our case belongs, was also recently considered in [2].

In this paper we derive the asymptotic behaviour for large $n$ of the
orthogonal polynomials associated to the exponential weight (1.5) in the
critical regime $t_{c}$. The zeros of the polynomials accumulate on the Szegö
curve that was first observed in relation to the zeros of the Taylor
polynomials of the exponential function [40]. For large $n$ we obtain the
expansion of the orthogonal polynomials in terms of a special solution of a
Painlevé IV transcendent. We show that such solution is pole free on the
negative semi-axis. We obtain the Fredholm determinant associated to such
solution and compute it numerically on the real line. Painlevé IV critical
asymptotic was observed for orthogonal polynomials with deformed Laguerre
weights [15]. Finally we remark that normal matrix models are strictly related
to the 2-dimensional Toda lattice, see, e.g., [41, 43]. It would be
interesting to study the critical case considered in this paper in this
perspective as done in [6] or [12, 13, 15] for the one-dimensional Toda
lattice.

###  1.1 Zeros of orthogonal polynomials: statement of the result

The zeros of the orthogonal polynomials (1.3) in the critical case $t=t_{c}$
concentrate on the curve $\hat{\mathcal{C}}$ that is described below. Let us
introduce the function

$\displaystyle\hat{\varphi}(\lambda)=\log\big{(}t_{c}-\lambda^{d}\big{)}+\dfrac{\lambda^{d}}{t_{c}}-\log
t_{c},$ (1.7)

and consider the level curve $\hat{\mathcal{C}}$

$\displaystyle\hat{\mathcal{C}}:=\left\\{\lambda\in\mathbb{C},\,\operatorname{Re}\hat{\varphi}(\lambda)=\log\dfrac{\big{|}t_{c}-\lambda^{d}\big{|}}{t_{c}}+\dfrac{\operatorname{Re}\lambda^{d}}{t_{c}}=0,\,\big{|}\lambda^{d}-t_{c}\big{|}\leq
t_{c}\right\\}.$ (1.8)

Observe that the level curve $\hat{\mathcal{C}}$ consists of a closed contour
contained in the set $D$ defined in (1.6) with $t=t_{c}$ because it is the
pullback under the map $z=1-\frac{\lambda^{d}}{t_{c}}$ of the celebrated
(closed) Szegö curve $\big{|}z{\rm e}^{1-z}\big{|}=1$ [40] (see Fig. 2).

Figure 2: The boundary of the domain $D$ that describes the support of the
eigenvalues distribution in black and the level curve $\hat{\mathcal{C}}$ in
red at the critical time $t_{c}$ and for $d=5$.

Let us define the measure $\hat{\nu}$, associated with this family of curves,
given by

$\displaystyle\mathrm{d}\hat{\nu}(\lambda)=\dfrac{1}{2\pi{\rm
i}d}\mathrm{d}\hat{\varphi}(\lambda),$ (1.9)

and supported on $\hat{\mathcal{C}}$.

######  Lemma 1.1.

The a-priori complex measure $\mathrm{d}\hat{\nu}$ in (1.9) is a probability
measure on the contour $\hat{\mathcal{C}}$ defined in (1.8).

When $t=t_{c}$, the zeros of the orthogonal polynomials (1.3) behave in the
following way.

######  Theorem 1.2.

The zeros of the polynomials $p_{n}(\lambda)$ defined in (1.3) for
$t=t_{c}=\sqrt{T/d}$, behave as follows

  * •

for $n=kd+d-1$, let $\omega={\rm e}^{\frac{2\pi{\rm i}}{d}}$. Then
$t^{\frac{1}{d}},\omega t^{\frac{1}{d}},\dots,\omega^{d-1}t^{\frac{1}{d}}$ are
zeros of the polynomials $p_{kd+d-1}$ with multiplicity $k$ and $\lambda=0$ is
a zero with multiplicity $d-1$.

  * •

for $n=kd+\ell$, $\ell=0,\dots,d-2$ the polynomial $p_{n}(\lambda)$ has a zero
in $\lambda=0$ with multiplicity $\ell$ and the remaining zeros in the limit
$n$, $N\to\infty$ such that

$\displaystyle N=\dfrac{n-\ell}{T},$

accumulate on the level curve $\hat{\mathcal{C}}$ defined in (1.8). The
measure $\hat{\nu}$ in (1.9) is the weak-star limit of the normalized zero
counting measure $\nu_{n}$ of the polynomials $p_{n}$ for $n=kd+\ell$,
$\ell=0,\dots,d-2$.

Our next result gives strong uniform asymptotics as $n\to\infty$ for the
polynomials $p_{n}(\lambda)$ in the complex plane. We consider the double
scaling limit

$\displaystyle t\to t_{c},\qquad k\to\infty,\qquad k=\dfrac{n-\ell}{d},$

in such a way that

$\displaystyle\lim_{k\to\infty,\,t\to
t_{c}}\sqrt{k}\left(\dfrac{t^{2}}{t_{c}^{2}}-1\right)=\mathcal{S},$

with $\mathcal{S}$ in compact subsets of the real line. In the description of
the asymptotic behaviour of the orthogonal polynomials, $p_{n}(\lambda)$, in
this double-scaling limit, the Painlevé IV transcendent (see Section 2) with
$\Theta_{0}=\Theta_{\infty}=\frac{\gamma}{2}$ and
$\gamma=\frac{d-\ell-1}{d}\in(0,1)$ plays a major role.

######  Theorem 1.3 (double scaling limit).

The polynomials $p_{n}(\lambda)$ with $n=kd+\ell$, $\ell=0,\dots,d-2$,
$\gamma=\frac{d-\ell-1}{d}\in(0,1)$, have the following asymptotic behaviour
when $n$, $N\to\infty$ in such a way that $NT=n-\ell$ and

$\displaystyle\lim_{k\to\infty,\,t\to
t_{c}}\sqrt{k}\left(\dfrac{t^{2}}{t_{c}^{2}}-1\right)=\mathcal{S},$

with $\mathcal{S}$ in compact subsets of the real line so that the solution
$Y(\mathcal{S})$ of the Painlevé IV equation (2.3) does not have poles. Below,
the function $Z=Z(\mathcal{S})$, $U=U(\mathcal{S})$ and the Hamiltonian
$H=H(\mathcal{S})$ are related to the Painlevé IV equation (2.3) by the
relations (2.2) and (2.7), respectively.

  * $(1)$

For $\lambda$ in compact subsets of the exterior of $\hat{\mathcal{C}}$ one
has

$\displaystyle
p_{n}(\lambda)=\lambda^{d-1}\big{(}\lambda^{d}-t_{c}\big{)}^{k-\gamma}\left(1-\dfrac{H(\mathcal{S})t_{c}}{\sqrt{k}\lambda^{d}}+\mathcal{O}\left(\frac{1}{k}\right)\right),$

with $H(\mathcal{S})$ the Hamiltonian (2.7) of the Painlevé IV equation (2.3).

  * $(2)$

For $\lambda$ in the region near $\hat{\mathcal{C}}$ and away from the point
$\lambda=0$ one has

$\displaystyle
p_{n}(\lambda)=\lambda^{d-1}\big{(}\lambda^{d}-t_{c}\big{)}^{k-\gamma}$
$\displaystyle\hphantom{p_{n}(\lambda)=}{}\times\left(1-\dfrac{H(\mathcal{S})t_{c}}{\sqrt{k}\lambda^{d}}+\dfrac{Z(\mathcal{S})}{U(\mathcal{S})}\dfrac{t_{c}{\rm
e}^{-k\hat{\varphi}(\lambda)}}{\lambda^{d}k^{\frac{1+\gamma}{2}}}\left(\dfrac{\lambda^{d}-t_{c}}{\lambda^{d}}\right)^{{\gamma}}+\mathcal{O}\left(\frac{1}{k}\right)\right),$
(1.10)

with $\hat{\varphi}(\lambda)$ defined in (1.7).

  * $(3)$

For $\lambda$ in compact subsets of the interior region of $\hat{\mathcal{C}}$
and away from the point $\lambda=0$ one has

$\displaystyle
p_{n}(\lambda)=\lambda^{\ell}\big{(}\lambda^{d}-t_{c}\big{)}^{k}\dfrac{{\rm
e}^{-k\hat{\varphi}(\lambda)}}{k^{\frac{1}{2}+\gamma}}\left(\dfrac{Z(\mathcal{S})}{U(\mathcal{S})}\dfrac{t_{c}}{\lambda^{d}}+\mathcal{O}\left(\frac{1}{k}\right)\right).$

  * $(4)$

In the neighbourhood of the point $\lambda=0$ one has

$\displaystyle
p_{n}(\lambda)=\lambda^{\ell}\big{(}\lambda^{d}-t_{c}\big{)}^{k-\gamma}\dfrac{{\rm
e}^{-\frac{k}{2}\hat{\varphi}(\lambda)}\lambda^{d\gamma}}{(-2k\hat{\varphi}(\lambda))^{\frac{\gamma}{4}}}$
$\displaystyle\hphantom{p_{n}(\lambda)=}{}\times\begin{cases}\Psi_{11}\big{(}\sqrt{-2k\hat{\varphi}(\lambda)};\mathcal{S}\big{)}+\mathcal{O}\left(\frac{1}{k}\right),&\sqrt{-\hat{\varphi}(\lambda)}\in\Omega_{\infty},\\\
-\Psi_{12}\big{(}\sqrt{-2k\hat{\varphi}(\lambda)};\mathcal{S}\big{)}+\mathcal{O}\left(\frac{1}{k}\right),&\sqrt{-\hat{\varphi}(\lambda)}\in\Omega_{0},\\\
\Psi_{11}\big{(}\sqrt{-2k\hat{\varphi}(\lambda)};\mathcal{S}\big{)}-\Psi_{12}\big{(}\sqrt{-2k\hat{\varphi}(\lambda)};\mathcal{S}\big{)}+\mathcal{O}\left(\frac{1}{k}\right),&\sqrt{-\hat{\varphi}(\lambda)}\in\Omega_{2},\end{cases}$

where $\Psi$ is the solution of the Painlevé IV Riemann–Hilbert Problem 2.1
with Stokes multipliers specified in (2.9) and (2.10) and the regions
$\Omega_{0}$, $\Omega_{\infty}$ and $\Omega_{2}$ are specified in Fig. 4. Here
$\Psi_{11}$ and $\Psi_{12}$ are the entries of the matrix $\Psi$.

We observe that, in compact subsets of the exterior of $\hat{\mathcal{C}}$,
there are no zeros of the polynomials $p_{n}(\lambda)$. The only possible
zeros are located in $\lambda=0$ and in the region where the second term in
parenthesis in the expression (1.10) is of order one. Since
$\operatorname{Re}\hat{\varphi}(\lambda)$ is negative inside
$\hat{\mathcal{C}}$ and positive outside $\hat{\mathcal{C}}$, it follows that
the possible zeros of $p_{n}(\lambda)$ lie inside $\hat{\mathcal{C}}$ and are
determined by the condition

$\displaystyle\operatorname{Re}\hat{\varphi}(\lambda)=-\dfrac{1+\gamma}{2}\dfrac{\log
k}{k}+\dfrac{1}{k}\log\left(\left|\dfrac{\lambda^{d}-t_{k}}{\lambda^{d}}\right|^{\gamma}\left|\dfrac{t_{c}Z(\mathcal{S})}{\lambda^{d}U(\mathcal{S})}\right|\right)+\dfrac{1}{k^{\frac{3}{2}}}\operatorname{Re}\left(\dfrac{t_{c}H(\mathcal{S})}{\lambda^{d}}\right)+\mathcal{O}\big{(}k^{-2}\big{)}.$

We conclude with the following proposition.

######  Proposition 1.4.

The support of the counting measure of the zeros of the polynomials
$p_{n}(\lambda)$ in the limit $n\to\infty$ outside an arbitrary small disk
surrounding the point $\lambda=0$ tends uniformly to the curve
$\hat{\mathcal{C}}$ defined in (1.8). The zeros are within a distance
$\mathcal{O}(1/k)$ from the curve defined by (1.8).

This manuscript is organised as follows. In Section 2 we present the Lax pair
and Riemann–Hilbert problem for the Painlevé IV equation and we show how to
associate for a particular value of the Stokes matrices the solution of the
Painlevé IV equation in terms of a Fredholm determinant. We compute such
solution numerically for several values of the parameter $\gamma$ that
parametrizes the asymptotic behaviour of the orthogonal polynomials. In
Section 3 we show how to associate to the orthogonal polynomials on the planar
domain a Riemann–Hilbert problem and in Section 4 we perform the asymptotic
analysis of this Riemann–Hilbert problem thus deriving the asymptotic
behaviour of the orthogonal polynomials $p_{n}(\lambda)$ as $n\to\infty$ in
the complex plane.

##  2 Painlevé IV and Fredholm determinant

The purpose of this section is to introduce the Painlevé IV equation and study
various aspects related to it that are connected with the work that we will be
doing in the following sections.

We will begin by introducing and deriving the general Painlevé IV equation
following the ideas outlined by Miwa, Jimbo and Ueno [30, 31]. In the sequence
of this, we will also discuss the associated Stokes’ phenomenon and the
monodromy problem (see Wasow [42]).

The second part of this section will be devoted to applying the general
setting of Painlevé IV in a particular case, in order to obtain a special
solution for it. This will correspond to a special Riemann–Hilbert problem
whose properties will be outlined. The reason and motivation for introducing
this solution is the fact that this Riemann–Hilbert problem will be seen to be
useful in the following section.

The third and final part of this section will be devoted to the study of the
Fredholm determinant, as it has been established by Its [28], Harnad [25],
Bertola and Cafasso [8], among others. This will be used in order to study the
$\tau$-function of the special solution of Painlevé IV that we introduced
before and analyze its validity with respect to the value of the parameter
$s$.

###  2.1 The general Painlevé IV

We recall the Lax pair formulation of the Painlevé IV equation as in the
original work of Miwa, Jimbo and Ueno [30, 31], with notation adapted to our
purposes.

Consider the matrix–valued function $\Psi$, associated with Painlevé IV, as
the fundamental joint solution of the system of linear differential equations,
the Lax pair, given by

$\displaystyle\frac{\partial}{\partial\lambda}\Psi=A(\lambda;s)\Psi,\qquad\frac{\partial}{\partial
s}\Psi=B(\lambda;s)\Psi,$ (2.1)

where

$\displaystyle
A(\lambda;s)=-\frac{1}{2}\left(\lambda+s\right)\sigma_{3}+\begin{pmatrix}0&\dfrac{Z}{U}\vspace{1mm}\\\
-U&0\end{pmatrix}+\frac{1}{\lambda}\begin{pmatrix}\Theta_{\infty}-Z&\dfrac{(\Theta_{\infty}-Z)^{2}-\Theta^{2}_{0}}{YU}\vspace{1mm}\\\
-UY&Z-\Theta_{\infty}\end{pmatrix},$ $\displaystyle
B(\lambda;s)=-\frac{1}{2}\lambda\sigma_{3}+\begin{pmatrix}0&\dfrac{Z}{U}\vspace{1mm}\\\
-U&0\end{pmatrix},$

where $\sigma_{3}$ is the third Pauli matrix, $Y$, $Z$ and $U$ are functions
of $s$ whose dependence we now derive, while $\Theta_{0}$ and
$\Theta_{\infty}$ are arbitrarily chosen constants. The compatibility of (2.1)
yields the zero-curvature equation

$\displaystyle[\partial_{\lambda}-A,\partial_{s}-B]\equiv 0,$

which, in turn, yields the following system of nonlinear ODEs

$\displaystyle U^{\prime}=U(Y-s),$ $\displaystyle
Z^{\prime}=ZY+\dfrac{-(Z-\Theta_{\infty})^{2}+\Theta^{2}_{0}}{Y},$
$\displaystyle Y^{\prime}=-Y^{2}+sY-2Z+2\Theta_{\infty}.$ (2.2)

In the equations above the prime means differentiation with respect to $s$. By
eliminating $U$, $Z$ from the above equations we obtain the classical form of
the fourth Painlevé equation

$\displaystyle
Y^{\prime\prime}=\frac{1}{2}\frac{(Y^{\prime})^{2}}{Y}+\frac{3}{2}Y^{3}-2sY^{2}+\left(1+\frac{s^{2}}{2}-2\Theta_{\infty}\right)Y-\frac{2\Theta^{2}_{0}}{Y}.$
(2.3)

####  2.1.1 Stokes’ phenomenon and the Monodromy problem

The review material of this section is mostly based on [21] and [30, 31]. The
solution of the first equation in (2.1) admits a formal sectorial expansion at
$\lambda=\infty$ of the form

$\displaystyle\Psi_{\rm
formal}(\lambda)=\left(\mathbf{1}+\frac{\Psi_{1}(s)}{\lambda}+\frac{\Psi_{2}(s)}{\lambda^{2}}+\mathcal{O}\left(\frac{1}{\lambda^{3}}\right)\right)\lambda^{\Theta_{\infty}\sigma_{3}}{\rm
e}^{-\theta(\lambda)\sigma_{3}},$
$\displaystyle\Psi_{1}(s)=\begin{bmatrix}H(s)&\dfrac{Z(s)}{U(s)}\vspace{1mm}\\\
U(s)&-H(s)\end{bmatrix},$ (2.4)
$\displaystyle\Psi_{2}(s)=\left[\begin{matrix}\dfrac{1}{2}\big{(}{H(s)^{2}+Z(s)-sH(s)}\big{)}&{\dfrac{Z(s)(Z(s)-\gamma-Y(s)s-H(s)Y(s))}{U(s)Y(s)}}\vspace{1mm}\\\
U(s)H(s)+U(s)Y(s)-sU(s)&\dfrac{1}{2}\big{(}{H(s)^{2}+Z(s)+sH(s)}\big{)}\end{matrix}\right],$
(2.5) $\displaystyle\theta(\lambda)=\frac{\lambda^{2}}{4}+\frac{s}{2}\lambda,$
(2.6)

where

$\displaystyle
H=\left(s+\frac{2\Theta_{\infty}}{Y}-Y\right)Z-\frac{\Theta_{\infty}^{2}-\Theta_{0}^{2}}{Y}-\frac{Z^{2}}{Y}.$
(2.7)

In each of the sectors $q=(0,{\rm I},{\rm II},{\rm III},{\rm IV},{\rm V})$ of
Fig. 3 there exists a unique solution $\Psi_{(q)}(\lambda)$ which is Poincaré
asymptotic to $\Psi_{\textrm{formal}}(\lambda)$, even though the latter is
only a non-convergent (in general) formal series. The matrices in Fig. 3 on
the rays $\Gamma_{q}$ are the Stokes’ matrices $\mathbf{S}_{q}$ (and the
formal monodromy on the negative axis) and relate the different solutions as
follows: $\Psi_{(q+1)}=\Psi_{q}{\mathbf{S}}_{q}$ ($q$ modulo $6$). The rays
$\Gamma_{1,5}$ can be chosen arbitrarily within the respective quadrants, or
more generally any smooth infinite arc within the quadrants that admits an
asymptotic direction.

Figure 3: Stokes’ phenomenon and the corresponding jumps in the
$\Psi(\lambda)$ Riemann–Hilbert problem.

Each of these solutions of Painlevé IV for each sector $q$,
$\Psi_{(q)}(\lambda)$, has an asymptotic expansion near $\lambda=0$ with the
following behaviour

$\displaystyle\Psi_{(q)}(\lambda)=G(\mathbf{1}+\mathcal{O}(\lambda))\lambda^{\Theta_{0}\sigma_{3}}C_{(q)},$

where the matrices $C_{(q)}$ defined by this equation are called “connection
matrices” and $G$ is a diagonalizing matrix for the $\frac{1}{\lambda}$
coefficient of $A(\lambda;s)$. From this equation and using an argument of
analytic continuation, we obtain the following monodromy relation

$\displaystyle C^{-1}_{(0)}{\rm e}^{2\pi{\rm
i}\Theta_{0}\sigma_{3}}C_{(0)}=\begin{pmatrix}1&-S_{-2}\\\
0&1\end{pmatrix}\begin{pmatrix}1&0\\\ -S_{1}&1\end{pmatrix}$
$\displaystyle\hphantom{C^{-1}_{(0)}{\rm e}^{2\pi{\rm
i}\Theta_{0}\sigma_{3}}C_{(0)}=}{}\times\begin{pmatrix}1&-S_{0}\\\
0&1\end{pmatrix}\begin{pmatrix}1&0\\\
-S_{-1}&1\end{pmatrix}\begin{pmatrix}1&-S_{-2}\\\ 0&1\end{pmatrix}{\rm
e}^{2\pi{\rm i}\Theta_{\infty}\sigma_{3}}.$ (2.8)

The remaining connection matrices $C_{(q)}$ are obtained from $C_{(0)}$ by
multiplying it by an appropriate sequence of the Stokes’ matrices; for example
$C_{\left(\mathrm{V}\right)}=C_{(0)}\left(\begin{smallmatrix}1&S_{-2}\\\
0&1\end{smallmatrix}\right)$.

The general solution to the isomonodromic equations (2.2) for given
$\Theta_{0}$, $\Theta_{\infty}$ is parametrized by the choices of parameters
$S_{-1},\dots$ and connection matrix $C_{(0)}$ that satisfy (2.8).

###  2.2 Special Riemann–Hilbert problem $\boldsymbol{\Psi(\lambda)}$ and
Painlevé IV

Several particular solutions of the Painlevé IV have been considered in the
literature, see, e.g., [7, 14, 29]. In our work we only need special values
for the Stokes’ parameters and $\Theta_{0,\infty}$:

$\displaystyle S_{1}=1,\qquad S_{0}=0,\qquad S_{-1}=-1,\qquad S_{-2}=-1,$
(2.9) $\displaystyle
C_{(\textrm{III})}=C_{(\textrm{IV})}=\mathbf{1},\qquad\Theta_{0}=\frac{\gamma}{2},\qquad\Theta_{\infty}=\frac{\gamma}{2}.$
(2.10)

The resulting Riemann–Hilbert problem is the one depicted in Fig. 4, where we
have more appropriately renamed $\Gamma_{1}$ and $\Gamma_{\infty}$. It is
formalized below:

Figure 4: Jumps in the $\Psi(\lambda)$ Riemann–Hilbert problem.

######  Riemann–Hilbert Problem 2.1.

The matrix $\Psi(\lambda;s)$ is analytic in
$\mathbb{C}\setminus(\Gamma_{1}\cup\Gamma_{\infty}\cup\mathbb{R}_{-})$ and
admits non-tangential boundary values. Moreover

  1. $1.$

The boundary values are related by $($see Fig. 4):

$\displaystyle\Psi_{+}(\lambda;s)=\Psi_{-}(\lambda;s)v_{\Psi},\qquad\lambda\in\Sigma_{\Psi},$
$\displaystyle v_{\Psi}=\begin{cases}\begin{pmatrix}1&-1\\\
0&1\end{pmatrix},&\lambda\in\Gamma_{1},\vspace{1mm}\\\ \begin{pmatrix}1&0\\\
1&1\end{pmatrix},&\lambda\in\Gamma_{\infty},\\\ {\rm e}^{\gamma\pi{\rm
i}\sigma_{3}},\quad&\lambda\in\mathbb{R}_{-},\end{cases}$

where $\Psi_{\pm}$ are the boundary values on the left and right of the
oriented contour $\Sigma_{\Psi}$.

  2. $2.$

Near $\lambda=\infty$ the matrix has the following sectorial behaviour

$\displaystyle\Psi(\lambda;s)=\left(\mathbf{1}+\mathcal{O}\left(\frac{1}{\lambda}\right)\right)\lambda^{\frac{\gamma}{2}\sigma_{3}}{\rm
e}^{-\theta(\lambda)\sigma_{3}},\qquad\lambda\rightarrow\infty.$

  3. $3.$

Near $\lambda=0$ in the region $\Omega_{\infty}$, $\Psi(\lambda;s)$ has the
behaviour

$\displaystyle\Psi(\lambda;s)=\mathcal{O}(1)\lambda^{\frac{\gamma}{2}\sigma_{3}}.$

We remind the reader that subcripts ± denote the boundary values from the left
($+$) or right ($-$) of an oriented countour. We will henceforth omit explicit
notation for the $s$ dependence of $\Psi$.

We now proceed to do a further modification of this
$\Psi(\lambda)$-Riemann–Hilbert problem, that will be used in the next section
for the study of the Fredholm determinant:

  * •

move the jump contour $\Gamma_{\infty}={\rm i}\mathbb{R}$ in the following way
${\rm i}\mathbb{R}\rightarrow{\rm i}\mathbb{R}+\varepsilon$, with
$\epsilon>0$;

  * •

join together (collapse) the two parts of the contour $\Gamma_{1}$ and
$\mathbb{R}_{-}$.

The use of the verb “move” means that we define a new matrix-valued function
in the region between the old and new position of the contours that
analytically continues $\Psi$ across the original contours. This is possible
because the jump matrices are all constants in $\lambda$. We still use the
symbol $\Psi$ for the resulting function. Finally, we introduce $M$ as follows

$\displaystyle M(\lambda):=\Psi(\lambda){\rm
e}^{\theta(\lambda)\sigma_{3}}\lambda^{-\frac{\gamma}{2}\sigma_{3}}.$

The powers of $\lambda$ are intended as the principal determination with
$\arg(\lambda)\in(-\pi,\pi)$. The matrix $M$ solves a Riemann–Hilbert problem
with jumps indicated in Fig. 5 and formalized below.

Figure 5: Jumps in the $M(\lambda)$ Riemann–Hilbert problem.

######  Riemann–Hilbert Problem 2.2.

The matrix $M(\lambda)$ is analytic in
$\mathbb{C}\setminus(\Gamma_{\infty}\cup\mathbb{R}_{-})$ and admits non-
tangential boundary values. Moreover

  1. $1.$

The boundary values $M_{\pm}$ are related by $($see Fig. 5):

$\displaystyle M_{+}(\lambda)=M_{-}(\lambda)v_{M},\qquad\lambda\in\Sigma_{M},$
$\displaystyle v_{M}=\begin{cases}\begin{pmatrix}1&2{\rm
i}\sin{(\pi\gamma)}|\lambda|^{\gamma}{\rm e}^{-2\theta(\lambda)}\\\
0&1\end{pmatrix},&\lambda\in\mathbb{R}_{-},\vspace{1mm}\\\
\begin{pmatrix}1&0\\\ \lambda^{-\gamma}{\rm
e}^{2\theta(\lambda)}&1\end{pmatrix},&\lambda\in\Gamma_{\infty}.\end{cases}$

  2. $2.$

Near $\lambda=\infty$ the matrix has the following sectorial behaviour:

$\displaystyle
M(\lambda)=\left(\mathbf{1}+\mathcal{O}\left(\frac{1}{\lambda}\right)\right),\qquad\lambda\rightarrow\infty.$

  3. $3.$

The matrix $M(\lambda)$ is bounded in a full neighbourhood of $\lambda=0$.

###  2.3 Fredholm determinant and the $\boldsymbol{\tau}$-function of Painlevé
IV

The Riemann–Hilbert Problem 2.2 is associated to an integral operator on
$L^{2}(\mathbb{R}_{-}\cup\Gamma_{\infty})$, where $\Gamma_{\infty}={\rm
i}{\mathbb{R}}+\epsilon$, and will be used in the next section in order to
compute the Fredholm determinant. We first introduce some general definitions.
Let $(X,\mathrm{d}\nu)$ be a measurable space and $K(x,y)\colon X\times
X\to{\mathbb{C}}$ a function. We define the operator $\mathcal{\hat{K}}\colon
L^{2}(X,\mathrm{d}\nu)\rightarrow L^{2}(X,\mathrm{d}\nu)$ as

$\displaystyle\big{(}\mathcal{\hat{K}}f\big{)}(x)=\int_{X}K(x,y)f(y)\mathrm{d}\nu(y).$

If $\iint_{X^{2}}|K(x,y)|^{2}{\rm d}x{\rm d}y<+\infty$, the operator
$\hat{\mathcal{K}}$ is a Hilbert–Schmidt operator and its spectrum is purely
discrete, the eigenvalues can accumulate only at $0$ and all the
multiplicities of the nonzero eigenvalues are finite (see, e.g., [39]).

######  Definition 2.3.

The $\tau$-function is defined as the following Fredholm determinant

$\displaystyle\tau(\rho)=\det{\big{(}\operatorname{Id}-\rho\mathcal{\hat{K}}\big{)}},$

where $\rho$ is a nonzero complex parameter.

From the properties of the Fredholm determinant we know that $\tau(\rho)=0$ if
and only if $\frac{1}{\rho}$ is an eigenvalue of $\mathcal{\hat{K}}$. With
this being said, since in our case we are interested in a different parametric
dependence of $\mathcal{\hat{K}}$ and not its spectral properties, we shall
set $\rho=1$ (see (2.13)).

####  2.3.1 Riemann–Hilbert problems and Fredholm determinants

Following [28] let us consider a set of contours $\Sigma\subset\mathbb{C}$ and
$f,g\colon\Sigma\to\mathbb{C}^{p}$ matrix–valued (smooth) functions that
satisfy the condition $f^{\rm T}(\lambda)g(\lambda)=0$. Then we define the
following scalar kernel

$\displaystyle K(\lambda,\mu):=\frac{f^{\rm T}(\lambda)g(\mu)}{\lambda-\mu}.$

We denote by $\hat{\mathcal{K}}\colon L^{2}(\Sigma,\mathbb{C}^{p})\rightarrow
L^{2}(\Sigma,\mathbb{C}^{p})$ the integral operator with kernel
$K(\lambda,\mu)$. Let $\mathcal{R}$ be the resolvent operator
$\mathcal{R}:=-\hat{\mathcal{K}}\circ(\operatorname{Id}-\hat{\mathcal{K}})^{-1}.$

######  Theorem 2.4 (see [28]).

The kernel, $R(\lambda,\mu)$, of the resolvent operator $\mathcal{R}$ is given
by

$\displaystyle R(\lambda,\mu):=\frac{f^{\rm T}(\lambda)A^{\rm
T}(\lambda)A^{-{\rm T}}(\mu)g(\mu)}{\lambda-\mu},$

where the $p\times p$ matrix $A(\lambda)$ solves the Riemann–Hilbert problem

$\displaystyle A_{+}(\lambda)=A_{-}(\lambda)\big{(}\mathbf{1}-2\pi{\rm
i}f(\lambda)g^{\rm T}(\lambda)\big{)},\qquad\lambda\in\Sigma,$ $\displaystyle
A(\lambda)=\mathbf{1}+\mathcal{O}\big{(}\lambda^{-1}\big{)},\qquad\lambda\rightarrow\infty.$

Moreover, the above Riemann–Hilbert problem admits a unique solution if and
only if the Fredholm determinant
$\det\big{(}{\mathrm{Id}}-\hat{\mathcal{K}}\big{)}$ is non-zero.

The case of interest to us, below, is with $p=2$.

####  2.3.2 Fredholm determinant and the Riemann–Hilbert problem
$\boldsymbol{M}$

The Riemann–Hilbert Problem 2.2 is a special example of the previous setup, as
we now see. We first define the following vectors

$\displaystyle f(z)=\begin{bmatrix}2{\rm
i}\sin(\pi\gamma)|z|^{\frac{{\gamma}}{2}}{\rm
e}^{-\theta(z)}\chi_{\mathbb{R}_{-}}(z)\\\ z^{-\frac{{\gamma}}{2}}{\rm
e}^{\theta(z)}\chi_{\infty}(z)\end{bmatrix},$ (2.11) $\displaystyle
g(z)=-\frac{1}{2\pi{\rm i}}\begin{bmatrix}z^{-\frac{{\gamma}}{2}}{\rm
e}^{\theta(z)}\chi_{\infty}(z)\\\ |z|^{\frac{{\gamma}}{2}}{\rm
e}^{-\theta(z)}\chi_{\mathbb{R}_{-}}(z)\end{bmatrix},$ (2.12)

where $\chi_{\infty}$ and $\chi_{\mathbb{R}_{-}}$ are the characteristic
functions of $\Gamma_{\infty}$ and $\mathbb{R}_{-}$, respectively. This leads
to the following jump matrices for Theorem 2.4

$\displaystyle\mathbf{1}-2\pi{\rm i}f(z)g^{\rm
T}(z)=\begin{cases}\begin{pmatrix}1&2{\rm i}\sin{(\pi\gamma)}|z|^{\gamma}{\rm
e}^{-2\theta(z)}\\\ 0&1\end{pmatrix},&z\in\mathbb{R}_{-},\vspace{1mm}\\\
\begin{pmatrix}1&0\\\ z^{-\gamma}{\rm
e}^{2\theta(z)}&1\end{pmatrix},&z\in\Gamma_{\infty},\end{cases}$

which match those in the Riemann–Hilbert Problem 2.2 ($\mathbf{1}$ denotes the
$2\times 2$ identity matrix).

Let $\Sigma=\Gamma_{\infty}\cup\mathbb{R}_{-}$, as a consequence of this, we
can now write the Riemann–Hilbert Problem (2.2) for the matrix $M$ in the
following way

$\displaystyle M_{+}(z)=M_{-}(z)\big{(}\mathbf{1}-2\pi{\rm i}f(z)g^{\rm
T}(z)\big{)},\qquad z\in\Sigma,$ $\displaystyle
M=\mathbf{1}+\mathcal{O}\left(\frac{1}{z}\right),\qquad z\rightarrow\infty.$

It should also be noted that, due to our construction, we have $f^{\rm
T}(z)g(z)\equiv 0$.

Let us consider the operator $\hat{\mathcal{K}}\colon L^{2}(\Sigma)\rightarrow
L^{2}(\Sigma)$ associated to the above Riemann–Hilbert problem. Due to the
decomposition of
$L^{2}(\Sigma)=L^{2}(\mathbb{R_{-}}\cup\Gamma_{\infty})=L^{2}(\mathbb{R_{-}})\oplus
L^{2}(\Gamma_{\infty})$ and $\varphi=\varphi_{0}\oplus\varphi_{1}$ with
$\varphi_{0}\in L^{2}(\mathbb{R}_{-})$ and $\varphi_{1}\in
L^{2}(\Gamma_{\infty})$, we can write the operator $\hat{{\cal K}}$ in the
form

$\displaystyle\big{(}\hat{\mathcal{K}}\varphi\big{)}(z)=\int_{\Sigma}\frac{f^{\rm
T}(z)g(w)}{z-w}\varphi(w)\mathrm{d}w$
$\displaystyle\hphantom{\big{(}\hat{\mathcal{K}}\varphi\big{)}(z)}{}=-\frac{2{\rm
i}\sin{(\pi\gamma)}}{2\pi{\rm
i}}\left(\int_{\Gamma_{\infty}}\frac{w^{-\frac{\gamma}{2}}{\rm
e}^{\theta(w)}\varphi_{1}(w)}{z-w}\mathrm{d}w\right)|z|^{\frac{\gamma}{2}}{\rm
e}^{-\theta(z)}\chi_{\mathbb{R}_{-}}(z)$
$\displaystyle\hphantom{\big{(}\hat{\mathcal{K}}\varphi\big{)}(z)=}{}-\frac{1}{2\pi{\rm
i}}\left(\int_{\mathbb{R}_{-}}\frac{|w|^{\frac{\gamma}{2}}{\rm
e}^{-\theta(w)}\varphi_{0}(w)}{z-w}\mathrm{d}w\right)z^{-\frac{\gamma}{2}}{\rm
e}^{\theta(z)}\chi_{\Gamma_{\infty}}(z).$ (2.13)

Due to this fact, the operator $\hat{\mathcal{K}}$ can now be written in block
form $\hat{\mathcal{K}}_{ij}$ with

$\displaystyle\hat{\mathcal{K}}_{00}\colon\ L^{2}(\mathbb{R}_{-})\rightarrow
L^{2}(\mathbb{R}_{-}),$ $\displaystyle\hat{\mathcal{K}}_{01}\colon\
L^{2}(\Gamma_{\infty})\rightarrow L^{2}(\mathbb{R}_{-}),$
$\displaystyle\hat{\mathcal{K}}_{10}\colon\ L^{2}(\mathbb{R}_{-})\rightarrow
L^{2}(\Gamma_{\infty}),$ $\displaystyle\hat{\mathcal{K}}_{11}\colon\
L^{2}(\Gamma_{\infty})\rightarrow L^{2}(\Gamma_{\infty}).$

In our case, due to the form of the vectors $f$ and $g$ in (2.11), (2.12), the
diagonal blocks are null operators:
$\hat{\mathcal{K}}_{00}=\hat{\mathcal{K}}_{11}=0$, while

$\displaystyle\hat{{\cal K}}_{01}\colon\ L^{2}(\Gamma_{\infty})\rightarrow
L^{2}(\mathbb{R}_{-}),\qquad\big{(}\hat{{\cal
K}}_{01}\varphi_{1}\big{)}(z)=\int_{\Gamma_{\infty}}K_{0,1}(z,w)\varphi_{1}(w)\mathrm{d}w,$
$\displaystyle K_{01}(z,w)=-\frac{2{\rm i}\sin{(\pi\gamma)}}{(2\pi{\rm
i})}\frac{w^{-\frac{\gamma}{2}}|z|^{\frac{\gamma}{2}}{\rm
e}^{\theta(w)-\theta(z)}}{z-w},$ $\displaystyle\hat{{\cal K}}_{10}\colon\
L^{2}(\mathbb{R}_{-})\rightarrow
L^{2}(\Gamma_{\infty}),\qquad\big{(}\hat{{\cal
K}}_{10}\varphi_{0}\big{)}(z)=\int_{R_{-}}K_{1,0}(z,w)\varphi_{0}(w)\mathrm{d}w,$
$\displaystyle K_{10}(z,w)=-\dfrac{K_{01}(w,z)}{2{\rm
i}\sin(\pi\gamma)}=-\frac{1}{(2\pi{\rm
i})}\frac{|w|^{\frac{\gamma}{2}}z^{-\frac{\gamma}{2}}{\rm
e}^{\theta(z)-\theta(w)}}{z-w}.$

Following the Definition 2.3 of the tau function, we can now compute it for
the case when $\hat{\mathcal{K}}$ is given by this operator

$\displaystyle\tau=\det{\big{(}\operatorname{Id}-\mathcal{\hat{K}}\big{)}}=\det{\left(\begin{bmatrix}\mathbf{1}_{0}&0\\\
0&\mathbf{1}_{1}\end{bmatrix}-\begin{bmatrix}0&\hat{{\cal K}}_{01}\\\
\hat{{\cal
K}}_{10}&0\end{bmatrix}\right)}=\det{\left(\begin{bmatrix}\mathbf{1}_{0}&\hat{{\cal
K}}_{01}\\\
0&\mathbf{1}_{1}\end{bmatrix}\begin{bmatrix}\mathbf{1}_{0}&-\hat{{\cal
K}}_{01}\\\ -\hat{{\cal K}}_{10}&\mathbf{1}_{1}\end{bmatrix}\right)}$
$\displaystyle\hphantom{\tau=\det{\big{(}\operatorname{Id}-\mathcal{\hat{K}}\big{)}}}{}=\det{\begin{bmatrix}\mathbf{1}_{0}-\hat{{\cal
K}}_{01}\hat{{\cal K}}_{10}&0\\\ -\hat{{\cal
K}}_{10}&\mathbf{1}_{1}\end{bmatrix}},$

where $\mathbf{1}_{0}=\operatorname{Id}_{L^{2}(\mathbb{R}_{-})}$,
$\mathbf{1}_{1}=\operatorname{Id}_{L^{2}\left(\Gamma_{\infty}\right)}$. We
thus see that the determinant of $\mathcal{\hat{K}}$ can be written as the
Fredholm determinant of the operator on $L^{2}({\mathbb{R}}_{-})$ given by
$\hat{\cal K}_{01}\hat{\cal K}_{10}$ and we are thus left with the equality

$\displaystyle\det{(\operatorname{Id}-\mathcal{\hat{K}})}=\det\big{[}\operatorname{Id}_{L^{2}(\mathbb{R}_{-})}-\hat{{\cal
K}}_{01}\hat{{\cal K}}_{10}\big{]}.$ (2.14)

We define $\hat{\cal{F}}:=\hat{{\cal K}}_{01}\hat{{\cal K}}_{10}\colon
L^{2}({\mathbb{R}}_{-})\to L^{2}({\mathbb{R}}_{-})$. The kernel of this
operator is is finally written as

$\displaystyle\big{(}\hat{\cal{F}}\varphi_{0}\big{)}(z):=\big{(}\hat{{\cal
K}}_{01}\hat{{\cal
K}}_{10}\varphi_{0}\big{)}(z)=\int_{\mathbb{R}_{-}}F(z,w)\varphi_{0}(w)\mathrm{d}w,$
(2.15) $\displaystyle
F(z,w)=\int_{\Gamma_{\infty}}K_{0,1}(z,x)K_{1,0}(x,w){\rm d}x$
$\displaystyle\hphantom{F(z,w)}{}=\frac{2{\rm
i}\sin{(\pi\gamma)}}{4\pi^{2}}|zw|^{\frac{\gamma}{2}}{\rm
e}^{-\theta(w)-\theta(z)}\int_{\Gamma_{\infty}}\frac{x^{-\gamma}{\rm
e}^{2\theta(x)}}{(z-x)(w-x)}\mathrm{d}x,$

where we recall that $\theta(x)=x^{2}/4+sx/2$ so that the operator
$\hat{\cal{F}}=\hat{\cal{F}}(\gamma,s)$ depends on the parameters $s$ and
$\gamma$. We observe that
$\det{\big{(}\operatorname{Id}-\mathcal{\hat{K}}\big{)}}$ is non zero as long
as the sup-norm of the operator $\hat{\cal{F}}$ is less than one. The values
of $s$ which satisfy this condition correspond to a regular solution of the
Painlevé IV equation. We also observe that the kernel $F(z,w)$ is real-valued
(and symmetric) and hence the tau function is real-valued.

####  2.3.3 Norm estimate for small $\boldsymbol{s}$

For $s$ negative and sufficiently large in absolute value, we can estimate the
norm of the operator $\hat{{\cal F}}=\hat{{\cal K}}_{01}\hat{{\cal K}}_{10}$
and guarantee that the determinant is non-zero. Introducing the quantities

$\displaystyle h(z):=\int_{\Gamma_{\infty}}{\rm i}\frac{x^{-\gamma}{\rm
e}^{\frac{x^{2}}{2}+{s}x}}{z-x}\mathrm{d}x,\qquad
a(z)=|z|^{\frac{\gamma}{2}}{\rm e}^{-\frac{z^{2}}{4}-\frac{s}{2}z},\qquad
c_{0}=\frac{\sin{(\pi\gamma)}}{2\pi},$ (2.16)

we can write the operator $\hat{\cal{F}}$ in the form

$\displaystyle\big{(}\hat{\cal{F}}\varphi_{0}\big{)}(z)=\int_{\mathbb{R}_{-}}F(z,w)\varphi_{0}(w)\mathrm{d}w=c_{0}a(z)\int_{\mathbb{R}_{-}}a(w)\varphi_{0}(w)\frac{h(z)-h(w)}{w-z}\frac{\mathrm{d}w}{\pi}.$

Next, we introduce the Hilbert transform
$(\mathcal{P}\phi)(z)=\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\int_{{\mathbb{R}}_{-}}\frac{\phi(w)}{z-w}\frac{{\rm
d}w}{\pi}$ which is a bounded operator from $L^{2}({\mathbb{R}}_{-})$ to
itself with norm one:
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathcal{P}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}=1$.
Here and below, we denote by

$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|A\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}=\sup_{f\in
L^{2}({\mathbb{R}}_{-}),\,f\not\equiv
0}\dfrac{\|Af\|_{L^{2}({\mathbb{R}}_{-})}}{\|f\|_{L^{2}({\mathbb{R}}_{-})}}.$

Given an essentially bounded function $f\in L^{\infty}({\mathbb{R}}_{-})$, we
denote by $M_{f}$ the corresponding multiplication operator and recall that
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|M_{f}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}=\|f\|_{\infty}$.
This leads to

$\displaystyle\big{\|}\big{(}\hat{\cal{F}}\varphi_{0}\big{)}(z)\big{\|}_{L^{2}}=\left\|\int_{\mathbb{R}_{-}}F(z,w)\varphi_{0}(w)\mathrm{d}w\right\|_{L^{2}}=|c_{0}|\|(M_{ah}\mathcal{P}M_{a}-M_{a}\mathcal{P}M_{ah})(\varphi_{0})\|_{L^{2}}$
$\displaystyle\hphantom{\big{\|}\big{(}\hat{\cal{F}}\varphi_{0}\big{)}(z)\big{\|}_{L^{2}}=\left\|\int_{\mathbb{R}_{-}}F(z,w)\varphi_{0}(w)\mathrm{d}w\right\|_{L^{2}}}{}\leq|c_{0}|\big{(}\|M_{ah}\mathcal{P}M_{a}\varphi_{0}\|_{L^{2}}+\|M_{a}\mathcal{P}M_{ah}\varphi_{0}\|_{L^{2}}\big{)}$
$\displaystyle\hphantom{\big{\|}\big{(}\hat{\cal{F}}\varphi_{0}\big{)}(z)\big{\|}_{L^{2}}=\left\|\int_{\mathbb{R}_{-}}F(z,w)\varphi_{0}(w)\mathrm{d}w\right\|_{L^{2}}}{}\leq|c_{0}|\|\varphi_{0}\|_{L^{2}}\big{(}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|M_{ah}\mathcal{P}M_{a}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}+{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|M_{a}\mathcal{P}M_{ah}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\big{)}.$

As a result of this procedure, we obtain the following estimate for the norm
of $\hat{\cal{F}}$

$\displaystyle\big{|}\big{|}\big{|}\hat{\cal{F}}\big{|}\big{|}\big{|}\leq\frac{\sin{(\pi\gamma)}}{\pi}\|a\|^{2}_{\infty}\|h\|_{\infty}.$

A calculus exercise shows that the function $a(z)$ in (2.16), satisfies
$\|a\|_{\infty}\leq 1$, for $z\in{\mathbb{R}}_{-}$, $\gamma\in(0,1)$, ${s}<0$.
We now estimate the values of $\|h\|_{\infty}$. By setting $x=\varepsilon+{\rm
i}t$, $t\in\mathbb{R}$ in the formula for $h(z)$ in (2.16) and $s<0$ and
$\varepsilon=\max(0,-{s})$ we have

$\displaystyle|h(z)|\leq\int_{{\rm
i}\mathbb{R}-{s}}|\mathrm{d}x|\frac{|x|^{-\gamma}{\rm
e}^{\operatorname{Re}\big{(}\frac{x^{2}}{2}+{s}x\big{)}}}{|z-x|}\qquad(|z-x|\geq-{s})$
$\displaystyle\hphantom{|h(z)|}{}\leq\frac{1}{|s|}\int_{{\rm
i}\mathbb{R}-{s}}|\mathrm{d}x||x|^{-\gamma}{\rm
e}^{\operatorname{Re}\big{(}\frac{x^{2}}{2}+{s}x\big{)}}\leq\frac{1}{|s|}\int_{\mathbb{R}}\mathrm{d}t|{-}{s}+{\rm
i}t|^{-\gamma}{\rm e}^{-\frac{t^{2}}{2}-\frac{{s}^{2}}{2}}$
$\displaystyle\hphantom{|h(z)|}{}\leq\frac{{\rm
e}^{-\frac{{s}^{2}}{2}}}{|s|}|s|^{-\gamma}\int_{\mathbb{R}}\mathrm{d}t{\rm
e}^{-\frac{t^{2}}{2}}=\frac{{\rm
e}^{-\frac{{s}^{2}}{2}}}{|s|^{1+\gamma}}\sqrt{2\pi}.$

The estimate for $\big{|}\big{|}\big{|}\hat{\cal{F}}\big{|}\big{|}\big{|}$
becomes

$\displaystyle\big{|}\big{|}\big{|}\hat{\cal{F}}\big{|}\big{|}\big{|}\leq\frac{\sin{(\pi\gamma)}}{\pi}\frac{{\rm
e}^{-\frac{{s}^{2}}{2}}}{|s|^{1+\gamma}}\sqrt{2\pi}.$ (2.17)

At this point, we need to see for which values of ${s}<0$ the norm of
$\hat{\cal{F}}$ is smaller than $1$, which guarantees that the determinant
(2.14) is not zero and hence that our Riemann–Hilbert Problem 2.2 admits a
solution. From (2.17) we have (recall $\gamma\in(0,1)$)

$\displaystyle\big{|}\big{|}\big{|}\hat{\cal{F}}\big{|}\big{|}\big{|}\leq\frac{\sin{(\pi\gamma)}}{\pi}\frac{\rm{e}^{-\frac{{s}^{2}}{2}}}{|s|^{1+\gamma}}\sqrt{2\pi}\leq\begin{cases}\dfrac{{\sqrt{2}}}{\sqrt{\pi}}\dfrac{{\rm
e}^{-\frac{s^{2}}{2}}}{s^{2}},&{s}\in[-1,0),\vspace{1mm}\\\
\dfrac{{\sqrt{2}}}{\sqrt{\pi}}\dfrac{{\rm
e}^{-\frac{s^{2}}{2}}}{|s|},&{s}<-1,\end{cases}$

and we can easily see that the norm $\hat{\cal{F}}$ is less than one for
${s}<{s}_{0}$ where ${s}_{0}\simeq-0.7701449782$. In summary, we have proven

######  Theorem 2.5.

The Riemann–Hilbert Problem 2.2 admits a solution for $s\in(-\infty,{s}_{0})$,
$s_{0}\simeq-0.7701449782$. In particular, the solution of the fourth Painlevé
equation (2.3) for our choice of monodromy data (2.9), (2.10) is pole-free
within that range.

It would be desirable to show that the solution of the Painlevé equation (2.3)
for our choice of monodromy data (2.9), (2.10) is pole-free on the whole real
line. However, the numerical analysis which will be performed in the following
section, shows that there is a discrete (in principle infinite) number of
values of $s$ for which the Fredholm determinant (2.14) vanishes and hence the
solution of (2.3) has poles.

###  2.4 Gaussian quadrature and numerics

The following is an exposition, for the benefit of the reader, of the main
idea of numerical evaluation of Fredholm determinants based on Nystrom’s
method, as explained by Bornemann [10]. The simple idea is to suitably
“discretize” the integral operator using appropriate Gaussian quadrature (see
[40]).

More specifically, the type of quadrature formulæ that we will use are the so-
called “Gauss–Hermite” quadratures, which are defined in the following way

######  Definition 2.6.

Considering an integral of the form

$\displaystyle\int^{\infty}_{-\infty}f(x){\rm e}^{-\Lambda x^{2}}\mathrm{d}x,$

the Gauss–Hermite quadrature is an approximation of the value of this integral
and is defined as

$\displaystyle\int^{\infty}_{-\infty}f(x){\rm e}^{-\Lambda
x^{2}}\mathrm{d}x\simeq\sum^{m}_{i=1}f\big{(}x^{(m)}_{i}\big{)}w^{(m)}_{i},$
(2.18)

where $m$ is the number of points used in this approximation: the points
$\big{\\{}x_{j}^{(m)}\big{\\}}_{j=1}^{m}$ are called the nodes and the
coefficients $\big{\\{}w_{j}^{(m)}\big{\\}}_{j=1}^{m}$ are called the weights
of the quadrature rule. The nodes $x^{(m)}_{j}$ are the roots of the $m$-th
Hermite polynomial $H_{m}\big{(}\sqrt{\Lambda}x\big{)}$ and the weights
$w^{(m)}_{i}$ are given by

$\displaystyle
w^{(m)}_{i}=\frac{2^{m-1}m!\sqrt{\pi}}{\sqrt{\Lambda}m^{2}\big{(}H_{m-1}\big{(}\sqrt{\Lambda}x^{(m)}_{i}\big{)}\big{)}^{2}}.$

We will have to adapt this definition so that it is suitable to be applied to
our case.

In order to use this technique for the computation of $F(x,y)$ in (2.15) we
rewrite it as

$\displaystyle F(x,y)={\rm e}^{\frac{-x^{2}-y^{2}}{4}}H(x,y),$ $\displaystyle
H(x,y):=\frac{\sin{(\pi\gamma)}}{2\pi^{2}}|x|^{\frac{\gamma}{2}}|y|^{\frac{\gamma}{2}}{\rm
e}^{-\frac{s}{2}(x+y)}\int_{\mathbb{R}}\mathrm{d}t\frac{({\rm
i}t+\varepsilon)^{-\gamma}{\rm e}^{-\frac{t^{2}}{2}}{\rm e}^{{\rm
i}t(s+\varepsilon)}}{(x-\varepsilon-{\rm i}t)(y-\varepsilon-{\rm i}t)}{\rm
e}^{\frac{\varepsilon^{2}}{2}+s\varepsilon}.$

Then Nystrom’s method states that the Fredholm determinant is approximated by

$\displaystyle\det{}_{L^{2}(\mathbb{R}_{-})}\big{[}\operatorname{Id}_{L^{2}(\mathbb{R}_{-})}-\hat{\cal{F}}\big{]}\simeq\det{}_{n\times
n}\left[\operatorname{Id}_{n}-\begin{bmatrix}H\big{(}x^{(2n)}_{i},x^{(2n)}_{j}\big{)}\sqrt{w^{(2n)}_{i}w^{(2n)}_{j}}\end{bmatrix}^{n}_{i,j=1}\right],$

where the nodes $x^{(2n)}_{i}$ and the weights $w^{(2n)}_{i}$ are chosen as in
(2.18) with $\Lambda=\frac{1}{2}$ and where we are selecting only half of the
nodes, namely $n$ out of $2n$ that lie on the negative real axis. Therefore,
the matrix $\mathbb{H}=[H_{jk}]_{j,k=1}^{n}$ that we are interested in
computing is given by

$\displaystyle\begin{bmatrix}H_{jk}\end{bmatrix}^{n}_{j,k=1}:=\begin{bmatrix}H\big{(}x^{(2n)}_{j},x^{(2n)}_{k}\big{)}\sqrt{w^{(2n)}_{j}w^{(2n)}_{k}}\end{bmatrix}^{n}_{j,k=1},$

where

$\displaystyle H_{jk}=\frac{{\rm
e}^{\frac{\varepsilon^{2}}{2}+s\varepsilon}\sin{(\pi\gamma)}}{2\pi^{2}}\big{|}x^{(2n)}_{j}x^{(2n)}_{k}\big{|}^{\frac{\gamma}{2}}{\rm
e}^{-\frac{s}{2}\big{(}x^{(2n)}_{j}+x^{(2n)}_{k}\big{)}}$
$\displaystyle\hphantom{H_{jk}=}{}\times\sqrt{w^{(2n)}_{j}w^{(2n)}_{k}}\int_{\mathbb{R}}\mathrm{d}t\frac{{\rm
e}^{-\frac{t^{2}}{2}}({\rm i}t+\varepsilon)^{-\gamma}{\rm e}^{{\rm
i}t(s+\varepsilon)}}{\big{(}x^{(2n)}_{j}-\varepsilon-{\rm
i}t\big{)}\big{(}x^{(2n)}_{k}-\varepsilon-{\rm i}t\big{)}}.$

In order to compute the above integral, we once again apply the Gauss–Hermite
quadrature

$\displaystyle\int_{\mathbb{R}}\mathrm{d}t\frac{{\rm
e}^{-\frac{t^{2}}{2}}({\rm i}t+\varepsilon)^{-\gamma}{\rm e}^{{\rm
i}t(s+\varepsilon)}}{\big{(}x^{(2n)}_{j}-\varepsilon-{\rm
i}t\big{)}\big{(}x^{(2n)}_{k}-\varepsilon-{\rm
i}t\big{)}}=\sum^{2n}_{\ell=1}\frac{w^{(2n)}_{\ell}\big{(}{\rm
i}x^{(2n)}_{\ell}+\varepsilon\big{)}^{-\gamma}{\rm e}^{{\rm
i}x^{(2n)}_{\ell}(s+\varepsilon)}}{\big{(}x^{(2n)}_{j}-\varepsilon-{\rm
i}x^{(2n)}_{\ell}\big{)}\big{(}x^{(2n)}_{k}-\varepsilon-{\rm
i}x^{(2n)}_{\ell}\big{)}}.$

Note that, in the above, we use all the $2n$ nodes because the integral is
over the whole ${\mathbb{R}}$. In principle, we could use a different number
of nodes for this last quadrature, but this is the way the code was actually
implemented for simplicity.

The matrix $\mathbb{H}=\left[H_{jk}\right]_{j,k=1}^{n}$ can now be written as

$\displaystyle\mathbb{H}=\begin{bmatrix}H_{jk}\end{bmatrix}^{n}_{j,k=1}=\frac{{\rm
e}^{\frac{\varepsilon^{2}}{2}+s\varepsilon}}{2\pi^{2}}\sin{(\pi\gamma)}\mathbb{A}.\mathbb{A}^{\rm
T},$

where $\mathbb{A}$ is the $n\times 2n$ matrix given by

$\displaystyle\mathbb{A}_{j\ell}=\left|x_{j}\right|^{\frac{\gamma}{2}}\sqrt{w_{j}}{\rm
e}^{-\frac{s}{2}x_{j}}\frac{({\rm
i}x_{\ell}+\varepsilon)^{-\frac{\gamma}{2}}{\rm e}^{{\rm
i}\frac{x_{\ell}}{2}(s+\varepsilon)}}{\left(x_{j}-\varepsilon-{\rm
i}x_{\ell}\right)}\sqrt{w_{\ell}},$

and $1\leq j\leq n$ and $1\leq\ell\leq 2n$.

The result of this approximation is expressed in terms of the parameters
$s\in{\mathbb{R}}$, $\gamma\in\left[0,1\right)$ and $n\in{\mathbb{N}}$, which
is the number of nodes in the Gaussian quadrature. Therefore, the final
approximation that we have for $\tau$ is

$\displaystyle\tau(s,\gamma,n)=\det{\left(\mathbf{1}_{n\times n}-\frac{{\rm
e}^{\frac{\varepsilon^{2}}{2}+s\varepsilon}}{2\pi^{2}}\sin{(\pi\gamma)}\mathbb{A}.\mathbb{A}^{\rm
T}\right)}.$

We take $\varepsilon=\max{[0,-s]}$ so as to optimize the distance from the
saddle-point of the function $\theta(x)=x^{2}/4+sx/2$ while respecting the
condition that the relative interior of the contours ${\mathbb{R}}_{-}$ and
$\Gamma_{\infty}$ do not intersect. It is worth pointing out that the value of
the Fredholm determinant (numerical considerations aside) is independent of
$\varepsilon$ thanks to Cauchy’s theorem. Various numerical results are
pictures in Figs. 6, 7, 8.

(a) $\gamma=0.1$.

(b) $\gamma=0.5$.

Figure 6: $\tau$-function with $n=30$ points.

(a) $\gamma=0.1$.

(b) $\gamma=0.5$.

Figure 7: $\tau$-function with $n=80$ points. Figure 8: Arctan of the
$\tau$-function for $\gamma=0.1$ and with $n=150$ with an extended range of
values for $s$ being considered.

##  3 Orthogonal polynomials and their Riemann–Hilbert problem

The orthogonality relations for the polynomials $p_{n}$ can be rewritten

$\displaystyle\int_{\mathbb{C}}p_{n}(\lambda)\bar{\lambda}^{jd+\ell}{\rm
e}^{-NW(\lambda)}\mathrm{d}A(\lambda)=0,\qquad j=0,\dots,k-1,$ $\displaystyle
n=kd+\ell,\qquad 0\leq\ell\leq d-1.$

As a consequence, the $n$-th monic orthogonal polynomial has a discrete
symmetry

$\displaystyle p_{n}\big{(}{\rm e}^{\frac{2\pi{\rm i}}{d}}\lambda\big{)}={\rm
e}^{\frac{2\pi{\rm i}n}{d}}p_{n}(\lambda)\quad\Rightarrow\quad
p_{n}(\lambda)=\lambda^{\ell}q^{(\ell)}_{k}\big{(}\lambda^{d}\big{)}.$

Using this definition, the initial sequence of orthogonal polynomials
$\\{p_{n}(\lambda)\\}^{\infty}_{n=0}$ can be split in $d$ sub-sequences, each
of which labelled by the remainder $\ell\equiv n$ mod $d$. Through a change of
coordinate $\lambda^{d}=u$ these sequences of monic orthogonal polynomials are
seen to satisfy the orthogonality relations

$\displaystyle\int_{\mathbb{C}}q_{k}^{(\ell)}(u)\bar{u}^{j}|u|^{-2\gamma}{\rm
e}^{-N(|u|^{2}-tu-t\bar{u})}\mathrm{d}A(u)=0,\qquad j=0,\dots,k-1,$ (3.1)
$\displaystyle\gamma:=\frac{d-\ell-1}{d}\in[0,1).$

With the further change of coordinate $u=-t(z-1)$, $z\in\mathbb{C}$, we
introduce the transformed monic polynomial

$\displaystyle\pi_{k}(z):=\frac{(-1)^{k}}{t^{k}}q^{(\ell)}_{k}(-t(z-1)),$
(3.2)

Figure 9: The contour $\Sigma$.

Next we transform the orthogonality on the plane to the orthogonality on a
contour as in [4] and [3]

######  Theorem 3.1 ([4, Theorem 2.1]).

For any polynomial $q(u)$ the following identity is satisfied

$\displaystyle\int_{\mathbb{C}}\\!q(u)\bar{u}^{j}|u|^{-2\gamma}{\rm
e}^{-N(|u|^{2}-tu-t\bar{u})}\mathrm{d}A(u)=\frac{\pi\Gamma(j-\gamma+1)}{N^{j-\gamma+1}}\frac{1}{2\pi{\rm
i}}\oint_{\tilde{\Sigma}}\\!q(u)\frac{{\rm
e}^{Ntu}}{(u-t)^{j+1}}\left(1-\frac{t}{u}\right)^{\gamma}\\!\mathrm{d}u,$

where $\gamma\in(0,1)$, $j$ is an arbitrary non-negative integer, and
$\tilde{\Sigma}$ is a positively oriented simple closed loop enclosing $u=0$
and $u=t$.

Consider now $q^{(\ell)}_{k}(u)$ to be the monic polynomial of degree $k$ with
the orthogonality relations (3.1). Then according to Theorem 3.1 the
transformed monic polynomial $\pi_{k}(z)$ of degree $k$ given by (3.2), is
characterised by the non-hermitian orthogonality relations

$\displaystyle\oint_{\Sigma}\pi_{k}(z)z^{j}\frac{{\rm
e}^{-Nt^{2}z}}{z^{k}}\left(\frac{z}{z-1}\right)^{\gamma}\mathrm{d}z=0,\qquad
j=0,1,\dots,k-1,$ (3.3)

where $\Sigma$ is a simple, positively oriented contour encircling $z=0$ and
$z=1$, as it can be seen in Fig. 9, and the function
$\big{(}\frac{z}{z-1}\big{)}^{\gamma}$ is analytic in
$\mathbb{C}\setminus[0,1]$ and tends to one for $|z|\to\infty$.

###  3.1 The Riemann–Hilbert problem

We now consider the polynomials $\pi_{k}(z)$ in the limit $k\rightarrow\infty$
and $N\rightarrow\infty$ in such a way that, for $n=kd+\ell$, one has

$\displaystyle T=\frac{n-\ell}{N}>0.$

We set the notation

$\displaystyle V(z)=\frac{z}{z_{0}}+\log{z},\qquad
z_{0}=z_{0}(t,N)=\frac{k}{N}\frac{1}{t^{2}},\qquad
z_{0}=\frac{t^{2}_{c}}{t^{2}},\qquad t^{2}_{c}=\frac{T}{d},$ $\displaystyle
w_{k}(z):={\rm e}^{-kV(z)}\left(\frac{z}{z-1}\right)^{\gamma}.$

The orthogonality relations (3.3) now read

$\displaystyle\oint_{\Sigma}\pi_{k}(z)z^{j}w_{k}(z)\mathrm{d}z=0,\qquad
j=0,1,\dots,k-1.$

When the limit $k\rightarrow\infty$ is taken, three different regimes arise

  * •

pre-critical case: $0<t<t_{c}$, leading to $z_{0}>1$,

  * •

critical case: $t=t_{c}$, leading to $z_{0}=1$,

  * •

post-critical case: $t>t_{c}$, leading to $z_{0}<1$.

The pre- and post-critical case were already analyzed in [4] and we are now
interested in analyzing the critical case. We begin by defining the so-called
complex moments as

$\displaystyle\nu_{j}:=\oint_{\Sigma}z^{j}w_{k}(z)\mathrm{d}z,$

where the dependency on $k$ is omitted in order to simplify notation, and use
this to introduce the auxiliary polynomial

$\displaystyle\Pi_{k-1}(z):=\frac{1}{\det[\nu_{i+j}]_{0\leq i,j\leq
k-1}}\det\left[\begin{matrix}\nu_{0}&\nu_{1}&\dots&\nu_{k-1}\\\
\nu_{1}&\nu_{2}&\dots&\nu_{k}\\\ \vdots&&&\vdots\\\
\nu_{k-2}&\dots&&\nu_{2k-3}\\\ 1&z&\dots&z^{k-1}\end{matrix}\right].$

It can be seen that this polynomial is not necessarily monic and its degree
may be less than $k-1$. In order to guarantee the existence of such a
polynomial, the determinant in the denominator must not vanish.

######  Proposition 3.2 ([4, Proposition 2.2]).

The determinant $\det[\nu_{i+j}]_{0\leq i,j\leq k-1}$ does not vanish and is
given by

$\displaystyle\det[\nu_{i+j}]_{0\leq i,j\leq
k-1}=(-1)^{k(k-1)/2}(2\mathrm{i})^{k}\left(\prod_{j=0}^{k-1}t^{2j+2\gamma-2}\frac{N^{j-\gamma+1}}{\Gamma(j-\gamma+1)}\right)$
$\displaystyle\hphantom{\det[\nu_{i+j}]_{0\leq i,j\leq
k-1}=}{}\times\det\left[\iint_{\mathbb{C}}z^{i}\bar{z}^{j}|z-1|^{-2\gamma}{\rm
e}^{-Nt^{2}|z|^{2}}\mathrm{d}A(z)\right].$

We can now reformulate the condition of orthogonality for the polynomials
$\pi_{k}(z)$ as a Riemann–Hilbert boundary value problem. To do this, we
define the matrix

$\displaystyle Y(z)=\begin{bmatrix}\pi_{k}(z)&\dfrac{1}{2\pi{\rm
i}}\displaystyle{\int_{\Sigma}}\dfrac{\pi_{k}(z^{\prime})}{z^{\prime}-z}w_{k}(z^{\prime})\mathrm{d}z^{\prime}\\\
-{2\pi{\rm
i}}\Pi_{k-1}(z)&-\displaystyle{\int_{\Sigma}}\dfrac{\Pi_{k-1}(z^{\prime})}{z^{\prime}-z}w_{k}(z^{\prime})\mathrm{d}z^{\prime}\end{bmatrix},$

which is the unique solution of following Riemann–Hilbert problem for
orthogonal polynomials [22].

Figure 10: Jump in the $Y(z)$ Riemann–Hilbert problem.

######  Riemann–Hilbert Problem 3.3.

The matrix $Y(z)$ is analytic in $\mathbb{C}\setminus\Sigma$, where $\Sigma$
is the oriented curve in Fig. 10. The limits $Y_{\pm}(z)$ exist and are
continuous along $\Sigma$. Moreover

  1. $1.$

Jump on $\Sigma$: The continuous boundary values $Y_{\pm}(z)$ are such that

$\displaystyle Y_{+}(z)=Y_{-}(z)\begin{pmatrix}1&w_{k}(z)\\\
0&1\end{pmatrix},\qquad z\in\Sigma.$

This jump can be seen in Fig. 10.

  2. $2.$

Behaviour at infinity: $Y(z)$ has the following behaviour as
$z\rightarrow\infty$

$\displaystyle
Y(z)=\left(\mathbf{1}+\mathcal{O}\left(\frac{1}{z}\right)\right)z^{k\sigma_{3}}.$

##  4 Riemann–Hilbert analysis

###  4.1 Transforming the Riemann–Hilbert problem

We will now begin by doing a simple transformation through the first
undressing step and, after that, use the steepest descent method of Deift–Zhou
[16] to proceed with further transformations and simplifications of our
problem.

####  4.1.1 First undressing step

The first undressing step consists in a simplification of the Riemann–Hilbert
Problem 3.3 that is done by defining a new matrix-valued function
$\tilde{Y}(z)$ as

$\displaystyle\tilde{Y}(z):=Y(z)\left(1-\frac{1}{z}\right)^{-\frac{\gamma}{2}\sigma_{3}},\qquad
z\in\mathbb{\setminus}(\Sigma\cup[0,1]).$

Clearly the matrix $\tilde{Y}$ now satisfies a new Riemann–Hilbert problem
whose simple formulation we leave to the reader. From the solution of
Riemann–Hilbert problem we can recover the orthogonal polynomials through

$\displaystyle\pi_{k}(z)=\tilde{Y}_{11}(z)\left(1-\frac{1}{z}\right)^{\frac{\gamma}{2}}.$

####  4.1.2 Transformation $\boldsymbol{\tilde{Y}\rightarrow{R}}$

We now define the so-called $g$-function; this is a scalar function that is
analytic off a contour $\Gamma_{r}$ (see Fig. 13), homotopically equivalent to
$\Sigma$ in $\mathbb{C}\setminus[0,1]$ and suitably chosen. Following [4] let

$\displaystyle
g(z)=\begin{cases}\dfrac{z}{z_{0}}+\ell,&z\in\operatornamewithlimits{Int}(\Gamma_{r}),\\\
\log{z},&z\in\operatornamewithlimits{Ext}(\Gamma_{r}),\end{cases}$ (4.1)
$\displaystyle\ell:=\log{r}-\dfrac{r}{z_{0}},\qquad r>0,$

where the logarithm is the principal determination and

$\displaystyle\Gamma_{r}=\big{\\{}z\in\mathbb{C}\colon\operatorname{Re}(\varphi(z;r))=0,\,|z|\leq
z_{0}\big{\\}},\qquad 0<r\leq z_{0},$ (4.2)

and

$\displaystyle\varphi(z;r)=\begin{cases}\log{z}-\dfrac{z}{z_{0}}-\log{r}+\dfrac{r}{z_{0}},&z\in\operatornamewithlimits{Int}(\Gamma_{r}),\vspace{1mm}\\\
\dfrac{z}{z_{0}}-\log{z}+\log{r}-\dfrac{r}{z_{0}},&z\in\operatornamewithlimits{Ext}(\Gamma_{r}),\end{cases}$

or equivalently,

$\displaystyle\varphi(z;r)=\log{r}-\dfrac{r}{z_{0}}+V-2g(z).$ (4.3)

Note that $g$ is such that

$\displaystyle g_{+}(z)+g_{-}(z)-\ell-V(z)\equiv
0,\qquad\forall\,z\in\Gamma_{r}.$ (4.4)

In order to establish the transformation $\tilde{Y}\rightarrow{R}$ we will now
deform the contour $\Sigma$ to the contour $\Gamma_{r}$. The choice of $r$
will be explained later. For the time being we define the matrix ${R}$ as

$\displaystyle{R}(z)={\rm e}^{-k\frac{\ell}{2}\sigma_{3}}\tilde{Y}(z){\rm
e}^{-kg(z)\sigma_{3}}{\rm e}^{k\frac{\ell}{2}\sigma_{3}},\qquad
z\in\mathbb{C}\setminus(\Gamma_{r}\cup[0,1]).$

The matrix ${R}(z)$ solves the following Riemann–Hilbert problem

######  Riemann–Hilbert Problem 4.1.

The matrix ${R}(z)$ is analytic in $\mathbb{C}\setminus(\Gamma_{r}\cup[0,1])$
and admits non-tangential boundary values that satisfy:

  1. $1.$

Jumps on $\Gamma_{r}$ and $[0,1]$ $($see Fig. 11):

$\displaystyle{R}_{+}(z)={R}_{-}(z)\begin{cases}\begin{pmatrix}{\rm
e}^{-k(g_{+}-g_{-})}&{\rm e}^{k(g_{+}+g_{-}-\ell-V(z))}\\\ 0&{\rm
e}^{k(g_{+}-g_{-})}\end{pmatrix},&z\in\Gamma,\\\ {\rm e}^{-\gamma\pi{\rm
i}\sigma_{3}},&z\in(0,1).\end{cases}$ (4.5)

Figure 11: Jumps for the Riemann–Hilbert problem of the matrix ${R}(z)$.

  2. $2.$

Large $z$ boundary behaviour:

$\displaystyle{R}(z)=\left(\mathbf{1}+\mathcal{O}\left(\frac{1}{z}\right)\right),\qquad
z\rightarrow\infty.$

  3. $3.$

Endpoint behaviour:

$\displaystyle{R}(z)=\mathcal{O}(1)z^{\frac{\gamma}{2}\sigma_{3}},\qquad
z\rightarrow 0,$
$\displaystyle{R}(z)=\mathcal{O}(1)(z-1)^{-\frac{\gamma}{2}\sigma_{3}},\qquad
z\to 1.$

The orhtogonal polynomials can be recovered from this Riemann–Hilbert problem
through

$\displaystyle\pi_{k}(z)={R}_{11}(z){\rm
e}^{kg(z)}\left(1-\frac{1}{z}\right)^{\frac{\gamma}{2}}.$ (4.6)

####  4.1.3 Transformation $\boldsymbol{{R}\rightarrow T}$

The jump matrix of ${R}(z)$ on $z\in\Gamma$ given by (4.5) can be factorized
in the following way

$\displaystyle\begin{pmatrix}{\rm e}^{-k\left(g_{+}-g_{-}\right)}&{\rm
e}^{k\left(g_{+}+g_{-}-\ell-V\right)}\\\ 0&{\rm
e}^{k\left(g_{+}-g_{-}\right)}\end{pmatrix}$
$\displaystyle\qquad{}=\begin{pmatrix}1&0\\\ {\rm
e}^{k\left(\ell+V-2g_{-}\right)}&1\end{pmatrix}\begin{pmatrix}0&{\rm
e}^{k\left(g_{+}+g_{-}-\ell-V\right)}\\\ -{\rm
e}^{k\left(g_{+}+g_{-}-\ell-V\right)}&0\end{pmatrix}\begin{pmatrix}1&0\\\ {\rm
e}^{k\left(\ell+V-2g_{+}\right)}&1\end{pmatrix}$
$\displaystyle\qquad{}=\begin{pmatrix}1&0\\\ {\rm
e}^{k\varphi(z)}&1\end{pmatrix}\begin{pmatrix}0&1\\\
-1&0\end{pmatrix}\begin{pmatrix}1&0\\\ {\rm e}^{k\varphi(z)}&1\end{pmatrix},$

where, in order to write the second line, we have used the definition of the
$\varphi(z)$ function defined in (4.3) and the function $g(z)$ defined in
(4.1) as well as the relation (4.4).

We will now consider three different loops $\Gamma_{i}$, $\Gamma_{r}$ and
$\Gamma_{e}$, so that the space is split into four different domains
$\Omega_{0}$, $\Omega_{1}$, $\Omega_{2}$ and $\Omega_{\infty}$, as it can be
seen in Fig. 12. $\Gamma_{i}$ is in the interior of $\Gamma_{r}$ and
$\Gamma_{e}$ is in the exterior. Using this, we will now define a new matrix-
valued function $T(z)$ in the following way

$\displaystyle T(z)=\begin{cases}{R}(z),&z\in\Omega_{0}\cup\Omega_{\infty},\\\
{R}(z)\begin{pmatrix}1&0\\\ -{\rm
e}^{k\varphi(z)}&1\end{pmatrix},&z\in\Omega_{1},\vspace{1mm}\\\
{R}(z)\begin{pmatrix}1&0\\\ {\rm
e}^{k\varphi(z)}&1\end{pmatrix},&z\in\Omega_{2}.\end{cases}$ (4.7)

This matrix $T(z)$ satisfies the following Riemann–Hilbert problem

Figure 12: Jumps in the $T(z)$ Riemann–Hilbert problem.

######  Riemann–Hilbert Problem 4.2.

The matrix $T(z)$ is analytic in
$\mathbb{C}\setminus(\Gamma_{i}\cup\Gamma_{r}\cup\Gamma_{e}\cup[0,1])$ and
admits non-tangential boundary values. Moreover

  1. $1.$

On $\Sigma_{T}=\Gamma_{i}\cup\Gamma_{r}\cup\Gamma_{e}\cup[0,1]$ the boundary
values satisfy

$\displaystyle T_{+}(z)=T_{-}(z)v_{T},\qquad z\in\Sigma_{T},$

where

$\displaystyle v_{T}=\begin{cases}\begin{pmatrix}1&0\\\ {\rm
e}^{k\varphi(z)}&1\end{pmatrix},&z\in\Gamma_{i},\vspace{1mm}\\\
\begin{pmatrix}0&1\\\ -1&0\end{pmatrix},&z\in\Gamma_{r},\vspace{1mm}\\\
\begin{pmatrix}1&0\\\ {\rm e}^{k\varphi(z)}&1\end{pmatrix},&z\in\Gamma_{e},\\\
{\rm e}^{-\gamma\pi{\rm i}\sigma_{3}},&z\in(0,1).\end{cases}$ (4.8)

  2. $2.$

Large $z$ boundary behaviour:

$\displaystyle
T(z)=\left(\mathbf{1}+\mathcal{O}\left(\frac{1}{z}\right)\right),\qquad
z\rightarrow\infty.$

  3. $3.$

Endpoint behaviour:

$\displaystyle T(z)=\mathcal{O}(1)z^{\frac{\gamma}{2}\sigma_{3}},\qquad
z\rightarrow 0,$ $\displaystyle
T(z)=\mathcal{O}(1)(z-1)^{-\frac{\gamma}{2}\sigma_{3}},\qquad z\to 1.$

###  4.2 Choice of contour

Since our analysis will be performed in the critical case when $z_{0}\to 1$,
the choice of the contour $\Gamma_{r}$ with $0<r\leq z_{0}=1$ is forced (Fig.
13).

Figure 13: The contours $\Gamma_{r}$ for $z_{0}=1$ and different values of
$0<r\leq 1$. Figure 14: Contours in the critical case $z_{0}=1$.

The corresponding contours $\Gamma_{i}$ and $\Gamma_{e}$ have to be chosen in
such a way that the jump matrix $v_{T}$ defined in (4.8) converges (as
$k\to\infty$ and hence also as $n\to\infty$) exponentially fast to a constant.
We make the choice of the deformed contours $\Gamma_{i}$ and $\Gamma_{e}$ as
in Fig. 14. This is the correct choice because the sign of
$\operatorname{Re}\varphi(z)$ on the contours $\Gamma_{i}$ and $\Gamma_{e}$ is
negative (see Fig. 15 where the region $\operatorname{Re}\varphi(z)<0$ is
plotted).

Therefore, the jump matrices on $\Gamma_{i}$ and $\Gamma_{e}$ converge
exponentially fast to the identity matrix in any $L^{p}$ norm except
$p=\infty$ because $\varphi(z=1)=0$.

Figure 15: Contours $\Gamma$ chosen for the critical case $z_{0}=r=1$. The
orange region corresponds to the part of the complex plane where
$\operatorname{Re}{(\varphi(z))}<0$ and the green contour is the region where
$\operatorname{Re}{(\varphi(z))}>0$.

With this choice of contours we arrive to the following theorem (whose proof
is a simple inspection and hence omitted)

######  Theorem 4.3.

The jump matrix $v_{T}$ converges exponentially fast as $k\to\infty$ to the
constant jump matrix

$\displaystyle v_{\infty}=\begin{cases}\begin{pmatrix}1&0\\\
0&1\end{pmatrix},&z\in(\Gamma_{e}\cup\Gamma_{i})\setminus\mathbb{D},\vspace{1mm}\\\
\begin{pmatrix}0&1\\\ -1&0\end{pmatrix},&z\in\Gamma_{1},\\\ {\rm
e}^{-\gamma\pi{\rm i}\sigma_{3}},&z\in(0,1),\end{cases}$ (4.9)

where $\mathbb{D}$ is a small disk surrounding the point $z=1$.

###  4.3 Approximate solutions to $\boldsymbol{T(z)}$

We are now ready to approximate the matrix $T(z)$ with two solutions, one
outside a neighbourhood of $z=1$ and one inside. We call exterior parametrix
the Riemann–Hilbert problem solved by the matrix $M(z)$ with jump
$v_{\infty}(z)$. We call local parametrix the solution $P(z)$ of the
Riemann–Hilbert problem obtained within a neighbourhood of $z=1$. These two
solutions are approximations of the exact solution, $T(z)$, in the limit
$k\to\infty$. In order to obtain the asymptotics of the orthogonal polynomials
$\pi_{k}(z)$, we need sub-leading corrections to the matrices ${N}(z)$ and
$P(z)$ and this will be accomplished by evaluating perturbatively the error
matrix $E(z)$, which is defined as

$\displaystyle
E(z)=\begin{cases}T(z){N}(z)^{-1},&z\in\mathbb{C}\setminus\mathbb{D},\\\
T(z)P(z)^{-1},&z\in\mathbb{D}.\end{cases}$

We will first construct the matrix ${N}(z)$ and then the matrix $P(z)$.

####  4.3.1 Exterior parametrix

The exterior parametrix, ${N}(z)$, is the piecewise analytic $2\times 2$
matrix

$\displaystyle{N}(z)=\left(1-\frac{1}{z}\right)^{\frac{\gamma}{2}\sigma_{3}}\begin{pmatrix}0&1\\\
-1&0\end{pmatrix}^{\chi_{L}},$

where $\chi_{L}$ is the characteristic function that is one on the left of the
contour $\hat{\Gamma}_{r=1}$ and is zero otherwise. By design, it satisfies
${N}(\infty)=\mathbf{1}$ and

$\displaystyle{N}_{+}(z)={N}_{-}(z)v_{\infty},\qquad
z\in{\Gamma_{1}\cup[0,1]},$

where $v_{\infty}$ has been defined in (4.9).

####  4.3.2 Local parametrix and double scaling limit

The local parametrix near the point $z=1$ is the solution, $P(z)$, of a matrix
Riemann–Hilbert problem that has the same jumps as $T(z)$ and that matches
${N}(z)$ on the boundary of a disc centered at $z=1$.

######  Riemann–Hilbert Problem 4.4.

The matrix $P(z)$ is analytic in $\mathbb{D}\setminus\Sigma_{P}$,
$\Sigma_{P}=\Gamma_{e}\cup\Gamma_{r=1}\cup\Gamma_{i}\cup[0,1]$ and admits non-
tangential boundary values. Moreover

  1. $1.$

Jumps on $\Sigma_{P}$:

$\displaystyle P_{+}(z)=P_{-}(z)v_{P}(z),\qquad z\in\Sigma_{P}\cap\mathbb{D},$

where

$\displaystyle v_{P}(z)=\begin{cases}\begin{pmatrix}1&0\\\ {\rm
e}^{k\varphi(z)}&1\end{pmatrix},&z\in\Gamma_{i}\cap\mathbb{D},\vspace{1mm}\\\
\begin{pmatrix}0&1\\\
-1&0\end{pmatrix},&z\in\Gamma_{r}\cap\mathbb{D},\vspace{1mm}\\\
\begin{pmatrix}1&0\\\ {\rm
e}^{k\varphi(z)}&1\end{pmatrix},&z\in\Gamma_{e}\cap\mathbb{D},\\\ {\rm
e}^{-\gamma\pi{\rm i}\sigma_{3}},&z\in(0,1)\cap\mathbb{D},\end{cases}$

which are the jumps that can be seen in Fig. 16.

Figure 16: Jumps of the Riemann–Hilbert problem for the matrix $P(z)$.

  2. $2.$

Behaviour at the boundary $\partial\mathbb{D}$:

$\displaystyle P(z)={N}(z)(\mathbf{1}+o(1)),\qquad\text{as}\quad
k\to\infty\quad\text{and}\quad z\in\partial\mathbb{D}.$ (4.10)

####  4.3.3 The double scaling limit

Let $\phi(z;z_{0})$ be defined by

$\displaystyle\phi(z;z_{0})=\left(\frac{z-1}{z_{0}}-\log{z}\right).$

We note that $\phi(z;z_{0})$ is the analytic continuation of
$\varphi(z;z_{0})$ (4.3) (setting $r=1$) from the exterior of
$\Gamma=\Gamma_{1}$. For $z_{0}\sim 1$ we see that it has a single critical
point in the neighbourhood of $z=1$; then

######  Proposition 4.5.

There exists a jointly analytic function $\zeta(z;z_{0})$ which is univalent
in a fixed neighbourhood of $z=1$, with $z_{0}$ in a neighbourhood of
$z_{0}=1$, and an analytic function $A(z_{0})$ near $z_{0}=1$ such that

$\displaystyle\phi(z;z_{0})=\frac{1}{2}\zeta^{2}(z;z_{0})+A(z_{0})\zeta(z;z_{0}),\qquad\zeta(1;z_{0})\equiv
0,\qquad A(1)=0.$ (4.11)

###### Proof.

The function $\phi(z;z_{0})$ has a critical point at $z=z_{0}$ with critical
value $\phi_{\rm cr}=\frac{z_{0}-1}{z_{0}}-\ln z_{0}$. Since we are interested
in the values $z_{0}\simeq 1$ we observe that $\phi_{\rm
cr}=-\frac{1}{2}(z_{0}-1)^{2}(1+\mathcal{O}(z_{0}-1))$. Hence we can write

$\displaystyle\phi_{\rm cr}=-\frac{A(z_{0})^{2}}{2}$

with $A(z_{0})$ an analytic series near $z_{0}=1$. Moreover

$\displaystyle\phi(z;z_{0})-\phi_{\rm
cr}=\frac{1}{2z_{0}^{2}}(z-z_{0})^{2}(1+\mathcal{O}(z-z_{0})).$

Define

$\displaystyle\frac{\zeta(z;z_{0})}{\sqrt{2}}:=\sqrt{\phi(z;z_{0})+\frac{A^{2}(z_{0})}{2}}-\frac{A(z_{0})}{\sqrt{2}}.$

This is a conformal map in a neighbourhood of $z=z_{0}$ of uniform radius of
convergence for $z_{0}\simeq 1$; moreover, from $\phi(1;z_{0})=0$ follows that
$\zeta(1;z_{0})\equiv 0$. The equation (4.11) is trivially satisfied. ∎

######  Remark 4.6.

A simple computation yields the following first few terms for the expansion
($\delta=z_{0}-1$)

$\displaystyle\zeta=\left(1+\frac{\delta}{3}+\frac{\delta^{2}}{4}\right)(z-1)-\frac{(z-1)^{2}}{3}\left(1-\frac{\delta}{4}+\frac{11\delta^{2}}{60}\right)+\mathcal{O}\big{(}(z-1)^{3},\delta^{3}\big{)},$
$\displaystyle
A(\delta)=-\delta+\frac{2}{3}\delta^{3}+\mathcal{O}\big{(}\delta^{3}\big{)}.$

From this it follows

$\displaystyle\sqrt{k}\zeta(z;z_{0})=\left(z-1-\frac{(z-1)^{2}}{3}\right)\sqrt{k}+\mathcal{O}(1),$

which manifests the fact that $\sqrt{k}\zeta$ homothetically expands the image
of the disk $\mathbb{D}$.

######  Definition 4.7 (double-scaling limit).

The double-scaling limit is defined by taking $k\rightarrow\infty$ and
$A\rightarrow 0$ (or equivalently $z_{0}\rightarrow 1$) so that

$\displaystyle\lim_{k\to\infty,\,z_{0}\to 1}\sqrt{k}A=\mathcal{S}$

with $\mathcal{S}$ in compact subsets of the complex plane. Equivalently, it
can be defined as

$\displaystyle A\sim\frac{\mathcal{S}}{\sqrt{k}}\qquad\hbox{or}\qquad
z_{0}\sim{1-\frac{\mathcal{S}}{\sqrt{k}}}.$ (4.12)

####  4.3.4 Model problem and the Painlevé IV equation

We will now recall the Riemann–Hilbert Problem 2.1 for the function
$\Psi(\lambda;s)$ associated to the Painlevé IV equation. For convenience, we
replace the contour $\Gamma_{1}$ with $\hat{\Gamma}_{i}$ and the contour
$\Gamma_{\infty}$ with $\hat{\Gamma}_{e}$ in order to be consistent with the
notation we are using in this section.

Given that $\Theta_{\infty}=\Theta_{0}=\frac{\gamma}{2}$, the equations (2.2)
and (2.7) reduce to a somewhat simpler form

$\displaystyle U^{\prime}=U(Y-s),\qquad
Z=\dfrac{1}{2}\big{(}sY+\gamma-Y^{\prime}-Y^{2}\big{)},\qquad
H=\left(s+\dfrac{\gamma}{Y}-Y\right)Z-\dfrac{Z^{2}}{Y},$ (4.13)

where $Y=Y(s)$ solves the Painlevé IV equation (2.3) in the form

$\displaystyle
Y^{\prime\prime}=\frac{1}{2}\frac{(Y^{\prime})^{2}}{Y}+\frac{3}{2}Y^{3}-2sY^{2}+\left(1+\frac{s^{2}}{2}-\gamma\right)Y-\frac{\gamma^{2}}{2Y}.$

For our purposes, we need to modify the Riemann–Hilbert problem for
$\Psi(\lambda)$ to be able to match the Riemann–Hilbert problem for $P(z)$.
Let us introduce the contour $\hat{\Gamma}_{r=1}$ on the $\lambda$ plane,
which is a contour between $\hat{\Gamma}_{i}$ and $\hat{\Gamma}_{e}$ and
passing through $\lambda=0$. Then, let us define

$\displaystyle\widehat{\Psi}(\lambda)=\Psi(\lambda){\rm
e}^{\theta\sigma_{3}}\begin{pmatrix}0&1\\\ -1&0\end{pmatrix}^{\chi_{L}},$

where $\theta$ is given by (2.6) and $\chi_{L}$ is the characteristic function
that is one on the left of the contour $\hat{\Gamma}_{r=1}$ and is zero
otherwise. It is now straightforward to check that the Riemann–Hilbert problem
for $\widehat{\Psi}(\lambda)$ is given by

######  Riemann–Hilbert Problem 4.8.

The matrix $\widehat{\Psi}(\lambda)\textrm{is analytic in
}\mathbb{C}\setminus\Sigma_{\widehat{\Psi}}$,
$\Sigma_{\widehat{\Psi}}=\hat{\Gamma}_{e}\cup\hat{\Gamma}_{i}\cup\hat{\Gamma}_{r=1}\cup\mathbb{R}^{-}$
and admits non-tangential boundary values. Moreover

  1. $1.$

Jumps on $\Sigma_{\widehat{\Psi}}$:

$\displaystyle\widehat{\Psi}_{+}(\lambda)=\widehat{\Psi}_{-}(\lambda)v_{\widehat{\Psi}}(\lambda),\qquad\zeta\in\Sigma_{\widehat{\Psi}},$

where

$\displaystyle v_{\widehat{\Psi}}(\lambda)=\begin{cases}\begin{pmatrix}1&0\\\
{\rm
e}^{-2\theta(\lambda)}&1\end{pmatrix},&z\in\hat{\Gamma}_{i},\vspace{1mm}\\\
\begin{pmatrix}0&1\\\
-1&0\end{pmatrix},&z\in\hat{\Gamma}_{r=1},\vspace{1mm}\\\
\begin{pmatrix}1&0\\\ {\rm
e}^{2\theta(\lambda)}&1\end{pmatrix},&z\in\hat{\Gamma}_{e},\\\ {\rm
e}^{-\gamma\pi{\rm i}\sigma_{3}},&z\in\mathbb{R}^{-}.\end{cases}$ (4.14)

  2. $2.$

Behaviour for $\lambda\to\infty$

$\displaystyle\widehat{\Psi}(\lambda)=\left(\mathbf{1}+\dfrac{\Psi_{1}}{\lambda}+\dfrac{\Psi_{2}}{\lambda^{2}}+\mathcal{O}\big{(}\lambda^{-3}\big{)}\right)\lambda^{\frac{\gamma}{2}\sigma_{3}}\begin{pmatrix}0&1\\\
-1&0\end{pmatrix}^{\chi_{L}}.$ (4.15)

Comparing the jump matrices for $P(z)$ and $\widehat{\Psi}(\lambda)$, we are
now ready to obtain the local parametrix $P(z)$ defined by the Riemann–Hilbert
Problem 4.4, which is given by

$\displaystyle P(z)={N}(z)\begin{pmatrix}0&-1\\\
1&0\end{pmatrix}^{\chi_{L}}\big{(}\sqrt{k}\zeta(z)\big{)}^{-\frac{\gamma}{2}\sigma_{3}}\widehat{\Psi}\big{(}\sqrt{k}\zeta(z);\sqrt{k}A\big{)},$
(4.16)

where $\zeta(z)=\zeta(z;z_{0})$ as defined in (4.5). We observe that the
product of the first three terms of (4.16) is holomorphic in the neighbourhood
of $z=1$ and therefore does not change the Riemann–Hilbert problem.
Furthermore, the first three terms have been inserted in order to have the
behaviour (4.10) for $z\in\partial\mathbb{D}$ in the limit $k\to\infty$ and
$A\to 0$ in such a way that

$\displaystyle\lim_{k\to\infty}\sqrt{k}A=\mathcal{S}.$

In the analysis above and below, we assume that $\mathcal{S}$ belongs to
compact sets where the solution of the Painlevé IV equation does not have
poles. Using (4.15), (2.4) and (2.5), we obtain in the limit $k\to\infty$

$\displaystyle P(z)={N}(z)\begin{pmatrix}0&-1\\\
1&0\end{pmatrix}^{\chi_{L}}\big{(}\sqrt{k}\zeta(z)\big{)}^{-\frac{\gamma}{2}\sigma_{3}}$
$\displaystyle\hphantom{P(z)=}{}\times\left(\mathbf{1}+\dfrac{1}{\sqrt{k}\zeta}\begin{bmatrix}H&\dfrac{Z}{U}\vspace{1mm}\\\
U&-H\end{bmatrix}+\mathcal{O}\big{(}k^{-1}\big{)}\right)\big{(}\sqrt{k}\zeta\big{)}^{\frac{\gamma}{2}\sigma_{3}}\begin{pmatrix}0&1\\\
-1&0\end{pmatrix}^{\chi_{L}}$
$\displaystyle\hphantom{P(z)}{}={N}(z)\\!\left(\mathbf{1}+\begin{pmatrix}0&-1\\\
1&0\end{pmatrix}^{\chi_{L}}\\!\dfrac{1}{\sqrt{k}\zeta}\\!\begin{bmatrix}H&\big{(}\sqrt{k}\zeta\big{)}^{-\gamma}\dfrac{Z}{U}\vspace{1mm}\\\
\big{(}\sqrt{k}\zeta\big{)}^{\gamma}U&-H\end{bmatrix}\\!\begin{pmatrix}0&1\\\
-1&0\end{pmatrix}^{\chi_{L}}\\!+\mathcal{O}\big{(}k^{-1+\gamma/2}\big{)}\right)\\!$
$\displaystyle\hphantom{P(z)}{}={N}(z)\big{(}\mathbf{1}+\mathcal{O}\big{(}k^{\frac{\gamma-1}{2}}\big{)}\big{)}.$
(4.17)

From the above expansion we can see that the subleading terms are not
uniformly small with respect to $\gamma$ as $k\to\infty$ since
$\gamma\in[0,1)$. Note that the source of the non-uniformity in the error
analysis is the element $(2,1)$ in (4.17). For this reason we need to
introduce an improved parametrix in the next section.

####  4.3.5 Improved parametrix

To construct an improved parametrix we define

$\displaystyle\widetilde{{N}}(z)=\left(\mathbf{1}+\frac{B\sigma_{-}}{z-1}\right){N}(z),\qquad\sigma_{-}:=\left[\begin{matrix}0&0\\\
1&0\end{matrix}\right],\qquad\sigma_{+}:=\left[\begin{matrix}0&1\\\
0&0\end{matrix}\right],$ (4.18)

where $B$ is to be determined. In the same way we define

$\displaystyle\widetilde{P}(z)=\widetilde{{N}}(z)\begin{pmatrix}0&-1\\\
1&0\end{pmatrix}^{\chi_{L}}\big{(}\sqrt{k}\zeta(z)\big{)}^{-\frac{\gamma}{2}\sigma_{3}}\left(\mathbf{1}-\dfrac{U\sigma_{-}}{\sqrt{k}\zeta(z)}\right)\widehat{\Psi}\big{(}\sqrt{k}\zeta(z);\sqrt{k}A\big{)},$
(4.19)

where $U$ is the $(2,1)$ entry of the subleading term of the expansion of
$\widehat{\Psi}(\lambda)$ for $\lambda\to\infty$. Now we have in the limit
$k\to\infty$

$\displaystyle\begin{split}&\widetilde{P}(z)=\widetilde{{N}}(z)\left(\mathbf{1}+\begin{bmatrix}0&-1\\\
1&0\end{bmatrix}^{\chi_{L}}\dfrac{1}{\sqrt{k}\zeta}\begin{bmatrix}H&(\sqrt{k}\zeta)^{-\gamma}\dfrac{Z}{U}\vspace{1mm}\\\
0&-H\end{bmatrix}\begin{bmatrix}0&1\\\
-1&0\end{bmatrix}^{\chi_{L}}+\mathcal{O}\big{(}k^{-1+\gamma/2}\big{)}\right)\\\
&\hphantom{\widetilde{P}(z)}{}=\widetilde{{N}}(z)\big{(}\mathbf{1}+\mathcal{O}\big{(}k^{-\frac{1}{2}}\big{)}\big{)}.\end{split}$

The improved parametrics $\widetilde{{N}}(z)$ and $\widetilde{P}(z)$ have the
same jump discontinuities as before, but $\widetilde{P}(z)$ might have poles
at $z=1$. The constant $B$ in (4.18) is then determined by the requirement
that $\widetilde{P}(z)$ is bounded at $z=1$. This gives the constant $B$ as

$\displaystyle B=Uk^{\frac{\gamma-1}{2}}.$

We are now ready to compute the error matrix $E(z)$.

####  4.3.6 Error matrix

The error matrix $E(z)$ is defined as

$\displaystyle
E(z)=\begin{cases}T(z)\widetilde{{N}}(z)^{-1},&z\in\mathbb{C}\setminus\mathbb{D},\\\
T(z)\widetilde{P}(z)^{-1},&z\in\mathbb{D},\end{cases}$ (4.20)

where the boundary of $\mathbb{D}$ is oriented clockwise. The matrix has also
jumps on $\widehat{\Gamma}_{i}\cup\widehat{\Gamma}_{e}\setminus\mathbb{D}$
which are exponentially close to the identity matrix and can be ignored for
the purposes of the analysis. Then, the matrix $E(z)$ satisfies the
Riemann–Hilbert problem

$\displaystyle E_{+}(z)=E_{-}(z)v_{E}(z),\qquad z\in\partial\mathbb{D},\qquad
v_{E}(z)=\widetilde{{N}}(z)\widetilde{P}^{-1}(z).$

In the double scaling limit $k\to\infty$, after using (4.14), (4.15), (2.4),
(2.5) and the last equation in (4.10), the jump matrix $v_{E}(z)$ takes the
form

$\displaystyle
v_{E}(z)=\left(\mathbf{1}+\dfrac{Uk^{\frac{\gamma-1}{2}}}{z-1}\sigma_{-}\right)\left(\dfrac{z-1}{\sqrt{k}z\zeta}\right)^{\frac{\gamma}{2}\sigma_{3}}$
$\displaystyle\hphantom{v_{E}(z)=}{}\times\left(\mathbf{1}-\dfrac{\begin{pmatrix}H&\dfrac{Z}{U}\vspace{1mm}\\\
0&-H\end{pmatrix}}{\sqrt{k}\zeta}+\dfrac{\begin{bmatrix}\dfrac{H^{2}+sH-Z}{2}&\dfrac{HZY+HY+ZY^{2}}{UY}\vspace{1mm}\\\
sU-
UY&\dfrac{H^{2}-sH+Z}{2}\end{bmatrix}}{k\zeta^{2}}+\mathcal{O}\big{(}k^{-\frac{3}{2}}\big{)}\right)$
$\displaystyle\hphantom{v_{E}(z)=}{}\times\left(\dfrac{z-1}{\sqrt{k}z\zeta}\right)^{-\frac{\gamma}{2}\sigma_{3}}\left(\mathbf{1}-\dfrac{Uk^{\frac{\gamma-1}{2}}}{z-1}\sigma_{-}\right)=\mathbf{1}+\dfrac{v_{E}^{(1)}}{\sqrt{k}}+\dfrac{v_{E}^{(2)}}{k^{\frac{1}{2}+\frac{\gamma}{2}}}+\dfrac{v_{E}^{(3)}}{k^{1-\frac{\gamma}{2}}}+\mathcal{O}\big{(}k^{-1}\big{)},$

where the expansion is for $z\in\partial\mathbb{D}$ (while $v_{E}$ is
exponentially close to the identity on $\Gamma_{E}\setminus\mathbb{D}$), and

$\displaystyle v_{E}^{(1)}=-\dfrac{H}{\zeta}\sigma_{3},\qquad
v_{E}^{(2)}=-\left(\dfrac{z-1}{z\zeta}\right)^{\gamma}\left(\dfrac{Z}{U\zeta}\right)\sigma_{+},$
$\displaystyle
v_{E}^{(3)}=\left(\left(\dfrac{z-1}{z\zeta}\right)^{-\gamma}\dfrac{(s-Y)U}{\zeta^{2}}-2\dfrac{HU}{(z-1)\zeta}\right)\sigma_{-},$

where $\sigma_{+}=\left(\begin{smallmatrix}0&1\\\
0&0\end{smallmatrix}\right)$. By the standard theory of small norm
Riemann–Hilbert problems (see for example [17, Chapter 7]), one has a similar
expansion for $E(z)$, namely

$\displaystyle
E(z)=\mathbf{1}+\dfrac{E^{(1)}}{\sqrt{k}}+\dfrac{E^{(2)}}{k^{\frac{1}{2}+\frac{\gamma}{2}}}+\dfrac{E^{(3)}}{k^{1-\frac{\gamma}{2}}}+\mathcal{O}\big{(}k^{-1}\big{)},$

so that

$\displaystyle E^{(j)}_{+}(z)=E^{(j)}_{-}(z)+v_{E}^{(j)}(z),\qquad
z\in\partial\mathbb{D},\qquad j=1,2,3.$

By solving the corresponding Riemann–Hilbert problem, we obtain, using the
Plemelj–Sokhtski formula

$\displaystyle E^{(j)}(z)=\dfrac{1}{2\pi{\rm
i}}\int_{\partial\mathbb{D}}\dfrac{v_{E}^{(i)}(\xi)}{\xi-z}\mathrm{d}\xi,\qquad
j=1,2,3,$

which gives

$\displaystyle
E^{(1)}(z)=-\dfrac{\mbox{Res}_{\xi=1}v_{E}^{(1)}(\xi)}{z-1}=\frac{H\sigma_{3}}{z-1},\qquad
z\in\mathbb{C}\setminus\mathbb{D},$ $\displaystyle
E^{(2)}(z)=-\dfrac{\mbox{Res}_{\xi=1}v_{E}^{(2)}(\xi)}{z-1}=\frac{(Z/U)\sigma_{+}}{z-1},\qquad
z\in\mathbb{C}\setminus\mathbb{D},$ $\displaystyle
E^{(3)}(z)=-\dfrac{\mbox{Res}_{\xi=1}v_{E}^{(3)}(\xi)}{(z-1)}-\dfrac{\mbox{Res}_{\xi=1}(\xi-1)v_{E}^{(2)}(\xi)}{(z-1)^{2}}$
$\displaystyle\hphantom{E^{(3)}(z)}{}=\left(\dfrac{2}{3}\dfrac{\gamma(2H+\gamma-s)\gamma+H+\gamma-s}{z-1}+\dfrac{U(2H+Y-s)}{z-1}\right)\sigma_{-},\qquad
z\in\mathbb{C}\setminus\mathbb{D},$ (4.21)

and

$\displaystyle
E^{(1)}(z)=v_{E}^{(1)}-\dfrac{\mbox{Res}_{\xi=1}v_{E}^{(1)}(\xi)}{z-1}=v_{E}^{(1)}+\frac{H\sigma_{3}}{z-1},\qquad
z\in\mathbb{D},$ $\displaystyle
E^{(2)}(z)=v_{E}^{(2)}-\dfrac{\mbox{Res}_{\xi=1}v_{E}^{(2)}(\xi)}{z-1}=v_{E}^{(2)}+\frac{(Z/U)\sigma_{+}}{z-1},\qquad
z\in\mathbb{D},$ (4.22) $\displaystyle
E^{(3)}(z)=v_{E}^{(3)}+\left(\dfrac{2}{3}\dfrac{\gamma(2H+\gamma-s)\gamma+H+\gamma-s}{z-1}+\dfrac{U(2H+Y-s)}{z-1}\right)\sigma_{-},\qquad
z\in\mathbb{D}.$

###  4.4 Asymptotics for the polynomials $\boldsymbol{\pi_{k}(z)}$ and proof
of Theorems 1.2 and 1.3

We are now ready to determine the asymptotic expansions for the orthogonal
polynomials

$\displaystyle\pi_{k}(z)={\rm
e}^{kg(z)}\left(1-\frac{1}{z}\right)^{\frac{\gamma}{2}}[U(z)]_{11}.$

Using (4.6), (4.7) and (4.20) we have

$\displaystyle\pi_{k}(z)={\rm
e}^{kg(z)}\left(1-\frac{1}{z}\right)^{\frac{\gamma}{2}}\begin{cases}\big{[}E(z)\widetilde{{N}}(z)\big{]}_{11},&z\in(\Omega_{\infty}\cup\Omega_{0})\setminus\mathbb{D},\vspace{1mm}\\\
\left[E(z)\widetilde{{N}}(z)\begin{pmatrix}1&0\\\ {\rm
e}^{k\varphi(z)}&1\end{pmatrix}\right]_{11},&z\in\Omega_{1}\setminus\mathbb{D},\vspace{1mm}\\\
\left[E(z)\widetilde{{N}}(z)\begin{pmatrix}1&0\\\ -{\rm
e}^{k\varphi(z)}&1\end{pmatrix}\right]_{11},&z\in\Omega_{2}\setminus\mathbb{D},\vspace{1mm}\\\
\big{[}E(z)\widetilde{P}(z)\big{]}_{11},&z\in(\Omega_{0}\cup\Omega_{\infty})\cap\mathbb{D},\vspace{1mm}\\\
\left[E(z)\widetilde{P}(z)\begin{pmatrix}1&0\\\ {\rm
e}^{k\varphi(z)}&1\end{pmatrix}\right]_{11},&z\in\Omega_{1}\cap\mathbb{D},\vspace{1mm}\\\
\left[E(z)\widetilde{P}(z)\begin{pmatrix}1&0\\\ -{\rm
e}^{k\varphi(z)}&1\end{pmatrix}\right]_{11},&z\in\Omega_{2}\cap\mathbb{D}.\end{cases}$

We want to analyze each distinct region. Using (4.21), (4.22), (4.19) and
(4.18) we obtain the following expressions:

The region $\boldsymbol{\Omega_{\infty}\setminus\mathbb{D}}$

$\displaystyle\pi_{k}(z)={\rm
e}^{kg(z)}\left(1-\frac{1}{z}\right)^{\gamma}\left(1+\dfrac{H}{\sqrt{k}(z-1)}+\mathcal{O}\left(\frac{1}{k}\right)\right)$
$\displaystyle\hphantom{\pi_{k}(z)}{}=z^{k}\left(1-\frac{1}{z}\right)^{\gamma}\left(1+\dfrac{H}{\sqrt{k}(z-1)}+\mathcal{O}\left(\frac{1}{k}\right)\right),$
(4.23)

with $H$ defined in (4.13).

The region $\boldsymbol{\Omega_{0}\setminus\mathbb{D}}$

$\displaystyle\pi_{k}(z)={\rm
e}^{kg(z)}\left(-\dfrac{Z}{U(z-1)k^{\frac{1}{2}+\gamma}}+\mathcal{O}\left(\frac{1}{k}\right)\right),$

with $Z$ and $U$ defined in (4.13).

The region $\boldsymbol{\Omega_{1}\setminus\mathbb{D}}$

$\displaystyle\pi_{k}(z)={\rm e}^{kg(z)}\left(\frac{z-1}{z}\right)^{\gamma}$
$\displaystyle\hphantom{\pi_{k}(z)=}{}\times\left({\rm
e}^{k\varphi(z)}\left(1+\dfrac{H}{\sqrt{k}(z-1)}\right)-\dfrac{Z}{U(z-1)k^{\frac{1+\gamma}{2}}}\left(\frac{z-1}{z}\right)^{-\gamma}+\mathcal{O}\left(\frac{1}{k}\right)\right),$
(4.24)

with $H$, $U$ and $Z$ defined in (4.13). In a similar way we can obtain the
expansion in the region $\Omega_{2}\setminus\mathbb{D}$.

The region $\boldsymbol{\Omega_{2}\setminus\mathbb{D}}$

$\displaystyle\pi_{k}(z)={\rm
e}^{kg(z)}\left(\frac{z-1}{z}\right)^{\gamma}\left(1+\dfrac{H}{\sqrt{k}(z-1)}-\dfrac{Z{\rm
e}^{k\varphi(z)}}{U(z-1)k^{\frac{1+\gamma}{2}}}\left(\frac{z-1}{z}\right)^{-\gamma}+\mathcal{O}\left(\frac{1}{k}\right)\right).$
(4.25)

The region $\boldsymbol{\mathbb{D}}$. In the region
$(\Omega_{0}\cup\Omega_{\infty})\cap\mathbb{D}$ we have

$\displaystyle\pi_{k}(z)={\rm
e}^{kg(z)}\left(\frac{z-1}{z}\right)^{\gamma}\left(\dfrac{\widehat{\Psi}_{11}\big{(}\sqrt{k}\zeta(z);\sqrt{k}A\big{)}}{k^{\frac{\gamma}{4}}\zeta(z)^{\frac{\gamma}{2}}}+\mathcal{O}\left(\frac{1}{k^{\frac{1}{2}+\frac{\gamma}{4}}}\right)\right),$

where $\widehat{\Psi}_{11}$ is the $11$ entry of the Painlevé isomonodromic
problem (4.8).

In the region $\Omega_{1}\cap\mathbb{D}$ and $\Omega_{2}\cap\mathbb{D}$ we
have

$\displaystyle\pi_{k}(z)={\rm e}^{kg(z)}\left(\frac{z-1}{z}\right)^{\gamma}$
$\displaystyle\hphantom{\pi_{k}(z)=}{}\times\left(\dfrac{\widehat{\Psi}_{11}\big{(}\sqrt{k}\zeta(z);\sqrt{k}A\big{)}\pm{\rm
e}^{k\varphi(z)}\widehat{\Psi}_{12}\big{(}\sqrt{k}\zeta(z);\sqrt{k}A\big{)}}{k^{\frac{\gamma}{4}}\zeta(z)^{\frac{\gamma}{2}}}+\mathcal{O}\left(\frac{1}{k^{\frac{1}{2}+\frac{\gamma}{4}}}\right)\right),$

where $\pm$ refers to the region $\Omega_{1}$ and $\Omega_{2}$, respectively,
and $\widehat{\Psi}_{12}$ is the $12$ entry of the solution Painlevé
isomonodromic problem (4.8). Making the change of variables
$z=1-\frac{\lambda^{d}}{t_{c}}$, the proof of Theorem 1.3 follows in a
straightforward way from the above expansions. With these expansions, we are
able to locate the zeros of the orthogonal polynomials.

Figure 17: Zeros of the polynomials $\pi_{k}(z)$ for $k=40$, $60$, $70$. The
values of the parameters are $z_{0}=1$ and $t=2$, $d=3$ and $\ell=0$. The plot
of the support of the limiting measure (Szegö curve) of the zeros of the
orthogonal polynomials $\pi_{k}(z)$ is in red.

######  Proposition 4.9.

The support of the counting measure of the zeros of the polynomials
$\pi_{k}(z)$ outside an arbitrary small disk $\mathbb{D}$ surrounding the
point $z=1$ tends uniformly to the curve $\Gamma_{r=1}$ defined in (4.2) for
$z_{0}=1$. The zeros are within a distance $o(1/k)$ from the curve defined by

$\displaystyle\log|z|-\frac{|z-1|}{|z_{0}|}=-\dfrac{1+\gamma}{2}\dfrac{\log
k}{k}+\dfrac{1}{k}\log\left(\left|\dfrac{z}{z-1}\right|^{\gamma}\left|\dfrac{Z(\mathcal{S})}{(z-1)U(\mathcal{S})}\right|\right),$
(4.26)

where we recall from (4.12) that $z_{0}=\frac{\sqrt{k}}{\sqrt{k}+\mathcal{S}}$
and that the function $Z=Z(\mathcal{S})$ and $U=U(\mathcal{S})$ are related to
the Painlevé IV equation via (4.13). The curves in (4.26) approach
$\Gamma_{r=1}$ at the rate $\mathcal{O}(\log k/k)$ and lie in
$\operatornamewithlimits{Int}(\Gamma_{r=1})$. The normalized counting measure
of the zeros of $\pi_{k}(z)$ converges to the probability measure $\nu$
defined in (1.9).

###### Proof.

Observing the asymptotic expansion (4.23) of $\pi_{k}(z)$ in
$\Omega_{\infty}\setminus{\mathbb{D}}$, it is clear that $\pi_{k}(z)$ does not
have any zeros in that region, since $z=0$ and $z=1$ do not belong to
$\Omega_{\infty}\setminus{\mathbb{D}}$. The same reasoning applies to the
region $\Omega_{0}\setminus{\mathbb{D}}$, where there are no zeros of
$\pi_{k}(z)$ for $k$ sufficiently large.

From the relations (4.24) and (4.25), one has that in
$\Omega_{1}\cup\Omega_{2}$ using the explicit expression of $g(z)$ defined in
(4.1)

$\displaystyle\pi_{k}(z)=z^{k}\left(\frac{z-1}{z}\right)^{\gamma}\left(1+\dfrac{H}{\sqrt{k}(z-1)}-\dfrac{{\rm
e}^{\pm
k\varphi(z)}Z}{U(z-1)k^{\frac{1+\gamma}{2}}}\left(\frac{z-1}{z}\right)^{-\gamma}+\mathcal{O}\left(\frac{1}{k}\right)\right),$

where $\pm$ refers to $\Omega_{2}$ and $\Omega_{1}$, respectively. The zeros
of $\pi_{k}(z)$ may only lie asymptotically where the expression

$\displaystyle 1+\dfrac{H}{\sqrt{k}(z-1)}=\dfrac{{\rm e}^{\pm
k\varphi(z)}Z}{U(z-1)k^{\frac{1+\gamma}{2}}}\left(\frac{z-1}{z}\right)^{-\gamma},\qquad
z\in\Omega_{1}\cup\Omega_{2}.$

Since $\Omega_{2}\cup\Omega_{1}\subset\\{\operatorname{Re}(\varphi)\leq 0\\}$,
it follows that the zeros of $\pi_{k}(z)$ may lie only in the region
$\Omega_{1}$ and such that $\operatorname{Re}\varphi(z)=\mathcal{O}(\log k/k)$
(where $z_{0}$ is the value given by the double scaling (4.12)). Taking the
logarithm of the modulus of the above equality, we obtain

$\displaystyle\operatorname{Re}\varphi(z)=-\dfrac{1+\gamma}{2}\dfrac{\log
k}{k}+\dfrac{1}{k}\log\left(\left|\dfrac{z}{z-1}\right|^{\gamma}\left|\dfrac{Z(\mathcal{S})}{(z-1)U(\mathcal{S})}\right|\right)+\dfrac{1}{k^{\frac{3}{2}}}\operatorname{Re}\left(\dfrac{H(\mathcal{S})}{z-1}\right)+\mathcal{O}\left(\frac{1}{k^{2}}\right).$

Namely, the zeros of the polynomials $\pi_{k}(z)$ lie on the curve given by
(4.26) with an error of order $\mathcal{O}\big{(}1/k^{2}\big{)}$. Such curves
converge to the curve $\Gamma_{r=1}$ defined in (4.2) with $z_{0}=1$ at a rate
$\mathcal{O}(\log k/k)$. ∎

The proof of Proposition 1.4 follows immediately from the proof of Proposition
4.9. The rest of the proof of Theorem 1.2 follows the steps obtained in [4]
and for this reason we omit it.

### Acknowledgements

The authors wish to thank the anonymous referees for their many suggestions
for improving this manuscript. T.G. acknowledges the support of the
H2020-MSCA-RISE-2017 PROJECT No. 778010 IPADEGAN. M.B. acknowledges the
support by the Natural Sciences and Engineering Research Council of Canada
(NSERC) grant RGPIN-2016-06660 and the FQRNT grant “Matrices Aléatoires,
Processus Stochastiques et Systèmes Intégrables” (2013–PR–166790).

## References

  * [1] Ameur Y., Kang N.-G., Makarov N., Wennman A., Scaling limits of random normal matrix processes at singular boundary points, arXiv:1510.08723.
  * [2] Ameur Y., Seo S.-M., On bulk singularities in the random normal matrix model, Constr. Approx. 47 (2018), 3–37, arXiv:1603.06761.
  * [3] Balogh F., Bertola M., Lee S.-Y., McLaughlin K.D.T.-R., Strong asymptotics of the orthogonal polynomials with respect to a measure supported on the plane, Comm. Pure Appl. Math. 68 (2015), 112–172, arXiv:1209.6366.
  * [4] Balogh F., Grava T., Merzi D., Orthogonal polynomials for a class of measures with discrete rotational symmetries in the complex plane, Constr. Approx. 46 (2017), 109–169, arXiv:1509.05331.
  * [5] Balogh F., Merzi D., Equilibrium measures for a class of potentials with discrete rotational symmetries, Constr. Approx. 42 (2015), 399–424, arXiv:1312.1483.
  * [6] Bambusi D., Maspero A., Birkhoff coordinates for the Toda lattice in the limit of infinitely many particles with an application to FPU, J. Funct. Anal. 270 (2016), 1818–1887, arXiv:1407.4315.
  * [7] Bassom A.P., Clarkson P.A., Hicks A.C., McLeod J.B., Integral equations and exact solutions for the fourth Painlevé equation, Proc. Roy. Soc. London Ser. A 437 (1992), 1–24.
  * [8] Bertola M., Cafasso M., Fredholm determinants and pole-free solutions to the noncommutative Painlevé II equation, Comm. Math. Phys. 309 (2012), 793–833, arXiv:1101.3997.
  * [9] Bleher P., Silva G., The mother body phase transition in the normal matrix model, arXiv:1601.05124.
  * [10] Bornemann F., On the numerical evaluation of Fredholm determinants, Math. Comp. 79 (2010), 871–915, arXiv:0804.2543.
  * [11] Chau L.-L., Zaboronsky O., On the structure of correlation functions in the normal matrix model, Comm. Math. Phys. 196 (1998), 203–247.
  * [12] Claeys T., Grava T., McLaughlin K.D.T.R., Asymptotics for the partition function in two-cut random matrix models, Comm. Math. Phys. 339 (2015), 513–587, arXiv:1410.7001.
  * [13] Claeys T., Vanlessen M., Universality of a double scaling limit near singular edge points in random matrix models, Comm. Math. Phys. 273 (2007), 499–532, math-ph/0607043.
  * [14] Dai D., Kuijlaars A.B.J., Painlevé IV asymptotics for orthogonal polynomials with respect to a modified Laguerre weight, Stud. Appl. Math. 122 (2009), 29–83, arXiv:0804.2564.
  * [15] Deift P., McLaughlin K.T.-R., A continuum limit of the Toda lattice, Mem. Amer. Math. Soc. 131 (1998), x+216 pages.
  * [16] Deift P., Zhou X., A steepest descent method for oscillatory Riemann–Hilbert problems. Asymptotics for the MKdV equation, Ann. of Math. 137 (1993), 295–368.
  * [17] Deift P.A., Orthogonal polynomials and random matrices: a Riemann–Hilbert approach, Courant Lecture Notes in Mathematics, Vol. 3, New York University, Courant Institute of Mathematical Sciences, New York, Amer, Math, Soc,, Providence, RI, 1999.
  * [18] Elbau P., Random normal matrices and polynomial curves, Ph.D. Thesis, ETH Zürich, 2007, arXiv:0707.0425.
  * [19] Elbau P., Felder G., Density of eigenvalues of random normal matrices, Comm. Math. Phys. 259 (2005), 433–450, math.QA/0406604.
  * [20] Etingof P., Ma X., Density of eigenvalues of random normal matrices with an arbitrary potential, and of generalized normal matrices, SIGMA 3 (2007), 048, 13 pages, math.CV/0612108.
  * [21] 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.
  * [22] Fokas A.S., Its A.R., Kitaev A.V., The isomonodromy approach to matrix models in $2$D quantum gravity, Comm. Math. Phys. 147 (1992), 395–430.
  * [23] Ginibre J., Statistical ensembles of complex, quaternion, and real matrices, J. Math. Phys. 6 (1965), 440–449.
  * [24] Gustafsson B., Putinar M., Saff E.B., Stylianopoulos N., Bergman polynomials on an archipelago: estimates, zeros and shape reconstruction, Adv. Math. 222 (2009), 1405–1460, arXiv:0811.1715.
  * [25] Harnad J., Its A.R., Integrable Fredholm operators and dual isomonodromic deformations, Comm. Math. Phys. 226 (2002), 497–530, solv-int/9706002.
  * [26] Hedenmalm H., Makarov N., Coulomb gas ensembles and Laplacian growth, Proc. Lond. Math. Soc. 106 (2013),, arXiv:1106.2971.
  * [27] Hedenmalm H., Wennman A., Planar orthogonal polynomials and boundary universality in the random normal matrix model, arXiv:1710.06493.
  * [28] Its A.R., Izergin A.G., Korepin V.E., Slavnov N.A., Differential equations for quantum correlation functions, Internat. J. Modern Phys. B 4 (1990), 1003–1037.
  * [29] Its A.R., Kapaev A.A., Connection formulae for the fourth Painlevé transcendent; Clarkson–McLeod solution, J. Phys. A: Math. Gen. 31 (1998), 4073–4113.
  * [30] Jimbo M., Miwa T., Ueno K., Monodromy preserving deformation of linear ordinary differential equations with rational coefficients. I. General theory and $\tau$-function, Phys. D 2 (1981), 306–352.
  * [31] Jimbo M., Miwa T., Monodromy preserving deformation of linear ordinary differential equations with rational coefficients. II, Phys. D 2 (1981), 407–448.
  * [32] Kuijlaars A.B.J., Tovbis A., The supercritical regime in the normal matrix model with cubic potential, Adv. Math. 283 (2015), 530–587, arXiv:1412.7597.
  * [33] Lee S.-Y., Makarov N.G., Topology of quadrature domains, J. Amer. Math. Soc. 29 (2016), 333–369, arXiv:1307.0487.
  * [34] Lee S.-Y., Teodorescu R., Wiegmann P., Shocks and finite-time singularities in Hele–Shaw flow, Phys. D 238 (2009), 1113–1128, arXiv:0811.0635.
  * [35] Mehta M.L., Random matrices, Pure and Applied Mathematics (Amsterdam), Vol. 142, 3rd ed., Elsevier/Academic Press, Amsterdam, 2004.
  * [36] Saff E.B., Stylianopoulos N., On the zeros of asymptotically extremal polynomial sequences in the plane, J. Approx. Theory 191 (2015), 118–127, arXiv:1409.0726.
  * [37] Saff E.B., Totik V., Logarithmic potentials with external fields, Grundlehren der Mathematischen Wissenschaften, Vol. 316, Springer-Verlag, Berlin, 1997.
  * [38] Sakai M., Regularity of a boundary having a Schwarz function, Acta Math. 166 (1991), 263–297.
  * [39] Simon B., Trace ideals and their applications, Mathematical Surveys and Monographs, Vol. 120, 2nd ed., Amer. Math. Soc., Providence, RI, 2005.
  * [40] Szegő G., Orthogonal polynomials, Colloquium Publications, Vol. 23, 4th ed., Amer. Math. Soc., Providence, R.I., 1975.
  * [41] Teo L.-P., Conformal mappings and dispersionless Toda hierarchy II: general string equations, Comm. Math. Phys. 297 (2010), 447–474, arXiv:0906.3565.
  * [42] Wasow W., Asymptotic expansions for ordinary differential equations, Dover Publications, Inc., New York, 1987.
  * [43] Wiegmann P.B., Zabrodin A., Conformal maps and integrable hierarchies, Comm. Math. Phys. 213 (2000), 523–538, hep-th/9909147.