Why do we use determinant for multivariate normal distribution? While learning statistics, I have a question why is the determinant used in the multivariate normal distribution. 
When I look for the answer on the internet, so far every answer I looked at was basically saying that it works, so we use that. 
But what I want is if there is a mathematical relation between multivariate normal distribution and determinant (volume factor of linear transformation or some other definition). 
There was one answer that by using determinant we can make the integral of the density over $R^{n}$ equal to $1$. This sounds nice, but if there is another intuition, please share it. 
 A: I don't think there's anything more behind the appearance of the determinant of the covariance matrix, $\ \Sigma\ $, say, beyond the fact that if $\ \det(\Sigma)\ne 0\ $ then $\ \int_{\mathbb{R}^n}e^{-\frac{1}{2}(x-\mu)^\top \Sigma^{-1}(x-\mu)}dx = (2\pi)^{\frac{n}{2}}\sqrt{\det(\Sigma)}\ $, so $\ N=(2\pi)^{-\frac{n}{2}}\det(\Sigma)^ {-\frac{1}{2}}\ $ is the normalising factor for which $\ N\int_{\mathbb{R}^n}e^{-\frac{1}{2}(x-\mu)^\top \Sigma^{-1}(x-\mu)}dx = 1\ $.
The integral $\ \int_{\mathbb{R}^n}e^{-\frac{1}{2}(x-\mu)^\top \Sigma^{-1}(x-\mu)}dx\ $ can be evaluated by change of variables from $\ x\ $ to $\ y= U(x-\mu)\ $, where $\ U\ $ is the orthogonal matrix that diagonalises $\ \Sigma\ $:
$$
U^\top\Sigma U=\pmatrix{\sigma_1&0&\dots&0\\
0&\sigma_2&\dots&0\\
\vdots&\vdots&\ddots&\vdots\\
0&0&\dots& \sigma_n}\ .
$$
Because $\ \Sigma\ $ is a positive definite matrix whenever $\ \det(\Sigma)\ne 0\ $, $\ U\ $ always exists , the eigenvalues $\ \sigma_1, \sigma_2, \dots,\sigma_n\ $ of $\ \Sigma\ $ are all positive, and  $\ \det(\Sigma)=\sigma_1\sigma_2\dots\sigma_n\ $.
Because orthogonal matrices preserve volume, the change of variables gives
\begin{align}
\int_{\mathbb{R}^n}e^{-\frac{1}{2}(x-\mu)^\top \Sigma^{-1}(x-\mu)}dx&=\int_{U^{-1}\mathbb{R}^n+\mu}e^{-\frac{1}{2} y^\top U\Sigma^{-1} U^{-1}y}dy\\
&= \int_{\mathbb{R}^n}e^{-\frac{1}{2} y^\top\left(U^\top\Sigma U\right)^{-1}y}dy\\
&=  \int_{\mathbb{R}^n}e^{-\frac{y_1^2}{2\sigma_1}- \frac{y_2^2}{2\sigma_2}-\dots-\frac{y_n^2}{2\sigma_n}}dy\\
&=\int_{-\infty}^\infty e^{-\frac{y_1^2}{2\sigma_1}}dy_1 \int_{-\infty}^\infty e^{-\frac{y_2^2}{2\sigma_2}}dy_2\dots\int_{-\infty}^\infty e^{-\frac{y_n^2}{2\sigma_n}}dy_n\\
&=\sqrt{2\pi\sigma_1}\sqrt{2\pi\sigma_2}\dots \sqrt{2\pi\sigma_n}\\
&= (2\pi)^\frac{n}{2}\sqrt{\det(\Sigma)}
\end{align}
