Matrix product ensembles of Hermite-type and the hyperbolic Harish-Chandra-Itzykson-Zuber integral

TL;DR

This paper studies the spectral properties of a Hermitised random matrix product with eigenvalues on the full real line, deriving explicit formulas for bi-orthogonal functions and correlation kernels, and introducing a new matrix transformation related to hyperbolic integrals.

Contribution

It introduces a new matrix transformation that maps polynomial ensembles onto themselves and connects it to the hyperbolic Harish-Chandra-Itzykson-Zuber integral, advancing understanding of matrix product spectral properties.

Findings

Eigenvalues form a bi-orthogonal ensemble

Asymptotic reduction to Hermite Muttalib-Borodin ensemble

Explicit correlation kernel involving Meijer G-functions

Abstract

We investigate spectral properties of a Hermitised random matrix product which, contrary to previous product ensembles, allows for eigenvalues on the full real line. We prove that the eigenvalues form a bi-orthogonal ensemble, which reduces asymptotically to the Hermite Muttalib-Borodin ensemble. Explicit expressions for the bi-orthogonal functions as well as the correlation kernel are provided. Scaling the latter near the origin gives a limiting kernel involving Meijer G-functions, and the functional form of the global density is calculated. As a part of this study, we introduce a new matrix transformation which maps the space of polynomial ensembles onto itself. This matrix transformation is closely related to the so-called hyperbolic Harish-Chandra-Itzykson-Zuber integral.

Tables

Table 1. Table 1: Summary of weight functions and support for the three canonical orthogonal ensembles in random matrix theory given by the joint distribution ( 1.6 ).

ensemble	weight	support
Jacobi	$w\left(\lambda \right)={\lambda }^a{\left(1-\lambda \right)}^b$	$\lambda \in (0,1)$
Laguerre	$w\left(\lambda \right)={\lambda }^a{e}^{-\lambda }$	$\lambda \in (0,\infty )$
Hermite	$w\left(\lambda \right)={\left|\lambda \right|}^a{e}^{-{\lambda }^2}$	$\lambda \in (-\infty ,\infty )$

Equations

$$W_M=G_M^{†}\cdots G_1^{†}H{G}_1\cdots G_M$$

$$\tilde{P}(x_1,\dots ,x_N)=\frac{1}{{\tilde{Z}}_N^{(M)}}\prod _{1\le j<k\le N}(x_k-x_j)(x_k^{2M+1}-x_j^{2M+1})\prod _{\ell =1}^N|x_{\ell }{|}^{\alpha }{e}^{-x_{\ell }^2},$$

$$a_1<a_2<\cdots <a_{n_0}<0<a_{n_0+1}<\cdots <a_n,$$

$$P^{n_{0}}_{n}(\{a_{j}\}_{j=1}^{n};\{x_{j}\}_{j=1}^{n})=\\ \prod_{l=1}^{n}\frac{1}{|a_{l}|}\frac{(x_{l}/a_{l})^{N-n}}{(N-l)!}\prod_{1\leq j<k\leq n}\frac{x_{k}-x_{j}}{a_{k}-a_{j}}\det\big{[}e^{-x_{i}/a_{j}}\big{]}_{i,j=1}^{n_{0}}\det\big{[}e^{-x_{i}/a_{j}}\big{]}_{i,j=n_{0}+1}^{n},$$

$$x_1<\cdots <x_{n_0}<0<x_{n_0+1}<\cdots <x_n.$$

$$P(d{x}_1,\dots ,d{x}_n)=\frac{1}{Z_n}{\Delta }_n(\{x\}{)}^2\prod _{k=1}^{n}w(x_k)d{x}_k,$$

$$\Delta_{n}(\{x\})=\det_{1\leq i,j\leq n}\big{[}x_{i}^{\,j-1}\big{]}=\prod_{1\leq i<j\leq n}(x_{j}-x_{i}).$$

$$P(d{x}_1,\dots ,d{x}_n)=\frac{1}{Z_n}{e}^{-\beta E(x_1,\dots ,x_n)}\prod _{k=1}^{n}d{x}_k$$

$$E(x_1,\dots ,x_n)=\frac{1}{2}\sum _{k=1}^{n}V(x_k)-\sum _{1\le i<j\le n}U(x_j,x_i).$$

$$G_M^{†}\cdots G_1^{†}{G}_1\cdots G_M.$$

$$P^{(M)}(x_{1},\ldots,x_{n})=\frac{1}{Z_{n}^{(M)}}\Delta_{n}(\{x\})\det_{1\leq i,j\leq n}\big{[}g_{j}^{(M)}(x_{i})\big{]},\qquad x_{i}>0\>(i=1,\dots,n)$$

$${\tilde{P}}^{(M)}(x_1,\dots ,x_n)=\frac{1}{{\tilde{Z}}_n^{(M)}}{\Delta }_n(\{x\}){\Delta }_n(\{x^M\})\prod _{k=1}^n{x}_k^a{e}^{-x_k},x_k>0(k=1,\dots ,n)$$

$$P(x_{1},\ldots,x_{n})=\frac{1}{Z_{n}}\prod_{j=1}^{n}w(x_{l})\prod_{1\leq j<k\leq n}\big{|}x_{k}-x_{j}\big{|}\,\big{|}\operatorname{sgn}(x_{k})|x_{k}|^{\theta}-\operatorname{sgn}(x_{j})|x_{j}|^{\theta}\big{|},$$

$${\int }_{U(N)/U(1{)}^N}{e}^{-\mathrm{Tr}AVB{V}^{-1}}(V^{-1}dV)={\pi }^{N(N-1)/2}\frac{\det [e^{-a_i{b}_j}{]}_{i,j=1}^N}{{\prod }_{1\le i<j\le N}(a_j-a_i)(b_j-b_i)},$$

$$P_A(\{a_k{\}}_{k=1}^n)=\frac{1}{Z_n}\prod _{1\le i<j\le n}(a_j-a_i)\det [w_j(a_i){]}_{i,j=1}^n,$$

$$\frac{1}{Z_{n}}\prod_{l=1}^{n}\frac{1}{(N-l)!}\prod_{1\leq j<k\leq n}(x_{k}-x_{j})\det_{1\leq i,j\leq n}\bigg{[}\int_{0}^{\infty}\frac{da\,e^{-a}}{a^{n-N+1}}\,w_{j}\big{(}\frac{x_{i}}{a}\big{)}\bigg{]}$$

$$\begin{gathered} {\int }_D{P}_A(\{a_k{\}}_{k=1}^n)P_n^{n_0}(\{a_j{\}}_{j=1}^n;\{x_j{\}}_{j=1}^n)d{a}_1\cdots d{a}_n=\frac{1}{Z_n}\prod _{1\le j<k\le n}(x_k-x_j) \\ \times {\int }_D\prod _{l=1}^n\frac{1}{|a_l|}\frac{(x_l/a_l{)}^{N-n}}{(N-l)!}\det \left[\begin{array}{cc}\{e^{-x_i/a_j}{\}}_{i,j=1}^{n_0}& 0\\ 0& \{e^{-x_i/a_j}{\}}_{i,j=n_0+1}^n\end{array}\right]\det [w_j(a_i){]}_{i,j=1}^{n}d{a}_1\cdots d{a}_n, \end{gathered}$$

$$\frac{1}{Z_{n}}\prod_{l=1}^{n}\frac{1}{(N-l)!}\prod_{1\leq j<k\leq n}(x_{k}-x_{j})\det\begin{bmatrix}\displaystyle\bigg{\{}-\int_{-\infty}^{0}\frac{da}{a}(x_{i}/a)^{N-n}e^{-x_{i}/a}w_{j}(a)\bigg{\}}_{\begin{subarray}{c}i=1,\ldots,n_{0}\\ j=1,\ldots,n\end{subarray}}\\ \displaystyle\bigg{\{}\int^{\infty}_{0}\frac{da}{a}(x_{i}/a)^{N-n}e^{-x_{i}/a}w_{j}(a)\bigg{\}}_{\begin{subarray}{c}i=n_{0}+1,\ldots,n\\ j=1,\ldots,n\end{subarray}}\end{bmatrix}.$$

$$\frac{1}{Z_{n}}\prod_{l=1}^{n}\frac{1}{(N-l)!}\prod_{1\leq j<k\leq n}(x_{k}-x_{j})\det_{1\leq i,j\leq n}\bigg{[}\int_{0}^{\infty}\frac{da}{a}\,\frac{e^{-|x_{i}|/a}}{(|{x_{i}}|/a)^{n-N}}\,w_{j}((\operatorname{sgn}x_{i})a)\bigg{]};$$

$$P^{(M)}(x_{1},\ldots,x_{n})=\frac{1}{Z^{(M)}_{n}}\prod_{1\leq i<j\leq n}(x_{j}-x_{i})\det_{1\leq i,j\leq n}\big{[}g_{j-1}^{(M)}(x_{i})\big{]},$$

$$g_j^{(0)}(x)=x^j{e}^{-x^2}\text{and}{g}_j^{(m)}(x)={\int }_0^{\infty }\frac{dy}{y}{y}^{{\nu }_m}{e}^{-y}{g}_j^{(m-1)}(x/y),m=1,\dots ,M$$

$$Z_n^{(M)}=2^{-n(n-1)/2}{\pi }^{n/2}\prod _{m=0}^M\prod _{j=1}^n\Gamma ({\nu }_m+j).$$

$$W_2=G_2^{†}{W}_1{G}_2\stackrel{d}{=}\left[\begin{array}{cc}{\tilde{G}}_2^{†}& g_2^{†}\end{array}\right]\left[\begin{array}{cc}{X}_1& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{\tilde{G}}_2\\ g_2\end{array}\right]={\tilde{G}}_2^{†}{X}_1{\tilde{G}}_2,$$

