Real eigenvalue statistics for products of asymmetric real Gaussian matrices

Introduction

where $\Delta(\{z_{p}\}_{p=1}^{m}):=\prod_{j<l}^{m}(z_{l}-z_{j})$ denotes the Vandermonde determinant and

(see also ). Integrating (1.1) over $\{\lambda_{l}\}\cup\{x_{j}+iy_{j}\}$ gives the probability $p_{N,k}$ that there are exactly $k$ real eigenvalues. The simplest case to compute is when $k=N$ and thus all eigenvalues are real, for which the probability was found to equal $2^{-N(N-1)/4}$.

A second example comes from quantum entanglement. Lakshminarayan considered the problem of quantifying when two-qubits $|\phi_{1}\rangle$ and $|\phi_{2}\rangle$ are an optimal pair, in the case that the states are chosen from a uniform distribution on the unit 3–sphere. The condition of being an optimal pair is known as particular inequalities for certain weighted inner products between the qubits. It was shown in that these can be interpreted as the condition for the probability that the random matrix product $X_{1}X_{2}$, with each $X_{i}$ a $2\times 2$ real Gaussian matrix, having all eigenvalues real, which was furthermore shown to be equal to $\pi/4$.

An intriguing effect was observed in the study , which seems to have escaped early notice. Thus, noting from the result of Edelman cited above that for a single real Gaussian $2\times 2$ random matrix the probability of all eigenvalues being real is equal to $2^{-1/2}$, while for a product of two independent Gaussian $2\times 2$ random matrices it is $\pi/4$, the fact that $2^{-1/2}<\pi/4$ led Lakshminarayan to investigate if the probability of all eigenvalues being real was an increasing function of the number of matrices in the product. Numerical simulation indicated that this is indeed the case, and further the probability that all eigenvalues are real tends to unity as the number of random matrices in the product tends to infinity. Evidence that this is also true for products of $d\times d$ real Gaussian matrices was given in , while the follow up work provided similar evidence for random matrices with independent non-Gaussian entries. A proof in the instance of the latter circumstance that the entries are all independent and identically distributed with a PDF containing an atom has recently been given in .

The appearance of the work coincided with the appearance of works containing other surprising advances relating to the eigenvalues of products of random matrices. Consider the random matrix product

where each $X_{i}$ is an $N\times N$ standard Gaussian matrix. In the case of complex entries, Akemann and Burda showed that the eigenvalues form a determinantal point process in the complex plane. This means that the $k$-point correlation function for the eigenvalues $\rho_{(k)}(z_{1},\dots,z_{k})$ is fully determined by a single function $K(w,z)$, referred to as the correlation kernel, according to

In the case of real entries, Forrester found a closed form expression for the probability that all eigenvalues are real.

To specify this latter result requires introducing the Meijer $G$-function

where $\gamma$ is an appropriate contour relating to the validity of the inverse Mellin transform formula. With $p_{N,k}^{P_{m}}$ denoting the probability that the random matrix product $P_{m}$ (1.3) has exactly $k$ real eigenvalues it was shown in that, for each $X_{i}$ a real Ginibre matrix, we have

A simple identity for the Meijer $G$-function — evident from the definition (1.5) — shows that $a_{j,k}$ is equal to the Meijer $G$-function occurring in . Moreover, these explicit formulas were used to prove that $p_{N,N}^{P_{m}}\to 1$ as $m\to\infty$. Extension to rectangular matrices where given in , while special arithmetic properties were shown to be present in the case $m=2$ .

A primary aim of the present paper is to extend this result to the calculation of $p_{N,k}^{P_{m}}$, for general $0\leq k\leq N$ with the same parity as $N$. The following theorem will be proved in section 3.3.

Consider the random matrix product (1.3), in which each $X_{i}$ is a real Ginibre matrix. Let

with $h_{j}=(2\sqrt{2\pi}\Gamma(2j+1))^{m}$, $a_{j,l}\ (l>0)$ given by (1.7) and $a_{j,0}=0$. For $N$ even, the probability $p_{N,2k}^{P_{m}}$ that exactly $2k$ eigenvalues are real is given by

and similarly for the $N$ odd case in terms of (1.10).

A formula closely related to Theorem 1 in the case $m=1$ was derived by Akemann and Kanzieper , and this working was soon after refined to obtain a formula equivalent to (1.9). Also for this case Forrester and Nagao gave a result more general than (1.9), applying to a real random matrix formed from a general linear combination of Gaussian symmetric and anti-symmetric matrices.

We now turn our attention to the other primary aim of our work. This relates to the statistical state formed by the eigenvalues of the product (1.3). In the complex case, it has been remarked that the statistical state is a determinantal point process. In the real case, it is known from the work of Ipsen and Kieburg that the eigenvalue correlations form instead a Pfaffian point process. Thus, considering for definiteness the real eigenvalues, one now has

A concern of the present paper is to compute the explicit form of the correlation kernel in (1.11) in the case of the real eigenvalues of (1.3) for real standard Gaussian matrices, and also for the case of the complex eigenvalues. In this paper, we will see that these correlation kernels possess many similarities with other results for product of random matrices. For example, the kernel for the Pfaffian point process specifying the scaled statistical state about the origin of the real eigenvalues of products of real Ginibre matrices is given in terms of Meijer $G$-functions. In the simplest case of the one point function $\rho_{(1)}^{r}(x)$ the resulting functional form is very succinct.

