Probing non-orthogonality of eigenvectors in non-Hermitian matrix models: diagrammatic approach

Using large $N$ arguments, we propose a scheme for calculating the two-point eigenvector correlation function for non-normal random matrices in the large $N$ limit. The setting generalizes the quaternionic extension of free probability to two-point functions. In the particular case of biunitarily invariant random matrices, we obtain a simple, general expression for the two-point eigenvector correlation function, which can be viewed as a further generalization of the single ring theorem. This construction has some striking similarities to the freeness of the second kind known for the Hermitian ensembles in large $N$. On the basis of several solved examples, we conjecture two kinds of microscopic universality of the eigenvectors - one in the bulk, and one at the rim. The form of the conjectured bulk universality agrees with the scaling limit found by Chalker and Mehlig [JT Chalker, B Mehlig, PRL, \textbf{81}, 3367 (1998)] in the case of the complex Ginibre ensemble.


Introduction
Non-normal operators are ubiquitous in physical models. Examples include hydrodynamics, open quantum systems, PT-symmetric Hamiltonians, Dirac operators in the presence of a chemical potential or finite angle θ. Non-normality is responsible for the transient dynamics, sensitivity of the spectrum to perturbations, pseudoresonant behavior and rapid growth of the perturbations of the system [1]. These effects are relevant in plasma physics [2], fluid mechanics [3], ecology [4,5], laser physics [6], atmospheric science [7], and magnetohydrodynamics [8], just to mention a few. Non-normality is common in dynamical systems as its simplest source is the asymmetry of coupling between components.

JHEP06(2018)152
Historically, most of the studied properties of non-normal random operators dealt with the eigenvalues. The eigenvalues of such operators are usually complex, requiring new calculational techniques, at the level of both macroscopic and microscopic correlations. Surprisingly, this quest for complex eigenvalues has eclipsed the study of eigenvectors, which are perhaps most distinctive features of non-normal operators. In particular, nonnormal operators have two sets of eigenvectors, left and right, which are non-orthogonal among themselves, but can be chosen to be bi-orthogonal, provided that the non-normal operator can be diagonalized at all.
One of the first attempts to develop a systematic understanding of the non-orthogonality of eigenvectors in non-Hermitian random matrices was made by Chalker and Mehlig [9,10]. Despite their study concentrated on the complex Ginibre ensemble, perhaps the simplest non-normal random operator, the results turned out quite non-trivial and revealed the possibilities of well-hidden universal properties of eigenvectors of nonnormal operators. Another connection of the properties of non-normal operators and their eigenvectors to free probability was established soon after [11], but the systematic study of this topic has not followed. Only very recently, the topic of eigenvectors of non-normal operators was picked back up. First, the transient growth driven by eigenvector nonorthogonality was proposed as a mechanism of amplification of neural signals in balanced neural networks [12][13][14] and giant amplification of noise crucial in the formation of Turing patterns [15][16][17]. Second, the non-orthogonality of eigenfunctions was related to the statistics of resonance width shifts in open quantum systems [18], which was soon confirmed experimentally [19]. Third, the essential role of eigenvectors in stochastic motion of eigenvalues was revealed [20][21][22]. Last but not least, the topic has triggered the attention of the mathematical community [23,24].
In this work we focus on statistical ensembles of complex non-Hermitian matrix models, the probability density of which is invariant under the action of the unitary group P (X) = P (U XU † ). We also assume that in the N → ∞ limit, at which we are working, the eigenvalues of X concentrate on a compact domain of a complex plane. Our results are valid for |z − w| of order 1. We will study one-point and two-point Green functions built out of left and right eigenvectors. Here we recall, that if a non-normal matrix X can be diagonalized by a similarity transformation X = SΛS −1 , it possesses two eigenvectors for each eigenvalue λ i : right |R i (a column in the matrix notation) and left L i | (row), satisfying the eigenequations These eigenvectors are not orthogonal L i |L j = δ ij = R i |R j , but normalized by the biorthogonality condition L i |R j = δ ij . They also satisfy the completeness relation N k=1 |R k L k | = 1. These two properties leave a freedom of rescaling each eigenvector by a non-zero complex number, |R i → c i |R i with L i | → L i | c −1 i . They also allow for multiplication by a unitary matrix |R i → U |R i and L i | → L i | U † . Upon the second transformation the new vectors are not the eigenvectors of the original matrix but of one given by the adjoint action of the unitary group X → U XU † , which suggests that a natural probability measure should assign these two matrices the same probability density JHEP06(2018)152 function (pdf). The simplest object, which is invariant under these transformations, is the matrix of overlaps O ij = L i |L j R j |R i [9,10].
To see how the eigenvector correlation functions appear naturally, let us consider an average 1 N Trf (X)g(X † ) , where f, g are two functions analytic in the spectrum of X and f (X) = f (X)P (X)dX denotes the average with respect to the pdf P (X). Taking f = g, we get the (normalized) Frobenius norm of a function of matrix. The 1/N normalization was taken to get a finite quantity in the N → ∞ limit. Using the spectral decomposition X = N k=1 |R k λ k L k | and inserting the identity, 1 = dµ(z)δ (2) (z − λ k ) twice, we obtain the expression 1 N Trf (X)g(X † ) = dµ(z)dµ(w)f (z)g(w)D(z, w), The two dimensional Dirac delta is understood as two deltas for real and imaginary parts δ (2) (z) = δ(Rez)δ(Imz), and the measure dµ(z) = dxdy for z = x + iy. D(z, w) introduced in [9,10] is the density of eigenvalues weighted by the invariant overlap of the corresponding eigenvectors. It is split into a regular and singular part D(z, w) =Õ 1 (z)δ (2) (1.4) A one-point function, defined this way, in the bulk and far from the rims of the complex spectra grows linearly with the size of a matrix. To have a finite limit in large N , one considers the scaled function O 1 (z) = 1 NÕ 1 (z). Throughout the paper we shall use only the 'untilded' function.
The one-point function O 1 plays an important role in scattering in open chaotic cavities [18,25] and random lasing [26,27], where the so-called Petermann factor [28] modifies the quantum-limited linewidth of a laser. It is also crucial in the description of the diffusion processes on matrices [21,22] and gives the expectation of the squared eigenvalue condition number [29], an important quantity governing the stability of eigenvalues [1,30]. The exact calculations are possible for Gaussian matrices [9,10,23], in the matrix model for open chaotic scattering [26,27,31] and for products of small Gaussian matrices [32]. For the Ginibre matrix the full distribution of the diagonal overlap is available and turns out to be heavy-tailed, as discovered by Bourgade and Dubach [24] with the use of probabilistic techniques, and investigated later using integrable structure and sypersymmetry by Fyodorov [33].
Despite that the overlap between eigenvectors are crucial in the description of the dynamic of the linear system [34] and in the decay laws in open quantum systems [35], the two-point function is much less known. The exact results are obtained only for the Ginibre JHEP06(2018)152 matrix [9,10,24] and for open chaotic scattering with a single channel coupling [31]. Even the asymptotic results are known only for Gaussian matrices [9,10] and the quantum scattering ensemble [36]. The aim of this paper is to extend the known asymptotic results and develop a diagrammatic technique for calculation of the two-point function in the large N limit.
The paper is organized as follows. In section 2 we briefly recall the cornerstones of diagrammatic calculus [11] for one-point Green's functions in the large non-Hermitian ensembles, to show an analogy between the formalism developed in this paper and the diagrammatic approach to one-point functions. This encapsulates both the mean spectral density and the one-point eigenvector correlation function O 1 . Appendix A shows concrete calculations within this formalism for the elliptic ensemble.
Section 3 contains the main body of the paper -a formalism for the calculation of O 2 in the large N limit. We extend the diagrammatic technique introduced by Chalker and Mehlig for Gaussian matrices to any probability distribution with unitary symmetry. Regularizing and linearizing the product of resolvents, we embed them into the quaternionic space. The analysis of planar Feynman diagrams leads us to the matrix Bethe-Salpeter equation (3.15), which relates the product of resolvents with the one-point Green's function and planar cumulants. The latter are encoded in their generating function -quaternionic R-transform, see (3.16). As a result, the two-point eigenvector correlation function is completely determined by the one-point functions encoding the spectral density and O 1 . This result resembles the Ambjørn-Jurkiewicz-Makeenko universality for Hermitian ensembles [37].
We also study the traced product of resolvents h(z,w) = 1 N Tr(z1−X) −1 (w1−X † ) −1 , which allows for the calculation of the average (1.2) as a Dunford-Taylor integral [38,39] where contoursγ,γ encircle all eigenvalues of X clockwise and counterclockwise, respectively. We derive the equation for h, expressing it in terms of quaternionic R-transform and traced resolvents, see (3.18) and (3.19). An important and still quite large subclass of non-Hermitian ensembles for which the main equations (3.15)(3.16) admit further simplifications consists of matrices, the pdf of which is invariant under the transformation by two independent unitary matrices U, V ∈ U(N ), i.e. P (X) = P (U XV † ), thus called the biunitarily invariant ensemble [40]. In this case we obtain a compact formula for the two-point eigenvector correlation function Here F is the radial cumulative distribution function (cdf), defined as F (r) = 2π r 0 ρ(s)sds, with ρ(s) the spectral density circularly symmetric on the complex plane. The one-point eigenvector function is related to F via [29]

