Prove or disprove the equality of these two integrals Let $\alpha$ be an arbitrary positive real number in:
$$
F_1 = \int_0^1 x^2 \left[ \int_{-1}^{+1} \frac{e^{-\alpha\sqrt{1+x^2+2xy}}(xy+1)}{(1+x^2+2xy)^{3/2}}dy\right]dx $$ $$
F_2 = \int_1^\infty x^2 \left[ \int_{-1}^{+1} \frac{e^{-\alpha\sqrt{1+x^2-2xy}}(xy-1)}{(1+x^2-2xy)^{3/2}}dy\right]dx
$$
Prove or disprove that $F_1 = F_2$.
Source of the problem are the equations (3) and (4) in
A Paradox of Newtonian Gravitation and Laplace’s Solution
by Amitabha Ghosh and Ujjal Dey.
They have done already numerical experiments that seem to confirm equality.
Quote: an analytical proof showing that F1 and F2 are exactly equal will be an interesting mathematical exercise.
And that's it. I have no idea how to proceed.
Progress so far
I promised myself not to do numerical experiments.
But the outcome of the inner integral - the one between square brackets - is indeed terrible.
So what else would be an option?
With MAPLE 8 some values in the publication can be reproduced.

for k from 0 to 10 do
alpha := k*0.1;
g(x,alpha) := int(exp(-alpha*sqrt(1+x^2+2*x*y))*(x*y+1)/(1+x^2+2*x*y)^(3/2),y=-1..1);
F1 := evalf(int(g(x,alpha)*x^2,x=0..1));
F2 := evalf(int(-g(x,alpha)*x^2,x=1..10^3)); 
end do;

The special case $\alpha = 0$ gives $\,F_1=\frac{2}{3}\,$ and $\,F_2=0\,$ exactly.
So it appears that nearby $\alpha = 0$ the integrals must be unequal.
But MAPLE keeps calculating endlessly for those low values and I had to manually stop the program.

Feynman trick.
$F_1$ and $F_2$ are both a function of $\alpha$ . We can take derivatives under the integral sign and see what happens near $\alpha=0$.
$$
\left.\frac{dF_1}{d\alpha}\right|_{\alpha=0} =
- \int_0^1 x^2 \left[ \int_{-1}^{+1} \frac{xy+1}{1+x^2+2xy}dy\right]dx = -\frac{1}{2} \\
\Longrightarrow \quad dF_1 = -\frac{1}{2}d\alpha
$$ $$
\left.\frac{dF_2}{d\alpha}\right|_{\alpha=0} =
- \int_1^\infty x^2 \left[ \int_{-1}^{+1} \frac{xy-1}{1+x^2-2xy}dy\right]dx = \infty \\
\Longrightarrow \quad dF_2 = \infty\,d\alpha
$$
Leading to the following heuristics.
$F_1(\alpha)$ is somewhat decreasing from $F_1(0)=2/3$ to lower values, but the increase in $F_2(0)=0$ is infinitely large at that place.
Based upon this, it's impossible to keep up appearances )-:
we can actually say nothing yet whether the outcome is $\,F_1(\alpha) = F_2(\alpha)\,$ for all $\,\alpha \gt 0\,$.
 A: Proof sketch:
\begin{align*}
\frac{{{\rm d}^3 F_1 }}{{{\rm d}\alpha ^3 }} & = -\int_0^1 {x^2 \int_{ - 1}^1 {{\rm e}^{ - \alpha \sqrt {1 + x^2  + 2xy} } (xy + 1){\rm d}y} \,{\rm d}x} \\ & = -\frac{{{\rm e}^{ - 2\alpha } (4\alpha ^4  + 14\alpha ^3  + 27\alpha ^2  + 30\alpha  + 15)}}{{\alpha ^6 }} - \frac{{3(\alpha ^2  - 5)}}{{\alpha ^6 }}
\end{align*}
and
\begin{align*}
\frac{{{\rm d}^3 F_2 }}{{{\rm d}\alpha ^3 }} & = -\int_1^{ + \infty } {x^2 \int_{ - 1}^1 {{\rm e}^{ - \alpha \sqrt {1 + x^2  - 2xy} } (xy - 1){\rm d}y} \,{\rm d}x} \\ &  = -\frac{{{\rm e}^{ - 2\alpha } (4\alpha ^4  + 14\alpha ^3  + 27\alpha ^2  + 30\alpha  + 15)}}{{\alpha ^6 }} - \frac{{3(\alpha ^2  - 5)}}{{\alpha ^6 }}.
\end{align*}
Thus,
$$
\frac{{{\rm d}^3 (F_1  - F_2 )}}{{{\rm d}\alpha ^3 }} = 0,
$$
i.e., $F_1-F_2$ is a polynomial of degree at most $2$. By changing the variable $y\to -y$ in the definition of $F_2$, we see that
$$
F_1  - F_2  = \int_0^{ + \infty } {x^2 \int_{ - 1}^1 {\frac{{{\rm e}^{ - \alpha \sqrt {1 + x^2  + 2xy} } ( xy+1)}}{{(1 + x^2  + 2xy)^{3/2} }}{\rm d}y}\, {\rm d}x} .
$$
The right-hand side tends to $0$ as $\alpha\to +\infty$, i.e., the polynomial $F_1-F_2$ must be the zero polynomial.
A: Define
$$ F(\alpha) = \int_0^{\infty} \int_0^{\pi} {\rm e}^{-\alpha r} \sin \theta \, \cos \theta \, {\rm d}\theta \, {\rm d}r $$
Integrating yields
$$ F(\alpha) = -\frac{1}{\alpha} \left. {\rm e}^{-\alpha r} \right|_{r=0}^\infty \; \cdot \; \frac{1}{2} \left. \sin^2 \theta \right|_{\theta=0}^\pi = 0 \quad \text{for all } \alpha>0, $$
but does not exist for $\alpha = 0$.
Map the function using $u = -r \cos \theta$, $v = r \sin \theta$ and the determinant of the Jacobian
$$ \left| J \right| = \left|
\begin{matrix}
  \frac{\partial u}{\partial \theta} = r \sin \theta & \frac{\partial u}{\partial r} = -\cos \theta \\
  \frac{\partial v}{\partial \theta} = r \cos \theta & \frac{\partial v}{\partial r} = \sin \theta
\end{matrix} \right| = r$$
This gives
$$ F(\alpha) =
\int_0^\infty \int_0^\pi  \frac{{\rm e}^{-\alpha r} \sin\theta \cos\theta}{r} \, r \, {\rm d}\theta \, {\rm d}r =
\int_0^\infty \int_{-\infty}^\infty  \frac{{\rm e}^{-\alpha \sqrt{u^2 +v^2}} uv}{(u^2 +v^2)^{3/2}} \, {\rm d}u \, {\rm d}v 
$$
The integration limits are determined from this map:
map from r, $\theta$ to u, v
The origin is shifted by one, $u \rightarrow u +1$, to give
$$ F(\alpha) =
\int_0^\infty \int_{-\infty}^\infty  \frac{{\rm e}^{-\alpha \sqrt{u^2 +2u +1 +v^2}} (u +1)v}{(u^2 +2u +1 +v^2)^{3/2}} \, {\rm d}u \, {\rm d}v 
$$
Map again using $y = -\cos\theta' = u/\sqrt{u^2 +v^2}$, $x = r' = \sqrt{u^2 +v^2}$, $u = xy$ and the determinant of the Jacobian
$$ \left| J \right| = \left|
\begin{matrix}
  \frac{\partial y}{\partial u} = 1/x -u^2/x^3 & \frac{\partial y}{\partial v} = -uv/x^3 \\
  \frac{\partial x}{\partial u} = u/x & \frac{\partial x}{\partial v} = v/x
\end{matrix} \right| = v/x^2$$
$$ F(\alpha) =
\int_0^\infty \int_{-\infty}^\infty  \frac{x^2}{v} \frac{{\rm e}^{-\alpha \sqrt{u^2 +2u +1 +v^2}} (u +1)v}{(u^2 +2u +1 +v^2)^{3/2}} \, \frac{v}{x^2} \, {\rm d}u \, {\rm d}v $$
$$ = \int_0^\infty x^2 \int_{-1}^1 \frac{{\rm e}^{-\alpha \sqrt{x^2 +2xy +1}} (xy +1)}{(x^2 +2xy +1)^{3/2}} \, {\rm d}y \, {\rm d}x 
$$
$x$ acts as a radius, $y$ is the fraction of $x$ along $u$.  The integration limits are determined from this map:
map from u, v to x, y
This is the same function given above by @Gary, so we have
$$ F(\alpha) = F_1 - F_2 = \int_0^{\infty} x^2 \int_{-1}^1 \frac{{\rm e}^{-\alpha \sqrt{1 +x^2 +2xy}} (xy +1)} {(1 +x^2 +2xy)^{3/2}} \, {\rm d}y \, {\rm d}x = 0 \quad \text{for all } \alpha>0, $$
but does not exist for $\alpha = 0$.  Therefore $F_1 = F_2$ for all $\alpha > 0$.
