Consider the one-parameter family of functions $f_\epsilon:\mathbb{R}^2\to \mathbb{C}$ given by $$f_\epsilon(x,y)=\frac{1}{(\alpha^2\epsilon^2+xy)^\beta},$$
where $\alpha\in \mathbb{R}$ and $\beta\in \mathbb{C}$ are fixed numbers. We can even assume that $\beta \in (0,+\infty)$ if it makes things better. I want to evaluate the expansion of $f_\epsilon(x,y)$ in powers of $\epsilon$ in the distributional sense. Indeed, the very similar case where we instead consider $g_\epsilon : \mathbb{R}^n\to \mathbb{C}$
$$g_\epsilon(x)=\frac{1}{(\alpha^2\epsilon^2+|x|^2)^\beta}$$
is well-understood: we integrate $g_\epsilon$ against a test function $\phi\in C_0^\infty(\mathbb{R}^n)$ and break the integration domain in two regions: $|x|<\zeta$ and $|x|>\zeta$. In the first region we change variables setting $x = \epsilon y$ and then we see that the $\epsilon$ expansion is just the Taylor expansion of the test function $\phi$ (which further reveals that the distributional expansion only makes sense in the subspace of test functions which are analytic at $x = 0$), whereas in the second region we can directly Taylor expand $g_\epsilon$ in powers of $\epsilon$.
Now, I tried doing the same for $f_\epsilon(x,y)$ but I'm confused on how to proceed. Now I have
$$\langle f_\epsilon,\phi\rangle=\int_{-\infty}^\infty dx \int_{-\infty}^\infty dy \frac{\phi(x,y)}{(\alpha^2\epsilon^2+xy)^\beta}.$$
Now break the $x$ integrals into $|x|<\zeta$ and $|x|>\zeta$ regions
$$\langle f_\epsilon,\phi\rangle=\int_{|x|<\zeta} dx \int_{-\infty}^\infty dy \frac{\phi(x,y)}{(\alpha^2\epsilon^2+xy)^\beta}+\int_{|x|>\zeta} dx \int_{-\infty}^\infty dy \frac{\phi(x,y)}{(\alpha^2\epsilon^2+xy)^\beta}.$$
In the first it looks like we want to set $x = \alpha^2\epsilon^2 u$ so that we can factor $\epsilon$
$$\langle f_\epsilon,\phi\rangle=\int_{|u|<\zeta/\epsilon^2} du \int_{-\infty}^\infty dy \frac{\phi(\epsilon^2 u,y)}{\alpha^{2\beta}\epsilon^{2\beta-2}(1+uy)^\beta}+\int_{|x|>\zeta} dx \int_{-\infty}^\infty dy \frac{\phi(x,y)}{(\alpha^2\epsilon^2+xy)^\beta}.$$
Now I guess I would Taylor expand $\phi(\epsilon^2 u,y) \simeq \phi(0,y)$ and use $\epsilon\to 0$ to also replace the $u$ integral by an integral over the whole real line. Then I should be able to evaluate the $u$ integral:
$$\int_{-\infty}^\infty \dfrac{du}{(1+uy)^\beta},$$
but this integral seems to give either zero or diverges depending on $\beta$. I'm confused. The thing I expect for $\lim_{\epsilon\to 0}f_\epsilon$ is a contribution proportional to $\frac{1}{y^\beta}\delta(x)$ and one proportional o $\frac{1}{x^\beta}\delta(y)$, but I don't see how this can actually happen.
So how do we tackle this distributional limit in the correct, rigorous way?
