Primes out of Thin Air

Left: von Mangoldt additive prime counting function generated using the explicit formula over the Riemann zeta zeros.

Right: Fourier sine transform of the additive von Mangoldt function showing coincidence with the zeta zeros.

The Riemann zeta zeros behave very similarly to a kind of Fourier transform of the primes. This led to the idea of performing a Fourier sine transform of the additive Mangoldt function which is the subject of the explicit prime formula. In his AMS paper Brian Conrey illustrates this process, which in a personal e-mail he attributes to Mike Rubinstein with Mike's quotes "I made the picture to illustrate that one could easily detect the zeros just from numerical data, looking at it on a logarithmic scale."

While writing my research review of the Riemann hypothesis and zeta fractal dynamics I began playing around with alternatives to the von Mangoldt function in the Fourier transform which might give rise to some other zeta-like zero distributions, I accidentally programmed in a simple step function and found it gave a highly smoothed transform, which still had irregular characteristics and a close correspondence to the zeta zeros, from which it was much easier to extract meaningful local minima than the original von mangoldt transform above.

Having got a very close correspondence, with only minor variations, to the zeta zeros, running in a consistent sequence for hundreds of values, I decided to explore what would happen if these values were chosen on the critical line and applied to the explicit formula generating prime counting from zeta zeros (complex value = rho):

To my fascination, the minima produced a meaningful integer-step pseudo-prime distribution, which showed both relative maxima at integer values, as would have perhaps been expected if the explicit formula were to produce any sort of meaningful result, which is itself surprising, but also, reflecting the zeta-approximate irregularities in the transform, it appears to have consistent semi-local peaks at primes and prime powers. So we have generated a distribution back out of the Fourier transform that intimates the primes out of thin air, in the sense that the prime counting function has been removed entirely from the Fourier transform.

Click to enlarge: (Top left) Fourier sin transform of integer step function (black) overlaid on that of the prime counting function (green) and actual zeta zeros (blue). (Top right) Step function minima distribution blue closely coincides with zeta zero distribution (green). (Lower right) Difference between the two. (Lower left) Applying the explicit formula to the step minima at real part 1/2 results in a function with integer local maxima with neighbourhood peaks at primes. (Inset) The step function and summed Mangoldt prime distribution compared.

Anyone familiar with the relationship between the Euler product over primes and the zeta function as a sum
will realize immediately that a similar Euler product over integers would correspond to a sum whose coefficients, instead of being 1, would be the rapidly increasing number of possible factorizations of n, unlikely to give rise to a convergent function with an analytic extension. This makes the fact that the explicit formula does give a meaningful distribution from the integer step function all the more surprising.

Of course the challenge is to you to find if there is a glitch in the calculation causing this phenomenon. I presume it is because the Fourier transform is performing a subtle version of prime sieving, just as the zeta function does with the Euler product. So you can verify it on the spot, I have included the Matlab m-file generating the process below.

function fintstep(scale,start,fin,its,xres)
% fintstep(10^6,1,50,490,600);
% plots Fourier sin transform of integer step distribution giving minima
% coinciding closely with zeta zeros
% 1/log(X) int_0^log(X) sin(ty) (stp(e^y)-e^y)/e^(y/2) dy
% parameters start = first prime, fin = last prime, its = how many zeta zeros to use xres = number of x points in graph
% sav will save the arrays z, x if needed prtits=1 print iterations
y=zeros(1,scale);
for j=2:scale;
y(j)=y(j-1)+mangoldt(j); %make the additive Mangoldt prime counting function
end
clf;
hold on
X=scale;
myspan=800;
xres=4000;
xx=linspace(1,myspan,xres);
stpval=10^4;
F=zeros(1,xres);
for t=1:xres
for yy=0:1/stpval:log(X)
F(t)=F(t)+sin(t*myspan/xres*yy)*(y(floor(exp(yy)))-exp(yy))/exp(yy/2);
end
end
F=F/log(X);
plot(xx(1:xres/10),F(1:xres/10),'g'); %first figure plots Fourier sin transform
%additive Mangoldt transform in green
y=zeros(1,scale);
for j=2:scale;
y(j)=y(j-1)+1; %make the integer step function
end
F=zeros(1,xres);
for t=1:xres
for yy=0:1/stpval:log(X)
F(t)=F(t)+sin(t*myspan/xres*yy)*(y(floor(exp(yy)))-exp(yy))/exp(yy/2);
end
end
F=F/log(X);
plot(xx(1:xres/10),F(1:xres/10),'k'); %integer step function transform in black
zros=textread('zeros1b.txt','%f',30);
for j=1:100 %also plot in blue vertical lines at the actual zeta zeros
if zros(j)>myspan/10 break;
else line([zros(j) zros(j)], [-500 500]);
end
end
hold off
x=[]; %making a vector out of the local minima of the Fourier transform
for k=2:length(F)-1;
if F(k)<=F(k-1)&&F(k)<=F(k+1) x=[x xx(k)]; end
end
xzeros=x(2:length(x)); %xzeros becomes the vector of minima
figure;
zros=textread('zeros1b.txt','%f',length(xzeros));
plot(zros(1:length(xzeros))-xzeros(1:length(xzeros))')
figure;
plot(xzeros); %plot step zeros in blue compared with zeta zeros in green
hold on
plot(zros,'g');
hold off
figure;
if its>length(zros)
its=length(zros);
end
x=linspace(start,fin,xres);
y=zeros(its,xres);
z=y;
for l=1:its
for j=1:xres
p=0.5+xzeros(l)*i;
v=x(j)^(p)/p;
y(l,j)=y(l,j)-v-conj(v); %pair off zeros above and below (conjugates)
if l>1
y(l,j)=y(l,j)+y(l-1,j); %sum the contributions from all the zeros
end
z(l,j)=x(j)+y(l,j)-log(2*pi)-log(1-x(j)^(-2))/2; %the explicit prime formula
end
plot(x,z(l,:),'k');
end
p=primes(fin);
for j=1:length(p)
if (p(j)<=fin) line([p(j) p(j)],[-10 80]); end
end

function o=mangoldt(n)
%von Mangoldt function %o=log(p) if n=p^k or 0 otherwise
f=factor(n);
if f(1)==f(length(f))
o=log(f(1));
else
o=0;
end