Eigenvector statistics of the product of Ginibre matrices

We develop a method to calculate left-right eigenvector correlations of the product of $m$ independent $N\times N$ complex Ginibre matrices. For illustration, we present explicit analytical results for the vector overlap for a couple of examples for small $m$ and $N$. We conjecture that the integrated overlap between left and right eigenvectors is given by the formula $O = 1 + (m/2)(N-1)$ and support this conjecture by analytical and numerical calculations. We derive an analytical expression for the limiting correlation density as $N\rightarrow \infty$ for the product of Ginibre matrices as well as for the product of elliptic matrices. In the latter case, we find that the correlation function is independent of the eccentricities of the elliptic laws.


I. INTRODUCTION
Products of random matrices have continuously attracted attention since the sixties [1][2][3][4][5]. They are of relevance in many fields of mathematics, physics and engineering including dynamical systems [2,6], disordered systems [7][8][9], statistical mechanics [10], quantum mechanics [11], quantum transport and mesoscopic systems [12,13], hidden Markov models [14], image processing [15], quantum chromodynamics [16], wireless telecommunication [17,18], quantitative finance [19][20][21] and many others [22]. Recently, an enormous progress has been made in the understanding of macroscopic  and microscopic  statistics of eigenvalues and singular values as well as of Lyapunov spectra for products of random matrices [67][68][69][70][71][72][73][74][75][76][77]. In contrast, not much has been learned about the eigenvector statistics of the products of random matrices so far. In this paper, we address this problem by considering a correlation function for eigenvectors of the product of Ginibre matrices. More precisely, we study the overlap between left and right eigenvectors for finite N and for N → ∞. In the first part of the paper, we adapt ideas developed in [78,79] to the product of random matrices by using the generalized Schur decomposition [45] for finite N , while in the second part we combine the generalized Green function method [81][82][83][84] with linearization (subordination) [10,28,37] to derive the limiting law for the overlap for N → ∞.

II. DEFINITIONS
Consider a diagonalizable matrix X over the field of complex numbers. Let {Λ α } be the eigenvalues of X. The corresponding left eigenvectors L α | and right eigenvectors |R α satisfy the relations Note that the Hermitian conjugate of the second equation has the form: X † |L α =Λ α |L α , where the symbol 'bar' denotes the complex conjugation of Λ α . The eigenvectors fulfill the bi-orthogonality and closure relations in the form The two relations are invariant with respect to the scale transformation with arbitrary non-zero coefficients c α 's. According to Refs. [78,79], an overlap of the left and right eigenvectors is defined in the following way By construction, the quantity O αβ is invariant with respect to the scale transformation given by Eq. (3) and consequently does not depend on the vector normalizations. If X is a random matrix, one defines averages over the ensemble where dµ(X) is the probability measure for the random matrix in question. The dependence of O αβ on X is suppressed in the notation. We use this notation throughout the paper also for other observables that depend on random matrices. The global diagonal overlap averaged over the ensemble is given by while the global off-diagonal one is expressed by the formula We are interested here in unitarily invariant random matrices for which the probability measure is invariant with respect to the similarity transformation X → U XU −1 , where U is a unitary matrix. In particular, this invariance implies that O αα = O 11 and O αβ = O 12 for any α and β. It follows that We can also define the local diagonal overlap density by the formula and the off-diagonal one by The symbol δ(z) denotes the Dirac delta function on the complex plane. Clearly, the diagonal global overlap is equal to the integrated overlap density given by Eq. (9), i.e.

III. PRODUCT OF GINIBRE MATRICES
Consider the product of m independent identically distributed N × N Ginibre random matrices [85] with complex entries. The probability measure factorizes and can be written as a product of measures for individual Ginibre matrices dµ(X) ≡ dµ(X 1 , X 2 , . . . , X m ) = dµ(X 1 )dµ(X 2 ) · · · dµ(X m ) , each of which is given by where σ is a scale parameter, and DX i = αβ dReX i,αβ dImX i,αβ . According to Eq. (9), the local diagonal overlap density can be calculated with respect to the measure dµ(X) in the following way where Λ α 's correspond to the eigenvalues of the product X (12). An analogous formula holds for the off-diagonal density. In the calculations we set σ = 1. One can easily transform the result to other values of σ (14) by using the formula which merely corresponds to the scale transformation of all Ginibre matrices X i −→ σX i in the product (12). Later, when discussing the limiting laws for N → ∞ we will choose σ = N −1/2 . This choice of the scale parameter σ will ensure the existence of the limiting eigenvalue density on a compact support being the unit disk in the complex plane.

