Impact of modelling foreground uncertainties

Impact of modelling foreground uncertainties on future CMB polarization satellite experiments

Carlos Hervías-Caimapo, Anna Bonaldi and Michael L. Brown
Jodrell Bank Centre for Astrophysics, School of Physics & Astronomy, University of Manchester, Oxford Road, Manchester M13 9PL, U.K.
SKA Organisation, Lower Withington Macclesfield, Cheshire SK11 9DL, U.K.
Accepted XXX. Received YYY; in original form ZZZ

We present an analysis of errors on the tensor-to-scalar ratio due to residual diffuse foregrounds. We use simulated observations of a CMB polarization satellite, the Cosmic Origins Explorer, using the specifications of the version proposed to ESA in 2010 (COrE). We construct a full pipeline from microwave sky maps to likelihood, using two models of diffuse Galactic foregrounds with different complexity, and assuming component separation with varying degrees of accuracy. Our pipeline uses a linear mixture (Generalized Least Squares) solution for component separation, and a hybrid approach for power spectrum estimation, with a Quadratic Maximum Likelihood estimator at low s and a pseudo- deconvolution at high s. In the likelihood for , we explore modelling foreground residuals as nuisance parameters. Our analysis aims at measuring the bias introduced in by mismodelling the foregrounds, and to determine what error is tolerable while still successfully detecting . We find that can be measured successfully even for a complex sky model and in the presence of foreground parameters error. However, the detection of is a lot more challenging, as inaccurate modelling of the foreground spectral properties may result in a biased measurement of . Once biases are eliminated, the total error on allows setting an upper limit rather than a detection, unless the uncertainties on the foreground spectral indices are very small, i.e. equal or better than 0.5% error for both dust and synchrotron. This emphasizes the need for pursuing research on component separation and foreground characterization in view of next-generation CMB polarization experiments.

cosmic background radiation – inflation – diffuse radiation – early Universe
pubyear: 2017pagerange: Impact of modelling foreground uncertainties on future CMB polarization satellite experimentsImpact of modelling foreground uncertainties on future CMB polarization satellite experiments

1 Introduction

The successful detection of primordial Cosmic Microwave Background (CMB) polarization -modes would confirm the inflationary paradigm, by probing the existence of gravitational waves that sets up the primordial tensor perturbations in the new born Universe during inflation. Also, this would allow us access into the energy scale of the very early Universe, GeV. See kamionkowski_2016; eco_paper_inflation for further details.

The astronomical community has put significant effort on the measurement of -modes. Planned future satellite and balloon experiments, such as CORE (core_paper), LiteBIRD (litebird_2014), PIXIE (pixie_2011), PRISM (prism_2014), LSPE (lspe_2012), and ground-based experiments, such as SPT (spt_2015), BICEP2-Keck (bicep2keck_2014), POLARBEAR (polarbear_2014), among others, aim at detecting the large-scale -mode polarization from the CMB in the near future. To accomplish this, the development of new detector technologies will allow an unprecedented high polarization sensitivity at microwave frequencies, capable of detecting , if indeed the final error is dominated by instrument noise.

It is worth pointing out that the signal could be much smaller, which would definitively test the limits of our instrumentation and abilities. Even if this is not the case, however, achieving the required sensitivity is not enough, because of the presence of bright diffuse Galactic and extra-galactic foregrounds that block our clean view into the CMB. Therefore, component separation techniques are developed to model and subtract these foregrounds, in order to obtain the cleanest possible CMB maps. The question is then, how accurate can we model and clean the foregrounds to the level required for measuring ?

Several forecasts of tensor-to-scalar ratio measurements including foreground residuals have been performed for different experiments (betoule_2009; armitage-caplan_2012; errard_2012; bonaldi_2014; remazeilles_2015; alonso_2016). In this work, we study how the error in the diffuse foregrounds component separation modelling propagates into the tensor-to-scalar ratio. Our approach is quite agnostic from the point of view of physical modelling of the Galactic emission, and it focuses on quantifying the bias on corresponding to some arbitrary modelling error levels. We also consider component separation and error mitigation techniques of different level of complexity.

