The rank 1 real Wishart spiked model

Introduction

In many of these applications, one has to deal with data in which both $N$ and $M$ are large, while the ratio $M/N$ is finite and non-zero. In this case, the sample covariance matrix $Y$ becomes a poor approximation for the covariance matrix $\Sigma$. However, since it is often reasonable to approximate the data by $N$-variate gaussian random variables, a comparison between the eigenvalue distribution of the sample covariance matrix and Wishart matrices with a known covariance matrix will give a good estimate for the spectrum of the true covariance matrix $\Sigma$. In particular, in applications to principle component analysis, one would like to study the asymptotic behavior of the largest eigenvalue of $Y$ as $N$, $M\rightarrow\infty$ with $M/N\rightarrow\gamma^{2}\geq 1$ fixed.

For many statistical data with large $N$ and $M$ and $M/N$ finite, it was noted in that the eigenvalue distribution of the sample covariance matrix is well-approximated by the Marchenko-Pastur law inside a bulk region (See ,, , , .)

where $\chi_{[b_{-},b_{+}]}$ is the characteristic function for the interval $[b_{-},b_{+}]$ and $b_{\pm}=(1\pm\gamma^{-1})^{2}$. However, outside of the bulk region, there are often a finite number of large eigenvalues at isolated locations. This behavior prompted the introduction of spiked models in , which are Wishart matrices with a covariance matrix $\Sigma$ such that only finitely many eigenvalues of $\Sigma$ are different from one. These non-trivial eigenvalues in $\Sigma$ will then be responsible for the spikes that appear in the eigenvalue distribution of the sample covariance matrix . The number of these non-trivial eigenvalues in $\Sigma$ is called the rank of the spiked model.

Of particular interest is a phase transition that arises in the largest eigenvalue distributions when the first of these spikes starts leaving the bulk region. This phenomenon was first studied in for the complex Wishart spiked model and then in for the rank 1 quarternionic Wishart spiked model. Let $1+\tau$ be the non-trivial eigenvalue in $\Sigma$, in both cases, it was shown that the phase transition in the largest eigenvalue distribution occurs when $\tau=\gamma^{-1}$. Their results are expressed in terms of the Hastings-McLeod solution of the Painlevé II equation, which is the unique solution to the Painlevé II equation

In particular, they have proved the following.

For $-1<\tau\leq\gamma^{-1}$, the largest eigenvalue distribution is given by

where $TW_{\beta}(\zeta)$ is the Tracy-Widom distribution.

where $erf(\zeta)$ is the error function.

The functions $TW_{\beta}$ are called the Tracy-Widom distributions in the literature and they give the largest eigenvalue distributions for the real ($\beta=1$), complex ($\beta=2$) and quarternionic ($\beta=4$) Wishart ensembles with $\Sigma=I$, as well as the largest eigenvalue distributions for a large class of random matrix models.

Note that in , the phase transition was in fact computed for spiked models of any finite rank. These previous results naturally divides the range of the non-trivial eigenvalue $1+\tau$ into 3 regimes, which are illustrated in Figure 1.

The subcritical regime: When $-1<\tau<\gamma^{-1}$, the perturbation in $\Sigma$ is not strong enough to form a spike in the eigenvalue distribution of $Y$ and the largest eigenvalue distribution in the Wishart matrix is not affected by this non-trivial eigenvalue. For real Wishart ensembles, this case was studied in and it was shown that the largest eigenvalue distribution remains the same as the case when $\Sigma=I$. The result of also applies to much more general sample covariance matrices that are not necessarily Gaussian.

The critical regime: When $1-\gamma^{-1}\tau=O\left(M^{-\frac{1}{3}}\right)$, the perturbation is just strong enough to form a spike and a phase transition occurs in the largest eigenvalue distribution. The result in for the complex phase transition in fact extends to the whole regime where the authors showed that the largest eigenvalue distribution is described by a function $F_{2}(\zeta,w)$ which equals $TW_{1}^{2}(\zeta)$ when $w=0$. The critical regime for the real case is the subject of this paper. It was also studied recently in using a completely different approach.

The super-critical regime: When $\tau>1$, the perturbation in $\Sigma$ is strong enough to form a spike and the largest eigenvalue is located at the spike instead of the edge of the bulk region. For rank 1 real Wishart spiked model, this was studied in . The results in has also been generalized to a large class of rank 1 spiked perturbed random matrix models in .

Despite having the most applications, the phase transition for real Wishart spiked model has not been solved until very recently . The main goal of this paper is to obtain the asymptotic largest eigenvalue distribution for the rank 1 real Wishart spiked model in the critical regime. In a recent paper , the asymptotic largest eigenvalue distribution for the rank 1 real Wishart ensemble was obtained by using a completely different approach to ours. In , the authors first use the Housefolder algorithm to reduce a Wishart matrix into tridiagonal form. Such tridiagonal matrix is then treated as a discrete random Schrödinger operator and by taking an appropriate scaling limit, the authors obtained a continuous random Schrödinger operator on the half-line. By doing so, the authors in bypass the problem of determining the eigenvalue j.p.d.f. for the real Wishart ensemble and obtain the largest eigenvalue distribution in the asymptotic limit. In , the largest eigenvalue distribution is characterized in several different ways, one of which is the solution of a PDE.

(Theorem 1.7 of ) Let $w=\left(\frac{N}{(1+\gamma^{-1})^{2}}\right)^{\frac{1}{3}}(1-\gamma\tau)\in(-\infty,\infty)$, then the following boundary value problem

has a unique solution $F_{\beta}(\zeta,w)$ and

where $\beta=1$, $2$ and $4$ for the real, complex and quarternionic Wishart ensembles respectively.

