# Calculating the density of nearest neighbours

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}][]/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!

• I like your question (+1) because of interesting topic even though I think you are asking about a trivial Mathematica problem. Nevertheless I have no time to answer. – Artes Jun 5 '14 at 10:31
• Thank you. I realise it is trivial in terms of Mathematica computation - sorry! I would really appreciate your assistance if you manage to get time at some point :) – martin Jun 5 '14 at 10:33
• I assume you've read Odlyzko's paper - all the (mathematical) details of how such a graph is created are in it. – ciao Jun 5 '14 at 10:46
• I suggest you look up the actual definition of the pair correlation function. Wikipedia: "The radial distribution function (or pair correlation function) is usually determined by calculating the distance between all particle pairs and binning them into a histogram" (emphasis mine). For example, just taking the first 1000 points (because I didn't want to wait a long time): zzd = Flatten@Table[Abs[zz[[i]] - zz[[j]]], {i, 2, 1000}, {j, 1, i - 1}]; Histogram[zzd, {0, 3, 0.2}] i.stack.imgur.com/Pc3Tw.png – Rahul Jun 5 '14 at 19:26
• The histogram is supposed to then be "normalized with respect to an ideal gas, where particle histograms are completely uncorrelated" so that the horizontal asymptote is at 1 instead of at about 140 in my plot, but I didn't get around to computing that. – Rahul Jun 5 '14 at 19:31

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}]]}] • works a treat! Note the scaling (where you have 2500 is approx $n/(2\pi^2)$ where $n$ is $\#$ of zeros in sample. – martin Oct 23 '15 at 0:11
• Works with 2M zeros! image – martin Oct 23 '15 at 0:15
• Nice. Thanks for the upvote, and the scaling tip. – KennyColnago Oct 23 '15 at 0:17