Put $x = e^t$ ($t \ge 0$). Because $\operatorname{li}(e^t) = \operatorname{Ei}(t)$, the original integral becomes
$$
I \;=\; \int_{0}^{\infty} (e^{t}-1)\,e^{-3t}\,
\operatorname{Ei}^2(t)\,\text{d}t = \int_{0}^{\infty} e^{-2t}\,
\operatorname{Ei}^2(t)\,\text{d}t - \int_{0}^{\infty} e^{-3t}\,
\operatorname{Ei}^2(t)\,\text{d}t.
$$
Therefore, we must evaluate the integral
$$
I(p) = \int_0^\infty e^{-pt} \operatorname{Ei}^2(t) \, dt, \quad \text{where } p \in \mathbb{N}.
$$
The exponential integral satisfies $\operatorname{Ei}(t) \sim \frac{e^t}{t}$ as $t \to \infty$, so the integrand behaves like
$$
e^{-pt} \cdot \left( \frac{e^t}{t} \right)^2 = \frac{e^{-(p-2)t}}{t^2}.
$$
Therefore, the integral converges if and only if $p \geq 2$.
We integrate by parts:
- Let $u = \operatorname{Ei}^2(t)$, $dv = e^{-pt}\,dt$,
- Then $du = \frac{2\,\operatorname{Ei}(t)\,e^t}{t}\,dt$, and $v = -\frac{e^{-pt}}{p}$.
This gives:
$$
I(p) = \frac{2}{p} \int_0^\infty \frac{e^{-(p-1)t}}{t} \operatorname{Ei}(t)\,dt. \tag{1}
$$
For $t > 0$, we use the representation:
$$
\operatorname{Ei}(t) = \gamma + \log(t) + \int_0^1 \frac{e^{tu} - 1}{u}\,du,
$$
where $\gamma$ is the Euler–Mascheroni constant.
Substitute this into $(1)$ and split the integral:
$$
I(p) = \frac{2}{p} \left( J_1 + J_2\right),
$$
where:
- $J_1 = \gamma \int_0^\infty \left(\frac{e^{-(p-1)t}}{t} + e^{-(p-1)t} \frac{\log(t)}{t}\right) \, dt = -\frac{\gamma}{p-1} \log(p-1) + \frac{1}{2(p-1)} \left[ \log^2(p-1) + \zeta(2) \right],$
- $J_2 = \int_0^1 \frac{du}{u} \int_0^\infty \frac{e^{-(p-1-u)t} - e^{-(p-1)t}}{t} \, dt$.
Use the known result:
$$
\int_0^\infty \frac{e^{-a t} - e^{-b t}}{t} \, dt = \log\left(\frac{b}{a}\right), \quad a, b > 0,
$$
to compute the inner integral:
$$
J_2 = \int_0^1 \frac{1}{u} \log\left( \frac{p - 1}{p - 1 - u} \right) du = \operatorname{Li}_2\left( \frac{1}{p - 1} \right).
$$
Putting everything together:
$$
I(p) = \frac{2}{p} \left( -\frac{\gamma}{p-1} \log(p-1) + \frac{1}{2(p-1)} \left[ \log^2(p-1) + \zeta(2) \right] + \operatorname{Li}_2\left( \frac{1}{p - 1} \right) \right).
$$
After simplifying:
$$
I(p) = \frac{1}{p} \left[ \zeta(2) + \log^2(p-1) + 2 \operatorname{Li}_2\left( \frac{1}{p-1} \right) \right], \quad p \ge 2.
$$
Then, the integral of the question is
$$I = \frac{1}{2}[\zeta(2) + \log^2(1) + 2\text{Li}_2(1)] - \frac{1}{3}[\zeta(2) + \log^2(2) + 2\text{Li}_2(1/2)] = \frac{5}{6}\zeta(2),$$
where we used that
$$2\text{Li}_2(1/2) = \zeta(2) - \log^2(2).$$
