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. 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$

$...$

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.