The distribution also has another characterization which requires more set up to explain. We will refer the readers to for more details. The above result in in fact also extends to the spiked perturbations of Gaussian ensembles in random matrix theory and to their general $\beta$ analogues, which are defined in .

On the other hand, the approach in this paper uses orthogonal polynomial techniques that are closer to those in and and we obtain a characterization of the largest eigenvalue distribution in the critical regime in terms of Painlevé transcendents. It will be very interesting see how these two representations of the largest eigenvalue distribution can be converted into one another. We will now state our results and outline the method used in the paper.

Statement of result

We will now explain the approach used in this paper.

Let $\lambda_{j}$ be the eigenvalues of the Wishart matrix, then the j.p.d.f. for the real Wishart ensemble is given by

where $g^{T}dg$ is the Haar measure on $O(N)$ and $Z_{M,N}$ is a normalization constant. The j.p.d.f. for complex and quarternionic Wishart ensembles are similar

A major difficulty in the asymptotic analysis of the real Wishart ensembles is to find a simple expression for the j.p.d.f. in terms of its eigenvalues. For complex Wishart ensembles, the integral over $U(N)$ in the expression of the j.p.d.f. can be evaluated using the Harish-Chandra Itzykson Zuber formula to obtain a compact determinantal formula for the j.p.d.f. in terms of the eigenvalues $\lambda_{j}$. This has led to the results in and results on the largest eigenvalue distributions for more general complex Wishart ensembles , . In the quarternonic case, although the Harish-Chandra and Itzykson Zuber formula does not apply, the integral over $Sp(N)$ in the j.p.d.f. can still be evaluated using Zonal polynomial expansion to obtain a determinantal formula for the j.p.d.f. in the rank 1 case . In the real case, however, neither of these methods apply and so far there has not been any simple formula for the j.p.d.f. that allows the asymptotic analysis in the real case. Our first result is a contour integral formula that would allow the asymptotic analysis of the rank 1 real Wishart spiked model.

Assuming $N$ is even. Let the non-trivial eigenvalue in the covariance matrix $\Sigma$ be $1+\tau$ and suppose $\tau\neq 0$. Then the j.p.d.f. of the eigenvalues in the rank 1 real Wishart spiked model with covariance matrix $\Sigma$ is given by

We will present two different proofs of this in the paper. The first one is a geometric proof which involves choosing a suitable set of coordinates on $O(N)$ and decompose the Haar measure into two parts so that the integral in (2.1) can be evaluated. This will be achieved in Sections 3.1 and 3.2. The second proof is an algebraic proof that uses the Zonal polynomial expansion to verify the formula in Theorem 2.1. This proof will be given in Appendix A where integral formulae of the form (2.2) for the complex and quarternionic Wishart ensembles will also be derived.

The integral formula derived here is very similar to a more general formula in , in which the matrix integral over $O(N)$ is given by

Shortly after the first part of this paper, which contains (2.2) appeared in the preprint server ArXiv, we learnt that D. Wang has also derived (2.2) independently and use it for the asymptotic analysis in the super-critical regime .

From the expression of the j.p.d.f., we see that the largest eigenvalue distribution is given by

and $\Gamma$ is a close contour that encloses the interval $[0,z]$.

We can analyze the integrand as in , , and . By an identity of Brujin , we can express the multiple integral as a Pfaffian.

where $r_{j}(x)$ is an arbitrary sequence of degree $j$ monic polynomials and $\left<f,g\right>_{1}$ is the skew product

where $\chi=\chi_{[z,\infty)}$ and $K$ is the operator whose kernel is given by

and $S_{1}(x,y)$ and $IS_{1}(x,y)$ are the kernels

Then we can write these down in terms of Laguerre polynomials.

Let $L_{k}$ be the degree $k$ monic Laguerre polynomial with respect to the weight $w_{0}(x)$

If $\left<L_{2k-1},L_{2k-2}\right>_{1}\neq 0$, then the skew orthogonal polynomials $\pi_{2k,1}$ and $\pi_{2k+1,1}$ both exist and $\pi_{2k,1}$ is unique while $\pi_{2k+1,1}$ is unique up to an addition of a multiple of $\pi_{2k,1}$. Moreover, we have $\left<L_{2k},L_{2k-1}\right>_{1}=0$ and the skew orthogonal polynomials are given by

Next, by representing skew orthogonal polynomials as multi-orthogonal polynomials and write them in terms of the solution of a Riemann-Hilbert problem as in , we can apply the results of and to express the kernel $S_{1}(x,y)$ as a finite rank perturbation of the Christoffel Darboux kernel of the Laguerre polynomials.

Suppose $\left<L_{N-1},L_{N-2}\right>_{1}\left<L_{N-3},L_{N-4}\right>_{1}\neq 0$. Let $S_{1}(x,y)$ be defined by (2.8) and choose the sequence of monic polynomials $r_{j}(x)$ such that $r_{j}(x)$ are arbitrary degree $j$ monic polynomials that are independent on $t$ and $r_{j}(x)=\pi_{j,1}(x)$ for $j=N-2,N-1$. Then we have

where $K_{2}(x,y)$ is the kernel of the Laguerre polynomials

The largest eigenvalue distribution of the rank 1 real Wishart ensemble can be written in the following integral form.

for some constant $c_{0}$ and $K$ is the operator with kernel given by (2.8). The integration contour $\Gamma$ is a close contour that encloses the interval $[0,z]$ in the anti-clockwise direction.

These results so far are valid for all $N$ and $M$ and are exact. As $N$, $M\rightarrow\infty$, the integral expression (2.13) can be used for the asymptotic analysis to obtain the largest eigenvalue distribution when the phase transition occur. In this paper, we will demonstrate the asymptotic analysis for the case where $M/N\rightarrow\gamma^{2}=1$, but the analysis can also be applied to the case of any finite $\gamma$.

