On Statistics of Bi-Orthogonal Eigenvectors in Real and Complex Ginibre Ensembles: Combining Partial Schur Decomposition with Supersymmetry

We suggest a method of studying the joint probability density (JPD) of an eigenvalue and the associated ‘non-orthogonality overlap factor’ (also known as the ‘eigenvalue condition number’) of the left and right eigenvectors for non-selfadjoint Gaussian random matrices of size 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}. First we derive the general finite N expression for the JPD of a real eigenvalue λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\lambda}$$\end{document} and the associated non-orthogonality factor in the real Ginibre ensemble, and then analyze its ‘bulk’ and ‘edge’ scaling limits. The ensuing distribution is maximally heavy-tailed, so that all integer moments beyond normalization are divergent. A similar calculation for a complex eigenvalue z and the associated non-orthogonality factor in the complex Ginibre ensemble is presented as well and yields a distribution with the finite first moment. Its ‘bulk’ scaling limit yields a distribution whose first moment reproduces the well-known result of Chalker and Mehlig (Phys Rev Lett 81(16):3367–3370, 1998), and we provide the ‘edge’ scaling distribution for this case as well. Our method involves evaluating the ensemble average of products and ratios of integer and half-integer powers of characteristic polynomials for Ginibre matrices, which we perform in the framework of a supersymmetry approach. Our paper complements recent studies by Bourgade and Dubach (The distribution of overlaps between eigenvectors of Ginibre matrices, 2018. arXiv:1801.01219).