$$g_{j}^{(M)}(x)=\frac{(\operatorname{sgn}x)^{j}}{4\pi i}\int_{c-i\infty}^{c+i\infty}ds\,|x|^{s}\,\Gamma\Big{(}\frac{j-s}{2}\Big{)}\prod_{m=1}^{M}\Gamma(\nu_{m}-s),$$

$$g_{j}^{(0)}(x)=x^{j}e^{-x^{2}}=\frac{(\operatorname{sgn}x)^{j}}{4\pi i}\int_{c-i\infty}^{c+i\infty}ds\,|x|^{s}\,\Gamma\Big{(}\frac{j-s}{2}\Big{)}.$$

$$g_{j}^{(M)}(x)=\int_{0}^{\infty}\frac{dy}{y}\,y^{\nu_{M}}e^{-y}\,\frac{(\operatorname{sgn}x)^{j}}{4\pi i}\int_{c-i\infty}^{c+i\infty}ds\,\Big{(}\frac{|x|}{y}\Big{)}^{s}\,\Gamma\Big{(}\frac{j-s}{2}\Big{)}\prod_{m=1}^{M-1}\Gamma(\nu_{m}-s).$$

$$g_{j}^{(M)}(x)=\frac{(\operatorname{sgn}x)^{j}}{4\pi i}\int_{c-i\infty}^{c+i\infty}ds\,|x|^{s}\,\Gamma\Big{(}\frac{j-s}{2}\Big{)}\prod_{m=1}^{M-1}\Gamma(\nu_{m}-s)\,\int_{0}^{\infty}\frac{dy}{y}\,y^{\nu_{M}-s}e^{-y}$$

$$g_{j}^{(M)}(x)=\frac{(\operatorname{sgn}x)^{j}}{2}H^{{M+1},{0}}_{{0},{M+1}}\bigg{(}\begin{matrix}-\\ (\nu_{1},1),\ldots,(\nu_{M},1),(\frac{j}{2},\frac{1}{2})\end{matrix}\,\bigg{|}\,|x|\bigg{)}$$

$$g_{j}^{(M)}(x)=(\operatorname{sgn}x)^{j}\prod_{m=1}^{M}\frac{2^{\nu_{m}-1}}{\sqrt{\pi}}G^{{2M+1},{0}}_{{0},{2M+1}}\bigg{(}\begin{matrix}-\\ \frac{\nu_{1}}{2},\frac{\nu_{1}+1}{2},\ldots,\frac{\nu_{M}}{2},\frac{\nu_{M}+1}{2},\frac{j}{2}\end{matrix}\,\bigg{|}\,\frac{x^{2}}{4^{M}}\bigg{)}.$$

$$G^{{q},{0}}_{{0},{q}}\bigg{(}\begin{matrix}-\\ b_{1},\ldots,b_{q}\end{matrix}\,\bigg{|}\,x\bigg{)}\sim\frac{1}{q^{1/2}}\Big{(}\frac{2\pi}{x^{1/q}}\Big{)}^{(q-1)/2}x^{(b_{1}+\cdots+b_{q})/{q}}e^{-qx^{1/q}}\big{(}1+O(x^{1/q})\big{)},$$

Full text

1 Introduction

1.1 Statement of the problem and summary of results

Let $H$ be a matrix from the Gaussian unitary ensemble (GUE) and let each $G_{i}$ $(i=1,\dots,M)$ denote a complex Ginibre matrix, i.e. a matrix with i.i.d. standard complex Gaussian entries. In this paper, we investigate the eigenvalues of the Hermitised product matrix

[TABLE]

under the assumption that all matrices, $H$ and $G_{i}$ ($i=1,\ldots,M$), are independent. We will see that the eigenvalues form a bi-orthogonal ensemble [11]. Furthermore, this ensemble is closely related (in a sense that will be specified in the next subsection) to the so-called Hermite Muttalib–Borodin ensemble [47, 11]. The latter is defined by the joint eigenvalue probability density function (PDF)

[TABLE]

where $\tilde{Z}_{N}^{(M)}$ is a normalisation constant and $\alpha$ is a non-negative constant.

It will transpire that the bi-orthogonal ensemble structure associated with the eigenvalue PDF of the product matrix (1.1) is a corollary of the following more basic result.

Theorem 1**.**

Let $G$ be an $n\times N$ ($n\leq N$) standard complex Gaussian matrix and let $A$ be an $n\times n$ Hermitian matrix with eigenvalues $a_{1},\ldots,a_{N}$. If the eigenvalues of $A$ are pairwise distinct and ordered as

[TABLE]

then the PDF of the non-zero eigenvalues of matrix $X=G^{\dagger}AG$ is given by

[TABLE]

where

[TABLE]

In particular, we see that $X$ has $n_{0}$ ($n-n_{0}$) negative (positive) eigenvalues, i.e. the same number as $A$. The remaining $N-n$ eigenvalues are all identically zero.

We remark that in Theorem 1 the case of $n<N$ remains unanswered although this is certainly of high interest; see [18, 2] and references therein for a relevant question.

The rest of this paper is organised as follows: In Section 2, we use Theorem 1 to find the PDF for the eigenvalues of the product (1.1) as a bi-orthogonal ensemble. Moreover, the explicit expression for the PDF is seen to reduce asymptotically to the functional form (1.2) specifying the Hermite Muttalib–Borodin ensemble. Explicit expressions for the bi-orthogonal functions are derived in Section 3. Analogous to the theory of Hermite polynomials (see e.g. (4.15) below), we will see that it is convenient to consider bi-orthogonal functions of even and odd degree separately. Section 4 provides reformulations of the bi-orthogonal functions and the correlation kernel as integral representations, which are more suited for asymptotic analysis. These integral representations can also be expressed in terms of Meijer $G$-functions, and we will see that they are closely related to known formulae stemming from the product ensemble of Laguerre type. The local scaling limit at the origin is derived and seen to be related to the Meijer $G$-kernel in Section 5. This result is also compared with the local scaling limit of the Hermite Muttalib–Borodin ensemble. Section 6 includes derivations of the global spectrum of the product (1.1) as well as the Hermite Muttalib–Borodin ensembles. Since our product ensemble reduces asymptotically to the Hermite Muttalib–Borodin ensemble, they are as expected seen to have the same global spectrum, which in turn is given in terms of the Fuss–Catalan density. Finally, Theorem 1 is proven in the appendix. This theorem is an important result by itself. For this reason, we provide three separate proofs each with their own merits.

1.2 First motivation: Muttalib–Borodin ensembles

Orthogonal polynomial ensembles are point processes on (a subset of) the real line with a joint distribution given by

[TABLE]

where $Z_{N}$ is a normalisation constant, $w(x)$ is a non-negative weight function, and $\Delta_{n}(\{x\})$ denotes the Vandermonde determinant,

[TABLE]

Like the corresponding moment problem, it is often useful to distinguish between models with support on a finite, semi-infinite, and double-infinite interval. The canonical examples are the Jacobi-, Laguerre-, and Hermite-ensembles summarised in Table 1. These ensembles are named according to the corresponding classical orthogonal polynomials. In fact, for $a\neq 0$ the latter would more appropriately be called the generalised Hermite ensemble.

In random matrix theory these three ensembles play a fundamental role as they appear as the distribution of the eigenvalues (or singular values) for the transfer (or truncated unitary) ensemble, the complex Wishart (or chiral) ensemble, and the Gaussian unitary ensemble, respectively; see e.g. [20].

A fundamental insight, which can be traced back to Wigner [55], is that the joint distribution (1.6) allows an interpretation as the equilibrium measure for a one-dimensional gas of pairwise repulsive point particles in a confining potential. More precisely, consider the Gibbs measure for a classical gas of $n$ point particles which are pairwise repulsive according to a two-point potential $U(x,y)$ and confined by a common one-point potential $V(x)$, i.e.

[TABLE]

with $\beta$ denoting the inverse temperature and $E$ the energy functional

[TABLE]

We see that at $\beta=2$ (sometimes referred to as the free fermion point) the Gibbs measure (1.8) is identical to (1.6) provided we set $V(\lambda)=-\log w(\lambda)$ and $U(x_{i},x_{j})=\log|x_{j}-x_{i}|$. In this way, the eigenvalues of random matrices relate to the Boltzmann factor of a simple statistical mechanical system with one- and two-body interactions only.

A recent development in random matrix theory is the study of exactly solvable product ensembles; see [3] for a review. As an example, let $G_{i}$ ($i=1,\ldots,M$) be independent complex Ginibre matrices (matrices whose entries are i.i.d. standard complex Gaussians) and consider the Hermitian product

[TABLE]

From the work of Akemann et al. [5, 4], we know that the explicit PDF for the eigenvalues of the matrix (1.10) is

[TABLE]

where $Z_{n}^{(M)}$ is a known normalisation constant and $g_{j}^{(M)}(x)$ ($j=1,\ldots,n$) are given by certain Meijer $G$-functions. Generally, such PDFs are known as polynomial ensembles [40].

It seems natural to ask whether these product ensembles also have (at least approximately) an interpretation as a Gibbs measure of the form (1.8) and (1.9). However, unlike the Vandermonde determinant, the determinant in (1.11) cannot be evaluated as a product (for $M\geq 2$). This prohibits a literal interpretation of the eigenvalues of (1.10) as a statistical mechanical system with only one- and two-body interactions. One could fear that this meant that there was no simple physical interpretation related to (1.11). However, if we consider (1.11) with each $x_{j}$ large, the Meijer $G$-functions can be replaced by their asymptotic approximation [19]. After a change of variables, the joint density (1.11) to leading order in the asymptotic expansion becomes [25]