In the asymptotic limit, the $t$ integral in (2.13) can be computed using steepest descent analysis. In fact, we shall see that when $t\neq b_{+}$, the integrand in (2.13) is of order

where $\rho(s)$ is given in (1.1). The saddle point $t_{saddle}$ of this equation does not belong to the bulk region $[b_{-},b_{+}]$ unless $\tau$ is at the critical value $\tau_{c}=\gamma^{-1}$. When $\tau<\tau_{c}$, the steepest descent contour $\Gamma$ can be deformed such that it does not intersect $[b_{-},b_{+}]$. In this case, $t$ will always be of finite distance from $[b_{-},b_{+}]$ and the factor $(t-x)^{-\frac{1}{2}}$ in $w(x)$ will have no effect on the asymptotics of the kernel $S_{1}$ at the edge point $b_{+}$. In this case, the kernel at the edge point will be given by the Airy kernel.

where $\xi_{1}=\left(x-b_{+}\right)\frac{\gamma M^{\frac{2}{3}}}{(1+\gamma)^{\frac{4}{3}}}$ and $\xi_{2}=\left(y-b_{+}\right)\frac{\gamma M^{\frac{2}{3}}}{(1+\gamma)^{\frac{4}{3}}}$ are finite.

In the critical case, however, the main contribution to the contour integral (2.13) comes from a small neighborhood of $t=b_{+}$ and the factor $(t-x)^{-\frac{1}{2}}$ in $w(x)$ will now significantly alter the behavior of the kernel $S_{1}$ and change it from the Airy kernel to into a more complicated kernel. This gives rise to the phase transition and a new distribution function. Our next result is the representation of this new distribution function in terms of the Hastings-McLeod solution of the Painlevé II equation. First let us define some functions that will appear in our formula.

Let $\psi(u,T)=(T-u)^{\frac{1}{2}}$ and $H_{j}(u,T)=Ai^{(j)}(u)\psi^{-1}(u,T)$ and let $\mathcal{S}_{10}$, $\mathcal{S}_{11}$ be

Define the function $\mathcal{S}(T)$ to be

where the contour of integration in the first term remains in the upper half plane and $\mathcal{S}_{11,\pm}$ are the boundary values of $\mathcal{S}_{11}$ as it approaches the real axis in the upper/lower half plane.

and define $\vec{h}_{j}$ to be the vector $\vec{h}_{j}=\left(0,0,0,\frac{\phi_{j}}{\psi}\right)^{T}$ for $j=0,1,2$ and $\vec{h}_{j}=0$ for $j=3$ and $j=4$, where $\phi_{0}$ is the Hastings-McLeod solution of Painlevé II (1.3) and

Let $\vec{v}_{j}$ be the vector that satisfies the linear system of ODEs with the following boundary condition

Then the largest eigenvalue distribution at the phase transition is given by

Suppose $N$ is even. Let $w=\left(\frac{N}{4}\right)^{\frac{1}{3}}(1-\tau)\in(-\infty,\infty)$ and let $\zeta=(z-4)\left(N/4\right)^{\frac{2}{3}}$, then as $N,M\rightarrow\infty$ such that $M/N\rightarrow\gamma^{2}=1$, the largest eigenvalue distribution at the phase transition is given by

The distribution in Theorem 2.4 is expressed in terms of integrals of the Airy function, together with the functions $v_{jk}$, which are solutions of linear ODEs with known boundary conditions. The coefficients of the ODEs satisfied by the $v_{jk}$ are given in terms of the Hastings-McLeod solutions and its derivatives and are therefore known functions.

The distribution in Theorem 2.4 can also be expressed in terms of an integral involving the solution of a Riemann-Hilbert problem. (see Appendix B) With the recent advancement in the numerical computation of Riemann-Hilbert problems , this representation may be useful for the numerical computation of the distribution in Theorem 2.4.

In obtaining the main result in Theorem 2.4, we have obtained some other results which may also be of interest to mathematicians and physicists working in random matrix theory.

One applications of the orthogonal polynomial approach developed in this paper is in the studies of random matrix with a rank 1 external source. Random matrix with external source are random matrix models on the space of real symmetric, Hermitian or Hermitian self-dual $N\times N$ matrices with the following probability measure

for some real-valued function $V(x)$ such that $e^{-V(x)}$ decays fast enough as $x\rightarrow\pm\infty$ and a real symmetric, Hermitian or Hermitian self-dual matrix $A$. The function $V(x)$ is usually called the potential while the matrix $A$ is known as the external source. The real Wishart ensemble can be thought of as a special case when $V(x)=\frac{x}{2}-\frac{M-N-1}{2M}\log x$. Random matrices with external source are first studied by Brézin and Hikami , and P. Zinn-Justin , as a model of systems with both random and deterministic parts. For Hermitian random matrices, the external source model can be studied using multi-orthogonal polynomial and Riemann-Hilbert techniques and there are many recent advancements in the asymptotic analysis of Hermitian external source models with non-Gaussian potential $V(x)$ , , , , , , . In , the contour integral formula (2.2) was derived independently and was used to study real symmetric random matrix with a rank 1 external source and a large class of potential $V(x)$. By using the linear statistics results of Johansson and a representation for the largest eigenvalue distribution that is different to (2.13), the largest eigenvalue distribution was obtained in the super-critical regime. The orthogonal polynomial approach developed in this paper can be used to extend the results in to the critical regime. In fact, for real random matrix with a rank 1 external source, the expression (2.13) for the largest eigenvalue distribution remains valid, although the kernels $S_{1}$ and $K$ will have to be modified according to the potential. For a polynomial potential $V(x)$, the analogue of (2.10), which expresses the kernel $S_{1}$ in terms of orthogonal polynomials is well-known . By using the asymptotics of orthogonal polynomials found in and , such representation can be used to compute the asymptotics of the kernel $S_{1}$, which can then be used in (2.13) to obtain the largest eigenvalue distribution.