This paper is organized as follows: In Section 2, we introduce the model we use to create simulated observations of the microwave sky by a representative future CMB satellite. In Section 3, we describe the pipeline we use to forecast the bias on the tensor-to-scalar ratio. In Section 4, we show the resulting bias for two different sky models, under different assumptions on component separation modelling complexity. Finally, in Section 5, we draw our conclusions.

Figure 1: Polarization intensity maps of the simulated sky at 105 GHz (left) and 555 GHz (right), for the sky model with variable spectral indices. Both maps are dominated by thermal dust emission. The maps of the model with spatially constant spectral indices look very similar.

2 Simulated observations

For our analysis, we use the specifications of the Cosmic Origins Explorer (COrE), as outlined in an earlier version of the experiment presented to ESA in 2010 (core_paper).

Recently the mission, renamed CORE, was formally submitted to ESA (CORE collaboration et al. in prep.) with different specifications (in particular more frequency channels, more tightly packed in the 60-600 GHz frequency range). However, in this work, we use the specifications from the earlier proposal to be representative of the capabilities of a future CMB polarization experiment. The frequencies, beam sizes, and sensitivities used in this work are listed in Table 1.

We perform our simulations using healpix (gorski_2005) maps with a resolution parameter of , corresponding to a pixel size of  arcmin. Some of the actual COrE bands have better resolution than the one allowed by such a pixel size, so we limit the band resolution to 7 arcmin in these frequencies, marked with in Table 1. This modification does not change our results appreciably because we focus on diffuse foreground components and primordial -modes, both dominant at low multipoles.

To simulate the full-sky observations of the microwave sky, we use the model presented in hervias_2016, based on the polarization results from the 2015 data release of Planck (planck_2015_10).

Band [GHz] 45 75 105 135 165 195 225 255 285 315 375 435 555 675 795
Beam FWHM [arcmin] 23.3 14.0 10.0 7.8
Noise [arcmin] 8.61 4.09 3.5 2.9 2.38 1.84 1.42 2.43 2.94 5.62 7.01 7.12 3.39 3.52 3.60
Table 1: COrE satellite specifications used in this work to simulate observations, taken from core_paper. As explained in the main text, the bands marked with a have better resolution than 7 arcmin, but have been simulated with a 7 arcmin pixel size ( Healpix maps) to limit the computational complexity of our analysis.

In this work, we consider three polarized sky components: CMB, thermal dust, and synchrotron. The main features of the sky model are as follows:

  • CMB It is a Gaussian realization of a theory power spectrum produced with CAMB (howlett_2012). The adopted cosmology is the following:  K, , , , ,  km/s/Mpc, , and . We include tensor perturbations with two different values of the tensor-to-scalar ratio and . Our simulation includes lensing -modes, generated from the CAMB power spectrum.

  • Thermal dust We use the dust polarization template described in hervias_2016, smoothed to . As a spectral law, we use a modified black body with a constant temperature of  K. For the spectral index, we use two models: one constant () and one spatially variable (based on the thermal dust spectral index map presented in planck_2015_10 and smoothed to ).

  • Synchrotron We use the syncrotron polarization template described in hervias_2016, smoothed to . We use a power law frequency scaling with a spectral index. Again we use either a constant (, as used in planck_2015_10), or a spatially variable (giardino_2002, having a resolution of ) spectral index.

Some maps for the sky model with spatially-variable indices are shown in Fig. 1.

For each model, we produce 100 sets of fits maps of the observed sky at each band. Each set has the frequency bands, resolution and white noise levels as specified in Table 1. Each set has a different CMB and white noise realization, but the same foreground components. We produce them with a healpix resolution parameter of . We also produce 100 low-resolution sets with . In this case, the modelled sky is produced with a resolution of across all bands, according with the larger size of the pixels. Although the beam size is not very well sampled by this pixel size, we have verified that, once both beam and pixel window function are deconvolved, the pipeline described in Sec. 3.2 yields an unbiased recovery of the CMB polarization power spectra. Since the seed used to create the CMB realization for a given set is the same always, the and CMB maps are the same realization, but with different resolution.

3 Methodology

In this section, we describe the various steps of our pipeline: component separation, power spectrum and likelihood estimation.

3.1 Component separation

