The timing residuals of millisecond pulsar arrays, as measured by IPTA-grade radio telescopes, contain sufficient Shannon entropy density (≥256 bits per observation epoch from a 3-pulsar array) to serve as a cryptographically secure key derivation source, with residual unpredictability arising from gravitational wave background noise and spin noise that is irreducible by any known forecasting method.

A minimum 3-pulsar array observed with IPTA-grade radio telescopes (sensitivity ≥ 10 μJy, timing precision ≤ 100 ns RMS) produces timing residuals per observation epoch that contain ≥ 256 bits of Shannon entropy, where that entropy is sourced from gravitational wave background (GWB) noise and pulsar spin noise components that are provably irreducible by any polynomial-time forecasting algorithm, making the residuals suitable as input to a cryptographically secure pseudorandom number generator (CSPRNG) or key derivation function (KDF) per NIST SP 800-90B standards.

PHASE 1 (Weeks 1–4): Archival data entropy audit using existing IPTA Data Release 2 (DR2) and NANOGrav 15-yr datasets PHASE 2 (Weeks 5–12): Noise decomposition and entropy attribution (GWB vs. spin noise vs. instrumental) PHASE 3 (Weeks 13–20): Adversarial predictability testing using state-of-the-art GP forecasting PHASE 4 (Weeks 21–28): NIST SP 800-90B compliance testing on extracted bitstreams PHASE 5 (Weeks 29–36): Prototype KDF implementation and end-to-end security evaluation

# Pulsar Entropy Extraction Pipeline (PEEP) — Architecture Sketch

# === MODULE 1: DATA INGESTION ===
class PTADataLoader:
    def load_ipta_dr2(self, pulsars: list[str]) -> dict[str, TimingResiduals]:
        # Load .tim and .par files via libstempo/TEMPO2
        # Returns: {pulsar_name: TimingResiduals(epochs, residuals_ns, uncertainties_ns)}
        pass
    
    def select_top_pulsars(self, n=3, criterion='timing_rms') -> list[str]:
        # Rank by RMS timing residual, return top-n
        # Target: J0437-4715 (~30ns), J1909-3744 (~50ns), J1713+0747 (~70ns)
        pass

# === MODULE 2: NOISE DECOMPOSITION ===
class EnterpriseNoiseModel:
    def fit_noise_model(self, residuals: TimingResiduals) -> NoiseComponents:
        # Enterprise PTMCMC: white noise + red noise + DM + GWB
        # Returns variance fractions: {gwb: 0.xx, spin: 0.xx, dm: 0.xx, white: 0.xx}
        # Runtime: ~48 CPU-hours per pulsar
        pass
    
    def extract_stochastic_residuals(self, residuals, noise_model) -> np.ndarray:
        # Subtract deterministic timing model
        # Return: residuals_stochastic [n_epochs] in nanoseconds
        pass

# === MODULE 3: ENTROPY ESTIMATION ===
class EntropyEstimator:
    def quantize_residuals(self, residuals_ns: np.ndarray, 
                           n_bits: int = 8) -> np.ndarray:
        # Uniform quantization to n_bits integers
        # Scale: [-5*sigma, +5*sigma] -> [0, 2^n_bits - 1]
        return np.clip(
            ((residuals_ns - residuals_ns.mean()) / 
             (10 * residuals_ns.std()) * 256 + 128).astype(int), 
            0, 255
        ).astype(np.uint8)
    
    def estimate_shannon_entropy(self, samples: np.ndarray) -> float:
        # H = -sum(p_i * log2(p_i)) in bits
        counts = np.bincount(samples, minlength=256)
        probs = counts[counts > 0] / len(samples)
        return -np.sum(probs * np.log2(probs))
    
    def estimate_min_entropy(self, samples: np.ndarray) -> float:
        # H_min = -log2(max(p_i))
        counts = np.bincount(samples, minlength=256)
        p_max = counts.max() / len(samples)
        return -np.log2(p_max)
    
    def nist_800_90b_assessment(self, samples: np.ndarray) -> dict:
        # Call NIST reference implementation via subprocess
        # Returns: {min_entropy_estimate, iid_result, restart_test_result}
        import subprocess
        result = subprocess.run(
            ['./ea_non_iid', '-v', 'samples.bin'],
            capture_output=True
        )
        return parse_nist_output(result.stdout)
    
    def compute_epoch_entropy(self, psr_residuals: dict[str, np.ndarray],
                               epoch_idx: int) -> float:
        # Per-epoch entropy: sum individual - mutual information
        H_individual = []
        for psr, res in psr_residuals.items():
            window = res[max(0, epoch_idx-10):epoch_idx+10]  # 20-epoch window
            H_individual.append(self.estimate_min_entropy(
                self.quantize_residuals(window)))
        
        # Mutual information via k-NN estimator (sklearn)
        MI_pairs = compute_pairwise_MI(psr_residuals, epoch_idx, k=5)
        
        H_total = sum(H_individual) - sum(MI_pairs)
        return H_total  # in bits

