I would like to know a rather precise asymptotics of the sum
$$ S(x) = S_{\alpha,\beta}(x) := \sum_{y \in \Bbb{Z}^d} \frac{1}{(|y| + 1)^{\alpha}(|x-y| + 1)^{\beta}}$$
as $|x| \to \infty$. Here, $\alpha, \beta > 0$. I suspect that
$$ S(x) = \Theta_{\alpha,\beta}( |x|^{-(\alpha+\beta-d)} ) $$
for each fixed $\alpha, \beta$ with $\alpha+\beta > d$. However, I would like to know about this asymptotics in a better precision, such as an expansion of the form
$$ S(x) = \frac{a_0}{|x|^{\alpha+\beta-d}} + \frac{a_1}{|x|^{\alpha+\beta+1-d}} + \cdots $$
as $|x| \to \infty$, with constants $a_0, a_1, \cdots$ depending only on $d$, $\alpha$ and $\beta$. Even a keyword would be appreciated.
Remark. It is easy to see that for some constants $c, C > 0$ depending only on the dimension $d$, we have
$$ c^{\alpha+\beta} I(x) \leq S(x) \leq C^{\alpha+\beta} I(x), $$
where $I$ is defined as
$$ I(x) = I_{\alpha,\beta}(x) := \int_{\Bbb{R}^d} \frac{1}{(|y|^2 + 1)^{\alpha/2}(|x-y|^2 + 1)^{\beta/2}} \, dy $$
This one is easier to deal with, and in fact we have
$$ I(x) = \frac{2\pi^{d/2}\Gamma(\frac{\alpha+\beta-d}{2})}{\Gamma(\frac{\alpha}{2})\Gamma(\frac{\beta}{2})} \int_{0}^{\pi/2} \frac{\cos^{\alpha-1}\theta\sin^{\beta-1}\theta}{(|x|^2 \cos^2\theta\sin^2\theta + 1)^{(\alpha+\beta-d)/2}} \, d\theta. $$
This might be used to give a good guess on the leading coefficient of the expansion, but I am not confident about this. Or, recognizing $S(x)$ as a convolution might be helpful in conjunction with Fourier analysis technique, but I am even less confident about this.