For singular values of products of complex Ginibre matrices, it is similarly the case that the kernel for the scaled determinantal point process in the neighbourhood of the origin can be expressed in terms of Meijer $G$-functions ; see also the recent review . Moreover, for fixed $N$, knowledge of the real-to-real eigenvalue correlations gives information about the moments of the distribution function for the probability that there are $k$ real eigenvalues. In particular, integration of the spectral density (one-point function) gives the expected number of real eigenvalues.

The rest of this paper is organised as follows. In section 2 we find the joint eigenvalue PDF for a Gaussian product matrix with a given number of real eigenvalues. In section 3 we introduce the generalised partition function and find the skew-orthogonal polynomials; we combine these results with the joint eigenvalues PDF to prove Theorem 1. Section 4 focuses on the real-to-real and the complex-to-complex eigenvalue correlations. In particularly, we study local and global scaling limits for the spectral densities and use the real spectral density to compute the expected number of real eigenvalues. The final section briefly sketches how all these results may be extended to products of rectangular matrices.

Joint probability density function

Our first task is to find the explicit functional form for the eigenvalue PDF of the random matrix product (1.3) in the case that each $X_{i}$ is an independent $N\times N$ standard real Gaussian matrix. With this specification the joint probability measure for $\{P_{m},X_{1},\dots,X_{m}\}$ is equal to

Actually this task, extended to the general bi-orthogonal invariant ensembles, has already been addressed by Ipsen and Kieburg . However the workings therein are not sufficient for all our purposes. In particular proportionality constants are ignored, meaning that it is not possible to proceed to derive the formulas of Theorem 1 for the probabilities $p_{N,k}^{P_{m}}$. These normalisation constants were included in the thesis but the PDF were given in terms of $2\times 2$ matrices, which is impractical for our purpose. Furthermore, the case that the working of — which is a generalisation of the strategies used in and in the cases $m=1$ and $m=2$ respectively — treats the real and complex eigenvalues on an equal footing, whereas we prefer to proceed in the way used in for $m=1$ which distinguishes the real and complex eigenvalues from the outset. Below we give a more practical formulation of the joint eigenvalue PDF.

referred to as the real (or one-point) weight function and let

referred to as the complex (or two-point) weight function.

Consider the product (1.3). Given that there are $k$ real eigenvalues ($k$ of the same parity as the matrix dimension $N$), the joint eigenvalue PDF is

with $Z_{N}$ given by (1.2) and $w_{r},w_{c}$ as above.

The starting point is to use a generalised real Schur decomposition to triangularise the matrices $\{X_{l}\}_{l}$ which appear in the product (1.3). Assuming that the product matrix (1.3) has $k$ real eigenvalues, the decompositions states that for invertible matrices (Gaussian matrices are invertible almost surely) we may write [30, Prop. A.26]

with $Q_{m+1}:=Q_{1}$. Here each $Q_{l}$ is a real orthogonal matrix in $O^{*}(N)/O^{*}(2)^{(N-k)/2}$ with $O^{*}(N)$ defined to be the set of matrices in $O(N)$ with the first entry in each column positive. Each $D_{l}$ is a (block) diagonal matrix with the first $k$ diagonal entries scalars $\{\lambda_{1}^{(l)},\dots,\lambda_{k}^{(l)}\}$ and the next $(N-k)/2$ block entries $2\times 2$ matrices $\{G_{s}^{(l)}\}_{s=k+1}^{(N+k)/2}$, while each $T_{l}$ is a strictly upper triangular matrix consisting of $N(N-1)/2-(N-k)/2$ independent Gaussian random variables.

The generalised Schur decomposition may be verified by applying an ordinary Schur decomposition on the product matrix (1.3) itself and then using (partial) QR decompositions on $\{Q_{l}X_{l}\}_{l=1,\ldots,m-1}$, recursively (see [30, Appendix A] for details). We stress that while it is possible to choose $m-1$ of the matrices $D_{l}$ in (2.6) to be strictly diagonal rather than block diagonal (due to the $m-1$ QR decompositions), we do not do so as it would complicate the derivation of the Jacobian.

For the following, it will be convenient to introduce the product $D:=D_{1}\cdots D_{m}$ which again is a block diagonal matrix. The first $k$ diagonal entries are scalars, $\{\lambda_{t}:=\lambda_{t}^{(1)}\cdots\lambda_{t}^{(l)}\}_{t=1}^{k}$, while the latter $(N-k)/2$ entries are $2\times 2$ matrices, $\{G_{s}:=G_{s}^{(1)}\cdots G^{(l)}_{s}\}_{s=k+1}^{(N+k)/2}$. With this notation, the Jacobian for the above given change of variables reads [30, Prop. A.26]

This notation is the same as used by Edelman [13, Eq. (6)]. More compactly, we may write

where the Vandermonde determinant is defined as in (1.1).

where we can integrate out the dependence on $\{T_{l}\}$ and $\{Q_{l}\}$ according to

