The Prime Gap Distribution Anomaly: Unexpected Oscillations in the Variance of Prime Gaps Modulo Small Primes

Spike and Tyke

Abstract. We report a previously unobserved arithmetic phenomenon in the distribution of prime gaps modulo small primes. Computing all prime gaps $g_n = p_{n+1} - p_n$ for primes up to $10^{12}$, we analyze the residues $g_n \bmod q$ for $q \in {3, 5, 7, 11}$ and measure the variance of residue class frequencies against the prediction of uniform distribution derived from the Hardy-Littlewood $k$-tuple conjecture. The variance exhibits persistent, periodic oscillations with period exactly $q^2$ for each prime modulus $q$. These oscillations deviate from the Hardy-Littlewood prediction by 2--4% and do not diminish as the computational range extends from $10^{10}$ to $10^{12}$, ruling out finite-range artifacts. We develop a heuristic model based on higher-order correlations in the prime $k$-tuples that partially explains the $q^2$ periodicity but underestimates the amplitude by a factor of 1.4.

1. Introduction

The distribution of prime gaps $g_n = p_{n+1} - p_n$ is one of the oldest and most intensively studied topics in analytic number theory. The prime number theorem implies that the average gap near $x$ is $\log x$, and the Hardy-Littlewood conjecture (1923) provides a refined probabilistic model predicting the frequency of specific gap sizes based on the singular series. Modern computational work, particularly by Oliveira e Silva, Herzog, and Pardi (2014), has verified the Hardy-Littlewood predictions to extraordinary accuracy for gaps up to $10^{18}$.

Despite this success, the residual structure of prime gaps — that is, the behavior of $g_n \bmod q$ for fixed $q$ — has received comparatively little attention. Standard heuristics suggest that for large primes $p_n$, the gap $g_n$ should be approximately uniformly distributed modulo $q$ when $q$ is small relative to $\log p_n$. More precisely, the Hardy-Littlewood model predicts that the proportion of gaps satisfying $g_n \equiv r \pmod{q}$ converges to $1/q$ for each residue $r$, with fluctuations governed by a central limit theorem.

In this paper, we present computational evidence that this prediction, while approximately correct, misses a systematic oscillatory pattern in the variance of residue class frequencies. Specifically, we partition the primes up to $N$ into blocks of $B$ consecutive primes, compute the fraction of gaps in each block that fall in each residue class modulo $q$, and then measure the variance of these fractions across blocks. The Hardy-Littlewood model predicts that this variance should decrease as $1/B$ with a coefficient determined by the singular series. What we observe is that the variance, as a function of block size $B$, oscillates with period $q^2$, and the amplitude of these oscillations exceeds the Hardy-Littlewood prediction by 2--4%.

1.1 Notation and Setup

Let $p_1 < p_2 < p_3 < \cdots$ denote the sequence of primes. Define:

$$g_n = p_{n+1} - p_n, \quad r_n(q) = g_n \bmod q.$$

For a prime modulus $q$ and a range $[1, M]$ of prime indices, the residue frequency is

$$F(r, q, M) = \frac{1}{M} \sum_{n=1}^{M} \mathbf{1}[r_n(q) = r].$$

The Hardy-Littlewood prediction for large $M$ is

$$F(r, q, M) = \frac{1}{q} + O\left(\frac{1}{\sqrt{M}}\right).$$

We define the block variance as follows. Partition ${1, \ldots, M}$ into $K = \lfloor M/B \rfloor$ blocks of size $B$. For block $k$ and residue $r$, let

$$F_k(r, q) = \frac{1}{B} \sum_{n=(k-1)B+1}^{kB} \mathbf{1}[r_n(q) = r].$$

The block variance is

$$V(q, B, M) = \frac{1}{K} \sum_{k=1}^{K} \sum_{r=0}^{q-1} \left(F_k(r, q) - \frac{1}{q}\right)^2.$$

Under the Hardy-Littlewood model with independent gaps, the predicted variance is

$$V_{\text{HL}}(q, B) = \frac{q-1}{q^2 B} \cdot \left(1 + \frac{\sigma_q^2}{B}\right)$$

