Condition numbers for real eigenvalues in real Elliptic Gaussian ensemble

We study the distribution of the eigenvalue condition numbers $\kappa_i=\sqrt{ (\bb{l}_i^* \bb{l}_i)(\bb{r}_i^* \bb{r}_i)}$ associated with real eigenvalues $\lambda_i$ of partially asymmetric $N\times N$ random matrices from the Gaussian elliptic ensemble. The large values of $\kappa_i$ signal about the non-orthogonality of (bi-orthogonal) set of left $\bb{l}_i$ and right $\bb{r}_i$ eigenvectors and enhanced sensitivity of the associated eigenvalues against perturbations of the matrix entries. We derive the general finite $N$ expression for the joint probability density (JPD) ${\cal P}_N(z,t)$ of $t_i=\kappa_i^2-1$ for $\lambda$ conditioned to have a value $z$ and investigate its several scaling regimes in the limit $N\to \infty$. When the degree of asymmetry is fixed as $N\to \infty$, the number of real eigenvalues is $O(\sqrt{N})$, and in the bulk of the real spectrum $t_i=O(N)$, while on approaching the spectral edges the non-orthogonality is weaker: $t_i=O(\sqrt{N})$. In both cases the corresponding JPDs, after appropriate rescaling, coincide with those found in the earlier studied case of fully asymmetric (Ginibre) matrices, see \cite{F2018}. A different regime of weak asymmetry arises when a finite fraction of $N$ eigenvalues remain real as $N\to \infty$. In such a regime eigenvectors are weakly non-orthogonal, $t=O(1)$, and we derive the associated JPD, finding that the characteristic tail ${\cal P}(z,t)\sim t^{-2}$ 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 nonnormal 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 spectra of such matrices numerically: keeping 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 always can be chosen 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 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 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 [39,40]. Note also that the Cauchy-Schwarz inequality implies κ ≥ 1, with κ = 1 only when X is normal.
It is natural to ask how well-conditioned is a 'typical' asymmetric matrix. 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 Gaussian numbers. This defines the standard real Ginibre ensemble which we denote Gin 1 . Note that the question is equally interesting for matrices whose entries are complex rather than real, defining complex Ginibre ensemble which we denote Gin 2 . Note that for such 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,33]. 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 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,33] 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 [41]. Remarkably, the function O 1 (z) can be efficiently studied within the formalism of the free probability [31] which recently allowed to extend the Chalker-Mehlig formulas to a broad class of invariant ensembles beyond the Gaussian case [3,35]. O 1 (z) is also known for finite size of the complex Ginibre matrix [8,33] 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 characterising rich properties of eigenvectors of nonnormal random matrices have been also studied in [4] and [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 [28,30,38], see also [15,32], analysis of spectral outliers in non-selfadjoint matrices [34], and, last but not least, the description of the Dyson Brownian motion for non-normal matrices [5,6,29]. 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 ) has been studied for a few special models different from Ginibre (both theoretically [18,24,25] and very recently experimentally [10,11]) and in the associated models of random lasing [36,37]. In the scattering context all eigenvalues are necessarily complex, and the lasing threshold is associated with an eigenvalue with the smallest imaginary part. For that special eigenvalue even the distribution of the overlap O ii has been studied [37].
In fact, already Chalker and Mehlig not only analysed 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 exactly solvable case of 2 × 2 matrices and numerical simulations for complex Ginibre case they predicted that for large overlaps the density will show 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 in [5] (where some information about O l =k was also provided) and by Fyodorov [20]. The latter paper also revealed that for real eigenvalues of a 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 (conditional) probability density function of the 'diagonal' (or 'self-overlap') 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. We will call it for brevity the joint probability density (JPD) of the two variables, t and z. As was shown in [5,20] the JPD P N (z, t) tends (after 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 the real Ginibre matrices (in which case the parameter z should be chosen real) and β = 2 to the complex Ginibre case. Recall that in the above the limiting spectral density of real eigenvalues for β = 1 is ρ(z) = 1 2 √ 2π for the interval |z| < 1, whereas the limiting spectral density of real eigenvalues for β = 2 is ρ(z) = π −1 inside the unit circle |z| < 1.
The limiting expression (4) naturally incorporates for complex Ginibre case 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 defined via O 1 (z) does not exist, its rescaled version,Õ 1 (z) = 1 2 √ 2π (1 − z 2 ), appears as a parameter in the inverse γ 1 distribution and defines therefore 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 nonvanishing 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 [20] 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 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 . The parameter τ ∈ [0, 1] tunes the degree of correlation, E(X ij X ji ) = τ for i = j, and (5) interpolates between the Real Ginibre Ensemble for τ = 0 and 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 [14,[21][22][23]26]. 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 JPD in eq. (3) suggested for Ginibre case in [20] 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 probability 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 more compact form when the rescaled variable q = (1 − τ )t is considered.
Theorem 2.1. Let X N be N × N random matrix with the probability density function given by (5). Let us define the rescaled and shifted eigenvalue condition number q = (1 − τ )(κ 2 − 1). The joint pdf (3) of real eigenvalue z and 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 P N = e z 2 Γ N + 1, , and the known result [20, eq. 2.5] is recovered.
Remark 2.3. The exact mean density of purely real eigenvalues ρ (r) N (z) for real Elliptic matrices of even size N is known due to Forrester and Nagao [17]. It is given by ρ (r) For N odd the density can be obtained using the method from [16]. 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. However, we checked that performing such integral numerically for moderate values of N gives indistinguishable results from the density of real eigenvalues, see Appendix A.
Being exact, the expression (6) 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 → √ N z) 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 behaviour for the Ginibre case τ = 0. Moreover, by recalling that the asymptotic density of real eigenvalues in 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 parametrize z and q as z exists and is equal to Note that this form is essentially the same as found for real Ginibre case in [20]. Finally, the ultimate goal of our study is to investigate the weak non-Hermiticity regime occuring 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 order of N eigenvalues are real, and their mean density asymptotically is given by [14,23] 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 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 ≥ 1 be fixed. Consider the limit P weak (z, t) = lim where A = (πρ sc (z)a) 2 and ρ sc (z) = 1 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,

