Muttalib--Borodin ensembles in random matrix theory --- realisations and correlation functions

Muttalib--Borodin ensembles are characterised by the pair interaction term in the eigenvalue probability density function being of the form $\prod_{1 \le j<k \le N}(\lambda_k - \lambda_j) (\lambda_k^\theta - \lambda_j^\theta)$. We study the Laguerre and Jacobi versions of this model --- so named by the form of the one-body interaction terms --- and show that for $\theta \in \mathbb Z^+$ they can be realised as the eigenvalue PDF of certain random matrices with Gaussian entries. For general $\theta>0$, realisations in terms of the eigenvalue PDF of ensembles involving triangular matrices are given. In the Laguerre case this is a recent result due to Cheliotis, although our derivation is different. We make use of a generalisation of a double contour integral formula for the correlation functions contained in a paper by Adler, van Moerbeke and Wang to analyse the global density (which we also analyse by studying characteristic polynomials), and the hard edge scaled correlation functions. For the global density functional equations for the corresponding resolvents are obtained; solving this gives the moments in terms of Fuss--Catalan numbers (Laguerre case --- a known result) and particular binomial coefficients (Jacobi case). For $\theta \in \mathbb Z^+$ the Laguerre and Jacobi cases are closely related to the squared singular values for products of $\theta$ standard Gaussian random matrices, and truncations of unitary matrices, respectively. At the hard edge the double contour integral formulas provide a double contour integral form of the scaled correlation kernel obtained by Borodin in terms of Wright's Bessel function.


Introduction
Recent studies in random matrix theory [14,29,13,10,20] have drawn renewed attention to the class of eigenvalue probability density functions (PDFs) proportional to These PDFs were proposed by Muttalib [38] in the context of a simplified model of the joint distribution of the transmission eigenvalues for disordered conductors in the metallic regime, and with no time reversal symmetry. The latter is known to have its exact form proportional to [7] N l=1 e −V (λ l ) 1≤j<k≤N (λ k − λ j ) arsinh 2 λ 1/2 k − arsinh 2 λ 1/2 j , λ l > 0, (1.2) where for large λ, V (λ) = N c arsinh 2 λ 1/2 (1 + O(N −1 )), with c = /L, denoting the mean free path length and L the length of the wire. In practice one has 1 c N .
Recalling that arsinh z = log(z + √ z 2 + 1) one sees (1.1) relates to (1.2) in the limit θ → 0 + when we have although this is still only an approximation to the corresponding factor in (1.2). Actually in [38] attention was restricted to θ a positive integer; on this point we remark that the change of variables λ → λ 1/θ (1.4) maps θ to 1/θ in (1.1) at the expense of altering V (λ). Our interest is two special cases of (1.1). The first is when This is referred to as the Laguerre weight, due to its appearance as the weight function in the orthogonality of the Laguerre polynomials in the theory of classical orthogonal polynomials. The choice (1.5), together with the choice of e −V (λ) as a Gaussian or Jacobi weight (for the latter see (1.6) below), was considered in some detail by Borodin [9]. Due to the significant advancement contained in [9], we will refer to the general class of PDFs (1.1) as Muttalib-Borodin ensembles, and the particular choice of weight (1.5) in (1.1) as the Laguerre Muttalib-Borodin ensemble.
In relation to the Laguerre Muttalib-Borodin ensemble, we first realise the special cases c, θ ∈ Z + as a particular class of complex Wishart matrices isolated in [1]. For general parameters we give a new derivation of a recent result of Cheliotis [13] which gives a realisation in terms of a particular class of random upper-triangular matrices, and we furthermore develop working in [1] to generalise this result (Section 2). The differential equation satisfied by the characteristic polynomial of the ensemble under the mapping (1.4) is studied, and we relate this to the resolvent and global density (Section 3). We use results contained in [1], and take inspiration from the recent work [32], to obtain a derivation of the global density directly from a double contour formula for the one-point function using the saddle point method (Section 4). Furthermore the double contour integral form of the correlation kernel given in [1], suitably generalised from integer to real parameters, is used to rederive the hard edge scaled limit known from [9] (Section 5).
Parallel to the analysis of the Laguerre Muttalib-Borodin ensemble, we also undertake an analogous program of study in relation to the Jacobi weight This substituted in (1.1) gives the Jacobi Muttalib-Borodin ensemble. Our realisation and double contour integral formula for the correlation kernel makes essential use of results contained in [1]. In relation to the global density, the resolvent is specified by a nonlinear equation which we solve using the Lagrange inversion formula to deduce that the moments of the global density are given in terms of particular binomial coefficients. A trigonometric parametrisation of the spectral variable is given which allows for the determination of an explicit functional form for the global density. The hard edge scaled limit gives the same double contour integral form as found for the Laguerre case, in keeping with the findings of [9].