JHEP06(2018)152
and r = |z|. The traced product of resolvents is shown to take a universal form where r out is the spectral radius. This result has been already obtained for the Ginibre ensemble [10] and, recently, for matrices with independent identically distributed (iid) entries [41]. Applications of the developed formalism are presented in section 4, where we consider an elliptic ensemble, some instances from the biunitarily invariant class: truncated unitary, induced Ginibre, the product of two Ginibres and their ratio. As the last example we consider a pseudohermitian matrix -a product of two shifted GUE matrices. In section 5, we discuss the consequences of our large N results on the microscopic regime. We conjecture, on the basis of the few examples solved in the literature and using our own results, that the two-point eigenvector correlation functions may exhibit universal bulk scaling, as what happens for the microscopic spectral two-pointers in Hermitian matrix models. More precisely, we conjecture that in generic complex non-Hermitian matrices for all points in the bulk at which the spectral density does not develop singularities there exists a limit lim N →∞ (1.10) Section 6 concludes the paper and points at some possible further developments.

Non-Hermitian random matrices
In non-Hermitian random matrix theory one is primarily interested in the distribution of the eigenvalues ρ(z) = 1 The 2-dimensional Dirac delta can be recovered using the relation ∂z 1 z = πδ (2) (z). Unfortunately, the average over the ensemble of the trace of the resolvent g(z) = 1 N Tr(z1 − X) −1 does not yield the correct result inside the spectrum, as one would naively expect. The reason for this failure is that differentiation and taking the ensemble average are not interchangeable. This phenomenon was termed the spontaneous breaking of holomorphic symmetry [42].
A way to circumvent this obstacle is to consider a regularization of the Dirac delta. In RMT one mostly considers the 2D Poisson kernel The expression on the right hand side provides a prescription for how the resolvent in the spectrum of X should be regularized. Having this hint in mind, one defines

JHEP06(2018)152
The spectral density can be now calculated via ρ(z,z) = 1 π lim |w|→0 ∂zg(z,z, w,w), which can be also understood as a Poisson law in 2D electrostatics, since is the (regularized) electrostatic potential of charges interacting via repulsive central force F (r) ∼ 1 r .

Linearization
Due to the quadratic expression in X in the denominator, the average in (2.2) is intractable when non-normal matrices are considered. To circumvent this problem one introduces the 2N × 2N matrix [42][43][44][45]  Then, one defines the 2 × 2 Green's function Its upper-left entry is exactly the desired function g (cf. (2.2)). Once Green's function is known, one gets four elements of G. The entry G11 is just the complex conjugate of G 11 , thus does not provide any additional information. The off-diagonal entries G 11 = −Ḡ1 1 in the large N limit give the one-point eigenvector correlation function [11]

