The inverse Riemann zeta function

The Riemann zeta function is classically defined by an infinite series

which is absolutely convergent $\Re(s)>1$, where $s=\sigma+it$ is a complex variable. The values for the first few special cases are:

and so on. For $s=1$, the series diverges asymptotically as $\gamma+\log(k)$, where $\gamma=0.5772156649\dots$ is the Euler-Mascheroni constant. The special values for even positive integer argument are generated by the Euler’s formula

for which the value is expressed as a rational multiple of $\pi^{2k}$ where the constants $B_{2k}$ are Bernoulli numbers defined such that $B_{0}=1$, $B_{1}=-\frac{1}{2}$, $B_{2}=\frac{1}{6}$, and so on. For odd positive integer argument, the values of $\zeta(s)$ converge to unique constants, which are not known to be expressed as a rational multiple of $\pi^{2k+1}$ as occurs in the even positive integer case. For $s=3$, the value is commonly known as Apéry’s constant, who proved its irrationality [2]. The key connection to prime numbers is by means of Euler’s infinite product formula

where $p_{1}=2$, $p_{2}=3$, $p_{3}=5$ and so on, denote the prime number sequence. These prime numbers can be recursively extracted from the Euler product by the Golomb’s formula [9]. Hence, if we define a partial Euler product up to the $n$th order as

for $n>1$ and $Q_{0}(s)=1$, then we obtain a recurrence relation for the $p_{n+1}$ prime

This leads to representation of primes by the following limit identities

and so on, whereby all the previous primes are used to excite the Riemann zeta function in a such as way as to extract the next prime. A detailed proof and numerical computation is shown in [11][12]. We will find that this recursive structure (when taken in the limit as $s\to\infty$) is a basis for the rest of this article, and will lead to formulas for trivial and non-trivial zeros, and the inverse Riemann zeta function.

Furthermore, the Riemann zeta series (1) induces a general Weierstrass factorization of the form

where it analytically extends the zeta function to the whole complex plane and reveals its full structure of the poles and zeros [1, p.807]. Only a simple pole exists at $s=1$, hence $\zeta(s)$ is convergent everywhere else in the complex plane, i.e., the set $\mathbb{C}\backslash 1$. Moreover, there are two kinds of zeros classified as the trivial zeros $\rho_{t}$ and non-trivial zeros $\rho_{nt}$. The first infinite product term of (8) encodes the factorization due to trivial zeros

for $n\geq 1$ which occur at negative even integers $-2,-4,-6\ldots$, where $n$ is index variable for the $n$th zero. The second infinite product term of (8) encodes the factorization due to non-trivial zeros, which are complex numbers of the form

and, as before, $n$ is the index variable for the $n$th zero (this convention is straightforward if the zeros are on the critical line). But in general, the real components of non-trivial zeros are known to be constrained to lie in a critical strip in a region between $0<\sigma_{n}<1$. It is also known that there is an infinity of zeros located on the critical line at $\sigma=\frac{1}{2}$, but it is not known whether there are any zeros off of the critical line, a problem of the Riemann hypothesis (RH), which proposes that all zeros should lie on the critical line. The first few zeros on the critical line at $\sigma_{n}=\frac{1}{2}$ have imaginary components $t_{1}=14.13472514...$, $t_{2}=21.02203964...$, $t_{3}=25.01085758...$, and so on, which were computed by an analytical recurrence formula as

as we have shown in [12][13] assuming (RH). Also, $\zeta(s,a)$ is the Hurwitz zeta function

which is a shifted version of (1) by an arbitrary parameter $a\neq-1,-2,\ldots$. Now, substituting the Weierstrass infinite product (8) into the Golomb’s formula (6) can be used to generate primes directly from non-trivial zeros. In later sections, we will show that non-trivial zeros can also be generated directly from primes.

Furthermore, the Riemann zeta function can have many points $s_{n}$ such that $w=\zeta(s_{n})$ can map to the same $w$ value. In Figure 1 we plot $\zeta(s)$ for real $s$ and $w$, and note that for $s>1$ the function is monotonically decreasing from $+\infty$ and tends $O(1)$ as $s\to\infty$, and for the domain $-2.7172628292\ldots<s<1$ it’s monotonically decreasing from $0.0091598901\ldots$ to $-\infty$, as also shown in Figure 2. And for $s<-2.7172628292\ldots$, it becomes oscillatory where there are many $s_{n}$ solutions. For example, the first two $s$ values