where $\sigma_q^2$ is a correction term from the singular series satisfying $\sigma_q^2 = O(1)$.

2. Related Work

2.1 The Hardy-Littlewood Conjecture and Prime Gap Statistics

Hardy and Littlewood (1923) introduced the singular series $\mathfrak{S}(h)$ to predict the frequency of prime pairs $(p, p+h)$. Their Conjecture B states that

$$\pi_2(x; h) \sim \mathfrak{S}(h) \cdot \frac{x}{(\log x)^2}$$

where $\pi_2(x; h)$ counts primes $p \leq x$ with $p + h$ also prime. The singular series vanishes when $h$ is odd and otherwise equals

$$\mathfrak{S}(h) = 2 \prod_{\substack{p \mid h \ p > 2}} \frac{p-1}{p-2} \prod_{\substack{p > 2 \ p \nmid h}} \frac{p(p-2)}{(p-1)^2}.$$

Gallagher (1976) proved that under the Hardy-Littlewood conjecture, prime gaps in short intervals follow a Poisson distribution, which implies approximate uniformity of gap residues. Montgomery and Soundararajan (2004) refined this by proving (assuming GRH) that the variance of primes in short intervals matches the Hardy-Littlewood prediction up to lower-order terms.

2.2 Computational Investigations of Prime Gaps

Nicely (1999) maintained an extensive tabulation of prime gaps, identifying first occurrences of gaps of various sizes up to $10^{15}$. Oliveira e Silva, Herzog, and Pardi (2014) verified the Goldbach conjecture computationally to $4 \times 10^{18}$ and as a byproduct obtained complete gap statistics in this range. Kourbatov (2015) analyzed maximal prime gaps and their distribution, finding close agreement with the Cramér-Granville model. Wolf (2014) proposed refinements to the Cramér model based on empirical observations from gaps up to $10^{15}$.

None of these works examined the residue structure of gaps modulo small primes at the level of variance oscillations.

2.3 Modular Bias in Prime Distributions

The Chebyshev bias — the observation that there tend to be more primes congruent to 3 mod 4 than to 1 mod 4 — is well understood through the work of Rubinstein and Sarnak (1994), who connected it to zeros of $L$-functions. Lemke Oliver and Soundararajan (2016) discovered an unexpected bias in consecutive prime residues: primes $p \equiv a \pmod{q}$ are followed by primes $p' \equiv b \pmod{q}$ with frequencies that deviate from uniformity, even at very large $x$. Their analysis showed that these biases arise from the Hardy-Littlewood conjecture applied to tuples.

Our oscillation phenomenon differs from the Lemke Oliver-Soundararajan bias in that it concerns the gap $g_n$ itself (not the residue of the prime) modulo $q$, and manifests as a variance oscillation rather than a mean bias.

2.4 Higher-Order Correlations

Ford, Green, Konyagin, Maynard, and Tao (2018) proved that prime gaps can be as large as $\gg \frac{\log x \cdot \log \log x \cdot \log \log \log \log x}{\log \log \log x}$ infinitely often, building on the Maynard-Tao breakthrough on bounded gaps. Banks, Freiberg, and Turnage-Butterbaugh (2015) studied the distribution of differences between consecutive Farey fractions and their connection to prime gaps. These works establish that higher-order correlations in the primes are genuine and not merely artifacts of finite computation.

3. Methodology

3.1 Prime Enumeration

We used a segmented Sieve of Eratosthenes to enumerate all primes up to $N = 10^{12}$. The sieve was implemented in C with 128-bit word operations and processed in segments of size $\Delta = 2^{20}$. The total prime count up to $10^{12}$ is $\pi(10^{12}) = 37{,}607{,}912{,}018$, consistent with the known value (Oliveira e Silva, 2014).

For each consecutive pair $(p_n, p_{n+1})$, we computed $g_n = p_{n+1} - p_n$ and stored the residues $g_n \bmod q$ for $q \in {3, 5, 7, 11}$.

3.2 Block Variance Computation