Quaternionic structure
Green's function can be conveniently written as This form of Green's function resembles its traditional form as a traced resolvent, but now its argument is a 2 × 2 matrix and the matrix X is duplicated to incorporate also X † . The matrix Q is a representation of a quaternion q = x + iy + ju + kv with the identification z = x + iy and w = v + iu [46]. The entries of G satisfy the same algebraic constraints as Q, therefore G is itself a quaternion and we refer to it as the quaternionic Green's function, similarly G is called the quaternionic resolvent.

Averages in large N
We are interested in calculations of the averages of some functions of X, e.g. f (X, X † ) , with respect to distributions invariant under the adjoint action of the unitary group P (X) = P (U XU † ). We parameterize them by V (X, X † ), often referred to as potential, has to be Hermitian and growing sufficiently fast at infinity. To simplify calculations, we split the potential into the Gaussian and the residual part. The Gaussian part can be conveniently parameterized with σ > 0 and τ ∈ [−1, 1] [47] Averages with respect to the Gaussian potential by the virtue of Wick's theorem reduce to products of pairwise expectations, called propagators 14) The exponent of the residual part of the potential is expanded into series, which produces additional terms, called vertices, to be averaged with respect to the Gaussian distribution.
To cope with the multitude of terms, we represent them as diagrams (see table 1 for an overview). This introduces a natural hierarchy of diagrams according to their scaling with the size of the matrix. The dominant contribution, which is of the order of 1, comes from planar diagrams (see figure 1). The subleading corrections can be classified by the genus of the surface at which they can be drawn without the intersection [48].

Moment expansion of the quaternionic resolvent
To calculate the average of the quaternionic resolvent, we write it as and expand it into the geometric series propagator Table 1. Diagrammatic representation of the basic expressions in the moment expansion of the resolvent. The propagator represents the averages with respect to the Gaussian potential (2.14). An exemplary vertex is drawn for the theory which contains the cubic interaction N g 3 TrX α X β X γ in the potential. A cumulant (dressed vertex) represents a sum over all connected diagrams connected to the baseline. Its structure in matrix indices (Latin letters) is the same as that of the vertex, because the propagators are the Kronecker deltas in this indices. The dashed line without arrows represent summation over Latin indices only. and perform averages in the large N limit, as described in the previous section. The expansion is valid, provided that ||Q −1 X || < 1, thus for z inside the spectrum of X, we need to keep w finite. If the spectrum is bounded, one can always find sufficiently large w, so that this series is absolutely convergent. For the calculations with z outside the spectrum one can safely set w = 0. It is convenient to introduce a notation, which incorporates the block structure of the duplicated matrices. We therefore endow each matrix with two sets of indices, writing for example G αβ,ij . The first two Greek indices, which we also refer to as quaternionic indices, enumerate blocks and take values 1 and1. The Latin ones, running from 1 to N enumerate matrices within each block. The space described by the Latin indices we call simply the matrix space. The block trace operation can be represented as a partial trace over the matrix space (see also table 1). Due to the fact, that the propagators are expressed in terms of Kronecker deltas, all averaged expressions have trivial matrix structure, e.g. G = G ⊗ 1, but we need this notation for the next section.
Among all diagrams contributing to G (see figure 2 for an example) we distinguish a class of one-line irreducible diagrams (1LI), i.e. the ones that cannot be split into two JHEP06(2018)152 . . Figure 2. Some exemplary planar diagrams in a model with a quartic term g 4 X 4 contributing to the Green's function. All diagrams (except for the first) are 1PI. parts, connected only through Q −1 . Let us denote as Σ a sum of all 1LI diagrams. This is a building block of the quaternionic resolvent, since any diagram can be decomposed into 1LI subdiagrams connected through Q −1 . Having the absolute convergence of the series, we rearrange terms, obtaining the Schwinger-Dyson equation (here we restrict it only to the quaternionic part) presented also diagrammatically in figure 3a). This is a geometric series, which can be summed and written in a closed form

Quaternionic R-transform
To find G, one needs to relate Σ to G. To this end, let us consider diagrams contributing to averages of traced strings of X's and X † 's such that all X's and X † 's are connected with each other. Their sum we call a cumulant (in field theory language it is known as a dressed vertex) and endow the respective average with a subscript c. We adopt a convenient notation for cumulants in which † is associated with the1 index and, trivially,

JHEP06(2018)152
lack of conjugation with 1. We also encode the length of the string. An example reads We also introduce the R-transform, which is a 2 × 2 matrix, defined through its expansion for small arguments This definition written in terms of indices may not seem to be intuitive, but in the matrix notation takes a clearer form which is the counterpart of (2.15). The matrix R is also a quaternion, so it is dubbed the quaternionic R-transform. Q is always associated with two consecutive indices in the cumulant and can be thought of as a transfer matrix. It is crucial for encoding all cumulants in the R-transform that matrices X and Q do not commute. Consider now any 1LI diagram. Due to its irreducibility it can be depicted as a certain cumulant connecting the first and last X and possibly some others in between. The subdiagrams disconnected from the cumulant can be in any form, which is already encoded in the quaternionic Green's function. This allows us to write the second Schwinger-Dyson equation relating Σ and G via the quaternionic R-transform (see also figure 3b)) Σ(Q) = R(G(Q)). (2.21) The knowledge of all cumulants allows us to solve the matrix model, since equations (2.17) and (2.21) can be brought to a single 2 × 2 matrix equation Once the averaging with respect to the ensemble was taken at the level of diagrams, we can safely remove the regularization and solve the above algebraic equation, setting first w = 0. We refer to [49,50] for more detailed calculations in the diagrammatic formalism. The construction presented in this section has been recently rigorously formalized in the framework of free probability [51].

