Planck CMB: from measured pixels to physical hypotheses

Author: Biswajit Jana · CMB Anomaly Lab

Make the blue-orange sky map from real Planck measurements, understand the equations, run Gaussian simulations, and test a primordial hypothesis with CAMB.

Project website · Source repository · Research report

Run from top to bottom in a Python CPU runtime. The first run downloads about 2 GB of official SMICA data; disk and RAM must accommodate it. No Google Drive access is required. Runtime and download speed vary. Saved outputs were actually executed in a local Python kernel and are included for inspection; opening in Colab creates your own execution environment.

Measurements: ESA / Planck Collaboration, distributed by NASA/IPAC IRSA. Code, explanation and project authorship: Biswajit Jana, with disclosed AI coding/drafting assistance. These are processed scientific products, not raw detector timelines. Statistical conclusions are conditional on the specified idealized models.

1. Set up the reproducible research code

The analysis code is pinned to a specific public commit. A local execution can set CMB_WORKSPACE to the existing checkout; fresh Colab sessions fetch the pinned revision.

import os, sys, subprocess, platform, json, importlib.metadata
from pathlib import Path
os.environ['OMP_NUM_THREADS']='2'
os.environ['OPENBLAS_NUM_THREADS']='2'
REVISION='47282e0e5f2b1c63f2bdbae87c213a453d133a8e'
WORK=Path(os.environ.get('CMB_WORKSPACE', str(Path.cwd()/'cmb-anomaly-lab'))).resolve()
if not (WORK/'src/cmb_lab').exists():
    subprocess.run(['git','clone','--quiet','https://github.com/Biswajit1999/cmb-anomaly-lab.git',str(WORK)],check=True)
    subprocess.run(['git','-C',str(WORK),'checkout','--quiet',REVISION],check=True)
actual=subprocess.check_output(['git','-C',str(WORK),'rev-parse','HEAD'],text=True).strip()
if actual!=REVISION:
    raise RuntimeError('Use a clean checkout of the pinned analysis revision or a fresh workspace.')
subprocess.run([sys.executable,'-m','pip','install','-q','-e',str(WORK)+'[maps,physics]'],check=True,capture_output=True)
print('Author: Biswajit Jana')
print('Python:',platform.python_version(),'| Analysis revision:',actual)
print('Kernel: Google Colab' if 'google.colab' in sys.modules else 'Local Python / Jupyter')

Author: Biswajit Jana
Python: 3.10.6 | Analysis revision: 47282e0e5f2b1c63f2bdbae87c213a453d133a8e
Local Python / Jupyter

2. Load the scientific libraries

SciencePlots controls figure styling. HEALPix describes equal-area pixels on the sphere. FITS stores the measured map and its metadata.

import numpy as np
import matplotlib.pyplot as plt
from IPython.display import display, Image, Markdown
from astropy.io import fits
import healpy as hp
from cmb_lab.data import fetch, sha256, load_spectra
from cmb_lab.statistics import dl_to_cl, s_kernel, s_stat, correlation_matrix, simulate, lower_tail
from cmb_lab.plots import style
style()
print({p:importlib.metadata.version(p) for p in ['numpy','scipy','healpy','astropy','SciencePlots','camb']})

{'numpy': '2.2.6', 'scipy': '1.15.3', 'healpy': '1.20.1', 'astropy': '6.1.7', 'SciencePlots': '2.2.2', 'camb': '2.0.5'}

3. Retrieve the actual Planck measurements

The downloader verifies cached files, resumes interrupted transfers and records provenance. A SHA-256 hash establishes file consistency, not a publisher signature.

manifest=fetch(WORK, include_map=True)
for name,record in manifest.items():
    print(name, '| bytes:',record['bytes'],'| SHA-256:',record['sha256'])
    print(record['url'])

Verified COM_PowerSpect_CMB-TT-full_R3.01.txt

Verified COM_PowerSpect_CMB-base-plikHM-TTTEEE-lowl-lowE-lensing-minimum-theory_R3.01.txt

Verified COM_CMB_IQU-smica_2048_R3.00_full.fits

