Boltzmann–Gibbs measures. If $E$ is a nice space such as $\mathbb{R}^d$, $\mathbb{Z}^d$, or $\mathbb{S}^{n}$, equipped with a reference measure denoted $\mathrm{d}x$, if $V:E\to(-\infty,+\infty]$ is a function, and if $\beta$ is a real, then the associated Boltzmann–Gibbs measure at inverse temperature $\beta$ with energy $V$ is $$\mathrm{d}\mu_\beta(x)=\frac{1}{Z_\beta}\exp(-\beta V(x))\mathrm{d}x\quad\text{provided that}\quad Z_\beta:=\int\exp(-\beta V(x))\mathrm{d}x<\infty.$$ It maximizes entropy at fixed average energy, and this explains its emergence. More precisely, let $\mu$ be a probability measure on $E$ with finite entropy and same average energy as $\mu_\beta$ :
$$\mathrm{S}(\mu):=-\int\frac{\mathrm{d}\mu}{\mathrm{d}x}\log\frac{\mathrm{d}\mu}{\mathrm{d}x}\mathrm{d}x<\infty\quad\text{and}\quad\int V\mathrm{d}\mu=\int V\mathrm{d}\mu_\beta<\infty.$$ Then, by the Jensen inequality for the strictly convex function $u\mapsto u\log(u)$, $$\mathrm{S}(\mu_\beta)-\mathrm{S}(\mu)=\int\frac{\mathrm{d}\mu}{\mathrm{d}\mu_\beta}\log\frac{\mathrm{d}\mu}{\mathrm{d}\mu_\beta}\mathrm{d}\mu_\beta\geq0\quad\text{with equality iff \mu=\mu_\beta.}$$ This is standard in statistical physics (Boltzmann) and statistical mechanics (Gibbs). It turns out that these objects are here and there in mathematics, and the interpretation in terms of temperature and energy is helpful for the intuition, regardless of the true connection to physics.

Random matrix ensembles. Wigner, Mehta, Dyson and others have defined Boltzmann–Gibbs measures on matrices, by taking for $E$ say the set of $n\times n$ Hermitian matrices, namely $$\mathrm{d}\mu(H)=\frac{1}{Z_\beta}\exp(-\beta\mathrm{Tr}(V(H)))\mathrm{d}H.$$ This law is unitary invariant namely invariant by conjugation with respect to all unitary matrices, due to the cyclicity of the trace and the spectrality of the energy. Under the spectral decomposition $H=UDU^*$, this gives a random unitary matrix $U$ which is uniform, independent of $D$, while the diagonal entries of $D=\mathrm{diag}(\lambda_1,\ldots,\lambda_n)$ have density
$$\propto\exp\Bigr(-\beta\sum_{i=1}^nV(\lambda_i)\Bigr)\prod_{i < j}|\lambda_i-\lambda_j|^2.$$ The term $\prod_{i < j}|\lambda_i-\lambda_j|^2$ is the Jacobian of the change of variable, and the power $2$ is algebraic : two reals are needed to encode a complex number. If we play with real symmetric matrices we would get $1$ instead of $2$, while quaternionic symplectic matrices would give $4$. We can view this law as a Coulomb gas, which is a Boltzmann–Gibbs measure
$$\propto\exp\Bigr(-\beta\sum_{i=1}^nV(\lambda_i)-2\sum_{i < j}\log\frac{1}{|\lambda_i-\lambda_j|}\Bigr)$$ at inverse temperature $\beta$ and energy $\sum_{i=1}^nV(\lambda_i)+\frac{2}{\beta}\sum_{i< j}\log\frac{1}{|\lambda_i-\lambda_j|}$.

