Condition Numbers for Real Eigenvalues in the Real Elliptic Gaussian Ensemble

We study the distribution of the eigenvalue condition numbers κi=(li∗li)(ri∗ri)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa _i=\sqrt{ ({\mathbf{l}}_i^* {\mathbf{l}}_i)({\mathbf{r}}_i^* {\mathbf{r}}_i)}$$\end{document} associated with real eigenvalues λi\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _i$$\end{document} of partially asymmetric N×N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\times N$$\end{document} random matrices from the real Elliptic Gaussian ensemble. The large values of κi\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa _i$$\end{document} signal the non-orthogonality of the (bi-orthogonal) set of left li\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf{l}}_i$$\end{document} and right ri\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf{r}}_i$$\end{document} eigenvectors and enhanced sensitivity of the associated eigenvalues against perturbations of the matrix entries. We derive the general finite N expression for the joint density function (JDF) PN(z,t)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal {P}}}_N(z,t)$$\end{document} of t=κi2-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=\kappa _i^2-1$$\end{document} and λi\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _i$$\end{document} taking value z, and investigate its several scaling regimes in the limit N→∞\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\rightarrow \infty $$\end{document}. When the degree of asymmetry is fixed as N→∞\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\rightarrow \infty $$\end{document}, the number of real eigenvalues is O(N)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {O}(\sqrt{N})$$\end{document}, and in the bulk of the real spectrum ti=O(N)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_i=\mathcal {O}(N)$$\end{document}, while on approaching the spectral edges the non-orthogonality is weaker: ti=O(N)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_i=\mathcal {O}(\sqrt{N})$$\end{document}. In both cases the corresponding JDFs, after appropriate rescaling, coincide with those found in the earlier studied case of fully asymmetric (Ginibre) matrices. A different regime of weak asymmetry arises when a finite fraction of N eigenvalues remain real as N→∞\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\rightarrow \infty $$\end{document}. In such a regime eigenvectors are weakly non-orthogonal, t=O(1)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=\mathcal {O}(1)$$\end{document}, and we derive the associated JDF, finding that the characteristic tail P(z,t)∼t-2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal {P}}}(z,t)\sim t^{-2}$$\end{document} survives for arbitrary weak asymmetry. As such, it is the most robust feature of the condition number density for real eigenvalues of asymmetric matrices.