[TABLE]

where $a$ is a known non-negative constant. This does correspond to the Boltzmann factor of a statistical mechanical system with one- and two-body interactions only.

A comparison between (1.11) and (1.12) can be done a posteriori. A connection between the two ensembles was first noted by Kuijlaars and Stivigny [40], who observed that the hard edge scaling limit of (1.12) found in [11] took the same functional form as the Meijer $G$-kernel found in the product ensemble [41], albeit with a different choice of parameters. Due to recent progress, even more is known about the scaling limits of both models, and their similarities. Thus it has been established that the two ensembles also share the same global spectral distribution [46, 13, 10, 51, 27]. Furthermore, in both cases the local correlations in the bulk and near the soft edge are given by the familiar sine and Airy process, respectively [44, 56].

The ensemble (1.12) had, in fact, appeared in earlier random matrix literature. It was first isolated by Muttalib [47], who suggested it as a naive approximation to the transmission eigenvalues in a problem about quantum transport. A feature of the new interaction is that bi-orthogonal polynomials (rather than orthogonal polynomials) are needed in the study of correlation functions. Such bi-orthogonal ensembles were considered in greater generality by Borodin [11], who devoted special attention to PDFs

[TABLE]

with $\theta>0$ and $w(x)$ representing one of the three classical weight functions from Table 1. Following [27], we will henceforth refer to these ensembles as the (Jacobi, Laguerre, Hermite) Muttalib–Borodin ensembles. We note that the awkward dependence of signs in the last factor in (1.13) disappears when the eigenvalues are non-negative (e.g. for Laguerre- and Jacobi-ensembles) and when $\theta$ is an odd integer as in (1.2).

At the time of their introduction, the Muttalib–Borodin ensembles had no obvious relation to any random matrix models defined in terms of PDFs on their entries (except for the trivial case $\theta=1$), and could merely be interpreted as a simple one-parameter generalisation of the classical ensembles. However, we now see that the Laguerre Muttalib–Borodin ensemble has a close connection to products of complex Gaussian random matrices (1.10) through the approximation (1.12).

Knowing that the Laguerre Muttalib–Borodin ensemble appears as an asymptotic approximation to the Gaussian product (1.10), it seems natural to ask the reverse question: Can we find product ensembles which reduce asymptotically to the Jacobi and Hermite Muttalib–Borodin ensembles? If this is possible, it would be reasonable to say we have completed a link between the Muttalib–Borodin ensembles with classical weights and the new family of product ensembles.

For the Jacobi Muttalib–Borodin ensemble a link to products of random matrices is provided by looking at the squared singular values of a product of truncated unitary matrices [35, 27]. In this paper, it is our aim to isolate a random matrix product structure for which the eigenvalue PDF reduces asymptotically to the functional form of the Hermite Muttalib–Borodin ensemble. This construction therefore completes the correspondence between product ensembles and the three Muttalib–Borodin ensembles with classical weights, i.e. Laguerre, Jacobi, Hermite. Furthermore, the relevant product ensemble provides by itself a new interesting class of integrable models, which unlike all previous product ensembles (see the review [3]) allows for negative eigenvalues.

As the product ensemble in question must allow for negative eigenvalues, it is no longer sufficient to investigate Wishart-type matrices like (1.10) which are positive-definite by construction. It turns out that the correct structure is the Hermitised product of a GUE matrix and $M$ complex Ginibre matrices given by (1.1). The case $M=1$ of (1.1) has previously been isolated in the recent paper of Kumar [42] as an example of a matrix ensemble which permits an explicit eigenvalue PDF.

1.3 Second motivation: hyperbolic Harish-Chandra–Itzykson–Zuber integrals

Another reason that the Hermitised random matrix product (1.1) is of particular interest is its relation to the so-called hyperbolic Harish-Chandra–Itzykson–Zuber (HCIZ) integral. By way of introduction on this point, we note that it is by now evident that the family of exactly solvable product ensembles is intimately linked to a family exactly solvable group integrals sometimes referred to as integrals of HCIZ type. For the study of products of Ginibre matrices (1.10) it was sufficient to know the familiar (and celebrated) HCIZ integral [30, 33]:

[TABLE]

where $(V^{-1}dV)$ denotes the Haar measure on the unitary quotient group $U(N)/U(1)^{N}$, while $A$ and $B$ are Hermitian $N\times N$ matrices with eigenvalues $a_{1}<\cdots<a_{N}$ and $b_{1}<\cdots<b_{N}$, respectively. However, for studies of products of spherical, truncated unitary, or coupled random matrices generalisations of the HCIZ integral are needed, see [35, 6, 43] for the two latter cases. We emphasise that the product of truncated unitary matrices considered by Kieburg et al. [35] required a previously unknown generalisation of the HCIZ integral. Likewise, our study of the Hermitised random matrix product (1.1) requires knowledge about the so-called hyperbolic HCIZ integral in which the integration on the left-hand side of (1.14) should be replaced with an integration over the pseudo-unitary group (see Section A.3 for details). The study of such hyperbolic group integrals was initiated by Fyodorov [28, 29]. An interesting feature of the hyperbolic HCIZ integral is that the integration over the pseudo-unitary is non-compact, which forces us to introduce some additional constraints on the Hermitian matrices $A$ and $B$ to ensure convergence; this is a difficulty which does not arise in other HCIZ-type integrals. Finally, we mention that HCIZ-type integrals have other applications in theoretical and mathematical physics beyond products of random matrices, e.g. the hyperbolic HCIZ integral was used to find the spectral properties of the Wilson–Dirac operator in lattice quantum chromodynamics [36]. Moreover, HCIZ-type integrals represent a rich area of mathematical research, for example within the study of Lie groups, harmonic analysis, combinatorics, and probability (e.g. matrix-valued Brownian motion); see e.g. the text [54].

2 Products of random matrices and Hermite Muttalib–Borodin ensembles

In this section, we establish that the eigenvalue PDF of the matrix product (1.1) is a polynomial ensemble and show that it reduces asymptotically to the Hermite Muttalib–Borodin ensemble (1.2).

As stated in the introduction, the eigenvalue PDF of (1.1) follows as a consequence of Theorem 1. The idea is simple: let $A$ be a random matrix from a polynomial ensemble, i.e. it has an eigenvalue PDF of the form

[TABLE]

where $a_{1}\leq a_{2}\leq\cdots\leq a_{n}$ are the (ordered) eigenvalues of $A$, $w_{j}$ ($j=1,\ldots,n$) is a family of weight functions, and $Z_{n}$ is a normalisation constant. Now, let $G$ be an $n\times N$ $(n\leq N)$ standard complex Gaussian matrix. Then Theorem 1 gives the eigenvalue PDF of $G^{\dagger}AG$. Moreover, it is seen that this new eigenvalue PDF is also a polynomial ensemble. In other words, Theorem 1 provides a map from the class of polynomial ensembles into itself. Thus, we may apply Theorem 1 recursively to construct hierarchies of polynomial ensembles. Let us make this statement more precise.

Lemma 2**.**

Let $G$ be an $n\times N$ $(n\leq N)$ standard complex Gaussian matrix, and let $A$ be a random matrix from a polynomial ensemble with eigenvalue PDF (2.1), independent of $G$. Then the PDF for the non-zero eigenvalues of the random matrix product $G^{\dagger}AG$ is equal to

[TABLE]

with the eigenvalues ordered $x_{1}\leq x_{2}\leq\cdots\leq x_{n}$.

Proof.

In order to use Theorem 1, we fix an $n_{0}\in\{0,1,\ldots,n\}$ and assume that the eigenvalues of $A$ are ordered as (1.3). Consequently, the non-zero eigenvalues of $G^{\dagger}AG$ can be ordered as (1.5) almost surely. It follows from the conditional eigenvalue PDF (1.4) and (2.1) that the eigenvalue PDF of $G^{\dagger}AG$ (up to $N-n$ eigenvalues which are identically zero) is given by

[TABLE]

where the domain of integration $D$ is given according to (1.3) and the eigenvalues $x_{i}$ ($i=1,\ldots,n$) are ordered according to (1.5). We note that the integral on the second line in (2.3) is a close cousin to the well-known Andreief integral [8, 12]. Upon expansion of the determinants, it is readily seen that (2.3) may rewritten as

[TABLE]

If we make a change of variables $a\mapsto-a$ in the first $n_{0}$ rows in the determinant in (2.4) then we get

[TABLE]

recall that $x_{i}$ ($i=1,\ldots,n$) is ordered according (1.5). Finally, if we make another change of variables $a\mapsto|x_{i}|/a$ in the $i$-th row, then (2.2) follows for any fixed $n_{0}$. Note that this result is independent of the choice of $n_{0}$, so we have the final result. ∎

Remark 3**.**

The study of maps from the space of polynomial ensembles onto itself is an interesting endeavour, since such maps give rise to new random matrix ensembles without destroying integrability. In fact, the study of such maps is presently an active area of research in random matrix theory [15, 38, 34, 39]. Lemma 2 provides a new transformation to this class of maps, which cannot be obtained directly from any of the previously established transformations. We note that Lemma 2 includes the transformation [40, Theorem 2.1] as a special case arising when the matrix $A$ is positive definite. A restriction of Lemma 2 is that $n\leq N$. Thus it is seen that the PDF (2.2) develops a singularity for $n>N$ indicating that the formula is no longer generally valid in this case, depending on the properties of $w_{j}$. It would be interesting to extend the above results to include the case $n>N$ more generally.

With Lemma 2 at hand, we are ready to write down the eigenvalue PDF for the product (1.1).

Theorem 4**.**

