I would like to thank M.N.C.E. for suggesting the use of the identity $$2\log(x)\log(y) = \log^{2}(x) + \log^{2}(y) - \log^{2} \left(\frac{x}{y} \right), $$ where $x$ and $y$ are positive real values.
As I stated here, we have
\begin{align} &\int_{0}^{\infty} [\operatorname{Ei}(-x)]^{4} \, dx \\ &= -4 \int_{0}^{\infty} [\operatorname{Ei}(-x)]^{3} e^{-x} \, dx \tag{1} \\ &= 4 \int_{0}^{\infty} \int_{1}^{\infty} \int_{1}^{\infty} \int_{1}^{\infty} \frac{1}{wyz} e^{-(w+y+z+1)x} \, dw \, dy \, dz \,dx \tag{2} \\& =4 \int_{1}^{\infty} \int_{1}^{\infty} \int_{1}^{\infty} \frac{1}{wyz} \int_{0}^{\infty} e^{-(w+y+z+1)x} \, dx \, dw \, dy \, dz \\ &= 4 \int_{1}^{\infty} \int_{1}^{\infty} \int_{1}^{\infty} \frac{1}{wyz} \frac{1}{w+y+z+1} \, dw \, dy \, dz \\ &= 4 \int_{1}^{\infty} \int_{1}^{\infty} \int_{1}^{\infty} \frac{1}{yz} \frac{1}{y+z+1} \left(\frac{1}{w} - \frac{1}{w+y+z+1} \right) \, dw \, dy \, dz \\ &= 4 \int_{1}^{\infty} \int_{1}^{\infty} \frac{1}{yz} \frac{1}{y+z+1} \log(y+z+2) \, dy \, dz \\ &= 8 \int_{1}^{\infty} \int_{1}^{z} \frac{1}{yz} \frac{1}{y+z+1} \log(y+z+2) \, dy \, dz \tag{3} \\ &= 8 \int_{2}^{\infty} \int_{u-1}^{u^{2}/4} \frac{1}{v} \frac{1}{u+1} \log(u+2) \frac{dv \, du}{\sqrt{u^{2}-4v}} \tag{4} \\& =16 \int_{2}^{\infty} \frac{\log(u+2)}{u+1} \int_{0}^{u-2} \frac{1}{u^{2}-t^{2}} \, dt \, du \tag{5} \\& =16 \int_{2}^{\infty} \frac{\log(u+2)}{u+1} \frac{1}{u} \, \operatorname{artanh} \left( \frac{u-2}{2}\right) \, du \\ &= 8 \int_{2}^{\infty} \frac{\log(u+2) \log(u-1)}{u(u+1)} \, du \\ &= 8 \Bigg(\int_{0}^{1/2} \frac{\log(1-u)\log(1+2u)}{1+u} \, du - \int_{0}^{1/2} \frac{\log(u) \log(1+2u)}{1+u} \, du \\ &- \int_{0}^{1/2} \frac{\log(u) \log (1-u) }{1+u} \, du + \int_{0}^{1/2} \frac{\log^{2}(u)}{1+u} \, du \Bigg). \tag{6} \end{align}
$(1)$ Integrate by parts.
$(2)$ $\operatorname{Ei}(-x) = - \int_{x}^{\infty} \frac{e^{-t}}{t} \, dt = - \int_{1}^{\infty} \frac{e^{-xu}}{u} \, \mathrm du $
$(3)$ The integrand is symmetric about the line $z=y$.
$(4)$ Make the change of variables $u= y+z$, $v=yz$.
$(5)$ Make the substitution $ t^{2}=u^2-4v$.
$(6)$ Replace $u$ with $\frac{1}{u}$.
Then using the identity I mentioned at the beginning of this answer, we have
$$ \begin{align} &\int_{0}^{\infty} [\text{Ei}(-x)]^{4} \, dx \\ &= -4\int_{0}^{1/2} \frac{\log^{2} \left(\frac{1-u}{1+2u} \right)}{1+u} \, du + 4 \int_{0}^{1/2} \frac{\log^{2} \left(\frac{u}{1+2u} \right)}{1+u} \, du + 4 \int_{0}^{1/2} \frac{\log^{2} \left(\frac{u}{1-u} \right)}{1+u} \, du \\& = -12 \int_{1/4}^{1} \frac{\log^{2} (w)}{(1+2w)(2+w)} \, dw +4 \int_{0}^{1/4} \frac{\log^{2} (y)}{(1-2y)(1-y)} \, dy + 4 \int_{0}^{1} \frac{\log^{2}(z)}{(1+2z)(1+z)} \, dz \\&= -6 \int_{1/4}^{1} \frac{\log^{2} (w)}{(1+2w)(1+\frac{w}{2})} \, dw +4 \int_{0}^{1/4} \frac{\log^{2} (y)}{(1-2y)(1-y)} \, dy + 4 \int_{0}^{1} \frac{\log^{2}(z)}{(1+2z)(1+z)} \, dz. \end{align}$$
EDIT:
Partial fraction decomposition followed by two applications of integration by parts shows that $ \begin{align} \int \frac{\log^{2}(x)}{(1+ax)(1+bx)} \, dx = \frac{1}{a-b} \Big(&2\big(\operatorname{Li}_{3}(-bx) - \operatorname{Li}_{3}(-ax) \big) + 2 \log (x) \big( \operatorname{Li}_{2}(-ax) - \operatorname{Li}_{2}(-bx)\big) \\ &+\log^{2}(x) \big(\ln(1+ax) - \log(1+bx) \big) \Big) + C. \end{align}$
This primitive can be used to determine the value of all three definite integrals above.
The final result won't immediately be in the form given by Cleo. Getting it in that form will require the use of polylogarithm identities.
I will open with the main theorem.
Let $s$ be a positive real number. Then the following equality holds.
$$\frac{\pi^2}{4}\int_0^{\infty} \dfrac{\tanh(x)\,\tanh(x s)}{x^2}\,dx = s \int_0^1 \ln\left(\frac{1-x}{1+x}\right) \ln\left(\frac{1-x^s}{1+x^s}\right) \,\frac{dx}{x}.\tag{1}$$
Proof.
First, loagrithmically differentiate the Weierstrass-form product of the hyperbolic cosine,
to obtain $\displaystyle \, \frac{1}{8x}\tanh(x)=\sum_{n=0}^{\infty} \frac1{\pi^2 (2n+1)^2+(2x)^2}\,.$
Also, since (elementarily) we have $\displaystyle \,\,\int_0^{\infty} \frac1{a^2+x^2}\,dx=\frac{\pi}{2a},\,\,$ it follows that
$\displaystyle \,\int_0^{\infty} \frac1{(a^2+x^2)(b^2+x^2)}\,dx=\frac1{b^2-a^2}\int_0^{\infty}\left(\frac1{a^2+x^2}-\frac1{b^2+x^2}\right)dx=\frac{\pi}{2}\,\frac1{ab(a+b)}.$
So,
$$\begin{align*}
\int_0^{\infty} \frac{\tanh(x)\,\tanh(x s)}{64 x^2 s}\,dx\\
&=\int_0^{\infty} \sum_{n,m=0}^{\infty} \dfrac1{(\pi^2(2n+1)^2+(2x)^2)(\pi^2(2m+1)^2+(2xs)^2)}\,\,dx\\
&=\frac1{2^2 (2s)^2} \sum_{n,m=0}^{\infty} \int_0^{\infty} \dfrac1{x^2+\left(\frac{\pi}{2}(2n+1)\right)^2}\,\dfrac1{x^2+\left(\frac{\pi}{2s}(2m+1)\right)^2}\,dx\\
&=\frac1{4\pi^2} \sum_{n,m=0}^{\infty} \dfrac1{(2n+1)\,(2m+1)\,\,((2n+1)+s(2m+1))}\\
&=\frac1{4\pi^2} \sum_{n,m=0}^{\infty} \dfrac1{(2n+1)(2m+1)} \int_0^1 x^{(2n+1)+s(2m+1)}\,\frac{dx}{x}\\
&=\frac1{16\pi^2} \int_0^1 \ln\left(\frac{1-x}{1+x}\right)\ln\left(\frac{1-x^s}{1+x^s}\right)\,\frac{dx}{x}.
\end{align*}$$
Example no. 1
Set $s=1$. With the substitution $x \mapsto \frac{1-x}{1+x},$ we have
$$\begin{align*}
\frac{\pi^2}{4}\int_0^{\infty} \frac{\tanh^2 x}{x^2}\,dx\\
&=\int_0^1 \ln^2\left(\frac{1-x}{1+x}\right)\,\frac{dx}{x}\\
&=2\int_0^1 \frac{\ln^2 x}{1-x^2}\,dx\\
&=2\int_0^1 \ln^2 x \sum_{n=0}^{\infty} x^{2n} \,dx\\
&=\sum_{n=0}^{\infty} \frac{4}{(2n+1)^3}=\frac{7}{2}\zeta(3).
\end{align*}$$
Therefore,
$$ \,\,\int_0^{\infty} \frac{\tanh^2 x}{x^2}\,dx=\frac{14\zeta(3)}{\pi^2}.$$
Example no. 2
Set $s=3$. Employing the same substitution, we have
$$\int_0^1 \ln\left(\frac{1-x}{1+x}\right)\ln\left(\frac{1-x^3}{1+x^3}\right)\,\frac{dx}{x}=2\int_0^1 \frac{\ln x}{1-x^2} \ln\left(\frac{x(x^2+3)}{3x^2+1}\right)dx$$
Now, I was able to obtain the following:
$$ f(a)=\int_0^1 \dfrac{\ln x \,\ln(a^2+x^2)}{1-x^2}dx=-\frac{\pi^2}{8}\ln(1+a^2)+\frac14 F\left(-\frac1{a^2}\right)-2\operatorname{Re} F\left(\frac{i}{a}\right),$$
where
$$ F(z)=\text{Li}_3(z)+2\text{Li}_3(1-z)-\ln(1-z)\text{Li}_2(1-z)-\frac{\pi^2}{6}\ln(1-z)-2\zeta(3).$$
It comes from the fact that whenever $|z|<1$, $\displaystyle \,\, F(z)=\sum_{n=1}^{\infty} \frac{H_{n}^{(2)}}{n}\,z^n.$
Using that notation,
$\displaystyle \,\,2\int_0^1 \frac{\ln x}{1-x^2}\ln\left(\frac{x(x^2+3)}{3x^2+1}\right)dx=\frac{7}{2}\zeta(3)+\frac{\pi^2}{4}\ln3+2f(\sqrt{3})-2f\left(\frac1{\sqrt{3}}\right).$
The trouble was simplifying the hideous, monstrous expression. After several hours of painful simplification by hand, I finally obtained
$$\int_0^1 \ln\left(\frac{1-x}{1+x}\right)\,\ln\left(\frac{1-x^3}{1+x^3}\right)\frac{dx}{x}=\frac{\pi^2}{18}\ln3-\frac{2\pi^2}{3}\ln2+\frac{8}{3}\ln^3 2-\frac{7}{2}\zeta(3)-2\text{Li}_3\left(\frac14\right)\\+16\operatorname{Re} \text{Li}_3(1-i\sqrt{3})-\frac{2\pi}{3}\operatorname{Im}\text{Li}_2\left(-\frac{i}{\sqrt{3}}\right).$$
(can it be simplified more?)
Or,
$$\int_0^{\infty} \frac{\tanh(x)\tanh(3x)}{x^2}\,dx=\frac23\ln3-8\ln2+\frac{32\ln^3 2}{\pi^2}-\frac{42\zeta(3)}{\pi^2}-\frac{24\text{Li}_3\left(\frac14\right)}{\pi^2}\\+\frac{192}{\pi^2}\operatorname{Re}\text{Li}_3(1-i\sqrt{3})-\frac{8}{\pi}\operatorname{Im}\text{Li}_2\left(-\frac{i}{\sqrt{3}}\right).$$
If we set $s=4$ and again substitute $x \mapsto \frac{1-x}{1+x},\,$ we find
$$\int_0^1 \ln\left(\frac{1-x}{1+x}\right)\,\ln\left(\frac{1-x^4}{1+x^4}\right)\frac{dx}{x}=2\int_0^1 \frac{\ln x}{1-x^2} \ln\left(\frac{4x(1+x^2)}{x^4+6x^2+1}\right)\,\frac{dx}{x}
\\=-\frac{\pi^2}{2}\ln2+\frac{7}{2}\zeta(3)+2f(1)-2f(\sqrt{2}+1)-2f(\sqrt{2}-1)$$
I am not courageous enough to even begin simplifying that.
It seems that a closed form may exist for each natural $s$ (and therfore, for each $1/n$ where $n$ is natural).
For example, under the subsitution $x \mapsto \frac{1-x}{1+x}$, the expression $\displaystyle \frac{1-x^5}{1+x^5}$ turns into $\displaystyle \frac{x(x^4+10x^2+5)}{5x^4+10x^2+1}$.
Factorizing $x^4+10x^2+5=(x^2+5+2\sqrt{5})(x^2+5-2\sqrt{5})$ and $\displaystyle 5x^4+10x^2+1=5\left(x^2+\frac1{5+2\sqrt{5}}\right)\left(x^2+\frac1{5-2\sqrt{5}}\right)$,
it follows that
$$ \int_0^1 \ln\left(\frac{1-x}{1+x}\right)\ln\left(\frac{1-x^5}{1+x^5}\right)\frac{dx}{x}=\frac72\zeta(3)+\frac{\pi^2}{4}\ln5+2f\left(\sqrt{5+2\sqrt{5}}\right)+2f\left(\sqrt{5-2\sqrt{5}}\right)-2f\left(\frac1{\sqrt{5+2\sqrt{5}}}\right)-2f\left(\frac1{\sqrt{5-2\sqrt{5}}}\right).$$
Similarly,
$$\int_0^1 \ln\left(\frac{1-x}{1+x}\right)\ln\left(\frac{1-x^6}{1+x^6}\right)\frac{dx}{x}=\frac72\zeta(3)-\frac{\pi^2}{4}\ln6+2f(\sqrt{3})+2f\left(\frac1{\sqrt{3}}\right)-2f(1)-2f(2+\sqrt{3})-2f(2-\sqrt{3}).$$
A pattern can be seen.
It seems that we can always factorize the expression that results by applying $x\mapsto \frac{1-x}{1+x}$ to $\frac{1-x^n}{1+x^n}$ into factors of the form $x$ and $x^2+r^2$ ($r \in \mathbb{C}$), implying that there is indeed a closed form in terms of logarithms and polylogarithms
for $\int_0^{\infty} \tanh(x)\tanh(xn)/x^2\,\,dx$.
In fact, the roots $r$ in the factorization seem to always be $\tan(q \pi)$ (up to multiplication by an imaginary unit), where q is a rational number (whose denominator is $n$ when $n$ is odd and $2n$ when it is even) . For example, $\displaystyle \sqrt{2}-1=\tan\left(\frac{\pi}{8}\right),\,\,\,\,\sqrt{5-2\sqrt{5}}=\tan\left(\frac{\pi}{5}\right),\,\,\,\,\,2-\sqrt{3}=\tan\left(\frac{\pi}{12}\right)\,\,$etc.
I am too lazy to write a general formula that works for every natural number.
$\Large \mathbf {EDIT}$
For the sake of completeness, I add below the general formula, which works for every natural number.
Theorem 2.$\,$
Let $n$ be a positive integer. Define the function $F$ by
$$ F(z)=\text{Li}_3(z) + 2\text{Li}_3(1-z) - \ln(1-z)\text{Li}_2(1-z) -\frac{\pi^2}{6}\ln(1-z)-2\zeta(3),$$
where we assume the principal value of the logarithm.
Now, define the function $g$ by
$$g(z)=F(-z^2)-8 \operatorname{Re} F(i z).$$
Then
$$
\int_0^1 \ln\left(\frac{1-x}{1+x}\right)\ln\left(\frac{1-x^n}{1+x^n}\right)\frac{dx}{x}
=\frac{7}{2}\zeta(3)+\frac12 \sum_{k=1}^{n-1} (-1)^k g\left(\cot\left(\frac{k \pi}{2 n}\right)\right) .$$
This result, with the help of theorem $(1)$, can easily be seen to establish a closed form
for $\mathcal A(n) = \int_0^{\infty} \frac{\tanh(x) \tanh(n x)}{x^2}dx $ for each natural $n$.
Best Answer
$$\begin{align}\int_0^\infty\operatorname{Ei}^3(-x)\,dx&=-3\operatorname{Li}_2\left(\frac14\right)-6\ln^22.\\\\\int_0^\infty\operatorname{Ei}^4(-x)\,dx&=24\operatorname{Li}_3\left(\frac14\right)-48\operatorname{Li}_2\left(\frac13\right)\ln2-13\,\zeta(3)\\&-32\ln^32+48\ln^22\,\ln3-24\ln2\,\ln^23+6\pi^2\ln2.\end{align}$$