Preliminaries
In order to investigate the 2-point eigenvector correlation function associated with the off-diagonal overlap, we follow the paradigm outlined in the previous section for calculations of Green's function. A naive approach, i.e. calculation of h(z 1 , spectrum of X, which we refer to as the holomorphic solution. Inside the spectrum, we regularize each resolvent, using the rule The weighted density is therefore given by In this paper we will calculate h by diagrammatic 1/N expansion in the planar limit. The singular part of D(z 1 , z 2 ) containing the Dirac delta is not accessible in perturbative calculations, so we get There exists a class of matrices which despite not being Hermitian have a real spectrum. A simple example is the product of two Hermitian matrices A, B, one of which (let us say A) is positive definite. The resulting matrix is not Hermitian, but isospectral to A 1/2 BA 1/2 , which must have real eigenvalues. The eigenvectors of AB are not orthogonal, which makes O 2 non-trivial. The realness of the spectrum facilitates calculations, as the knowledge of the traced resolvent is sufficient. By the virtue of the Sochocki-Plemelj formula valid for real x we can write 5) and the two-point function can be calculated from the holomorphic function via where h(±, ±) = lim and signs are uncorrelated.

Linearization
The expression for the regularized product of resolvents (3.2) contains two quadratic nonlinearities. We overcome this difficulty, by using 2N × 2N matrices Q = Q ⊗ 1, P = P ⊗ 1 and X , where

JHEP06(2018)152
As a natural generalization of the quaternionic resolvent to two-point functions, we define the average of the Kronecker product of two quaternionic resolvents Such an object is quite unusual in Random Matrix Theory. A similar construction was used by Brezin and Zee for the calculation of the connected 2-point density in Hermitian models [52], but there one deals only with matrix indices. To the best of our knowledge the quaternionic approach to two-point functions for non-Hermitian matrices is considered for the first time, thus we will discuss it in more detail.
H is a 4N 2 × 4N 2 matrix with a very specific block structure. To keep track of it, we endow H with 8 indices. The upper ones refer to the first matrix in the Kronecker product, while the lower ones to the second. As in the case of the quaternionic Green's function, Greek indices, taking values in {1,1}, enumerate blocks, while Latin indices ranging in {1, . . . , N } denote elements within each block. In the index notation, its components read (note the transpose of the second matrix) With the same assumptions as for one-point functions, the resolvents are then expanded into the power series and taking the expectation produces diagrams. The flow of Latin (matrix) indices in the diagrams follows the lines in the double line notation. The propagators are symmetric, thus the direction does not matter. The flow of quaternionic (Greek) indices is governed by their order in the expansion of the resolvent. Since the quaternion matrices Q and P are not symmetric, the direction of the line representing Q −1 matters and is depicted by an arrow. We draw diagrams in such a way that the terms in the expansion of the resolvents are in two rows, hereafter called baselines, with the first resolvent above. The quaternionic indices flow from left to right in the upper baseline and in the opposite direction below. There are two ways of contracting matrix indices, 1 thus we define two functions which correspond to contractions presented in figure 4a). It will become clear later that K encodes correlations of eigenvectors and L of eigenvalues. These two possible contractions define two different classes of planar diagrams. The admissible diagrams have to be drawn in the region of the plane bounded by baselines and dashed lines depicting contractions. The diagrams contributing to K are of the ladder type (see figure 5), while the class of planar diagrams contributing to L, termed wheel diagrams, is broader, as it admits for  [52,53] and K is given as a sum of ladder diagrams. b) An example of a diagram which contributes to L but is subleading in the calculation of K. c) An example of a diagram appearing during the calculation of L, which despite its planarity is subleading. circumventing one of the baselines if the points on the baseline are connected through propagators and vertices, see figure 4b). Not all planar diagrams contribute equally to L. Diagrams in which a propagator connects two sides of a vertex and encircles a baseline is subleading, see figure 4c). In this section we concentrate on the ladder diagrams.

Ladder diagrams
In this section we are interested in the calculation of K. The contraction of matrix indices in H, which leads to K, is in fact a summation of all N 4 elements within each 4 × 4 block.
To make the notation of Greek indices even more explicit, we write the entries of K An important consequence of this construction is that the K 11 11 element is exactly the desired function h (3.2) for the calculation of the eigenvector correlation function.
Let us define K αβ,ij µν,kl the sum of all ladder diagrams contributing to K (before we contract indices). A vertex can connect two points on a baseline (a side rail of the ladder), dressing the part of the rail. There are also vertices connecting two baselines, which give rise to the rungs of the ladder. If we denote Γ αβ,ij µν,kl a sum of all connected subdiagrams JHEP06(2018)152 which connect two rails, one can express K in terms of Γ as a geometric series, presented in figure 6, which can be written in a closed form (a sum over repeating indices is implicit) This relation, shown diagramatically in figure 7 and known as the matrix Bethe-Salpeter equation, is the counterpart of the Schwinger-Dyson equation for the two-point function, with Γ the counterpart of the self-energy.
A direct analysis of planar diagrams yields Γ αβ,ij µν,kl = 1 N B αβ µν δ i k δ j l , where B is of order 1, see figure 8. Using the matrix structure of Γ, we trace out the matrix indices and find the equation for K, which in the matrix notation reads K(Q, P ) = G(Q) ⊗ G T (P ) + G(Q) ⊗ G T (P ) B(Q, P )K(Q, P ). (3.15) We now turn our attention to the rungs. Any diagram contributing to Γ can be decomposed as a certain cumulant of length n ≥ 2, the first k legs of which are attached to the upper rail, while the last legs are connected to the lower rail. The part of the rail between the legs of the cumulant gets dressed to the quaternionic Green's function G(Q) above and G(P ) below. The space between k-th and (k + 1)-th legs is left unfilled, because the quaternonic indices at the end of rails are not contracted. This decomposition of Γ is depicted in figure 9. As Γ is completely determined by the planar cumulants, B αβ µν can be calculated from the quaternionic R transform (2.19). The rule is simple and goes as follows.
Consider the expansion of R αµ in Q (2.19) and differentiate it with respect to Q βν . As a result for some 0 < k < n − 1 the k-th quaternion Q µ k µ k+1 will be replaced by δ µ k β δ νµ k+1 . Then all Q µ l µ l+1 's from the l.h.s. of the removed Q (i.e. for l < k) are replaced by G µ l µ l+1 (Q) and all Q µ l µ l+1 on the right (l > k) by G T µ l µ l+1 (P ). Then the sum over all possible positions (i.e. k's), where Q has been removed, is performed. JHEP06(2018)152 g 4 Figure 8. An example of a diagram contributing to Γ. It contributes to the second-to-last diagram in figure 9. Since the matrix indices follow the solid lines and propagators are given by Kronecker deltas, Γ αβ,ij µν,kl = 1 N B αβ µν δ i k δ j l , allowing for the calculation of K. B can be therefore expressed in terms of cumulants as a power series 16) where all σ i and ρ j take values in {1,1} and for k = 1 or l = 1 G σ k σ k+1 reduces to Kronecker delta. An application of this procedure to the quantum scattering ensemble is presented in appendix B.
We remark that the additivity of the quaternionic R-transform under the addition of unitarily invariant non-Hermitian matrices implies additivity of B.