Let $\nu_{0}=0,\nu_{1},\ldots,\nu_{M}$ be non-negative integers. Suppose that $H$ is an $n\times n$ GUE matrix and $G_{1},\ldots,G_{M}$ are independent standard complex Gaussian matrices where $G_{m}$ is of size $(\nu_{m-1}+n)\times(\nu_{m}+n)$. Then the joint PDF for the non-zero eigenvalues of the matrix (1.1) is given by

[TABLE]

where the weight functions $g_{j}^{(M)}$ are defined recursively by

[TABLE]

and the normalisation constant is

[TABLE]

Proof.

First, let us consider the simplest situation, that is a product of square matrices, i.e. $\nu_{1}=\cdots=\nu_{M}=0$. As the eigenvalue PDF of an $n\times n$ GUE matrix is given by (2.6) with $M=0$, the theorem follows immediately by applying Lemma 2 $M$ times.

We need to be a little more careful when the case of rectangular matrices is considered. The $M=1$ case of the theorem is again an immediate consequence of Lemma 2, which gives us the eigenvalues of $W_{1}=G_{1}^{\dagger}HG_{1}$. However, in order to apply Lemma 2 a second time and find the non-zero eigenvalues of $W_{2}=G_{2}^{\dagger}W_{1}G_{2}$, we have to take into account that $W_{1}$ has a zero eigenvalue with multiplicity $\nu_{1}$. To proceed, we can use the same idea as in [32]. The unitary invariance of Gaussian matrices tells us that $W_{1}\stackrel{{\scriptstyle d}}{{=}}U^{\dagger}W_{1}U$ and $G_{2}\stackrel{{\scriptstyle d}}{{=}}VG_{2}$ for any $U,V\in U(n+\nu_{1})$. It thus follows

[TABLE]

where $X_{1}=\operatorname*{diag}(x_{1},\ldots,x_{n})$ is an $n\times n$ diagonal matrix distributed according to (2.6) with $M=1$, while $\tilde{G}_{2}$ and $g_{2}$ are standard Gaussian matrices of size $n\times(n+\nu_{2})$ and $\nu_{1}\times(n+\nu_{2})$, respectively. Now, Lemma 2 can be applied to the right-hand side in (2.9), which gives us the PDF of the non-zero eigenvalues of $W_{2}$. Repeating this procedure completes the proof. ∎

Remark 5**.**

We note that the case $M=1$ of Theorem 4 is in agreement with the result stated by Kumar [42, Eq. (46) and (47)]. However, the derivation therein is incomplete due to the reliance on the HCIZ integral (1.14), rather than its hyperbolic variant (A.38) below.

There are many other representations for the weight functions in Theorem 4 beyond the recursive definition (2.7). As usual, it is particularly useful for analytic purposes to write the weight functions in their contour integral representation.

Lemma 6**.**

We have

[TABLE]

where $c<0$ is a negative constant.

Proof.

By means of the residue theorem, it is seen that

[TABLE]

Now, assume that $g_{j}^{(M-1)}(x)$ is given by (2.10). From the recursive formula (2.7), we have

[TABLE]

It is a straightforward exercise, considering the asymptotic decay, to show that with $c<0$ the order of the integrals may be interchanged. Thus we have

[TABLE]

and the lemma follows by induction. ∎

As already mentioned, functions with a contour integral representation like (2.10) have certain properties which are useful for analytical purposes. In fact, many of these properties may be found in the literature if we first recognise the contour integral as a Fox $H$-function

[TABLE]

or as a Meijer $G$-function

[TABLE]

We refer to the book [45] for an extensive review of these functions; the Fox $H$- and Meijer $G$-functions are defined by [45, Def. 1.1] and [45, Def. 1.5], respectively.

As discussed in Section 1.2, one of our goals is to find a ‘classical gas’ approximation for (2.6). For this purpose, we can use the asymptotic result [19]

[TABLE]

for $x\to\infty$ (recall that a typical eigenvalue grows with $n$). This immediately allows us to find a Muttalib–Borodin ensemble (1.13), which approximates the product ensemble (2.6). However, for notational simplicity, it is convenient to first make a change of variables

[TABLE]

for $i=1,\ldots,n$. Using the asymptotic formula (2.16) and making a change of variable (2.17), we find the approximate PDF

[TABLE]

where $\tilde{Z}^{(M)}$ is a new normalisation constant and $\alpha=\sum_{m=1}^{M}(2\nu_{m}+1)$. We recognise (2.18) as the Hermite Muttalib–Borodin ensemble (1.2). Note that the approximation (2.16) is valid for large $x$ and that the absolute value of a typical eigenvalue grows with the matrix dimension $n$. Thus, one might suspect agreement between the two models in the large-$n$ limit except for local correlations near the origin. We will return to a comparison between the two models in Section 5 and Section 6.

It is worth noting that the following exact relation holds

[TABLE]

for integer $q$. This may be proven by writing the Meijer $G$-function on the left-hand side as its integral representation and then using Gauss’ multiplication formula for the gamma functions. The exact relation (2.19) tells us that we can choose the parameters $b_{k}$ $(k=1,\ldots,q)$ in (2.16) such that all subleading terms in the expansion vanish, and (2.18) is exact.

3 Biorthogonality and correlations

Generally polynomial ensembles describe determinantal point processes. The correlation kernel may be written as

[TABLE]

where $p_{k}(x)$ and $\phi_{k}(x)$ are bi-orthogonal functions with normalisation $h_{k}$, i.e.

[TABLE]

For both (2.6) and (2.18) the $p_{k}(x)$ will be a monic polynomial of degree $k$, while

[TABLE]

for the product ensemble (2.6) and

[TABLE]

for the Hermite Muttalib–Borodin ensemble (2.18). In the latter case, the bi-orthogonal structure is already known [37, 14, 11, 22]. The main purpose of this section is to determine the bi-orthogonal functions — and consequently the kernel (3.1) — for the product ensemble (2.6).

3.1 The oddness of being even

For $M=0$, the bi-orthogonal functions are $p_{n}(x)=\tilde{H}_{n}(x)$ and $\phi_{n}(x)=\tilde{H}_{n}(x)e^{-x^{2}}$ with $\tilde{H}_{n}(x)$ denoting the Hermite polynomials in monic normalisation. We recall that the $n$-th Hermite polynomials is an even (odd) function when $n$ is even (odd); this is due to the reflection symmetry of the Gaussian weight about the origin. A similar phenomenon is present for our product generalisation. In order to see this, we use an alternative form of the kernel (3.1). We have [17, 11]

[TABLE]

where $B_{n}^{(M)}=(b_{i,j}^{(M)})_{i,j=0}^{n-1}$ is the $n$-th bi-moment matrix constructed from the bi-moments

[TABLE]

Simon [53] refers the inverse moment matrix representation (3.5) as the ABC (Aitken–Berg–Collar) theorem.

In the following, it will be useful to extend the concept of odd and even moments to bi-moments. We say that the bi-moments are odd (even) when $k+\ell$ is odd (even). Now, using that $g_{\ell}^{(M)}(-x)=(-1)^{\ell}g_{\ell}^{(M)}(x)$, we see that the bi-moments satisfy

[TABLE]

which implies that all odd moments are equal to zero. In other words, the entries in the bi-moment matrix vanishes in a chequerboard pattern. Thus, by reordering rows and columns we may write the bi-moment matrix in a block diagonal form $B_{n}^{(M)}\mapsto\operatorname*{diag}(B_{n}^{(M,\text{even})},B_{n}^{(M,\text{odd})})$ with $B_{n}^{(M,\text{even})}=(b_{2k,2\ell}^{(M)})_{k,\ell}$ and $B_{n}^{(M,\text{odd})}=(b_{2k+1,2\ell+1}^{(M)})_{k,\ell}$. Using this reordering in the sum (3.5), we see that the kernel splits into two parts

[TABLE]

where

[TABLE]

Here, the latter equality in both (3.9) and (3.10) follow from comparison with (3.1). Finally, we note that

[TABLE]

Thus, in the following we can restrict our attention to kernels with an even subscript.

3.2 Bi-orthogonal functions

We are now ready to write down the bi-orthogonal functions for our product ensemble. As explained in the previous subsection, it is convenient to consider functions of odd and even degree separately.

Proposition 7**.**

The ensemble defined by Theorem 4 is bi-orthogonalised by

[TABLE]

with notation as above (recall that $\nu_{0}=0$).

There are several different approaches to prove Proposition 7. Here, we will present a method which emphasizes the relation to the Hermite polynomials (see [31, Prop. 3.5] for the same method applied to the product ensemble of Laguerre type). In order to use this approach, it is convenient to first calculate the bi-moments.

Lemma 8**.**

The bi-moments are given by

[TABLE]

for $k+\ell$ even and zero otherwise. The bi-moment determinant is

[TABLE]

Proof.

We insert the contour integral representation of the weight functions (2.10) into the expression for the bi-moments (3.6), then we see that the even moments are

[TABLE]

The integrals in this expression can be recognised as a combination of a Mellin and an inverse Mellin transform, which yields (3.13).

In order to evaluate the bi-moment determinant (3.14), we note that

[TABLE]

This completes the proof, since the $M=0$ case is the well-known Hermite (or GUE) case. ∎

Proof of Proposition 7.

The bi-orthogonal functions may be expressed by means of their bi-moments exactly as orthogonal polynomials through their moments. Thus, we have

[TABLE]

Furthermore, we have $h_{n}=D_{n}^{(M)}/D^{(M)}_{n-1}$.

The constants $h_{n}$ are immediate from the above and (3.14). Thus, it remains only to find the bi-orthogonal functions. To do so, we first note that

[TABLE]

This observation is important, since we know that the monic Hermite polynomials (with respect to the weight $e^{-x^{2}}$) are given by