Introduction
Let x be a N -component column vector, real or complex. We will use x T = (x 1 , . . . , x N ) to denote the corresponding transposed row vector (and similar notation for matrices), and x * = (x 1 , . . . , x N ) for the Hermitian conjugate, with bar standing for complex conjugation. The inner product of two such vectors will be denoted as x * 1 x 2 = N i=1 x 1i x 2i . Let G be a N × N matrix which we assume to be non-selfadjoint and not normal: G * = G, G * G = GG * . We will further assume that all N eigenvalues λ a , a = 1, . . . , N of this matrix, which are in general complex numbers, have multiplicity one. Then the matrix is diagonalizable by a similarity transformation: G = S S −1 where = diag (λ 1 , . . . , λ N ) and S is in general non-unitary: S * = S −1 . The associated right eigenvectors defined by G x Ra = λ a x Ra are columns of the matrix S, whereas their left counterparts satisfying x * La G = λ a x * La form the rows of S −1 , and generically x Ra = x La . The sets x * La and x Ra of left and right eigenvectors can always be chosen to satisfy the bi-orthonormality condition x * La x Rb = δ ab for a, b = 1, . . . , N , but non-unitarity of S implies that x * Rb x Ra = δ ab and similarly x * La x Lb = δ ab . Then the simplest informative object characterizing the eigenvector non-orthogonality is the so-called 'overlap matrix' O ab = (x * La x Lb )(x * Rb x Ra ). In particular, the real diagonal entries O aa are known in the literature on numerical analysis as eigenvalue condition numbers and characterize sensitivity of eigenvalues λ a to perturbation of entries of G, see e.g. [48]. Namely, consider a family of matrices G(α) = G + αV , with V being an arbitrary matrix whose 2-norm is fixed as ||V || 2 = 1, whereas α is a real parameter controlling the magnitude of the perturbation. Denote, for a given V , the eigenvalues of G(α) as λ a (α) and considerλ a (α) = dλ a dα . A standard calculation using bi-orthonormality shows thaṫ λ a (0) = x * La V x Ra and therefore λ a (0) ≤ |x La | ||V || 2 |x Ra | = O 1/2 aa showing indeed that O aa controls the speed of change of eigenvalues under perturbation. As for some classes of non-normal matrices O aa 1, their eigenvalues could be much more sensitive to perturbations in comparison with their normal counterparts.
If the matrix G is random, it makes sense to be interested in statistics of O ab . This line of research originated from the influential papers by Chalker and Mehlig [7,37] who were the first to evaluate asymptotically, for large N 1, the lowest moments of the form where δ(z − λ a ) stands for the appropriate Dirac delta-distribution (so that e.g. the empirical density of eigenvalues at a (in general, complex) point z is given by N a=1 δ(z− λ a )). The brackets . . . Gin 2 denote here the expectation with respect to the probability measure on G known as the complex Ginibre ensemble, which we denote in this paper as Gin 2 to reflect that ensembles with complex entries are usually characterized by the Dyson index β = 2, see below. The probability measure on G with real entries known as the real Ginibre ensemble will be denoted correspondingly with Gin 1 .
For β = 2 Chalker and Mehlig were able to extract the leading asymptotic behaviour of O(z) and O(z 1 , z 2 ) in the N 1 limit. In particular, they found that O(z) ≈ N 2 (1−|z| 2 ) inside the 'Ginibre circle' characterized by the asymptotic mean eigenvalue density N a=1 δ(z − λ a ) Gin 2 ≈ N π 2 for |z| 2 < 1 and zero otherwise. This suggests that typically one should expect O aa ∼ N for eigenvalues inside the circle, which is parametrically larger than O aa = 1 typical for normal matrices.
In the last decades there was steady growth of interest in understanding properties of non-orthogonal random eigenvectors in theoretical physics, see [1,[4][5][6]13,22,27,31,33,38,41,45], with emphasis on calculating the Chalker-Mehlig correlators (1.1) and related objects beyond the framework of the complex Ginibre ensemble. One motivation comes from the abovementioned relevance of eigenvector correlations for describing the motion of complex eigenvalues under perturbations of the ensemble, see e.g. [39], and associated Dysonian dynamics, see e.g. [5] and Appendix A of [3]. Note that the non-orthogonality factors reflect non-normality of the matrix, which in the context of dynamical systems is known to give rise to a long transient behaviour, see a general discussion in [26,47]. In a related setting, non-symmetric matrices appear very naturally via linearization around an equilibrium in a complicated nonlinear dynamical system [24,36], and the non-orthogonality factors then control transients in a relaxation towards equilibrium [29]. Non-orthogonality also plays some role in analysis of spectral outliers in non-selfadjoint matrices, see e.g. [40] and references therein. Another strong motivation comes from the field of quantum chaotic scattering, where non-selfadjoint random matrices of special type (different from the Ginibre ensembles) play a prominent role, see e.g. [12,18,21,35,43] for the background information. The corresponding non-orthogonality overlap matrix O ab shows up in various scattering observables, such as e.g. decay laws [44], 'Petermann factors' describing excess noise in open laser resonators [45], as well as in sensitivity of the resonance widths to small perturbations [22,31]. Unfortunately, main progress in understanding properties of the bi-orthogonal eigenvectors for such ensembles relied on treating non-Hermiticity perturbatively in a small parameter, whereas non-perturbative results are scarce [13,38,45].
In the Mathematics community, a systematic rigorous research in this direction seems to have started only recently [49]. In a very recent development, Bourgade and Dubach [3] demonstrated a possibility to find the law of the random variable O aa for the complex Ginibre ensemble, asymptotically for large N , and provided a valuable information about the off-diagonal correlations between the two different eigenvectors at various scales of eigenvalue separation (the so-called 'microscopic' vs. 'mesoscopic' scales). That work motivated the present paper, where we use a rather different approach to consider the following object (1.2) interpreted as the (conditional) probability density of the 'diagonal' (or 'self-overlap') non-orthogonality factor t = O aa − 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. Note that this JPD is normalized in such a way that the integral R + P(t, z) dt = a δ(z − λ a ) := ρ N (z) coincides with the mean density ρ N (z) of eigenvalues around point z in the complex plane.
Naturally, the JPD function P(t, z) can be defined for a general random matrix G and, in particular, may be used to quantify the statistics of eigenvalue sensitivity parameters for such matrices. To give an example, consider again the family of matrices G(α) = G + αV , but choose the perturbation V to be a random matrix independent of G. For simplicity one may take V to be proportional to a random complex Ginibre matrix, and normalized in such a way that its entries V i j are i.i.d. mean zero complex numbers with the variance V i j V kl Gin 2 = 1 N δ ik δ jl . Then the eigenvalue sensitivity to such a perturbation is given byλ a (0) = x * La V x Ra and for a fixed G becomes a complex Gaussian variable with mean zero and variance |λ a (0)| 2 V = 1 N O aa . Define now the probability density π(w, z) of the eigenvalue sensitivity at a point z of the complex plane via π(w, z) = a δ w −λ a (0) δ(z − λ a ) G,V where the ensemble averaging goes both over G and over V . Since the complex Gaussian variable w =λ a (0) has the density π(w) = 1 π |λ a (0)| 2 V e −|w| 2 / |λ a (0)| 2 V with respect to the Lebesgue measure d(Im w)d(Re w) = 1 2 dwdw and recalling O aa = 1 + t we immediately see that relating the statistics of the eigenvalue sensitivity in that case to the knowledge of P(t, z).
In this paper we concentrate on finding explicit expressions for P(t, z) for Ginibre matrices, both real and complex. We first consider in Sect. 3 the case Gin 1 of real Ginibre matrices with β = 1. To this end it is useful to recall that real-valued matrices may have either purely real eigenvalues or pairs of complex conjugate eigenvalues. As the result, ρ N (z) for real Ginibre ensemble necessarily has the form ρ N (z) = ρ (c) where the non-singular part ρ (c) N (z) describes the mean density of complex eigenvalues, whereas ρ (r ) N (x) describes the mean density of purely real eigenvalues, so that b a ρ (r ) N (λ)dλ stands for the mean number of real eigenvalues in an interval [a, b] of the real axis. As a consequence, the introduced JPD P(t, z) inherits the same structure The organization of the paper is as follows. A summary of the main results and discussion of possible directions for the future work is presented in the Sect. 2. We start our consideration with demonstrating in Sect. 3 a way to evaluate P (r ) (t, λ), which describes non-orthogonality factor for eigenvectors associated with a real eigenvalue λ of the real Ginibre ensemble. First, we reduce the problem of finding P (r ) (t, λ) to a problem of evaluating certain ratios of determinants of random real-symmetric matrices with block structure, which as one may eventually see are intimately related to a deformed version of the so-called real chiral ensemble. Technical calculations within a framework of the supersymmetry approach which proves to be an efficient technical tool for dealing with such ratios of determinants are presented in Sect. 3.3. Our approach yields exact and explicit formula for any size N , which is then amenable to extracting the appropriate 'bulk' and 'edge' scaling limits as N → ∞. The problem of evaluating P (c) (t, z) for real Ginibre matrices remains presently outstanding, and we hope to be able to address it in a future publication.
In the next section, we apply essentially the same method for evaluating P (c) (t, z) in the complex Ginibre ensemble Gin 2 , i.e. β = 2. The computations and results become somewhat more technically involved, and considerably simplify only for the special case |z| = 0. General case is treated again by the supersymmetry approach outlined in Sect. 4.1.2. Eventually, we present an explicit finite-N expression for any z, and then extract the corresponding 'bulk' and 'edge' scaling limits.