1.2 Orthogonal ensembles

A large part of this paper involves the analysis of an orthogonal ensemble with weight (2.4) and in doing so, we have obtained some new results for orthogonal ensembles.

in the limit $N\rightarrow\infty$, where $Z_{n,\beta}(w)$ is the partition functions of the ensembles (See remark 2.4 of and remark 1.5 of )

In order to proof the universality in orthogonal and symplectic ensembles using the method in , and , one needs to show that $\lim_{N\rightarrow\infty}\det Q\neq 0$. While the leading order asymptotics of these partition functions can be found using the estimates in , their combined contributions to $\det Q$ cancel in the leading order and hence higher order terms in the asymptotics of $Z_{n,\beta}(w)$ are needed to show that $\lim_{N\rightarrow\infty}\det Q\neq 0$. At the moment, asymptotics for the sub-leading order terms of partition functions are only available for $\beta=2$. By using a differential identity for $\log Z_{N,2}$, the authors in computed the sub-leading order terms in $Z_{N,2}$ for the potential $V(x)=\frac{x^{4}}{4}+\frac{tx^{2}}{2}$ as $N\rightarrow\infty$. A combination of the method in and the differential identity (2.16) may provide a way to compute the sub-leading order terms in the partition function $Z_{N,1}$ and help extend the universality results in orthogonal and symplectic ensembles to more general potentials.

Another interesting observation is Corollary 4.1, in which we showed that $\left<L_{k},L_{k-1}\right>_{1}=0$ whenever $k$ is even. This turns out to be a very useful identity in the analysis of the phase transition when the double scaling limit $t-4=T\left(4/N\right)^{\frac{2}{3}}$ for the ensemble with weight (2.4) has to be considered. When analyzing this double scaling limit, the identity $\left<L_{k},L_{k-1}\right>_{1}=0$ leads to cancelation in the leading order terms of the kernel $S_{1}$. This enables us to show that $S_{1}$ is of order $N^{\frac{2}{3}}$ instead of $N^{\frac{4}{3}}$, which is essential for the scaled limit of the determinant $\det_{2}\left(I-\chi K\chi\right)$ to exist. We believe this type of identity will also be useful in the analysis of other double scaling limits in orthogonal ensembles.

Throughout the paper, we shall assume that $N$ is even and that $M-N>0$.

Contour integral formula for the j.p.d.f.

In this section we will prove the integral formula for the j.p.d.f. in Theorem 2.1.

In this section, we will find a convenient set of coordinate on $O(N)$ to evaluate the integral

that appears in the expression of the j.p.d.f. (2.1). As both $\Sigma^{-1}$ and $Y$ are symmetric matrices, they can be diagonalized by matrices in $O(N)$. We can therefore replace both $\Sigma^{-1}$ and $Y$ by the diagonal matrices $\Sigma_{d}^{-1}$ and $\Lambda_{d}$.

The group $O(N)$ has two connected components, $SO(N)$ and $O_{-}(N)$ that consists of orthogonal matrices that have determinant $1$ and $-1$ respectively. Let $T$ be the matrix

then the left multiplication by $T$ defines an diffeomorphism from $O_{-}(N)$ to $SO(N)$. In particular, we can write the integral over $O(N)$ in (2.1) as

Note that $g^{T}dg$ is also the Haar measure on $SO(N)$.

Let $g_{ij}$ be the entries of $g\in SO(N)$. Then the integral $I$ can be written as

We will now find an expression of the Haar measure and use it to compute the integral $I(\Sigma,\Lambda)$.

First let us define a set of coordinates on $SO(N)$ that is convenient for our purpose. We will then express the Haar measure on $SO(N)$ in terms of these coordinates.

An element $g\in SO(n)$ can be written in the following form

In particular, the matrix $V$ with entries $v_{ij}$ for $1\leq i,j\leq N-1$ is in $SO(N-1)$. The following then gives a set of coordinates on $SO(N)$.

In the above equation, $\vec{g}_{N}$ is identified with the coordinates $\phi_{j}$ in (3.2), while the matrix $V$ is identified with the coordinates in $SO(N-1)$ that correspond to $V$. In terms of these coordinates, the left action of an element $S\in SO(N)$ on $g$ is given by the following.

From the fact that $G(S\vec{g}_{N})SG(\vec{g}_{N})^{-1}$ is an orthogonal matrix, it is easy to check that $\vec{s}=0$ and $s_{N}=\pm 1$. To determine $s_{N}$, let us consider the action of $G(S\vec{g}_{N})SG(\vec{g}_{N})^{-1}$ on $(0,0,\ldots,1)^{T}$. We have

We will now write the Haar measure on $SO(N)$ in terms of the coordinates (3.4). These coordinates give a local diffeomorphism between $SO(N)$ and $S^{N-1}\times SO(N-1)$ as $\vec{g}_{N}\in S^{N-1}$ and $V\in SO(N-1)$. Let $dX$ be a measure on $S^{N-1}$ that is invariant under the action of $SO(N)$ and $V^{T}dV$ be the Haar measure on $SO(N-1)$, then the following measure

as $dX$ is invariant under the action of $S$. Therefore if we can find a measure on $S^{N-1}$ that is invariant under the action of $SO(N)$, then $dX\wedge V^{T}dV$ will give us a left invariant measure on $SO(N)$. Since the left invariant measure on a compact group is also right invariant, this will give us the Haar measure on $SO(N)$. As the metric on $S^{N-1}$ is invariant under the action of $SO(N)$, it is clear that the volume form on $S^{N-1}$ is invariant under the action of $SO(N)$. Let $dX$ be the volume form on $S^{N-1}$, then from (3.7), we see that the measure $dX\wedge V^{T}dV$ is invariant under the action of $SO(N)$.