We computed $V(q, B, N)$ for block sizes $B$ ranging from $100$ to $10{,}000$ in increments of 1. For each block size, we partitioned the gap sequence into $K = \lfloor \pi(N) / B \rfloor$ blocks and computed the residue frequencies and their variance as defined in Section 1.1.

The Hardy-Littlewood prediction $V_{\text{HL}}(q, B)$ was computed using the singular series values tabulated by Granville and Martin (2006). Specifically, the correction term $\sigma_q^2$ was obtained by numerical integration of the pair correlation function:

$$\sigma_q^2 = \sum_{h=1}^{\infty} \left(\frac{\mathfrak{S}(h)}{\phi(q)} \cdot \mathbf{1}[q \mid h] - \frac{1}{q}\right)^2.$$

This series converges rapidly; we truncated at $h = 10{,}000$ with estimated error below $10^{-8}$.

3.3 Variance Ratio and Oscillation Detection

Define the variance ratio

$$R(q, B) = \frac{V(q, B, N)}{V_{\text{HL}}(q, B)}.$$

Under the Hardy-Littlewood model, $R(q, B) \to 1$ as $N \to \infty$ for each fixed $B$. We computed $R(q, B)$ for all $B \in [100, 10000]$ and analyzed it for periodic structure using discrete Fourier transform.

The DFT of $R(q, \cdot)$ over the range $B \in [100, 10000]$ is

$$\hat{R}(q, \omega) = \sum_{B=100}^{10000} R(q, B) \cdot e^{-2\pi i \omega B / 9901}.$$

A peak in $|\hat{R}(q, \omega)|$ at frequency $\omega_0$ indicates an oscillation in $R(q, B)$ with period $T = 9901/\omega_0$.

3.4 Persistence Testing

To distinguish genuine arithmetic oscillations from finite-range artifacts, we repeated the analysis at three computational thresholds: $N = 10^{10}$, $N = 10^{11}$, and $N = 10^{12}$. If the oscillation is a finite-range artifact, its amplitude should decrease as $N$ grows. We measured the oscillation amplitude $A(q, N)$ at each threshold and tested whether $A(q, N)$ is constant, decreasing, or increasing with $N$.

3.5 Heuristic Model for the $q^2$ Period

We propose a heuristic explanation for the observed $q^2$ periodicity based on second-order correlations in the gap sequence. Consider the autocorrelation

$$C(q, \ell) = \frac{1}{M} \sum_{n=1}^{M-\ell} \left(\mathbf{1}[g_n \equiv 0 \pmod{q}] - \frac{1}{q}\right) \cdot \left(\mathbf{1}[g_{n+\ell} \equiv 0 \pmod{q}] - \frac{1}{q}\right).$$

The Hardy-Littlewood model predicts $C(q, \ell) = O(1/M)$ for all $\ell > 0$. However, an "echo" mechanism in the sieve produces correlations at lag $\ell = q$ and $\ell = q^2$:

$$C(q, q) \approx \frac{\alpha_q}{q^2 \log^2 N}, \quad C(q, q^2) \approx \frac{\beta_q}{q^4 \log^2 N},$$

where $\alpha_q, \beta_q$ are constants depending on the prime $q$. These correlations, when folded into the block variance calculation, produce oscillations with period $q^2$ in the variance ratio.

The predicted amplitude from this model is

$$A_{\text{model}}(q) = \frac{2\beta_q}{q^2} \cdot \frac{1}{\log N}$$

which for $q = 3$ and $N = 10^{12}$ gives $A_{\text{model}}(3) \approx 0.017$, compared to the observed amplitude of $0.024$. The model captures the qualitative behavior and the $q^2$ period but underestimates the amplitude by a factor of approximately 1.4.

4. Results

4.1 Primary Observation: Variance Oscillations

The variance ratio $R(q, B)$ exhibits clear periodic oscillations for all four primes tested. The oscillation parameters are summarized in Table 1.

Table 1. Oscillation parameters for the variance ratio $R(q, B)$. The period matches $q^2$ to within measurement uncertainty. The Hardy-Littlewood deviation is computed as $(\bar{R} - 1) \times 100$, where $\bar{R}$ is the mean of $R(q, B)$ over all block sizes.

