Evaluate $\int_{0}^{\infty}x\left ( \operatorname{Ci}(x)^2 +\operatorname{si}(x)^2\right )\operatorname{Ci}(x)\text{d}x$ 
Evaluate
$$
I_1=\int_{0}^{\infty}x\left ( \operatorname{Ci}(x)^2
+\operatorname{si}(x)^2\right )\operatorname{Ci}(x)\text{d}x.
$$
Where $\operatorname{Ci}(x)=-\int_{x}^{\infty}\frac{\cos(t)}{t}\text{d}t,\operatorname{si}(x)=-\int_{x}^{\infty}\frac{\sin(t)}{t}\text{d}t$ are cosine integral function and (modified) sine integral function respectively.

By integrating numerically, one gives $I_1=0$. My evaluation is a bit complex, so I would like to see a canonical one. Thanks for replying.

A related version:
$$
I_2=\int_{0}^{\infty}x\left ( \operatorname{Ci}(x)^2
+\operatorname{si}(x)^2\right )\operatorname{si}(x)\text{d}x=\frac\pi2\left(1-2\ln(2)\right).
$$

A sophisticated version:
$$
I_3=\int_{0}^{\infty}\left ( \operatorname{Ci}(x)^2+\operatorname{si}(x)^2\right )^2\left ( 
 \operatorname{Ci}(x)^2+\operatorname{si}(x)^2+3\pi\operatorname{si}(x)\right ) \text{d}x
=6\pi^3\operatorname{Li}_2\left ( \frac{1}{4}  \right ) +12\pi^3\ln(2)^2.
$$
Where $\operatorname{Li}_2(z)=\sum_{k=1}^{\infty}\frac{z^k}{k^2}$ is the dilogarithm.
 A: We can make use of the following result:$$\boxed{\operatorname{Ci}^2(x)+\operatorname{si}^2(x)=\int_0^\infty \frac{e^{-xy}\ln(1+y^2)}{y}dy}$$

$$I_1=\int_0^\infty x\left(\operatorname{Ci}^2(x)+\operatorname{si}^2(x)\right)\operatorname{Ci}(x)dx$$
$$=\int_0^\infty \int_0^\infty \frac{xe^{-xy}\ln(1+y^2)\operatorname{Ci}(x)}{y}dydx$$
$$=\int_0^\infty \frac{\ln(1+y^2)}{y}\left(\int_0^\infty x e^{-xy}\operatorname{Ci}(x)dx\right)dy$$
$$=\int_0^\infty \frac{\ln(1+y^2)}{y}\left(-\frac{d}{dy}\int_0^\infty e^{-xy}\operatorname{Ci}(x)dx\right)dy$$
$$=\frac12\int_0^\infty \frac{\ln(1+y^2)}{y}\frac{d}{dy}\left(\frac{\ln(1+y^2)}{y}\right)dy=\frac14\left(\frac{\ln(1+y^2)}{y}\right)^2\bigg|_0^\infty =0$$
Above the Laplace transform of the cosine integral was used.

Similarly, for $I_2$ we have:
$$I_2=\int_0^\infty \frac{\ln(1+y^2)}{y}\left(\int_0^\infty x e^{-xy}\operatorname{si}(x)dx\right)dy$$
$$=\int_0^\infty \frac{\ln(1+y^2)}{y}\left(-\frac{d}{dy}\int_0^\infty e^{-xy}\operatorname{si}(x)dx\right)dy$$
$$=\int_0^\infty \frac{\ln(1+y^2)}{y}\frac{d}{dy}\left(\frac{\arctan y}{y}\right)dy$$
$$=\int_0^\infty \left(\frac{\ln(1+y^2)}{y^2(1+y^2)}-\color{blue}{\frac{\arctan y\ln(1+y^2)}{y^3}}\right)dy$$
$$\overset{\color{blue}{IBP}}=\int_0^\infty \left(\frac12\frac{\ln(1+y^2)}{y^2}-\frac12\frac{\ln(1+y^2)}{1+y^2}-\color{red}{\frac{\arctan y}{y(1+y^2)}}\right)dy$$
$$\overset{\color{red}{IBP}}=\int_0^\infty \left(\frac12\underbrace{\frac{\ln(1+y^2)}{y^2}}_{\pi}-\underbrace{\frac{\ln(1+y^2)}{1+y^2}}_{\pi\ln 2}+\underbrace{\frac{\ln y}{1+y^2}}_{0}\right)dy=\frac{\pi}{2}-\pi\ln 2$$
This time the Laplace transform of the sine integral was used.

Proof for the mentioned result
Since $2\cos x= e^{ix}+e^{-ix}$, we can write:
$$\operatorname{Ci}(x)=-\int_x^\infty \frac{\cos t}{t}dt=-\frac12\int_x^\infty \frac{e^{it}}{t}dt-\frac12\int_x^\infty \frac{e^{-it}}{t}dt$$
$$=-\frac12\int_x^\infty \int_0^\infty e^{it}e^{-ty}dydt-\frac12\int_x^\infty \int_0^\infty e^{-it}e^{-ty}dydt$$
$$=-\frac12 \int_0^\infty\int_x^\infty e^{-(y-i)t}dtdy-\frac12\int_0^\infty\int_x^\infty e^{-(y+i)t}dtdy$$
$$=-\frac12 e^{ix} \int_0^\infty \frac{e^{-xy}}{y-i}dy-\frac12 e^{-ix} \int_0^\infty \frac{e^{-xy}}{y+i}dy$$
By the same approach, using $2i\sin x= e^{ix}-e^{-ix}$, we arrive at:
$$\operatorname{si}(x)=-\frac1{2i} e^{ix} \int_0^\infty \frac{e^{-xy}}{y-i}dy+\frac1{2i} e^{-ix} \int_0^\infty \frac{e^{-xy}}{y+i}dy$$
Finally, writting $2(a^2+b^2)$ as $(a+b)^2+(a-b)^2$ yields:
$$\small \operatorname{Ci}^2(x)+\operatorname{si}^2(x)=\int_0^\infty \frac{e^{-xy}}{y-i}dy\int_0^\infty \frac{e^{-xz}}{z+i}dz\overset{y+z\to y}=\int_0^\infty \int_z^\infty \frac{e^{-xy}}{(y-z-i)(z+i)}dydz$$
$$\small =\int_0^\infty e^{-xy} \int_0^y \frac{1}{(y-z-i)(z+i)}dzdy=\int_0^\infty \frac{e^{-xy}\ln(1+y^2)}{y}dy$$
