The McCarthy multimatrix ensemble of random matrices. For all integers $n\geq1$ and $d\geq1$, let $\mathbb{M}_{n,d}$ be the set of $d$-tuples $(M_1,\ldots,M_d)$ of $n\times n$ Hermitian matrices such that
$M_pM_q=M_qM_p\quad\text{for all 1\leq p,q\leq d.}$ We equip this hypersurface with the trace of the Gaussian distribution, namely
$(M_1,\ldots,M_d) \mapsto\frac{1}{Z_{n,d}} \mathrm{e}^{-\sum_{k=1}^d\mathrm{Tr}(M_k^2)}\mathrm{d}M$ where $\mathrm{d}M$ is the trace on $\mathbb{M}_{n,d}$ of the (product) Lebesgue measure on $d$-uples of $n\times n$ Hermitian matrices, and $Z_{n,d}$ is the normalizing constant. Since the Hermitian matrices $M_1,\ldots,M_d$ commute, they are diagonalizable in the same orthonormal basis, namely there exists a single $n\times n$ unitary matrix $U$ and real $n\times n$ diagonal matrices $D_1,\ldots,D_d$ carrying the eigenvalues of $M_1,\ldots,M_d$ respectively such that
$U(M_1,\ldots,M_d)U^* =(D_1,\ldots,D_d).$ The eigenvectors couple the eigenvalues : $\lambda_i\in\mathbb{R}^d$ is an eigenvalue of $(M_1,\ldots,M_d)$ when there exists $u\in\mathbb{R}^n$, $u\neq0$, such that $M_pu=\lambda_{i,p}u$ for all $1\leq p\leq n$. The computation of the Jacobbian of the spectral change of variable provides the following remarkable fact: the joint law of the eigenvalues $\lambda_1,\ldots,\lambda_n$ of $(M_1,\ldots,M_d)$ has probability density function
$(\lambda_1,\ldots,\lambda_n)\in(\mathbb{R}^d)^n \mapsto\frac{1}{Z_{n,d}} \mathrm{e}^{-\sum_{i=1}^n\|\lambda_i\|^2} \prod_{1\leq i< j\leq n}\|\lambda_i-\lambda_j\|^2 \mathrm{d}\lambda$ where $\mathrm{d}\lambda$ is the Lebesgue measure on $(\mathbb{R}^d)^n$, $\|x\|^2:=x_1^2+\cdots+x_d^2$ is the squared Euclidean norm of $\mathbb{R}^d$, and $Z_{n,d}$ is the normalizing constant. When $d=1$, we recover the usual formula for the Gaussian Unitary Ensemble (GUE), studied notably by Freeman Dyson. Thus the model generalizes the GUE into a “MacCarthy Unitary Ensemble” (MUE). When $d=2$, we obtain the formula for the spectrum of the complex Ginibre ensemble, providing a novel interpretation of this formula in terms of spectrum of two commuting Hermitian matrices!

The log-gas picture and the analogue of the Wigner theorem. The law above is a log-gas of $n$ particles on $\mathbb{R}^d$ with quadratic external field, namely
$\mathrm{e}^{-\sum_{i=1}^n\|\lambda_i\|^2} \prod_{1\leq i< j\leq n}\|\lambda_i-\lambda_j\|^2 =\exp\Bigr(n\int V\mathrm{d}\mu_n+\iint_{\neq}W\mathrm{d}\mu_n^{\otimes 2}\Bigr)$ where
$V(x):=\|x\|^2,\quad W(x,y):=\log\frac{1}{\|x-y\|},\quad \mu_n:=\frac{1}{n}\sum_{i=1}^n\delta_{\lambda_i}.$ What is the behavior of the empirical measure $\mu_n$ under this law? By computing the correlation functions, or by using the Laplace method, we get that after scaling by $n^{-1/2}$, the empirical measure in high dimension $n$ tends to the equilibrium measure
$\mu_{\mathrm{eq}} =\arg\min_\mu\Bigr(\int V\mathrm{d}\mu+\iint W\mathrm{d}\mu^{\otimes 2}\Bigr)$ where the minimum runs over the set of probability measures on $\mathbb{R}^d$. It turns out that these equibrium measures were already computed, quite recently for $d\geq3$, namely
$\mu_{\mathrm{eq}} = \begin{cases} \frac{1}{\pi} \sqrt{2-x^2} \mathbf{1}_{|x|\leq\sqrt{2}} & \text{if d=1}\\ \frac{1}{\pi} \mathbf{1}_{\|x\|\leq1} & \text{if d=2}\\ \frac{3}{2\pi^2}\frac{1}{\sqrt{\frac{2}{3}-\|x\|^2}}\mathbf{1}_{\|x\|\leq\sqrt{\frac{2}{3}}} &\text{if d=3}\\ \sigma_{S^{d-1}(\frac{1}{\sqrt{2}})} & \text{if d\geq 4} \end{cases}$ where $\sigma_{S^{d-1}(R)}$ is the uniform distribution on the sphere $S^{d-1}:=\{x\in\mathbb{R}^d:\|x\|=R\}$ of radius $R$. Also $\mu_{\mathrm{eq}}$ is radially symmetric and its one-dimensional projections are semi-circle distributions. We have an analogue or generalization of the Wigner theorem for the McCarthy multimatrix Ensemble. The Wigner theorem for GUE corresponds to $d=1$.

Note. Lydia Giacomin is studying these questions as a side project during her PhD.