Traced product of resolvents
In the holomorphic domain outside the spectrum the situation simplifies considerably, because one can set |w| → 0 at the very beginning of calculations. Green's function is Due to such a structure, B is also diagonal with components where we assume the standard convention g 1 (z) = g(z) and g1(z) =ḡ(z). A matrix equation (3.15) splits into decoupled scalar equations with the explicit solution for the component of our interest (3.18)

JHEP06(2018)152
The desired component of B obtained from (3.17) reads Despite the fact that the mapping between cumulants and the R-transform is not one to one [49], some cumulants can be uniquely determined from the knowledge of R(Q). The cumulants c in the expansion of R 11 (Q).
One can easily see that there are no other cumulants contributing to this term.
All cumulants contributing to R 11 have at least one X † following X in the string, therefore R 11 is divisible by Q 11 . Let us defineR 11 = R 11 /Q 11 , which is regular at 0. The considered cumulants are the only ones in which X is followed by X † exactly once. To exclude all other possibilities in the expansion ofR 11 , we set Q 11 = 0 = Q1 1 inR 11 (Q). To reproduce (3.19) fromR 11 one also needs to replace Q 11 by g(z 1 ) and Q11 byḡ(z 2 ). Finally, B 11 11 =R 11 (diag(g(z 1 ),ḡ(z 2 ))) . (3.20)

Biunitarily invariant ensembles
In this subsection we consider a class of ensembles, the pdf of which is invariant under multiplication by two unitary matrices, i.e. P (X) = P (U XV † ). In the large N limit the spectral density, which is rotationally invariant, is supported on either a disc or an annulus, a phenomenon termed 'the single ring theorem' [54,55]. The enhanced symmetry allows one to relate the distribution of eigenvalues and singular values both in the N → ∞ limit [56] and for finite N [40]. More precisely, the radial cumulative distribution function F (r) = 2π r 0 ρ(s)sds is given by the solution of the simple functional equation S XX † (F (r)−1) = 1 r 2 , where S XX † is the Voiculescu S-transform of the density of squared singular values [56]. Recently, this result was extended to the one-point eigenvector correlation function, which is determined solely by F [29,49] Such simple results in the large N limit are possible because of the exceptionally simple structure of cumulants. The only non-zero planar cumulants are the alternating ones [49], α n = c

JHEP06(2018)152
A direct application of formula (3.16) leads to The components of B can be expressed through auxiliary functions with A being the determining sequence. We remark that the average over the ensemble has been already taken at the level of Feynman diagrams and at this moment, we can safely remove the regularization. There are further simplifications for the biunitarily invariant matrices [49] G 11 G1 1 A(G 11 G1 1 ) = F (r) − 1, G 11 G1 1 = −πO 1 (r). (3.30) Having calculated B and knowing Green's function, we determine K 11 11 from (3.15) and, after algebraic manipulations, we get a compact formula for the 2-point eigenvector correlation function from (3.4) The quaternionic R-transform of biunitarily invariant ensembles takes a remarkably simple form [49], in particular R 11 (Q) = Q 11 A(Q 11 Q1 1 ). Moreover, due to the rotational symmetry of the spectrum, g(z) = 1/z. According to (3.18), the traced product of resolvents is given by Interestingly, A(0) = r 2 out , where r out is the external radius of the spectrum. This result shows a high level of universality, since for any two functions f, g analytic in the spectrum the expectation in the N → ∞ limit

JHEP06(2018)152
is given by the same formula, irrespectively of the specific biunitarily invariant ensemble. The only parameter -spectral radius r out -can be set to 1 by rescaling the matrix. This result, appearing naturally in the language of cumulants, from the point of the spectral decomposition, X = k |R k λ k L k |, is far from being obvious and may explain the simplicity of formula (3.31).

Elliptic ensemble
As the first example of application of this formalism, we take the elliptic ensemble. Due to the fact that only the second cumulants do not vanish, the sum in (3.16) reduces to a single term and B is diagonal, B ell = diag(σ 2 τ, σ 2 , σ 2 , σ 2 τ ). However, the equations (3.15) do not decouple, because Green's functions are not diagonal in the non-holomorphic regime. Denoting for j = 1, 2 Green's function of the elliptic ensemble in the non-holomorphic regime (see appendix A), we find K, solving (3.15) Then we focus on the component K 11 11 and differentiate it twice, according to (3.4), obtaining This result was derived for the first time by Chalker and Mehlig [10]. 2 For the Ginibre Ensemble (σ = 1, τ = 0) it reduces to For completeness, we remark that the holomorphic part of the two point function, calculated from (3.18), reads h(z 1 ,z 2 ) = 4 (4.5)

Biunitarily invariant ensembles
We consider some examples where the two-point function can be easily calculated. This list is by no means exhaustive. In fact, biunitary invariance is preserved under addition and multiplication, thus one can easily generate new ensembles. We do not present results for products of the ensembles considered below, solely due to the fact that the expressions for O 2 (z 1 , z 2 ) become lengthy.
• Ginibre. As a cross-check of correctness of our formula, let us first consider the Ginibre ensemble. Its spectral density is uniform on the unit disk, therefore F (r) =
• Truncated Unitary [59]. Let us consider a (N + L) × (N + L) random unitary matrix with a pdf given by the Haar measure on U(N + L) and remove its last L rows and columns. The radial cdf of the remaining square matrix in the limit N, L → ∞, with κ = L N fixed, reads F (r) = κ r 2 1−r 2 for r < (1 + κ) −1/2 and 1 otherwise [60]. Therefore the two-point eigenvector function reads (4.8) • Spherical Ensemble. Consider the product Y = X 1 X −1 2 , where X 1 and X 2 are Ginibre matrices. Its radial cdf reads F (r) = r 2 1+r 2 and its spectrum is unbounded [61]. This ensemble is beyond the assumptions made for the derivation of (3.31). Nevertheless, motivated by the successful application of these methods for the one-point correlation function in this ensemble [29], we assume the correctness of our formulas and calculate the two-point function (4.9) • Product of two Ginibre. We consider a matrix Y = X 1 X 2 , where X 1 and X 2 are Ginibre matrices. The radial cdf of Y is F (r) = min(r, 1), thus (4.10)

