Calculating the density of nearest neighbours

Question

I am trying to plot this

which is a numerical simulation of the Montgomery-Odlyzko law for the nontrivial 1st $10^5$ zeros of the Riemann zeta function $ζ(s)$. The solid line is given by 1-(Sin[π x]/(πx))^2.

The following is what I can't seem to interpret:

The blue symbols indicate the pair correlation function of normalized spacings $δn=(γ_{n+1}-γ_n)\log(γ_n/2π)/2π$ between two consecutive nontrivial zeros $1/2+iγ_n$ and $1/2+iγ_{n+1},\ n=1\dots10^5)$ of the Riemann zeta function $ζ(s)$.

from this page. I have looked at various ways of calculating the density, using Histogram, NearestFunction, FindClusters, etc. but am getting nothing like the image above. I realise that this is largely due to a conceptual misunderstanding, but I was hoping someone could point me in the right direction.

Update

plotted with

DeleteCases[Flatten[Table[Abs[z1 - z2] Log[z2/(2 π)]/(2 π), 
{z1, Take[zz,1000]}, {z2, Take[zz,1000]}]], 0];
Show[ListPlot[HistogramList[%, {0, 3, 0.1}][[2]]/200, 
DataRange -> 3], Plot[1 - (Sin[π u]/(π u))^2, {u, 0, 3}]]

using First $10^5$ zeros as zz, thanks to Rahul Narain's comments below. Unfortunately, my computer won't calculate for $10^5$ zeros, but link left for those that will.

... I am now beginning to appreciate the computational power that would have to be used to generate an image like this or this!

Answer 1

I like your questions since we seem to be doing similar things.

Note that your code correctly finds the normalized spacings between all pairs of zeros. This is what the pair correlation function is, but you only want to plot the normalized spacings up to a limit, usually 3. Values returned by your code, MartinPairs[t] where list t contains the imaginary parts of the zeta zeros, are mostly considerably larger than this upper limit.

MartinPairs[t_List] :=
   DeleteCases[Flatten[
      Table[Abs[z1 - z2] Log[z2/(2 Pi)]/(2 Pi), {z1, t}, {z2, t}]],
      0.];

It is much faster to use Differences[t,1,k] to return the first differences of the input list of zeros t with step increment k.

NormedZetaZeroSpacing[t_List, k_Integer] := 
   Differences[t, 1, k] * Log[Drop[t, -k]/(2 Pi)]/(2 Pi)

Years ago, I constructed and stored a list of the first 40000 zeta zeros with Mathematica. The normalized spacings of these zeros, with steps k=1 to k=6, are shown together in the following plot. The calculation took about 5 s on my machine. For comparison, using MartinPairs[t] to find all pair spacings between only the first 1000 zeros took about 7 s.

Histogram[Table[NormedZetaZeroSpacing[t, k], {k, 1, 6}], {0, 5, 0.05},
   Frame -> True, FrameLabel -> {"Zero Spacing", "Number"}]

The k=6 histogram barely shows in the bottom right corner of the above plot.

The match to the theoretical function is shown below.

Histogram[
   Flatten[Table[NormedZetaZeroSpacing[t, k], {k, 1, 8}]], {0, 5, 0.05},
   Frame -> True, FrameLabel -> {"Normalized Spacing  x", "Number"},
   PlotLabel -> "1-(Sin[\[Pi] x]/(\[Pi] x)\!\(\*SuperscriptBox[\()\), \(2\)]\)",
   Epilog -> {Thick, Red,
      Line[Table[{x, 2500 (1 - (Sin[\[Pi] x]/(\[Pi] x))^2)}, {x, 0.1, 5, 0.1}]]}]