Introduction
A (real-valued) square matrix X is asymmetric if it is different from its transpose X T , and non-normal if XX T = X T X. Generically, asymmetric matrices are non-normal, and their eigenvalues are much more sensitive to the perturbations of the matrix entries than for their symmetric (hence selfadjoint and normal) counterparts. It is well known that non-normality may raise serious issues when calculating the spectra of such matrices numerically: keeping a fixed precision of calculations might not be sufficient, as some eigenvalues can be 'ill-conditioned'.
To be more specific, we assume that X can be diagonalized (which for random matrices happens with probability one). Then to each eigenvalue λ i , real or complex (in the latter case being always accompanied by its complex conjugate partner λ i ) correspond two sets of eigenvectors, left l i and right r i which can always be chosen to be bi-orthogonal: l * i r j = δ ij , where l * i := l T i stands for Hermitian conjugation. The corresponding eigenproblems are Xr i = λ i r i and X T l i = λ i l i . Consider now a matrix X = X + P , where the second term represents an error one makes by storing the matrix entries with a finite precision, with > 0 controlling the magnitude of the error and P reflecting the matrix structure of the perturbation. In the first order perturbation theory in parameter eigenvalues are shifted by The latter quantity, κ i = (l * i l i )(r * i r i ), shows that the sensitivity of eigenvalues is essentially controlled by the non-orthogonality of the corresponding left and right eigenvectors. Correspondingly, in the numerical analysis context κ i is called the eigenvalue condition number of the eigenvalue λ i [41,42]. Note also that the Cauchy-Schwarz inequality implies κ ≥ 1, with κ = 1 only when X is normal.
It is natural to ask how well-conditioned a 'typical' asymmetric matrix is. This question can be most meaningfully answered in the context of Random Matrix Theory (RMT), by defining 'typical' as randomly chosen according to a probability measure specified by a particular choice of the ensemble. The simplest yet nontrivial choice is to assume that all entries are mean-zero, independent, identically distributed (i.i.d.) Gaussian numbers. This defines the standard real Ginibre ensemble which we denote Gin 1 . Note that the question is equally interesting for matrices with entries that are complex rather than real, defining the complex Ginibre ensemble, which we denote Gin 2 . Note that for such an ensemble eigenvalues λ i are purely complex with probability one.
It is the latter ensemble for which the study of the eigenvalue condition numbers has been initiated two decades ago by Chalker and Mehlig [8,35]. More precisely, Chalker and Mehlig introduced a matrix of inner products O ij = (l * i l j )(r * j r i ), which they called 'eigenvector overlaps'. The diagonal elements of that matrix are simply the squared eigenvalue condition numbers. They further associated with the diagonal elements of the overlap matrix the following single-point correlation function: where the angular brackets stand for the expectation with respect to the probability measure associated with the complex Ginibre ensemble, and δ(z − λ k ) is the Dirac delta mass at the eigenvalue λ k , so that the empirical density of eigenvalues in the complex plane z reads ρ(z) = 1 is the mean spectral density around z [3]. It turned out that in the bulk of the spectrum of the complex Ginibre ensemble the magnitude of a typical diagonal overlap O ii grows linearly with the size of the matrix N , so one needs to consider a rescaled objectÕ 1 (z) = 1 N O 1 (z) to obtain a non-trivial limit. In their influential papers [8,35] Chalker and Mehlig used the 'formal' perturbation theory expansion to evaluate asymptotically, for N 1, both the diagonal overlap O 1 (z) and its more general off-diagonal counterpart The first mathematically rigorous verification of the Chalker and Mehlig result for the diagonal overlap has been done in [43]. Remarkably, the function O 1 (z) can be efficiently studied within the formalism of free probability [33] which recently allowed to extend the Chalker-Mehlig formulas to a broad class of invariant ensembles beyond the Gaussian case [3,37]. O 1 (z) is also known for a finite size of the complex Ginibre matrix [8,35] and products of small Ginibre matrices [7]. It has been recently shown that for complex Ginibre matrices the one-and two-point functions conditioned on an arbitrary number of eigenvalues are related to determinantal point processes [2]. Various features characterizing the rich properties of eigenvectors of nonnormal random matrices have been also studied in [4,9].
Here it is necessary to mention that the interest in statistical properties of the overlap matrix O kl and related objects extends much beyond the issues of eigenvalue stability under perturbation and is driven by numerous applications in theoretical and mathematical physics. In particular, non-orthogonality governs transient dynamics in complex systems [30,32,40] (see also [16,34]), analysis of spectral outliers in non-selfadjoint matrices [36], and, last but not least, the description of the Dyson Brownian motion for non-normal matrices [5,6,31]. Another steady source of interest in the statistics of eigenvector overlaps is due to its role in chaotic wave scattering. In that context O 1 (z) and O 2 (z 1 , z 2 ) have been studied for a few special models different from Ginibre (both theoretically [19,26,27] and very recently experimentally [10,11]) and in the associated models of random lasing [38,39]. In the scattering context all eigenvalues are necessarily complex, and the lasing threshold is associated with the eigenvalue with the smallest imaginary part. For that special eigenvalue even the distribution of the overlap O ii has been studied [39].
In fact, Chalker and Mehlig not only analyzed O 1 (z), but also put forward a conjecture on the tail for the density P (O ii ) of the distribution of diagonal overlaps O ii . Namely, based on the exactly solvable case of 2 × 2 matrices and numerical simulations for the complex Ginibre case they predicted that for large overlaps the density will exhibit a tail P ii , making all the positive integer moments beyond O 1 (z) divergent. This conjecture has been settled only recently with two different methods, by Bourgade and Dubach [5] (where some information about O l =k was also provided) and by Fyodorov [21]. The latter paper also revealed that for real eigenvalues of real Ginibre matrices Gin 1 the diagonal overlaps O ii are distributed with even heavier tail: ii , making even the mean of the overlap divergent. To address the above distributions it is convenient to introduce the following natural generalization of the equation (2): interpreted as the joint density function (JDF) of the 'diagonal' (or 'selfoverlap') non-orthogonality factor t = O ii − 1 for the right and left eigenvectors corresponding to eigenvalues in the vicinity of a point z = x + iy in the complex plane. The summation runs over all real eigenvalues for the real ensembles and over all eigenvalues for complex ensembles. As such, it is not a probability density, because it is normalized to the mean total number of (real) eigenvalues. As was shown in [5,13,21] the JDF P N (z, t) tends (after an appropriate rescaling of the variables z and t with the size N ) to the inverse gamma distribution as N 1: Here parameter β = 1 corresponds to the real eigenvalues of real Ginibre matrices (in which case the parameter z should be chosen to be real), β = 2 to the complex Ginibre case and β = 4 to the quaternionic Ginibre case.
For complex Ginibre the limiting expression (4) naturally incorporates the Chalker-Mehlig result. In the formula aboveÕ(z) = π −1 (1 − |z| 2 ), which is the large N limit of the rescaled one-point correlation function. Interestingly, despite the fact that for β = 1 the mean value O 1 (z) defined via (2) does not exist, the closely related rescaled combination,Õ 1 (z) = 1 2 √ 2π (1 − z 2 ), appears as a parameter in the inverse γ 1 distribution and therefore defines the typical value of the diagonal overlap. Further calculations in a few non-Gaussian rotationally invariant matrix ensembles (in particular, associated with 'truncations' of unitary matrices) done very recently in [12] suggest that (4) might exhibit a certain degree of universality. Note that the statistics of O ii for complex eigenvalues of real Ginibre matrices remains an outstanding problem, though it would be natural to expect that also in that case, for a fixed z with a non-vanishing imaginary part, the limit should be the same as for the complex Ginibre case.
Returning to the original question of eigenvalue condition numbers for real-valued matrices, the above results in particular imply that in contrast to well-conditioned eigenvalues of symmetric matrices with κ = 1 the typical condition numbers in fully asymmetric random matrices grow with matrix size as √ N [21] and show strong fluctuations. One of natural questions is then to ask how those properties evolve for matrices with a controlled degree of asymmetry in their entries. The aim of this work is to answer this question. To this end we consider matrices with i.i.d. real Gaussian entries, such that the entries X ij and X ji are correlated. The joint pdf for the elements of this ensemble, known in the literature either as the real partly symmetric Ginibre ensemble, or, alternatively, as the real Elliptic Gaussian ensemble, is given by Here dX = N i,j=1 dX ij is the flat Lebesgue measure over all matrix elements and the normalization constant reads (5) interpolates between the Real Ginibre Ensemble for τ = 0 and an ensemble of real symmetric matrices (Gaussian Orthogonal Ensemble) for τ = 1. In particularly, it is well known that for large sizes N 1 a nontrivial scaling regime of weak non-Hermiticity arises as long as the product N (1−τ ) is kept of the order of unity [15,[23][24][25]28]. It is this regime when non-normality gradually develops, and the condition numbers κ i start growing away from the value κ i = 1. Our considerations allow us to address this regime in a quantitative way.