Pseudohermitian matrix
Let us consider the product X = AB of two Hermitian matrices A, B. The product is not Hermitian, X † = BA = X, but if one of the matrices, let us say A, is positive definite, X is isospectral to the Hermitian matrix A 1/2 BA 1/2 , thus X, despite its non-Hermiticity, has a real spectrum. The diagonalising matrix is, however, not unitary, resulting in nonorthogonality of eigenvectors. Such matrices can be toy-models for more complicated physical system described by Hamiltonians which are not Hermitian but possesses parity-time (PT) symmetry [62]. The most interesting models have a parameter which controls how far the system is from breaking of the symmetry. At a critical value, called the exceptional point, two real eigenvalues coalesce and move away in the imaginary direction, spontaneously breaking the PT-symmetry.
As an example we consider the matrix X = (2+G 1 )(2+G 2 ), where G i 's are independent matrices drawn from the Gaussian Unitary Ensemble, the spectral density of which in the large N limit is the Wigner semicircle, ρ GUE (x) = 1 2π √ 4 − x 2 , supported on the interval [ −2, 2]. This model has an exceptional point at x = 0.
The components of the quaternionic R-transform of X read [63] (4.11) The other two components are given by the exchange of indices 1 ↔1. Inserting them into (2.22) and focusing only on the holomorpic solution (|w| = 0), we arrive at the equation for Green's function 4 (1 − g(z)) 2 + 1 g(z) = z. We choose a branch which gives the asymptotic behavior g(z) ∼ 1/z for large z. The spectrum is supported on a single interval [0, z + ], with z + = 1 2 (11 + 5 √ 5). The Green's function infinitely close to the spectrum reads where a = 1 + 10x + x 2 and r = 1 The imaginary part of Green's function yields the spectral density, calculated in [64]. The traced product of resolvents according to (3.18) satisfies the equation where g(z 1 ) and g(z 2 ) are the solutions of (4.13) with 1/z asymptotic behavior. The two-point function is calculated from (3.6) and its cross-sections are juxtaposed with the numerical simulations in figure 10, showing an excellent agreement.

Towards microscopic universality of eigenvectors
Random matrices show the phenomenon of universality at certain regions of the spectra. In the case of Hermitian ensembles, such universalities appear in the bulk (the so-called sine kernel) and at the edges of the spectra (Airy, Bessel, Pearcey, etc.). For a given generic Hermitian ensemble represented by N × N matrices H, one of the tools for investigating the existence of universalities are the multi-trace correlation functions (5.1) The subscript c denotes the connected part. Such objects were studied extensively using various techniques including loop equations [37], Coulomb gas analogy [65] and Feynman diagrams [52,53]. They were put into a formal mathematical formulation of the higher order freeness [66][67][68]. When the eigenvalues occupy a single interval, they obey the Ambjørn-Jurkiewicz-Makeenko universality [37]. The divergences of the double-trace correlation function signal the breakdown of the 1/N expansion and the need to resum the whole series and rescale its arguments. Different universal limits are manifested as different types of singularities.
A natural generalization of the two-point double-trace function to the non-Hermitian setting is the connected average of two copies of the electrostatic potential (2.5) introduced in [42], where Gaussian models were also considered. As the quaternionic Green's function, encoding all expectations of the traces can be obtained from the potential (see (2.10)), the function above generates all covariances of traces  Figure 11. Hierarchy of wheel diagrams contributing to the two-point double-trace correlation function (5.2), which in turn corresponds to the contraction of indices in figure 4 leading to L. The combinatorial factor 1/m corresponds to rotational symmetry and prevents overcounting the diagrams.
being a natural extension of the second order freeness for large non-Hermitian matrices.
As we mentioned earlier, the indices in the product of a resolvent can be contracted in two ways, see figure 4. One of them gives access to the eigenvector correlation function, while the second one yields F . More precisely, F (Q, P ) = TrL.
Since we consider connected expectation, we may write ln det(Q − X ) = ln det(1 − X Q −1 ) and use the identity ln det = Tr ln. Then, logarithms are expanded in power series, ln(1 + z) = ∞ k=1 z k k , which allows for convenient calculation of Feynman diagrams. Due to the presence of traces, the baselines from (X Q −1 ) k are now drawn as two concentric rings. 3 The dominant diagrams are the planar ones in which vertices and propagators are drawn between the two rings, but propagators connecting vertices do not encircle the inner ring (as in figure 4), see also [52,53]. The diagrams have an additional symmetry, namely rotating each ring leads to a new admissible diagram contributing equally. The resulting symmetry factors exactly cancel coefficients in the expansion of logarithms.
Each diagram can be decomposed into m segments in which X 's from two rings are connected through propagators and vertices. Segments are connected to each other through rings. As a result, each diagram looks like a wheel with m spokes. It turns out that the sum of all diagrams contributing to the spoke is exactly the rung, Γ, from the ladder diagrams in section 3.3. The X 's on ring, which are not part of a spoke can be connected with each other through propagators and vertices in any way, thus contributing to the Green's function. The general structure of such diagrams is presented in figure 11.
The wheel diagrams with m spokes have an additional symmetry, namely they can be rotated by an angle 2π/m. In order not to overcount the diagrams in the sum, we must include 1/m factor. Finally, we get

(5.4)
This means that the result, derived in [42] and used for deducing the existence of the edge universality for the spectral density [69], holds for the entire class of non-Hermitian models.

