Quadratic Estimator Module

Quadratic Estimator for Cosmic Birefringence Module

This module implements quadratic estimator (QE) methods for reconstructing cosmic birefringence rotation angle fields from CMB polarization observations.

The module provides: - Wiener filtering of E and B modes using C⁻¹ filters - Quadratic estimator reconstruction of rotation angle alms - N0 and RDN0 bias estimation - Mean field subtraction - Likelihood analysis for cosmic birefringence amplitude

Classes

FilterEB

Implements C⁻¹ Wiener filtering for E and B mode polarization using the curvedsky library. Handles component-separated maps with noise and beam deconvolution.

QE

Quadratic estimator class for cosmic birefringence reconstruction. Computes rotation angle power spectra with various bias correction methods (analytical N0, realization-dependent N0, mean field).

Simulation Index Structure (sim_config is required):

From sim_config={‘set1’: 400, ‘reuse_last’: 100} and nsims_mf=100: - stat_index: 0-299 (set1 - nsims_mf) - Statistics and OCL computation - mf_index: 300-399 (last nsims_mf from set1) - Mean field simulations - vary_index: 300-399 (last reuse_last from set1) - Varying alpha (CMB mode=’vary’) - const_index: 400-499 (reuse_last sims) - Constant alpha (CMB mode=’constant’) - null_index: 500-599 (reuse_last sims) - Null alpha (CMB mode=’null’)

Note: When nsims_mf equals reuse_last, mf_index and vary_index are identical.

Bias Estimation: - MCN0(‘stat’): Uses stat_index (0-299) - MCN0(‘vary’): Uses vary_index (300-399) - MCN0(‘const’): Uses const_index (400-499) - N1 = MCN0(‘const’) - MCN0(‘vary’) - Nlens: Uses null_index (500-599) and MCN0(‘vary’)

AcbLikelihood

Likelihood analysis class for constraining the cosmic birefringence amplitude A_CB from reconstructed rotation angle power spectra.

Algorithm

The quadratic estimator exploits EB correlations induced by cosmic birefringence:

1. Apply C⁻¹ Wiener filter to observed E and B maps

2. Compute quadratic combination: qlm ∝ Ẽlm B̃*lm

3. Apply normalization from Fisher information

4. Estimate and subtract biases (N0, mean field)

5. Compare to theoretical prediction for A_CB

Requirements

Requires the curvedsky library for efficient curved-sky operations. Install from: https://github.com/toshiyan/cmblensplus

Example

from cobi.quest import FilterEB, QE from cobi.simulation import LATsky, Mask

# Set up sky and mask sky = LATsky(libdir, nside=512) mask = Mask(libdir, nside=512, mask_str=’LAT’)

# Initialize Wiener filter filter_eb = FilterEB(

sky=sky, mask=mask, lmax=3000, fwhm=2.0, # arcmin sht_backend=’ducc0’

)

# Create quadratic estimator qe = QE(

filter=filter_eb, lmin=30, lmax=300, recon_lmax=100, norm_nsim=100

)

# Reconstruct rotation angle for simulation 0 qlm = qe.qlm(idx=0)

# Get bias-corrected power spectrum cl_aa = qe.qcl(idx=0, rm_bias=True, rdn0=True)

# Run likelihood analysis likelihood = AcbLikelihood(qe, lmin=2, lmax=50) samples = likelihood.samples(nwalkers=100, nsamples=2000)

Notes

The curvedsky library must be installed to use this module. MPI parallelization is supported for RDN0 computation.

class FilterEB(sky, mask, lmax, fwhm=2, sht_backend='healpy')

Bases: object

Parameters:

• sky (LATsky)

• mask (Mask)

• lmax (int)

• fwhm (float)

• sht_backend (str)

__init__(sky, mask, lmax, fwhm=2, sht_backend='healpy')

Parameters:

• sky (LATsky)

• mask (Mask)

• lmax (int)

