  We evaluate, in the large-$N$ limit, the complete probability distribution
$\mathcal{P}(A,m)$ of the values $A$ of the sum $\sum_{i=1}^{N} |\lambda_i|^m$,
where $\lambda_i$ ($i=1,2,\dots, N$) are the eigenvalues of a Gaussian random
matrix, and $m$ is a positive real number. Combining the Coulomb gas method
with numerical simulations using a matrix variant of the Wang-Landau algorithm,
we found that, in the limit of $N\to \infty$, the rate function of
$\mathcal{P}(A,m)$ exhibits phase transitions of different characters. The
phase diagram of the system on the $(A,m)$ plane is surprisingly rich, as it
includes three regions: (i) a region with a single-interval support of the
optimal spectrum of eigenvalues, (ii) a region emerging for $m<2$ where the
optimal spectrum splits into two separate intervals, and (iii) a region
emerging for $m>2$ where the maximum or minimum eigenvalue ``evaporates" from
the rest of eigenvalues and dominates the statistics of $A$. The phase
transition between regions (i) and (iii) is of second order. Analytical
arguments and numerical simulations strongly suggest that the phase transition
between regions (i) and (ii) is of (in general) fractional order $p=1+1/|m-1|$,
where $0<m<2$. The transition becomes of infinite order in the special case of
$m=1$, where we provide a more complete analytical and numerical description.
Remarkably, the transition between regions (i) and (ii) for $m\leq 1$ and the
transition between regions (i) and (iii) for $m>2$ occur at the ground state of
the Coulomb gas which corresponds to the Wigner's semicircular distribution.