IV. CALCULATIONS OF THE OVERLAP FOR FINITE N
In order to calculate the global left-right vector overlap, defined by Eq. (4), for the product of Ginibre matrices (12), we will change the parametrization of the matrices X i 's using the generalized Schur decomposition [45] for i = 1, . . . , m, where U i are unitary matrices from the unitary group U (N ), and τ i are upper triangular matrices of size N × N . We use a cyclic indexing U i ≡ U m+i , in particular U 0 ≡ U m . Sometimes it is convenient to express each τ i as a sum of a diagonal matrix λ i and a strictly upper triangular one t i , namely In this representation, the product X is unitarily equivalent to a matrix T , that is The matrix T has also an upper triangular form The diagonal elements of T are given by and the off-diagonal ones by Any instance of τ i,αα with two identical Greek indices can be replaced by λ i,α and of τ i,αβ with two different Greek indices by t i,αβ in the last formula. One can also express the integration measure in terms of U 's, λ's and t's. Since one is interested in invariant observables, the U 's can be integrated out. For the scale parameter σ = 1 one gets [45] dµ(λ, t) where the normalization factor Z is given by the formula and the Vandermonde determinant ∆ (Λ) for the product X = X 1 X 2 · · · X m has the form The square of the determinant in Eq. (23) comes from the Jacobian of the transformation (17). The next step is to express the observables in terms of t s and λ s. For example, to calculate the diagonal overlap density [cf. Eq. (15)], we have to find O 11 = O 11 (t, λ) and to integrate over t's and λ's with the Dirac delta constraint while for the global overlap O = dµ(λ, t)O 11 (t, λ). The measure dµ(λ, t) (23) factorizes dµ(λ, t) = dµ(λ)dµ(t). One can first integrate over t's. This is a Gaussian integral and can be easily performed. After this integration, only the dependence on λ's is left where dµ(t) is a normalized Gaussian measure equal to the t-dependent piece of dµ(λ, t) (23). The last step is to integrate over λ's with the measure given by Eq. (23) where as before Λ α 's stand for Λ α = λ 1,α λ 2,α · · · λ m,α . We will do this below. First we have to find the function O 11 (t, λ). This can be done as follows. We choose the basis in which the product matrix X is equal to T . Such a basis exists since the two matrices are unitarily equivalent. In this basis, the first right eigenvector |R 1 is represented as a column vector with '1' in the position 1 and zeros elsewhere: The vector is written here as transpose of a row vector to save space. Denote the elements of the first left eigenvector L 1 | = (B 1 , B 2 , . . .). The eigenvalue equation L 1 |T = L 1 |Λ 1 leads to the following recursion relation for B β 's [78,79] The recursion is initiated by B 1 = 1 as follows from the bi-orthogonality relation (2). One finds The element O 11 of the overlap matrix is related to B's as and B's depend on t's and λ's through T 's and Λ's. Combining Eqs. (30), (31) with Eq. (26) we obtain an explicit form of the integral over t's and λ's which can be done. We will give a couple of examples below.

V. EXAMPLES
Let us first illustrate the calculations for N = 2, m = 2 and σ = 1 -that is for the product of two 2 × 2 Ginibre matrices. Firstly, we express T 12 in terms of t's and λ's as follows This gives T 12 = λ 1,1 t 2,12 + t 1,12 λ 2,2 and Λ α = λ 1,α λ 2,α for α = 1, 2. Thus we have According to Eq. (27), the integration over t's leads to the following result Now we have to compute the integral over λ's given by Eq. (28), namely We first integrate over the λ's that do not appear in the Dirac delta, that is λ 1,2 and λ 2,2 . These integrals are in general of the Gaussian type combined with a power function, i.e. d 2 z|z| 2k exp (−|z| 2 ) = πk!. As a result of the integration, we obtain Now we integrate over λ 2,1 . We use the scaling property of the Dirac delta δ(a(z − z 0 )) = (1/|a| 2 )δ(z − z 0 ) to get The integral over λ 1,1 can be conveniently done in polar coordinates, yielding where K ν denotes the modified Bessel function of the second kind. The global overlap is The overlap density depends on the modulus |z|. It is convenient to represent this quantity as a radial function in the variable r = |z|, Clearly O rad (r)dr is equal to the overlap density integrated over the annulus r ≤ |z| ≤ r + dr. In our case we have O rad (r) = 2r(2 + r 2 )K 0 (2r) + 2r 2 K 1 (2r) .