Partial Schur decomposition
Let λ be a real eigenvalue of a N × N real matrix X N . Then it is well-known, see e.g. [13], that the matrix X N can be decomposed as where w is a column vector with N − 1 components and X N −1 is a matrix of size (N − 1) × (N − 1). The matrix O is known in the literature as the Hausholder reflection matrix. Note that although the left/right eigenvectors ofX N corresponding to λ are different from those of X N , the inner products (hence, the eigenvalue condition numbers) are the same. Parameterizing these eigenvectors as we immediately obtain for the associated condition number κ 2 Demanding thatl λ is the left eigenvector ofX N leads us to the relation b = (λ − X T N −1 ) −1 w. As a consequence [20] The Lebesgue measure on X N can be decomposed as dX N =C N | det(λ − X N −1 )|dλdwdX N −1 dO, with the known proportionality constantC N . 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 JPD 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 O yields the volume of the space of Hausholder transformations The integral over w is Gaussian and can be easily performed, giving the factor .
(22) 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)-(10).  (24) into (20) we see that L(z, p) is already represented as a Laplace transform 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 right-hand 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 auxilliary 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 where a ⊗ b T stands for the matrix with entries a i b j . The Berezin Gaussian integration yields Pfaffian of the matrix M , evaluating which explicitly gives We then see that the resulting integrand 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 [19] The Jacobian of this change of variables is (detQ) N −3 2 , while the integration over redundant angular variables yields the factor C (o) [27]. 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 [20] it convenient to change variables into r = (detQ) 1/2 , and parameterize the integration region by Q 2 = r 2 +Q 2 Q 1 . The change of measure reads dQ 1 dQ 2 dQ = 2 dQ 1 Q 1 rdrdQ. After rescaling Q 1 → uQ 1 , we have Pf(M). (33) Noticing that TrA = (1 + τ )Q, TrA 2 = (1 + τ 2 )Q 2 + 2τ Q 1 Q 2 , TrAA T = (1 + τ 2 )Q 1 Q 2 + 2τ Q 2 and det A = −τ det Q, the Pfaffian Pf(M ) can be expressed as 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 Cristoffel-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 elliptic ensemble [1]. All other integrals over a, b, c, f in (33) 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 (38) where R N and T N are defined by (9) and (10). Note also that R N (z) = 1 2(N +1) . It is convenient to exploit the structure of (38) and exponent in (33) and rescale further Q → Q 1+τ and, similarly, Q 1 → Q 1 1+τ . Recalling that u 2 = 2p(1 − τ 2 ), one then arrives at The remaining integration over Q is Gaussian, while the one over r is of the type . The integral over Q 1 formally looks like logarithmically divergent. To see the cancellation of the divergent part one should exploit a non-trivial identity (40) which is verified in the Appendix B. After further algebraic manipulations with the help of Mathematica 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 (37) and using the integral representation for Hermite polynomials in (7) we obtain The sum is evaluated using N k=0 An analogous procedure applied to T N gives Lemma 3.5.

Bulk scaling
Let us give the proof of the Corollary 2.4.
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 (46) 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 yielding the required asymptotic formula: The same type of reasoning applied to (44) gives Finally, the asymptotics is obtained from (47) using the fact that R N (z) = 1 2(N +1) dP N +1 (z) dz . Upon inserting this asymptotics into (6) and rescaling q → N q 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 (47)-(50) is close to unity by a distance O(N −1/2 ), the correponding asymptotics need to be evaluated with higher accuracy. Such regime is known as the edge scaling, which features in the Corollary 2.5 which we now prove.
Proof. In the proof we choose the vicinity of z = 1 + τ . Correspondingly, in (45) we now scale z = 1+τ + w √ N , where w is of order 1. The transition from (45) to (46) 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 helps to rewrite the integral (46) 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 widths O(N −1/2 ). Parametrizing in such a vicinity u = 1 + v √ N one then arrives at the leading -order asymptotic 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( we obtain that in such regime asymptotically R N ∼ √ N P N . From the asymptotics of (50) 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 (43) and (44) and following the analogous reasoning as above we obtain which is of the same order as √ N P N , as expected. To get the correct asymptotic at the edge, we rescale q → q √ N in (6). It is now clear that the first term in square bracket in (6) is subleading and 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 the Corollary 2.5 follows.

Weak nonhermiticity
Our final goal is to provide the proof of Corollary 2.6.
Proof. We start again from (45) keeping z fixed and N −independent like before in the bulk case, but for the weak nonhermiticity regime replace τ = 1 − a 2 2N . It is then immediately obvious that 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 (45) 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) = 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 square bracket are of the same order. Direct use of the asymptotic forms (59) and (60) leads to (16).

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 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 Forrester-Nagao result, 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 Fig. 2.  Figure 2: Histograms of the density of real eigenvalues for the elliptic ensemble with τ = 0.9 obtained by direct diagonalization of 10 6 matrices of size N = 3 (left) and 2 · 10 5 matrices of size N = 10 (right). Blue solid lines present 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.

Appendix B Proof of the identity (40)
We shall prove (40) 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 (40) 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 lhs of (40) as The the expression in square brackets is zero by the induction assumption. Verification that the second line equals 0 relies on the substitution of (66) and (67) and consecutive use of the three term recursion He N +1 (x) = xHe N (x) − N He N −1 (x).