Primes and the Zeta Function

Math 242 Modern Computational Mathematics

The Riemann Zeta Function

In the homework for today we encountered the Riemann zeta function. One of the most important functions in number theory, the zeta function is defined for s>1 by:

ζ(s)=(∑_(n=1)^∞)(1/(n^s))

The Mathematica function Zeta[s] returns values ζ(s).

Your Turn: Evaluate ζ(s) for positive and negative integers s.

Zeta[2]

N[Zeta[2]]

Negative even integers: trivial zeros

Table[Zeta[k],{k,-20,-2,2}]

Negative odd integers: values alternate positive and negative

Table[Zeta[k],{k,-19,-1,2}]

Positive odd integers:

Table[Zeta[k],{k,1,20,2}]//N

Positive even integers :

Table[Zeta[k],{k,2,20,2}]

Visualizing ζ(s)

The Riemann zeta function extends to the entire complex plane, except for the point s=1, by the process of analytic continuation. However, it is tough to visualize graphs of complex-valued functions. Since the domain and range each have two real dimensions, plotting the graph of a complex-valued function would require four dimensions! Since we don’t have four spatial dimensions, we’ll have to do something else.

Domain Coloring

One common approach for visualizing a function f:ℂ→ℂ is to color each point z of the complex plane according to the value of f(z). For example, if f(z)=r (e^(i θ)), we could color z such that θ determines the hue of the color and r determines the brightness of the color.

The next plot shows the complex plane colored according to this scheme. First, we define a “paint” function, and then we produce the plot.

paint[z_]:=Module[{color, shade},
color=Hue[Rescale[Arg[-z],{-Pi,Pi}]];
shade=Exp[-Abs[z]]/1.3;
Darker[color,shade]
]

ParametricPlot[{x,y},{x,-5,5},{y,-5,5},ColorFunction->Function[{x,y},paint[x+y I]],ColorFunctionScaling->False,PlotPoints->{100,100},PlotRangePadding->0,ImageSize->250,Axes->False,FrameLabel->{"real","imaginary"}]

Note that in this plot, the positive real numbers are red (■) and negative real numbers are cyan (■). Numbers on the positive imaginary axis are green (■), and numbers on the negative imaginary axis are purple (■). Darker colors mean that the number is closer to zero (that is, smaller in absolute value). Points very close to the origin are almost black.

Domain coloring the Riemann zeta function

Now we will color the complex plane by the values of the Riemann zeta function.

ParametricPlot[{x,y},{x,-20,20},{y,-30,30},ColorFunction->Function[{x,y},paint[Zeta[x+y I]]],ColorFunctionScaling->False,PlotPoints->{200,300},PlotRangePadding->0,ImageSize->300,Axes->False,FrameLabel->{"real","imaginary"}] // Quiet

From the plot, we can see the trivial zeros on the negative real axis. We can also see the “pole” of the function, where the function blows up to infinity, at (1,0). These features are easier to see if we zoom in.

In the following plot, observe the small black regions at negative even integers. These are the “trivial” zeros of the Riemann zeta function.

ParametricPlot[{x,y},{x,-24,4},{y,-4,4},ColorFunction->Function[{x,y},paint[Zeta[x+y I]]],ColorFunctionScaling->False,PlotPoints->{280,80},PlotRangePadding->0,ImageSize->500,Axes->False,FrameLabel->{"real","imaginary"}] // Quiet

Here is a close-up plot of the pole at (1,0). There is no black dot here, because the values of the function get huge rather than going to zero. Hence, the point (1,0) is a pole, not a zero.

ParametricPlot[{x,y},{x,0.5,1.5},{y,-0.5,0.5},ColorFunction->Function[{x,y},paint[Zeta[x+y I]]],ColorFunctionScaling->False,PlotPoints->{100,100},PlotRangePadding->0,ImageSize->300,Axes->False,FrameLabel->{"real","imaginary"}] // Quiet

The Riemann hypothesis says that the non-trivial zeros of ζ(s) all have real part (1/2). This is considered to be one of the greatest unsolved problems in mathematics. More than 1 trillion nontrivial zeros of ζ(s) have been computed numerically, and they all have real part (1/2). The hypothesis appears to be true, but no one has found a proof.

In our domain coloring, we can see some nontrivial zeros on the line Re(s)=(1/2). Again, zoom in:

ParametricPlot[{x,y},{x,-4,6},{y,-2,28},ColorFunction->Function[{x,y},paint[Zeta[x+y I]]],ColorFunctionScaling->False,PlotPoints->{100,300},PlotRangePadding->0,ImageSize->200,Axes->False,FrameLabel->{"real","imaginary"}] // Quiet

Nontrivial Zeros of ζ(s)

Computing the nontrivial zeros of ζ(s) is challenging. (Shall we say it’s “nontrivial”?) Fortunately, Mathematica provides a built-in function, ZetaZero[k], that gives the (k^th) nontrivial zero of the Riemann zeta function. For example:

N[ZetaZero[1]]

As expected, the real part of this zero is 0.5. The imaginary part is about 14.1347. Does this agree with the location of a zero in the domain coloring plot above?

We can also ask Mathematica for a bunch of nontrivial zeros:

ZetaZero /@Range[10] //N

Do these zeros agree with what you observe in the domain coloring plot?

