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 *)
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.