[TABLE]

It follows that the expressions for the bi-orthogonal function $p_{n}(x)$ and $\phi_{n}(x)$ can be found using the known expressions for the Hermite polynomials and then making substitutions

[TABLE]

respectively. We recall that

[TABLE]

which makes it a straightforward exercise to verify the proposition. ∎

Remark 9**.**

We recall that the bi-orthogonal functions also can be obtained from the characteristic polynomial using that

[TABLE]

with $\langle\cdots\rangle$ denoting the matrix average and $z\in\mathbb{C}\setminus\mathbb{R}$. The first relation allows for an alternative method to calculate $p_{N}(x)$; see e.g. [23].

4 Integral representations and correlation kernels

The explicit expressions for the bi-orthogonal functions given by Proposition 7 allow us to write down an explicit form for the correlation kernel by insertion in (3.1). However, this formulation of the kernel is not optimal for asymptotic analysis. For this reason, in this section we will provide integral representations of the bi-orthogonal functions as well as the kernel.

Proposition 10**.**

The bi-orthogonal functions given by Proposition 7 have integral representations

[TABLE]

where the contour $\Sigma$ encloses the integers $0,1,\ldots,n$ in the negative direction and $-\frac{1}{2}<c<0$. We recall that $p_{2n}(x)$ and $\phi_{2n}(x)$ are even functions, while $p_{2n+1}(x)$ and $\phi_{2n+1}(x)$ are odd functions.

Proof.

The integrands in (4.1) and (4.2) has $n+1$ simple poles located at $0,1,\ldots,n$, thus the series representations in Proposition 7 follow upon a straightforward application of the residue theorem.

In order to find the integral representation of the bi-orthogonal functions $\phi_{n}(x)$, we first note that the weight functions can be written as

[TABLE]

this is easily seen starting from the recursive definition (2.7). Now, using (4.5) in the expression for $\phi_{n}(x)$ (cf. Proposition 7), we see that

[TABLE]

In other words, the bi-orthogonal functions $\phi_{n}(x)$ are an integral transform of the Hermite polynomials with respect to a Meijer $G$-function as integral kernel. The Hermite polynomial can itself be expressed as a Meijer $G$- or Fox $H$-function (see [45, Sec. 1.8.1.]) and the remaining integral is well-known from the literature [45, Sec. 2.3.]. ∎

Proposition 11**.**

Integral representations of the even and odd kernels are

[TABLE]

with $-1<c<0$ and the contour $\Sigma$ is chosen such that it encircles $0,1,\ldots,2n-1$ in the negative direction with $\textup{Re}\{s\}>c$ for all $s\in\Sigma$.

Proof.

As the proofs for the odd and even kernels are almost identical, we provide only the proof for the even case. The odd case is easily verified by the reader.

It follows the definition of the even kernel (3.9) together with contour integral representation of the bi-orthogonal functions from Proposition 10, that

[TABLE]

Following similar steps as in [41], we note the sum allows a telescopic evaluation. This gives

[TABLE]

Here, the integrand on the second line is zero as it has no poles encircled by the contour $\Sigma$ and, thus

[TABLE]

Finally, the proposition follows by a change of variables $s\mapsto s/2$ and $t\mapsto t/2$. ∎

The above integral representations for the bi-orthogonal functions and the kernel are probably the most convenient form for further asymptotic analysis, as we will see in Section 5. However, it is also often helpful to express these formulae in terms of special functions as (for example) it allows for use of pre-defined mathematical software. Furthermore, such reformulations often guide us to recognise patterns which otherwise would have been left unseen.

The integral representations for the bi-orthogonal functions given by Proposition 10 can also be recognised as several different types of special function; this includes generalised hypergeometric, Meijer $G$-, and Fox $H$-functions. Here, we will restrict ourselves to their Meijer $G$-function formulation.

Let us first consider the bi-orthogonal polynomials which may be written as

[TABLE]

It is worth comparing these polynomials with the polynomial found in the study of the Laguerre-like matrix product (1.10). Akemann et al. [4] found that in this case the bi-orthogonal polynomial is given by

[TABLE]

It is clear that the two families of polynomials are related as

[TABLE]

with

[TABLE]

respectively. This is a generalisation of the relation between Hermite and Laguerre polynomials. Recall that

[TABLE]

where $\tilde{H}_{n}(x)$ and $\tilde{L}_{n}^{(\alpha)}(x)$ denote the Hermite and Laguerre polynomials in monic normalisation.

Likewise, the (non-polynomial) bi-orthogonal functions may be written as

[TABLE]

Again, we want to compare to the formula in [4] which this time reads

[TABLE]

Evidently, we have the following relations

[TABLE]

with (4.14) as before. Yet again, this is a generalisation of the relation between Hermite and Laguerre polynomials. In the simplest case the relations (4.19) reduces to

[TABLE]

where $w_{\text{H}}(x)=e^{-x^{2}}$ and $w_{\text{L}}^{(\alpha)}(x)=x^{\alpha}e^{-x}$ are the Hermite and Laguerre weight functions.

It is, of course, well-known that there are relations between ensembles with reflection symmetry about the origin and ensembles on the half-line (albeit explicit formulae may be elusive). A general description of such relations in the Muttalib–Borodin ensembles can be found in [22].

5 Scaling limits at the origin in product and Muttalib–Borodin ensembles

With the integral representations of the correlation kernels established by Proposition 11, we can turn to a study of asymptotic properties. Perhaps the most interesting scaling regime is that of the local correlations near the origin, referred to as the hard edge when the eigenvalues are strictly positive. For other product ensembles [41, 21, 35, 24], it has been observed that correlations at the hard edge is determined by the so-called Meijer $G$-kernel, which generalises the more familiar Bessel kernel. Below, we will see that the Meijer $G$-kernel appears once again, but this time involving a sum.

Theorem 12**.**

Let $K_{n}(x,y)=K_{n}^{\textup{even}}(x,y)+K_{n}^{\textup{odd}}(x,y)$ with the even and odd kernels given by Proposition 11. For $x,y\in\mathbb{R}\setminus\{0\}$ and $\nu_{1}\ldots,\nu_{M}$ fixed, the microscopic limit near the origin is

[TABLE]

with

[TABLE]

where $-1<c<-1/2$ and $\Sigma$ encloses the non-negative half-line in the negative direction starting and ending at $+\infty$ such that $\textup{Re}\{s\}>c$ for all $s\in\Sigma$.

Proof.

We only consider the even kernel in Proposition 11 since the odd case is very similar. After rescaling we rewrite the integral representation of the even kernel in Proposition 11 as

[TABLE]

with

[TABLE]

For any fixed $t\in c+i\mathbb{R}$ and $s\in\Sigma$, using [49, eq. 5.11.13] we see

[TABLE]

Formally, substituting (5.6) in (5.4) gives (5.2). To proceed rigorously, we need to verify a condition for the exchange of limit and integration. For this purpose, we will proceed to find two dominated functions respectively corresponding to $1/|f_{n}(t)|$ and $|f_{n}(s)|$.

First, using [49, eq. 5.11.13] we have for sufficiently large $n$

[TABLE]

Second, we require an upper bound for $|f_{n}(s)|$. Noting the asymptotic expansion, that as $z\rightarrow\infty$ in the sector $|\mathrm{arg}(z)|\leq\pi-\delta$ (with $0<\delta<\pi$)

[TABLE]

it is easy to see that for a given $y_{0}>0$ we can choose the contour $\Sigma=\Sigma_{l}\cup\Sigma_{r}$ with

[TABLE]

Thus, we get from (5.8) and the boundedness of $\Gamma(-s/2)$ over $\Sigma_{l}$ that for large $n$ there exists a constant $C_{1}=C_{1}(y_{0})>0$ such that

[TABLE]

In order to estimate $f_{n}(s)$ with $s\in\Sigma_{r}$, we use the integral representation

[TABLE]

where $\mathcal{C}_{0}$ is a counter-clockwise path which begins and ends at $1$ and encircles the origin once; see e.g. [49, eq. 5.12.10]. Note that we choose $(-u)^{-1-s/2}=e^{-(1+s/2)\log(-u)}$ with $-\pi<\mathrm{arg}(-u)<\pi$. Change $u$ by $u/n$ and deform the resulting contour into the path which starts from $n$, proceeds along the (upper) real axis to 1, describes a circle of radius one counter-clock round the origin and returns to $n$ along the (lower) real axis. That is,

[TABLE]

Let $s=v\pm iy_{0},v>1$. On the unit circle of the $u$-integral above write $-u=e^{i\theta}$. Then we easily obtain for $n\geq 1$

[TABLE]

On the upper and lower real axis, we have

[TABLE]

Using the simple fact $|\sin\frac{\pi s}{2}|\geq|\sinh\frac{\pi}{2}\mathrm{Im}(s)|$, combination of (5.13) and (5.14) shows that there exists a constant $C_{2}=C_{2}(y_{0})>0$ such that

[TABLE]

Together with (5.10) this gives us a bound $C>0$, that is, for large $n$

[TABLE]

Finally, use (5.8) and the asymptotic formula that as $y\rightarrow\pm\infty$

[TABLE]

with bounded real value of $x$ (see [49, eq. 5.11.9]), it is easy to conclude that the function of variables $s$ and $t$

[TABLE]

is integrable along the chosen contours, whenever $-1<c<-1/2$. Here we emphasize that the assumption $-1/2\leq c<0$ does not ensure the convergence in the special case $M=0$ while for $M\geq 1$ it can be relaxed to $-1<c<0$ as in Proposition 11.

With this, combing (5.7) and (5.16), we have indeed justified the interchange of limit and integrals for every $M$ by the dominated convergence theorem, which completes the proof. ∎

