When trying to approximate a distribution of random vectors $D$ by using multivariate Gaussian, what properties must we ensure that $D$ has? I.e., what distributions can be estimated by multivariate Gaussian fitting?
In a multivariate Gaussian distribution different coordinates can be related only in a linear way. Run regression of $X(i)$ on $X(j)$ and $X(j)^2$. If the coefficient of $X(j)^2$ is statistically significant, then the multivariate distribution is not Gaussian. A scatterplot matrix allows one to detect any aberrations informally.
Also, it goes without saying that the marginal distribution of each coordinate $X(i)$ must be "bell-shaped". In other words, it should pass Gaussian goodness-of-fit tests. Such tests are Kolmogorov-Smirnov, Shapiro-Wilk, Jarque-Bera, ...
Multivariate Gaussian distribution is not allowed to be degenerate. In particular, the covariance matrix of the coordinates must have positive determinant. For example, $Y$ can be a Gaussian variable but $(Y, -Y)$ will not be a Gaussian vector since the coordinates are linearly related.