Characteristic polynomials of products of Wigner matrices: finite-N results and Lyapunov universality

We compute the average characteristic polynomial of the hermitised product of $M$ real or complex Wigner matrices of size $N\times N$ and the average of the characteristic polynomial of a product of $M$ such Wigner matrices times the characteristic polynomial of the conjugate matrix. Surprisingly, the results agree with that of the product of $M$ real or complex Ginibre matrices at finite-$N$, which have i.i.d. Gaussian entries. For the latter the average characteristic polynomial yields the orthogonal polynomial for the singular values of the product matrix, whereas the product of the two characteristic polynomials involves the kernel of complex eigenvalues. This extends the result of Forrester and Gamburd for one characteristic polynomial of a single random matrix and only depends on the first two moments. In the limit $M\to\infty$ at fixed $N$ we determine the locations of the zeros of a single characteristic polynomial, rescaled as Lyapunov exponents by taking the logarithm of the $M$th root. The position of the $j$th zero agrees asymptotically for large-$j$ with the position of the $j$th Lyapunov exponent for products of Gaussian random matrices, hinting at the universality of the latter.


Introduction and main results
Characteristic polynomials represent one of the central building blocks when studying the spectral statistics of random matrices. For example, in invariant ensembles the Heine-formula directly relates the expectation value of a single characteristic polynomial over a random matrix of size N × N to the orthogonal polynomial of degree N . Invariant ensembles represent determinantal point processes, and the corresponding kernel of orthogonal polynomials follows from the expectation value of two characteristic polynomials [31,9] at finite-N . This statement extends to non-Hermitian ensembles as well [6]. In applications of random matrices characteristic polynomials also play a key role, e.g. in comparison to moments and correlations of the Riemann ζ-function [21,23], or in the theory of strong interactions in the presence of baryon chemical potential [27].
One of the central questions in random matrix theory is that of universality, that is the independence of the distribution of matrix elements in asymptotic regimes such as the limit of large matrix size. Two main classes of deformations of the classical ensembles with independent Gaussian distribution of matrix elements exist: Wigner ensembles, where the independence is kept and invariance is dropped, allowing for more general distributions than Gaussian, and the invariant ensembles, where independence is dropped while keeping invariance under orthogonal or unitary transformations. This introduces a dependence among matrix elements, typically through a potential.
Given that invariant ensembles represent determinantal point processes at finite-N , the knowledge of the kernel at finite-N allows a direct asymptotic analysis of the marginals or k-point correlation functions of the matrix eigenvalues. Here, sophisticated techniques as the Riemann-Hilbert method have been developed. The universality of products and ratios of characteristic polynomials has been directly addressed as well, yielding a generating functional for the kernel, both for invariant [16] and Hermitian Wigner ensembles [17], see also [26,1] for recent work using supersymmetry. We refer to [28,7] for most concise expressions for averages of products and ratios of characteristic polynomials at finite-N , and to [18] for the supersymmetric perspective on that.
In (non-Hermitian) Wigner ensembles, however, such determinantal structures seem to be completely absent at finite matrix size. In consequence powerful probabilistic tools have been developed by several groups, in oder to prove universality in the various scaling regimes, for Hermitian and non-Hermitian random matrices, cf. [12,29], respectively and references therein. How can we understand this broad universality? Are there perhaps also objects that show a similar structure as for Gaussian ensembles at finite-N ? Indeed it was shown by Forrester and Gamburd [13], that the expectation value of a single characteristic polynomial of Hermitian Wigner matrices of size N × N agrees with that of the corresponding Gaussian ensemble. It is given by the Hermite polynomial for the Gaussian Unitary Ensemble (GUE) and by the Laguerre polynomial for the complex Wishart ensemble (also called chiral GUE or Laguerre unitary ensemble). The same polynomials are obtained for real symmetric Wigner matrices [13]. In this short article we will extend the list of such examples of an exact agreement at finite-N to products of M non-Hermitian Wigner matrices, both for singular values and complex eigenvalues of the product matrix. Furthermore, given the combinatorial meaning of averaged characteristic polynomials pointed out by Forrester and Gamburd in the above examples for Hermitian Wigner matrices, we expect such interpretations to exist for our examples as well.
Consider M independent Gaussian random matrices G 1 , . . . , G M of size N × N . Each matrix G k has independent matrix elements g (k) i,j with identical normal distribution, with zero mean and variance σ 2 m,n ] = δ i,m δ j,n δ k,l σ 2 k . The squared singular values of G k are called Wishart ensembles, whereas the complex eigenvalues of G k are called Ginibre ensembles. In [3] it was shown for complex matrices with unit variances σ k = 1 that the squared singular values of the product matrix G 1 · · · G M form a determinantal point process, representing an example for a polynomial ensemble. The corresponding kernel of biorthogonal functions was explicitly determined in [3] using Gram-Schmidt orthogonalisation 1 . As the Heine formula trivially extends to polynomial ensembles, the following holds for the orthogonal polynomials: For M = 1 this relation is well known to hold for a single Wishart ensemble, where the polynomial reduces to the Laguerre polynomial p normalisation. We find that the same relation extends to the product of real or complex independent non-Hermitian Wigner matrices, where we allow to choose a different variance for each matrix. i,j of every matrix X k are independent as well, having arbitrary real or complex distributions with zero mean and variance σ k > 0, i.e., (  (1.5) in terms of a Hermite polynomial H N (x) orthognal with respect to e −x 2 , in monic normalisation. It is tempting to expect that Theorem 1.1 implies similar combinatorial consequences as [13] for M = 1, in particular as the moment generating function of products of random matrices relates to Fuss-Catalan numbers, see e.g. [24]. The fact that the two equations in Theorem 1.1 agree extend the observed asymptotic commutativity of (rectangular) matrices G k for N → ∞ described first in [10], and later on for finite-N in [20] in a weak sense.
We turn to the complex eigenvalues and corresponding characteristic polynomials. Due to independence of matrices and matrix elements, for both products of Ginibre and non-Hermitian Wigner matrices the expectation value of a single characteristic polynomial is trivial, For products of independent complex respectively real Ginibre matrices it was shown in [4] respectively [14] that the complex eigenvalues of the product matrix G 1 · · · G M form a determinantal respectively Pfaffian point process, with a rotationally invariant weight function. Thus the monomials z N are the orthogonal polynomials respectively the even subset of skew orthogonal polynomials as well, albeit trivial ones. The determinantal point process is of orthogonal polynomial type. Thus the kernel is given by a sum over orthonormalised polynomials, containing nontrivial information about the weight through their (squared) norms h k . A similar statement holds for the kernel of skew orthogonal polynomials for complex eigenvalues of the Pfaffian point process [14]. In orthogonal polynomial ensembles it is known, both for real [31,9] and complex eigenvalues [6], that the kernel can be expressed by an averaged product of two characteristic polynomials. The same holds true for the kernel of skew orthogonal polynomials in Pfaffian process including the real Ginibre ensemble [5]. Combined with the result for the kernel K (M ) N of products of M independent complex Ginibre matrices [4], extended to non homogeneous variances in [2], we have the following statement For products of M real Ginibre matrices the expectation value (1.7) yields precisely the anti-symmetric kernel K [5,14]. We will show that the same expectation value is obtained for products of independent non-Hermitian Wigner matrices.
It is simple to understand why the above results do not easily extend to more products of characteristic polynomials: In all averages (1.3), (1.4) and (1.9) every matrix X k and its adjoint X * k appear exactly once. This also explains the absence of higher moments of these matrices. In principle, the agreement between Gaussian and non-Hermitian Wigner matrices at finite-N could thus be extended to the expectation of any polynomial that shares this property. Let us emphasise that the proofs of Theorem 1.1 and 1.2 are purely algebraic and constructive. They directly yield the explicit combinatorial result for non-Hermitian Wigner matrices, without recurring to independent calculations for Gaussian matrix elements. For simplicity we have restricted ourselves to square matrices, see e.g. [3] for the generalisation of (1.2) to products of rectangular complex Gaussian matrices. We expect that this agreement holds for products of rectangular non-Hermitian Wigner matrices as well.
Since the above identities between Gaussian and non-Hermitian Wigner ensembles already hold for finite-N , the universality of these expectations in various large-N limits is guaranteed. Let us emphasise, however, that this does not imply an identity for all k-point singular value or complex eigenvalue correlation functions at finite-N , as then (non-Hermitian) Wigner ensembles do not possess any determinantal or Pfaffian structure. For singular values of complex matrices, the derivation of the kernel within polynomial ensembles requires to evaluate the expected ratio of two characteristic polynomials, see [11] for more details. While for complex eigenvalues (1.7) indeed establishes the (skew-)kernel, the k-point correlation functions also depend on the weight function multiplying this (skew-)kernel, which contributes non-trivially in the large-N limit.
In order to derive a non-trivial universality statement for products of non-Hermitian Wigner matrices based on the above findings, we consider a growing number of factors, choosing M → ∞, while keeping N fixed. Consider the zeros of the average in increasing order. These are all non-negative as is shown in [24]. We wish to compare these zeros to the limiting Lyapunov exponents of the product matrix (X 1 · · · X M ) * (X 1 · · · X M ). They are defined in terms of the ordered non-negative eigenvalues (or squared singular values) λ (1.11) The Lyapunov exponents are obtained in the almost sure limit (1.12) We refer to [30] for the vast literature, including existence for Gaussian and other random matrices. In the same scaling as in (1.11) we obtain the following for the zeros. (1.13) In the case that all matrices X j are real or complex Ginibre matrices, labelled by β = 1, 2 respectively, the corresponding Lyapunov exponents are explicitly known [25,15] µ j = 1 2 Ψ (βj/2) + 1 2 log 2σ 2 /β , j = 1, . . . , N, (1.14) where Ψ denotes the Digamma function. For large j, applying the standard large argument asymptotic for the Digamma-function, we get Ψ(j) = log(j) + O (1/j) .
(1. 15) In consequence, for large orders j ≤ N , Lyapunov exponents and zeros of the averaged characteristic polynomial (which is universal) agree up to an error O(1/j). This strongly suggests that the jth Lyapunov exponents for products of non-Hermitian Wigner matrices also become universal.

Finite-N identities for non-Hermitian Wigner and Ginibre matrices
The proofs of Theorems 1.1 and 1.2 use the independence and simple linear algebra, including the Cauchy-Binet formula. For convenience of notation, let us denote the set of all subsets of {1, . . . , N } with exactly r elements by K r,N . For K = {k 1 < . . . < k r }, L = { 1 < . . . < r } ∈ K r,N and a matrix X of size N × N we write the determinant of the corresponding r × r sub-matrix as follows: Thus the sub-matrix on the right-hand side is obtained from X by choosing the rows with indices 1 ≤ k 1 < . . . < k r ≤ N and then the columns with indices 1 ≤ 1 < . . . < r ≤ N .
We begin by introducing the following Lemma about expectations of two determinants of different sub-matrices of equal size of the same matrix and its adjoint. Characteristic polynomials of products of non-Hermitian Wigner matrices Lemma 2.1. For K,K, L,L ∈ K r,N and j = 1 . . . , N we have .

(2.4)
First, we assume that K = L orK =L. In writing both determinants as sums according to Leibniz' rule and expanding the product, every summand is a product of entries of X j without repetition. That means, due to K = L orK =L there will be at least one mismatch in the first or second index pairs. Using the independence of all entries, every summand vanishes in expectation, the second case in (2.2). Now let us assume K = L andK =L. In this case we can write (2.4) as where we used that all summands with permutations π = σ vanish. This shows Eq.
Notice that on the left-hand side the sub-matrix of X * j X j is of size r × r, whereas the matrices multiplied on the right-hand side are of sizes r × N and N × r, respectively. For a set of indices 1 ≤ ν 1 < . . . < ν r ≤ N we write V = {ν 1 , . . . , ν r }, and an application of the Cauchy-Binet formula gives (2.5) By part one of the lemma, Eq. (2.2), we know that for every summand with set of indices K, V , and V, L we obtain a non-vanishing contribution r! σ 2r j for K = L only, and thus Now we can turn to the derivation of the explicit expression for the average of a single characteristic polynomial.
Proof of Theorem 1.1. We start by deriving the expression given in the first statement of Eq. (1.3). Therefore, we expand the characteristic polynomial into powers of x by expressing the corresponding coefficients in terms of the principal minors ECP 26 (2021), paper 30.
The minors of the product on the right-hand side can be expanded by means of an iterated application of the Cauchy-Binet formula (2.5). This way we obtain We reorder the factors in the product of the 2M determinants on the right-hand side by pairing matrices X * j and X j and take the expectation to obtain where we also used the independence of the matrices in order to distribute the expectation over pairwise matching factors. Now we are able to apply Lemma 2.1. Thus, we first see that a summand vanishes as soon as one of the conditions K (j) = K (2M +2−j) , j = 2, . . . , M , is not satisfied. Evaluating the expectations explicitly we obtain which now leads to the claimed expression.
To derive the second average characteristic polynomial in Eq. (1.4) we can proceed analogously. We first expand the polynomial in terms of principal minors, which in turn can be expanded by means of the Cauchy-Binet formula in the following way We can use the independence of the matrices X j again in order to distribute the expectations over the factors, and observe that by Lemma 2.1 we only have contributions from index sets with K (1) = K (2) = · · · = K (M ) . Hence, we obtain Let us turn to the average over two characteristic polynomials, related to the complex eigenvalues of the product matrix.
Proof of Theorem 1.2. We first deal with one of the determinants in Eq. (1.9), i.e., we have The minors of the product on the right-hand side can be expanded by means of an iterated application of the Cauchy-Binet formula (2.5). This way we obtain and similarly for det wI N − (X 1 · · · X M ) * . We thus have for their average using the independence of the matrices X 1 , . . . , X M . If ν = µ, then for every index j = 1, . . . , N −1 we find that in one of the two matrices K (j+1) K (j) X j and L (j) L (j+1) X * j there are entries which do not appear in the other matrix. Thus, the expectation of the determinant of these two matrices vanishes in this case, and using Lemma 2.1, we obtain,

Large-M asymptotic of the zeros and Lyapunov spectrum
For simplicity, we first study the case with unit variance σ k = 1 for k = 1, . . . , M . We are interested in the behaviour for large M and fixed dimensions N of the suitably rescaled zeros of (1.10) To this end, we study the behaviour of the quantities (z  This rescaling introduces many additional zeros in the complex plane, however, this happens in a regular way and later we will be interested in comparing the positive zeros only to the incremental Lyapunov exponents.
First, we study the asymptotic behaviour of the polynomials P   Proof. We note that the sequence |z| 2k k! , k = 0, . . . , N , is unimodal with a unique maximum This suggests to consider the index ν z , so that this sequence strictly increases up to the index ν z , and strictly decreases afterwards. First, let us consider |z| 2 ≤ 1 − , then we have as M → ∞, where C N is some positive constant depending on N only. Next, let us consider the case ν + ≤ |z| 2 ≤ ν + 1 − for some ν ∈ {1, . . . , N − 1}. Then the sequence |z| 2k k! , k = 0, . . . , N , attains its unique maximum at k = ν. Moreover, for k = ν we can estimate the following quotients for k < ν by as M → ∞, and for k > ν in a similar way by The sum can be estimated by The right-hand side equals Next, we show that every point of these circles indeed is a limit point of the zeros, from which we can deduce that the zeros converge weakly to the uniform distribution on the union of these circles. Due to the specific rescaling of the polynomials we have for all integer values of Hence, it is sufficient to show that every point 1, √ 2, . . . , √ N is a limit point of the zeros. To this end we study the behaviour of the polynomials P (M ) N in the neighbourhood of these points.

Proposition 3.2.
Let us consider a fixed ν ∈ {1, . . . , N } and a fixed N . Then we have for some q ∈ (0, 1) as M → ∞, uniformly in w on compact subsets of the complex plane.
Using the same estimates as in the proof of the statement (3.4) it is not difficult to see that the latter sum in fact is of order O q M , as M → ∞, uniformly in w on compact subsets of C, which gives the statement in (3.6).
From the statement (3.6) we can infer that we have