Real
. Then the joint probability density P (r ) (t, λ) of the variables t = O λ − 1 and the eigenvalue λ is given for any N ≥ 2 by In what follows we will frequently omit the index λ in eigenvectors to lighten the notations, simply writing x T L or x R .
Remark 2.2. The above expression is a generalization of the exact mean density of purely real eigenvalues ρ (r ) N (λ) for real Ginibre matrices of size N explicit expression for which is known due to Edelman et al. [9], see also [10]: is the incomplete -function. The expression (2.3) can be most easily recovered from (2.2) as ρ (r ) N (λ) = R + P(t, λ) dt after rewriting (2.2) in an equivalent form as and then introducing u = |λ| t 1+t as the integration variable, cf. (3.31) below.
The expressions (2.2) or (2.5) show that the random 'overlap' variable t is maximally heavy-tailed, so that all positive integer moments E[t μ ], μ ≥ 1 for a fixed λ are divergent due to the 'fat tail' P(t, λ) ∼ t −2 as t → ∞, and only the normalization μ = 0 is finite and yields the density (2.3).
Being exact, the expression (2.2) can be further analyzed in interesting scaling limits as N → ∞. In fact, we find the form (2.5) most suitable for such an analysis. In particular, by rescaling λ = √ N x, t = N s (which is standard to call the bulk scaling limit), then considering x, s as fixed when N → ∞ and exploiting the appropriate asymptotic behaviour of the incomplete -function: and P bulk (x) being the limiting mean density of real eigenvalues within the bulk of the spectrum of the real Ginibre ensemble, which is known to be uniform inside its support.
Another natural edge scaling limit arises in the vicinity of the edge of the support of limiting spectral measure for real eigenvalues, that is for λ = √ N + δ, with δ < ∞ being fixed. It is easy to understand that the variable t needs to be rescaled in this regime as t = √ N σ , keeping σ fixed. A straightforward calculation using the well-known asymptotics (2.10) In particular, integrating the above over σ gives This expression is in full agreement with one for the limiting mean density of real eigenvalues at the edge of the spectrum of the real Ginibre ensemble, see [11].