The latter is equal to ${{\rm vol}\,O^{*}(N)}/{({\rm vol}\,O^{*}(2))^{(N-k)/2}}$.

Using all the above results, it follows that, for a given $k$, the joint probability measure for the eigenvalues is equal to

We have, at this point, not yet explicitly introduced the constraint that the eigenvalues of each $G_{s}$ are not real and thus are consequently a complex conjugate pair. For this reason, we have a similarity with [30, Prop. 4.26].

In order to explicitly impose our constraint that the product matrix has exactly $k$ real eigenvalues, we suppose an orthogonal similarity transformation has been used to bring each matrix $G_{i}$ into the form

with $b,c>0$. The eigenvalues are then $x\pm iy$ with $y^{2}=bc$, and we know too (see e.g. [16, Proof of Prop. 15.10.1 and Prop. 15.10.2]) that changing variables from the elements of $G_{i}$ to $\{x,y,\delta,\theta\}$, where $\theta$ parametrises the orthogonal similarity transformation and $\delta=b-c$ introduces the Jacobian

The fact that the integrand in (2.4) is invariant under real orthogonal transformations allow us to simplify further. Firstly, we may integrate out $\theta$, which contributes with an extra factor of $\pi$. Secondly, we may replace the matrix $G_{i}$ by the diagonal matrix of its singular values, $\mu_{+}$ and $\mu_{-}$ say. In terms of the variables $x,y,\delta$ it is straightforward to compute that the singular values are given by (2.3). Combining these results completes the proof. ∎

Due to the relatively involved expression for the two-point weight (2.3), it might be beneficial to briefly expand on the simplest cases, $m=1$ and $m=2$, where explicit expressions are known.

For $m=1$, the joint PDF (2.5) must, of course, reduce to the classical result (1.1). Inspection of (2.2) and (2.4) shows that the integration therein are immediate for $m=1$ due to the delta functions. In the real case we then read off that $w_{r}(\lambda)=e^{-\lambda^{2}/2}$. In the complex case, substituting in (2.3) gives

where the second equality first appeared in , albeit out by a factor of 2 as remarked in . Substituting these evaluations in (2.5) indeed reproduces (1.1).

Returning now to the case $m=2$, the Meijer $G$-function (2.2) is a modified Bessel function,

To simplify (2.3) requires simplifying (2.4). For this purpose, and without yet restricting $m$, we introduce $2\times 2$ real matrices $\{M^{(l)}=G^{(l)}\cdots G^{(1)}\}_{l=1,\ldots,m}$ and set $M^{(0)}$ equal to the $2\times 2$ identity matrix. We note that

which allows the integration over $M^{(m)}$ to be carried out in (2.4) using the delta function, showing that

This is the two-by-two matrix version of [30, Eq. (2.20)]. A further change variables $A^{(l)}=M^{(l)}M^{(l)T}$ for each $l=1,\dots,m-1$ shows

where the integration is over positive-definite real symmetric matrices $A^{(l)}$, $l=1,\dots,m-1$. In the case $m=2$ we can also express the integral in terms of modified Bessel functions.

Write for the $2\times 2$ positive definite matrix $A$

Using the notation $h=b_{1}b_{2}-c^{2}$ for the determinant, and expressing this equation as a delta function constraint allows us to write

The working now is elementary. We first integral over $c$, change variables $h\mapsto b_{1}b_{2}h$, $w\mapsto w/b_{1}b_{2}$, and integrate over $w$, then $b_{1}$ and $b_{2}$, using the fact that

The last step is to change variables $s=1/\sqrt{h}$. ∎

Alternative expressions for $I(\mu_{+},\mu_{-})$ are known. One, which involves not the $K_{0}$ Bessel function but rather the $I_{0}$ Bessel function is based on changing variables to the eigenvalues and eigenvectors of $A$, $B$, and using the matrix integration formula for the integral over Haar measure of the $2\times 2$ orthogonal group restricted to matrices with elements in the first entry of each column positive,

Another, which is based on working similar to that used in the proof of Lemma 4, but starting from (2.12) rather than (2.13) tells us that

There is some advantage in the form (2.14), due to its functional dependence on $\mu_{+}^{2}+\mu_{-}^{2}$ and $\mu_{+}\mu_{-}$, which according to (2.3) are given in terms of $\delta,x,y$ by

Recalling the definition of $I(\mu_{+},\mu_{-})$ in (2.13), this tells us that

Substituting in (2.3) and using the integral in (2.10) to integrate over $\delta$ we obtain

In the following, we will see that it is possible to calculate the probability finding exactly $k$ eigenvalues without such explicit knowledge of the two-point weight function (2.3). Here, we make note of them to make contact with the existing literature and as a reference for a comment in section 4.1.

Finally, we note that an important difference compared to the result presented in [30, Prop. 4.26] is the shift from the two-by-two matrix weight function (2.4) to (2.3) which will be essential in the remaining sections.

Generalised partition function, skew-orthogonal polynomials and proof of Theorem 1

Let us denote the joint PDF (2.5) by $\mathcal{Q}(P_{m})$, and define the generalised partition function for $k$ real and $(N-k)/2$ complex conjugate pairs of eigenvalues by

We have that with $u=v=1$ the generalised partition function (3.1) is the probability of finding $k$ real eigenvalues and $(N-k)/2$ complex conjugate pairs of eigenvalues. Functional differentiation of