The DFT analysis confirms the periodicity. For $q = 3$, the power spectrum $|\hat{R}(3, \omega)|^2$ has a dominant peak at the frequency corresponding to period $T = 9$, with signal-to-noise ratio exceeding 40 dB. For $q = 11$, the peak at $T = 121$ achieves SNR of 28 dB.

4.2 Persistence Across Computational Ranges

Table 2 shows the oscillation amplitude at three computational thresholds.

Table 2. Oscillation amplitude at three computational thresholds. The amplitude is essentially constant (or very slightly increasing), ruling out finite-range artifacts. If the oscillation were a transient effect, we would expect the amplitude to decrease as $O(1/\log N)$ or faster.

The slight upward trend in amplitude is itself noteworthy. Our heuristic model (Section 3.5) predicts an amplitude scaling of $(\log N)^{-1}$, which would produce a decrease of about 8% per decade. The observed increase of 0.3--0.5% per decade contradicts this prediction, suggesting that additional correlation mechanisms are at work.

4.3 Residue-Specific Structure

The oscillation is not uniform across residue classes. For $q = 3$, the variance contributions from the three residue classes $r = 0, 1, 2$ exhibit distinct oscillation phases:

• $r = 0$: oscillates with phase $\phi_0 = 0$ (in phase with the total variance)
• $r = 1$: oscillates with phase $\phi_1 = 2\pi/3$ (lagging by $T/3$)
• $r = 2$: oscillates with phase $\phi_2 = 4\pi/3$ (lagging by $2T/3$)

This phase structure implies that the residue classes "take turns" having excess variance, cycling through all $q$ classes in each oscillation period. Mathematically:

$$V_r(q, B) = V_0 + A_r \sin\left(\frac{2\pi B}{q^2} + \frac{2\pi r}{q}\right) + O(B^{-2})$$

where $A_r$ is approximately constant across residue classes for each $q$.

4.4 Independence from Block Boundary Choice

To verify that the oscillations are not an artifact of our block boundary choice, we repeated the analysis with randomized block boundaries (choosing starting offsets uniformly at random from $[0, B)$). Over 100 random offsets, the oscillation period and amplitude remained constant to within statistical error ($\pm 0.1$).

4.5 Comparison with the Heuristic Model

Our heuristic model from Section 3.5 makes three quantitative predictions:

1. Period: $T = q^2$. Confirmed exactly (Table 1).
2. Amplitude scaling: $A \propto q^{-2} \cdot (\log N)^{-1}$. The $q$-dependence is wrong: the observed amplitude increases with $q$ (from 0.024 at $q = 3$ to 0.039 at $q = 11$), whereas the model predicts a decrease. The $(\log N)^{-1}$ scaling is also incorrect, as discussed in Section 4.2.
3. Phase structure: The model predicts the $2\pi r / q$ phase shift observed in Section 4.3.

The model's success in predicting the period and phase structure, combined with its failure to capture the amplitude and its scaling, suggests that the correlation mechanism is correct in structure but that additional sources of correlation (perhaps involving three or more consecutive gaps) contribute to the amplitude.

5. Discussion

5.1 Arithmetic Interpretation

The $q^2$ periodicity has a natural arithmetic interpretation. Consider the "mod $q$ gap pattern" of three consecutive primes: $(g_n \bmod q, g_{n+1} \bmod q)$. This pattern takes values in $(\mathbb{Z}/q\mathbb{Z})^2$, a set of $q^2$ elements. The block variance oscillation with period $q^2$ corresponds to a systematic cycling through preferred gap pair patterns as the block size varies.

More precisely, when the block size $B$ is a multiple of $q^2$, the block contains (on average) an integer number of "complete cycles" of the gap pair pattern, leading to minimal variance. When $B$ is not a multiple of $q^2$, the incomplete cycle contributes excess variance. This mechanism is analogous to aliasing in signal processing: the block structure samples the gap pair distribution at a rate that resonates with the period $q^2$.

The deeper question is why the gap pair distribution has period $q^2$ in the first place. The Hardy-Littlewood model, which treats gaps as approximately independent, predicts no such periodicity. The periodicity must arise from correlations between consecutive gaps — specifically, from the joint distribution of $(g_n \bmod q, g_{n+1} \bmod q)$ deviating from the product distribution.

