With the substitution $t = zu^{2}$, we have
$$ \mathrm{Ei}(-z) = - \int_{z}^{\infty} \frac{e^{-t}}{t} \, dt = -2 \int_{1}^{\infty} \frac{e^{-zu^{2}}}{u} \, du. $$
Thus for $\Re s > 0$, we have
\begin{align*}
\int_{0}^{\infty} z^{s-1} \mathrm{Ei}(-z)^{2} \, dz
&= 4 \int_{0}^{\infty}\int_{1}^{\infty}\int_{1}^{\infty} \frac{z^{s-1} e^{-z(x^{2}+y^{2})}}{xy} \, dxdydz \\
&= 4 \Gamma(s) \int_{1}^{\infty}\int_{1}^{\infty} \frac{dxdy}{xy(x^{2}+y^{2})^{s}} \\
&= 4 \Gamma(s) \int_{0}^{1}\int_{0}^{1} \frac{x^{2s-1}y^{2s-1}}{(x^{2}+y^{2})^{s}} \, dxdy \tag{1} \\
&= 8 \Gamma(s) \iint\limits_{0 \leq y \leq x \leq 1} \frac{x^{2s-1}y^{2s-1}}{(x^{2}+y^{2})^{s}} \, dxdy, \tag{2}
\end{align*}
where we utilized the substitution $(x, y) \mapsto (x^{-1}, y^{-1})$ at $(1)$ and exploited the symmetry of the domain of integration $\{ x \geq 1, y \geq 1 \}$ along $y = x$ to get $(2)$.
Now with the polar coordinate transform,
\begin{align*}
\int_{0}^{\infty} z^{s-1} \mathrm{Ei}(-z)^{2} \, dz
&= 8 \Gamma(s) \int_{0}^{\frac{\pi}{4}}\int_{0}^{\sec\theta} r^{2s-1}\cos^{2s-1}\theta\sin^{2s-1}\theta \, drd\theta \\
&= \frac{4 \Gamma(s)}{s} \int_{0}^{\frac{\pi}{4}} \frac{\sin^{2s-1}\theta}{\cos\theta}\,d\theta \\
&= \frac{2 \Gamma(s)}{s} \int_{0}^{\frac{\pi}{4}} \frac{\sin^{2s-2}\theta}{1 - \sin^{2}\theta} \cdot 2\cos\theta\sin\theta\,d\theta \\
&= \frac{2 \Gamma(s)}{s} \int_{0}^{\frac{1}{2}} \frac{u^{s-1}}{1 - u} \, du \tag{3} \\
&= \frac{2 \Gamma(s)}{s} \beta\left(\frac{1}{2} ; s, 0 \right)
\end{align*}
where $\beta(x ; a, b)$ is the incomplete beta function and we used the substitution $u = \sin^{2}\theta$ at $(3)$. Now the rest is just a bit of calculus, yielding the desired result.
A slightly general result is as follows:
\begin{align*}
\int_{0}^{\infty} z^{s-1} \Gamma(p, z) \Gamma(q, z) \, dz = \frac{\Gamma(s+p+q)}{s} \left\{ \beta\left( \frac{1}{2}; s+p, q\right) + \beta\left(\frac{1}{2}; s+q, p\right) \right\},
\end{align*}
where $\Gamma(s, z) = \int_{z}^{\infty} t^{s-1}e^{-t} \, dt$ is the incomplete gamma function.
