A Neat Identity Involving Zeta Zeroes

The result found by OP turns out to be quite generic; it holds for a wide range of sequences $\rho_n$ and not just the zeros of the $\zeta$-function. A precise formulation of what is observed is the following:

$$\lim_{N\to\infty} \frac{1}{N}\sum_{n=1}^N \text{Ei}(\rho_n) = \pi i \tag{1}$$

The sum above is the Cesaro mean of the sequence $a_n = \text{Ei}(\rho_n)$. If a sequence $a_n$ converges to $a$ then the Cesaro mean of $a_n$ also converges to $a$ (see e.g. this question) so $(1)$ would follow if we could show that $\lim_{n\to\infty}\text{Ei}(\rho_n) = \pi i$.

This is indeed true and the reason for this is as follows:

1. We can write $\text{Ei}(z)$ on the form $\text{Ei}(z) = -\Gamma(0,z) + \pi i$ where the incomplete $\Gamma$-function has the asymptotical expansion $\Gamma(0, z) = \frac{e^{-z}}{z} + \mathcal{O}(z^{-2})$. This implies that $\Gamma(0, z_n) \to 0$ for any sequence $z_n = x_n + i y_n$ where $y_n\to \infty$ and where $x_n$ is bounded.
2. The non-trivial zeros of the $\zeta$-function satisfy the conditions on the sequence above since they all lie in the strip $0<\Re[z]< 1$ in the complex plane and $\lim_{n\to\infty}\Im[\rho_n] = \infty$. This means that $\lim_{n\to\infty}\Gamma(0,\rho_n) = 0$ so $\lim_{n\to\infty} \text{Ei}(\rho_n) = \pi i$ and $(1)$ follows.

Note that the conditions in point 1. above are far from being very restrictive so there are infinitely many sequences for which $(1)$ holds (random examples are $\rho_n = 7 + n^2i$ and $\rho_n = \frac{2n}{1+n} + n\log(n)i$).

This is not a proof, just some tips that can may be leads to a proof

(I really hope someone manage to find a proof for this result)

You can consider the other form for the exponential integral:

$$Ei(z)=\gamma +\ln(x)+\sum_{k=1}^\infty\frac {x^k}{k\cdot k!},$$

where $\gamma$ is the Euler-Mascheroni's constant.

Call the $n$-th non trivial zero

$$p_n=:\frac 12 +i \alpha_n.$$

So you get

$$Ei(p_n)=\gamma +\ln\left(\frac 12 +i \alpha_n\right)+\sum_{k=1}^\infty\frac {\left(\frac 12 +i \alpha_n\right)^k}{k\cdot k!}.$$

Use the principal determination of the complex logarithm: $\ln(x+iy)=\ln(x)+i\arg(x+iy)$, where

$$\arg(x+iy)=\arctan\left(\frac yx\right).$$

So

$$Ei(p_n)=\gamma +\ln\left(\frac 12\right) +i\arctan(2\alpha_n)+\sum_{k=1}^\infty\frac {\left(\frac 12 +i \alpha_n\right)^k}{k\cdot k!}.$$

So

$$\mathfrak I(Ei)(p_n)=\arctan(2\alpha_n)+\sum_{k=1}^\infty\left(\sum_{\substack{m=0 \\ m\equiv 1 (4)}}^{k}\frac {\binom km\alpha_n^m}{2^{k-m}k\cdot k!}-\sum_{\substack{m=0 \\ m\equiv 3 (4)}}^{k}\frac {\binom km\alpha_n^m}{2^{k-m}k\cdot k!}\right),$$

using

$$\binom km=\frac{k!}{m!(k-m)!}$$

you get

$$\mathfrak I(Ei)(p_n)=\arctan(2\alpha_n)+\sum_{k=1}^\infty\left(\sum_{\substack{m=0 \\ m\equiv 1 (4)}}^{k}\frac {\alpha_n^m}{2^{k-m}k\cdot m!(k-m)!}-\sum_{\substack{m=0 \\ m\equiv 3 (4)}}^{k}\frac {\alpha_n^m}{2^{k-m}k\cdot m!(k-m)!}\right).$$

The key is to use now the approximation:

$$\alpha_n\approx \frac{2\pi n}{\ln(n)}$$

to put $\pi$ in the place.

To conclude, you have to prove that

$$\lim_{x\to \infty}\frac 1x\sum_{n=1}^x\mathfrak I(Ei)(p_n)=\pi.$$

For the first term, you can say that

$$\frac 1x\sum_{n=1}^x \arctan(2\alpha_n)=\frac 1x\left(\text{a finite number of terms}\right)+\frac 1x(x-K)\left(\text{something close to $\frac\pi2$}\right),$$

so

$$\lim_{x\to \infty}\frac 1x\sum_{n=1}^x \arctan(2\alpha_n)=\frac\pi 2.$$

Now, you need to deal with the second term, you can use Stirling formula

$$y!\sim \sqrt{2\pi y}\frac{y^y}{e^y}$$

to get

$$\frac {\alpha_n^m}{2^{k-m}k\cdot m!(k-m)!}=\frac{(2\pi n)^m e^m e^{k-m}}{\ln(n)^m 2^{k-m}k\sqrt{2\pi m}m^m \sqrt{2\pi(k-m)}(k-m)^{k-m}}$$

$$=k\left(\frac{2\pi n}{m\ln(n)}\right)^m\left(\frac e{2(k-m)}\right)^{k-m}\sqrt{2\pi k}.$$