VI. CONJECTURE
The calculations of the global density are slightly easier because there is no Dirac delta δ(z − Λ 1 ) in the integrand. They are particularly simple for N = 2. In this case and after inserting this into Eq. (33) and integrating the t's, one obtains Each integral over λ is either of the form |z| 2 exp (−|z| 2 )d 2 z = π or exp (−|z| 2 )d 2 z = π, so all together the integration over λ's gives the factor π 2m which cancels the pre-factor π −2m yielding Now, consider the case m = 1 for any N . This case was discussed in Ref. [79]. As follows from the discussion presented in this paper, one can cast the overlap into the form of the following multidimensional integral where . What remains to do is to compute this integral. We do this in Appendix A, where we show that the integral yields The results given by Eqs. (56) and (58) suggest that O grows linearly with m and N , hence it is tempting to conjecture that for any m and N the global overlap is given by the formula The result given by Eq. (53) is in agreement with this formula and Monte Carlo simulations fully corroborate this conjecture as shown in Fig. 4.

VII. LARGE N LIMIT
We now consider the limit N → ∞. We set the width parameter σ 2 = 1/N in the measure (14). The limit N → ∞ has to be taken carefully since we expect O N (z) to grow with N as it results from Eq. (59). In order to explicitly indicate the size dependence of O(z) on N here we exceptionally added the subscript N to O(z) = O N (z), which is implicit in the remaining part of the paper. It is convenient to define the growth rate of the overlap density as It depends on N but is expected to approach a N -independent function o(z): In the calculations, we shall use the method [28] that was previously employed to calculate the limiting eigenvalue density The method is based on the generalized Green function [81][82][83] For clarity, the symbol 'hat' is reserved for matrices with a superimposed block structure. By defining the block-trace Tr b as a matrix of traces of individual blocks The elements of this matrix are related to each other, g 22 (z) =ḡ 11 (z) and g 21 (z) = −ḡ 12 (z) [44], so we have In the large N limit, the eigenvalue density is related to the diagonal element [81][82][83] and the growth rate of the overlap to the off-diagonal one [84]   For large N , the leading contribution to the overlap grows linearly with N : Equations (67) and (68) are general and can be applied to any random matrix provided the Green function g(z) can be calculated. So the goal is now to calculate the Green function for the problem at hand. To this end, we use the planar diagrams enumeration technique [86][87][88].

VIII. DYSON-SCHWINGER EQUATIONS
The enumeration of planar Feynman diagrams is a method to derive the large N limit for matrix models [86][87][88]. The method is based on a field-theoretical representation of multidimensional integrals in terms of Feynman diagrams. One is interested in calculating the Green function where Q is a constant matrix and X is the random matrix that is averaged over. Matrix indices are denoted by A and B in the last equation. In this approach, the Green function plays the role of generating function for connected two-point Feynman diagrams. The contributions from non-planar diagrams are suppressed at least as 1/N in the large N limit, so for N → ∞ only planar diagrams survive in the counting. One can write a set of equations that relate the Green function G AB to a generating function Σ AB for one-line irreducible diagrams. Such equations are known in the field-theoretical literature as Dyson-Schwinger equations. Here, we are interested only in Gaussian random matrices. In this case, the Dyson-Schwinger equations assume a simple form in the planar limit N → ∞ [28] where P AB,CD represents the propagator The matrix Q AB and the propagator P AB,CD are inputs to be injected into these equations, while G AB and Σ AB are unknown functions to be determined for the given inputs. In other words, one has first to specify what Q and P are, and then, using these equations, one can find the Green function G, from there g and finally the eigenvalue density ρ(z) [cf. Eq. (67)] and the overlap growth rate o(z) [cf. Eq. (68)].

