<figure><img src="image_1.png" /></figure>

HAL
open science

# On the number of real eigenvalues of a product of truncated orthogonal random matrices

Alex Little, Francesco Mezzadri, Nick Simm

<figure><img src="image_2.png" /></figure>

## To cite this version:

Alex Little, Francesco Mezzadri, Nick Simm. On the number of real eigenvalues of a product of truncated orthogonal random matrices. Electronic Journal of Probability, 2022, 27 (none), 10.1214/21-EJP732. hal-04872039

HAL Id: hal-04872039

https://hal.science/hal-04872039v1

Submitted on 7 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
<figure><img src="image_4.png" /></figure>

Electron. J. Probab. 27 (2022), article no. 5, 1–32.
ISSN: 1083-6489 https://doi.org/10.1214/21-EJP732

# On the number of real eigenvalues of a product of truncated orthogonal random matrices*

Alex Little†

Francesco Mezzadri†

Nick Simm‡

## Abstract

Let $O$ be chosen uniformly at random from the group of $(N + L) \times (N + L)$ orthogonal matrices. Denote by $\tilde{O}$ the upper-left $N \times N$ corner of $O$, which we refer to as a truncation of $O$. In this paper we prove two conjectures of Forrester, Ipsen and Kumar (2020) on the number of real eigenvalues $N_R^{(m)}$ of the product matrix $\tilde{O}_1 \dots \tilde{O}_m$, where the matrices $\{\tilde{O}_j\}_{j=1}^m$ are independent copies of $\tilde{O}$. When $L$ grows in proportion to $N$, we prove that

$$E(N_R^{(m)}) = \sqrt{\frac{2mL}{\pi}} \operatorname{arctanh} \left( \sqrt{\frac{N}{N+L}} \right) + O(1), \quad N \to \infty.$$

We also prove the conjectured form of the limiting real eigenvalue distribution of the product matrix. Finally, we consider the opposite regime where $L$ is fixed with respect to $N$, known as the regime of weak non-orthogonality. In this case each matrix in the product is very close to an orthogonal matrix. We show that $E(N_R^{(m)}) \sim c_{L,m} \log(N)$ as $N \to \infty$ and compute the constant $c_{L,m}$ explicitly. These results generalise the known results in the one matrix case due to Khoruzhenko, Sommers and Życzkowski (2010).

**Keywords:** products of random matrices; real eigenvalues; truncated orthogonal matrices.

**MSC2020 subject classifications:** 15B52, 60B20.

Submitted to EJP on February 25, 2021, final version accepted on December 9, 2021.

<footnote>*A. L. would like to gratefully acknowledge the support of the UK Engineering and Physical Sciences Research Council (EPSRC) DTP (grant number EP/N509619/1). F. M. is grateful for support from a University Research Fellowship of the University of Bristol. N. S. gratefully acknowledges support of the Royal Society University Research Fellowship ‘Random matrix theory and log-correlated Gaussian fields’, reference URF\R1\180707.</footnote>

<footnote>†School of Mathematics, University of Bristol, BS8 1UG, United Kingdom.</footnote>

<footnote>E-mail: al17344@bristol.ac.uk, F.Mezzadri@bristol.ac.uk</footnote>

<footnote>‡Department of Mathematics, University of Sussex, BN1 9RH, United Kingdom.</footnote>

<footnote>E-mail: n.j.simm@sussex.ac.uk</footnote>
# 1 Introduction and main results

Research on products of random matrices started in 1960 with the work of Furstenberg and Kesten [22]. Early investigations concerned the properties and computation of the Lyapunov exponents of products of random matrices when the dimension $N$ is fixed and the number of factors $m$ tends to infinity, see e.g. [8, 12] and references therein. More recently there were several new developments related to the opposite regime, namely when the number of factors $m$ is finite and the matrix dimension $N$ tends to infinity, or is also finite [5]. Early progress in this direction came with the work of Burda et al. [10], who computed the spectral density for a product consisting of independent non-Hermitian matrices whose elements are i.i.d. standard complex normal random variables, the *complex Ginibre ensemble*, in the limit $N \to \infty$ and for finite $m$. Subsequently, Akemann and Burda [2] discovered that the determinantal structure of the eigenvalue point process known to hold for a single complex Ginibre matrix continues to hold for products of such matrices. However, the correlation kernel that arises for products is more complicated and has to be expressed in terms of Meijer G-functions. The determinantal machinery allowed them to analyse the point process of eigenvalues in various microscopic regimes for fixed $m$. A similar theory has been developed to analyse products of truncated unitary random matrices [4, 1, 35]. Recently these developments have been extended to study the question of double scaling limits as both $m$ and $N$ tend to infinity simultaneously [3, 33, 34].

The focus of this paper will be on products of real random matrices. One of the most well-studied examples consists of matrices whose elements are real i.i.d. standard normal random variables, known as the *real Ginibre ensemble*. A distinguishing feature of the real case is that there is a non-zero probability of finding purely real eigenvalues and therefore the point process of eigenvalues consists of two correlated components of purely real and complex conjugated points. This makes the study of real random matrices considerably more challenging than their complex counterparts. The real Ginibre ensemble was introduced in 1965 by Ginibre [25]. After the original paper there was little progress until the work of Lehmann and Sommers [32] who obtained the joint probability density function of all real and complex eigenvalues, later discovered independently by Edelman [14]. These results paved the way for the calculation of correlation functions: it was then discovered that the eigenvalues of the real Ginibre ensemble form a *Pfaffian point process* [21, 40, 7]. Statistics of the real eigenvalues in particular have been shown to arise in many diverse areas, notably in connection to annihilating Brownian motions [43, 42, 37, 16], random dynamical systems [23] and recently to an inverse scattering solution of the Zakharov-Shabat system [6]. Real eigenvalues of products have been applied to physical problems [31, 26] and recently appeared in connection to random Markov matrices [27]. The analogue of the results discussed above for products of complex matrices also hold here: products of independent real Ginibre matrices are again described by a Pfaffian point process [28, 18].

Much less is known about truncated orthogonal random matrices and their products. Their study began in the single matrix case with the work of Khoruzhenko et al. [29] who showed that the eigenvalues form a Pfaffian point process with an explicit correlation kernel. Unlike the real Ginibre ensemble, there is a qualitatively distinct asymptotic regime known as *weak non-orthogonality*, defined by truncating only a finite number of rows and columns from the original orthogonal matrix. This asymptotic regime has shown up in connection to the zeros of random Kac polynomials [17]. More recently, it appeared in [36, 24] where it is related to the persistence exponent of the 2d diffusion equation with random initial conditions. As with the real Ginibre ensemble, the Pfaffian structure of a single truncated orthogonal random matrix has been extended to products
of such matrices [28]. Using this integrable structure, the probability that all eigenvalues of a product of truncated orthogonal random matrices are real was calculated in [20]. Simplified expressions for the correlation kernel of the Pfaffian point process were obtained in the recent work [19].

We now discuss the main results of the paper. Let $O$ denote a random orthogonal matrix of size $(N+L) \times (N+L)$ sampled uniformly with respect to the Haar measure on the orthogonal group and let $\tilde{O}$ denote the upper-left $N \times N$ truncation of $O$. We are interested in the product matrix $X^{(m)} := \tilde{O}_1 \dots \tilde{O}_m$ where $\{\tilde{O}_j\}_{j=1}^m$ are independent copies of $\tilde{O}$. In our first result we shall make the assumption that $L$ grows in proportion to $N$. Let $L := L_N \to \infty$ be a sequence of positive integers such that

$$ \gamma := \frac{L_N}{N} \to \tilde{\gamma}, \quad N \to \infty, \qquad (1.1) $$

where $\tilde{\gamma} > 0$. Throughout the paper we will suppress the $N$-dependence of $L$ and $\gamma$ for notational convenience. We denote by $N_R^{(m)}$ the total number of real eigenvalues of $X^{(m)}$.

**Theorem 1.1.** Let $L := L_N$ be such that hypothesis (1.1) holds and define $\alpha := \frac{1}{1+\gamma}$. Then for any fixed $m \in \mathbb{N}$, we have

$$ \mathbb{E}(N_R^{(m)}) = \sqrt{\frac{2m\gamma N}{\pi}} \operatorname{arctanh}(\sqrt{\alpha}) + O(1), \quad N \to \infty. \qquad (1.2) $$

Next we consider the individual real eigenvalues $\lambda_1, \dots, \lambda_n$ of $X^{(m)}$, where $n = N_R^{(m)}$. The averaged global density of real eigenvalues is defined by

$$ \rho_N(x) := \mathbb{E} \left[ \sum_{j=1}^{n} \delta(x - \lambda_j) \right], \qquad (1.3) $$

and we denote its appropriately normalised version $\tilde{\rho}_N(x) = \frac{1}{\mathbb{E}(N_R^{(m)})} \rho_N(x)$.

**Theorem 1.2.** Let $L := L_N$ be such that (1.1) holds and define $\tilde{\alpha} = \frac{1}{1+\tilde{\gamma}}$. Then for any bounded continuous test function $h$ and fixed $m \in \mathbb{N}$, we have

$$ \lim_{N \to \infty} \int_{-1}^{1} h(x) \tilde{\rho}_N(x) dx = \int_{-1}^{1} h(x) \rho(x) dx, \qquad (1.4) $$

where

$$ \rho(x) = \frac{1}{2m \operatorname{arctanh}(\sqrt{\tilde{\alpha}})} \frac{1}{|x|^{1-\frac{1}{m}} (1 - |x|^{\frac{2}{m}})} \mathbb{1}_{x \in (-\tilde{\alpha}^{\frac{m}{2}}, \tilde{\alpha}^{\frac{m}{2}})}. \qquad (1.5) $$

Furthermore, for any fixed $x \in (-1, 1) \setminus \{-\tilde{\alpha}^{\frac{m}{2}}, 0, \tilde{\alpha}^{\frac{m}{2}}\}$, we have the pointwise limit $\lim_{N \to \infty} \tilde{\rho}_N(x) = \rho(x)$.

Theorems 1.1 and 1.2 prove two conjectures of Forrester, Ipsen and Kumar [19, Conjectures 4.1 and 4.2] and generalise the $m = 1$ asymptotic results of [29]. Theorem 1.1 goes further by providing a precise error bound on the remainder, see Remark A.2 for further discussion on this point. Results of the type (1.2), showing that on average there are of order $\sqrt{N}$ real eigenvalues go back to the work of Edelman et al. [15] in the case of a single real Ginibre matrix, proving an earlier conjecture in [32]. Owing to the mentioned developments linking products of real random matrices to Pfaffian point processes in [28, 18], this was extended to products of independent real Ginibre matrices in [39]. The ‘$\sqrt{N}$-law’ is expected to hold for a broad class of real random matrices [13]. On this point, for $m = 1$ Tao and Vu [41] used moment matching methods to extend the $\sqrt{N}$ estimate in [15] to a class of i.i.d. real random matrices. The fluctuations of $N_R^{(1)}$
about the $\sqrt{N}$-law have been investigated and are known to be Gaussian for the real Ginibre ensemble [38, 30].

The limiting density (1.5) has a simple description in terms of its single matrix $m = 1$ version. If $A$ is the random variable whose density is given by the $m = 1$ case of (1.5), then the symmetrised power $A^m B$ gives the density (1.5) for any $m \ge 1$, where $B$ is an independent Bernoulli random variable on $\{-1, 1\}$. This type of result is familiar in the study of free probability [11] which is very effective at computing the complex spectrum of a product matrix like $X^{(m)}$ in terms of the spectra of the individual factors. However, methods of free probability do not seem to be applicable to the problems solved in this paper regarding real eigenvalues. The main technical challenges in the proofs of Theorems 1.1 and 1.2 arise in the evaluation of integrals (2.6) and (2.7) in Corollary 2.2, whose asymptotics are not uniform due to singularities of the integrands near the origin. This problem is common when dealing with products of random matrices and is discussed in more detail in Section 2.1. It is likely that the technique we develop is useful for other problems involving products of random matrices. For example, our method applies equally well to real Ginibre matrices and this is sketched out in Appendix A.

Finally we consider the regime of weak non-orthogonality, defined by fixing $L$ and letting $N \to \infty$. In this case one is only truncating a finite number of rows and columns, and each matrix in the product is very close to an orthogonal matrix. The $\sqrt{N}$-law discussed above does not hold for these matrices. We find

**Theorem 1.3.** Let $L$ and $m$ be fixed positive integers. Then as $N \to \infty$ we have

$$E(N_{\mathbb{R}}^{(m)}) \sim \frac{1}{B\left(\frac{mL}{2}, \frac{1}{2}\right)} \log(N), \quad (1.6)$$

where $B(a, b) = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}$ is the beta function.

The case $m = 1$ of Theorem 1.3 was obtained in the work of Khoruzhenko, Sommers and Życzkowski [29]. We are not aware of the general $m$ case given by Theorem 1.3 appearing in the literature, conjecturally or otherwise.

This paper is structured as follows. In Section 2 we begin by recalling the results of [19] regarding the Pfaffian structure associated with products of truncated orthogonal random matrices. Then we give an overview of the proof of Theorem 1.1, postponing the full details until Section 3. Then in Section 4 we give the proof of Theorem 1.2 based on the results obtained in Section 3. In Section 5 we prove Theorem 1.3. Section 6 is devoted specifically to the asymptotic analysis of multiple integrals of Laplace type that are needed throughout Sections 2, 3 and 4. Finally, Appendix A includes a discussion about how the approach of the present paper is easily adapted to the case of products of real Ginibre matrices.

## 2 Strategy of the proof and leading order asymptotics

The goal of this section is to discuss the key steps in the proof of Theorem 1.1 and give a quick derivation of the leading asymptotics in (1.2). The starting point of the proof is the fact that that the real eigenvalues of the product matrix $X^{(m)} = \tilde{O}_1 \dots \tilde{O}_m$ form a Pfaffian point process [28]. In the recent work [19] this has been simplified and explicit formulae for the correlation kernel of the Pfaffian point process were identified.

**Theorem 2.1** (Forrester, Ipsen and Kumar [19]). Define the function

$$S(x_1, x_2) := \int_{-1}^{1} dy (x_1 - y) \text{sgn}(x_2 - y) w_L(x_1) w_L(y) f_{N-2,L}(x_1 y), \quad (2.1)$$
where the weight function is

$$w_L(x) = \left( \frac{L}{2B(\frac{L}{2}, \frac{1}{2})} \right)^{\frac{m}{2}} \int_{[-1,1]^m} \delta(x - y_1 \dots y_m) \prod_{i=1}^{m} (1 - y_i^2)^{\frac{L}{2}-1} dy_i, \quad (2.2)$$

and

$$f_{N-2,L}(x) = \sum_{n=0}^{N-2} \binom{L+n}{n}^m x^n. \quad (2.3)$$

Then the real eigenvalues of $X^{(m)}$ form a Pfaffian point process with correlation functions given by

$$\rho_{(k)}(x_1, \dots, x_k) = \mathrm{Pf} [K(x_i, x_j)]_{i,j=1}^k, \quad (2.4)$$

where the correlation kernel is of the form

$$K(x_i, x_j) = \begin{pmatrix} D(x_i, x_j) & S(x_i, x_j) \\ -S(x_j, x_i) & \tilde{I}(x_i, x_j) \end{pmatrix}. \quad (2.5)$$

In (2.5), $D$ and $\tilde{I}$ are certain anti-symmetric functions of $x_i$ and $x_j$ that are not required here; see [19] for further details.

In particular, we have

**Corollary 2.2.** The one-point density of real eigenvalues $\rho_N(x)$ in (1.3) is obtained by taking $x_1 = x_2 = x$ in (2.1), so that $\rho_N(x) = S(x, x)$ and therefore

$$\rho_N(x) = \int_{-1}^{1} dy |x - y| w_L(x) w_L(y) f_{N-2,L}(xy). \quad (2.6)$$

Furthermore, the expected number of real eigenvalues of $X^{(m)}$ is the total integral of $\rho_N(x)$, given by

$$E(N_{\mathrm{R}}^{(m)}) = 2 \int_{0}^{1} dx \int_{-1}^{1} dy |x - y| w_{L}(x) w_{L}(y) f_{N-2,L}(xy), \quad (2.7)$$

where we used that $\rho_N(x)$ is an even function of $x$.

**Remark 2.3.** It is well known that (2.2) has a probabilistic interpretation in terms of products of independent scalar random variables. Namely, up to a constant of proportionality, $w_L(x)$ is equal to the probability density function for a product of $m$ independent Beta $(\frac{L}{2}, \frac{L}{2})$ random variables. Functions of the form (2.2) can also be represented in terms of Meijer G-functions, see [4, 19] for further details.

## 2.1 Preliminary estimates