For a comparison with other known results, it is useful to rewrite the hard edge correlation function of Theorem 12 in terms of Meijer $G$-functions. Using that

[TABLE]

we see that the even and odd kernel can be written as

[TABLE]

respectively. We recall that the so-called Meijer $G$-kernel is given by [41]

[TABLE]

with $x,y>0$. We note that this kernel is single-sided ($x,y\in\mathbb{R}_{+}$) while the kernel from Theorem 12 is double-sided ($x,y\in\mathbb{R}\setminus\{0\}$). However, it is also evident that our new kernel may be re-expressed in terms of the Meijer $G$-kernel. We have

[TABLE]

with

[TABLE]

respectively. Thus, the random product matrix (1.1) provides yet another appearance of the Meijer $G$-kernel; albeit this time in a double-sided version. For graphical representation of the Meijer $G$-kernel we refer to [31, Fig. 3.2], which shows plots of the local density (i.e. the kernel with $x=y$) for different values of $M$.

A double-side hard edge scaling limit near the origin is also present in the Hermite Muttalib–Borodin ensemble. In this case the kernel is found to be [11]

[TABLE]

where

[TABLE]

with $J_{a,b}(x)$ denoting Wright’s Bessel function. In the case relevant to us (1.2), we also have $\theta=2M+1$. Furthermore, it is known from [40] that the kernel (5.26) is a Meijer $G$-kernel whenever $\theta$ is a positive integer. In particular, we have

[TABLE]

where the Meijer $G$-kernel on the right-hand side has indices

[TABLE]

and as always $\nu_{0}=0$. It follows from (5.25) and (5.26) that the hard edge correlations for the Hermite Muttalib–Borodin ensemble with appropriately chosen parameters may be expressed in terms of the Meijer $G$-kernel in a similar fashion as done for the product ensemble above. We note that the choice of variables in (5.26) should be compared to the change of variables (2.17) performed in the derivation of the asymptotic reduction (2.18).

It is worth verifying consistency of the simplest scenario of $M=0$. When $M=0$ our matrix ensemble (1.1) reduces to the GUE, hence the kernel given by Theorem 12 must reduce to the sine kernel for $M=0$. To see this, we use

[TABLE]

It follows that

[TABLE]

which upon insertion into (5.1) indeed reproduces the sine kernel.

In the end of this section, let us emphasize that there also exists a contour integral representation of the limiting kernel in Theorem 12 which combines the odd and even into a single formula.

Proposition 13**.**

With the same notation as in Theorem 12, the limiting kernel at the origin can be rewritten as

[TABLE]

where the kernel on the right-hand side is defined as

[TABLE]

with $C_{R}$ denoting a path in the right-half plane from $-i$ to $i$.

Proof.

Using Euler’s reflection formula and duplication formula for the gamma function, we see that

[TABLE]

where

[TABLE]

In order to proceed, we will consider the cases $xy<0$ and $xy>0$ separately. For $xy<0$, it is seen that

[TABLE]

Now (5.32) can be obtained using the integral representation

[TABLE]

with the contour $C_{R}$ as above, together with the definition of Meijer $G$-function. For $xy>0$, we note that

[TABLE]

The $s$-variable integrand in the second part has no pole within the contour $\Sigma$. Thus, the problem reduces to the proven situation. ∎

The simplest non-trivial case is $M=1$. Here, we get

[TABLE]

with the modified Bessel functions $I_{\nu}$ and $K_{\nu}$, which follows immediately from the fact that

[TABLE]

6 Global spectra in product and Muttalib–Borodin ensembles

The study of the scaling limit at the origin in the previous section introduces a scale in which the average spacing between eigenvalues is of order unity. A very different, but still well-defined, limiting process is the so-called global scaling regime. In this regime the average spacing between eigenvalues tends to zero in such way that the spectral density tends to a quantity $\rho(x)$ with compact support $I\subset\mathbb{R}$ and $\int_{I}\rho(x)dx=1$. Here $\rho(x)$ is referred to as the global density. Throughout this section, the indices $\nu_{1}\ldots,\nu_{M}$ are kept fixed.

For the Laguerre Muttalib–Borodin ensemble specified by the density (1.12) the global scaling limit corresponds to a change of variables $x_{j}\mapsto nx_{j}$. Introducing the further change of variables $x_{j}\mapsto Mx_{j}^{M}$, the global density is known to be the so-called Fuss–Catalan density with parameter $M$ [27]. It can be specified by the moment sequence

[TABLE]

These are the Fuss–Catalan numbers (the Catalan numbers are the case $M=1$).

Now, consider the product of $M$ standard complex Gaussian random matrices. Consistent with the discussion in Section 1.2, the corresponding global density is again the Fuss–Catalan density with parameter $M$ [46, 7, 10, 48].

It is known that the Fuss–Catalan density, $\rho^{(M)}_{\text{FC}}(x)$ say, can also be characterised as the minimiser of the energy functional

[TABLE]

with $L=(M+1)^{M+1}/M^{M}$; see [16, 23, 25]. Note that the energy functional (6.2) relates to (1.12) through the aforementioned change of variables. Similarly, the energy functional corresponding to (1.2) is

[TABLE]

with $\theta=2M+1$. We note that changing variables $x^{2}\mapsto x$ and $y^{2}\mapsto y$, then setting $\tilde{\rho}(x)=x\rho(x^{2})$ reduces (6.3) to (6.2) with $L=\tilde{L}^{2}$. Thus, the minimiser in (6.3) is given in terms of the Fuss-Catalan density

[TABLE]

and is symmetric about the origin.

As an illustration, let us consider the simplest case, $M=1$. The Fuss–Catalan density becomes the celebrated Marčenko–Pastur density,

[TABLE]

The formula (6.4) then gives the standard result (see e.g. [50]) that the energy functional

[TABLE]

is minimised by

[TABLE]

which is Wigner’s semi-circle law.

It has been demonstrated in Section 2, that the energy function implicit in (6.3) underlies the eigenvalue distribution of the random matrix product (1.1). Thus, we can anticipate that after appropriate scaling the global density for the product ensembles is given by (6.4). A direct proof of this can obtained through a number of different strategies. We consider first a method based on the characteristic polynomial.

In terms of the global scaled variables, the key equation relating the averaged characteristic polynomial to the global is the asymptotic formula [27]

[TABLE]

where $\tilde{G}_{M}(z)$ is the Stieltjes transform of the global spectral density,

[TABLE]

Following the strategy first used in [23], the formula (6.8) leads to a characterisation of the Stieltjes transform (6.9), upon realising that the characteristic polynomial satisfy a linear differential equation.

Proposition 14**.**

Consider a matrix product (1.1) with even matrix dimension, $n=2N$. Let

[TABLE]

denote the characteristic polynomial. Then (6.10) is a solution to the $(2M+2)$-th differential equation,

[TABLE]

with asymptotic boundary condition $f(z)\sim z^{n}$ for $|z|\to\infty$.

Proof.

The characteristic polynomial (6.10) is identical to the bi-orthogonal polynomial $p_{2N}(z)$. As shown earlier, this polynomial is proportional to a Meijer $G$-function (4.11). It is well-known that such Meijer $G$-functions satisfy the differential equation (6.11). The asymptotic boundary condition follows trivially, since $f(z)$ is a monic polynomial. ∎

Changing variables $z\mapsto n^{M+\frac{1}{2}}\hat{z}/\sqrt{2}$ in (6.11) and using that [23]

[TABLE]

to leading order in $n$, we see that for large $n$ the differential equation (6.11) reduces to the algebraic equation (see e.g. [9] for $M=1$)

[TABLE]

with asymptotic condition $\tilde{G}_{M}(z)\sim 1/z$ as $|z|\to\infty$. This equation is to be compared to the algebraic equation satisfied by the Stieltjes transform of the Fuss–Catalan density,

[TABLE]

see e.g. [23]. With $z\mapsto z^{2}$ and $M\mapsto 2M+1$ and setting $\tilde{G}_{M}(z)=zG_{2M+1}(z^{2})$, we see that (6.14) reduces to (6.13). This prescription is equivalent to  (6.4), thus verifying this formula as the evaluation of the global density.

The same result can also be obtained using free probability techniques. To see this, we need some additional notation. Let $a$ be a non-commutative random variable with distribution $d\mu(x)=\rho(x)dx$. The Stieltjes transform $G_{a}(z)$ of the variable $a$ is defined analogous to (6.9). The $S$-transform is defined as

[TABLE]

Now assume that $a$ and $b$ are two freely independent non-commutative random variables and that the Stieltjes transform $G_{b}(z)$ satisfies a functional equation $P(z,G_{b}(z))=0$. It is known [48] that under these conditions the Stieltjes transform $G_{ab}(z)$ of the product $ab$ satisfies

[TABLE]

Moreover, we know that if $a$ is given by the free normal distribution (i.e. Wigner’s semi-circle) and $b$ is given by the free Poisson distribution (i.e. Marčenko–Pastur), then

[TABLE]

We can now use that the limiting distributions for the GUE and the Wishart ensemble are the free normal and the free Poisson, respectively. Thus, using (6.16) $M$ times, we see that our product (1.1) indeed gives rise to the the functional equation (6.13).

It is also possible to construct a parametrisation of the global density in terms of elementary functions based on the polynomial equation (6.13). With

[TABLE]

we have

[TABLE]

see e.g. [23]. We remark that it follows that the singularity at the origin blows up like

[TABLE]

as $x_{0}\to 0$.

7 Conclusion and outlook