To perform the component separation, we rely on the linear mixture model, stated as follows. The intensity of each foreground at a frequency band and in a line of sight can be expressed as , where is the corresponding assumed spectral law and would correspond to the template map of each foreground at a fixed arbitrary frequency. Then, the observed intensity is


where is the instrumental beam depending on the frequency channel , denotes convolution and is the instrumental noise. If the resolution of all frequency channels is the same, for each line of sight it is possible to rewrite the previous equation in matrix notation,


where is the mixing matrix, with dimensions (number of components) times (number of spectral bands). The vector now contains all the components convolved by the frequency-constant beam and contains all the data maps .

If the linear mixture models holds, it is possible to obtain an estimate of the components with a suitable linear mixture of the frequency channels, . If we know the mixing matrix, one possible solution is the Generalized Least Square solution (GLS), given by the matrix , where is the covariance matrix of the instrumental noise. This solution is unbiased in recovering , but retains a noise contribution. However, it minimizes the variance of the error when the sky signal is deterministic (delabrouille_2009). In practice, an estimate of the mixing matrix is typically calculated by parametrizing the spectral laws of the CMB and foreground components and by estimating the relevant parameters from the data. In this paper, we skip such estimation: we assume some error on the spectral parameters describing the true mixing matrix and propagate them through the full pipeline.

As stated above, one important assumption of the linear mixture model (at least when applied in pixel domain) is that the instrumental beam does not depend on frequency. This is not true in general, nor it is for COrE, as shown by Table 1. To overcome this problem, we pre-processed all maps by smoothing them with a Gaussian beam, thus equalizing their resolution to 23.3 arcmin (that of the lowest frequency channel, for the high-resolution sets ) or (the resolution sampled by the maps for the low-resolution sets).

3.2 Power spectra estimation

We estimated the polarization power spectra from the CMB maps with a hybrid approach: using a Quadratic Maximum Likelihood (QML) estimator at low () multipoles and a pseudo- estimator at the remaining intermediate and high multipoles. The QML estimator is optimal at low multipoles, and it is able to recover the reionization bump at . However, it gets very computationally demanding very quickly with increasing resolution. The pseudo- estimator is appropriate for high multipoles, which are unobtainable for the QML estimator, where it can recover the first acoustic peak at . This hybrid approach has been shown to be nearly optimal in the whole range and at the same time computationally feasible (e.g. Efstathiou2004; efstathiou_2006). The simulated observations at are used for estimating the pseudo- power spectra, while the low-resolution maps with are used for the QML estimator.

The QML method we use is based on tegmark_1997; tegmark_2001, see also Efstathiou_2004_2; Gruppuso2009. It works on pixel space, constructing an estimator based on the covariance matrices of the data. This method gives minimal error bars but it is very computationally demanding, since it requires operations of order , where is the number of pixels outside the mask.

The pseudo- deconvolution method we use is described in brown_2005 and brown_2009, which extended to polarization the technique proposed by hivon2002. This method uses a fast spherical harmonic transform to estimate the pseudo- spectra on the masked sky, and corrects them for the effect of the sky cut, noise and filtering with a deconvolution process. The output power spectrum is binned with bandpass window functions , and needs to be compared to a binned theory power spectrum


in the likelihood for .

Figure 2: Default Galactic mask used for the power spectrum estimation, retaining a fraction of the sky . Top: mask used for the QML power spectrum estimation; bottom: apodized mask used for the pseudo- power spectrum estimation.

3.2.1 Galactic mask

To exclude the foreground residual contamination due to the Galactic emission, we estimate the power spectrum outside a Galactic mask. The default mask is constructed using the dust and synchrotron polarization templates from planck_2015_10, smoothing them to a FWHM of , and masking every pixel with an intensity of  K or higher. We repeat this procedure for both and maps, and dust and synchrotron. We combine all of them to produce the final mask. For the pseudo- power spectrum estimation, it is beneficial to use an apodized mask, because sharp edges make the deconvolution kernel more complicated. Therefore, we apodize the mask by using the function


where is the distance between the pixel of interest and the closest masked pixel (with value 0), and is the distance scale of apodization (the scale in which the function goes from 1 to 0, in our case). The resulting apodized mask is shown in Fig. 2, bottom. The sky fraction retained is . The version of this mask, needed for the QML estimator, is not apodized, and has been constructed by rounding the mask and degrade to . This mask is shown in Fig. 2, top.