Complex Ginibre ensemble.
For the case of complex Ginibre ensemble Gin 2 the corresponding joint probability density P (c) (t, z) of the non-orthogonality variable t = O z −1 and the associated complex eigenvalue z can be found in explicit form for finite N as well, but turns out to be given by a much more cumbersome expression in comparison with the real Ginibre case.

12)
Let z be a complex eigenvalue of G. Then the joint probability density P (c) (t, z) of the self-overlap non-orthogonality variable t = O z − 1 and the associated complex eigenvalue z for N ≥ 2 is given by are functions of |z| 2 explicitly defined via the relations to the incomplete -function as depend only on |z| 2 but not on the variable t.
Remark 2.4. The expression (2.13) is a generalization of the well-known mean density of eigenvalues of the complex Ginibre ensemble, see e.g. [34]: The latter can indeed be recovered as ρ (c) The joint density P (c) (t, z) at fixed z decays at large arguments t 1 as P (c) (t, z) ∼ t −3 as was already anticipated by Mehlig and Chalker on the basis of N = 2 example and informal eigenvalue repulsion arguments [37]. In contrast to the real Ginibre case such density does have the finite first moment. Only for the special value z = 0 the above joint density significantly simplifies and is given by Despite the relative complexity of (2.13), its bulk rescaling limit z = √ N w, t = N s, with w, s being fixed when N → ∞, can be straightforwardly extracted. To this end, it is convenient to use the following integral representations for d 1 are easily amenable to the standard asymptotic analysis by the Laplace method. In this way we find in the 'bulk' scaling limit for |w| 2 < 1 the following w-independent asymptotic behaviour: and then via (2.14) and (2.15) we further find We then see that in the 'bulk' limit the first term D (N ) 1 in the brackets of (2.13) is dominant in comparison with the other two, and taking the limit and zero otherwise. This expression agrees with results obtained by Bourgade and Dubach [3] in a different approach to the problem. Its first moment is precisely 1 π (1−|w| 2 ) inside the bulk of the spectrum, in agreement with the expression by Chalker and Mehlig. Finally, one also can extract the corresponding edge asymptotics by replacing |z| = √ N +δ and t = √ N σ and performing the limit N → ∞. With a help of the Mathematica package 1 one then finds that and we denoted = 1 − 2σ δ. In particular, one can check that integrating over σ yields a well-known formula for the mean edge density of complex eigenvalues: One also can see that the 'bulk' (2.24) and 'edge' (2.25) asymptotics match by replacing in the latter δ = 1 2 √ N (|w| 2 − 1) and σ = √ N s and letting N → ∞ for fixed |w| < 1 and s, checking that bulk (s, w).

