After removal of all deterministic components via high-fidelity pulsar timing models, residual pulsar timing noise constitutes a physically irreducible stochastic process — providing a universally accessible, publicly verifiable source of certified true randomness for cryptographic applications requiring only a radio telescope and an open pulsar ephemeris.

Residual timing noise in millisecond pulsars (MSPs), after subtraction of all deterministic signal components (spin-down, proper motion, parallax, Shapiro delay, dispersion measure variations, binary orbital parameters) using state-of-the-art timing models (TEMPO2/PINT), constitutes a statistically irreducible stochastic process whose entropy is sourced from quantum-mechanical and astrophysical processes (turbulent interstellar medium, magnetospheric plasma fluctuations, superfluid vortex dynamics) that are computationally infeasible to predict or reproduce. Specifically: (1) residuals from ≥5 MSPs pass all NIST SP 800-22 and TestU01 randomness test suites at p>0.01 after whitening; (2) no adversary with access to all public ephemerides can predict future residuals better than chance (Shannon entropy H ≥ 0.99 bits/bit); (3) the system can deliver ≥1 bit/s of certified randomness per pulsar observable with a standard radio telescope (dish ≥25m diameter, receiver sensitivity SEFD ≤2000 Jy).

Phase 1 — Data Acquisition (Days 1–90): Obtain archival ToA datasets from NANOGrav 15-year data release, EPTA DR2, and PPTA DR3 for 10 MSPs. Select 5 MSPs with best timing precision (RMS < 500 ns): PSR J0437−4715, PSR J1909−3744, PSR J1713+0747, PSR J0030+0451, PSR J1744−1134.

Phase 2 — Deterministic Subtraction (Days 30–120): Apply TEMPO2 and PINT timing models with full parameter sets. Perform F-test to confirm no additional deterministic parameters improve fit at p < 0.01. Compute whitened residuals using Cholesky decomposition of the noise covariance matrix.

Phase 3 — Randomness Characterization (Days 90–180): Convert whitened residuals to bit strings via von Neumann extractor and hash-based extractor (SHA-3). Apply NIST SP 800-22 (15 tests), TestU01 BigCrush (106 tests), and Dieharder (18 tests) suites. Compute min-entropy via compression-based estimator (NIST SP 800-90B).

Phase 4 — Adversarial Prediction Challenge (Days 120–210): Train LSTM, Transformer, and Gaussian Process models on 80% of residual data; test on held-out 20%. Compute prediction accuracy, R², and compare to null model (mean prediction).

Phase 5 — Multi-Observatory Cross-Validation (Days 180–270): Compare simultaneous residuals from NANOGrav and EPTA for overlapping pulsars. Compute cross-correlations and test for shared structure.

Phase 6 — Throughput Benchmarking (Days 240–300): Measure bits/second output rate for a 25m dish observing PSR J0437−4715 at 1.4 GHz. Estimate practical randomness generation rate.

# Pulsar Timing Residual TRNG — Implementation Architecture
# ============================================================

# LAYER 1: DATA INGESTION
class PulsarDataIngester:
    def __init__(self, pulsar_id: str, pta_source: str = "NANOGrav15"):
        self.pulsar_id = pulsar_id  # e.g., "J1713+0747"
        self.pta_source = pta_source
        
    def load_toa_file(self, path: str) -> np.ndarray:
        """Load Time-of-Arrival data from .tim file format"""
        # Returns: array of (MJD, ToA_uncertainty_us, telescope_code)
        return tempo2_reader.parse_tim(path)
    
    def load_ephemeris(self, path: str) -> dict:
        """Load .par timing model file"""
        # Returns: dict of {parameter_name: (value, uncertainty)}
        return tempo2_reader.parse_par(path)

# LAYER 2: DETERMINISTIC SUBTRACTION
class TimingModelSubtractor:
    def __init__(self, ephemeris: dict, solar_system_ephem: str = "DE440"):
        self.model = PINT_TimingModel(ephemeris, solar_system_ephem)
        
    def compute_residuals(self, toas: np.ndarray) -> np.ndarray:
        """
        Subtract full deterministic model from ToAs.
        Returns raw residuals in microseconds.
        """
        # Apply: barycentric correction, Shapiro delay, DM correction,
        # spin-down model, binary orbital model (if applicable)
        return self.model.residuals(toas)
    
    def f_test_model_completeness(self, residuals: np.ndarray, 
                                   candidate_params: list) -> dict:
        """
        F-test to verify no additional deterministic parameters
        improve fit at p < 0.01.
        Returns: {param: p_value} for all candidates
        """
        results = {}
        baseline_chi2 = np.sum(residuals**2 / residuals_uncertainty**2)
        for param in candidate_params:
            extended_model = self.model.add_parameter(param)
            new_residuals = extended_model.residuals(toas)
            new_chi2 = np.sum(new_residuals**2 / residuals_uncertainty**2)
            delta_chi2 = baseline_chi2 - new_chi2
            p_value = chi2.sf(delta_chi2, df=1)
            results[param] = p_value
        return results