map to the same $w$ value as shown in Figure 2. It is usually customary to report that $\zeta(-3)=\frac{1}{120}$, but in fact is the second solution $s_{2}$, the first solution, or the principal solution $s_{1}$, is the value $-2.47270347\ldots$.

In this article, we seek to develop a formula for an inverse Riemann zeta function $s=\zeta^{-1}(w)$. In general, the existence of an inverse is established by an inverse function theorem as shown in [14, p.135], where for any holomorphic function $f(z)$, an inverse exists when $f^{\prime}(z)\neq 0$ in the $z$-domain. Hence, for

there is an inverse function

which implies that

and

for some real and complex domains $w$ and $s$, except at those points where $\zeta^{\prime}(s)=0$ in the $s$-domain. The inverse zeta is a multi-valued function with infinite number of branches, and $s_{n}$ are the multi-valued solutions for which $w=\zeta(s_{n})$. The presented method can also recursively access these multi-valued solutions of $s_{n}=\zeta^{-1}(w)$, but as we will find, the computational precision becomes very high, and so we will primarily focus on the principal solution $s_{1}$, which as we will find is very easy to compute and will cover almost the entire complex plane. Additionally, in [8] discusses additional theoretical basis behind solutions to (14), also known as $a$-points, which we refer to as $s_{n}$.

We develop a recursive formula for an inverse Riemann zeta function as:

where $s_{n}$ is the value for which $w=\zeta(s_{n})$ for some branch $\lambda$ of the $m$th root; it will recursively generate solutions $s_{n}$ for a given $w$-domain. The formula for the principal branch is

and is valid for all complex $w$-domain except in a small strip region $\Re(w)\in(j_{1},1)\cup\{\Im(w)=0\}$ that we determined experimentally, and $j_{1}=0.00915989\ldots$ is a constant. As we will show in more details in later sections, this formula can easily invert the Basel problem

or the Apéry’s constant

and essentially the entire complex domain $w\in\mathbb{C}$ with an exception of a strip region $\Re(w)\in(j_{1},1)\cup\{\Im(w)=0\}$ where there lie (possibly) an infinite number of branch non-isolated singularities.

The way in which we will arrive at the presented formula is by connecting two simple Theorems. The first Theorem 1, as presented in Section 2, is the $m$th log-derivative formula for obtaining a generalized zeta series over the zeros of a function. Such a method appears in the literature from time to time and can be traced back to Euler, who, according to a reference in [21, p.500], used it to devise a means of computing several zeros of the Bessel function of the first kind by solving a system of equations generated by the log-derivative formula. In the works of Voros, Lehmer, and Matsuoka, it has been used to find a closed-form formula for the secondary zeta functions [15][16][19][20]. Then, the second Theorem 2, as given in Section 3, develops a method for finding a zeta recurrence formula for the $n$th+1 term of a generalized zeta series, i.e., all terms of a generalized zeta series must be known in order to generate the $n$th+1 term, as we have shown in our previous work [12][13]. In Section 4, we connect these two theorems and find formulas for trivial and non-trivial zeros of the Riemann zeta function and explore their properties in greater detail. We then develop a formula for an inverse zeta and also explore some its identities. We empirically observe that there are infinitely many branch singularities of the inverse zeta that are spread out along a narrow strip $(j_{1},1)$ forming branch singularity attractor. And in the last part, we briefly extend the presented results to find zeros (and inverses) of many common functions, such as the gamma function, the Bessel function of the first kind, the trigonometric functions, Lambert-W function, and any entire function or finite/infinite degree polynomial or rational function (provided that the function fits the constraints of this method).

Throughout this article, we numerically compute the presented formulas to high precision in PARI/GP software package [17], as it is an excellent platform for performing arbitrary precision computations, and show that these formulas do indeed converge to the inverse of the Riemann zeta function to high precision. We note that when running the script in PARI, the precision has to be set very high (we generally set precision to 1000 decimal places). Also, the Wolfram Mathematica software package [22] was instrumental in developing this article.

In this section, we outline the zeta $m$th log-derivative method. We recall the argument principle, that for an analytic function $f(z)$ we have

which equals the number of zeros $N_{z}$ minus the number of poles $N_{p}$ (counting multiplicity) which are enclosed in a simple contour $\Omega$. But now, instead working with the number of zeros or poles, the aim of the $m$th log-derivative method is to find a generalized zeta series over the zeros and poles of the function $f(z)$ in question. Hence if $Z=\{z_{1},z_{2},z_{3},\ldots,z_{N_{z}}\}$ is a set of all zeros of $f(z)$ and $P=\{p_{1},p_{2},p_{3},\ldots,p_{N_{p}}\}$ is a set of all poles of $f(z)$ in the whole complex plane, then the generalized zeta series are

