Send me back to Mathematics Stack Exchange if this question is not research level.
At Terence Tao's blog post there is the expression:
$$\sum\limits_{n \leq X} \Lambda(n)\Lambda(n+h) \ \ \ \ \ \ \ \ \ \ \ \ \ (1)$$
It has been proven that the von Mangoldt function is:
$$\Lambda(n)=\lim\limits_{s \rightarrow 1} \zeta(s)\sum\limits_{d|n} \frac{\mu(d)}{d^{(s-1)}}$$
which for each $n$ has an associated Dirichlet series.
Therefore I made a program that looks at the resulting matrix of Dirichlet series of the convolutions:
$$\Lambda(n)\Lambda(n+h) = \left(\lim\limits_{s \rightarrow 1} \zeta(s)\sum\limits_{d|n} \frac{\mu(d)}{d^{(s-1)}}\right)\left( \lim\limits_{s \rightarrow 1} \zeta(s)\sum\limits_{d|n+h} \frac{\mu(d)}{d^{(s-1)}}\right)$$
Example: For $h=2$ this is a matrix starting:
$$M=\left( \begin{array}{cccccc} 1 & 1 & 1 & 1 & 1 & 1 & \dots \\ 2 & -2 & 2 & -2 & 2 & -2 & \\ -1 & 2 & -1 & -1 & 2 & -1 & \\ 3 & -1 & 3 & -1 & 3 & -1 & \\ 2 & 2 & -3 & 2 & -3 & 2 & \\ -2 & -4 & -2 & 2 & 4 & 2 & \\ 2 & 2 & 2 & 2 & -5 & 2 & \\ 4 & 0 & 4 & 0 & 4 & 0 & \\ \vdots & & & & & & \ddots \end{array} \right)$$
It appears that the rows in matrix $M$ are periodic and it therefore suffices to look at the row sums of the lower triangular matrix:
$$R=\left( \begin{array}{cccccc} 1 & & & & & & \dots \\ 2 & -2 & & & & & \\ -1 & 2 & -1 & & & & \\ 3 & -1 & 3 & -1 & & & \\ 2 & 2 & -3 & 2 & -3 & & \\ -2 & -4 & -2 & 2 & 4 & 2 & \\ \vdots & & & & & & \ddots \end{array} \right)$$
and compare them for different $h$.
In the program below I have compared $h=2$ to $h=2,3,4,5,6,7,8,...,32$ and I notice that the row sums of $R$ are the same for $h$ a power of $2$ although the entries in matrices $M$ and $R$ are different for different $h$.
(*start*)
Clear[ConvolvedMangoldt];
Monitor[Table[h = 2;
nn = 140;
A = Table[Table[If[Mod[n, k] == 0, 1, 0], {k, 1, nn}], {n, 1, nn}];
B = Table[
Table[If[Mod[k, n] == 0, MoebiusMu[n]*n, 0], {k, 1, nn}], {n, 1,
nn}];
MatrixForm[T = (A.B)];
TableForm[ConvolvedMangoldt = Table[a = T[[All, kk]];
F1 =
Table[Table[If[Mod[n, k] == 0, a[[n/k]], 0], {k, 1, nn}], {n, 1,
nn}];
b = T[[All, kk + h]];
F2 =
Table[Table[If[Mod[n, k] == 0, b[[n/k]], 0], {k, 1, nn}], {n, 1,
nn}];
(F1.F2)[[All, 1]], {kk, 1, nn - h}]];
Flatten[
Table[Table[ConvolvedMangoldt[[n - k + 1, k]], {k, 1, n}], {n, 1,
nn - h}]];
TableForm[
Table[Table[ConvolvedMangoldt[[n - k + 1, k]], {k, 1, n}], {n, 1,
nn - h}]];
aa = Table[
Sum[ConvolvedMangoldt[[k, n]], {k, 1, n}], {n, 1, nn - h}];
h = hh;
nn = 140;
A = Table[Table[If[Mod[n, k] == 0, 1, 0], {k, 1, nn}], {n, 1, nn}];
B = Table[
Table[If[Mod[k, n] == 0, MoebiusMu[n]*n, 0], {k, 1, nn}], {n, 1,
nn}];
MatrixForm[T = (A.B)];
TableForm[ConvolvedMangoldt = Table[a = T[[All, kk]];
F1 =
Table[Table[If[Mod[n, k] == 0, a[[n/k]], 0], {k, 1, nn}], {n, 1,
nn}];
b = T[[All, kk + h]];
F2 =
Table[Table[If[Mod[n, k] == 0, b[[n/k]], 0], {k, 1, nn}], {n, 1,
nn}];
(F1.F2)[[All, 1]], {kk, 1, nn - h}]];
Flatten[
Table[Table[ConvolvedMangoldt[[n - k + 1, k]], {k, 1, n}], {n, 1,
nn - h}]];
TableForm[
Table[Table[ConvolvedMangoldt[[n - k + 1, k]], {k, 1, n}], {n, 1,
nn - h}]];
bb = Table[
Sum[ConvolvedMangoldt[[k, n]], {k, 1, n}], {n, 1, nn - h}];
{" hh = ", hh,
aa[[Range[100]]] - bb[[Range[100]]]}, {hh, 2, 32}], hh]
(*end*)
I have put the output here at Pastebin: https://pastebin.com/WiU9GSm0
The program explains better my comparison, but I don't have the intelligence to write out in LaTex what the program compares.
But the question is the same anyways. Is it known that powers of $h$ a prime gap, have the same density?
The program is written in Mathematica.
Edit 18.2.2018:
An improved version of the program where one can vary the prime gap parameter "p".
(*start*)
p = 6
Do[
nn = 140;
h = p; (* h = 2 gives twin primes *)
hh = hhh; (* hh = 4 gives cousin primes *)
TableForm[
A = Table[Table[If[Mod[n, k] == 0, 1, 0], {k, 1, nn}], {n, 1, nn}]];
TableForm[
B = Table[
Table[If[Mod[k, n] == 0, MoebiusMu[n]*n, 0], {k, 1, nn}], {n, 1,
nn}]];
TableForm[T = (A.B)];
(*Print["h = ",h]*)
TableForm[ConvolvedMangoldt = Table[
a = T[[All, nnn]];
b = T[[All, nnn + h]];
Table[
Sum[If[Mod[n, k] == 0, a[[n/k]]*b[[k]], 0], {k, 1, nn - h}], {n,
1, nn - h}], {nnn, 1, nn - h}]];
TableForm[CMT = Transpose[ConvolvedMangoldt]];
TableForm[
AA = Table[
Sum[CMT[[n, k]], {k, 1, n}], {n, 1, Length[CMT] - (hh - h)}]];
(*TableForm[AAA=Table[Sort[Table[CMT[[n,k]],{k,1,n}]],{n,1,Length[\
CMT]-(hh-h)}]];*)
(*Print["h = ",hh]*)
h = hh;
TableForm[ConvolvedMangoldt = Table[
a = T[[All, nnn]];
b = T[[All, nnn + h]];
Table[
Sum[If[Mod[n, k] == 0, a[[n/k]]*b[[k]], 0], {k, 1, nn - h}], {n,
1, nn - h}], {nnn, 1, nn - h}]];
TableForm[CMT = Transpose[ConvolvedMangoldt]];
TableForm[
BB = Table[Sum[CMT[[n, k]], {k, 1, n}], {n, 1, Length[CMT]}]];
(*TableForm[BBB=Table[Sort[Table[CMT[[n,k]],{k,1,n}]],{n,1,Length[\
CMT]}]];*)
Print["h = ", hh, " ", Total[BB - AA]]
(*Print["The difference between the sorted rows are zero"]*)
(*BBB-AAA;*)
, {hhh, p, 100}]
(*end*)
Edit 19.2.2018:
Anyways,
$$\Lambda(n)\Lambda(n+h) = \left(\lim\limits_{s \rightarrow 1} \zeta(s)\sum\limits_{d|n} \frac{\mu(d)}{d^{(s-1)}}\right)\left( \lim\limits_{s \rightarrow 1} \zeta(s)\sum\limits_{d|n+h} \frac{\mu(d)}{d^{(s-1)}}\right)$$
where $n$ is unconventionally the column index.
For $h=2$, the convolution above is associated with the matrix of Dirichlet series:
$$M_{h=2}=\left(\begin{array}{lllllllll} 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & \cdots\\ 2 & -2 & 2 & -2 & 2 & -2 & 2 & -2 \\ -1 & 2 & -1 & -1 & 2 & -1 & -1 & 2 \\ 3 & -1 & 3 & -1 & 3 & -1 & 3 & -1 \\ 2 & 2 & -3 & 2 & -3 & 2 & 2 & -3 \\ -2 & -4 & -2 & 2 & 4 & 2 & -2 & -4 \\ 2 & 2 & 2 & 2 & -5 & 2 & -5 & 2 \\ 4 & 0 & 4 & 0 & 4 & 0 & 4 & 0 \\ \vdots & & & & & & & &\ddots \\ \end{array}\right)$$
$R_{h=2}$ is equal to $M_{h=2}$ for $\text{ row index } \geq \text{ column index}$:
$$R_{h=2}=\left(\begin{array}{lllllllll} 1 & \text{} & \text{} & \text{} & \text{} & \text{} & \text{} & \text{} &\cdots \\ 2 & -2 & \text{} & \text{} & \text{} & \text{} & \text{} & \text{} \\ -1 & 2 & -1 & \text{} & \text{} & \text{} & \text{} & \text{} \\ 3 & -1 & 3 & -1 & \text{} & \text{} & \text{} & \text{} \\ 2 & 2 & -3 & 2 & -3 & \text{} & \text{} & \text{} \\ -2 & -4 & -2 & 2 & 4 & 2 & \text{} & \text{} \\ 2 & 2 & 2 & 2 & -5 & 2 & -5 & \text{} \\ 4 & 0 & 4 & 0 & 4 & 0 & 4 & 0 \\ \vdots & & & & & & & & \ddots \end{array}\right)$$
Sorting the entries in the rows of $R_{h=2}$ we get:
$$R_{h=2,\text{ sorted rows}}=\left(\begin{array}{lllllllll} 1 & \text{} & \text{} & \text{} & \text{} & \text{} & \text{} & \text{} & \cdots\\ -2 & 2 & \text{} & \text{} & \text{} & \text{} & \text{} & \text{} \\ -1 & -1 & 2 & \text{} & \text{} & \text{} & \text{} & \text{} \\ -1 & -1 & 3 & 3 & \text{} & \text{} & \text{} & \text{} \\ -3 & -3 & 2 & 2 & 2 & \text{} & \text{} & \text{} \\ -4 & -2 & -2 & 2 & 2 & 4 & \text{} & \text{} \\ -5 & -5 & 2 & 2 & 2 & 2 & 2 & \text{} \\ 0 & 0 & 0 & 0 & 4 & 4 & 4 & 4 \\ \vdots & & & & & & & & \ddots \\ \end{array}\right)$$
Again:
$$\Lambda(n)\Lambda(n+h) = \left(\lim\limits_{s \rightarrow 1} \zeta(s)\sum\limits_{d|n} \frac{\mu(d)}{d^{(s-1)}}\right)\left( \lim\limits_{s \rightarrow 1} \zeta(s)\sum\limits_{d|n+h} \frac{\mu(d)}{d^{(s-1)}}\right)$$
where $n$ is unconventionally the column index.
For $h=4$, the convolution above is associated with the matrix of Dirichlet series:
$$M_{h=4}\left(\begin{array}{lllllllll} 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & \cdots\\ 2 & -2 & 2 & -2 & 2 & -2 & 2 & -2 \\ 2 & -1 & -1 & 2 & -1 & -1 & 2 & -1 \\ 3 & -1 & 3 & -1 & 3 & -1 & 3 & -1 \\ -3 & 2 & 2 & 2 & -3 & -3 & 2 & 2 \\ 4 & 2 & -2 & -4 & -2 & 2 & 4 & 2 \\ 2 & 2 & -5 & 2 & 2 & 2 & -5 & 2 \\ 4 & 0 & 4 & 0 & 4 & 0 & 4 & 0 \\ \vdots & & & & & & & & \ddots \end{array}\right)$$
$R_{h=4}$ is equal to $M_{h=4}$ for $\text{ row index } \geq \text{ column index}$:
$$R_{h=4}=\left(\begin{array}{lllllllll} 1 & \text{} & \text{} & \text{} & \text{} & \text{} & \text{} & \text{} & \cdots \\ 2 & -2 & \text{} & \text{} & \text{} & \text{} & \text{} & \text{} \\ 2 & -1 & -1 & \text{} & \text{} & \text{} & \text{} & \text{} \\ 3 & -1 & 3 & -1 & \text{} & \text{} & \text{} & \text{} \\ -3 & 2 & 2 & 2 & -3 & \text{} & \text{} & \text{} \\ 4 & 2 & -2 & -4 & -2 & 2 & \text{} & \text{} \\ 2 & 2 & -5 & 2 & 2 & 2 & -5 & \text{} \\ 4 & 0 & 4 & 0 & 4 & 0 & 4 & 0 \\ \vdots & & & & & & & & \ddots \end{array}\right)$$
Sorting the entries in the rows of $R_{h=4}$ we get:
$$R_{h=4,\text{ sorted rows}}=\left(\begin{array}{lllllllll} 1 & \text{} & \text{} & \text{} & \text{} & \text{} & \text{} & \text{} & \cdots \\ -2 & 2 & \text{} & \text{} & \text{} & \text{} & \text{} & \text{} \\ -1 & -1 & 2 & \text{} & \text{} & \text{} & \text{} & \text{} \\ -1 & -1 & 3 & 3 & \text{} & \text{} & \text{} & \text{} \\ -3 & -3 & 2 & 2 & 2 & \text{} & \text{} & \text{} \\ -4 & -2 & -2 & 2 & 2 & 4 & \text{} & \text{} \\ -5 & -5 & 2 & 2 & 2 & 2 & 2 & \text{} \\ 0 & 0 & 0 & 0 & 4 & 4 & 4 & 4 \\ \vdots & & & & & & & & \ddots \end{array}\right)$$
which looks very much as:
$$R_{h=2,\text{ sorted rows}}=\left(\begin{array}{lllllllll} 1 & \text{} & \text{} & \text{} & \text{} & \text{} & \text{} & \text{} & \cdots\\ -2 & 2 & \text{} & \text{} & \text{} & \text{} & \text{} & \text{} \\ -1 & -1 & 2 & \text{} & \text{} & \text{} & \text{} & \text{} \\ -1 & -1 & 3 & 3 & \text{} & \text{} & \text{} & \text{} \\ -3 & -3 & 2 & 2 & 2 & \text{} & \text{} & \text{} \\ -4 & -2 & -2 & 2 & 2 & 4 & \text{} & \text{} \\ -5 & -5 & 2 & 2 & 2 & 2 & 2 & \text{} \\ 0 & 0 & 0 & 0 & 4 & 4 & 4 & 4 \\ \vdots & & & & & & & & \ddots \\ \end{array}\right)$$
So proving that the asymptotic density of the twin primes $(h=2)$ is equal to asymptotic density of the cousin primes $(h=4)$, appears to be equivalent to proving:
$$R_{h=2,\text{ sorted rows}}=R_{h=4,\text{ sorted rows}}$$
Edit 28.2.2018:
In an attempt to verify $R_{h=2}=R_{h=4}$ I wrote this program that compares rows in $R_{h=2}$ and $R_{h=4}$. The shifts are found in the output from the program. Squares of primes is a problem in the comparison, but inspection shows that they so far match entrywise between $R_{h=2}$ and $R_{h=4}$.
(*start*)
nn = 230;
h = 2; (* h = 2 gives twin primes *)
hh = 4; (* hh = 4 gives cousin primes *)
(* TableForm[A=Table[Table[If[Mod[n,k]==0,1,0],{k,1,nn}],{n,1,nn}]];
TableForm[B=Table[Table[If[Mod[k,n]==0,MoebiusMu[n]*n,0],{k,1,nn}],{n,\
1,nn}]];
TableForm[T=(A.B)]; *)
Clear[a, n, d];
a[n_] := If[n < 1, 0,
Sum[d MoebiusMu@d, {d, Divisors[n]}]] (*Michael Somos,Jul 18 2011*)
TableForm[T = Table[Table[a[GCD[n, k]], {k, 1, nn}], {n, 1, nn}]];
nnn = 12;
TableForm[AA = Table[
Divisors[n];
Clear[a, b];
TableForm[c1 = Table[a = T[[Reverse[Divisors[n]], k]];
b = T[[Divisors[n], k + h]];
a*b, {k, 1, n}]];
TableForm[A = Reverse[Transpose[c1]]];
TableForm[c1 = Table[a = T[[Reverse[Divisors[n]], k]];
b = T[[Divisors[n], k + hh]];
a*b, {k, 1, 2*n - 1}]];
TableForm[B = Reverse[Transpose[c1]]];
TableForm[Table[A[[1]] -
B[[1]][[n + Range[Length[A[[1]]]]]], {n,
0, -n + Length[B[[1]]]}]];
Table[Min[Flatten[Position[Table[A[[k]] -
B[[k]][[n + Range[Length[A[[k]]]]]], {n,
0, -n + Length[B[[k]]]}],
Range[n]*0]]], {k, 1, Length[Divisors[n]]}], {n, 1, nnn}]];
Table[AA[[n]][[1]], {n, 1, nnn}]
AA = Table[Reverse[AA[[n]]], {n, 1, nnn}]
TableForm[
BB = Table[
Reverse[Accumulate[
Reverse[Table[If[Mod[n, k] == 0, 1, 0], {k, 1, nnn}]]]]*
Table[If[Mod[n, k] == 0, 1, 0], {k, 1, nnn}], {n, 1, nnn}]];
TableForm[
CC = Table[
Table[If[Mod[n, k] == 0, AA[[n, BB[[n, k]]]], 0], {k, 1,
nnn}], {n, 1, nnn}]];
MatrixForm[ReplaceAll[CC, Infinity -> 0]]
(*end*)
Shifts between rows in $R_{h=2}$ and $R_{h=4}$ $$=\left( \begin{array}{cccccccccccc} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 2 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 4 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 5 & 5 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 6 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 2 & 0 & 3 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 9 & 9 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 10 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 5 & 2 & 1 & 5 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 \end{array} \right)$$