IX. SINGLE GINIBRE MATRIX
In this section, we review the calculations [82,84] for a single Ginibre matrix [85]. In the next section, we will then show how to generalize the method to the product of Ginibre matrices [28].
As mentioned before, first one has to identify the matrix Q and to calculate the propagator P AB,CD . The Green function (62) reads and where The symbol ⊗ denotes the Kronecker product. The blocks of the matrix X can be identified with the Ginibre matrix and its Hermitian conjugate: X 11 = X, X 22 = X † and X 12 = X 21 = 0, respectively. In order to calculate the propagator, we recall that the two-point correlations for the Ginibre matrix (14) with σ 2 = 1/N are and  2N × 2N . This block structure is also inherited by the propagators. Using the identification X 11 ↔ X, X 22 ↔ X † , along with Eq. (76) and (77), we see that the propagator factorizes into the inter-block part (in Greek indices) and intra-block part (in Latin indices) The only non-trivial elements of the inter-block part are p 11,22 = p 22,11 = 1. All remaining elements vanish: p αβ,γδ = 0. Since both the propagator (78) and the matrix Q AB = q αβ δ ab are proportional to the Kronecker deltas in Latin indices, this implies that the matrices G and Σ, being the solution of the Dyson-Schwinger equations (70), also are proportional to the Kronecker delta in the intra-block indices Alternatively, one can write G = g ⊗ 1 and Σ = σ ⊗ 1. Therefore, one can reduce the Dyson-Schwinger equation (70) to equations for inter-block elements (in Greek indices) In the second equation, we used that p 11,22 = p 22,11 = 1 and p αβ,γδ = 0 for other combinations of indices. The limit N → ∞ has already been taken in these equations, since they count contributions of planar diagrams. Now we can take the limit → 0 [cf. Eq. (65)]. This merely corresponds to setting = 0. Eliminating the {σ αβ }, we get Setting g = g 11 =ḡ 22 and γ = g 12 = −ḡ 21 we obtain The solution reads and g(z) =z, |γ(z)| = 1 − |z| 2 for |z| ≤ 1 .
The solution for γ(z) inside the unit circle is given up to the phase, but this is sufficient for our purposes since the correlations density o(z) given by Eq. (68) depends only on the modulus of γ(z). Using Eqs. (67) and (68), one eventually finds: and where χ D is an indicator function for the unit disk, χ D (z) = 1 for |z| ≤ 1 and χ D (z) = 0 for |z| > 1.