In this paper, we have shown that it is possible to construct a Hermitised random matrix product for which the eigenvalues form a determinantal point process on the entire real line with an explicit kernel. This is a fundamental new contribution to the study of random matrix product ensembles, since all previous exactly solvable models of this type have had eigenvalues restricted to the positive half-line. Furthermore, we have argued that this Hermitised product ensemble can be considered a natural generalisation of the classical Hermite ensemble (i.e. GUE) in similar way as the squared singular values of matrix products with Gaussian matrices [5, 4] and truncated unitary matrices [35] can be considered generalisations of the Laguerre and Jacobi ensembles, respectively. To this point, we have shown that the joint eigenvalue PDF reduces asymptotically to the Muttalib–Borodin ensemble of Hermite type.

On another front, we have shown that the local scaling limit near the origin is described by a two-sided generalisation of the so-called Meijer $G$-kernel [41]. This two-sided kernel reduces to the sine kernel in the simplest case. We have also seen that the global density can be found explicitly and that it is expressed in terms of the so-called Fuss–Catalan distribution in a simple manner. Our result relies on an explicit double contour integral formulation of the correlation kernel (Proposition 11). It is worth stressing that we could make full use of this double contour integral formulation and give an analytical proof of the global density. In fact, following almost exactly the same steps introduced in [44], it can be proven that the sine-kernel arises in the bulk and that Airy-kernel arises at the soft edge; cf. the proof of Theorems 1.1 and Theorem 1.3 as well as Remark 2 in [44]. The full details are beyond the scope of this paper, so let us only mention that a basic starting point is to approximate the integrand by elementary functions and rewrite the kernel, say the even part, as

[TABLE]

where the phase function is given by

[TABLE]

Hence, the saddle point equation $g^{\prime}(z)=0$ is exactly expressed through the equation (6.13). We stress that the above parametrisation representation plays a key role in the proof of the sine-kernel via the steepest decent method.

Finally, we emphasize that our construction of a Hermitised random product ensemble is based on a matrix transformation which maps the space of polynomial ensembles onto itself (Theorem 1). This type of matrix transformation are important since they preserve exact solvability. Our proof of Theorem 1 is applicable to the Hermitised product ensemble multiplied by a Gaussian matrix, crucially with the help of the hyperbolic HCIZ integral over the pseudo-unitary group. However, it would be interesting to see whether this could be extended to the product ensemble multiplied by other types of random matrices, say, truncated unitary matrices. For this, a possible way is to first extend the matrix integral formula stated in [35, Theorem 2.3] from the unitary group to the pseudo-unitary case, and then perform the same steps as in Section A.3. This will be an interesting and challenging problem for us.

Acknowledgements

We thank S. Kumar for useful discussions on his paper and M. Kieburg for comments on a first draft. We acknowledge support by the Australian Research Council through grant DP170102028 (PJF), the ARC Centre of Excellence for Mathematical and Statistical Frontiers (PJF,JRI), and partially by ERC Advanced Grant #338804, the National Natural Science Foundation of China #11771417, the Youth Innovation Promotion Association CAS #2017491, the Fundamental Research Funds for the Central Universities #WK0010450002, Anhui Provincial Natural Science Foundation #1708085QA03 (DZL). D.-Z. Liu is particularly grateful to László Erdős for funding his one-year stay at IST Austria.

Appendix A Three different proofs of Theorem 1

This appendix contains three separate proofs of Theorem 1. Each proof has its own merits and provides a different perspective on the matrix transformation.

The first proof is based on the method of additive rank-one deformations. A benefit of this method is that it avoids group integrals of HCIZ type, and therefore might be a more suitable starting point for generalisations to studies of real or quaternion matrices. The second proof is based on a theorem by Forrester and Rains [26] that gives the eigenvalue density of a Hermitised matrix product by means of an inverse double-sided Laplace transform. This idea is closely related to the spherical transforms used in the context of other matrix transformations [34, 39]. The third and final proof uses a generalisation of the HCIZ integral previously studied by Fyodorov [28, 29].

The three proofs will be given in the three subsections below. However, before we start it is worth noting the following reduction result related to Theorem 1.

Remark 15**.**

Suppose $n=N$, and consider the limit $a_{n_{0}+1}\to 0^{+}$. Recalling the ordering (1.5), inspection of (1.4) shows that its leading contribution comes from a Laplace expansion via the top left entry of the second determinant. This entry is in turn significant only for $x_{n_{0}+1}\to 0^{+}$, telling us that to leading order (1.4) with $n=N$ and $a_{n_{0}+1}\to 0^{+}$ is equal to

[TABLE]

After relabelling, (A.1) is equivalent to (1.4) with $n=N-1$ times a Dirac delta function corresponding to an eigenvalue at zero. Repeating this limiting procedure a total of $N-n$ times shows that (1.4) in the case $n=N$ reduces to the general $n$ case.

A.1 First proof: Recursive structure using additive rank-one deformations

In this section, we prove Theorem 1 by induction. The induction step will be constructed using the method of rank-one deformations. For the reader’s convenience, we start by providing an outline of the proof; details will follow in the steps i), ii), and iii) below.

We first need some additional notation. Let $G^{(p)}$ denote the $p\times N$ matrix consisting of the first $p$ rows of the $n\times N$ complex Gaussian matrix $G$ and let $A^{(p)}=\operatorname*{diag}(a_{1},\dots,a_{p})$. Define

[TABLE]

We see that $X^{(p)}$ is an $N\times N$ matrix and that $X^{(n)}=X$. Moreover, $X^{(p)}$ has rank (less than or equal to) $p$ and we therefore know that it has (at most) $p$ non-zero eigenvalues; we will denote these eigenvalues $\lambda_{k}^{(p)}$ ($k=1,\ldots,p$). The crucial observation is that $X^{(p)}$ is an additive rank-one deformation of $X^{(p-1)}$ for $p>1$. More precisely, we have

[TABLE]

where $\vec{x}$ is an $N\times 1$ column vector with standard complex Gaussian entries. We will see below that if the eigenvalues of $X^{(p-1)}$ are known, then the rank-one deformation (A.3) can be used to find the conditional PDF for the eigenvalues of $X^{(p)}$. Let us denote this conditional PDF by

[TABLE]

It is clear that if the PDF $P^{n_{0}}_{1}(a_{1};\lambda_{1}^{(1)})$ is known, then $P^{n_{0}}_{p}(\{a_{j}\};\{\lambda_{j}^{(p)}\})$ can be constructed recursively using

[TABLE]

for $p=1,\ldots,n$ and a suitable integration domain $D$. Thus, our proof can be divided into three steps:

• i)

Use the additive rank-one deformation (A.3) to find the conditional PDF (A.4).

• ii)

Use the conditional PDF (A.4) together with the recursion (A.5) to show that if $P^{n_{0}}_{p-1}(\{a_{j}\};\{\lambda_{j}^{(p-1)}\})$ is given by (1.4) with $n=p-1$, then $P^{n_{0}}_{p}(\{a_{j}\};\{\lambda_{j}^{(p)}\})$ is given by (1.4) with $n=p$.

• iii)

Show that $P^{n_{0}}_{p=1}(a_{1};\lambda_{1}^{(1)})$ is given by (1.4) with $n=1$.

We will look at these three steps separately below.

i)

We want to find the eigenvalues of $X^{(p)}$ using the rank-one deformation (A.3) assuming that the eigenvalues of $X^{(p-1)}$ are known. The matrix $X^{(p-1)}$ has (at most) $p-1$ non-zero eigenvalues. We assume that these eigenvalue are pairwise distinct and ordered as

[TABLE]

i.e. $X^{(p-1)}$ has $p-n_{0}-1$ positive eigenvalues if $p>n_{0}$ and no positive eigenvalues if $p\leq n_{0}$.

If $p<N$ then the matrix $X^{(p)}$ must have an eigenvalue equal to zero with multiplicity $N-p$ (or higher). The remaining $p$ eigenvalues will be random variables which are non-zero and have multiplicity one (almost surely), since the deformation (A.3) includes a Gaussian vector $x$. Thus, we know from [26] that the eigenvalues of $X^{(p)}$ are given as solutions to the secular equation

[TABLE]

where, with $\Gamma[\alpha,\beta]$ denoting a gamma-distibuted variable with shape parameter $\alpha$ and rate parameter $\beta$, each $q_{j}$ is a random variable given by

[TABLE]

Furthermore, the eigenvalues of $X^{(p)}$ must be interlaced with the eigenvalues $X^{(p-1)}$, i.e. interlaced with $\{0\}\cup\{\lambda_{k}^{(p-1)}\}_{k=1}^{p-1}$. This interlacing may be verified by sketching the plot of the secular equation (A.7) as a function of $\lambda$. Moreover, we note that whether the interlacing starts from the left or from the right depends on whether $a_{p}$ is negative or positive, or equivalently on whether $p\leq n_{0}$ or $p>n_{0}$, cf. (1.4). For $p\leq n_{0}$ we have the interlacing

[TABLE]

while for $p>n_{0}$ we have the interlacing

[TABLE]

Subject to these interlacings, we read off from [26, Cor. 3] that the corresponding conditional PDF (A.4) is given by

[TABLE]

with indices $1\leq i,j\leq p$ and $1\leq k,\ell\leq p-1$.

ii)

We can now turn to the recursive formula (A.5). The probability density (1.4) with $n=p-1$ is given by

[TABLE]

Furthermore, we know from step i) that the conditional PDF is given by (A.11) and that the integration domain is given by either (A.9) or (A.10) depending on whether $p\leq n_{0}$ or $p>n_{0}$. Considering the case $p\leq n_{0}$, we have the recursion

[TABLE]

We note that there is only one determinant since $p<n_{0}$. Let us focus on the integral on the last line in (A.13). We see that

[TABLE]