Statement of the Main Results
It turns out that the method of evaluating the JDF in (3) suggested for the Ginibre case in [21] works for the Elliptic ensemble as well, though actual calculations turn out to be significantly more involved. Relegating technical detail to the rest of the paper, in this section we present our main findings.
Our main theorem gives the joint density function of the eigenvalue λ i and the shifted overlap t = O ii −1 for elliptic matrices of a given size N distributed according to (5). It takes a more compact form when the rescaled variable Theorem 2.1. Let X N be an N × N random matrix with the probability density function given by (5). Let us define the rescaled and shifted eigenvalue condition . The joint density (3) of a real eigenvalue z and the associated squared condition number expressed via the variable q is given by where the functions: as Remark 2.2. Note that for τ = 0 these quantities simplify to N (z) for the real Elliptic ensemble of even size N is known thanks to Forrester and Nagao [18]. It is given by ρ (r) For odd N the density can be obtained using the method from [17]. Our expression (6) by its very definition must reproduce the Forrester-Nagao result after integration over the variable t. Performing such an integration analytically is, however, a challenging task which we managed to complete for N = 2, 3, 4. Nevertheless, we checked that performing such an integral numerically for moderate values of N gives results that are indistinguishable from the density of real eigenvalues, see "Appendix A'.
As the expression (6) is exact, the joint density in the original variable, P N (z, t), can be further analyzed in interesting scaling limits as N → ∞. The first of such limits is the so-called bulk scaling corresponding to the eigenvalues inside the limiting support of the spectrum which (after appropriate rescaling z → √ Nz) for a fixed z and 0 ≤ τ < 1 represents an ellipse in the complex plane (hence the name for the ensemble), centered at the origin, with semi-axis 1 − τ along imaginary axis and 1 + τ along the real axis. Since we are dealing only with real eigenvalues, we restrict ourselves to real z such that |z| < 1 + τ , where the following asymptotics holds: This asymptotics shows that the typical value of the diagonal overlap t = O ii − 1 in this regime is always of the order N as N 1, similarly to the behavior for the Ginibre case τ = 0. Moreover, by recalling that the asymptotic density of real eigenvalues in the elliptic case is ρ , we see that (13) is exactly of the form (4) for β = 1.
When approaching the boundary |z| = 1+τ of the eigenvalue support the typical diagonal overlapÕ 1 (z) tends to zero, and in the appropriate scaling vicinity of the boundary it becomes parametrically weaker, as the variable t in such a regime becomes of the order √ N : Corollary 2.5. (edge scaling) Take a fixed 0 ≤ τ < 1 and parameterize z and