• fwhm (float)

• sht_backend (str)

property Bl

convolved_EB(idx, split=0)

convolve the component separated map with the beam

Parameters:

idx (int : index of the simulation)

NL(idx, split=0)

array manipulation of noise spectra obtained by ILC weight for the filtering process

QU(idx, split=0)

deconvolve the beam from the QU map

Parameters:

idx (int : index of the simulation)

cinv_EB(idx, split=0, test=False)

C inv Filter for the component separated maps

Parameters:

• idx (int : index of the simulation)

• test (bool : if True, run the filter for 10 iterations (not cached))

check_file_exist(nsims=600, split=0)

plot_cinv(idx, split=0, lmin=2, lmax=3000)

plot the cinv filtered Cls for a given idx

Parameters:

idx (int : index of the simulation)

class QE(filter, lmin, lmax, recon_lmax, nsims_mf=100, nlb=10, lmax_bin=1024)

Bases: object

Parameters:

• filter (FilterEB)

• lmin (int)

• lmax (int)

• recon_lmax (int)

__init__(filter, lmin, lmax, recon_lmax, nsims_mf=100, nlb=10, lmax_bin=1024)

Parameters:

• filter (FilterEB)

• lmin (int)

• lmax (int)

• recon_lmax (int)

property cl_aa

qlm(idx)

check_file_exist()

qcl(idx, n0=None, mf=False, nlens=False, n1=False, binned=False)

mean_field()

mean_field_cl()

RDN0(idx)

RDN0_mpi(idx)

MPI-parallel version of RDN0(…) using mpi4py.

Parallelizes the 100 Monte-Carlo iterations; safe for any number of ranks. Only rank 0 touches the on-disk cache.

N0_sim(idx, which='vary', debug=False)

Calculate the N0 bias from the Reconstructed potential using filtered Fields with different CMB fields. If E modes is from ith simulation then B modes is from (i+1)th simulation

idx: int : index of the N0 which: str : ‘stat’, ‘vary’, or ‘const’ to select which index range to use debug: bool : if True, print the indices used for computation and return None

Index ranges and wrapping: - ‘stat’: stat_index, wraps within stat range - ‘vary’: vary_index, wraps within vary range - ‘const’: const_index, wraps within const range

Requires sim_config to be set, or manually set the corresponding index array.

MCN0(which='vary')

Monte Carlo average of N0_sim over specified index range

which: str‘stat’, ‘vary’, or ‘const’ to select which index range to use

‘stat’ uses stat_index ‘vary’ uses vary_index ‘const’ uses const_index

Requires sim_config to be set, or manually set the corresponding index arrays. Note: vary_index overlaps with mf_index only when nsims_mf equals reuse_last.

N1(binned=False)

N1 bias: difference between MCN0 for constant and varying alpha N1 = MCN0(‘const’) - MCN0(‘vary’)

Nlens(MCN0=True, binned=False)

Lensing bias: average qcl on null_index minus MCN0(‘vary’) Nlens = <qcl(null_alpha)> - MCN0(‘vary’)

RDN0_stat()

RDN0 for all stat_index simulations

qcl_stat(n0=None, mf=False, nlens=False, n1=False, binned=True)

class CrossQE(filter, lmin, lmax, recon_lmax, nsims_mf=100, nlb=10, lmax_bin=1024)

Bases: object

Parameters:

• filter (FilterEB)

• lmin (int)

• lmax (int)

• recon_lmax (int)

__init__(filter, lmin, lmax, recon_lmax, nsims_mf=100, nlb=10, lmax_bin=1024)

Parameters:

• filter (FilterEB)

• lmin (int)

• lmax (int)

• recon_lmax (int)

property cl_aa

precompute_split_ocl(splits=(1, 2, 3, 4))

Computes split-level total spectra: oclE[s] = ClEE_lensed + <N_EE_split_s> oclB[s] = ClBB_lensed + <N_BB_split_s> Saves dict {s: (oclE, oclB)}.

