Deriving the sampling distribution of MLE for Normal distribution Let $X_1,\ldots,X_n$ be an observed random sample from $N_p(\mu, \Sigma)$. 
I know that the MLE of $\Sigma$ is $\frac{1}{n} \sum_i^n(X_i -\bar X)(X_i -\bar X)^T$, which is biased. 
We define $S = \frac{1}{n-1}\sum_i^n(X_i-\bar X)(X_i-\bar X)^T$ which is unbiased. 
I am proving that $(n-1)S \sim W_p(\Sigma, n-1)$ where $W$ represents the Wishart distribution. 
I let $X$ be the $n\times p$ matrix with $X_i^T$ as its i-th row. If $H = (I-\textbf{11}^T/n)$, then $HX$ returns the column-centered data matrix, and hence $S=X^THX/(n-1)$. Furthermore, $H=H^T$ and $\operatorname{Tr}(H)=n-1$, and hence $H$ is a $n\times n$ orthogonal projection matrix with rank $n-1$. It therefore has the spectral decomposition $H=\sum_{j=1}^{n-1} u_j u_j^T$, where the $u_j$s are orthonormal. Using this spectral decomposition, I get $$(n-1)S=\sum_{j=1}^{n-1}(X^Tu_j)(X^Tu_j)^T$$. 
If I can prove that the vectors $(X^Tu_j)_j$ are jointly MVN, and are $\text{i.i.d. } N_p(0,\Sigma)$, then the Wishart distribution is shown. Showing the expectation is 0 and variance is $\Sigma$ however is what I am stuck on. I wonder whether it has something to do with the column-centering matrix which subtracts the mean?   
 A: Disclaimer: This is not my proof, but rather an understanding of one of Lukács theorems.
First transform $X_k$ using orthogonal matrix $V$ so that it contains entries equal to $\frac{1}{ \sqrt{n}}$ in the last row. This means that the sum of entries in first $n-1$ rows are all zeros. Let's call the transformed vectors $Y_k$, i.e.
\begin{equation}
 Y_i = \sum_{k=1}^n v_{ik} X_k~~~\forall ~i=1, \ldots, n, 
\end{equation}
where $v_{ij}$ is the $(i,j)^{th}$ entry of $V.$ Equivalently, we can write
\begin{equation}
 Y = XV^\top.
\end{equation}
We can write $S$ as follows
\begin{equation}
 S 
 =
 \sum_{k=1}^n \left(X_k - \bar{X}\right)\left(X_k - \bar{X}\right)^\top
 =
 \sum_{k=1}^n X_kX_k^\top - n \bar{X}\bar{X}^\top.
\end{equation}
Using some straightforward math, the above could be shown to be
\begin{equation}
 S 
 =\sum_{k=1}^n Y_kY_k^\top- n \left(\frac{Y_n}{\sqrt{n}}\right)
 \left(\frac{Y_n}{\sqrt{n}}\right)^\top
 =
 \sum_{k=1}^{n-1} Y_kY_k^\top
\end{equation}
i.e. a sum of $n-1$ outer-product. Now, we can see that $S$ is the dyadic sum of $n-1$ independent normal vectors of mean $0$ and variance $\Sigma$. Therefore, by using the Wishart-distribution definition, we have that
\begin{equation}
 S \sim \mathcal{W}_p(n-1, \Sigma). 
\end{equation}
