We can actually use $I_{1}$ and $I_{2}$ to find closed-form expressions for $$\int_{0}^{\infty} x \operatorname{Ci}(x)^{3} \, \mathrm dx$$ and
$$\int_{0}^{\infty} x \operatorname{si}(x)^{3} \, \mathrm dx.$$
Let $E_{1}(z)$ be the exponential integral function defined in the right half-plane by the integral $$E_{1}(z) = \int_{1}^{\infty} \frac{e^{-zt}}{t} \, \mathrm dt, \quad \left( \Re(z) \ge 0, z \ne 0 \right). $$
For positive values of $x$, $E_{1}(x) = - \operatorname{E}_{i}(-x)$.
Similar to the my answer here, if we integrate the function $zE_{1}(z)^{3}$ around an infinitely-large closed quarter circle in the first quadrant of the complex plane, we get $$ \begin{align} \oint zE_{1}(z)^{3} dz &=\int_{0}^{\infty} x E_{1}(x)^{3} \, \mathrm dx - \int_{0}^{\infty} (ix) E_{1}(ix)^{3} \, i \, \mathrm dx \\ &= \small\int_{0}^{\infty} x E_{1}(x)^{3} \, \mathrm dx + \int_{0}^{\infty}x\left(-\operatorname{Ci}(x)^{3}+3i \operatorname{Ci}(x)^{2}\operatorname{si}(x) +3 \operatorname{Ci}(x) \operatorname{si}(x)^{2} -i \operatorname{si}(x)^{3}\right) \, \mathrm dx \\ &=0. \end{align}$$
Equating the real parts on both side of the equation, we have $$\int_{0}^{\infty} x E_{1}(x)^{3} \, \mathrm dx+ \int_{0}^{\infty} x\left(-\operatorname{Ci}(x)^{3}+3\operatorname{Ci}(x) \operatorname{si}(x)^{2} \right) \, \mathrm dx =0.$$
And combining the above equation with $I_{1}$, we get $$\int_{0}^{\infty}x \operatorname{Ci}(x)^{3} \, \mathrm dx = \frac{1}{4} \int_{0}^{\infty}x E_{1}(x)^{3} \, \mathrm dx. $$
(If you ask WolframAlpha to approximate the value of $\int_{0}^{a} x \left(\operatorname{Ci}(x)^{3}+ \frac{1}{4} \operatorname{E}_{i}(-x)^{3} \right) \, \mathrm dx $ for increasing positive values of $a$, it returns results that appear to be approaching zero. But if you ask it to approximate $\int_{0}^{\infty} x \left(\operatorname{Ci}(x)^{3}+ \frac{1}{4} \operatorname{E}_{i}(-x)^{3} \right) \, \mathrm dx $, it strangely returns an approximation that is nowhere near zero. This could simply be the result of WolframAlpha needing more computation time.)
To evaluate $$\int_{0}^{\infty}xE_{1}(x)^{3}\, \mathrm dx, $$ we can do what I did here to evaluate $\int_{0}^{\infty}E_{1}(x)^{3}\, \mathrm dx $.
Specifically, we have
$\begin{align} \int_{0}^{\infty}xE_{1}(x)^{3}\, \mathrm dx &\overset{(1)}{=} \frac{3}{2} \int_{0}^{\infty} x e^{-x}E_{1}(x)^{2} \, \mathrm dx \\ &= \frac{3}{2} \int_{0}^{\infty}x e^{-x} \int_{1}^{\infty} \int_{1}^{\infty} \frac{e^{-xy}}{y} \frac{e^{-xz}}{z} \, \mathrm dy \, \mathrm dz \, \mathrm dx \\ &= \frac{3}{2} \int_{1}^{\infty} \int_{1}^{\infty} \frac{1}{yz} \int_{0}^{\infty} x e^{-(1+y+z)x} \, \mathrm dx \, \mathrm dy \, \mathrm dz \\ &= \frac{3}{2} \int_{1}^{\infty} \int_{1}^{\infty} \frac{1}{yz} \frac{1}{(1+y+z)^{2}} \, \mathrm dy \, \mathrm dz \\ &= \frac{3}{2} \int_{1}^{\infty} \frac{1}{z} \int_{1}^{\infty}\left(\frac{1}{(1+z)^2} \left(\frac{1}{y}-\frac{1}{1+y+z} \right) -\frac{1}{1+z} \frac{1}{(1+y+z)^{2}}\right) \, \mathrm dy \, \mathrm dz \\ &= \frac{3}{2} \int_{1}^{\infty} \frac{1}{z} \left(\frac{\ln(2+z)}{(1+z)^{2}} -\frac{1}{(1+z)(z+2)} \right) \, \mathrm dz \\ &= \frac{3}{2} \left( \int_{1}^{\infty} \frac{1}{z} \frac{\ln(2+z)}{(1+z)^{2}} \, \mathrm dz - \ln(2) + \frac{\ln(3)}{2} \right) \\ &\overset{(2)}{=} \frac{3}{2} \left(\int_{0}^{1} \frac{u \ln(1+2u)}{(1+u)^{2}} \, \mathrm du - \int_{0}^{1} \frac{u \ln(u)}{(1+u)^{2}} \, \mathrm du - \ln(2) + \frac{\ln(3)}{2} \right) \\ &=\frac{3}{2} \left(\int_{0}^{1} \frac{u \ln(1+2u)}{(1+u)^{2}} \, \mathrm du + \frac{\pi^{2}}{12} -2 \ln (2) + \frac{\ln(3)}{2} \right) \\ &\overset{(3)}{=} \small \frac{3}{2} \left(\frac{\ln(3)}{2} + \ln(2) \ln(3) -2 \int_{0}^{1} \frac{\mathrm du}{(1+2u)(1+u)} -2 \int_{0}^{1} \frac{\ln(1+u)}{1+2u} \, \mathrm du + \frac{\pi^{2}}{12} -2 \ln (2) + \frac{\ln(3)}{2} \right) \\ &=\frac{3}{2} \left(\ln(2) \ln(3) -\ln(3) -2 \int_{0}^{1} \frac{\ln(1+u)}{1+2u} \, \mathrm du + \frac{\pi^{2}}{12} \right) \end{align}$
where $$ \begin{align} \int_{0}^{1} \frac{\ln(1+u)}{1+2u} \, \mathrm du &= \frac{1}{2} \int_{1}^{3} \frac{\ln (1+w) - \ln(2)}{w} \, \mathrm dw \\ &= \frac{1}{2} \left( \operatorname{Li}_{2}(-1)- \operatorname{Li}_{2}(-3) -\ln(2) \ln(3) \right) \\ & \overset{(4)}{=} \frac{\pi^{2}}{24} + \frac{\ln^{2}(3)}{4}+\frac{1}{2} \, \operatorname{Li}_{2}\left(-\frac{1}{3} \right) - \frac{\ln(2) \ln(3)}{2}. \end{align}$$
Therefore,
$$ \begin{align} \int_{0}^{\infty}xE_{1}(x)^{3}\, \mathrm dx &= -\int_{0}^{\infty} x \operatorname{E}_{i}(-x)^{3} \mathrm dx \\ &= \frac{3}{2} \left(2 \ln(2) \ln(3) -\ln(3) -\frac{\ln^{2}(3)}{2} -\operatorname{Li}_{2} \left(-\frac{1}{3} \right) \right) \\ &= 0.1949195673... \end{align}$$
and $$\int_{0}^{\infty} x \operatorname{Ci}(x)^{3} \, \mathrm dx = \frac{3}{8} \left(2 \ln(2) \ln(3) -\ln(3) -\frac{\ln^{2}(3)}{2} -\operatorname{Li}_{2} \left(-\frac{1}{3} \right) \right). $$
$(1)$ Integrate by parts.
$(2)$ Let $z= \frac{1}{u}$.
$(3)$ Integrate by parts.
$(4)$ Use the dilogarithm inversion formula.
If we equate the imaginary parts on both sides of the equation $$\small\int_{0}^{\infty} x E_{1}(x)^{3} \, \mathrm dx + \int_{0}^{\infty}x\left(-\operatorname{Ci}(x)^{3}+3i \operatorname{Ci}(x)^{2}\operatorname{si}(x) +3 \operatorname{Ci}(x) \operatorname{si}(x)^{2} -i \operatorname{si}(x)^{3}\right) \, \mathrm dx =0, $$ we have $$\int_{0}^{\infty} \left(3x\operatorname{Ci}(x)^{2}\operatorname{si}(x) - x \operatorname{si}(x)^{3} \right) \, \mathrm dx =0. $$
And combining the above equation with $I_{2}$, we get $$\int_{0}^{\infty} x \operatorname{si}(x)^{3} \, \mathrm dx = \frac{3 \pi}{8} \left(1-2 \ln(2) \right).$$