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

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.


Introduction
Double Wishart random matrices, defined as W = (W 1 + W 2 ) −1 W 1 , with W 1 and W 2 n × 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 W, either the complete set or simply the extreme largest/smallest eigenvalues [1]. For   [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 → ∞ with α 1 = m 1 − n and α 2 = m 2 − n fixed [9,12,13]. Consequently, under hard-edge asymptotics, it follows that the extreme eigenvalue distributions of W can be expressed 25 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 W was shown to be expressible in terms of a Bessel hypergeometric function of a α 1 -dimensional matrix argument. 30 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 W under hard-edge asymptotics. Remarkably, in the context of the LUE , this asymptotic distribution was shown to admit a very 35 simple representation in terms of a finite-dimensional determinant involving Bessel functions [16]. Such a representation should also apply for the extreme eigenvalues of 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 → ∞ and either α 1 → ∞ or α 2 → ∞. Under this regime, the Jacobi kernel was shown to converge to the 2. Fredholm determinant [5,6] 3. Gauss hypergeometric function of α 1 or α 2 -dimensional argument [9] 4. Painlevé VI [7] 5. Polynomial of degree nα 1 or nα 2 involving partitions [9,10] Airy kernel [9]; thus, the extreme eigenvalue distributions of 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 45 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 50 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 β-ensemble, which subsumes the real (β = 1), complex (β = 2) and symplectic (β = 4) cases, and provide a first-order correction proportional to 55 n −2 for β = 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 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 W, and 60 therefore, for the F model and the classical JUE, under hard-edge asymptotics. By exploiting new exact representations of the extreme eigenvalue distributions of W, we perform an asymptotic analysis under the hard-edge scaling regime. In the process, we unveil a striking connection of these exact distributions, classi-cally 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 W can be expressed in terms of α 1 -and α 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 W, giving insights 70 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.

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 α and β, admits [20, eq. (8.960.1)] Jacobi polynomials are orthogonal with respect to the weight w( where δ kl denotes the Kronecker delta function, which equals 1 if l = k and 0 otherwise.

75
The Legendre polynomial of degree l admits [20, eq. (8.910.2)] Legendre polynomials are particular cases of Jacobi polynomials when α = β = 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)] for z ∈ C.

Models
Let X 1 ∈ C m1×n (m 1 ≥ n) and X 2 ∈ C m2×n (m 2 ≥ n) be independent complex Gaussian matrices with independent columns that have the same covariance matrix Σ. Then, W 1 = X † 1 X 1 and W 2 = X † 2 X 2 are n × n complex Wishart matrices with m 1 and m 2 degrees of freedom respectively; i.e., W 1 ∼ CW n (m 1 , Σ) and W 2 ∼ CW n (m 2 , Σ). Define α 1 = m 1 − n, α 2 = m 2 − n. The joint probability density function (JPDF) of the eigenvalues of is proportional to [6] n k=1 φ α1 with 1 ≥ φ 1 > . . . > φ n ≥ 0. This JPDF does not depend on Σ since the eigenvalues of W do not change [6]. The ensemble of n × n matrices of the form (6) is said to have the multivariate complex beta distribution with parameters α 1 and 80 α 2 [22]. The JPDF of its eigenvalues is related to that of other well-known ensembles as follows.
In this work, we first study the extreme eigenvalue distributions of 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 → ∞ with α 1 and α 2 fixed. By simply applying the corresponding transformations aforementioned, our results for W can immediately be rephrased for the JUE and F models.

90
Before presenting our main results, we make note of the following: Hence, from the smallest eigenvalue distribution of W, we can deduce that of the largest eigenvalue by simply applying the transformation ξ → 1 − ξ and interchanging α 1 with α 2 [22].