Let $dX$ be the volume form on $S^{N-1}$ given by

in terms of the coordinates $\phi_{1},\ldots,\phi_{N-1}$ in (3.2) and (3.4), then the Haar measure on $SO(N)$ is equal to a constant multiple of

where $V^{T}dV$ is the Haar measure on $SO(N-1)$ in terms of the coordinates (3.4).

We can now compute the integral $I(\Sigma,\Lambda)$.

2 Integral formula

By using the expression of the Haar measure derived in the last section, we can now write the integral $I(\Sigma,\Lambda)$ as

for some constant $C$, where the $N-1$ sphere $S^{N-1}$ in the above formula is defined by $\sum_{j=1}^{N}g_{jN}^{2}=1$ and $dX$ is the volume form on it. If we let $g_{jN}=x_{j}$, then the above can be written as

Therefore in terms of polar coordinates, we have

To compute the integral $I(\Sigma,\Lambda)$, we use a method in the studies of random pure quantum systems (see, e.g. ). The idea is to consider the Laplace transform of the function $I(\Sigma,\Lambda,t)$ defined by

then $I(\Sigma,\Lambda,1)=I(\Sigma,\Lambda)$. The Laplace transform of $I(\Sigma,\Lambda,t)$ in the variable $t$ is given by

Taking the inverse Laplace transform, we obtain an integral expression for $I(\Sigma,\Lambda)$.

where $\Gamma$ is a contour that encloses all the points $\frac{\tau M}{2(1+\tau)}\lambda_{1},\ldots,\frac{\tau M}{2(1+\tau)}\lambda_{N}$ that is oriented in the counter-clockwise direction. Rescaling the variable $s$ to $s=Mt$, we obtain

This then give us an integral expression for the j.p.d.f. in Theorem 2.1.

In this section we will prove the representation of the kernel $S_{1}$ in Theorem 2.2. We will do so by using the multi-orthogonal polynomial representation of skew orthogonal polynomials in and then apply the Christoffel-Darboux formula for multi-orthogonal polynomials in to write the kernel $S_{1}$ as a finite sum of multi-orthogonal polynomials. We then simplify this sum further by using a result in . This gives a new proof to a well-known result of Widom .

As explain in the introduction, we need to find the skew orthogonal polynomials with the weight (2.4). Let us consider the skew orthogonal polynomials with respect to the weight $w$ in (2.4). We shall use the ideas in to express the skew orthogonal polynomials in terms of a linear combinations of Laguerre polynomials.

Let $Q_{j}(x)$ to be the degree $j+2$ polynomial

Then as we assume $M-N>0$, it is easy to see that

Note that $w_{0}(x)$ is not the square of $w(x)$. The fact that $w_{0}(x)$ is the weight for the Laguerre polynomials allows us to express the skew orthogonal polynomials for the weight (2.4) in terms of Laguerre polynomials.

In particular, this implies that the conditions (2.9) is equivalent to the following conditions

and the exactly same conditions for $\pi_{2k+1,1}(x)$. In particular, the second condition implies the skew orthogonal polynomials can be written as

where $L_{j}(x)$ are the degree $j$ monic Laguerre polynomials that are orthogonal with respect to the weight $w_{0}(x)$.

The constants $\gamma_{k,j}$ are to be determined from the first condition in (4.3). We will now show that if $\left<L_{2k-1},L_{2k-2}\right>_{1}\neq 0$, then the skew orthogonal polynomials $\pi_{2k,1}$ and $\pi_{2k+1,1}$ exist and that $\pi_{2k,1}$ is unique. First let us show that the first condition in (4.3) is equivalent to

To do this, we will first define a map $\varrho_{N}$ from the span of $L_{2k-1}$ and $L_{2k-2}$ to the span of $y$ and $1$.

Let $P(x)$ be a polynomial of degree $m$. Then we can write the polynomial $P(x)$ as

where $q(x)$ is a polynomial of degree $m-2$ and $R(x)$ is a polynomial of degree less than or equal to 1. By writing down the system of linear equations satisfied by the coefficients of $q(x)$ and $R(x)$, we see that the polynomials $q(x)$ and $R(x)$ are uniquely defined for any given $P(x)$. In particular, the map $f:P(x)\mapsto R(x)$ is a well-defined linear map from the space of polynomial to the space of polynomials of degrees less than or equal to $1$. Let $\varrho_{k}$ be the following restriction of this map.

For any polynomial $P(x)$, let $f$ be the map that maps $P(x)$ to $R(x)$ in (4.5). Then the map $\varrho_{k}$ is the restriction of $f$ to the linear subspace spanned by the orthogonal polynomials $L_{k},L_{k-1}$.

If $\left<L_{k},L_{k-1}\right>_{1}\neq 0$, then the map $\varrho_{k}$ is invertible.

Suppose there is exists non-zero constants $a_{1}$ and $a_{2}$ such that

for some polynomial $q(x)$ of degree $k-2$, then by taking the skew product $\left<\right>_{1}$ of this polynomial with $L_{k}$, we obtain

As $q(x)$ is of degree $k-2$. Since $\left<L_{k-1},L_{k}\right>_{1}\neq 0$, this shows that $a_{2}=0$. By taking the skew product with $L_{k-1}$, we conclude that $a_{1}=0$ and hence the map $\varrho_{k}$ has a trivial kernel. ∎

If $k$ is even, then $\left<L_{k},L_{k-1}\right>_{1}=0$.

Let $q(x)$ be a polynomial of degree $k-2$ that satisfies the following conditions