precompute_pair_norms(pairs=((1, 2), (3, 4), (1, 3), (2, 4), (1, 4), (2, 3)))

Precomputes normalization arrays for each split pair.

Returns dict: norms[(i,j)] = {“EiBj”: norm_ij, “EjBi”: norm_ji} where norm_ij normalizes Q(E^i, B^j).

qlm_pair(idx, si, sj)

Parameters:

• si (int)

• sj (int)

mean_field_sim(idx, splits=(1, 2, 3, 4))

Cross-only lensing power estimate: average over cross-spectra of reconstructions from disjoint split pairs. For 4 splits, uses the 3 disjoint pairings.

mean_field()

mean_field_cl()

qcl_cross_only(idx, splits=(1, 2, 3, 4))

Cross-only lensing power estimate: average over cross-spectra of reconstructions from disjoint split pairs. For 4 splits, uses the 3 disjoint pairings.

Nlens(binned=False)

Lensing bias: average qcl on null_index minus MCN0(‘vary’) Nlens = <qcl(null_alpha)> - MCN0(‘vary’)

N1(binned=False)

N1 bias: difference between MCN0 for constant and varying alpha N1 = MCN0(‘const’) - MCN0(‘vary’)

qcl_stat(n0=None, n1=False, mf=False, nlens=False, binned=False)

Compute statistics of qcl over stat_sims with various bias corrections.

Parameters:

• n0 (str or None) – N0 bias correction: None (no correction), ‘mcn0’, or ‘rdn0’

• mf (bool) – If True, subtract mean field bias

• nlens (bool) – If True, subtract lensing bias

• binned (bool) – If True, return binned spectra

Returns:

cl – Array of shape (nsims, nlmax+1) or (nsims, nbins) if binned

Return type:

array

N0_sim(idx, which='vary', debug=False)

Calculate the N0 bias from the Reconstructed potential using filtered Fields with different CMB fields. If E modes is from ith simulation then B modes is from (i+1)th simulation

idx: int : index of the N0 which: str : ‘stat’, ‘vary’, or ‘const’ to select which index range to use debug: bool : if True, print the indices used for computation and return None

Index ranges and wrapping: - ‘stat’: stat_index, wraps within stat range - ‘vary’: vary_index, wraps within vary range - ‘const’: const_index, wraps within const range

Requires sim_config to be set, or manually set the corresponding index array.

MCN0(which='vary')

RDN0(idx, navg=100)

Realization-dependent N0 (RDN0) compatible with qcl_cross_only() (symmetric estimator).

Uses the same Eq.(43)-style structure: RDN0 = Cx(ds,ds) + Cx(ds,sd) + Cx(sd,ds) + Cx(sd,sd) - Cx(ss0,ss0) - Cx(ss0,s0s) where Cx is the cross-only spectrum operator (12×34+13×24+14×23)/3.

Returns: array C_L[0..recon_lmax]

Parameters:

• idx (int)

• navg (int)

RDN0_mpi(idx, navg=100)

MPI-parallel cross-only RDN0 for SYM estimator (compatible with qcl_cross_only()).

Implements the Eq.(43)-style combination (same algebra as your rot version), but with SYM estimator: RDN0 = Cx(ds,ds)+Cx(ds,sd)+Cx(sd,ds)+Cx(sd,sd) - Cx(ss0,ss0) - Cx(ss0,s0s)

Parallelizes navg Monte-Carlo draws over (s, s0) pairs. Only rank 0 reads/writes the cache.

Parameters:

• idx (int)

• navg (int)

RDN0_stat(array=False)

RDN0 for all stat_index simulations

class CrossQEv1(filter, lmin, lmax, recon_lmax, nsims_mf=100, nlb=10, lmax_bin=1024)

Bases: object

Cross-only quadratic estimator for cosmic birefringence, implementing Eq. 38 of Madhavacheril et al. (2020).