Beta-ensembles. For a real $\beta\geq0$ and $V:\mathbb{R}\to\mathbb{R}$, the $\beta$-ensemble is the set of $n\times n$ Hermitian matrices of the form $H=UDU^*$ where $U$ and $D$ are independent with $U$ unitary uniform and $D=\mathrm{diag}(\lambda_1,\ldots,\lambda_n)$ with diagonal entries of density
$$\propto\exp\Bigr(-\beta\sum_{i=1}^nV(\lambda_i)\Bigr)\prod_{i< j}|\lambda_i-\lambda_j|^\beta$$ which is nothing else but the Boltzmann–Gibbs measure (more precisely the Coulomb gas) $$\propto\exp\Bigr(-\beta\sum_{i=1}^nV(\lambda_i)-\beta\sum_{i< j}\log\frac{1}{|\lambda_i-\lambda_j|}\Bigr)$$ at inverse temperature $\beta$ and energy $\sum_{i=1}^nV(\lambda_i)+\sum_{i< j}\log\frac{1}{|\lambda_i-\lambda_j|}$. The law of $H$ is not the Boltzmann–Gibbs measure $\propto\exp(-\mathrm{Tr}(V(H))$ introduced previously, except when $\beta=2$.

Beware that the term beta ensemble is used for the Boltzmann–Gibbs measure on $\mathbb{R}^n$ as well as for the unitary invariant random Hermitian matrix for which it describes the spectrum.

When $V$ is homogeneous, in particular when $V(x)=x^2$, then we can cancel the $\beta$ in front of $V$ by scaling, while the one in front of the interaction is a rigid shape parameter.

Fishy? Since the spectral decomposition $H=UDU^*$ is not linear, putting a Boltzmann–Gibbs measure on $H$ or on $D$ does not produce the same objet. Nothing more.

Gaussians. A famous theorem due to Maxwell states that a random vector of $\mathbb{R}^n$, $n\geq2$, has independent components and has a rotationnally invariant law if and only if it follows an isotropic Gaussian law. Following Mehta, there exists a matrix analogue which states that an $n\times n$ say Hermitian random matrix has independent centered entries and has a unitary invariant law if and only if it follows a Boltzmann–Gibbs measure of the form $$\propto\exp(-\beta\mathrm{Tr}(H^2).$$ Note that $V(x)=x^2$ is $2$-homogeneous, the temperature $1/\beta$ is a Gaussian variance.

On the other hand, in the rigid cases where $V$ is Hermite/Laguerre/Jacobi, and following Dumitriu and Edelman, there exists a random matrix model, which is tridiagonal, with independent entries, for which the spectrum is the $\beta$ ensemble, but it is not unitary invariant. This means that for $\beta$ ensembles with $\beta\neq2$, we have to choose between unitary invariance without independent entries or independent entries without unitary invariance.

Is $\beta=2$ special? The classical random matrix ensembles, see for instance Mehta book, which are Gaussians on symmetric, Hermitian, or symplectic matrices, correspond to physical symmetries in the Hamiltonian. This allows to see the special value $\beta=2$ as being associated to something special and rigid. On another direction, it is believed that the high dimensional local behavior of the spectrum of beta ensembles has a threshold phenomenon at $\beta=2$.

Orthogonal polynomials. The beta ensemble as a Coulomb gas on $\mathbb{R}^n$ for arbitrary beta in the rigid Hermite/Laguerre/Jacobi choice for $V$ also emerges from combinatorics, and the associated multivariate orthogonal polynomials are then the solution of a Dyson ODE : this was explored by Jack, Macdonald, Lassalle, and others. In dimension $n=1$, the Hermite/Laguerre/Jacobi trilogy was also identified by Bakry as the unique possible cases for which the orthogonal polynomials of the invariant measure are the eigenfunctions of the diffusion Markov process. This trilogy is also at the heart of the Dumitriu-Edelman tridiagonal random matrix model for beta ensembles.

Physical relevance of beta ensembles. The modelling of the energy levels of heavy nuclei lead to put a Boltzmann–Gibbs measure on a Hamiltonian, and then to the classical Gaussian ensembles of random matrix theory, for which the spectrum is a Coulomb gas with $\beta\in\{1,2,4\}$ (moreover this beta has nothing to do with the $\beta$ of the Boltzmann–Gibbs measure on matrices, as explained above). The three values of $\beta$ in the Coulomb gas describing the eigenvalues are algebraic, they correspond to the nature of the matrix entries, related to the symmetries of the Hamiltonian. It is the same story for multivariate statistics, in which the empirical covariance matrices lead to random matrix models associated with the same three algebraic values of $\beta$. In both cases, the restriction on $\beta$ comes from the fact that we put a randomness on the matrix entries and we consider then the law of the spectrum.

In contrast, the beta ensembles with arbitray beta seem artificial for classical nuclear physics and classical multivariate statistics. However they appear in other contexts notably as Coulomb gases, for instance the Coulomb gas on $\mathbb{C}^n$ at inverse temperature $\beta\in\{2,4,6,8,\ldots\}$ with $V(z)=|z|^2$ appears from the wave function of the Laughlin fractional quantum Hall effect. This gas also appears in the description of Ginzburg–Landau vortices, and there is a conjecture about the convergence to a rigid structure when $\beta\to\infty$ (crystallization). It turns out that when $\beta=2$ it also describes the spectrum of complex square random matrices with iid complex Gaussian entries (complex Ginibre ensemble), but this random matrix ensemble is not normal and does not coincide with the random normal matrix ensemble that we can also construct from the same spectral gas, also it is not a random matrix beta ensemble in the sense introduced above.

All this suggests that the relevance of the beta ensembles for arbitrary beta may come directly from a specific modelling producing directly a Coulomb gas, while the random matrix cases impose an algebraic value of $\beta$ : $1$, $2$, or $4$. However, the random matrix models associated to the beta ensembles with arbitrary beta such as the Dumitriu-Edelman random matrices may help to study the initial modelling that produces a Coulomb gas with arbitrary beta.