95
Theorem 1. The cumulative distribution function of the smallest eigenvalue φ n of W admits F φn (ξ) = g α1,α2 (ξ) (9) and that of the largest eigenvalue φ 1 admits with E γ (y) the (α + β) × γ matrix with entries The entries in (12) admit the explicit representation Theorem 1 reveals a tight connection between the distributions of the extreme eigenvalues of W and Legendre polynomials, which are simple particular cases of the Jacobi and Gegenbauer polynomials [24]. The derived exact expressions involve (α 1 + α 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 α 1 and α 2 are kept fixed. This 100 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 Wronskiantype (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 105 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 α 1 = 0 or α 2 = 0.
while when α 2 = 0, while when α 1 = 0, 6 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 110 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 α 1 -dimensional determinant, while that for the largest eigenvalue involves an α 2 -dimensional determinant; in both cases the entries involve relatively complicated linear combinations of derivatives of Jacobi polynomials. The result is as follows: The cumulative distribution function of φ n also admits F φn (ξ) = h α1,α2 (ξ) (18) and that of φ 1 also admits where with G the α × α matrix with entries Since the kth 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 It is noteworthy that the simplified special cases (14) and (16) are also easily recoverable from this 115 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.

Asymptotic extreme eigenvalue distributions of W 120
We consider the hard-edge scaling limit, for which n grows large, with α 1 and α 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 7 smallest eigenvalue of the LUE, which has been shown to admit a remarkably simple form involving a finitedimensional determinant whose entries are Bessel functions [16]. Using this correspondence, along with the simple mapping between the eigenvalues of W and the JUE, given in Section 1.2, an analogous determinant expression should be obtained for the extreme eigenvalue distributions of 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 W. This is accomplished by 130 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 135 α 1 and α 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 while, for the standard deviation, σ φn = n(n + α 2 ) (1 + n(n + α 2 )) 2 (2 + n(n + α 2 )) It is therefore natural to scale φ 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 we arrive at the following: Theorem 2. For fixed α 1 and α 2 , and for x ≥ 0.
8 Contrasting this result with Theorem 1, where the exact extreme eigenvalue distributions of W were given in terms of (α 1 + α 2 )-dimensional determinants, the asymptotic distributions in Theorem 2 involve α 1 -or α 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 α 1 or α 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.

145
If one now considers the JUE, for whichφ k = 2φ k − 1 (see Section 1.2), we easily establish that and therefore These asymptotic distributions coincide precisely with the smallest eigenvalue distribution of the LUE under similar hard-edge scaling (with suitable parameterization of α 1 and α 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 150 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 α 1 or α 2 , as we present below. A generalized formula for arbitrary α 1 , α 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] f (α) Our main result is the following:
Proposition 2 reveals that, for the cases of α 1 ,α 2 considered, the first-order correction of the extreme eigenvalue distributions of W is proportional to the density (30) which, similar to the asymptotic distributions of Theorem 2, is given as an α 1 -or α 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 (28), Proposition 2 can be immediately rephrased for the JUE. Focusing in particular on the smallest eigenvalue, and for the cases of α 1 ,α 2 considered in the proposition, 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 α, the distribution of the smallest eigenvalue for large (but finite) n is given by [25,Theorem 4 which coincides with that of the JUE when α 1 = α and α 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 α 1 , α 2 considered.

160
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 α 1 and α 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 165 validates Proposition 2, and checks numerically whether the stated first-order corrections may continue to hold beyond the cases of α 1 , α 2 considered in Proposition 2.
We first computed the empirical density of the smallest eigenvalue of W for the cases α 1 = 2 and ∞ (x), scaled 1 by ne x , along with the theoretical first-order correction to the asymptotic density, obtained from (31) in Proposition 2 (and correspondingly scaled) as which gives, for α 1 = 2 and arbitrary α 2 , 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. suggesting that Proposition 2 may hold in general. This is formally conjectured as follows: 175 Conjecture 1. For arbitrary α 1 and α 2 , 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 α 1 = 0 or α 2 = 0 when respectively considering the largest or the smallest eigenvalue of the JUE. When both α 1 and α 2 are non-zero, the 180 suggested universality of the first-order corrections under hard-edge scaling does not persist, due to the scaling factor (α 1 + α 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, 185 and remains an interesting topic for future investigation.

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.

Additional Definitions
The Legendre polynomial P n (x) is alternatively defined by the Rodrigues' formula [21, eq. (22.11.5)] The associated Legendre polynomial P −m n (x) is defined by the Rodrigues' formula where n ≥ m.
The shifted Legendre polynomial of degree n is defined by
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.
Proof. First, we prove (54) by following the strategy of [26,Section IV]. Let n > m. Using (39), we rewrite 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., or, equivalently, Let we have Applying (40), we obtain Now, we are ready to make n large. Since for k = 0, 1, 2, . . ., and we obtain that where we have only presented the (k + 2)th term of (62) with Therefore, which leads to the result (54) with the help of (5).
The second term of the expansion (54) depends on a modified Bessel function of one order less than 205 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 knowledge 2 .
Corollary 3. For fixed m ≥ 0 and x > 0, For fixed m ∈ Z and x ≥ 0, Remark 2. Although Corollary 3 presents asymptotic results for Legendre polynomials of degree n when 210 n → ∞, they are also valid for Legendre polynomials of degree n + c when n → ∞, with fixed c ∈ Z, as shown in Lemma 3. We will use this in the proof of Theorem 2.

Proof of Theorem 1
We make use of the following result: . . .
where K is a normalization constant and π n is the nth order monic polynomial orthogonal with respect to 215 the weight function w(y) in the integration interval. We start by writing where 0 ≤ ξ < 1 and .
with C the normalization constant of the eigenvalue JPDF in (7). Since the integrand is symmetric in x 1 , . . . , x n , we may write After the multiple change of variables y i = ( Rearranging the expression yields where T α1,α2 n (β, γ) = 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, . . . , α 1 + α 2 , are not distinct. To proceed, first recognize that if t l for all l were distinct, then where C is a normalization constant and P ν (x) is the ν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 .

Proof of Proposition 1
From [9, eq. (3.16)] and [31], where α 1,2 = α 1 + α 2 and 2 F 1 is the Gauss hypergeometric function of scalar argument. Our desired expression in (92) can be evaluated from (93) by taking limits as where the Gauss hypergeometric function of scalar argument can be expressed in terms of Jacobi polynomials as [21, eq. (15.4.6)] To evaluate these limits, we apply [30, Lemma 2] and we have Finally, we obtain the result by applying the Leibniz rule to the entries of the determinant.

Proof of Theorem 2 225
As before, we prove the result (26) for the smallest eigenvalue, with the result (27) then following from Remark 1.
First consider the case α 2 = 0. Applying (40) to the entries of E α1 (y) and some algebraic simplifications, we obtain If one takes n large and applies Corollary 3 to the entries of the numerator determinant of (97), given in (12), we obtain while for the denominator, with (13), we obtain 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 230 20 E α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 , 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 235 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 α 1 = 3, in order to shed light on the set of manipulations required to prove the general result.
When y = x , by virtue of Corollary 3, we identify for large n 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 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(E 3 (y)) = y −2 det(Θ n (y)) with 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.

22
Note that the entries of the third row of Θ (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 , 1 − y 2 = −4x/n 2 . With this in mind, by virtue of Corollary 3, we now have where in contrast to (101), as alluded earlier, the order in n has been successively reduced in the second 245 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 At this point, by virtue of Corollary 3, the columns ofΘ (2) n (x ) are linearly independent when n → ∞. We then rewrite (97) for α 1 = 3 and α 2 = 0 as and take the n → ∞ limit. Specifically, applying Corollary 3 to the entries ofΘ (2) n (x ), and recalling that 1 − y 2 = −4x/n 2 when y = x , we obtain that lim n→∞ det(Θ (2) n (x )) = det L (0) (x) with 23 To make all complex constants vanish, we divide the jth column by (−ι) 3−j for j = 1, . . . , 3 and then multiply the ith row by (−ι) 3−i for i = 1, . . . , 3, so that det(L (0) (x)) = det(L (0) (x)) with Notice that we can simplify the determinant by subtracting the second row scaled by (4x) −1/2 from the third row, so that det(L (0) (x)) = det(L (1) (x)) with We then evaluate the n → ∞ limit of (114) as where we have used the limit definition of the exponential. Since I m (0) = 0 for all m ∈ Z + and I 0 (0) = 1, explicit computation of the denominator determinant in (118) gives 1. Finally, noting that since I m (z) = I −m (z) for all m ∈ Z, we obtain the result (26) for α 1 = 3 and α 2 = 0.

5.2.
Proof for arbitrary α 1 and α 2 = 0 For arbitrary α 1 and α 2 = 0, we use the same approach as in the case α 1 = 3 and α 2 = 0. First, we make row operations to E α1 (y) to successively reduce the order in n of entries of the same column, similar 255 to (112). Then, we perform some manipulations to facilitate the application of Corollary 3. Finally, after taking n → ∞, we perform some row operations to simplify the entries of the limiting determinant.
In each iteration, for specified k values, we successively scale the kth last row of E α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 kth last row. In the first iteration, we perform those row operations for k = 1, 2 . . . , α 1 − 1. This does not alter the first row. Let τ = n + α 1 − 1. This procedure yields det(E α1 (y)) = y −α1+1 det(Θ (1) n (y)) with 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.
with fixed c (i) k ∈ Z + , and c (i) 0 = 1, for all k, i ∈ 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 ≥ i + 2, since they can be written in terms of derivatives of Legendre polynomials thanks to the Rodrigues' formula (39).

260
Now, we divide the jth column by (−ι) α1−j for j = 1, . . . , α 1 and multiply the ith row by (−ι) α1−i for i = 1, . . . , α 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 ∈ 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 (α 1 − 1)/2 iterations, so that we can write the n → ∞ limit of (138) as where we have used the limit definition of the exponential and Considering that I k (0) = 0 for k ∈ Z + and I 0 (0) = 1, explicit computation of the denominator in (143) gives 1. Hence, we have proved the result for arbitrary α 1 and α 2 = 0.

Extension for arbitrary α 2
Now consider the case α 2 = 0. Substituting (13) and (40) for the entries of the determinants E α (y) in (11) along with some algebraic simplifications, we obtain with where A (0) (y) is a (α 1 + α 2 ) × α 1 matrix with entries and B (0) is a (α 1 + α 2 ) × α 2 matrix with entries In the following, we apply a set of iterative operations to Ξ (0) n (y) to reduce the case α 2 = 0 to the case α 2 = 0 when n → ∞. In the first iteration, we take advantage of the fact that the entries of the (α 1 + 1)th column of Ξ (0) n (y) are all ones to reduce the dimension of the determinant by one. We successively subtract the (k + 1)th row of Ξ (0) n (y) from the kth row, for k = 1, . . . , α 1 + α 2 − 1, to make zero all the entries of the (α 1 + 1)th column of Ξ (0) n (y) except that of the last row. We then simplify the modified entries beyond the (α 1 + 1)th column with the help of Corollary 2 and the second line of (13), i.e., After this set of operations, the entries of that (α 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(Ξ and In the second iteration, we successively subtract the (k + 1)th row of Ξ (1) n (y) scaled by (n + k)/(n + k + 1) from the kth row, for k = 1, . . . , α 1 + α 2 − 2. Then, the entries of the (α 1 + 1)th column of Ξ (1) n (y) become all zeros except for that of the last row, which remains unchanged. We expand along this column to obtain where A (2) (y) is a (α 1 + α 2 − 2) × α 1 matrix with entries and B (2) is a (α 1 + α 2 − 2) × (α 2 − 2) matrix with entries 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 (α 1 + 1)th column of Ξ (r) n (y) except for that of the last row. After a total of α 2 iterations, we rewrite (145) as with s which is O(1). Therefore, We also have This yields which concludes the proof.

Proof of Proposition 2 265
We here prove the Proposition result for the different cases considered. We will only prove the result for the smallest eigenvalue distribution of W since, once this is established, the result for the largest eigenvalue follows immediately from Remark 1.

270
The case α 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 along with the large-n behaviour (including finite-n corrections) of Legendre polynomials (Lemma 3) and that of the coefficients s (α2) k (n).
Note that these polynomial ratios can be expanded as 1 + O (1/n) since, for i ∈ N, 2(n + i) + 1 2(n + i + 1) + 1 From these expansions, we can see that and we generally write where a k is the total number of aggregated polynomial terms; indeed, we see that Note also that evaluating (182) at y = 1+x/n 2 1−x/n 2 = 1 (i.e., at x = 0) yields Therefore, we evaluate (166) (equivalently (165)) as where the first equality follows from (182), (186) and (164). The Proposition result (for α 1 = 1 and α 2 = 3) is obtained after noting that From the second equality of (187), observe that, when considering the asymptotic distribution to order O 1 n , the quantityb (3) has no effect (i.e., it drops out in the analysis). This will also occur in the general case of arbitrary α 2 , where we will only need to determineā (α2) andc (α2) , which depend only on a (α2) k . Indeed, following the same steps of the previous example, we find that Ξ In light of (188), the proof will be complete ifc holds. Let us now prove this equality.
From the previous example, we see that a (α2) k equals the number of aggregated polynomial terms of degree n + k, after α 2 iterations, andā (α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, 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 Ξ and, using this recursion, we can writē where we have used the facts that a Thanks to (196), we now prove (191) by induction. For α 1 = 1 and α 2 = 1, where we identify that s 1 (n) = 1 in light of (165). Then, a 1 = 1, which givesc (1) = 4 andā (1) =2. Thus,c (α2) /ā (α2) = 1 + α 2 holds for α 2 = 1.
6.3. Case α 1 = 2 and α 2 = 0 285 For this case, the scaled smallest eigenvalue distribution is given by ( Performing the same row operations as in Section 5.1, we rewrite (200) in the form of (114), i.e.
Noticing that 1 det   1 0 we develop the numerator determinant in (202) and we obtain the result after some manipulations.
When applying Lemma 3 to (209), we obtaiñ We then apply some row operations to (210) so that the determinant remains unaltered. Specifically, we scale the second row by 2 √ x/n and we subtract it from the first row. Then, we scale the first row by 4 √ x/n and we subtract it from the second row. Therefore, we evaluate (208) as where (164) has been used. Since we have the result when developing the numerator determinant in (211), performing some simplifications and applying (188). For this case, the scaled smallest eigenvalue distribution is given by Performing the same row operations as in Section 5.3, we rewrite (213) in the form of (156), i.e., where Ξ (2) n (y) is a 2 × 2 matrix with entries Ξ (2) n (y) = d j−1 d j−1 P n+i−1 (y) + P n+i + n + i n + i + 1 (P n+i (y) + P n+i+1 (y)) .
We then apply Lemma 3 to the entries ofΞ (2) n (y), and expand the ratio of polynomials in n to obtain, after aggregating terms, 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.