Extreme Eigenvalue Distributions of Jacobi Ensembles: New Exact Representations, Asymptotics and Finite Size Corrections

TL;DR

This paper develops new exact formulas and asymptotic analyses for the extreme eigenvalue distributions of Jacobi ensembles, providing finite-size corrections and insights into universality, using elementary algebraic methods and properties of Legendre polynomials.

Contribution

It introduces new exact determinant formulas and asymptotic results for eigenvalues of Jacobi ensembles, with finite-size corrections, using elementary methods and Legendre polynomial properties.

Findings

Derived exact distribution formulas involving matrix determinants and Legendre polynomials.

Established large-$n$ asymptotics for extreme eigenvalues under hard-edge scaling.

Provided finite-$n$ corrections that enhance understanding of universality in eigenvalue distributions.

Abstract

Let $\mathbf{W}_1$ and $\mathbf{W}_2$ be independent $n\times n$ complex central Wishart matrices with $m_1$ and $m_2$ degrees of freedom respectively. This paper is concerned with the extreme eigenvalue distributions of double-Wishart matrices $(\mathbf{W}_1+\mathbf{W}_2)^{-1}\mathbf{W}_1$, which are analogous to those of F matrices ${\bf W}_1 {\bf W}_2^{-1}$ and those of the Jacobi unitary ensemble (JUE). Defining $\alpha_1=m_1-n$ and $\alpha_2=m_2-n$, we derive new exact distribution formulas in terms of $(\alpha_1+\alpha_2)$-dimensional matrix determinants, with elements involving derivatives of Legendre polynomials. This provides a convenient exact representation, while facilitating a direct large-$n$ analysis with $\alpha_1$ and $\alpha_2$ fixed (i.e., under the so-called "hard-edge" scaling limit); the analysis is based on new asymptotic properties of Legendre polynomials and their relation with Bessel functions that are here established. Specifically, we present limiting formulas for the smallest and largest eigenvalue distributions as $n \to \infty$ in terms of $\alpha_1$- and $\alpha_2$-dimensional determinants respectively, which agrees with expectations from known universality results involving the JUE and the Laguerre unitary ensemble (LUE). We also derive finite-$n$ corrections for the asymptotic extreme eigenvalue distributions under hard-edge scaling, giving new insights on universality by comparing with corresponding correction terms derived recently for the LUE. Our derivations are based on elementary algebraic manipulations, differing from existing results on double-Wishart and related models which often involve Fredholm determinants, Painlev\'e differential equations, or hypergeometric functions of matrix arguments.

Tables

Table 1. Table 1 : Previous exact results for the extreme eigenvalue distributions of complex Jacobi ensembles.

1. $n$-dimensional determinant	[8]
2. Fredholm determinant	[5, 6]
3. Gauss hypergeometric function of ${\alpha }_1$ or ${\alpha }_2$-dimensional argument	[9]
4. Painlevé VI	[7]
5. Polynomial of degree $n{\alpha }_1$ or $n{\alpha }_2$ involving partitions	[9, 10]

Equations

$$P_l^{(\alpha ,\beta )}(x)=\frac{1}{2^l}\sum _{k=0}^l\left(\genfrac{}{}{0ex}{}{l+\alpha }{k}\right)\left(\genfrac{}{}{0ex}{}{l+\beta }{l-k}\right)(x-1{)}^{l-k}(x+1{)}^l.$$

$${\int }_{-1}^1(1-x{)}^{\alpha }(1+x{)}^{\beta }{P}_k^{(\alpha ,\beta )}(x)P_l^{(\alpha ,\beta )}(x)dx=\frac{2^{\alpha +\beta +1}}{2l+\alpha +\beta +1}\frac{(l+\alpha )!(l+\beta )!}{l!(l+\alpha +\beta )!}{\delta }_{kl}$$

$$P_l(x)=\sum _{k=0}^l\left(\genfrac{}{}{0ex}{}{l}{k}\right)\left(\genfrac{}{}{0ex}{}{l+k}{k}\right){\left(\frac{x-1}{2}\right)}^k.$$

$$P_l^p(x)=(-1{)}^p(1-x^2{)}^{p/2}\frac{d^p{P}_l(x)}{d{x}^p}.$$

$$I_l(z)=\frac{z^l}{2^l}\sum _{k=0}^{\infty }\frac{z^{2k}}{2^{2k}k!(l+k)!}$$

$$\mathbf{W}=({\mathbf{W}}_1+{\mathbf{W}}_2{)}^{-1}{\mathbf{W}}_1$$

$$\prod _{k=1}^n{\varphi }_k^{{\alpha }_1}(1-{\varphi }_k{)}^{{\alpha }_2}\prod _{i<j}^n({\varphi }_i-{\varphi }_j{)}^2$$

$$\prod _{k=1}^n(1+{\tilde{\varphi }}_k{)}^{{\alpha }_1}(1-{\tilde{\varphi }}_k{)}^{{\alpha }_2}\prod _{i<j}^n({\tilde{\varphi }}_i-{\tilde{\varphi }}_j{)}^2$$

$$F_{{\varphi }_n}(\xi )=g_{{\alpha }_1,{\alpha }_2}(\xi )$$

$$F_{{\varphi }_1}(\xi )=1-g_{{\alpha }_2,{\alpha }_1}(1-\xi )$$

$$g_{\alpha ,\beta }(\xi )=1-(1-\xi {)}^{n^2+n\alpha +n\beta +\alpha \beta }\frac{\det \left[{\mathbf{E}}_{\alpha }\left(-\frac{1+\xi }{1-\xi }\right){\mathbf{E}}_{\beta }(1)\right]}{\det \left[{\mathbf{E}}_{\alpha }(-1){\mathbf{E}}_{\beta }(1)\right]}$$

$$[{\mathbf{E}}_{\gamma }(y){]}_{ij}=\frac{d^{j-1}{P}_{n+i-1}\left(y\right)}{d{y}^{j-1}}.$$