Figure 3: Left: polarization intensity maps of the foreground residuals (reconstructed-true CMB, top) compared to the thermal dust (middle) and synchrotron (bottom) maps reconstructed by the component separation. Right: full-sky power spectrum of foregrounds residuals compared to the full-sky power spectrum of the reconstructed thermal dust and synchrotron foregrounds. Notice the similar shape between the residuals and the thermal dust.

3.3 Cosmological parameters likelihood

We calculate the likelihood on the power spectra averaged over the 100 realizations of simulated observations, where we varied both the CMB and noise realizations. This effectively eliminates the cosmic variance bias, and only leaves the foreground residuals bias, which is of our interest. We define a standard Gaussian likelihood to calculate the posterior distribution of the tensor-to-scalar ratio. We define the as


where is the measured -mode bandpower at bin , is the binned -mode theory spectrum and is the inverse of the binned signal+noise covariance matrix. We construct from the low-multipole and high-multipole analysis, by joining at (with no overlap) the results from the QML and the pseudo- estimators.

The theory power spectrum is binned using the bandpass window functions at high multiple range and using a top hat function centered at each bin in the low multipole range. The theory power spectrum is calculated as


where is the primordial (scalar+tensor perturbations) power spectrum at a given , and is the weak gravitational lensed power spectrum, which we assume as known.

The covariance matrix is calculated using the 100 realizations signal+noise complete runs of the pipeline (including the component separation). Therefore, it accounts for cosmic and noise variance but also foreground residuals effects.

3.3.1 Modelling foreground residuals with nuisance parameters

The likelihood presented in equation (5) assumes that the measured power spectra contain only CMB and noise. In reality, there are also some foreground residuals, due to non-perfect component separation. We are now going to extend this likelihood to explicitly model a foreground residual contribution


where and are models for the power spectra of synchrotron and dust residuals, respectively, and the amplitudes and are two free nuisance parameters that can be varied, together with , and finally marginalized over. The need for adding one or both such extra parameters can be checked by seeing whether they improve the fit, by means of the reduced value.

In practice, a way to derive the foreground residual template models and is to assume that they are proportional to the dust and synchrotron power spectra. In our case, these can be computed from the foreground maps which are output of the component separation, as exemplified by Fig. 3.

In the analysis that follows, and are the binned power spectra of the thermal dust and synchrotron, respectively, reconstructed by the component separation. We process these maps through the same procedure we use for the reconstructed CMB, that is, the power spectra estimation with pseudo- for the high-resolution map and the QML estimator for the low-resolution map, under the same conditions.

Simulation run Sky model Component separation model Reference
Simple model Spatially constant , Constant s, with a 1,2,3% error 4.1; Table 3
Complex model Spatially variable , Constant s (the average of the true variable maps) 4.2.1; Table 4, top
Variable s, with a global error of 1% and 0.5% 4.2.2, Table 4, bottom
Table 2: Summary of the different runs performed in this work.

4 Results

We run the component separation pipeline, described in Section 3, for the two sky models, with constant and spatially-variable spectral indices. The summary of all the runs performed in this paper, together with descriptions and the referenced section where the results appear, is shown in Table 2.

Figure 4: Reconstructed power spectrum for the simulation with , constant foreground spectral indices and a % estimation error on both and . The reconstructed CMB (on 10 frequency bands, with  GHz) is biased (black circles). The modelled foreground residuals are shown as the diamonds and stars. The reconstructed CMB minus the modelled foregrounds residuals is shown as the green triangles, which can be compared with the theory power spectrum, shown as the grey curve. The multi-parameter likelihood yields
Figure 5: Left: tensor-to-scalar ratio likelihoods for the model with constant spectral indices and . The dashed grey curve shows the likelihood for the perfect knowledge of foreground spectral indices, centered in the correct . The dashed curves show the likelihood when an error of % is made on both spectral indices, using the entire frequency range (blue) and limited ( GHz) one (red). The solid curves show the results in the same cases when the multi-parameter likelihood is used. Right: same as the left panel, for the model with constant spectral indices and . We do not show the 15 bands, 1 parameter likelihood case with % error, since it is extremely biased, measuring .
Figure 6: Tensor-to-scalar ratio bias (estimated – true , top) and error (, bottom) for different cases of fixed constant errors on both spectral indices, for the simulation with constant spectral indices and . In the case of multi-parameter likelihoods, the empty symbols are the ones using only , and the filled symbols are the ones using both and .
Figure 7: Same as Fig. 6 for the simulation with constant spectral indices and . If the measured value is less than 2  away from , instead we plot the 95% upper limit, with an arrow down symbol.