observed | bytes: 170547 | SHA-256: ccf3113604020536f6f13ccf51680a7316ad0f32da558eee7f625e613bdd5522
https://irsa.ipac.caltech.edu/data/Planck/release_3/ancillary-data/cosmoparams/COM_PowerSpect_CMB-TT-full_R3.01.txt
theory | bytes: 205647 | SHA-256: 50eecd4422e22026070bdacf9fa8a7e7ef3badf4fa50fffea040df669715cf51
https://irsa.ipac.caltech.edu/data/Planck/release_3/ancillary-data/cosmoparams/COM_PowerSpect_CMB-base-plikHM-TTTEEE-lowl-lowE-lensing-minimum-theory_R3.01.txt
smica | bytes: 2013312960 | SHA-256: 60952c645eb33d151905ddf5837477e15ca02a3261feca4ceca3c3feece4f9ac
https://irsa.ipac.caltech.edu/data/Planck/release_3/all-sky-maps/maps/component-maps/cmb/COM_CMB_IQU-smica_2048_R3.00_full.fits

4. Inspect the FITS header

Before plotting, check the units and pixel ordering. K_CMB becomes microkelvin by multiplying by $10^6$. Never interpret NESTED pixels as RING pixels.

map_record=manifest['smica']
with fits.open(WORK/map_record['path'],memmap=True) as h:
    print({key:h[1].header[key] for key in ['NSIDE','ORDERING','COORDSYS']})
    print('Temperature column:',h[1].columns['I_STOKES'].unit)
    print('Columns:',h[1].columns.names)

{'NSIDE': 2048, 'ORDERING': 'NESTED', 'COORDSYS': 'GALACTIC'}
Temperature column: K_CMB
Columns: ['I_STOKES', 'Q_STOKES', 'U_STOKES', 'TMASK', 'PMASK', 'I_STOKES_INP', 'Q_STOKES_INP', 'U_STOKES_INP', 'TMASKINP', 'PMASKINP']

5. Define the measured-map renderer

The blue-cream-orange palette is inspired by Planck outreach imagery. It is our rendering of PR3 values, not a copied ESA image. The full-sky display includes pixels excluded from the confidence-masked statistical analysis. The renderer below is included explicitly so its operations are inspectable.

"""Render measured PR3 SMICA in an ESA-inspired blue/cream/orange palette.

Author: Biswajit Jana. The palette is our visualization, not an ESA image copy.
Full-sky display does not imply foreground-clean confidence in every pixel.
"""
from pathlib import Path
import json
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
import numpy as np
import healpy as hp
from astropy.io import fits
from cmb_lab.data import sha256
from cmb_lab.plots import style


def read_display(root, nside=512):
    record=json.loads((root/'data/manifest.json').read_text())['smica']
    path=root/record['path']
    if sha256(path)!=record['sha256']:
        raise ValueError('Measured FITS file checksum mismatch')
    with fits.open(path,memmap=True) as h:
        if h[1].columns['I_STOKES'].unit!='K_CMB':
            raise ValueError('Unexpected temperature unit')
        order=h[1].header['ORDERING']
        sky=hp.ud_grade(h[1].data['I_STOKES'].ravel(),nside,order_in=order,order_out='RING')*1e6
        fraction=hp.ud_grade(h[1].data['TMASK'].ravel().astype(float),nside,order_in=order,order_out='RING')
    return sky,fraction,record


