  Random hyperspherical harmonics are Gaussian Laplace eigenfunctions on the
unit $d$-dimensional sphere ($d\ge 2$). We study the convergence in Total
Variation distance for their nonlinear statistics in the high energy limit,
i.e., for diverging sequences of Laplace eigenvalues. Our approach takes
advantage of a recent result by Bally, Caramellino and Poly (2020): combining
the Central Limit Theorem in Wasserstein distance obtained by Marinucci and
Rossi (2015) for Hermite-rank $2$ functionals with new results on the
asymptotic behavior of their Malliavin-Sobolev norms, we are able to establish
second order Gaussian fluctuations in this stronger probability metric as soon
as the functional is regular enough. Our argument requires some novel estimates
on moments of products of Gegenbauer polynomials that may be of independent
interest, which we prove via the link between graph theory and diagram
formulas.