4.1 Sky model with constant spectral indices (Simple model)

We use the model described in Section 2, where the foregrounds have spatially constant spectral indices ( and ) and we run the component separation assuming fixed errors on these spectral indices. These are , 2 and 3% errors on both and . As a reference, we also examined the case of perfect knowledge on the foregrounds, that is, 0% error on the spectral indices.

We consider two cases for the GLS reconstruction: a linear mixture of all the 15 COrE frequency bands, and one of only the lowest 10 bands, having  GHz. This is motivated by the fact that high-frequency bands are strongly contaminated by thermal dust, so including them in the CMB reconstruction increases the dust residuals for a given error on the dust spectral index. The drawback is an increase in the noise level, that needs to be weighted against the reduction in the foreground residuals. In any case, it is worth pointing out that the whole frequency range should be used in order to estimate the spectral indices before the GLS reconstruction, as this strategy in general achieves the smallest errors on and .

For all the assumed error cases, we calculate both the likelihood of equation (6), where the only parameter is , and the multi-parameter likelihood of equation (7), where we include either one or both of the foreground parameters and , depending on what is achieving the lowest reduced value.

As an example of the multi-parameter likelihood method, we show in Fig. 4 the power spectrum for the case with , % error on both spectral indices and using only the frequency bands  GHz. The reconstructed CMB (shown as the black circles) contains extra power because of the foreground residuals. However, we are able to model the foreground residuals, shown as diamonds for the thermal dust and as stars for the synchrotron. The 3-parameter model yields and unbiased value of despite the large foreground residuals present. The 1-parameter model gave the highly biased result of for the same case.

In Fig. 5 we show some example tensor-to-scalar ratio likelihoods for the simulations with (left) and (right). The grey curves show the perfect knowledge component separation using all the 15 frequency bands, which always yields an unbiased result ( and ). All the other curves assume a % error on both spectral indices. The blue curves correspond to a CMB reconstruction using 15 frequency bands, and a 1-parameter (dashed) and multi-parameter (solid) likelihood. The red curves are the same for a CMB reconstruction using only the first 10 frequency bands.

For both the and cases, the multi-parameter likelihood on the 10 frequency bands case allow removing the large bias corresponding to the % spectral index error. However, there is a degradation in the measured error . For the case, this does not allow a detection over 2  any more, but only corresponds to an upper limit.

We show the summary of all the results for the simulated observation with in top half of Table 3 and in Fig. 6. In the top panel, we show the measured bias (estimated minus true ); in the bottom panel, we show the width of the likelihood, . As expected, when we assume perfect knowledge of the foregrounds, the result is unbiased. However, when we introduce some error in the component separation, the likelihood is biased towards higher values of . If we adopt the multi-parameter likelihood instead of the 1-parameter one, the bias is either reduced or removed.

Limiting the frequency bands used in the CMB solution to  GHz is also effective in reducing the bias. In fact, a large fraction of the foreground residuals is introduced by the high frequency bands that are strongly dominated by thermal dust. By comparing the red stars (multi-parameter likelihood, full frequency range) to the green circles (1-parameter likelihood, limited frequency range), we see that they give similar biases for the same error on the spectral parameters. However, the green circles have smaller values than the red stars, this showing that, in this case, it is preferable to limit the bands used in the component separation than to introduce a multi-parameter likelihood. Even so, in some cases, when the bias is large (e.g. % spectral index error), both approaches must be used at the same time (shown by the yellow triangles).

The summary of all the cases for the simulations with is reported in the bottom half of Table 3 and shown in Fig. 7. It follows the same scheme from Fig. 6, for the same assumed component separation error cases. Getting an unbiased result is much more difficult in this case, due to the small value of . In particular, for errors in the spectral indices larger than %, we always need both the multi-parameter likelihood and the limited frequency range. We note that, for the perfect knowledge case, the value of is only just below the value allowing a 2  detection. Therefore, the multi-parameter likelihood increases and only allows for a 95% upper limit.

