How to evaluate $$F(\Lambda)=\int_0^{\infty}\int_0^{\infty}\left(\frac{e^{-\Lambda x}-e^{-\Lambda y}}{x-y}\right)^2 \,dx\,dy$$ where $\Lambda$ is a positive real number?
I tried evaluating the innermost integral first,
$\displaystyle\int_{0}^{\infty}\left(\frac{e^{-\Lambda x}-e^{-\Lambda y}}{x-y}\right)^2 dx$
and found the solution is terms of the exponential integral $E_1(u)$ and I am not sure if that would help at all.
First, the change of variables $u=\Lambda x$, $v=\Lambda y$, shows that $F(\Lambda)$ does not depend on $\Lambda>0$ hence one can simply compute $F(1)$.
Second, decomposing the domain $x>0$, $y>0$, into its $x>y$ and $y<x$ parts and using the symmetry of the integrand, one sees that $F(1)=2I$, where $$I=\iint_{0<x<y}\left(\frac{e^{-x}-e^{-y}}{x-y}\right)^2dxdy$$ Third, the change of variable $y=x+z$ shows that $$I=\iint_{x>0,z>0}\left(\frac{e^{-x}(1-e^{-z})}z\right)^2dxdz$$ and now the domain of integration is a product and the integrand can be factored hence $I=J\cdot K$ where $$J=\int_0^\infty e^{-2x}dx=\frac12$$ and $$K=\int_0^\infty\left(\frac{1-e^{-z}}z\right)^2dz$$ Fourth, to estimate $K$, one can use the identity $$1-e^{-z}=\int_0^ze^{-u}du$$ which allows to rewrite $K$ as $$K=\int_0^\infty\frac{dz}{z^2}\int_0^ze^{-u}du\int_0^ze^{-v}dv=\iint_{(0,\infty)^2}e^{-u}e^{-v}dudv\int_{\max\{u,v\}}^\infty\frac{dz}{z^2}$$ that is, $$K=\iint_{(0,\infty)^2}e^{-u}e^{-v}\frac{dudv}{\max\{u,v\}}$$ Using symmetry again, $K=2L$, where $$L=\iint_{0<u<v}e^{-u}e^{-v}\frac{dudv}v=\int_0^\infty(1-e^{-v})e^{-v}\frac{dv}v$$ and finally, $$L=\int_0^\infty\frac{e^{-v}-e^{-2v}}vdv$$ Thus, $L$ is a Frullani integral hence $$L=\log\left(\frac21\right)$$ and $$F(\Lambda)=2\cdot\frac12\cdot2\cdot\log2=2\cdot\log2$$ ...unless I missed some factor $2$ en route. :-)
By symmetry, for any $\Lambda>0$ we have
$$\begin{eqnarray*} F(\Lambda)=\iint_{(0,+\infty)}\left(\frac{e^{-\Lambda x}-e^{-\Lambda y}}{x-y}\right)^2\,dx\,dy&=&2\int_{0}^{+\infty}\int_{0}^{x}\left(\frac{e^{-\Lambda x}-e^{-\Lambda y}}{x-y}\right)^2\,dy\,dx \\ &\stackrel{y\mapsto xz}{=}&2\int_{0}^{+\infty}\int_{0}^{1}\frac{\left(e^{-\Lambda x}-e^{-\Lambda xz}\right)^2}{x(1-z)^2}\,dz\,dx \\ &\stackrel{x\mapsto w/\Lambda}{=}&2\int_{0}^{1}\int_{0}^{+\infty}\frac{\left(e^{-w}-e^{-wz}\right)^2}{w(1-z)^2}\,dw\,dz \\ &\stackrel{\text{Frullani}}{=}&2\int_{0}^{1}\frac{2\log(1+z)-\log(4z)}{(1-z)^2}\,dz\end{eqnarray*}$$ hence $F(\Lambda)$ is constant and it equals $$\begin{eqnarray*} F(\Lambda)&=&2\int_{0}^{1}\frac{2\log\left(1-\frac{u}{2}\right)-\log(1-u)}{u^2}\,du\\&=&2\int_{0}^{1}\sum_{n\geq 2}\left(\frac{1}{n}-\frac{2}{n 2^{n}}\right)u^{n-2}\,du\\&=& 2\sum_{n\geq 2}\frac{1}{n(n-1)}\left(1-\frac{2}{2^n}\right)=\color{blue}{2\log 2}. \end{eqnarray*}$$