# LAYER 3: NOISE CHARACTERIZATION (Bayesian)
class BayesianNoiseCharacterizer:
    def __init__(self, pulsar_id: str):
        self.pta = enterprise.PTA([pulsar_id])
        
    def fit_noise_model(self, residuals: np.ndarray, 
                        toas: np.ndarray) -> dict:
        """
        Fit white noise (EFAC, EQUAD) and red noise (power law)
        using MultiNest sampler.
        Returns: posterior samples for all noise hyperparameters
        """
        sampler = PTMCMCSampler(self.pta, Niter=500000, 
                                 resume=False)
        chain = sampler.sample()
        return {
            'EFAC': chain['efac_posterior'],
            'EQUAD': chain['equad_posterior'],
            'red_noise_amplitude': chain['rn_amp_posterior'],
            'red_noise_spectral_index': chain['rn_gamma_posterior'],
            'DM_noise_amplitude': chain['dm_amp_posterior']
        }
    
    def compute_noise_covariance(self, posteriors: dict, 
                                  toas: np.ndarray) -> np.ndarray:
        """
        Construct full noise covariance matrix C from posterior medians.
        Shape: (N_toa, N_toa)
        """
        C_white = np.diag((posteriors['EFAC'] * toa_errors)**2 + 
                           posteriors['EQUAD']**2)
        C_red = red_noise_covariance(toas, 
                                      posteriors['red_noise_amplitude'],
                                      posteriors['red_noise_spectral_index'])
        C_dm = dm_noise_covariance(toas, posteriors['DM_noise_amplitude'])
        return C_white + C_red + C_dm

# LAYER 4: WHITENING
class ResidualWhitener:
    def whiten(self, residuals: np.ndarray, 
               C: np.ndarray) -> np.ndarray:
        """
        Cholesky whitening: r_w = L^{-1} r where C = L L^T
        Returns: whitened residuals (should be iid N(0,1))
        """
        L = np.linalg.cholesky(C)
        r_whitened = np.linalg.solve(L, residuals)
        return r_whitened
    
    def verify_whitening(self, r_w: np.ndarray) -> dict:
        """
        Verify whitened residuals have flat PSD and unit variance.
        """
        psd_freqs, psd_power = lombscargle(toas, r_w)
        max_peak_sigma = (max(psd_power) - mean(psd_power)) / std(psd_power)
        return {
            'variance': np.var(r_w),  # Should be ~1.0
            'max_psd_peak_sigma': max_peak_sigma,  # Should be < 3.0
            'normality_p': scipy.stats.normaltest(r_w).pvalue  # Should be > 0.05
        }

# LAYER 5: RANDOMNESS EXTRACTION
class RandomnessExtractor:
    def von_neumann_extract(self, bits: np.ndarray) -> np.ndarray:
        """
        Von Neumann extractor on consecutive bit pairs.
        (0,1) -> 0; (1,0) -> 1; (0,0),(1,1) -> discard
        Efficiency: ~50% for fair coins
        """
        output = []
        for i in range(0, len(bits)-1, 2):
            if bits[i] == 0 and bits[i+1] == 1:
                output.append(0)
            elif bits[i] == 1 and bits[i+1] == 0:
                output.append(1)
        return np.array(output)
    
    def hash_extract(self, residuals: np.ndarray, 
                     block_size: int = 256) -> bytes:
        """
        SHA3-256 hash extractor for stronger entropy concentration.
        Input: 256 quantized residual values (8-bit each = 2048 bits)
        Output: 256 bits of extracted randomness
        """
        output_bits = []
        for i in range(0, len(residuals) - block_size, block_size):
            block = residuals[i:i+block_size]
            quantized = np.int8(block * 127).tobytes()
            hash_output = hashlib.sha3_256(quantized).digest()
            output_bits.extend(hash_output)
        return bytes(output_bits)

# LAYER 6: STATISTICAL TESTING
class RandomnessTestSuite:
    def run_nist_sp800_22(self, bit_sequence: bytes) -> dict:
        """Run all 15 NIST SP 800-22 tests. Returns p-values."""
        return nist_tests.run_all(bit_sequence, n_bits=10**6)
    
    def run_testu01_bigcrush(self, bit_sequence: bytes) -> dict:
        """Run TestU01 BigCrush (106 tests). Returns pass/fail."""
        return testu01.bigcrush(bit_sequence)
    
    def compute_min_entropy(self, bit_sequence: bytes) -> float:
        """
        NIST SP 800-90B min-entropy estimation.
        Returns: H_min in bits/bit (range 0-1)
        """
        estimators = [
            nist_90b.compression_estimator(bit_sequence),
            nist_90b.collision_estimator(bit_sequence),
            nist_90b.markov_estimator(bit_sequence),
            nist_90b.most_common_value_estimator(bit_sequence)
        ]
        return min(estimators)  # Conservative: take minimum

# LAYER 7: ADVERSARIAL PREDICTION CHALLENGE
class AdversarialPredictor:
    def train_lstm(self, residuals: np.ndarray, 
                   train_frac: float = 0.8) -> dict:
        """3-layer LSTM, 256 hidden units, trained on 80% of data."""
        model = LSTM(layers=3, hidden=256, dropout=0.2)
        train_size = int(len(residuals) * train_frac)
        model.fit(residuals[:train_size], epochs=100, 
                  batch_size=32, patience=10)
        predictions = model.predict(residuals[train_size:])
        return self._evaluate(residuals[train_size:], predictions)
    
    def _evaluate(self, true: np.ndarray, 
                  pred: np.ndarray) -> dict:
        sign_accuracy = np.mean(np.sign(true) == np.sign(pred))
        r2 = 1 - np.sum((true-pred)**2) / np.sum((true-np.mean(true))**2)
        return {'sign_accuracy': sign_accuracy, 'r2': r2,
                'mse': np.mean((true-pred)**2)}

# LAYER 8: MAIN PIPELINE
def run_pulsar_trng_validation(pulsar_ids: list) -> dict:
    results = {}
    for pid in pulsar_ids:
        ingester = PulsarDataIngester(pid)
        toas = ingester.load_toa_file(f"data/{pid}.tim")
        ephem = ingester.load_ephemeris(f"data/{pid}.par")
        
        subtractor = TimingModelSubtractor(ephem)
        residuals = subtractor.compute_residuals(