# 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. Jun 5, 2014 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 :) Jun 5, 2014 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, 2014 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
– user484
Jun 5, 2014 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.
– user484
Jun 5, 2014 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. Oct 23, 2015 at 0:11
• Works with 2M zeros! image Oct 23, 2015 at 0:15
• Nice. Thanks for the upvote, and the scaling tip. Oct 23, 2015 at 0:17