value , Using all 15 bands Using 10 bands  GHz
1-parameter multi-parameter 1-parameter multi-parameter
bias [] bias [] bias [] bias []
0% (0.1)
+1% (7.0) (1.1) (1.3)
-1% (6.6) (0.0) (1.0)
+2% (20) (3.6) (5.0)
-2% (18) (2.5) (3.7)
+3% (32) (6.2) (9.3) (0.5)
-3% (31) (4.6) (7.0)
0% 3.7
+1% (8.3) 8.0 (2.6) 5.3
-1% (7.3) 9.0 (2.3) 5.4
+2% (23) 10 (2.0) 15 (6.6) 6.7
-2% (21) 11 (2.6) 14 (6.4) 7.5
+3% (37) 13 (4.1) 18 (11) 7.9
-3% (34) 13 (5.4) 17 (11) 9.4
Table 3: Measured tensor-to-scalar ratio biases and values for runs on the simulation with spatially constant spectral indices. The component separation is modelled using the true spectral indices with a small error. In the bias columns, the bias expressed as number of is shown in parenthesis. The values with are 95% upper limits.

4.2 Sky model with spatially variable spectral indices (Complex model)

Now, we consider a more realistic model of the sky, where the spectral indices of the foreground components are spatially variable, as explained in Section 2. We simulate the measurement of the tensor-to-scalar ratio for two levels of component separation modelling complexity, as detailed in the two following sub-sections.

Figure 8: Histograms of the dust (left) and synchrotron (right) spectral index residuals (true – average spectral index) outside the Galactic mask. The true spectral indices are the ones in the sky model from hervias_2016. The true average is calculated on a map, while the true average is calculated on a map. The histograms are normalized so that they integrate to 1. The standard deviations of the distribution of spectral index residuals are equivalent to a 1.7% error for thermal dust and to a 3.5% error for synchrotron.

4.2.1 Modelling the component separation with spatially constant and

The first approach we adopt is to model the spatially variable spectral indices as a constant value across the sky. We set this value to the average of the true and maps outside the Galactic mask of Fig. 2. As such, these values ( and ) better represent most of the pixels used for the analysis. The histograms of the spectral index residuals, defined as the difference between the true indices (at a given pixel) and the true average value, for all the pixels in the sky outside the mask are shown in Fig. 8. The standard deviation of these residuals are 0.0253 for and 0.1074 for , which corresponds to a 1.7 % and 3.5 % error, respectively. These errors are qualitatively similar to the ones we considered in Section 4.1.

Figure 9: Same as Fig. 2, but now showing the optimized Galactic mask for the runs described in Section 4.2.1 and labelled as best in the top half of Table 4. The mask has .

We start by considering the 1-parameter and multi-parameter likelihood on the CMB reconstructed using the first 10 frequency bands. This was the best-performing case in the previous (constant spectral index) exercise. The results are reported in the top half of Table 4 as the base case. As we can see, this case is not good enough any more, because the bias on is still significant. We therefore proceeded to optimize the analysis to reduce the bias on , and obtained the results quoted in Table 4 as the best case. The modifications we introduced are:

  1. Secondary dust component: we use the idea presented in stolyarov_2005, which is that a component with spatially variable spectral dependence can be modelled as a series of components with constant spectral dependence, each one corresponding to a term in a Taylor expansion. In this way, the first-order thermal dust is the usual grey-body spectral law, and we add a second-order thermal dust component, whose spectrum is the derivative of that spectral law with respect to . We also explored the possibility to add a second synchrotron component to account for the variability of the synchrotron spectral index, but this achieved no significant improvement in our case.

  2. Optimization of the Galactic mask: we produced a new mask specifically optimized for -modes and tailored to the component separation approach we used. Specifically, we produced an estimate of foreground errors on the CMB map with a Monte Carlo (MC) approach, and then excluded all pixels for which this map was over some threshold. To derive the error map, we repeated the GLS CMB reconstruction 100 times by varying randomly the assumed spectral indices with Gaussian distributions (having and  K). We transformed each MC output CMB from to and finally computed the standard deviation of the 100 MC maps for each pixel. The optimized Galactic mask used in the best case analysis is shown in Fig. 9. This mask has , which is very similar to that of the mask used before.