$$[{\mathbf{E}}_{\gamma }(y){]}_{ij}=\{\begin{array}{ll}(-1{)}^{j-1}{\left(1-y^2\right)}^{-\frac{j-1}{2}}{P}_{n+i-1}^{j-1}(y),& y\ne 1,y\ne -1\\ 2^{1-j}(n+i-j+1{)}_{2j-2}/(j-1)!,& y=1\\ (-1{)}^{n+i+j}[{\mathbf{E}}_{\gamma }(1){]}_{ij},& y=-1.\end{array}$$

$$F_{{\varphi }_n}(\xi )=1-(1-\xi {)}^{n^2+n{\alpha }_2},$$

$$F_{{\varphi }_n}(\xi )=1-K_{{\alpha }_1}(1-\xi {)}^{n^2+n{\alpha }_1}\det \left({\mathbf{E}}_{{\alpha }_1}\left(\frac{1+\xi }{1-\xi }\right)\right)$$

$$F_{{\varphi }_1}(\xi )={\xi }^{n^2+n{\alpha }_1},$$

$$F_{{\varphi }_1}(\xi )=K_{{\alpha }_2}{\xi }^{n^2+n{\alpha }_2}\det \left({\mathbf{E}}_{{\alpha }_2}\left(\frac{2}{\xi }-1\right)\right).$$

$$F_{{\varphi }_n}(\xi )=h_{{\alpha }_1,{\alpha }_2}(\xi )$$

$$F_{{\varphi }_1}(\xi )=1-h_{{\alpha }_2,{\alpha }_1}(1-\xi )$$

$$h_{\alpha ,\beta }(\xi )=1-\prod _{k=1}^{\alpha }\frac{(n+k-1)!(\alpha -k)!}{(\alpha +n-1)!(k-1)!}(1-\xi {)}^{n^2+n\alpha +n\beta }\det (\mathbf{G})$$

$$\begin{split}[\mathbf{G}]_{ij}=\ &\sum_{k=0}^{\alpha-j}{\alpha-j\choose k}(-2)^{k}(j-i+k+1)_{\alpha-j-k}\left(-\frac{\xi}{1-\xi}\right)^{j-i+k}\frac{d^{k}}{dy^{k}}P_{n+i-1}^{(\alpha-i,\beta-i+1)}\left(y\right)\Big{|}_{y=\frac{1+\xi}{1-\xi}}.\end{split}$$

$$\begin{array}{rl}[\mathbf{G}{]}_{ij}=& \sum _{k=0}^{\alpha -j}\left(\genfrac{}{}{0ex}{}{\alpha -j}{k}\right)(-1{)}^k(j-i+k+1{)}_{\alpha -j-k}(n+\alpha +\beta -i+1{)}_k\\ & \times {\left(-\frac{\xi }{1-\xi }\right)}^{j-i+k}{P}_{n+i-k-1}^{(\alpha +k-i,\beta +k-i+1)}\left(\frac{1+\xi }{1-\xi }\right).\end{array}$$

$$\begin{array}{rl}\mathbb{E}[{\varphi }_n]=& \frac{1}{n^2+n{\alpha }_2+1}\\ =& \frac{1}{n^2}+o\left(\frac{1}{n^2}\right)\end{array}$$

$$\begin{array}{rl}{\sigma }_{{\varphi }_n}=& \sqrt{\frac{n(n+{\alpha }_2)}{(1+n(n+{\alpha }_2){)}^2(2+n(n+{\alpha }_2))}}\\ =& \frac{1}{n^2}+o\left(\frac{1}{n^2}\right).\end{array}$$

$$F_{\infty }^{(\alpha )}(x)=1-e^{-x}\det {\left[I_{i-j}(\sqrt{4x})\right]}_{i,j=1,\dots ,\alpha },$$

$$\underset{n\to \infty }{\lim }{F}_{n^2{\varphi }_n}\left(x\right)=\underset{n\to \infty }{\lim }{F}_{{\varphi }_n}\left(\frac{x}{n^2}\right)=F_{\infty }^{({\alpha }_1)}(x)$$

$$\underset{n\to \infty }{\lim }{F}_{n^2(1-{\varphi }_1)}(x)=1-\underset{n\to \infty }{\lim }{F}_{{\varphi }_1}\left(1-\frac{x}{n^2}\right)=F_{\infty }^{({\alpha }_2)}(x)$$

$$F_{n^2({\tilde{\varphi }}_n+1)}(x)=F_{n^2{\varphi }_n}(x/2),$$

$$F_{n^2(1-{\tilde{\varphi }}_1)}(x)=F_{n^2(1-{\varphi }_n)}(x/2),$$

$$\underset{n\to \infty }{\lim }{F}_{n^2({\tilde{\varphi }}_n+1)}(x)=F_{\infty }^{({\alpha }_1)}(x/2),$$

$$\underset{n\to \infty }{\lim }{F}_{n^2(1-{\tilde{\varphi }}_1)}(x)=F_{\infty }^{({\alpha }_2)}(x/2).$$

Full text

††journal: Journal of Multivariate Analysis

\xpatchcmd\MaketitleBox

\xpatchcmd\MaketitleBox

1 Introduction

Double Wishart random matrices, defined as ${\mathbf{W}=(\mathbf{W}_{1}+\mathbf{W}_{2})^{-1}\mathbf{W}_{1}}$, with $\mathbf{W}_{1}$ and $\mathbf{W}_{2}$ $n\times n$ Wishart with $m_{1}$ and $m_{2}$ degrees of freedom respectively, are an important class of random matrix models. They find application in multivariate analysis of variance (MANOVA), where corresponding test statistics involve the eigenvalues of $\mathbf{W}$, either the complete set or simply the extreme largest/smallest eigenvalues [1]. For linear hypothesis testing, the natural “null hypothesis” considers $\mathbf{W}_{1}$ and $\mathbf{W}_{2}$ independent, central, having identical covariance matrix. The eigenvalues of $\mathbf{W}$ in this case are intimately connected with those of classical Jacobi ensembles and those of Fisher (or F) matrices $\mathbf{F}=\mathbf{W}_{1}\mathbf{W}_{2}^{-1}$, by appropriate variable transformations. Here, we present new results for the extreme eigenvalues of ${\bf W}$ for the case of $\mathbf{W}_{1}$ and $\mathbf{W}_{2}$ being complex Wishart, hence yielding analogous results for the classical Jacobi unitary ensemble (JUE) and complex F model. In addition to their use in statistical testing, the extreme eigenvalues of such complex models arise in multi-antenna communication systems with co-channel interference [2] and in quantum conductance in mesoscopic physics [3, 4].

Our results further contribute to a large amount of prior work on the extreme eigenvalues of double Wishart models (equivalently, JUE/F models). Exact expressions for the extreme eigenvalue distributions of ${\bf W}$ have been given in terms of Fredholm determinants [5, 6], or through equivalent representations in terms of solutions to Painlevé differential equations [5, 7]. Other exact results have been given in terms of $n$-dimensional determinants [8], Gauss hypergeometric functions [9], and polynomial expansions involving combinatorial sums [9, 10]. These results are summarized in Table 1.

The asymptotic distributions of the extreme eigenvalues of $\mathbf{W}$ have also been studied as $n$, $m_{1}$, and $m_{2}$ become large. Particularly noteworthy are results obtained by taking asymptotics on the Fredholm determinant representation [5, 6], a determinant expansion involving the so-called Jacobi kernel [6]. This kernel has been shown to converge to the well-known Bessel kernel [11], when appropriately scaled under the “hard-edge” scaling regime, $n\rightarrow\infty$ with $\alpha_{1}=m_{1}-n$ and $\alpha_{2}=m_{2}-n$ fixed [9, 12, 13]. Consequently, under hard-edge asymptotics, it follows that the extreme eigenvalue distributions of $\mathbf{W}$ can be expressed in terms of a Fredholm determinant involving the Bessel kernel, which has also been shown to admit an equivalent integral form involving the solution of a Painlevé III differential equation [14]. Along a different line, hard-edge asymptotics were evaluated directly in [9], based on the exact hypergeometric function representation, where the smallest eigenvalue distribution of $\mathbf{W}$ was shown to be expressible in terms of a Bessel hypergeometric function of a $\alpha_{1}$-dimensional matrix argument.

It is important to note that an analogous Fredholm determinant representation involving the Bessel kernel has also been established for the smallest eigenvalue distribution of the Laguerre unitary ensemble (LUE) under a similar hard-edge scaling limit [15], suggesting a form of universality among the behavior of the smallest eigenvalue of the LUE and that of the extreme eigenvalues of ${\bf W}$ under hard-edge asymptotics. Remarkably, in the context of the LUE , this asymptotic distribution was shown to admit a very simple representation in terms of a finite-dimensional determinant involving Bessel functions [16]. Such a representation should also apply for the extreme eigenvalues of ${\bf W}$, though this has not been explicitly shown.

In turn, there exists another asymptotic regime, referred to as the “soft-edge” scaling regime, for which $n\rightarrow\infty$ and either $\alpha_{1}\rightarrow\infty$ or $\alpha_{2}\rightarrow\infty$. Under this regime, the Jacobi kernel was shown to converge to the Airy kernel [9]; thus, the extreme eigenvalue distributions of $\mathbf{W}$ can be expressed in terms of a Fredholm determinant involving the Airy kernel. Alternatively, an integral form has been established, involving the solution of a Painlevé II differential equation [17], which is the widely recognized Tracy-Widom law. Also noteworthy is the fact that an analogous representation involving the Airy kernel has been established for the smallest eigenvalue of the LUE [5], so that the said form of universality between the extreme eigenvalues of the JUE and the LUE holds more generally, not only under hard-edge, but also under soft-edge asymptotics. Such universality has been suggested to hold even more generally, for other statistics beyond the extreme eigenvalue distributions [18].

Despite the extensive literature regarding the JUE and LUE, results are scarce when one considers departure from universality. In particular, Edelman, Guionnet and Péché have recently conjectured a first-order correction proportional to $n^{-1}$ for the smallest eigenvalue distribution of the LUE under the hard-edge scaling limit. Independent proofs for this correction have been provided by Perret and Schehr [12] and Bornemann [19]. In a very recent unpublished manuscript [13], Forrester and Trinh studied the optimal scaling for the smallest eigenvalue distribution of the Laguerre $\beta$-ensemble, which subsumes the real ($\beta=1$), complex ($\beta=2$) and symplectic ($\beta=4$) cases, and provide a first-order correction proportional to $n^{-2}$ for $\beta=2$ in integral form, involving a solution to a second-order differential equation. However, to the best of our knowledge, finite-$n$ corrections for the extreme eigenvalue distributions of $\mathbf{W}$ (or the JUE) are not available thus far. A question remains as to whether the universality between LUE and JUE persists when considering finite-$n$ corrections?

A goal of this paper is to provide finite-$n$ corrections for the extreme eigenvalue distributions of $\mathbf{W}$, and therefore, for the F model and the classical JUE, under hard-edge asymptotics. By exploiting new exact representations of the extreme eigenvalue distributions of $\mathbf{W}$, we perform an asymptotic analysis under the hard-edge scaling regime. In the process, we unveil a striking connection of these exact distributions, classically associated with the Jacobi polynomials, with the simpler Legendre polynomials. This new connection allows us to firstly give an explicit proof that shows that the extreme eigenvalue distributions of $\mathbf{W}$ can be expressed in terms of $\alpha_{1}$- and $\alpha_{2}$-dimensional determinants involving Bessel functions, without resorting to study correlation kernels. The proof is a direct one, which takes $n$ large in our new exact formulas, and can be viewed as the “double Wishart analogue” of a similar proof provided for the LUE in [16], but it now boils down to manipulating Legendre polynomials instead of Laguerre polynomials. Secondly, following similar manipulations, we provide finite-$n$ corrections for the extreme eigenvalue distributions of $\mathbf{W}$, giving insights on the universality for the JUE and LUE at the left edge of the spectrum support. To this end, we derive new asymptotic results for the Legendre and associated Legendre polynomials, which are non-standard and may be of independent interest.

1.1 Basic Definitions

The exact results involve Jacobi, Legendre and associated Legendre polynomials. These are defined as follows. The Jacobi polynomial of degree $l$, and parameters $\alpha$ and $\beta$, admits [20, eq. (8.960.1)]

[TABLE]

Jacobi polynomials are orthogonal with respect to the weight $w(x)=(1-x)^{\alpha}(1+x)^{\beta}$ in the interval $[-1,1]$, i.e. [21, eq. (22.2.1)]

[TABLE]

where $\delta_{kl}$ denotes the Kronecker delta function, which equals $1$ if $l=k$ and 0 otherwise.

The Legendre polynomial of degree $l$ admits [20, eq. (8.910.2)]

[TABLE]

Legendre polynomials are particular cases of Jacobi polynomials when $\alpha=\beta=0$. They are then orthogonal with respect to the weight $w(x)=1$ in the interval $[-1,1]$. The associated Legendre polynomial of degree $l$ and order $p$ is defined by [20, eq. (8.810)]

[TABLE]

Our asymptotic results involve the $l$th order modified Bessel function of the first kind, defined by [20, eq. (8.406.3)]

[TABLE]

for $z\in\mathbb{C}$.

1.2 Models

Let $\mathbf{X}_{1}\in\mathbb{C}^{m_{1}\times n}$ ($m_{1}\geq n$) and $\mathbf{X}_{2}\in\mathbb{C}^{m_{2}\times n}$ ($m_{2}\geq n$) be independent complex Gaussian matrices with independent columns that have the same covariance matrix $\mathbf{\Sigma}$. Then, $\mathbf{W}_{1}=\mathbf{X}_{1}^{\dagger}\mathbf{X}_{1}$ and $\mathbf{W}_{2}=\mathbf{X}_{2}^{\dagger}\mathbf{X}_{2}$ are $n\times n$ complex Wishart matrices with $m_{1}$ and $m_{2}$ degrees of freedom respectively; i.e., $\mathbf{W}_{1}\sim\mathcal{CW}_{n}(m_{1},\mathbf{\Sigma})$ and $\mathbf{W}_{2}\sim\mathcal{CW}_{n}(m_{2},\mathbf{\Sigma})$. Define $\alpha_{1}=m_{1}-n$, $\alpha_{2}=m_{2}-n$. The joint probability density function (JPDF) of the eigenvalues of

[TABLE]

is proportional to [6]

[TABLE]

with $1\geq{\phi}_{1}>\ldots>{\phi}_{n}\geq 0$. This JPDF does not depend on $\mathbf{\Sigma}$ since the eigenvalues of $\mathbf{W}$ do not change under the joint transformation $\mathbf{W}_{1}\mapsto\mathbf{\Sigma}^{-1/2}\mathbf{W}_{1}\mathbf{\Sigma}^{-1/2}$, $\mathbf{W}_{2}\mapsto\mathbf{\Sigma}^{-1/2}\mathbf{W}_{2}\mathbf{\Sigma}^{-1/2}$ [6]. The ensemble of $n\times n$ matrices of the form (6) is said to have the multivariate complex beta distribution with parameters $\alpha_{1}$ and $\alpha_{2}$ [22]. The JPDF of its eigenvalues is related to that of other well-known ensembles as follows.

The first is the classical JUE, the ensemble of $n\times n$ matrices with eigenvalue JPDF proportional to

[TABLE]

with $1\geq\widetilde{\phi}_{1}>\ldots>\widetilde{\phi}_{n}\geq-1$. From (7), one obtains (8) by performing the transformation $\phi_{k}=(1+\widetilde{\phi}_{k})/2$, $k=1,\ldots,n$ [6]. The second is the set of random matrices $\mathbf{F}=\mathbf{W}_{1}\mathbf{W}_{2}^{-1}$, commonly referred to as the complex F model [22]. The JPDF of the eigenvalues $\infty>\hat{\phi}_{1}>\ldots>\hat{\phi}_{n}\geq 0$ of $\mathbf{F}$ is obtained from (7) by performing the transformation $\phi_{k}=\hat{\phi}_{k}/(1+\hat{\phi}_{k})$, $k=1,\ldots,n$ [23]. It is said that $\mathbf{F}$ and $\mathbf{W}$ are “matrix analogue” [6].

In this work, we first study the extreme eigenvalue distributions of $\mathbf{W}$, providing new exact determinant expressions which are then leveraged to present asymptotic results and finite-$n$ corrections under hard-edge scaling, i.e., for $n\rightarrow\infty$ with $\alpha_{1}$ and $\alpha_{2}$ fixed. By simply applying the corresponding transformations aforementioned, our results for $\mathbf{W}$ can immediately be rephrased for the JUE and F models.

Before presenting our main results, we make note of the following:

Remark 1**.**

One sees that $1-\phi_{n}$ is the largest eigenvalue of $\mathbf{I}-(\mathbf{W}_{1}+\mathbf{W}_{2})^{-1}\mathbf{W}_{1}=(\mathbf{W}_{1}+\mathbf{W}_{2})^{-1}\mathbf{W}_{2}$. Hence, from the smallest eigenvalue distribution of $\mathbf{W}$, we can deduce that of the largest eigenvalue by simply applying the transformation $\xi\mapsto 1-\xi$ and interchanging $\alpha_{1}$ with $\alpha_{2}$ [22].

1.3 Exact extreme eigenvalue distributions of $\mathbf{W}$

Theorem 1**.**

The cumulative distribution function of the smallest eigenvalue $\phi_{n}$ of $\mathbf{W}$ admits

[TABLE]

and that of the largest eigenvalue $\phi_{1}$ admits

[TABLE]

where

[TABLE]

with $\mathbf{E}_{\gamma}(y)$ the $(\alpha+\beta)\times\gamma$ matrix with entries

[TABLE]

The entries in (12) admit the explicit representation

[TABLE]

Theorem 1 reveals a tight connection between the distributions of the extreme eigenvalues of $\mathbf{W}$ and Legendre polynomials, which are simple particular cases of the Jacobi and Gegenbauer polynomials [24]. The derived exact expressions involve $(\alpha_{1}+\alpha_{2})$-dimensional determinants, whose entries are given exclusively in terms of derivatives of these Legendre polynomials. This has some interesting implications. First, the dimensionality of the determinants in Theorem 1 does not explode when $n$ grows large, if $\alpha_{1}$ and $\alpha_{2}$ are kept fixed. This allows for efficient computation of the extreme eigenvalue distributions in such cases. Moreover, the simple structure of the block matrices inside the determinant in (12) is analogous to a determinant representation derived previously for the smallest eigenvalue distribution of the LUE [16], which is said to be of Wronskian-type (in that case, the successive derivatives were with respect to Laguerre polynomials). Due to this analogy, despite the derivation being more challenging, we will show that we can employ similar manipulations to those presented in [16] to study the large-$n$ behavior of the extreme eigenvalue distributions.

The results of Theorem 1 reduce to simplified forms when either $\alpha_{1}=0$ or $\alpha_{2}=0$.

Corollary 1**.**

When $\alpha_{1}=0$,

[TABLE]

while when $\alpha_{2}=0$,

[TABLE]

where $K_{\alpha}=\prod_{k=0}^{\alpha-1}\frac{2^{k-1}}{(2n+2\alpha-2k)_{k}}$, with $K_{0}=1$.

Similarly, when $\alpha_{2}=0$,

[TABLE]

while when $\alpha_{1}=0$,

[TABLE]

It also turns out that by manipulating known extreme eigenvalue distribution results which were expressed in terms of a Gauss hypergeometric function of a matrix argument (see Table 1), one can also obtain an equivalent expression involving a smaller size determinant, albeit with more complicated entries. Specifically, the expression for the distribution of the smallest eigenvalue involves an $\alpha_{1}$-dimensional determinant, while that for the largest eigenvalue involves an $\alpha_{2}$-dimensional determinant; in both cases the entries involve relatively complicated linear combinations of derivatives of Jacobi polynomials. The result is as follows:

Proposition 1**.**

The cumulative distribution function of $\phi_{n}$ also admits

[TABLE]

and that of $\phi_{1}$ also admits

[TABLE]

where

[TABLE]

with $\mathbf{G}$ the $\alpha\times\alpha$ matrix with entries

[TABLE]

Since the $k$th derivative of a Jacobi polynomial can be expressed in terms of another Jacobi polynomial [20, eq. (8.961.4)], the entries in (21) admit the explicit representation

[TABLE]

It is noteworthy that the simplified special cases (14) and (16) are also easily recoverable from this proposition; however directly recovering the simplified forms (15) and (17) does not appear straightforward. Moreover, due to its simplified structure and dependence on Legendre polynomials (which, recall, are simplified cases of Jacobi polynomials), the results in Theorem 1 are more amenable to direct asymptotic analysis than those given in Proposition 1. We now pursue such asymptotic analysis.

1.4 Asymptotic extreme eigenvalue distributions of $\mathbf{W}$

We consider the hard-edge scaling limit, for which $n$ grows large, with $\alpha_{1}$ and $\alpha_{2}$ fixed. Results under this scaling have been considered previously, where the eigenvalue correlation kernel of the JUE (with appropriate centering and scaling) has been shown to coincide with that of the LUE in this asymptotic limit [11]. This implies that the extreme eigenvalue distributions of the JUE should coincide with the distribution of the smallest eigenvalue of the LUE, which has been shown to admit a remarkably simple form involving a finite-dimensional determinant whose entries are Bessel functions [16]. Using this correspondence, along with the simple mapping between the eigenvalues of ${\bf W}$ and the JUE, given in Section 1.2, an analogous determinant expression should be obtained for the extreme eigenvalue distributions of ${\bf W}$.

Here we provide a direct proof of this result, without resorting to a study of correlation kernels etc., by simply taking $n$ large in our exact formulas for the extreme eigenvalues of $\mathbf{W}$. This is accomplished by deriving asymptotic expansions of Legendre polynomials that are non-standard, and may be of independent interest. In principle, this direct approach is the “double Wishart analogue” (or “JUE analogue”) of a similar direct proof provided for the LUE in [16], which exploited asymptotic properties of Laguerre polynomials. Our derivation, while being based on elementary operations, also enables explicit computation of the large (but finite) $n$ correction terms to the asymptotic distribution. We present this for some particular cases of $\alpha_{1}$ and $\alpha_{2}$.

To guide our asymptotic analysis, it is insightful to first study the scaling of the mean and standard deviation of the smallest eigenvalue, using the simple representation in (14). Specifically, explicit computation of the mean yields

[TABLE]

while, for the standard deviation,

[TABLE]

It is therefore natural to scale $\phi_{n}$ by $n^{2}$ to study its asymptotic distribution (see also [11]). Recall also that the asymptotic distribution of the largest eigenvalue can be deduced from that of the smallest one, as indicated in Remark 1. With this in mind, defining

[TABLE]

we arrive at the following:

Theorem 2**.**

For fixed $\alpha_{1}$ and $\alpha_{2}$,

[TABLE]

and

[TABLE]

for $x\geq 0$.

Contrasting this result with Theorem 1, where the exact extreme eigenvalue distributions of $\mathbf{W}$ were given in terms of $(\alpha_{1}+\alpha_{2})$-dimensional determinants, the asymptotic distributions in Theorem 2 involve $\alpha_{1}$- or $\alpha_{2}$-dimensional determinants, as in Proposition 1. However, contrary to Proposition 1, other than in defining the determinant size, there is no further dependence on either $\alpha_{1}$ or $\alpha_{2}$ in the asymptotic expression. Hence, for both the largest and smallest eigenvalue distributions, the dependence on one of the alphas is fully washed out when taking hard-edge asymptotics, while the dependence on the other is simply to determine the dimensionality of the matrix determinant.

If one now considers the JUE, for which $\tilde{\phi}_{k}=2\phi_{k}-1$ (see Section 1.2), we easily establish that

[TABLE]

and therefore

[TABLE]

These asymptotic distributions coincide precisely with the smallest eigenvalue distribution of the LUE under similar hard-edge scaling (with suitable parameterization of $\alpha_{1}$ and $\alpha_{2}$), suggesting a form of “universality” under the hard-edge scaling limit. While this is aligned with previous results relating the hard-edge scaling of the JUE and LUE [11], an open question is whether such correspondence persists when considering first-order correction terms to the asymptotic distribution? Recently, these correction terms were computed explicitly for the LUE [25], though for the JUE (or the double Wishart model) we are unaware of any corresponding results. The explicit exact eigenvalue distribution in Theorem 1 lends itself to this analysis, at least for specific values of $\alpha_{1}$ or $\alpha_{2}$, as we present below. A generalized formula for arbitrary $\alpha_{1}$, $\alpha_{2}$ may also be possible, although we have been unable to establish a generalized proof at this point.

We first recall that the density corresponding to the asymptotic distribution (25) admits [16]

[TABLE]

Our main result is the following:

Proposition 2**.**

For $\alpha_{1}=0,1$ and arbitrary (but fixed) $\alpha_{2}$,

[TABLE]

This also holds for $\alpha_{1}=2$ and $\alpha_{2}=0,1,2$.

Similarly, for arbitrary $\alpha_{1}$ and $\alpha_{2}=0,1$,

[TABLE]

which also holds for $\alpha_{1}=0,1,2$ and $\alpha_{2}=2$.

Proposition 2 reveals that, for the cases of $\alpha_{1}$,$\alpha_{2}$ considered, the first-order correction of the extreme eigenvalue distributions of $\mathbf{W}$ is proportional to the density (30) which, similar to the asymptotic distributions of Theorem 2, is given as an $\alpha_{1}$- or $\alpha_{2}$-dimensional determinant. Interestingly, the alpha parameter which was washed out in those asymptotic expressions appears when considering finite-$n$ corrections, as a scaling factor in the first-order correction term of (31)-(32). From the equivalence (1.4), Proposition 2 can be immediately rephrased for the JUE. Focusing in particular on the smallest eigenvalue, and for the cases of $\alpha_{1}$,$\alpha_{2}$ considered in the proposition,

[TABLE]

which bears a strong analogy with a recent corresponding result for the LUE, conjectured in [25] and proved in [12, 19]; specifically, for the LUE with fixed parameter $\alpha$, the distribution of the smallest eigenvalue for large (but finite) $n$ is given by [25, Theorem 4.2]

[TABLE]

which coincides with that of the JUE when $\alpha_{1}=\alpha$ and $\alpha_{2}=0$. Therefore, Proposition 2 shows that the correspondence under the hard-edge scaling between the extreme eigenvalue distributions of the JUE and LUE still holds for finite-$n$ corrections to first order, at least for the specific values of $\alpha_{1},\alpha_{2}$ considered.

A natural question is whether the result of Proposition 2 and, therefore, the suggested universality of the first-order corrections under hard-edge scaling is still valid for arbitrary $\alpha_{1}$ and $\alpha_{2}$. The proof of the general case is particularly challenging, due to the overwhelming number of terms that appear in the iterative procedure to reduce the dimensions of the involved determinants (see Section 6). Although we have not been able to establish such proof, in the following we present some numerical results which both validates Proposition 2, and checks numerically whether the stated first-order corrections may continue to hold beyond the cases of $\alpha_{1}$, $\alpha_{2}$ considered in Proposition 2.

We first computed the empirical density of the smallest eigenvalue of $\mathbf{W}$ for the cases $\alpha_{1}=2$ and $\alpha_{2}=0,1$. This was computed from 50 million realizations of $\mathbf{W}$ for $n=20$ and $n=100$; the simulations took $11$ hours for $n=20$ and $37$ hours for $n=100$ on a $12$-core computer. In Figs. 1(a) and 1(b) we show the empirical correction, computed as the difference between the empirical density and the theoretical asymptotic density $f_{\infty}^{(\alpha_{1})}(x)=\frac{d}{dx}F_{\infty}^{(\alpha_{1})}(x)$, scaled111We find it convenient to scale the correction term by $ne^{x}$, as opposed to simply by $n$, to cancel the exponential factor that appears from the derivative in (35), which would otherwise dominate the behavior of the correction term, rendering numerical validations visually less clear. by $ne^{x}$, along with the theoretical first-order correction to the asymptotic density, obtained from (31) in Proposition 2 (and correspondingly scaled) as

[TABLE]

which gives, for $\alpha_{1}=2$ and arbitrary $\alpha_{2}$,

[TABLE]

As expected, the simulated correction approaches the theoretical first-order correction as $n$ increases, since the contribution of higher-order terms in the simulated correction decreases. The agreement between the simulated and theoretically-predicted correction is already evident at $n=100$.

To further evaluate whether the theoretical first-order correction holds beyond the cases of Proposition 2, we again computed the empirical density and compared the empirical correction with the theoretical one, just as in Figs. 1(a) and 1(b), but now for $\alpha_{1}=3$ and $\alpha_{2}=0,2$. The results are presented in Figs. 1(c) and 1(d), which again show an excellent agreement between simulated and theoretically-predicted corrections, suggesting that Proposition 2 may hold in general. This is formally conjectured as follows:

Conjecture 1**.**

For arbitrary $\alpha_{1}$ and $\alpha_{2}$,

[TABLE]

With this, we equivalently conjecture that the first-order corrections to the asymptotic distribution for the extreme eigenvalues of the JUE are indeed equivalent to those for the smallest eigenvalue of the LUE in general, upon suitable JUE-LUE parametrization; recall that the LUE is parametrized by a single alpha, so that for the equivalence with the JUE to hold, we must either have $\alpha_{1}=0$ or $\alpha_{2}=0$ when respectively considering the largest or the smallest eigenvalue of the JUE. When both $\alpha_{1}$ and $\alpha_{2}$ are non-zero, the suggested universality of the first-order corrections under hard-edge scaling does not persist, due to the scaling factor $(\alpha_{1}+\alpha_{2})$ in the correction terms given in Conjecture 1.

A further interesting question is whether the correspondence between the LUE and JUE, and the suggested universality under the hard-edge scaling, still hold for second-order (or higher-order) correction terms, upon suitable JUE-LUE parametrization. An answer to this question requires substantial further analysis, and remains an interesting topic for future investigation.

2 Legendre Polynomials and Bessel Functions

Our analysis relies heavily on properties of Legendre polynomials and their asymptotic connection to Bessel functions. In this section, we summarize the properties needed for the proofs of Theorem 1, Theorem 2 and Proposition 2.

2.1 Additional Definitions

The Legendre polynomial $P_{n}(x)$ is alternatively defined by the Rodrigues’ formula [21, eq. (22.11.5)]

[TABLE]

with $P_{n}(1)=1$, where it is clear that [21, eq. (22.4.6)]

[TABLE]

The associated Legendre polynomial $P_{n}^{-m}(x)$ is defined by the Rodrigues’ formula

[TABLE]

where $n\geq m$.

The shifted Legendre polynomial of degree $n$ is defined by

[TABLE]

Shifted Legendre polynomials are orthogonal with respect to $1$ in the interval $[0,1]$, i.e., [21, eq. (22.2.11)]

[TABLE]

2.2 Identities

Lemma 1**.**

For $m\geq-(n+1)$ and $n\geq 0$

[TABLE]

Proof 1**.**

To prove such result, we manipulate recurrence properties of associated Legendre polynomials. We start with [20, eq. (8.731.1)]

[TABLE]

and [20, eq. (8.731.1(1))]

[TABLE]

Applying the Rodrigues’ formula (41) and the chain rule, we rewrite (45) and (46) as

[TABLE]

and

[TABLE]

respectively. Replacing $n$ with $n+1$ in (48), multiplying by $x$ and then subtracting (47), we obtain the result after some simplifications. ∎

Corollary 2**.**

For $m\geq 0$,

[TABLE]

Proof 2**.**

The proof is straightforward from Lemma 1 and the Rodrigues’ formula (39).

Corollary 2 is a special case of Lemma 1 and will be key to give insight on the proof of Theorem 2 in Section 5.1.

Lemma 2**.**

For $-n\leq m\leq n$,

[TABLE]

Proof 3**.**

This follows from (41) and [20, eq. (8.752.2)]

[TABLE]

∎

2.3 Asymptotics

In [25, Lemma 4.1.], Edelman et al. provided the following two-term asymptotic expansion for Laguerre polynomials

[TABLE]

which extends the result [16, eq. (3.29)], with $L_{l}^{(p)}(x)$ the associated Laguerre polynomial of degree $l$ and order $p$. Here, we provide an analogous property for derivatives of Legendre polynomials and associated Legendre polynomials. Our results extend the classical result by Laurent [26, Section IV]

[TABLE]

Lemma 3**.**

For fixed $m\geq 0$, $c\in\mathbb{Z}$ and $x>0$,

[TABLE]

For fixed $m,c\in\mathbb{Z}$ and $x\geq 0$,

[TABLE]

where $\iota^{2}=-1$.

Proof 4**.**

First, we prove (54) by following the strategy of [26, Section IV]. Let $n>m$. Using (39), we rewrite

[TABLE]

Applying Leibnitz formula for the $(n+c+m)$-times differentiation of the product of functions $f(y)=(y+1)^{n+c}$ and $g(y)=(y-1)^{n+c}$, i.e.,

[TABLE]

yields

[TABLE]

or, equivalently,

[TABLE]

Let

[TABLE]

we have

[TABLE]

Applying (40), we obtain

[TABLE]

Now, we are ready to make $n$ large. Since

[TABLE]

[TABLE]

for $k=0,1,2,\ldots$, and

[TABLE]

we obtain that

[TABLE]

where we have only presented the $(k+2)$th term of (62) with

[TABLE]

Therefore,

[TABLE]

which leads to the result (54) with the help of (5).

Using (4), we have

[TABLE]

where

[TABLE]

Combining the limiting properties of (68) and (70), we obtain the result (55).∎

The second term of the expansion (54) depends on a modified Bessel function of one order less than that of the leading term. This is in contrast to the second term of the expansion (52), which is given by a modified Bessel function of two order less than that of the leading term.

As a by-product of Lemma 3, we also give the following Corollary, which presents results that have not been reported elsewhere, to the best of our knowledge222A similar result to (73) was presented in [27, p. 156 eq. (3)] without proof. Although that result relates Bessel functions with associated Legendre polynomials when their arguments lie outside $[-1,1]$ as $n$ grows, it involves a different argument, omits the complex constant $\iota^{m}$ and is not valid for the whole range of values indicated in [27]. That result should instead read

$\lim_{n\rightarrow\infty}n^{m}P_{n}^{-m}\left(\cosh\left(\frac{x}{n}\right)\right)=\iota^{m}I_{m}(x).$

(71) ..

Corollary 3**.**

For fixed $m\geq 0$ and $x>0$,

[TABLE]

For fixed $m\in\mathbb{Z}$ and $x\geq 0$,

[TABLE]

Remark 2**.**

Although Corollary 3 presents asymptotic results for Legendre polynomials of degree $n$ when ${n\rightarrow\infty}$, they are also valid for Legendre polynomials of degree $n+c$ when $n\rightarrow\infty$, with fixed $c\in\mathbb{Z}$, as shown in Lemma 3. We will use this in the proof of Theorem 2.

3 Proof of Theorem 1

We make use of the following result:

Lemma 4**.**

Let $w(y)$ be a non-negative function with all its moments finite and $t_{l}\in\mathbb{R}$ for all $l$, with $t_{l}\neq t_{k}$ for all $l\neq k$. Then [28, eq. (22.4.11)],

[TABLE]

where $K$ is a normalization constant and $\pi_{n}$ is the $n$th order monic polynomial orthogonal with respect to the weight function $w(y)$ in the integration interval.

First, we prove the result (9) for the smallest eigenvalue. For the most part, the proof follows the strategy of [29, Appendix A], which considered the smallest eigenvalue distribution of the non-central complex Wishart model with rank-$1$ mean.

We start by writing

[TABLE]

where $0\leq\xi<1$ and

[TABLE]

with $C$ the normalization constant of the eigenvalue JPDF in (7). Since the integrand is symmetric in $x_{1},\ldots,x_{n}$, we may write

[TABLE]

After the multiple change of variables $y_{i}=(x_{i}-\xi)/(1-\xi)$, for $1\leq i\leq n$, we have

[TABLE]

Rearranging the expression yields

[TABLE]

where

[TABLE]

with

[TABLE]

Note that this is of the same form as (74) in Lemma 4; however, we cannot apply the lemma directly since $t_{l}$, $l=1,\ldots,\alpha_{1}+\alpha_{2}$, are not distinct. To proceed, first recognize that if $t_{l}$ for all $l$ were distinct, then

[TABLE]

where $\widetilde{C}$ is a normalization constant and $\widetilde{P}_{\nu}(x)$ is the $\nu$th order polynomial orthogonal with respect to $1$ in $[0,1]$. This is precisely the shifted Legendre polynomial, defined in Section 2. Our desired integral in (80) can be evaluated from (82) by taking limits as

[TABLE]

To evaluate these limits, we apply [30, Lemma 2] and we have

[TABLE]

where $\mathbf{\Upsilon}(\beta,\gamma)$ is a $(\alpha_{1}+\alpha_{2})\times(\alpha_{1}+\alpha_{2})$ matrix defined by

[TABLE]

with $\mathbf{D}_{\alpha}(y)$ a $(\alpha_{1}+\alpha_{2})\times\alpha$ matrix with entries

[TABLE]

and

[TABLE]

The product $\prod_{i=1}^{\alpha-1}i!$ is taken as $1$ when $\alpha=0$.

Using (84), (79) and (75), we obtain

[TABLE]

Since $F_{\phi_{n}}(0)=0$,

[TABLE]

With the help of (42) and the chain rule, we write the derivatives of the shifted Legendre polynomials in (86) in terms of standard Legendre ones as

[TABLE]

which gives the result for the smallest eigenvalue distribution in (9). The result (10) follows from (9) by applying the transformation in Remark 1.

4 Proof of Proposition 1

From [9, eq. (3.16)] and [31],

[TABLE]

where333As mentioned in [9], one can write $h_{\alpha_{1},\alpha_{2}}(\xi)$ in terms of a polynomial in $\xi/(1-\xi)$. However, this polynomial is difficult to compute since it involves a sum over all partitions $\kappa$ of $k\in\{0,1,\ldots,n\alpha_{1}\}$ into no more than $\alpha_{1}$ parts.

[TABLE]

with ${}_{2}\widetilde{F}_{1}$ the $\alpha_{1}$-dimensional complex Gauss hypergeometric function. First recognize that if $s_{j}$, $j=1,\ldots,\alpha_{1}$, were distinct in (92), then [32, eq. (2.9)]

[TABLE]

where $\alpha_{1,2}=\alpha_{1}+\alpha_{2}$ and ${}_{2}\mathcal{F}_{1}$ is the Gauss hypergeometric function of scalar argument. Our desired expression in (92) can be evaluated from (93) by taking limits as

[TABLE]

where the Gauss hypergeometric function of scalar argument can be expressed in terms of Jacobi polynomials as [21, eq. (15.4.6)]

[TABLE]

To evaluate these limits, we apply [30, Lemma 2] and we have

[TABLE]

Finally, we obtain the result by applying the Leibniz rule to the entries of the determinant.

5 Proof of Theorem 2

As before, we prove the result (26) for the smallest eigenvalue, with the result (27) then following from Remark 1.

First consider the case $\alpha_{2}=0$. Applying (40) to the entries of $\mathbf{E}_{\alpha_{1}}(y)$ and some algebraic simplifications, we obtain

[TABLE]

where $x^{\star}=(1+x/n^{2})/(1-x/n^{2})$.

If one takes $n$ large and applies Corollary 3 to the entries of the numerator determinant of (97), given in (12), we obtain

[TABLE]

while for the denominator, with (13), we obtain

[TABLE]

Hence, replacing the entries of both determinants with their leading order terms gives clearly a $0/0$ indetermination in (97). To circumvent this, we iteratively make a set of manipulations to the determinant of $\mathbf{E}_{\alpha_{1}}(y)$ by using properties of Legendre polynomials (Lemmas 1 and 2 and Corollary 2), following a similar method as in [16] for the Laguerre case. In particular, we iteratively make row operations and use the recurrence property in Corollary 2 to modify the derivative orders of the entries within a specific column, so that when $y=x^{\star}$, by virtue of Corollary 3, they approach a different Bessel function in the limit, which avoids the $0/0$ indetermination. However, the recurrence property in Corollary 2 for the Legendre case presents a certain range of validity, which will prevent its application for some entries. Also, contrary to the recurrence property of Laguerre polynomials used in [16], having the constant $(n-m+1)$ on the left-hand side of Corollary 2 makes the derivation more cumbersome. We first demonstrate the result for the case $\alpha_{1}=3$, in order to shed light on the set of manipulations required to prove the general result.

5.1 An illustrative case: $\alpha_{1}=3$ and $\alpha_{2}=0$

This corresponds to the simplest case that demonstrates the challenge posed by the recurrence property of Corollary 2 or, more generally, Lemma 1. For $\alpha_{1}=3$ and $\alpha_{2}=0$, we have

[TABLE]

When $y=x^{\star}$, by virtue of Corollary 3, we identify for large $n$

[TABLE]

We apply a set of iterative operations which will successively decrease the order (in $n$) from one row to the next in (101). This will make use of the recurrence properties of Legendre polynomials in Lemma 1 and Corollary 2. Although we could apply the more general Lemma 1 instead, Corollary 2 will be useful to illustrate the purpose of each iteration. In the first iteration, to facilitate the application of Corollary 2, we scale the third row of $\mathbf{E}_{3}(y)$ by $y$ and then subtract the second row. We then scale the second row by $y$ and then subtract the first row. Note that this does not alter the first row. This procedure yields $\det(\mathbf{E}_{3}(y))=y^{-2}\det(\mathbf{\Theta}_{n}^{(1)}(y))$ with

[TABLE]

We can now apply Corollary 2 to the modified entries of the second and third columns. For the modified entries of the first column, we use their Rodrigues’ formula representation (39) and then employ Lemma 1. This leads to

[TABLE]

concluding the first iteration. In the second iteration, we repeat the same manipulations, but this time we only scale the third row by $y$ and subtract the second row. This gives $\det(\mathbf{E}_{3}(y))=y^{-3}\det(\mathbf{\Theta}_{n}^{(2)}(y))$ with

[TABLE]

where

[TABLE]

We then employ Lemma 1 and Corollary 2 to rewrite the entries $a_{1,j}(y)$. Specifically, we rewrite $a_{1,1}(y)$ as

[TABLE]

where the first term of the last line followed from Lemma 1. The entries $a_{1,2}(y)$ and $a_{1,3}(y)$ are handled similarly, by employing Lemma 1 and Corollary 2 respectively, giving

[TABLE]

At this point, we apply Lemma 2 to the entries below the main diagonal of $\mathbf{\Theta}^{(2)}_{n}(y)$ to obtain

[TABLE]

Note that the entries of the third row of $\mathbf{\Theta}^{(2)}_{n}(y)$ have an additional term with respect to the previous rows as a result of the manipulations. This is due to the fact that the recurrence formula of Legendre polynomials in Lemma 1 (or Corollary 2) has the constant factor $(n-m+1)$ in its left-hand side. We will need to consider such additional terms when taking limits.

When $y=x^{\star}$, $1-y^{2}=-4x/n^{2}$. With this in mind, by virtue of Corollary 3, we now have

[TABLE]

where in contrast to (101), as alluded earlier, the order in $n$ has been successively reduced in the second and third rows. This is a consequence of iteratively applying Corollary 2 or, more generally, Lemma 1. Effectively, when applying Corollary 2 to the modified upper-triangular entries, the order in $n$ is reduced by one. This can be seen from the reduction in the derivative order of the Legendre polynomials, which reduces the order in $n$ by two (Corollary 3), and from the factor $(n-m+1)$ (left-hand side of Corollary 2), which increases that order by one. The same effect occurs when applying Lemma 1 to the lower-triangular entries, even though Lemma 1 does not explicitly show this order reduction.

Next, we perform some manipulations in order to apply the limiting results of Corollary 3, while making the entries all of order $1$. We divide the $j$th column by $n^{j-1}(1-y^{2})^{(3-j)/2}$ for $j=1,\ldots,3$, and multiply the $i$th row by $(1-y^{2})^{(3-i)/2}$ for $i=1,\ldots,3$. This gives $\det(\mathbf{E}_{3}(y))=(y/n)^{-3}\det(\mathbf{\tilde{\Theta}}_{n}^{(2)}(y))$ with

[TABLE]

At this point, by virtue of Corollary 3, the columns of $\mathbf{\tilde{\Theta}}^{(2)}_{n}(x^{\star})$ are linearly independent when $n\rightarrow\infty$. We then rewrite (97) for $\alpha_{1}=3$ and $\alpha_{2}=0$ as

[TABLE]

and take the $n\rightarrow\infty$ limit. Specifically, applying Corollary 3 to the entries of $\mathbf{\tilde{\Theta}}^{(2)}_{n}(x^{\star})$, and recalling that $1-y^{2}=-4x/n^{2}$ when $y=x^{\star}$, we obtain that $\lim_{n\rightarrow\infty}\det({\mathbf{\tilde{\Theta}}}^{(2)}_{n}(x^{\star}))=\det\left(\mathbf{L}^{(0)}(x)\right)$ with

[TABLE]

To make all complex constants vanish, we divide the $j$th column by $(-\iota)^{3-j}$ for $j=1,\ldots,3$ and then multiply the $i$th row by $(-\iota)^{3-i}$ for $i=1,\ldots,3$, so that $\det(\mathbf{L}^{(0)}(x))=\det(\mathbf{\tilde{L}}^{(0)}(x))$ with

[TABLE]

Notice that we can simplify the determinant by subtracting the second row scaled by $(4x)^{-1/2}$ from the third row, so that $\det(\mathbf{L}^{(0)}(x))=\det(\mathbf{\tilde{L}}^{(1)}(x))$ with

[TABLE]

We then evaluate the $n\rightarrow\infty$ limit of (114) as

[TABLE]

where we have used the limit definition of the exponential. Since $I_{m}(0)=0$ for all $m\in\mathbb{Z}^{+}$ and $I_{0}(0)=1$, explicit computation of the denominator determinant in (118) gives $1$. Finally, noting that

[TABLE]

since $I_{m}(z)=I_{-m}(z)$ for all $m\in\mathbb{Z}$, we obtain the result (26) for $\alpha_{1}=3$ and $\alpha_{2}=0$.

5.2 Proof for arbitrary $\alpha_{1}$ and $\alpha_{2}=0$

For arbitrary $\alpha_{1}$ and $\alpha_{2}=0$, we use the same approach as in the case $\alpha_{1}=3$ and $\alpha_{2}=0$. First, we make row operations to $\mathbf{E}_{\alpha_{1}}(y)$ to successively reduce the order in $n$ of entries of the same column, similar to (112). Then, we perform some manipulations to facilitate the application of Corollary 3. Finally, after taking $n\rightarrow\infty$, we perform some row operations to simplify the entries of the limiting determinant.

In each iteration, for specified $k$ values, we successively scale the $k$th last row of $\mathbf{E}_{\alpha_{1}}(y)$ by $y$ and then subtract the $(k-1)$th last row to facilitate the application of the Legendre recurrence properties to the entries of the $k$th last row. In the first iteration, we perform those row operations for $k=1,2\ldots,\alpha_{1}-1$. This does not alter the first row. Let $\tau=n+\alpha_{1}-1$. This procedure yields $\det(\mathbf{E}_{\alpha_{1}}(y))=y^{-\alpha_{1}+1}\det(\mathbf{\Theta}^{(1)}_{n}(y))$ with

[TABLE]

where the modified entries of the first column followed from the Rodrigues’ formula (39) and Lemma 1, and the rest of modified entries followed from Corollary 2.

For the sake of notational simplicity, we unify the remaining iterative operations by only employing Lemma 1 to simplify the modified entries. As noted in the previous subsection, Corollary 2 allowed to better illustrate the purpose of each iteration. Here, using Corollary 2 produces cumbersome notation, and we will resort to the more general Lemma 1. Then, we first apply the Rodrigues’ formula (39) to the modified entries beyond the first column of $\mathbf{\Theta}_{n}^{(1)}(y)$ to obtain

[TABLE]

and we repeat the same row operations, but this time for ${k=1,\ldots,\alpha_{1}-2}$, so that $\det(\mathbf{E}_{\alpha_{1}}(y))=y^{-2\alpha_{1}+3}\det(\mathbf{\Theta}^{(2)}_{n}(y))$ with

[TABLE]

where

[TABLE]

with the first term of (123) following from Lemma 1. Like in the case $\alpha_{1}=3$ and $\alpha_{2}=0$, the second iteration gives an additional term for the entries below the second row. In the third iteration, we repeat the same procedure, but this time for $k=1,\ldots,\alpha_{1}-3$, obtaining $\det(\mathbf{E}_{\alpha_{1}}(y))=y^{-3\alpha_{1}+6}\det(\mathbf{\Theta}^{(3)}_{n}(y))$ with

[TABLE]

where

[TABLE]

In the fourth iteration, we repeat the same steps, but this time for $k=1,\ldots,\alpha_{1}-4$, to obtain $\det(\mathbf{E}_{\alpha_{1}}(y))=y^{-4\alpha_{1}+10}\det(\mathbf{\Theta}^{(4)}_{n}(y))$ with

[TABLE]

where

[TABLE]

Observe that an additional term appears every two rows. After a total of $\alpha_{1}-1$ iterations, we obtain $\det(\mathbf{E}_{\alpha_{1}}(y))=y^{-\alpha_{1}(\alpha_{1}-1)/2}\det(\mathbf{\Theta}^{(\alpha_{1}-1)}_{n}(y))$ with

[TABLE]

where

[TABLE]

with fixed $c_{k}^{(i)}\in\mathbb{Z}^{+}$, and $c_{0}^{(i)}=1$, for all $k,i\in\mathbb{Z}$. We then apply Lemma 2 to the entries below the main diagonal, except for the terms generated by the sum in the second line of (129) when $j+k\geq i+2$, since they can be written in terms of derivatives of Legendre polynomials thanks to the Rodrigues’ formula (39). For the main diagonal and the upper-triangular entries, we also apply (39), so that we obtain

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

As in the case $\alpha_{1}=3$ and $\alpha_{2}=0$, we further manipulate the entries to facilitate the application of Corollary 3 and to make the entries of order $1$. We divide the $j$th column by $n^{j-1}(1-y^{2})^{(\alpha_{1}-j)/2}$ for $j=1,\ldots,\alpha_{1}$, while multiplying the $i$th row by $(1-y^{2})^{(\alpha_{1}-i)/2}$ for $i=1,\ldots,\alpha_{1}$. This gives $\det(\mathbf{E}_{\alpha_{1}}(y))=(y/n)^{-\alpha_{1}(\alpha_{1}-1)/2}\det({\mathbf{\tilde{\Theta}}}^{(\alpha_{1}-1)}_{n}(y))$ with

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

Recall $x^{\star}=(1+x/n^{2})/(1-x/n^{2})$. We can now rewrite (97) as

[TABLE]

and take the $n\rightarrow\infty$ limit. Specifically, applying Corollary 3 to the entries of $\mathbf{\tilde{\Theta}}^{(\alpha_{1}-1)}_{n}(x^{\star})$, and recalling that $1-y^{2}=-4x/n^{2}$, we obtain $\lim_{n\rightarrow\infty}\det({\mathbf{\tilde{\Theta}}}^{(\alpha_{1}-1)}_{n}(x^{\star}))=\det(\mathbf{L}^{(0)}(x))$ with

[TABLE]

where

[TABLE]

since $I_{m}(z)=I_{-m}(z)$ for all $m\in\mathbb{N}$.

Now, we divide the $j$th column by $(-\iota)^{\alpha_{1}-j}$ for $j=1,\ldots,\alpha_{1}$ and multiply the $i$th row by $(-\iota)^{\alpha_{1}-i}$ for $i=1,\ldots,\alpha_{1}$, so that all complex constants vanish, just as in the previous subsection. We also perform row operations to get rid of the sums in (140)-(142), recalling that $I_{m}(z)=I_{-m}(z)$ for all $m\in\mathbb{N}$. In the first iteration, we manipulate the third and fourth rows. We subtract the second row scaled by $(4x)^{-1/2}$ from the third row. Then, we subtract the third row scaled by $(4x)^{-1/2}$ from the fourth row. In the second iteration, we manipulate the fifth and the sixth rows similarly. We repeat this procedure for a total of $\lfloor(\alpha_{1}-1)/2\rfloor$ iterations, so that we can write the $n\rightarrow\infty$ limit of (138) as

[TABLE]

where we have used the limit definition of the exponential and

[TABLE]

Considering that $I_{k}(0)=0$ for $k\in\mathbb{Z}^{+}$ and $I_{0}(0)=1$, explicit computation of the denominator in (143) gives $1$. Hence, we have proved the result for arbitrary $\alpha_{1}$ and $\alpha_{2}=0$.

5.3 Extension for arbitrary $\alpha_{2}$

Now consider the case $\alpha_{2}\neq 0$. Substituting (13) and (40) for the entries of the determinants $\mathbf{E}_{\alpha}(y)$ in (11) along with some algebraic simplifications, we obtain

[TABLE]

with

[TABLE]

where $\mathbf{A}^{(0)}(y)$ is a $(\alpha_{1}+\alpha_{2})\times\alpha_{1}$ matrix with entries

[TABLE]

and $\mathbf{B}^{(0)}$ is a $(\alpha_{1}+\alpha_{2})\times\alpha_{2}$ matrix with entries

[TABLE]

In the following, we apply a set of iterative operations to $\mathbf{\Xi}^{(0)}_{n}(y)$ to reduce the case $\alpha_{2}\neq 0$ to the case $\alpha_{2}=0$ when $n\rightarrow\infty$. In the first iteration, we take advantage of the fact that the entries of the $(\alpha_{1}+1)$th column of $\mathbf{\Xi}^{(0)}_{n}(y)$ are all ones to reduce the dimension of the determinant by one. We successively subtract the $(k+1)$th row of $\mathbf{\Xi}^{(0)}_{n}(y)$ from the $k$th row, for $k=1,\ldots,\alpha_{1}+\alpha_{2}-1$, to make zero all the entries of the $(\alpha_{1}+1)$th column of $\mathbf{\Xi}^{(0)}_{n}(y)$ except that of the last row. We then simplify the modified entries beyond the $(\alpha_{1}+1)$th column with the help of Corollary 2 and the second line of (13), i.e.,

[TABLE]

After this set of operations, the entries of that $(\alpha_{1}+1)$th column become all zero, except for that in the last row, which is not altered. We expand the determinant along this column to obtain $\det(\mathbf{\Xi}^{(0)}_{n}(y))=2^{-\alpha_{2}+2}\alpha_{2}^{-1}(\alpha_{2}-1)^{-1}\det(\mathbf{\Xi}^{(1)}_{n}(y))$ with

[TABLE]

where $\mathbf{A}^{(1)}(y)$ is a $(\alpha_{1}+\alpha_{2}-1)\times\alpha_{1}$ matrix with entries

[TABLE]

and $\mathbf{B}^{(1)}$ is a $(\alpha_{1}+\alpha_{2}-1)\times(\alpha_{2}-1)$ matrix with entries

[TABLE]

In the second iteration, we successively subtract the $(k+1)$th row of $\mathbf{\Xi}^{(1)}_{n}(y)$ scaled by $(n+k)/(n+k+1)$ from the $k$th row, for $k=1,\ldots,\alpha_{1}+\alpha_{2}-2$. Then, the entries of the $(\alpha_{1}+1)$th column of $\mathbf{\Xi}^{(1)}_{n}(y)$ become all zeros except for that of the last row, which remains unchanged. We expand along this column to obtain $\det(\mathbf{\Xi}^{(0)}_{n}(y))=2^{-\alpha_{2}+2}((n+\alpha_{1}+\alpha_{2}-1)\alpha_{2})^{-1}(\alpha_{2}-1)^{-1}\det(\mathbf{\Xi}^{(2)}_{n}(y))$ with

[TABLE]

where $\mathbf{A}^{(2)}(y)$ is a $(\alpha_{1}+\alpha_{2}-2)\times\alpha_{1}$ matrix with entries

[TABLE]

and $\mathbf{B}^{(2)}$ is a $(\alpha_{1}+\alpha_{2}-2)\times(\alpha_{2}-2)$ matrix with entries

[TABLE]

In the following iterations, we repeat the same steps as in the second iteration, where we modify the scalings of the rows appropriately to make zero all the entries of the $(\alpha_{1}+1)$th column of $\mathbf{\Xi}^{(r)}_{n}(y)$ except for that of the last row. After a total of $\alpha_{2}$ iterations, we rewrite (145) as

[TABLE]

where $\mathbf{\Xi}^{(\alpha_{2})}_{n}(y)$ is a $\alpha_{1}\times\alpha_{1}$ matrix with entries

[TABLE]

with $s_{k}^{(\alpha_{2})}(n)$ a ratio of polynomials in $n$ of the same order that does not depend on $j$. We do not need to explicitly define $s_{k}^{(\alpha_{2})}(n)$ to complete the proof, since $s_{k}^{(\alpha_{2})}(n)$ does not depend on $j$ and $\lim_{n\rightarrow\infty}\frac{d^{j-1}}{dy^{j-1}}P_{n+k+i-1}(y)$ does not depend on $k$ by virtue of Remark 2. Additionally, $s_{k}^{(\alpha_{2})}(n)=\mathcal{O}(1)$ for all $k$, since $s_{k}^{(\alpha_{2})}(n)$ is the result of multiplying and dividing entries of the $(\alpha_{1}+1)$th column of $\mathbf{\Xi}^{(r)}_{n}(y)$, which have the same order in $n$ for all $r$, and are exactly the same for $\mathbf{\Xi}^{(r)}_{n}(x^{\star})$ and $\mathbf{\Xi}^{(r)}_{n}(1)$. Recalling Corollary 3 and Remark 2, we then have

[TABLE]

where

[TABLE]

which is $\mathcal{O}(1)$. Therefore,

[TABLE]

We also have

[TABLE]

This yields

[TABLE]

which concludes the proof.

6 Proof of Proposition 2

We here prove the Proposition result for the different cases considered. We will only prove the result for the smallest eigenvalue distribution of $\mathbf{W}$ since, once this is established, the result for the largest eigenvalue follows immediately from Remark 1.

6.1 Case $\alpha_{1}=0$ and arbitrary $\alpha_{2}$

The case $\alpha_{1}=0$ is straightforward. From Corollary 1, we have

[TABLE]

and the result follows upon noting that, for fixed $\beta,\gamma$ (independent of $n$),

[TABLE]

6.2 Case $\alpha_{1}=1$ and arbitrary $\alpha_{2}$

The case $\alpha_{1}=1$ is more complicated, and for this we apply some results from Section 5.3; for consistency, we will use the same notation as in that section. Specifically, we will use (156) to write the distribution of the scaled smallest eigenvalue as

[TABLE]

along with the large-$n$ behaviour (including finite-$n$ corrections) of Legendre polynomials (Lemma 3) and that of the coefficients $s_{k}^{(\alpha_{2})}(n)$.

Since computing $s_{k}^{(\alpha_{2})}$ for arbitrary $\alpha_{2}$ is not easy, let us first explicitly compute them for the example case $\alpha_{1}=1$ and $\alpha_{2}=3$, in order to shed light on the properties required to prove the general result. To obtain (165), we start from (145) where, for $\alpha_{1}=1$ and $\alpha_{2}=3$,

[TABLE]

with

[TABLE]

Following the steps described in Section 5.3, we apply a set of iterative operations to $\det(\mathbf{\Xi}^{(0)}_{n}(y))$, reducing successively the matrix dimension. In the first iteration, we obtain $\det(\mathbf{\Xi}^{(0)}_{n}(y))=12^{-1}\det(\mathbf{\Xi}^{(1)}_{n}(y))$ with

[TABLE]

in the second iteration, we obtain $\det(\mathbf{\Xi}^{(0)}_{n}(y))=(12(n+3))^{-1}\det(\mathbf{\Xi}^{(2)}_{n}(y))$ with

[TABLE]

while after the last iteration, we obtain $\det(\mathbf{\Xi}^{(0)}_{n}(y))=(12(n+3)(n+2)(2n+5))^{-1}\det(\mathbf{\Xi}^{(3)}_{n}(y))$ with

[TABLE]

Observing each element of the first column of matrices $\mathbf{\Xi}_{n}^{(t)}(y)$, $t=0,1,2,3$, it is clear that the number of Legendre polynomial terms is doubled in each iteration, if one does not aggregate polynomials of the same degree. After aggregating polynomials of the same degree, we write

[TABLE]

in agreement with (165), where

[TABLE]

Note that these polynomial ratios can be expanded as $1+\mathcal{O}\left(1/n\right)$ since, for $i\in\mathbb{N}$,

[TABLE]

From these expansions, we can see that

[TABLE]

and we generally write

[TABLE]

where $a_{k}^{(3)}$ is equal to the number of aggregated polynomial terms of degree $n+k$. It then follows that $\sum_{k=0}^{3}a_{k}^{(3)}$ is the total number of aggregated polynomial terms; indeed, we see that $\sum_{k=0}^{3}a_{k}^{(3)}=2^{3}$, consistent with the fact that the number of terms is doubled in each iteration.

Applying Lemma 3 to (171), we obtain

[TABLE]

where

[TABLE]

Note also that evaluating (182) at $y=\frac{1+x/n^{2}}{1-x/n^{2}}=1$ (i.e., at $x=0$) yields

[TABLE]

Therefore, we evaluate (166) (equivalently (165)) as

[TABLE]

where the first equality follows from (182), (186) and (164). The Proposition result (for $\alpha_{1}=1$ and $\alpha_{2}=3$) is obtained after noting that

[TABLE]

From the second equality of (187), observe that, when considering the asymptotic distribution to order $\mathcal{O}\left(\frac{1}{n}\right)$, the quantity $\bar{b}^{(3)}$ has no effect (i.e., it drops out in the analysis). This will also occur in the general case of arbitrary $\alpha_{2}$, where we will only need to determine $\bar{a}^{(\alpha_{2})}$ and $\bar{c}^{(\alpha_{2})}$, which depend only on $a_{k}^{(\alpha_{2})}$.

Indeed, following the same steps of the previous example, we find that $\mathbf{\Xi}_{n}^{(\alpha_{2})}(y)=\sum_{k=0}^{\alpha_{2}}s_{k}^{(\alpha_{2})}(n)P_{n+k}(y)$ with $s_{k}^{(\alpha_{2})}(n)=a_{k}^{(\alpha_{2})}-b_{k}^{(\alpha_{2})}/n+\mathcal{O}(1/n^{2})$ and, using Lemma 3 and (164) in (165) we have, for arbitrary $\alpha_{2}$,

[TABLE]

where

[TABLE]

In light of (188), the proof will be complete if

[TABLE]

holds. Let us now prove this equality.

From the previous example, we see that $a_{k}^{(\alpha_{2})}$ equals the number of aggregated polynomial terms of degree $n+k$, after $\alpha_{2}$ iterations, and $\bar{a}^{(\alpha_{2})}$ coincides with the total number of aggregated polynomial terms, irrespective of their degree (i.e., terms with the same degree are counted as different terms). We have also seen that the number of terms doubles after each iteration and, therefore,

[TABLE]

Furthermore, we have seen that, in each iteration, the additional polynomial terms have degrees increased by one (with respect to those terms in the previous iteration). This can be seen in the $(1,1)$ element of the matrices $\mathbf{\Xi}_{n}^{(t)}(y)$, $t=0,1,\ldots,\alpha_{2}$; see (167)-(170) in the $\alpha_{2}=3$ example. From this, it becomes clear that

[TABLE]

and, using this recursion, we can write

[TABLE]

where we have used the facts that $a_{-1}^{(t)}=0$ and $a_{t+1}^{(t)}=0$. From this, along with (190) and (192), it is then clear that

[TABLE]

Thanks to (196), we now prove (191) by induction. For $\alpha_{1}=1$ and $\alpha_{2}=1$,

[TABLE]

where we identify that $s_{0}^{(1)}(n)=s_{1}^{(1)}(n)=1$ in light of (165). Then, $a_{0}^{(1)}=a_{1}^{(1)}=1$, which gives $\bar{c}^{(1)}=4$ and $\bar{a}^{(1)}$=2. Thus, $\bar{c}^{(\alpha_{2})}/\bar{a}^{(\alpha_{2})}=1+\alpha_{2}$ holds for $\alpha_{2}=1$.

For $\alpha_{1}=1$ and $\alpha_{2}=t$, suppose that

[TABLE]

Then, for $\alpha_{2}=t+1$,

[TABLE]

where (196) has been used. This proves (191) and the result of Proposition 2 is obtained after applying (191) and (188) to (189).

6.3 Case $\alpha_{1}=2$ and $\alpha_{2}=0$

For this case, the scaled smallest eigenvalue distribution is given by

[TABLE]

Performing the same row operations as in Section 5.1, we rewrite (200) in the form of (114), i.e.

[TABLE]

By virtue of Lemma 3 and (164), we obtain

[TABLE]

since

[TABLE]

Noticing that

[TABLE]

we develop the numerator determinant in (202) and we obtain the result after some manipulations.

6.4 Case $\alpha_{1}=2$ and $\alpha_{2}=1$

For this case, the scaled smallest eigenvalue distribution is given by

[TABLE]

Performing the same row operations as in Section 5.3, we rewrite (205) in the form of (156), i.e.,

[TABLE]

where

[TABLE]

Now, we perform the manipulations described in Section 5.1, and we write

[TABLE]

where

[TABLE]

When applying Lemma 3 to (209), we obtain

[TABLE]

We then apply some row operations to (210) so that the determinant remains unaltered. Specifically, we scale the second row by $2\sqrt{x}/n$ and we subtract it from the first row. Then, we scale the first row by $4\sqrt{x}/{n}$ and we subtract it from the second row. Therefore, we evaluate (208) as

[TABLE]

where (164) has been used. Since

[TABLE]

we have the result when developing the numerator determinant in (211), performing some simplifications and applying (188).

6.5 Case $\alpha_{1}=\alpha_{2}=2$

For this case, the scaled smallest eigenvalue distribution is given by

[TABLE]

Performing the same row operations as in Section 5.3, we rewrite (213) in the form of (156), i.e.,

[TABLE]

where $\mathbf{\Xi}^{(2)}_{n}(y)$ is a $2\times 2$ matrix with entries

[TABLE]

As in Section 5.1, we perform some row operations to $\mathbf{\Xi}^{(2)}_{n}(y)$ to facilitate the application of Lemma 3. Specifically, we scale the second row of $\mathbf{\Xi}^{(2)}_{n}(y)$ by $y$ and we subtract the first row. We also divide the $j$th column by $n^{j-1}(1-y^{2})^{(2-j)/2}$ for $j=1,2$, and multiply the $i$th row by $(1-y^{2})^{(2-i)/2}$ for $i=1,2$. We then rewrite (214) as

[TABLE]

where

[TABLE]

with

[TABLE]

We then apply Lemma 3 to the entries of $\mathbf{\tilde{\Xi}}^{(2)}_{n}(y)$, and expand the ratio of polynomials in $n$ to obtain, after aggregating terms,

[TABLE]

Using these asymptotic expansions for the determinant entries in (216), we then compute those determinants, make some simplifications and obtain the result with the help of (164) and (188), similarly as in the previous cases.