Exploring Cosmic Origins with CORE: B-mode Component Separation

TL;DR

This paper assesses the ability of the CORE satellite mission to detect primordial B-mode polarization in the CMB, demonstrating that with effective foreground control, it can measure tensor-to-scalar ratios down to about 5×10⁻³, despite challenges from lensing and foreground residuals.

Contribution

It provides a detailed simulation-based analysis of foreground separation and bias estimation for the CORE mission's B-mode detection capabilities, highlighting the importance of foreground control and delensing.

Findings

CORE can measure r ≈ 5×10⁻³ with minimal bias under ideal conditions.

Foreground residuals, especially polarized point sources, limit detection sensitivity at low r.

Delensing improves the significance of B-mode detection by reducing lensing contamination.

Abstract

We demonstrate that, for the baseline design of the CORE satellite mission, the polarized foregrounds can be controlled at the level required to allow the detection of the primordial cosmic microwave background (CMB) $B$-mode polarization with the desired accuracy at both reionization and recombination scales, for tensor-to-scalar ratio values of ${r\gtrsim 5\times 10^{-3}}$. We consider detailed sky simulations based on state-of-the-art CMB observations that consist of CMB polarization with $\tau=0.055$ and tensor-to-scalar values ranging from $r=10^{-2}$ to $10^{-3}$, Galactic synchrotron, and thermal dust polarization with variable spectral indices over the sky, polarized anomalous microwave emission, polarized infrared and radio sources, and gravitational lensing effects. Using both parametric and blind approaches, we perform full component separation and likelihood analysis of the simulations, allowing us to quantify both uncertainties and biases on the reconstructed primordial $B$-modes. Under the assumption of perfect control of lensing effects, CORE would measure an unbiased estimate of $r=\left(5 \pm 0.4\right)\times 10^{-3}$ after foreground cleaning. In the presence of both gravitational lensing effects and astrophysical foregrounds, the significance of the detection is lowered, with CORE achieving a $4\sigma$-measurement of $r=5\times 10^{-3}$ after foreground cleaning and $60$% delensing. For lower tensor-to-scalar ratios ($r=10^{-3}$) the overall uncertainty on $r$ is dominated by foreground residuals, not by the 40% residual of lensing cosmic variance. Moreover, the residual contribution of unprocessed polarized point-sources can be the dominant foreground contamination to primordial B-modes at this $r$ level, even on relatively large angular scales, $\ell \sim 50$. Finally, we report two sources of potential bias for the detection of the primordial $B$-modes.[abridged]

Tables

Table 1. Table 1 : Summary of simulated sky components.

Component	Emission law	Template
CMB	Blackbody derivative	$r={10}^{-2}$ (simulation #1)
		$r=5\times {10}^{-3}$ (simulation #2)
		$r={10}^{-3}$ (simulation #3)
		$r={10}^{-3}$, with lensing (simulation #4)
		$r={10}^{-3}$, with 40% lensing (simulation #5)
Synchrotron	Power-law $\nu {}^{{\beta }_s}$	WMAP $23$ GHz polarization maps
	Non-uniform $⟨{\beta }_s⟩=-3$	(Miville-Deschênes et al. 2008)
Thermal dust	Modified blackbody	Planck GNILC $353$ GHz map
	$\nu {}_{}{}^{{\beta }_d}B_{\nu }^{}\left(T_d\right)$	(Planck Collaboration XLVIII 2016)
	Non-uniform $⟨{\beta }_d⟩=1.6$	$[Q_{\nu },U_{\nu }]=f_d{g}_d{I}_{\nu }^{\mathrm{𝙶𝙽𝙸𝙻𝙲}}[\cos \left(2{\gamma }_d\right),\sin \left(2{\gamma }_d\right)]$
	Non-uniform $⟨T_d⟩=19.4$ K	$f_d=$ 15%, $⟨f_d{g}_d⟩=$ 5%
		$g_d$ and ${\gamma }_d$ coherent with synchrotron polarization
AME	Cold Neutral Medium	Thermal dust map rescaled by 0.91 K/K
		at $23$  GHz. Same polarization angles as
		thermal dust. Uniform 1% polarization fraction.
Point-sources	Four power-laws	Radio source surveys at $4.85$, $1.4$, $0.843$ GHz
		2.7% to 4.8% mean polarization fraction
	Modified blackbodies+free-free	IRAS ultra-compact H ii regions
	Modified blackbodies	IRAS infrared sources
		1% mean polarization fraction

Table 2. Table 2 : Set of simulations. Checkmarks indicate which components are included.

	CMB	dust	synchrotron	AME	sources	lensing
Simulation #1	$r={10}^{-2}$	$✓$	$✓$	-	-	-
Simulation #2	$r=5\times {10}^{-3}$	$✓$	$✓$	-	-	-
Simulation #2-bis	$r=2.5\times {10}^{-3}$	$✓$	$✓$	-	-	-
Simulation #3	$r={10}^{-3}$	$✓$	$✓$	-	-	-
Simulation #4	$r={10}^{-3}$	$✓$	$✓$	$✓$	$✓$	$✓$
Simulation #5	$r={10}^{-3}$	$✓$	$✓$	$✓$	$✓$	(40%)

Table 3. Table 3 : Instrumental specifications for the CORE  mission. For the purposes of generating reasonably sized simulations at N side = 2048 subscript 𝑁 side 2048 N_{\rm side}=2048 , the high-frequency observations ( ≥ 340 absent 340 \geq 340  GHz) have been simulated at 4 ′ superscript 4 ′ 4^{\prime} beam resolution instead of their native instrumental beam resolution.

Frequency	Beam	$Q$ and $U$ noise RMS
$\left[\mathrm{GHz}\right]$	$\left[\mathrm{arcmin}\right]$	$\left[\mu K.\mathrm{arcmin}\right]$
60	17.87	10.6
70	15.39	10.0
80	13.52	9.6
90	12.08	7.3
100	10.92	7.1
115	9.56	7.0
130	8.51	5.5
145	7.68	5.1
160	7.01	5.2
175	6.45	5.1
195	5.84	4.9
220	5.23	5.4
255	4.57	7.9
295	3.99	10.5
340	3.49 (4.0)	15.7
390	3.06 (4.0)	31.1
450	2.65 (4.0)	64.9
520	2.29 (4.0)	164.8
600	1.98 (4.0)	506.7

Table 4. Table 4 : Summary table of the tensor-to-scalar ratio estimates derived from simulated observations of the sky by CORE . First column : component separation methods. Second column : minimum multipole used in the likelihood estimation. Third column : maximum multipole used in the likelihood estimation. Fourth column : recovered value of the tensor-to-scalar ratio. Fifth column : uncertainty on r 𝑟 r . Sixth column : bias, defined as the difference between the recovered and input values of r 𝑟 r . For the cases where r = 10 − 3 𝑟 superscript 10 3 r=10^{-3} , the values in parenthesis correspond to estimates after correction for residual foregrounds in the likelihood (see text).

Simulation #1	$𝒓={\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}}$, dust, synchrotron
	${\ell }_{\min }$	${\ell }_{\max }$	$r\left[{10}^{-3}\right]$	$\sigma \left(r\right)\left[{10}^{-3}\right]$	$\left|r-r_{in}\right|/\sigma \left(r\right)$
Commander	2	47	$9.7$	$0.9$	$0.3$
NILC	48	349	$11.9$	$2.5$	$0.8$
SMICA	48	600	$11.7$	$1.0$	$1.7$
Commander $+$ NILC	2	349	$10.0$	$0.9$	$0.0$
Commander $+$ SMICA	2	600	$10.6$	$0.7$	$0.9$
Simulation #2	$𝒓=\mathrm{𝟓}\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟑}}$, dust, synchrotron
	${\ell }_{\min }$	${\ell }_{\max }$	$r\left[{10}^{-3}\right]$	$\sigma \left(r\right){10}^{-3}$	$\left|r-r_{in}\right|/\sigma \left(r\right)$
Commander	2	47	$4.6$	$0.4$	$1.0$
NILC	48	349	$7.2$	$2.5$	$0.9$
SMICA	48	600	$6.1$	$0.9$	$1.2$
Commander $+$ NILC	2	349	$4.7$	$0.4$	$0.7$
Commander $+$ SMICA	2	600	$4.9$	$0.4$	$0.2$
Simulation #3	$𝒓={\mathrm{𝟏𝟎}}^{-\mathrm{𝟑}}$, dust, synchrotron
	${\ell }_{\min }$	${\ell }_{\max }$	$r\left[{10}^{-3}\right]$	$\sigma \left(r\right)\left[{10}^{-3}\right]$	$\left|r-r_{in}\right|/\sigma \left(r\right)$
Commander	2	47	$1.3\left(1.0\right)$	$0.1\left(0.1\right)$	$3.0\left(0.0\right)$
NILC	48	349	$3.7$	$2.2$	$1.2$
SMICA	48	600	$2.1$	$0.8$	$1.4$
Commander $+$ NILC	2	349	$1.3\left(1.0\right)$	$0.1\left(0.1\right)$	$3.0\left(0.0\right)$
Commander $+$ SMICA	2	600	$1.3\left(1.0\right)$	$0.1\left(0.1\right)$	$3.0\left(0.0\right)$
Simulation #4	$𝒓={\mathrm{𝟏𝟎}}^{-\mathrm{𝟑}}$, dust, synchrotron,
	AME, sources, lensing
	${\ell }_{\min }$	${\ell }_{\max }$	$r\left[{10}^{-3}\right]$	$\sigma \left(r\right)\left[{10}^{-3}\right]$	$\left|r-r_{in}\right|/\sigma \left(r\right)$
Commander	2	47	$1.3\left(1.0\right)$	$0.5\left(0.6\right)$	$0.6\left(0.0\right)$
NILC	48	349	$9.1$	$3.7$	$2.2$
SMICA	48	600	$2.4$	$1.4$	$1.0$
Commander $+$ NILC	2	349	$1.4\left(1.2\right)$	$0.5\left(0.6\right)$	$0.8\left(0.3\right)$
Commander $+$ SMICA	2	600	$1.4\left(1.2\right)$	$0.5\left(0.6\right)$	$0.8\left(0.3\right)$
Simulation #5	$𝒓={\mathrm{𝟏𝟎}}^{-\mathrm{𝟑}}$, dust, synchrotron,
	AME, sources, 40% lensing
	${\ell }_{\min }$	${\ell }_{\max }$	$r\left[{10}^{-3}\right]$	$\sigma \left(r\right)\left[{10}^{-3}\right]$	$\left|r-r_{in}\right|/\sigma \left(r\right)$
Commander	2	47	$1.5\left(0.9\right)$	$0.5\left(0.5\right)$	$1.0\left(0.2\right)$
NILC	48	349	$8.5$	$3.4$	$2.2$
SMICA	48	600	$2.3$	$1.0$	$1.3$
Commander $+$ NILC	2	349	$1.6\left(1.1\right)$	$0.5\left(0.5\right)$	$1.2\left(0.2\right)$
Commander $+$ SMICA	2	600	$1.7\left(1.2\right)$	$0.4\left(0.5\right)$	$1.7\left(0.4\right)$

Equations

$$C_{\ell }^{BB}=C_{\ell }^{tensor}(r,\tau =0.055)+A_{lens}{C}_{\ell }^{lensing},$$

$$Q_{\nu }^{sync}$$

$$U_{\nu }^{sync}$$

$$Q_{\nu }^{dust}$$

$$U_{\nu }^{dust}$$

$$I_{\nu }^{\mathtt{GNILC}}={\tau }_{353}{\left(\frac{\nu }{353\mathrm{GHz}}\right)}^{{\beta }_d}{B}_{\nu }(T_d).$$

$$Q_{\nu }^{AME}$$

$$U_{\nu }^{AME}$$

$$I_{23\mathrm{GHz}}^{AME}$$

$$n^{FM}=\frac{n^{HM1}+n^{HM2}}{2},\text{mbox}with⟨n^{HM1},n^{HM2}⟩=0.$$

$$C_{\ell }^{thEE}(\tau )=\{\begin{array}{ll}{\left(\frac{\tau }{0.05}\right)}^2{C}_{\ell }^{EE}(\tau =0.05)& 2\le \ell \le 12\\ \frac{e^{-2\tau }}{e^{-2\times 0.05}}{C}_{\ell }^{EE}(\tau =0.05)& \ell >12\end{array},$$

$$P\left(\tau \right)\approx \frac{1}{N_G}\sum _{i=1}^{N_G}\mathcal{L}\left[{\hat{C}}_{\ell }^{EE(i)}|C_{\ell }^{thEE}\left(\tau \right)\right],$$

$$-2\ln \mathcal{L}\left[{\hat{C}}_{\ell }^{(i)}|C_{\ell }^{th}\right]=\sum _{\ell }(2\ell +1)\left[\ln \left(\frac{C_{\ell }^{th}}{{\hat{C}}_{\ell }^{(i)}}\right)+\frac{{\hat{C}}_{\ell }^{(i)}}{C_{\ell }^{th}}-1\right].$$

$$C_{\ell }^{thBB}(r,A_{lens})=\frac{r}{0.01}{C}_{\ell }^{tensor}(r=0.01)+A_{lens}{C}_{\ell }^{lensing}(r=0),$$

$$P\left(r,A_{lens}\right)\approx \frac{1}{N_G}\sum _{i=1}^{N_G}\mathcal{L}\left[{\hat{C}}_{\ell }^{BB(i)}|C_{\ell }^{thBB}\left(r,A_{lens}\right)\right],$$

$$C_q=\sum _{\ell \in [{\ell }_{\min }(q),{\ell }_{\max }(q)]}\frac{(2\ell +1)}{{\sum }_{\ell \in [{\ell }_{\min }(q),{\ell }_{\max }(q)]}(2\ell +1)}{C}_{\ell }.$$

$$f_l=\exp \left[-\frac{1}{2}{\left(\frac{l-l_{centre}}{100}\right)}^2\right],$$

$$\mathcal{L}(p)$$

$${\chi }^2(p)$$

$${\left(\frac{\nu }{{\nu }_0}\right)}^{{\beta }_d+C\log \left(\frac{\nu }{{\nu }_0}\right)}$$

$$C_{\ell }={\int }_0^{S_d}\frac{dN}{dS}{S}^{2}dS,$$

$$\mathbf{m}(\nu ,p)$$

$$+$$

$$\mathbf{s}=\left({\mathbf{s}}^{cmb},{\mathbf{s}}^{dust},{\mathbf{s}}^{sync}\right),\mathbf{\beta }=\left({\beta }_d,T_d,{\beta }_s\right),C_{\ell }=\left\{C_{\ell }^{EE},C_{\ell }^{BB}\right\}=⟨|{\mathbf{s}}_{\ell m}^{cmb}{|}^2⟩.$$

$$\displaystyle P\left(\boldsymbol{s},\boldsymbol{\beta},C_{\ell}\big{|}\boldsymbol{d}\right)\propto\mathcal{L}\left(\boldsymbol{d}\big{|}\boldsymbol{s},\boldsymbol{\beta},C_{\ell}\right)P\left(\boldsymbol{s},\boldsymbol{\beta},C_{\ell}\right),$$

$$\displaystyle\mathcal{L}\left(\boldsymbol{d}\big{|}\boldsymbol{s},\boldsymbol{\beta},C_{\ell}\right)\propto\exp\left\{{-1\over 2}\sum_{\nu}\left(\boldsymbol{d}(\nu,p)-\boldsymbol{m}(\nu,p)\right)^{t}{\boldsymbol{\rm N}}^{-1}\left(\boldsymbol{d}(\nu,p)-\boldsymbol{m}(\nu,p)\right)\right\},$$

$${\chi }^2(p)=\sum _{\nu }{\left(\frac{\mathbf{d}(\nu ,p)-\mathbf{m}(\nu ,p)}{{\sigma }_{\nu }(p)}\right)}^2,$$

$${\hat{\mathbf{s}}}^{(i+1)}$$

$${\hat{\mathbf{\beta }}}^{(i+1)}$$

$${\hat{C}}_{\ell }^{(i+1)}$$

$$\displaystyle P\left(\boldsymbol{s}\big{|}C_{\ell},\boldsymbol{d}\right)$$

$$\propto$$

$$\hat{\mathbf{s}}={\left({\mathrm{S}}^{-1}+{\mathrm{N}}^{-1}\right)}^{-1}{\mathrm{N}}^{-1}\mathbf{d}.$$

Full text

1 Introduction

The standard model of cosmology is based on the inflationary paradigm (Starobinsky 1980; Guth 1981; Linde 1982, 1983; Albrecht and Steinhardt 1982), yet direct observational evidence of an inflationary epoch remains elusive. Measurements of the cosmic microwave background (CMB) provide the cleanest experimental approach to address this issue, in particular since primordial gravitational waves, generated during an inflationary phase in the early Universe, induce a specific signature in its polarization properties.

The CMB polarization signal can be divided into even and odd parity $E$- and $B$-modes (Kamionkowski et al. 1997). The former are generated by both scalar and tensor perturbations, the latter by tensor modes only. $E$-mode polarization has been detected at high significance, as shown in studies of stacked fields centred on temperature hot and cold spots around which radial and tangential polarization patterns can be observed (Komatsu et al. 2011; Planck Collaboration XVI 2016). However, primordial $B$-modes, arising only from tensor perturbations which are intrinsically weaker than the $E$-mode generating scalar perturbations, are yet to be discovered.

We quantify constraints on $B$-modes in terms of the ratio, $r$, of the tensor fluctuations (gravitational waves) to scalar (density) fluctuations, evaluated at a given spatial wavenumber. The $B$-mode power spectrum has a peak at the horizon scale at recombination ($\ell\sim 90$) with an amplitude proportional to this value. Reionisation then introduces an additional peak at low-$\ell$ ($\ell\sim 10$) with an amplitude that depends on the optical depth of the Universe, $\tau$. Recent results from Planck (Planck Collaboration XLVI 2016) have determined a value for $\tau$ of $0.055\pm 0.009$, a decrease from the WMAP9 result of $0.089\pm 0.014$ (Hinshaw et al. 2013) and the previous Planck result of $0.078\pm 0.019$ (Planck Collaboration XIII 2016) obtained when combining a low-$\ell$ likelihood based on the 70 GHz polarization data and a high-$\ell$ temperature-based likelihood. This will have some implications for the possibility of detection of the primordial $B$-modes. Since the primordial $B$-mode power spectrum scales as $r\times\tau^{2}$ at the reionization scales $\ell\sim 10$, the current $15$% uncertainty on $\tau=0.055$ translates into a 30% uncertainty, and a possible shift, on the amplitude of the reionization bump of $B$-modes, and hence on $r$. Nevertheless, the determination of the $B$-mode power spectrum will provide a powerful probe of the physics of inflation. The current upper limit on the tensor-to-scalar ratio from the BICEP2 and Keck Array experiments (BICEP2 Collaboration et al. 2016) is $r\sim 0.07$.

Current measurements of the Galactic foreground emission (Planck Collaboration XXX 2016; Krachmalnicoff et al. 2016; Choi and Page 2015) imply that primordial $B$-modes will be sub-dominant relative to foregrounds on all angular scales and over all observational frequencies in the microwave regime. The detection of $B$-modes must, therefore, be regarded as a component separation problem. This issue has been addressed previously in the literature in the context of dedicated $B$-mode satellite experiments (Armitage-Caplan et al. 2012; Baccigalupi et al. 2004b; Betoule et al. 2009; Bonaldi and Ricciardi 2011; Dunkley et al. 2009a; Katayama and Komatsu 2011; Remazeilles et al. 2016; Errard et al. 2016; Stompor et al. 2016; Hervías-Caimapo et al. 2017). However, some of the conclusions are open to question due, in some cases, to the assumption of simplified foreground emission properties or the adoption of a higher $\tau$ value.

In this paper, one of a series of publications dedicated to the preparation of a future post-Planck CMB space mission, the Cosmic Origins Explorer (CORE, Delabrouille et al. 2017), we focus on evaluating the accuracy with which CORE can measure $r$ in the presence of foregrounds. Closely related papers include the companion papers on inflation (CORE Collaboration et al. 2016) and on cosmological parameters (Di Valentino et al. 2016). Other obstacles to the observation of $B$-modes, due to instrumental noise and parasitic systematic signal contributions, or the effects of gravitational lensing, are addressed in detail in Natoli et al. 2017 and Challinor et al. 2017 , respectively. Nevertheless, we do consider relevant issues related to gravitational lensing, which mixes the $E$ and $B$ polarization modes and creates a lensed $B$-mode spectrum peaking at $\ell\approx 1000$, in the context of determining $r$. Additionally, given the utilisation of the polarization $E$-modes for delensing purposes, we also briefly present the corresponding spectra, and infer the quality of the $E$-mode reconstruction by fitting the $\tau$ parameter.

Component separation will be the most critical step for measuring the primordial CMB $B$-mode signal at a level of $r\sim 10^{-3}$. A common approach to estimating the constraints on $r$ that a given experiment might achieve is via the Fisher forecasting formalism. However, the predicted uncertainties are usually optimistic, and do not capture potential biases in the recovery of the tensor-to-scalar ratio. Therefore we directly perform full component separation analysis on simulated CORE sky maps for several values of $r$, and include challenging simulations of foreground emission including contributions from synchrotron, dust, anomalous microwave emission (AME) and radio and infrared point sources. We then adopt an approach close in spirit to the analysis of actual real-world data. Specifically, various component separation approaches are applied to the simulated data, and Galactic diffuse and point source masks are inferred directly from the analysis, before evaluating $r$ via a likelihood method. This paper can be regarded as a follow-up to the comprehensive tretament in Leach et al. (2008) of component separation issues for intensity observations, but focussed instead on polarization.

The paper is organised as follows. In Sect. 2, we provide an overview of the important foregrounds that must be addressed when searching for primordial $B$-modes. In Sect. 3, we produce and describe challenging sky simulations for CORE. In Sect. 4, we perform a full component separation analysis for CORE, as it would be for real data analysis, on the sky simulations: point-source detection and pre-processing; component separation with parametric, blind and semi-blind methods; and likelihood estimation of the tensor-to-scalar ratio. Note that details of the component separation methods are given in Appendix A. Then, in Sect. 5, we present a hybrid likelihood analysis of simulations with $r=10^{-3}$, explicitly combining results from multiple component separation approaches. Section 6 highlights several important issues to be addressed in order to improve component separation approaches to future data sets. We conclude in Sect. 7, including a comparison to a forecasting approach for $r$ described in more detail in Appendix B.

2 Complexity of foregrounds

It has by now been established that the primordial $B$-mode CMB signal cannot be measured without correction for foreground emission. Here we provide a synthesis of the current understanding about the nature of such foregrounds. Since lensing induced $B$-modes are an effective foreground to the primordial $B$-mode signal, we also include a brief overview of their nature.

2.1 Diffuse Galactic emission

Our picture of the Galactic emission components in the microwave frequency range largely originates in the WMAP and Planck observations of the microwave sky from 23 to 857 GHz. The total intensity sky maps are consistent with an overall picture of the Galactic foreground that comprises four components (for a recent review, see e.g. Delabrouille and Cardoso 2007): synchrotron emission from relativistic cosmic ray electrons, free-free (thermal bremstrahlung) emission in the diffuse ionised medium, thermal (vibrational) emission from dust heated by the interstellar radiation field, and finally an anomalous microwave emission (AME) component strongly correlated spatially with the thermal dust emission but that exhibits a rising spectrum towards lower frequencies. The latter has been associated with rotational modes of excitation of small dust grains (so-called ’spinning dust’). A time-variable contribution on large angular scales from interplanetary dust (zodiacal light emission) has also been detected by Planck (Planck Collaboration XIV 2014), which may lead to systematic leakage from temperature to polarization at a level that might be non-negligible for very sensitive $B$-mode experiment, depending on the specifics of the scanning strategy.

In contrast to the situation for intensity where the foreground emission dominates over only 20% of the sky, the polarized flux at 20 GHz exceeds the level of CMB polarization over the full sky, and reveals the presence of large coherent emission features. Analysis of the WMAP and Planck data has demonstrated that the polarized Galactic emission is well-described by a simple two component model of the interstellar medium comprising synchrotron radiation and thermal dust emission. However, this picture is likely to become more complex as the sky is measured with increasing accuracy.

Synchrotron emission is produced by cosmic-ray electrons spiralling in the Galactic magnetic field. The measured synchrotron emission is dependent on the density of the relativistic electrons along a given line-of-sight, and approximately to the square of the plane-of-sky magnetic field component, and can be strongly polarized in the direction perpendicular to the Galactic magnetic field. It constitutes by far the most important component of the polarized foreground at low frequencies ($<50$ GHz). In detail, the observed polarized emission is seen to arise mainly in a narrow Galactic plane and well-defined filamentary structures – the loops and spurs well-known in total intensity measurements – that can extend over 100 degrees across the sky and be polarized at a level of $\sim$40% (see Vidal et al. 2015). However, away from these features, the polarization fraction remains relatively low, corresponding to values of less than $\sim$15% at high latitudes.

The synchrotron spectral energy distribution (SED) is typically modelled as a power law, often with some form of spectral curvature that is relevant for observations at microwave frequencies (Kogut et al. 2007). There is no precise determination of the spatial variation of the synchrotron spectral index, either in intensity or polarization. In the former case, this is, in part, due to the difficulty of separating the emission from free-free, AME and CMB, which may dominate the integrated emission in the 20$-$100 GHz range. It is also a consequence of the fact that the fidelity of current low-frequency data (e.g. Haslam et al. 1982; Reich et al. 2001) is not as good as for CMB data in general, although new measurements from experiments such as C-BASS (Irfan et al. 2015) should improve this situation. This uncertainty will have obvious implications for our attempts to model the diffuse emission in this paper. A further important observation is that the polarized synchrotron and dust contributions are spatially correlated (Page et al. 2007; Planck Collaboration XXII 2015) on large angular scales.

Polarized dust emission results from non-spherical grains that adopt a preferential orientation with the Galactic magnetic field and then emit thermal radiation along their longest axis. This will be perpendicular to the Galactic magnetic field and so the observed thermal dust emission is polarized in the same direction as the synchrotron emission. Studies of the Planck data have yielded a wealth of new knowledge about the nature of this emission over the entire sky.

Planck Collaboration XIX (2015) have shown that over a very small fraction of the sky the polarization fraction may reach $\sim$25%, although values are typically smaller ($\sim$12% at high latitudes) yet strongly variable and scale-dependent. This can be interpreted in terms of the structure of the turbulent part of the line-of-sight magnetic field and an associated depolarization effect.

There is no single theoretical emission law for dust, which is composed of many different populations of particles of matter. However, on average, an SED can be fitted to the observational data, generally in the form of a modified blackbody spectrum. Planck Collaboration XXII (2015) have determined the mean SED of dust emission in both intensity and polarization from WMAP and Planck data that is spatially correlated with the Planck 353 GHz emission. This is well fitted by a mean dust temperature of 19.6 K, and an opacity spectral index of $1.59\pm 0.02$ for polarization, slightly lower than that measured for total intensity. This modest evidence for a different frequency dependence in intensity and polarization may be connected to the variation of alignment efficiencies for various types of dust grains. It should also be noted that Meisner and Finkbeiner (2015); Finkbeiner et al. (1999) have demonstrated that a two-component dust model, with independent spectral indices and dust temperatures for the cold and hot components, provides a marginally better fit in intensity when combining the Planck and DIRBE data.

There is no precise estimate to date of the spatial variation of the polarized dust spectral index. However, Planck Collaboration L (2017) provides evidence for significant variations of the dust polarization SED at high Galactic latitude, larger than those measured for dust intensity. Moreover, Planck Collaboration L (2017) further demonstrates that the polarized dust emission may decorrelate across frequencies, because of the co-addition of different dust component spectra along the line-of-sight, in which case the use of fixed spectral indices across frequencies might be inadequate. Such variations can lead to an erroneous detection of primordial $B$-modes if not properly taken into account in any component separation analysis.

Although the synchrotron and thermal dust emission are clearly the dominant contributors to the diffuse polarized Galactic foreground emission, uncertainties in the current data may still allow other components, more evident in intensity measurements, to contribute at fainter levels.

If AME is solely due to spinning dust particles (Draine and Lazarian 1998), then we expect it to have a very low polarization percentage, $\lesssim$1%, with a specific level determined by the alignment efficiency of small grains in the interstellar magnetic field, and a polarization fraction that decreases with increasing frequency. Recent theoretical work (Draine and Hensley 2016) predicts that dissipative processes suppress the alignment of small grains contributing to the AME through rotational emission, such that negligible polarization ($\sim 10^{-6}$) will be observed at frequencies above 10 GHz, although such low levels still need to be confirmed empirically. However, AME might arise from other physical mechanisms. For example, Draine and Hensley (2013); Hoang and Lazarian (2016) have suggested that part of the observed AME emission may be due to magnetic dipolar emission, which would lead to a contribution with polarization perpendicular to that of thermal dust. Recent measurements place upper limits on the AME polarization at the few per cent level in individual clouds (Dickinson et al. 2011; López-Caraballo et al. 2011; Rubiño-Martín et al. 2012). Moreover, Planck Collaboration XXV (2016) have determined a 2$\sigma$ upper limit of 1.6% in the Perseus region, although setting such limits in other areas of the sky was hindered by significant synchrotron contamination. More recently, Génova-Santos et al. (2017) obtained the most stringent upper limit on AME polarization of $<0.22\%$ from the W43 molecular complex. Macellari et al. (2011) have analysed the AME polarization in diffuse regions of the sky, obtaining an upper limit of 5% for the diffuse AME polarization. Although it appears to be low, it is difficult to infer the true intrinsic AME polarization due to various potential depolarization effects such as averaging of polarization along the line-of-sight or within the telescope beam. Nevertheless, while the level at which the AME is polarized appears to be low, a failure to account for it could bias the measurement of $r$ in $B$-mode searches (Remazeilles et al. 2016).

Conversely, the contribution from free-free emission at high latitudes and away from bright H ii regions is negligible in polarization (Macellari et al. 2011), as expected theoretically from the randomness of Coulomb interactions in H ii regions.

In summary, diffuse Galactic foregrounds are potentially less complex for polarization studies because only a subset of the Galactic foreground emissions seen in intensity are significantly polarized. However, the component separation problem for $B$-mode polarization is more challenging because the CMB $B$-mode signal is itself intrinsically weak, especially if $r\sim 10^{-3}$ and $\tau=0.055$, compared to the foreground minimum at $\sim 70$ GHz (Fig. 1), so that a more precise understanding of the polarized foreground properties is required. It is worth noting that, contrary to expectations, an experiment with restricted low frequency coverage, e.g., without detector bands $<150$ GHz, still cannot avoid synchrotron contamination to $B$-modes: Fig. 1 demonstrates that at frequencies in excess of 200 GHz the synchrotron foreground has a similar spectral shape and amplitude to the primordial CMB $B$-mode signal for $r\lesssim 10^{-2}$. This may prevent multi-frequency component separation methods from disentangling the CMB and synchrotron $B$-modes in the absence of low-frequency observations acting as lever arms. Therefore, wide frequency coverage is essential to allow the accurate measurement of the primordial CMB $B$-modes, as will be provided by a CMB $B$-mode satellite mission like CORE.

2.2 Point sources

The two significant contributors at mm and sub-mm wavelengths are radio sources and dusty star-forming galaxies. Our understanding of both populations in the CORE spectral range has greatly improved in recent years, primarily thanks to Planck’s all sky surveys Planck Collaboration XXVI (2016a) and to the much deeper surveys over limited sky areas carried out by the Herschel satellite and by ground-based facilities such as the South Pole Telescope (SPT; Mocanu et al. 2013), the Atacama Cosmology Telescope (ACT; Marsden et al. 2014) and SCUBA-2 on the James Clerk Maxwell Telescope (JCMT; Geach et al. 2017).

The dominant radio source populations in the CORE frequency range are the compact flat- and inverted- spectrum ones, primarily blazars (BL Lac objects and flat-spectrum radio quasars). Observations with the Australia Telescope Compact Array (ATCA) in the frequency range between 4.5 and 40 GHz of 3 complete samples of such sources (for a total of 464 objects), carried out almost simultaneously with the first two Planck surveys, have shown that the spectra of most objects steepen above $\simeq 30\,$GHz, consistent with synchrotron emission becoming optically thin Massardi et al. (2016). The median high-frequency ($\nu\geq 70\,$GHz) slope was found to be in the range 0.6\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}\alpha\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}0.7 ($S_{\nu}\propto\nu^{-\alpha}$). However, individual sources show a broad variety of spectral shapes: flat, steep, upturning, peaked, inverted, downturning; (see also Planck Collaboration XIV 2011; Planck Collaboration XV 2011; Planck Collaboration XLV 2016). This complexity greatly complicates the removal of the point source contamination from CMB maps.

Extended, steep-spectrum radio sources are minor contributors at mm and sub-mm wavelengths. Nevertheless WMAP and Planck surveys have detected a few tens of these sources (López-Caniego et al. 2007; Massardi et al. 2009; Gold et al. 2011; Planck Collaboration XIII 2011; Planck Collaboration XXVI 2016a). A small fraction of them were resolved by Planck, in spite of its large beam and will also be resolved by CORE  complicating their removal from the CMB maps.

The local population of dusty star-forming galaxies was characterized by the InfraRed Astronomy Satellite (IRAS; Neugebauer et al. 1984). IRAS detected, in addition to relatively quiescent galaxies like the Milky Way, Luminous InfraRed Galaxies (LIRGs) with star-formation rates of tens to hundred $M_{\odot}$/yr and infrared luminosities in the range with $10^{11}\,L_{\odot}<L_{\rm IR}<10^{12}\,L_{\odot}$ and UltraLuminous Infrared Galaxies (ULIRGs) with $L_{\rm IR}>10^{12}\,L_{\odot}$ and up to $\geq 10^{13}\,L_{\odot}$, and star-formation rates of up to thousands $M_{\odot}$/yr. The dust emission of these galaxies is reasonably well described by a grey-body spectrum, which peaks at rest-frame wavelengths $\sim 100\,\mu$m. At mm and sub-mm wavelengths such a spectrum is approximated by $S(\nu)\propto\nu^{2+\beta}$ where $\beta$ is the dust emissivity index, which typically takes values in the range $1.5\leq\beta\leq 2$.

As mentioned above, the average frequency spectra of the two populations are widely different: the radio emission declines with increasing frequency while the dust emission steeply increases. This makes the crossover frequency between radio and dust emission components only weakly dependent on their relative intensities. Moreover, dust temperatures tend to be higher for distant high luminosity sources, partially compensating for the effect of redshift. As a consequence there is an abrupt change in the populations of bright sources above and below $\sim 1\,$mm: radio sources dominate at longer wavelengths, while in the sub-mm region dusty galaxies take over.

The wavelength at which the contribution of extragalactic sources is minimum is therefore shorter than that of minimum Galactic emission ($\simeq 5\,$mm, see Fig. 1). The power spectra of extragalactic sources are also very different than those of Galactic foregrounds. The Galactic dust power spectrum scales approximately as $C_{\ell}\propto\ell^{-2.7}$ or $\ell^{-2.8}$ for $\ell>110$ (Planck Collaboration XXX 2014) and the Galactic synchrotron power spectrum is similarly steep (La Porta et al. 2008). The point source power spectrum is much flatter. It is the sum of two components: Poisson fluctuations with $C_{\ell}=\hbox{constant}$ and clustering. However, the contribution of clustering to the angular power spectrum of radio sources is strongly diluted by the broadness of their luminosity function which mixes up, at any flux density level, sources distributed over a broad redshift range. As a consequence, the power spectrum is dominated by the Poisson term.

On the contrary, the power spectrum of dusty galaxies making up the Cosmic Infrared Background (CIB) is dominated by clustering for \ell\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}2000 (Planck Collaboration XXX 2014; De Zotti et al. 2015), while the Poisson contribution takes over on smaller scales (higher multipoles). Although the clustering power spectrum deviates from a simple power law, a reasonably good approximation is $C_{\ell}\propto\ell^{-1.2}$ (Planck Collaboration XXX 2014). The flatter point source power spectra compared to diffuse Galactic emissions imply that extragalactic sources are the main contaminants of CMB maps on small angular scales. This happens already for \ell\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}200 for \nu\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}100\,GHz, where the dominant population are radio sources and for \ell\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}1000–2000 at higher frequencies, where dusty galaxies dominate.

