On statistics of bi-orthogonal eigenvectors in real and complex Ginibre ensembles: combining partial Schur decomposition with supersymmetry

Introduction

Let $\mathbf{x}$ be a $N-$ component column vector, real or complex. We will use $\mathbf{x}^{T}=(x_{1},\ldots,x_{N})$ to denote the corresponding transposed row vector ( and similar notation for matrices), and $\mathbf{x}^{*}=(\overline{x}_{1},\ldots,\overline{x}_{N})$ for the Hermitian conjugate, with bar standing for complex conjugation. The inner product of two such vectors will be denoted as $\mathbf{x}_{1}^{*}\mathbf{x}_{2}=\sum_{i=1}^{N}\overline{x}_{1i}x_{2i}$.

Let $G$ be a $N\times N$ matrix which we assume to be non-selfadjoint and not normal : $G^{*}\neq G,\,G^{*}G\neq GG^{*}$. We will further assume that all $N$ eigenvalues $\lambda_{a},\,a=1,\ldots,N$ of this matrix, which are in general complex numbers, have multiplicity one. Then the matrix is diagonalizable by a similarity transformation: $G=S\Lambda S^{-1}$ where $\Lambda=\mbox{diag}\left(\lambda_{1},\ldots,\lambda_{N}\right)$ and $S$ is in general non-unitary: $S^{*}\neq S^{-1}$. The associated right eigenvectors defined by $G\,\mathbf{x}_{Ra}=\lambda_{a}\mathbf{x}_{Ra}$ are columns of the matrix $S$, whereas their left counterparts satisfying $\mathbf{x}_{La}^{*}G=\lambda_{a}\mathbf{x}^{*}_{La}$ form the rows of $S^{-1}$, and generically $\mathbf{x}_{Ra}\neq\mathbf{x}_{La}$. The sets $\mathbf{x}^{*}_{La}$ and $\mathbf{x}_{Ra}$ of left and right eigenvectors can always be chosen to satisfy the bi-orthonormality condition $\mathbf{x}^{*}_{La}\mathbf{x}_{Rb}=\delta_{ab}$ for $a,b=1,\ldots,N$, but non-unitarity of $S$ implies that $\mathbf{x}^{*}_{Rb}\mathbf{x}_{Ra}\neq\delta_{ab}$ and similarly $\mathbf{x}^{*}_{La}\mathbf{x}_{Lb}\neq\delta_{ab}$. Then the simplest informative object characterizing the eigenvector non-orthogonality is the so-called ’overlap matrix’ $\mathcal{O}_{ab}=(\mathbf{x}^{*}_{La}\mathbf{x}_{Lb})(\mathbf{x}^{*}_{Rb}\mathbf{x}_{Ra})$. In particular, the real diagonal entries $\mathcal{O}_{aa}$ are known in the literature on numerical analysis as eigenvalue condition numbers and characterize sensitivity of eigenvalues $\lambda_{a}$ to perturbation of entries of $G$, see e.g. . Namely, consider a family of matrices $G(\alpha)=G+\alpha V$, with $V$ being an arbitrary matrix whose $2-$norm is fixed as $||V||_{2}=1$, whereas $\alpha$ is a real parameter controlling the magnitude of the perturbation. Denote, for a given $V$, the eigenvalues of $G(\alpha)$ as $\lambda_{a}(\alpha)$ and consider $\dot{\lambda_{a}}(\alpha)=\frac{d\lambda_{a}}{d\alpha}$. A standard calculation using bi-orthonormality shows that $\dot{\lambda}_{a}(0)=\mathbf{x}^{*}_{La}V\mathbf{x}_{Ra}$ and therefore $\left|\dot{\lambda}_{a}(0)\right|\leq|\mathbf{x}_{La}|\,||V||_{2}\,|\mathbf{x}_{Ra}|=\mathcal{O}_{aa}^{1/2}$ showing indeed that $\mathcal{O}_{aa}$ controls the speed of change of eigenvalues under perturbation. As for some classes of non-normal matrices $\mathcal{O}_{aa}\gg 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 $\mathcal{O}_{ab}$. This line of research originated from the influential papers by Chalker and Mehlig who were the first to evaluate asymptotically, for large $N\gg 1$, the lowest moments of the form