Discussion of the method and open problems.
Our approach consists of two steps. In the first step we show that the partial Schur decomposition of Ginibre matrices employed in works [9,10] allows one to represent the JPD's P (r ) (t, λ) and P (c) (t, z), Laplacetransformed with respect to the variable t, in terms of the following object: where L = 0, 1, 2, . . . is an integer, p ≥ 0, the parameter β = 2 stands for the complex Ginibre ensemble and β = 1 for the real Ginibre one (in the latter case z = λ is real), and I N standing for the N × N identity matrix. In fact, the goals of the present paper require evaluation of (2.26) only for L = 2, but it is interesting to consider a more general problem, see below.
Note that for β = 2 we deal here with expectation values involving integer powers of characteristic polynomials for non-selfadjoint matrices G in both numerator and denominator. Studying similar objects for self-adjoint random matrices has a long history, see e.g. [2,19] for a background discussion and further references. At the same time, for β = 1 a half-integer power in the denominator is involved. To deal with the latter challenge we employ one of very few techniques available in that case, the so-called supersymmetry approach, see [32,51] for concise introductions and also [23,50] for earlier computations involving half-integer powers of characteristic polynomials for real symmetric Gaussian random matrices. We find it convenient to use a (rigorous) variant of the approach proposed originally in [14] and the final expression for D (2) N ,1 (λ, p) is given in (3.11) or (3.27). As a by-product of the same calculation one also finds for L = 0: The same procedure works, with due modifications, for the complex case β = 2, with the computational challenge now coming not from the half-integer power in the denominator, but from the higher integer power of the determinant in the numerator. The actual calculation is very straightforward for L = 0, becomes slightly more involved for L = 1, and for the case of actual interest L = 2 produces much more cumbersome expressions, see Sect. 4.1.2 for the derivation. In the end we have to resort to symbolic computer manipulations to deal with the ensuing integrals. Here we simply quote the results for L = 0 and L = 1 for the sake of completeness: In particular, by a direct integration one can check that D Given the complexity of arising expressions for β = 2, it is worth to give a different perspective on the problem. To that end we note that by introducing the matrices W = z − G our main object for β = 2 case, namely D (L) N ,2 (z, p), can be formally rewritten as with integration going over complex N × N matrices W with the weight function P L ,z (W, W * ) depending on an integer parameter L = 0, 1, 2, . . ., and on the complex parameter z: The right-hand side of (2.30) can be obviously interpreted as the mean inverse characteristic polynomial of the matrix W * W averaged over this 'ensemble' 3 closely related (though not identical for L > 0) to a limiting case of versions of the chiral ensemble with a 'source' considered in [8,46]. Namely, let us consider a more general version of (2.31): (2.32) where the 'source' matrix A is a fixed N × N complex matrix with the singular values (i.e. eigenvalues of A * A) being (in general, distinct) non-negative real numbers a 1 , a 2 , . . . , a N . Obviously our previous choice corresponded to all a i equal to |z| 2 . Note that the n-point correlation functions of eigenvalue densities for such type of a chiral ensemble (with a Hermitean source A) were derived in [46], but their knowledge is not sufficient for our purposes. Some information for the mean inverse characteristic polynomial for the chiral ensemble with a 'source' similar to (2.32) was given in the framework of the method of multiple orthogonal polynomials in [8]. In a separate paper [25] we are providing the full analysis of the problem for β = 2 for any integer positive L and N by deriving the following representation (see Proposition 3.9 in [25]): where we defined the following function of ρ = |z| 2 and τ = t/(1 + t): are Laguerre polynomials. The equivalence with (2.28)-(2.29) for L = 0, 1 can be straightforwardly verified (see the Appendix A of [25]). One can further perform the asymptotic analysis of (2.33)-(2.34) for N 1 and extract the bulk scaling asymptotics relevant for the present paper in a more transparent and systematic way than is provided by the supersymmetry approach in β = 2 case. Nevertheless, supersymmetric treatment has its own merits: the method is robust and is expected to be generalizable to more general ensembles of non-selfadjoint random matrices lacking the full invariance of the Ginibre ensembles.
The approach suggested in the present paper can be certainly adjusted for addressing overlaps of left/right eigenvectors corresponding to complex eigenvalues of real Ginibre ensemble, although in this way one encounters a few challenging technical problem not yet fully resolved. One can also envisage extensions addressing overlaps of two different eigenvectors, as well as posing similar questions for other types of non-Hermitian matrices, including those with quaternion structure for β = 4, those relevant in the theory of chaotic scattering and those relevant in the Quantum Chromodynamics context [42]. We hope to be able to answer some of these questions in future publications.

Eigenvectors of real matrices via partial Schur decomposition.
Let λ be a real eigenvalue of a matrix G (N ) with real entries, and denote the associated real left and right eigenvectors as x T L and x R . Then, as is well-known, see e.g. [9,10], it is always possible to represent the matrix G (N ) as for some real (N − 1)-component column vector w and a real matrix G (N −1) of size (N −1)×(N −1), whereas the matrix P is real symmetric and orthogonal: P T = P and P 2 = I N . As the left/right eigenvectors of the matrixG (N ) = λ w T 0 G (N −1) are given in terms of x T L and x R byx T L = x T L P andx R = Px R , we see that the corresponding selfoverlaps remain invariant: O λ =Õ λ . Hence one can usex T L andx R for our calculation. Moreover, it is immediately clear from the form ofG (N ) thatx R = e 1 : = (1, 0, . . . , 0) reducing the study of the self-overlap to the statistics of the norm b T λ b λ . To this end we find that the vector b T λ can be readily expressed in terms of the resolvent ofG (N −1) as

