Assume, $\mu=0$, to ease notation.
It is easily seen that $R=EXX'$ has two important properties: first, it is symmetric
(obviously); second it is positive semi-definite, i.e.\ the quadratic form
$u'Ru=\sum_k\sum_\ell R_{k\ell}u_k u_\ell$ is non-negative for all values of
the variables, positive or negative. The reason for the latter is that
$u'Ru $ is the expectation of a non-negative quantity:
$u'Ru=u'EXX'u=E(X'u)'(X'u) = E(X'u)^2 \ge 0$. Standard linear algebra
shows that $R$ has exactly $d$ (counting multiplicities) non-negative eigenvalues.
(In fact, they are strictly positive,
due to $\det(R) \ne 0$). Furthermore, the eigenvectors can be chosen to be orthonormal.
Letting $U$ be the matrix whose columns are these $d$ orthonormal eigenvectors
and $\Lambda$ the diagonal matrix with the eigenvalues in its diagonal,
we have $RU = U \Lambda$, from the very definition of the eigenvectors.
Now $U'=U^{-1}$, hence $R=U \Lambda U' = U \Lambda^{1/2} \Lambda^{1/2} U' =
P P'$, with $P=U \Lambda^{1/2}$. The matrix $P$ is non-singular and is called
the square root of $R$. We now define new random variables $Z$ by 
\[
X = P Z\,.
\]
The thing to observe is that the covariance matrix of $Z$ is
$\cov(Z) = EZZ' = EP^{-1}XX'P'^{-1}=P^{-1}R'P'^{-1}=P^{-1}PP'P'^{-1}=I$,
the identity matrix. Thus $Ee^{i\theta'Z} = e^{-\theta'\theta/2}=
\prod_{j=1}^d e^{-\theta_j^2/2} = \prod_{j=1}^d E e^{i\theta_jZ_j}$,
implying that the components of $Z$ are independent standard normal.
Hence the density $g$ of $Z$ is product:
$g(z) = \prod_j \frac{1}{\sqrt{2\pi}} e^{-z_j^2/2}$. 
Now $X=PZ$, hence its density $f$ is computed easily
by $f(x) = g(P^{-1}x)/|\det(P)|$, which yields the desired formula.
\QED

If $R$ has determinant zero then it is possible to ``reduce the dimension''
of the random vector $X$ so that the density exists.
In fact,
\begin{lemma}
The support of $X$ is the range of its covariance matrix $R$.
\end{lemma}
{\em Proof.}
Suppose that $R$ has rank $r$. Then it has $d-r$ eigenvalues at zero,
so that the matrix $\Lambda$ consists of a $d-r$ size block of zeros and 
the remaining non-zero eigenvalues. Hence now $R = PP'$,
where $P$ is a $d\times r$ matrix with rank $r$. 
We try again to find $Z$ so that $X=PZ$, where $Z$ is an $r$-dimensional random vector
with density.
If we manage to do this we will finish, since the range of $P$ is
the range of $R$. Observing that $P'P$ is an $r\times r$ non-singular matrix,
we pre-multiply by it to get $P'PZ = P'X$, hence $Z=(P'P)^{-1}P'X$.
So if we define $Z$ this way, we see that $EZZ'=
(P'P)^{-1}P'RP(P'P)^{-1} =
(P'P)^{-1}P'PP'P(P'P)^{-1}=I$, i.e.\ $Z$ is a collection of $r$ independent
standard normal variables.
The formula for $Z$ is the formula that solves a full-rank over-estimated
linear system. We usually write $Z = P^+ X$ and call $P^+$ the pseudo-inverse of $P$.
It remains to show that every linear function of $X$ is a linear function of $Z$.
Let $F=\lambda'X$, $G=\lambda'PZ$ with $Z=P^+X$.
Consider $F-G=\lambda'(X-PZ)$. Then $E(F-G)^2 = E\lambda'(X-PZ)(X'-Z'P')\lambda
=E\lambda'(XX'-XZ'P'-PZX'+PZZ'P')\lambda$.
But $EZX'=(P'P)^{-1}P'EXX'=(P'P)^{-1}P'PP'=P'$.
So $E(F-G)^2 = \lambda'(PP'-PP'-PP'+PP')\lambda=0$.
\QED