and

The number of zeros or poles may be finite or infinite, but in the latter case, the convergence of a generalized zeta series to the right-half side of the line $\Re{(s)}>\mu$ may depend on the distribution of its terms. If the number of zeros and poles is finite, then one could count them by evaluating $Z(0)-P(0)=N_{z}-N_{p}$, which reduces to the argument principle (22). The argument principle may be modified further by introducing another function $h(z)$ as such

The proof is a slight modification of a standard proof of (22) using the Residue theorem, and if we let $h(z)=z^{-s}$, then we have

but the contour has to be deformed as to not encircle the origin. However, in the following sections, we will not pursue (26) as the integral makes the study and computation more difficult, as well as being dependent on the contour. There is a simpler variation of (26) which is better suited for our study. The main formula of the zeta $m$th log-derivative method is:

We also remark if the number of zeros or poles is infinite, then the Weierstrass factorization product naturally imposes that $Z(m)$ and $P(m)$ are absolutely convergent, which is very critical for next section.

Additionally, we can obtain an integral representation using the Cauchy integral formula applied to coefficients of Taylor expansion (34) as

where $\Omega_{0}$ is a simple contour encircling the origin, but unlike (26), is not enclosing zeros or poles. In order to apply (34) and (35) successfully, one has to judiciously choose $g(z)$ as to remove it from $f(z)$ so that the $m$th log-differentiation will not produce unwanted artifacts due to $g(z)$, as well as the zero or pole at origin.

We now outline a method to extract the terms of a generalized zeta series by means of a recurrence formula satisfied by the terms of such series. Hence, if the terms of a generalized zeta series are to be zeros of a function, then such a method effectively gives a recurrence formula satisfied by the zeros, which in turn can be used to compute the zeros. And similarly, the same holds if the terms of a generalized zeta series are poles of a function. In fact, any quantities represented by the terms of a generalized zeta series can be recursively found. In [12], we developed a formula for the $n$th+1 prime based on the prime zeta function. But for the purpose of this paper, let us consider the generalized zeta series over zeros $z_{n}$ of a function

and let us also assume that the zeros are positive, real, and ordered from smallest to largest such that $0<z_{1}<z_{2}<z_{3}<\ldots<z_{n}$, then the asymptotic relationship holds

as $s\to\infty$. To illustrate how fast the terms decay, let us take $z_{1}=2$ and $z_{2}$=3. Then for $s=10$ we compute

where we roughly see an order of magnitude difference. But for $s=100$ we compute

where we see a difference by $17$ orders of magnitude. Hence, as $s\to\infty$, then

as the $z^{-s}_{n+1}$ term completely vanishes in relation to $z^{-s}_{n}$, and so the $z^{-s}_{n}$ term dominates the limit. As a result, we write

From this we have

as $s\to\infty$ where the lowest order term dominates, and we refer to it as the principal term, or in the case where the zeta series are considered to be zeros of a function, the principal zero. This rapid decay of higher order zeta terms (41) opens a possibility for a recursive root extraction as shown by Theorem 2 next.

Let us next give an example of how Theorem 1 and Theorem 2 is applied to find a formula for zeros for any meromorphic function. For simplicity, suppose that we wish to find zeros of the function

We know in advance that the zeros are just integer multiples: $z_{n}=\pm n$ (for any non-zero integer $n=1,2,3\ldots$) and zero at origin is already canceled. But now, if we apply the $m$th log-derivative formula to (52), then we get a generalized zeta series of over the zeros as

which is equal to the zeta series over all zeros (including the negative ones) as such

From this we can deduce that even values become double of a half the side of zeros (which in this example is equivalent to $\zeta(s)$) as

and the odd values cancel

In view of this, the Euler’s formula (3) is

is an example of a closed-form representation of $\zeta(2k)$ that is not involving zeros directly, i.e., the positive integers. But such a formula may not always be available, so we will just use the $m$th log-derivative formula (53) as the main closed-form representation of $Z(s)$. Also, on a side note, the Bernoulli numbers are expansions coefficients of another function

where they can be similarly obtained by $m$th differentiation:

Hence in essence, the Euler’s formula (57) is just a transformed version of (53), but it just happens that Bernoulli numbers are rational constants. And so, when putting this together, we obtain a full solution to (52) as a recurrence relation

