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).

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)

convolve the component separated map with the beam

Parameters:

idx (int : index of the simulation)

NL(idx)

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

QU(idx)

deconvolve the beam from the QU map

Parameters:

idx (int : index of the simulation)

cinv_EB(idx, 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)

check_file_exist()

plot_cinv(idx, 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=2, lmax_bin=100)

Bases: object

Parameters:

• filter (FilterEB)

• lmin (int)

• lmax (int)

• recon_lmax (int)

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

Parameters:

• filter (FilterEB)

• lmin (int)

• lmax (int)

• recon_lmax (int)

property cl_aa

qlm(idx)

check_file_exist()

qcl(idx, rm_bias=False, rdn0=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')

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 : ‘vary’ or ‘cons’ to select which alpha_config index to use

MCN0(which='vary')

Monte Carlo average of N0_sim over N0 indexes

which: str : ‘vary’ or ‘cons’ to select which alpha_config index to use

N1()

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

Nlens(MCN0=False)

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

qcl_stat(rm_bias=True, rdn0=False, binned=True)

class AcbLikelihood(qelib, lmin=2, lmax=50, rdn0=False, simple_chi=False)

Bases: object

Parameters:

qelib (QE)

__init__(qelib, lmin=2, lmax=50, rdn0=False, simple_chi=False)

Parameters:

qelib (QE)

theory(Acb)

plot(Acb)

ln_prior(Acb)

ln_likelihood(Acb)

ln_posterior(Acb)

sampler(nwalkers=32, nsamples=1000)

samples(nwalkers=100, nsamples=2000, discard=300)

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)

convolve the component separated map with the beam

Parameters:

idx (int : index of the simulation)

NL(idx)

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

QU(idx)

deconvolve the beam from the QU map

Parameters:

idx (int : index of the simulation)

cinv_EB(idx, 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)

check_file_exist()

plot_cinv(idx, 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=2, lmax_bin=100)

Parameters:

• filter (FilterEB)

• lmin (int)

• lmax (int)

• recon_lmax (int)

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

Parameters:

• filter (FilterEB)

• lmin (int)

• lmax (int)

• recon_lmax (int)

property cl_aa

qlm(idx)

check_file_exist()

qcl(idx, rm_bias=False, rdn0=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')

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 : ‘vary’ or ‘cons’ to select which alpha_config index to use

MCN0(which='vary')

Monte Carlo average of N0_sim over N0 indexes

which: str : ‘vary’ or ‘cons’ to select which alpha_config index to use

N1()

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

Nlens(MCN0=False)

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

qcl_stat(rm_bias=True, rdn0=False, binned=True)

Likelihood Class

class AcbLikelihood(qelib, lmin=2, lmax=50, rdn0=False, simple_chi=False)

Parameters:

qelib (QE)

__init__(qelib, lmin=2, lmax=50, rdn0=False, simple_chi=False)

Parameters:

qelib (QE)

theory(Acb)

plot(Acb)

ln_prior(Acb)

ln_likelihood(Acb)

ln_posterior(Acb)

sampler(nwalkers=32, nsamples=1000)

samples(nwalkers=100, nsamples=2000, discard=300)