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?
On
$\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}
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}}. $$