Some brief numerics gives the following error distributions for the fit of 30 samples created from a standard normal distribution then fit to a univariate Gaussian.
The quartiles are indicated. It is assumed that this level of variation is desired in the multi-dimensional case.
I don't have the time to beat up MatLab to get the total result, so I will share my "rule of thumb". The 30 is provided as a rule of thumb, or heuristic so it is assumed that heuristics are not unacceptable.
My heuristic is to use Pascal's triangle multiplied by the univariate case.
If I am using 2d data then I go to the 2nd row and sum it to get 2x the number of samples, or 60 samples. For 3d data I go to the 3rd row and sum it to get 4x the number of samples or 120 samples. For 5d data I go to the 5th row and sum it to get 16x the number of samples, or 480 samples.
Best of luck.
EDIT:
It was intuitive, but everything has to be defended in math. I can't just take leaps from formulation of polynomial forms from Finite Elements with experience to get a ballpark.
The equation for the sum of the $ k^{th}$ row of Pascal's triangle is $ 2^k$.
My idea for the approach here is to equate the AIC of a higher-dimensional distribution with more samples to a reduced dimensional distribution with fewer samples.
The Akaike Information Criterion (AIC) is defined as $ AIC = n \log( \frac {RSS}{n}) + 2*k$ where $ RSS$ is residual sum of squares, $ n$ is sample count, and $ k$ is parameter count for the model.
$ AIC_1 = AIC_2$
$ n_1 \log(\frac {RSS_1}{n_1}) +2k_1 = n_2 \log(\frac {RSS_2}{n_2}) +2k_2$
For each dimension that we eliminate this means that the mean loses a row and the covariance loses both a row and a column. We can state this as
$ k \left( d\right)= d^2+d$.
of
$k \left( d+1 \right) - k \left( d\right) = 2 d + 2$
Assuming the error per sample point is constant relates the residual sum of squares to the sample count, and the term in the logarithm stays constant. The difference in sample count becomes a scaling constant.
so we have:
$ n_1 A +2(k_2+2d+2) = n_2 A +2k_2 $
solving for the increase in samples with dimension gives:
$ n_2- n_1 = (2(k_2+2d+2) - 2k_2) A^{-1} = (4 d+4 ) \cdot A^{-1} $
So what is the scaling function? Lets assume that for a 2-dimensional multivariate Gaussian the number of samples required is 15 per parameter. There are 2 means and 4 elements of the covariance therefore 6 parameters or 90 samples. The difference is 60 samples, the value of $ A^{-1} = 5$.
At this point I would say that the heuristic starts a little low but ends up being about 2x the number of samples required. Its range of best utility, in my personal opinion, is around 4 dimensions or so.
If I have to edit this answer again, I'm going to delete it.