X. PRODUCT OF TWO GINIBRE MATRICES
In this section, we generalize the approach from the previous section to the product of two Ginibre matrices [28]. The integration measure for the product X = X 1 X 2 of independent Ginibre matrices X 1 and X 2 is the product of individual integration measures dµ(X 1 )dµ(X 2 ) given by Eq. (14). According to Eq. (76) the only non-vanishing two-point correlations are The Green function (62) for the product reads This form is difficult to handle because the product of Gaussian matrices X 1 X 2 is not Gaussian. One can however linearize the problem by considering a block matrix of dimensions 2N × 2N which is Gaussian. We call it root matrix because its square, reproduces two copies of the product, X 1 X 2 and X 2 X 1 . The two copies have identical eigenvalues. The Green function for the root matrix is which is actually a 4N × 4N block matrix and In this representation, the resolvent (92) has the standard form in which R is linear in the random matrices X's. Indexing blocks of R by R αβ , with α = 1, . . . , 4 and β = 1, . . . , 4, we have R 12 = X 1 , R 21 = X 2 , R 34 = X † 2 , R 43 = X † 1 . As follows from Eq. (87), the block R 12 is correlated with R 43 and R 21 with R 34 , so the propagator has the following non-zero elements, p 12,43 = p 43,12 = p 21,34 = p 34,21 = 1. All other elements of p αβ,γδ = 0. The situation is completely analogous to that discussed in the previous section, except that now the problem has dimensions 4 × 4 in inter-block indices. The intra-block correlations are the same as before -that is they are proportional to (1/N )δ ad δ bc -so the solution has the diagonal form proportional to the Kronecker delta in Latin indices (79). The Dyson-Schwinger equations (70) for the inter-block elements of the Green function of the root matrix read for → 0    g 11 g 12 g 13 g 14 g 21 g 22 g 23 g 24 g 31 g 32 g 33 g 34 g 41 g 42 g 43 g 44 In the second equation, we used the propagator structure: p 12,43 = p 43,12 = p 21,34 = p 34,21 = 1 and p αβ,γδ = 0 otherwise. Inserting {σ αβ } into the first equation, we get    g 11 g 12 g 13 g 14 g 21 g 22 g 23 g 24 g 31 g 32 g 33 g 34 g 41 g 42 g 43 g 44 It is convenient to solve this equation by defining matricesg andσ unitarily equivalent to g and σ:g = P gP −1 andσ = P σP −1 where The effect of the similarity transformation is equivalent to permutation of indices of the corresponding matrices: g αβ =g π(α)π(β) σ αβ =σ π(α)π(β) with π : (1, 2, 3, 4) → (1, 3, 2, 4). After this transformation, Eq. (98) is equivalent to The matrixg is a block matrix made of 2 × 2 blocks. The off-diagonal blocks are zero while the diagonal ones fulfill the following equations and The two equations admit only a symmetric solution being a solution of The last equation is exactly the same as for a single Ginibre matrix (81), so the solution eventually reads where g and γ are given by Eqs. (83) and (84). If we permute indices back to the original order g = P −1g P , we find    g 11 g 12 g 13 g 14 g 21 g 22 g 23 g 24 g 31 g 32 g 33 g 34 g 41 g 42 g 43 g 44 We see that the Green function for the root matrix consists of two identical blocks equal to the Green function of a single Ginibre matrix. In other words, the Green function of the root matrix behaves exactly as a pair of copies of the Green function of a single Ginibre matrix. The eigenvalue density and the growth rate of correlations between left and right eigenvectors of this matrix are given by Eqs. (85) and (86) as and Note that the size of the root matrix is 2N × 2N , so the leading term of the overlap behaves for large N as From these expressions, one may derive the corresponding expressions for R 2 , which are directly related to the product X 1 X 2 as follows from Eq. (90). The eigenvalues of R 2 are related to those of R as λ = λ 2 R , so one can find the densities by the change of variables z = w 2 : ρ(z)d 2 z = ρ R (w)d 2 w and O(z)d 2 z = O R (w)d 2 w. This gives and respectively. The result for the eigenvalue density ρ(z) was first found in [28]. The overlap O(z) is a new result. The product X 1 X 2 is of size N × N , so the growth rate is obtained by dividing O(z) by N , The radial profile is obtained from the last expression by setting r = |z| and multiplying the result by 2πr [cf. Eq. (41)]. This gives a triangle law where χ I is an indicator function for the interval [0, 1]: χ I (r) = 1 for r ∈ [0, 1] and χ I (r) = 0 otherwise. This prediction is compared to Monte Carlo data for N = 100 in Fig. 5.
As one can see in the figure, there are deviations from the limiting law for finite N . The radial profile drops to zero at the origin and develops a tail going beyond the support of the limiting profile for large r. We study the N -dependence of these effects in Fig. 6. We see that the gap at the origin closes in a way characteristic of the hard edge behavior, while the tail at the edge of the support gets shorter and falls off quicker as N increases. The behavior at the origin can be probably related to the microscopic behavior of the gap probabilities, which are driven by the Bessel kernel and were first studied in the context of QCD [89]. More generally, for the product of m matrices the behavior at the origin is controlled by the hypergeometric kernel [45]. In turn, the tail behavior at the soft edge is described by the error-function type of corrections [28,45].

XI. PRODUCT OF ELLIPTIC GAUSSIAN MATRICES
For completeness, we also consider the product of elliptic matrices defined by the measure [90] dµ(X) As before, we set σ 2 = 1/N and scale it with N while taking the limit N → ∞. The parameter κ belongs to the range [−1, 1]. It is related to the ellipse eccentricity. For κ = 0, (114) reproduces the Ginibre measure. Generically, the support of the eigenvalue density of matrices generated according to the measure given by Eq. (114) is elliptic. When κ approaches 1 (or −1), the support flattens and in the limit κ → ±1 gets completely squeezed to an interval of the real (or imaginary) axis. The corresponding matrix becomes Hermitian (or anti-Hermitian). The two-point correlations for the elliptic ensemble (114) are and Consider the product X = X 1 X 2 of two elliptic matrices X 1 and X 2 with different eccentricity parameters κ 1 and κ 2 . As in the previous section, we construct the root matrix (94), which is a 4N × 4N matrix. The propagator for the root matrix elements is where p αβ,γδ has now more nonzero elements. In addition to Inserting the {σ αβ } into Eq. (96) and permuting indices as in the previous section, we get This equation is much more complicated than that for the product of Ginibre matrices (105), because the two off-diagonal blocks on the right-hand side are non-zero. However, making the ansatz that the off-diagonal blocks of the solution vanish forces the two remaining blocks to satisfy the very same equation as for the product of Ginibre matrices (100), hence the solution is the same as before. This solution is independent of the eccentricity parameters κ 1 and κ 2 and moreover it is always spherically symmetric, even though the two matrices in the product are elliptic. To summarize, in the large N limit the eigenvalue density and the left-right eigenvector correlations for the product of two elliptic matrices are spherically symmetric (110) [28] and the eigenvector correlations are identical as for the product of Ginibre matrices (113). This prediction is compared to Monte Carlo data for N = 100 in Fig. 7. We see that it also follows the triangle law as for the product of Ginibre matrices. The finite N data exhibit however stronger finite size effects as compared to those for the product of two Ginibre matrices which manifest as a stronger deviation from the limiting density for small values of r. Compare Figs. 5 and 7. More generally, the limiting profile for N → ∞ is independent of κ 1 and κ 2 , while the finite size corrections do depend on the eccentricities. We checked numerically that the overlap density for the product of elliptic matrices is isotropic (circularly invariant).

