While I was trying to write a PDE methods-centric answer to this question I ran into the problem in the title. The original problem was (with new notation) to calculate
$$\int_{\Bbb{R}^3\times\Bbb{R}^3} dx\:dx'\:\frac{e^{-2a(x+x')}}{|\mathbf{x-x'}|}$$
As a commenter noted in the original post, this quantity is homogeneous in $a$ with degree $-5$, so we will remove $a$ from consideration entirely. The integral also is locally integrable everywhere. Denoting
$$f(\mathbf{x}) = \int_{\Bbb{R}^3}dx'\:\frac{e^{-2x'}}{|\mathbf{x-x'}|}$$
I noticed that
$$\Delta f = -4\pi e^{-2x}$$
So we were being asked to compute
$$-\frac{1}{4\pi}\int_{\Bbb{R}^3}dx\:f\Delta f$$
This was a nice form to attempt integration by parts. Choosing $\nabla f$ to be purely radial, we obtain that
$$\nabla f = \frac{-4\pi}{x^2}\hat{\mathbf{x}}\int dx\:x^2e^{-2x} = \pi e^{-2x}\left(2+\frac{2}{x}+\frac{1}{x^2}\right)\hat{\mathbf{x}}$$
Already we can see there is some sort of pole at $0$, so I tried to remove a ball of size $\epsilon$ and continue integrating by parts
$$-\frac{1}{4\pi}\int_{\Bbb{R}^3}dx\:f\Delta f = \lim_{\epsilon\to0^+}-\frac{1}{4\pi}\int_{\Bbb{R}^3-B_\epsilon(0)}dx\:f\Delta f$$
$$= \lim_{\epsilon\to0^+}+\frac{1}{4\pi}\int_{\partial B_\epsilon(0)}d\sigma\:f\left(\nabla f \cdot \hat{\mathbf{x}}\right)-\frac{1}{4\pi}\int_{\Bbb{R}^3-B_\epsilon(0)}dx\:|\nabla f|^2$$
Computing the boundary integral is easier since the integrand purely radial $$\lim_{\epsilon\to0^+}\frac{1}{4\pi}f(\epsilon) \left(\nabla f(\epsilon)\cdot\hat{x}\right)\int_{S^2}\epsilon^2\:d\Omega = \lim_{\epsilon\to0^+}\pi f(\epsilon) e^{-2\epsilon}\left(2\epsilon^2+2\epsilon+1\right) = 2\pi^2$$
wherein lies the problem. Usually when the other integral has a singularity in it, the boundary integral has some sort of equal and oppositely valued singularity that would cancel out and give me the regularized value of the original, well-defined integral. Where did I go wrong, why does the singularity not cancel out?
For completeness, computing the volume integral we obtain
$$\lim_{\epsilon\to0^+}\pi^2\int_{\epsilon}^\infty dx \:e^{-4x}\left(2x+2+\frac{1}{x}\right)^2 = \frac{21\pi^2}{8} + \lim_{\epsilon\to0^+}\frac{\pi^2 e^{-4\epsilon}}{\epsilon}$$
Now one interesting thing is that if I was off by a factor of epsilon in the boundary integral, then we get the resultant expression
$$4\pi^2 -\frac{21\pi^2}{8} + \lim_{\epsilon\to0^+}\frac{\pi\left(f(\epsilon) e^{-2\epsilon}-\pi e^{-4\epsilon}\right)}{\epsilon}$$
EDIT: I found my mistake, I will write up an answer later.