The most extensive study of the polarization properties of extragalactic radio sources at high radio frequencies was carried out in Massardi et al. (2013). These authors obtained polarization data for 180 extragalactic sources extracted from the Australia Telescope 20-GHz (AT20G) survey catalogue and observed with the Australia Telescope Compact Array (ATCA) during a dedicated, high-sensitivity run ($\sigma_{\rm p}\simeq 1\,$mJy). Complementing their data with polarization information for seven extended sources from the 9–yr Wilkinson Microwave Anisotropy Probe (WMAP) co-added maps at 23 GHz, they obtained a roughly $99\%$ complete sample of extragalactic sources brighter than $S_{20\rm GHz}=500\,$mJy at the selection epoch. The distribution of polarization degrees was found to be well described by a log-normal function with mean of $2.14\%$ and dispersion of $0.90\%$. Higher frequency surveys indicate that the distribution does not change appreciably, at least up to $\sim 40\,$GHz (cf. Battye et al. (2011) and Galluzzi et al. (2016)). The log-normal distribution of the polarization fractions of the radio sources with a mean of $\sim 3$% has now been confirmed up to $353$ GHz by Bonavera et al. (2017).

In the case of star-forming galaxies, the polarized emission in the CORE frequency range is dominated by dust at wavelengths \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}3\,mm. At longer wavelengths the synchrotron emission takes over; but at these wavelengths the extragalactic sky is dominated by radio sources, also in polarization.

Polarization properties of dusty galaxies as a whole at (sub-)mm wavelengths are almost completely unexplored. The only available information has come from SCUPOL, the polarimeter for SCUBA on the James Clerk Maxwell Telescope, that has provided polarization measurements at $850\,\mu$m for only two galaxies, M 82 (GreavesHolland2002) and M 87 (Matthews et al. 2009). However the global polarization degree has been published only for M 82 and is $\Pi=0.4\%$. Integrating the Planck dust polarization maps over a $20^{\circ}$ wide band centred on the Galactic plane, De Zotti et al. (2016) found an average value of the Stokes $Q$ parameter of about 2.7%. We may then expect a similar value for spiral galaxies seen edge-on. For a galaxy seen with an inclination angle $\theta$ the polarization degree is reduced by a factor $\cos(\theta)$. If all galaxies are about as polarized as ours, the average polarization fraction for unresolved galaxies, averaged over all possible orientations, should be about half of 2.7%, i.e. around 1.4%.

2.3 Lensing $B$-modes

