Let us consider the vacuum energy of a scalar field in $d+1$-dimensional spacetime. We have the integral $$I=\int\frac{d^d k}{(2\pi)^d}\frac{E_k}{2},$$ where $E_k=\sqrt{k^2+m^2}$ and $k$ is a $d$-dimensional vector. Now this integral is apparently divergent. We will employ the dimensional regularzation. We have a useful formula (see arXiv:1701.01554 , page 20) $$\Phi(m,d,A)=\int\frac{d^dk}{(2\pi)^d}\frac{1}{(k^2+m^2)^A}=\frac{1}{(4\pi)^{d/2}}\frac{\Gamma(A-d/2)}{\Gamma(A)}\frac{1}{(m^2)^{A-d/2}}.$$
If we consider a 4-dimensional spacetime, then we take $d=3-2\epsilon$ and the integral has a contribution proportional to $1/\epsilon$ which comes from $\Gamma(-1/2-3/2+\epsilon)=\Gamma(-2+\epsilon)$, indicating the divergence for $d=3$. Now, if we take $d=4-2\epsilon$, we have $$\Gamma(-5/2+\epsilon)=\frac{\Gamma(1/2+\epsilon)}{(-5/2+\epsilon)(-3/2+\epsilon)(-1/2+\epsilon)}=\frac{\Gamma(1/2)+\Gamma(1/2)'\epsilon}{(-5/2+\epsilon)(-3/2+\epsilon)(-1/2+\epsilon)}.$$ The above expression, however, has no pole for $\epsilon\rightarrow 0$. Why isn't the integral divergent? Shouldn't one expect worse behavior for convergence with larger $d$? Is there anything wrong with the above mathematics or a physical reason to explain it?