with integration domain on the right-hand side on the first line given by (A.9). The first equality in (A.1) follows by shifting the integration inside the determinant, while the second equality follows by a standard row manipulation. Performing the integral within the determinant on the last line (A.1), we see that

[TABLE]

where the last inequality can be understood by applying elementary row operations of adding multiples of the first row to the rows below so as to get zero entries in the final column, expect for the first entry, then Laplace expanding by that entry.

Finally using this evaluation of the integral from (A.13) verifies the recursion for $p\leq n_{0}$. The verification for $p>n_{0}$ follows the same lines.

iii)

It only remains to show that $P^{n_{0}}_{p=1}(a_{1};\lambda_{1}^{(1)})$ is given by (1.4) with $n=1$. There are two cases $n_{0}=0$ and $n_{0}=1$, but they are both immediate since we are considering scalars.

A.2 Second proof: Limit of inverse Laplace transform expression

This second proof of our main theorem starts by looking at the eigenvalues of a more general matrix

[TABLE]

where $G$ be an $n\times N$ complex Gaussian random matrix, while $A$ and $B$ are Hermitian matrices. Due to unitary invariance we can in fact choose $A$ and $B$ to be diagonal, say $A=\operatorname*{diag}(a_{1},\ldots,a_{n})$ and $B=\operatorname*{diag}(b_{1},\ldots,b_{N})$, without loss of generality.

It is evident that this eigenvalue problem reduces to that of Theorem 1, when

[TABLE]

This observation is crucial, since it was shown by Forrester and Rains [26, Thm. 6] that (assuming the eigenvalues of $A$ and $B$ are pairwise distinct) the eigenvalue PDF for the matrix (A.16) can be written as

[TABLE]

where $\operatorname{ePDF}(M)$ is a short-hand notation for the eigenvalue density for a random matrix $M$ and $\mathcal{L}^{-1}$ denotes the inverse two-sided Laplace transform,

[TABLE]

In some sense, this result is more general than the one we are trying to prove, but it is also far less explicit and therefore less useful for our purposes. Thus, the strategy to prove Theorem 1 presented in this subsection is to show that given (1.3) then (A.18) reduces to (1.4) in the limit (A.17). Here (1.4) refers to the PDF with the $N$ non-zero eigenvalues. However, we know from Remark 15 that the $n=N$ case of (1.4) implies the general $n$ case. It is therefore sufficient for us to set $n=N$ in (A.18) and to show that in the limit (A.17) the case $n=N$ of (1.4) appears.

For this purpose, we begin by noting that with $f(s)=\det(\mathbb{I}+As)^{-1}$, the limit $\tau\to 0^{+}$ in (A.19) can be taken inside the integral, telling us that

[TABLE]

Our evaluation of this contour integral follows the standard procedure of first closing the contour and then using the residue theorem. Due to the inequalities (1.5), we see that the contour must be closed in the positive half-plane for $i=1,\ldots,n_{0}$, which according to the inequalities (1.4) picks up contributions from $n_{0}$ simple poles located at $-1/a_{1},\ldots,-1/a_{n_{0}}$. For $i=n_{0}+1,\ldots,N$ the contour must be closed in the negative half-plane picking up contributions from the remaining $N-n_{0}$ poles located at $-1/a_{n_{0}+1},\ldots,-1/a_{N}$. After substituting this straightforward evaluation of the contour integral into the last determinant in (A.18), the limit (A.17) may be found by successive use of L’Hôpital’s rule. This yields

[TABLE]

The latter determinant in this expression may be simplified considerably by noting that the only factor inside determinant which depends on both index $i$ (through $k_{i}$) and index $j$ is $a_{k_{i}}^{N-j}$, while the only factor depending on both index $i$ and index $k_{i}$ is $e^{-x_{i}/a_{k_{i}}}$. Thus, upon expansion and reordering of products we see that

[TABLE]

Thus, the eigenvalue PDF becomes

[TABLE]

In order to evaluate the sums on the second line in (A.23), we note that the indices satisfy $k_{1},\dots,k_{n_{0}}<k_{n_{0}+1},\dots,k_{N}$ and that the determinant is anti-symmetric in both $\{k_{1},\ldots,k_{n_{0}}\}$ and $\{k_{n_{0}+1},\ldots,k_{N}\}$. This allow us to write

[TABLE]

The first determinant on the right-hand side in (A.24) and the first determinant in (A.23) are both Vandermonde determinants, so the eigenvalue PDF (A.23) becomes

[TABLE]

which we recognise as the desired statement (1.3) with $n=N$.

A.3 Third proof: Matrix integral over the pseudo-unitary group

For this third proof of Theorem 1 we will again restrict our attention to the case $n=N$. The Gaussian matrix $G$ specified in Theorem 1 has distribution

[TABLE]

where $(dG)$ is the Lebesgue measure on the space of complex $N\times N$ matrices. It is a standard result from random matrix theory that the positive semi-definite Hermitian matrix $\tilde{W}=GG^{\dagger}$ is distributed according to

[TABLE]

where $(d\tilde{W})$ is the Lebesgue measure on the space of Hermitian matrices subject to the constraint that $\tilde{W}$ is positive semi-definite. This may be seen by decomposing the matrix $G$ using a polar decomposition, i.e. $G=U\tilde{W}^{1/2}$ with $U$ unitary. Making this change of variables and integrating over the unitary degrees of freedom contribute an extra factor, $2^{-N}{\rm vol}\,U(N)$, to the normalisation; see e.g. [20].

The proof presented in this section is based on an integration formula for the pseudo-unitary group with a pseudo-metric tensor $\eta$ determined according to the number of positive (negative) eigenvalues of the matrix $A$. More precisely, with $A$ as specified in Theorem 1, we define

[TABLE]

such that $A=A_{+}\eta=\eta A_{+}$. In fact, we have the more general relation $A=A_{+}^{p}\eta A_{+}^{q}$ with $p+q=1$, since the matrix $A$ is assumed to be non-singular.

Next, introduce the matrix $Z=A_{+}^{1/2}\tilde{W}A_{+}^{1/2}$ where $\tilde{W}$ is an $N\times N$ matrix distributed according to (A.27). Since both $Z$ and $\tilde{W}$ are complex Hermitian random matrices, we have $(dZ)=(\det A_{+})^{N}(d\tilde{W})$, see e.g. [20, Eq. (1.35)]. Thus, using that $A_{+}$ is invertible we know from (A.27) that the distribution of $Z$ is equal to

[TABLE]

In terms of $Z$ we may define $\tilde{Z}=\eta Z=\eta A_{+}^{1/2}GG^{\dagger}A_{+}^{1/2}$. This matrix is important, since its eigenvalues are identical to those of

[TABLE]

Here, the right-hand side is the matrix of interest for Theorem 1 and the equality is a simple consequence of the definition (A.28). Moreover, since $(d\tilde{Z})=(dZ)$ (the action of $\eta$ on $Z$ is only to change the sign of some of the entries of $Z$), we read off from (A.29) that the distribution of $\tilde{Z}$ is equal to

[TABLE]

where it is further required that $\eta\tilde{Z}$ is positive semi-definite.

It is a known result that if $Z$ is a positive definite matrix, then the matrix $\eta Z$ has exactly $n_{0}$ negative eigenvalues and $N-n_{0}$ positive eigenvalues [52]. Furthermore, the matrix $\tilde{Z}$ can be diagonalised using a pseudo-unitary similarity transformation, i.e. there exists a matrix $V\in U(\eta)$ such that

[TABLE]

where $L=\operatorname*{diag}(x_{1},\dots,x_{N})$ is a real diagonal matrix. We recall that the pseudo-unitary group is defined as

[TABLE]

The Jacobian associated to the change of variables (A.32) is [52]

[TABLE]

where $(V^{-1}dV)$ is the Haar measure on the pseudo-unitary group. For the measures (A.33) to be in one-to-one correspondence, it is necessary to restrict the overall phase of each eigenvector, or equivalently require that $V\in U(\eta)/U(1)^{N}$.

Substituting (A.33) into (A.31) shows that the eigenvalue PDF of $\tilde{Z}$ (or equivalently of $G^{\dagger}AG$) is equal to

[TABLE]

We note that for $n_{0}=0$ and $n_{0}=N$, the signature $\eta$ becomes proportional to the identity and the group integral in (A.34) reduces to the well-known HCIZ integral (1.14). The generalisation of the HCIZ integral from an integral over the unitary group to an integral over the pseudo-unitary group (i.e. $0<n_{0}<N$) has been studied by Fyodorov [28, 29], who showed that

[TABLE]

with $K_{N,n_{0}}$ denoting an undetermined proportionality constant, and $A=\operatorname*{diag}(a_{1},\ldots,a_{N})$ and $B=\operatorname*{diag}(b_{1},\ldots,b_{N})$ denoting diagonal matrices subject to the constraints

[TABLE]

The constraints (A.36) must be included to ensure convergence of the group integral on the right-hand side in (A.35) for an $0<n_{0}<N$. This is necessary since the pseudo-unitary group is non-compact except for $n_{0}=0$ or $n_{0}=N$ in which case the aforementioned constraints may be ignored.

Now, using the integration formula (A.35) to evaluate the group integral in (A.34), we see that

[TABLE]

which agrees with the $n=N$ case of (1.4) provided that $K_{N,n_{0}}=\pi^{N(N-1)/2}$. We note that the proportionality constant is independent of $n_{0}$. Moreover, the cases $n_{0}=0$ and $n_{0}=N$ (where the group integral is over $U(N)/U(1)^{N}$) are consistent with the known proportionality constant from the HCIZ integral; recall that our choice of measure is not normalised to unity rather we have ${\rm vol}\,U(N)/U(1)^{N}=\pi^{N(N-1)/2}/\prod_{j=1}^{N}\Gamma(j)$.