q) exists and is equal to
Note that this form is essentially the same as the one found for the real Ginibre case in [21].
Finally, the ultimate goal of our study is to investigate the weak non-Hermiticity regime occurring for N → ∞ when τ approaches unity with the rate O(N −1 ), so that the parameter 2N (1 − τ ) = a 2 is fixed. Such parameter therefore controls the deviation from the fully symmetric limit. In this regime of 'almost symmetric' matrices already a finite fraction of eigenvalues of the order of N are real, and asymptotically their mean density is given by [15,25] where ρ sc (z) = 1 2π √ 4 − z 2 is the standard Wigner semicircle density characterizing real symmetric GOE matrices, and A = (πρ sc (z)a) 2 . As anticipated, such a regime turns out to be not only 'weakly non-Hermitian', but also 'weakly non-normal' as the typical value of the diagonal overlap t = O ii − 1 turns out to be of the order of unity in the bulk of the spectrum, namely Corollary 2.6. Let |z| < 2 and t ≥ 0 be fixed. Consider the limit P weak (z, t) = lim where A = (πρ sc (z)a) 2 and ρ sc (z) = 1 2π √ 4 − z 2 . See also Fig. 1 for a plot of this distribution. Remark 2.7. After integration by parts one can rewrite the above as From this form it is easy to check that ∞ 0 P weak (z, t)dt agrees with the mean density (15), as expected.
We thus conclude that the characteristic tail P weak (z, t) ∼ t −2 is the most robust feature of the condition number density for real eigenvalues of asymmetric matrices, as it survives in the regime of arbitrary weak asymmetry as long as a > 0 (Fig. 1).

Derivation of the Main Results
We briefly outline an adaptation of the method of evaluating the JDF in (3) following [21] with necessary modifications.

Partial Schur Decomposition
Let λ be a real eigenvalue of an N × N real matrix X N . Then it is well known, see, e.g., [14], that the matrix X N can be decomposed as The Lebesgue measure on X N can be decomposed as dX where dS(v) denotes the uniform measure on the (N − 1)-hemisphere [14] (see also supplementary material of [22]). It turns out to be more technically convenient to concentrate on evaluating a characteristic function L(z, p) = δ(z − λ)e −pb T b N representing the Laplace transform of the JDF P(z, t). Hereafter by . . . N we denote the expected value with respect to the probability measure (5) for matrices X N of size N .
Proof. Substituting the decomposition (18) together with the associated decomposition of the Lebesgue measure into the probability measure of the elliptic ensemble (5) one can easily see that the ensemble average in (20) amounts to performing the following integral: The integral over v yields the volume of the (N − 1)-hemisphere dS(v) = π N/2 Γ( N 2 ) [14]. The integral over w is Gaussian and can be easily performed, giving the factor Taking all the multiplicative numerical constants into account and the factor | det(z − X N −1 )| from the Jacobian, we arrive at (20).