Assuming $k$ is even and let $q(x)$ be a polynomial that satisfies (4.6). By taking the inner product $\left<\right>_{2}$ with $x^{j}$, we see that there exists non-zero constants $a_{1}$ and $a_{2}$ such that

Therefore by Lemma 4.1, we see that if $k$ is even, we will have $\left<L_{k},L_{k-1}\right>_{1}=0$. ∎

Lemma 4.1 shows that if $\left<L_{i},L_{i-1}\right>_{1}\neq 0$, then there exist two independent polynomials $R_{0}(y)$ and $R_{1}(y)$ in the span of $y$ and $1$ such that $R_{j}(y)=\varrho_{i}(L_{i-j})$. Then we have

In particular, the skew product of $R_{j}(y)$ with $L_{i-l}$, $l<2$ is given by

As $q_{j}$ is a polynomial of degree less than or equal to $i-2$ and $l<2$, the first term on the right hand side is zero. Therefore we have

We can now show that the skew orthogonal polynomials $\pi_{2k,1}$ and $\pi_{2k+1,1}$ exist if $\left<L_{2k-1},L_{2k-2}\right>_{1}\neq 0$.

If $\left<L_{2k-1},L_{2k-2}\right>_{1}\neq 0$, then the skew orthogonal polynomials $\pi_{2k,1}$ and $\pi_{2k+1,1}$ both exist and $\pi_{2k,1}$ is unique while $\pi_{2k+1,1}$ is unique up to an addition of a multiple of $\pi_{2k,1}$. Moreover, we have $\left<L_{2k},L_{2k-1}\right>_{1}=0$ and the skew orthogonal polynomials are given by

for $k\geq 2$, where $c$ is an arbitrary constant.

Let $\pi_{2k,1}$ and $\pi_{2k+1,1}$ be polynomials defined by

for some constants $\gamma_{j,k}$. If we can show that $\left<\pi_{2k-l,1},y^{j}\right>_{1}=0$ for $j=0,1$ and $l=-1,0$, then $\pi_{2k-l,1}$ will be the skew orthogonal polynomial. Let $R_{0}$ and $R_{1}$ be the images of $L_{2k-1}$ and $L_{2k-2}$ under the map $\varrho_{2k-1}$. Then by the assumption in the Proposition, they are independent in the span of $y$ and $1$. Therefore the conditions $\left<\pi_{2k-l,1},y^{j}\right>_{1}=0$ are equivalent to $\left<\pi_{2k-l,1},R_{j}(y)\right>_{1}=0$. By taking $i=2k-1$ in (4.7), we see that this is equivalent to $\left<\pi_{2k-l,1},L_{2k-1-j}\right>_{1}=0$. This implies

which exist and are unique as $\left<L_{2k-1},L_{2k-2}\right>_{1}\neq 0$. This determines $\pi_{2k,1}$ uniquely. By Corollary 4.1, we have $\left<L_{2k},L_{2k-1}\right>_{1}=0$ and hence $\gamma_{2k,2}=0$. Similarly, the coefficients for $\pi_{2k+1,1}$ are

Again, these coefficients exist and are unique. However, as $\left<\pi_{2k,1},\pi_{2k,1}\right>_{1}=0$ and $\left<\pi_{2k,1},y^{j}\right>_{1}=0$ for $j=0,\ldots,2k-1$, adding any multiple of $\pi_{2k,1}$ to $\pi_{2k+1,1}$ will not change the orthogonality conditions $\left<\pi_{2k+1,1},y^{j}\right>_{1}=0$ that is satisfied by $\pi_{2k+1,1}$ and hence $\pi_{2k+1,1}$ is only determined up to the addition of a multiple of $\pi_{2k,1}$. ∎

2 The Christoffel Darboux formula for the kernel

In , skew orthogonal polynomials were interpreted as multi-orthogonal polynomials and represented as the solution of a Riemann-Hilbert problem. This representation allows us to use the results in to derive a Christoffel-Darboux formula for the kernel (2.8) in terms of the Riemann-Hilbert problem.

Let us recall the definitions of multi-orthogonal polynomials. First let the weights $w_{0}$, $w_{1}$ and $w_{2}$ be

Note that the weights $w_{l}(x)$ are defined with the polynomials $L_{N-3}$ and $L_{N-4}$ instead of $L_{N-1}$ and $L_{N-2}$. This is because the construction below involves the polynomial $\pi_{N-2,1}$ as well as $\pi_{N,1}$. By taking $i=N-3$ in (4.7), we see that the orthogonality conditions for $\pi_{N,1}$ is also equivalent to

provided $\left<L_{N-3},L_{N-4}\right>_{1}$ is also non-zero.

Then the multi-orthognoal polynomials of type II $P^{II}_{N,l}(x)$, , , , are polynomials of degree $N-1$ such that

More accurately, these are in fact the multi-orthogonal polynomials with indices $(N-2,\vec{e}_{l})$, where $\vec{e}_{1}=(0,1)$ and $\vec{e}_{2}=(1,0)$.

We will now define the multi-orthogonal polynomials of type I. Let $P^{I}_{N,l}(x)$ be a function of the following form

where $B_{N,l}(x)$ is a polynomial of degree $N-3$ and $A_{N,k,l}$ is independent on $x$. Moreover, let $P^{I}_{N,l}(x)$ satisfy

Then the polynomials $B_{N,l}(x)$ and $\delta_{kl}x+A_{N,k,l}$ are multi-orthogonal polynomials of type I with indices $(N,2,1)$ for $l=1$ and $(N,1,2)$ for $l=2$. We will now show that $P^{II}_{N,1}$ and $P^{II}_{N,2}$ exist and are unique if both $\left<L_{N-1},L_{N-2}\right>_{1}$ and $\left<L_{N-3},L_{N-4}\right>_{1}$ are non-zero.