Requires 4 splits with independent noise. Uses the O(m^2) algorithm from Section 3.2 of the paper.

Parameters:

• lmin (int)

• lmax (int)

• recon_lmax (int)

__init__(filter, lmin, lmax, recon_lmax, nsims_mf=100, nlb=10, lmax_bin=1024)

Parameters:

• lmin (int)

• lmax (int)

• recon_lmax (int)

property cl_aa

four_split_phi(idx)

CHANGED: Completely rewritten to follow the so-lenspipe implementation.

Computes all the intermediate quantities needed for Eq. 38:

• phi^(ij) for all i,j (including diagonals)

• phi_hat = (1/m^2) sum_{ij} phi^(ij) (full coadd reconstruction)

• phi^x = phi_hat - (1/m^2) sum_i phi^(ii) (cross coadd, no diags)

• phi^(i) = (1/m) sum_j phi^(ij) (per-split marginals)

• phi^(i)x = phi^(i) - (1/m) phi^(ii) (per-split cross)

Returns a dict with all quantities needed by split_phi_to_cl.

The original code only computed phi for disjoint pairs and averaged them at the map level, which does not implement Eq. 38.

split_phi_to_cl(xy, uv=None)

CHANGED: Completely new method implementing Eq. 38.

Computes the cross-only power spectrum:

C_L^x(XY, UV) = 1/[m(m-1)(m-2)(m-3)] * [

m^4 * C_L(phi_X, phi_X’)

• 4*m^2 * sum_i C_L(phi_ix, phi_ix’)

• 4 * sum_{i<j} C_L(phi_ij, phi_ij’)

]

This is the O(m^2) algorithm from Section 3.2 of the paper. Follows the so-lenspipe split_phi_to_cl function exactly.

Parameters:

• xy – dict from four_split_phi (first reconstruction, e.g. data)

• uv – dict from four_split_phi (second reconstruction, for cross). If None, compute auto-spectrum.

qlm(idx)

CHANGED: This now returns the cross-coadd map phi^x for mean-field etc. Note: this is biased by (1 - 1/m), but that’s fine for mean-field estimation where we average and subtract.

qcl(idx, n0=None, mf=False, nlens=False, n1=False, binned=False)

CHANGED: The power spectrum is now computed via split_phi_to_cl (Eq. 38) rather than by taking the auto-spectrum of an averaged map.

mean_field_sim(idx)

CHANGED: Mean field per-sim is now just the cross-coadd map phi^x. The original averaged phi over 3 disjoint pairings at the map level, which is not the correct cross-only mean field.

mean_field()

Mean field: average of phi^x over mf_index sims.

mean_field_cl()

N0_sim(idx, which='vary')

CHANGED: N0 for the cross-only estimator.

Uses Eq. 38 applied to reconstructions where E and B come from DIFFERENT simulations (to get disconnected/Gaussian piece only).

For each split pair (i,j), the symmetrized estimator is:

phi^(ij) = (1/2) * norm * [Q(E_i^sim1, B_j^sim2) + Q(E_j^sim1, B_i^sim2)]

Note: for i != j in the data case, we had 4 terms for full symmetrization over both split indices AND E/B leg assignment. For N0 with two sims, the symmetrization is: swap which sim provides E vs B for each split.

MCN0(which='vary')

Monte Carlo N0: average N0_sim over the relevant index range.

N1(binned=False)

N1 = MCN0(‘const’) - MCN0(‘vary’)

Nlens(MCN0=True, binned=False)

Lensing bias from null sims.

RDN0(idx)

CHANGED: Realization-dependent N0 for the cross-only estimator.

Implements Eq. 43 of the paper: apply the C^x algorithm (Eq. 38) to each of the 4-point data/sim combinations from the standard RDN0 formula (Eq. 17), then average over sim realizations.

For each MC iteration with sims s1, s2:

RDN0 = <C^x(d,s1) + C^x(s1,d) - C^x(s1,s2) - C^x(s2,s1)>

where C^x(a,b) means: E legs from ‘a’, B legs from ‘b’.

qcl_stat(n0=None, mf=False, nlens=False, n1=False, binned=True)

Power spectra for all stat_index simulations.

Filtering Class

class FilterEB(sky, mask, lmax, fwhm=2, sht_backend='healpy')

Parameters:

• sky (LATsky)

• mask (Mask)

• lmax (int)

• fwhm (float)

• sht_backend (str)

__init__(sky, mask, lmax, fwhm=2, sht_backend='healpy')

Parameters:

• sky (LATsky)

• mask (Mask)

• lmax (int)

• fwhm (float)

• sht_backend (str)

property Bl

convolved_EB(idx, split=0)

convolve the component separated map with the beam

Parameters:

idx (int : index of the simulation)

NL(idx, split=0)

array manipulation of noise spectra obtained by ILC weight for the filtering process

QU(idx, split=0)

deconvolve the beam from the QU map

Parameters:

idx (int : index of the simulation)

cinv_EB(idx, split=0, test=False)

C inv Filter for the component separated maps

Parameters:

• idx (int : index of the simulation)

• test (bool : if True, run the filter for 10 iterations (not cached))

check_file_exist(nsims=600, split=0)

plot_cinv(idx, split=0, lmin=2, lmax=3000)

plot the cinv filtered Cls for a given idx

Parameters:

idx (int : index of the simulation)

Quadratic Estimator Class

class QE(filter, lmin, lmax, recon_lmax, nsims_mf=100, nlb=10, lmax_bin=1024)

Parameters:

• filter (FilterEB)

• lmin (int)

• lmax (int)

• recon_lmax (int)

__init__(filter, lmin, lmax, recon_lmax, nsims_mf=100, nlb=10, lmax_bin=1024)

Parameters:

• filter (FilterEB)

• lmin (int)

• lmax (int)

• recon_lmax (int)

property cl_aa

qlm(idx)

check_file_exist()

qcl(idx, n0=None, mf=False, nlens=False, n1=False, binned=False)

mean_field()

mean_field_cl()

RDN0(idx)

RDN0_mpi(idx)

MPI-parallel version of RDN0(…) using mpi4py.

Parallelizes the 100 Monte-Carlo iterations; safe for any number of ranks. Only rank 0 touches the on-disk cache.

N0_sim(idx, which='vary', debug=False)

Calculate the N0 bias from the Reconstructed potential using filtered Fields with different CMB fields. If E modes is from ith simulation then B modes is from (i+1)th simulation

idx: int : index of the N0 which: str : ‘stat’, ‘vary’, or ‘const’ to select which index range to use debug: bool : if True, print the indices used for computation and return None

Index ranges and wrapping: - ‘stat’: stat_index, wraps within stat range - ‘vary’: vary_index, wraps within vary range - ‘const’: const_index, wraps within const range

Requires sim_config to be set, or manually set the corresponding index array.

MCN0(which='vary')

Monte Carlo average of N0_sim over specified index range

which: str‘stat’, ‘vary’, or ‘const’ to select which index range to use

‘stat’ uses stat_index ‘vary’ uses vary_index ‘const’ uses const_index

Requires sim_config to be set, or manually set the corresponding index arrays. Note: vary_index overlaps with mf_index only when nsims_mf equals reuse_last.

N1(binned=False)

N1 bias: difference between MCN0 for constant and varying alpha N1 = MCN0(‘const’) - MCN0(‘vary’)

Nlens(MCN0=True, binned=False)

Lensing bias: average qcl on null_index minus MCN0(‘vary’) Nlens = <qcl(null_alpha)> - MCN0(‘vary’)

RDN0_stat()

RDN0 for all stat_index simulations

qcl_stat(n0=None, mf=False, nlens=False, n1=False, binned=True)

Likelihood Class