I expand my comment into a full answer, providing more detail.
Fix $r \ge 1$ and a modulus $q\ge 1$. Let us define more generally, for any given $\mathbf{a}=(a_1,\ldots,a_r) \in (\mathbb{Z}/q\mathbb{Z})^r$,
$$\pi(x; q,\mathbf{a}) = \# \{ p_n \le x: p_{n+i-1} \equiv a_i \bmod q, (i=1,\ldots,r)\}.$$
That is, we are counting tuples $(p_n,p_{n+1},\ldots,p_{n+r-1})$ of $r$ consecutive primes with $p_n \le x$, such that $p_{n+i-1} \bmod q$ is prescribed.
For $r=1$, the limit of $\pi(x;q,a)/\pi(x)$ as $x \to \infty$ is described by Dirichlet's Theorem on primes in arithmetic progression. It tells us that $\lim_{x \to \infty} \pi(x;q,a)/\pi(x)$ is equal to $1/\phi(q)$ if $\gcd(a,q)=1$ and is equal to $0$ otherwise.
In a recent work of Lemke Oliver and Soundararajan, ''Unexpected biases in the distribution of consecutive primes'' (Proc. Natl. Acad. Sci. USA 113 (2016), no. 31, E4446–E4454), the authors used a modified Hardy-Littlewood Conjecture to understand $\pi(x;q,\mathbf{a})$, down to lower order terms, when $k\ge 2$. Obviously, $\pi(x;q,\mathbf{a})$ is $O(1)$ unless each $a_i$ is coprime to $q$, which we shall assume from now on.
Their Main Conjecture states that
$$\pi(x;q,\mathbf{a})= \frac{\mathrm{li}(x)}{\phi(q)^r}\left( 1+c_1(q;\mathbf{a})\frac{\log \log x}{\log x} + \frac{c_2(q;\mathbf{a})}{\log x} + O\left( (\log x)^{-7/4}\right)\right)$$
for certain arithmetic constants $c_1(q;\mathbf{a})$ and $c_2(q;\mathbf{a})$ which, in a sense, are the most interesting part of their work. Here $\mathrm{li}(x)$ is the logarithm integral. In particular,
$$\frac{\pi(x;q,\mathbf{a})}{\pi(x)}$$
is equidistributed across 'reduced' (i.e. coprime to $q$) vectors $\mathbf{a}$. This a very natural conjecture. As the authors write,
Any model based on the randomness of the primes would suggest strongly that every permissible pattern of $r$ consecutive primes appears roughly equally often.
Now let's turn to the gaps you study. One can express your $P(k,m,a)$ as
$$P(k, m,a) = \lim_{x \to \infty} \sum_{\substack{\mathbf{a} \in (\mathbb{Z}/2m\mathbb{Z})^{k+1} \\ a_{k+1}-a_1 \equiv 2a \bmod {2m}\\ \gcd(a_i,2m)=1} } \frac{\pi(x;2m,\mathbf{a})}{\pi(x)},$$
which, by the above conjecture, is simply
$$P(k, m,a) = \frac{1}{\phi(2m)^{k+1}} \sum_{\substack{\mathbf{a} \in (\mathbb{Z}/2m\mathbb{Z})^{k+1} \\ a_{k+1}-a_1 \equiv 2a \bmod {2m}\\ \gcd(a_i,2m)=1} }1.$$
Note that this expression for $P(k,m,a)$ makes no reference to primes. Moreover, the main condition, $a_{k+1}-a_1 \equiv 2a \bmod{2m}$, does not depend on $a_2,\ldots,a_k$, so that $P(k,m,a)$ may be simplified as
$$P(k,m,a) = P(m,a)$$
where
$$P(m,a):=\frac{1}{\phi(2m)^{2}} \sum_{\substack{\mathbf{a} \in (\mathbb{Z}/2m\mathbb{Z})^{2} \\ a_{2}-a_1 \equiv 2a \bmod {2m}\\ \gcd(a_i,2m)=1} }1$$
is independent of $k$ (explaining your observation).
It remains to evaluate $P(m,a)$ which is a simple quantity as it does not involve primes whatsoever and instead only involves 'local' number theory. It is an exercise in elementary number theory.
First, the Chinese Remainder Theorem explains the multiplicative property that you found (your equation (38)), since any solution to the equation $a_2 - a_1 \equiv 2a \bmod{2m}$ satisfying $(a_1,2m)=(a_2,2m)=1$ can be reduced to a solution to $a_{2,p_i^{\alpha_i}}-a_{1,p_i^{\alpha_i}} \equiv 2a \bmod {2p_i^{\alpha_i}}$ for each prime power $p_i^{\alpha_i}$ dividing $m$, where $(a_{1,p_i^{\alpha_i}},2p_i^{\alpha_i})=(a_{2,p_i^{\alpha_i}},2p_i^{\alpha_i})=1$ (and solutions modulo $2p_i^{\alpha_i}$ can be lifted to a solution modulo $2m$).
We've reduced matters to prime powers, so we are just counting $(a_1,a_2) \in (\mathbb{Z}/2p^{k}\mathbb{Z})^2$ with $(a_1a_2, 2p^k)=1$ and $a_2-a_1 \equiv 2a \bmod {2p^k}$. Let us suppose first that $p=2$. In this case, we count $x \bmod {2^{k+1}}$ with $x$ and $x+2a$ being odd. This is equivalent simply to counting odd $x \bmod 2^{k+1}$, of which there are $2^k$ in total, hence $P(2^{k},a) = 2^k/\phi(2^{k+1})^2=1/2^k$, as you've observed in your equation (35).
If $p>2$, we count $x \bmod {2p^k}$ with $x$ and $x+2a$ coprime to $2p$. Again using the Chinese remainder theorem (seeing that $x$ must be $1 \bmod 2$ and that $\phi(2)=1$), we may simply count $x \bmod p^k$ with $\gcd(x(x+2a),p)=1$, and then divide by $\phi(p^k)^2 = p^{2k-2}(p-1)^2$. Such $x$ are easy to parameterize, since the condition $\gcd(x(x+2a),p)=1$ simply says that the least significant digit of $x$ in base-$p$ is not $0$ or $-2a \bmod p$. More precisely, these are $x$ with $x \not\equiv 0,-2a \bmod p$. If $p \nmid a$, there are $p^{k-1}(p-2)$ such $x$s, and otherwise there are $p^{k-1}(p-1)$. This explains your formula for $P(p^{k},a)$ in equations (36) and (37).
Although current technology cannot yet establish the equidistribution result mentioned above, it is very easy to prove a local analogue of your result. Namely, fix an integer $M$ (say, a primorial) and consider the integers that are coprime to $M$. You can list them as $p'_i$ and define an analogue of your $d_n(k)$ and of your $\pi(x;k,d)$ and $P(x;k,m,a)$ for these numbers, at least when your $2m$ divide $M$, and you'll be able to rigorously prove similar formulas for the analogues of $P(m,a)$ in this setting.