def render(root, sky, fraction, record):
    style()
    cmap=LinearSegmentedColormap.from_list('cmb_blue_orange',
        ['#083478','#147cc1','#69c8e8','#f4f3ca','#ffd365','#ed972f','#a74714'],N=256)
    out=root/'results/maps'; out.mkdir(parents=True,exist_ok=True)
    for masked,name,title in [(False,'planck_fullsky','Planck 2018 | measured SMICA temperature'),
                              (True,'planck_blue_orange_masked','Planck 2018 | confidence-masked temperature')]:
        field=hp.ma(sky)
        field.mask=(fraction<.95) if masked else (~np.isfinite(sky))
        fig=plt.figure(figsize=(12,6.2))
        hp.mollview(field,fig=fig.number,min=-300,max=300,cmap=cmap,badcolor='#d7dfe5',
                    title=title,unit=r'$\Delta T$ [$\mu$K]',coord='G',xsize=1600)
        note='Grey: temperature confidence mask.' if masked else 'Full-sky display; confidence exclusions are shown in the companion map.'
        fig.text(.5,-.075,'Biswajit Jana | CMB Anomaly Lab\nESA / Planck Collaboration measurements | Galactic coordinates\n'+note,
                 ha='center',fontsize=9)
        fig.savefig(out/f'{name}.png',dpi=180,bbox_inches='tight')
        fig.savefig(out/f'{name}.pdf',bbox_inches='tight')
        plt.close(fig)
    provenance={'author':'Biswajit Jana','source_sha256':record['sha256'],'source_url':record['url'],
        'display_nside':hp.npix2nside(len(sky)),'units':'microkelvin','limits':[-300,300],
        'coordinate_system':'Galactic','palette':'custom blue-cream-orange, inspired by Planck outreach images',
        'processing':'native pixel averaging for visualization only; not the map-statistics estimator',
        'full_sky_is_unmasked_display':True,'confidence_mask_min_fraction':.95}
    (out/'blue_orange_provenance.json').write_text(json.dumps(provenance,indent=2)+'\n')
    return out

6. Read and reduce pixels for display

Pixel averaging to NSIDE=512 is a display operation. The quantitative map experiment later extracts harmonics at native resolution before smoothing.

sky, mask_fraction, verified_record=read_display(WORK)
print('Display pixels:',len(sky),'| NSIDE:',hp.npix2nside(len(sky)))
print('Finite measured temperatures:',bool(np.isfinite(sky).all()))
print('Display confidence fraction:',float(np.mean(mask_fraction>=.95)))

Display pixels: 3145728 | NSIDE: 512
Finite measured temperatures: True
Display confidence fraction: 0.8399804433186849

7. Make the full-sky blue-orange map

The Mollweide projection shows Galactic coordinates. The colour scale is fixed at ±300 microkelvin.

map_outputs=render(WORK,sky,mask_fraction,verified_record)
display(Image(filename=str(map_outputs/'planck_fullsky.png')))

/home/biswajit/.venvs/cmb-anomaly-lab/lib/python3.10/site-packages/healpy/visufunc.py:225: UserWarning: Ignoring specified arguments in this call because figure with num: 1 already exists
  f = pylab.figure(fig, figsize=(8.5, 5.4))

'created' timestamp seems very low; regarding as unix timestamp

'modified' timestamp seems very low; regarding as unix timestamp

'created' timestamp seems very low; regarding as unix timestamp

'modified' timestamp seems very low; regarding as unix timestamp

'created' timestamp seems very low; regarding as unix timestamp

'modified' timestamp seems very low; regarding as unix timestamp

'created' timestamp seems very low; regarding as unix timestamp

'modified' timestamp seems very low; regarding as unix timestamp

8. Show which pixels are excluded

Grey regions fail the temperature confidence threshold. The full-sky display and the confidence-masked view must not be confused with independent observations.

display(Image(filename=str(map_outputs/'planck_blue_orange_masked.png')))
del sky,mask_fraction

9. Convert the measured angular power spectrum

$$T(\hat n)=\sum_{\ell m}a_{\ell m}Y_{\ell m}(\hat n),\qquad D_\ell=\frac{\ell(\ell+1)}{2\pi}C_\ell.$$ Published spectrum estimates and the best-fit theoretical spectrum are separate inputs. The estimates are not a substitute for the full observational likelihood.

observed,theory_table=load_spectra(WORK,64)
ell=observed[:,0].astype(int)
c_observed=dl_to_cl(ell,observed[:,1])
c_theory=dl_to_cl(ell,theory_table[:,1])
print('Multipoles:',ell[0],'to',ell[-1])
print('Measured C2 =',c_observed[0],'microkelvin squared')

Multipoles: 2 to 64
Measured C2 = 236.55669082755543 microkelvin squared

10. Reconstruct correlations and integrate large angles

$$C(x)=\sum_\ell\frac{2\ell+1}{4\pi}C_\ell P_\ell(x),\qquad S_{1/2}=\int_{-1}^{1/2}C(x)^2dx=\mathbf C^TK\mathbf C.$$ This is the spectrum-based statistic. The cut-sky map estimator later is different.