where the sum is restricted to $k$ of the same parity of $N$ allows the correlation functions to be computed; see e.g. [16, §15.10].

Independent of the specific functional form of $w_{r}$ and $w_{c}$ in (2.5), an observation of Sinclair tells us that due to the product of difference $\Delta$, the method of integration over alternative variables implies that $Z_{k,(N-k)/2}[u,v]$ can be written as a Pfaffian. The details of the necessary working can be found in e.g. [16, Prop. 15.10.3, $N$ even] and [41, §4.3.1 ($N$ even) and §4.3.2 $N$ odd)]. We report the final result only.

Let $\{p_{l-1}(x)\}_{l=1,\dots,N}$ be a set of monic polynomials, with $p_{l-1}(x)$ of degree $l-1$. Let

with $[\xi^{k}]f(\xi)$ defined as in Theorem 1 and $Z_{k,N}$ given by (1.2).

Skew orthogonal polynomials

The matrix $[\zeta\alpha_{j,l}+\beta_{j,l}]$ is antisymmetric. For $\zeta=1$ and $u=v=1$, it is possible to choose the monic polynomials $\{p_{l-1}(x)\}_{l=1,\dots,N}$ so that this anti-symmetric matrix is block diagonal, with the blocks $2\times 2$ anti-symmetric matrices

$j=1,\dots,N/2$, $N$ even, and $j=1,\dots,(N-1)/2$, $N$ odd, with the last diagonal entry 0 in this latter case. In fact, from a theoretical perspective this is also true for general $\zeta$, however our method for these polynomials (given below) is only valid if $\zeta=1$. The use of skew-orthogonal polynomials is standard in random matrix theory; see e.g. [16, Ch. 6]. A Gram–Schmidt procedure shows that the construction of such polynomials is always possible, and that they are unique up to the mapping

This mapping, for $\gamma_{2m}$ an arbitrary constant, leaves the skew-orthogonality property unchanged.

With $\alpha_{j,l}$ and $\beta_{j,l}$ specified by (5) define the skew-product

In the case $m=1$, when the underlying point eigenvalues PDF is given by (1.1), the corresponding skew-orthogonal polynomials were first determined by Forrester and Nagao . They were found to be

The case $m=2$ has been considered by Akemann and collaborators . In fact these authors determined the skew-orthogonal polynomials for a more general model, in which the matrices $X_{1}$ and $X_{2}$ in (1.3) with $m=2$ are a general linear combination of Gaussian symmetric and anti-symmetric matrices, as already noted below Theorem 1. Specialising to the case that $X_{1}$ and $X_{2}$ are both standard Gaussian matrices, we read off the skew-orthogonal polynomials

Comparing the skew-orthogonal polynomials (3.8) and (3.10) as well as the normalisations (3.9) and (3.11) a simple pattern seems apparent: the coefficient in the odd skew-orthogonal polynomials as well as the normalisations constants are raised to powers of $m$, $m=1$ and $m=2$, respectively. This pattern indeed persists in the general case.

form a skew-orthogonal set with normalisation

In the method used to find the skew-orthogonal polynomials was to first establish that with $m=1$

which in turn made essential use of knowledge of the explicit functional form of $w_{r}(x)$ and $w_{c}(x,y)$. Some of the details of the working are given in [16, Proof of Prop. 15.10.4]. Soon after Sommers noted that knowledge of the functional form of the averages of the product of two characteristic polynomials $C_{N}(z)=\prod_{j=1}^{k}(z-\lambda_{j})\prod_{s=k+1}^{(N+k)/2}(z-(x_{s}+iy_{s}))(z-(x_{s}-iy_{s}))$, summed over $k$ contains sufficient information to fully determine the skew-orthogonal polynomials. Subsequently, Akemann, Kieburg and Phillips [6, Eqns. (4.6)–(4.7)] gave the explicit matrix averages formulas

for the skew-orthogonal polynomials, where in the present setting the average over $G$ is over the $m$ standard Gaussian matrices $X_{1},\dots,X_{m}$ of size $2n\times 2n$.

Forrester gave a systematic way to compute averages of the form (3.2) in the cases that $G$ is drawn from an ensemble invariant under real orthogonal transformations. This method was based on the use of zonal polynomials, and was built on ideas contained in . Here we will show that elementary methods suffice to evaluate (3.2).

According to (1.3), $P_{m}$ is the product of $m$ independent standard Gaussian matrices $X_{1},\dots,X_{m}$. Moreover, from the rule for matrix multiplication, and this specification of the $X_{i}$, we see that elements taken from distinct rows $j_{1},\dots,j_{r}$ and columns $k_{1},\dots,k_{r}$, each $j_{\mu}\neq k_{\nu}$ are uncorrelated, so that with $P_{m}=[y_{jk}]_{j,k=1,\dots,N}$

From the definition of a determinant we have

where $\varepsilon(\sigma)$ denotes the parity of the permutation $\sigma$ and $\delta_{l,\sigma(l)}$ denotes the Kronecker delta. Averaging over $X_{1},\dots,X_{m}$ using (3.16) shows that the only non-zero term comes from the identity permutation and furthermore this average is equal to $z^{2n}$. This establishes $p_{2j}(z)$ in (3.12).