Partial Schur decomposition of the Real Ginibre Ensemble and overlap statistics.
In this section we show how to reduce the calculation of the Laplace transform L λ ( p) := ∞ 0 e − pt P (r ) (t, λ) dt of the JPD P (r ) (t, λ) defined in (1.2) to evaluating the ensemble average for the ratio of certain determinants, see (3.8-3.9).
For the ensemble of N × N square matrices G (N ) ≡ G ∈ M N (R) with independent identically distributed real matrix element G j,k ∼ N (0, 1) we use the notation Gin 1 and denote by the angular brackets · · · Gin 1 the expectation of any function F: R N ×N → C with respect to the associated probability distribution, though we will frequently omit the corresponding subscript to lighten the formulas.
Remark 3.1. In a similar way one defines the complex Ginibre ensemble Gin 2 G j,k = g (1) j,k + i · g (2) j,k , with i.i.d. g (·) j,k ∼ N (0, 1/2) , as well as the so-called quaternion Gin 4 ensemble which is however not considered in the present work.
Assigning the Dyson's index β = 1, 2, 4 4 one can write for all three ensembles the Joint Probability Density (JPD) with respect to the flat Lebesgue measure in the form Tr GG * . (3.4) where G * = G T stands for the Hermitian conjugation, and the bar for the complex conjugation. In this section we will concentrate in detail on the real case β = 1; similar treatment of β = 2 case will be briefly described in the last section. We start with exploiting the following Proposition 3.2 (see Lemma 3.1, Lemma 3.2 and Theorem 3.1 in [9]). Assume that the JPD of G is given by (3.4) with β = 1, and apply the transformation (3.1). Then the joint probability density of elements of the matrixG is given by with some normalization constant C 1,N .
The above j.p.d. can be used to calculate the Laplace transform of the probability density for our main object of interest, the random variable b T λ b λ which for a given value of λ is given by (3.3), or equivalently the characteristic function As the integral over w is Gaussian and p > 0 it can be readily performed yielding the factor where we have used det λ . Combining all the factors we finally see that the characteristic function in question is proportional to the ensemble average of the ratio of determinants, cf. (2.26) for β = 1, which we also may present in an equivalent, but different form convenient for further evaluation: 1 (λ, p), (3.8) where as will be found below (see the Corollary 3.4) N = 2 N /2 (N /2) and Here the ensemble average is performed over the j.p.d. (3.4) of real Ginibre matrices G of the reduced size (N − 1) × (N − 1). The problem of averaging the ratio of determinants in the above expression can be efficiently solved in the framework of the supersymmetry approach. The main steps of the corresponding procedure are presented in the following section. Interestingly, when implementing such an approach inverting the Laplace transform comes as a part of the procedure. In this way one recovers first (2.5) which by straightforward algebraic manipulations can be shown to be equivalent to (2.2).

Supersymmetry approach to the ratio of determinants and proof of Theorem 2.1.
In this section we evaluate that ensemble average for real Ginibre matrices G of size N × N , with the main object of interest being Our goal is to verify the following with the constant C N = Proof of the Theorem 2.1. The Prop. 3.3 when combined with (3.8) immediately provides the proof of (2.5), hence of the Theorem 2.1. Namely, to arrive at (2.5) one replaces N → N − 1 in (3.11), and substitutes it into (3.8). Noting that the result assumes the form of a Laplace transform in variable t makes its inversion trivial, and we recover the JPD P (r ) (t, λ) of the random variables t = b T λ b λ and the real eigenvalue λ as is given in (2.5).