kernel=s_kernel(ell)
s_observed=float(s_stat(c_observed,kernel))
theta=np.linspace(0,180,361)
basis=correlation_matrix(ell,np.cos(np.deg2rad(theta)))
fig,ax=plt.subplots(figsize=(8,4))
ax.plot(theta,basis@c_observed,label='Published TT reconstruction')
ax.plot(theta,basis@c_theory,label='Theory mean')
ax.axvline(60,color='gray',ls='--'); ax.axhline(0,color='gray',lw=.5)
ax.set(xlabel='Angular separation [degrees]',ylabel=r'$C(\theta)$ [$\mu$K$^2$]')
ax.legend(); fig.savefig(WORK/'results/maps/notebook_correlation.png',dpi=180,bbox_inches='tight'); plt.close(fig)
display(Image(filename=str(WORK/'results/maps/notebook_correlation.png')))
print('Measured spectral S1/2 =',s_observed,'microkelvin to the fourth')

Measured spectral S1/2 = 6767.531524807367 microkelvin to the fourth

11. A transparent Monte Carlo experiment

$$\widehat C_\ell=C_\ell^{\rm theory}\frac{\chi^2_{2\ell+1}}{2\ell+1}.$$ Generate 10,000 pedagogical realizations with their own seed. The main study below uses larger independent reference/calibration ensembles. A p-value describes a tail under a model, not the probability that the model is true.

rng=np.random.default_rng(731445)
draws=simulate(ell,c_theory,10000,rng)
s_draws=s_stat(draws,kernel)
print('Unconditional:',lower_tail(s_draws,s_observed))
conditioned=draws.copy(); conditioned[:,0]=c_observed[0]
print('Conditioned on measured C2:',lower_tail(s_stat(conditioned,kernel),s_observed))

Unconditional: {'p': 0.07949205079492051, 'hits': 794, 'n': 10000, 'ci95': [0.07417405807409827, 0.084872408362534]}
Conditioned on measured C2: {'p': 0.48425157484251574, 'hits': 4842, 'n': 10000, 'ci95': [0.47436087667790416, 0.49404832980086755]}

12. Run the full spectral research protocol

Two independent sets of 100,000 spectra calibrate fixed statistics and declared searches. The pipeline exports figures, tables, configuration, counts and source hashes.

from cmb_lab.pipeline import run as run_spectra
spectral_config=json.loads((WORK/'config/analysis.json').read_text())
spectral_results=run_spectra(WORK,spectral_config,WORK/'results/baseline')
display(Markdown((WORK/'results/baseline/SUMMARY.md').read_text()))

Generating independent reference and calibration ensembles

Conditioning on the observed quadrupole

Calibrating multipole and angular search penalties

'created' timestamp seems very low; regarding as unix timestamp

'modified' timestamp seems very low; regarding as unix timestamp

'created' timestamp seems very low; regarding as unix timestamp

'modified' timestamp seems very low; regarding as unix timestamp

'created' timestamp seems very low; regarding as unix timestamp

'modified' timestamp seems very low; regarding as unix timestamp

'created' timestamp seems very low; regarding as unix timestamp

'modified' timestamp seems very low; regarding as unix timestamp

'created' timestamp seems very low; regarding as unix timestamp

'modified' timestamp seems very low; regarding as unix timestamp

'created' timestamp seems very low; regarding as unix timestamp

'modified' timestamp seems very low; regarding as unix timestamp

'created' timestamp seems very low; regarding as unix timestamp

'modified' timestamp seems very low; regarding as unix timestamp

'created' timestamp seems very low; regarding as unix timestamp

'modified' timestamp seems very low; regarding as unix timestamp

'created' timestamp seems very low; regarding as unix timestamp

'modified' timestamp seems very low; regarding as unix timestamp

'created' timestamp seems very low; regarding as unix timestamp

'modified' timestamp seems very low; regarding as unix timestamp

'created' timestamp seems very low; regarding as unix timestamp

'modified' timestamp seems very low; regarding as unix timestamp

{
  "quadrupole": 0.04698953010469895,
  "s_half": 0.07520924790752093,
  "parity_2_29": 0.01010989890101099,
  "variance_2_29": 0.1035289647103529
}

Results saved to /mnt/c/Users/biswa/Documents/GitHub/cmb-anomaly-lab/results/baseline

Computed results