Large-scale structures induce gravitational lensing in the CMB which mixes the $E$ and $B$ polarization modes (Blanchard and Schneider 1987; Cole and Efstathiou 1989; Bernardeau 1997; Zaldarriaga and Seljak 1998; Benabed et al. 2001; Challinor and Lewis 2005). The lensing $B$-mode power spectrum approximates that of white noise on large angular scales, peaks at $\ell\sim 1000$, and for $r\sim 0.07$, the current upper limit from the BICEP2 and Keck Array experiments (BICEP2 Collaboration et al. 2016), its amplitude is always larger than the primordial signal for scales smaller than the reionisation bump. Such a signal therefore acts as an effective foreground in the search for primordial $B$-modes. The recovery of the $B$-mode polarization may be attempted through a process called ‘delensing’. This requires an unlensed estimate of the $E$-mode signal and of the lensing potential (Hu and Okamoto 2002; Hirata and Seljak 2003). The latter can be derived from the CMB itself (Carron et al. 2017), or from alternative measures of large-scale structure, e.g., the CIB (Simard et al. 2015; Sherwin and Schmittfull 2015). In this context, it is worth noting that Larsen et al. (2016) have recently provided the first demonstration on Planck temperature data of CIB delensing, supporting its utility for lensing removal from high precision $B$-mode measurements.

In a companion paper (Challinor et al. 2017 , in prep.) it is shown that, for a CORE-like experiment, 60% of the lensing effect will be removed. We describe in a later section our approach to this signal and its treatment, and assess its impact on component separation and derived results on $r$.

3 Sky simulations

We produce detailed simulations of the polarized emission of the sky by using a modified version of the publicly released Planck Sky Model (PSM version 1.7.8, Delabrouille et al. 2013).111In our modified version of the PSM version 1.7.8, we have added the options to generate polarization for spinning dust and polarization for thermal dust with a modified blackbody spectrum as parametrized by the Planck GNILC dust model (Planck Collaboration XLVIII 2016). These additional models will be included in future releases of the PSM. The simulation is more challenging than has generally been considered in the literature to date.

The sky simulation consists of: (i) CMB $E$- and $B$-mode polarization with a low optical depth to reionization, $\tau=0.055$, and tensor-to-scalar ratios spanning the range $r=10^{-2}$ down to $10^{-3}$, including or not gravitational lensing effects; (ii) polarized synchrotron radiation with a power-law spectrum and variable spectral index over the sky; (iii) polarized thermal dust radiation with a modified blackbody spectrum and variable spectral index and temperature over the sky; (iv) polarized anomalous microwave emission (AME); and (v) infrared and radio polarized point-sources. The main characteristics of these components are summarized in Table 1. We analyse a set of 5 simulations, spanning different values of the tensor-to-scalar ratio and different amounts of gravitational lensing effects. The specific content of each simulation is given in Table 2. All of the simulated maps are provided in HEALPix format (Górski et al. 2005),222http://healpix.sourceforge.net with a pixel size defined by the $N_{\rm side}$ parameter, here set to a value of $2048$. A lower resolution set of simulations generated directly at $N_{\rm side}=16$ are also provided for the Commander analysis (see Sects. 4.2 and 6.2).

3.1 CMB

By using the Boltzmann solver CAMB (Lewis et al. 2000), we generate both lensed and unlensed $E$- and $B$-mode CMB angular power spectra from a $\Lambda$CDM+$r$ cosmology, an optical depth to reionization, $\tau=0.055$, motivated by the latest Planck results (Planck Collaboration XLVI 2016) and a tensor-to-scalar ratio varying from $r=10^{-2}$ down to $r=10^{-3}$, which is the ambitious detection goal for the CORE space mission. The other $\Lambda$CDM cosmological parameters are set to the Planck best-fit values from Planck Collaboration I (2016). The CMB $B$-mode angular power spectrum, $C_{\ell}^{BB}$, generated by CAMB therefore is the combination of a pure tensor power spectrum and a lens-induced power spectrum:

[TABLE]

where $r$ is either set to $10^{-2}$, $5\times 10^{-3}$, $2.5\times 10^{-3}$, or $10^{-3}$, and $A_{lens}$ is either set to 0 (unlensed), 1 (no delensing), or 0.4 (60% delensing). The delensed case is idealized – when analyzing observations of the real sky, delensing would be applied post-component separation (Carron et al. 2017), but such a treatment is beyond the scope of this work. Instead, we assume a scale-independent delensing efficiency, and simply rescale the full lensing contribution. The application of component separation methodologies to such sky realizations does not affect the effectivness of the foreground removal, and allows the generation of signal covariance matrices for the likelihood analysis described later in Sect. 4.3.

Gaussian random realizations of the CMB Stokes $Q$ and $U$ polarization components can then be simulated with the appropriate $E$- and $B$-mode power spectra using the HEALPix routine synfast. The corresponding $E$- and $B$-mode maps are then generated from the spherical harmonic transforms of the $Q$ and $U$ components computed by the HEALPix anafast routine, where the resulting $a_{\ell m}^{E}$, and $a_{\ell m}^{B}$ (pseudo)scalar coefficients are transformed to full-sky maps using synfast. In this paper, we consider only a single CMB realization per simulation.

The lensed CMB $Q$ polarization map with $r=10^{-3}$, smoothed to two degree resolution (FWHM) for illustrative purposes, is shown in the top left panel of Fig. 2. The CMB polarization $Q$ and $U$ maps are scaled across the CORE frequency channels through the derivative of a blackbody spectrum that is achromatic in thermodynamic temperature units. The component of interest, i.e. the primordial CMB $B$-mode polarization map, is shown for $r=10^{-3}$ and $\tau=0.055$ in the top left panel of Fig. 3, while the lensed CMB $B$-mode polarization map is shown in the top right panel of Fig. 3. Note that gravitational lensing effects add significant small-scale noise to the anisotropies of the primordial CMB $B$-mode polarization, for which the bulk of the cosmological signal is on the degree scale and larger.

3.2 Synchrotron

The Galactic synchrotron radiation is simulated by extrapolating the WMAP 23 GHz polarization maps, $Q_{23\,GHz}$ and $U_{23\,GHz}$, to CORE frequencies through a power-law frequency dependence

[TABLE]

with an average spectral index, $\langle\beta_{s}\rangle=-3$, and including spatial variations over the sky. The variable spectral index map was estimated by fitting a power-law to the Haslam et al. 408 MHz map and a WMAP 23 GHz synchrotron map derived using polarization data (Miville-Deschênes et al. 2008). The Stokes $Q$ map of the synchrotron polarization component at $23$ GHz and the synchrotron spectral index map are shown in Fig. 2. The synchrotron $B$-mode map at $60$ GHz is shown in the bottom left panel of Fig. 3.

To date there is still no concensus in the literature as to an optimal estimate of the synchrotron spectral indices (see, e.g., Dickinson et al. 2009, to justify the choice of templates). The characterization of the synchrotron spectral indices is problematic due to three main reasons: the difficulty in separating synchrotron from free-free emission (and AME at higher frequencies) in intensity; uncertainties in modelling the spectral shape, which is not well described by a single power-law over the wide frequency range considered; and the quality of the low-frequency 408 MHz data arising from significant variations of the calibration with angular scale, that can result in artificial variations in inferred spectral index maps.

For our simulations we elected to use the spectral index map estimated in each pixel by Miville-Deschênes et al. (2008) from the WMAP 23 GHz polarization map and the Haslam et al. 408 MHz intensity data, as currently implemented in the PSM software. This template has a representative mean value of $-3$ which is close to typical values observed in the literature at CMB frequencies (Davies et al. 1996; Planck Collaboration X 2016; Kogut et al. 2007; Dickinson et al. 2009; Miville-Deschênes et al. 2008; Bennett et al. 2013). The exact choice of the spectral index template is not critical for the simulations, as long as a reliable, physically-motivated SED (here a power-law with a mean spectral index of $-3$) is used to scale the synchrotron emission across frequencies.

3.3 Thermal dust

In this study, we focus on the spectral variations in the dust emission over the sky as the main complexity of the dust foreground, and postpone the inclusion of other potentially important effects, such as the frequency decorrelation333The spectral index of the Galactic dust emission varies along a given line-of-sight and with frequency due to the emission from multiple components, so that the actual spectrum is not a power-law. noted in Sect. 2, for future investigations. We therefore consider only a single modified blackbody dust component in the simulations. However, we note that mismodelling the dust emission in a parametric component separation method, e.g., by parametrising the emission with a single modified blackbody when two are required to fit the data accurately, can strongly bias the estimate of the tensor-to-scalar ratio (Remazeilles et al. 2016). The accurate characterization of the spectral properties of dust is essential for $B$-mode foreground studies and active research is being pursued in this field.

Here, the polarization maps of the Galactic thermal dust radiation are generated from the intensity map of the Planck GNILC 2016 dust model (Planck Collaboration XLVIII 2016), from which the CIB fluctuations have been removed.

[TABLE]

where the dust intensity map $I^{\tt GNILC}_{\nu}$ is scaled to the CORE frequencies through a modified blackbody spectrum,

[TABLE]

Here, $\tau_{353}$, $\beta_{d}$, and $T_{d}$ are respectively the Planck GNILC dust optical depth map at $353$ GHz, the Planck GNILC dust emissivity map, and the Planck GNILC dust temperature map that were derived in Planck Collaboration XLVIII (2016). The dust emissivity and temperature are both variable over the sky with average values, $\langle\beta_{d}\rangle=1.6$ and $\langle T_{d}\rangle=19.4$ K, respectively. The maps of dust emissivity and dust temperature are shown in Fig. 2. The dust $B$-mode map at $600$ GHz is shown in the bottom right panel of Fig. 3.

We make the assumption that the dust polarization angle map, $\gamma_{d}$, and the geometric depolarization map, $g_{d}$, due to the specific magnetic field configuration (Miville-Deschênes et al. 2008), are coherent with those of the polarized synchrotron model. The dust polarization fraction $f_{d}$ is set to 15 % on the sky, which, after modulation with the geometric depolarization factor, gives an overall polarization fraction fraction of $f_{d}g_{d}\sim 5$ % on average444In fact, the dust polarization fraction adopted in the public version of the PSM corresponds to about half of the currently accepted value. , with spatial variations over the sky. The Stokes $Q$ polarization map of thermal dust is shown in the right panel of the second row of Fig. 2.

3.4 AME

For sensitive CMB experiments, AME, even with a low polarization fraction, may be a relevant low-frequency foreground to the primordial CMB $B$-modes, especially for low values of the tensor-to-scalar ratio (Remazeilles et al. 2016). Therefore, we include in the sky simulation an AME component with a uniform $\pi=1$ % polarization fraction over the sky.

[TABLE]

where the AME intensity map, $I^{AME}_{\nu}$, is the Planck thermal dust intensity map at 353 GHz, $I^{\tt GNILC}_{353}$ (Planck Collaboration XLVIII 2016), but rescaled by a factor $0.91$ K/K at $23$ GHz using the correlation coefficient between AME and thermal dust measured by Planck Collaboration XXV (2016), and extrapolated to CORE frequencies from the $23$ GHz value by assuming a Cold Neutral Medium (CNM) for modelling the emission law (Ali-Haïmoud et al. 2009; Draine and Lazarian 1998). Because of the correlation between AME and thermal dust, we choose the AME polarization angles, $\gamma_{d}$, to be identical to those of thermal dust. The Stokes $Q$ polarization map for AME at $60$ GHz is shown in the bottom left panel of Fig. 2.

3.5 Point-sources

While polarized compact extragalactic sources are a negligible foreground for CMB $B$-modes on very large angular scales near the reionization peak ($\ell\lesssim 12$), they are expected to be the dominant foreground for $r=10^{-3}$ once delensing has been applied to the data, from the recombination peak to smaller angular scales ($\ell>47$, Curto et al. 2013). Therefore, we include in the sky simulation both radio and infrared extragalactic sources in polarization.

Radio sources are taken from radio surveys at $4.85$, $1.4$, and $0.843$ GHz (Delabrouille et al. 2013), and extrapolated to CORE frequencies assuming four kinds of power-laws for assigning the radio sources to either a steep- or flat-spectrum class. The polarization degree of the radio sources is randomly selected from the observed sample of flat and steep radio sources (Ricci et al. 2004), so that the polarization fraction of the radio sources in our PSM simulation is 2.7% on average for flat sources and 4.8% on average for steep sources. The $Q$ polarization map of strong radio sources at $60$ GHz is shown in the top right panel of Fig. 2.

Strong and faint infrared sources are taken from the IRAS point-source catalogue (Beichman et al. 1988; Moshir et al. 1992) and extrapolated to CORE frequencies by assuming modified blackbody spectra (see Delabrouille et al. (2013) for more details). The polarization fraction of the infrared sources is distributed around an average value of 1% through a chi-square distribution. We also include ultra compact H ii regions extracted from IRAS, which we extrapolate to CORE frequencies by assuming both modified blackbody spectra and power-law spectra due to free-free emission.

3.6 CORE instrumental specifications

The proposed CORE space mission can observe the polarized sky emission in 19 frequency bands ranging from $60$ to $600$ GHz. The goal is that this wide frequency coverage will provide lever arms that allow non-trivial foregrounds at low and high frequencies to be modelled adequately. The large number of frequency channels is also essential for component separation when facing multiple degrees of freedom for foregrounds, e.g., decorrelation effects that may result from multi-layer dust emission (Planck Collaboration L 2017), spectral index curvature, or as-yet undiscovered foregrounds.

The instrumental specifications of the CORE space mission are summarized in Table 3. The optical performance of CORE will allow high-resolution observations with FWHM $<10^{\prime}$ over the primary CMB frequency channels and a few arcminute resolution at high frequencies. The high resolution of the CORE observations will play an important role in the correction of the primordial CMB $B$-modes for gravitational lensing effects. CORE will have unprecedented sensitivity with detector noise levels of order $\sim 5\,\mu$K.arcmin in polarization observations at the primary CMB frequencies.

We convolve the component maps of our sky simulation with Gaussian beams with the FWHMs given in Table 3 for each frequency channel. Note that we have limited the high-frequency observations ($\geq 340$ GHz) to $4^{\prime}$ beam resolution instead of the native instrumental beam resolution in order to avoid an oversized data set. This does not impact the results of this study in which we are interested in detecting CMB $B$-modes on large angular scales. Sky maps at each frequency are obtained by co-adding the beam-convolved component maps.

For the purposes of this work it is sufficient to assume that the noise is Gaussian and white, and uncorrelated between the $Q$ and $U$ Stokes parameters. We also ignore any variation of the noise properties across the pixeized sky as would be introduced by a realistic scanning strategy. We do, however, consider the division of the data into ‘half-mission’ surveys. These are also idealized, splitting the full-mission data into two equal parts. We then simulate two distinct realizations of white noise for the half-mission $Q$ and $U$ maps at a given frequency by using the noise RMS values listed in Table 3 multiplied by a factor of $\sqrt{2}$. The two sets of noise realizations are then co-added to the same sky realisation to generate two half-mission surveys which have uncorrelated noise properties. The resulting two sets of observation maps in the $19$ frequency bands are referred to as the CORE half-mission 1 (HM1) and half-mission survey 2 (HM2) surveys. The corresponding full-mission survey (FM) is formed simply by adding to the sky maps at each frequency the full-mission noise maps, $n^{FM}$, that are related to the half-mission noise maps, $n^{HM1}$ and $n^{HM2}$, by:

[TABLE]

The component separation is performed on the full-mission simulation. The appropriately cleaned half-mission simulations can be used to compute the CMB power spectrum free from noise bias via cross-spectral estimators.

A companion paper (Natoli et al. 2017 ) considers in detail more realistic instrumental simulations, addressing topics including the presence of 1/f noise, asymmetric beams, realistic scanning strategies, temperature to polarization leakage, and bandpass mismatch. A more comprehensive study of the CORE mission and its capabilities will need to apply component separation methods to such simulations to assess the impact of these effects on the fidelity of the component separation. Note that Dick et al. (2010) have previously demonstrated the highly detrimental effect of calibration errors on various classes of component separation algorithms.

4 Component separation and likelihood analysis

Several component separation approaches have been proposed in the literature for cleaning foregrounds in CMB temperature and polarization maps, although the focus has predominantly been on intensity (Leach et al. 2008). For $B$-mode detection, several blind and parametric methods have been proposed.

Blind methods (e.g., Tegmark and Efstathiou 1996; Delabrouille et al. 2003; Baccigalupi et al. 2004a; Delabrouille et al. 2009; Betoule et al. 2009; Basak and Delabrouille 2012; Fernández-Cobos et al. 2016) exploit minimal prior information on the foregrounds, therefore they are not prone to systematic errors due to mis-modelling of the foreground properties. One of the drawbacks of such approaches is that foreground error estimation is difficult and must usually rely on Monte-Carlo simulations or other bootstrapping techniques.

A complementary strategy is provided by parametric methods (e.g, Eriksen et al. 2008; Stompor et al. 2009a; Ricciardi et al. 2010; Remazeilles et al. 2016), that explicitly model the foreground properties by means of a set of parameters that are fitted to the data. Such methods provide an easy way to characterise and propagate foreground errors, but their effectiveness depends on the consistency of the foreground model adopted.

Which of these two strategies will yield the best results with future experiments ultimately depends on the intrinsic complexity of the true sky in polarization and our ability to model it. This, and the need for cross-checks, continues to motivate the community to develop multiple independent approaches.

In this work, we use the Commander (Eriksen et al. 2008), and the NILC (Delabrouille et al. 2009; Basak and Delabrouille 2013), and SMICA methods (Delabrouille et al. 2003; Cardoso et al. 2008), as representative of the parametric and of the blind approaches, respectively. As detailed in Appendix A.1 the Commander code is currently limited by computational resources to the analysis of low-resolution maps and by the precision with which foregrounds can be modelled. However, it is ideally suited for studies on large angular scales and specifically the reionization peak ($2\leq\ell\leq 47$). Conversely, while NILC (Appendix A.2) and SMICA (Appendix A.3) can process full resolution maps, their ability to clean the data from foregrounds at low $\ell$ are limited by the minimum variance that can be reached depending on the number of modes and available frequency channels. In terms of $B$-mode detection, these methods more naturally target the recombination bump ($\ell\sim 100$). Since the two types of approaches are complementary, they can be considered to form low-multipole and high-multipoles analysis pipelines respectively.

The fidelity of the $B$-mode component separation is assessed by evaluating the tensor-to-scalar ratio $r$ through an appropriate cosmological likelihood function. Since the $E$-mode spectrum is an important input to delensing methods as applied to the $B$-mode spectrum, we also quantify its accuracy by fitting the optical depth to reionization, $\tau$ using a likelihood formalism. However, we fix all other cosmological parameters to their input values, thus the derived errors will not be representative of those determined from a combined analysis of the spectra derived from both temperature and polarization.

4.1 Point-source detection and masking

Point-sources are expected to constitute the dominant astrophysical foreground for CMB $B$-modes on scales beyond the recombination peak (Curto et al. 2013). The detection and subsequent masking or removal of such sources is therefore an important pre-processing step to be applied to the data before the application of component separation algorithms that predominantly target diffuse Galactic foreground emission.

Two different procedures have been used to create the temperature and polarization point-source masks. First, a blind source detection pipeline is run on the 19 CORE intensity maps, producing catalogues of sources above a given signal-to-noise ratio in each of the simulated maps. Then, a different pipeline is run on the polarization $Q$ and $U$ maps using as input the positions of the sources detected in intensity. In both cases, intensity and polarization, the single frequency masks are built using a hole radius of about 1.25 times the FWHM of each channel. All the intensity masks are then combined into a single intensity union mask. Similarly, all the polarization masks are combined together into a single polarization union mask. At the end of the process we have 40 masks, 19 intensity masks, 19 polarization masks and the two combined union masks.