Ratio of Determinants
The problem has been therefore reduced to the calculation of the expectation for the ratio of two random determinants which is evaluated as where P N , R N and T N are defined in (8) (24) into (20) we see that L(z, p) is already represented as a Laplace transform, but of the rescaled variable t(1 − τ ), and (6) follows.
The proof of Theorem 3.2 proceeds via employing the supersymmetry approach to ratios of determinants.
Proof. Let χ, ρ, θ, η denote N -component vectors in anticommuting (Grassmann) variables. This allows us to rewrite the determinant in the numerator as a standard Berezin Gaussian integral The inverse square root of the determinant of a symmetric positive definite matrix A can be represented as a standard Gaussian integral. Namely, introducing N -component real vectors s 1 , s 2 we can write where we denoted u 2 = 2p (1 − τ 2 ). This provides a representation of the righthand side in (24) in the form The identity e −TrXA N = e 1 2 Tr(AA T +τA 2 ) allows us to perform the ensemble average. This in turn produces terms that are quartic in Grassmann variables, which we further bilinearize by employing a few auxiliary Gaussian integrals, the step known as the Hubbard-Stratonovich transformation: where we use the notation d 2 z = dx dy for z = x + iy.
Applying these transformations converts all integrations over anticommuting variables into a Gaussian Berezin integral which we can write as dχdρdθdηe − 1 2 ξ T Mξ , where ξ T = (χ T η T θ T ρ T ) and the antisymmetric matrix M is given by Here we denoted g = z + if √ τ , h = z + ic √ τ for brevity, and introduced the rank-two matrix A = s 1 s T 2 + τ s 2 s T 1 . The Berezin Gaussian integration yields a Pfaffian of the matrix M . The Pfaffian is evaluated using the fact that through the Schur decomposition A is equivalent toÃ ⊕ 0 N −2 , where 0 N −2 denotes the matrix of size N − 2 filled with zeros, andÃ As a consequence, Pf(M ) = Pf(F ) N −2 Pf(G), where the matrices of size 4 and 8 read explicitly The Pfaffian reads We notice that Therefore, the resulting expression depends on the vectors s 1 and s 2 only via their scalar products, so it is convenient to parameterize integrals by the entries of the associated Gram matrix [20, Section 2] [29, Theorem 1a in Appendix The Jacobian of this change of variables is (detQ) N −3 2 , while the integration over redundant angular variables yields the factor C (o) [29]. The range of integration is restricted by the non-negativity conditions Q 1 ≥ 0, Q 2 ≥ 0, detQ = Q 1 Q 2 − Q 2 ≥ 0. Following [21] it is convenient to change variables into r = (detQ) 1/2 and parameterize the integration region by Q 2 = r 2 +Q 2 Q1 . The change of measure reads dQ 1 dQ 2 dQ = 2 dQ1 Q1 rdrdQ. After rescaling Q 1 → uQ 1 , we have After the change of variables (34), Pf(M ) can be expressed as Vol. 22 (2021) Condition Numbers for Real Eigenvalues 321 The integrals over a, b, c, f are performed in the following way. Let us denote Expanding the expression in the bracket and using the binomial theorem twice, we obtain where He m (x) = (2π) −1/2 ∞ −∞ e − y 2 2 (x+iy) m dy are the monic Hermite polynomials. The internal sum can be performed via the Christoffel-Darboux formula, finally yielding Note that P N can be interpreted as the expectation of the squared characteristic polynomial det(z − X)(z − X T ) N and in this capacity has been already studied for the Real Elliptic ensemble [1]. All other integrals over a, b, c, f in (35) are performed in a similar way. After exploiting the three term recurrence for Hermite polynomials He N +1 (x) = xHe N (x) − N He N −1 (x), the integrals are evaluated to where R N and T N are defined by (9) and (10). Note also that R N (z) = . It is convenient to exploit the structure of (40) and exponent in (35) and rescale further Q → Q 1+τ and, similarly, Q 1 → Q1 1+τ . Recalling that u 2 = 2p(1 − τ 2 ), one then arrives at The integrations over Q and r are elementary. The integral over Q is Gaussian, while the one over r is of the type . The remaining integral over Q 1 formally looks logarithmically divergent. To see the cancellation of the divergent part, one should exploit a non-trivial identity which is verified in "Appendix B". After elementary but tedious algebraic manipulations with the help of Mathematica in rearranging terms, we finally obtain

Asymptotic Analysis
As P N , R N and T N are the building blocks of the determinant, we consider here their large-N asymptotics. First, we find convenient integral representations, which should allow the use of the Laplace method. For this we start from (39) and, using the integral representation for Hermite polynomials in (7) with both signs, we obtain The sum is evaluated using N k=0 An analogous procedure applied to T N gives Lemma 3.5.
For N 1 the integral over p can be most straightforwardly evaluated by the Laplace method, yielding that the leading contribution to P N (z √ N ) can be written as For large N the q-integral above can be performed (for a fixed, N -independent real value of z) by the standard saddle point method, with the saddle point position given by q = q * := − iz √ 2 1+τ yielding the required asymptotic formula: The same type of reasoning applied to (46) gives Finally, the asymptotics is obtained from (49) using the fact that R N (z) = 1 2(N +1) dPN+1(z) dz . Upon inserting this asymptotics into (6) and rescaling q → Nq it is clear that only the second to last term in the square bracket provides the leading order contribution, which happens to be As a consequence, the joint pdf reads Changing variables to t = q 1−τ , one immediately recovers Corollary 2.4.

Edge Scaling.
When z is tuned to values parametrically close to z = ±(1 + τ ) where the step-function argument in equations (49)-(52) is close to unity with a distance O(N −1/2 ), the corresponding asymptotics need to be evaluated with higher accuracy. Such a regime is known as the edge scaling, which features in Corollary 2.5 which we now prove.
Proof. In the proof we choose the vicinity of z = 1+τ . Correspondingly, in (47) we now scale z = 1 + τ + w √ N , where w is of order 1. The transition from (47) to (48) remains the same as before. Now we use the integral representation of the incomplete gamma function Γ (N, x) = x N ∞ 1 u N −1 e −ux du, which helps to rewrite the integral (48) in the form An inspection shows that whereas the q−integration is dominated by the contribution from the saddle point q = − √ 2i, the last u-integral is dominated by the vicinity of u = 1 of the width O(N −1/2 ). Parameterizing in such a vicinity u = 1 + v √ N one then arrives at the leading-order asymptotics After the change of variables u = 1+τ 1−τ (v + 2w 1+τ ) and the use of Stirling's approximation Γ (N + 1) ∼ √ 2πN N + 1 2 e −N , we obtain The last integral is related to the complementary error function erfc(x) = we obtain that in such a regime asymptotically R N ∼ √ NP N . From the asymptotics of (52) one expects that the leading order contributions from (N + 1)P N and T N cancel. Therefore, one needs to work with the appropriate integral representation. Combining (45) and (46) and following the analogous reasoning as above we obtain which is of the same order as √ NP N , as expected. To get the correct asymptotics at the edge, we rescale q → q √ N in (6). It is now clear that the first term in the square bracket in (6) is subleading and the contribution of other terms is of the same order. The asymptotics of and we obtain After denoting w = δ τ √ 1 − τ 2 and q = σ √ 1 − τ 2 , the statement of Corollary 2.5 follows.

Weak Non-Hermiticity.
Our final goal is to provide the proof of Corollary 2.6.
Proof. We start again from (47) keeping z fixed and N −independent like before in the bulk case, but for the weak non-Hermiticity regime we replace τ = 1 − a 2 2N . It is then immediately obvious that the p-integral is no longer dominated by the small vicinity p ∼ N −1/2 , but rather by the integration range of order unity. Then a quick inspection shows that for extracting the leading asymptotics in the large N limit one can effectively replace (47) by Performing the integral over q by the saddle point method, we see that the range of integration over p is given by |p| < . After a few simple changes of variables and straightforward manipulations one arrives at The asymptotic behavior of R N simply follows from the relation R N (z) = (62) Note that in (6) we used the rescaled quantity q = (1 − τ )t. Therefore, for the correct asymptotics, we need to rescale q → 2N a 2 t. This shows that all terms in the square bracket are of the same order. Direct use of the asymptotic forms (61) and (62) leads to (16).
Open Access. This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/ licenses/by/4.0/.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Appendix A: Density of Real Eigenvalues for Moderate Matrix Size
For N = 2 the joint pdf (6) reads The substitution t 2 = q q+1 allows one to calculate the integral. After integration by parts, we obtain  Figure 2. Histograms of the density of real eigenvalues for the Real Elliptic ensemble with τ = 0.9 obtained by the direct diagonalization of 10 6 matrices of size N = 3 (left) and 2 · 10 5 matrices of size N = 10 (right). Blue solid lines represent the formula obtained by analytical (left) and numerical (right) integration of P(z, q) (6). Formulas are rescaled so that the density is normalized to 1 which agrees with (11)-(12) when we substitute N = 2. This way, with the help of Mathematica software, we were also able to perform integration for N = 3, 4.
For N = 4 we again see agreement with the Forrester-Nagao result [18], while for N = 3 we compared the results of integration with the numerical diagonalizations of random matrices, see Fig. 2. For moderate matrix sizes, where the symbolic calculations were not possible, we numerically integrated (6) and compared with numerical diagonalization, observing good agreement, see

Appendix B: Proof of Identity (42)
We shall prove (42) by induction.
Proof. The first step is trivial as this identity can be verified by substituting Hermite polynomials for low N . Let us assume that (42) holds for N −1. Using formulas (8)- (10) it is easy to find the recurrence relations where, for simplicity, we omitted the argument z √ τ of Hermite polynomials. These recursions allow us to rewrite the left-hand side of (42) as The expression in square brackets is zero by the induction assumption. Verification that the second line equals 0 relies on the substitution of (68) and (69) and the consecutive use of the three term recursion He N +1 (x) = xHe N (x) − N He N −1 (x).