Realisations and extensions
2.1. The Laguerre upper-triangular ensemble. We use the term complex Wishart matrix to refer to a random matrix of the form W † W with W containing complex Gaussian entries with mean and standard deviation to be specified. Furthermore, for any M × N matrix A, we denote A m×n as the m × n matrix consisting of the upper-left m×n block of A. We are interested in a particular complex Wishart matrix, due to Adler, van Moerbeke and Wang [1], which is parametrised by non-negative integers α 1 , . . . , α N satisfying to have entries otherwise. All nonzero entries of X are therefore standard complex Gaussians. Moreover, the condition (2.1) implies all entries on and above the diagonal of X are nonzero. Due to its relationship to a certain family of special functions by the same name, the ensemble of matrices X † M ×n X M ×n , (n = 1, . . . , M ) was referred to in [1] as the multiple Laguerre ensemble. Our first result is to specify an uppertriangular matrix obtained from X by a sequence of Householder transformations. Such transformations were introduced into random matrix theory in [45], [43], and furthermore underpin the construction of β-ensembles as formulated in [17].
(2.2) We say that Y is a random matrix in the upper-triangular ensemble, and that Y † Y belongs to the Laguerre upper-triangular ensemble.
Recall that a complex Householder reflection matrix acting on the left of the M × N matrix X has the form where † denotes the operation of complex conjugate and transpose and u is a M × 1 complex column vector with the property that u † · u = 1. This latter requirement implies U † U = I M , so U is unitary. Geometrically U corresponds to a reflection in the complex hyperplane orthogonal to u † .
To prove the proposition, we construct a sequence of M × N random matrices X (0) = X, X (1) , . . . , X (N ) and a sequence of M × M random Householder reflection matrices U (1) , . . . , U (N ) inductively. The matrix U (l) is determined by X (l−1) and X (l+1) = U (l+1) X (l) (l = 0, 1, . . . , N − 1), where X (l) satisfies: (1) the (j, k) entry is zero for k < j ≤ M if 1 ≤ k ≤ l, or k + α k < j ≤ M if l < k ≤ N ; (2) the diagonal (k, k) entry is real positive and its square is in Γ[α k + 1, 1] distribution if k ≤ l; (3) all other entries are in complex standard Gaussian distribution. By the method of construction to be detailed below, the block of X (N ) consisting of the first N rows is such that X j,l+1 with j = l +1, . . . , l +1+α l+1 are in independent complex standard Gaussian distribution. We denote the (1 + α l+1 )-dimensional vector We construct the Householder reflection matrix U (l+1) using u (l+1) , which is a concatenation of an l-dimensional zero vector, an ( if x (l+1) = 0, and is defined simply as (1, 0, . . . , 0) T if x (l+1) = 0. From the definition of U (l+1) , and recalling X (l+1) = U (l+1) X (l) , it is clear that X (l+1) and X (l) are identical in the upper block consisting of the first l rows and the lower block consisting of the last (M − l − 1 − α l+1 ) rows. In the middle block consisting of the remaining 1 + α l+1 rows, the left part consisting of the middle part of the left-most l columns of X (l) are zeros, so they remain zero in X (l+1) . The entries of X (l) in the right part of the middle block consisting of the right-most (N − l − 1) columns are in independent complex standard Gaussian distribution, and they are all independent of u (l+1) . So after the left multiplication by U (l+1) , the entries of X (l+1) in that part of the middle block are also in independent complex standard Gaussian distribution by the rotational invariance of random Gaussian vectors. The (l + 1)-th column of the middle block becomes ( x (l+1) , 0, . . . , 0) T , so its first entry is positive and the square of the first entry is in Γ[α k +1, 1] distribution, while all other entries are zero. Hence each X (l+1) satisfies the required properties, so finishing the proof by induction.
Remark 2.2. Since we are interested in the spectral properties of X † M ×n X M ×n in the multiple Laguerre ensemble of [1], and by Proposition 2.1 they are identical to Y † N ×n Y N ×n with α 1 , . . . , α N restricted to some subset to their domain, we can think of the upper-triangular ensemble as a generalisation of the multiple Laguerre ensemble.
The upper-triangular matrix Y distributed as in (2.2) with has recently been shown by Cheliotis [13] to have the PDF for its squared singular values given by the Laguerre Muttalib-Borodin ensemble (1.1) with weight (1.5). For this to relate to our construction from complex Gaussian matrices according to Proposition 2.1 we must have θ and c non-negative integers. Thus in this circumstance we have identified a realisation of this ensemble as the eigenvalue PDF of a Wishart matrix.
Consider the matrix X as defined below (2.1), and with θ, c ∈ Z ≥0 , let α j be specified as in (2.3). We have that the eigenvalue PDF of the Wishart matrix X † X is given by the Laguerre Muttalib-Borodin ensemble (1.1) with weight (1.5).
Additional details of the spectral properties of X † X, or equivalently according to Proposition 2.1, of Y † Y with the parameters α j of integer values, beyond the joint distribution of the eigenvalues have been given in [1]. In particular, we can read off from the multiple Laguerre part of [1, Thm. 1] the explicit functional form of the conditional distribution of the eigenvalues of Y † N ×n Y N ×n , given the eigenvalues Below we show that the results can be generalised to arbitrary upper-triangular random matrices Y with real-valued α j > −1. The proofs are similar to those in [1] and we mainly emphasis the differences. Proposition 2.4. Let the N × N random matrix Y be in the upper-triangular ensemble with diagonal entries specified by (2.2). Denote by {λ 1 , . . . , λ n } and {µ 1 , . . . , µ n−1 } the eigenvalues of Y † N ×n Y N ×n and Y † N ×(n−1) Y N ×(n−1) respectively in descending order. The conditional PDF of {λ 1 , . . . , λ n }, with {µ 1 , . . . , µ n−1 } fixed and distinct, is equal to , (2.4) subject to the interlacing constraint Actually we can prove a slightly stronger result: . . , µ n−1 } in descending order. Define the n × n matrix B = (b i,j ) by letting: (1) the (n − 1) × (n − 1) upper-triangular block equal A; (2) the bottom row has all entries but the rightmost one equal 0; (3) all entries of the rightmost row be independent random variables, such that the (n, n)-entry b n,n is real positive with |b n,n | 2 d = Γ[α n +1, 1], and all other entries in the row are in standard complex normal distribution. Then the eigenvalues of B † B, denoted by {λ 1 , . . . , λ n } in descending order, satisfies the interlacing constraint (2.5) and have the distribution given by p n,n−1 ({λ j } n j=1 , {µ j } n−1 j=1 ) in (2.4).

Proof of Lemma 2.5.
A key point is that where y = (b 1,n , b 2,n , . . . , b n,n ) T . The distribution of b k,n implies that Moreover there exists an (n − 1) × (n − 1) unitary matrix V depending on A = B (n−1)×(n−1) such that where 1 .
Using the property that the multiplication of a unitary matrix and a vector of independent standard complex Gaussians yields another vector of independent standard complex Gaussians, we have that the vector z = (z 1 , . . . , z n ) T defined as has the properties that all its components are independent, z 1 , . . . , z n−1 are in standard complex normal distribution, and . . , n − 1), and z n = b n,n .
Conjugating both sides of (2.6) as in (2.8), we have that where ePDF denotes the eigenvalue PDF. It is standard computation to show that the eigenvalue equation for the matrix on the RHS of (2.9) is given by The distribution of the roots of this rational function, with residues distributed according to (2.7), and thus the conditional PDF of {λ 1 , . . . , λ n }, can now be read off as a special case of [22,Cor. 3], and (2.4) with interlacing constraint (2.5) follows.
Knowledge of Proposition 2.4 allows us to rederive the result of Cheliotis [13] noted below (2.3). We will require the use of a particular multiple integral evaluation.
Lemma 2.6. Let R λ denote the region (2.5) for (µ 1 , . . . , µ n−1 ), and suppose α k = α n (k = 1, . . . , n − 1). We have Since the dependence on µ j is entirely in row j, the integration over λ j+1 > µ j > λ j can be done row-by-row. Furthermore, by adding row 1, . . . , j − 1 to row j the integration can be taken to be over λ j+1 > µ j > λ 1 . Applying this operation to (2.12) gives as can be seen by subtracting the first row from each of the next rows in the determinant on the RHS, then expanding by the final column to obtain the LHS. Thus (2.11) now follows.
(2. 16) In the special case the α j are given by (2.3) this reduces to The joint probability distribution function of λ (1) , λ (2) , . . . , λ (N ) is equal to where C n is defined in (2.16), and 1 µ λ is the indicator function of the region that satisfies inequality (2.5). In case that some α i are identical, we understand the formulas in the limiting sense with l'Hôpital's rule.
Proof of Part (a). We prove the case that α i are distinct, and the general result follows by analytical continuation.
(2.20) Substituting for the integrand, then making use of Lemma 2.6 shows that the LHS is equal to Comparison with (2.15) and recalling (2.16) shows that this is precisely the RHS.
To deduce (2.17) from (2.15) we make use of the Vandermonde determinant identity Proof of Part (b). By Lemma 2.5, we have that the eigenvalues λ (1) , λ (2) , . . . , λ (N ) constitute an inhomogeneous Markov chain with the transition probability density function from time n − 1 to time n given by (2.4). Thus the joint distribution function of λ (n) (n = 1, . . . , N ) is obtained by multiplying (2.4) repeatedly. The argument is the same as the proof of [1, Cor. 1] and we omit the details.
The proof of the proposition is by a standard argument for determinantal processes based on the joint probability density function (2.19). In [1,Thm. 3(c)], the proposition for non-negative integer α i under condition (2.1) is proved. Since the proof does not use these additional conditions, it is also a complete proof to Proposition 2.10.