Suppose both $\left<L_{N-1},L_{N-2}\right>_{1}$ and $\left<L_{N-3},L_{N-4}\right>_{1}$ are non-zero, then the polynomials $P^{II}_{N,1}$ and $P^{II}_{N,2}$ exist and are unique.

Then by the first condition in (4.9), we see that

As the matrix with entries $c_{ij}$ is invertible, we see that the linear equations (4.11) has a unique solution in the linear span of $L_{N-1}$ and $L_{N-2}$ if and only if $\left<L_{N-1},L_{N-2}\right>_{1}\neq 0$. ∎

As we shall see, existence and uniqueness of $P^{II}_{N,l}$ would imply that the multi-orthogonal polynomials of type I also exist and are unique. As in , the multi-orthogonal polynomial together with the skew orthogonal polynomials form the solution of a Riemann-Hilbert problem. Let $Y(x)$ be the matrix

Then by using the orthogonality conditions of skew orthogonal polynomials and multi-orthogonal polynomials, together with the jump discontinuity of the Cauchy transform, one can check that $Y(x)$ satisfies the following Riemann-Hilbert problem.

Similarly, multi-orthogonal polynomials of type I can also be arranged to satisfy a Riemann-Hilbert problem. Let $\epsilon$ be the operator

First note that the functions $\psi_{j}(x)=\epsilon(\pi_{j,1}w)$ for $j\geq 2$ can be express in the form of (4.10). By Lemma 4.1, we can write $\pi_{j,1}(x)$ as

where $B_{j}(x)$ is a polynomial of degree $j-2$.

Let $X(z)$ be the matrix value function defined by

Then by using the orthogonality and the the jump discontinuity of the Cauchy transform, it is easy to check that $X^{-T}(z)$ and $Y(z)$ satisfy the same Riemann-Hilbert problem and hence the multi-orthogonal polynomials of type I also exist and are unique.

We will now show that the kernel $S_{1}(x,y)$ given by (2.8) can be expressed in terms of the matrix $Y(z)$.

Suppose $\left<L_{N-3},L_{N-4}\right>_{1}\left<L_{N-1},L_{N-2}\right>_{1}\neq 0$ and let the kernel $S_{1}(x,y)$ be

Then the kernel $S_{1}(x,y)$ exists and is equal to

As in , let us now expand the functions $xr_{j}(x)$ and $x\epsilon(r_{j}w)$.

Then the coefficients $c_{jk}$ and $d_{jk}$ for $j,k=0,\ldots,N-1$ are given by

Therefore if we let $\mathcal{C}$ be the matrix with entries $c_{jk}$, $j,k=0,\ldots N-1$ and $\mathcal{D}$ be the matrix with entries $d_{jk}$ for $j,k=0,\ldots N-1$, then we have

where $r(x)$ is the column vector with components $r_{k}(x)$. Now by (4.24), we see that

which are equal as $\mu^{T}=-\mu$. Now from the form of $\mu$ in (4.22), we see that $\mu_{N-1,j}=\delta_{j,N-2}h_{N-1,1}^{-1}$. Hence we have

Let us now consider the second term in (4.25). As in (4.17) we can write $\epsilon(r_{k}w)w$ as

where $q_{k}(x)$ is of degree $k-2$. Then from the form of $P_{N,j}^{I}$ in (4.10) and the orthogonality condition (4.9), we see that the coefficients $d_{k,N+l}$, $l=1,2$ are given by

For $k\neq N-1$, the polynomial $q_{k}$ is of degree less than or equal to $N-4$, while for $k=N-1$, $B_{N}$ in (4.17) is a polynomial of degree $N-2$ and hence the coefficient $d_{k,N}$ is zero unless $k=N-1$. For $k=N-1$, it is given by the leading coefficient of $q_{N-1}$ divided by the leading coefficient of $B_{N}$. Since

we see that both the leading coefficient of $q_{N-1}(x)$ and $B_{N}$ is $\frac{2}{M}$. Hence $d_{k,N}$ is $\delta_{k,N-1}$. This gives us

To express the second term in (4.25) in terms of the multi-orthogonal polynomials, let us now express $P^{II}_{N,l}$ in terms of the polynomials $r_{k}$. Let us write $P^{II}_{N,l}=\sum_{j=0}^{N-1}a_{j}r_{j}(x)$. Then we have

Hence $P^{II}_{N,l}(x)$ can be written as

Therefore the second term in (4.25) is given by

By using the fact that $Y^{-1}(y)=X^{T}(y)$ and the expressions of the matrix $Y$ (4.12) and $X$ (4.18), together with (4.13), we see that this is the same as (4.21). ∎

3 The kernel in terms of Laguerre polynomials

We will now use a result in to further simplify the expression of the kernel $S_{1}(x,y)$ so that its asymptotics can be computed using the asymptotics of Laguerre polynomials. Let us recall the set up in . First let $Y(x)$ be a matrix satisfying the Riemann-Hilbert problem

and let $\mathcal{K}_{1}(x,y)$ be the kernel given by

Let $\pi_{j,2}(x)$ be the monic orthogonal polynomials with respect to the weight $w_{0}(x)$. Let $\mathcal{K}_{0}(x,y)$ be the following kernel.

We can now apply Proposition 4.3 to our case. In our case, the vectors $\pi(z)$ and $v(z)$ are given by

while the matrix $\mathcal{W}$ is given by

By corollary 4.1, we see that $\left<L_{N-2},L_{N-3}\right>_{1}=0$ and hence the determinant of $\mathcal{W}$ is

Suppose $\det\mathcal{W}\neq 0$ and that the multi-orthogonal polynomials $P^{II}_{N,l}$ exist. Let us now consider the vector $u(z)$. By the Christoffel-Darboux formula, we have

