{'fraction_non_alphanumeric': 0.043140399878321946, 'fraction_numerical': 0.028760266585547967, 'mean_word_length': 4.1315453384418905, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Simple foreground cleaning algorithm for detecting primordial B-mode polarization of the cosmic microwave background

We reconsider the pixel-based,"template"polarized foreground removal method within the context of a next-generation, low-noise, low-resolution (0.5 degree FWHM) space-borne experiment measuring the cosmological B-mode polarization signal in the cosmic microwave background (CMB). This method was put forward by the Wilkinson Microwave Anisotropy Probe (WMAP) team and further studied by Efstathiou et al. We need at least 3 frequency channels: one is used for extracting the CMB signal, whereas the other two are used to estimate the spatial distribution of the polarized dust and synchrotron emission. No external template maps are used. We extract the tensor-to-scalar ratio (r) from simulated sky maps consisting of CMB, noise (2 micro K arcmin), and a foreground model, and find that, even for the simplest 3-frequency configuration with 60, 100, and 240 GHz, the residual bias in r is as small as Delta r~0.002. This bias is dominated by the residual synchrotron emission due to spatial variations of the synchrotron spectral index. With an extended mask with fsky=0.5, the bias is reduced further down to<0.001.


INTRODUCTION
Why study the B-mode polarization of the cosmic microwave background (CMB)? Detection of the primordial gravitational waves generated during inflation would give us a direct insight into the physical condition of the universe when the energy scale was close to the grand unification scale, ∼ 10 16 GeV (see Liddle & Lyth 2009, for a recent review and references therein). While a direct detection of the primordial gravitational waves using, e.g., laser interferometers, seems not possible with the present-day technology, an indirect detection using the B-mode polarization of the CMB (Seljak & Zaldarriaga 1997;Kamionkowski et al. 1997) may be possible in the near future (most optimistically, within a few years), provided that the energy scale of inflation at which the observed gravitational waves were generated was indeed as high as the grand unification scale.
We often characterize the amplitude of gravitational waves (also known as tensor perturbations) using the socalled "tensor-to-scalar ratio," which is conventionally defined as where h + k and h × k are the Fourier transform of the amplitudes of two linear polarization states of gravitational waves, and R k is the primordial curvature perturbation, which is a scalar perturbation (hence the name, "tensor-to-scalar ratio"). It is R k that seeded the ob-served structure in the universe, as well as the dominant component of the observed CMB temperature anisotropy (see Weinberg 2008, for a recent review and references therein).
The dominant, scalar part of the temperature anisotropy generates radial and tangential polarization patterns around hot and cold spots (Coulson et al. 1994). This is called the E-mode polarization, and has been detected with high statistical significance (Brown et al. 2009;Chiang et al. 2010;Larson et al. 2010;Komatsu et al. 2010;QUIET 2010). However, the B-mode polarization, which cannot be generated by the scalar perturbations but can be generated by the tensor perturbations, has not been found yet. The current 95% upper limit on the tensor-to-scalar ratio is r < 0.24, which mainly comes from the upper limit on the tensor contribution to the temperature anisotropy on large angular scales (Komatsu et al. 2010).
Given the upper limit on r, one can calculate the expected level of the B-mode power spectrum (see Figure 1). For r = 0.24, the B-mode power spectrum is smaller than the E-mode power spectrum by a factor of 10 at the first bump (created by electrons at z 10). At the second bump (created by electrons at z ≃ 1090), the B-mode power spectrum is smaller than the E-mode power spectrum by a factor of 50. It is the smallness of the B-mode power spectrum that makes the detection of this signal challenging.
There are three sources of noise for B-mode detection: (1) Detector noise; (2) Galactic foreground emission; and (3) Gravitational lensing. In this paper, we shall focus on the Galactic foreground. We use a mapbased method for reducing the Galactic foreground, and study how the residual foreground limits a measurement of the primordial B-mode polarization. The foreground reduction technique we use is motivated by the "template  (Larson et al. 2010), the QUaD final data (Brown et al. 2009), the BICEP two-year data (Chiang et al. 2010), and the QUIET 43 GHz data (QUIET 2010), respectively. The upper solid line shows the scalar E-mode power spectrum of the WMAP seven-year best-fit model. The dashed lines show the primordial B-mode power spectra with the tensor-to-scalar ratio of r = 0.24, which corresponds to the current 95% upper limit (Komatsu et al. 2010), as well as of r = 0.03 and 0.003. These lines are linearly proportional to r. The dotted line shows the secondary B-mode power spectrum expected to be generated by the weak gravitational lensing effect converting E modes to B modes (Zaldarriaga & Seljak 1998). This line is fixed (by the WMAP seven-year best-fit model) and acts as noise for the primordial B-mode detection. The lensing contribution becomes comparable to the primordial bump at l = 10 and 100 for r = 0.003 and 0.03, respectively. cleaning method" used by the WMAP team (Page et al. 2007;Gold et al. 2009Gold et al. , 2010. This method was further investigated by Efstathiou et al. (2009) in the context of the Planck mission. We shall study this technique in the context of a next-generation, low-noise, low-resolution (0.5 degree FWHM) space-borne experiment.
There is a large body of literature on the issue of polarized foreground cleaning for the B-mode detection. Our method is one specific (and relatively simpler) example. For the other methods in the literature, see review articles Fraisse et al. 2008) and references therein.
This paper is organized as follows. In Section 2, we show how the detector noise and the lensing noise influence the statistical errors on r. In Section 3, we describe our method for estimating r in the presence of the Galactic foreground and the dominant scalar E-mode polarization. In Section 4, we describe our simulation including CMB, detector noise, and foreground. In Section 5, we present the main results of this paper. We conclude in Section 6.

DETECTOR NOISE AND LENSING NOISE
Before we study the effect of the foreground, we show how the detector noise and the lensing noise influence our ability to detect r. The detector noise enters into the likelihood of r via the noise power spectrum, N BB l .
Assuming white noise, we write the noise power spectrum as Fig. 2.-B-mode polarization signal and noise power spectra. The dashed lines show the primordial B-mode power spectra with the tensor-to-scalar ratio of r = 0.01 and 0.001, while the dotted line shows the secondary B-mode power spectrum from gravitational lensing. We also show the noise power spectra (Equation (2) where w −1/2 p is the noise in Stokes parameters Q or U per pixel whose solid angle, Ω pix , gives Ω pix = 1 arcmin. This quantity is useful because one can compare various experiments on the same scale.
Current and future experiments use many (of order 10 3 − 10 4 ) detectors to reduce the noise equivalent temperature (NET) down to a few µK arcmin level. Is this sufficient for detecting primordial B modes? For comparison, the expected sensitivity of Planck combining 70, 100, and 143 GHz is w −1/2 p = 63 µK arcmin (see, e.g., Appendix A of Zaldarriaga et al. 2008;Planck Blue Book 2005).
In Figure 2, we compare the noise power spectra for w −1/2 p = 2 and 10 µK arcmin to the primordial and lensing B modes. For r = 10 −3 and the 10 µK arcmin noise, only a few modes (l = 2, 3, and 4) are above noise. For the 2 µK arcmin noise, the noise power spectrum is below the lensing B-mode power spectrum, and thus noise is no longer the limiting factor (unless we "de-lens" maps and remove the lensing noise). How would this influence our ability to detect r?
To see this, let us calculate the likelihood of r for a given noise level. For simplicity, we assume that we cover the full sky and the noise per pixel is homogeneous. 3 Then, one can write down the probability distribution function of the measured B-mode power spectrum,Ĉ BB l , for a given value of r as (e.g., Equation (  -Effect of detector noise and gravitational lensing on the likelihood for r (foreground is not included). The input value of r is r input = 10 −3 . The values on the horizontal axis have been multiplied by 10 3 . In the top left, top right, bottom left, and bottom right panels, we sum the multipoles up to lmax = 2, 5, 10, and 100, respectively. For the detector noise level, we use 2 µK arcmin for the solid and dashed lines, and 10 µK arcmin for the dotted lines. The solid and dotted lines include the gravitational lensing contribution to the total noise, while the dashed lines do not. Even if we set the detector noise to be zero, the solid lines do not change very much: for r = 10 −3 , the gravitational lensing effect prevents us from measuring r accurately beyond l ∼ 10. Note that a single multipole, l = 2, is enough for us to detect r = 10 −3 if the detector noise is smaller than 10 µK arcmin.
where c GW l is the primordial B-mode power spectrum from gravitational waves with r = 1, and c L l is the secondary B mode from gravitational lensing. We then use Bayes' theorem to calculate the likelihood for r as L(r|Ĉ BB (4) Figure 3 shows the likelihood of r for the input value of r input = 10 −3 and l max = 2, 5, 10, and 100. One useful number to keep in mind is that a single multipole, l = 2, is sufficient for detecting r = 10 −3 , if the noise is smaller than 10 µK arcmin. However, the precision on r does not improve beyond l = 5. This is apparent also in Figure 2: the noise power spectrum exceeds the signal at l ≥ 5.
We can improve the precision further if we lower the noise level to, say, 2 µK arcmin. Even so, the gravitational lensing prevents us from improving on the precision beyond l ∼ 10 if r = 10 −3 . (If there were no lensing in the universe, we would be able to continue to improve on the precision, as indicated by the dashed lines.) In fact, 2 µK arcmin is essentially the same as zero detector noise, as the lensing term dominates the error budget. Again, this is apparent in Figure 2.
Of course, these results are overly optimistic, as the error would be dominated by the foreground rather than by the detector noise. Nevertheless, it is still useful to know what would be possible when we ignore the foreground.
To quantify the precision on r, it is convenient to use the variance, σ 2 r , given by the second moment of the like- Here, σr is the square-root of the second-order moment of the likelihood function given by Equation (5). (The foreground is not included. The full sky coverage is assumed.) From the top to bottom lines, we show r input = 0.001 with and without the lensing noise, and r input = 0.01 with and without the lensing noise, respectively. For the instrumental noise level, we have used 2 µK arcmin. Note that σr/r input ∼ 1 at lmax = 3 does not mean that we do not detect r; on the contrary, we detect r with high significance even at lmax = 2. Rather, it just means that the likelihood for r is highly non-Gaussian and has a long tail toward large values of r (see the top left panel of Figure 3). In other words, we detect r with high statistical significance, but the value of r is not determined very well. lihood: Here, we have assumed that the likelihood is normalized such that ∞ 0 drL(r) = 1. One should be careful when interpreting this quantity. For l max = 2, σ r would be greater than the input value, r input = 10 −3 ; however, this does not mean that we cannot detect r. This just means that the distribution is highly non-Gaussian and has a long tail toward large values of r (see the top left panel of Figure 3). For large values of l max , e.g., l max 10, the distribution of r becomes approximately a Gaussian, and thus the value of σ r may be interpreted as the size of the usual 1σ error bar. Figure 4 shows the fractional error, σ r /r input , on the determination of the value of r as a function of l max . First, as one may expect from Figure 3, the fractional error for r input = 10 −3 saturates at l max ∼ 10 and does not improve further due to the lensing noise. For this case, while we can detect r with high statistical significance, we can determine the actual value of r to only ∼ 20%. For r input = 10 −2 , we can determine the value of r to ∼ 4% at l max = 200 (beyond which the fractional error no longer improves due to the lensing noise).
This study gives us an estimate of statistical errors on the measured values of r. On the other hand, the Galactic foreground gives us systematic errors (and bias). Now we shall turn to the foreground issue, which is the main subject of this paper.

Motivation
The basic idea behind our methodology is simple: we have (at least) 3 polarized components on the sky that we know and have been detected by the WMAP: CMB, synchrotron emission, and thermal dust emission. As the synchrotron dominates at lower frequencies and the dust at higher frequencies, we use one map at a low frequency and another map at a high frequency as the foreground "templates." We put the quotation marks here because these maps also contain the CMB. No external template maps are used in our method.
The WMAP team has applied this method for modeling the synchrotron: they used the lowest frequency (K-band, 23 GHz) map as a template, fitted it to the higher frequency maps (Ka, Q, V, and W bands), and subtracted from those maps. One can write this operation as where Q ′ and U ′ are the template-cleaned Stokes Q and U maps, respectively, and α S is the best-fit synchrotron coefficient for a given frequency ν. The denominator accounts for the fact that the K-band map also contains the CMB signal. However, the WMAP team had to rely on an external map for modeling the dust emission, as the highest frequency, the W band (94 GHz), was not high enough for being a good template of the polarized dust emission. This issue would probably be resolved by the Planck satellite, which has higher frequency channels such as 217 and 353 GHz. Efstathiou et al. (2009) have studied this by using a simulated Planck 217 GHz or 353 GHz map as a template for dust, and a simulated 30 GHz map as a template for synchrotron. They find that this simple method removes the foreground efficiently, bringing the bias in r down to a few times 10 −3 , which is much smaller than the expected statistical uncertainty on r from Planck, σ r = O(10 −2 ).
The goal of this paper is to put this method in the context of a next-generation, low noise (2 µK arcmin) polarization satellite experiment, and see if this method yields a promising result for measuring r ∼ 10 −3 (which is easy to detect in the absence of foreground, as we just saw in Section 2).

"Template" cleaning method
Our methodology is similar to that given in Section 4.2 of Efstathiou et al. (2009).
The main parameter that we wish to extract from data is the tensor-to-scalar ratio, r. (We do not vary the tensor tilt, n t .) The foreground coefficients, α, are nuisance parameters that we wish to marginalize over. The foreground coefficients may be spatially varying.
Another nuisance parameter (for detecting B modes) is the amplitude of the scalar E-mode power spectrum, which is by far the dominant source of CMB polarization. The signal power spectra are thus given as where c l denotes the power spectra with s = 1 and r = 1. The fiducial value of s is s = 1. We shall maximize the following likelihood function for estimating r, s, and α i : where is a template-cleaned map. This is a generalization of Equation (6) for a multi-component case. In this paper, i takes on "S" and "D" for synchrotron and dust, respectively, unless noted otherwise. For definiteness, we shall choose: These choices are somewhat arbitrary, but our preliminary optimization study indicates that this is a good configuration for achieving a smaller bias in r. A fuller optimization study, including more frequency channels, would require a more detailed specification of a given experiment (e.g., how many detectors one can fit in a given focal place; how low the detector noise can be as a function of frequencies), which is beyond the scope of this paper, but will be presented elsewhere. The covariance matrix in pixel space, C, for Stokes Q and U maps is given as (11) where c is the signal covariance matrix calculated from the theoretical power spectra, c l , (see Appendix A) and the noise matrices, N 1 and N 2 , are a noise covariance of a smoothed map (which is not diagonal) before the template cleaning is applied, and a small artificial diagonal noise matrix for a matrix regularization, respectively (see Section 4.1 for details).
For simplicity and clarity, we have ignored noise in template maps. For, if we assume that all three channels are similar in detector noise level, it is a good approximation, as α D ∼ 0.08 and α S ∼ 0.25, and the fractional contribution of the template noise to the covariance matrix is given by α 2 i , i.e., 6% effect in the derived error bars. Note that this is equivalent to ignoring P in Section 4.2 of Efstathiou et al. (2009).

CMB and detector noise
For CMB, we first generate the scalar and tensor polarization power spectra using the CAMB code (Lewis et al. 2000) with and without lensing contributions. We then generate Stokes Q and U maps at the Healpix resolution of N side = 128. The signal map has been smoothed with a 30 ′ beam (FWHM), representing a low-angularresolution CMB polarization satellite experiment targeting the primordial B modes.
To this smoothed signal map, we add random Gaussian noise given by σ 0 / N obs (n) per pixel in the direction of n. Here, σ 0 is related to noise w −1/2 p as σ 0 = π 10800 where N pix = 12(128) 2 = 196608 is the total number of pixels at N side = 128, and N obs is the number of observations per pixel. We adopt N obs from the "EPIC low-cost" (EPIC-LC) design (Bock et al. 2008). The noise is highest on the ecliptic plane and lowest on the ecliptic poles, similar to the N obs pattern of the WMAP. Note that the absolute value of N obs will cancel out in σ 0 / N obs (n) if we use the above formula: only the spatial distribution is taken from N obs , and the overall noise level is set by the assumed value of w −1/2 p . We shall use w −1/2 p = 2 µK arcmin for the rest of this paper. For this low noise configuration, the results are not sensitive to the details of the N obs pattern.
As we described at the end of Section 3, noise in template maps (at ν S = 60 GHz and ν D = 240 GHz) makes only a small contribution to the final covariance matrix. Therefore, for simplicity we add noise only to our CMB channel at 100 GHz. 4 We then apply an additional Gaussian smoothing to this signal-plus-noise map with 9.16 degrees (FWHM), which is 2.5 times the pixel size at N side = 16, and resample the smoothed map to N side = 16. Finally, as the smoothed map at N side = 16 is dominated by the scalar E-mode signal at all angular scales supported by the map resolution, the covariance matrix of this map is singular. In order to regularize the covariance matrix, we add an artificial, homogeneous white noise of 0.2 µK arcmin such that the map becomes noise dominated at the Nyquist frequency, l max = 3N nside − 1 = 47. The polarized synchrotron and dust emission are mod- Fig. 6.-B-mode polarization signal and foreground power spectra. The dashed lines show the primordial B-mode power spectra with the tensor-to-scalar ratio of r = 0.01 and 0.001, while the dotted line shows the secondary B-mode power spectrum from gravitational lensing. We also show the B-mode power spectra of the dust (dot-dashed line) and synchrotron (triple-dot-dashed line) emission at 100 GHz outside the WMAP P06 mask (N side = 128). (The measured foreground power spectra have been divided by f sky = 0.733 to approximately correct for the mask, as well as by the pixel window function at N side = 128.) The total (synch+dust) power spectrum is measured from the total map, which is slightly larger than the sum of the synchrotron and dust power spectra, as these two foreground components are spatially correlated. Note that the original PSM v1.6.2 dust map has the average polarization fraction of 1.5%, but we have multiplied the dust map by a factor of 3 to approximate a more recent dust template map adopted by the Planck collaboration.  .
Here,Q PSM andŨ PSM are the PSM Stokes parameters in units of antenna temperature, and g(ν) ≡ (e x − 1) 2 /(x 2 e x ) where x = hν/k B T CMB = ν/56.780 GHz converts the antenna temperature to thermodynamic temperature. (Q and U are in units of thermodynamic temperature.) For synchrotron, the position-dependent spectral index, β S (n), is calculated from the Haslam 408 MHz map (Haslam et al. 1981) and the three-year WMAP temperature map at 23 GHz (Page et al. 2007). The template maps at 30 GHz are taken from Miville-Deschênes et al. (2008).
For dust, the position-dependent spectral index, β D (n), as well as the unpolarized intensity map are taken from Model 8 of Finkbeiner et al. (1999). The polarization angles of dust approximately follow those of the synchrotron maps. The original PSM dust map has the average polarization fraction of 1.5% over the full sky, but we will multiply this map by a factor of 3 to approximate a more recent dust map used by the Planck collaboration.
Top panels of Figure 5 show the amplitude (P = Q 2 + U 2 ) of polarization intensity of synchrotron and dust at 100 GHz, while the bottom panels show the spectral indices, β S and β D . After adding the above foreground maps (smoothed with a 9.16-degree beam at N side = 128 and degraded to N side=16 ) to the CMB-plusnoise map, we mask the simulated sky by the WMAP P06 mask (f sky = 73%) (Page et al. 2007).
The norm of the pixel vector, [Q,U ], is 2259 × 2, where 2259 is the number of pixels outside the P06 mask. In order to mask the covariance matrix, we use the technique described in Appendix D of Page et al. (2007): we compute an inverse of 6144 × 6144 matrix and reduce it to 4518 × 4518 matrix using Equation (D7) of Page et al. (2007). (Note that there is a typo in this equation: D should be replaced by D −1 .) In Figure 6, we show the B-mode power spectra measured from the PSM (N side = 128) at 100 GHz outside the P06 mask. The total foreground power spectrum has l(l + 1)C BB l /(2π) ≈ 10 −1 µK 2 at l 10, which is 250 and 2500 times larger than the primordial B-mode spectra with r = 0.01 and 0.001, respectively. The problem seems formidable; however, as we show below, the simple cleaning method can reduce the foreground-induced bias in r to ∆r ≈ 0.002(< 0.001) with the P06(extended) mask.

Fixing the scalar E-mode amplitude
Before we use our full likelihood function given by Equation (9), let us first try a simpler version and show that it actually fails.
For the moment (only within this subsection), we fix the amplitude of the scalar E modes, i.e., s = 1, and consider cleaning dust using a map at 240 GHz. (Synchrotron will not be discussed in this subsection.) Our model is thus (16) As we described at the end of Section 3, we ignore noise at 240 GHz. We then fit the 240 GHz map to the 100 GHz map: Minimizing with respect to α D gives the following least-square solution: As the polarization signal is dominated by scalar E modes, we can set r = 0 when computing the covariance matrix C in this equation. (In practice, we used r input .) Finally, we maximize the likelihood given in Equation (9) with respect to r, with s = 1 and α D given by the above least-square solution.
The left panel of Figure 7 shows the values of r and α D obtained from many random realizations of noise and CMB skies. (The input tensor-to-scalar ratio is r input = 0.003.) There is a clear correlation between r and α D , indicating a failure of this algorithm. This correlation is caused by a chance correlation between foreground and the dominant scalar E modes (Chiang et al. 2008;Efstathiou et al. 2009). The correlation disappears -Correlation between the tensor-to-scalar ratio, r, and the amplitude of dust, α D . Left: The amplitude of the scalar E modes is held fixed at s = 1. Right: The amplitude of the scalar E modes is treated as a nuisance parameter and marginalized over. The input tensor-to-scalar ratio is r input = 0.003 and, for this figure only, the original PSM dust map (with an average polarization of ∼ 1.5%) is used. when we set C EE l = 0. This result motivates our treating the amplitude of scalar modes as a nuisance parameter.
The right panel of Figure 7 shows the results when s is treated as a nuisance parameter and marginalized over. For this, we have maximized the likelihood given by Equation (9) by varying r, s, and α D simultaneously. The correlation between r and α D has disappeared.
How well was dust cleaned? We have repeated this onecomponent foreground cleaning test for various values of r input from 0.001 to 0.1. The results are shown in Table 1: in all cases, the method recovers r successfully.

Cleaning synchrotron in multi-region
We are now ready to include synchrotron. Our model is  (240) (21) It turns out cleaning synchrotron is more challenging than cleaning dust, as the spatial distribution of synchrotron tends to be more extended above the Galactic plane than that of dust (see the top panels of Figure 5). We start by adding a mock synchrotron model (MSM) map to the PSM dust map. The MSM map has the same synchrotron polarization intensity across the sky as PSM at 30 GHz, but has a spatially invariant spectral index of β = −3.0 (where β = −3.0 is the average of spatially varying spectral index of PSM). With MSM plus PSM dust, r is recoverd successfully; mean(r) = 0.0012 and 0.0031 for r input = 0.001 and 0.003. (see the second and third columns of Table 2) Even more problematic is the spatial variation of the synchrotron spectral index (see the bottom left panel of Figure 5), which causes a mismatch between a template map at 60 GHz and the actual synchrotron distribution at 100 GHz. When we use a single synchrotron coefficient, α S , for the whole sky for the PSM model even without dust (= synchrotron only), we find a bias in r of order ∆r ≈ 0.002: mean(r) = 0.0028 and 0.0120 for r input = 0.001 and 0.01, respectively (see the fourth and fifth columns of Table 2).
One way to mitigate this issue would be to extend the Galactic mask (Efstathiou et al. 2009). In addition, one may give up using a single synchrotron amplitude for the whole sky, and use multiple amplitudes depending on the locations on the sky. 5 In this paper, we (Method I) continue to use the P06 mask, but divide the sky using the Healpix map with N side = 2, as shown in Figure 8a and (Method II) extend the mask to f sky = 50% and divide the sky using the Healpix map with  N side = 1, as shown in Figure 8b.
We give the details of our definition of the extended mask in Appendix B. In short, we choose the threshold polarization intensity values at 60 and 240 GHz above which the pixels are masked, such that we retain 50% of the sky.
While we would probably do a better job at cleaning synchrotron if we divide the sky according to our knowledge of the polarized synchrotron measured by WMAP, for this paper we prefer to explore a simpler algorithm and see how well we can recover r.
Each region I will be cleaned as (c.f., Equation (10)) (22) Note that we still use a single amplitude for dust on the whole sky. Similarly, the covariance matrix is given by (c.f., Equation (11)) where C IJ , c IJ and N IJ denote a block of matrices for pixels within regions I and J. The free parameters in the maximization are r, s, α D , and α i S (i=1. . . 12 (48)). In principle we wish to maximize the full likelihood function with respect to these parameters; however, in practice, this process is too time consuming to do brute-force, as varying each of these 15(51) parameters requires re-inverting a 4518 × 4518 matrix. Therefore, we make one approximation: we fix α D and α S in the covariance matrix (Equation (23)) at the bestfit values, α 0 D and α i,0 S . This is a good approximation as long as the noise term is sub-dominant compared to the dominant scalar E-mode signal, which is always the case for our low-noise configuration. 6 With this approximation, can be maximized with respect to r, s and α i where i runs from 1 (dust) to 13(49) (synchrotron for 2 to 13(49)). We use the MINUIT package (James 1988) for the maximization.
In the fourth and fifth columns of Table 2, we show the recovered values of r for the synchrotron-only cases, in order to see if dividing the sky into 48 regions helps to 6 As we have shown in Section 5.1, the foreground amplitudes and the dominant scalar E modes are covariant. Therefore, in order to find the best-fit αs without running the full likelihood, we had to "cheat" and measure αs in maps that do not contain the CMB signal or noise. Of course, we cannot do this in real life and thus we will have to come up with an efficient numerical algorithm for maximizing the full likelihood without this approximation. We believe that this is doable, so this will not be a limiting factor for our method.  reduce the bias in r that we have just seen. We find that the bias has reduced, but not by much: ∆r = 0.0018 → 0.0014, 0.0019 → 0.0016, and 0.0020 → 0.0015 for r input = 0.001, 0.003, and 0.01, respectively. This is probably due to the division not being tailored to match the distribution of synchrotron emission. While we keep this simple division and do not pursue a more complex division in this paper, we shall come back to this issue in the future work.

5.3.
Recovering r Now, we recover r from the full dust-plus-synchrotron cases. In Figure 9, we show the distribution of r and s for all of 68 realizations that we have run with r input = 0.003 using Method I. In Table 3, we show the recovered values of r in the second and fourth columns. Comparing them to the input values, r input , in the first column, we conclude that our Method I recovers r with a foregroundinduced bias of ∆r ≈ 0.002, which is consistent with the bias we have just seen from the synchrotron-only cases.

-3
10 -2 10 -1  -Recovered values of r (mean(r)) and error bars (std(r)) as a function of r input . The green points with error bars show the recovered r from the dust-only results (Table 1); the red (cyan) points with error bars show the synchrotron-plusdust results using Method I(II) ( Table 3) A systematic bias of ∆r ≈ 0.002(0.0006) is seen for the synchrotron-plus-dust results using Method I(II), which can be described by r recovered = r input + (0.0018 ± 0.0004)((0.0006 ± 0.0004)) (red(cyan) dotted line). We do not detect an offset for the dust-only results: (0.0000 ± 0.0003) (green dashed line) With Method II we recover r with a much smaller bias of ∆r ≈ 0.0006. We visualize our results in Figure 10.
The bias in r is important, but the uncertainty in the recovered r is equally important. In Figure 4, we have shown the predicted fractional errors on the determination of r, σ r /r input , for idealistic full-sky, foregroundfree cases. How would they look when synchrotron and dust are included and cleaned with our method? In Figure 11, we show the same figure but with the predictions made for the simulated CMB-plus-noise maps at N side = 16 as described in Section 2, and scaled to the P06 mask. We also show σ r /r input where σ r is extracted from the simulations. First, when both the foreground and lensing noise are ignored, the simulation and the analytical prediction are in a good agreement (to within 10%) for r input = 0.001 (see the star symbol in Figure 11). When the foreground is included, however, the error increases. For r input = 0.01, the foreground cleaning increases the error by about 60%. We see larger discrepancies between the foreground-free predictions and the foreground-cleaned results for r input = 0.001: the foreground-cleaned error is a factor of two larger than the foreground-free prediction; thus, the increase in the error due to foreground cleaning can be substantial when r is as small as 10 −3 .
Further optimizations could be done, given the details of a given experiment, and we intend to explore this issue within the context of some specific experimental designs.  Figure 4, but the foreground-free predictions are made for simulated maps at N side = 16 (as described in Section 2) and divided by f sky where f sky = 0.735 for the P06 mask at N side = 16. The filled circles use the error bars, σr, derived from our simulations with synchrotron and dust foreground cleaning (denoted as "std(r)" in Table 3). The denominators, r input , are not the values recovered from the simulations, but the input values fixed at either 0.01 or 0.001. Each filled circle (foreground-cleaned result from simulations) should be compared to the corresponding line with the same color (foreground-free, analytical prediction). The star shows the foreground-free and lensing-free result from the simulation with r input = 0.001, which agrees with the analytical prediction to within 10%.
Another improvement can be made by using N side = 32 (or 64) map, so that the Nyquist frequency is close to (or beyond) the second bump of the B-mode spectrum and more information is used. Such analysis, however, would take 2 6 = 64 (4 6 = 4096) times more computation.

CONCLUSIONS
In this paper, we have studied the pixel-based foreground cleaning method within the context of a nextgeneration, low-noise CMB polarization satellite. This method was originally applied to polarized data by the WMAP team (Page et al. 2007;Gold et al. 2009Gold et al. , 2010, and further investigated by Efstathiou et al. (2009) in the context of Planck.
Despite the simplicity of the method (namely, we have maps at 3 different frequencies, two of which are used for removing the synchrotron and dust emission), we are able to recover the input tensor-to-scalar ratio with only a small bias, ∆r ≈ 0.002(< 0.001) for the P06(extended) mask, which is dominated by the residual synchrotron emission. Further improvements should be straightforward: one can tune the Galactic mask, and divide the synchrotron fitting regions according to the actual distribution of the synchrotron spectral index in the Galaxy (rather than using the regular division shown in Figure 8). One may also increase the number of frequencies for measuring the spatial distribution of the synchrotron spectral index, provided that we have enough space on the focal plane. These will be investigated in the context of specific experimental designs such as LiteBIRD 7 , and presented elsewhere.
Our study suggests that a detection of the primordial B-mode polarization at the level of r ≈ 10 −3 should be possible with carefully optimized mask and α regions. Note that our statistical error and systematic bias becomes comparable with f sky = 50% mask case. However, let us mention one important caveat in our analysis. While our knowledge of the distribution and properties of the polarized synchrotron is fairly secure thanks to the WMAP data, our knowledge of the polarized dust emission, especially the spatial variation of the dust spectral index, is still highly limited. Therefore, the estimated bias in r that we have presented in this paper cannot be too accurate. Fortunately, Planck will soon provide us with maps of the polarized dust emission with the unprecedented sensitivity; thus, we intend to revisit this issue once the Planck data become available.
We thank J. Dunkley for providing us with the Planck Sky Model maps v1.6.2, and T. Matsumura for providing us with the map of N obs . We acknowledge use of the HEALPix (Gorski et al. 2005), CAMB (Lewis et al. 2000), and MINUIT (James 1988) packages. This work was supported by MEXT KAKENHI 21111002 and 22111510.