Generated from the official PR3 TT estimates; all p-values below use an idealized full-sky Gaussian reference. They are not official Planck likelihood results.

Diagnostic	Observed	Lower-tail p	95% Monte Carlo interval
quadrupole	236.557	0.04699	0.04568 - 0.04831
s_half	6767.53	0.07521	0.07357 - 0.07685
parity_2_29	0.762881	0.01011	0.00949 - 0.01074
variance_2_29	2074.52	0.10353	0.10164 - 0.10542

Four-test minimum-p family calibration: 0.03792. All declared fixed/scan tests: 0.03013. Parity cutoff scan (10-29), corrected: 0.01327. Angular threshold scan (45-90 degrees), corrected: 0.11577.

These corrections cover only the declared finite families, not all historical choices or all follow-up diagnostics.

Holding the quadrupole fixed

• s_half: conditional p = 0.48423.
• parity_2_29: conditional p = 0.02989.
• variance_2_29: conditional p = 0.40245.

Conditioning is a diagnostic. It removes the probability cost of the low quadrupole, assumes its published estimate is exact, and does not establish a physical cause.

Resolution of the experiment

Reference spectra: 100,000; independent calibration spectra: 100,000. The conditional ensemble reuses the calibration skies with C2 replaced. Zero observed events never imply zero probability. See results.json for exact counts, confidence intervals, source hashes, and software versions.

13. Inspect joint behaviour

Conditional distributions expose shared quadrupole information; they do not explain its physical origin.

display(Image(filename=str(WORK/'results/baseline/figures/null_distributions.png')))
display(Image(filename=str(WORK/'results/baseline/figures/dependencies.png')))

14. Define the matched map experiment

On retained pixels fit and remove monopole and dipole with $P=I-X(X^TX)^{-1}X^T$. The same beam, pixel window, confidence mask and projection are applied to every Gaussian sky. The finite-band correlation divides the masked-field correlation by the mask-pair correlation. This is a CMB-only null, without foreground/noise realizations.

import runpy
map_module=runpy.run_path(str(WORK/'scripts/analyze_maps.py'))
map_config=json.loads((WORK/'config/map_analysis.json').read_text())
print(json.dumps(map_config,indent=2))

{
  "seed": 20261102,
  "n_reference": 4096,
  "n_calibration": 4096,
  "input_lmax": 128,
  "nside": 32,
  "correlation_lmax": 95,
  "smoothing_fwhm_deg": 5.0,
  "mask_min_fraction": 0.95,
  "axis_nside": 2,
  "robustness": [
    {
      "name": "stricter_mask",
      "mask_min_fraction": 0.99,
      "smoothing_fwhm_deg": 5.0
    },
    {
      "name": "broader_smoothing",
      "mask_min_fraction": 0.95,
      "smoothing_fwhm_deg": 8.0
    }
  ],
  "status": "Retrospective exploratory protocol; fixed before executing this map ensemble. Four diagnostic types and two robustness variants form one 12-column calibrated family."
}

15. Execute the measured-sky ensemble

Run 4,096 reference and 4,096 independent calibration skies under three processing variants. The hemisphere statistic repeats the full 24-axis direction search on each sky.

map_module['run'](WORK,map_config,WORK/'results/masked_sky')
display(Markdown((WORK/'results/masked_sky/SUMMARY.md').read_text()))

Verifying native SMICA input and extracting low harmonics...

Batch 1: 512/4096 skies; 23s

Batch 1: 1024/4096 skies; 27s

Batch 1: 1536/4096 skies; 30s

Batch 1: 2048/4096 skies; 34s

Batch 1: 2560/4096 skies; 37s

Batch 1: 3072/4096 skies; 41s

Batch 1: 3584/4096 skies; 44s

Batch 1: 4096/4096 skies; 48s

Batch 2: 512/4096 skies; 52s

Batch 2: 1024/4096 skies; 56s

Batch 2: 1536/4096 skies; 60s

Batch 2: 2048/4096 skies; 64s

Batch 2: 2560/4096 skies; 68s

Batch 2: 3072/4096 skies; 73s

Batch 2: 3584/4096 skies; 77s

Batch 2: 4096/4096 skies; 81s

'created' timestamp seems very low; regarding as unix timestamp