where $h_{j,0}=\left<L_{j},L_{j}\right>_{2}$. We shall show that $L_{N-3}$ and $L_{N-4}$ can be written in the following form

for some polynomial $q_{l}(x)$ of degree $N-1$. By Lemma 4.1, we see that if $\left<L_{N-3},L_{N-4}\right>_{1}\neq 0$, then the map $\varrho_{N-3}$ in Definition 4.1 is invertible. Therefore if the restriction of the map $f$ in Definition 4.1 is also invertible on the span of $\pi_{N+1,1}$ and $\pi_{N,1}$, we will be able to write $L_{N-3}$ and $L_{N-4}$ in the form of (4.31).

Let $f$ be the map in Definition 4.1 and let $\varrho_{\pi}$ be its restriction to the span of $\pi_{N+1,1}$ and $\pi_{N,1}$. Then $\varrho_{\pi}$ is invertible.

Suppose there exist $a_{1}$ and $a_{2}$ such that

for some polynomial $q(x)$ of degree at most $N-1$. Then we have

As the degree of $q(x)$ is at most $N-1$, this is only possible if $q(x)=0$. ∎

The composition of $\varrho_{N-3}$ and $\varrho_{\pi}^{-1}$ will therefore give us a representation of $L_{N-3}$ and $L_{N-4}$ in the form of (4.31). By using this representation and the fact that $q_{l}(x)$ is of degree at most $N-1$, we see that

where $Q$ is the matrix with entries $Q_{i,j}$. We will now determine the constants $Q_{i,j}$.

Let $\pi_{N,1}$ and $\pi_{N+1,1}$ be the monic skew orthogonal polynomial with respect to the weight $w(x)$ and choose $\pi_{N+1,1}$ so that the constant $c$ in (4.8) is zero. Then the vector $u(y)$ in (4.30) is given by

where $Q$ is the matrix whose entries $Q_{i,j}$ are given by

First let us compute the leading order coefficients of the polynomial $q_{l}(x)$ in (4.31). Let $q_{l}(x)=q_{l,N-1}x^{N-1}+q_{l,N-2}x^{N-2}+O(x^{N-3})$, then we have

On the other hand, by orthogonality, we have

Let us now compute $Q_{i,2}$. By taking the skew product, we have

By using the expansion (4.8) of the skew orthogonal polynomials in terms of $L_{N-k}$, we have

By substituting (4.36) into this, we obtain

Note that the matrix $\det(Q\mathcal{W}^{-1})=M^{2}/(4h_{N-1,0}h_{N-2,0})$ and hence $Q\mathcal{W}^{-1}$ is invertible. From Lemma 4.4 and (4.30), we obtain Theorem 2.2.

We will now further simplify the expression for the correlation kernel $S_{1}(x,y)$. First we make the following observation. From the definition (4.19) of the kernel $S_{1}$, it is easy to check that the following identity is true.

We can use this simple observation in the expression (2.10) that we obtained in the last section. First let us write the Christoffel Darboux kernel terms of the solution of a Riemann-Hilbert problem, we have

where $\sigma_{3}$ is the Pauli matrix $\sigma_{3}=\operatorname{diag}(1,-1)$, $\kappa_{0,n}=-\frac{2\pi i}{h_{0,n}}$ and $C(f)$ is the Cauchy transform.

It is well-known that if the logarithmic derivative of $w_{0}$ is rational, then the matrix $Z(x)$ in (4.39) satisfies a linear system of ODE

for some rational function $A(x)$. By using (4.4) and the recurrence relation of the $L_{n}$

Let us now consider the derivative $\frac{\partial K_{2}(x,y)}{\partial y}$. We have, from (4.38)

Similarly, $\frac{\partial}{\partial x}K_{2}(y,x)$ is given by

The sum of these two derivatives is therefore given by

By using the recurrence relations of the Laguerre polynomials, we obtain

Let $K_{1}(x,y)$ be the correction kernel in (2.10).

Then by using $\frac{\partial S_{1}(y,x)}{\partial x}+\frac{\partial S_{1}(x,y)}{\partial y}=0$ and the formula for $h_{0,n}$,

In particular, the correction term $K_{1}(x,y)$ in the kernel can be written as

Derivative of the partition function

where $S_{1}(x,y)$ is the kernel given in (4.19).

Computing the individual determinants using the Laplace formula, we obtain

As $r_{j}(x)$ are independent on $t$ for $j<N-2$, the derivative $\partial_{t}D_{ij}$ is given by

for $i,j<N-2$. For either $i$ or $j$ equal to $N-2$ or $N-1$, we have

Note that by orthogonality, the last two terms in the above expression are zero as $\partial_{t}\pi_{N-1,1}$ is of degree $N-2$ and $\partial_{t}\pi_{N-2,1}$ is of degree $N-3$. Applying the same argument to $\partial_{t}D_{i,N-2}$, we obtain

From this, (5.3) and the expression of the kernel in (4.19), we obtain (5.1). ∎

Asymptotic analysis

We will now use the asymptotics of the Laguerre polynomials to obtain asymptotic expression for the kernel $S_{1}(x,y)$. We will demonstrate the asymptotic analysis with the assumption that $M/N\rightarrow\gamma=1$. The same analysis also be applied for general $\gamma$.

Let $n=N-j$, where $j=O(1)$. Let us define the function $\varphi(x)$ to be

Similarly, the map $f_{n}$ in a neighborhood of $4$ is defined to be

As in , these maps are conformal inside the small discs around $b_{\pm,n}$, provided $\delta$ is sufficiently small. In particular, they behave as follows as near $b_{\pm,n}$. (Recall that $b_{-,n}=0$ and $b_{+,n}=4$ when $\gamma=1$.)

where $J_{\alpha}$ is the Bessel function and $\zeta$ is the following