# Eigenvalues of $H_1 H_2 H_1$, where $H_1$, $H_2$ independent $\mathit{GUE}$

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.

\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.