where a $2m$ limit value ensures that it is even, and so, an additional factor of $\frac{1}{2}$ is added due to (54). The principal zero is:

and a numerical computation for $m=20$ yields

which is accurate to 13 digits after the decimal place, and the script to compute it in PARI is presented in Listing 1. The key aspect of the script is the derivnum function for computing the $m$th derivative very accurately which will be very useful for the rest of this article. The next zero is recursively found as

but we must know the first zero in advance in order to compute the next zero. A numerical computation for $m=20$ yields

which is accurate to 8 digits after the decimal place. Then, the next zero in the sequence is

but we must know the first two zeros in advance in order to compute the next zero. A numerical computation for $m=20$ yields

which is accurate to 6 digits after the decimal place, and so on. We see that the accuracy becomes lesser for higher zeros, and so the limit variable $m$ has to be increased to get better accuracy. One can then continue this process and extract the $z_{n+1}$ zero.

As a second example, we find a formula for zeros of the Bessel function of the first kind

for all $\nu>-1$ real orders. We denote the zeros of (67) as $x_{\nu,n}$ where $n$ is the index variable for the $n$th zero for $\nu$ order of the Bessel function. The Weierstrass product representation of (67) is

for $\nu>-1$ involving the zeros directly [1, p.370]. To find the roots of this function, we apply Theorem 1 to obtain a generalized zeta series over the Bessel zeros

which is essentially the $2m$th derivative of $\log[J_{\nu}(x)/x^{\nu}]$ evaluated at $x=0$, where it is taken in a limiting sense $x\to 0$ as to avoid division by zero. Here we cancel the zero $x^{\nu}$ at origin as shown in previous section. This will prevent any artifacts of $x^{\nu}$ in (68) from being generated by the log-differentiation. Furthermore, the Bessel function is even; hence we consider the $2m$ values.

The first few special values of $Z_{\nu}(2m)$ for even orders are:

and so on. These generated values are rational functions of the Bessel order $\nu>-1$ for even orders, and this implies that if $\nu>-1$ is rational, so is $Z_{\nu}(2m)$. These formulas for $Z_{\nu}(2m)$ were generated using another recurrence relation found Sneddon in [18] as an alternative to (69), which is given in Appendix B. In Table 1, we give the values of $Z_{\nu}(2m)$ for different $m$ and $\nu$.

And now, by applying Theorem 2, we obtain a full recurrence formula satisfied by the Bessel zeros:

To verify (71) numerically, we compute the principal zero using (69), since it is more efficient than (70), for the limit variable $m=250$ and $\nu=0$ which results in

which is accurate to 181 decimal places (we are just showing the first 30 digits). In [21, p.500-503], Rayleigh-Cayley generated values for $Z_{\nu}(2m)$ as shown in (70) and extended Euler’s original work and computed the smallest Bessel zero using this method in papers dating back to year 1874. Moving on, the next zero is found the same way, we compute

for $m=250$, and is accurate to $99$ decimal places, but in order to ensure convergence, the first zero $x_{0,1}$ has to be known to high enough precision (usually much higher than can be efficiently computed using this method as we did above). Henceforth, as a numerical experiment, we take $x_{0,1}$ that was already pre-computed to high enough precision using more efficient means to 1000 decimal places using the standard equation solver found on mathematical software packages (such root finding algorithms are very effective, but must assume an initial condition), rather than taking the zero computed above with less accuracy. And similarly, the third zero is computed as

which is accurate to $68$ decimal places and that it was assumed $x_{0,1}$ and $x_{0,2}$ was already pre-computed to high enough precision (1000 decimal places using the standard equation solver) in order to ensure convergence. Hence in general, one can continue on and keep removing all the known zeros up to the $n$th order in order to compute the $n$th+1 zero. In numerical computations, the key here is that the accuracy of the previous zeros must be much higher than the next zero in order to ensure convergence, i.e., $x_{\nu,n}^{-s}\gg x_{\nu,n+1}^{-s}$, and also one cannot use the same limit variable to compute the next zero based on the previous zero as it will cause a self-cancelation in the formula. Numerically, there is a fine balance as to how many accurate digits are available and the magnitude of the limit variable $m$ used to compute the next zero. We also note that this method is not numerically an efficient method to compute zeros, but it allows to have a true closed-form representation of the zeros, and also that one does not need to make an initial guess for the zero, as is generally the case for many root finding algorithms.