$$\text{Compute} \qquad \int_{0}^{\infty}e^{-2x}\operatorname{Ei}^{2}(x)dx $$ $$\text{where} \qquad\operatorname{Ei}x=-\int_{-x}^\infty\frac{e^{-t}}tdt.\tag1$$
My attempt :
$\operatorname{Ei}^{2}( x)=\int_{-x}^\infty\int_{-x}^\infty\frac{e^{-(t+y)}}{ty}dydt.\tag2$
We obtain:
$\displaystyle\int_{0}^\infty\operatorname{Ei}^{2}(x)e^{-2x}dx=\int_0^\infty\int_{-x}^\infty\int_{-x}^\infty\frac{e^{-(2+t+y)}}{ty}dydtdx.\tag3$
Now how I can complete this work ?
For positive real values of $x$, the exponential integral $\operatorname{Ei}{\left(x\right)}$ is defined as the Cauchy principal value integral
$$\operatorname{Ei}{\left(x\right)}:=-\mathcal{P}\int_{-x}^{\infty}\mathrm{d}t\,\frac{\exp{\left(-t\right)}}{t}=-\lim_{\varepsilon\to0^{+}}\left[\int_{-x}^{-\varepsilon}\mathrm{d}t\,\frac{\exp{\left(-t\right)}}{t}+\int_{+\varepsilon}^{\infty}\mathrm{d}t\,\frac{\exp{\left(-t\right)}}{t}\right];~~~\small{x\in\mathbb{R}_{>0}}.$$
Suppose $x\in\mathbb{R}_{>0}$. Starting from the Cauchy principal value definition of $\operatorname{Ei}{\left(x\right)}$, an alternative representation in terms of an ordinary convergent improper integral can be found:
$$\begin{align} \operatorname{Ei}{\left(x\right)} &=-\lim_{\varepsilon\to0^{+}}\left[\int_{-x}^{-\varepsilon}\mathrm{d}t\,\frac{\exp{\left(-t\right)}}{t}+\int_{+\varepsilon}^{\infty}\mathrm{d}t\,\frac{\exp{\left(-t\right)}}{t}\right]\\ &=\gamma+\ln{\left(x\right)}+\int_{0}^{1}\mathrm{d}t\,\frac{\exp{\left(xt\right)}-1}{t},\\ \end{align}$$
where $\gamma$ here denotes the Euler-Mascheroni constant, whose value for our purposes here can be defined by the definite integral
$$\gamma:=-\int_{0}^{\infty}\mathrm{d}x\,\exp{\left(-x\right)}\ln{\left(x\right)}.$$
The following improper integrals can be shown to have a value expressed in terms of $\gamma$:
$$-z\int_{0}^{\infty}\mathrm{d}x\,\exp{\left(-zx\right)}\ln{\left(x\right)}=\gamma+\ln{\left(z\right)};~~~\small{z\in\mathbb{R}_{>0}},$$
and
$$\int_{0}^{\infty}\mathrm{d}x\,\exp{\left(-x\right)}\ln^{2}{\left(x\right)}=\gamma^{2}+\zeta{\left(2\right)}.$$
Now, define the function $\mathcal{I}:\left[2,\infty\right)\rightarrow\mathbb{R}$ via the definite integral
$$\begin{align} \mathcal{I}{\left(a\right)} &:=\int_{0}^{\infty}\mathrm{d}x\,\exp{\left(-ax\right)}\left[\operatorname{Ei}{\left(x\right)}\right]^{2}.\\ \end{align}$$
Given $x\in\mathbb{R}_{>0}$, the square of the exponential integral can be expanded as
$$\begin{align} \left[\operatorname{Ei}{\left(x\right)}\right]^{2} &=\left[\gamma+\ln{\left(x\right)}+\int_{0}^{1}\mathrm{d}t\,\frac{\exp{\left(xt\right)}-1}{t}\right]^{2}\\ &=\left[\gamma+\ln{\left(x\right)}\right]^{2}+2\left[\gamma+\ln{\left(x\right)}\right]\int_{0}^{1}\mathrm{d}t\,\frac{\exp{\left(xt\right)}-1}{t}+\left[\int_{0}^{1}\mathrm{d}t\,\frac{\exp{\left(xt\right)}-1}{t}\right]^{2}\\ &=\gamma^{2}+2\gamma\ln{\left(x\right)}+\ln^{2}{\left(x\right)}\\ &~~~~~+2\gamma\int_{0}^{1}\mathrm{d}t\,\frac{\exp{\left(xt\right)}-1}{t}+2\ln{\left(x\right)}\int_{0}^{1}\mathrm{d}t\,\frac{\exp{\left(xt\right)}-1}{t}\\ &~~~~~+\int_{0}^{1}\mathrm{d}t\,\frac{\exp{\left(xt\right)}-1}{t}\int_{0}^{1}\mathrm{d}u\,\frac{\exp{\left(xu\right)}-1}{u}\\ &=\gamma^{2}+2\gamma\ln{\left(x\right)}+\ln^{2}{\left(x\right)}\\ &~~~~~+2\gamma\int_{0}^{1}\mathrm{d}t\,\frac{\exp{\left(xt\right)}-1}{t}+2\int_{0}^{1}\mathrm{d}t\,\frac{\left[\exp{\left(xt\right)}-1\right]\ln{\left(x\right)}}{t}\\ &~~~~~+\int_{0}^{1}\mathrm{d}t\int_{0}^{1}\mathrm{d}u\,\frac{\left[\exp{\left(xt\right)}-1\right]\left[\exp{\left(xu\right)}-1\right]}{tu}.\\ \end{align}$$
Assuming $a\in\left(2,\infty\right)$, we then find
$$\begin{align} \mathcal{I}{\left(a\right)} &=\int_{0}^{\infty}\mathrm{d}x\,\exp{\left(-ax\right)}\left[\operatorname{Ei}{\left(x\right)}\right]^{2}\\ &=\int_{0}^{\infty}\mathrm{d}x\,\exp{\left(-ax\right)}\bigg{[}\gamma^{2}+2\gamma\ln{\left(x\right)}+\ln^{2}{\left(x\right)}\\ &~~~~~+2\gamma\int_{0}^{1}\mathrm{d}t\,\frac{\exp{\left(xt\right)}-1}{t}+2\int_{0}^{1}\mathrm{d}t\,\frac{\left[\exp{\left(xt\right)}-1\right]\ln{\left(x\right)}}{t}\\ &~~~~~+\int_{0}^{1}\mathrm{d}t\int_{0}^{1}\mathrm{d}u\,\frac{\left[\exp{\left(xt\right)}-1\right]\left[\exp{\left(xu\right)}-1\right]}{tu}\bigg{]}\\ &=\gamma^{2}\int_{0}^{\infty}\mathrm{d}x\,\exp{\left(-ax\right)}+2\gamma\int_{0}^{\infty}\mathrm{d}x\,\exp{\left(-ax\right)}\ln{\left(x\right)}+\int_{0}^{\infty}\mathrm{d}x\,\exp{\left(-ax\right)}\ln^{2}{\left(x\right)}\\ &~~~~~+2\gamma\int_{0}^{\infty}\mathrm{d}x\,\exp{\left(-ax\right)}\int_{0}^{1}\mathrm{d}t\,\frac{\exp{\left(xt\right)}-1}{t}\\ &~~~~~+2\int_{0}^{\infty}\mathrm{d}x\,\exp{\left(-ax\right)}\int_{0}^{1}\mathrm{d}t\,\frac{\left[\exp{\left(xt\right)}-1\right]\ln{\left(x\right)}}{t}\\ &~~~~~+\int_{0}^{\infty}\mathrm{d}x\,\exp{\left(-ax\right)}\int_{0}^{1}\mathrm{d}t\int_{0}^{1}\mathrm{d}u\,\frac{\left[\exp{\left(xt\right)}-1\right]\left[\exp{\left(xu\right)}-1\right]}{tu}\\ &=\frac{\gamma^{2}}{a}\int_{0}^{\infty}\mathrm{d}y\,\exp{\left(-y\right)}+\frac{2\gamma}{a}\int_{0}^{\infty}\mathrm{d}y\,\exp{\left(-y\right)}\ln{\left(\frac{y}{a}\right)}\\ &~~~~~+\frac{1}{a}\int_{0}^{\infty}\mathrm{d}y\,\exp{\left(-y\right)}\ln^{2}{\left(\frac{y}{a}\right)};~~~\small{\left[x=\frac{y}{a}\right]}\\ &~~~~~+2\gamma\int_{0}^{\infty}\mathrm{d}x\int_{0}^{1}\mathrm{d}t\,\frac{1}{t}\left[\exp{\left(-(a-t)x\right)}-\exp{\left(-ax\right)}\right]\\ &~~~~~+2\int_{0}^{\infty}\mathrm{d}x\int_{0}^{1}\mathrm{d}t\,\frac{1}{t}\left[\exp{\left(-(a-t)x\right)}\ln{\left(x\right)}-\exp{\left(-ax\right)}\ln{\left(x\right)}\right]\\ &~~~~~+\int_{0}^{\infty}\mathrm{d}x\int_{0}^{1}\mathrm{d}t\int_{0}^{1}\mathrm{d}u\,\frac{1}{tu}\left[e^{-\left(a-t-u\right)x}-e^{-\left(a-t\right)x}-e^{-\left(a-u\right)x}+e^{-ax}\right]\\ &=\frac{\gamma^{2}}{a}\int_{0}^{\infty}\mathrm{d}y\,\exp{\left(-y\right)}+\frac{2\gamma}{a}\int_{0}^{\infty}\mathrm{d}y\,\exp{\left(-y\right)}\left[\ln{\left(y\right)}-\ln{\left(a\right)}\right]\\ &~~~~~+\frac{1}{a}\int_{0}^{\infty}\mathrm{d}y\,\exp{\left(-y\right)}\left[\ln^{2}{\left(a\right)}-2\ln{\left(a\right)}\ln{\left(y\right)}+\ln^{2}{\left(y\right)}\right]\\ &~~~~~+2\gamma\int_{0}^{1}\mathrm{d}t\int_{0}^{\infty}\mathrm{d}x\,\frac{1}{t}\left[\exp{\left(-(a-t)x\right)}-\exp{\left(-ax\right)}\right]\\ &~~~~~+2\int_{0}^{1}\mathrm{d}t\int_{0}^{\infty}\mathrm{d}x\,\frac{1}{t}\left[\exp{\left(-(a-t)x\right)}\ln{\left(x\right)}-\exp{\left(-ax\right)}\ln{\left(x\right)}\right]\\ &~~~~~+\int_{0}^{1}\mathrm{d}t\int_{0}^{1}\mathrm{d}u\int_{0}^{\infty}\mathrm{d}x\,\frac{1}{tu}\left[e^{-\left(a-t-u\right)x}-e^{-\left(a-t\right)x}-e^{-\left(a-u\right)x}+e^{-ax}\right]\\ &=\frac{\gamma^{2}}{a}\int_{0}^{\infty}\mathrm{d}y\,\exp{\left(-y\right)}+\frac{2\gamma}{a}\int_{0}^{\infty}\mathrm{d}y\,\exp{\left(-y\right)}\ln{\left(y\right)}-\frac{2\gamma}{a}\int_{0}^{\infty}\mathrm{d}y\,\exp{\left(-y\right)}\ln{\left(a\right)}\\ &~~~~~\small{+\frac{1}{a}\int_{0}^{\infty}\mathrm{d}y\,\exp{\left(-y\right)}\ln^{2}{\left(a\right)}-2\frac{1}{a}\int_{0}^{\infty}\mathrm{d}y\,\exp{\left(-y\right)}\ln{\left(a\right)}\ln{\left(y\right)}+\frac{1}{a}\int_{0}^{\infty}\mathrm{d}y\,\exp{\left(-y\right)}\ln^{2}{\left(y\right)}}\\ &~~~~~+2\gamma\int_{0}^{1}\mathrm{d}t\,\frac{1}{t}\left[\frac{1}{a-t}-\frac{1}{a}\right]+2\int_{0}^{1}\mathrm{d}t\,\frac{1}{t}\left[-\frac{\gamma+\ln{\left(a-t\right)}}{a-t}+\frac{\gamma+\ln{\left(a\right)}}{a}\right]\\ &~~~~~+\int_{0}^{1}\mathrm{d}t\int_{0}^{1}\mathrm{d}u\,\frac{1}{tu}\left[\frac{1}{a-t-u}-\frac{1}{a-t}-\frac{1}{a-u}+\frac{1}{a}\right].\\ \end{align}$$
Continuing,
$$\begin{align} \mathcal{I}{\left(a\right)} &=\frac{\left[\gamma-\ln{\left(a\right)}\right]^{2}}{a}\int_{0}^{\infty}\mathrm{d}y\,\exp{\left(-y\right)}+\frac{2\left[\gamma-\ln{\left(a\right)}\right]}{a}\int_{0}^{\infty}\mathrm{d}y\,\exp{\left(-y\right)}\ln{\left(y\right)}\\ &~~~~~+\frac{1}{a}\int_{0}^{\infty}\mathrm{d}y\,\exp{\left(-y\right)}\ln^{2}{\left(y\right)}+2\int_{0}^{1}\mathrm{d}t\,\frac{1}{t}\left[-\frac{\ln{\left(a-t\right)}}{a-t}+\frac{\ln{\left(a\right)}}{a}\right]\\ &~~~~~+\int_{0}^{1}\mathrm{d}t\,\frac{1}{t}\int_{0}^{1}\mathrm{d}u\,\left[\frac{1}{\left(a-t\right)\left(a-t-u\right)}-\frac{1}{a\left(a-u\right)}\right]\\ &=\frac{\left[\gamma-\ln{\left(a\right)}\right]^{2}}{a}+\frac{2\left[\gamma-\ln{\left(a\right)}\right]\left(-\gamma\right)}{a}+\frac{\gamma^{2}+\zeta{\left(2\right)}}{a}\\ &~~~~~+2\int_{0}^{1}\mathrm{d}t\,\left[-\frac{\ln{\left(a-t\right)}}{\left(a-t\right)t}+\frac{\ln{\left(a\right)}}{at}\right]\\ &~~~~~+\int_{0}^{1}\mathrm{d}t\,\frac{1}{t}\,\left[\int_{0}^{1}\mathrm{d}u\,\frac{1}{\left(a-t\right)\left(a-t-u\right)}-\int_{0}^{1}\mathrm{d}u\,\frac{1}{a\left(a-u\right)}\right]\\ &=\frac{\zeta{\left(2\right)}+\ln^{2}{\left(a\right)}}{a}\\ &~~~~~-2\int_{0}^{1}\mathrm{d}t\,\left[\frac{\ln{\left(a-t\right)}}{a\left(a-t\right)}+\frac{\ln{\left(a-t\right)}}{at}-\frac{\ln{\left(a\right)}}{at}\right]\\ &~~~~~+\int_{0}^{1}\mathrm{d}t\,\frac{1}{t}\,\left[\frac{\ln{\left(a-t\right)}-\ln{\left(a-t-1\right)}}{\left(a-t\right)}-\frac{\ln{\left(a\right)}-\ln{\left(a-1\right)}}{a}\right]\\ &=\frac{\zeta{\left(2\right)}+\ln^{2}{\left(a\right)}}{a}\\ &~~~~~-\frac{2}{a}\int_{0}^{1}\mathrm{d}t\,\left[\frac{\ln{\left(a-t\right)}}{\left(a-t\right)}+\frac{\ln{\left(a-t\right)}}{t}-\frac{\ln{\left(a\right)}}{t}\right]\\ &~~~~~+\int_{0}^{1}\mathrm{d}t\,\left[\frac{\ln{\left(a-t\right)}-\ln{\left(a-1-t\right)}}{\left(a-t\right)t}-\frac{\ln{\left(a\right)}-\ln{\left(a-1\right)}}{at}\right]\\ &=\frac{\zeta{\left(2\right)}+\ln^{2}{\left(a\right)}}{a}-\frac{1}{a}\int_{0}^{1}\mathrm{d}t\,\left[\frac{2\ln{\left(a-t\right)}}{\left(a-t\right)}+\frac{2\ln{\left(\frac{a-t}{a}\right)}}{t}\right]\\ &~~~~~+\frac{1}{a}\int_{0}^{1}\mathrm{d}t\,\left[\frac{\ln{\left(a-t\right)}}{\left(a-t\right)}-\frac{\ln{\left(a-1-t\right)}}{\left(a-t\right)}+\frac{\ln{\left(\frac{a-t}{a}\right)}}{t}-\frac{\ln{\left(\frac{a-1-t}{a-1}\right)}}{t}\right]\\ &=\frac{\zeta{\left(2\right)}+\ln^{2}{\left(a\right)}}{a}\\ &~~~~~-\frac{1}{a}\int_{0}^{1}\mathrm{d}t\,\frac{\ln{\left(a-t\right)}}{\left(a-t\right)}-\frac{1}{a}\int_{0}^{1}\mathrm{d}t\,\frac{\ln{\left(a-1-t\right)}}{\left(a-t\right)}\\ &~~~~~-\frac{1}{a}\int_{0}^{1}\mathrm{d}t\,\frac{\ln{\left(\frac{a-t}{a}\right)}}{t}-\frac{1}{a}\int_{0}^{1}\mathrm{d}t\,\frac{\ln{\left(\frac{a-1-t}{a-1}\right)}}{t}\\ &=\frac{\zeta{\left(2\right)}+\ln^{2}{\left(a\right)}}{a}\\ &~~~~~-\frac{1}{a}\left[-\frac12\ln^{2}{\left(a-1\right)}+\frac12\ln^{2}{\left(a\right)}\right]\\ &~~~~~-\frac{1}{a}\int_{0}^{1}\mathrm{d}t\,\frac{\ln{\left(a-1-t\right)}}{\left(a-t\right)}\\ &~~~~~-\frac{1}{a}\left[-\operatorname{Li}_{2}{\left(\frac{1}{a}\right)}\right]-\frac{1}{a}\left[-\operatorname{Li}_{2}{\left(\frac{1}{a-1}\right)}\right]\\ &=\frac{\zeta{\left(2\right)}+\ln^{2}{\left(a\right)}}{a}+\frac{\ln^{2}{\left(a-1\right)}-\ln^{2}{\left(a\right)}}{2a}+\frac{1}{a}\operatorname{Li}_{2}{\left(\frac{1}{a}\right)}+\frac{1}{a}\operatorname{Li}_{2}{\left(\frac{1}{a-1}\right)}\\ &~~~~~-\frac{1}{a}\int_{a-1}^{a}\mathrm{d}u\,\frac{\ln{\left(u-1\right)}}{u};~~~\small{\left[t=a-u\right]}\\ &=\frac{\zeta{\left(2\right)}}{a}+\frac{\ln^{2}{\left(a-1\right)}+\ln^{2}{\left(a\right)}}{2a}+\frac{1}{a}\operatorname{Li}_{2}{\left(\frac{1}{a}\right)}+\frac{1}{a}\operatorname{Li}_{2}{\left(\frac{1}{a-1}\right)}\\ &~~~~~-\frac{1}{a}\int_{\frac{1}{a}}^{\frac{1}{a-1}}\mathrm{d}v\,\frac{\ln{\left(\frac{1-v}{v}\right)}}{v};~~~\small{\left[u=\frac{1}{v}\right]}\\ &=\frac{\zeta{\left(2\right)}}{a}+\frac{\ln^{2}{\left(a-1\right)}+\ln^{2}{\left(a\right)}}{2a}+\frac{1}{a}\operatorname{Li}_{2}{\left(\frac{1}{a}\right)}+\frac{1}{a}\operatorname{Li}_{2}{\left(\frac{1}{a-1}\right)}\\ &~~~~~+\frac{1}{a}\int_{\frac{1}{a}}^{\frac{1}{a-1}}\mathrm{d}v\,\frac{(-1)\ln{\left(1-v\right)}}{v}+\frac{1}{a}\int_{\frac{1}{a}}^{\frac{1}{a-1}}\mathrm{d}v\,\frac{\ln{\left(v\right)}}{v}\\ &=\frac{\zeta{\left(2\right)}}{a}+\frac{\ln^{2}{\left(a-1\right)}+\ln^{2}{\left(a\right)}}{2a}+\frac{1}{a}\operatorname{Li}_{2}{\left(\frac{1}{a}\right)}+\frac{1}{a}\operatorname{Li}_{2}{\left(\frac{1}{a-1}\right)}\\ &~~~~~+\frac{1}{a}\left[\operatorname{Li}_{2}{\left(\frac{1}{a-1}\right)}-\operatorname{Li}_{2}{\left(\frac{1}{a}\right)}\right]+\frac{1}{2a}\left[\ln^{2}{\left(\frac{1}{a-1}\right)}-\ln^{2}{\left(\frac{1}{a}\right)}\right]\\ &=\frac{\zeta{\left(2\right)}}{a}+\frac{\ln^{2}{\left(a-1\right)}}{a}+\frac{2\operatorname{Li}_{2}{\left(\frac{1}{a-1}\right)}}{a}.\\ \end{align}$$
Then,
$$\begin{align} \mathcal{I}{\left(2\right)} &=\lim_{a\to2^{+}}\mathcal{I}{\left(a\right)}\\ &=\lim_{a\to2^{+}}\left[\frac{\zeta{\left(2\right)}}{a}+\frac{\ln^{2}{\left(a-1\right)}}{a}+\frac{2\operatorname{Li}_{2}{\left(\frac{1}{a-1}\right)}}{a}\right]\\ &=\frac{\zeta{\left(2\right)}}{2}+\frac{\ln^{2}{\left(1\right)}}{2}+\frac{2\operatorname{Li}_{2}{\left(1\right)}}{2}\\ &=\frac{3}{2}\zeta{\left(2\right)}\\ &=\frac{\pi^{2}}{4}.\blacksquare\\ \end{align}$$