where $\delta(z-\lambda_{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 $\sum_{a=1}^{N}\delta(z-\lambda_{a})$). The brackets $\left\langle\ldots\right\rangle_{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 $\beta=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 $\beta=2$ Chalker and Mehlig were able to extract the leading asymptotic behaviour of $\mathcal{O}(z)$ and $\mathcal{O}(z_{1},z_{2})$ in the $N\gg 1$ limit. In particular, they found that $\mathcal{O}(z)\approx N^{2}(1-|z|^{2})$ inside the’ Ginibre circle’ characterized by the asymptotic mean eigenvalue density $\left\langle\sum_{a=1}^{N}\delta(z-\lambda_{a})\right\rangle_{Gin_{2}}\approx\frac{N}{\pi^{2}}$ for $|z|^{2}<1$ and zero otherwise. This suggests that typically one should expect $\mathcal{O}_{aa}\sim N$ for eigenvalues inside the circle, which is parametrically larger than $\mathcal{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 , 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. , and associated Dysonian dynamics, see e.g. and Appendix A of . 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 . In a related setting non-symmetric matrices appear very naturally via linearization around an equilibrium in a complicated nonlinear dynamical system , and the non-orthogonality factors then control transients in a relaxation towards equilibrium . Non-orthogonality also plays some role in analysis of spectral outliers in non-selfadjoint matrices, see e.g. 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. for the background information. The corresponding non-orthogonality overlap matrix $\mathcal{O}_{ab}$ shows up in various scattering observables, such as e.g. decay laws , ’Petermann factors’ describing excess noise in open laser resonators , as well as in sensitivity of the resonance widths to small perturbations . 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 .

In the Mathematics community a systematic rigorous research in this direction seems to have started only recently . In a very recent development Bourgade and Dubach demonstrated a possibility to find the law of the random variable $\mathcal{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

Naturally, the JPD function $\mathcal{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(\alpha)=G+\alpha 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_{ij}$ are i.i.d. mean zero complex numbers with the variance $\left\langle\overline{V}_{ij}V_{kl}\right\rangle_{Gin_{2}}=\frac{1}{N}\delta_{ik}\delta_{jl}$. Then the eigenvalue sensitivity to such a perturbation is given by $\dot{\lambda}_{a}(0)=\mathbf{x}^{*}_{La}V\mathbf{x}_{Ra}$ and for a fixed $G$ becomes a complex Gaussian variable with mean zero and variance $\left\langle|\dot{\lambda}_{a}(0)|^{2}\right\rangle_{V}=\frac{1}{N}\mathcal{O}_{aa}$. Define now the probability density $\pi(w,z)$ of the eigenvalue sensitivity at a point $z$ of the complex plane via $\pi(w,z)=\left\langle\sum_{a}\delta\left(w-\dot{\lambda}_{a}(0)\right)\delta(z-\lambda_{a})\right\rangle_{G,V}$ where the ensemble averaging goes both over $G$ and over $V$. Since the complex Gaussian variable $w=\dot{\lambda}_{a}(0)$ has the density $\pi(w)=\frac{1}{\pi\left\langle|\dot{\lambda}_{a}(0)|^{2}\right\rangle_{V}}e^{-|w|^{2}/\left\langle|\dot{\lambda}_{a}(0)|^{2}\right\rangle_{V}}$ with respect to the Lebesgue measure $d(Im\,w)d(Re\,w)=\frac{1}{2}dwd\overline{w}$ and recalling $\mathcal{O}_{aa}=1+t$ we immediately see that

relating the statistics of the eigenvalue sensitivity in that case to the knowledge of $\mathcal{P}(t,z)$.

In this paper we concentrate on finding explicit expressions for $\mathcal{P}(t,z)$ for Ginibre matrices, both real and complex. We first consider in Section 3 the case $Gin_{1}$ of real Ginibre matrices with $\beta=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, $\rho_{N}(z)$ for real Ginibre ensemble necessarily has the form $\rho_{N}(z)=\rho^{(c)}_{N}(z)+\delta\left(y\right)\rho^{(r)}_{N}(x)$, where the non-singular part $\rho^{(c)}_{N}(z)$ describes the mean density of complex eigenvalues, whereas $\rho^{(r)}_{N}(x)$ describes the mean density of purely real eigenvalues, so that $\int_{a}^{b}\rho^{(r)}_{N}(\lambda)d\lambda$ stands for the mean number of real eigenvalues in an interval $[a,b]$ of the real axis. As a consequence, the introduced JPD $\mathcal{P}(t,z)$ inherits the same structure $\mathcal{P}(t,z)=\mathcal{P}^{(c)}(t,z)+\delta(y)\mathcal{P}^{(r)}(t,x)$.

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 Section 2. We start our consideration with demonstrating in Section 3 a way to evaluate $\mathcal{P}^{(r)}(t,\lambda)$, which describes non-orthogonality factor for eigenvectors associated with a real eigenvalue $\lambda$ of the real Ginibre ensemble. First, we reduce the problem of finding $\mathcal{P}^{(r)}(t,\lambda)$ 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 Section 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\to\infty$. The problem of evaluating $\mathcal{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 4 we apply essentially the same method for evaluating $\mathcal{P}^{(c)}(t,z)$ in the complex Ginibre ensemble $Gin_{2}$, i.e. $\beta=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 Section 4.4. Eventually, we present an explicit finite-$N$ expression for any $z$, and then extract the corresponding ’bulk’ and ’edge’ scaling limits.

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 to the author’s attention, Gernot Akemann for pointing out and Peter Forrester for mentioning . Jacek Grela and Eugene Strahov are acknowledged for their collaboration on the associated analysis of Eq.(2.30) using different methods . 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”.

Discussion of the main results

Note that the left and right eigenvectors of real-valued matrices corresponding to real eigenvalues $\lambda$ can be chosen real as well. Hence we may write $\mathbf{x}_{\lambda,L}^{T}$ instead of $\mathbf{x}_{\lambda,L}^{*}$.

In what follows we will frequently omit the index $\lambda$ in eigenvectors to lighten the notations, simply writing $\mathbf{x}^{T}_{L}$ or $\mathbf{x}_{R}$.

The above expression is a generalization of the exact mean density of purely real eigenvalues $\rho^{(r)}_{N}(\lambda)$ for real Ginibre matrices of size $N$ explicit expression for which is known due to Edelman, Kostlan and Schub , see also :

and then introducing $u=|\lambda|\sqrt{\frac{t}{1+t}}$ as the integration variable, cf. (3.31) below.

Being exact, the expression (2.2) can be further analyzed in interesting scaling limits as $N\to\infty$. In fact, we find the form (2.5) most suitable for such an analysis. In particular, by rescaling $\lambda=\sqrt{N}x,\,\,t=Ns$ (which is standard to call the bulk scaling limit), then considering $x,s$ as fixed when $N\to\infty$ and exploiting the appropriate asymptotic behaviour of the incomplete $\Gamma$-function:

one easily finds that $\lim_{N\to\infty}N\mathcal{P}^{(r)}(Ns\,,\sqrt{N}x)=\mathcal{P}^{(r)}_{bulk}(s,x)$ where

and $\rho^{(r)}_{bulk}(x)=0$ otherwise, in full agreement with $\rho^{(r)}_{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 $\lambda=\sqrt{N}+\delta$, with $\delta<\infty$ being fixed. It is easy to understand that the variable $t$ needs to be rescaled in this regime as $t=\sqrt{N}\sigma$, keeping $\sigma$ fixed. A straightforward calculation using the well-known asymptotics

then yields $\lim_{N\to\infty}\sqrt{N}\mathcal{P}(\sqrt{N}\sigma,\,\,\sqrt{N}+\delta)=\mathcal{P}^{(r)}_{edge}(\sigma,\delta)$ where

In particular, integrating the above over $\sigma$ gives

where $\mbox{erf}(\delta)=\frac{2}{\sqrt{\pi}}\int_{0}^{\delta}e^{-u^{2}}\,du$. 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 .

2 Complex Ginibre ensemble

For the case of complex Ginibre ensemble $Gin_{2}$ the corresponding joint probability density $\mathcal{P}^{(c)}(t,z)$ of the non-orthogonality variable $t={\cal 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.

Let $z$ be a complex eigenvalue of $G$. Then the joint probability density $\mathcal{P}^{(c)}(t,z)$ of the self-overlap non-orthogonality variable $t={\cal O}_{z}-1$ and the associated complex eigenvalue $z$ for $N\geq 2$ is given by

where $D^{(N)}_{1}$ and $D^{(N)}_{2}$ are defined as

and $d_{1}^{(N)}$, $d_{2}^{(N)}$ are functions of $|z|^{2}$ explicitly defined via the relations to the incomplete $\Gamma-$function as

Note that $d^{(N)}_{1},d^{(N)}_{2}$ and $D^{(N)}_{1},D^{(N)}_{2}$ depend only on $|z|^{2}$ but not on the variable $t$.

The expression (2.13) is a generalization of the well-known mean density of eigenvalues of the complex Ginibre ensemble, see e.g. :

The joint density $\mathcal{P}^{(c)}(t,z)$ at fixed $z$ decays at large arguments $t\gg 1$ as $\mathcal{P}^{(c)}(t,z)\sim t^{-3}$ as was already anticipated by Mehlig and Chalker on the basis of $N=2$ example and informal eigenvalue repulsion arguments . 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=\sqrt{N}w,\,t=Ns$, with $w,s$ being fixed when $N\to\infty$, can be straightforwardly extracted. To this end, it is convenient to use the following integral representations for $d_{1}^{(N)}$ and $d_{2}^{(N)}$ (following from combining (2.16) and (2.17) with (2.4) and appropriate rescaling):

which for large $N\gg 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:

where $a_{N}\sim b_{N}$ means $\lim_{N\to\infty}a_{N}/b_{N}=1$. This implies $d_{1}^{(N-1)}\sim\frac{1}{N^{2}}\,d_{1}^{(N)},\quad d_{2}^{(N-1)}\sim\frac{2}{N}\,d_{1}^{(N)}$ 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 $\lim_{N\to\infty}N\mathcal{P}^{(c)}(t=Ns,z=\sqrt{N}w)=\mathcal{P}^{(c)}_{bulk}(s,w)$ one finds

and zero otherwise. This expression agrees with results obtained by P. Bourgade and G. Dubach in a different approach to the problem. Its first moment is precisely $\frac{1}{\pi}(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|=\sqrt{N}+\delta$ and $t=\sqrt{N}\sigma$ and performing the limit $N\to\infty$. With a help of the Mathematica packageThe author is grateful to J. Grela for his help with utilizing Wolfram Mathematica for that purpose. one then finds that

and we denoted $\Delta=1-2\sigma\delta$. In particular, one can check that integrating over $\sigma$ 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 $\delta=\frac{1}{2}\sqrt{N}(|w|^{2}-1)$ and $\sigma=\sqrt{N}s$ and letting $N\to\infty$ for fixed $|w|<1$ and $s$, checking that

3 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 and allows one to represent the JPD’s $\mathcal{P}^{(r)}(t,\lambda)$ and $\mathcal{P}^{(c)}(t,z)$, Laplace-transformed with respect to the variable $t$, in terms of the following object:

where $L=0,1,2,\ldots$ is an integer, $p\geq 0$, the parameter $\beta=2$ stands for the complex Ginibre ensemble and $\beta=1$ for the real Ginibre one (in the latter case $z=\lambda$ is real), and $I_{N}$ standing for the $N\times 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 $\beta=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. for a background discussion and further references. At the same time, for $\beta=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 for concise introductions and also 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 and the final expression for $\mathcal{D}_{N,1}^{(2)}(\lambda,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 $\beta=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 straighforward 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 Sec. 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:

and In fact the $L=1$ case was considered by a different variant of the supersymmetry approach in , though the result was not presented in the form 2.29.

In particular, by a direct integration one can check that $\mathcal{D}_{N}^{(L=1)}(z,p=0)=1$, in agreement with the definition (2.26).

Given the complexity of arising expressions for $\beta=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 $\beta=2$ case, namely $\mathcal{D}_{N,2}^{(L)}(z,p)$, can be formally rewritten as

with integration going over complex $N\times N$ matrices $W$ with the weight function $P_{L,z}\left(W,W^{*}\right)$ depending on an integer parameter $L=0,1,2,\ldots$, 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’ Formally the weight defined in (2.31) is not a probability measure for any $L\neq 0$ as it is not normalized to unity, but we disregard such difference for our goals. closely related (though not identical for $L>0$) to a limiting case of versions of the chiral ensemble with a ’source’ considered in and . Namely, let us consider a more general version of (2.31):

where the ’source’ matrix $A$ is a fixed $N\times 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},\ldots,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 , 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 . In a separate paper we are providing the full analysis of the problem for $\beta=2$ for any integer positive $L$ and $N$ by deriving the following representation (see Proposition 3.9 in ):

where we defined the following function of $\rho=|z|^{2}$ and $\tau=t/(1+t)$:

where ${\bf L}_{k}(x)$ are Laguerre polynomials. The equivalence with (2.28) - (2.29) for $L=0,1$ can be straightforwardly verified (see the Appendix A of ). One can further perform the asymptotic analysis of (2.33) - (2.34) for $N\gg 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 $\beta=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 $\beta=4$, those relevant in the theory of chaotic scattering and those relevant in the Quantum Chromodynamics context. We hope to be able to answer some of these questions in future publications.

Proof of Theorem 2.1

Let $\lambda$ be a real eigenvalue of a matrix $G^{(N)}$ with real entries, and denote the associated real left and right eigenvectors as $\mathbf{x}^{T}_{L}$ and $\mathbf{x}_{R}$. Then, as is well-known, see e.g. , it is always possible to represent the matrix $G^{(N)}$ as

2 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 $\mathcal{L}_{\lambda}(p):=\int_{0}^{\infty}e^{-pt}\mathcal{P}^{(r)}(t,\lambda)\,dt$ of the JPD $\mathcal{P}^{(r)}(t,\lambda)$ defined in (1.2) to evaluating the ensemble average for the ratio of certain determinants, see (3.8-3.9).

In a similar way one defines the complex Ginibre ensemble $Gin_{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 $\beta=1,2,4$ In the literature one frequently uses the notation $GinOE$ emphasizing an orthogonal symmetry of the distribution, and correspondingly $GinUE$ and $GinSE$ for complex and quaternion real versions of Ginibre ensemble with $\beta=2$ and $\beta=4$, correspondingly. one can write for all three ensembles the Joint Probability Density (JPD) with respect to the flat Lebesgue measure in the form

where $G^{*}=\overline{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 $\beta=1$; similar treatment of $\beta=2$ case will be briefly described in the last section.

with some normalization constant ${\cal 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 $\mathbf{b}_{\lambda}^{T}\mathbf{b}_{\lambda}$ which for a given value of $\lambda$ 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 ${\,\rm det}\>\left(\lambda\,I_{N-1}-G^{(N-1)}\right)^{T}={\,\rm det}\>\left(\lambda\,I_{N-1}-G^{(N-1)}\right)$. 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 $\beta=1$, which we also may present in an equivalent, but different form convenient for further evaluation:

where as will be found below (see the Corollary 3.4) $\mathcal{N}=2^{N/2}\Gamma(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)\times(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).

3 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\times N$, with the main object of interest being

with the constant $C_{N}=\frac{1}{2^{N/2}\Gamma\left(\frac{N}{2}\right)}$.

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\to 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 $\mathcal{P}^{(r)}(t,\lambda)$ of the random variables $t=\mathbf{b}_{\lambda}^{T}\mathbf{b}_{\lambda}$ and the real eigenvalue $\lambda$ as is given in (2.5). Proof of the Proposition 3.3:

Let $\Psi_{1},\,\Psi_{2},\,\Phi_{1},\,\Phi_{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 real-symmetric matrix $A$ and any positive $\epsilon>0$:

where the integration goes over the vector $\mathbf{S}$ with $N$ real commuting components. This allows to represent the denominator in (3.10) as a Gaussian integral over two such vectors $\mathbf{S}_{1,2}$ :

where $\propto$ here and below stands for (temporaly) 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

etc, where $M=\mathbf{a}\otimes\mathbf{b}^{T}$ stands for the matrix with entries $M_{ij}=a_{i}b_{j}$, one can easily perform the averaging over the real Ginibre matrices by using the identity

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 $\left(\Phi_{1}^{T}\Phi_{2}\right)\left(\Psi_{1}^{T}\Psi_{2}\right)$. 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 $\mathbf{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,\overline{q}$ is now easy to perform via using the polar coordinates:

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 in such a situation is to pass from the pair of vectors $(\mathbf{S}_{1},\mathbf{S}_{2})$ to the matrix $\hat{Q}$ as a new integration variable. Such change is non-singular for $N\geq 2$ and incurs a Jacobian factor proportional to ${\,\rm det}\>{\hat{Q}}^{(N-3)/2}$ (see the Appendix D of ). This finally brings $\mathcal{D}^{(2)}_{N,1}(\lambda,p)$ to the form

The next step requires employing a convenient parametrization of the integration domain defined via the inequalities $Q_{1}\geq 0,\,Q_{2}\geq 0,\,-\infty<Q<\infty$ and $Q_{1}Q_{2}\geq Q^{2}$ which ensure that $\hat{Q}=\begin{pmatrix}Q_{1}&Q\\ Q&Q_{2}\end{pmatrix}$ 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\geq 0,\,-\infty<\theta<\infty$ chosen in such a way that $r=\left({\,\rm det}\>\hat{Q}\right)^{1/2},\,\,\theta=\frac{1}{2}\ln{\left(Q_{1}/Q_{2}\right)}$. By explicitly writing $Q_{1}=e^{\theta}\sqrt{r^{2}+Q^{2}},\,Q_{2}=e^{-\theta}\sqrt{r^{2}+Q^{2}}$ and evaluating the associated Jacobian we get in those coordinates $d\hat{Q}:=dQ_{1}dQ_{2}dQ=2\,r\,dr\,dQ\,d\theta$. 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}\geq 0$ and $-\infty<Q<\infty$ as new coordinates, complemented with $r=\left({\,\rm det}\>\hat{Q}\right)^{1/2}\geq 0$, and expressing the remaining entry as $Q_{2}=\frac{r^{2}+Q^{2}}{Q_{1}}\geq 0$. This finally gives

and further changing $Q_{1}\to\sqrt{2p}Q_{1}$ brings (3.20) to the form

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 , 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}\to t$ the equation (3.23) can be rewritten as

Now by using the relation $\Gamma(N+1,\lambda^{2})=e^{-\lambda^{2}}\lambda^{2N}+N\Gamma(N,\lambda^{2})$ one can see that

for some real constant $C_{N}$. Simple manipulations with incomplete $\Gamma-$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\to\infty$ where according to the definition (3.10)

On the other hand, performing $p\to\infty$ limit in (3.27) gives after a simple calculation

The definition of the left-hand side implies that the coefficient in front of $\lambda^{2N}$ must be equal to unity, giving finally

The normalization constant $\mathcal{N}$ in Eq.(3.8) is given by $\mathcal{N}=2^{N/2}\Gamma(N/2)$.

To establish the value of the constant $\mathcal{N}$ we consider the limit $p\to 0$ in both sides of (3.8). By the very definition of the Laplace transform $\mathcal{L}_{\lambda}(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=|\lambda|\sqrt{\frac{t}{1+t}}$ as new integration variable. One gets in this way:

To get $\mathcal{D}_{N-1,1}^{(2)}(\lambda,0)$ featuring in the right-hand side of (3.8) replace in the above $N\to N-1$ and use $\Gamma(N)=\frac{2^{N-1}}{\sqrt{\pi}}\Gamma\left(\frac{N}{2}\right)\Gamma\left(\frac{N+1}{2}\right)$. Multiplying with $e^{-\lambda^{2}/2}$ and comparing with the left-hand side gives the value for the constant $\mathcal{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 non-degenerate eigenvalues, and assuming it has an eigenvalue $z$ (in general, complex) it can be represented as (see e.g. Sec. 6 of )

Now we exploit the analogue of Prop. (3.2)

with some normalization constant ${\cal C}_{2,N}$.

Using the above j.p.d. to calculate the Laplace transform of the probability density for the random variable $\mathbf{b}_{z}^{*}\mathbf{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 $\beta=2$,

with averaging performed over the j.p.d. (3.4) of complex Ginibre matrices $G$ of the size $N\times N$. Note, that the value of the constant normalization factor in (4.4) is found aposteriori 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 $\mathcal{L}_{z}(0)$) and the corresponding limit in the right-hand side of (4.4).

In the general case evaluating $\mathcal{D}_{N,2}^{(L)}(z,p)$ for integer $L$ can be done essentially by the same supersymmetry method which was used in section (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 or, in a different way, in . The corresponding calculation is sketched in the first part of the next section. The supersymmetry approach for $z\neq 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.

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 , 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\sqrt{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 $\mathcal{D}_{N,2}^{(L=2)}(z=0,p)$ as

Now the integrals over $R_{1}$ and $R_{2}$ can be readily evaluated, with the result being simply

Further employing a well-known integral representation for the Bessel function of the second kind

which after substituting to the Laplace transform (4.4) is equivalent to (2.13).

1.2 Evaluation of (4.5) for L=2𝐿2L=2 and |z|≠0𝑧0|z|\neq 0 by supersymmetry approach

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.

One starts with using two copies of the set of four anticommuting vectors, namely $\Psi_{A1},\,\Psi_{A2},\,\Phi_{A1},\,\Phi_{A2}$ and $\Psi_{B1}\,\Psi_{B2},\,\Phi_{B1},\,\Phi_{B2}$, to represent separately two determinants in the numerator via Gaussian integrals, see (3.12). At the same time, one needs two commuting vectors ${\bf S}_{1},{\bf S}_{2}$ with $N$ complex-valued components to represent the denominator:

The ensemble averaging is performed by exploiting $\beta=2$ analogue of (3.15)

Performing the average one collects all terms in the exponential which are quartic in anticommuting variables, e.g. $\left(\Psi_{A1}^{T}\Phi_{A1}\right)\left(\Psi_{A2}^{T}\Phi_{A2}\right)$, etc.. The corresponding exponential factor can be then represented via a matrix version of the Hubbard-Stratonovich transformation generalizing (3.16):

with $\hat{Q}_{F},\hat{Q}_{F}^{*}$ being a pair of general $2\times 2$ complex conjugate matrices. This trick allows to integrate out the vectors with anticommuting component completely. The analogue of (3.17) takes the form

At the next step one can simplify the above expression by employing the singular value decomposition $\hat{Q}_{F}=U\mbox{diag}(\sqrt{R_{F1}},\sqrt{R_{F2}})V^{*},\hat{Q}_{F}=V\mbox{diag}(\sqrt{R_{F1}},\sqrt{R_{F2}})U^{*}$ with unitary $U,V$ and $R_{F1},R_{F2}\geq 0$ and replacing the integration over complex vectors $\textbf{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\geq 2$ as

The Hermitian matrix $\hat{Q}\geq 0$ can be parametrized very similarly to (3.21). Namely, writing for the complex variable $Q=\rho\,e^{i\phi}$ and using $r=\left({\,\rm det}\>{\hat{Q}}\right)^{1/2}\geq 0$ together with $Q_{1}\geq 0$ as the coordinates, the domain of integration is parametrized by matrices

In this way we arrive at an analogue of (3.23):

The integrals over $\rho$, $r$, $dR_{F1}$ and $dR_{F2}$ can be performed. Namely, denoting $r^{2}=R_{B}$ and $\rho^{2}=R$ we introduce a function, cf. (3.24),

After further denoting $e^{|z|^{2}}\Gamma(N,\lambda^{2}):=\gamma\left(N,|z|^{2}\right)$ we then define the functions

and renaming $Q_{1}\to t$ we notice that (4.17) can be represented, after restoring the normalization constants, as

Substituting (4.19) to the above and taking the limit, the expression can be further simplified with the help of symbolic manipulations using Wolfram Mathematica and is fianlly represented as (4.9). In particular, it is easy to see that $D^{(N+1)}_{1}|_{|z|=0}=N(N+1)d_{1}^{(N+1)}|_{|z|=0}=N!(N+1)!$ so that (4.9) at $|z|=0$ indeed reproduces (4.8). ∎

References