Proof of the Proposition 3.3:
Proof. Let 1 , 2 , 1 , 2 be four column vectors with N anticommuting components each. Using the standard rules of Berezin integration one represents the numerator in the ratio (3.10) as a Gaussian integral Now we further use a form of the standard Gaussian integral well-defined for any realsymmetric matrix A and any positive > 0: , (3.13) where the integration goes over the vector S with N real commuting components. This allows to represent the denominator in (3.10) as a Gaussian integral over two such vectors S 1,2 : where ∝ here and below stands for (temporally) ignored multiplicative constants (in general, N -dependent) whose product will be restored in the very end of the procedure. After substituting the above representations to (3.10) and rearranging in the exponent as After the ensemble average is performed, there exists only one term in the exponential in the integrand which is quartic in anticommuting variables, and it is of the form The corresponding exponential factor is then represented as: where the formula above represents the simplest instance of what is generally known as the Hubbard-Stratonovich transformation. After such a representation is employed, it allows to perform the (by now, Gaussian) integration over the anticommuting vectors explicitly, and reduce the whole expression to the integral over the two vectors S 1,2 and over a single complex variable q: A straightforward calculation shows that the determinant in the above expression is equal to and we see that the integration over q, q is now easy to perform via using the polar coordinates: (N + 1, λ 2 ). (3.18) As to the remaining integrations, one may notice that the integrand depends only on the entries of a positive semidefinite real symmetric matrix A useful trick suggested in [14] in such a situation is to pass from the pair of vectors (S 1 , S 2 ) to the matrixQ as a new integration variable. Such change is non-singular for N ≥ 2 and incurs a Jacobian factor proportional to detQ (N −3)/2 (see the Appendix D of [16]). This finally brings D (2) N ,1 (λ, p) to the form (3.20) The next step requires employing a convenient parametrization of the integration domain defined via the inequalities is a real symmetric positive semidefinite matrix.
First, it is easy to see that such domain can be parametrized by expressing the diagonal entries Q 1 and Q 2 in terms of two real coordinates r ≥ 0, −∞ < θ < ∞ chosen in such a way that r = detQ 1/2 , θ = 1 2 ln (Q 1 /Q 2 ). By explicitly writing Q 1 = e θ r 2 + Q 2 , Q 2 = e −θ r 2 + Q 2 and evaluating the associated Jacobian we get in those coordinates dQ := d Q 1 d Q 2 d Q = 2 r dr d Q dθ . Although calculation in that parametrization is already quite convenient, it turns out that it becomes even shorter if one parametrizes the same domain in a related, but slightly less obvious way using instead the matrix entries Q 1 ≥ 0 and −∞ < Q < ∞ as new coordinates, complemented with r = detQ 1/2 ≥ 0, and expressing the remaining entry as Q 2 = r 2 +Q 2 Q 1 ≥ 0. This finally giveŝ where all integrals are well-defined and convergent for p > 0; in particular, the latter one can be evaluated explicitly in terms of the Bessel function of second kind as ( see [28], p. 363) In principle, one can demonstrate existence of a chain of integral identities which allows to perform the remaining integrations in (3.22) explicitly without changing the order of integrations. This way leads however to quite cumbersome intermediate formulas, and we proceed instead by changing the order in (3.22) (which can be justified by Fubini's theorem) to which allows to perform the integrals over Q and r much more efficiently. Namely, introduce the function Then it is easy to see that after renaming Q 1 → t the equation (3.23) can be rewritten as Now by using the relation (N + 1, λ 2 ) = e −λ 2 λ 2N + N (N , λ 2 ) one can see that for some real constant C N . Simple manipulations with incomplete -function (2.4) show that this is equivalent to (3.11). To establish the value for the constant one can use, for example, the limit p → ∞ where according to the definition (3.10) On the other hand, performing p → ∞ limit in (3.27) gives after a simple calculation lim p→∞ The definition of the left-hand side implies that the coefficient in front of λ 2N must be equal to unity, giving finally Proof. To establish the value of the constant N we consider the limit p → 0 in both sides of (3.8). By the very definition of the Laplace transform L λ ( p) its value at p = 0 must be equal to the mean density of real eigenvalues for the real Ginibre ensemble given in (2.3). On the other hand, for p = 0 the integration over t in (3.11) can be easily performed introducing u = |λ| t 1+t as new integration variable. One gets in this way: To get D (2) N −1,1 (λ, 0) featuring in the right-hand side of (3.8) replace in the above N → N − 1 and use (N ) Multiplying with e −λ 2 /2 and comparing with the left-hand side gives the value for the constant N .

Proof of Theorem 2.3
Our approach to complex Ginibre matrices follows essentially the same steps as for the real case, with very little modifications, and we only briefly indicate necessary changes. Similarly to (3.1), suppose that a complex-valued matrix G (N ) has only nondegenerate eigenvalues, and assuming it has an eigenvalue z (in general, complex) it can be represented as (see e.g. Sec. 6 of [20]) Now we exploit the analogue of Prop. (3.2) Proposition 4.1 (see Appendix B of [20]). Assume that the j.p.d. of G is given by (3.4) with β = 2, and apply the transformation (4.1). Then the joint probability density of elements of the matrixG is given by with some normalization constant C 2,N .
Using the above j.p.d. to calculate the Laplace transform of the probability density for the random variable b * z b z for a fixed value of z we arrive after standard manipulations at representing it as the expectation of the ratios of the determinants of the form where we introduced the notation, cf. (2.26) for β = 2, with averaging performed over the j.p.d. (3.4) of complex Ginibre matrices G of the size N × N . Note, that the value of the constant normalization factor in (4.4) is found a posteriori exactly in the same way as in the real case, by comparing the known expression for the mean density of complex eigenvalues (2.18) (coinciding with L z (0)) and the corresponding limit in the right-hand side of (4.4).
In the general case evaluating D (L) N ,2 (z, p) for integer L can be done essentially by the same supersymmetry method which was used in Sect. 3.3, with obvious necessary modifications imposed by symmetries. In particular, presence of higher powers of the determinants in the numerator of (4.5) necessitates to use L sets of anticommuting vectors for their representation, making the resulting integral representation in our version of the supersymmetry method significantly more cumbersome than in the real case. In the most relevant case L = 2 and a special choice of the spectral parameter z = 0 the expected value of the ratio featuring in the right-hand side of (4.5) can be relatively easily extracted as a special limiting case of a more general object evaluated in [17] or, in a different way, in [15]. The corresponding calculation is sketched in the first part of the next section. The supersymmetry approach for z = 0 works along exactly the same general lines as in the real case, but is somewhat more involved technically. The corresponding calculation is outlined in the second part of the next section. (4.5) for L = 2 . 4.1.1. Evaluation of (4.5) for L = 2 and |z| = 0. In the special case |z| = 0 an integral representation of the averaged ratio of determinants featuring in (4.4) which we find most convenient for our purposes was derived in [15], see Eq. (29) there. Actually, our object arises as a particular case n f = 2, n b = 1 of that formula, identifying x b = 2 √ p and considering a special limit X f = 0 [the latter limit is highly degenerate, and it is easier to perform it directly in eq. (25), and then rederive (29)]. In this way we arrive at representing D

Evaluation of
N ,2 (z = 0, p) as Now the integrals over R 1 and R 2 can be readily evaluated, with the result being simply where we introduced the notation Further employing a well-known integral representation for the Bessel function of the second kind shows that We finally conclude that  (4.5) for L = 2 and |z| = 0 by supersymmetry approach. Our goal is to verify the following Proposition 4.2.

Evaluation of
where the entering quantities were defined in equations (2.16)-(2.15).
The proof is very similar to the real case, and is outlined below.
Proof. One starts with using two copies of the set of four anticommuting vectors, namely A1 , A2 , A1 , A2 and B1 B2 , B1 , B2 , to represent separately two determinants in the numerator via Gaussian integrals, see (3.12). At the same time, one needs two commuting vectors S 1 , S 2 with N complex-valued components to represent the denominator: The ensemble averaging is performed by exploiting β = 2 analogue of (3.  T  A1 A1   T  A1 B1  T  B1 A1   T  B1 B1   T  A2 A2   T  A2 B2  T  B2 A2   T withQ F ,Q * F being a pair of general 2×2 complex conjugate matrices. This trick allows to integrate out the vectors with anticommuting component completely. The analogue of (3.17) takes the form × det Q F ⊗ I N i zI N + S 1 ⊗ S * 2 ⊗ I 2 i z I N + S 2 ⊗ S * 1 ⊗ I 2Q * F ⊗ I N (4.13) At the next step one can simplify the above expression by employing the singular value )U * with unitary U, V and R F1 , R F2 ≥ 0 and replacing the integration over complex vectors S 1,2 with one over the Hermitian positive semidefinite matrix [cf. (3.19)]: After straightforward algebraic manipulations this allows to represent (4.13) for N ≥ 2 as The Hermitian matrixQ ≥ 0 can be parametrized very similarly to (3.21). Namely, writing for the complex variable Q = ρ e iφ and using r = detQ 1/2 ≥ 0 together with Q 1 ≥ 0 as the coordinates, the domain of integration is parametrized by matriceŝ Q = Q 1 ρ e iφ ρ e −iφ r 2 +ρ 2 Q 1 , dQ = 2 d Q 1 Q 1 r dr ρ dρdφ; 0 ≤ φ ≤ 2π, 0 ≤ Q 1 , ρ, r < ∞.
Acknowledgements. The author is most grateful to Paul Bourgade and Guillaume Dubach for generously communicating their unpublished results at an early stage which stimulated his own research on the topic. Ramis Movassagh is acknowledged for an interesting discussion and bringing reference [39] to the author's attention, Gernot Akemann for pointing out [46] and Peter Forrester for mentioning [8]. Jacek Grela and Eugene Strahov are acknowledged for their collaboration on the associated analysis of Eq. (2.30) using different methods [25]. The present paper was started when preparing a lecture course for PCMI Summer School 2017, and essentially completed during the Beg Rohu Summer School 2017. The author would like to thank the organizers and participants of the schools for creating a stimulating atmosphere, and for the financial support of his participation in the events, in particular from the NSF Grant DMS:1441467. The research at King's College London was supported by EPSRC Grant EP/N009436/1 "The many faces of random characteristic polynomials".
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.