'modified' timestamp seems very low; regarding as unix timestamp

'created' timestamp seems very low; regarding as unix timestamp

'modified' timestamp seems very low; regarding as unix timestamp

'created' timestamp seems very low; regarding as unix timestamp

'modified' timestamp seems very low; regarding as unix timestamp

'created' timestamp seems very low; regarding as unix timestamp

'modified' timestamp seems very low; regarding as unix timestamp

'created' timestamp seems very low; regarding as unix timestamp

'modified' timestamp seems very low; regarding as unix timestamp

'created' timestamp seems very low; regarding as unix timestamp

'modified' timestamp seems very low; regarding as unix timestamp

/home/biswajit/.venvs/cmb-anomaly-lab/lib/python3.10/site-packages/healpy/visufunc.py:225: UserWarning: Ignoring specified arguments in this call because figure with num: 1 already exists
  f = pylab.figure(fig, figsize=(8.5, 5.4))

'created' timestamp seems very low; regarding as unix timestamp

'modified' timestamp seems very low; regarding as unix timestamp

{
  "settings": {
    "name": "baseline",
    "mask_min_fraction": 0.95,
    "smoothing_fwhm_deg": 5.0
  },
  "f_sky": 0.8118489583333334,
  "statistics": {
    "cut_sky_s_half": {
      "observed": 2036.8056017035017,
      "tail": "lower",
      "p": 0.008054674151818403,
      "hits": 32,
      "n": 4096,
      "ci95": [
        0.005349704647585767,
        0.011011168204278781
      ]
    },
    "masked_variance": {
      "observed": 1590.5224597111323,
      "tail": "lower",
      "p": 0.04369050524774225,
      "hits": 178,
      "n": 4096,
      "ci95": [
        0.03741887596103993,
        0.050156704628523015
      ]
    },
    "pseudo_parity": {
      "observed": 0.7624026736607924,
      "tail": "lower",
      "p": 0.018550158652672687,
      "hits": 75,
      "n": 4096,
      "ci95": [
        0.01442895004305957,
        0.022899003396203203
      ]
    },
    "hemisphere_max_contrast": {
      "observed": 0.31084859602326814,
      "tail": "upper",
      "p": 0.08054674151818404,
      "hits": 329,
      "n": 4096,
      "ci95": [
        0.07217830716617371,
        0.08907056200303796
      ]
    }
  },
  "family": {
    "p": 0.03148645350256285,
    "hits": 128,
    "n": 4096,
    "ci95": [
      0.026136084358227417,
      0.03704599498741056
    ],
    "min_local_p": 0.007566512082011228,
    "n_tests": 4
  },
  "variance_closure": {
    "analytic": 2087.657814809686,
    "simulation_mean": 2089.5817243425627,
    "standard_errors": 0.3555707893665821
  },
  "quadrature_relative_change": 4.50195436485501e-13,
  "axis_count": 24,
  "maximizing_axis_galactic_lon_lat_deg": [
    67.5,
    41.81031489577859
  ]
}

Measured SMICA map analysis

Matched finite-band map estimators under an ideal Gaussian CMB-only null; not Planck end-to-end likelihood significances.

Variant	Statistic	Observed	Tail p	95% MC interval
baseline	cut_sky_s_half	2036.81	0.00805	0.00535 - 0.01101
baseline	masked_variance	1590.52	0.04369	0.03742 - 0.05016
baseline	pseudo_parity	0.762403	0.01855	0.01443 - 0.02290
baseline	hemisphere_max_contrast	0.310849	0.08055	0.07218 - 0.08907
stricter_mask	cut_sky_s_half	1985.11	0.00659	0.00415 - 0.00929
stricter_mask	masked_variance	1582.49	0.04174	0.03560 - 0.04807
stricter_mask	pseudo_parity	0.774671	0.02441	0.01969 - 0.02935
stricter_mask	hemisphere_max_contrast	0.312029	0.08055	0.07218 - 0.08907
broader_smoothing	cut_sky_s_half	1767.85	0.00610	0.00376 - 0.00871
broader_smoothing	masked_variance	1217.25	0.05370	0.04678 - 0.06080
broader_smoothing	pseudo_parity	0.729427	0.02075	0.01639 - 0.02533
broader_smoothing	hemisphere_max_contrast	0.355647	0.09617	0.08710 - 0.10538

