Calculating a ``density of states'' Let $n>0$ and $T > n$ be integers and let ${\bf c}$ be a symmetric, real and positive definite matrix. We define the ''density of states'' of a sample covariance matrix. In other words we have:
\begin{equation}
\omega^{(n,T)}({\bf c}) := \int\limits_{{\mathcal R}^{n T}} \delta\left( {\bf c} - x \cdot x^T \right) d^{n T} x
\end{equation}
where $x$ is a $n \times T$ real matrix.  The ''matrix delta function'' is defined as a product of delta functions over all distinct elements of the matrix. 
The question is to show that the quantity above equals:
\begin{equation}
\omega^{(n,T)}({\bf c}) := \left(\det({\bf c})\right)^{\frac{T-n-1}{2}} \cdot \frac{1}{2^n} \cdot \prod\limits_{j=1}^n S_{T-j}
\end{equation}
where $S_{T-1}$ is a surface of a unit sphere in $T$ dimensions. In my old question Multidimensional integral involving delta functions  I have calculated (with help of another user) this quantity for $n=2$.  

Firstly, note that the integrand is a product of $n$ delta functions each one related to one diagonal element of the matrix ${\bf c}$ and a product of another $\binom{n}{2}$ delta functions each one related to one of the cross-diagonal elements above the diagonal. Using the scaling property $\delta( a x) = 1/|a| \delta(x)$  we can always rescale the integration variable  so that the ''diagonal delta functions''  are centered at values equal to unity and the $(i,j)$ cross-diagonal delta function is centered  at $c_{i,j}/\sqrt{c_{i,i} c_{j,j}}$. Therefore we have:
\begin{equation}
\omega^{(n,T)}({\bf c}) = \frac{(\prod\limits_{i=1}^n c_{i,i})^{-1}}{\prod\limits_{1\le i < j \le n} \sqrt{ c_{i,i} c_{j,j}} } \cdot
\int\limits_{{\mathcal R}^{n T}} \prod\limits_{j=1}^n \delta\left(1 - ||\vec{x}_j||^2\right) \cdot \prod\limits_{1\le i < j \le n} \delta\left(\frac{c_{i,j}}{\sqrt{c_{i,i} c_{j,j}}}  - \vec{x}_i \vec{x}_j\right) d^{n T} x \cdot \left(\prod\limits_{i=1}^n \sqrt{c_{i,i}}\right)^T
\end{equation}
Here the prefactor on the right hand side comes from ''taking out'' the elements of the matrix ${\bf c}$ from all the delta's and the last factor on the right hand side comes form rescaling the integration variables, ie dividing them by square roots of the diagonal elements. After simplifying the expression we get:
\begin{equation}
\omega^{(n,T)}({\bf c}) = 
\left(\prod\limits_{j=1}^n c_{j,j}\right)^{\frac{T-n-1}{2}} \cdot
\int\limits_{{\mathcal R}^{n T}} \prod\limits_{j=1}^n \delta\left(1 - ||\vec{x}_j||^2\right) \cdot \prod\limits_{1\le i < j \le n} \delta\left(\frac{c_{i,j}}{\sqrt{c_{i,i} c_{j,j}}}  - \vec{x}_i \vec{x}_j\right) d^{n T} x
\end{equation}
Before proceeding further we note that the integral in question is ''rotationally invariant''. Indeed, if we write ${\bf c} = O\cdot {\mathcal D} \cdot O^T$  where $O \cdot O^T = O^T \cdot O = 1$ and ${\mathcal D}$ is a diagonal matrix then we have:
\begin{equation}
\omega^{(n,T)}({\bf c}) = 
\int\limits_{{\mathcal R}^{n T}} \delta\left( O\cdot {\mathcal D} \cdot O^T- O y y^T O^T\right) \cdot d^{n T} y \cdot \left|\det{O}\right|^T
\end{equation}
Here we simply changed variables by setting $x = O \cdot y$ where $y$ is a $n \times T$  matrix. The last term on the right hand side is the Jacobian of that transformation. Since from orthogonality $\left|\det(O)\right| = 1$ we have:
\begin{eqnarray}
&&\omega^{(n,T)}({\bf c}) = 
\int\limits_{{\mathcal R}^{n T}} \delta\left( O\cdot \left( {\mathcal D} -  y y^T \right) \cdot O^T\right) \cdot d^{n T} y \cdot \left|\det{O}\right|^T \\
&&= 
\frac{1}{\left|\det(O \cdot O^T)\right|}
\int\limits_{{\mathcal R}^{n T}} \delta\left(  \left( {\mathcal D} -  y y^T \right) \right) \cdot d^{n T} y  = \int\limits_{{\mathcal R}^{n T}} \delta\left(  \left( {\mathcal D} -  y y^T \right) \right) \cdot d^{n T} y
\end{eqnarray}
Therefore, we can without loss of generality assume that the cross-diagonal elements of the matrix ${\bf c}$ are zero. So we have:
\begin{equation}
\omega^{(n,T)}({\bf c}) = 
\left(\det(\bf c)\right)^{\frac{T-n-1}{2}} \cdot
\underbrace{\int\limits_{{\mathcal R}^{n T}} 
\prod\limits_{j=1}^n \delta\left(1 - ||\vec{x}_j||^2\right) \cdot
\prod\limits_{1\le i < j \le n} \delta\left(\vec{x}_i \vec{x}_j\right) d^{n T}x}_{I_{n,T}}
\end{equation}
Now, let us do the integral over $\vec{x}_1$. This clearly produces $1/2 S_{T-1}$. Now, in integrating over $\vec{x}_2,\cdots,\vec{x}_n$ we choose such a coordinate system that the azimuthal axis  coincides with $\vec{x}_1$. Now, the first $n-1$ ''cross-diagonal delta's'' enforce that all the vectors $\vec{x}_2,\cdots,\vec{x}_n$ are perpendicular to $\vec{x}_1$. In other words in the coordinate system that we have chosen  we have $|| \vec{x}_j||^2 = x_{j,n}^2 + ||\vec{y}_j||^2$  where each $x_{j,n}=0$ and the $\vec{y}_j$ vector is a projection onto a space perpendicular to $\vec{x}_1$.  Taking all that into account the integral reads:
\begin{equation}
I_{n,T} = \frac{1}{2} S_{T-1} \int\limits_{{\mathcal R}^{(n-1)(T-1)}} \prod\limits_{j=2}^n \delta\left(1 - ||\vec{y}_j||^2\right) \cdot
\prod\limits_{2\le i < j \le n} \delta\left(\vec{y}_i \vec{y}_j\right) d^{(n-1) (T-1)}y
= \frac{1}{2} S_{T-1} I_{n-1,T-1}
\end{equation}
To reiterate, the first diagonal delta produced the prefactor on the right hand side and the first $n-1$ cross-diagonal deltas , in our coordinate system, changed the scalar products of the $x$-vectors into scalar products of their projections onto the perpendicular space. Now by iterating the last equality on the right hand side we readily get the answer.