Here the final equality follows by noting that $y_{l,l}$ consists of a sum of $(2n)^{m-1}$ terms which are monomials in the elements of the $X_{i}$, and due to (3.16) the only terms that survives this averaging after squaring are the $(2n)^{m-1}$ perfect squares, which contribute unity. Substituting this result into (3.2) establishes $p_{2j+1}(z)$ in (3.12).

It remains to establish (3.13). On this point, we first note that from the meaning of $Z_{k,(N-k)/2}[u,v]|_{u=v=1}$ as the probability that there are exactly $k$ real eigenvalues, it follows that $Z_{N}[u,v]|_{u,v=1}=1$, where $Z_{N}[u,v]$ is specified by (3.2). On the other hand, it follows from (3.5) that for $N,k$ even

Setting $u=v=1$, and using the skew-orthogonal polynomials, the RHS can be evaluated to give

Examination of the above proof shows that invariance of a single matrix entry under the reflection $y_{jk}\mapsto-y_{jk}$ implies

Probability of k real eigenvalues

It has already been remarked below the definition of the generalised partition function (3.1) that with $u=v=1$ this quantity can be interpreted as the probability $p_{N,k}^{P_{m}}$ that for the ensemble of matrices specified by (1.3), with each $X_{i}$ therein an $N\times N$ real standard Gaussian, there are exactly $k$ real eigenvalues. This assumes $k$ and $N$ have the same parity; if not the probability is zero. According to Proposition 5 these probabilities can be written as Pfaffians. Let us suppose the polynomials therein are furthermore even (odd) when there degree is even (odd). We then know, by the symmetry of the integrands, that each $(\zeta\alpha_{j,l}+\beta_{j,l})|_{u=v=1}=0$ unless the parity of $j$ and $l$ is opposite. Furthermore making use of the fact that the $(\zeta\alpha_{j,l}+\beta_{j,l})|_{u=v=1}$ is anti-symmetric in $j,l$ allows the Pfaffian to be written as a determinant of half the size, telling us that for $N,k$ even

We are now well placed to establish (1.9) and (1.10).

We choose the polynomials $\{p_{j}(x)\}$ as the skew-orthogonal polynomials (3.12) so we have

with the explicit value of $u_{j-1}$ being given by (3.13). Thus we have been able to eliminate the dependence on $\beta_{j,k}$, which from the definition (5) involves the weight $w_{c}(x,y)$ — a quantity which from (2.3) is not in general known in terms of explicit special functions. The remaining quantity $\alpha_{2j-1,2l}$ is specified by (5), and the weight therein $w_{r}(x)$ is given as a Meijer $G$-function according to (2.2). In fact this very same quantity, up to a proportionality has appeared in the earlier study [18, Proposition 3] and we read off the evaluation

where we use the definition (1.7) with $a_{j,0}:=0$ since the lowest order odd skew-orthogonal polynomial is a monomial. Substituting (3.22) in (3.21), and substituting the result in turn in (3.19) we obtain after minor manipulation the formula (1.9).

To deduce (1.10) we require the additional evaluation, also contained in [18, Proposition 3], $\mu_{2j-1}=\Gamma(j-{1/2})^{m}$, and similarly substitute in (3.20). ∎

For $m=1$ the probabilities $p_{N,k}^{P_{1}}$ have been known since the late nineties and they are all of the form $r+\sqrt{2}s$ where $r$ and $s$ are rational numbers . Tabulations can be found in [13, Table 5] and [5, Table 2]. Recently, an evaluation of the Meijer $G$-function

as a summation over a linear combination of $\{{}_{2}F_{1}(\mu+a,\mu+b;\mu+c;1-z)\}_{\mu=0}^{n}$ has been given by Kumar , and this was used to show

which allows us to get explicit expressions for the probabilities $p_{N,k}^{P_{2}}$ (i.e. $m=2$). Note in particular that this is of the form $\pi^{2}$ times a rational number, a feature which was conjectured in . Substituting in (1.9) and (1.10) in the case $m=2$ makes the structure of the probabilities explicit for $p_{N,k}^{P_{2}}$. These are all polynomials of degree $\lfloor N/2\rfloor$ in $\pi$ with rational coefficients; probabilities for low values of $N$ are tabulated in Table 1. It is worth noting that similar probabilities for the real spherical and the truncated orthogonal ensembles are also given as polynomials in $\pi$ and $1/\pi$; see and references therein for an extensive summary.

Beyond the cases $m=1$ and $m=2$, evaluation formulas for the Meijer $G$-function in (1.7) are challenging. In addition to the contour integral representation (1.5), we may also write the Meijer $G$-function as an $m$-fold integral on the real line,

which may be checked to agree with (3.23) for $m=2$. Such $m$-fold integral representations give a relation to product of random scalars. However, explicit expressions in terms of elementary functions remain unknown for $m\geq 3$.