To prove Theorems 1.1 and 1.2, the idea is to obtain asymptotics of the functions $w_L(x)$ and $f_{N-2,L}(xy)$ separately, and then insert the results into the explicit formulae (2.6) and (2.7). Although this would appear to be a straightforward approach, particular care must be taken with the estimates due to parts of the integral (2.7) where the required asymptotics are not uniform. In particular, when $m > 1$ both $w_L(x)$ and $f_{N-2,L}(xy)$ have singular terms in their asymptotics near the origin, while for $m \ge 1$ the function $f_{N-2,L}(xy)$ undergoes a transition in its asymptotic behaviour near the hyperbola $xy = \alpha^m$.

The purpose of this subsection is simply to detail our asymptotic results for $w_L(x)$ and $f_{N-2,L}(xy)$ that are precise enough to prove Theorems 1.1 and 1.2.
**Remark 2.4.** Throughout the paper, quantities c and C will always denote absolute positive constants that only depend on fixed parameters such as m or κ (see below) and do not depend on N, L or any variables of integration. Furthermore we deem the precise value of these constants as unimportant and caution that their value may change from line to line.

**Proposition 2.5.** Fix a small constant $\kappa > 0$ and a large constant $M > 0$. Then the weight function satisfies the following estimate uniformly on $|x| \in [ML^{-\frac{m}{2}}, 1 - \kappa]$,

$$w_L(x) = d_{L,m} (1 - |x|^{2/m})^{\frac{mL}{2}-1} |x|^{\frac{1}{m}-1} \left( 1 + O \left( \frac{1}{|x|^{\frac{2}{m}L}} \right) \right), \quad L \to \infty, \quad (2.8)$$

where $d_{L,m}$ is defined by

$$d_{L,m} = \frac{1}{2} \sqrt{\frac{L}{\pi m}} \left( \frac{2\pi}{B(\frac{L}{2}, \frac{1}{2})} \right)^{\frac{m}{2}} = \frac{1}{2} \sqrt{\frac{L}{\pi m}} (2\pi L)^{\frac{m}{4}} \left( 1 + O \left( \frac{1}{\sqrt{L}} \right) \right), \quad L \to \infty. \quad (2.9)$$

Furthermore, for $|x| \in [ML^{-\frac{m}{2}}, 1]$ we have the crude estimate

$$w_L(x) \le N^c (1 - |x|^{\frac{2}{m}})^{\frac{mL}{2}}. \quad (2.10)$$

To estimate $f_{N-2,L}(x)$ in (2.3) the idea will be to replace it with the infinite series

$$f_{\infty,L}(x) := \sum_{n=0}^{\infty} \binom{L+n}{n}^m x^n. \quad (2.11)$$

Note that (2.11) converges absolutely inside the unit disc. We will show that the sum (2.3) naturally splits into two contributions, the first coming from $f_{\infty,L}(x)$ restricted to a specific interval, and the second an error term coming from the tail of the sum.

**Proposition 2.6.** Under the hypotheses of Theorem 1.1 the following holds. We have the following estimate as $N \to \infty$, uniformly on $x \in [-1, 1] \setminus ((\alpha - \omega)^m, (\alpha + \omega)^m)$,

$$f_{N-2,L}(x) = f_{\infty,L}(x) \mathbb{1}_{[-(\alpha+\omega)^m, x] \subset (\alpha-\omega)^m} + \frac{x^{N-1}}{x-\alpha^m} e_{N,m} \left( 1 + O \left( \frac{1}{\sqrt{N\omega}} \right) \right), \quad (2.12)$$

where $\omega$ may depend on N, and

$$e_{N,m} := \frac{\gamma^{-m\gamma N - \frac{m}{2}} (1 + \gamma)^{mN(1+\gamma) - \frac{3m}{2}}}{(2\pi N)^{\frac{m}{2}}}. \quad (2.13)$$

In practice we will take $\omega = N^{-\frac{1}{2}}$, such that the big-O term in (2.12) is $O(1)$.

The function $f_{\infty,L}(x)$ can then be analysed separately in the limit $L \to \infty$.

**Proposition 2.7.** Fix a small constant $\kappa > 0$ and a large constant $M > 0$. Then the function $f_{\infty,L}(x)$ satisfies the following estimate as $L \to \infty$, uniformly on $x \in [ML^{-m}, 1 - \kappa]$,

$$f_{\infty,L}(x) = (1 - x^{\frac{1}{m}})^{-mL-1} L^{-\frac{m-1}{2}} \frac{1}{(2\pi)^{\frac{m-1}{2}}} \frac{1}{\sqrt{m}} x^{-\frac{m-1}{2m}} \left( 1 + O \left( \frac{1}{x^{\frac{1}{m}L}} \right) \right). \quad (2.14)$$

Furthermore, on the interval $x \in [0, 1)$ we have the crude estimate

$$f_{\infty,L}(x) \le (1 - x^{\frac{1}{m}})^{-m(L+1)}, \quad (2.15)$$

while for $x \in (-1, 0]$ we have

$$|f_{\infty,L}(x)| \le (1 - |x|^{\frac{1}{m}})^{-m(L+1)} e^{-cL|x|^{\frac{1}{m}}}. \quad (2.16)$$
The proofs of Propositions 2.5, 2.6 and 2.7 are given in full in Section 6. As we shall see, all the quantities involved have convenient integral representations that are amenable to the Laplace or saddle point method for the asymptotic analysis of multiple integrals. Although this is relatively standard, there is a threshold beyond which the Laplace method breaks down, for example if $0 < x < L^{-\frac{m}{2}-\epsilon}$ for any $\epsilon > 0$ in Proposition 2.5.

**Remark 2.8.** The decomposition (2.12) was inspired by analogous results for the truncated exponential function given in [9]. At the level of pointwise asymptotics, the leading terms in (2.8) and (2.14) were also derived in [4]. However, our situation requires more precise estimates and uniform error bounds as given in the above three Propositions. This is due to the specific integral form in (2.1) that is used to construct the correlation kernel for products of truncated orthogonal matrices.

## 2.2 Proof of Theorem 1.1

A straightforward check shows that after formally substituting the asymptotics (2.8) and (2.14) into (2.7), there is a critical point at $y = x$ in the positive quadrant. Then keeping $xy < (\alpha - N^{-\frac{1}{2}})^m$ in accordance with (2.12), our main results are obtained by a final saddle point computation.

To make this approach rigorous we need to remove some regions of the integration in (2.7) where the mentioned asymptotics do not hold or otherwise give negligible contributions. These include a small microscopic layer of width $MN^{-\frac{m}{2}}$ near the origin, a small region near $\pm 1$ of fixed size $\kappa > 0$ and the entire negative y quadrant, see Lemmas 3.1, 3.3 and 3.4 for precise statements. The last of these Lemmas also shows that we can always neglect the second term in (2.12). This results in the following approximation of (2.7):

**Lemma 2.9.** For any $\kappa > 0$ sufficiently small and $M > 0$ a fixed large constant, we have $E(N_R) = I_+ + O(1)$ as $N \to \infty$, where

$$I_+ = 2m^2 \int_{N^{-1/2}M'}^{1-\kappa'} du \int_{N^{-1/2}M'}^{1-\kappa'} dv |v^m - u^m| w_L(u^m) w_L(v^m) \times f_{\infty,L}((uv)^m) (uv)^{m-1} \mathbb{1}_{uv<\alpha-N^{-1/2}}, \quad (2.17)$$

and $M' = M^{\frac{1}{m}}$, $1 - \kappa' = (1 - \kappa)^{\frac{1}{m}}$.

Note that (2.17) is the integral (2.7) with the domain of integration restricted to a certain subset of the positive quadrant, within which we have made the change of variables $x = u^m$ and $y = v^m$. Throughout the paper we will frequently work in the $u$ and $v$ coordinates instead of $x$ and $y$ coordinates as this turns out to greatly simplify the presentation.

The advantage of the approximation (2.17) is that the asymptotics (2.8) and (2.14) are now applicable and have uniform error bounds. As we explained, we may expect the main contributions to come from a small $\delta$-neighbourhood of the line $v = u$ for some fixed $\delta > 0$. This motivates the introduction of the critical region:

