$$ \int _\infty ^\infty\int _\infty ^\infty \frac {(z_1 -z_2)^2}{(k+(z_1 -z_2)^2)^2} \exp(-z_1^2 -z_2^2)dz_1 dz_2$$
I have been tried substituting $z_1-z_2$ by y, then using some substitution like $x=\frac {\exp(-y)^2}{(1+y^2)}$, but things have not worked out. Is there a better way?
By letting $\frac{z_1+z_2}{\sqrt{2}}=u, \frac{z_1-z_2}{\sqrt{2}}=v$ the given integral takes the form $$ \iint_{\mathbb{R}^2}\frac{2v^2}{(k+2v^2)^2}e^{-u^2-v^2}\,du\,dv =\sqrt{\pi}\int_{0}^{+\infty}\frac{4v}{(k+2v^2)^2}\cdot ve^{-v^2}dv$$ and by integration by parts the RHS equals $$ \sqrt{\pi}\int_{0}^{+\infty}\frac{1-2v^2}{k+2v^2}e^{-v^2}\,dv =\frac{\pi}{2k}-\sqrt{\pi}\left(1+\frac{1}{k}\right)\int_{0}^{+\infty}\frac{2v^2}{k+2v^2}e^{-v^2}\,dv$$ where the last term is non-elementary (it depends on the error function), but by the Cauchy-Schwarz inequality $$\int_{0}^{+\infty}\frac{2v^2}{k+2v^2}e^{-v^2}\,dv\leq\sqrt{\int_{0}^{+\infty}\left(\frac{2v}{k+2v^2}\right)^2\,dv\int_{0}^{+\infty}\left(ve^{-v^2}\right)^2\,dv}=\frac{\pi^{3/4}}{2^{5/2}k^{1/4}}.$$