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](https://oeis.org/A002293) of the Catalan numbers. These numbers arose in a different context (see [Another generalization of parity of Catalan numbers](https://mathoverflow.net/questions/409735/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](https://mathematica.stackexchange.com/a/274712/69425) (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 *)
g=1/Sqrt[2*Pi]*InverseFourierTransform[f,k,x];

(* check normalization, gives 1 *)
NIntegrate[g,{x,-Infinity,Infinity}]

(* plot *)
Plot[g,{x,-4,4},PlotStyle->Thick]

```

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.



Eigenvalues of $H_1 H_2 H_1$, where $H_1$, $H_2$ independent $\mathit{GUE}$
A follow-up answer by [Bob Hanlon](https://mathematica.stackexchange.com/a/274735/69425) 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.