4.1.1 Intensity point-source analysis

The detection of compact sources in the 19 simulated intensity CORE frequency maps is attempted on a frequency-by-frequency basis by searching for peaks above a given signal-to-noise ratio (SNR) in maps that have been optimally filtered using the Mexican Hat Wavelet 2 (hereafter MHW2 López-Caniego et al. 2006; González-Nuevo et al. 2006). Such a wavelet is employed as a filter because it simultaneously removes both small scale noise fluctuations and large scale structures in the vicinity of the compact sources, thereby improving the detection process. The analysis is performed on flat patches of the sky, generated from full sky maps in HEALPix format. Specifically, a given map is sub-divided into a sufficiently large number of overlapping projected patches so as to effectively cover the full sky. Each of these patches is then filtered with the MHW2 taking into account the local statistics of the background in the vicinity of the source in order to optimize the shape the filter. After filtering, candidate sources are identified and repetitions from overlapping regions removed, retaining those sources with the highest SNR. The process is repeated twice, first, following this blind detection procedure on the full sky, then, second, repeating the analysis with the patches centred on each of the previously detected sources. In this second iteration, the sources are characterized again and those with a SNR > 5 are kept. As a result, a catalogue of sources is produced with information about the position, amplitude and SNR of each object. This procedure is identical to that used to produce the Planck compact source catalogues (Planck Collaboration XXVIII 2014; Planck Collaboration XXVI 2016b).

In order to produce the corresponding compact source masks, a simple masking procedure has been followed, masking all the pixels within a radius proportional to the beam FWHM of each simulated map. This number depends on the amplitude of the source allowing for slightly larger masked regions when the sources are very bright and slightly smaller masked regions for weaker sources. Typically, the radius of the masked holes is of the order of three times the beam ($3\sigma$), where the beam is defined as $\sigma\approx FWHM/\sqrt{8\ln 2}$.

4.1.2 Polarization point-source analysis

The detection of compact sources in the CORE polarization $Q$ and $U$ maps is performed in a non-blind fashion by assessing the significance of the polarized signal at the position of the sources detected in the intensity maps. The procedure is a two-step process as proposed in López-Caniego et al. (2009). First, a maximum likelihood filter is applied on the $Q$ and $U$ maps (Argüeso et al. 2009). Second, a map of polarization intensity, $P=\sqrt{Q^{2}+U^{2}}$, is built from the filtered maps of $Q$ and $U$, and the polarized flux density at the position of the source is calculated. Then, analyzing the background statistics of the $P$ map, a $99.90$,% significance threshold is calculated. If the flux density at the position of the source is above this significance threshold, we consider that the signal is from the source and not from the background. Figure 4 shows the union mask for the polarized compact sources that we have detected in the nineteen CORE frequency bands.

4.1.3 Point-source pre-processing

For some foreground cleaning methods, e.g. NILC (see Appendix A.2) and SMICA (see Appendix A.3), the removal of polarized point-sources from the sky maps prior to component separation is required to optimize the reconstruction of the CMB $B$-mode power spectrum. In the case of NILC, the presence of any bright point-source increases the value of the local frequency-frequency covariance matrix in the neighbourhood of the source since the set of pixels from which the covariance is computed is more extended than the source itself. The NILC weights then ultimately adjust themselves to suppress the point-source power as much as possible at the expense of a lower control of the diffuse foreground signal around the source. By pre-processing the point-sources in the CORE sky maps, the NILC weights computed from the pre-processed maps readjust themselves to better suppress the contributions from diffuse foregrounds and noise. In the case of SMICA, the performance of foreground cleaning relies on the correct assumption of the dimension of the foreground subspace (i.e., the number of independent foreground degrees of freedom or the rank of the foreground covariance matrix). Appropriate masking of the diffuse Galactic foreground emission must ensure that the dimension of the foregrounds is close to constant over the sky, as assumed by SMICA. However, the presence of point-sources locally increases the effective dimension of the foregrounds, which are then no longer uniform over the sky. By pre-processing the point-sources in the CORE sky maps, we avoid any local increase of the dimension of the foreground subspace, therefore optimizing the foreground removal by SMICA.

Many techniques have been proposed in the literature to remove point-sources from sky maps, either relying on a direct fit of the source profiles which are then subtracted, or on the restoration of signal to the masked source pixels by extrapolating the background signal determined from neighbouring pixels (often referred to as ’inpainting’). Different approaches have been developed for restoring missing data, including sparse representation of the data in a wavelet frame (e.g., Abrial et al. 2008), constrained Gaussian realisations (e.g., Bucher and Louis 2012), and minimum curvature spline surface inpainting (e.g., Remazeilles et al. 2015).

In this paper, a shortcut has been adopted to subtract the brightest polarized sources from the CORE sky maps. Specifically, the point-source component maps input to the PSM simulation at each CORE frequency are masked using the source masks produced in Sect. 4.1.2. These masked source maps are then coadded to the other diffuse foreground component maps in order to produce the simulated CORE sky maps. The resultant simulated skies contain no contribution from those sources that have been detected in polarization as would result from a perfect source subtraction approach, although a background of undetected polarized sources remains.

4.2 Bayesian parametric fitting at low multipoles using Commander

We have applied the Commander algorithm (Eriksen et al. 2008) to the simulated multi-frequency CORE $Q$ and $U$ sky maps described in Sect. 3 in order to separate the individual components of polarized emission. Commander performs a parametric fit to the data, on a pixel by pixel basis, in a Bayesian framework by using Markov Chain Monte Carlo (MCMC) Gibbs sampling (Wandelt et al. 2004). For details of the methodology, we refer the reader to Sect. A.1. Since the fit is performed for each pixel and requires a large number of Gibbs samples in the MCMC, the algorithm has a large computational cost. This implies that the method is better suited for the analysis of large angular scales and low-resolution maps, where the primordial CMB $B$-mode signal is most significant. The Commander component separation method has previously been applied successfully to Planck data to recover the CMB and foreground contributions at low multipoles (Planck Collaboration X 2016).

We follow the methodology described in Remazeilles et al. (2016) for the analysis of $B$-modes. The simulated data, $\boldsymbol{d}(\nu,p)$, consist of both Stokes $Q_{\nu}(p)$ and $U_{\nu}(p)$ polarization maps at CORE frequencies, $\nu$, that were simulated directly at the HEALPix $N_{\rm side}=16$ resolution. This avoids issues connected to the downgrading of sky maps and its impact on foreground spectral behaviour which has consequences for parametric component separation algorithms. Specifically, using lower resolution simulations ensures that any residuals in the data after component separation are not due to a mismatch between the assumed spectral model and that induced by the downgrading. We refer the reader to Sect. 6.2 for a more detailed discussion. Therefore, we use Commander to reconstruct the CMB $B$-mode power spectrum on scales corresponding to reionization at low multipoles $2\leq\ell\leq 47$.

In order to partially mitigate the impact of diffuse foreground contamination on the Bayesian fit by Commander, we have applied a Galactic mask to the CORE frequency maps so as to avoid the brightest regions of the sky. The mask has been constructed as follows. The polarization intensity, $P=\sqrt{Q^{2}+U^{2}}$, of the CORE 600 GHz map, which is dominated by thermal dust emission, was computed and extrapolated to $70$  GHz (the foreground minimum) through a modified blackbody spectrum with $\beta_{d}=1.6$ and $T_{d}=19.4$ K. A dust mask was then defined by thresholding the extrapolated $P$ map with respect to a simulated CMB $P$ map at $70$ GHz. The regions of the sky where the extrapolated $P$ map amplitude exceeds ten standard deviations of the simulated CMB $P$ map signal at $70$ GHz are masked. Similarly, we have used the CORE 60 GHz map, which is the map with most synchrotron emission, to compute a corresponding $P$ map that is then extrapolated to $70$  GHz using a power-law spectrum with $\beta_{s}=-3$. A similar thresholding procedure defines a synchrotron mask. The combination of these two masks then mitigates the brightest synchrotron and dust contamination. Figure 5 shows the Galactic mask (grey), superimposed on the polarization intensity map of the combined CORE $60$ and $600$ GHz data extrapolated to $70$ GHz (background colours). The mask leaves $f_{\rm sky}=51$% of the sky usable. It should also be noted that, despite masking the sky maps, the CMB map reconstructed by Commander is still a full-sky map because of the inference of the $C_{\ell}$ to the map and vice-versa through the Gibbs sampling iteration scheme described by Eq. A.9. Of course, the application of a mask is not mandatory for component separation, and foreground parameters can be computed over the full sky if required.

The Commander component separation outputs consist of a set of MCMC $C_{\ell}$ samples characterized by the posterior distribution Eq. (A.20). As a consequence, the Commander CMB $C_{\ell}$ samples are particularly suited to a Blackwell-Rao type estimation of the cosmological parameters $\tau$ and $r$. Note that all other cosmological parameters are set to the Planck 2015 best-fit values as used in the simulations.

Given that the theoretical CMB $E$-mode power spectrum scales as (e.g., Planck Collaboration XI 2016)

[TABLE]

then the posterior distribution of the optical depth to reionization, $\tau$, can be computed with the Blackwell-Rao estimator using the MCMC Gibbs samples, ${\widehat{C}}_{\ell}^{EE\,(i)}$, of the reconstructed CMB $E$-mode power spectra:

[TABLE]

where the sum runs over $N_{G}$ Gibbs samples and the log-likelihood function reads as

[TABLE]

If we then consider the theoretical CMB $B$-mode power spectrum as the combination of the tensor modes and the lensing modes,

[TABLE]

then the posterior distribution of the tensor-to-scalar ratio, $r$, and the amplitude of the lensing $B$-modes, $A_{lens}$, can be computed with the Blackwell-Rao estimator using the MCMC Gibbs samples, ${\widehat{C}}_{\ell}^{BB\,(i)}$, of the reconstructed CMB $B$-mode power spectra.

[TABLE]

where the log-likelihood function reads as in Eq. 4.3. We specify the quantity ${\mathcal{D}_{\ell}\equiv\ell(\ell+1)C_{\ell}/2\pi}$ for plotting purposes. We also define ${\Delta\mathcal{D}_{\ell}/\sigma_{\ell}=\left(\widehat{\mathcal{D}}_{\ell}-\mathcal{D}^{th}_{\ell}\right)/\sigma_{\ell}}$ as the difference between the estimated power spectrum, $\widehat{\mathcal{D}}_{\ell}$, and that of the input theoretical model sky, $\mathcal{D}^{th}_{\ell}$, normalized by the error bar, $\sigma_{\ell}$, on $\widehat{\mathcal{D}}_{\ell}$.

Throughout the paper, the reconstructed angular power spectra are binned over specific multipole ranges in order to minimize the correlations between angular scales. We adopt the following weighting in the binning procedure,

[TABLE]

Five multipole bins are defined for Commander: $\ell\in[2,6]$, $[7,12]$, $[13,20]$, $[21,31]$, and $[32,47]$.

In Fig. 6, we first show the reconstruction of the CMB $E$-mode power spectrum after component separation. The reconstruction of the CMB $E$-mode signal is easier than in the case of the primordial CMB $B$-mode signal since the latter is sub-dominant to foreground signals over all angular scales. The accurate reconstruction of the CMB $E$-modes fulfils several functions: it serves as a validation criterion for the fidelity of the component separation algorithm; E-modes can be used as an input to delensing algorithms (Challinor et al. 2017 , in prep.) operating via quadratic estimators (Hu and Okamoto 2002; Hirata and Seljak 2003); and finally it provides important information on cosmological parameters, in particular the optical depth to reionization (Di Valentino et al. 2016).

By applying the Blackwell-Rao estimator Eqs. (4.2)-(4.3) to the Commander reconstructed $E$-mode power spectrum over the multipole range $2\leq\ell\leq 12$, we determine the posterior distribution for the optical depth to reionization (shown in the right panel of Fig. 6), indicating an unbiased estimate of $\tau=0.0522\pm 0.0044$. This is a more than $11\sigma$ measurement of $\tau$ by CORE to be compared to the $6\sigma$ measurement from the latest Planck results (Planck Collaboration XLVI 2016).

The left panels of Fig. 7 present the reconstructed CMB $B$-mode power spectra as determined by Commander for different values of the tensor-to-scalar ratio, in the absence of gravitational lensing and when only synchrotron and dust foregrounds are present. These correspond to simulations #1 and #2, as described in Tables 1 and 2, plus an additional simulation #2-bis with $r=2.5\times 10^{-3}$. In these three cases, Commander successfully recovers the primordial CMB $B$-mode signal from CORE simulations over the multipole range $2\leq\ell\leq 47$ and provides unbiased measurements of the tensor-to-scalar ratio, finding $r=\left(0.97\pm 0.09\right)\times 10^{-2}$, $r=\left(4.6\pm 0.4\right)\times 10^{-3}$, and $r=\left(2.3\pm 0.3\right)\times 10^{-3}$ respectively (as seen in the right panels of Fig. 7). These correspond to detections of $r$ at a level exceeding $8\sigma$.

