I have the following limit to take
$$\lim_{\epsilon\to0}\int\mathrm{d}x\mathrm{d}y\,\frac{g_\epsilon(x)g_\epsilon(y)}{4\pi|x-y|}e^{-\sqrt{\lambda}|x-y|} \overset{?}{=} \int\mathrm{d}x\mathrm{d}y\,\lim_{\epsilon\to0}\frac{g_\epsilon(x)g_\epsilon(y)}{4\pi|x-y|}e^{-\sqrt{\lambda}|x-y|}\tag{1} $$
where $$g_\epsilon(x) = \epsilon^{-3}g(\epsilon^{-1}x)\qquad g\in L^1(\mathbb{R}^3)\cap L^2(\mathbb{R}^3)\qquad\int\mathrm{d}x\,g(x)=1$$ is a $\delta$-approximating function.
I would like the result $(1)$ to hold since, under the appropriate change of variables, I would get $$\lim_{\epsilon\to0}\int\mathrm{d}x\mathrm{d}y\,\frac{g_\epsilon(x)g_\epsilon(y)}{4\pi|x-y|}e^{-\sqrt{\lambda}|x-y|} = \lim_{\epsilon\to 0}\int\mathrm{d}x\mathrm{d}y\,\frac{g(x)g(y)}{4\pi\epsilon|x-y|}e^{-\sqrt{\lambda}\epsilon|x-y|} \\ =\int\mathrm{d}x\mathrm{d}y\,\frac{g(x)g(y)}{4\pi\epsilon|x-y|}\left(1-\epsilon\sqrt{\lambda}|x-y|+o(\epsilon^2)\right)$$
How would one justify such a passage of the limit under the integral sign so to make the Tayor expansion of the exponential possible?
Just leave behind the fact that the result is divergent since this calculation has to do with a specific rinormalization condition.