Figure 10: Same as Fig. 5 for the simulation with spatially varying spectral indices and component separation assuming spatially constant spectral indices, for (left) and (right).

The results for this run are reported in the top of Table 4 and shown in Fig.  10. The left panel shows the measured likelihoods for the simulated observations with . The base case measurements (the blue curves) are biased, even when foreground residuals are modelled in the likelihood. However, the improvements we introduced to the analysis allow for an unbiased detection (shown by the solid red curve). This detection is only 2  away from .

Fig. 10, right, shows the same result likelihoods for . In this case, the results are all biased, even with the improvements in the best case (the red curves). This means that, for such a small value of , the systematic error we commit by neglecting the spatial variability of the spectral indices is too big to be compensated. In this case, the spatial variability needs to be modelled directly in the component separation, as we do in the next subsection.

value Case Measured values
1-parameter multi-parameter
bias bias
Base 0.00257
Best 0.00390
Base 0.00115
value , global error
0,0 %
1,1 %
1,0.5 %
0.5,1 %
Table 4: Measured tensor-to-scalar ratio biases and values for runs on the simulation with variable spectral indices. On the top half, we show the results for modelling with constant spectral indices in the component separation. In the bottom half, we show the results for modelling with the true spatially variable spectral indices with a small level of error in the component separation. In the bias columns, the bias expressed as number of is shown in parenthesis. The values with are 95% upper limits.

4.2.2 Modelling the component separation with spatially variable spectral indices

For the model with , which did not give an unbiased result in the previous subsection, we model the spatial variability of the spectral indices directly into the component separation. Indeed, most component separation methods are able to perform a local estimation of the foreground spectral properties, either pixel-by-pixel (e.g. commander, eriksen_2008, MIRAMARE, stompor2009), on sky patches (e.g. CCA, ricciardi2010) or by means of other kind of spatial localization (e.g., NILC delabrouille_2009_nilc; basak2013). As a drawback, estimation errors might be larger on a local estimation than on a global one, especially where foregrounds are weaker.

That is, in lines of sights where the (polarized) intensity is stronger, the error in the determination of spectral properties would in general be smaller, since there is a higher signal-to-noise ratio.

In order to model such error properties in our analysis, we investigate the spatial correlation of errors in and made with the commander algorithm. These simulated observations were produced with the Planck Sky Model (PSM, PSM2013) for the “Exploring Cosmic Origins with CORE” foregrounds paper (CORE collaboration et al. in prep.), and using a Galactic mask with . Although the instrumental specifications and the sky model are slightly different from what we use here, this allows us to derive the basic error properties that are needed for our modelling.

We find that the error on the synchrotron spectral index, , is consistent with a random distribution and it is not significantly correlated with the synchrotron polarized intensity. This is because of the frequency coverage of CORE, which does not include many synchrotron dominated channels, which means that the synchrotron estimation error is essentially noise-dominated. We then modelled the synchrotron spectral index as having a Gaussian distribution with standard deviation:


where is the error percentage we assume in our analysis and is the average synchrotron spectral index outside the Galactic mask.

We find that the thermal dust error, , is clearly anti-correlated with the polarized dust intensity. To model this property, we binned the pixels outside the Galactic mask into ranges of polarized intensity having roughly the same number of data points and fitted Gaussian density functions to the error distribution in each bin. The standard deviation as a function of is well modelled by a power-law:


where is the normalization corresponding to outside the Galactic mask, is the slope of the anti-correlation of with , and is the assumed error percentage.

To simulate the estimation of spatially variable spectral indices on our study, we start from the true input spectral indices maps, and add random error maps consistent with the error characterization of eqns. (8) and (9). We assume error levels of () and ().

We generated 100 realizations of both and at . We run our component separation pipeline on each of the 100 simulation sets, having a different CMB, noise and random spectral index error realization, on the  GHz limited frequency range. We used the spectral index maps directly for the low-resolution pipeline, and we upgraded them to for the high-resolution pipeline.

