How can a multivariate normal be derived from a 1D normal using orthogonal matrices?

I've seen a very hand-wavy explanation of how a multivariate normal can be derived from the 1D case, but I wan't really able to find it done rigorously anywhere. Here's an outline of how it went:

• consider two independent random variables $X$ and $Y$ with a gaussian distribution
• because they're independent $P(X, Y) = P(X)P(Y)$
• we can directly write their covariance matrix (since it's diagonal), and we can also write as an outer product
• a correlation between $X$ and $Y$ can be viewed as rotating the joint distribution along the origin
• orthogonal matrix multiplication corresponds to rotation, which means we can (not sure how) find an orthogonal matrix and multiply the covariance matrix with it, which somehow gives the formula for a multivariate normal with an arbitrary covariance matrix

I don't remember how the determinant of the covariance matrix is derived.

Is this just a hand-wavy derivation that is mostly based on intution, or can this be done in a rigorous way to arrive at the actual multivariate normal?

A covariance matrix is symmetric, and hence can be diagonalized by an orthogonal transform. That is, for a desired covariance matrix $C$, there is an orthogonal matrix $O$ such that $O^\top CO=D$ is diagonal. Find this transform and scale individual independent random variables with a standard normal distribution by the square roots of the diagonal elements (which are positive since $C$ is positive semidefinite). The covariance matrix of this vector of scaled variables is $D$. Then multiply the vector by $O$, which makes the covariance matrix of the new vector $ODO^\top=C$.