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.


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!


Posted 2014-06-05T09:57:17.280

Reputation: 7 587

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 – 2014-06-05T10:31:05.397

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 – 2014-06-05T10:33:41.217

I assume you've read Odlyzko's paper - all the (mathematical) details of how such a graph is created are in it. – ciao – 2014-06-05T10:46:17.453

@rasher - do you have a link? – eldo – 2014-06-05T10:53:20.993

@eldo: No, I have it on JSTOR, read it long ago. If you have JSTOR (I'm sure it's elsewhere too), worth a look. – ciao – 2014-06-05T10:55:30.527

In the paper, it refers to the integral of the pair correlation function. Would one then have to find the derivative for each pair? I am sorry, but I am at a bit of a loss :( – martin – 2014-06-05T11:21:53.283

Using your above link and having a look at the image (you show it without its text) it can be seen that 10^5 zeroes were used to compute it. You overlooked the exponentiation operator. – eldo – 2014-06-05T11:32:29.020

Ah - many thanks! – martin – 2014-06-05T11:37:00.463

Would using $10^3$ yield much? – martin – 2014-06-05T11:39:03.790

@martin - better than just 100. Another question is the somehow vague meaning of "normalized spacings". If you could produce a y-axis like in the image you probably would be on the right track. – eldo – 2014-06-05T11:44:41.757

As I understand it, normalised spacings are given by $δn=(γ_{n+1}-γ_n)\log(γ_n/2π)/2π$, no? – martin – 2014-06-05T11:46:40.420

@martin - yes, they are. I overlooked this definition in the text. Sorry. – eldo – 2014-06-05T12:07:10.913

as well studied as the topic is i'd strongly suspect you can readily find tabulated values for the zeros to work with. – george2079 – 2014-06-05T12:34:20.703

@ george2079, That would be great - it would certainly make calculations faster, but I don't know what to do with them once I have them :/ – martin – 2014-06-05T13:09:15.930


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}]

– None – 2014-06-05T19:26:26.767

Thank you :) ... I admit it has had me rather confused! – martin – 2014-06-05T19:28:23.080

Please feel free to add answer & I will accept (though if you wouldn't mind skipping the first sentence, I would be eternally grateful ;)). – martin – 2014-06-05T19:30:45.023

1The 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. – None – 2014-06-05T19:31:15.947

Yes, I figured that :) – martin – 2014-06-05T19:31:49.077

@martin: You can get a text table of the first 100K zeros here, while there, go to home page, follow publications link, you can get a copy of the original paper if you've not already.

– ciao – 2014-06-06T08:27:03.140

Yes, thank you. - The limitations are with my computer now - but I understand the concept. Can't plot with 10^5 on my machine - :/ – martin – 2014-06-06T08:29:31.573



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] :=
      Table[Abs[z1 - z2] Log[z2/(2 Pi)]/(2 Pi), {z1, t}, {z2, t}]],

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"}]

zeta zero correlation histograms

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

The match to the theoretical function is shown below.

   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}]]}]

zeta zero pair correlation


Posted 2014-06-05T09:57:17.280

Reputation: 14 269

works a treat! Note the scaling (where you have 2500 is approx $n/(2\pi^2)$ where $n$ is $#$ of zeros in sample. – martin – 2015-10-23T00:11:06.463

Works with 2M zeros! image

– martin – 2015-10-23T00:15:23.610

1Nice. Thanks for the upvote, and the scaling tip. – KennyColnago – 2015-10-23T00:17:36.293