In this question., a very inefficient, yet rigorous analytic approach for finding the next prime was established. I wondered whether a similar approach could exist to find the next non-trivial zero ($\rho_k$) of the the Riemann $\zeta$-function.
The series:
$$\sigma_r = \sum_{k=1}^{\infty} \left( \frac{1}{\rho_k^r} + \frac{1}{(1-\rho_k)^r}\right) \quad r \in \mathbb{N}$$
has a couple of related closed forms, see f.i. Lehmer's groundwork on this, or using the coefficients of the Taylor expansion of the Riemann $\xi$-function or applying the Stieltjes constants $\gamma_x$. EDIT: e.g.:
$\sigma_1= 1 + \frac{\gamma}{2}- \frac{\ln(4\pi)}{2}$
$\sigma_2= 1 + \gamma^2- \frac34\zeta(2)+2\gamma_1$
$\sigma_3= 1 + \gamma^3- \frac78\zeta(3)+3\gamma\gamma_1 + \frac32\gamma_2$
$...$ $\displaystyle\sigma_r=1+\left(\frac{1}{2^r}-1 \right )\zeta(r)+\frac{\gamma\,\gamma_{r-2}}{\Gamma(r-1)}+\frac{r\gamma_{r-1}}{\Gamma(r)}-\sum_{j=1}^{r-2}\frac{\gamma_{j-1}}{\Gamma(j)}\left( 1+\left(\frac{1}{2^{r-j}}-1\right)\zeta(r-j)-\sigma_{r-j}\right)$ for $r>1$, see here.
Now define the following function to 'recover' a $\rho$ from these closed forms:
$$f(r,N,x)= \sigma_r -\left( \sum_{k=1}^{N} \left( \frac{1}{\rho_k^r} + \frac{1}{(1-\rho_k)^r}\right)+ \left( \frac{1}{x^r} + \frac{1}{(1-x)^r}\right)\right)$$
where $r, N \in \mathbb{N}$ and $x$ is the unknown next non-trivial zero ($\rho_{N+1}$).
In Maple:
for N from 0 to 6 do N, fsolve(f(33, N, x), x = 0 + 12*I .. 1 + 42*I, complex) end do;
yields this encouraging list for $\rho_{N+1}$ for $r=33$:
0, 0.50000000000000000000 + 14.134724467544674288 I
1, 0.50000000000000000000 + 21.020287719482273773 I
2, 0.50000000000000000000 + 25.009821593878642551 I
3, 0.50000000000000000000 + 30.365139393045860165 I
4, 0.50000000000000000000 + 32.923423713014869349 I
5, 0.50000000000000000000 + 37.517148219902925704 I
6, 0.50000000000000000000 + 40.750950425647022339 I
Accuracy improves for higher $r$, let's test $N=2$ and derive $\rho_3$ for increasing $r$ :
for r from 5 by 5 to 35 do r, fsolve(f(r, 2, x), x = 0 + 12*I .. 1 + 42*I, complex) end do;
5, 0.50000000000000000000 + 22.755757318210846210 I
10, 0.50000000000000000000 + 24.477994717471646165 I
15, 0.50000000000000000000 + 24.922438734940328122 I
20, 0.50000000000000000000 + 24.979676684194649394 I
25, 0.50000000000000000000 + 25.004076475496886667 I
30, 0.50000000000000000000 + 25.008087373823152633 I
35, 0.50000000000000000000 + 25.010200234483713493 I
Computations quickly require higher precision since $\sigma_r$ becomes very small.
Now, contrary to the primes that are integers and have a minimal distance of $2$ between the odd primes, we are now dealing with probably irrational numbers and no known minimal distance between them (nothing stops Lehmer-pairs from becoming infinitely small).
Q: Could the above approach of analytically 'recovering' the next zero be made rigorous up to a fixed $n$-digits accuracy?
Added observation:
Applying the same approach to the relation between the Keiper-Li constants and their series expression with (powered) $\rho$s :
$$\lambda_n = \sum_{k=1}^{\infty} \left(\left(1-\left(1- \frac{1}{\rho_k}\right)^n\right) +\left(1-\left(1- \frac{1}{1-\rho_k}\right)^n\right)\right) \quad n \in \mathbb{N}$$
fails to 'recover' the non-trivial zeros. So, $\sigma_r$ seems to be a special case for this method.