CORRECTED:
We may get the limit for $\operatorname{Ei}(\sigma+it)$ as $t\to \infty$ by using following asymptotic series :
$$\tag{1} \operatorname{Ei}(z)\sim \operatorname{sgn}(\Im(z))i\pi+\frac {e^z}z\left(1+\frac {1!}{z}+\frac {2!}{z^2}+\frac {3!}{z^3}+\cdots\right)$$
(the terms at the right may be obtained by repetitive integration by parts of $\frac {e^z}{z^n}$ getting each time a greater power $n$ at the denominator, for more about this see $(12)$ in Pegoraro and Slusallek's paper of 2011 'On the Evaluation of the Complex-Valued
Exponential Integral' where the Mathematica conventions are used)
If we apply this to $z=\sigma+it$ as $t\to \pm\infty$ we get
$$\tag{2} \operatorname{Ei}(\sigma+it)\sim \operatorname{sgn}(t)i\pi+\frac {e^{\sigma+it}}{\sigma+it}$$
Since $\left|\frac {e^{\sigma+it}}{\sigma+it} \right|\sim \frac {e^{\sigma}}{|t|}$ as $t\to \pm\infty$ we get $\ \operatorname{Ei}(\sigma+it)\sim \operatorname{sgn}(t)i\pi\ $ and your result.
But you wanted a proof without series expansion at $\infty$ so let's try to use directly the integral definition of $\operatorname{Ei}$ and 'suppose' that :
$$\tag{3} \overline{\operatorname{Ei}}(s):=\int_{-\infty}^s \frac {e^z}z\,dz$$
Let's study $\overline{\operatorname{Ei}}(it)$ as $t \to +\infty$ first (case $\sigma=0$) and use following contour at the top right of $0$ supposing that $R=t$ (for $t\to -\infty$ simply consider the symmetric contour!) :
$\overline{\operatorname{Ei}}$ is of course analytic inside this contour so that any 'internal' path from $(-R)$ to $(it)$ should give the same integral.
Using the direct upper path we get :
$$\int_{-R}^{iR} \frac {e^z}z\,dz = \int_{\pi}^{\frac {\pi}2} \frac {e^{R e^{i\phi}}}{Re^{i\phi}}\,iRe^{i\phi}\,d\phi=i\int_{\pi}^{\frac {\pi}2} e^{R e^{i\phi}}\,d\phi$$
Since $\cos(\phi)<0$ for $\phi \in (\frac {\pi}2,\pi)\ $ and $\ \left|\int_{-R}^{iR} \frac {e^z}z\,dz\right| \le \int_{\frac {\pi}2}^{\pi} e^{R \cos(\phi)}\,d\phi\ $ we conclude that the integral from $(-R)$ to $(iR)$ admits the limit $0$ as $R\to \infty$ (with an error of order $\frac 1R$).
Of course the other path should give the same result and :
$$\overline{\operatorname{Ei}}(it)-\overline{\operatorname{Ei}}(-R)=\int_{-R}^{-r} \frac {e^z}z\,dz + i\int_{\pi}^{\frac {\pi}2} e^{r e^{i\phi}}\,d\phi+\int_r^t \frac {e^{iu}}u\,du$$
(allowing us to compare the exponential integral of real and pure imaginary values but that's another subject...)
Find the limit of $\overline{\operatorname{Ei}}(\sigma+it)$ as $t\to \infty$ is not more difficult since :
$$\overline{\operatorname{Ei}}(\sigma+it)-\overline{\operatorname{Ei}}(it)=\int_0^{\sigma} \frac {e^{u+it}}{u+it}\,du$$
So that, for $\sigma > 0$, the absolute value of the difference will be bounded by $\frac{\sigma e^{\sigma}}{|t|}$ and go to $0$ as $t\to \infty$ (consider a rectangle with a larger value of $t$ say $T$ that we bring to $\infty$ would be better but this should be enough to get the idea).
This is all nice and well and seems to prove that $\overline{\operatorname{Ei}}(\sigma+it)$ goes to $0$ as $t\to \infty$ but where is the $\pi$ if not in the sky ?? :-)
The $\pi$ is merely an artifact of the cut chosen for $\operatorname{Ei}$ (on the negative real axis) instead of the cut on the positive real axis implied by the integral expression I used in $\overline{\operatorname{Ei}}$ (that's the reason of the over-line everywhere). To go from the cut on the positive axis ($\overline{\operatorname{Ei}}$) to the cut at the negative axis ($\operatorname{Ei}$) is not really difficult : add $\pi i$ to the result if the imaginary part of the parameter is positive and subtract $\pi i$ if the imaginary part is negative (I updated my answer to this older thread with graphics to make things clearer ; $\overline{\operatorname{Ei}}(z)$ corresponds there to the definition given by Lebedev and is equal to $-\operatorname{E_1}(-z)$ ).
This study showed that the limit of $\overline{\operatorname{Ei}}(z)$ was $0$ when $z$ went to $\infty$ in the left half-plane or along vertical lines in the right half-plane (it should go to $\infty$ in the remaining cases).
Concerning your second question I still don't know what to think about it (or even if this specific approach is useful at all as $t\to \infty$). Isn't $\operatorname{li}\left(x^{1-s}\right)$ decreasing while the error term is increasing ?
By the way let's observe that a more precise bound was given by Schoenfeld if you suppose the R.H. and $x\ge 2657$ : $\displaystyle |\pi(x)-\operatorname{li}(x)| < \frac{\sqrt{x}\ln(x)}{8\pi}$
I tried some numerical experiments and found that for $x$ fixed and growing values of $t$ $P_x(t)$ was oscillating 'randomly' without noticeable modification. In fact it seemed scale invariant relatively to $t$ and this for small or large values of $\sigma$ (so that it doesn't behave as $\operatorname{li}$ nor as the error term!).
It could perhaps be interesting to know the algorithm used by Wolfram Alpha to compute the 'prime zeta function' ($P_x(s)$ as $x \to \infty$) in the complex plane since they could evaluate a finite sum combined with an extrapolation (see Neat examples for the case $\sigma =\frac 12$). For real parts of order $0.4$ or more if not on the imaginary line I got this plot of PrimeZetaP[0.4+t*i] :