With an explicit method for calculating the probability of finding $k$ real eigenvalues, it seems natural to ask for different types of number statistics. A prime example would be the expected number of real eigenvalues. Albeit such expectation values may be calculated using Theorem 1, we will see in section 4.2 that the spectral density for the eigenvalues can be used to obtain a more efficient formula. The interest in real eigenvalue statistics, of course, extends beyond the expected number of real eigenvalues. Another common question is to ask for extreme value statistics, i.e. the probability that there are abnormally many (or few) real eigenvalues. As mentioned in the introduction, the probability that all eigenvalues are real has already be studied in , which led to the remarkable conclusion that this probability tends to unity for $m\to\infty$. It is more challenging to ask for the probability of finding only a few real eigenvalues in the large-$N$ limit, say the probability that an even dimensional product matrix has no real eigenvalues.

A step in this direction was taken in , where using the relation to the Brownian annihilation process $A+A\to\varnothing$, the first two terms of the large $s$ asymptotics of the probability that there are no real eigenvalues in an interval of size $s$ near the origin for $N\to\infty$ real Ginibre ($m=1$) was computed. It was realized Kanzieper et al. that heuristic at least this result implies for large $N$

with $\zeta(x)$ denotes the Riemann zeta function and

and moreover these authors gave a rigorous proof of the leading term. It is not known how to generalize the workings of , which are based on Theorem 1, beyond $m=1$. However, our Theorem 1 at least allows us to establish numerical estimates, e.g. fitting $aN^{1/2}+bN^{0}+cN^{-1/2}$ to $\log p_{N,0}^{P_{2}}$ for $N=50,52,\ldots,120$ suggest that

for $N$ even. We note that $1.474>\zeta(3/2)/\sqrt{2\pi}\approx 1.042$, which is in the agreement with the expectation that $p_{N,0}^{P_{m}}$ decreases with increasing $m$.

Correlation functions

The Pfaffian formulae of Proposition 5 for the generalised partition function, combined with the simplification inherent in the use of skew-orthogonal polynomials, $\{p_{j}(x)\}$, allow the $k$-point correlation to be expressed in the form (1.11) with entries given in terms of $\{p_{j}(x)\}$. While (1.11) refers to the real-to-real eigenvalue correlations, this same structure remains true for the general correlation functions. In fact, the entries of the correlation kernel also have the same structure; see e.g. [41, §4.5 and §4.6].

In this notation, the entries of the correlation kernel (1.11) in the case of the correlation between real eigenvalues only, or the correlation between complex conjugate pairs of eigenvalues are given by

For $N$ odd these expressions require modification; see e.g. , [41, §4.6]. For efficiency of presentation, we will restrict attention to the $N$ even case.

Our main interest in section 4.1 and 4.2 will be spectral densities (one-point correlation functions) and quantities derivable from these. For this reason, we focus on the complex-to-complex and the real-to-real eigenvalue correlations, but real-to-complex correlations can be treated in a similar manner.

We see from (4.3) and (4) that in the case of the correlation between complex eigenvalues, up to factors involving $w_{c}(x,y)$ all the quantities are polynomials, and are related by

Thus it suffices to consider $S(w,z)$, where $w=u+iv$ and $z=x+iy$. For this, (4.3) and (4) tell us that

Upon use of the skew-orthogonal polynomials given by Proposition 9 this simplifies to

We are typically interested in either a global scaling regime (where the eigenvalues are concentrated within a region with compact support) or local scaling regimes (where the eigenvalue interspacing is order unity). For simplicity, let us focus on the one-point function (i.e. the spectral density) which for complex eigenvalues is given by $\rho_{(1)}^{c}(z)=S(z,z)$.

The global scaling regime for the spectral density is known from free probability ,

where $\chi(A)=1$ if $A$ is true, otherwise. This holds because the full spectral density (i.e. including complex as well as real eigenvalues) is dominated by the complex spectrum in the global scaling regime. We note that there also exists a global scaling regime for the real spectrum, albeit sub-dominant. We will return to this limit in section 4.2.

On the local scale, the region near the origin is of greatest interest since it gives rise to new types scalings (i.e. different than the ordinary Ginibre case). The local density near the origin is given by

We recall from section 2 that the weight function $w_{c}(x,y)$ has an explicit and concise expression for $m=1,2$ but not for $m>2$. We note that if $m=1,2$ then the Meijer $G$-function in (4.14) evaluates as

with the latter being a modified Bessel function. Combining this with the weight functions from section 2 reproduces known formulae for the density (the $m=2$ case was given in ).

Real eigenvalues

In this section we focus on the part of the spectrum which is located on the real axis. Similarly to the complex spectrum described above, all correlations may be expressed in terms of the pre-kernel $S(x,y)$. We see from (4.3) and (4) that

which produce the correlation functions by insertion in (1.11). We note that the relations between the pre-kernels (4.16) are more complicated for the real-to-real correlations than for the complex-to-complex correlations where the pre-kernels are related according to (4.11). On the other hand, the weight functions are simpler in the real case (2.2) than in the complex case (2.3).

Using (4.3) and (4) and the skew-orthogonal polynomials (Proposition 9), we write the pre-kernel as

For $m=1$ (i.e. the ordinary Ginibre ensemble), the sum may be rewritten as an incomplete gamma function times an exponential and the integral over $v$ can be performed, which yields