When the tensor-to-scalar ratio falls to the $r=10^{-3}$ level, the reconstruction of the primordial CMB $B$-mode signal by Commander at low-$\ell$ becomes more challenging. Figure 8 presents results for such a case in the presence of dust and synchrotron foregrounds with variable spectral indices, but in the absence of lensing (simulation #3). In particular, the reconstructed CMB $B$-mode power spectrum (top left panel) shows evidence of residual foreground contamination on scales $\ell\leq 12$. As a consequence, the Blackwell-Rao estimate of $r$ (top right panel) is biased by $3\sigma$. In the absence of lensing the uncertainty on $r$ after foreground cleaning is $\sigma(r=10^{-3})=10^{-4}$ (see Table 4 in Sect. 5), therefore the sensitivity target of CORE, $\sigma(r)=3\times 10^{-4}$, is achieved by Commander, even when including the bias in the uncertainty.

In order to provide some intuition on the origin of the bias in the estimated primordial $B$-modes for $r=10^{-3}$, we plot in Fig. 9 the recovered distributions of the foreground spectral indices by Commander (red) with respect to the input distributions (black). While the dust spectral indices and dust temperatures are recovered with the desired accuracy, the recovered distribution of the synchrotron spectral index presents slight deviations, $\Delta\beta_{s}\sim 0.02$, with respect to the input distribution. Those small discrepancies might still be significant for $B$-modes at level of sensitivity of $r=10^{-3}$. Scaling the synchrotron map at $145$ GHz through either the recovered or the input spectral index distribution actually yields to a difference of $B$-mode power of order $\sim 2\times 10^{-5}$ $\mu{\rm K}^{2}$ on reionization scales, as shown in the bottom right panel of Fig. 9. Therefore, errors on the synchrotron spectral index of order $\Delta\beta_{s}\sim 0.02$ are sufficient to add an excess of $B$-mode power on top of the primordial signal at $r=10^{-3}$, but still below $r=5\times 10^{-3}$. Of course, the fit is multidimensional so that the exact shape of the excess power might be different, but this calculation at least provides the order of magnitude of those effects at reionization scales. Interestingly, the means and standard deviations of the recovered and the input distributions of the synchrotron spectral index are consistent, but the skewness and kurtosis are not because of the non-Gaussian nature of the distribution, which yields to imperfect recovery of the synchrotron spectral index over the sky. For accurate Bayesian fitting at the level of sensitivity of $r=10^{-3}$, we may think of using more informative priors on the foreground spectral indices instead of the Gaussian priors adopted in this analysis. This will be investigated in a future work. We also anticipate that ongoing low-frequency surveys, e.g. C-BASS (Irfan et al. 2015), will provide tighter constraints on the synchrotron spectral index at a precision smaller than $\Delta\beta_{s}\sim 0.02$, which will help the parametric component separation. Paradoxically, using chi-square statistics to evaluate the agreement between the modelled total sky emission and the simulated data indicates that the fit is statistically adequate. However, the foreground residuals present in the reconstructed spectrum and the bias on $r$ indicate that the recovered CMB component is not perfectly fitted. As discussed in Remazeilles et al. (2016), the absence of chi-square evidence for incorrect foreground modelling arises from a lack of frequency coverage. With more channels at frequencies below $60$ GHz, any mismatch $\Delta\beta_{s}\sim 0.02$ between the input and modelled foregrounds would be indicated by the chi-square statistics, thus avoiding either biased or false detections of the tensor-to-scalar ratio. However, this assumes that the spectral index above 60 GHz is the same, an assumption that is not guaranteed to hold at the required level of precision.

Residual $B$-mode foreground contamination at the reionization peak could, in principle, be corrected for in the likelihood estimation of $r$. After introducing a parametrized phenomelogical model with a specific spectral shape for the residual foregrounds, then $r$ and the foreground nuisance parameter, $A_{\rm FG}$ (the amplitude of the power spectrum of residual foregrounds), could be fitted simultaneously to the Commander reconstructed CMB $B$-mode power spectrum. We defer an attempt to implement such a strategy to Sect. 5.

The presence of gravitational lensing effects, which transform $E$-modes into $B$-modes, is an additional complication for the measurement of the primordial signal. In simulation #4, we attempt to analyse a simulation more closely matching reality by including both the lensing $B$-mode contribution and additional foreground contributions from AME, and point-sources. The results are shown in the bottom panels of Fig. 8. Note that the reconstruction of the lensed CMB $B$-mode power spectrum is unbiased even at the reionization peak. We then apply the Blackwell-Rao estimator Eqs. (4.3)-(4.5) to derive the posterior distributions for $r$ by fixing $A_{lens}=1$. Although treating the lensing contribution effectively as a nuisance term does allow the bias on $r$ to be accommodated, it does not correct for the additional cosmic variance introduced by the lensing signal. The effect of this is to reduce the significance of the measurement of $r$, as indicated in the right hand panel of the plot. Specifically, we obtain $r=\left(1.3\pm 0.5\right)\times 10^{-3}$ when fixing the amplitude of the lensing $B$-modes to $A_{lens}=1$. Although the information on $A_{lens}$ derived from $T$- and $E$-modes, together with the possible lensing potential reconstruction from $E$-modes, would help this analysis, the inclusion of additional modes at higher $\ell$ obtained from appropriate component separation approaches such as NILC or SMICA would clearly be beneficial. This is discussed in Sect. 5. In addition, actual delensing of the CMB map with the quadratic estimators proposed in Hu and Okamoto (2002); Hirata and Seljak (2003) will reduce the cosmic variance on $r$ (Challinor et al. 2017 , in prep.). With this in mind, Sect. 5 presents revised estimates of $r$ after imposing an effective $60$% delensing of the data, in addition to foreground cleaning (simulation #5).

Of course, the Commander results also include maps of the reconstructed physical components. Figure 10 compares the the input CMB $Q$ map (top left panel) with the reconstructed CMB $Q$ map (top right panel), and presents their difference in the bottom left panel. The bottom right panel shows a map of the Commander $\chi^{2}$ statistic given by Eq. (A.8) divided by the number of input frequency channels, which measures the mismatch between the fitted model and the data in each pixel. In principle, the $\chi^{2}$ map can be useful for constructing a posteriori an optimal mask to then be employed for parametric fitting, by rejecting any pixels where the $\chi^{2}$ value is too large. We have not adopted this approach here. We also note that the inspection of foreground residuals in the reconstructed CMB $Q$ and $U$ maps is not strongly informative with respect to $B$-modes, since they are dominated by the $E$-mode contribution. However, it can highlight a significant failure in component seperation, likely driven by incorrect assumptions in the foreground models adopted by Commander, also indicated by the $\chi^{2}$ map.

In Fig 11, we show the reconstruction by Commander of the foreground spectral indices for thermal dust and synchrotron for simulation #4. The inputs and Commander estimates can be compared between the left and right columns. While the dust spectral indices (top panels) are accurately recovered over the (unmasked) sky, with a similar fidelity for the dust temperature (middle panels), it is evident that the reconstructed synchrotron spectral indices (bottom right panel) are noisy. CORE does not have enough sensitivity to measure the synchrotron parameters in the $60-600$ GHz frequency range, and would benefit from additional low-frequency information, either from external measurements provided by other instruments such as C-BASS (Irfan et al. 2015) or QUIJOTE (Rubiño-Martín et al. 2012), or by the addition of a low-frequency channel below $60$ GHz (see the discussion in Sect. 6.4).

4.3 Needlet Internal Linear Combination at high multipoles

The NILC component separation method (Delabrouille et al. 2009; Basak and Delabrouille 2012, 2013; Planck Collaboration IX 2016) has been applied to the CORE simulations to determine the power spectrum of the polarized CMB signal at high-$\ell$. For details of the method, we refer the reader to Sect. A.2.

The simulated sky maps are first convolved or de-convolved in harmonic space to put them all at the same angular resolution prior to the application of the NILC algorithm. Here we pick that of the smallest beam, i.e. $4^{\prime}$ (see Table 3). Since the method as currently implemented is applicable to scalar fields on the sphere, sky maps of the $E$ and $B$ modes are constructed from the input Stokes parameters, $Q$ and $U$, on the full sky (i.e., without masking). The NILC weights used to combine the multi-frequency input data in order to estimate the CMB signal are then computed from the full-mission $E$ an $B$ sky maps. The derived weights are also applied to the half-mission maps, which are later used for both power spectrum and noise estimation.

Figure 12 shows for the simulation #4 the full sky average of the NILC weights for each frequency channel and needlet band (as specified in Sect. A.2). The overall features of the fullsky average of the NILC weights at different frequencies and needlet bands do not notably change when AME, point-sources, and lensing effects are included or not in the simulation. The needlet weights are mostly determined by the galactic contamination, which dominates on large angular scales, and by the noise level, which dominates on small angular scales. It is apparent that most of the contributions to the reconstructed CMB polarization maps come from those channels with frequencies around $130$ GHz and $255$ GHz. On large angular scales (the first needlet bands), the $130$ GHz channel gets a more significant weighting than the $255$ GHz channel. However, on small angular scales (the last needlet bands), the situation is reversed due to the higher angular resolution of the $255$ GHz channel. Although the weights of the remainder of the frequency channels are relatively low, they are important for removing the Galactic foreground contamination.

However, as shown by Eq. A.24, the reconstructed CMB $E$- and $B$-mode maps cannot be completely free from contamination by residual foregrounds and noise. Figure 13 clearly shows residual foreground contamination along a narrow strip of the Galactic plane for E modes, while the B-mode reconstructed map is dominated by residuals of galactic contamination. Therefore, for further analysis, a set of conservative masks are derived from the full-mission NILC CMB maps. These have been generated using the following procedures.

For $E$-modes, the NILC CMB map is filtered in harmonic space through the Gaussian window

[TABLE]

peaking at $l_{centre}=2275$ which corresponds to the multipole where the sum of the power of the CMB and noise is minimum. The resulting map is then squared and smoothed with a Gaussian beam of FWHM=$8^{\circ}$. The variance map obtained in this way is then corrected for the noise contribution by subtracting a noise variance map obtained using the same procedure as applied to the half-difference of the NILC half-mission maps. The $E$-mode confidence mask shown in Fig. 14 (left panel) is then obtained by thresholding the noise-corrected variance map. The NILC CMB $B$-mode map is dominated by noise at almost all angular scales. Since most of the power of the CMB $B$-modes is concentrated on large angular scales, the NILC CMB map is filtered through a Gaussian window function of FWHM=$30^{\prime}$. The resulting map is then squared and smoothed with a Gaussian beam of FWHM=$20^{\circ}$. The resulting variance map is corrected for the noise contribution as before. The confidence mask for $B$-modes shown in the right panel of Fig. 14 is then obtained by thresholding the noise-corrected variance map.

For the computation of the angular power spectrum of the high resolution NILC CMB maps, we use a pseudo-$C_{\ell}$ estimator (Hivon et al. 2002; Chon et al. 2004; Szapudi et al. 2005). This method is computationally much faster than maximum likelihood and provides optimal results at intermediate to high $\ell$s. In order to compute the covariance on our measurement of angular power spectrum, we have followed the method described in Tristram et al. (2005). More details are given in Appendix A.2. The impact of instrumental noise residuals on the measurement of angular power spectra is avoided by evaluating cross-spectra between the NILC half-mission maps. Then, in order to reduce effects related to the sharp edges of the masks on the measurement of the power spectra, the masks are apodized through a cosine transition with length $1^{\circ}$ for $E$-modes and $2^{\circ}$ for $B$-modes, the apodization length for $B$-mode masks being larger because of the weakness of the signal compared to $E$-modes. Given that the effective sky fraction of an apodized mask is $f_{\rm sky}=\sum_{i}w^{2}_{i}\Omega_{i}/(4\pi)$, where $w_{i}$ is the value of the mask in pixel $i$ and $\Omega_{i}$ is the solid angle of the pixel, the resulting $E$- and $B$-more masks have retained sky fractions of 80% and 60% respectively. The computed pseudo-$C_{\ell}$ spectra are then used as input to a likelihood-based estimation of the cosmological parameters.

Unlike in the case for Commander, where a non-Gaussian likelihood has been used via the Blackwell-Rao approximation to calculate the posterior distribution of either the optical depth to reionization, $\tau$, or the tensor-to-scalar ratio, $r$, here we adopt a simple Gaussian likelihood. Given an estimated and binned CMB power spectrum $\widehat{C}^{XX}_{b}$, and a model CMB power spectrum $C^{th\,XX}_{\ell}(p)$, which, in principle, can depend on any arbitrary number of parameters $p$, we define the likelihood as

[TABLE]

where $b,b^{\prime}$ runs through the corresponding multipole bins, $C^{th\,XX}_{b}(p)$ corresponds to the binned model power spectrum, and $\left(\widehat{\Xi}^{XX}\right)^{-1}_{bb^{\prime}}$ is the inverse of the covariance matrix. The covariance matrix is computed for a fiducial model, defined by the input tensor-to-scalar ratio.

When considering the likelihood for $\tau$, then $XX$=$EE$ and the theoretical CMB $E$-mode power spectrum is given by Eq. 4.1. In the case of $r$ and $A_{\rm lens}$, then $XX$=$BB$ and the model $B$-mode power spectrum is defined as in Eq. 4.4. However, in either case, it may be necessary to consider the addition of various nuisance terms to account for diffuse foreground residuals at low-$\ell$ or point-source residuals at high-$\ell$. These effects, $i$, will then be specified by a template spectrum, $C_{\ell}^{BB,i}$, accompanied by the corresponding nuisance amplitude, $A_{i}$, used to mitigate the effect (Hervías-Caimapo et al. 2017).

The $A_{\rm lens}$ parameter may not be well constrained by the $B$-mode spectra alone. However, tight priors on it can be derived from other CMB power spectra, together with the lensing power spectrum. Indeed, the uncertainty on $A_{\rm lens}$ for CORE has been shown to be 0.013 from a joint analysis of the temperature and $E$-mode spectra, improving to 0.012 when combined with the lensing power spectrum, (see Di Valentino et al. 2016). This justifies the choice below to fix $A_{\rm lens}$ at the fiducial value for a $B$-mode only likelihood.

The left panel of Fig. 15 shows the $E$-mode angular power spectrum derived from the NILC CMB reconstruction. The reconstruction of the CMB $E$-modes is of good quality over a large range of angular scales up to $\ell=3000$, with only three outliers biased by more than $3\sigma$. The resulting estimate of the optical depth to reionization is $\tau=0.054\pm 0.0022$, therefore a more than $24\sigma$ measurement of $\tau$ by CORE when accounting for multipoles $20\leq\ell\leq 359$. When we add more multipoles to the likelihood, we measure unbiased $\tau$ estimates with higher significance, although in real-sky data analysis other $\Lambda CDM$ parameters should be jointly fitted with $\tau$ on small angular scales. Our $\tau$ estimates are just used as a validation criterion of the accurate CMB $E$-mode reconstruction, not as forecasts on the measurement of $\tau$ by CORE. Detailed forecasts on the measurement of $\tau$ by CORE are given in the accompanying paper on cosmological parameters (Di Valentino et al. 2016), where it has been found that $\tau$ can be detected at $\sim 29\sigma$ significance when using multipoles up to $\ell=3000$ and marginalizing over the other $\Lambda CDM$ cosmological parameters, but neglecting the foreground contamination. Using more multipoles and negelecting foregrounds obviously improves the significance of the detection, while marginalizing over other cosmological parameters decreases it.

As in the case of Commander, the accurate reconstruction of the CMB $E$-mode power spectrum serves as a validation of the component separation method, while guaranteeing a good proxy for potential delensing of CMB $B$-modes.

The left panels of Fig. 16 present the reconstructed primordial CMB $B$-mode angular power spectrum by NILC for two values of the tensor-to-scalar ratio, $r=10^{-2}$ (top) and $r=5\times 10^{-3}$ (bottom) in the absence of gravitational lensing and when only synchrotron and dust foregrounds are present. The binning of the reconstructed $B$-mode angular power spectra again follows the weighting defined in Eq. 4.6. The five multipole bins defined for NILC are $\ell\in[2,14]$, $[15,47]$, $[48,124]$, $[125,224]$, $[225,349]$. While the reconstruction of the reionization peak of the primordial CMB $B$-mode power spectrum is evidently biased by more than $3\sigma$ in both cases, the recombination peak is accurately recovered by NILC. The corresponding estimates of the tensor-to-scalar ratio seen in Fig. 16 (right panels) are therefore based on the three NILC multipole bins spanning the range $48\leq\ell\leq 349$. In the two cases examined, we find $r=\left(1.19\pm 0.25\right)\times 10^{-2}$ and $r=\left(7.2\pm 2.4\right)\times 10^{-3}$ respectively.

When the tensor-to-scalar ratio falls to $r=10^{-3}$, the reconstruction of primordial CMB $B$-modes becomes more problematic. Figure 17 shows the $B$-mode power spectra derived from the NILC CMB sky maps both when lensing effects are either excluded (top) or included (bottom). There is a significant bias on many angular scales, and clear consequences for the estimated tensor-to-scalar ratio. We have estimated the residual foreground contamination in the NILC maps by applying the NILC weights that were used for the reconstruction of the CMB $B$-mode map to the individual foreground components. As shown in Fig. 17, the residual Galactic foregrounds, thermal dust (solid magenta line) and synchrotron (solid green line), dominate the cosmological signal at $r=10^{-3}$ on large angular scales ($2\leq\ell\leq 50$), while residual extra-galactic point-sources (solid yellow line) are the main contaminant on the recombination peak and angular scales $\ell>50$.

The origin of the bias is the fact that the ILC (and its needlet variant, the NILC) is not a power spectrum estimation tool, but a tool to make minimum variance maps of a component of interest. The reconstructed map is contaminated by a mixture of foreground residuals and of noise. Residual foregrounds are seen in Fig. 17 to be always larger than the primordial $B$ modes. Unlike the noise, which does not contribute to the estimated power spectrum as it is uncorrelated between the half-mission maps, foreground emission residuals bias the estimated power spectrum, with no natural prescription to de-bias. However, the level of this bias is quite informative as to what is the level of foreground residuals in the reconstructed maps.

Two ways to monitor such residuals with NILC can be implemented. One could modify the way that weights are computed, to reduce foreground contamination rather than noise in the maps, e.g. following the idea of the MILCA algorithm (Hurier et al. 2013). One can also at least avoid false detection of primordial B-modes by computing the power spectrum using variable masks, to check for a dependence of the measured power spectrum on the level of foreground emission. If this test shows that there indeed is a dependence, it is a clear indication that the measurement is foreground-dominated.

4.4 Spectral Matching Independent Component Analysis at high multipoles

As a third independent component separation method, we apply the SMICA algorithm (Delabrouille et al. 2003; Cardoso et al. 2008) to the CORE simulations. This method, based on matching the data to a model in the spectral domain, allows the joint estimation of the spatial power spectra of the components and of the noise, and of their mixing coefficients. For details of the methodology, we refer the reader to Sect. A.3. In this paper, no attempt is made to reconstruct maps corresponding to the derived best-fit spectra.

Although the method works directly with $E$- and $B$-mode spectra, the analysis utilises a mask constructed in a similar manner to that of Commander (Sect. 4.2). Polarization intensity maps at 60 and 600 GHz are smoothed with a beam of 10 degrees FWHM and extrapolated to the foreground minimum frequency of 70 GHz assuming that the maps are dominated by synchrotron and thermal dust emission respectively. Each of the extrapolated maps is then thresholded with respect to the smoothed 70 GHz simulated $P$ map to form two independent masks, which are then combined. The resulting mask is then apodized with an apodization length of 4 degrees, generating an effective usable sky fraction of $f_{\mathrm{sky}}=0.40$, as shown in Fig. 18.

Figure 19 shows the reconstruction of the CMB $E$-mode power spectrum by SMICA for a simulation with optical depth to reionization $\tau=0.055$, and dust and synchrotron foregrounds. The recovered CMB $E$-mode power spectrum is consistent with the input theory spectrum at the $2\sigma$ level over the majority of angular scales up to $\ell=3000$, although there are quite a few outliers showing deviations exceeding $3\sigma$ on scales $\ell>200$. The reconstructed spectrum is an improvement on measurements by Planck (Planck Collaboration XI 2016) in terms of the absolute deviation from the theory spectrum. However, the error-normalized deviations are more significant since the uncertainties obtained with CORE are smaller, and indicate that further optimization in the SMICA approach should be explored. The resulting estimate of the optical depth to reionization is $\tau=0.0558\pm 0.0082$, a measurement of $\tau$ by CORE at a level exceeding $6\sigma$ when including multipoles over the range $20\leq\ell\leq 199$ (and holding all other cosmological parameters fixed at their input vaues). If we consider a multipole range with more bins, $20\leq\ell\leq 600$, we measure $\tau=0.077$. This result is biased by the bins above $\ell\sim 500$.

The left panels of Fig. 20 show the reconstructed primordial CMB $B$-mode power spectra determined by SMICA from simulations #1 and #2 for different values of the tensor-to-scalar ratio, in the absence of gravitational lensing and when only synchrotron and dust foregrounds are present. The binning of the reconstructed $B$-mode angular power spectra follows the weighting scheme defined in Eq. 4.6. The 18 multipole bins defined for SMICA are $\ell\in[5,12]$, $[13,19]$, $[20,27]$, $[28,34]$, $[35,47]$, $[48,57]$, $[58,70]$, $[71,86]$, $[87,105]$, $[106,130]$, $[131,159]$, $[160,196]$, $[197,241]$, $[242,296]$, $[297,364]$, $[365,447]$, $[448,549]$, $[550,600]$. In both cases, the recombination peak, between $\ell\sim 60$ and $\ell\sim 200$, is recovered by SMICA within the $1\sigma$ error bound. The right panels in these figures show the estimated tensor-to-scalar ratio obtained from the reconstructed CMB $B$-mode power spectra on angular scales $48\leq\ell\leq 600$. We recover values of $r=\left(1.2\pm 0.1\right)\times 10^{-2}$ and $\left(6.1\pm 0.9\right)\times 10^{-3}$ for the fiducial values of $r=10^{-2}$ and $r=5\times 10^{-3}$ respectively. Using only those multipole scales corresponding to the recombination peak and above, these results are consistent with detections of the tensor-to-scalar ratio at a level of $5\sigma$ or above.

However, as we have seen in previous sections, for a value of $r=10^{-3}$, the reconstructed CMB $B$-mode power spectrum can show evidence of foreground residuals and bias relative to the input theory spectrum. The top left panel of Fig. 21 presents results for simulation #3, and indicates that an excess of power is seen on most angular scales for $\ell\leq 400$ due to foreground residuals. The corresponding estimate of the tensor-to-scalar ratio (top right panel), $r=\left(2.1\pm 0.8\right)\times 10^{-3}$, is not formally biased, albeit more than $1\sigma$ high, but the error has increased so that only an upper limit on $r$ can be claimed.

In the presence of gravitational lensing effects (simulation #4), SMICA performs an accurate reconstruction of the lensed CMB $B$-mode power spectrum when $r=10^{-3}$ over a large multipole range $20<\ell<600$, and is consistent with the fiducial lensed CMB $B$-mode power spectrum within $1\sigma$ (see the bottom left panel of Fig. 21). The tensor-to-scalar ratio $r=10^{-3}$ is recovered without significant bias (bottom right panel of Fig. 21), but with a large uncertainty, when either fitting or fixing $A_{lens}$ in the likelihood, with an estimate of $r=\left(2.4\pm 1.4\right)\times 10^{-3}$.

Currently, it appears that the measurement of the primordial CMB $B$-mode power spectrum with SMICA based on CORE observations can not be achieved at the required level of accuracy and precision when $r=10^{-3}$, yet is possible for $r=5\times 10^{-3}$ in the absence of lensing effects. Evidently, further optimization of the methodology is desirable.

5 CORE results on CMB $B$-mode measurements

In this section, we attempt to determine improved estimates of the the tensor-to-scalar ratio by adopting a hybrid likelihood approach, following the work of Planck (Planck Collaboration XV 2014; Planck Collaboration XI 2016) and WMAP (Dunkley et al. 2009b). Specifically, we combine the Blackwell-Rao estimate of the CMB $B$-mode power spectrum from Commander on scales $2\leq\ell\leq 47$ (including the reionisation peak) with a pseudo-$C_{\ell}$ estimate from either SMICA or NILC on scales spanning the recombination peak and above ($\ell>47$). Specifically, we take the product of the posterior on $r$ from Commander, obtained via the Blackwell-Rao approximation and a non-Gaussian likelihood as presented in Sect. 4.2, with that from either NILC or SMICA, obtained using the Gaussian likelihood formalism described in Sect. 4.3, to obtain the resulting joint posterior. SMICA and NILC provide their own estimation of the covariance matrix (see e.g. Appendix A.2), to use in the Gaussian likelihood approach for $\ell>47$.

5.1 Component separation results on $\boldsymbol{r=10^{-2}}$ and $\boldsymbol{r=5\times 10^{-3}}$

Figures 22 and 23 present the results of the joint reconstruction of the primordial $B$-mode power spectra (left panels) and the joint likelihood fits (right panels) to simulations #1 and #2. For these values of the tensor-to-scalar ratio, both the reionization and the recombination peaks of the primordial $B$-modes are accurately recovered by CORE between $\ell=2$ and $\ell=200$, with neither of the the combined power spectra deviating from the input theory spectrum by more than $2\sigma$.

The right panels of Figs. 22 and 23 indicate that observations of a CMB sky with $r=10^{-2}$ by the current CORE configuration yield estimates of the tensor-to-scalar ratio of $r=\left(1.06\pm 0.07\right)\times 10^{-2}$ and $r=\left(1.00\pm 0.09\right)\times 10^{-2}$ by the hybrid Commander-SMICA or Commander-NILC analyses, respectively. For an input value of $r=5\times 10^{-3}$, the corresponding estimates are $r=\left(4.9\pm 0.4\right)\times 10^{-3}$ and $r=\left(4.7\pm 0.4\right)\times 10^{-3}$, respectively. These results correspond to high significance detections of the tensor-to-scalar ratio, at a level exceeding $10\sigma$ for unlensed realisations of the CMB sky after correcting for foregrounds consisting of synchrotron and thermal dust emission which have spectral variations over the sky.

5.2 Component separation results on $\boldsymbol{r=10^{-3}}$

We next consider the critical case where the tensor-to-scalar ratio is of order $r=10^{-3}$, initially excluding lensing effects (simulation #3). Figure 24 summarizes the results. As we have seen previously in Sect. 4.2, there is clear evidence of residual foregrounds over the reionization scales probed by Commander. As a consequence, the recovered estimate of the tensor-to-scalar ratio for both combinations of component separation methods is $r=\left(1.3\pm 0.1\right)\times 10^{-3}$. However, in order to mitigate the potential bias due to the foreground residuals, we implement a strategy outlined in Sect. 4.2, and fit a parametric nuisance model. The actual shape of the power spectrum of foreground residuals is not known a priori, and may differ from that of the input foregrounds. This is particularly plausible for those component separation methods which combine the frequency maps in harmonic space, such as NILC, which computes different weights over different multipole ranges and regions, depending on the relative ratio between foregrounds and noise contamination. A conservative first approximation, then, is to assume a power-law shape for the power spectrum of the foreground residuals at low multipoles, and consider that this, up to an arbitrary amplitude to be fitted, can be described by the Commander thermal dust map estimate, presented in Fig. 25. For the frequency coverage considered in the CORE simulations, thermal dust emission dominates, and therefore most of the foreground residuals will originate therein. The best-fitting spectral index of the power-law is determined by a direct fit to the $B$-mode power spectrum of this thermal dust map, calculated with a Quadratic Maximum Likelihood (QML) estimator. The Blackwell-Rao likelihood is then modified in order to fit simultaneously for both $r$ and the amplitude of the foreground residuals, $A_{\rm FG}$, over the Commander multipole range. This results in an effective debiasing of the tensor-to-scalar ratio, as shown in the right panels of Fig. 24.

Indeed, in this case, we find that the joint estimates of the tensor-to-scalar ratio estimated from Commander and either SMICA or NILC are $r=\left(1.02\pm 0.10\right)\times 10^{-3}$ and $r=\left(1.00\pm 0.10\right)\times 10^{-3}$ respectively, corresponding to $10\sigma$ detections of cosmological $B$-modes.

5.3 Component separation results on $\boldsymbol{r=10^{-3}}$ with lensing and delensing

Finally, we consider challenging simulations where both gravitational lensing effects and complex foregrounds comprising synchrotron, thermal dust emission, AME, strong and faint point sources have been included. Figure. 26 presents the results for simulation #4. In order to account for the lensing contribution which will otherwise bias the estimate of $r$, we fix $A_{\rm lens}=1$ and fit for $r$ only. The posterior distribution of the tensor-to-scalar ratio from the joint $B$-mode power spectra analyses is shown in the right panels of Fig. 26. The measured value of $r$ is $\left(1.4\pm 0.47\right)\times 10^{-3}$ and $\left(1.4\pm 0.49\right)\times 10^{-3}$ for the Commander-SMICA and Commander-NILC joint likelihoods, respectively. As is the case for simulation #3, foreground residuals on the reionization scales are a source of bias at low-$\ell$. If the nuisance parameter, $A_{\rm FG}$, is again fitted in the likelihood, the estimated tensor-to-scalar ratios are $\left(1.2\pm 0.6\right)\times 10^{-3}$ and $\left(1.2\pm 0.6\right)\times 10^{-3}$ for Commander-SMICA and Commander-NILC, respectively. Note that that the presence of the lensing signal leads to increased uncertainty in the measurement due to the related cosmic variance contribution, and reduces the measurements of $r$ effectively to upper limits.

Consequently, we consider the impact of delensing on the determination of the tensor-to-scalar ratio. In principle, CORE should allow the reduction of the lensing contribution to the measured $B$-mode power spectrum by a factor typically of order 60% (Challinor et al. 2017 , in prep.). In this paper, we adopt a shortcut to mimic such a delensing analysis. Specifically, we analyse simulation #5, which is identical to simulation #4 except that only 40% of the pure $B$-mode lensing signal is included in the CMB map. In order to properly account for residuals from delensing in the likelihood estimation of the tensor-to-scalar ratio, the $A_{lens}$ parameter is fixed to a value of $0.40$ while fitting for $r$ only. Figure 27 presents the results of the analyses, indicating best-fit values for the tensor-to-scalar ratio of $r=\left(1.70\pm 0.45\right)\times 10^{-3}$ and $r=\left(1.60\pm 0.49\right)\times 10^{-3}$ for the joint likelihood with Commander-SMICA and Commander-NILC, respectively. As before, fitting for foreground residuals at low-$\ell$ improves the consistency with the input value of $r$, resulting in estimates of $r=\left(1.20\pm 0.46\right)\times 10^{-3}$ and $r=\left(1.10\pm 0.51\right)\times 10^{-3}$. Note that there is no clear improvement in the error term, $\sigma(r=10^{-3})$, due to the delensing procedure. This is because the uncertainty is still dominated by foreground residuals. We discuss this important issue further in Sect. 6.1.

6 Discussion

6.1 Concerning B-mode delensing versus foreground cleaning

It is apparent from Table 4 that the overall uncertainty on $r=10^{-3}$ is not reduced much after de-lensing, either with Commander alone, or with Commander in combination with SMICA or NILC. The uncertainty on the tensor-to-scalar ratio remains at a level $\sigma(r=10^{-3})\sim 0.5\times 10^{-3}$ whether $60$% of the lensing is subtracted (simulation #5) or not (simulation #4). This means that the error on $r$ is not dominated by residual lensing. It is not dominated by noise in CMB channels either, since noise is lower than the lensing contamination. Hence, the dominant source of error can be either foreground uncertainty, propagating to the final power spectrum error, or cosmic variance, or a combination of both. Cosmic variance can dominate the error if most of the sensitivity to $B$ modes comes from the reionization bump.

We also note, however, that when de-lensing is performed, the error on $r$ from the SMICA fit alone decreases substantially, from $1.4\times 10^{-3}$ for simulation #4 to $1.0\times 10^{-3}$ for simulation #5. As our implementation of SMICA in this paper constrains $r$ in the range $48\leq\ell\leq 600$, this decrease shows that de-lensing has an impact on the detectability of the recombination bump. The overall error on $r$, however, is not impacted much because the lowest harmonic modes, measured by Commander in this paper, dominate the final sensitivity in this simulation.

To clarify further the impact of delensing on $\sigma(r)$, we apply Commander to additional specific simulations, as follows. In the case where there are no foregrounds at all in the observed maps, a run of Commander finds that $\sigma(r=10^{-3})$ decreases from $0.3\times 10^{-3}$ without delensing to $0.2\times 10^{-3}$ after $60$% delensing of the CMB $B$-mode signal (left panel of Fig. 28, corresponding to analyses where no foregrounds are included). There also is a notable difference in $\sigma(r=10^{-3})$ between the cases with and without foreground emission. It hence appears that when foreground emission is included, the global error budget is dominated by foreground-induced uncertainties.

The right panel of Fig. 28 presents results from the analysis of three simplified versions of simulation #5. Firstly, we consider a simulated sky including $40$% lensing residuals but with no point-source contamination. In this case (denoted by the solid red line), the uncertainty on $r$ is reduced to a level $\sigma(r=10^{-3})=0.4\times 10^{-3}$ (20% lower than the case with point sources). When Commander has to deal only with galactic foregrounds, it performs better than when the additional confusion due to polarized sources is included.

If in addition an idealized improvement in modelling the synchrotron foreground emission is imposed on this simulation by fixing the synchrotron spectral indices to their input values rather than fitting them, then the uncertainty on $r$ is reduced to the level $\sigma(r=10^{-3})=0.3\times 10^{-3}$. Finally, when the dust temperatures are fixed to their input values instead of fitting them, then in addition to reducing the uncertainty to the level of ${\sigma(r=10^{-3})=0.3\times 10^{-3}}$ we also suppress the shift on the best fit value of $r$ previously seen in analyses of the reionization scales.

Hence, when Commander is given information about the real frequency scaling in each pixel, the code can interpolate better to CMB frequencies the foreground templates observed at 60 GHz and 600 GHz. The final performance of the mission is as good as if there was no foreground contamination at all.

We note that combining the Commander low-$\ell$ solution with information from higher multipoles determined either by NILC or SMICA does not seem to help much in reducing the uncertainty on $r$, as compared with what is achieved by Commander alone. This is attributable to the fact that the two blind methods, as implemented in this work, provide information in a range of $\ell$ where the primordial B-modes are subdominant anyway. But one must keep in mind that the simulations analyzed here have been tuned to match the Commander model assumptions. The impact of this is discussed in the next subsection.

Possibly, further optimization of the methods can help to better detect the recombination peak, using as a prior improved understanding of the foregrounds (including the point-source contribution at high-$\ell$). Detecting the recombination bump in addition to the reionization bump is important to confirm the inflationary origin of any detected B-mode excess, as witnessed by the bias in Commander data points at low $\ell$. De-lensing is an important tool to achieve this.

6.2 Concerning pixelization-related effects on foreground parametrization

The emission laws of thermal dust and synchrotron are described, at least in part, by power-law spectral indices that vary over the sky. Even if we consider that, for any given line-of-sight, an effective spectral index can be defined (ignoring the fact that it may include emission with varying spectral properties along the line-of-sight itself), maps from any CMB experiment are pixelized at a finite spatial resolution and observed through a finite optical beam response. Thus a given pixel (and beam) encompasses many line-of-sights with different spectral indices. The coaddition of the different power-laws with different spectral indices within the finite pixel or beam solid angle can perturb the effective emission law away from the idealized power-law dependence, and result in an effective, unphysical, curvature of the spectral indices. This can represent a difficulty for parametric component separation methods faced with the challenge of correctly modelling foreground contamination to measure the primordial CMB $B$-modes. A related discussion on these issues for intensity can be found in Chluba et al. (2017).

In particular, we have experienced this issue when utilising the Bayesian parametric fitting method Commander to fit the spectral properties of foreground emission in low-resolution maps, at ${N_{\rm side}=16}$, which have been derived from simulated $N_{\rm side}=2048$ data sets, within which the foreground spectral dependencies vary from high resolution pixel-to-pixel.

We demonstrate the effect as follows. In Fig. 29, we explicitly fit both a spectral index, $\beta_{d}$, and a curvature term, $C$, i.e.,

[TABLE]

to simulated CORE dust frequency maps in which the simulated emission law is a perfect power-law with zero curvature ($C=0$) in the original simulation at $N_{\rm side}=2048$. We adopt a simple chi-square fitting procedure to illustrate the issue. The fit is either performed on the $N_{\rm side}=2048$ frequency maps degraded to $N_{\rm side}=16$ (top panels of Fig. 29) or on direct $N_{\rm side}=16$ simulations with one spectral index per low resolution pixel (bottom panels of Fig. 29). Clearly, degrading the pixelized sky maps generates a non-zero unphysical curvature in the dust spectral index (top right panel). Although the effective curvature of the dust spectral index caused by degrading the pixelized sky maps is low ($C\sim 0.05$ at maximum), omitting this systematic curvature in the fitted model can have a non-negligible impact on low $r=10^{-3}$ $B$-mode signal detection, and this can potentially bias the detection of CMB $B$-modes.

This issue is not specific to degrading sky maps to lower pixel size. It is more general and due to the finite-sized pixelization of the data and is a consequence of the foreground spectral indices varying from one line-of-sight to another with a characteristic scale smaller than the adopted pixel size. The scales over which the effect may be important are those related to the detailed physics of the interstellar medium, and the connection in relation to polarized dust and synchrotron emission with the magnetic field properties of the Galaxy.

6.3 Concerning point-source processing for $B$-mode studies

The impact of point source contamination on $B$-mode analyses is an issue that requires careful appraisal for both current and future CMB experiments that aim to detect the tensor-to-scalar ratio at levels of $10^{-3}$ and below. In general, the compact source contribution does not impact the large angular scales (near the reionization peak), but can play an important role on intermediate and small angular scales where the lensing-induced $B$-mode signal is present (Curto et al. 2013). Moreover, the consequences of point-source contamination of the $E$-mode signal to be used in delensing needs to be assessed. In fact, a careful investigation of the power spectrum of radio sources in polarization was carried out by Tucci and Toffolatti (2012) concluding that radio sources should not be a strong contaminant of the CMB $E$-mode polarization at frequencies above $\simeq 70\,$GHz. However, careful treatment is required for CMB $B$-mode polarization observations if the tensor-to-scalar ratio, $r$, is $\ll 10^{-2}$. If the lensing $B$-mode contribution is removed, the contamination by point sources remains the main nuisance. The problem can be addressed in several different ways.

Firstly, it is important to improve our ability to detect sources at the map level (in both intensity and polarization) using blind and non-blind techniques. The current approach selects those sources detected in intensity to mask out any possible related contamination in the $Q$ and $U$ maps. This is a very conservative approach that, given the typical compact source detection levels in CMB experiments, could result in the masking out of a significant fraction of the sky, particularly at lower frequencies where the beam sizes are large. Moreover, heavily masked maps can have implications for the fidelity of diffuse component separation methods as applied to the data.

Alternatively, an attempt can be made to assess whether polarized sources are present in the data given the list of detections in intensity. Such a search is made by studying a map of polarized intensity, $P$, and searching at all of the candidate positions. However, the statistical properties of the $P$ map can not be approximated by a Gaussian distribution, following instead a Rayleigh distribution. Thus, the usual thresholding approach based on variance-levels of the signal are not applicable, and the assessment must define significance levels that account for the local statistics in the vicinity of the putative source location. This is the method followed in this paper, and previously by Planck (Planck Collaboration XXVIII 2014; Planck Collaboration XXVI 2016b). Such an approach is much less conservative than masking all the sources detected in intensity, but an important caveat applies: it is a non-blind approach and sources with flux densities just below the detection level in intensity but with very large polarization fractions can be missed. Therefore, blind searches in the $P$ maps would be desirable, and perhaps even extended to $E$ and $B$ maps given that some component separation methods work directly with these quantities.

Secondly, we can seek to improve the source detection process using as a proxy deep catalogues of compact sources detected with high resolution telescopes or interferometers in the microwave and sub-mm frequency range to either mask or enable non-blind searches in the CMB data sets. We can motivate this as follows.

As mentioned in Sect. 2.2, unresolved extragalactic sources contribute both shot noise (Poisson) and clustered power to the CMB maps. There is a substantial difference between the contributions of the two main populations, radio sources and dusty galaxies. The former, that dominate the point source temperature fluctuations up to 100–150 GHz, have a low clustering power, so that their power spectrum is simply white noise ($C_{\ell}=\hbox{constant}$). The amplitude of the Poisson power spectrum is

[TABLE]

where $dN(S)/dS$ are the differential number counts per steradian of sources weaker that the detection limit $S_{d}$. The slope $\beta$ of the radio source $dN(S)/dS$ ($dN(S)/dS\propto S^{-\beta}$) is always $<3$, implying that the maximum contribution comes from sources just below the detection limit. Thus the radio source contamination can be mitigated, at least in some sky areas, by masking sources detected by large area ground based surveys, such as the Australia Telescope 20 GHz survey, covering the full southern sky down to flux densities in the range 40–100 mJy. Over areas of $\sim 1,000\,\hbox{deg}^{2}$ one can take advantage of the deep SPT (at 95 and $\simeq 150$ GHz) and ACT (at $\simeq 150$ GHz) surveys.

Moreover, in the future, significant progress will be made possible by the Stage-IV CMB (CMB-S4) surveys (Abazajian et al. 2016) which aim to map most of the sky to sensitivities of $\sim 1\,\mu$K$\cdot$arcmin. The counts based on the combination of the SPT and ACT surveys with those from Planck (Planck Collaboration XIII 2011; Planck Collaboration VII 2013), that extend them to bright flux densities, tightly constrain the amplitude of the power spectrum, making possible a correction of the CMB contamination at the power spectrum level.

However, it should also be noted that a large fraction of the radio sources that can be detect in CMB experiments are variable. The fact that a source was detected with a high resolution instrument at any given time does not imply that it is contributing to the signal as measured by a CMB experiment measured at a different epoch. Conversely, some sources might be much brighter in the space mission data than when observed from the ground. One way of addressing this problem is to plan follow-up observations of sources with high resolution instruments in such a way that these observations are performed near simultaneously with the CMB observations. Such a strategy was already undertaken by Planck, and involved the simultaneous observation of radio sources using the Australian Compact Array. Something similar could be envisaged for future CMB experiments, although ground based observatories are often heavily over-subscribed. This has led some CMB projects to invest in the contruction of dedicated instruments to perform simultaneous follow-up observations of compact sources, e.g., the Very Small Array at Teide Observatory.

The situation is substantially different in the case of unresolved dusty galaxies. The main contribution to their power spectrum comes from faint, high surface density, high-$z$ galaxies. In fact, the amplitude of the Poisson power spectrum can be determined quite accurately exploiting Planck plus ground based surveys: the SCUBA-2 counts at 350 GHz (Hsu et al. 2016; Geach et al. 2017, and references therein), SPT and ACT counts at $\simeq 150$ GHz and $\simeq 200$ GHz, extended to $\mu$Jy levels by ALMA surveys (Carniani et al. 2015; Aravena et al. 2016, and references therein). But their power spectrum is dominated by clustering for \ell\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}2000 (Planck Collaboration XXX 2014; De Zotti et al. 2015). Planck has allowed measurements of the power spectrum of unresolved point sources in the multipole range dominated by clustering at 217 GHz, 353 GHz, and 545 GHz, with some preliminary estimate also at 100 GHz (Planck Collaboration XVIII 2011; Planck Collaboration XXX 2014; Mak et al. 2016). Herschel (Viero et al. 2013), ACT (Dunkley et al. 2013) and SPT (Reichardt et al. 2012) measurements have extended the determinations to smaller angular scales. It should be mentioned that despite the low level of polarization of individual dusty galaxies, the actual polarization of the CIB emission does not average out in variance since this measures fluctuations in power.

As in the case of the temperature power spectrum, if indeed dusty galaxies do not contribute much, the radio source contamination can be mitigated exploiting deep, large area surveys. Furthermore, the counts in polarization that CORE will provide should allow accurate modelling of the contamination at the power spectrum level.

6.4 Further improvements

The work done for this paper shows the difficulty of reaching precision and accuracy for the detection of inflationary $B$ modes if $r\simeq$ few $\times 10^{-3}$. A big risk for a future space mission is to detect an excess of $B$-mode power on large scales, without the capability to ascertain its cosmological origin. Such a false detection would reproduce the uncertainty that followed the original BICEP2 announcement with ground-based observations.

To avoid this, the future space mission must be designed not only for the capability to reduce foreground contamination to allow for the detection of primordial $B$-modes. It must also be designed to provide means of evaluating the level of residual foreground emission in $B$-mode maps and spectra.

To this effect, the capability of comparing the results obtained on different patches of sky, with different band subsets, with different component separation methods is key.

Component separation to detect $r\simeq 10^{-3}$ has been shown to be hard. We now identify paths for further improvement of the component separation, and for the characterization of errors for the three main methods used in this work.

Complementary observations:

If the model assumed by Commander is accurate enough to precisely hold over a large range of frequencies, then additional observations might help. We anticipate that by the CORE launch-date external surveys of polarized sources will have accurately measured their spectral properties, and that better constraints on the spectral index of the diffuse synchrotron polarization emission over the full sky will have been achieved by C-BASS (Irfan et al. 2015) and QUIJOTE (Rubiño-Martín et al. 2012) and follow-up surveys. CORE would be the most sensitive experiment to precisely map dust emission in the relevant frequency range, although an experiment such as PIXIE (Kogut et al. 2016), with 600 absolutely-calibrated frequency channels from 30 to 6000 GHz, could bring additional information on large angular scales to better model large-scale galactic foreground emission.

Commander:

In this simulation, to alleviate the issue discussed in Sec. 6.2, we have applied Commander on maps generated directly at $N_{\rm side}=16$. The analysis of real sky data will require to find a way to mitigate the impact of variations of emission laws within sky pixels and along the line-of-sight. This difficulty, which has been bypassed in the present work, will have to be addressed in the future.

As a first step, we must investigate ways to extend the Commander parametric fit for $B$-modes to smaller pixels and larger multipole ranges, e.g. $2\leq\ell\leq 200$, by optimizing the algorithm in some way to relax the actual large computational costs. We must also find a way to validate within the analysis the parameterization that is assumed for each of the foreground emissions.

SMICA:

The SMICA method relies on the effective dimension of the foreground subspace in polarization (Sect. A.3), which is a strong prior assumption, to extract the cosmological signal. In the current SMICA algorithm, the dimension of the foreground subspace is fixed to an ad-hoc value, although this is likely to vary both over the sky and the angular scales probed, depending on the relevance of each of the foregrounds with respect to the noise. A possible improvement would be to follow the approach adopted by the GNILC component separation method, developed for intensity in Remazeilles et al. (2011b) and Planck Collaboration XLVIII (2016). Here, any prior assumption on the number of foregrounds is relaxed and the effective dimension of the foreground subspace is estimated locally both over the sky and over physical scales by using wavelet decomposition of the sky maps and thresholding with respect to the local foreground signal-to-noise ratio.

NILC:

In the present paper, a standard NILC algorithm has been used. There has been no attempt to optimize spectral windows, nor to optimize the rejection of foreground contamination, rather than minimizing the total variance of the reconstructed maps. Such improvements may prove useful to reduce the foreground contamination in the reconstructed maps, and possibly to find ways of evaluating residual emission. Further investigations are left to future work on this topic.

7 Conclusions

In this paper, we have considered whether the CORE satellite mission, as proposed to ESA’s M5 call, allows the removal of astrophysical foregrounds at the accuracy required to measure the power spectrum of primordial CMB $B$-modes on both the reionization and recombination scales for inflationary models with tensor-to-scalar ratios over a range of interesting values, an optical depth to reionization, $\tau=0.055$ (Planck Collaboration XLVI 2016) and the other $\Lambda$CDM model cosmological parameters defined by the Planck best-fit values (Planck Collaboration I 2016).

In particular, three independent component separation methods of either a parametric (Commander) or blind (SMICA and NILC) nature were applied to a set of simulations consistent with observations of the sky by CORE. These simulations include foregrounds covering a range of complexity as modelled by the PSM, but at least including contributions due to synchrotron and dust with variable spectral indices, with additional components due to AME and point-sources considered in the most challenging cases. The impact of gravitational lensing effects has also been tested by comparing results either without the addition of a lensing $B$-mode contribution, or including them at both the predicted level or after an idealized 60% delensing has been applied.

We have determined that, even at the current stage of development more than one decade before the potential launch of the CORE satellite, each of the three component separation methods are independently able to reconstruct the cosmological $B$-mode signal, without bias, for tensor-to-scalar values $r\gtrsim 5\times 10^{-3}$. Specifically, in the absence (or perfect control) of lensing CORE allows the determination of a tensor-to-scalar ratio at the $10^{-2}$ or $5\times 10^{-3}$ level at a significance exceeding $10\sigma$. However, when lensing effects are included, the significance of detection at the latter amplitude is lowered to about $4\sigma$, with, for example, the Commander estimate yielding a value of $r=\left(5.4\pm 1.5\right)\times 10^{-3}$. Nevertheless, one of the methods – Commander – operating in the low-$\ell$ regime only suggests that a value of $r$ as low as ${r=2.5\times 10^{-3}}$ can be determined without bias and with $\sim 8\sigma$ significance in the absence of lensing. Therefore, it appears that, at the very least, observations by CORE would allow constraints to be imposed on the energy-scale of inflation for the Starobinsky $R^{2}$ inflation model, for which $r\sim 4.2\times 10^{-3}$ (Starobinsky 1980; Mukhanov and Chibisov 1981; Starobinskij 1983).

Importantly, however, we have demonstrated that when the tensor-to-scalar ratio is decreased to a critical value of $r\sim 10^{-3}$, then the spectral complexity of the foregrounds is sufficient to bias the measurement of the primordial CMB $B$-mode power spectrum, both on reionization scales ($\ell<12$) in the case of parametric component separation techniques, and over all angular scales when utilizing blind methods, even in the absence of lensing effects. What is evident is that, in order to measure $r=10^{-3}$, foreground residuals after component separation have to be modelled and marginalized over in the likelihood estimation, in which case CORE is able to correct for the bias and measure the tensor-to-scalar ratio with an accuracy $\sigma(r=10^{-3})\sim 2\times 10^{-4}$, in the absence of lensing effects. However, the presence of lensing $B$-modes, together with the addition of further foreground contributions due to AME and point-sources, decreases the sensitivity to $r$ to $\sigma(r=10^{-3})\sim 5\times 10^{-4}$. We note that Alonso et al. (2017) have claimed similar uncertainties on $r$ for a high-resolution ground-based CMB-Stage IV experiment, although the authors do not include a contribution from extra-galactic sources. In addition, and perhaps unexpectedly, even if delensing of the $B$-mode signal is achieved at the 60% level, the uncertainty on the tensor-to-scalar ratio remains at this level, despite the fact that the cosmic variance contribution from the lensing $B$-modes must have been reduced by a corresponding amount. The uncertainty is therefore dominated by foreground residuals rather than lensing cosmic variance. In this case, the use of external information on the spectral index variation of the diffuse synchrotron emission and the polarized point-source contribution should increase the significance of the detection by CORE to $\sigma(r=10^{-3})\sim 3\times 10^{-4}$, after foreground cleaning and 60% delensing (Sect. 6.1).

We have identified different sources of potential bias in the reconstruction of the primordial CMB $B$-mode signal: (i) incorrect foreground modelling; (ii) inadequate frequency range or lack of frequency channels sampling it; and (iii) systematic bias due to averaging foreground spectral indices through pixelization. While these sources of bias likely have no significant impact on the reconstruction of CMB temperature and $E$-mode polarization, they become critical for the reconstruction of the faint primordial CMB $B$-mode polarization especially if $r\sim 10^{-3}$. Despite the absence of spectral mismatch between the foreground model and the data, and a perfect scaling of the foregrounds in each pixel area through direct $N_{\rm side}=16$ simulations, we have found that the spectral variations of the synchrotron spectral index are not accurately reconstructed, even though the fit of the total sky emission in polarization can be considered satisfactory on the basis of low chi-square statistics. This, ultimately, generates some bias on the estimate of a tensor-to-scalar ratio at the level of $r=10^{-3}$. However, fixing the synchrotron spectral indices and the dust temperatures to their input values and fitting the dust spectral index only with Commander helps to remove the bias on the reconstructed CMB $B$-mode power spectrum on reionization scales, and reduces the overall uncertainty on $r$ after foreground cleaning and delensing (right panel of Fig. 28). Therefore, CORE could benefit from the addition of low-frequency channels below $60$ GHz. This could be achieved in part by including information from full-sky polarization surveys at lower frequencies, such as C-BASS (Irfan et al. 2015) and QUIJOTE (Rubiño-Martín et al. 2012). A counterpoint to this, however, is that we would then need accurate modelling of diffuse AME polarization, which can no longer be ignored as a contaminant to low-$r$ CMB $B$-modes at such frequencies.

The high resolution of the CORE frequency channels also allows the determination of the CMB $E$-mode polarization signal up to relatively high multipoles. This, in principle, will then allow the accurate reconstruction of the lensing potential, $\phi$, using well-known quadratic lensing estimators (Hu and Okamoto 2002) as applied to the reconstructed CMB $E$-mode polarization map. Subsequently, proper delensing of the CMB $B$-mode power spectra could reduce the lensing contribution to the cosmic variance by $\sim 60$%, enabling the detection of the recombination bump for a level of $r$ that is out of reach for lower-resolution CMB polarization experiments such as LiteBIRD (Matsumura et al. 2016) or PIXIE (Kogut et al. 2016).

It should be noted that systematic effects (calibration errors, beam asymmetries, bandpass mismatch, etc) will have an impact on component separation results and on the accuracy of foreground removal. Dick et al. (2010) demonstrated that small calibration errors (less than $1$%) can perturb the weighting of the frequency maps implemented in the NILC component separation method, resulting in a partial cancellation of the reconstructed CMB signal, especially for the high signal-to-noise measurements. Beam asymmetries can also be regarded as an imperfect calibration for NILC, and also affect the inferred NILC weighting similarly. Given the high sensitivity required to measure the primordial $B$-modes at a level of $r\sim 10^{-3}$, an accurate calibration of the instrument becomes essential. Detector bandpass mismatch can also impact the results of parametric fitting methods, such as Commander, due to the modication of the foreground spectra integrated over the bandpass from one detector to another. The control of such systematics effects is discussed in a companion paper (Natoli et al. 2017 ). However, ideally the effects of foregrounds and systematics actually should be controlled collectively rather than independently. This will be investigated in future work.

In summary, we have performed a detailed component separation and likelihood analysis, going beyond a simple Fisher approach, to reconstruct the exact shape of the primordial CMB $B$-mode power spectrum as might be measured by CORE in order to study both sensitivities and biases to the tensor-to-scalar ratio measurements after foreground cleaning. In Sect. B, we have compared our component separation results with a forecasting method that performs an effective statistical averaging of CMB and noise realisations and foreground uncertainties, and we have found consistency between both approaches when analysing data with a tensor-to-scalar ratio of $r=10^{-3}$ in the presence of foregrounds and lensing. Although additional refinements are required, the component separation results seem promising and we can speculate that a detection of $B$-mode polarisation corresponding to $r\sim 10^{-3}$ will be achievable by a CORE-like mission with a potential launch-date around 2030.

Acknowledgments

The research leading to these results has received funding from the ERC Starting Consolidator Grant (no. 307209). Some of the results in this paper have been derived using the HEALPix package (Górski et al. 2005). We acknowledge the use of the PSM package (Delabrouille et al. 2013) , developed by the Planck working group on component separation, for making the simulations used in this work. We acknowledge the use of the Ulysses cluster at SISSA. This research was partially supported by the RADIOFOREGROUNDS project, funded by the European Commission’s H2020 Research Infrastructures under the Grant Agreement 687312, and the INDARK INFN Initiative. GDZ acknowledges support by ASI/INAF agreement no. 2014-024-R.1. R.F.-C., E.M.-G., and P.V. acknowledge support from the Spanish MINECO project ESP2015-70646-C2-1-R (cofinanced with EU FEDER funds), Consolider-Ingenio 2010 project CSD2010-00064 and from the CSIC “Proyecto Intramural Especial” project 201550E091. JGN acknowledges financial support from the Spanish MINECO for a ’Ramon y Cajal’ fellowship (RYC-2013-13256) and the I+D 2015 project AYA2015-65887-P (MINECO/FEDER). CJM is supported by an FCT Research Professorship, contract reference IF/00064/2012, funded by FCT/MCTES (Portugal) and POPH/FSE.

Appendix A Component Separation Methods

A.1 Bayesian Parametric Fitting

A physical model is fit to a set of observations within a Bayesian parametric framework using the Commander algorithm (Eriksen et al. 2008). The sky model, $\boldsymbol{m}(\nu,p)$, as fitted to the CORE sky maps, $\boldsymbol{d}(\nu,p)$, at each frequency, $\nu$, and for each pixel, $p$, is parametrized by

[TABLE]

where $\boldsymbol{s}^{cmb}(p)$ is the amplitude of the CMB ${Q,U}$ polarization anisotropies, $\boldsymbol{s}^{sync}(p)$ is the amplitude of the polarized synchrotron ${Q,U}$ radiation, $\boldsymbol{s}^{dust}(p)$ is the amplitude of the polarized thermal dust ${Q,U}$ radiation, $\boldsymbol{n}(\nu,p)$ is the instrumental noise for CORE in the ${Q,U}$ Stokes parameters, $a(\nu)$ is the frequency spectrum of the CMB, $\beta_{s}(p)$ is the synchrotron spectral index, $\beta_{d}(p)$ is the thermal dust spectral index, and $T_{d}(p)$ is the dust temperature. Here, we fix $a(\nu)$ to be constant and equal to unity in thermodynamic temperature units. In addition to the amplitudes and spectral indices of the components, the $E$- and $B$-mode angular power spectra $C_{\ell}=\left\{C_{\ell}^{EE},C_{\ell}^{BB}\right\}$ of the CMB are fitted in a self-consistent way with the amplitude and spectral index maps using the Gibbs sampling scheme detailed later in the section. Considering both the $Q$ and $U$ Stokes parameters in the analysis, there are $2\times 19=38$ channels from CORE, while the parametric model consists of $2\times 3+3+2=11$ parameters in total to be fitted to the data,

[TABLE]

As shown in Remazeilles et al. (2016), the performance of the Commander component separation method relies on the accurate parametrization of the spectral properties of the foregrounds, and small foreground modelling errors can strongly bias an estimate of a low tensor-to-scalar ratio. In this work, we fit a power-law for synchrotron and a single modified blackbody for thermal dust, thereby avoiding any mismatch between the data and the model, in order to characterize the uncertainty on the tensor-to-scalar ratio that is only due to foreground residuals after component separation. It should be noted that we have verified that omitting AME in the parametric model has no impact for Commander on the estimate of the tensor-to-scalar ratio, so that we do not attempt to fit for its contribution. We have determined that this absence of impact on the results is due to the frequency coverage of CORE. The CNM spectrum of AME steepens so sharply at frequencies $\gtrsim 50$ GHz that AME polarization becomes insignificant in the lowest frequency band, $60$ GHz, of the instrument.

Component separation is achieved by computing the joint CMB-foreground posterior distribution, which, according to Bayes’ theorem, is given by

[TABLE]

where $P\left(\boldsymbol{s},\boldsymbol{\beta},C_{\ell}\right)$ is the prior distribution of the parameters. Here, we do not set any prior on $C_{\ell}$ and on the amplitudes of the components but we assume Gaussian priors for the foreground spectral parameters, $P\left(\beta_{s}\right)\sim\mathcal{N}(-3,0.1)$, $P\left(\beta_{d}\right)\sim\mathcal{N}(1.6,0.3)$, and $P\left(T_{d})\sim\mathcal{N}(19.4\,{\rm K},1.5\,{\rm K}\right)$. The Gaussian priors are multiplied by a Jeffreys ignorance prior in order to suppress the prior volume of the likelihood space for non-linear parameters (Jeffreys 1946; Eriksen et al. 2008). By assuming the noise to be Gaussian and uncorrelated across frequencies, the likelihood is then

[TABLE]

where ${\boldsymbol{\rm N}}$ is the noise covariance matrix in pixel space.

The Bayesian parametric fitting approach enables the end-to-end propagation of the foreground uncertainties to the reconstructed CMB power spectrum and then to an estimate of the tensor-to-scalar ratio. The algorithm also provides a map of the chi-square goodness-of-fit, measuring the mismatch between the model and the data in each pixel,

[TABLE]

where $\sigma_{\nu}(p)$ is the noise pixel $p$. The chi-square map provides useful feedback on the fidelity of the foreground modelling as well as a possible criterion to improve masking a posteriori, by discarding the pixels where the chi-square statistic is too high in successive iterations of the component separation. However, in this paper we use a simpler masking approach, as described in Sect. 4.2.

Direct computation of the joint CMB-foreground posterior distribution Eq. (A.6) is not feasible as it would require the distribution to be mapped out over a multidimensional grid, the size of which grows exponentially with the number of parameters. A more tractable approach is to map out the joint CMB-foreground posterior distribution by Gibbs sampling (Wandelt et al. 2004; Eriksen et al. 2004b, 2008), i.e., each parameter is sampled alternately from iterative conditional probabilities, according to the following Gibbs sampling scheme,

[TABLE]

where $(i)$ denotes the samples at the $i^{th}$ iteration in the Markov chain. It has been demonstrated mathematically that as the number of iterations approaches infinity then the Gibbs scheme converges to the sampling with the full joint CMB-foreground posterior distribution, after some initial burn-in (Wandelt et al. 2004). Here, we use $2000$ iterations in the Markov chain but discard the first $500$ Gibbs samples (burn-in) to ensure that we have reached convergence.

The amplitudes of the components are sampled by the conditional Gaussian distribution

[TABLE]

where ${\boldsymbol{\rm S}}$ and ${\boldsymbol{\rm N}}$ are the CMB and noise covariance matrices in pixel space, and $\widehat{\boldsymbol{s}}$ is the Wiener filter solution:

[TABLE]

In practice, the Wiener estimate is computed by solving the following equation with the conjugate gradient method,

[TABLE]

where $w_{0},w_{1}\sim\mathcal{N}(0,1)$. By assuming the CMB fluctuations to be Gaussian and isotropic, the CMB power spectrum is thus sampled from the conditional distribution,

[TABLE]

which corresponds to an Inverse-Gamma distribution when interpreted as a function of $C_{\ell}$.555In the case of polarization, the conditional distribution actually corresponds to a mix of an Inverse-Wishart distribution (the correlated T and E-mode parts) and an Inverse-Gamma distribution (the B-mode part). By writing the pseudo-power spectrum of the CMB map as

[TABLE]

the $C_{\ell}$ sampling distribution can be recast as

[TABLE]

For the sampling of the foreground spectral indices, $\boldsymbol{\beta}$, we refer to Eriksen et al. (2008) for more details.

The posterior distribution of the CMB $E$- and $B$-mode power spectra can then be estimated by using the Blackwell-Rao approximation (Chu et al. 2005),

[TABLE]

where the sum runs over $N_{G}$ Gibbs samples $\widehat{C}_{\ell}^{(i)}$, and P\left(C_{\ell}\big{|}\widehat{C}_{\ell}^{(i)}\right) is given by Eq. (A.19). The Blackwell-Rao estimate becomes an exact approximation of the posterior $C_{\ell}$ distribution as the number of Gibbs samples increases (Chu et al. 2005; Eriksen et al. 2008).

In summary, the Commander algorithm enables component separation, power spectrum estimation, and cosmological parameter estimation (see Sect. 4.2) in a self-consistent way using MCMC Gibbs sampling and the Blackwell-Rao approximation.

Note that in Sect. 5, where we integrate low- and high-$\ell$ results for an optimal estimation of $r$ based on a larger range of multipoles, we multiply the posteriors for the cosmological parameters obtained from Commander, using a Blackwell-Rao approximation and non-Gaussian likelihood, with those obtained from either NILC or SMICA, using a Gaussian likelihood. This is because the latter are not suited to the Blackwell-Rao approach.

A.2 Internal Linear Combination in Needlet space

As an alternative approach for the removal of foregrounds present in the sky maps, we implement a blind Internal Linear Combination (ILC) method (Tegmark and Efstathiou 1996; Bennett et al. 2003; Tegmark et al. 2003; Eriksen et al. 2004a). This method only attempts to reconstruct the CMB signal, without using any prior information about foregrounds. It is based on two specific assumptions. Firstly, that the CMB is frequency independent in thermodynamic unit, and secondly, that the CMB is uncorrelated with foreground signals. The ILC method then estimates the CMB, ${\widehat{S}}$, as a weighted linear combination of the set of input multi-frequency sky maps such that the variance of the estimate is minimum, with unit response to the flat CMB frequency spectrum,

[TABLE]

where $X$ is the vector of frequency maps, $a$ the constant frequency spectrum of the CMB signal $S$, $F$ the total foreground signal, $N$ the instrumental noise for the different frequency channels, and ${\widehat{R}}$ the frequency-frequency covariance matrix. The first condition guarantees minimum contamination by foregrounds and instrumental noise whereas the second condition guarantees that the CMB signal is conserved without bias. The presence of the foregrounds induces correlated errors across frequencies, so that the ILC weights adjust themselves to minimise the foreground residuals present in the the weighted linear combination. However, in reality the weights result from a trade-off between minimising the foregrounds and minimising the instrumental noise contribution in the reconstructed CMB map.

However, there are also some drawbacks to the method. As discussed by a number of authors (Saha et al. 2006; Souradeep et al. 2006; Hinshaw et al. 2007; Saha et al. 2008; van de Weygaert and Schaap 2009; Delabrouille et al. 2009; Saha 2011), the component of interest (CMB) and the foreground signals must be uncorrelated for proper ILC performance. On finite data sets, this can only be approximately true, and empirical correlations between the CMB and foregrounds generate a bias in the reconstructed CMB on large angular scales. In addition, as shown by (Dick et al. 2010), the ILC method tends to amplify calibration errors such that the reconstructed CMB map exhibits significantly lower variance than the true sky with strong suppresssion of CMB features, in particular in the high signal-to-noise regime.

The ILC method can be straightforwardly implemented in either real (pixel) space or in harmonic space. Thus, sets of ILC weights can either be computed for different regions of the sky or for different angular scales, respectively, which allows for variations of the data covariance matrix in either space. However, the ILC in harmonic space does not take into account the fact that noise can be a significant source of CMB measurement error in at high Galactic latitude, while foreground signals are more important at low Galactic latitude. Conversely, the ILC in pixel space does not take into account the fact that the noise dominates on small angular scales, while diffuse Galactic foreground emission dominates on large angular scales.

In order to overcome this problem, we implement the ILC on a frame of spherical wavelets called needlets (Narcowich et al. 2006), a component-separation approach that we now refer to as the Needlet Internal Linear Combination (NILC) method. This technique has already been applied broadly in CMB data analysis (Delabrouille et al. 2009; Remazeilles et al. 2011a, b; Basak and Delabrouille 2012, 2013; Remazeilles et al. 2013). The needlets enable localized filtering in both pixel space and harmonic space because they have compact support in the harmonic domain, while still being very well localized in the pixel domain (Narcowich et al. 2006; Marinucci et al. 2008). The needlet decomposition allows the ILC weights to vary both smoothly on large angular scales and rapidly on small angular scales, which is not possible by sub-dividing the sky into different areas prior to any processing.

The needlet windows, $h^{j}_{l}$, in harmonic space that we use in our analysis are defined as follows,

[TABLE]

In terms of $h_{l}^{j}$, the spherical needlets are defined as

[TABLE]

where the $\{\xi_{jk}\}$ denote a set of cubature points on the sphere for scale $j$. In practice, we identify these points with the pixel centres of the HEALPix pixelization scheme (Górski et al. 2005). Each index $k$ corresponds to a particular HEALPix pixel, at a resolution parameter $N_{\mathrm{side}}(j)$ specific to that scale $j$. The cubature weights $\lambda_{jk}$ are inversely proportional to the number $N_{j}$ of pixels used for the needlet decomposition, i.e., $\lambda_{jk}=\frac{\displaystyle{4\pi}}{\displaystyle{N_{j}}}$.

Given a set of needlet functions, any sky map of a spin-[math] field $X(\hat{n})$ (such as the CMB temperature anisotropy, or the $E$- and $B$-modes of CMB polarization) on the sphere can be expressed as

[TABLE]

where the needlet coefficients, $\beta^{X}_{jk}$, of the sky map are denoted as

[TABLE]

For each scale $j$, the ILC filter has a compact support between the multipoles $l^{j}_{min}$ and $l^{j}_{max}$ with a peak at $l^{j}_{peak}$ (Fig. 30) The needlet coefficients, $\beta^{X}_{jk}$, are computed on the HEALPix grid points, $\xi_{jk}$, with a resolution parameter, $N_{\mathrm{side}}$  equal to the smallest power of $2$ larger than $l^{j}_{max}/2$.

Once an estimate of the CMB sky has been evaluated, the corresponding CMB angular power spectrum can be evaluated using a a pseudo-$C_{\ell}$ estimator (Hivon et al. 2002; Chon et al. 2004; Szapudi et al. 2005). This method is computationally much faster than maximum likelihood and provides optimal results at intermediate to high $\ell$s. In order to compute the covariance on our measurement of angular power spectrum, we have followed the method described in Tristram et al. (2005).

In this paper, although the NILC weights are computed from full-mission sky maps, the impact of instrumental noise residuals on the measurement of the angular power spectra is avoided by producing independent CMB maps for the two half-mission data sets. These maps are obtained by applying the full mission NILC weights to the half-mission sky maps. Each data point of the angular power spectra is then obtained from the average of all possible cross-half-mission angular power spectra.

Suppose we have $N$ measurements (one per half-mission) of CMB fields such that the residuals of noise in these measurements are statistically independent from each other. The estimator of angular power spectra with $(2l+1)$ weights for binning is given by

[TABLE]

where $l_{min}^{b}$ and $l_{max}^{b}$ are respectively the minimum and maximum values of the multipole for the b-th bin, and $\widehat{D}^{XY}_{l}$ is the average of all possible cross-half-mission spectra for the recovered CMB fields after applying the mask under consideration:

[TABLE]

The corresponding coupling matrix $M^{XY}_{ll^{\prime}}$, in terms of the angular power spectrum $W^{XY}_{l}$ of the mask, is given by

[TABLE]

The covariance $\Xi^{XY,XY}_{ll^{\prime}}$ of $\widehat{C}^{XY}_{l}$ is by definition:

[TABLE]

where $\Sigma^{XY,XY}_{l_{1}l_{2}}$ is the covariance of masked angular spectra $\widehat{D}^{XY}_{l}$:

[TABLE]

The coupling matrix $M^{(2)}_{ll^{\prime}}$, in terms of the angular power spectrum, $W^{XY,X^{\prime}Y^{\prime}}$, of the product of the masks for the fields $X$ and $Y$ is given by

[TABLE]

[TABLE]

The covariance of the estimator of the binned angular power spectra, $\widehat{C}^{XY}_{b}$, is then given by

[TABLE]

A.3 Spectral Matching

The use of the Spectral Matching Independent Component Analysis (SMICA) method (Delabrouille et al. 2003) assumes that the total sky emission observed across a set of frequency bands, in any pixel $p$, or any harmonic mode $(\ell,m)$ is a noisy sum of components, that can be written, in the most general form, as,

[TABLE]

where the sum runs over all components of the sky emission, and where the contribution of each component to the total observed sky emission is fully described by a vector $\mathbf{c}$ that represents its emission in the set of all observed frequency bands. For each pixel or harmonic mode, all of $\mathbf{x}$, the various components $\mathbf{c}$ and $\mathbf{n}$ are vectors of dimension $n_{\nu}$.

While no assumption has been made so far about the properties of the components, such a decomposition is of particular interest when the components are independent. Indeed, together with the assumption of independence between the sky components and noise, this guarantees that the multivariate power spectrum of $\mathbf{x}$ is the sum of the multivariate power spectra of the components and of the noise, with no cross terms.

[TABLE]

where the sum runs over independent sky components, and for each $\ell$, each of ${\boldsymbol{\rm X}}_{\ell}$, ${\boldsymbol{\rm N}}_{\ell}$ and of all of the ${\boldsymbol{\rm C}}_{\ell}$ is an $n_{\nu}\times n_{\nu}$ covariance matrix.

The independence of components depends on their physical origin. CMB emission is independent from Galactic foreground emission, however, the spinning dust emission in our Galaxy can not be assumed to be independent from the thermal dust emission.

We now turn to modeling further the emission of each physical component, $\mathbf{c}$, as a linear mixture of $n\leq n_{\nu}$ independent, unphysical, templates.

[TABLE]

where ${\boldsymbol{\rm A}}$ is an $n_{\nu}\times n$ matrix, and $\mathbf{s}$ a set of $n$ templates that describe component $\mathbf{c}$ (hereafter denoted as the sources of component $\mathbf{c}$). Such a decomposition is always possible independently of the nature and of the physical properties of the component: in the worst case scenario, the emission of a component can be modeled using as sources the $n=n_{\nu}$ templates of emission of the component in all the bands of observation. The matrix ${\boldsymbol{\rm A}}$ is then the identity matrix. However, this decomposition into sources is more interesting when $n<n_{\nu}$. We denote as the dimension of a component the minimum number of sources required to represent with sufficient accuracy666So that the difference between the real emission of the component and its model with $n<n_{\nu}$ sources is not detectable given the observational noise. its emission in all of the observed bands. The extreme case of $n=1$ corresponds to a component that scales perfectly with frequency, a model valid for the CMB emission. For example, Galactic ISM emission, considered as one single component, typically requires $n>1$ templates that may or may not be identified with emission from specific processes such as free-free or synchrotron.

The multivariate power spectrum ${\boldsymbol{\rm C}}_{\ell}$ of an $n$–dimensional component can be written as

[TABLE]

where ${\boldsymbol{\rm S}}_{\ell}$ is the multivariate power spectrum of the sources of component (for each $\ell$, an $n\times n$ matrix). The diagonal elements of ${\boldsymbol{\rm S}}_{\ell}$ represent the elements of the auto-spectra and the off-diagonal elements the cross-spectra of the sources $\mathbf{s}$. For a one-dimensional component, ${\boldsymbol{\rm S}}_{\ell}$ is for each $\ell$ a $1\times 1$ matrix, that represents the value of the power spectrum of that particular component for multipole $\ell$.

The idea of the SMICA method is to adjust a model of the covariance matrices of the various independent components and of the noise to achieve the best match, in the maximum likelihood sense, to the observed covariance matrices of the data. Observed covariances are estimated from the observations using spectral band averages of the form

[TABLE]

where $q$ indexes the spectral band, and ${{\mathcal{N}}_{q}}$ is the number of modes $(\ell,m)$ in the band. In practice, one must take into account the fact that the original maps are in general at a different resolution, and that a mask must be applied to select only regions where the Galactic emission is low and simple enough to be accurately modeled as a component of sufficiently low dimensionality for the spectral fit to be possible.

The model that we adjust comprises in general three sky emission components, and noise

• •

CMB, modeled as a 1-dimensional component and specified by its frequency dependence;

• •

a multi-dimensional ‘catch-all’ component that models the emission of the Galactic ISM as well as the diffuse background of extra-galactic sources;

• •

an additional 1-dimensional component that models the emission of strong extra-galactic sources;

• •

a noise component for which the covariance matrix is fixed to its exact value.

For the data, at each frequency a $B$-mode map is generated from the $Q$ and $U$ Stokes parameter data. A spherical harmonic transform is computed up to $\ell_{\mathrm{max}}=6000$, and then the corresponding $a^{B}_{\ell m}$ are convolved with a Gaussian beam of 20 arcminutes FWHM, before synthesising the $B$-mode map. These maps are then masked before we finally compute the weighted covariance matrices that SMICA will adjust to the model.

The statistical properties of the noise are assumed to be known a priori (i.e., the term ${\boldsymbol{\rm N}}_{\ell}$ in Eq. A.46 is known). Alternatively, a parametric model can be assumed for a joint estimation of sky model and noise model parameters. Generically, the set of model parameters $\theta$ that we assume comprises the set of band-average auto and cross spectra of all sources of sky emission, and the elements of the ‘mixing matrix’ ${\boldsymbol{\rm A}}_{\rm fg}$ of the ‘catch all’ component of foreground sky emission. The adjustment criterion is to maximize the likelihood of the model given the observed band-averaged covariance of the observations. The best fit estimated parameter set $\widehat{\theta}$ is obtained as :

[TABLE]

where

[TABLE]

and where $D(\cdot,\cdot)$ is a measure of divergence between two positive $n\times n$ matrices defined by

[TABLE]

We refer the reader to the SMICA publications (Delabrouille et al. 2003; Cardoso et al. 2008) for further details about the method.

Appendix B Forecasts for the CORE component separation problem

We describe in this section a complementary approach to the component separation methods applied on single CMB and noise realizations as presented in the main sections of this article. The framework, named xForecast (Stompor et al. (2016)), is an extension of the CMB4cast777Publicly accessible at http://portal.nersc.gov/project/mp107/index.html method described in Errard et al. (2011); Errard and Stompor (2012); Errard et al. (2016).

xForecast optimizes a CMB- and noise-averaged spectral likelihood, therefore providing an estimate of the ensemble-averaged spectral parameters, the statistical and systematic foreground residuals, and the likelihood for the tensor-to-scalar ratio, $r$. Whereas CMB4cast allows the estimation of the impact of statistical residuals on $\sigma(r)$, xForecast derives the statistically-meaningful performance of CORE in the specific case of a parametric component separation approach, in particular wih respect to the bias on the estimation of $r$.

B.1 Formalism

Similarly to the formalism used by Commander (see Sect. 4.2), we use the parametric maximum-likelihood approach as introduced in, e.g., Brandt et al. (1994); Eriksen et al. (2006); Stompor et al. (2009b). At a sky pixel $p$, the measured amplitudes at all frequencies are concatenated in a data vector $d$, such that

[TABLE]

where

• •

$\mathbf{A}$ is the mixing matrix, which contains the scaling laws of all sky components (CMB, foregrounds). Under the parametric formalism, we assume that the mixing matrix $\mathbf{A}$ can be parametrized by a set of spectral parameters $\beta$:

[TABLE]

• •

$s^{\rm true}(p)$ contains the true amplitudes of the sky signals, scaled at a reference frequency;

• •

$n(\nu,p)$ is the instrumental noise, assumed white in this analysis.

We will not write the $\nu$ argument in most of the equations below. Given Eq. B.1, the component separation and cosmological analysis is performed in three steps:

• •

the estimation of the mixing matrix or, equivalently, the estimation of the spectral parameters. This is achieved through the optimization of a spectral likelihood, $\mathcal{L}_{spec}(\beta)$, as detailed in Stompor et al. (2009b). In order to estimate the statistically-averaged performance of the component separation for a given instrumental configuration, Stompor et al. (2016) propose a spectral likelihood averaged over a statistical ensemble of the noise realizations,

[TABLE]

where frequency-frequency $\mathbf{N}$ is the noise covariance matrix and $\mathbf{\hat{d}}$ is the noiseless sky signal i.e., $\mathbf{A}(p)\,s(p)$. In Eq. B.3, the dependence on the spectral parameters is confined to the projection operator, $\mathbf{P}(p)$,

[TABLE]

$\langle{\cal L}_{spec}\rangle$ can be maximized very efficiently numerically, given that the number of unknown spectral parameters is usually limited and that one can capitalize on the analytical derivatives of the likelihood. xForecast, in the same manner as CMB4cast, based on Errard et al. (2011), uses a semi-analytical expression for the covariance of the error bars on spectral indices, $\mathbf{\Sigma}(\beta)$. This gives a computationally efficient way of estimating the statistical foregrounds residuals.

• •

the “inversion” of Eq. B.1 with the estimated $\mathbf{A}$, in order to disentangle sky components and obtain estimates of the sky signals $\tilde{s}$,

[TABLE]

From Eq. B.6, it should be evident that the noise variance, $\sigma_{\rm CMB}$, associated with the recovered CMB map is given by

[TABLE]

The noiseless foreground residuals are then given by

[TABLE]

Eq. B.8 can be rewritten and specialized for the CMB component residual,

[TABLE]

which does not contain CMB signal. $\mathbf{F}$ is a foreground matrix, and the $k$th column defines the total foreground contribution to the $k$th frequency channel. Stompor et al. (2016) perform a Taylor expansion of the residuals with respect to the scaling parameters around the maximum-likelihood values, $\tilde{\beta}$,

[TABLE]

Evaluating $\langle\mathbf{r}^{\rm cmb\ \dagger}_{\ell m}\mathbf{r}^{\rm cmb}_{\ell^{\prime}m^{\prime}}\rangle_{m,m^{\prime}}$ in harmonic space leads to an analytical expression for residuals power spectrum which can be found in Stompor et al. (2016). Schematically,

[TABLE]

The first term corresponds to the statistical residuals, generated by the finite error bar on the spectral indices, $\mathbf{\Sigma}$. Note that this term is the only one used in the CMB4cast framework (Errard et al. 2011, 2016). The second term corresponds to the systematic residuals, sourced by the mismatch between the fitted mixing matrix, $\mathbf{A}(\tilde{\beta})$, and the true mixing matrix used to generate the sky simulations.

• •

Optimization of a CMB+noise-averaged cosmological likelihood. Stompor et al. (2016) start from a standard Gaussian likelihood, which accounts only for noise, CMB signal and statistical foreground residuals in the recovered CMB map, $\tilde{s}(p)$. Assuming that they all are Gaussian with the total covariance given by $\mathbf{C}$, one can write

[TABLE]

where $\mathbf{a}$ is a harmonic representation of the map obtained after the component separation procedure. In addition to the CMB signal, it can include the noise as well as the statistical and systematic residuals. Similarly to Errard et al. (2011), this latter contributionis ignored in the assumed data covariance matrix, $\mathbf{C}$. The cosmological likelihood averaged over the instrumental noise and CMB signal realizations is given by (Stompor et al. (2016)),

[TABLE]

where $\mathbf{E}\equiv\langle\mathbf{a}\mathbf{a}^{t}\rangle$ is the correlation matrix of the data.

B.2 Application of xForecast to CORE simulations

In the exercise of cleaning the foregrounds in the CORE simulated sky maps, we simply consider two diffuse polarized astrophysical foregrounds: dust and synchrotron. Similarly to Eq. A.1, these are assumed to follow a grey-body and power-law spectra as,

[TABLE]

where the reference frequency $\nu_{\rm ref}=150$ GHz. We consider a modified grey-body emission law for the dust,

[TABLE]

These expressions are used to build the mixing matrix, $\mathbf{A}$, involved in the spectral likelihood, cf. Eq. B.3.

After the numerical optimization of the spectral likelihood, Eq. B.3, we evaluate the level of statistical and systematic residuals, Eq. B.13, and finally look at the cosmological likelihood, Eq. B.15, for three values of the tensor-to-scalar ratio, $r=10^{-2}$, $r=5\times 10^{-3}$ and $r=10^{-3}$. xForecast hence gives complementary results to the ones presented in the previous sections as it gives the ensemble-averaged level of foregrounds residuals as well as the averaged distribution of recovered $r$.

Results obtained in the case of the CORE baseline after foreground cleaning, but without any delensing, are shown in Figs. 31 and 32. It should be noted that the bias on $r$ is strongly related to both the complexity of the simulated foregrounds and the fidelity of the chosen parametrization, Eqs. B.16, B.17. In addition, the uncertainty on the measurement of $r$ depends primarily on the noise variance, the level of the statistical foreground residuals (sourced by the uncertainty on the estimation of the spectral parameters), and the lensing variance. The latter can potentially be improved on by delensing the cleaned $B$-modes map.