Twelve-column family calibration: 0.03100. This covers these variants, not historical statistic selection.

The hemisphere statistic includes the complete finite direction search on every realization. The pseudo-parity statistic is a mask-dependent diagnostic, not a full-sky power estimate.

16. Inspect the correlation and robustness

The bands are pointwise simulation ranges. The numerical uncertainty intervals do not cover all observational systematics.

display(Image(filename=str(WORK/'results/masked_sky/figures/cut_sky_correlation.png')))
display(Image(filename=str(WORK/'results/masked_sky/figures/map_robustness.png')))

17. Propagate a physical power-spectrum hypothesis

$$\mathcal P_{\mathcal R}(k)=A_s(k/k_*)^{n_s-1}[1-Ae^{-(k/k_c)^2}],\qquad C_\ell^{XY}=4\pi\int d\ln k\,\mathcal P_{\mathcal R}(k)\Delta_\ell^X\Delta_\ell^Y.$$ CAMB computes TT, TE and EE together. The code searches 161 distinct grid models and repeats that search on 100,000 null spectra. The low-ell likelihood expression is an ideal full-sky proxy on published TT estimates, not the official Planck likelihood.

subprocess.run([sys.executable,str(WORK/'scripts/analyze_primordial.py')],check=True)
physical_results=json.loads((WORK/'results/primordial/results.json').read_text())
print('Search-calibrated p:',physical_results['search_calibrated_p'])
display(Image(filename=str(WORK/'results/primordial/figures/primordial_search.png')))

Computing CAMB scalar transfer functions...

'created' timestamp seems very low; regarding as unix timestamp
'modified' timestamp seems very low; regarding as unix timestamp
'created' timestamp seems very low; regarding as unix timestamp
'modified' timestamp seems very low; regarding as unix timestamp
'created' timestamp seems very low; regarding as unix timestamp
'modified' timestamp seems very low; regarding as unix timestamp

'created' timestamp seems very low; regarding as unix timestamp
'modified' timestamp seems very low; regarding as unix timestamp
'created' timestamp seems very low; regarding as unix timestamp
'modified' timestamp seems very low; regarding as unix timestamp

{
  "config": {
    "seed": 20261104,
    "n_calibration": 100000,
    "fit_lmin": 2,
    "fit_lmax": 29,
    "prediction_lmax": 128,
    "cutoff_min_mpc_inverse": 3e-05,
    "cutoff_max_mpc_inverse": 0.0008,
    "cutoff_count": 16,
    "amplitude_steps": 11,
    "cosmology": {
      "H0": 67.36,
      "ombh2": 0.02237,
      "omch2": 0.12,
      "mnu": 0.06,
      "omk": 0.0,
      "tau": 0.0544,
      "As": 2.1e-09,
      "ns": 0.9649
    },
    "status": "Retrospective phenomenological model test. Fixed background near Planck 2018; ideal full-sky TT likelihood proxy, not the official likelihood or an inflationary potential reconstruction."
  },
  "created_utc": "2026-11-04T09:48:12.792196+00:00",
  "best_amplitude": 1.0,
  "best_cutoff_mpc_inverse": 0.0002677731800328679,
  "best_twice_log_likelihood_improvement": 1.0264184618785537,
  "search_calibrated_p": {
    "p": 0.2649973500264997,
    "hits": 26499,
    "n": 100000,
    "ci95": [
      0.26225651822213886,
      0.26773711837085806
    ]
  },
  "grid_entries": 176,
  "distinct_models": 161,
  "baseline_max_fractional_TT_difference_from_official_ell2_29": 0.0049486395484374945,
  "accuracy_check": {
    "accuracy_boost": 2,
    "best_model_improvement": 1.0260204212755628,
    "improvement_absolute_difference": 0.00039804060299086785,
    "max_TT_fractional_change_ell2_29": 0.0007834546801963604
  },
  "interpretation": "Exploratory fixed-background, unlensed linear scalar transfer model. Published TT points enter an ideal full-sky likelihood proxy. Search-calibrated frequentist p, not a Bayes factor or posterior. EE/TE are predictions only.",
  "packages": {
    "camb": "2.0.5",
    "numpy": "2.2.6",
    "scipy": "1.15.3",
    "SciencePlots": "2.2.2"
  },
  "input_sha256": {
    "observed": "ccf3113604020536f6f13ccf51680a7316ad0f32da558eee7f625e613bdd5522",
    "theory": "50eecd4422e22026070bdacf9fa8a7e7ef3badf4fa50fffea040df669715cf51"
  },
  "source_sha256": {
    "scripts/analyze_primordial.py": "11de4a160469e3f0f5cbdd002e9626b11645704b2a3f3a8fce393bbe5f0fc557",
    "src/cmb_lab/primordial.py": "b2b46e44a973bb350239b9b55070a66255c7ea0d532337bd0d54119b04ad4a0d"
  }
}