XII. PRODUCT OF M GINIBRE MATRICES
We now proceed analogously as in Sec. X, where we discussed the product of two Ginibre matrices in the large N limit. The integration measure for the product X = X 1 X 2 · · · X m of m independent Ginibre matrices X 1 , X 2 , . . . , X m is the product dµ(X 1 )dµ(X 2 ) · · · dµ(X m ) of the individual integration measures given by Eq. (14). In turn, the two-point correlations are given by Eq. (76) for µ, ν = 1, . . . , m and a, b, c, d = 1, . . . , N . As in Sec. X, instead of directly applying the Green function technique to the product X 1 X 2 · · · X m , we apply it to the root matrix R being a block matrix of dimensions mN × mN The m-th power of the root matrix reproduces m cyclic copies of the product X 1 X 2 · · · X m , which all have identical eigenvalues. The Green function for the root matrix is a 2mN × 2mN block matrix and The resolvent given by Eq. (125) has the standard form with R being linear in X's. We index blocks of R by Greek letters R αβ , with α, β = 1, . . . , 2m. We have the following equivalence R α,[α]+1 ≡ X α and R m+[α]+1,m+α ≡ X † α for α = 1, . . . , m and [α] = α modulo m. All other blocks are zero. As follows from Eq. (122), we see that the only non-zero two-point correlations are for α = 1, . . . , m. Thus the propagator has the form and p αβ,γδ = 0 otherwise. The situation is analogous to that discussed in Sec. X, except that now there are 2m × 2m blocks. The intra-block correlations are the same as before (1/N )δ ad δ bc , so the solution is given as before as Kronecker product with the Kronecker delta in the intra-block indices G AB = g αβ δ ab [cf. Eq. (79)]. The first Dyson-Schwinger equation (70) The second Dyson-Schwinger equation (70) yields for α = 1, . . . , m, and σ αβ = 0 for all other elements of the matrix σ.
Any integer x can be decomposed uniquely as x = y + 2mk where y ∈ {1, 2, . . . , 2m} and k is an integer. The function (x) selects y from this decomposition. In particular (x) = x and (2m + 1) = 1, (2m + 2) = 2, (0) = 2m, (−1) = 2m − (134) The matrixg can be viewed as a block matrix made of 2 × 2 blocks. The off-diagonal blocks are zero and the diagonal ones fulfill the following equations for α = 1, . . . , m. Making the ansatz that the solution should be symmetric -that isg 2α−1,2α−1 = g,g 2α,2α =ḡ,g 2α−1,2α = γ andg 2α,2α−1 = −γ for all α = 1, . . . , m, the last equations reduce to a single one which is identical as that for a single Ginibre matrix (82). Hence, the solution for γ and g is given by Eqs. (83) and (84). This ansatz is equivalent to the one we used for m = 2 and merely means that the solution should not break the symmetry between different cyclic permutations of Ginibre matrices in the product. Inserting the solution intog we find where g and γ are given by Eqs. (83) and (84). Permuting indices back to the original order g = PgP −1 Hence, we see that the Green function of the root matrix behaves as m copies of the Green function of a single Ginibre matrix. The eigenvalue density and the growth rate of correlations between left and right eigenvectors of this matrix are identical as Eqs. (85) and (86), namely and The leading term of the overlap is therefore The eigenvalues λ of R m are related to those of R as λ = λ m R , so by changing variables as z = w m we can find the corresponding distributions for R m : ρ(z)d 2 z = ρ R (w)d 2 w and O(z)d 2 z = O R (w)d 2 w. This gives and respectively. Thus for large N the growth rate of the overlap for the product X 1 X 2 · · · X m is The radial profile defined by Eq. (41) is where as before χ D is the indicator function for the unit disk |z| ≤ 1 and χ I for the interval [0, 1]. While finalizing the manuscript, we learned that this result was derived independently in [91] with the aid of an extension of the Haagerup-Larsen theorem [92,93]. The integrated growth rate is which means that for large N the overlap grows as O ∼ mN/2 in agreement with Eq. (59). In Fig. 8, we plot the expression given by Eq. (145) for m = 4 and compare it to Monte Carlo data for N = 100. So far we have discussed the product of m Ginibre matrices. We could repeat the whole discussion from this section for the product of elliptic matrices (114) with arbitrary eccentricity parameters κ 1 , κ 2 , . . . , κ m . We would then arrive at an equation for g like Eq. (134) except that the matrix on the right hand side would now have non-zero non-diagonal 2 × 2 blocks. These blocks would be made of elements of non-diagonal blocks ofg multiplied in some way by κ's. Adopting the ansatz from Section XI that all off-diagonal blocks ofg are equal zero we would reduce this equation to Eq. (134) and get the same result as for the product of m Ginibre matrices.

