Given $H_1$ and $H_2$ i.i.d. $\mathit{GUE}$ matrices, what is the single eigenvalue distribution of $H_1 H_2 H_1$ in the large $N$ limit? This matrix is Hermitian, and so its eigenvalues are still real.

As some background, I'm practicing moment methods to find the distribution of single eigenvalues of random $N\times N$ matrices in the large $N$ limit.

The typical example is for a matrix $H$ from the $\mathit{GUE}$. Diagrammatic methods find that moments of the eigenvalues behave asymptotically at large $N$ as $\mathbb{E}[\lambda^{2n}] = \mathbb{E}[\frac{1}{N} \operatorname{Tr}[H^{2n}]] \sim c_n$, where the $c_n = \frac{\binom{2n}{n}}{n+1}$ are the Catalan numbers.

The next step is that one recognizes these as the moments of the Wigner semicircle distribution $f_{\lambda}(x) = \frac{1}{2\pi} \sqrt{4 - x^2}$ supported on $[-2,2]$. Alternatively, one can use these moments to calculate the Stieltjes transform, $R(z) = \mathbb{E}[\frac{1}{z - \lambda}]$, by way of the generating function for the Catalan numbers. The inverse Stieltjes–Perron formula gives $f_{\lambda}(x) = \lim_{\epsilon \to 0^+} \frac{R(x+i\epsilon) - R(x-i\epsilon)}{-2\pi i}$, where the difference is across a branch cut. This gives an algorithmic way for identifying the probability distribution.

As a comparison to the problem at hand, notice that if one wanted the eigenvalue distribution of $H_1^3$, one would find from the Wigner semi-circle the distribution $f_{\lambda}(x) = \frac{1}{6\pi} \frac{\sqrt{4 - x^{2/3}}}{x^{2/3}}$ supported on $[-8,8]$. Numerically, I instead find that the distribution of eigenvalues of $H_1 H_2 H_1$ is supported on roughly $[-3, 3]$, and looks similar to but not quite a rescaling of the above density for $H_1^3$.

Running through diagrammatic arguments applied now to the case of $H_1 H_2 H_1$, I calculate by hand the moments $\mathbb{E}[\lambda^2] = 1$, $\mathbb{E}[\lambda^4] = 4$, $\mathbb{E}[\lambda^4] = 22$. These moments (and a couple higher moments estimated numerically) suggest $\mathbb{E}[\lambda^{2n}] = t_n$, where $t_n = \frac{ \binom{4n}{n}}{3n+1}$ are a generalization OEIS A002293 of the Catalan numbers. These numbers arose in a different context (see Another generalization of parity of Catalan numbers) in another problem of mine.

However, I do not recognize the probability distribution giving these moments. I find that the Stieltjes transform $R(z)$ is a root of the equation $z^2 R(z)^4 - z R(z) +1 =0$, but the resulting quartic roots appear both complicated, and, more critically, with a complicated branch cut structure.

As an aside, I would be satisfied by an answer identifying the probability distribution of $\lambda$ for which $\mathbb{E}[\lambda^{2n}] = t_n$, regardless of whether it uses random matrix theory techniques.

Update: Empirically, the divergence of the density of eigenvalues of $H_1 H_2 H_1$ at small argument $x$ goes asymptotically proportional to $|x|^{-1/2}$, rather than the $|x|^{-2/3}$ seen with $H^3$, so my comparison above might not give qualitative intuition about the solution.

With an answer from Mathematica Stack Exchange user293787 (code copied below), I also have the following Mathematica code that generates the density above in terms of hypergeometric functions. The main idea is to write the Fourier transform of the density as $\mathbb{E}[e^{i k X}] = \sum_{n=0}^{\infty} \frac{(-1)^n k^{2n}}{(2n)!} \frac{\binom{4n}{n}}{3n+1}$, which equals ${}_{2}F_{3}(\{1/4, 3/4\}, \{2/3, 1, 4/3\}, -((64 k^2)/27))$. Mathematica evaluates the inverse Fourier transform in terms of hypergeometric functions,

f=HypergeometricPFQ[{1/4,3/4},{2/3,1,4/3},-64 k^2/27];

(* this returns a symbolic result *)

(* check normalization, gives 1 *)

(* plot *)

The symbolic result is quite complicated and omitted in the code above. However, it appears from this solution that the exact bounds for the support are not $[-3, 3]$ but instead the close $[-\frac{16}{3 \sqrt{3}}, \frac{16}{3 \sqrt{3}}]$. My hope is that the result could be further simplified into a more (for me) intuitive form.


1 Answer 1


A follow-up answer by Bob Hanlon on Mathematica Stack Exchange further simplifies Mathematica's output of the probability density, found via the Fourier transform technique detailed in the update in my question:

$$ \begin{aligned} f_{\lambda}(x) =& \,\frac{1}{32 \pi^3} \left( 3^{3/4} \Gamma\left(\frac{3}{4}\right)^2 \Gamma\left(\frac{5}{12}\right) \Gamma\left(\frac{13}{12}\right) \left(32|x|^{-1/2} {}_{3}F_2\left(\frac{-1}{12}, \frac{1}{4}, \frac{7}{12}; \frac{1}{2}, \frac{3}{4}; \frac{27 x^2}{256}\right) \\ \vphantom{\int_1^2} - |x|^{1/2}{}_{3}F_2\left(\frac{5}{12}, \frac{3}{4}, \frac{13}{12}; \frac{5}{4}, \frac{3}{2}; \frac{27 x^2}{256}\right)\right) - 8 \pi^2 {}_{3}F_2\left(\frac{1}{6}, \frac{1}{2}, \frac{5}{6}; \frac{3}{4}, \frac{5}{4}; \frac{27 x^2}{256}\right)\right) \end{aligned} $$

for $x$ in $[-\frac{16}{3 \sqrt{3}}, \frac{16}{3 \sqrt{3}}]$, and $0$ otherwise.

If you see a nice way to interpret or further simplify this result, feel free to add that as a comment or an answer. At the least, the prefactor of the $x^{-1/2}$ divergence is clear, though I'm also curious about the behavior close to the edges.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.