For one further visualization, we can plot the real and imaginary parts of the Riemann Zeta function ζ(s) along the line Re(s)=(1/2).

Plot[{Re[Zeta[1/2+I*t]],Im[Zeta[1/2+I*t]]},{t,0,50},PlotLegends->{"real","imaginary"},AspectRatio->1/3]

The zeros of ζ(s) occur where its real and imaginary parts are simultaneously equal to zero. Confirm that the location of zeros given by the plot above agree with the values returned by ZetaZero. To do this, make a list of the imaginary parts of the nontrivial zeros:

Table[Im[ZetaZero[k]],{k,1,10}]//N

Zeta Zeros and Primes

A fascinating pattern emerges if we add up cosine functions whose “frequencies” are the imaginary parts of the nontrivial zeros. Let (Subscript[θ, k]) be the imaginary part of the kth nontrivial zero. That is, (Subscript[θ, 1])=14.143..., (Subscript[θ, 2])=21.022..., and so on. Then consider the function:

f(x)=-(∑_(k=1)^∞)cos(ln(x)(Subscript[θ, k]))

Your Turn: Write a Mathematica function that computes the partial sum of f(x) using a specified number of terms. You may use the built-in Sum function if you like.

f[x_,numTerms_]:= -1.0*Sum[Cos[Log[x]*Im[ZetaZero[k]]],{k,1,numTerms}]

f[4,20]

Your Turn: What does the graph of f(x) look like? Plot partial sums of f(x) on an interval such as 2<=s<=30. What do you observe?

Plot[f[x,200],{x,1,30},PlotRange->{-10,50}]

An approximation to the prime counting function

Now let’s return to the prime counting function π(x), which counts the number of primes less than or equal to a real number x. We will see one way in which the Riemann zeta function ζ(x) is connected to the prime counting function π(x).

First we define Riemann’s R-function R(x), which involves the Riemann zeta function.

R(x)=1+(∑_(k=1)^∞)((((ln x))^k)/(k k! ζ(k+1)))

It turns out that R(x) is a very good approximation to the prime-counting function. To see this, we would like to plot R(x) and π(x). Since R(x) is an infinite sum, we will instead use a partial sum with the first m terms; the next line implements this:

Your Turn: Write a Mathematica function that computes the partial sum of R(x) using a specified number of terms.

riemannR[x_,numTerms_]:=1+Sum[((Log[x]^k)/(k*k!*Zeta[k+1])),{k,1,numTerms}]

riemannR[6,20]//N

Now plot the prime counting function π(x) and our approximation consisting of the first 20 terms of R(x):

max=100;
Show[
Plot[ PrimePi[x], {x,0,max} ],
Plot[ riemannR[x,20],{x,1,max},PlotStyle->Red]
]

A very cool approximation

The second Chebyshev function ψ(x) is a variant on the prime counting function. It counts not only primes, but prime powers. The function is:

ψ(x)=(∑_((p^k)≤ x))ln(p)=(∑_(p≤x))⌊(Subscript[log, p])(x)⌋ln(p)

The first sum above is over all prime powers not exceeding x; the second sum is over all primes not exceeding x.

Your Turn: Write a Mathematica function that computes ψ(x).

ψ[x_]:=Sum[ If[PrimeQ[p],Floor[((Log[x])/(Log[p]))]*Log[p],0],{p,2,x} ]

ψ[6]

Your Turn: Now make a plot of ψ(x). Carefully observe where the graph of ψ(x) has vertical steps.

psiPlot=Plot[ψ[x],{x,0,50}]

Notice that the second Chebyshev function is not continuous. However, it can be approximated by a continuous function that involves the nontrivial zeros of the Riemann zeta function! The approximation function is:

(Subscript[ψ, 0])(x)=x-ln(2π)-(1/2)ln(1-(x^(-2)))-(∑_ρ)((x^ρ)/ρ)

Here, the sum is over all nontrivial zeros ρ of the Riemann zeta function. Since we can’t sum over infinitely many zeros, we will define (Subscript[ψ, 0]) in Mathematica using a finite sum, with a parameter m that determines how many zeros we use. Note that we are using ZetaZero to obtain values of the zeros.

ψ0[x_,m_]:= x-Log[2Pi]-Log[1-x^(-2)]/2-Sum[2*Re[(x^ZetaZero[k])/ZetaZero[k]],{k,1,m}]

Here is a plot, with (Subscript[ψ, 0]) evaluated using 20 zeros of the Riemann zeta function:

Plot[ψ0[x,20],{x,0,50},PlotStyle->Purple]

The function (Subscript[ψ, 0]) converges remarkably to ψ as we increase the number of terms in the sum. Here is an interactive window with a slider that controls the number of terms in the sum. (Unfortunately the response time when moving the slider might be slow because plotting our function (Subscript[ψ, 0]) is computationally intensive!)

Manipulate[ Show[psiPlot, Plot[ψ0[x,m],{x,0,50},PlotStyle->Purple]], {m,0,60,1}]

Move the m slider to the right and watch the approximation (in purple) get closer to the second Chebyshev function ψ (in orange). Incredible! As we include more and more zeros of the Riemann zeta function, the continuous function (Subscript[ψ, 0]) converges to the prime-counting function variant ψ!