JHEP06(2018)152
The two-point single-trace correlation functions encoding correlations between eigenvectors also have their counterpart in the Hermitian case, but because of realness of the spectrum and orthogonality of eigenvectors it trivially reduces to the one-point Green's function thus not attracting attention. Eigenvectors of non-normal matrices are no longer trivial, making such correlation functions meaningful quantities. In the spirit of the above analysis, one is tempted to ask if we can probe hypothetical eigenvector universality using similar tools. We would like to stress that, even in the case of the simplest Ginibre ensemble, the direct analysis of the eigenvector correlation functions is very hard. Whereas the finite N expression for the one point function is known [10,23], the only known non-perturbative results for the calculation of the two-point eigenvector correlation function are given implicitly [9,10] as where the matrix h is pentadiagonal with entries given by There is, however, a different possibility of inferring the existence of universality. Spectra of non-normal matrices are intimately linked with the properties of their eigenvectors. The completeness relation While the right hand side is of order 1, the one-point correlator gives a contribution of order N , thus there has to be a counterterm from the integral. As the region of integration is in fact compact in the large N limit, the divergence can stem only from the region when w is close to z. The exact calculations in this regime are not accessible within the diagrammatic approach, but below we give a qualitative argument that the microscopic scaling is responsible for the cancellation of divergences.
In RMT the microsopic universality can be probed on the scale of the typical distance between eigenvalues. Demanding that in the disk of radius δz centered at z we expect one eigenvalue, leads us to the scaling where u ∼ 1. We notice that in all examples presented in section 4 the two-point function can be expressed as

JHEP06(2018)152
with the same behavior in the denominator. The microscoping scaling (5.9) inserted in the denominator produces a term (N ρ(z)) 2 , while the Jacobian of the change of variables reduces the power to one, giving the desired behavior in N . Moreover, by the explicit evaluation of derivatives in (3.31) and the application of de l'Hospital's rule twice we get for biunitarily invariant ensembles P (z, z) = O(z) ρ(z) , which cancels densities, eventually leaving only the one-point function, which produces the desired counterterm. We hypothesize that this phenomenon is universal across all non-Hermitian ensembles in the bulk.
Motivated by the ubiquitousness of the |z − w| −4 divergence in the bulk we state the conjecture that in generic non-Hermitian matrices with complex entries for all points in the bulk at which the spectral density does not develop singularities, there exists a microscopic limit (5.12) The function Φ was calculated in [9,10] by evaluating O 2 (0, z) from the exact result (5.6) and taking the scaling limit z = ω √ N . It is presented in figure 12a) and compared with the evaluation of the exact formula.
Interestingly, performing an analogous reasoning for the Ginibre ensemble (in which P (z, w) = 1 − zw) when z is at the edge of the spectrum leads us to a different conclusion. When two arguments get close, the two-point function diverges, but it does so as |z − w| −3 , because P also vanishes, reducing the exponent. This suggests that at the edge the twopoint function scales as N 3/2 instead of N 2 in the bulk. This stays in agreement with the sum rule (5.8), since O 1 at the edge scales as N 1/2 [23]. The limiting scaling function is not available to us, hence we evaluate numerically (5.6) and show the results in figure 12b). The divergence of the two-point function at the origin for the product of two Ginibre matrices (4.10) also suggests the existence of a different scaling there.

Summary
Using the methods of the quaternion formalism [70] for non-Hermitian random matrices, we have proposed the explicit calculational scheme for the two-point eigenvector correlation function (1.4). First, we have checked that our formalism reproduces all known examples in the literature, i.e. the complex Ginibre ensemble, an elliptic ensemble and the open chaotic scattering ensemble. Second, we considered two subclasses of non-normal random matrices: the pseudo-hermitian and the biunitarily invariant ensembles, which in the large N limit are described by the R-diagonal operators from free probability [71,72]. In both cases we got new results for the two-point eigenvector correlation functions. In the case of the bi-unitarily invariant ensembles, the two-point function O 2 (z, w) has a particularly simple form. It is expressed solely as a function of the radial cumulative distribution function F (r) and the one-point eigenvector correlation function O 1 (r). Recently, it was proven [29] that for biunitarily invariant ensmbles, O 1 (r) can be expressed in terms of F (r) exclusively, which can be viewed as an extension of the single ring (Haagerup-Larsen) theorem [54,56]. Combining this result with our formalism, we arrive at the conclusion, that the two-point eigenvector correlation function for general biunitarily invariant ensembles in the large N limit depends functionally solely on the spectral density. Such a situation resembles the Ambjørn-Jurkiewicz-Makeenko (known also as the Brezin-Zee) universality in the case of Hermitian random matrix models, where the two-point spectral Green's function depends solely on the one-point Green's function, irrespectively on the specific ensemble. Mathematical formulation of such a construction is known as the second order freeness [66]. We are therefore tempted to speculate that, by combining second order freeness and freeness with amalgamation [73], the notion of the non-orthogonality of eigenvectors can be extended into a broader context of operators in von Neumann algebras. Indeed, an equation similar to (3.15) has recently appeared in the description of fluctuations of Gaussian block matrices [74]. Moreover, the diagrammatic calculations of the traced product of resolvents resemble the partition structure diagrams introduced in [75].

JHEP06(2018)152
The similarity of our result to AJM (BZ) universality has further consequences. In the case of the ABJ (BZ) universality, the singular points of the correlation functions identify the regions of the spectra where microscopic universality takes place. This includes both the cases of the bulk and edge universality. We are therefore inclined to apply a similar argument to our result, searching for the microscopic eigenvector universalities. An additional argument for the microscopic universality comes from a constraint on eigenvector correlation function (5.8), as originallly noted by Walters and Starr. The sum rules originating from this constraint strongly suggest the universal form of the microscopic two-point eigenvector correlations in analogy to a similar phenomenon for the sum rules of Dirac Euclidean operators found by Leutwyler and Smilga [76]. The latter lead the Stony Brook group to the discovery of the universal Bessel kernels for chiral random matrix models [77,78].

