Hi Terry. I am going to be somewhat lazy and sketch Here is a construction sketch of an argument that could should yield a presentable bound (meaning the estimate $\epsilon_d$ 1/\varepsilon_d = O(d^{3/2})$. That is polynomial in $d$) without explicitly calculating anythingnot your full conjecture, but it is getting there. I haven't checked every detail of the more delicate second half of the argument, but it should work.
First, let me rephrase the question with only polynomial overhead. The contrapositive of your statement is that if the integer lattice is a lattice packing of $K^*$, then it is a lattice covering of $\epsilon_d \varepsilon_d K$ (up to a factor of 2). You can replace $K^*$ by the largest inscribed ellipsoid $E$, and then John's theorem says that $\sqrt{d} \cdot K \supseteq E^*$. Then, after a linear transformation, we can say that $E = E^*$ is the round unit ball $B = B_1(0)$, i.e., the $\ell^2$ unit ball. The contrapositive hypothesis is that lattice $\Lambda$ is a lattice packing of $B$, i.e., a unit sphere packing. You want to bound the sphere covering radius of $\Lambda^*$.
By taking a Fourier transform on $L^2(\mathbb{R}^d)$, what you know is that $\Lambda^*$ is a "1-design" in the sense of Delsarte. (This is with a scaled Fourier transform so that none of the geometric lengths have factors of $\pi$.) If $f$ is a sufficiently regular function whose Fourier transform is supported on the unit ball in $\mathbb{R}^d$, then the integral of $f$ on $\mathbb{R}^d$ (or average or some regularized integral) equals its sum (or regularized sum ) on $\Lambda^*$. I have a paper on $t$-designs [arXiv:math/0405366] which suggests a method that could give you a covering radius, although in my paper it was the analytically easier case of a compact domain. The idea is to find an $f(x) = f(||x||_2)$ whose integral is positive, yet which is non-positive for $||x||_2 > c$, and whose Fourier transform satisfies the support condition. Then $\Lambda^*$ must have a lattice point in the ball $B_c(0)$, and indeed in $B_c(p)$ for any $p$. I call this the "positive island" method.
If
In one stage of the construction argument in my papercan be generalized , I made a positive island function on the manifold $\mathbb{C}P^d$ of the form $P(z)/(a-z)$, where $P$ is a Jacobi polynomial (with indices suppressed) and $a$ is its last zero. This was for $t$-designs on $\mathbb{C}P^n$ (and later the simplex, using the moment map). There is a similar formula for ordinary spherical $t$-designs on $S^d$. To make this expression relevant to your casequestion, then you would let can take the limit as the degree $f(x)$ look something like t \to \infty$ and the positive island shrinks to a point. In this limit, the geometry of the manifold becomes approximately Euclidean and so approximates your question. The island function $J(||x||_2)^2/(c^2-||x||^2_2)$, P(z)/(a-z)$ limits to the function$$f(x) = \frac{J(||x||_2)^2}{c^2-||x||^2_2}$$on $\mathbb{R}^d$, where $J(r)$ is a suitable hyperspherical Bessel function and $c$ is its first zero.
So, although at for roundabout reasons, what should happen is that the moment integral of this $f(x)$ vanishes and its Fourier transform has the right support property. If you perturb $f$ slightly, you can make its integral positive. I imagine that there is a direct argument for the properties of this $f$, but I did not work at it. However, I did check a few cases numerically with Maple and it seems like to work. For instance, the integral on $\mathbb{R}^3$ of$$f(x) = \frac{(\sin x)^2}{x^2(\pi^2-x^2)}$$vanishes.
Now, the first zero of the first hyperspherical Bessel function in $\mathbb{R}^d$ is the same as the first zero $j_{(d-2)/2,1}$ of an ordinary Bessel function. I believe that this number is $O(d)$. This would yield the estimate $1/\varepsilon_d = O(d^{3/2})$, since you also get a naive formulafactor of $O(d^{1/2})$ from John's theorem.

