I did a preliminary feasibility analysis of our methods and it appears possible that one may be able to tighten our $n^\epsilon$ loss to something more like $\exp( \sqrt{n} )$ in the Gaussian case, but this is still well short of what you want. The main obstacle is potential coupling between permanents of minors, which we were not able to fully avoid.
Here's the heuristic calculation. Suppose that a Gaussian $k \times k$ permanent has some distribution $P_k$. Then a Gaussian $k+1 \times k+1$ permanent $P_{k+1}$, by cofactor expansion, looks like
$$ P_{k+1} = \sum_{i=1}^{k+1} (-1)^i x_i P_k^{(i)}$$
where $x_i$ are iid Gaussians and the $P_k^{(i)}$ are copies of $P_k$ corresponding to various $k \times k$ minors of the $k+1 \times k+1$ matrix.
As the entries of the $k+1 \times k+1$ matrix are iid, the $x_i$ are independent of the $P_k^{(i)}$, and so for fixed values of $P_k^{(i)}$, we see that $P_{k+1}$ is distributed normally with mean zero and variance $\sum_{i=1}^{k+1} |P_k^{(i)}|^2$, which we can rewrite as
$$ P_{k+1} = (\sum_{i=1}^{k+1} |P_k^{(i)}|^2)^{1/2} \cdot N_{k+1}\qquad (1)$$
where $N_{k+1}$ is a standard normal random variable (independent of the $P_k^{(i)}$).
Now we come up against the key problem: the $P_k^{(i)}$ are identically distributed, but are not jointly independent, because there is a huge amount of overlap between the $k \times k$ minors. So, while heuristically one expects concentration of measure to kick in and make $(\sum_{i=1}^{k+1} |P_k^{(i)}|^2)^{1/2}$ more concentrated than any one of the $P_k^{(i)}$, we don't know how to prevent huge correlations from happening. In the worst case, all the $P_k^{(i)}$ are perfectly correlated to each other, and then (1) could become something more like
$$ P_{k+1} = (k+1)^{1/2} |P_k| \cdot N_{k+1}.$$
This multiplicative normal process would lead to $P_n$ to concentrate between $\sqrt{n!} \exp(-O(\sqrt{n}))$ and $\sqrt{n!} \exp(O(\sqrt{n}))$, as can be seen by taking logs and applying the central limit theorem.
But this worst case can't actually be the truth - among other things, it contradicts the second moment calculation. So there should be some way to prevent this correlation. Unfortunately, my paper with Van completely evades this issue - we try to get as far as we can just from the obvious fact that disjoint minors are independent from each other. This is why our bounds are significantly worse than $\exp(O(\sqrt{n}))$ from the truth.
As you say, the situation is much better for the determinant of a Gaussian iid matrix. Here, we can use the base-times-height formula to express the determinant as a product $\prod_{i=1}^n \hbox{dist}(X_i,V_i)$, where $X_i$ is the $i^{th}$ row and $V_i$ is the span of the first $i-1$ rows. With everything being Gaussian, $\hbox{dist}(X_i,V_i)^2$ is a chi-squared distribution, which is a martingale in $i$, and as a consequence one can get the determinant within $\exp(O(\sqrt{\log n}))$ of $\sqrt{(n-1)!}$, which would give you what you want. Unfortunately, there is nothing like the base-times-height formula for the permanent...
Finally, I am fairly certain that at the level of $n^{-\epsilon n}$ losses, one can replicate our paper in the Bernoulli case to the Gaussian case. I don't think the $n^{-\epsilon n}$ loss gets much better in that case, but the $\frac{1}{n^{0.1}}$ bound should improve substantially.
For every constant dimension $d$, (I.e. when $C = [0,1]^d$), the answer is asymptotically $\varepsilon_d(m) = \Theta\left(\frac{\log m}{m}\right)^{1/d}$. In the $\Theta$ notation I’m hiding constant factors that depend on $d$. This follows from the coupon collectors problem: let us partition a $[0,1]^d$ cube into $\varepsilon^{-d}$ sub-cubes of side-length $\varepsilon$, I.e. $[0,1]^d = \left( \bigcup_i [i \varepsilon, (i + 1) \varepsilon]\right)^d$.
Now, let us draw $m$ points uniformly at random from $[0, 1]^d$, and let’s keep track of which sub-cubes are being hit.
By coupon collectors, as soon as $m \gg \varepsilon^{-d} \log(\varepsilon^{-d})$, with high probability we will hit each sub-cube, and when $m \ll \varepsilon^{-d} \log(\varepsilon^{-d})$ with high probability we will miss at least one sub-cube.
Now it is enough to show the following two elementary geometric fact:
Any subset $S\subset [0, 1]^d$ which misses at least one sub-cube, has dispersion at least $\varepsilon/2$ (take the center of a missed cube as a witness $x$).
Any subset $S \subset [0,1]^d$ which hits every sub-cube has dispersion at most $\sqrt{d}\varepsilon$ (the diameter of a cube of side length $\varepsilon$).
Related: I recommend the great answer by Iosif Pinelis with much more precise calculation of the exact constant for the case where $d=1$.
Best Answer
The packing number of the cube with L2 balls of radius $\lambda$ is at least $ \lambda^{-n} C^n (n/2)! \approx \lambda^{-n} C^n n^{n/2} $. By the argument you had before you get the upper bound in the question is tight.