Search-calibrated p: {'p': 0.2649973500264997, 'hits': 26499, 'n': 100000, 'ci95': [0.26225651822213886, 0.26773711837085806]}

18. Examine polarization predictions

$$\mathbb E[\widehat C_\ell^{EE}\mid\widehat C_\ell^{TT}]=C_\ell^{EE}-\frac{(C_\ell^{TE})^2}{C_\ell^{TT}}+\left(\frac{C_\ell^{TE}}{C_\ell^{TT}}\right)^2\widehat C_\ell^{TT}.$$ These EE/TE curves are theoretical predictions. This notebook does not analyse a measured polarization likelihood.

display(Image(filename=str(WORK/'results/primordial/figures/transfer_predictions.png')))
print('No compelling preference for this suppression family:',physical_results['search_calibrated_p']['p']>.05)

No compelling preference for this suppression family: True

19. Save an execution record and downloadable outputs

The archive contains generated results and provenance, excluding intermediate simulation arrays. It excludes the 2 GB input map. The original file remains retrievable from the recorded archive URL.

from datetime import datetime, timezone
import zipfile
execution={'author':'Biswajit Jana','executed_utc':datetime.now(timezone.utc).isoformat(),
    'python':platform.python_version(),'analysis_revision':REVISION,
    'environment':'Google Colab' if 'google.colab' in sys.modules else 'Local Python kernel',
    'real_data_sha256':{k:r['sha256'] for k,r in manifest.items()},
    'spectral_samples':200000,'map_samples':8192,'primordial_bootstrap_samples':100000}
(WORK/'results/notebook_execution.json').write_text(json.dumps(execution,indent=2)+'\n')
archive=WORK/'tmp/cmb-notebook-results.zip'; archive.parent.mkdir(exist_ok=True)
with zipfile.ZipFile(archive,'w',zipfile.ZIP_DEFLATED) as z:
    for path in (WORK/'results').rglob('*'):
        if path.is_file() and path.suffix in {'.json','.csv','.png','.pdf','.md'} and 'quick' not in path.parts:
            z.write(path,path.relative_to(WORK))
print(json.dumps(execution,indent=2))
print('Saved: tmp/cmb-notebook-results.zip; use the Colab Files sidebar to download.')

{
  "author": "Biswajit Jana",
  "executed_utc": "2026-11-04T09:48:18.731246+00:00",
  "python": "3.10.6",
  "analysis_revision": "47282e0e5f2b1c63f2bdbae87c213a453d133a8e",
  "environment": "Local Python kernel",
  "real_data_sha256": {
    "observed": "ccf3113604020536f6f13ccf51680a7316ad0f32da558eee7f625e613bdd5522",
    "theory": "50eecd4422e22026070bdacf9fa8a7e7ef3badf4fa50fffea040df669715cf51",
    "smica": "60952c645eb33d151905ddf5837477e15ca02a3261feca4ceca3c3feece4f9ac"
  },
  "spectral_samples": 200000,
  "map_samples": 8192,
  "primordial_bootstrap_samples": 100000
}
Saved: tmp/cmb-notebook-results.zip; use the Colab Files sidebar to download.

Conclusions and credit

The results quantify estimator-dependent large-angle rarity, shared low-multipole information and a weak improvement for one primordial alternative. They do not establish a new physical origin for CMB anomalies.

Author: Biswajit Jana. Please cite the project, Planck measurements, HEALPix/healpy, CAMB and SciencePlots. The project report contains full references and limitations. The annotated literature guide distinguishes reading depth and prior art.