This formulation of the pre-kernel is extremely useful in the study of large-$N$ asymptotics. Unfortunately there are no direct generalisation of this result to $m\geq 2$, which makes asymptotic analysis more challenging. However, it is possible to perform the integral over $v$ in (4.17) for arbitrary $m$. To do so, we rewrite (4.17) as

Now, standard identities for the Meijer $G$-function give

The quantity $\alpha_{j,k}$ is precisely the same quantity appearing in the study , which evaluates to $\alpha_{2j-1,2k}=2^{(j+k-1/2)m}a_{j,k}$ with $a_{j,k}$ given by the Meijer $G$-function (1.7). In the case where the first index of $\alpha_{j,k}$ is even and the second index odd, we use the anti-symmetric property $\alpha_{j,k}=-\alpha_{k,j}$. This gives

where $\lceil\cdot\rceil$ and $\lfloor\cdot\rfloor$ denote the ceiling and floor function, respectively. We recall that the formulae above assume that $N$ is even (for odd $N$ the expression (4.22) is altered by the addition of unity). As already mentioned, an evaluation of $a_{j,l}$ in terms of arithmetic constants is only known for $m=1,2$; consequently the same holds for (4.22). The $m=1$ case is known since the mid nineties , while the $m=2$ case is evaluated using (3.23); the results for small $N$ are tabulated in Table 2. As anticipated, Table 2 reveals that the expected value of real eigenvalues are consistingly larger for $m=2$ than for $m=1$. For $m>2$ a computation of the expectation value (4.22) requires numerical evaluation of the Meijer $G$-functions. The expected number of real eigenvalues can, of course, also be obtained using the probabilities given by Theorem 1. In fact, for $m=2$ and small $N$ the expected number of real eigenvalues follows immediately from Table 1, e.g. for $N=4$ we see that

Let us return to the pre-kernel (4.17) and consider large-$N$ asymptotics for the real spectral density. Similarly to section 4.1 we focus on the local density near the origin and the global density. Using (4.15), it is immediately seen that the local scaling regime near the origin gives (1.13) announced in Theorem 2. Compared to the same result for the complex density (4.14), the real density has the advantage that the weight function $w_{r}(x)$ has a known expression as a Meijer $G$-function (2.2) for all $m$ while $w_{c}(x,y)$ does not. We note that for $m=1$ the Meijer $G$-functions in (1.13) are all simple exponentials; this allows integration over $v$ and confirms that the local spectral density is constant for $m=1$. Moreover, for $m=1$ the corresponding $k$-point correlation takes on the explicit form

as obtained in . We remark that it has been argued by Beenakker and co-workers that the statistical state implied by (4.23) is realised by the level crossings of so-called Majorana zero modes for a disordered semiconducting wire at a Josephson junction, in a weak magnetic field. And this same correlation kernel appears in the seemingly unrelated problem of the annihilation process $A+A\to\emptyset$ in the limit $t\to\infty$ .

A study of the global scaling regime for the real spectrum is more challenging. Unlike the complex spectral density (section 4.1), we have no help from free probability. A qualified guess for this spectral density might be obtained by looking at the $m$-th power of a real Ginibre matrix rather than at the product of $m$ independent matrices. It is immediate that the $m$-th power and the $m$-th product share the same complex macroscopic spectral density, thus assuming that this extends to the real spectrum we expect that

where $\chi(A)$ is defined as in (4.13). For $m=1$ the density (4.24) is well-known ; a verification follows from (4.18) using known asymptotics for the incomplete gamma functions. Moreover, we see that the real spectrum (4.24) develops a non-integrable singularity at the origin when $m$ tends to infinity similarly to (4.13) as we would expect. For $m\geq 2$ we have no rigorous derivation of (4.24) but the form (4.24) is supported by (i) a heuristic saddle point analysis and (ii) numerical data.

Let us first look at the saddle point analysis, which takes (4.17) as the starting point. The first step is to introduce an approximation for the sum in (4.17). We know from [3, Appendix C] that

for $\lvert x\rvert<1$ while exponentially suppressed in $N$ for $\lvert x\rvert>1$. An approximation for the weight function is known from the literature on special functions , and we have

We insert these approximations into (4.17) and want to evaluate the integral over $v$ using a saddle point approximation. Note that there are two maxima of the integrand symmetrically distributed around $v=x$ (the integrand is equal to zero at $v=x$). These two maxima tend to $x$ from left or right, respectively, as $N$ tends to infinity. Thus, we will use an ansatz $v_{*}=x\pm f(x)$ for our saddle points where $f(x)$ is sub-dominant in $N$. With this ansatz and expanding to lowest order, the saddle points are found to be

Evaluation at either of these saddle points yields the conjectured form (4.24) up to a normalisation.

Finally, let us compare the density (4.24) with a simulation of the random matrix product. Figure 1 shows the visual similarity between the density (4.24) for $m=2$ and numerical data stemming from a simulation of $1\,000$ matrix products with $N=1024$. It should be noted that convergence is expected to be exponentially fast in the bulk but considerably slower near the edges. Similar numerical tests have been performed for $m=3,4,5$ and it has been verified that the difference between the analytic formula (4.24) and the numerical data decreases with increasing $N$. Furthermore, we expect that the real global density (4.24) is universal in the sense that the Gaussian entries may be replaced by other independent variables under suitable assumptions on their moments. This type of universality is known to hold for the complex spectra and the expectation that such results extend the real spectra is strengthend by numerical comparison generated from random sign ($\pm 1$) matrices. Although it seems a very natural problem, this type of universality for the real global spectrum has received little attention in the literature; this is true even for the classical Ginibre ensemble ($m=1$).