Figure 11: Tensor-to-scalar ratio likelihoods for a complex model (with spatially variable spectral indices) and modelling the component separation as spatially variable spectral indices. All the runs are made with the optimized mask shown in Fig. 9 and limiting the frequency coverage to  GHz. The case with 1% global error has a 1% error (standard deviation of error in pixels outside the Galactic mask) modelled with a spatially uniform random Gaussian error for and with a spatially correlated (following equation 9) random Gaussian error for . The case with 0.5% global error is analogous.

The likelihoods for the reconstructed power spectrum averaged over the 100 realizations for both the 1% and 0.5% global error cases, along with the case assuming perfect knowledge on the spectral indices, are shown in Fig. 11 and reported in the bottom part of Table 4. As usual, assuming a perfect knowledge (spectral indices errors equal to 0) yields an unbiased result, . With a 1% global error on both spectral indices, we measure a bias on of , and a for the 1-parameter likelihood. The multi-parameter likelihood yields a very small bias of , but with a 95% upper limit of .

With 0.5% global error, we still measure a small bias, of and an error of for the 1-parameter likelihood. Fitting with the foregrounds nuisance parameters is not well motivated, since the marginalized likelihoods for and are consistent with 0.

We have verified that the residual systematic error on is due in a greater proportion by thermal dust than to synchrotron residual contamination. In fact, the 1-parameter likelihood yields , and , respectively, if we consider a 0.5% global error only on and and we leave the other foreground at 1% global error.

5 Conclusions

We have performed an analysis on the tensor-to-scalar ratio bias produced by the mis-modelling of foreground spectral parameters, taking into account a realistic model of the sky and a full data analysis pipeline. We have considered two sky models: a very simple one, where the foregrounds (synchrotron and dust) have constant frequency spectra across the sky, and a more complex one, where the spectral dependence is spatially-varying. We modelled component separation strategies and likelihood estimations of increasing complexity. The main results of our analysis can be summarized as follows.

For , we obtain an unbiased estimation of for all simulations considered. The requirements on the accuracy of foregrounds modelling for component separation purposes are not too stringent (for example modelling spatially-varying foreground spectral indices as spatially-constant still gives a successful measurement).

Depending on the error level on the synchrotron and thermal dust spectral indices (from 1% to 3%), the best results may require exploiting a limited set of “cleaner” frequency maps to reconstruct the CMB. The use of all channels is anyway recommended to obtain an accurate estimation of the foreground spectral indices. We achieve significant improvements by explicitly modelling the synchrotron and dust foreground residuals in the likelihood, and marginalizing over foreground amplitude nuisance parameters. Furthermore, an important role is played by the Galactic mask, that needs to be optimized for the component separation method used and for -mode detection.

For and a simple sky model, using a suitable mask, modelling foreground residuals in the likelihood and limiting the frequency range for CMB reconstruction always yields an unbiased value. The error on often does not allow a detection but just an upper limit; this result is however conservative because the Gaussian likelihood we adopted is not optimal in the low-multipole regime.

When increasing the complexity of the sky, large modelling errors, such as approximating a spatially-varying spectral index with a constant, are not allowed in this case, as they give raise to biases on that are too large to be corrected for (at the likelihood level). Modelling the spatial variability of the spectral indices at the component separation level is required. For a global error of 0.5–1% in and , we obtain an unbiased detection/upper limit on . We show that the foreground residuals biasing the measurement of are due in greater proportion to thermal dust emission than synchrotron emission, due to the particular frequency coverage of COrE.

Such level of accuracy in the determination of foreground spectral parameters is very challenging, and motivates further research on component separation and foreground characterization. However, our analysis does not take into account polarization ancillary data that will become available, such as C-BASS (cbass_paper). Our method could also be used to further optimize the instrumental specifications of future CMB B-modes experiments.


We thank M. Remazeilles for his kind help with the commander code. We thank the CORE collaboration for allowing us to use data to calculate spectral parameter errors. CHC acknowledges the funding from Becas Chile/CONICYT. AB and MLB acknowledge support from the European Research Council under the EC FP7 grant number 280127. MLB also acknowledges support from an STFC Advanced/Halliday fellowship. We gratefully acknowledge the anonymous referee for useful suggestions that led to the improvement of this paper.


Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description