Introduction
The posting continues Clem, Walt and Sergi's work on small sky patch analysis. Questions I want to answer here are, how strong is the dust foreground generated by PySM in the context of \(r\)? Does it match the reality? In order to get a clear conclusion, I made \(r\)-equivalent maps from PySM and Planck 353GHz full-sky maps, and put them together for a comparison.
PySM Simulation and Planck Data
PySM (arXiv:1608.02841) is a publicly available numerical code used for the simulation of Galactic emission in intensity and polarization at microwave frequencies - Clem has a CMB-S4 posting checking dust decorrelation in these models. In this analysis there are 3 PySM simulation maps involved.
- d1_150: 150GHz simulation map generated by PySM d1 model (as in arXiv:1608.02841) with tophat bandpass.
- d4_150: 150GHz simulation map generated by PySM d4 model (as in arXiv:1608.02841) with tophat bandpass.
- d7_150: 150GHz simulation map generated by PySM d7 model (Hensley/Draine model) with tophat bandpass.
On the data side there are 5 Planck maps available.
- Pl353 Det Set Split: Planck 353GHz detector-set maps (nside = 2048). The detectors at this frequency are divided into two groups.
- Pl353 Year Split: Planck 353GHz yearly maps (nside = 2048). The data from odd/even years (year1+3/year2+4) of observations are used for the two maps.
- Pl353 Half-ring Split: Planck 353GHz half-ring maps (nside = 2048). The data from each pointing period is divided in halves.
- Pl353 Half-mission Split: Planck 353GHz half-mission maps (nside = 2048). The data from the first/second half-mission are used for the two maps.
- PIP XXX: a redigitalized \(r\)-equivalent map directly passed from figure 8 of PIP XXX paper (arXiv:1409.5738). According to the description of PIP XXX it was derived from Pl353 det set split maps.
Computation
This analysis is computed by numpy, scipy and healpy - I made Python code on top of Walt's groundwork. Here is the recipe (same as in PIP XXX) of how I compute an \(r\)-equivalent map from a full-sky input map.
- Work in Galactic coordinates (\(l\), \(b\)).
- Downgrade the input map to nside = 512 whenever it is larger than that.
- Define the resolution of the output map by a value of nside. In this analysis nside = 8 is used as in PIP XXX.
- On a single pixel of the output map
- Create a mask centering on this pixel - the mask chosen is an 11.3° radius circular disk (i.e. 400 sq. deg., ~1% of sky) tapered with a FWHM=2° Gaussian as in PIP XXX. It yields \(f_{\text{sky}}^{\text{eff}} = 0.0080\) as stated in PIP XXX.
- Subtract mean T/Q/U (of the area covered by the mask) from the input map. Then apply the mask to it.
- Run anafast on this masked input map to get \({{C_l}^{\prime}}^{TT,EE,BB}\)
- Planck data maps: for each input map, as it composes map A and map B coming from the split of data, there are three ways to make \({C_l}^{XX}\): \(C_l^{AA}, C_l^{BB}, C_l^{AB}\). In this analysis \(C_l^{AB}\) is chosen to remove noise bias appearing in auto spectrum.
- PySM input maps: as simulation input maps should have no noise, proceed to use the auto spectrum.
- Multiply factor \(l(l+1)/2\pi\) to \(C_l^{\prime}\) to get \(D_l^{\prime}\).
(Note: an old version of this recipe stated that \(C_l^{\prime}\) is first binned before it is transformed to \(D_l^{\prime}\). It was incorrect. The figures shown in this posting have always been produced by the results coming from binning of \(D_l^{\prime}\).) - \({D_l}^{\prime}\) is binned using top-hat binning on the intervals defined between \(l = 40,70,110,160,220,290,370\). This also yields bin centers \(l_c = 55, 90, 135, 190, 255, 330\).
- \(D_l^\prime\) is corrected for partial sky patch size (i.e. divided by the square of the integral of the apodization mask.) in order to scale it back to \(D_l\). No E/B purification.
- Removal of \(\Lambda\text{CDM}\) contribution
- Planck data maps: \(D_l^{XX} \rightarrow D_l^{XX} - {D_l}_{\Lambda\text{CDM}}^{XX}\), where \(XX = TT, EE\). As in PIP XXX only T and E are removed - the argument is that CMB B-mode power is negligible with respect to the dust polarization at 353GHz even at \(l=500\).
- PySM input maps: it seems like the input maps contain no CMB contribution at all, so this step is skipped.
- Error in \(D_l\)
- Planck data maps: calculate \(\Delta D_l\) by Xpol formula. The details are well-documented in equation (29) - (32) of Tristram et al., arXiv:astro-ph/0405575. Sample variance is quadratically subtracted from it to make \(\Delta D_l\) entirely due to instrument noise.
- PySM input maps: simulation has no noise, so the sample variance is taken as the \(D_l\) error.
- Fit \( (l_c, {D_{l_c}}^{BB}, \Delta {D_{l_c}}^{BB}) \) to the model function \({D_l}^{BB} = A^{BB}\big(\frac{l}{80}\big)^{\alpha^{BB}+2}\) where \(\alpha^{BB} = -2.42\) to get amplitude \(A^{BB}_\text{fit}\) and the corresponding \(\Delta A^{BB}_\text{fit}\).
- Frequency scaling
- Planck data maps: scale \(A^{BB}_\text{fit}\) obtained from 353GHz to 150GHz by using factor \(0.0395^2\) as in PIP XXX. \(\Delta A^{BB}_\text{fit}\) is scaled as well with uncertainty in the scaling factor added in quadrature.
- PySM input maps: the input maps are already at 150GHz - no scaling is needed.
- As PIP XXX states \(A^{BB} = 6.71 \times 10^{-2} \mu\text{K}^{2}_{\text{CMB}}\) when \(r=1\) and \(l=80\) (at 150GHz), and as this amplitude increases linearly with \(r\), calculate the equivalent \(r_d\) by \(r_d = A^{BB}_\text{fit} \, / \, (6.71 \times 10^{-2} \mu\text{K}^{2}_{\text{CMB}})\). So this is the \(r_d\) value of a single pixel. A corresponding \(\sigma(r_d)\) can be derived as well from \(\Delta A^{BB}_\text{fit}\).
Comparison
Figure 4.1 presents the smallfield \(|r|\)-equivalent maps in Galactic coordinates and the perspective of orthographic projection. I cut the highly saturated rim of each map to show \(|b| > 35^\circ\) only. The graticule used in these plots is defined by \((\Delta l, \Delta b) = (45^\circ, 20^\circ)\). Patches with negative \(r\) are marked by white crosses while the center of BICEP field (\(l=316.1^\circ,b=-58.3^\circ\)) is marked by a circular dot with color showing the magnitude coming from the measurement of BK14 - \(A_\text{dust@353GHz} = 4.300^{+1.200}_{-1.100} \mu\text{K}^2\) or \(r_d = 0.238^{+ 0.066}_{-0.061}\).
Note: The conversion is done by exploiting the modified blackbody law with \(\beta_d = 1.59\) and \(T_d = 19.6\text{K}\). We can scale \(A_\text{dust@353GHz}\) to \(A_\text{dust@150GHz}\) by multiplying the factor \(((150\text{GHz})^{\beta_d}B(\nu = 150\text{GHz}, T = T_d) \big/ (353\text{GHz})^{\beta_d}B(\nu = 353\text{GHz}, T = T_d))^2 = (0.0609)^2 \).
It is clear that d4_150 is an overestimate. d1_150 and d7_150 is closer to the reality but they still have higher \(r\). Moreover PySM models in general make the sky look more even - they cannot reproduce the spatial variation.
(Note: Pl353 Det Set Split is supposed to be the same as PIP XXX but they don't match. Walt has found this before and I again get the same conclusion.)
Updated on 20180503: I reran the analysis with the mask area increased to 1000 sq. deg. (~2.4% of sky). Now there is an extra option in all figures showing r/σ(r) from a larger area. Also the circular dot now has a color reflecting \(A_\text{dust}\) value from BK14.
- All negative r values disappear; 4 data set split r/σ(r) maps basically converge - det set split still shows a smaller \(r\) but they in general agree with each other. This is due to the increase in S/N ratio. Moreover adjacent patches become more correlated as common area is larger.
- Both PySM and Planck maps yield higher r. However the conclusion that PySM models overestimate \(r\) still holds.
- In the 400 sq. deg. case, \(A_\text{dust}\) from BK's measurement is basically the same as Planck's estimate; in the 1000 sq. deg. case, it is a bit smaller.
In figure 4.2 I make histograms comparing Planck data vs. PySM d1/d4/d7 models. As before I drop all data in \(|b| \leq 35^\circ\), since the emission directly coming from Galactic plane is out of our intereset. One can see Planck data cluster around smaller \(r\).
I also make the following graphs in order to show how \(r_d\) varies over Galactic latitude \(b\). In figure 4.3, "r_d vs. b" presents every \(r_d\) value at a particular \(b\), while "median(r_d) vs. b" takes median over Galactic longitude \(l\). The vertical dashed black line indicates \(b\) of BICEP field center and the black dot stands for BK14 measurement.
Products
Here are the nside = 8 HEALPix \(r\)-equivalent maps. In each FITS file it contains 4 columns: \(\big( r_d, \sigma(r_d), r_E, \sigma(r_E) \big)\). The latter two columns contain \(r\) and \(\sigma(r)\) derived from E-modes - under the assumption \(\alpha_{EE} = \alpha_{BB}\), I used \(C_l^{EE}\) in fitting to get amplitude \(A^{EE}_{\text{fit}}\), and calculated \(A^{BB}_{\text{fit}}\) by the relation \(A_{EE}/A_{BB} = 2\). They are relatively experimental so I did not show them in figure 4.1 .
- PySM \(r\)-equivalent maps:
- Planck \(r\)-equivalent maps:
And here are the nside = 8 HEALPix \(A^{BB}_\text{fit}\)-equivalent maps at 353GHz. They are in the unit of \( \mu\text{K}^{2}_{\text{CMB}} \).
- Planck \(A^{BB}_\text{fit}\)-equivalent maps:
Related postings and papers
- 2017 June 27: Checking dust decorrelation in models d1/d4/d7 and hipdt (Clem P.)
- The Python Sky Model: software for simulating the Galactic microwave sky (arXiv:1608.02841)
- Planck intermediate results. XXX. The angular power spectrum of polarized dust emission at intermediate and high Galactic latitudes (arXiv:1409.5738)
- Xspect, estimation of the angular power spectrum by computing cross-power spectra with analytical error bars (arXiv:1409.5738)