{'fraction_non_alphanumeric': 0.03937153800218273, 'fraction_numerical': 0.012272717912814283, 'mean_word_length': 4.247892803112168, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Free products of large random matrices - a short review of recent developments

We review methods to calculate eigenvalue distributions of products of large random matrices. We discuss a generalization of the law of free multiplication to non-Hermitian matrices and give a couple of examples illustrating how to use these methods in practice. In particular we calculate eigenvalue densities of products of Gaussian Hermitian and non-Hermitian matrices including combinations of GUE and Ginibre matrices.


I. INTRODUCTION
Products of random matrices arise in many areas of research and engineering [1]. In studies of dynamical systems one uses products of random matrices to derive limiting laws for Lyapunov exponents [2]; in information theory and wireless telecommunication -to calculate channel capacities for serial MIMO (multiple-input multiple-output) transmission [3]. With the help of products of random matrices one investigates unitary evolution [4], matrix diffusion [5] or phase transitions in quantum Yang-Mills theories [6]. One encounters products of random matrices in studies of quantum entanglement [7], financial engineering [8], and many other fields of research [9,10].
The problem of random matrix multiplication has attracted attention since the sixties [11,12]. Originally the studies were mostly concentrated on deriving limiting laws for products of infinitely many finite random matrices [10] with some exceptions where also limiting laws for infinite matrices were formulated [2]. Recently the focus has shifted to products of infinitely large matrices due to the discovery of a correspondence between infinitely large invariant random matrices and free random variables [13,14] which have been studied in the framework of free probability [15] -a new branch of probability theory that has been intensively developed since the eighties [16,17]. Many results derived in this theory have been applied to large random matrices in the limit of infinite matrix size N → ∞.
In this paper we review some of these results, and give practical recipes to calculate eigenvalue densities of products of large random matrices. Instead of giving formal proofs we gradually build up intuition passing from matrix addition through multiplication of Hermitian random matrices to multiplication of non-Hermitian matrices. Free addition of large random matrices is a counterpart of addition of independent random variables in classical probability. We start with addition because it is much simpler than multiplication but it already has essential ingredients needed to understand a special role of invariant random matrices played in the formulation of limiting laws for eigenvalue distributions in the large N limit. The key feature of invariant random matrices is that their eigenvectors are not correlated with the eigenvalues, and that they are uniformly distributed. The last statement means that the unitary matrix constructed from the normalized eigenvectors is uniformly distributed on the unitary group. In effect, randomness of such matrices is entirely shaped by the probability distribution of the eigenvalues. The sum of independent invariant random matrices is an invariant random matrix. In the limit N → ∞ the eigenvalue density of the sum depends only on eigenvalue densities of the matrices added. This law is very similar to the law of addition in classical probability where the probability distribution of the sum of independent random variables can be calculated from the probability distributions of these variables. Similarities to classical probability go far beyond that. One can write a dictionary between classical and free probability. In this dictionary, random variable corresponds to invariant random matrix for N → ∞, probability density function to eigenvalue density, characteristic function to Green's function, cumulant generating function (logarithm of the characteristic function) to R transform, Mellin transform of the probability density function to S transform, independence to freeness, etc. Using this correspondence one can find a full analogy between free addition and free multiplication of large invariant matrices and addition and multiplication of independent random variables.
It is worth emphasizing that the laws of free addition and free multiplication can be independently derived using field theoretical methods, including perturbation theory of invariant matrix models [18,19], without referring to free probability. Perturbation theory is a way of graphical enumeration of planar Feynman diagrams that arise in the large N limit [20]. A great advantage of the field theoretical approach is that it also works in the case of non-Hermitian matrices and it allows to derive a non-Hermitian version of the multiplication law [21]. To the best of our knowledge this result is beyond the reach of free probability so far. We shortly recall this result here.
The paper is organized as follows. In section II we discuss eigenvalues of a sum of random matrices. In particular we define a slightly modified version of matrix addition -such that it is independent of the eigenvectors of the matrices added. This brings us to the class of invariant random matrices that are discussed in section III. In section IV we recall the definition of the Green's function, the moment generating function and the R transform. These functions are used in the derivation of limiting laws for invariant matrices in the large N limit. We employ these functions in the formulation of the law of free addition. In section V we discuss the law of free multiplication for invariant Hermitian matrices. In section VI we discuss consequences of this law for multiplication of infinitely large isotropic non-Hermitian matrices. Finally in section VII we show how to calculate eigenvalue densities for products of large independent non-Hermitian matrices. We illustrate calculations with examples. We conclude the paper with a short summary in section VIII.

II. EIGENVALUES OF THE SUM OF RANDOM HERMITIAN MATRICES
Let a and b be N × N Hermitian matrices. We are interested in the eigenvalues of the sum a + b. To calculate them it is not sufficient to know the eigenvalues of a and b. One would also need to know their eigenvectors. Actually it would be enough to know relative positions of the eigenvectors of a and b: for example the decomposition of beigenvectors in the eigenbasis of a. One can, however, define a non-standard addition of matrices that depends only on the eigenvalues of individual matrices in the sum and not on their eigenvectors. This new addition is defined in a probabilistic way as a sum where u, v are independent Haar unitary matrices and u † denotes the Hermitian conjugate of u. Let us recall that Haar unitary matrix is a random matrix distributed according to the uniform measure on the unitary group U (N ). The matrix a b is not a concrete matrix but it is a random matrix or equivalently an ensemble of matrices with a certain probability measure. The eigenvalues of a b are random variables. One is interested in the joint probability distribution of these random variables. The marginal distribution of the joint probability distribution is referred to as the eigenvalue distribution. One should note that the transformations a → u † au and b → v † bv preserve the eigenvalues of the individual terms, but uniformly randomize their eigenvectors. One can therefore expect that the eigenvalue distribution of a b depends only on eigenvalues of a and b. Moreover, since the sum u † au + v † bv can be written as u † (a + w † bw)u, where w = vu −1 is also a Haar unitary matrix, the eigenvalue distribution of the matrix a + w † bw averaged over w is the same as the eigenvalue density of a b. In other words it is sufficient to randomize the eigenvectors of one matrix against the other. One can consider a more general version of the problem in which a and b are independent random matrices instead of being deterministic ones. In this case the question is whether one can determine the eigenvalue distribution of a b from the eigenvalue distributions of individual matrices a and b. The answer is affirmative for N → ∞. In this limit the eigenvalue density ρ a b (λ) of the sum a b can indeed be determined solely from the eigenvalue densities ρ a (λ) and ρ b (λ). For finite N the answer is more complicated since it depends also on eigenvalue correlations. The correlations can be neglected only in the limit N → ∞. The addition (1) is called free addition when N → ∞ and a and b are independent. In the next section we shall see that the addition (1) is equivalent to the standard matrix addition a + b for invariant random matrices.

III. INVARIANT HERMITIAN MATRICES AND FREE MATRICES
Random matrices are formally defined by matrix ensembles equipped with probability measures. Let us denote the measure of an N × N Hermitian random matrix x by µ(x). Random matrix x is called invariant if the measure is invariant under the transformation x → U † xU , where U is an arbitrary unitary matrix: µ(x) = µ(U † xU ). Let us give two standard examples.
One can construct an invariant random matrix from a diagonal random matrix Λ = diag(λ 1 , . . . , λ N ) and Haar unitary random matrix u a = u † Λu . ( The invariance of the measure for a: µ(a) = µ(U † aU ) is a direct consequence of the invariance of the Haar measure µ Haar (u) = µ Haar (uU ) under the right multiplication by a unitary matrix U . One can choose the diagonal elements λ i 's of Λ to be independent, identically distributed real random variables. In such a case the joint eigenvalue density factorizes and ρ N (λ 1 , . . . , λ N ) = p(λ 1 ) . . . p(λ N ) where p(λ) is the probability density function for the random variables representing diagonal elements of Λ. The factorization of the joint eigenvalue density reflects the independence of the eigenvalues. The eigenvalue density is ρ(λ) = p(λ).
In physical applications one frequently encounters invariant measures of the form [22,23] dµ(a) ∝ e −N TrV (a) da , where V (a) is a polynomial or a power series in a and da = 1≤i≤N da ii 1≤i<j≤N dRea ij dIma ij , denotes flat measure for N × N Hermitian matrices. The proportionality sign in (3) means that the right hand side is displayed up to a constant. This implicit constant ensures the probabilistic normalization dµ(a) = 1. The joint probability eigenvalue distribution is given by a well-known Dyson formula [22,24] which can be interpreted as a partition function that describes a gas of N particles in one dimension. We again skipped a normalization constant that ensures R N ρ N (λ 1 , . . . , λ N )dλ 1 . . . dλ N = 1. In this picture the i-th eigenvalue λ i corresponds to the position of the i-th particle in one-dimensional space. Particles are in the potential V (λ). They repel each other with a logarithmic potential depending of the distance |λ i − λ j |. The repulsion makes the eigenvalues interact. As a result they are correlated. The eigenvalue density can be calculated by integrating all but one eigenvalue from the joint eigenvalue distribution ρ(λ) = ρ(λ, λ 2 , . . . , λ N )dλ 2 . . . dλ N . In the limit N → ∞ one can relate the eigenvalue density to the potential by an integral equation [18] V obtained by applying the saddle point approximation to (4). P.V. stands for principal value. One can use this equation to calculate the eigenvalue density for the given potential or to find the potential corresponding to a given eigenvalue density. In particular one can choose the potential to reproduce the same eigenvalue density as in the first ensemble (2) where the eigenvalue density is just equal to the probability density function for the diagonal elements of Λ. The two random matrix ensembles (2) and (3) have completely different statistics of eigenvalues even if they have identical eigenvalue densities. Let us come back to the general discussion of invariant random matrices. The standard addition a + b of invariant random matrices is equivalent to free addition a b (1). This follows from the invariance of the measure, which implies that the matrix a is identically distributed as u † au. The same holds for b and v † bv. For a and b being invariant random matrices, the two sums a + b and a b have identical probability measures. They represent the same random matrix ensemble. Therefore from here on we shall not distinguish between a + b and a b while discussing invariant random matrices. The sum of invariant random matrices is also an invariant random matrix. An important property of the sum a + b of invariant matrices is that the eigenvalue density depends in the limit N → ∞ solely on the eigenvalue densities of a and b, as follows from the discussion in section II. Independent invariant random matrices for N → ∞ are called free. There is one-to-one correspondence between them and the free random variables known from free probability [15].
To summarize this section, the standard addition of invariant random matrices is equivalent to free addition (1). In the limit N → ∞ independent invariant random matrices are called free and the eigenvalue density ρ a+b (λ) of the free sum can be calculated from eigenvalues ρ a (λ) and ρ b (λ) of a and b. The resulting density is independent of the detailed statistics of eigenvalues. For example, it does not depend on whether eigenvalues are correlated or not or whether they are generated using the ensemble of type (2), (3) or any other. The only thing that matters is the eigenvalue densities of the individual matrices in the sum. In the next section we briefly summarize the algorithm to calculate ρ a+b (λ).

IV. ADDITION OF FREE MATRICES
For N → ∞ the eigenvalue density of an invariant random matrix a can be described by an eigenvalue density ρ a (λ) having support I on the real axis. Typically I is a finite interval. In calculations it is convenient to use the Stieltjes transform of the density that is a complex function defined on the complex plane outside I. This function is often called Green's function. The density can be reconstructed from the Green's function G a (z) by calculating the imaginary part of the Green's function close to the real axis The power expansion G a (z) at infinity generates moments of the eigenvalue distribution These moments are equivalent to the trace moments calculated with respect to the probability measure µ(a) for the entire matrix a. The zeroth moment is fixed m a,0 = 1 by the normalization of the measure. One often defines another moment generating function as a power series at zero It is related to the Green's function The eigenvalue density ρ a (λ), the moment generating function φ a (z), and the Green's function G a (z) contain roughly speaking the same information about the eigenvalue distribution of the invariant matrix a and one can reproduce these functions from each other. The idea is to calculate these functions for a + b when the corresponding functions for a and b are given. Let us first calculate the moments of a + b The measure factorizes, since a and b are independent by assumption. For n = 1, 2, . . . the last equation gives In the second equation we replaced 1 N Trab by m a,1 m b,1 . The equality 1  (14) can be written in a compact way using free cumulants which are certain (specific) combinations of moments [15]: for n = 1, 2, . . .. The concept of free cumulants is analogous to the standard cumulants in classical probability. Let us recall that cumulants for a real random variable are generated as coefficients of the power series expansion of the logarithm of the characteristic function. The logarithm of the characteristic function for a sum of independent random variables is equal to the sum of logarithms of the characteristic functions of these random variables. Therefore cumulants for the sum of independent variables are additive, exactly as free cumulants (15). Free cumulants are, however, generated in a different way. The generating function for free cumulants is called R transform [15,16]. Free cumulants are equal to coefficients of the power series expansion Note that there is a shift between the power of z and the index of the corresponding cumulant. The R transform is related to the Green's function which generates moments of the eigenvalue distribution. The relation reads .
We give it without a proof [15]. One can find a diagrammatic interpretation of this equation [21] using field theoretical methods for planar graphs enumeration [18,19]. Using equation (17) one can find G(z) if R(z) is known. One can also invert equation (17) G R(z) This form is useful when one wants to determine the R transform for the given Green's function. The R transform for the sum of free matrices a + b is a sum of R transforms for a and b [15,16] exactly as the logarithm of the characteristic function for the sum of independent real random variables in classical probability. The additivity of R transform reflects the invariance and independence of the measures for a and b in (13). For N → ∞ these two properties correspond to freeness. Explicit relations between free cumulants κ n and moments m n can be found by inserting the power series R(z) = κ 1 + κ 2 z + . . . and G(z) = 1/z + m 1 /z 2 + . . . into equation (18). The first three free cumulants are identical as the standard ones The R transform was first introduced in free probability by Voiculescu [16]. Later it was observed that free cumulants are related to a combinatorial problem of non-crossing lines on a plane [14,25]. For this reason free cumulants are sometimes called non-crossing cumulants. For completeness we mention that equation (18) can be deduced [21,26] using field theoretical methods of planar graph enumeration [18]. The Green's function and the R transform correspond to sums over certain classes of diagrams, and equations like (17) -to Dyson-Schwinger equations [18]. The combinatorics for planar Feynman diagrams for matrix models is equivalent to the combinatorial approach to freeness [14,25].
We close this section by giving a standard example. Recall that in classical probability the normal random variable with the probability density function p(x) = 1 √ 2π e −x 2 /2 has all cumulants equal zero except κ 2 = 1. What is the corresponding random matrix eigenvalue density in free probability that has all free cumulants equal zero except the second one κ 2 = 1? We want to find the eigenvalue density corresponding to the R transform equal R(z) = z (16). Using (17) we first find an equation for the Green's function G(z) = 1/(z − G(z)). It can be solved for G(z): . We select the branch of the square root to match the asymptotic behavior G(z) ∼ 1/z for large z. The eigenvalue density is (7) and zero otherwise. It is the Wigner semicircle law. Odd moments of this distribution are equal zero while even moments m 2k = 2 −2 dλρ(λ)λ 2k = 1 k+1 2k k are equal to the Catalan numbers. One can use equation (5) to find that the corresponding potential is quadratic V (x) = x 2 /2, so the invariant measure (3) is Gaussian dµ(a) ∝ e −N Tra 2 /2 . Let us slightly modify the R transform and consider R(z) = α + σz. Without repeating the calculations one can easily argue on general grounds that the corresponding eigenvalue density is shifted as compared to the previous case by the first cumulant κ 1 = α and rescaled by the second one κ 2 = σ leading to ρ(λ) = 1 Let us now calculate the eigenvalue density of the sum of two independent large (N → ∞) invariant matrices a and b that have eigenvalue densities given by the Wigner semicircle law with different shift and scale parameters. Denote these parameters by α a , σ a for a and α b , σ b for b, respectively. The R transforms are R a (z) = α a + σ a z and R b (z) = α b + σ b z. To find the eigenvalue density of the sum c = a + b we first calculate R transform that is equal to R c (z) = α c + σ c z with α c = α a + α b and σ c = σ a + σ b . The R transform is exactly of the same form as for a and b. Therefore the eigenvalue density for c is also given by the Wigner semicircle law. The situation is analogous to that in classical probability where the sum of two random Gaussian variables is known to be a Gaussian random variable.
One says that the semicircle law is stable with respect to free addition in the same way as the Gaussian law is stable with respect to addition in classical probability. Actually there exist a bijection between stable laws in classical and free probability that allows one to classify all stable laws in free probability [27]. Just as a curiosity we mention that the only distribution that is stable with respect to addition (random variables) and free addition (random matrices) is the Cauchy distribution.

V. PRODUCTS OF LARGE INVARIANT MATRICES
As we have seen the law of addition for independent large invariant matrices can be concisely expressed in terms of the R transform (19). Before discussing the law of multiplication it is convenient to recall the corresponding law in classical probability. Let x and y be independent non-negative real random variables. What is the distribution of the random variable xy? A way to derive this distribution is via the Mellin transform of the probability density function. The Mellin transform M xy (s) for the product xy is the product of Mellin transforms M xy (s) = M x (s)M y (s) for x and y. Clearly this multiplication law is analogous to the addition law discussed in the previous section (19). It remains to find the corresponding law for products of large invariant matrices. Let a and b be independent invariant matrices. For the moment we also assume a and b to be positive semi-definite. Later we will weaken this restriction. The first trivial observation is that even if a and b are Hermitian, the product ab is generically non-Hermitian, so its eigenvalues are complex. Therefore, it is convenient to introduce a slightly modified version of the multiplication [28] We use dot in the notation to distinguish a · b from the standard matrix multiplication ab. Square root of a positive semi-definite matrix is defined as The product a · b of Hermitian matrices is Hermitian. The product a · b of invariant random matrices is an invariant random matrix. Moreover, the probability measure µ a·b for the product a · b is identical as for b · a: µ a·b = µ b·a so this product is in a sense commutative. One should also note that the trace moments for the product a · b (20) are identical as for the standard product: 1 N Tr(a · b) k = 1 N Tr(ab) k . The law of free multiplication was formulated by Voiculescu in free probability [15,17]. Later it was shown that it applies also to multiplication (20) of invariant random matrices in the limit N → ∞ [28]. The multiplication law is expressed in terms of S transform where S a (z) is a complex function. The S transform S a (z) can be constructed from the moment generating function φ a (z) (11) or actually the inverse function of φ a that is usually denoted by χ a : χ a (φ a (z)) = φ a (χ a (z)) = z. The relation reads [17] S a (z) = z + 1 z χ a (z) .
The S transform can be viewed as a series in z that is obtained by inserting the inverse series in place of φ a (z) in equation (22). The coefficients of this series are related to the moments m a,k . The coefficients are singular for m a,1 = 0, in which case the S transform is ill-defined. The S transform is related to the R transform as follows [15,21] S(z) = 1 R(zS(z)) .
If one knows the R transform then the last equation can be used to calculate the S transform. One can invert the last equation [21] R(z) = 1 S(zR(z)) , which allows one to calculate the R transform if the S transform is known. Using these relations one can show that the multiplication law (21) can be alternatively written in the R transform form [21]: This set of equations looks more complicated than (21) but in some situations it is more advantageous. For example, it can be used when the first moments vanish, that is when the S transforms are ill-defined. Moreover, this form can be generalized to the case of non-Hermitian matrices, as we discuss in section VII. The two forms of the multiplication law are in a sense complementary. The original S transform form of (21) has a very simple and appealing structure that reflects important properties of free products, for example that the multiplication a · b (20) of invariant random matrices is commutative and associative: S (a·b)·c (z) = S c·(a·b) (z). This follows from the fact that multiplication on the right hand side of (21) is just multiplication of complex numbers which is commutative and associative. Let us give an example to illustrate how these formal relations work in practice. Consider the product a · b of two independent identical Wishart matrices [29]. The Wishart matrix can be thought of as a square of a Gaussian matrix from the GUE ensemble [22]. The matrices a and b have the eigenvalue densities ρ a (λ) = ρ b (λ) = ρ(λ) for λ ∈ [0, 4] and ρ(λ) = 0 for λ outside this interval. The Green's function (6) corresponding to this density is and the moment-generating function (12) φ The inverse function of φ fulfills the equation z = 1 2χ(z) 1 − 1 − 4χ(z) − 1 which can be solved for χ(z) Thus the S transform is (22) We now want to find the S transform for the product S a·b (z) from S a (z) = S b (z) = S(z). Using the multiplication law (21) we obtain We now proceed to find the eigenvalue density ρ a·b (λ) that corresponds to the S transform given above. First we write the equation for the inverse function of the moment generating function χ a·b (z) = z/(1 + z) 3 using (22). We invert it for φ a·b (z): It is a cubic equation for φ a·b (z). Replacing φ a·b (z) with the Green's function (12) we have This equation has three solutions given by the Cardano formulas. We select the one which asymptotically behaves as 1/z for large z. Using (7) we find the density ρ a·b (λ) = 2 1/3 3 1/2 12π for λ ∈ 0, 27 4 and zero otherwise. We close this section with a few general comments. The multiplication law (21) or (26) can be derived using solely perturbation expansion for invariant matrix models in the limit N → ∞ without referring to free probability. We refer the interested reader to [21]. The second comment concerns the multiplication law for matrices which are not positive semi-definite. In the discussion above we restricted the class of matrices to positive semi-definite ones because we wanted to use √ a in the definition of the product a · b (20). For (invariant) Hermitian matrices the product a · b is (invariant) Hermitian. One can also consider the multiplication law for the standard product ab of Hermitian matrices. In this case a and b do not have to be positive semi-definite. The matrix ab is non-Hermitian and thus has complex eigenvalues. Its trace moments 1 N Tr(ab) k are, however, real. Since the S transform depends only on the moments, all calculations are identical as discussed before for a · b except that the final result ρ ab (λ) cannot be interpreted as an eigenvalue density function, but merely as a function to calculate the trace moments of the matrix ab: The eigenvalue density has in this case generically a two-dimensional support on the complex plane. For example, the product of two independent Hermitian matrices from the GUE ensemble is a non-Hermitian matrix. Each GUE matrix has real eigenvalues with the eigenvalue distribution given by the Wigner semicircle law. The R transform R a (z) = R b (z) = z (we consider standardized GUE with the mean zero and unit variance). Using the R transform form of the multiplication law (26) one finds that R ab (z) = 0 and thus G ab (z) = 1/z. This means that all the moments m ab,n = 0. It does not mean, however, that eigenvalue density of the product is zero. To calculate the eigenvalue density in this case one has to use methods that allow one to treat non-Hermitian matrices. We discuss them in the subsequent sections. Here we only mention that the eigenvalue density for the product is where z is a complex variable. This result was first derived in [30] with the help of field theoretical methods. Note that the result is invariant under rotations around the origin of the complex plane.

VI. PRODUCT OF ISOTROPIC RANDOM MATRICES
The extension from Hermitian to non-Hermitian random matrices is similar to the extension from real to complex random variables in classical probability. In this section we employ this similarity to define isotropic random matrices [31] being probably the simplest extension of invariant random matrices to non-Hermitian ones.
Consider a complex random variable with a probability measure that is invariant under rotation around the origin of the complex plane µ(z) = µ(ze iα ). The probability density is a function of the distance from the origin dµ(z) = d 2 zρ (|z|). Such a random variable z = re iφ can be constructed as a product of a real non-negative random variable r and a uniformly distributed random phase φ ∈ [0, 2π). When one knows the probability density function for the radial variable r, one can easily derive the probability density for z on the complex plane.
Isotropic random matrices are constructed in a way that imitates this construction. A random matrix A is called isotropic if the probability measure µ(A) is invariant under the left and right multiplication by unitary matrices U, V : µ(A) = µ(U † AV ). In Math and Engineering community such random matrices are called Bi-Unitarily Invariant (BUI) and in general they can be rectangular. Here we concentrate on square matrices and use the name isotropic in reference to them.
The simplest way of constructing isotropic random matrices is to build them from real semi-positive diagonal random matrices a and Haar unitary matrices u, v independent of a (free of a for N → ∞) Clearly, the diagonal matrix a corresponds to the radial variable r of z = re iφ and the Haar unitary matrices v and u to the uniform phase φ. There are two of them to ensure the invariance under multiplication by a unitary matrix on both sides. Alternatively one could consider a one-sided invariance, in which case the last equation would take the form A = au. Most of the results given below hold in either case. Another way to construct such matrices is to define them by a probability measure of the form [32] dµ where dA = ij dImA ij dReA ij is flat measure for N × N complex matrices and V (x) is a polynomial or a power series in x. There are many other ways of defining isotropic matrices. An interesting property of isotropic random matrices is that independently of the way they are defined, in the limit N → ∞ there is a one-to-one correspondence between the eigenvalue density of the matrix A and the eigenvalue density of the Hermitian matrix A † A associated with A [33]. By construction the matrix A † A is an invariant positive semi-definite random matrix. We denote the square root of this matrix by a = √ A † A. Note that the eigenvalues of a correspond to singular values of A. We can use the whole arsenal of methods discussed in the previous section to calculate the eigenvalue distribution of such matrices or their products.
The correspondence between the eigenvalue distribution of a large isotropic matrix and the distribution of its singular values was found by Feinberg and Zee [32] who applied field theoretical perturbative methods to solve the matrix model (39) and independently by Haagerup and Larsen [33] who used methods of free probability. Below we quote the result in the form given by Haagerup and Larsen [33]. In the terminology of free probability isotropic matrices in the limit N → ∞ correspond to variables that are called R-diagonal [34]. Here we prefer to call them large isotropic matrices. The invariance of the measure implies that the eigenvalue density of an isotropic matrix is invariant under rotation around the origin of the complex plane. For an isotropic matrix A it is convenient to define radial cumulative density function as the probability that a random eigenvalue of the isotropic matrix lies on the complex plane within the distance x from the origin: F A (x) = Prob(|λ A | ≤ x). The eigenvalue density is related to the radial cumulative function as for x > 0. This equation does not fix the value of the radial cumulative function for x = 0. It is fixed by an additional equation where F a (x) = x 0 ρ a (λ)dλ is the cumulative eigenvalue density function for a. This equation merely means that the number of zero eigenvalues of A is equal to the number of zero eigenvalues of a. The value F a (0) is generically zero except the situation when there is a point measure (Dirac delta) at zero in the eigenvalue density ρ a (λ), as for instance in the distribution with λ ± = (1 ± √ r) 2 , 0 < r < 1 and the support being {0} ∪ [λ − , λ + ]. This distribution naturally arises in the context of random matrices and it is called anti-Wishart distribution in physics literature and free Poissonian in math and engineering literature. The first term on the right hand side means that the fraction of zero eigenvalues of the corresponding random matrix a is (1 − r). The corresponding isotropic random matrix A = au 'inherits' all zero eigenvalues from a and thus the corresponding fraction of zero eigenvalues of A is also 1 − r.
Equations (41) and (42) provide an implicit relation between ρ A (z) and ρ a (λ) that allows one to calculate one from another.
A few comments are in order. Equation (41) was first published in physical literature in a slightly different (less transparent) albeit equivalent form [32]. In this paper it was also observed that for N → ∞ the support of the eigenvalue density of an isotropic random matrix has the shape of a single ring on the complex plane and that the radii of the ring can be derived from the eigenvalue density of a. The external radius is equal to R e = ( ∞ 0 ρ a (λ)λ 2 dλ) 1/2 and the internal one to R i = ( ∞ 0 ρ a (λ)λ −2 dλ) −1/2 . The ring reduces to a disk when R i = 0. In the paper [33] the discussion was formalized and all statements were proven rigorously using concepts of free probability. In addition to that the S transform form of equation (41) and the discussion of the point measure for zero eigenvalues p 0 = 0 (42) were given. From the discussion [33] it follows that in general there are three different types of solutions for F A . The ring solution that corresponds to F A (x) growing monotonically from 0 to 1 for x ∈ [R i , R e ] and for R i > 0, the disk solution that corresponds to F A (x) growing monotonically from 0 to 1 for x ∈ [0, R e ] and the disk solution with a finite fraction of zero eigenvalues p 0 > 0 that corresponds to F A (x) growing monotonically from p 0 to 1 for x ∈ [0, R e ].
Let us illustrate how the Haagerup-Larsen theorem (41) works in practice by giving an example. Consider a random matrix A from the Ginibre ensemble [36,37]. It is an N ×N random matrix whose elements are independent identically distributed, normally distributed, centered complex random variables. The probability measure for A is where dA is flat measure for N × N complex matrices. Obviously A is an isotropic random matrix. The limiting eigenvalue density for this matrix for N → ∞ is given by the uniform distribution on the unit disk [36,37] ρ A (z) = 1 π , |z| ≤ 1 .
It is zero outside the disk. The cumulative density Inserting it into equation (41) we obtain an equation for the S transform Next we find (22) χ a 2 (z) = z (z+1) 2 , the moment generating function (the inverse function for χ a 2 ) and the Green's function (12) Finally we calculate the eigenvalue density of a 2 (7) and correspondingly of a This is the density of singular values of the Ginibre matrix. We see that it is given by the quarter-circle law.
Let us now discuss the product of n independent isotropic matrices A 1 , A 2 , . . ., A n where A i 's are independent isotropic matrices [30,35]. The goal is to calculate the eigenvalue density and the singular value density in the limit N → ∞ having the knowledge of eigenvalue densities for the matrices in the product. Since the product itself is an isotropic matrix we can use equation (41) to write The S transform for p 2 = P P † can be calculated using the multiplication law (21) Since the right-hand side of this equation is a product of complex numbers, it does not depend on the order of multiplication. If instead of P we considered a product of matrices permuted with a permutation π P π = A π(1) A π(2) . . . A π(n) , we would obtain the S transform S p 2 π (z) identical to S p 2 (z) because multiplication on the right hand side of (53) is commutative. Therefore F P (x) = F Pπ (x) are equal for any permutation π. In other words, the product of isotropic matrices is spectrally commutative in the limit N → ∞.
Isotropic random matrices possess a striking property that has no counterpart in classical probability. Consider the product P = A 1 . . . A n of independent identically distributed isotropic random matrices. This can be realized by picking up n independent matrices from the same ensemble and multiplying them. It turns out that in the limit N → ∞ such a product has exactly the same eigenvalue density as the n-th power Q = A n of one matrix picked up at random from this ensemble [31]. Actually it was already observed in [33] that singular values of Q and P are identically distributed. This observation can be trivially extended to eigenvalues if one additionally assumes that p 0 = 0 (42). This is a sort of spectral self-averaging of multiplication in the limit N → ∞ because a single matrix from the ensemble is sufficiently representative in the n-fold product to give the same result as n independent random realizations. The proof is short so we recall it here in the form given in [31] where it was tacitly assumed that p 0 = 0. All constituents needed for this proof can be found already in [33].
We assume that p 0 = F A (0) = F a (0) = 0 (42). The S transform for the product of independent identically distributed matrices is the n-th power of the S transform for a single one S p 2 (z) = S n a 2 (x) as follows from equation (53) in which all S a 2 i (z) are replaced by S a 2 (z). Inserting S n a 2 (x) in place of S p 2 (z) in equation (52) and taking the n-th root of both sides we get Comparing the last equation to equation (41) we see that F P (x) = F A (x 1/n ). We can now argue that the radial cumulative density F Q (x) for the power Q is given by the same expression. Indeed, using directly the definition of radial cumulative density we have for x > 0, and for F Q (x) = F P (x) = 0 for x = 0 and hence F P (x) = F Q (x) = F A x 1/n . This completes the proof. If F A (0) = p 0 > 0 then the proof breaks down at x = 0. Recall that p 0 is the probability of drawing at random a zeroeigenvalue from the spectrum of random matrix A. The corresponding probability for the product of n independent A matrices is F P (0) = 1 − (1 − p 0 ) n , while for the n-th power F Q (0) = 1 − p 0 , so we see that F P (0) = F Q (0) for p 0 > 0.
In particular, one can use the product-power equivalence to calculate the eigenvalue density of the product P = A 1 . . . A n of n independent Ginibre matrices in the limit N → ∞. We know that F A (x) = x 2 for x ∈ [0, 1] for one Ginibre matrix (46) thus for the product of n The eigenvalue density for the product is (40) ρ P (z) = 1 nπ |z| and zero outside the unit circle. This calculation reproduces in a very simple way the result that was originally obtained by different methods [30]. This result can be generalized also to products of rectangular Gaussian matrices [38] or products [40,41] of Wigner matrices [39] with non-Gaussian entries. One can also find the distribution of singular values for the product of Ginibre matrices [7,35]. Inserting radial cumulative density function (57) to the Haagerup-Larsen formula (41) one can find the Green's function for p 2 = P † P [35] z n G n+1 and calculate the density function ρ p 2 (λ) by using (7). For n = 1 it gives the Wishart distribution (27); for n = 2 the distribution given in equation (35). One can find an explicit expression for any n [42]. The last step is to change variables λ → λ 2 to go from ρ p 2 (λ) to ρ p (λ).
To summarize, large isotropic random matrices are exceptional in many respects. There exists a one-to-one correspondence between the densities of eigenvalues and singular values [33]. Multiplication of independent isotropic matrices is spectrally commutative and self-averaging in the large N limit [31,35]. One can also show that if the eigenvalue density for an isotropic matrix has a power-law singularity at zero ρ A (z) ∼ |z| −s with 0 < s < 2, then the density of singular values behaves at zero as ρ p (λ) ∼ λ − s 4−s [35].

VII. GENERAL MULTIPLICATION LAW
So far we have discussed invariant and isotropic random matrices. We are now going to discuss the more general case of products of non-Hermitian matrices. The multiplication law for non-Hermitian matrices was derived in reference [21] by field theoretical methods.
To build some intuition let us discuss Gaussian non-Hermitian random matrices. One can construct such matrices from independent identically distributed Hermitian Gaussian matrices from the GUE ensemble. The probability measure for GUE matrices is dµ(a) ∼ e −N Tra 2 /2 da, where da is flat measure for Hermitian matrices. Let us take two independent such matrices a and b. The following combination is a complex random matrix; a gives rise to Hermitian sector of X and ib to the anti-Hermitian one. The resulting matrix X is a Ginibre random matrix [36,37] with the uniform eigenvalue distribution on the unit disk (45). One can also define elliptic random matrices [37,43] X = cos(α)a + i sin(α)b being linear combinations of Hermitian and anti-Hermitian random matrices with different coefficients controlled by a mixing parameter α. For α = 0 the random matrix X reduces to the Hermitian matrix a, for α = π/4 to the Ginibre matrix and for α = π/2 to the anti-Hermitian matrix ib. The eigenvalue density of the matrix X is constant on the support being an ellipse whose oblateness depends on α. The probability measure for matrix X can be calculated from the measure dµ(a, b) = dµ(a)dµ(b) ∝ e −N Tr(a 2 +b 2 )/2 dadb that factorizes since a and b are independent. Changing variables to X = cos(α)a + i sin(α)b and X † = cos(α)a − i sin(α)b one finds the probability measure for X [44,45] dµ(X) where τ = cos(2α). For this measure one can calculate two-point correlation functions X ab X cd = dµ(X)X ab X cd .
There are four such functions : that can be compactly written in the matrix form The matrix on the right hand side of the last equation defines the correlation structure of the model [21]. It plays an important role, as we shall see below. The parameter τ takes values between −1 and 1. For Ginibre matrix it is 0, and for GUE it is 1. Without going into details we quote the law of non-Hermitian matrix multiplication and refer the interested reader to [21]. We present it in the R transform form (26). The main difference as compared to the Hermitian case (26) is that the R transform is not a complex function but a quaternionic one. We denote quaternions by calligraphic letters. The R transform maps quaternions G onto quaternions G → R(G). The multiplication law reads Here G A = G A (z), G B = G B (z) and G AB = G AB (z) correspond to generalized quaternionic Green's functions for A, B and AB. They have the structure: where b(z) and g(z) are complex functions. We suppressed the indices of the corresponding matrix (A, B or AB) in the last equation. The density is related to the upper left element g(z) of the generalized quaternionic Green's function The superscripts L and R refer to the left and right adjoint actions of the following unitary matrix on matrices in the square brackets The diagonal elements of U are constructed from z being the argument of the Green's functions (66). The generalized R transform is related to the Green's function by an equation analogous to (17) This law looks complicated but in practice it reduces to a set of algebraic equations for g(z) from which one can calculate the eigenvalue density by taking the derivative with respect toz (67). Before we give an example let us mention an electrostatic analogy [45]. In this analogy eigenvalues of a random matrix represent positions of charges in two-dimensions in the given potential, the function g(z) is interpreted as the electric field g = (g x , g y ) with components g x = Reg(z) and g y = −Img(z) and equation (67) as the Gauss law in two dimensions [45]. It takes a more familiar form when written in the components of z = x + iy: We now illustrate how to use the multiplication law (65) to calculate the eigenvalue density of the product. We consider as an example the multiplication of Gaussian matrices constructed from the centered elliptic matrix X (62) by shifting and rescaling The scale parameter σ is a real number and the shift parameter q is a complex number. The quaternionic R transform for A is There are only two terms on the right-hand side because the matrix is Gaussian and thus all higher order cumulants vanish (16). The second term on the right-hand side is directly related to the matrix on the right hand side of equation (64) that defines the two-point correlation structure of the model [21,30]. The simplest example is q = 0, σ = 1 and τ = 0 (73). It corresponds to the Ginibre matrix. Let us calculate the eigenvalue density for the product of two such matrices. The first equation in the multiplication law (65) gives We used a shorthand notation c = (z/z) 1/4 . Because A and B are identically distributed, we search a symmetric solution b A (z) = b B (z) and g A (z) = g B (z).  X2) where X1 and X2 are independent standardized Ginibre matrices, each having the eigenvalue density uniformly distributed on the unit circle. The result is derived analytically [21]. The horizontal axes correspond to the real and imaginary axes of the complex plane and the vertical one to the value of the eigenvalue density ρ(z). The eigenvalue density has a support surrounded by a contour given by the Pascal's limaçon. (b) Solid line represents the Pascal's limaçon. Dots correspond to eigenvalues collected for ten matrices of size 1000 × 1000 generated by Monte-Carlo simulations. A few eigenvalues lie outside the limiting contour but this effect is expected by finite size analysis. The number of outliers decreases when the matrix size increases.
FIG. 2. Solid line represents the contour of the eigenvalue density for the product (1 + X)(1 + H) where X is a standardized Ginibre matrix and H is an independent standardized Hermitian GUE matrix. The curves forming the contour are calculated analytically [46] using the multiplication law (65). The extrapolation of these curves outside the support is shown in red dots. Black dots correspond to eigenvalues collected for ten matrices of size 1000 × 1000 generated by Monte-Carlo simulations.
law (65) is conveniently written in the R transform form (26). The R transform is a quaternionic function. One can use this law to find eigenvalue densities for products of large non-Hermitian random matrices. Isotropic matrices are a special class of non-Hermitian random matrices. These matrices have two exceptional features as far as matrix multiplication is concerned. Isotropic random matrices are spectrally commutative and self-averaging in the limit N → ∞. The first notion means that the eigenvalue distribution of the product of such matrices does not depend on the order of matrix multiplication. For example, the matrix abcd has the same eigenvalue density as bdac. The second one means that the product of n independent isotropic identically distributed matrices has the same eigenvalue distribution as the n-th power of a single matrix.
We finish by mentioning that recently a progress has been also made in the understanding of microscopic properties of eigenvalues statistics for products of random matrices. For example, an explicit form of the joint probability distribution for eigenvalues and singular values of the product of Ginibre matrices has been explicitly found [47,48]. Using these results one can determine not only the eigenvalue density but also joint eigenvalue distribution and eigenvalue correlations for products of large random matrices. One can also infer finite size corrections to the limiting laws discussed in this paper.