I arrived at the following result
$$\tag{1}\int^\infty_0 z^{p-1} E^2(z)\,dz=\frac{\Gamma(p)}{p}\int^1_0\frac{_2F_1(p,p;p+1;-\frac{1}{z})}{z}\,dz$$
where the exponential integral $E(z)$ is defined as
$$E(z)=\int^\infty_z \frac{e^{-t}}{t}\,dt$$
I have two questions
[1] Does (1) hold for all $p>0$ ?
[2] Is there a way to simplify or solve the integral on the right ?
Let \begin{align} E_{1}(x) = \int_{x}^{\infty} \frac{e^{-t}}{t} \, dt \end{align} then the following two identities can be seen to be \begin{align} \int_{0}^{\infty} x^{n} E_{1}(ax) E_{1}(bx) \, dx &= - \frac{n!}{n+1} \left[ \frac{1}{a^{n+1}} \left\{ \ln\left(\frac{b}{a+b}\right) + \sum_{m=1}^{n} \frac{1}{m} \left( \frac{a}{a+b} \right)^{m} \right\} \right. \\ & \hspace{20mm} \left. + \frac{1}{b^{n+1}} \left\{ \ln\left(\frac{a}{a+b}\right) + \sum_{m=1}^{n} \frac{1}{m} \left( \frac{b}{a+b} \right)^{m} \right\} \right]. \end{align} When $a=b=1$ this becomes \begin{align} \int_{0}^{\infty} x^{p} E_{1}^{2}(x) \, dx &= \frac{\Gamma(p+1)}{2^{p} \, (p+1)} \, \sum_{m=0}^{\infty} \frac{1}{2^{m} (m+p+1)} . \end{align}