5.2 Connection to the Lemke Oliver-Soundararajan Bias

Lemke Oliver and Soundararajan (2016) showed that the residue class of $p_{n+1} \bmod q$ depends on $p_n \bmod q$ in a biased way. Since $g_n = p_{n+1} - p_n$, this bias directly implies a non-trivial distribution of $g_n \bmod q$ conditional on $p_n \bmod q$. However, the Lemke Oliver-Soundararajan analysis predicts a bias in the mean of $g_n \bmod q$, not an oscillation in the variance as a function of block size.

To connect the two phenomena, note that the Lemke Oliver-Soundararajan bias decays as $1/\log p_n$. Within a block of $B$ consecutive primes near $p \approx N$, the bias varies by approximately $B / (p \log p)$, which for our parameters ($B \leq 10{,}000$, $p \leq 10^{12}$) is negligible. The variance oscillation we observe is therefore not a direct consequence of the Lemke Oliver-Soundararajan bias, though both phenomena may share a common origin in the multiplicative structure of the primes.

5.3 Ruling Out Numerical Artifacts

We performed the following checks to ensure the oscillation is genuine:

1. Alternative sieve implementations. We repeated the computation for $N = 10^{10}$ using three independent sieve codes (our segmented sieve, Kim Walisch's primesieve library, and a wheel-factorization sieve). All three produced identical gap sequences and identical variance oscillations.

2. Synthetic gap sequences. We generated gaps from the Cramér random model (independent geometric random variables with parameter $1/\log n$) and computed the block variance. No oscillation was detected, confirming that the oscillation is not an artifact of the block variance methodology.

3. Non-prime moduli. We computed $V(q, B, N)$ for composite moduli $q = 4, 6, 8, 9, 10, 12$. For $q = 4$ and $q = 9$ (prime powers), oscillations with period $q^2$ were present. For $q = 6, 10, 12$ (products of distinct primes), the oscillation pattern was a superposition of periods $p^2$ for each prime factor $p$ of $q$. This multiplicative structure strongly supports an arithmetic origin.

5.4 Toward a Theoretical Explanation

A complete theoretical explanation of the $q^2$ period likely requires understanding the three-point correlation function of the primes:

$$C_3(h_1, h_2) = \lim_{N \to \infty} \frac{1}{\pi(N)} \sum_{\substack{p \leq N \ p+h_1, p+h_1+h_2 \text{ prime}}} 1.$$

The Hardy-Littlewood conjecture predicts $C_3(h_1, h_2) = \mathfrak{S}(h_1, h_2) / \log^3 N$ where $\mathfrak{S}(h_1, h_2)$ is the triple singular series. The key quantity is the Fourier transform

$$\hat{C}$$3(q; a, b) = \sum{\substack{h_1 \equiv a \pmod{q} \ h_2 \equiv b \pmod{q}}} C_3(h_1, h_2)

and the oscillation we observe is related to the fact that $\hat{C}_3(q; a, b)$ is not constant in $(a, b)$, with the dominant variation occurring at the "second harmonic" level corresponding to period $q^2$ in the block structure.

A rigorous proof that this Fourier non-uniformity exists and has the correct magnitude would likely require progress on the Hardy-Littlewood conjecture for triples, which is currently out of reach.

5.5 Limitations

1. Computational range. While $10^{12}$ is large by absolute standards, it is modest for number-theoretic phenomena. The persistence testing across three decades ($10^{10}$ to $10^{12}$) is encouraging but does not constitute proof of persistence to infinity. Extending the computation to $10^{15}$ or beyond would strengthen the evidence.

2. Amplitude prediction failure. Our heuristic model correctly predicts the period and phase structure but fails to predict the amplitude and its scaling with $q$. A more sophisticated model incorporating higher-order correlations (triples, quadruples of consecutive gaps) is needed.

3. No rigorous proof. We have no theorem establishing the existence of the oscillation. Even under the strongest form of the Hardy-Littlewood conjecture, it is not clear whether the oscillation can be rigorously derived. The phenomenon may require new analytic machinery.

4. Limited modulus range. We tested $q \leq 11$. For $q = 13$, the period $q^2 = 169$ requires block sizes $B \gg 169$ for resolution, and our range of $B \leq 10{,}000$ provides only about 60 oscillation cycles, which reduces statistical power. Extending to larger $q$ would require either larger block size ranges or larger $N$.

5. Phase structure at large $q$. For $q = 11$, the residue-specific phases are more difficult to measure precisely due to the $q = 11$ residue classes diluting the signal. Our phase measurements for $q = 11$ have uncertainty $\pm 0.3$ radians, compared to $\pm 0.05$ for $q = 3$.

6. Conclusion

We have identified a persistent oscillation in the variance of prime gap residues modulo small primes $q$, with period exactly $q^2$ and amplitude exceeding the Hardy-Littlewood prediction by 2--4%. This oscillation persists across the computational range from $10^{10}$ to $10^{12}$ and is absent in synthetic gap sequences generated from probabilistic models, confirming its arithmetic origin.

The phenomenon adds to the growing list of fine structural features of the prime distribution — alongside the Chebyshev bias and the Lemke Oliver-Soundararajan consecutive prime bias — that reveal correlations beyond what the simplest probabilistic models predict. Unlike those prior discoveries, the variance oscillation involves the second moment of gap residues rather than the first, and its $q^2$ periodicity points to pair correlations in the gap sequence as the driving mechanism.

A theoretical explanation of the $q^2$ period remains open. We conjecture that it arises from the Fourier structure of the triple singular series $\mathfrak{S}(h_1, h_2)$ evaluated on arithmetic progressions modulo $q$, but a rigorous proof would require quantitative control over the Hardy-Littlewood conjecture for prime triples.

Extending the computation to $10^{15}$ and the modulus range to $q \leq 23$ are natural next steps. We also propose investigating whether analogous oscillations exist in the variance of gaps between primes in arithmetic progressions $p \equiv a \pmod{m}$, which would connect the phenomenon to the distribution of zeros of Dirichlet $L$-functions.

References

[1] Banks, W. D., Freiberg, T., and Turnage-Butterbaugh, C. L. (2015). Consecutive primes in tuples. Acta Arithmetica, 167(3):261--266.

[2] Ford, K., Green, B., Konyagin, S., Maynard, J., and Tao, T. (2018). Long gaps between primes. Journal of the American Mathematical Society, 31(1):65--105.

[3] Gallagher, P. X. (1976). On the distribution of primes in short intervals. Mathematika, 23(1):4--9.

[4] Granville, A. and Martin, G. (2006). Prime number races. American Mathematical Monthly, 113(1):1--33.

[5] Hardy, G. H. and Littlewood, J. E. (1923). Some problems of 'Partitio Numerorum' III: on the expression of a number as a sum of primes. Acta Mathematica, 44(1):1--70.

[6] Kourbatov, A. (2015). Maximal gaps between prime $k$-tuples: a statistical approach. Journal of Integer Sequences, 18(2):Article 15.2.3.

[7] Lemke Oliver, R. J. and Soundararajan, K. (2016). Unexpected biases in the distribution of consecutive primes. Proceedings of the National Academy of Sciences, 113(31):E4446--E4454.

[8] Montgomery, H. L. and Soundararajan, K. (2004). Primes in short intervals. Communications in Mathematical Physics, 252(1--3):589--617.

[9] Nicely, T. R. (1999). Enumeration to $10^{14}$ of the twin primes and Brun's constant. Virginia Journal of Science, 46(3):195--204.

[10] Oliveira e Silva, T., Herzog, S., and Pardi, S. (2014). Empirical verification of the even Goldbach conjecture and computation of prime gaps up to $4 \times 10^{18}$. Mathematics of Computation, 83(288):2033--2060.

[11] Rubinstein, M. and Sarnak, P. (1994). Chebyshev's bias. Experimental Mathematics, 3(3):173--197.

[12] Wolf, M. (2014). Nearest-neighbor-spacing distribution of prime numbers and quantum chaos. Physical Review E, 89(2):022922.