JHEP06(2018)152
Our analysis, as well as explicit examples for the biunitarily invariant ensembles calculated in section 4.2, point at the generic shape of such universality, coming from the ubiquitous factor |z − w| 4 in the denominator. Explicitly, O 2 (z, w) = − 1 π 2 P (z,w) |z−w| 4 , where lim z→w P (z, w) = O 1 (z)/ρ(z) yields the Petermann factor. Such unique behavior is responsible for the crucial cancellation of the divergent terms in the leading order in N in the sum rules. The identification of this mechanism leads us to predict the existence of the universal microscopic scaling of the eigenvector correlation function Φ(|ω|). Such a limit was obtained in the special case of the Ginibre ensemble [9,10]. We conjecture that this universality extends to at least biunitarily invariant random ensembles.
Interestingly, the sum rule (5.8) leads also to interesting predictions at the edge. It is well known that correlations of eigenvalues of non-Hermitian matrices exhibit universal behavior at the edge, given by the error function. Our large N results for the eigenvector correlations show that the leading singularity weakens at the edge, |z − w| 4 → |z − w| 3 , leading to N 3/2 scaling of the two-point correlation function. The numerical evaluation of the implicit exact result (5.6) confirms this hypothesis, but the analytic form of the scaling function is not yet available, even in the case of the simplest, complex Ginibre ensemble.
Our results are only one step towards understanding the statistical properties of nonnormal random operators and give rise to new questions. The matrix of overlaps O ij is the simplest invariant object. It is natural to ask what kind of non-trivial higher order invariants can be built out of eigenvectors. This problem is even more cumbersome in the light of recent results [24,33], because the distribution of the diagonal overlap O ii is heavy tailed and some objects, for instance O 2 ii , do not exist. For the real Ginibre ensemble the situation is even more hopeless, since at the real axis the one-point function O 1 does not exist! While one expects the existence of certain correlation functions involving local averages of distinct eigenvectors, it is unclear whether their mathematical structure simplifies as it does for spectral statistics, which form determinantal point processes. Even though an event with two or more eigenvalues lying close to each other is unlikely to happen due to the eigenvalue repulsion, correlations between their eigenvectors do not decay, as can be seen from the microscopic scaling of O 2 . It is therefore very appealing to study microscopic eigenvector correlations involving more than two eigenvalues.
Although the real eigenvalues and corresponding eigenvectors of the real Ginibre ensemble are beyond the scope of perturbative techniques, we expect that the results for the two-point function remain unchanged for the eigenvectors associated with complex eigenvalues of the real Ginibre. Despite that the eigenvector overlaps are heavy-tailed, the traces of powers of X and its conjugate transpose are localized around their mean value [23]. Such a big cancellation is possible due to the sum rule originating from the completeness relation. Based on this fact, we expect that the formula for the traced product of resolvents (3.18) holds also for matrices with real entries.
The issue of a hypothetical microscopic eigenvector universality for generic non-Hermitian ensembles is also of primary importance, since unraveling the unknown microscopic eigenvector correlations may give hope in the case of notorious sign problems by giving an insight into the properties of the Dirac operator in Euclidean QCD at non-zero chemical potential.

JHEP06(2018)152
Note added. After completing this manuscript, we became aware of a recent work by Bourgade and Dubach [24], which tackles the issue of eigenvector correlations in the complex Ginibre ensemble in the bulk using different probabilistic techniques. They found the full probability of the diagonal overlap as an inverse gamma distribution and also studied the first two moments of the off-diagonal overlap. Moreover they proved that the result for the macroscopic two-point function (4.4) extends to mesoscopic scales.

JHEP06(2018)152
This is the holomorphic part, valid outside the spectrum and we have to choose the branch of the solution with a minus sign for correct asymptotic behavior at infinity g(z) ∼ 1/z. In the holomorphic domain, the off-diagonal elements of Green's function vanish, because the one-point eigenvector correlation function is trivially zero as there are no eigenvalues there.
Considering the non-trivial solution of (A.3) and inserting it into the equations for 11 and11 components, we obtain a system of two linear equations σ 2 τ G 11 + σ 2 G11 = z, .
The spectral density is calculated according to (2.3): ρ(z,z) = 1 π ∂zG 11 = 1 πσ 2 (1 − τ 2 ) . (A.6) One can also calculate G 11 and get the following formula for the one-point eigenvector correlation function from (2.9) (A.7) The boundary of the spectrum can be calculated in two ways: by requiring that the holomorphic and non-holomorphic solutions match at the boundary or by imposing vanishing of the one-point eigenvector correlation function. Both methods give which is the equation for the ellipse with semi-axes σ(1 + τ ) and σ(1 − τ ), hence the name of the ensemble.

B Quantum scattering ensemble
Let us see how the procedure for determining the rung of the ladder works in practice. We consider the quantum scattering ensemble [79] given by where H is a N × N complex matrix with Gaussian entries of zero mean and variance |H kl | 2 = N −1 δ kl and Γ = M a=1 V a (V a ) † . The components of N -dimensional vectors V a are complex Gaussians with variance V a kV b l = N −1 δ kl δ ab . The two-point eigenvector correlation function in the limit M, N → ∞ with M/N = m fixed (planar limit) was studied by Mehlig and Santer [36]. We show how this result can be rederived within this formalism in a simpler way.

JHEP06(2018)152
Γ is the complex Wishart matrix [80] multiplied by m. The planar cumulants of the Wishart matrix are stored in the Voiculecu's R-transform from free probability, which reads R Γ (z) = m 1−z . The considered matrix X is non-Hermitian, therefore we need its quaternionic R-transform. Using the embedding of the complex R-transform into the quaternionic structure [70], we get R Γ (Q) = m(1 2 − Q) −1 . The Gaussian matrix is a particular instance of the elliptic ensemble corresponding to τ = 1, therefore R H (Q) = Q. Further, Γ is rescaled by a complex number iγ. The quaternionic R-transform of such a rescaled matrix is obtained from the relation [70] R iγΓ (Q) = gR Γ (Qg), where g = diag(iγ, −iγ). As the R-transform is additive under addition of two matrices, we have R X (Q) = Q + mg(1 − Qg) −1 , which we then expand into a power series Then we perform the procedure with acting derivatives on the quaternionic R-transform and substituting the argument, as described in section 3.3. After summing up the resulting series, we get Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.