$$A := \left\{ (u, v) \in [M'N^{-\frac{1}{2}}, 1]^2 : \{|u - v| < \delta\} \wedge \{uv < \alpha - N^{-\frac{1}{2}}\} \right\}. \quad (2.18)$$

Combined with Lemma 2.9, the following Lemma immediately yields the proof of Theorem 1.1.

**Lemma 2.10.** Let $I_{+,A}$ denote the contribution of the integral (2.17) from the set $A$ in (2.18). Then $I_+ = I_{+,A} + O(e^{-cN})$ and

$$I_{+,A} = \sqrt{\frac{2m\gamma N}{\pi}} \operatorname{arctanh}(\sqrt{\alpha}) + O(1), \quad N \to \infty. \quad (2.19)$$
For simplicity, in this section we will just show the following estimate, contingent on Lemma 2.9,

$$I_{+,A} \sim \sqrt{\frac{2m\tilde{\gamma}N}{\pi}} \arctanh(\sqrt{\tilde{\alpha}}), \quad N \to \infty, \qquad (2.20)$$

where $\tilde{\alpha} := \frac{1}{1+\tilde{\gamma}}$, leaving the full proof of Lemmas 2.9 and 2.10 to Section 3. Inserting (2.8) and (2.14) into (2.17) leads to the approximation

$$I_{+,A} = C_{m,L} \int_{N^{-1/2}M'}^{\sqrt{\alpha-N^{-1/2}}} du \int_{u}^{\min(u+\delta, \frac{\alpha-N-1/2}{u})} dv \times \frac{v^m - u^m}{(uv)^{\frac{m-1}{2}}} T_L(u, v) \left(1 + O\left(\frac{1}{u^2L}\right)\right), \qquad (2.21)$$

where

$$T_L(u, v) = \left[ \frac{(1-u^2)(1-v^2)}{(1-uv)^2} \right]^{\frac{m\gamma N}{2}} (1-u^2)^{-1} (1-v^2)^{-1} (1-uv)^{-1}. \qquad (2.22)$$

In (2.21) we used the symmetry of the set A to integrate only over the region $\{v > u\}$. The pre-factor $C_{m,L}$ is the result of collecting all the individual pre-factors from (2.14) and (2.8), together with a factor $4m^2$ from (2.17) and the symmetry of A. Using the standard asymptotics of the beta function, we have

$$C_{m,L} \sim \tilde{\gamma}N \sqrt{\frac{2m\tilde{\gamma}N}{\pi}}, \quad N \to \infty. \qquad (2.23)$$

The contribution of the term $O(\frac{1}{u^2L})$ in (2.21) will be analysed more precisely in the full proof of Lemma 2.10. For now observe that since $u > M'N^{-1/2}$ its contribution goes to zero when $N \to \infty$ followed by $M \to \infty$; this is sufficient to neglect it at leading order.

Proceeding now with the saddle point asymptotics, we see that the function (2.22) has a critical point at $v = u$. This motivates the change of coordinates $v \to v/\sqrt{N} + u$ in (2.21) and using (2.23) we obtain

$$I_{+,A} \sim \gamma \sqrt{\frac{2m\gamma N}{\pi}} \int_{N^{-1/2}M'}^{\sqrt{\alpha-N^{-1/2}}} du \int_{0}^{\sqrt{N}\min(\delta, \frac{\alpha-N-1/2}{u}-u)} dv \sqrt{N}Q(u, v/\sqrt{N}) \times T_L(u, u+v/\sqrt{N}), \qquad (2.24)$$

where

$$Q(u, v) = \frac{(u + v)^m - u^m}{(u(u + v))^{\frac{m-1}{2}}}. \qquad (2.25)$$

Regarding the integrand of (2.24), we have the easily derived pointwise limits

$$\lim_{N \to \infty} T_L \left( u, u + v/\sqrt{N} \right) = \frac{\exp \left( -\frac{m\tilde{\gamma}v^2}{2(1-u^2)^2} \right)}{(1-u^2)^3}, \qquad (2.26)$$

$$\lim_{N \to \infty} \sqrt{N} Q(u, v/\sqrt{N}) = mv. \qquad (2.27)$$

By (B.1) we have the uniform bound $T_L(u, u + v/\sqrt{N}) \le Ce^{-cv^2}$. Regarding the function $Q(u, v)$ in (2.25), we use the Taylor expansion (B.4) and the fact that $u > M'N^{-1/2}$ to see that $\sqrt{N}Q(u, v/\sqrt{N})$ is uniformly bounded by a polynomial in $v$ with finite degree. Hence
we can apply the dominated convergence theorem and pass the limit inside the integrals (2.24), leading to

$$I_{+,A} \sim \tilde{\gamma} \sqrt{\frac{2m\tilde{\gamma}N}{\pi}} \int_{0}^{\sqrt{\tilde{\alpha}}} du \int_{0}^{\infty} dv \, mv \frac{\exp\left(-\frac{m\tilde{\gamma}v^2}{2(1-u^2)^2}\right)}{(1-u^2)^3} \quad (2.28)$$

$$= \sqrt{\frac{2m\tilde{\gamma}N}{\pi}} \int_{0}^{\sqrt{\tilde{\alpha}}} du \frac{1}{1-u^2} \quad (2.29)$$

$$= \sqrt{\frac{2m\tilde{\gamma}N}{\pi}} \operatorname{arctanh}(\sqrt{\tilde{\alpha}}), \quad (2.30)$$

which is the desired leading order result for the number of real eigenvalues.

### 3 Proof of Lemmas 2.9 and 2.10

We now provide the key estimates that allow us to prove Lemmas 2.9 and 2.10. For a subset $E \subset [0, 1] \times [-1, 1]$, we will denote $I_E$ as the integral (2.7) with the domain of integration restricted to $E$. We begin with the proof of Lemma 2.9. This will be the immediate consequence of the following smaller Lemmas that follow now.

**Lemma 3.1.** Let $M > 0$ be a large fixed constant and define

$$E_1 := \left\{ (x, y) \in [0, 1] \times [-1, 1] : \{x < MN^{-\frac{m}{2}}\} \lor \{|y| < MN^{-\frac{m}{2}}\} \right\}. \quad (3.1)$$

Then $I_{E_1} = O(1)$ as $N \to \infty$.

*Proof.* From the definition (2.2), we have

$$w_L(x) = 2^{m-1} \left( \frac{L}{2B(\frac{L}{2}, \frac{1}{2})} \right)^{\frac{m}{2}} \int_{A_x} \left( 1 - x^2 \prod_{i=1}^{m-1} y_i^{-2} \right)^{\frac{L}{2}-1} \prod_{i=1}^{m-1} \frac{(1-y_i^2)^{\frac{L}{2}-1}}{y_i} d\vec{y}, \quad (3.2)$$

where $A_x = \{\vec{y} \in (0,1)^{m-1} : |x| < \prod_{i=1}^{m-1} |y_i| < 1\}$. Now using the inequality $1 - x^2 \le e^{-x^2}$,

$$w_L(x) \le CL^{\frac{3m}{4}} \int_{[0,\infty)^{m-1}} \exp\left(-\left(\frac{L}{2}-1\right)\left(\frac{x^2}{y_1^2 \cdots y_{m-1}^2} + \sum_{j=1}^{m-1} y_j^2\right)\right) \frac{d\vec{y}}{y_1 \cdots y_{m-1}}, \quad (3.3)$$

and changing variables $y_j \to y_j x^{\frac{1}{m}}$ we see that the last integral is precisely the weight function (A.2) of the real Ginibre ensemble, see the equivalent formula (A.8). We thus have the bound $w_L(x) \le CL^{\frac{3m}{4}} w_{\text{Gin}}((L/2)^{\frac{m}{2}}x)$ where we used that $\frac{L}{2} - 1 > \frac{L}{4}$. Furthermore, we bound the sum $|f_{N-2,L}(xy)|$ in terms of the analogous Ginibre sum (A.3) noting that

$$\sum_{n=0}^{N-2} \binom{L+n}{n}^m |xy|^n \le \sum_{n=0}^{N-2} \frac{(L+N)^{nm} |xy|^n}{(n!)^m} = \sum_{n=0}^{N-2} \frac{|(1+\gamma)^m N^m xy|^n}{(n!)^m}, \quad (3.4)$$

so that by definition of (A.3) we have

$$|f_{N-2,L}(xy)| \le f_{N-2}((1+\gamma)^m N^m |xy|), \quad (3.5)$$

Then applying the change of variables $x \to xL^{-m/2}$ and $y \to yL^{-m/2}$ exactly cancels the large pre-factors coming from the weights and gives the bound,

$$\begin{align}
I_{E_1} &\le C \int_0^{M'} dx \int_0^\infty dy (y+M) w_{\text{Gin}}(x) w_{\text{Gin}}(y) f_{N-2}(c|xy|) \tag{3.6} \\
&\le C \int_0^M w_{\text{Gin}}(x) F(cx) dx \nonumber
\end{align}$$
where we used $|x - y| < |y| + M$ (the sign of $x$ and $y$ is not relevant here) and defined

$$F(x) := \sum_{j=0}^{\infty} \frac{x^j}{(j!)^m} \int_{0}^{\infty} dy w_{\text{Gin}}(y) y^{j+1} + M \sum_{j=0}^{\infty} \frac{x^j}{(j!)^m} \int_{0}^{\infty} dy w_{\text{Gin}}(y) y^j. \quad (3.7)$$

To conclude that $I_{E_1} = O(1)$ it suffices to check that $F(x)$ is bounded on the compact set $[0, M]$. To compute $F(x)$ we use the moment formula

$$\int_{0}^{\infty} dy y^{j+1} w_{\text{Gin}}(y) = \left(2^{\frac{j}{2}} \Gamma(j/2 + 1)\right)^{m}. \quad (3.8)$$

This follows directly from the definition of the Ginibre weight in (A.2). Therefore

$$F(x) = \sum_{j=0}^{\infty} \frac{x^j}{(j!)^m} \left(2^{\frac{j}{2}} \Gamma(j/2 + 1)\right)^m + M \sum_{j=0}^{\infty} \frac{x^j}{(j!)^m} \left(2^{\frac{j-1}{2}} \Gamma(j/2 + 1/2)\right)^m. \quad (3.9)$$

It is straightforward to check that the radius of convergence of these power series is infinite. Consequently $F(x)$ defines an entire function of $x$ (in particular it is continuous and bounded on compact sets) and this implies that (3.6) is bounded. Therefore $I_{E_1} = O(1)$ as required. □

Now in the other regions outside $E_1$ in the positive quadrant, the contribution to (2.7) is exponentially small on any set which is uniformly bounded away from the main diagonal.

**Lemma 3.2.** Fix a large constant $M > 0$. Uniformly on the domain $(x, y) \in (MN^{-\frac{m}{2}}, 1)^2$, we have the bound

$$|x - y| w_L(x) w_L(y) f_{N-2,L}(xy) \le N^c e^{-c' N(u-v)^2} \quad (3.10)$$

where $u = x^{\frac{1}{m}}$, $v = y^{\frac{1}{m}}$. Consequently, if the domain $E \subset (MN^{-\frac{m}{2}}, 1)^2$ is strictly bounded away from the main diagonal $x = y$ by some fixed $\delta > 0$ independent of $N$, then we have the exponential decay $I_E = O(e^{-cN})$ as $N \to \infty$.

**Proof.** This follows immediately from putting together the individual bounds (2.10), $f_{N-2,L}(xy) \le f_{\infty,L}(xy)$, (2.15) and (B.1), together with the fact that $|x - y| \le 1$. □

**Lemma 3.3.** Consider the following small neighbourhood of the hyperbola $xy = \alpha^m$:

$$E_2 := \left\{ (x, y) \in [0, 1]^2 : \left\{ (\alpha - N^{-\frac{1}{2}})^m < xy < (\alpha + N^{-\frac{1}{2}})^m \right\} \right\}. \quad (3.11)$$

Then we have $I_{E_2} = O(1)$ as $N \to \infty$.

**Proof.** By symmetry we consider just the part of $E_2$ where $y > x$. We will make use of the bound $f_{N-2,L}(xy) \le f_{\infty,L}(xy)$ throughout the proof. Let us denote $\alpha_{\pm}^m := (\alpha \pm N^{-\frac{1}{2}})^m$ and split the integration domain as the disjoint union $E_2 \cap \{y > x\} = \bigcup_{j=1}^3 E_{2,j}$ where

$$E_{2,1} = \left\{ \alpha_{-}^{m} < x < t^{m}, \quad \frac{\alpha_{-}^{m}}{x} < y < \frac{\alpha_{+}^{m}}{x} \right\}, \quad (3.12)$$

$$E_{2,2} = \left\{ t^m < x < \alpha_{-}^{\frac{m}{2}}, \quad \frac{\alpha_{-}^{m}}{x} < y < \frac{\alpha_{+}^{m}}{x} \right\}, \quad (3.13)$$

$$E_{2,3} = \left\{ \alpha_{-}^{\frac{m}{2}} < x < \alpha_{+}^{\frac{m}{2}}, \quad x < y < \frac{\alpha_{+}^{m}}{x} \right\}, \quad (3.14)$$
and $t$ is any fixed constant such that $\alpha_- < t < \sqrt{\alpha_-}$. The main contribution will come from a small neighbourhood of the point $x = \alpha_-^{\frac{m}{2}}$, while the contribution from $E_{2,1}$ is exponentially small due to Lemma 3.2. On $E_{2,2}$ and $E_{2,3}$ the asymptotics (2.8) and (2.14) can be applied and we mimic the steps that led to (2.21). Also shifting $v \to v+u$, applying (B.1) and Lemma B.2, we find

$$I_{E_{2,2}} \le CN^{\frac{3}{2}} \int_t^{\sqrt{\alpha_{-}}} du \int_{\frac{\alpha_{-}}{u}-u}^{\infty} dv \, v e^{-cNv^2}, \quad I_{E_{2,3}} \le CN^{\frac{3}{2}} \int_{\sqrt{\alpha_{-}}}^{\sqrt{\alpha_{+}}} du \int_{0}^{\infty} dv \, v e^{-cNv^2}. \quad (3.15)$$

The integrals over $v$ are now explicit and we get

$$I_{E_{2,2}} \le CN^{\frac{1}{2}} \int_t^{\sqrt{\alpha_{-}}} du \, e^{-cN(\frac{\alpha_{-}}{u}-u)^2} = O(1), \quad N \to \infty, \quad (3.16)$$

where the $O(1)$ bound follows from a standard Laplace approximation near the critical point $u = \sqrt{\alpha_{-}}$. Evaluating the second integral in (3.15) we obtain $I_{E_{2,3}} \le CN^{\frac{1}{2}}(\sqrt{\alpha_{+}} - \sqrt{\alpha_{-}}) = O(1)$ as $N \to \infty$. This completes the proof of the Lemma. $\square$

We now show that the second term in (2.12) does not contribute to the leading order asymptotics. Note that for the purposes of proving Lemma 2.9 we do not need to show this on the union $E' = E_1 \cup E_2$ of the two negligible sets of Lemmas 3.1 and 3.3.

**Lemma 3.4.** Consider the contribution of the second term in (2.12) to the integral (2.7) on the set $D := [0, 1 - \kappa] \times [-1 + \kappa, 1 - \kappa] \setminus E'$, taking $\omega = N^{-\frac{1}{2}}$ in Proposition 2.6. Then

$$J := 2 \int_D dx dy |x - y| w_L(x) w_L(y) e_{N,m} \frac{(xy)^{N-1}}{xy - \alpha^m} = O(1), \quad N \to \infty. \quad (3.17)$$

Furthermore, consider the following thin layer near $\pm 1$ of width $\kappa > 0$,

$$E_3 := \{ (x, y) \in [0, 1] \times [-1, 1] : \{x > 1 - \kappa\} \cup \{|y| > 1 - \kappa\} \} \setminus E', \quad (3.18)$$

and finally the remaining part of the negative-$y$ quadrant

$$E_4 := ([0, 1] \times [-1, 0]) \setminus (E' \cup E_3). \quad (3.19)$$

Then for $\kappa > 0$ sufficiently small we have $I_{E_3} = O(e^{-cN})$ and $I_{E_4} = O(1)$ as $N \to \infty$.

**Proof.** We begin with the estimate (3.17). Write $J = J_+ + J_-$ where $J_+$ (resp. $J_-$) denotes the contribution to $J$ from $y > 0$ (resp. $y < 0$). First consider $J_+$, on which we bound the term $|(xy - \alpha^m)^{-1}| \le C\sqrt{N}$ due to the fact that the interval $(\alpha - N^{-\frac{1}{2}})^m < xy < (\alpha + N^{-\frac{1}{2}})^m$ lies outside the integration domain $D$. Inserting the asymptotics of the weights (2.8) and extending the domain of integration to $[0, 1]^2$, we arrive at

$$J_+ \sim CN^{\frac{m+3}{2}} e_{N,m} \int_{\sqrt{\alpha}-\epsilon}^{\sqrt{\alpha}+\epsilon} du \int_{\sqrt{\alpha}-\epsilon}^{\sqrt{\alpha}+\epsilon} dv |u^m - v^m| \times (1-u^2)^{\frac{m\gamma N}{2}-1} (1-v^2)^{\frac{m\gamma N}{2}-1} (uv)^{m(N-1)} \quad (3.20)$$

To obtain (3.20) we have noted that $\sqrt{\alpha}$ is a critical point of the integrand in both $u$ and $v$ variables, so that contributions from outside the small $\epsilon$-neighbourhood of $\sqrt{\alpha}$ are exponentially suppressed. Applying the standard techniques of the Laplace method, we obtain two factors of $N^{-\frac{1}{2}}$ coming from Taylor expanding near the critical point, and an extra factor of $N^{-\frac{1}{2}}$ from the term $|u^m - v^m|$. Then the leading contribution from evaluating the integrand of (3.20) at both critical points gives an upper bound of order

$$O(N^{\frac{m}{2}} e_{N,m} (1 - \alpha)^{m\gamma N} \alpha^{mN}) = O(1), \quad (3.21)$$
where to deduce the $O(1)$ bound, we inserted the explicit form of $e_{N,m}$ defined in (2.13), keeping in mind that $\alpha := \frac{1}{1+\gamma}$.

The procedure for $J_-$ is exactly the same as for $J_+$ except that since now $xy$ is negative, we can bound $|xy - \alpha^m|^{-1}$ by an absolute constant. For the same reason the term $|x - y|$ in (2.7) is of no use in the bounds and we use the simple fact that $|x - y| \le 2$. This results in an upper bound

$$ J_- \le CN^{\frac{m}{2}+1} e_{N,m} \int_{\sqrt{\alpha}-\epsilon}^{\sqrt{\alpha}+\epsilon} du \int_{\sqrt{\alpha}-\epsilon}^{\sqrt{\alpha}+\epsilon} dv \times (1-u^2)^{\frac{m\gamma N}{2}-1} (1-v^2)^{\frac{m\gamma N}{2}-1} (uv)^{m(N-1)}, \quad (3.22) $$

so that $J_- = O(1)$ as $N \to \infty$, by the same reasoning. This completes the proof of (3.17).

Regarding the set $E_3$ in (3.18), in the above bounds for $J_+$ and $J_-$ we saw that the second term in (2.12) can be neglected outside $E'$ and gives an $O(1)$ contribution. In fact repeating these estimates restricted to the set $E_3$ gives a contribution of order $O(e^{-cN})$ because $E_3$ lies outside the saddle point region discussed below (3.20). For this to be true we must choose $\kappa > 0$ small enough such that $1 - \kappa > \sqrt{\alpha}$ for all $N$ sufficiently large, which is guaranteed by hypothesis (1.1). Thus it remains to show the same estimate for the first term in (2.12). For $y > 0$, the indicator function in (2.12) combined with Lemma 3.2 gives the required exponential decay because $E_3 \cap \{xy < (\alpha - N^{-\frac{1}{2}})^m\}$ is uniformly bounded away from the main diagonal. When $y < 0$ we similarly note that $E_3 \cap \{xy > -(\alpha + N^{-\frac{1}{2}})^m\}$ is uniformly bounded away from the anti-diagonal $x = -y$. Then $I_{E_3 \cap \{y < 0\}} = O(e^{-cN})$ follows as in Lemma 3.2, using (2.16) in place of (2.15).

Finally we show that $E_4$ in (3.19) is negligible. By (3.17) we can ignore the second term in (2.12) and deal only with the first term $f_{\infty,L}(xy)$. When $y < 0$ we change variables with $x = u^m$ and $y = -(v^m)$. Making use of the bound (2.16), we insert the asymptotics of the weights (2.8). The combination of the main terms in (2.8) and (2.16) is bounded using (B.1), namely for any $\delta > 0$ we have

$$ \left( \frac{(1-u^2)(1-v^2)}{(1-uv)^2} \right)^{\frac{mL}{2}} e^{-cLuv} \le e^{-\frac{mL}{2}(u-v)^2-cLuv} \\ = e^{-\delta L(u^2+v^2)-L(\frac{m}{2}-\delta)(u-v)^2-L(c-2\delta)uv} \\ \le e^{-\delta L(u^2+v^2)}, \quad (3.23) $$

where we took $0 < \delta < \min(\frac{m}{2}, \frac{c}{2})$ in the final inequality. This gives

$$ I_{E_4} \le CL^{\frac{m}{2}+1} \int_0^1 du \int_0^1 dv (u^m + v^m) e^{-\delta L(u^2+v^2)} + J = O(1), \quad (3.24) $$

and this completes the proof of Lemma 3.4. $\square$

*Proof of Lemma 2.9.* This follows immediately from combining the three Lemmas 3.1, 3.3, 3.4 and noting the presence of the indicator function $xy < (\alpha - N^{-\frac{1}{2}})^m$ in Proposition 2.6. $\square$

*Proof of Lemma 2.10.* We start from the expression (2.21), which after the shift $v \to v+u$ we can write as

$$ I_{+,A} = C_{m,L} \int_{N^{-1/2}M'}^{\sqrt{\alpha - N^{-\frac{1}{2}}}} du \int_{0}^{\min(\delta, \frac{\alpha - N^{-\frac{1}{2}}}{u} - u)} dv e^{\frac{m\gamma N}{2}\phi(u,v)} \\ \times F(u,v) \left(1 + O\left(\frac{1}{u^2L}\right)\right), \quad (3.25) $$
where

$$ \phi(u, v) = \log \left( \frac{(1-u^2)(1-(v+u)^2)}{(1-u^2(v+u)^2)} \right), \qquad (3.26) $$

and in terms of (2.25),

$$ F(u, v) = Q(u, v) (1 - u^2)^{-1} (1 - (v+u)^2)^{-1} (1 - u(v+u))^{-1}. \qquad (3.27) $$

The pre-factor outside (2.21) satisfies

$$ C_{m,L} = \gamma N \sqrt{\frac{2m\gamma N}{\pi}} \left(1 + O\left(\frac{1}{\sqrt{N}}\right)\right), \quad N \to \infty, \qquad (3.28) $$

where we made use of (2.9). Let us denote $I_{+,A}^{(1)}$ the contribution of (3.25) where we discard the error terms $O(N^{-1/2})$ in (3.28) and $O(1/u^2)$ in (3.25); we will explain at the end of the proof why $I_{+,A} = I_{+,A}^{(1)} + O(1)$ as $N \to \infty$. We begin by Taylor expanding $F(u, v)$ and the action $\phi(u, v)$ near $v=0$. A careful estimation of the remainder, see the bound (B.4), shows that uniformly on $u \in [0, 1-\kappa)$, we have

$$ F(u, v) = \frac{mv}{(1 - u^2)^3} + O(v^2) + \sum_{j=0}^{m-1} O\left(\frac{v^{3+j}}{u^{2+j}}\right), \quad 0 < v < \delta. \qquad (3.29) $$

Note that the singular powers of $u$ in (3.29) should be handled carefully. We will show below that they contribute at most $O(1)$ as $N \to \infty$, due to the fact that the $u$-integration starts at $N^{-1/2}M'$, meanwhile the high powers of $v$ give successively smaller contributions to the Laplace asymptotics. In particular it is important that the summation in (3.29) starts at $v^3/u^2$ and not $v^2/u$, as the latter would give an estimate of order $O(\log(N))$ instead of $O(1)$.

For the action we have $\phi(u, v) = -\frac{v^2}{(1-u^2)^2} + E_3(u, v)$ where $E_3(u, v)$ is a remainder satisfying $E_3(u, v) = O(v^3)$ uniformly on $u \in [0, 1-\kappa)$. Let us denote $\xi = \frac{m\gamma N}{2} E_3(u, v)$ and $F_0 = \frac{mv}{(1-u^2)^3}$ in what follows. Then we make the decomposition

$$ e^{\xi} F = F_0 + (F - F_0)e^{\xi} + (e^{\xi} - 1)F_0 \qquad (3.30) $$

and thus write $I_{+,A}^{(1)}$ as the sum of three separate contributions coming from each term in (3.30): $I_{+,A}^{(1)} = J_1 + J_2 + J_3$.

For the leading term $J_1$, the integral over $v$ is explicit and we obtain

$$ J_1 = \sqrt{\frac{2m\gamma N}{\pi}} \int_{N^{-1/2}M'}^{\sqrt{\alpha-N^{-1/2}}} du \frac{1}{1-u^2} \qquad (3.31) \\ - \sqrt{\frac{2m\gamma N}{\pi}} \int_{N^{-1/2}M'}^{\sqrt{\alpha-N^{-1/2}}} du \frac{e^{-\frac{m\gamma N}{2}(\min(\delta, \frac{\alpha-N^{-1/2}}{u}-u))^2}}{1-u^2}. $$

In the second integral of (3.31), define by $u_{-}^*$ the unique positive solution of the equation

$$ \frac{\alpha - N^{-1/2}}{u} = u + \delta. \qquad (3.32) $$

For $u < u_{-}^*$ the minimum in (3.31) is $\delta$ and this yields an exponentially small contribution.

For $u > u_{-}^*$, the minimum is $\frac{\alpha-N^{-1/2}}{u} - u$, which gives an integral of the form (3.16), so that the second term in (3.31) is $O(1)$. Meanwhile the first integral in (3.31) clearly gives the correct leading term up to errors of order $O(1)$ and so

$$ J_1 = \sqrt{\frac{2m\gamma N}{\pi}} \operatorname{arctanh}(\sqrt{\alpha}) + O(1), \quad N \to \infty. \qquad (3.33) $$
For $J_2$ we use (B.1) and (3.29). The contribution from the error terms in (3.29) can then be upper bounded by integrals of the form

$$ \sum_{j=0}^{m-1} N^{\frac{3}{2}} \int_{N^{-\frac{1}{2}} M'}^{1} du \frac{1}{u^{2+j}} \int_{0}^{\infty} dv \, v^{3+j} e^{-cNv^2} dv = O(1), \quad N \to \infty \quad (3.34) $$

and

$$ N^{\frac{3}{2}} \int_{N^{-\frac{1}{2}} M'}^{1} du \int_{0}^{\infty} dv \, v^2 e^{-cNv^2} dv = O(1), \quad N \to \infty. \quad (3.35) $$

Therefore $J_2 = O(1)$. For $J_3$ we use $|e^\xi - 1| \le |\xi|e^{|\xi|} \le CNv^3e^{cN\delta v^2}$. Hence, choosing $\delta > 0$ small enough, the contribution of this error term gives rise to an integral of the form

$$ N^{\frac{5}{2}} \int_{N^{-\frac{1}{2}} M'}^{1} du \int_{0}^{\infty} dv \, v^4 e^{-cNv^2} dv = O(1), \quad N \to \infty, \quad (3.36) $$

and hence $J_3 = O(1)$. Finally, if we had included the error term $O(\frac{1}{Lu^2})$, we would gain an additional power of $N^{-\frac{1}{2}}$ in all the estimates, due to the fact that $\int_{N^{-\frac{1}{2}} M'}^{1} \frac{1}{Lu^2} du = O(N^{-\frac{1}{2}})$ and similarly in (3.34). The same is true had we included the error term of order $N^{-\frac{1}{2}}$ in (3.28). Therefore these error terms can only contribute at most $O(1)$ to $I_{+,A}$ and this completes the proof of (2.19). The fact that $I_+ = I_{+,A} + O(e^{-cN})$ is an immediate consequence of Lemma 3.2. This completes the proof of Lemma 2.10. $\square$

# 4 Convergence of the eigenvalue density

The goal of this section is to prove Theorem 1.2. Recall that

$$ \rho_N(x) = w_L(x) \int_{-1}^{1} dy |x - y| w_L(y) f_{N-2,L}(xy). \quad (4.1) $$

*Proof of (1.4).* Since $h$ is bounded and the integrand of (2.6) is positive, by Lemmas 2.9 and 2.10 we immediately find that

$$ \int_{-1}^{1} h(x) \tilde{\rho}_N(x) dx = \frac{I_{-,A}(h) + I_{+,A}(h)}{\mathbb{E}(N_{\mathrm{R}}^{(m)})} + O(N^{-\frac{1}{2}}), \quad N \to \infty, \quad (4.2) $$

where the analogue of (2.21) is

$$ I_{\pm,A}(h) = \frac{1}{2} C_{m,L} \int_{N^{-\frac{1}{2}} M'}^{\sqrt{\alpha - N^{-\frac{1}{2}}}} du h(\pm u^m) \int_{u}^{\min(u+\delta, \frac{\alpha - N^{-\frac{1}{2}}}{u})} dv G_L(u, v) E_L^{(1)}(u, v) \\ + \frac{1}{2} C_{m,L} \int_{N^{-\frac{1}{2}} M'}^{\sqrt{\alpha - N^{-\frac{1}{2}}}} du h(\pm u^m) \int_{\max(N^{-\frac{1}{2}} M', u-\delta)}^{u} dv G_L(v, u) E_N^{(2)}(u, v), \quad (4.3) $$

and

$$ G_L(u, v) = \frac{v^m - u^m}{(uv)^{\frac{m-1}{2}}} T_L(u, v). \quad (4.4) $$

Here $E_L^{(1)}(u, v) = 1 + O((u^2L)^{-1})$ and $E_L^{(2)}(u, v) = 1 + O((v^2L)^{-1})$. Finally the limit (1.4) follows from changing coordinates $v \to u + \frac{v}{\sqrt{N}}$ (similarly to (2.24)) and again using dominated convergence with the pointwise limits (2.26) and (2.27). $\square$

To obtain the pointwise limit of $\tilde{\rho}_N(x)$ requires some further estimates, but they follow a similar pattern to the proofs of Section 3. For the remainder of this section we assume that $x$ is fixed and satisfies $0 < x < 1$. As in Lemma 3.1, we start by showing that the small region near $y = 0$ can be removed from the integration in (4.1).
**Lemma 4.1.** For any $x \in (0, 1)$ we have

$$w_L(x) \int_{-MN^{-\frac{m}{2}}}^{MN^{-\frac{m}{2}}} dy |x-y| w_L(y) f_{N-2,L}(xy) = O(e^{-cN}), \quad N \to \infty. \quad (4.5)$$

*Proof.* We start with the case that $M_2 N^{-m} < y < MN^{-\frac{m}{2}}$ for a large fixed $M_2 > 0$. We bound $|x-y| \le 2$ and $|f_{N-2,L}(xy)| \le f_{\infty,L}(|xy|)$. In this sense the sign of $y$ is irrelevant and we assume $y > 0$. Then we mimic the proof of Lemma 3.2 noting that

$$\left| w_L(x) \int_{M_2 N^{-m}}^{M N^{-\frac{m}{2}}} dy w_L(y) f_{N-2,L}(xy) \right| \le N^c \int_0^{M N^{-\frac{1}{2}}} dv w_L(v^m) e^{-N(u-v)^2} \quad (4.6)$$

where $x = u^m$. Since $x > 0$ is fixed, so is $u > 0$, and we have that $u - v$ is uniformly bounded away from zero and (4.5) follows on this range. When $0 < y < M_2 N^{-m}$ we use (3.4) so that, with definitions (A.3) and (A.12), we have

$$f_{N-2,L}(|xy|) \le f_{N-2}(cL^m|xy|) \le f_{\infty}(cL^m|xy|) < C, \quad (4.7)$$

where the last bound follows from the fact that for $x > 0$ fixed, $|xy|L^m$ is bounded and $f_{\infty}$ has infinite radius of convergence. Using estimates (4.7) and (2.10) on the left-hand side of (4.5) completes the proof. $\square$

**Lemma 4.2.** For any $x \in (0, 1)$ with $x \neq (\tilde{\alpha})^{\frac{m}{2}}$, we have the following estimates as $N \to \infty$,

$$w_L(x) \int_{|y|>MN^{-\frac{m}{2}}} dy |x-y| w_L(y) e_{N,m} \frac{(xy)^{N-1}}{xy-\alpha^m} 1_{|y|<1} = O(e^{-cN}), \quad (4.8)$$

$$w_L(x) \int_{y<-MN^{-\frac{m}{2}}} dy |x-y| w_L(y) f_{\infty,L}(xy) 1_{|y|<1} = O(e^{-cN}). \quad (4.9)$$

Furthermore, for any $x > (\tilde{\alpha})^{\frac{m}{2}}$, we have $\rho_N(x) = O(e^{-cN})$. For $0 < x < (\tilde{\alpha})^{\frac{m}{2}}$, with $x = u^m$, we have

$$\rho_N(x) = \frac{m}{x^{1-\frac{1}{m}}} \int_{u-\epsilon}^{u+\epsilon} dv |u^m - v^m| w_L(v^m) w_L(u^m) f_{\infty,L}((uv)^m) (uv)^{m-1} + O(e^{-cN}). \quad (4.10)$$

*Proof.* Since the order of magnitude on the right-hand side of these bounds is at the exponentially small scale, it suffices to make use of the crude estimates (2.10), (2.15) and (B.1) (or (2.16) if $y < 0$). For example, the estimate (4.8) follows by repeating the first steps in the proof of Lemma 3.4, noting that since $x$ is strictly away from the critical point located at $x = (\alpha)^{\frac{m}{2}}$ we only obtain exponentially small contributions, as explained below (3.20). For (4.9) we repeat the estimates in (3.23). Because the integral over $u$ is absent and $u = x^{\frac{1}{m}} > 0$ is fixed we again obtain exponentially small contributions for (4.9). Next, when $x > (\tilde{\alpha})^{\frac{m}{2}}$, the indicator function in (2.12) shows that there exists $\delta > 0$ independent of $N$ such that $y < \tilde{\alpha}^{\frac{m}{2}} - \delta$. Hence $|y - x|$ is uniformly bounded away from zero and $\rho_N(x) = O(e^{-cN})$ follows as in Lemma 3.2. The estimate (4.10) follows similarly. $\square$

**Corollary 4.3.** For any $x$ with $0 < x < (\tilde{\alpha})^{\frac{m}{2}}$, we have

$$\lim_{N \to \infty} \frac{\rho_N(x)}{E(N_R^{(m)})} = \frac{1}{2m \arctan h(\sqrt{\tilde{\alpha}})} \frac{1}{x^{1-\frac{1}{m}}(1-x^{\frac{2}{m}})}. \quad (4.11)$$

*Proof.* This follows from the estimate (4.10) and repeating the steps in Section 2.2, starting at (2.17) and leading to (2.29). $\square$
## 5 Regime of weak non-orthogonality

In this section we will consider the opposite case of (1.1), that is where $L$ is a fixed positive integer, while the size of the original orthogonal matrices from which we truncate tends to infinity with $N$. Our goal is to prove Theorem 1.3 which gives the $\log(N)$ growth of the expected number of real eigenvalues (and the explicit multiplicative constant), as shown for the single matrix $m = 1$ case in [29]. It is reasonable to expect that the real eigenvalue density in this regime also behaves similarly to [29], namely that the real eigenvalues of the product matrix $X^{(m)}$ lie in a neighbourhood of the points $\pm 1$ of size $\frac{1}{N}$ and that after appropriate rescaling by this factor their density converges to a well-defined limit. Proving this rigorously is beyond the scope of the present paper and here we focus exclusively on the calculation of $E(N_{\mathbb{R}}^{(m)})$.

We will employ a slightly different approach based on earlier work of the third author [39] and on an alternative exact formula for the expected number of real eigenvalues in [19]. Specifically, it is shown there that the double integral in (2.7) can be computed as,

$$E(N_{\mathbb{R}}^{(m)}) = 2q_{L,m} \sum_{j=0}^{N-2} \binom{L+j}{L}^m (-1)^j g_j, \quad q_{L,m} = \left( \frac{L\Gamma(\frac{L}{2})\Gamma(\frac{L+1}{2})}{2\sqrt{\pi}} \right)^m, \quad (5.1)$$

where $g_j$ is a particular instance of a Meijer G-function,

$$g_j = G_{2m+1,2m+1}^{m+1,m} \left( \begin{array}{c} \frac{1}{2}, \dots, \frac{1}{2}; \frac{L}{2} + j + 1, \dots, \frac{L}{2} + j + 1, \lceil \frac{j}{2} \rceil + 1 \\ \lceil \frac{j}{2} \rceil, j + 1, \dots, j + 1; \frac{1-L}{2}, \dots, \frac{1-L}{2} \end{array} \middle| 1 \right). \quad (5.2)$$

We refer to [19] and references therein for the precise definition of the Meijer G-function and its basic properties. What is needed here is the following contour integral representation which follows directly from the definition,

$$g_j = \frac{1}{2\pi i} \int_{\gamma} \left( \frac{\Gamma(\frac{1}{2} + s)\Gamma(j + 1 - s)}{\Gamma(\frac{L+1}{2} + s)\Gamma(\frac{L}{2} + 1 + j - s)} \right)^m \frac{1}{\lceil \frac{j}{2} \rceil - s} ds, \quad (5.3)$$

where we may take the contour as the imaginary axis, $\gamma = i\mathbb{R}$. The Meijer G-function coefficients can alternatively be written as a real integral as follows.

**Lemma 5.1.**

$$g_j = \frac{1}{\Gamma(\frac{L}{2})^{2m}} \int_{[0,1]^2} \prod_{\ell=1}^{m} (1-t_\ell)^{\frac{L}{2}-1} (1-r_\ell)^{\frac{L}{2}-1} t_\ell^{\lceil \frac{j}{2} \rceil - \frac{1}{2}} r_\ell^{j-\lceil \frac{j}{2} \rceil} \mathbb{1}_{r_1 \dots r_m > t_1 \dots t_m} dt \, d\vec{r}. \quad (5.4)$$

*Proof.* We have the integral representations

$$\left( \frac{\Gamma(\frac{1}{2} + s)}{\Gamma(\frac{1+L}{2} + s)} \right)^m = \frac{1}{\Gamma(\frac{L}{2})^m} \int_{[0,1]^m} d\vec{t} \prod_{\ell=1}^m t_\ell^{s-\frac{1}{2}} (1-t_\ell)^{\frac{L}{2}-1}, \\ \left( \frac{\Gamma(j+1-s)}{\Gamma(\frac{L}{2}+1+j-s)} \right)^m = \frac{1}{\Gamma(\frac{L}{2})^m} \int_{[0,1]^m} d\vec{r} \prod_{\ell=1}^m r_\ell^{j-s} (1-r_\ell)^{\frac{L}{2}-1}. \quad (5.5)$$

Then inserting (5.5) into (5.3) and interchanging the order of integration, the proof is complete after noting the inverse Mellin transform

$$\frac{1}{2\pi i} \int_{\gamma} \left( \prod_{\ell=1}^{m} \frac{t_{\ell}}{r_{\ell}} \right)^s \frac{ds}{\lceil \frac{j}{2} \rceil - s} = \left( \prod_{\ell=1}^{m} \frac{t_{\ell}}{r_{\ell}} \right)^{\lceil \frac{j}{2} \rceil} \mathbb{1}_{r_1 \dots r_m > t_1 \dots t_m}. \quad (5.6)$$
**Lemma 5.2.** We have the following estimates for even and odd indices:

$$j^{mL} g_{2j} = j^{mL} g_j^{\text{sym}} + \frac{A_{L,m}}{j} + o\left(\frac{1}{j}\right), \quad j \to \infty \qquad (5.7)$$

$$j^{mL} g_{2j+1} = j^{mL} g_j^{\text{sym}} - \frac{A_{L,m}}{j} + o\left(\frac{1}{j}\right), \quad j \to \infty. \qquad (5.8)$$

where

$$g_j^{\text{sym}} = \frac{1}{2} \left( \frac{\Gamma(j+1)}{\Gamma(j+1+\frac{L}{2})} \right)^{2m} \qquad (5.8)$$

and

$$A_{L,m} = \frac{1}{\Gamma\left(\frac{L}{2}\right)^{2m}} \int_{\mathbb{R}_+^{2m}} d\vec{t} d\vec{r} \prod_{\ell=1}^{m} t_{\ell}^{\frac{L}{2}-1} r_{\ell}^{\frac{L}{2}-1} e^{-t_{\ell}-r_{\ell}} \frac{t_1 + \dots + t_m}{2} \mathbb{1}_{r_1+\dots+r_m < t_1+\dots+t_m}. \qquad (5.9)$$

*Proof.* We first prove the expansion of $g_{2j}$ in (5.7). In the integral representation (5.4), we replace the term $\prod_{\ell=1}^{m} t_{\ell}^{-\frac{1}{2}}$ with $1 + (\prod_{\ell=1}^{m} t_{\ell}^{-\frac{1}{2}} - 1)$, considering the second term as a perturbation. The contribution of the main term in (5.4) leads to an integral with permutation invariance between the $r_{\ell}$ and $t_{\ell}$ variables. Hence the indicator function can be dropped after multiplying by $\frac{1}{2}$. The integrals are then explicit and this gives (5.8). In the perturbation term, we change coordinates $t_{\ell} \to 1 - \frac{t_{\ell}}{j}$ and $r_{\ell} \to 1 - \frac{r_{\ell}}{j}$ and apply the dominated convergence theorem to take the limit $j \to \infty$ inside the integrals, using that

$$\lim_{j \to \infty} j \left( \prod_{\ell=1}^{m} \left(1 - \frac{t_{\ell}}{j}\right)^{-\frac{1}{2}} - 1 \right) = \frac{t_1 + \dots + t_m}{2}. \qquad (5.10)$$

When $j$ is odd the procedure is exactly the same, except the term that breaks the symmetry is $1 + (\prod_{\ell=1}^{m} t_{\ell}^{-\frac{1}{2}} - 1)$. Then we take the same limit in (5.10) except the exponent is $+\frac{1}{2}$ and this results in an overall minus sign on the right-hand side of (5.10). This completes the proof of the expansions (5.7). $\square$

**Lemma 5.3.** For any $L, m$ fixed, we have as $N \to \infty$,

$$E(N_{\mathbb{R}}^{(m)}) \sim \log(N) 2q_{L,m} \frac{2^{mL}}{(L!)^m} \left(2A_{L,m} - \frac{mL}{4}\right). \qquad (5.11)$$

*Proof.* We decompose (5.1) as a sum over even and odd indices via $E(N_{\mathbb{R}}) = T_{\text{even}} - T_{\text{odd}}$ and use Lemma 5.2. We have

$$T_{\text{even}} = 2q_{L,m} \sum_{j=0}^{\lfloor \frac{N-2}{2} \rfloor} \binom{L+2j}{L}^m g_{2j} \qquad (5.12)$$

$$= 2q_{L,m} \sum_{j=j_0}^{\lfloor \frac{N-2}{2} \rfloor} \frac{2^{mL}}{(L!)^m} \left(1 + \frac{mL(L+1)}{4j}\right) \left(\frac{1}{2} + \frac{\beta}{j} + \frac{A_{L,m}}{j} + o\left(\frac{1}{j}\right)\right) \qquad (5.13)$$

for some constant $\beta$ whose precise value is unimportant as we shall see below. To derive this we have used the asymptotic expansion of a binomial coefficient and a similar expansion to compute asymptotics of the coefficient $g_j^{(\text{sym})}$ in (5.8). Similarly, we have

$$T_{\text{odd}} = 2q_{L,m} \sum_{j=j_0}^{\lfloor \frac{N-2}{2} \rfloor} \frac{2^{mL}}{(L!)^m} \left(1 + \frac{mL(L+3)}{4j}\right) \left(\frac{1}{2} + \frac{\beta}{j} - \frac{A_{L,m}}{j} + o\left(\frac{1}{j}\right)\right). \qquad (5.14)$$

Then directly computing the difference $T_{\text{even}} - T_{\text{odd}}$ from (5.13) and (5.14) results in several cancellations so that only terms proportional to $\frac{1}{j}$ remain. Collecting the terms results in (5.11). $\square$
*Proof of Theorem 1.3.* In Appendix C we compute the integrals $A_{L,m}$ in (5.9) explicitly. We find that

$$A_{L,m} = \frac{Lm}{8} + \frac{1}{2} \frac{\Gamma(mL)}{\Gamma\left(\frac{mL}{2}\right)^2} 2^{-mL}. \qquad (5.15)$$

Inserting (5.15) into (5.11) and using the duplication formula for the Gamma function completes the proof. □

# 6 Laplace asymptotics

The purpose of this section is to prove the three Propositions 2.5, 2.6 and 2.7.

*Proof of Proposition 2.5.* After the scaling $y_j \to y_j x^{\frac{1}{m}}$ the integral in (2.2) takes the following form:

$$I := \int_{A_x} e^{K\phi(\vec{y})} h(\vec{y}) d\vec{y} \qquad (6.1)$$

where $K = \frac{L}{2} - 1$, $\vec{y} = (y_1, \dots, y_{m-1})$, $d\vec{y} = dy_1 \dots dy_{m-1}$, $h(\vec{y}) = (y_1 \dots y_{m-1})^{-1}$,

$$A_x = \left\{ \vec{y} \in [0, x^{-\frac{1}{m}}]^{m-1} : x^{\frac{1}{m}} \le \prod_{i=1}^{m-1} y_i \le x^{-(1-\frac{1}{m})} \right\}, \qquad (6.2)$$

and the action, which depends implicitly on $x$ throughout the proof, is

$$\phi(\vec{y}) = \log \left( 1 - \frac{x^{\frac{2}{m}}}{y_1^2 \dots y_{m-1}^2} \right) + \sum_{j=1}^{m-1} \log(1 - x^{\frac{2}{m}} y_j^2). \qquad (6.3)$$

Although the asymptotic analysis of (6.1) is standard, special care has to be taken if $x$ depends on $K$ such that $x \to 0$ as $K \to \infty$. One reason is that $\phi(\vec{y})$ vanishes as $x \to 0$ and this occurs at a rate $x^{\frac{2}{m}}$. As we shall see, provided that the combined quantity $\eta := K x^{\frac{2}{m}}$ is large, the Laplace method is still applicable, for which we now give the details.

It is straightforward to check that the only critical point of (6.3) is $\vec{y} = (1, \dots, 1)$. Thus for $\epsilon > 0$ small, we split $A_x$ into the region $E := [1 - \epsilon, 1 + \epsilon]^{m-1}$ and its complement $A_x \setminus E$ and denote the corresponding integrals $I_E$ and $I_{E^c}$ respectively. Note that because $x < 1 - \kappa$, no matter how small we choose $\kappa > 0$, we can make $\epsilon > 0$ small enough such that $E$ is contained strictly inside $A_x$. For $I_E$ we will Taylor expand the action at the point $\vec{y} = \vec{1}$ up to the fourth derivative:

$$\phi(\vec{y}) = \phi(\vec{1}) + \frac{1}{2}(\vec{y} - \vec{1})^T H(\vec{y} - \vec{1}) + E_3(\vec{y} - \vec{1}) + R_4(\vec{y}) \qquad (6.4)$$

where $H$ is the Hessian matrix of $\phi$ at $\vec{1}$, $E_3(\vec{y} - \vec{1})$ is the third order term in the Taylor expansion of $\phi$, and $R_4$ is the remainder term that can be bounded in terms of fourth order partial derivatives of $\phi$.

We will denote in what follows $\xi := K(E_3(\vec{y} - \vec{1}) + R_4(\vec{y}))$. A direct computation of the third and fourth derivatives of $\phi$ shows that they are uniformly bounded on the set $E$ by a factor $x^{\frac{2}{m}}$ times an absolute constant depending only on $m$ and $\kappa$. Therefore we have the following bounds for $\vec{y} \in E$,

$$\xi(\vec{y}) = K E_3(\vec{y}) + O(\eta |\vec{y} - \vec{1}|^4), \qquad (6.5)$$

$$\xi(\vec{y}) = O(\eta |\vec{y} - \vec{1}|^3), \qquad (6.6)$$

$$\xi(\vec{y}) = O(\eta \epsilon |\vec{y} - \vec{1}|^2). \qquad (6.7)$$
The Hessian matrix $H$ can be computed explicitly:

$$H_{ij} = - \frac{4x_{\tilde{m}}^{\tilde{2}}}{(1 - x_{\tilde{m}}^{\tilde{2}})^2} \times \begin{cases} 2, & i = j \\ 1, & i \neq j \end{cases} \qquad (6.8)$$

where $i, j = 1, \dots, m-1$. We have

$$\det(-H) = \frac{4^{m-1} x_{\tilde{m}}^{\frac{2(m-1)}{m}}}{(1 - x_{\tilde{m}}^{\frac{2}{m}})^{2(m-1)}} m. \qquad (6.9)$$

To approximate $I_E$, we note that on $E$ the functions $\xi$ and $h$ are close to zero and 1 respectively, and make use of the trivial identity

$$e^{\xi}h = 1 + (h-1) + (h-1)\xi + \xi + (e^{\xi} - 1 - \xi)h. \qquad (6.10)$$

Then we can decompose $I_E = e^{K\phi(\tilde{1})}(I_1 - I_2 + I_3 + I_4 + I_5 + I_6)$, where after the shift $y_j \to y_j + 1$ we have

$$I_1 = \int_{\mathbb{R}^{m-1}} e^{\frac{1}{2} K \tilde{y}^T H \tilde{y}} d\tilde{y}, \qquad I_2 = \int_{\mathbb{R}^{m-1} \setminus [-\epsilon, \epsilon]^{m-1}} e^{\frac{1}{2} K \tilde{y}^T H \tilde{y}} d\tilde{y}, \qquad (6.11)$$

$$I_3 = \int_{[-\epsilon, \epsilon]^{m-1}} e^{\frac{1}{2} K \tilde{y}^T H \tilde{y}} (\tilde{h}(\tilde{y}) - 1) d\tilde{y}, \qquad I_4 = \int_{[-\epsilon, \epsilon]^{m-1}} e^{\frac{1}{2} K \tilde{y}^T H \tilde{y}} (\tilde{h}(\tilde{y}) - 1) \tilde{\xi}(\tilde{y}) d\tilde{y}, \qquad (6.12)$$

$$I_5 = \int_{[-\epsilon, \epsilon]^{m-1}} e^{\frac{1}{2} K \tilde{y}^T H \tilde{y}} \tilde{\xi}(\tilde{y}) d\tilde{y}, \qquad (6.13)$$

and

$$I_6 = \int_{[-\epsilon, \epsilon]^{m-1}} e^{\frac{1}{2} K \tilde{y}^T H \tilde{y}} \left( e^{\tilde{\xi}(\tilde{y})} - (1 + \tilde{\xi}(\tilde{y})) \right) \tilde{h}(\tilde{y}) d\tilde{y}, \qquad (6.14)$$

where $\tilde{h}(\tilde{y}) = h(\tilde{y} + \tilde{1})$ and $\tilde{\xi}(\tilde{y}) = \xi(\tilde{y} + \tilde{1})$. The integral $I_1$ gives the leading order term:

$$\begin{aligned} e^{K\phi(\tilde{1})} I_1 &= (2\pi)^{\frac{m-1}{2}} \frac{e^{K\phi(\tilde{1})}}{K^{\frac{m-1}{2}} \det(-H)^{\frac{1}{2}}} \qquad (6.15) \\ &= \frac{1}{\sqrt{m}} (2\pi)^{\frac{m-1}{2}} 2^{-m+1} (1 - x_{\tilde{m}}^{\frac{2}{m}})^{Km+m-1} \eta^{-\frac{m-1}{2}}, \end{aligned}$$

where we remind that $\eta := Kx_{\tilde{m}}^{\tilde{2}}$. We must show that the relative error produced by the other terms $\{I_j\}_{j=2}^6$ is no larger than the error term claimed in (2.8), namely that $I_j/I_1 = O(\eta^{-1})$ for each $j = 2, 3, 4, 5, 6$. The smallest eigenvalue of $-H$ will be helpful to prove this and we denote it $\Gamma_1$. From (6.8) we have that there is an absolute constant $c_m > 0$ depending only on $m$ such that

$$\Gamma_1 = c_m \frac{x_{\tilde{m}}^{\frac{2}{m}}}{(1 - x_{\tilde{m}}^{\frac{2}{m}})^2}. \qquad (6.16)$$

Then a standard diagonalisation of $H$ followed by Cauchy-Schwarz gives the bound

$$e^{\frac{1}{2} K \tilde{y}^T H \tilde{y}} \le e^{-\frac{1}{2} K \Gamma_1 |\tilde{y}|^2} \le e^{-c\eta |\tilde{y}|^2}, \qquad (6.17)$$

for some constant $c > 0$. Applying this to $I_2$, we upper bound by extending all integrations to $\mathbb{R}$ except $y_1$. These $m-2$ integrations are explicit Gaussian integrals and yield

$$I_2 \le C\eta^{-\frac{m-2}{2}} \int_{\epsilon}^{\infty} e^{-c\eta y_1^2} dy_1 \le C e^{-c\eta}, \qquad (6.18)$$
which implies $I_2/I_1 = O(e^{-c\eta})$.

The estimation of $I_3, I_4, I_5$ and $I_6$ follows from Taylor expanding $\tilde{h}(\vec{y})$ and applying the estimates on $\xi$ in (6.5), (6.6) and (6.7). The degree of the monomials in $\vec{y}$ occurring in these expansions will give a corresponding power of $\eta^{-\frac{1}{2}}$ after integration. For example, in $I_3$ we Taylor expand $\tilde{h}(\vec{y})$ near 0 to quadratic order, the linear terms produced giving zero in the integration by symmetry. In the error term of quadratic order, we upper bound by replacing the integration with $\mathbb{R}^{m-1}$, use (6.17) and evaluate the resulting Gaussian integral giving $I_3/I_1 = O(\eta^{-1})$.

In $I_4$ and $I_5$ we apply the same approach, using $|\tilde{h}(\vec{y}) - 1||\xi(y)| = O(\eta|\vec{y}|^4)$, which upon integration shows that $I_4/I_1 = O(\eta^{-1})$. Then in $I_5$ the third order terms contained in $\tilde{\xi}(\vec{y})$ vanish in the integration by symmetry. The fourth order remainder term gives a contribution of order $O(\eta|\vec{y}|^4)$. Thus $I_5/I_1 = O(\eta^{-1})$. In $I_6$ we apply the bound $|e^{\tilde{\xi}} - 1 - \tilde{\xi}| \le \frac{|\tilde{\xi}|^2}{2}e^{\tilde{\xi}} = O(\eta^2|\vec{y}|^6e^{\eta\epsilon C|\vec{y}|^2})$, where in the last estimate we applied (6.6) and (6.7). This implies that

$$I_6 \le C \int_{\mathbb{R}^{m-1}} e^{-\eta|\vec{y}|^2(c-C\epsilon)} \eta^2 |\vec{y}|^6 d\vec{y}, \quad (6.19)$$

and choosing $\epsilon > 0$ small enough to ensure $c - C\epsilon > 0$, we also have $I_6/I_1 = O(\eta^{-1})$.

Finally we treat the contribution from the complementary region $E^c = A_x \setminus [1 - \epsilon, 1 + \epsilon]^{m-1}$ to the integral (6.1). Without loss of generality suppose that $|y_1 - 1| > \epsilon$. On this region we will use the bound $1 + x \le e^x$, which implies that

$$\left( \left( \frac{1 - x_{\mathrm{m}}^{2}(y_1^2 \dots y_{m-1}^2)^{-1}}{1 - x_{\mathrm{m}}^{2}} \right) \prod_{j=1}^{m-1} \left( \frac{1 - x_{\mathrm{m}}^{2} y_j^2}{1 - x_{\mathrm{m}}^{2}} \right) \right)^K \le \exp \left( -K \frac{x_{\mathrm{m}}^{2}}{1 - x_{\mathrm{m}}^{2}} \psi(\vec{y}) \right), \quad (6.20)$$

where the action $\psi(\vec{y})$ is independent of $x$ and is given by

$$\psi(\vec{y}) = \sum_{j=1}^{m-1} y_j^2 + \frac{1}{y_1^2 \dots y_{m-1}^2} - m. \quad (6.21)$$

The function $\psi(\vec{y})$ has a unique global minimum at $\vec{y} = \vec{1}$ and $\psi(\vec{1}) = 0$. The integration domain avoids the point $\vec{1}$ and so there exists a constant $c_\epsilon > 0$ such that $\psi(\vec{y}) > c_\epsilon$ for all $\vec{y} \in E^c$. We thus obtain the bound

$$(1 - x_{\mathrm{m}}^{2})^{-Km} I_{E^{c}} \le \int_{\mathbb{R}^{m-2}} \int_{1+\epsilon}^{\infty} \exp(-c\eta \psi(\vec{y})) \frac{d\vec{y}}{\vec{y}} \quad (6.22)$$

$$\le C e^{-c\eta}. \quad (6.23)$$

This implies that $I_{E^c}/I_1 = O(e^{-c\eta})$. Repeating the estimate (6.22) on the full integration interval gives the crude estimate (2.10).

Therefore the biggest error terms came from $I_3, I_4, I_5$ and $I_6$ and had a relative error with $I_1$ of order $O(\eta^{-1})$, as required. $\square$

**Proof of Proposition 2.6.** We recall from the hypothesis of Theorem 1.1 that $L := L_N$ is a sequence of positive integers such that $\gamma := L_N/N \to \tilde{\gamma}$ as $N \to \infty$ with $\tilde{\gamma} > 0$. We also recall that $\alpha := \frac{1}{1+\tilde{\gamma}}$. It will be important in what follows that due to the hypothesis on $\gamma$, we can find $N_0 \in \mathbb{N}$ and $\epsilon > 0$ independent of $N$ such that $\alpha < 1 - \epsilon$ for all $N > N_0$.

We start from the integral representation

$$\binom{L+n}{n} = \frac{1}{2\pi i} \oint_C (1-z)^{-L-1} z^{-n-1} dz, \quad (6.24)$$
where $C$ is any loop enclosing the origin and properly contained in a compact subset of the open unit disc. Introducing $z^{(m)} := z_1 \dots z_m$ and $d\vec{z} = dz_1 \dots dz_m$, this implies that

$$f_{N-2,L}(x) = \frac{1}{(2\pi i)^m} \oint_{C^m} \left( \prod_{j=1}^{m} (1-z_j)^{-L-1} \right) \sum_{n=0}^{N-2} \left( \frac{x}{z^{(m)}} \right)^n \frac{d\vec{z}}{z^{(m)}} \quad (6.25)$$

$$= \frac{1}{(2\pi i)^m} \oint_{C^m} \left( \prod_{j=1}^{m} (1-z_j)^{-L-1} \right) \frac{d\vec{z}}{z^{(m)} - x} \quad (6.26)$$

$$- \frac{x^{N-1}}{(2\pi i)^m} \oint_{C^m} \left( \prod_{j=1}^{m} (1-z_j)^{-L-1} z_j^{-N+1} \right) \frac{d\vec{z}}{z^{(m)} - x}. \quad (6.27)$$

Let us start with the exterior region defined by $(\alpha + \omega)^m < |x| \le 1$ with $\alpha = \frac{1}{1+\gamma}$. We now choose $C$ to be the circle around the origin of radius $\alpha$, noting that our hypothesis on $L_N$ implies that for $N > N_0$ the contour $C$ lies in a compact subset of the unit disc and is bounded away from the unit circle by some fixed $\delta > 0$. This is fundamental and without this assumption the asymptotics of $f_{N-2,L}$ and Theorem 1.1 can take a different form, indeed recall Theorem 1.3.

The integrand in (6.26) is analytic inside the unit circle except at the point $z_m = x/(z_1 \dots z_{m-1})$, which by the choice of contour and condition on $x$ lies outside $C$. Therefore by Cauchy's theorem we have

$$f_{N-2,L}(x) = - \frac{x^{N-1}}{(2\pi i)^m} \oint_{C^m} \left( \prod_{j=1}^{m} (1-z_j)^{-\gamma N-1} z_j^{-N+1} \right) \frac{1}{z^{(m)} - x} d\vec{z} \quad (6.28)$$

$$= - \frac{x^{N-1}}{(2\pi i)^m} \oint_{C^m} \prod_{j=1}^{m} e^{-N\phi(z_j)} h_x(\vec{z}) d\vec{z}, \quad (6.29)$$

where $\phi(z) := \gamma \log(1-z) + \log(z)$ and

$$h_x(\vec{z}) := \left( \prod_{j=1}^{m} \frac{z_j}{1-z_j} \right) \frac{1}{z^{(m)} - x}. \quad (6.30)$$

The integral (6.29) can now be analysed with the Laplace method. We set $z_j = \alpha e^{i\theta_j}$ with $\theta_j \in [-\pi, \pi]$ for each $j = 1, \dots, m$. A simple computation shows that $\theta = 0$ is a critical point of the function $\theta \to \phi(\alpha e^{i\theta})$ on $[-\pi, \pi]$. Now consider the region $E := [-\epsilon, \epsilon]^m$ and denote the contribution to the integral in (6.29) coming from the set $E$ as $I_E$ and similarly denote the contribution from the complementary region $E^c := [-\pi, \pi]^m \setminus E$ as $I_{E^c}$. On $E$ we Taylor expand the action

$$\phi(\alpha e^{i\theta}) = \gamma \log \left( \frac{\gamma}{1+\gamma} \right) + \log \left( \frac{1}{1+\gamma} \right) + \frac{1}{2} \frac{1+\gamma}{\gamma} \theta^2 + O(\theta^3), \quad |\theta| < \epsilon. \quad (6.31)$$

Regarding the function $h_x(\vec{z})$ in (6.30), we have

$$\frac{1}{z^{(m)} - x} = \frac{1}{\alpha^m - x} \left( 1 + \alpha^m \frac{1 - e^{i\theta^{(m)}}}{\alpha^m e^{i\theta^{(m)}} - x} \right), \quad (6.32)$$

where $\theta^{(m)} = \theta_1 + \dots + \theta_m$. As the distance between $\alpha^m$ and $x$ is at least $(\alpha + \omega)^m - \alpha^m \ge m\alpha^{m-1}\omega$, the second term in the brackets in (6.32) is $O(|\theta^{(m)}|\omega^{-1})$. The linear dependence on $\theta_j$'s will contribute an additional factor $N^{-\frac{1}{2}}$ in the Laplace method, so
that this correction will contribute at order $O(\frac{1}{\sqrt{N\omega}})$. Taylor expanding the other terms in (6.30), we see that

$$h_x(\vec{z}) d\vec{z} = i^m \frac{(\gamma(1+\gamma))^{-m}}{\alpha^m - x} (1 + O(\vec{\theta}\omega^{-1})) d\vec{\theta}. \quad (6.33)$$

Using (6.33) and (6.31) we obtain the approximation

$$I_E = \frac{e_{N,m}}{\alpha^m - x} \left( \frac{N(1+\gamma)}{\pi\gamma} \right)^{\frac{m}{2}} \left( \int_{-\epsilon}^{\epsilon} e^{-N\frac{1+\gamma}{\gamma}\theta^2} d\theta \right)^m \left( 1 + O\left( \frac{1}{\sqrt{N\omega}} \right) \right) \quad (6.34)$$

$$= \frac{e_{N,m}}{\alpha^m - x} \left( 1 + O(e^{-cN}) + O\left(\frac{1}{\sqrt{N\omega}}\right) \right). \quad (6.35)$$

On the other hand, the integral $I_{E^c}$ is exponentially small relative to $I_E$. To see this, note that for $z = \alpha e^{i\theta}$ the action satisfies $\mathrm{Re}(\phi(z)) = \gamma \log|1 - \alpha e^{i\theta}| + \log(\alpha)$ and that this function is monotonically increasing in $\theta$ with the minimum obtained at $\theta = 0$ with $\mathrm{Re}(\phi(\alpha)) = \phi(\alpha)$. The monotonicity follows immediately from the formula $|1 - \alpha e^{i\theta}|^2 = (1 - \alpha)^2 + 2\alpha(1 - \cos(\theta))$ and corresponding monotonicity of the log and square root functions. On the set $E^c$ there exists at least one $\ell = 1, \dots, m$ such that $\theta_\ell$ is bounded away from the critical point at 0, i.e. $|\theta_\ell| > \epsilon$, hence from the mentioned properties of $\mathrm{Re}(\phi(z))$ there exists a constant $c_\epsilon > 0$ such that $\mathrm{Re}(\phi(\alpha e^{i\theta_\ell})) - \phi(\alpha) > c_\epsilon$. The latter inequality and (6.32) imply that $I_{E^c}/I_E = O(\frac{e^{-c_\epsilon N}}{\omega})$. This completes the proof of the Lemma in the exterior region $(\alpha + \omega)^m < |x| < 1 - \kappa$.

For the region $-(\alpha + \omega)^m < x < (\alpha - \omega)^m$ we follow the steps leading to (6.27) but we use the big circle

$$C_R = \left\{ -R\alpha + (R+1)\alpha e^{i\theta} : -\pi < \theta \le \pi \right\}. \quad (6.36)$$

In this case we will use the following facts proved in Lemma B.4: Choosing $R$ large enough, in the integration over $z_m$ the pole at $x/(z_1 \dots z_{m-1})$ is strictly inside $C_R$. Furthermore, there exist absolute constants $\delta_0 > 0$ and $c > 0$ such that

$$\inf_{x \in [-(\alpha + \delta_0)^m, (\alpha - \omega)^m]} \inf_{\theta_k \in [-\pi, \pi], k=1, \dots, m} |z^{(m)} - x| > c\omega > 0. \quad (6.37)$$

Hence by Cauchy's integral formula we get

$$f_{N-2,L}(x) = f_{\infty,L}(x) - \frac{x^{N-1}}{(2\pi i)^m} \oint_{C_R^m} \left( \prod_{j=1}^{m} (1-z_j)^{-L-1} z_j^{-N+1} \right) \frac{1}{z^{(m)} - x} d\vec{z}, \quad (6.38)$$

where

$$f_{\infty,L}(x) = \frac{1}{(2\pi i)^{m-1}} \oint_{C_R^{m-1}} \left(1 - \frac{x}{z_1 \dots z_{m-1}}\right)^{-L-1} \prod_{j=1}^{m-1} (1-z_j)^{-L-1} \frac{d\vec{z}}{z_1 \dots z_{m-1}}, \quad (6.39)$$

and the identity (6.39) follows from repeating the steps that led to (6.27) with $N = \infty$ and again applying Cauchy's integral formula. The estimation of the second term in (6.38) can be obtained by repeating the same steps that led to (6.35), using (6.32) and (6.37). In the complementary region $E^c$ we have $\mathrm{Re}(\phi(z)) = \gamma \log|1 + R\alpha - (R+1)\alpha e^{i\theta}| + \log|-R\alpha + (R+1)\alpha e^{i\theta}|$ which again is an increasing function of $\theta$ with a unique minimum at $\theta = 0$. This is easy to see from the formulas

$$|-R\alpha + (R+1)\alpha e^{i\theta}|^2 = \alpha^2(1 + 2R(R+1)(1 - \cos(\theta))), \quad (6.40)$$

$$|1 + R\alpha - (R+1)\alpha e^{i\theta}|^2 = (1 - \alpha)^2 + 2(R\alpha(1 + (\alpha+1)R) + \alpha)(1 - \cos(\theta)), \quad (6.41)$$
so that we again have a constant $c_\epsilon > 0$ such that $I_{E^c}/I_E = O\left(\frac{e^{-c_\epsilon N}}{\omega}\right)$. This concludes the proof of (2.12) when $-(\alpha + \delta_0)^m < x < (\alpha - \omega)^m$. $\square$

**Proof of Proposition 2.7.** In what follows, without loss of generality, we may take $m \ge 2$ since for $m = 1$ the proposition is immediate. Consider the integral representation for $f_{\infty,L}(x)$ given in (6.39) with contour C given by the circle of radius $x^\frac{1}{m}$ centered at the origin. For each $k = 1, \dots, m-1$ we parameterise $z_k = x^\frac{1}{m}e^{i\theta_k}$ with $\theta_k \in [-\pi, \pi]$. This gives

$$f_{\infty,L}(x) = \frac{1}{(2\pi)^{m-1}} \int_{[-\pi, \pi]^{m-1}} e^{-(L+1)\phi(\vec{\theta})} d\vec{\theta}, \quad (6.42)$$

where $\vec{\theta} = (\theta_1, \dots, \theta_{m-1})$, $d\vec{\theta} = d\theta_1 \dots d\theta_{m-1}$ and the action is given by

$$\phi(\vec{\theta}) = \log \left( 1 - x^\frac{1}{m} e^{-i(\theta_1 + \dots + \theta_{m-1})} \right) + \sum_{j=1}^{m-1} \log \left( 1 - x^\frac{1}{m} e^{i\theta_j} \right). \quad (6.43)$$

It is straightforward to show that $\vec{\theta} = \vec{0}$ is a critical point of $\phi(\vec{\theta})$. As we shall demonstrate below, the main contribution to the integral (6.42) comes from a small neighbourhood of this point that we denote $E = [-\epsilon, \epsilon]^{m-1}$ for $\epsilon > 0$ sufficiently small. The Hessian matrix of $\phi$ at $\vec{0}$ is given by

$$H_{ij} = \frac{x_m^{\frac{1}{m}}}{(1 - x_m^{\frac{1}{m}})^2} \times \begin{cases} 2, & i = j \\ 1, & i \neq j \end{cases} \quad (6.44)$$

where $1 \le i, j \le m-1$. On the set $E$, the considerations follow the proof of Theorem 2.2 identically, taking now $\eta := Lx_m^{\frac{1}{m}}$, and we omit the details. We restrict ourselves to showing that the contribution from the complementary region $E^c := [-\pi, \pi]^{m-1} \setminus [-\epsilon, \epsilon]^{m-1}$ is uniformly exponentially subleading and to proving the bounds (2.15) and (2.16). By Lemma B.3 we have

$$|I_{E^c}| = \left| \frac{1}{(2\pi)^{m-1}} \int_{E^c} \left(1 - x_m^{\frac{1}{m}} e^{-i \sum_{k=1}^{m-1} \theta_k}\right)^{-L-1} \prod_{j=1}^{m-1} d\theta_j \left(1 - x_m^{\frac{1}{m}} e^{i \theta_j}\right)^{-L-1} \right| \quad (6.45)$$

$$\le \left(1 - x_m^{\frac{1}{m}}\right)^{-m(L+1)} \frac{1}{(2\pi)^{m-1}} \int_{E^c} e^{-\frac{L}{2\pi^2} x_m^{\frac{1}{m}} \sum_{j=1}^{m-1} \theta_j^2} d\vec{\theta} \quad (6.46)$$

$$\le \left(1 - x_m^{\frac{1}{m}}\right)^{-m(L+1)} e^{-\frac{Lx_m^{\frac{1}{m}} \epsilon^2}{2\pi^2}}, \quad (6.47)$$

which shows that $I_{E^c}/I_E = O(e^{-c\eta})$ for some $c > 0$.

The inequality (2.15) for $x \in [0, 1)$ follows from

$$f_{\infty,L}(x) = \sum_{n=0}^{\infty} \left( \binom{L+n}{n} x_m^{\frac{n}{m}} \right)^m \le \left( \sum_{n=0}^{\infty} \binom{L+n}{n} x_m^{\frac{n}{m}} \right)^m = \frac{1}{\left(1 - x_m^{\frac{1}{m}}\right)^{m(L+1)}} \quad (6.48)$$

To prove (2.16), let $x \in [0, 1)$ and consider the integral representation for $f_{\infty,L}(-x)$ in (6.39) with contour C given by the circle of radius $x^\frac{1}{m}$ centered at the origin. Setting $z_k = x^\frac{1}{m}e^{i\theta_k}$ with $\theta_k \in [-\pi, \pi]$ for each $k = 1, \dots, m-1$ we split the domain $[-\pi, \pi]^{m-1}$ into disjoint subsets

$$A_{\ell} = \left\{ \vec{\theta} \in [-\pi, \pi]^{m-1} : 2\pi\ell \le \sum_{k=1}^{m-1} \theta_k < 2\pi(\ell + 1) \right\}, \quad (6.49)$$
where $\ell = -\lfloor \frac{m-1}{2} \rfloor - 1, \dots, \lfloor \frac{m-1}{2} \rfloor$. On the set $A_\ell$ we have

$$ \theta_{\ell} := \sum_{k=1}^{m-1} \theta_k - \pi - 2\ell\pi \in [-\pi, \pi]. \qquad (6.50) $$

Applying Lemma B.3, with angle (6.50) we find

$$ \begin{aligned} |f_{\infty,L}(-x)| & \le \sum_{\ell=-\lfloor \frac{m-1}{2} \rfloor-1}^{\lfloor \frac{m-1}{2} \rfloor} \frac{1}{(2\pi)^{m-1}} \int_{A_\ell} |1-x^\frac{1}{m}e^{-i\theta_\ell}|^{-L-1} \prod_{k=1}^{m-1} |1-x^\frac{1}{m}e^{i\theta_k}|^{-L-1} d\vec{\theta} \\ & \le \frac{1}{(1-x^\frac{1}{m})^{m(L+1)}} \frac{1}{(2\pi)^{m-1}} \sum_{\ell=-\lfloor \frac{m-1}{2} \rfloor-1}^{\lfloor \frac{m-1}{2} \rfloor} \int_{A_\ell} e^{-\frac{Lx^\frac{1}{m}}{2\pi^2}} \Psi(\vec{\theta}) d\vec{\theta}, \end{aligned} $$

where $\Psi(\vec{\theta}) = (\sum_{k=1}^{m-1} \theta_k - \pi - 2\pi\ell)^2 + \sum_{k=1}^{m-1} \theta_k^2$. The function $\Psi(\vec{\theta})$ is minimised at $\vec{\theta}^* = (\theta^*, \dots, \theta^*)$ where $\theta^* = \frac{\pi+2\pi\ell}{m}$. Since $\Psi(\vec{\theta}^*) = \frac{1}{m}(\pi + 2\pi\ell)^2 \ge \frac{\pi^2}{m}$ we obtain

$$ |f_{\infty,L}(-x)| \le \frac{1}{(1-x^\frac{1}{m})^{m(L+1)}} e^{-\frac{Lx^\frac{1}{m}}{2m}}. \qquad (6.51) \qquad \square $$

## A The real Ginibre ensemble

Although we focus in this paper on the case of truncated orthogonal matrices, our approach applies almost without change to a similar model defined as follows. Let $X^{(m)} = G_1 \dots G_m$ denote the product of $m$ independent real Ginibre matrices of size $N \times N$, that is each matrix in the product is independent and consists of independent standard (real) Gaussian random variables. Like the truncated orthogonal matrices, they are an integrable model whose eigenvalues (including the real ones) form a Pfaffian point process [28, 18].

The analogues of exact formulae (2.7) and (2.6) are known from the work of Forrester and Ipsen [18]. It is assumed in [18] that $N$ is even in this context, and we make the same assumption in what follows. There it is shown that

$$ E(N_{\mathbb{R}}^{(m)}) = \frac{N^{\frac{3m}{2}}}{(2\sqrt{2\pi})^m} \int_{\mathbb{R}} dx \int_{\mathbb{R}} dy |x-y| w_{\text{Gin}}(N^{m/2}x) w_{\text{Gin}}(N^{m/2}y) f_{N-2}(N^m xy) \quad (A.1) $$

where

$$ w_{\text{Gin}}(N^{m/2}x) = \int_{\mathbb{R}^m} \exp \left( -\frac{1}{2} \sum_{j=1}^m \lambda_j^2 \right) \delta(xN^{m/2} - \lambda_1 \dots \lambda_m) d\lambda_1 \dots d\lambda_m, \quad (A.2) $$

and

$$ f_{N-2}(N^m x) = \sum_{n=0}^{N-2} \frac{(N^m x)^n}{(n!)^m}. \qquad (A.3) $$

Using this formula, we can derive analogues of Theorems 1.1 and 1.2 with an identical approach used in the present paper.

**Theorem A.1.** Let $N_{\mathbb{R}}^{(m)}$ denote the number of real eigenvalues of $X^{(m)} = G_1 \dots G_m$. For any fixed $m \in \mathbb{N}$, we have

$$ E(N_{\mathbb{R}}^{(m)}) = \sqrt{\frac{2Nm}{\pi}} + O(1), \quad N \to \infty. \qquad (A.4) $$
The leading term in Theorem A.1 was first proved in [39] using an approach similar to the one used here in Section 5.

**Remark A.2.** When $m = 1$ a precise asymptotic expansion is known for $E(N_R)$, see [15], and this shows that the $O(1)$ estimate in (A.4) is optimal. Similarly, we believe that our $O(1)$ estimates in (A.4) and (1.2) are optimal for any fixed $m \ge 1$. In [39] the estimate on the remainder in (A.4) is at the less precise order $O(\log(N))$, hence Theorem A.1 gives an improvement on this result.

Now again from [18] we have that the appropriately scaled and normalised density of real eigenvalues of the product matrix $N^{-\frac{m}{2}}G_1...G_m$ is given by

$$ \tilde{\rho}_N(x) = \frac{1}{E(N_R^{(m)})} \frac{N^{\frac{3m}{2}}}{(2\sqrt{2\pi})^m} \int_R dy |x-y| w_{\text{Gin}}(N^{\frac{m}{2}}x) w_{\text{Gin}}(N^{\frac{m}{2}}y) f_{N-2}(N^m x y). \quad (A.5) $$

**Theorem A.3.** For any bounded continuous test function $h$, we have

$$ \lim_{N \to \infty} \int_{\mathbb{R}} h(x) \tilde{\rho}_N(x) dx = \int_{\mathbb{R}} h(x) \rho(x) dx, \quad (A.6) $$

where

$$ \rho(x) = \frac{1}{2m x^{1-\frac{1}{m}}} \mathbb{1}_{x \in (-1,1)}. \quad (A.7) $$

Furthermore, for any fixed $x \in \mathbb{R} \setminus \{-1, 0, 1\}$, we have the pointwise limit $\lim_{N \to \infty} \tilde{\rho}_N(x) = \rho(x)$.

The weak convergence (A.6) was also established in [39], proving a conjecture of Forrester and Ipsen in [18].

Let us now explain the small changes required to prove Theorems A.1 and A.3 using the techniques of the present paper. First note that the definition of the weight (A.2) can be written in the equivalent form

$$ w_{\text{Gin}}(N^{\frac{m}{2}} x) = 2^{m-1} \int_{[0, \infty)^{m-1}} \exp \left( -\frac{N|x|^{\frac{2}{m}}}{2} \left( \frac{1}{u_1^2 \dots u_{m-1}^2} + \sum_{j=1}^{m-1} u_j^2 \right) \right) \frac{d\vec{u}}{u^{(m-1)}}, \quad (A.8) $$

where $u^{(m-1)} = u_1 \dots u_{m-1}$. Then a saddle point analysis as in the proof of Proposition 2.5 gives the following result.

**Proposition A.4.** Fix a large constant $M > 0$. Then we have the following asymptotic estimate as $N \to \infty$, uniformly on $|x| \in [MN^{-\frac{m}{2}}, \infty)$,

$$ w_{\text{Gin}}(N^{\frac{m}{2}} x) = N^{-\frac{m-1}{2}} e^{-\frac{N m}{2} |x|^{\frac{2}{m}}} \frac{(4\pi)^{\frac{m-1}{2}}}{\sqrt{m}} |x|^{-\frac{m-1}{m}} \left( 1 + O\left( \frac{1}{N|x|^{\frac{2}{m}}} \right) \right). \quad (A.9) $$

In fact the proof of Proposition A.4 is more straightforward than for Proposition 2.5, because it is already immediate from (A.8) that $\eta := N|x|^{\frac{2}{m}}$ is the appropriate large parameter in the Laplace asymptotics.

Regarding $f_{N-2}(N^m x)$ in (A.3), we mimic the proof of Proposition 2.6 almost identically. For instance, instead of (6.24), we use the integral representation

$$ \frac{N^n}{n!} = \frac{1}{2\pi i} \oint_C e^{Nz} z^{-n-1} dz, \quad (A.10) $$

and follow all the steps identically, the main difference being that the saddle point is now located at $z = 1$ instead of $z = \alpha$. This leads to
**Proposition A.5.** As $N \to \infty$ we have the following estimate uniformly on $x \in \mathbb{R} \setminus ((1-\omega)^m, (1+\omega)^m)$,

$$f_{N-2}(N^m x) = f_{\infty}(N^m x) \mathbb{1}_{-(1+\omega)^m < x < (1-\omega)^m} + \frac{x^{N-1}}{x-1} e'_{N,m} \left( 1 + O \left( \frac{1}{\sqrt{N\omega}} \right) \right), \quad (A.11)$$

where we take $\omega = N^{-\frac{1}{2}}$ and $e'_{N,m} = e^{mN}(2\pi N)^{-\frac{m}{2}}$.

Now regarding the function,

$$f_{\infty}(N^m x) = \sum_{n=0}^{\infty} \frac{(N^m x)^n}{(n!)^m}, \qquad (A.12)$$

we insert the integral representation (A.10) and find that for $x > 0$, we have

$$f_{\infty}(N^m x) = \frac{1}{(2\pi i)^{m-1}} \oint_{C^{m-1}} \exp \left( N x^{\frac{1}{m}} \left( \sum_{k=1}^{m-1} z_k + \frac{1}{z_1 \dots z_{m-1}} \right) \right) \frac{dz}{z^{(m-1)}}, \quad (A.13)$$

where $z^{(m-1)} = z_1 \dots z_{m-1}$. The representation (A.13) is well-suited for applying the Laplace method with large parameter $\eta := N x^{\frac{1}{m}}$, in complete analogy with the proof of Proposition 2.7.

**Proposition A.6.** Fix a large constant $M > 0$. Then we have the following estimate as $N \to \infty$, uniformly on $x \in (MN^{-m}, \infty)$,

$$f_{\infty}(N^m x) = (2\pi)^{-\frac{m-1}{2}} x^{-\frac{m-1}{2m}} e^{Nmx} N^{-\frac{m-1}{2}} \frac{1}{\sqrt{m}} \left( 1 + O \left( \frac{1}{x^{\frac{1}{m}}N} \right) \right). \quad (A.14)$$

**Remark A.7.** Analogously to Remark 2.8, at the pointwise level the asymptotics (A.9) and (A.14) were derived in the context of products of independent complex Ginibre matrices in [2].

**Remark A.8.** We make a remark on the $m = 1$ case of Proposition A.5. This corresponds to studying asymptotics of a truncated exponential function and has been considered by many authors. At the level of precision required in this paper, this appeared in [9], but more precise asymptotics are known. The proof given in [9] directly motivated our proof of Propositions 2.6 and A.5 that hold for any $m \ge 1$.

Now with Propositions A.4, A.5 and A.6, the proofs in Sections 2.2 and 3 go through in complete analogy and we omit the details. This is how we prove Theorems A.1 and A.3.

## B Miscellaneous bounds

**Lemma B.1.** For any $(u, v) \in [0, 1]^2$, we have the inequality

$$\frac{(1-u^2)(1-v^2)}{(1-uv)^2} \le e^{-(u-v)^2}. \qquad (B.1)$$

*Proof.* Using the bound $1+x \le e^x$ for all $x \in \mathbb{R}$ we have

$$\frac{(1-u^2)(1-v^2)}{(1-uv)^2} = 1 - \frac{(u-v)^2}{(1-uv)^2} \le e^{-(\frac{u-v}{1-uv})^2} \le e^{-(u-v)^2}, \quad (B.2)$$

where in the last inequality we used $\frac{1}{1-uv} \ge 1$ for all $(u, v) \in [0, 1]^2$. $\square$
**Lemma B.2.** Consider the function $Q : [0, \infty)^2 \to [0, \infty)$ defined by

$$Q(u, v) = \frac{(u + v)^m - u^m}{(u(u + v))^{\frac{m-1}{2}}}. \qquad (B.3)$$

Then for $\delta > 0$ sufficiently small and $M > 0$ arbitrary, we have for $(u, v) \in [0, M] \times [0, \delta]$ the uniform bound

$$Q(u, v) = mv + \sum_{j=0}^{m-1} O \left( \frac{v^{3+j}}{u^{2+j}} \right). \qquad (B.4)$$

Furthermore, for any $\epsilon > 0$, there is a constant $C_{\epsilon,M} > 0$ independent of $u$ such that on the domain $(u, v) \in [\epsilon, M] \times [0, M]$ we have $Q(u, v) \le C_{\epsilon,M}v$.

*Proof.* We have

$$Q(u, v) = u \frac{\left(1 + \frac{v}{u}\right)^m - 1}{\left(1 + \frac{v}{u}\right)^{\frac{m-1}{2}}} = \left(1 + \frac{v}{u}\right)^{-\frac{m-1}{2}} \sum_{j=0}^{m-1} \binom{m}{j+1} \frac{v^{j+1}}{u^j}. \qquad (B.5)$$

This implies that $Q(u, v) \le C_{\epsilon,M}v$ provided $u > \epsilon$. Taylor expanding near $v = 0$, we have

$$(1 + \frac{v}{u})^{-\frac{m-1}{2}} = 1 - \frac{m-1}{2} \frac{v}{u} + O\left(\frac{v^2}{u^2}\right), \quad 0 < v < \delta. \qquad (B.6)$$

To obtain the uniform big-O term in (B.6) note that

$$\frac{d^2}{dv^2} \left(1 + \frac{v}{u}\right)^{-\frac{m-1}{2}} = \frac{\left(1 + \frac{v}{u}\right)^{-\frac{m-1}{2}}}{4(u + v)^2} (m^2 - 1) \le \frac{m^2 - 1}{4u^2}. \qquad (B.7)$$

Then inserting (B.6) into (B.5) the term proportional to $\frac{v^2}{u}$ cancels and we obtain (B.4). $\square$

**Lemma B.3.** For $\theta \in [-\pi, \pi]$ and $0 \le x < 1$, we have the inequality

$$\frac{1}{\left|1 - x^{\frac{1}{m}} e^{i\theta}\right|} \le \frac{1}{1 - x^{\frac{1}{m}}} e^{-\frac{x^{\frac{1}{m}}}{2\pi^2} \theta^2}. \qquad (B.8)$$

*Proof.* The identity

$$\frac{\left(1 - x^{\frac{1}{m}}\right)^2}{\left|1 - x^{\frac{1}{m}} e^{i\theta}\right|^2} = \frac{1}{1 + \mathcal{A}(1 - \cos(\theta))}, \qquad \mathcal{A} := \frac{2x^{\frac{1}{m}}}{\left(1 - x^{\frac{1}{m}}\right)^2}, \qquad (B.9)$$

and the inequalities $\frac{1}{1+y} \le e^{-\frac{y}{1+y}}$ for $y \ge 0$ and $(1 - \cos(\theta)) \le 2$ show that

$$\frac{\left(1 - x^{\frac{1}{m}}\right)^2}{\left|1 - x^{\frac{1}{m}} e^{i\theta}\right|^2} \le e^{-\frac{\mathcal{A}}{1+2\mathcal{A}}(1-\cos(\theta))}. \qquad (B.10)$$

To conclude, we apply Jordan's inequality $1 - \cos(\theta) \ge \frac{2}{\pi^2} \theta^2$ for $\theta \in [-\pi, \pi]$ and the following bound to the right-hand side of (B.10),

$$\frac{\mathcal{A}}{1 + 2\mathcal{A}} = \frac{2x^{\frac{1}{m}}}{(1 - x^{\frac{1}{m}})^2 + 4x^{\frac{1}{m}}} \ge \frac{1}{2} x^{\frac{1}{m}}. \qquad (B.11)$$

$\square$
**Lemma B.4.** Let the circles $C_R$ be defined by

$$z_k = -R\alpha + (R+1)\alpha e^{i\theta_k}, \quad -\pi < \theta_k \le \pi, \qquad (B.12)$$

and define $z^{(m)} := z_1 \dots z_m$. Consider the parameters $\omega$ and $\alpha$ as in Proposition 2.6. Then for $R > 0$ large enough, there exist absolute constants $c > 0$ and $\delta_0 > 0$ independent of $N$ such that

$$\inf_{x \in [-(\alpha + \delta_0)^m, (\alpha - \omega)^m]} \inf_{\theta_k \in [-\pi, \pi], k=1, \dots, m} |z^{(m)} - x| > c\omega. \qquad (B.13)$$

Furthermore, for any $x \in [-(\alpha + \delta_0)^m, (\alpha - \omega)^m]$ and for each $k = 1, \dots, m$, we have that $\frac{xz_k}{z^{(m)}}$ belongs to the interior of $C_R$.

*Proof.* We have

$$|z_k|^2 = 2\alpha^2 R^2 (1 - \cos(\theta_k)) + 2R\alpha^2 (1 - \cos(\theta_k)) + \alpha^2, \qquad (B.14)$$

and clearly each $|z_k| \ge \alpha$. When $x > 0$ we have

$$|z^{(m)} - x| \ge ||z^{(m)}| - |x|| \ge \alpha^m - (\alpha - \omega)^m \ge c\omega. \qquad (B.15)$$

We also have

$$\left| \frac{xz_{\ell}}{z^{(m)}} \right| \le \frac{(\alpha - \omega)^m}{\alpha^{m-1}} \le \alpha(1 - \omega/\alpha)^m < \alpha, \qquad (B.16)$$

and therefore $\frac{xz_{\ell}}{z^{(m)}}$ belongs to the interior of $C_R$, as $C_R$ strictly includes any disc of radius smaller than $\alpha$.

If $x \le 0$, in particular when $x$ is close to $-\alpha^m$ we have to check that there is not too much winding in the product $z^{(m)}$ that might allow $\phi := \text{Arg}(z^{(m)}) = \pi$. We start by supposing that for at least one $k = 1, \dots, m$, we have $|\theta_k| \ge \epsilon/R$. Then for $R$ large enough and any fixed $\epsilon > 0$

$$1 - \cos(\theta_k) \ge 1 - \cos(\epsilon/R) \ge \left(\frac{\epsilon}{R}\right)^2 \frac{1}{\pi}. \qquad (B.17)$$

Hence, for this $k$, inserting (B.17) into (B.14) gives $|z_k|^2 \ge \alpha^2(1 + \epsilon^2/2)$. Using this and $|z_j| \ge \alpha$ shows that $|z^{(m)}| \ge \alpha^m(1 + \epsilon^2/8)$. Furthermore $|x| < \alpha^m(1 + c_m\delta_0)$ for some $c_m > 0$ depending only on $m$. Then

$$|z_1 \dots z_m - x| \ge ||z_1| \dots |z_m| - |x|| \ge \alpha^m (\epsilon^2/8 - c_m \delta_0) > c\omega, \qquad (B.18)$$

which is valid for any $\epsilon > \sqrt{8/c_m\delta_0 + \omega}$. Thus we see that by choosing $\delta_0$ small enough we can allow arbitrarily small $\epsilon > 0$. The same bound shows that for any $\ell = 1, \dots, m$, we have

$$|xz_{\ell}/z^{(m)}| \le \alpha \frac{1 + c_m \delta_0}{1 + \frac{\epsilon^2}{8}} < \alpha, \qquad (B.19)$$

which implies that $xz^{\ell}/z^{(m)}$ is strictly in the interior of $C_R$. This bounds the region where at least one of the $|\theta_k| \ge \epsilon/R$. The complement of this region is where $|\theta_j| < \epsilon/R$ for all $j = 1, \dots, m$. We claim that on this region

$$|\text{Arg}(z_1 \dots z_m)| = |\text{Arg}(z_1) + \dots + \text{Arg}(z_m)| \le 2\epsilon m. \qquad (B.20)$$

To prove this we compute the argument of each $z_j$,

$$\text{Arg}(z_j) = \tan^{-1} \left( \frac{\sin(\theta_j)}{-\frac{R}{R+1} + \cos(\theta_j)} \right). \qquad (B.21)$$
By symmetry it suffices to assume $\theta_j \ge 0$. Using $\sin(x) \le x$, we have for the numerator $\sin(\theta_j) \le \theta_j \le \epsilon/R$. The cosine inequality $\cos(x) \ge 1 - 2x/\pi$ gives for the denominator $-\frac{R}{R+1} + \cos(\theta_j) \ge \frac{1}{R+1} - 2\epsilon/(\pi R)$ and clearly $\frac{R}{R+1} - 2\epsilon/\pi > \frac{1}{2}$. These bounds imply that

$$ \frac{\sin(\theta_j)}{-\frac{R}{R+1} + \cos(\theta_j)} \le 2\epsilon. \qquad (B.22) $$

Now by monotonicity of $\tan^{-1}$ we have $\mathrm{Arg}(z_j) \le 2\epsilon$ and this implies (B.20). We thus have

$$ |z^{(m)} - x|^2 = |z|^2 + x^2 - 2|z|x \cos(\phi) \ge \alpha^{2m} > \omega, \qquad (B.23) $$

which holds for example if $\phi := \mathrm{Arg}(z^{(m)}) < \pi/2$ (note that (B.20) implies $\phi < 2\epsilon m$). We also have

$$ \begin{gathered} |\mathrm{Arg}(z_{\ell}x/z^{(m)}) - \pi| \le 2m\epsilon, \\ \left|\frac{z_{\ell}x}{z^{(m)}}\right| \le \alpha(1 + \delta_0)^m, \end{gathered} \qquad (B.24) $$

and choosing $\delta_0$ and $\epsilon$ small enough implies the strict inclusion of $z_{\ell}x/z^{(m)}$ inside the contour $C_R$. $\square$

## C Exact calculation of a multiple integral

In this section we compute the integral from (5.9) explicitly, namely the integral

$$ A_{L,m} = \frac{1}{\Gamma\left(\frac{L}{2}\right)^{2m}} \int_{\mathbb{R}_+^{2m}} d\vec{t} d\vec{r} \prod_{\ell=1}^{m} t_{\ell}^{\frac{L}{2}-1} r_{\ell}^{\frac{L}{2}-1} e^{-t_{\ell}-r_{\ell}} \frac{t_1 + \dots + t_m}{2} \mathbb{1}_{r_1+\dots+r_m < t_1+\dots+t_m}. \qquad (C.1) $$

**Lemma C.1.** We have

$$ A_{L,m} = \frac{mL}{8} + \frac{1}{2} \frac{\Gamma(mL)}{\Gamma\left(\frac{mL}{2}\right)^2} 2^{-mL}. \qquad (C.2) $$

*Proof.* By permutation invariance it is sufficient to consider the following

$$ B_{L,m} := \frac{1}{\Gamma\left(\frac{L}{2}\right)^{2m}} \int_{\mathbb{R}_+^{2m}} d\vec{t} d\vec{r} \prod_{\ell=1}^{m} t_{\ell}^{\frac{L}{2}-1} r_{\ell}^{\frac{L}{2}-1} e^{-t_{\ell}-r_{\ell}} t_m \mathbb{1}_{r_1+\dots+r_m < t_1+\dots+t_m}, \qquad (C.3) $$

so that clearly $A_{L,m} = \frac{m}{2}B_{L,m}$. We would also prefer to replace the $t_m$ with $r_m$. Note that if we replaced $t_m$ with $t_m + r_m$, then by symmetry between the $r_\ell$ and $t_\ell$ variables we can drop the indicator function after multiplying by a factor $\frac{1}{2}$. The remaining integrals are explicit and we have

$$ \frac{1}{\Gamma\left(\frac{L}{2}\right)^{2m}} \int_{\mathbb{R}_+^{2m}} d\vec{t} d\vec{r} \prod_{\ell=1}^{m} t_{\ell}^{\frac{L}{2}-1} r_{\ell}^{\frac{L}{2}-1} e^{-t_{\ell}-r_{\ell}} (t_m + r_m) \mathbb{1}_{r_1+\dots+r_m < t_1+\dots+t_m} = \frac{L}{2}. \qquad (C.4) $$

We thus have $B_{L,m} = \frac{L}{2} - F_{L,m}$ where

$$ F_{L,m} := \frac{1}{\Gamma\left(\frac{L}{2}\right)^{2m}} \int_{\mathbb{R}_+^{2m}} d\vec{t} d\vec{r} \prod_{\ell=1}^{m} t_{\ell}^{\frac{L}{2}-1} r_{\ell}^{\frac{L}{2}-1} e^{-t_{\ell}-r_{\ell}} r_m \mathbb{1}_{r_1+\dots+r_m < t_1+\dots+t_m}. \qquad (C.5) $$

To compute $F_{L,m}$, we will integrate by parts in the variable $r_m$, differentiating the factor $r_m^{\frac{L}{2}}$ and integrating $e^{-r_m}$. The second term in the integration by parts is thus completely explicit, again by symmetry, and gives a contribution $\frac{L}{4}$. Thus we have $F_{L,m} = \frac{L}{4} - D_{L,m}$
where $D_{L,m}$ is the boundary term of the integration by parts. To calculate this, we parameterise the integration such that $\vec{t} \in \mathbb{R}_+^m$ are unconstrained, while $\vec{r}$ satisfy the constraints $K_{m+1-\ell} := \{0 < r_{m+1-\ell} < t_1 + \dots + t_m - (r_1 + \dots + r_{m-\ell})\}$ where $\ell = 1, \dots, m$. Setting $r_m = t_1 + \dots + t_m - (r_1 + \dots + r_{m-1})$ results in the boundary term,

$$D_{L,m} = \frac{1}{\Gamma\left(\frac{L}{2}\right)^{2m}} \int_{\mathbb{R}_+^m} d\vec{t} \prod_{\ell=1}^{m} t_{\ell}^{\frac{L}{2}-1} e^{-2t_{\ell}} \prod_{\ell=1}^{m-1} \int_{K_{\ell}} dr_{\ell} r_{\ell}^{\frac{L}{2}-1} (t^{(m)} - r^{(m-1)})^{\frac{L}{2}-1}, \quad (C.6)$$

where $t^{(m)} = t_1 + \dots + t_m$ and $r^{(m-1)} = r_1 + \dots + r_{m-1}$. Then we integrate over each $r_\ell$ variable, starting with $r_{m-1}$, repeatedly using the formula

$$\int_0^u dr_\ell r_\ell^{a-1} (u-r_\ell)^{b-1} = u^{a+b-1} \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}, \quad (C.7)$$

for each integral. We obtain

$$D_{L,m} = \frac{1}{\Gamma\left(\frac{L}{2}\right)^{2m}} \left( \prod_{\ell=1}^{m-1} \frac{\Gamma\left(\frac{\ell L}{2} + 1\right) \Gamma\left(\frac{L}{2}\right)}{\Gamma\left(\frac{\ell L}{2} + \frac{L}{2} + 1\right)} \right) \int_{\mathbb{R}_+^m} d\vec{t} \prod_{\ell=1}^{m} t_{\ell}^{\frac{L}{2}-1} e^{-2t_{\ell}} (t_1 + \dots + t_m)^{\frac{mL}{2}}. \quad (C.8)$$

Now to compute the integrals over the $t_\ell$ variables we substitute the $t_m$ variable with $u = t_1 + \dots + t_m$ and parameterise such that $u \in \mathbb{R}_+$, while $t_1, \dots, t_{m-1}$ satisfy the constraints $H_{m-j} = \{0 < t_{m-j} < u - (t_1 + \dots + t_{m-j-1})\}$ where $j = 1, \dots, m-1$. The integrals over $H_{m-j}$ are computed again using (C.7) repeatedly and we find

$$D_{L,m} = \frac{1}{\Gamma\left(\frac{L}{2}\right)^{2m}} \left( \prod_{\ell=1}^{m-1} \frac{\Gamma\left(\frac{\ell L}{2} + 1\right) \Gamma\left(\frac{L}{2}\right)}{\Gamma\left(\frac{\ell L}{2} + \frac{L}{2} + 1\right)} \frac{\Gamma\left(\frac{\ell L}{2}\right) \Gamma\left(\frac{L}{2}\right)}{\Gamma\left(\frac{\ell L}{2} + \frac{L}{2}\right)} \right) \int_0^\infty du \, u^{mL-1} e^{-2u} \quad (C.9)$$

$$= \frac{1}{\Gamma\left(\frac{L}{2}\right)^{2m}} \left( \prod_{\ell=1}^{m-1} \frac{\Gamma\left(\frac{\ell L}{2} + 1\right) \Gamma\left(\frac{L}{2}\right)}{\Gamma\left(\frac{\ell L}{2} + \frac{L}{2} + 1\right)} \frac{\Gamma\left(\frac{\ell L}{2}\right) \Gamma\left(\frac{L}{2}\right)}{\Gamma\left(\frac{\ell L}{2} + \frac{L}{2}\right)} \right) 2^{-mL} \Gamma(mL) \quad (C.10)$$

$$= \frac{\Gamma(mL)}{m\Gamma(\frac{mL}{2})^2} 2^{-mL}, \quad (C.11)$$

where we exploited the cancellation of successive terms in the products (C.10). Substituting (C.11) in the formula $A_{L,m} = \frac{mL}{8} + \frac{m}{2} D_{L,m}$ concludes the proof of the Lemma. $\square$

## References

[1] K. Adhikari, N. Kishore Reddy, T. Ram Reddy, and K. Saha. Determinantal point processes in the plane from products of random matrices. Ann. Inst. Henri Poincaré Probab. Stat., 52(1):16–46, 2016. MR3449292

[2] G. Akemann and Z. Burda. Universal microscopic correlation functions for products of independent Ginibre matrices. J. Phys. A, 45(46):465201, 18, 2012. MR2993423

[3] G. Akemann, Z. Burda, and M. Kieburg. Universality of local spectral statistics of products of random matrices. Phys. Rev. E., 102:052134, 2020. MR4189356

[4] G. Akemann, Z. Burda, M. Kieburg, and T. Nagao. Universal microscopic correlation functions for products of truncated unitary matrices. J. Phys. A, 47(25):255202, 2014. MR3224113

[5] G. Akemann and J. R. Ipsen. Recent exact and asymptotic results for products of independent random matrices. Acta. Phys. Pol. B, 46(9):1747–1784, 2015. MR3403839

[6] J. Baik and T. Bothner. The largest real eigenvalue in the real Ginibre ensemble and its relation to the Zakharov-Shabat system. Ann. Appl. Probab., 30(1):460–501, 2020. MR4068316

[7] A. Borodin and C. D. Sinclair. The Ginibre ensemble of real random matrices and its scaling limits. Comm. Math. Phys., 291(1):177–224, 2009. MR2530159
[8] P. Bougerol and J. Lacroix. *Products of random matrices with applications to Schrödinger operators*, volume 8 of *Progress in probability and statistics (P. Huber and M. Rosenblatt, eds.)*. Birkhäuser, Boston, 1985. MR0886674

[9] R. Boyer and W. M. Y. Goh. On the zero attractor of the Euler polynomials. *Adv. in Appl. Math.*, 38(1):97-132, 2007. MR2288197

[10] Z. Burda, R. A. Janik, and B. Waclaw. Spectrum of the product of independent random Gaussian matrices. *Phys. Rev. E* (3), 81(4):041132, 12, 2010. MR2736204

[11] Z. Burda, M. Nowak, and A. Swiech. New spectral relations between products and powers of isotropic random matrices. *Phys. Rev. E.*, 86:061137, 2012.

[12] A. Crisanti, G. Paladin, and A. Vulpiani. Products of random matrices. In *Random Matrices and their applications*, volume 104 of *Springer series in Solid-State Sciences*. Springer-Verlag, Berlin, Heidelberg, 1993. MR1278483

[13] L. C. G. del Molino, K. Pakdaman, and J. Touboul. Real eigenvalues of non-symmetric random matrices: Transitions and Universality. eprint = arXiv:1605.00623[math-ph].

[14] A. Edelman. The probability that a random real Gaussian matrix has $k$ real eigenvalues, related distributions, and the circular law. *J. Multivariate Anal.*, 60(2):203-232, 1997. MR1437734

[15] A. Edelman, E. Kostlan, and M. Shub. How many eigenvalues of a random matrix are real? *J. Amer. Math. Soc.*, 7(1):247-267, 1994. MR1231689

[16] W. FitzGerald, R. Tribe, and O. Zaboronski. Sharp asymptotics for Fredholm Pfaffians related to interacting particle systems and random matrices. *Electron. J. Probab.*, 25:1-15, 2020. MR4161126

[17] P. J. Forrester. The limiting Kac random polynomial and truncated random orthogonal matrices. *J. Stat. Mech: Theo. Exper.*, P12018, 2010.

[18] P. J. Forrester and J. R. Ipsen. Real eigenvalue statistics for products of asymmetric real Gaussian matrices. *Linear Algebra and its Applications*, 510:259-290, 2016. MR3551633

[19] P. J. Forrester, J. R. Ipsen, and S. Kumar. How many eigenvalues of a product of truncated orthogonal matrices are real? *Exp. Math.*, 29(3):276-290, 2020. MR4134827

[20] P. J. Forrester and S. Kumar. The probability that all eigenvalues are real for products of truncated real orthogonal random matrices. *J. Theoret. Probab.*, 31(4):2056-2071, 2018. MR3866607

[21] P. J. Forrester and T. Nagao. Eigenvalue Statistics of the Real Ginibre Ensemble. *Phys. Rev. Lett.*, 99:050603, 2007.

[22] H. Furstenberg and H. Kesten. Products of random matrices. *Ann. Math. Statist.*, 31:457-469, 1960. MR0121828

[23] Y. V. Fyodorov and B. A. Khoruzhenko. A Nonlinear Analogue of May-Wigner Instability Transition. *Proc. Natl. Acad. Sci. USA*, 113:6827-6832, 2016. MR3521630

[24] M. Gebert and M. Poplavskyi. On pure complex spectrum for truncations of random orthogonal matrices and Kac polynomials. eprint = arXiv:1905.03154.

[25] J. Ginibre. Statistical ensembles of complex, quaternion, and real matrices. *J. Mathematical Phys.*, 6:440-449, 1965. MR0173726

[26] S. Hameed, K. Jain, and A. Lakshminarayan. Real eigenvalues of non-Gaussian random matrices and their products. *J. Phys. A*, 48(38):385204, 26, 2015. MR3400913

[27] G. C. P. Innocentini and M. Novaes. Time-inhomogeneous random Markov chains. *J. Stat. Mech. Theory Exp.*, (10):103202, 14, 2018. MR3886008

[28] J. Ipsen and M. Kieburg. Weak commutation relations and eigenvalue statistics for products of rectangular random matrices. *Phys. Rev. E.*, 89:032106, 2014.

[29] B. A. Khoruzhenko, H.-J. Sommers, and K. Życzkowski. Truncations of random orthogonal matrices. *Phys. Rev. E* (3), 82(4):040106, 4, 2010. MR2788023

[30] P. Kopel. Linear statistics of non-Hermitian matrices matching the real or complex Ginibre ensemble to four moments. eprint = arXiv:1510.02987.
[31] A. Lakshminarayan. On the number of real eigenvalues of products of random matrices and an application to quantum entanglement. *J. Phys. A: Math. Theor.*, 46:152003, 2013. MR3043878

[32] N. Lehmann and H.-J. Sommers. Eigenvalue statistics of random real matrices. *Phys. Rev. Lett.*, 67:941-944, 1991. MR1121461

[33] D.-Z. Liu, D. Wang, and Y. Wang. Lyapunov exponent, universality and phase transition for products of random matrices. eprint = arXiv:1810.00433.

[34] D.-Z. Liu and Y. Wang. Phase transitions for infinite products of large non-Hermitian random matrices. eprint = arXiv:1912.11910.

[35] D.-Z. Liu and Y. Wang. Universality for products of random matrices I: Ginibre and truncated unitary cases. *Int. Math. Res. Not.* IMRN, (11):3473-3524, 2016. MR3556416

[36] M. Poplavskyi and G. Schehr. Exact persistence exponent for the 2d-diffusion equation and related Kac polynomials. *Phys. Rev. Lett.*, 121:150601, 2018.

[37] M. Poplavskyi, R. Tribe, and O. Zaboronski. On the distribution of the largest real eigenvalue for the real Ginibre ensemble. *Ann. Appl. Probab.*, 27(3):1395-1413, 2017. MR3678474

[38] N. Simm. Central limit theorems for the real eigenvalues of large Gaussian random matrices. *Random Matrices Theory Appl.*, 6(1):1750002, 18, 2017. MR3612267

[39] N. Simm. On the real spectrum of a product of Gaussian matrices. *Electron. Commun. Probab.*, 22:1-11, 2017. MR3685239

[40] H.-J. Sommers and W. Wieczorek. General eigenvalue correlations for the Ginibre ensemble. *J. Phys. A: Math. Theor.*, 41(40), 2008. MR2439268

[41] T. Tao and V. Vu. Random matrices: universality of local spectral statistics of non-Hermitian matrices. *Ann. Probab.*, 43(2):782-874, 2015. MR3306005

[42] R. Tribe, S. K. Yip, and O. Zaboronski. One dimensional annihilating and coalescing particle systems as extended Pfaffian point processes. *Electron. Commun. Probab.*, 17:no. 40, 7, 2012. MR2981896

[43] R. Tribe and O. Zaboronski. Pfaffian formulae for one dimensional coalescing and annihilating systems. *Electron. J. Probab.*, 16:no. 76, 2080-2103, 2011. MR2851057

**Acknowledgments.** We are grateful to an anonymous referee who pointed out to us references [14, 32] and whose comments improved the presentation of the paper.