2.2.
The Jacobi upper-triangular ensemble. Adler, van Moerbeke and Wang [1] considered the joint eigenvalue PDF of a sequence of random matrices where X M ×n is the left M × n sub-block of the M × N matrix X specified by (2.2), while A M ×n is the top M × n sub-block of the M × N matrix A which has all elements independently distributed as standard complex Gaussians. Here it is required that M ≥ n and M ≥ n. Due to its relationship to particular multiple orthogonal polynomials, this was referred to as the Jacobi-Piñeiro ensemble. We consider in this section a Jacobi-type counterpart of the random matrix ensembles Y † N ×n Y N ×n in Section 2.1. To this end, we use the N ×N random matrix Y specified in Section 2.1 above Proposition 2.1, and denote the N ×N random matrix Z, which is also in the same upper-triangular ensemble as Y . Let Z = [z j,k ] j,k=1,...,N be an upper-triangular N × N random matrix with all upper-triangular entries independent, the diagonal entries be real positive with distributions depending on 24) and all entries strictly above the diagonal be complex and in standard complex Gaussian distribution. Then define for all n = 1, . . . , N 25) where Y N ×n (resp. Z N ×n ) is the left N × n sub-block of the N × N matrix Y (resp. Z). Note that bothJ n and J n depend on parameters α 1 , . . . , α n > −1, J n also depends on β 1 , . . . , β n > −1 and forJ n it is further assumed that α 1 , . . . , α n are non-negative integers satisfying inequality (2.1). We say that the matrices J n form the Jacobi upper-triangular ensemble. Parallel to Proposition 2.1, we have Proof. The proof is analogous to that of Proposition 2.1, and we divide it into two steps. As a bridge betweenJ n and J n , for all n = 1, . . . , N we definê Recall the sequence of random matrices X (0) = X, X (1) , . . . , X (N ) and random Householder matrices U (1) , . . . , U (N ) constructed in the proof of Proposition 2.1. We have that for all n = 1, . . . , N , (X M ×n are independent of A, and they have the same joint distribution as Y † N ×n Y N ×n . So the joint distribution ofĴ n is the same as that of On the other hand, X (N ) is a function of X and (X is identical toJ n given that X and A are the same. Thus we have showed that the joint distribution ofJ n is the same as that ofĴ n . Next, since A is also in the upper-triangular ensemble, we use the same algorithm in the proof of Proposition 2.1 to construct a sequence of random matrices A (0) = A, A (1) , . . . , A (N ) and random Householder matrices V (1) , . . . , V (N ) , such that A (n+1) = V (n+1) A (n) and A (n) satisfy the same conditions as X (n) but with the parameters M, α 1 , . . . , α N replaced by M , M ×n is independent of Y , and they have the same joint distribution as Z † N ×n Z N ×n . So the joint distribution of J n is the same as that of

On the other hand, A (N ) is a function of A and (A
is identical toĴ n given that Y and A are identical. Thus we show that the joint distribution ofĴ n is the same as that of J n .
The following proposition is the Jacobi counterpart of Proposition 2.4, and its proof is analogous to that of the Jacobi-Piñeiro part of [1, Thm. 1].
Proposition 2.12. Let the N × N random matrices Y and Z be in the uppertriangular ensemble with diagonal entries specified by (2.2) and (2.24), and thus with parameters α 1 , . . . , α N and β 1 , . . . , β N respectively. Denote by {λ 1 , . . . , λ n } and {µ 1 , . . . , µ n−1 } the eigenvalues of J n and J n−1 respectively in descending order, where J n and J n−1 are defined in (2.25). The conditional PDF of {λ 1 , . . . , λ n }, with {µ 1 , . . . , µ n−1 } fixed and distinct, is equal to subject to the interlacing constraint Analogous to Proposition 2.4, we actually prove a slightly stronger result: . . , µ n−1 } in descending order. Define the n × n matrix C n (resp. B n ) by letting: (1) the (n − 1) × (n − 1) upper-triangular block equal C n−1 (resp. B n−1 ); (2) the bottom row has all entries but the rightmost one equal 0; (3) all entries of the rightmost column be independent random variables, such that the (n, n)-entry c n,n (resp. b n,n ) is real positive with |c n,n | 2 d = Γ[α n + 1, 1] (resp. |b n,n | 2 d = Γ[β n + 1, 1]), and all other entries in the column are in standard complex normal distribution. Then the eigenvalues of (C † n C n )(B † n B n + C † n C n ) −1 , denoted by {λ 1 , . . . , λ n } in descending order, satisfy the interlacing constraint (2.28) and have distribution given Proof of Lemma 2.13. It is more convenient to consider for * = n or n − 1, and their eigenvalues {λ 1 , . . . ,λ n } for T n and R n , and {μ 1 , . . . ,μ n−1 } for T n−1 and R n−1 , assumed in ascending order, noting that T * and R * have the same eigenvalues. The relations betweenλ i ,μ i and λ i , µ i arẽ It is straightforward to check that the proposition is equivalent to the statement that the PDF of {λ 1 , . . . ,λ n } is equal tõ , (2.29) subject to the interlacing constraint For notational simplicity, we denote y = (c 1,n , c 2,n , . . . , c n−1,n ) T , z = (b 1,n , b 2,n , . . . , b n−1,n ) T , η = c n,n and ζ = b n,n . Then Taking the singular value decomposition to C n−1 B −1 n−1 , we have (n−1)-dimensional unitary matrices U, V such that Analogous to (2.10), the eigenvalue equation for the matrix R n = Q † Q is given by where the w k are components of w. Note that |w 1 | 2 , . . . , |w n−1 | 2 , η 2 , ζ 2 are independent, and their distribution functions are positive with densities Comparing ( This result can be used to derive the analogue of Corollary 2.8 in the Jacobi case.
(a) Then the eigenvalues of the random matrix J n , denoted by λ (n) = {λ n } in descending order, is equal to In the special case that α j is given by (2.3) this reduces to where C n is defined in (2.34), and 1 µ λ is the indicator function of the region that satisfies inequality (2.30). In case that some α i are identical, we understand the formulas in the limiting sense with l'Hôpital's rule.
Proof of Part (a). We prove the case that the α i are distinct, and the general result follows by analytical continuation. Recall the conditional PDF p n, It is a classical result that this combination of random variables is distributed according to the beta distribution B[α 1 + 1, β 1 + 1] and thus has for its PDF in agreement with (2.33) with n = 1. Proceeding by induction, let us now assume that (2.33) is valid in the case n − 1. Our task is to check the validity of which is analogous to (2.20), but the integral domain R J λ for {µ 1 , . . . , µ n−1 } is defined by (2.28), and the p n , p n−1 and p n,n−1 are defined differently. Substituting for the integrand, then making use of Lemma 2.6 shows that the LHS is equal to Remark 2.15. The assumption (2.32) that β k are in arithmetic progression with common difference −1 is crucial in the application of Lemma 2.6. By the symmetry of the model, it is also possible to let β k be arbitrary and α k in arithmetic progression with common difference −1.
The proof of the proposition is by a standard argument of determinantal process based on the joint probability density function (2.37). In [1,Thm. 3(c)], the proposition for non-negative integer β and non-negative integer α i under condition (2.1) is proved. Since the proof does not use these additional conditions, it is also a complete proof to Proposition 2.10. Note that in [1], only case (1) of the contours occurs.
It is possible to use Corollary 2.8(a) to give an alternative derivation of Corollary 2.14(a), in the case that if β 1 , . . . , β N satisfies (2.32) with β ∈ Z + . In this case we note that the eigenvalue PDF of J n is the same as the eigenvalue PDF ofĴ n defined in (2.26), where the height of the random matrix A is M = β + N . We also need a recent result due to Kuijlaars and Stivigny [29]. Below we give the derivation without the tedious calculation of the normalisation constant of the PDF. . Let the matrix W be an n × n random matrix such that W † W has an eigenvalue PDF proportional to ..,n . For ν ≥ 0, let G be an (n + ν) × n random matrix whose entries are in independent standard complex Gaussian distribution. The squared singular values of GW , or equivalently the eigenvalues of (GW ) † GW , have PDF proportional to In the application of Proposition 2.17, we let W = Y −1 n×n and G = A M ×n , such that ν = M − n = β n = β + N − n. In Corollary 2.8 we obtained the eigenvalue PDF of Y † n×n Y n×n . From this, by the change of variables λ j → 1/x j , we have that the eigenvalue PDF of and so we can take Hence we deduce that the eigenvalue PDF of (GW n×n , which is the same as the eigenvalue PDF of The eigenvalues {λ j } ofĴ n and the eigenvalues {σ j } ofŜ n are related by  [29,20] have revealed an intimate relationship between random matrix products and the Laguerre Muttalib-Borodin ensemble. To explain this requires the introduction of a family of integer sequences -the Fuss-Catalan numbers -parametrised by s ∈ R + and specified by For general s > 0 these are known to be moments of a PDF -the Fuss-Catalan density -with compact support [0, L], L > 0 [6,36], and they uniquely define the PDF.
Consider first a product of s N × N matrices X 1 , . . . , X s with each containing independent, identically distributed zero mean, unit standard deviation random variables. Alternatively, for one such matrix X 1 say, consider the power X s 1 . In either case, ask for the limiting spectral density of the squared singular values after dividing by N s -what results is precisely the Fuss-Catalan density with parameter s [3,6,40].
Consider now the Laguerre Muttalib-Borodin ensemble defined by (1.1) and be the eigenvalues in the N -dimensional ensemble. As N → ∞, by standard techniques we have that the empirical distribution of the scaled eigenvalues N −1 λ converges in distribution to a limiting probability distribution, also known as the equilibrium measure of the model; see [16,Sec. 6.4]. The equilibrium measure is characterised as the minimum of a variation problem; see Claeys and Romano [14,Eq. (1.22)]. By interpreting the recent results of [14], Forrester and Liu [20] have identified the global density (i.e. density scaled by an appropriate power of N to have compact support) for the Laguerre Muttalib-Borodin ensemble in terms of the Fuss-Catalan density for general s ≥ 1.
Proposition 3.1. Suppose θ ≥ 1 (this is for technical reasons in the working of [14]; in [20] it is commented that the same result is expected to hold for all θ > 0 and in fact this has recently been established in [21]). After changing variables x In particular, this shows a relationship between the product of s random matrices and the Laguerre Muttalib-Borodin ensemble with θ = s (see also [29] and the appendix in [21]). Here we will demonstrate the relationship in a different way, by considering the characteristic polynomial of the latter after the change of variables (1.4). The proof of Proposition 3.1 via Claeys and Romano's approach is valid for all θ ≥ 1 at least, but depends on the construction of a mapping J(s) [14,Eq. (1.25)], so it is unclear if it can be applied to the Jacobi case. On the other hand, the approach presented below is valid only for θ ∈ Z + , but it does generalise to the Jacobi case.
Crucial for the alternative proof is knowledge of certain biorthogonal polynomials associated with (1.1) where V (x) is given in (1.5). Thus for given j = 0, 1, 2, . . . let p j (x), q j (x) be monic polynomials of degree j, and suppose these polynomials where the positivity of h j follows from the positivity of the integral over ( Since Proposition 3.1 requires the change of variables (1.4), we see that q N (x) is equal to the corresponding averaged characteristic polynomial. For the Laguerre weight (1.5), the explicit form of q N (x) is known from a result of Konhauser [28]. A relation between q N (x) and generalised hypergeometric functions is observed in [44]. We will use the standard notation p F q to denote the generalized hypergeometric function defined by a series as presented in e.g. [34,Sec. 3]. .

7)
where dµ L denotes the equilibrium measure, J is the corresponding support, dμ L is the measure transformed from dµ L by the change of variablex = (x/θ) θ , andJ is the support of dμ L .
The proof of Lemma 3.3 is similar to [16,Lem. 6.77]. Note that in [16] it is required that the function φ, corresponding to the function log(z − (x/θ) θ ), is a bounded function in x, while our log(z − (x/θ) θ ) is not. But since the growth of the function as x → ∞ is mild, the argument there can be applied. Comparing the left-hand side of (3.7) with the formula (3.3) for q N , we have that is the limiting resolvent (or equivalently Stieltjes transform) of the measure dμ.
Below we show thatG L (z) satisfies a polynomial equation which uniquely characterises the Fuss-Catalan distribution.
Proposition 3.4. Transform the Laguerre Muttalib-Borodin ensemble according to (1.4). DefineG L (z) by (3.9) so that it is equal to the resolvent corresponding to the global scaled density, scaling x = (N θ) θ z. For θ ∈ Z + we have that Proof. First we note that the convergence (3.8) and (3.9) yields that as N → ∞, where o(1) is an analytic function in z that vanishes uniformly for z in any compact set of C \ [0, ∞), and C is a constant. The asymptotic formulas below are thus justified. To leading order in N , after substituting for {a i }, {b j } as implied by (3.5), we see that (3.6) reads Now change variables x = (N θ) θ z, and let q (k) N (x) denote the k-th derivative of q N (x). Working contained in [20, above Prop. 5.1], establishes that under the validity of (3.8), (note that (3.8) itself is the case k = 1; the general k case follows by expressing higher derivatives in terms of the logarithmic derivatives). Using this in (3.11) gives (3.10). (ii) The spectral density is the global scaled limit of the one-point correlation function. Let λ be the eigenvalues in the Ndimensional Laguerre Muttalib-Borodin ensemble. Then the one-point correlation function is (we suppress the N -dependence in the notation ρ (1) ) (3.14) where p N is the PDF implied by (1.1) with the Laguerre weight (1.5). The scaled spectral density, as N → ∞, converges to the equilibrium measure in distribution: where dµ L is the same as in (3.7). (iii) The resolventG L (z) defined in (3.9) can be expressed by the limiting spectral density/equilibrium measure of the Laguerre Muttalib-Borodin ensemble asG There is another viewpoint on deducing the global spectral density from knowledge of (3.5). This makes use of a recent result of Hardy [25], which subject to a mild technical condition [25,Eq. (1.12)] states that characteristic polynomials coming from a class of determinantal point processes including the multiple orthogonal polynomial ensemble under present discussion (see [25,Sec. 1.4]), the limiting density of zeros equals the limiting global spectral density. On the other hand, the limiting density of zeros for the generalised hypergeometric function 1 F θ in (3.5) has been shown by Neuschel [39] to be given by the Fuss-Catalan density with parameter s = θ. (Strictly speaking the result of [39] assumes the bottom line of parameters in (3.5) to be positive integers. Since the leading asymptotics are independent of these parameters, it is expected that this assumption in [39] is not necessary.) 3.2. The Jacobi Muttalib-Borodin ensemble. As with the Laguerre weight (1.5), for the Muttalib-Borodin ensemble with Jacobi weight (1.6), we also let λ be the eigenvalues in the N -dimensional ensemble. As N → ∞, by standard techniques, we have that the empirical distribution of the eigenvalues converges in distribution to a limiting probability distribution, also known as the equilibrium measure of this model. The equilibrium measure can be characterized as the minimum of a variation problem analogous to that for the Laguerre case, and from this, in the recent work [21] the moments of the corresponding density have been given in terms of certain binomial coefficients (see Proposition 3.10 below). Analogous to the Laguerre ensemble, there is a corresponding system of biorthogonal polynomials. For j = 0, 1, 2, . . . , let p j (x), q j (x) be monic polynomials of degree j, and with weight V (x) defined in (1.6) suppose these polynomials have the biorthogonal property analogous to (3.2), h j > 0. (3.17) Let P N,θ denote the PDF specified by (1.1) with V specified by (1.6) and 0 < λ l < 1. Then (3.3) also holds in the Jacobi ensemble with corresponding different meanings of p N , q N and P N,θ . Below we state algebraic results on the biorthogonal polynomials analogous to Proposition 3.2. Again we will focus on the polynomials q j (x).
Also we have the analogue of Lemma 3.3.

20)
where dµ J denote the equilibrium measure of the Jacobi ensemble and dμ J denote the measure transformed from dµ J by the change of variablex = x θ . Note that both dµ J and dμ J have support [0, 1].
The proof of Lemma 3.7 is similar to [16,Lem. 6.77], and we omit the details. From Lemma 3.7, we derive the counterpart of (3.8), that is the limiting resolvent of the measure dμ J . Below we show thatG J (z) satisfies a polynomial equation that is similar to (3.10) characterising the Fuss-Catalan distribution.
Proposition 3.8. Transform the Jacobi Muttalib-Borodin ensemble according to (1.4). DefineG J (z) by (3.22) so that it is equal to the resolvent corresponding to the global density. For θ ∈ Z + we have that Proof. Applying the identity (3.6) to (3.19) with j = N , we have, to leading order in N , while the analogue of (3.12) is (ii) Generally we expect the global density associated with the weight (1.6) to be independent of c 1 and c 2 provided those parameters are themselves independent of N . On the other hand, the change of variables (1.4) shows that for finite N the density ρ (1) (x), defined by (3.14) with p N the PDF implied by (1.1) with the Jacobi weight (1.6), has the functional property Proof. Let 1/z = Lt and X = zG J (z). Then by (3.31), X is a power series in 1/z and also a power series in t. Simple manipulation of (3.23) shows which is of the form (3.32) with a = 0. Applying (3.33) with f (ζ) = ζ and recalling (3.31) shows that (3.36) Using the binomial theorem to expand the two main factors on the right-hand side in (3.36) into power series in z, then combining the coefficients appropriately to form a single power series shows (3.37) The sum can be recognised as a polynomial example of a particular 2 F 1 Gaussian hypergeometric function, allowing us to writẽ The functional equation for the gamma function shows Also, in general, it is a simple exercise to verify from the series definition of 2 F 1 that In the case of the 2 F 1 in (3.38) we have 2b − c + 2 + (a − b − 1)z = 0, b + 2 = c. Thus on the RHS, only the second term contributes, and furthermore it can be simplified from the fact that 2 F 1 (a, b + 2; b + 2; z) = (1 − z) −a , implying the result

Substituting (3.39) and (3.40) into (3.38) we obtain (3.34).
Remark 3.11. (i) With 1/z = Lt, the large z expansion of zG J (z) is thus seen to be a special case of the function This is intimately related to Lambert's solution of the trinomial equation x = q + x m in a power series in q [24]. In mathematical physics, there are applications of (3.41), and its multivariable analogue, in the theory of anyons [4,5]. (ii) As z → 0, we read off from (3.23) that (zG J (z)) θ+1 ∼ −z(1/θ) θ . Since the corresponding global density is given in terms of the resolvent by dμ The corresponding leading term behaviour for x → 0 + agrees with that implied by (3.42).
In addition to the above remarks, we draw attention to a relationship between the Jacobi Muttalib-Borodin ensemble and products of truncations of Haar distributed unitary matrices. Let U j be a Haar distributed unitary random matrix of size m j × m j and let T j be the corresponding (n + ν j ) × (n + ν j−1 ), with ν j ≥ 0 and ν 0 = 0, and m j ≥ 2n upper left block such that ≥ 2n. Let G j denote a standard complex Gaussian matrix of size (n + ν j ) × (n + ν j−1 ). According to the recent work [23], in the limit n → ∞ the singular values squared of the random matrix product has a density such that its moments J r,s,1 (n) are given by J r,s,a (n) = a n a r (1 + a) s n P (αn−1,βn−1) n−1 where α n = rn + r + 1, β n = −(r + 1 − s)n − (r + 2 − s) and P (α,β) k (x) denotes the Jacobi polynomial. It is also required that r > s and thus at least one Gaussian matrix in the product. The immediate relevance of this work is due to the fact that the moments (3.44) are shown in [23] to be such that w = azG(z), where G(z) denotes the corresponding resolvent, satisfy w r+1 − z(w − a)(w + 1) s = 0. (3.45) This is the same as our equation (3.23) with w/a = zG J , r = s = a = θ. Indeed for these parameters we can check that (3.44) reduces to (3.34), up to the form of the scale factor L. Subsequent to [23], the work [26] has considered integrability and exactly solvable features of the random matrix product (3.43) with r = s. In particular, it has been shown that the corresponding characteristic polynomial for the squared singular values is given by [26,Eq. (2.31 where G 0,s+1 s+1,s+1 is a particular Meijer G-function and instead of the constraint m j ≥ 2n as required in [23], it is required m 1 ≥ 2n 1 + µ 1 and m j ≥ n + ν j + 1 (j = 2, . . . , n). According to the strategy introduced in [20], and applied in the derivation of Propositions 3.4 and 3.8 above, the significance of this result is that the Meijer G-function satisfies a linear differential equation (see e.g. [34,Sec. 5.8]) allowing us to deduce a polynomial equation for the corresponding resolvent. Specifically, this procedure gives (3.47) where α j = lim n→∞ m j /n, β j = lim n→∞ ν j /n. We see that (3.47) reduces to (3.34) in the case s = θ, β j = 0 (j = 1, . . . , s), α j = 1 + 1/θ (j = 1, . . . , s). Thus we learn that the global spectral density for the Jacobi Muttalib-Borodin ensemble in the case θ = s ∈ Z + is the same as the global spectral density for the squared singular values of the random matrix product T s · · · T 1 where each T j is an n×n sub-block of a Haar distributed unitary matrix of size (n + p) × (n + p) where lim n→∞ p/n = 1/s. Choosing instead β j = 0 (j = 1, . . . , s) and α j = 1 + 1/a we see that (3.47) is the same equation as (3.45) with w/a = zG in the latter. Thus J s,s,a (n) has an interpretation as the moments of the spectral density of the squared singular values of the random matrix product T s · · · T 1 where each T j is an n × n sub-block of a Haar distributed unitary matrix of size (n + p) × (n + p) where lim n→∞ p/n = 1/a for general a > 0. The case a → 0 exhibits further structure. Then the underlying unitary matrices are of infinite size, and after rescaling the sub-blocks have a Gaussian distribution. More specifically, an n×n sub-block U n of an (n+p)× (n + p) Haar distributed unitary matrix has a distribution proportional to (det(I n − U † n U n )) p−n [47]. It follows that for p → ∞, and with n fixed, the distribution of G n = √ pU n is proportional to e −G † n Gn . Hence G n is a standard complex Gaussian matrix. For a → 0 we see from (3.44) that These are Fuss-Catalan numbers (3.1), which we know are moments of global spectral limit for the squared singular values of a product of s standard complex Gaussian matrices.

The global density -saddle point method
It has already been remarked that the Borodin-Muttalib ensemble (1.1) is intimately related to biorthogonal polynomials of one variable as specified by (3.2). Moreover, as noted in [38] (1.1) is an example of a determinantal point process.
This means that the general k-point correlation function ρ (k) (x 1 , . . . , x k ) is fully determined by a function K(x, y) = K N (x, y), independent of k and referred to as the correlation kernel, according to the formula (4.1) Moreover, the correlation kernel is expressed in terms of the biorthogonal polynomials p l , q l defined and corresponding normalisation h l in (3.2) according to Since we focus on the classical Laguerre and Jacobi cases, we denote the correlation kernel by K L (x, y) for V defined in (1.5) and K J (x, y) for V defined in (1.6).
In [9], an indirect way to transform the summation in (4.2) was devised for the classical Laguerre and Jacobi weights. This transformed summation enabled the computation of the so called hard edge scaled limit. This refers to the limit N → ∞ with x, y scaled so that in the neighbourhood of the origin the spacing between eigenvalues is of order unity. Our present interest is in the global limit of the one-point function. The transformed summation of [9] is not suited for that purpose. Fortunately Propositions 2.10 and 2.16 provide the double contour integral formulas for both the Laguerre and Jacobi cases of the Borodin-Muttalib ensemble, which is suited. In Section 2 we showed that they are spectrally equivalent to the upper-triangular ensemble (defined in (2.2)) and its Jacobi counterpart (defined in (2.24)) respectively, but with restricted values of α i and β i .
where α j = θ(j − 1) + c is specified in (2.3), and the contours Σ and Γ α are specified in Proposition 2.10. The factor h(x)/h(y) in (4.3) cancels out of the determinant (4.1); this factor account for the different meaning of K(n 1 , x; n 2 , y) adopted in Proposition 2.10 in consistency with [1,Eq. (21)], which has the factor e −(V (x)+V (y))/2 in (4.2) replaced by e −x y c , see also [1, Eq. (108)] and derivations above it. Thus we have We recall from (2.21) that the upper-triangular ensemble is well defined for θ = 0 and c > −1 and it is the θ → 0 + limit of the Laguerre Muttalib-Borodin ensemble. We call it the θ = 0 case of the ensemble. The double contour integral formula (4.3) is still valid in this case, where α j = c for all j = 1, . . . , N . Note that this case is degenerate in the sense that the biorthogonal polynomials p l and q l are not well defined by (3.2).
Below we consider the limiting global density in two cases, first for θ > 0 and next for θ = 0. We give details of the derivation in the former case, while point out the differences of the argument in the latter. 4.1.1. The θ > 0 case. We seek to use (4.3) to compute the density function of the the transformed limiting counting measure dμ in (3.7). By definition the density function is given by (3.15) and the required change of variables isx = (x/θ) θ (recall the text below (3.7)), so (4.2) gives and is the global density appearing in (3.9). The scale factor θ is chosen with the benefit of hindsight as it allows (3.10) to be reclaimed without the need for rescaling. Under the assumption θ ∈ Z + , we showed in Section 3.1 that the resolventG L (z) defined in (3.9) corresponding to dμ L , satisfies the identity (3.10) characterizing the Fuss-Catalan density. Consequentlỹ is the Fuss-Catalan distribution supported on I = (0, (1+θ) 1+θ θ −θ ). The parametrization by Biane and independently Neuschel [8,39] gives that for x ∈ I, there is a unique ϕ ∈ (0, π/(θ + 1)) such that and which allows for the simple functional form The main result obtained in the subsection is: The remaining part of this section is devoted to the proof of Proposition 4.1. Our analysis, which is based on the method of steepest descents, is guided by a very similar calculation carried out recently by Liu, Wang and Zhang [32] in relation to the global density for the squared singular values of a product M standard complex Gaussian matrices. That these two computations should be closely related is not surprising upon recalling from Section 3.1 that the global density in the case θ = M ∈ Z + coincides with the global density for the squared singular values of a product M standard complex Gaussian matrices. As well as guiding our overall strategy, [32] will be referred to for the proof of some technical bounds, which we omit below.
It is clear that the Hankel like contour Σ can be deformed to an infinite contour that is from −i · ∞ to i · ∞, as long as it keeps the poles −1, −2, . . . of z to its left and Γ α to its right. Actually this deformation has more freedom. By the residue theorem, if Σ is deformed into the infinite vertical contour to the right of Γ α , the double contour integral in (4.10) remains the same. Moreover, we can split Γ α into two positively oriented contours that jointly enclose all the poles α 1 , . . . , α N , and let Σ be an infinite vertical contour passing between them. See [32, Eq. (2.8)] for the explicit computation in a similar case. Specifically, deform Σ from the Hankel contour into the upward vertical contour To express the shape of the deformed contour Γ α , we define first the contoursΓ andΓ where > 0. Thus definẽ and for a small enough > 0, (see Then we define Γ α , depending on two small positive constants and . Here we take the notational convention that if C is a contour and r > 0, then rC is the contour consisting of {z | z/r ∈ C} and with the same orientation. In terms of this notation and Note that Γ 1 ∪ Γ 3 and Γ 2 ∪ Γ 4 are disjoint closed contours, and we assume that Γ ǫ Figure  1. The schematic shape of Γ for the Laguerre Muttalib-Borodin ensemble (θ > 0). Σ Γ α Figure  2. The schematic shapes of Σ and Γ α for both the Laguerre and Jacobi Muttalib Borodin ensembles (θ > 0). they are both oriented counterclockwise. Here we choose small enough so that N u L ± ± lie between two poles of the integral α j and α j+1 , and then Γ 1 ∪ Γ 3 and Γ 2 ∪ Γ 4 jointly cover all the poles α 1 , . . . , α N . Our contours Σ and Γ α are identical to the contours C and Σ respectively defined in [32, Sec. 2.2] and ourΓ andΓ are identical toΣ andΣ in [32,Sec. 3.1] respectively, if θ = M ∈ Z + . See Figure 2 for the shapes of Σ and Γ α .
Our choice of the contours Σ andΓ α allows u L ± to be identified as maximum points of F (z; x) for z ∈ Σ and minimum points of F (z; x) for z ∈ Γ curved , which is key to the subsequent asymptotic analysis. To be precise, we state the result as follows, where we denote D r (z) = {w ∈ C | |w − z| < r}. Lemma 4.2. There exists δ > 0 such that for all N large enough, The most important ingredient of the proof of Lemma 4.2 is the estimate of the leading termF 1 (z; x), which is identical toF (z; x) in [32] if θ = M . The estimate ofF 1 (z; x) is stated in [32,Lem. 3.1 and 3.2], where the results and the proofs hold for general θ > 0. Then Lemma 4.2 is proved analogously to the proof of [32,Lem. 2.1] in [32,Sec. 3.2]. We omit the detail.
Remark 4.3. If N u L ± happens to be identical to a pole α j , we replace N u L ± into N u L ± + 1/2 in formulas (4.22), (4.26) and (4.27). All later arguments are valid with notational changes.
Taking the limit → 0, we havẽ where, with p.v. denoted the principal value integral, and (4.37) To evaluate I 1 , we define Our strategy is to consider the Cauchy principal integral (4.36) on Σ + local × Γ + local and Σ − local × Γ − local first, and then to show that the remaining part of the integral is negligible in the asymptotic analysis.
Make the change of variables z = N θu L + + N where C + is defined in (4.21). This gives p. v.
Note that the O(N −1/5 ) term in the integrand on the right-hand side of (4.41) is uniform and analytic in N N 3/5 (N θu L + ). Using (4.30) of (4.31) in Lemma 4.2, we find that there exists c 1 , c 2 > 0 such that for z ∈ Σ + local and w ∈ Γ + local , (4.42) Hence standard application of the saddle point method yields Analogous reasoning gives    Remark 4.5. The asymptotic analysis above can also prove the bulk local universality of the limiting distribution of the eigenvalues, as in [32]. The result, which is unsurprisingly the sine universality, is also the same as in [32]. Since the local universality is not our focus in this paper, we omit further discussion on it. In the case θ = 2 this problem, along with Airy kernel universality at the soft edge, was solved some time ago [33]. Very recently, these universalities have been established for general θ > 0 by Zhang [46] according to the method of [32]. In [46], an equivalent form of (4.3) is also derived.
(4.52) Changing variables z = N u, w = N v and using Stirling's formula shows that for large N , as long as arg u, arg v ∈ (−π + , π − ), we see that the stationary points of H 1 (u; x) occur when This equation is solved by the Lambert W function [42,Sec. 4.13] -recall too Remark 3.5(i) -and the two solutions are complex conjugates Then using steepest descent arguments like in the θ > 0 case, we have Proposition 4.6. As θ → 0, the limiting global density of the Muttalib-Borodin ensemble is .
The proof is omitted since it is similar to the θ > 0 case. Note that the relation (4.49) in Remark 4.4 still holds. The result (4.55) was derived in [13, Cor. 1] by a determination of the moments, with the latter also known from [18].
, (4.56) where α j = θ(j − 1) + c 1 as specified in (2.3) and c = c 1 . Analogous to (4.3) we have that the factorh(x)/h(y) in (4.56) cancels out of the determinant (4.1), and the reasoning leading to (4.4) tells us that Here we assume θ > 0. The θ = 0 case of the Jacobi upper-triangular ensemble, that is, α 1 = · · · = α N = c 1 , is well defined and can be thought as the θ → 0 + limit of the Jacobi Muttalib-Borodin ensemble, but we omit it in this paper. Our aim is to use a saddle point analysis to computẽ which is the limiting density of the eigenvalues under the change of variables x → x θ . The required workings is structurally identical to that just given to derive Proposition 4.1; only a brief sketch will be given below.
Recall that in the special case that θ ∈ Z + , the resolventG J (z) defined in (3.22), satisfies the identity (3.23), and then the global densityρ J (x) satisfies Below we show that the global densityρ J (x) has an explicit formula given as follows.
Thus we have, parallel to Proposition 4.1 (4.62) This is analogous to the resolvent for the Fuss-Catalan distribution satisfying the identity z(zG(z) − 1) = (zG(z)) θ+1 hold. Of interest is to extend the characterisation (4.62) to general θ > 0. Such a characterisation was established in [23] for the measures corresponding to the moments (3.44) with r > s, whereas as remarked below (3.45) the moments (3.34) deduced from (4.62) require r = s.
The method of the proof to Proposition 4.7 is the same as the proof of Proposition 4.1 for the Laguerre case. So we only give a sketch of the proof and point out the main differences.
Sketch of the proof of Proposition 4.7. First we consider the case that c 2 ∈ Z. Making the replacements x → x 1/θ , y → y 1/θ in (4.56) gives Changing variables in the integrand according to (4.12), then use of Stirling's formula and (4.14) shows that for large N where, after the change of variables of z, w into u, v as in (4.12), and the validity of the error bound O(N −1 ) is the same as in (4.15), i.e., all z satisfying (4.13) and arg z = arg u ∈ (−π + , π − ). Substituting in (4.63) gives This is structurally identical to (4.18), and is analysed accordingly. For the sake of steepest-descent analysis, we need to deform the contours Σ and Γ α . Since we assume that c 2 ∈ Z, Γ α is a finite contour similar to the Γ α in (4.10). It is clear that Σ can be deformed to an infinite contour that is from −i · ∞ to i · ∞,as long as it keeps the poles −1, . . . , −(β + N ) and Γ α to the left. Then we see that we can split Γ α into two positively oriented contours that jointly enclose all the poles α 1 , . . . , α N , and let Σ be an infinite vertical contour passing between them. Note that the precise shapes of Σ and Γ α in the Laguerre case depend on the computation of the critical points ofF 1 (u; x) that is the counterpart of ourĜ 1 (u; x), see (4.20) and (4.23). Below we explain the analogous construction of Σ and Γ α in the Jacobi case.
A difference in detail is now that so the stationary points now occur for u such that Analogous to the situation with (4.19), we observe that this is precisely the equation (3.23) after the identification z → x, zG J (z) → u in the latter. It is straightforward to verify that if x ∈ (0, 1) is parametrised by ϕ ∈ (0, π/(θ + 1)) in (4.60), then with v(ϕ) defined in (4.61) are two solutions to (4.69). As x runs over (0, 1), the locus of u J + (resp. u J − ) is a curve in C + (resp. C − ) whose two ends are on the real line. Thus we let Σ be the upward vertical contour through N θu J ± , analogous to (4.22) in the Laguerre case. To define Γ α , we first definẽ analogous toΓ in (4.23), and then deform it to the desired Γ α parallel to the deformation carried out in the Laguerre case. We omit the detail, but only point out that the schematic shapes of Σ and Γ α is also shown in Figure 2.
By saddle-point analysis analogous to that in the Laguerre case, we derive, as expected in Remark 4.4, and hence finish the proof.
In the case that c 2 / ∈ Z, the contour Γ α in (4.56) is infinite. We express it as the combination of Γ α ∪ Γ α . Here Γ α is an infinite, Hankel like contour starting at −∞ − i , running parallel to the negative real axis, looping around the poles w = −(c 2 + N + k) with k ∈ Z + , and finishing at −∞ + i after again running parallel to the negative real axis. Γ α is a finite positive oriented contour enclosing α 1 , . . . , α N . Then by argument of contour deformation used in the c 2 ∈ Z case and the Laguerre case, we can deform Σ into a contour going from −i · ∞ to i · ∞ between Γ α and Γ α . We writẽ In the computation ofK J, (x), we further deform the contours Γ α into two parts as the deformation of Γ α in the c 2 ∈ Z case, and let Σ be an vertical contour between the two parts. See Figure 3 for the deformed shapes of Γ α , Γ α and Σ. By the argument in the c 2 ∈ Z case we have that On the other hand, by estimating the factors of the integrand Thus we prove the proposition by combining the two limits above. Figure 3. The deformed shape of contours Σ and Γ α for the computation of the global density with θ > 0 and c 2 / ∈ Z. 5.1.1. The θ > 0 case. The meaning of the hard edge scaling has already been noted below (4.3). It is the N → ∞ limit with x, y scaled in the neighbourhood of the origin so that the spacing between eigenvalues is of order unity.

(5.2)
Proof. First we use the argument in Section 4.1, more specifically the paragraph above (4.22), to justify that the double contour integral formula (4.3) still holds with Σ deformed into a contour from e −(π/2+δ)i · ∞ to e (π/2+δ)i · ∞ and going on the left of Γ α . Thus analogous to (4.22) we can take Σ to be Σ δ −1/2 in (4.3). Furthermore it is clear that the closed contour Γ α can be deformed to the infinite Hankel loop contour Γ 0 . Substituting (2.3) and making use of the functional and reflection formula for the gamma function shows that Taking the limit N → ∞ inside the integral (a step which is justified by estimating the integrand and using dominated convergence; see [30,Sec. 5.2] for the details in a very similar setting) we thus see that x −z−1 y w z − w Γ(z + 1)Γ((z − c)/θ + 1) sin π(z − c)/θ Γ(w + 1)Γ((w − c)/θ + 1) sin π(w − c)/θ .
Borodin has previously given a different formula for the hard edge scaled limit in (5.1).
(ii) Here we use a nearly vertical contour Σ δ −1/2 to be in consistent with [30]. It is also possible to deform Σ into a Hankel-like form, which is symmetric to Γ 0 . With the help of this symmetry, and under the change of variables z → −(z + (c + 1)/θ), w → −(w + (c + 1)/θ) in (5.9) the contours Σ and Γ 0 interchange. Making use of the reflection formula for the gamma function then allows us to deduce that 1 θ x 1/θ−1 K (α,θ) (x 1/θ , y 1/θ ) = x y α K (α ,1/θ) (y, x), (5. 10) which as noted in [29] also follows from the original form (5.4). (iii) When our article was almost complete, we received a preprint from Zhang [46], containing amongst other things an independent analysis of the hard edge scaling of the Laguerre Muttalib-Borodin ensemble.

5.1.2.
The θ = 0 case. We could take the limit of the double contour integral formula (4.3) with α 1 = · · · = α N = c to derive the limiting correlation kernel near 0 for the θ = 0 case of the Laguerre Muttalib-Borodin ensemble, but here we introduce an alternative approach. Note that the joint PDF of the θ = 0 Laguerre Muttalib-Borodin ensemble is given in (2.21). Changing variables log λ j → µ j , we have that the joint PDF for µ 1 ≥ µ 2 ≥ · · · ≥ µ N is proportional to Then the results in [15] yields that the correlation functions for the smallest variables µ N , µ N −1 , . . . converge to the correlation functions in the Tracy-Widom distribution, upon proper scaling. (Note that in [15] there is a technical assumption that V (x) → +∞ faster than any linear equation as x → ±∞. But it is not hard to see that in the −∞ direction this requirement can be relax to a linear growth to +∞.) 5.2. The Jacobi Muttalib-Borodin ensemble. The hard edge scaled limit of the Jacobi Muttalib-Borodin ensemble is the same as that for the Laguerre Muttalib-Borodin ensemble, although the specific scale is different [9]. (5.11) Proof. First we consider the case that c 2 ∈ Z. Then the contour Γ α in the integral (4.56) is a finite contour enclosing α 1 , . . . , α N . By the argument similar to that in Section 4.2, we can deform the contour Σ in (4.56) into a vertical contour going from e −(π/2+δ)i · ∞ to e (π/2+δ)i · ∞ on the left of Γ α . Furthermore we can deform the closed contour Γ α into the infinite contour Γ 0 . Thus the double integral formula (4.56) still holds with Σ and Γ α replaced by Σ δ −1/2 and Γ 0 respectively. Comparing the integrands of K L (4.3) and K J (4.56) we see that the only difference is that in the latter there is an additional factor Γ(w + c 2 + N + 1)/Γ(z + c 2 + N + 1). For z and w fixed and N → ∞, this simplifies to Γ(w + c 2 + N + 1) Γ(z + c 2 + N + 1) = N (w−z) 1 + O(1/N ) , so the only difference in the present case is an additional factor of N (w−z) . The additional factor is accounted for by the different scalings: x → N −1/θ x for the Laguerre case, and x → N −1−1/θ x for the Jacobi case.