Double integral with strange Change of Variables I am trying to compute the following integral,
$$\iint_{\mathbb{R}^2} \left(\frac{1-e^{-xy}}{xy}\right)^2 e^{-x^2-y^2}dxdy$$
First I tried substituting $x=r\cos{\theta}, y=r\sin{\theta}$ but it didn't really give me anything.
For the second try, I tried $u=x^2+y^2, v=xy$ but after computing the Jacobian, it got really messy, I couldn't really continue further.
Would there be a way to change the variables so that I would be able to compute this integral? 
 A: Unleash the infinite power of symmetry! The value of the integral is twice the value of the integral over the region $0\leq y\leq x$, and by setting $y=x t, dy = x\,dt$ we get
$$ I = 2 \iint_{(0,+\infty)^2} x\left(\frac{1-e^{-tx^2}}{tx^2}\right)^2 e^{-x^2-t^2 x^2}\,dx\,dt\stackrel{x\mapsto\sqrt{s}}{=} \iint_{(0,+\infty)^2}\left(\frac{1-e^{-ts}}{ts}\right)^2 e^{-s-t^2 s}\,ds\,dt$$
By Frullani's theorem we have $\int_{0}^{+\infty}\frac{e^{-au}-e^{-bu}}{u}\,du = \log\frac{b}{a}$ for any $a,b>0$, hence by integrating with respect to $s$ we have
$$ I = \int_{0}^{+\infty}\left[(1+t^2)\log\frac{1+t^2}{(1+t)^2}-2(1+t+t^2)\log\frac{1+t+t^2}{(1+t)^2}\right]\frac{dt}{t^2} $$
and by integration by parts
$$ I = \int_{0}^{+\infty}\left[2t\log\frac{1+t^2}{(1+t)^2}-2(1+2t)\log\frac{1+t+t^2}{(1+t)^2}\right]\frac{dt}{t}. $$
Now $1+t+t^2=\Phi_3(t)=\frac{1-t^3}{1-t}$ and both $\log(1\pm t^k)$ and $\frac{1}{t}\log(1\pm t^k)$ have manageable primitives in terms of $\log$ and $\text{Li}_2$. The final outcome is
$$ \iint_{(0,+\infty)^2} \left(\frac{1-e^{-xy}}{xy}\right)^2 e^{-x^2-y^2}dxdy = \color{blue}{\pi -\frac{2 \pi }{\sqrt{3}}+\frac{\pi ^2}{9}}. $$
A: $\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,}
 \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace}
 \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack}
 \newcommand{\dd}{\mathrm{d}}
 \newcommand{\ds}[1]{\displaystyle{#1}}
 \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,}
 \newcommand{\ic}{\mathrm{i}}
 \newcommand{\mc}[1]{\mathcal{#1}}
 \newcommand{\mrm}[1]{\mathrm{#1}}
 \newcommand{\pars}[1]{\left(\,{#1}\,\right)}
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,}
 \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}}
 \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$
\begin{align}
&\bbox[10px,#ffd]{\ds{\iint_{\mathbb{R}^2}
\pars{1 - \expo{-xy} \over xy}^{2}\expo{-x^{2} - y^{2}}\dd x\,\dd y}}
\\[5mm] = &\
\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\
\overbrace{\pars{\int_{0}^{1}\expo{-xya}\dd a}}
^{\ds{\expo{-xy} - 1 \over -xy}}\
\overbrace{\pars{\int_{0}^{1}\expo{-xyb}\dd b}}
^{\ds{\expo{-xy} - 1 \over -xy}}
\expo{-x^{2} - y^{2}}\dd x\,\dd y
\\[5mm] = &\
\int_{0}^{1}\int_{0}^{1}\
\overbrace{\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}
\exp\pars{-x^{2} - \bracks{a + b}xy - y^{2}}\dd x\,\dd y}
^{\ds{2\pi \over \root{4 - \pars{a + b}^{2}}}}\
\,\dd a\,\dd b
\\[5mm] = &\
2\pi\int_{0}^{1}\int_{0}^{1}{\dd a\,\dd b \over \root{4 - \pars{a + b}^{2}}} =
\bbx{{4 \over 3}\,\pi\pars{3 - 3\root{3} + \pi}} \approx 3.9603
\end{align}