Rectangular matrices

A generalisation to the case of rectangular matrices is also available and we briefly treat it here. The main idea when dealing with a product of random matrices is to reformulate problem as a product of square random matrices with the same eigenvalue properties; this is possible due to a general reduction procedure presented in (see also [30, Prop. 2.4]). After this reformulation, the approach is similar to the previous sections because Proposition 5 as well as the formulae (4.3) and (4) are completely general. Due to this similarity we will only sketch the main ideas here.

where each matrix $X_{i}$ has dimensions $(N+\nu_{i-1})\times(N+\nu_{i})$ with $\{\nu_{i}\}_{i=0,\ldots,m}$ denoting non-negative integers such that $\nu_{0}=\nu_{m}=0$. Here the constraint is introduced to ensure that the product matrix is square and has $N$ non-trivial eigenvalues. We note that if $\nu_{0}=\nu_{m}>0$ but $\nu_{j}=0$ for some $0<j<m$ (i.e. the smallest matrix dimension is still $N$) then there will be $\nu_{0}$ eigenvalues which are trivially equal to zero (and therefore real) but the joint PDF otherwise remains the same except for an obvious change in normalisation. Consequently, all formulae given below may effortlessly be extended to the $\nu_{0}>0$ case if desired.

The generalisation of the probabilities (1.6) with (1.7) for a purely real spectrum have already appeared in the thesis [30, Prop. 4.29]. They are given by

These formulae allow us to make some straightforward generalisations of the exact expressions presented by Kumar in the $m=2$ case. Following , we have

The next step would be to rewrite gamma functions with a non-integer argument using Gauss’ duplication formula. The right-hand side of (5.4) evaluates as $r\pi^{2}$ for even $\nu$ and $r\pi$ for odd $\nu$ where $r$ denotes some rational constant (depending on both $N$ and $\nu$). This difference in the power of $\pi$ for even and odd $\nu$ has a remarkable consequence: for even $\nu$ the probabilities (5.2) are given as a rational number times $\pi^{\lfloor N/2\rfloor}$ but for odd $\nu$ these constants are simple rational constants (i.e. there is no powers of $\pi$). The probabilities of a purely real spectrum are tabulated in Table 3 for small values of $N$ and $\nu$.

As we have seen in previous sections, to extend the probabilities for a purely real spectrum (5.2) to the probabilities $P_{N,k}^{P_{m}^{\nu}}$ we need a formula for the joint PDF of the eigenvalues and a formula for the skew-orthogonal polynomials, i.e. generalisations of Theorem 3 and Proposition 6. Given such generalisations the rest of the results presented in previous sections may be extended as well due to the generality of Proposition 5.

Given a Gaussian product matrix (5.1) of dimension $N$ with $k$ real eigenvalues, $\{\lambda_{l}\}_{l=1}^{k}$, and $(N-k)/2$ complex conjugate pairs of a eigenvalues, $\{x_{j}\pm iy_{j}\}_{j=1}^{(N-k)/2}$, the joint PDF for these eigenvalues is given by

The proof follows the same lines as the proof of Theorem 3. We use generalised real Schur decomposition to get an expression for the joint PDF in terms of real eigenvalues and $2\times 2$ matrices, see [30, Prop. 4.26]. Finally, changing variables in this expression from the general $2\times 2$ matrix $G$ to a matrix (2.9) using an orthogonal similarity transformation and introducing the singular values, $\mu_{\pm}$, completes the proof. ∎

For the skew-product (3.7) defined in accordance with the joint PDF given by Proposition 8, the polynomials

form a skew-orthogonal set with normalisation

For a product square matrices, we found the skew-orthogonal polynomials by exploiting that elements taken of different rows and columns are uncorrelated. This property is still true for rectangular matrices, thus skew-orthogonal polynomials (5.10) are obtained following the exact same steps. Likewise for the normalisation (5.11) where we evaluate the generalised partition function (3.1) at $u=v=1$ and use (5.9). ∎

With these two propositions established, it is straightforward to extend the rest of our results from square to rectangular matrices. In particularly, we have that the probability of finding exactly $2k$ eigenvalues are real is given by

for $N$ even, while the probability of finding $2k+1$ real eigenvalues is

Moreover, the local densities at the origin is given by

for the real eigenvalues. This generalises (4.14) and (1.13), respectively. The generalised formulae (5.15) and (5.16) follows from the derivations in Section 4.1 and 4.2 now using the weights and polynomials from Proposition 8 and 9. The global densities remains unaltered as long as $\{\nu\}$ are kept fixed in the large-$N$ limit.

Acknowledgements

We would like to thank Mario Kieburg and Oleg Zaboronski comments on this manuscript. Remark 7 on page 7 was given to us by Mario Kieburg. The work of PJF was supported by the Australian Research Council grant DP140102613, and that of JRI by the ARC Centre of Excellence for Mathematical and Statistical Frontiers.

References