Given $H_1$ and $H_2$ i.i.d. $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$x$N$ matrices in the large $N$ limit.
The typical example is for a matrix $H$ from the $GUE$. Diagrammatic methods find that moments of the eigenvalues behave asymptotically at large $N$ as $\mathbb{E}[\lambda^{2n}] = \mathbb{E}[\frac{1}{N} 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 of the Catalan numbers. These numbers arose in a different context 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.