A Bayesian ILC method for CMB B-mode posterior estimation and reconstruction of primordial gravity wave signal
DD RAFT VERSION O CTOBER
1, 2020
Preprint typeset using L A TEX style emulateapj v. 12/16/11
A BAYESIAN ILC METHOD FOR CMB B-MODE POSTERIOR ESTIMATION AND RECONSTRUCTION OFPRIMORDIAL GRAVITY WAVE SIGNAL. S ARVESH K UMAR Y ADAV , R AJIB S AHA Draft version October 1, 2020
ABSTRACTThe Cosmic Microwave Background (CMB) radiation B mode polarization signal contains the unique sig-nature of primordial metric perturbations produced during the inflation. The separation of the weak CMBB-mode signal from strong foreground contamination in observed maps is a complex task, and proposed newgeneration low noise satellite missions compete with the weak signal level of this gravitational background. Inthis article, for the first time, we employ a foreground model-independent internal linear combination (ILC)method to reconstruct the CMB B mode signal using simulated observations over large angular scales of thesky of 6 frequency bands of future generation CMB mission Probe of Inflation and Cosmic Origins (PICO).We estimate the joint CMB B mode posterior density following the interleaving Gibbs steps of B mode angularpower spectrum and cleaned map samples using the ILC method. We extend and improve the earlier reportedBayesian ILC method to analyze weak CMB B mode reconstruction by introducing noise bias corrections attwo stages during the ILC weight estimation. By performing 200 Monte Carlo simulations of the Bayesian ILCmethod, we find that our method can reconstruct the CMB signals and the joint posterior density accuratelyover large angular scales of the sky. We estimate Blackwell-Rao statistics of the marginal density of CMB Bmode angular power spectrum and use them to estimate the joint density of scalar to tensor ratio r and a lensingpower spectrum amplitude A lens . Using 200 Monte Carlo simulations of the delensing approach, we find thatour method can achieve an unbiased detection of the primordial gravitational wave signal r with more than 8 σ significance for levels of r (cid:62) . Subject headings:
CMB B Modes, CMB B Mode Delensing, CMB Component Separation INTRODUCTION
The discovery of the CMB in the second half of the lastcentury led to an era of precision cosmology, which resultedin the demise of many cosmological models and a few’s sur-vival. The current observations have put very stringent con-straints (Planck Collaboration et al. 2018b,a) on the variousinflationary models (Martin et al. 2014), within the very suc-cessful inflationary paradigm for the origin of the primordialperturbations. The final Planck 2018 release (Planck Col-laboration et al. 2018b), ruled out the perfect scale invari-ance for the spectral index of scalar perturbations at 8.4 σ and the running and the running of the running (Chung et al.2003), of the spectral index have been negated with 95 % CL,consistent with the simplest slow-roll dynamics for the infla-ton, and the spatial curvature Ω K is − . + . − . at 95 % CL.BICEP2/Keck Array (BICEP2 Collaboration et al. 2016) to-gether with Planck 2018 strongly disfavors monomial mod-els with V ( φ ) ∝ φ p , p >
1, natural inflation, and low scaleSUSY models. The observations have also established thatthe primordial perturbations are adiabatic to a very high de-gree (Planck Collaboration et al. 2018b), and the primor-dial power spectrum does not deviate from a pure power-law(Planck Collaboration et al. 2018b). Further the Planck likeli-hood together with the B-mode polarization likelihood of theBICEP2-Keck Array puts a stringent 95 % CL upper limit of r . < .
056 corresponding to the energy scale of inflationof V / < . × GeV with 95 % CL bound.To further constraint the cosmological models and probe theenergy scale of the inflation and the existence of primordialgravitational waves, new generation of CMB space missions Physics Department, Indian Institute of Science Education and Re-search Bhopal, Bhopal, M.P, 462066, India. such as PICO (Hanany et al. 2019) COrE (Delabrouille et al.2018), LiteBIRD (Matsumura et al. 2014) and PIXIE (Kogutet al. 2011) have been proposed to detect the primordial CMBB mode on large angular scales at a level of r < − . TheCMB B-mode signal given the current bound will have typ-ical RMS fluctuations < . µ K , which are extremely weak,with the strong polarized Galactic foregrounds and instrumen-tal systematic making their detection and reconstruction ex-tremely difficult. To further add to difficulties, gravitationallensing introduces spurious cosmic variance from the lens-induced B modes or the lensing B modes. The lensing B-modes are due to the conversion of CMB E-mode to B-modedue to weak gravitational lenses along the line of sight (Seljak& Hirata 2004). Not only this biases the amplitude r but alsothe cosmic variance of the primordial CMB B-mode powerspectrum.Over the recent years, various studies dealing with the fore-grounds minimization in the context of CMB B-mode sky hasbeen undertaken (Baccigalpi et al. 2004; Betoule et al. 2009;Dunkley et al. 2009; Bonaldi & Ricciardi 2011; Katayama& Komatsu 2011; Armitage-Caplan et al. 2012; Errard et al.2016; Remazeilles et al. 2016, 2018b; Hervías-Caimapo et al.2017). Recently, new methods were proposed to investigatethe joint posterior density of CMB signal and correspondingtheoretical angular power spectrum on large scales, for CMBtemperature (Sudevan & Saha 2018a,b) and CMB E-mode(Purkayastha et al. 2020; Purkayastha et al. 2020) polariza-tion.In this article, we extend (Sudevan & Saha 2018b;Purkayastha et al. 2020) and develop a new non-parametricmethod by also taking care of detector noise to reconstructclean CMB B-mode signal and corresponding theoretical an-gular power spectrum. We perform a joint analysis of CMB a r X i v : . [ a s t r o - ph . C O ] S e p B-mode signal and its angular power spectrum posterior den-sity without considering any foreground model. Unlike thepolarized CMB, we do not have an accurate enough modelfor the polarized galactic foreground and the exact num-ber of independent polarized foregrounds is not known (Re-mazeilles et al. 2016). Our non-parametric method avoid ef-fects (Armitage-Caplan et al. 2012) due to inaccurate polar-ized Galactic foreground models. Our method also providesthe best fit estimates of both, CMB B-mode map and it’s the-oretical angular power spectrum along with their confidenceinterval regions. We apply our Bayesian ILC method to re-cover the weak CMB B-mode signal from the simulated fore-ground and noise-contaminated 6 PICO frequency channels.We look into the performance of our Bayesian ILC methodfollowing the Gibb’s procedure to reconstruct the primordialCMB B-mode signal and its theoretical angular power spec-trum. We also perform correction for lensing bias in CMBB-mode power spectrum without removing the lensing cos-mic variance contribution to B-modes. We use the samplesof the theoretical CMB B-mode power spectrum generatedat each Gibb’s step of our method and simultaneously fit theamplitude of the primordial B-mode power spectrum parame-ter and the amplitude of lensing B-mode power spectrum in aBayesian framework (Remazeilles et al. 2018a). The methodremoves the lensing bias on the posterior distribution of r andenables us to detect r with more than 8 σ significance for CMBB-mode satellite mission like PICO for r ≥ . r . Finally, in Section 6, we discuss and conclude. FORMALISM
This section discusses the formalism used to estimate thejoint posterior density of the CMB signal and its theoreti-cal angular power spectrum given the observed data. Weadopt and improve formalism used in this work for compo-nent separation as in (Sudevan & Saha 2018b; Purkayasthaet al. 2020) so that it applies to the weak CMB B-mode sig-nal. The method not only gives us the best-fit CMB B-modemap and it is best-fit power spectrum, but it also providesMCMC Gibbs samples for sky power spectrum ˆ σ B (cid:96) and the-oretical power spectrum C B (cid:96) , which we utilize to delens theangular power spectrum and to estimate of tensor-to-scalarratio r . Data Model
Given observations of CMB B-mode signal S at n differ-ent frequencies in thermodynamic temperature units, we canwrite for an observed i th frequency map X i , X i = S + F i + N i (1) where F i is the net foreground contribution from all the fore-ground components at the i th frequency channel and N i isthe corresponding detector noise. Each of the above bold-faced quantity is a column vector of size N pix represent-ing a HEALPix map where N pix = 12 N side , N side being thepixel resolution parameter, having common beam and pixelresolution. Let D denote the observed data set i.e. D = { X , X , ..., X n } . CMB Posterior Estimation
Given the observed data, D , P ( S , C B (cid:96) | D ) represents the jointdensity of CMB B-mode map, S , and the theoretical CMB Bmode angular power spectrum, C B (cid:96) . As it is difficult to obtainthe P ( S , C B (cid:96) | D ) analytically we evaluate it by drawing samplesfrom it. If we can sample from the conditional distributions P ( S | C B (cid:96) , D ) and P ( C B (cid:96) | S , D ) then utilizing Gibbs sampling ap-proach Rubin (1992), which says that samples ( S i , C B i (cid:96) ) canbe drawn from the joint distribution P ( S , C B (cid:96) | D ) by iteratingthe following symbolic sampling equations: S i + ← P ( S | C B i (cid:96) , D ) (2) C B i + (cid:96) ← P ( C B (cid:96) | S i + ) (3)The symbol “ ← ” implies that a sample of corresponding vari-ables is drawn from the distribution on the right-hand side.Once the initial burn-in period is over, the samples will con-verge to being drawn from the required joint distribution. Sampling CMB Signal
We use foreground model-independent method to drawsamples of S given the CMB B-mode theory C B (cid:96) and D . Wemodify the global ILC method described in (Sudevan & Saha2018; Purkayastha et al. 2020) to improve separating the weakCMB B-mode signal given the detector noise model. Let usassume that the mean corresponding to each frequency map X i , as discussed in (2.1), has already been subtracted. Thecleaned CMB B mode map S can be obtained by linear com-bination of n input maps X i , with weight factor w i , i.e., S = n (cid:88) i =1 w i X i . (4)Since the spectral distribution of CMB photons is a blackbodyto an excellent approximation, the CMB anisotropy signal S (in thermodynamic temperature units) is independent of thefrequency channel. In order to avoid multiplicative bias inamplitudes of CMB anisotropies the sum of weights is con-straint to unity i.e., (cid:80) ni =1 w i = 1. As discussed in (Sudevan &Saha 2018) instead of minimizing the clean map variance S T S we minimize σ = S T C † S (5)where C represents the CMB B-mode theoretical covariancematrix and † denotes the Moore-Penrose generalized inversePenrose (1955). Using Equation (4) in Equation (5) we write, σ = WAW T (6) Hierarchical Equal Area Isolatitude Pixellization of sphere, e.g., seeGórski et al. (2005) new ILC method for weak CMB B-mode Signal 3where W = ( w , w , w , ..., w n ) is a 1 × n weight row vectorand A is an n × n matrix with it element A i j given by A i j = X Ti C † X j (7)The weights that minimize the variance given by Equa-tion (6) subject to the above constraint is obtained followingLagrange’s multiplier approach (Saha et al. 2008; Tegmark &Efstathiou 1996; Saha et al. 2006; Tegmark et al. 2003) and isgiven by W = eA † eA † e T (8)where e = (1 , , ...,
1) is the 1 × n CMB shape vector in thermo-dynamic temperature units and A † is the Moore-Penrose gen-eralized inverse of the matrix A . Computing a dense matrix[ C ] N pix × N pix at every Gibbs iteration is computationally costly,hence we switch to the harmonic space where Equation (7) issimpler to compute, A i j = l max (cid:88) (cid:96) =2 (2 (cid:96) + σ i j (cid:96) C B (cid:96) , (9)where (cid:96) max denotes the maximum multipole used in he analy-sis, σ i j (cid:96) denotes the angular cross power spectrum between X i and X j channel maps and C B (cid:96) represents the beam and pixelsmoothed CMB BB theoretical power spectrum i.e., C B (cid:96) = C B (cid:48) (cid:96) B (cid:96) P (cid:96) (10)where C B (cid:48) (cid:96) does not have any smoothing effect, and B (cid:96) and P (cid:96) are respectively the polarization beam and polarizationpixel window functions. The internal linear combinationmethod for component separation performs well only in a lownoise environment, and since the CMB B-mode signal is evenweaker than CMB E mode by order of magnitude, to mini-mize the residual noise bias in the output CMB-B mode mapand power spectrum we subtract the noise auto-power initiallyfrom the input frequency cross power spectrum, A i j = l max (cid:88) (cid:96) =2 (2 (cid:96) +
1) 1 C B (cid:96) ( σ i j (cid:96) − δ i j σ N , i (cid:96) ) . (11)where σ N , i (cid:96) is noise auto power corresponding to detector at i th frequency. We use matrix A , the component for whichare given by Equation (11), in Equation (8) to obtain the rowvector W . We use the weights obtained, to sample the fore-ground minimized CMB B-mode signal S , by linearly com-bining the input channel maps X i at every Gibb’s step follow-ing the Equation (4). Sampling C B (cid:96) The signal sample S can be represented mathematically interms of spherical harmonics, S ( θ, φ ) = ∞ (cid:88) (cid:96) =2 (cid:96) (cid:88) m = − (cid:96) s (cid:96) m Y (cid:96) m ( θ, φ ) , (12)then the realization-specific power spectrum is given by ˆ σ B (cid:48) (cid:96) ≡ (cid:96) + (cid:96) (cid:88) m = − (cid:96) | s (cid:96) m | . (13) Frequency Beam FWHM Q and U noise RMS(GHz) (arcmin) ( µ K CMB arcmin )90 9.5 2.09108 7.9 1.70129 7.4 1.53155 6.2 1.28186 4.3 3.54268 3.2 2.63TABLE 1PICO
FREQUENCY MAPS USED IN THIS WORK ALONG WITH THEIRINSTRUMENTAL SPECIFICATIONS . In order to minimize the noise bias in the sampled theory C B (cid:96) we further subtract the weighted noise power from ˆ σ B (cid:48) (cid:96) to ob-tain, ˆ σ B (cid:96) = ˆ σ B (cid:48) (cid:96) − n (cid:88) i =1 w i σ N , i (cid:96) (14)Since the power spectrum only depends on the signal S through ˆ σ B (cid:96) , and not its phases, therefore to draw samples of C B (cid:96) given S , we sample from P ( C B (cid:96) | ˆ σ B (cid:96) ). The conditional den-sity P ( C B (cid:96) | ˆ σ B (cid:96) ) can be written Sudevan & Saha (2018b) as, P ( C B (cid:96) | ˆ σ B (cid:96) ) = (cid:18) C B (cid:96) (cid:19) (2 (cid:96) + / exp (cid:20) − ˆ σ B (cid:96) (2 (cid:96) + C B (cid:96) (cid:21) , (15)where the variable x = ˆ C B (cid:96) (2 (cid:96) + / C B (cid:96) is a χ distributed ran-dom variable having 2 (cid:96) − C B (cid:96) using Equation (14) we need to draw firstx from the χ distribution of 2 (cid:96) − (cid:96) − S we have estimates of ˆ σ B (cid:96) , we then obtain C B (cid:96) using C B (cid:96) = ˆ σ B (cid:96) (2 (cid:96) + / x . FREQUENCY MAPS
In this work, we simulate the foreground and noise-contaminated CMB B modes at 6 CMB dominating the leastnoisy frequency bands of proposed satellite mission PICO, inthe frequency range 90 GHz to 268 GHz. We list the fre-quency channel maps along with their instrumental specifica-tions in Table (1).
CMB B-mode Signal
We simulate lensed CMB Q and U Stoke parameter mapsfrom the lensed CMB B-mode angular power spectra gen-erated by the Boltzmann solver CAMB (Lewis et al. 2000).Since we perform our analysis on large scales ( (cid:96) ≤ λ CDM + r cosmol-ogy with optical depth to reionization τ = 0 .
055 (Planck Col-laboration et al. 2016c), A lens = 1 and other cosmological pa-rameters set to the Planck 2015 best-fit values (Planck Col-laboration et al. 2016a). Foreground B-mode Signal
We generate foreground maps at all the six PICO frequen-cies used in this work corresponding to synchrotron and ther-mal dust; two major CMB polarized foreground contribu-tors. To generate them, we follow the procedure similar to (cid:96) ( (cid:96) + ) C l BB (cid:96) / π , µ K Multipole, (cid:96) Cl
BBr =0 . + Cl BBlens Cl BBr =0 . + Cl BBlens F IG . 1.— The above plot shows the noise power spectrum corresponding to each of the six PICO frequencies used in this work, along with the lensed B-modepower spectrum for r = 0 .
01 and r = 0 .
05. From the figure, we find that the CMB B-mode theoretical angular power spectrum is well above the noise power forthe 6 PICO frequency bands used in this work. (Remazeilles et al. 2018a). We generate Q and U maps ateach frequency and use them to obtain corresponding B-modemaps for both the foregrounds.To generate the polarized Galactic synchrotron Stokesmaps, we extrapolate the Wilkinson Microwave AnisotropyProbe (WMAP) 23 GHz (Page et al. 2007) (Bennett et al.2013) stokes maps Q and U to the six PICO frequenciesthrough a power-law frequency dependence: Q sync ν ( p ) = Q ( p ) (cid:18) ν GHz (cid:19) β s (16) U sync ν ( p ) = U ( p ) (cid:18) ν GHz (cid:19) β s (17)We use a constant spectral index β s = − p is the pixel index.To simulate the Galactic polarized thermal dust Stokesmaps, we extrapolate the generalized needlet ILC (GNILC)Planck 353 GHz thermal dust optical depth map (Planck Col-laboration et al. 2016d) to the relevant PICO frequencies: Q dust ν ( p ) = f d g d ( p ) I GNILC ν ( p ) cos(2 γ d ( p )) (18) U dust ν ( p ) = f d g d ( p ) I GNILC ν ( p ) sin(2 γ d ( p )) (19)where f d is the pixel independent intrinsic dust polarizationfraction which depends on the level of depolarization alongthe line of sight, following (Delabrouille et al. 2013; Re-mazeilles et al. 2018a) we take it to be 0.15, g d is the pixeldependent geometric depolarization factor which we com-pute using the 3D Galactic magnetic field and 3D distributionalong the line of sight. To compute polarization angle γ d (De-labrouille et al. 2013) at each pixel, we use WMAP 23 GHz map after smoothing with Gaussian beam of 3 ◦ , γ d ( p ) = 12 tan − (cid:18) − U ( p ) Q ( p ) (cid:19) . (20)We compute the depolarization factor g d using WMAP 23GHz, and the residual monopole subtracted 408 MHz Haslamsynchrotron template ( I (0 . ), extrapolated to 23 GHz as-suming a constant spectral index of -3.0. To compute it,we smooth the extrapolated map to gauss beam of 3 ◦ at N side = 512 and use, g d ( p ) = (cid:113) Q ( p ) + U ( p ) f s I (0 . ( p )(23 . / . − . , (21)where for the spectral index used in above equation, the syn-chrotron polarization fraction f s = 0 .
75. The I GNILC ν is theGNILC dust intensity map free from the cosmic infraredbackground at the frequency ν and is given by the modifiedblackbody spectrum: I GNILC ν ( p ) = τ GNILC ( p ) (cid:18) ν GHz (cid:19) β d B ν ( T d ) (22)where τ GNILC is the Planck GNLIC dust optical depth at 353GHz, the dust emissivity β d = 1 . T d = 19 . K is the dusttemperature. B ν ( T d ) is the Planck function at thermal dusttemperature T d given by: B ν ( T d ) = 2 h c (cid:16) h ν B T d (cid:17) − Q and U maps to get B-mode synchrotron and thermaldust foreground maps at each of the PICO frequencies usedin this article at N side = 16. We smooth the obtained maps bypolarized Gaussian beam of FWHM 9 ◦ . new ILC method for weak CMB B-mode Signal 5 Detector Noise Simulations
We simulate Gaussian, isotropic, and pixel-pixel uncorre-lated random realizations of detector Q and U noise maps forthe six PICO (Young et al. 2018) frequency bands used in thiswork. We present detector specifications for each of the bandsin Table (1). We further assume that Q and U noise maps arepixel uncorrelated i.e. (cid:104) Q i ( p ) U i ( p (cid:48) ) (cid:105) = 0 . (24)We further assume that the pixel noise variances σ Q i for Q and σ U i for U maps at a frequency ν i are identical and given by σ Q i = σ U i = (cid:0) c ∆ Q i (cid:1) / ( ∆Ω ) (25)where ∆ Q i is the noise RMS in arcminute for Q i map, c isthe conversion factor from arcminute to radian and ∆Ω is thesolid angle subtended by single-pixel at N side =16 . We bringboth Q and U noise maps to the same beam resolution at N side = 16 by multiplying the ratio of a polarized Gaussianbeam of FWHM 9 ◦ and the polarized beam is given in Table(1) for corresponding frequency channel. We finally convertthe noise Stokes maps obtained to full sky B-mode noise mapat each of the frequencies. In Figure (1) we show the PICOdetector model noise power corresponding to the six chan-nels along with the lensed CMB B-mode theoretical powerfor r = 0 .
01 and r = 0 .
05. We can see that the noise power forall the used frequency is well below the B-mode signal powerspectrum at all multipoles used in this work. Finally, we ob-tain the simulated PICO input noisy foreground contaminatedCMB B-mode maps by combining all the three componentsfor each of the six frequencies using equation (1). METHODOLOGY
We implement our model-independent method on the sim-ulated foreground and noise-contaminated B-mode maps ob-tained above after removing the monopole and dipole compo-nents from each of them. We smooth theoretical C B (cid:96) obtainedfrom CAMB using τ = 0 . A lens = 1 and Planck 2015 best-fit values by Gauss beam of 9 ◦ and polarization pixel windowfunction corresponding to N side = 16 as in equation (10). Inorder to obtain sampled CMB B-mode sky S given the ob-served data set D and sampled theory C B i (cid:96) following the sym-bolic sampling equation (2) we first obtain A matrix follow-ing equation (11) and use it to obtain weights using equation(8). We use weights obtain in the last step to combine theinput foreground linearly, and noise-contaminated CMB B-mode maps to obtain the cleaned CMB B-mode map. Weobtain sky power ˆ σ B (cid:48) (cid:96) from the above cleaned CMB B-modemap. We use equation (14) to obtain weighted noise powersubtracted sky power ˆ σ B (cid:96) which is used to sample the theory C B (cid:96) . We use ten independent chains; each chain consists of10000 Gibbs steps. We discard the initial 50 samples for theburn-in period in each chain. In total we obtain 99500 sam-ples of C B (cid:96) and S . We perform this analysis on cases with 0 . .
05 tensor-to-scalar ratios.Using the samples obtained after applying our method, weforecast the proposed CMB space mission PICO’s ability toconstrain r in the presence of realistic lensing and foregroundcontributions. We simulate 200 different noise and fore-ground contaminated Gaussian random CMB B-mode realiza-tions as described in section (3) and apply Gibb’s ILC methodto obtain 99500 samples of C B (cid:96) and S for each of them. We usea set of samples { C B i (cid:96) } to compute the posterior distribution of the tensor-to-scalar ratio, r , and the amplitude of lensing, A lens , using the Blackwell-Rao estimator (Chu et al. 2005) foreach of the 200 cases. We use sampled C B i (cid:96) to obtain the best-fit value of the power spectrum for all the 200 simulationsand use them to study bias in the recovered power spectrum.We also obtain a mean map and study reconstruction error inrecovered CMB B-mode maps using our method. RESULTS
This section presents results obtained after applying ourmethod on the simulated foreground and noise-contaminatedCMB B-mode map at 6 frequency channels of proposed futureCMB mission PICO, with fiducial tensor-to-scalar ratios 0 . .
01. We present our method’s performance to recon-struct the CMB B-mode map, CMB B-mode angular powerspectrum, and power spectrum delensing in the following.
Cleaned Maps
In this subsection, we present the performance of ourmethod to reconstruct the CMB B-mode maps. In the Figure(2) and (3) we show pixel standard deviation maps obtainedusing 200 CMB signal reconstruction following our methodfor r = 0 .
05 and r = 0 .
01 respectively. The second last mapat the right bottom corner, labeled MEAN, in both the fig-ures shows the mean of all the 200 standard deviation mapsobtained from the 200 simulations. The last map at the rightbottom corner, labeled STDEV, of both the figures shows thestandard deviation maps obtained using the 200 standard de-viation maps from the 200 simulations. From the mean, stan-dard deviation map for both the cases, we find reconstructionbias along the galactic plane is ≤ − µ K . From the stan-dard deviation maps obtained using the 200 standard devia-tion maps, we find small variation of order ≤ − µ K in pixelreconstruction error from one simulation to another for bothcases of tensor-to-scalar ratios. In Figure (4), we show themean of 200 difference maps for both r values. We find amean map using 99500 samples of the map from a given sim-ulation, and subtract the input map to obtain the differencemap corresponding to the simulation. From the mean differ-ence maps, we find that the mean absolute pixel reconstruc-tion error is ≤ − µ K for both values of r, which indicatesaccurate signal reconstruction using our method. Angular Power Spectrum
In this subsection, we present our method’s performance toreconstruct the CMB B-mode angular power spectrum. Wepresent normalized densities of the Gibb’s samples of CMBtheoretical angular power spectrum from multipole 2 to 31,along with the input angular power spectrum (vertical blackdashed line) and best-fit angular power spectrum (vertical reddashed line) for a randomly chosen simulation seed 1, with r = 0 .
05 in Figure (5). In the figure the position of most ofthe histogram peeks agree well with the input sky angularpower spectrum. The deviation of the input angular powerspectrum from the peeks in some of the histograms is due topresence of detector noise in contaminated CMB frequencychannel maps. The plots in the Figure (5) confirms the ex-pected behavior of the C (cid:96) histograms at both low and highmultipoles. For the tensor-to-scalar ratio of 0.05, we presentin the top panel of Figure (6) mean over 200 simulations of in-put angular power spectrum and best-fit angular power spec-trum along with corresponding standard deviations to quan-tify the reconstruction error in CMB B-mode angular power F IG . 2.— Above figure shows standard deviation maps for seeds 1 to 23 obtained using CMB signal reconstruction method discussed in this article. The lastmap shown at the right bottom corner of the figure, labeled as STDEV, shows standard deviation map obtained using the 200 standard deviation maps from the200 simulations. The second last map shown in the last row of the figure, labeled as MEAN, shows the mean of all the 200 standard deviation maps obtainedfrom the 200 simulations. In all of the above maps, the standard deviation’s maximum value is well within the 10 µ K , which indicates accurate CMB signalreconstruction. The unit is in µ K thermodynamic.F IG . 3.— The figure shows 25 standard deviation maps for r = 0 .
01. The figure shows standard deviation maps for seeds 1 to 23 obtained using CMB signalreconstruction method discussed in this article. The last map shown at the right bottom corner of the figure, labeled as STDEV, shows standard deviation mapobtained using the 200 standard deviation maps from the 200 simulations. The second last map shown in the last row of the figure, labeled as MEAN, shows themean of all the 200 standard deviation maps obtained from the 200 simulations. All other maps in the above figure are standard deviation map for randomly chosenseeds out of the 200 simulations. In all these maps, the standard deviation’s maximum value is well within the 10 − µ K , which indicates accurate reconstructionsof the CMB signal. The unit is in µ K thermodynamic. new ILC method for weak CMB B-mode Signal 7 F IG . 4.— In the above figure, we show the mean of 200 difference maps for r values of 0.01 and 0.05 in the top and bottom panels, respectively. We finda mean map using 99500 map samples from a given simulation, and subtractthe input map to obtain the difference map corresponding to the simulation.From the above figure, we find that the mean over simulations of absolutepixel reconstruction error is ≤ − µ K for both values of r , indicating accu-rate map reconstruction using our method. The unit is in µ K thermodynamic. spectrum. We also plot in the bottom panel of the Figure (6),the mean over 200 simulations of difference between best-fitand input angular power spectrum along with correspondingstandard deviations to further quantify the reconstruction er-ror in recovered CMB B-mode angular power spectrum. Fromthe upper panel in the Figure (6) we find that the mean in-put and the mean best-fit power spectrum agree very well for r = 0 .
05. From the bottom panel of the Figure (6) we findthat the mean over simulations of absolute power reconstruc-tion error at each multipole is < × − µ K . Similarly wepresent normalized densities of Gibb’s samples of the theo-retical angular power spectrum from multipole 2 to 31, alongwith the input angular power spectrum (vertical black dashedline) and best-fit angular power spectrum (vertical red dashedline), for simulation seed 1, with r = 0 .
01 in Figure (7). Thebest-fit theoretical angular power spectrum estimate well theinput power spectrum. The Figure (7), confirms the expectedbehaviour of the angular power spectrum histograms at bothlow and high multipoles. In the Figure (8) for r = 0 .
01, wepresent in the top panel, mean over 200 simulations of inputand best-fit angular power spectrum, in the bottom panel, themean over 200 simulations of difference between the best-fitangular power spectrum and input angular power spectrumalong with corresponding standard deviations. From the up-per panel in Figure (8), we find that the mean best-fit powerspectrum has more power than the mean input power spec-trum at multipoles <
7. From the bottom panel of the Figure(8) we find that the mean over simulations of absolute power reconstruction error at each multipole is < × − µ K for r = 0 .
01. We plot in the Figure (9) the fractional bias ∆ C f b (cid:96) inrecovered angular power spectrum calculated using 200 best-fit angular power spectrum and input angular power spectrum,defined as: ∆ C f b (cid:96) = (cid:68) C best − fit (cid:96) (cid:69) − (cid:68) C input (cid:96) (cid:69)(cid:68) C input (cid:96) (cid:69) From the plot in the Figure (9), we find 2% to 3% morebias in the reconstructed power spectrum for simulated CMBB-mode maps with tensor-to-scalar ratio 0 .
01 than 0 .
05, in-dicating that our method does not have significant bias evenwhen r = 0 .
01. This shows that our method performs very wellin reconstructing the CMB B-mode angular power spectrumfor both the cases.
Reconstructing r
Using set of Gibbs samples C B i (cid:96) and Blackwell-Rao Estima-tor (Chu et al. 2005) we in a self-consistent Bayesian frame-work compute the joint posterior distribution P ( r , A lens ) of thetensor-to-scalar ratio r and the amplitude of lensing, A lens . Toestimate the cosmological parameter r and A lens we maximizethe likelihood L ( C B (cid:96) | C B th (cid:96) ) = exp (cid:32) − (cid:96) max (cid:88) (cid:96) =2 (2 (cid:96) + (cid:20) ln (cid:18) C B th (cid:96) C B (cid:96) (cid:19) + C B (cid:96) C B th (cid:96) − (cid:21)(cid:33) (26)where the model theoretical CMB B-mode power spectrum C B th (cid:96) is given by linear sum C B th (cid:96) = (cid:16) r . (cid:17) C T (cid:96) ( r = 0 . + A lens C L (cid:96) ( r = 0) (27)where C T (cid:96) ( r = 0 .
05) is the tensor B-mode power spectrum, fora tensor-to-scalar ratio r = 0 .
05, and C L (cid:96) ( r = 0) is the lensing-induced B-mode power spectrum. In order to estimate thejoint posterior distribution of r and A lens , P ( r , A lens ) we varyboth r and A lens and make use of the Blackwell-Rao approxi-mation P ( r , A lens ) ≈ N N (cid:88) i =1 L ( C B i (cid:96) | C B th (cid:96) ( r , A lens )) P prior ( A lens ) (28)where N is the total number of Gibbs samples of C B i (cid:96) used.For large N the Blackwell-Rao estimate becomes an exactapproximation of P ( r , A lens ) (Chu et al. 2005). In this work,we do not put any prior on A lens so that P prior ( A lens ) is a con-stant.In the following, we discuss our Bayesian method’s per-formance to delens the CMB B-mode angular power spec-trum and hence minimize the lensing contribution to the re-covered distribution of the tensor-to-scalar ratio P ( r ). Forfiducial tensor-to-scalar ratio 0.05 we plot in the Figure (10)normalized joint 2-D Blackwell-Rao posterior density esti-mates P ( r , A lens ) along with the normalized posterior distribu-tion P ( r ). We get the P ( r ) by slicing the joint 2-D Blackwell-Rao posterior density P ( r , A lens ), for maximum likelihood of A lens , using set of Gibb’s samples { C B , i (cid:96) } for each of the 200different simulations. Since the true value 0.05 is within 1 σ of the normalized posterior, we conclude that our method per-forms well in reconstructing the angular power spectrum for . . . .
81 0 3 6 9 00 . . . .
81 0 5 00 . . . .
81 1 2 00 . . . .
81 1 2 3 00 . . . .
81 1 2 00 . . . .
81 0 . . . . .
81 0 . . . . .
81 0 . . . . . .
81 0 . . . . .
81 0 . . . . . .
81 0 . . . . .
81 0 . . . . .
81 0 . . . . . . .
81 0 . . . . .
81 0 . . . . .
81 0 . . . . .
81 0 . . . . .
81 0 . . . . .
81 1 2 00 . . . .
81 1 00 . . . .
81 0 . . . . .
81 0 . . . . . . . . . . . . . . . . . . . . . . . .
81 0 . . . . . .
81 1 00 . . . .
81 1 1 . . . . .
81 1 2 00 . . . .
81 1 2 3 (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) F IG . 5.— In the above figure, we show normalized densities of the sampled CMB theoretical angular power spectrum obtained by Gibbs sampling for 2 to 31multipoles for simulation seed 1, with r = 0 .
05. The horizontal axis for each subplot represents (cid:96) ( (cid:96) + C (cid:96) / π in the unit of 10 − µ K . The above histogram givesus the best estimates of the theoretical C (cid:96) (vertical black dashed line) given the data. The vertical red dashed line is the value corresponding to input sky C (cid:96) .The position of most of the histogram peeks agree well with the input sky C (cid:96) . The deviation of the input C (cid:96) from the peeks in some of the histograms is due topresence of the residual detector noise along with the CMB. The above plots confirms the expected behavior of the C (cid:96) histograms at both low and high (cid:96) . -0.0003-0.0002-0.000100.00010.00020.0003 (cid:96) ( (cid:96) + ) C (cid:96) / π , µ K mean input clstddevmean best(cid:28)t cl (cid:96) ( (cid:96) + ) ∆ C (cid:96) / π , µ K Multipole, (cid:96) mean di(cid:27)rence clstddev F IG . 6.— For the tensor-to-scalar ratio r = 0 .
05, we plot the mean input angular power spectrum, mean best-fit angular power spectrum, and correspondingstandard deviation for 200 simulations in the top panel figure. We show the mean of 200 differences angular power spectrum and associated errors in the bottompanel figure. To get the difference angular power spectrum for a given simulation, we subtract the corresponding input power spectrum from the best-fit angularpower spectrum. We find that the mean input and the mean best-fit power spectrum agree very well from the top panel figure. In the bottom panel figure, themean over simulations of absolute power reconstruction error at each multipole is < × − µ K . The small absolute reconstruction error at each multipoleshows that our method performs very well in reconstructing the angular power spectrum. new ILC method for weak CMB B-mode Signal 9 . . . .
81 0 10 20 30 00 . . . .
81 10 20 00 . . . .
81 2 4 00 . . . .
81 3 6 00 . . . .
81 2 4 00 . . . .
81 1 200 . . . .
81 1 2 00 . . . .
81 1 00 . . . .
81 1 2 3 00 . . . .
81 1 2 00 . . . .
81 1 2 3 00 . . . .
81 1 2 300 . . . .
81 1 2 3 00 . . . .
81 1 2 00 . . . .
81 1 2 00 . . . .
81 1 2 00 . . . .
81 1 2 3 00 . . . .
81 1 2 300 . . . .
81 2 4 6 00 . . . .
81 2 4 00 . . . .
81 1 2 3 00 . . . .
81 2 4 00 . . . .
81 2 4 6 00 . . . .
81 2 4 600 . . . .
81 2 4 6 00 . . . .
81 2 4 6 00 . . . .
81 2 4 6 00 . . . .
81 2 4 00 . . . .
81 3 6 9 00 . . . .
81 6 12 (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) (cid:96) F IG . 7.— The figure shows normalized densities of the sampled CMB theoretical angular power spectrum obtained by Gibbs sampling for 2 to 31 multipolesfor simulation seed 1, with r = 0 .
01. The horizontal axis for each subplot represents (cid:96) ( (cid:96) + C (cid:96) / π in the unit of 10 − µ K . The vertical black dashed line is thebest-fit value, and the vertical red dashed line is the value corresponding to input C (cid:96) . The position of most of the histogram peeks agree well with the input sky C (cid:96) . The deviation of the input C (cid:96) from the peeks in some of the histograms is due to presence of the residual detector noise along with the CMB. The above plotsconfirms the expected behavior of the C (cid:96) histograms at both low and high (cid:96) . -0.0001-5e-0505e-050.0001 (cid:96) ( (cid:96) + ) C (cid:96) / π , µ K mean input clstddevmean best(cid:28)t cl (cid:96) ( (cid:96) + ) ∆ C (cid:96) / π , µ K Multipole, (cid:96) mean di(cid:27)rence clstddev F IG . 8.— In the top panel figure, we show the mean input angular power spectrum, mean best-fit angular power spectrum, and corresponding error bars. Weshow the mean of 200 difference angular power spectrum in the bottom panel figure and corresponding error bars. To get the difference angular power spectrumfor a given simulation, we subtract the corresponding input power spectrum from the best-fit angular power spectrum. We find the mean best-fit power spectrumfrom the top panel figure is slightly more than the mean input power spectrum at multipoles <
7. In the bottom panel figure, the mean over simulations of absolutepower reconstruction error at each multipole is < × − µ K . The small absolute reconstruction error at each multipole shows that our method performs verywell in reconstructing the angular power spectrum, and there is no significant bias in angular power spectrum reconstruction for r = 0 . -0.04-0.03-0.02-0.0100.010.020.030.040.050.06 r = 0 . -0.025-0.02-0.015-0.01-0.00500.0050.010.0150.020.025 r = 0 . (cid:0) (cid:10) C b e s t − f i t (cid:96) (cid:11) − (cid:10) C i n p u t (cid:96) (cid:11) (cid:1) / (cid:10) C i n p u t (cid:96) (cid:11) fractional biasstddevMultipole, (cid:96) fractional biasstddev F IG . 9.— In the top panel of this figure we show the fractional bias corresponding to reconstructed angular power spectrum at each multipole for fiducialtensor-to-scalar ratio of 0.01. In the bottom panel we present the same bias measure for fiducial tensor-to-scalar ratio of 0.05. From the above two plots we findthat we have 2% to 3% more positive bias for r = 0 .
01 than r = 0 .
05, indicating that our method do not have significant bias even for r = 0 .
01. This shows thatour method performs very well in reconstructing the CMB B-mode angular power spectrum for both the cases. r = 0 .
05 hence delensing the angular power spectrum. Sim-ilarly we plot normalized joint 2-D Blackwell-Rao posteriordensity P ( r , A lens ) in left panel and the posterior distribution P ( r ) in the right panel of the Figure (10) for fiducial tensor-to-scalar ratio 0.01. We find that the true value 0.01 is within 1 σ of the normalized posterior P ( r ), establishing that our methodalso performs well in reconstructing the angular power spec-trum for r = 0 .
01 and delensing the angular power spectrum.To show convergence of the posterior P ( r | D ) we show theproduct of the 200 posteriors in Figure (12) and Figure (13) asshaded gray band for fiducial tensor-to-scalar ratio 0.05 and0.01 respectively. Since the fiducial value of r in both thecases is within the corresponding gray band’s width, we con-clude that chains converge, and our method is correct. We alsoshow posterior for some randomly chosen simulations, whichwe normalize arbitrarily in each of the plots to fit on axes. DISCUSSIONS & CONCLUSION
We develop a new foreground model-independent approachto measure CMB B-mode signal and angular power spec-trum using simulated observations for proposed future gen-eration, PICO satellite mission in this work. Our new non-parametric method is useful since spatial and spectral varia-tions of the polarized foreground component may not be ac-curately known. In this article, we extend and improve theearlier reported Bayesian ILC method, to reconstruct weakCMB B-mode signals by introducing noise bias correctionsat two stages during the ILC weight estimation. We exten-sively test our new method’s performance to reconstruct the CMB B-mode sky signal and angular power spectrum and ob-tain the joint distribution of tensor to scalar ratio and lensingamplitude for two different values of r . The proposed futuregeneration CMB B-mode mission, like PICO, can break thepower spectrum degeneracy between primordial B-mode andlensing B-mode by detecting the reionization bump. Utilizingthis advantage, we further perform the delensing of the recov-ered CMB B-mode angular power spectrum. We use Gibb’ssamples of the CMB B-mode theoretical angular power spec-trum obtained from our method to separate the primordial andlensing B-mode contribution to the recovered distribution ofthe tensor-to-scalar ratio in a Bayesian manner and have quan-tified the performance of our method at two different tensor-to-scalar ratios.We summaries the findings of our method in the following: 1. From the mean, standard deviation maps for r = 0 . .
01, we find reconstruction bias along the galacticplane is ≤ − µ K . From the standard deviation oversimulations of standard deviation, we find small vari-ation of order ≤ − µ K in pixel reconstruction errorfrom one simulation to other for both 0 .
05 and 0 . ≤ − µ K ) for both the cases. In light of theabove, we conclude that our method accurately recon-structs the simulated primordial CMB B-mode sky forboth r = 0 .
05 and r = 0 .
01 cases. new ILC method for weak CMB B-mode Signal 11 A lens r . . . . P ( r ) r Fitting { r = 0 . , A lens } F IG . 10.— In the left panel figure, we plot normalized joint 2-D Blackwell-Rao posterior density estimates for A lens and r from 200 different simulations for thecase with fiducial tensor-to-scalar ratio 0.05. The color box shows the range of normalized likelihood. In the right panel figure, we show the posterior distributionof r , which we get by slicing the joint 2-D Blackwell-Rao posterior density on the left, for maximum likelihood of A lens . The vertical black dashed line representsthe fiducial tensor-to-scalar ratio of the input CMB B-mode maps, and the vertical red dashed line represents the mode of the r posterior. We can see in the rightpanel figure that the fiducial r value is within the 1 σ of the posterior, which indicates that we can achieve significant delensing of P ( r ) using our method. A lens r . . . . P ( r ) r Fitting { r = 0 . , A lens } F IG . 11.— The left panel figure shows normalized joint 2-D Blackwell-Rao posterior density estimates for A lens and r from 200 different simulations for thecase with fiducial tensor-to-scalar ratio 0.01. The color box shows the range of normalized likelihood. In the right panel figure, we show the posterior distributionof r , which we get by slicing the joint 2-D Blackwell-Rao posterior density on the left, for maximum likelihood of A lens . The vertical dashed line represents thefiducial tensor-to-scalar ratio of the input CMB B-mode maps, and the vertical red dashed line represents the r corresponding to the maximum likelihood for theposterior distribution. We can see in the right panel figure that the fiducial r value is within the 1 σ of the posterior, which indicates that we achieve significantdelensing of P ( r ) using our method. .
04 0 .
045 0 .
05 0 .
055 0 . r = 0 . P ( r | D ) r F IG . 12.— Blackwell-Rao posteriors from each of randomly selected 12different simulations for tensor-to-scalar ratio 0.05. The gray band is theproduct of the posteriors from all the 200 simulations with the tensor-to-scalarratio 0.05. We arbitrarily normalize the posteriors to fit on these axes. Thedark black dashed vertical line represents the simulated fiducial tensor-to-scalar ratio. We expect that the gray band covers the fiducial value of tensor-to-scalar ratio 0.05 to within its width, as is indeed the case showing that theposterior P ( r | D ) converges.
2. We find the mean input power spectrum and the meanbest-fit power spectrum agree very well for r = 0 . r = 0 .
01, there is slightly more meanpower in the reconstructed power spectrum at multi-poles <
7. Using fractional bias to quantify this positivebias, we find it to be only 2% to 3% more for the case .
008 0 .
009 0 .
01 0 .
011 0 .
012 0 . r = 0 . P ( r | D ) r F IG . 13.— In the above plot, we show Blackwell-Rao posteriors from ran-domly selected 15 simulations for the fiducial tensor-to-scalar ratio 0.01. Thegray band is the product of the posteriors from all the 200 simulations for thetensor-to-scalar ratio 0.01. We arbitrarily normalize the posteriors to fit onthese axes. The dark black dashed vertical line represents the simulated fidu-cial tensor-to-scalar ratio. We expect that the gray band covers the fiducialvalue of the tensor-to-scalar ratio 0.01 to within its width, as is indeed thecase. This shows that the posterior P ( r | D ) converges for r = 0 . with r = 0 .
01 than r = 0 .
05, which indicates that ourmethod does not have significant bias even for r = 0 . r and A lens , we find the fiducial tensor-to-scalar values to be within 1 σ of the recovered dis-tribution P ( r ) for both the cases. This shows that thesamples of the Gibb’s CMB B-mode theoretical angu-lar power spectrum, obtained using our method, givesan unbiased estimate of P ( r ) for both the cases. Thepower spectrum delensing method used in this articlecannot remove the lensing B-mode cosmic variance in-duced by E-modes. However, using our new method onthe foreground and noise-contaminated six PICO CMBB-mode channels, we can detect r with more than 9 σ and 8 σ significance if the true values of r were 0 . .
01 respectively at low resolutions without using any A lens prior.5. Our method does not explicitly require foreground model. Thus any error due to incorrect foreground B-mode model, does not bias our results.6. Our method is computationally fast, efficient, and ac-curate in delensing and detecting significant unbiaseddetection for levels of r ≥ − . ACKNOWLEDGMENTS