XIII. CONCLUSIONS
In this paper, we have studied macroscopic and microscopic eigenvector statistics of the product of Ginibre matrices. We have developed analytical methods to calculate the left-right eigenvector overlap for finite N and in the limit N → ∞. The overlap is not only an interesting object from the mathematical point of view but is also of interest for physical problems. In the physics literature, it is known as Petermann factor and is for example used as a measure of non-orthogonality of cavity modes in chaotic scattering [94,95]. The off-diagonal overlap has been recently used as a sensitive indicator of non-orthogonality occurring in open systems due to perturbations resulting from shifts of resonance widths [96,97]. It plays also an important role in the description of Dysonian diffusion for non-Hermitian random matrices [98,99].
There are many open problems and potential generalizations of the studies presented in this paper. For example, one may try to extend the studies of the microscopic eigenvector statistics to products of truncated unitary matrices [100], which can also be mapped onto a determinantal point process [51] via generalized Schur decomposition [45]. A great challenge is to determine the microscopic eigenvalue and eigenvector statistics for products of elliptic matrices or to find any non-trivial solvable example of products of random matrices having non-spherical measures.
We have considered complex random matrices here. It would also be interesting to study overlaps for products of real and quaternionic matrices. They are much more challenging since in these cases the microscopic correlations are driven by Pfaffian point processes rather than determinantal ones. The real and quaternionic ensembles have additional scaling regimes near the real axis, which introduce an additional complication. Moreover, the Schur decomposition, which is at the heart of the method used in this paper, cannot be applied in a straightforward way to real matrices since generically they are not orthogonally similar to upper triangular ones. On the other hand, we believe that the limiting laws for N → ∞ are identical for real and complex ensembles since the underlying Dyson-Schwinger equations are identical in the planar limit (N → ∞).
Concerning the large N limit and macroscopic statistics, it would be interesting to generalize the calculations of the overlap to polynomials of random matrices [37,40,42] and to go beyond isotropic (R-diagonal) matrices [91,93], as well as to better understand the overlap in terms of the quaternionic formalism [44], and finally to calculate the off-diagonal elements of the overlap (7) using the Bethe-Salpeter equation [79].

XIV. ACKNOWLEDGMENTS
We would like to thank Romuald Janik for many interesting discussions. P.V. acknowledges the stimulating research environment provided by the EPSRC Centre for Doctoral Training in Cross-Disciplinary Approaches to Non-Equilibrium Systems (CANES, Grant No. EP/L015854/1).
Appendix A: Calculation of the integral (57) In this Appendix, we detail the calculation of the integral given by Eq. (57) where we have renamed the Vandermonde determinant on N complex variables as ∆ N (λ) for convenience. We can rewrite this as which can be more compactly expressed as using the complex version of the Andréief identity [101]. The integral over z yields I jk (λ) = π (|λ| 2 + 1)(k − 1)! + k! δ j,k − πk!λδ j−1,k − π(k − 1)!λδ j+1,k .