# === MODULE 4: ADVERSARIAL TESTING ===
class AdversarialForecaster:
    def gp_forecast(self, residuals: np.ndarray, 
                    train_frac: float = 0.8) -> np.ndarray:
        from sklearn.gaussian_process import GaussianProcessRegressor
        from sklearn.gaussian_process.kernels import RBF, Matern, WhiteKernel
        kernel = RBF() + Matern(nu=1.5) + WhiteKernel()
        gpr = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=10)
        n_train = int(len(residuals) * train_frac)
        t = np.arange(len(residuals)).reshape(-1, 1)
        gpr.fit(t[:n_train], residuals[:n_train])
        predictions = gpr.predict(t[n_train:])
        return residuals[n_train:] - predictions  # prediction residuals
    
    def lstm_forecast(self, residuals: np.ndarray) -> np.ndarray:
        # PyTorch LSTM: 2 layers, 128 hidden units, seq_len=50
        # Returns: prediction residuals on held-out 20%
        pass
    
    def entropy_after_prediction(self, prediction_residuals: np.ndarray) -> float:
        return EntropyEstimator().estimate_min_entropy(
            EntropyEstimator().quantize_residuals(prediction_residuals))

# === MODULE 5: KDF PROTOTYPE ===
class PulsarKDF:
    def extract_bits(self, residuals_ns: np.ndarray, 
                     method: str = 'hash') -> bytes:
        if method == 'hash':
            # SHA3-256 of 64-byte blocks of quantized residuals
            import hashlib
            quantized = EntropyEstimator().quantize_residuals(residuals_ns)
            blocks = [quantized[i:i+64].tobytes() 
                      for i in range(0, len(quantized)-64, 64)]
            return b''.join(hashlib.sha3_256(b).digest() for b in blocks)
        elif method == 'lsb':
            # Least significant bits of quantized residuals
            quantized = EntropyEstimator().quantize_residuals(residuals_ns)
            return np.packbits(quantized & 0x01).tobytes()
    
    def derive_key(self, pulsar_residuals: dict[str, np.ndarray],
                   epoch_idx: int, key_length: int = 32) -> bytes:
        # HKDF (RFC 5869) with SHA3-256
        import hkdf
        ikm = b''.join(
            self.extract_bits(res[epoch_idx-20:epoch_idx], method='hash')
            for res in pulsar_residuals.values()
        )
        salt = f"pulsar-entropy-epoch-{epoch_idx}".encode()
        info = b"cryptographic-key-v1"
        return hkdf.hkdf_expand(
            hkdf.hkdf_extract(salt, ikm), info, key_length)

# === MODULE 6: MAIN PIPELINE ===
def run_validation_pipeline():
    loader = PTADataLoader()
    pulsars = loader.select_top_pulsars(n=3)
    data = loader.load_ipta_dr2(pulsars)
    
    noise_model = EnterpriseNoiseModel()
    stochastic_residuals = {}
    for psr in pulsars:
        nm = noise_model.fit_noise_model(data[psr])
        stochastic_residuals[psr] = noise_model.extract_stochastic_residuals(
            data[psr], nm)
    
    estimator = EntropyEstimator()
    epoch_entropies = []
    for epoch in range(50, len(list(stochastic_residuals.values())[0])):
        H = estimator.compute_epoch_entropy(stochastic_residuals, epoch)
        epoch_entropies.append(H)
    
    # Report: fraction of epochs with H >= 256 bits
    success_rate = np.mean(np.array(epoch_entropies) >= 256)
    print(f"Epochs with H_min >= 256 bits: {success_rate:.1%}")
    print(f"Median H_min per epoch: {np.median(epoch_entropies):.1f} bits")
    
    # NIST compliance
    all_residuals = np.concatenate(list(stochastic_residuals.values()))
    nist_result = estimator.nist_800_90b_assessment(
        estimator.quantize_residuals(all_residuals))
    print(f"NIST 800-90B min-entropy: {nist_result['min_entropy']:.3f} bits/sample")
    
    return epoch_entropies, nist_result

# === ABORT LOGIC ===
def check_abort_conditions(epoch_entropies_so_far: list) -> bool:
    if len(epoch_entropies_so_far) < 20:
        return False
    median_H = np.median(epoch_entropies_so_far)
    if median_H < 64:  # Less than 25% of threshold
        print("ABORT: Median entropy < 64 bits. Hypothesis likely false.")
        return True
    return False