Galactic foreground constraints on primordial $B$-mode detection for ground-based experiments

Contamination by polarized foregrounds is one of the biggest challenges for future polarized cosmic microwave background (CMB) surveys and the potential detection of primordial $B$-modes. Future experiments, such as Simons Observatory (SO) and CMB-S4, will aim at very deep observations in relatively small ($f_{\rm sky} \sim 0.1$) areas of the sky. In this work, we investigate the forecasted performance, as a function of the survey field location on the sky, for regions over the full sky, balancing between polarized foreground avoidance and foreground component separation modeling needs. To do this, we simulate observations by an SO-like experiment and measure the error bar on the detection of the tensor-to-scalar ratio, $\sigma(r)$, with a pipeline that includes a parametric component separation method, the Correlated Component Analysis, and the use of the Fisher information matrix. We forecast the performance over 192 survey areas covering the full sky and also for optimized low-foreground regions. We find that modeling the spectral energy distribution of foregrounds is the most important factor, and any mismatch will result in residuals and bias in the primordial $B$-modes. At these noise levels, $\sigma(r)$ is not especially sensitive to the level of foreground contamination, provided the survey targets the least-contaminated regions of the sky close to the Galactic poles.


INTRODUCTION
In the next decade, the potential detection of primordial polarized B-modes in the cosmic microwave background (CMB) radiation could be one of the most important milestones in modern cosmology. B-modes are a key prediction of inflation. They are predicted to arise due to a background of primordial gravitational waves created just after the Big Bang. The amplitude of these primordial gravitational waves could vary by a few orders of magnitude, from so small that is impossible to detect, to large enough to be detected in the immediate future (Knox & Song 2002). The best current limits constrain the tensor-to-scalar ratio to r 0.002 < 0.056 (BICEP2 Collaboration et al. 2018;Planck Collaboration et al. 2020a). In practice, the experimental detection of primordial B-modes is extremely challenging, and the cosmological community is investing significant effort in experiments specifically designed to meet this challenge, such as satellite experiments like LiteBIRD (Hazumi et al. 2019;Sugai et al. 2020), as well as ground-based telescopes, such as the Simons Observatory (SO; Lee et al. 2019) and CMB-S4 (Abazajian et al. 2019).
The polarized foregrounds are perhaps the biggest obstacle we need to overcome for such a detection because they block our otherwise unobstructed view of the CMB. The data analysis techniques used to model and mitigate contamination by foregrounds are termed component separation. We have methods to remove most of the contaminating signal, but small residuals left over by systematics and modeling mismatch are a source of bias for the polarization B-modes. These residuals could be comparable to the primordial B-mode signal we aim to measure. Several works in the recent literature have aimed at forecasting the performance of different satellite and ground-based experiments for measuring the primordial B-modes, with an emphasis on the foreground modeling (e.g. Betoule et al. 2009;Bonaldi & Ricciardi 2011;Katayama & Komatsu 2011;Armitage-Caplan et al. 2012;Errard & Stompor 2012;Remazeilles et al. 2016Remazeilles et al. , 2018bAlonso et al. 2017;Chluba et al. 2017a;Hervías-Caimapo et al. 2017; Thorne et al. 2019;Remazeilles et al. 2021).
Ground-based experiments, unlike the more expensive satellite experiments, do not cover the full sky. Instead, they must choose an area to target their observations arXiv:2105.06311v2 [astro-ph.CO] 12 Nov 2021 (e.g. Stevens et al. 2018;Abazajian et al. 2019) and balance between focusing on a small region to increase the map sensitivity and covering the largest possible region to reduce sample variance in the primordial B-mode signal (which peaks on large angular scales). Moreover, ideally we would like the complete absence of any foreground contamination when observing the CMB. This is not possible, so ground-based experiments often choose the areas where the foregrounds are the weakest, which in general correspond to areas near the south Galactic pole, observed from the southern hemisphere in Chile and Antarctica. No experiment so far has surveyed the north Galactic pole in its entirety, which may be just as good. However, for measuring primordial B-modes, avoiding the foregrounds is certainly not sufficient. Even component separation, which helps in the cleaning of foregrounds, is not sufficient. The selection of an optimal survey area for measuring the primordial B-modes (and r) needs to consider both factors: on one hand we wish the foregrounds to be as weak as possible (the foreground avoidance approach) while on the other, we require the foregrounds to have sufficient signal-to-noise ratio (S/N) so that component separation algorithms can measure their properties and thus clean them.
In this work, we investigate this issue. We analyze this balance and try to optimize the selection of the survey field in order to obtain the best measurement of the tensor-to-scalar ratio. We use SO as our example survey, but this paper should not be interpreted as an official forecast for this experiment, but rather a guide on the relative performance between different regions of the sky.
In Section 2, we describe the model and simulated observations we use, with a configuration based on the SO experiment. In Section 3, we describe our pipeline in detail. In Section 4, we present our results for the forecasted error on the detection of the tensor-to-scalar ratio, σ(r). In Section 5, we discuss a few specific issues, regarding the use of extra high-frequency channels to monitor the thermal dust emission. Finally, in Section 6, we draw our conclusions.

SIMULATIONS
We simulate observations by modeling the SO Small Aperture Telescopes (SATs) survey as our example experiment, which will target an f sky ∼ 0.1 area of the sky at degree resolution, with the intention of measuring the primordial BB angular power spectrum and constraining the value of r Lee et al. 2019). We refer to the area of this survey as the SO-SAT mask.
The simulated observations of the microwave Galactic foregrounds are created with the pysm3 code , smoothed with a Gaussian beam of appropriate width in order to model the SO-SAT instrument. The frequency channels, beam widths, and sensitivities used for the simulations are listed in Table 1. The maps are created at N side = 512. The model we use includes polarized thermal dust and synchrotron emission. Their spectral energy distribution (SED) is a modified black-Body (MBB) for thermal dust, given by in antenna temperature units, where the spectral index β dust and the dust temperature T dust are the free parameters, and h and k are the Planck and Boltzmann constant, respectively. For synchrotron, the SED is a power law with an index β syn close to 3, The modeled dust and synchrotron emission correspond to models d1 and s1, respectively, in Thorne et al. (2017). The synchrotron polarization template is based on the 9 yr WMAP maps (Bennett et al. 2013). The synchrotron SED is the power law from eq. 2 with a spatially variable index, which is taken from 'model 4' of Miville-Deschênes et al. (2008). The thermal dust polarization is constructed from the Planck 2015 data release (Planck Collaboration et al. 2016a), using the model fitted with the commander component separation code. The polarization template at the anchor frequency of 353 GHz, as well as the spatially variable maps of β dust and T dust are used to build this dust model.
The CMB realizations are created with a fiducial power spectra created with camb (Howlett et al. 2012) with the full lensing BB signal, an amplitude of primordial scalar perturbations A s = 2×10 −9 , an index for the primordial scalar perturbation power spectra n s = 0.965 and r = 0. Otherwise, we use a standard ΛCDM cosmology with the Hubble constant H 0 = 67.5 km/s/Mpc, baryonic matter density parameter Ω b h 2 = 0.022, cold dark matter density parameter Ω c h 2 = 0.122, sum of neutrino masses m ν = 0.06 eV, and reionization optical depth τ = 0.06.
We simulate the instrument noise with the public SO noise generator 1 . This includes the use of atmospheric 1/f noise curves with their optimistic knee , listed in Table 1. We create two sets of noise realization with the baseline and goal levels. Further details are provided in Ade et al. (2019). We use only the auto N spectra curves, not the cross-spectra, so the noise is uncorrelated between frequency channels in our simulations. The white-noise level sensitivities, along with the 1/f parameters, are listed in Table 1.
We also consider the Cerro Chajnantor Atacama Telescope-prime (CCAT-prime) project in the Fred Young Submillimeter Telescope (FYST) (Choi et al. 2020), a new submillimeter 6 m telescope scheduled to start observations roughly at the same time as SO. One of the science objectives of the experiment is to contribute to CMB experiments with high-frequency observations and thermal dust monitoring. The synergy between the CCAT-prime project and SO/CMB-S4 will allow for better constraining of the thermal dust foreground that hinders the potential detection of CMB Bmodes.
We simulate the observations in the 350, 410, and 850 GHz frequency channels from the Prime-Cam instrument in the FYST, using the same sky model. We use the noise power spectra curves to generate 1/f noise simulations as described in Choi et al. (2020). The details of the 1/f parameters and white-noise levels are listed in Table 1.

METHODOLOGY
Our full pipeline consists of running the Correlated Component Analysis (CCA) estimation over a given region of the sky on the simulated observations and subsequently reconstructing the angular power spectra of the three components that are included: CMB, thermal dust, and synchrotron emission. We apply this pipeline to 200 Monte Carlo iterations to estimate the covariance matrix of the B-mode angular power spectrum of the CMB. Finally, we estimate a σ(r) error using the Fisher information matrix.

Component separation and estimation of the BB power spectrum
The component separation method used for our analysis is the CCA (Bonaldi et al. 2006;Ricciardi et al. 2010). This is a parametric foreground fitting algorithm that exploits a generalized least-squares method applied to second-order statistics. In its harmonic domain version, which we use in this work, it estimates both the frequency behavior of all the foregrounds, expressed as a function of a few parameters, and the auto-and crosspower spectra of the components. A parametric method is particularly useful when investigating the correlation between the goodness of foreground fitting and the intensity of the foregrounds, which is the goal of this analysis.
As with many other component separation techniques, CCA is a linear method. It assumes that the observed signal is a linear combination of the true signal from the different components in the sky plus some noise. By making a few other assumptions (that spectral and spatial features on the components are not correlated, and that the instrumental beam is constant within each passband), the signal can be written for each line of sight r as x(r) = [B * Hs](r) + n(r), where x and n are vectors containing the observed signal and noise, respectively, B is a matrix containing the beam per frequency channel, s is the vector containing the true signal of each component, and the symbol * denotes convolution. H is the mixing matrix that contains the frequency spectra (SED) of all the components (three in our case: CMB, synchrotron, and thermal dust). If we transform to harmonic space, eq. 3 becomesx =BHs +ñ, wherex,s, andñ are the harmonic transforms of x, s and n, respectively, andB is the harmonic transform of B.
The relation between the covariances is given bỹ where the dagger superscript denotes the adjoint matrix, and the angular power spectra C for the observations x, for the true signal s and noise n. C x and C n are estimated from the data and noise properties,B is also known, while C s and H are the unknowns.
To reduce the number of unknowns, the mixing matrix is expressed as a function of a few parameters H = H(p). We use the blackbody spectrum for the CMB (with no free parameters); eq. 2 for synchrotron, with the spectral index β syn as a free parameter; and eq. 1 for thermal dust, where the spectral index β dust and the temperature T dust are, in general, both free parameters.
Finally, the data and source spectra are binned in multipoles to further reduce the number of unknowns. Then, in harmonic space, our model is the following equation: where d V are the noise-bias-subtracted auto-and crosspower spectra of the observations and arranged into a vector over the multipole binsˆ . H kB is a blockdiagonal matrix that contains one block per multipole bin. Each block is a matrix given by where ⊗ is the Kronecker product (hence the k subindex). This matrix, just like the mixing matrix, depends on the parameters p. c V is the power spectra of the different source components in the model, again arranged into a vector over multipole bins, and V represents the residuals of the power spectra, between the true sources in the sky and the modeled source components. Ricciardi et al. (2010) shows that the mixing matrix parameters, p, and the binned power spectra of the components, c V , can be found by minimizing the functional where N B is a block-diagonal matrix, containing one block per multipole bin. Each block is the matrix N (ˆ ), which corresponds to the auto and cross terms (over the frequency channels) of the Gaussian covariance matrix for the noise power spectra.
For each value of the mixing matrix parameters p, the estimated source angular power spectrac V are obtained as a suitable linear mixture of the data power spectra, depending on the estimated mixing matrix H kB (p) and the noise covariance matrix N B , given bȳ The solution for both the mixing matrix and the source spectra yields the minimum value of eq. 7. The minimization in CCA is performed with a simulated annealing over the foreground spectral parameters.
We have modified the CCA method slightly as compared to the one presented in Ricciardi et al. (2010). The estimation is done over the polarization spectra EE and BB, as opposed to working over the spectrum of any random scalar field, as before. We also used a more general noise covariance N (ˆ ) than the white-noise-only approximation used in Ricciardi et al. (2010). The namaster software  has been used to estimate the angular power spectra of the data and the noise as well as the noise covariance matrix (whose estimation adapts the analytical estimate for a Gaussian covariance proposed in Efstathiou (2004); Couchot et al. (2017)). We set the cross-power spectra to zero because we do not simulate noise cross-correlation among the frequency channels.

Estimation of the r uncertainty
To translate the error on the BB power spectrum to an uncertainty on r, σ(r), we use the Fisher information matrix. This is given as where N b is the number of multipole bins, C is the covariance matrix over multipole bins of the observable, and ∂cˆ /∂p i is the binned derivative of the fiducial model spectra with respect to the parameter p i . The covariance matrix of the parameters is given by the inverse of F ij . Therefore, the uncertainty of parameter i (marginalized over the rest of the parameters) is given by the square root of the ii element of the inverse of the Fisher information matrix. The off-diagonal elements will measure the correlation between parameters i and j. In our case, the observable is the estimated CMB BB power spectrum C BB and the parameter we are interested in is the tensor to scalar ratio r. To correct for residual foreground contamination in the CMB BB power spectrum, we introduce a second parameter, A fore , which represents the residual foreground amplitude. We do not consider the foreground spectral parameters explicitly in the Fisher matrix, but they are included in the estimation of the Cˆ ,ˆ covariance via the Monte Carlo simulation. Therefore, our fiducial model for the BB power spectrum is the following: The covariance matrix of the observable C is easily calculated from 200 Monte Carlo iterations over the simulations. The derivative of the fiducial BB power spectra with respect to r and A fore is calculated numerically with respect to mean values r = 0 and A fore = 1.
To estimate C BB,foregrounds , we produced companion simulations with only the foregrounds, without CMB and noise. This dataset is mixed together with the same weights calculated by CCA, given in eq. 8. Finally, the power spectra C BB,foregrounds is given by the average over the 200 Monte Carlo simulations. This is labeled method A. We realize this procedure is not feasible with real observations, because pure foregrounds maps are not available. We also consider a second method, which is very pessimistic: Because thermal dust is the dominant foreground (as will be shown in Sec. 5.1), we produce 200 Gaussian realizations following a power-law angular power spectrum fitted in Planck Collaboration et al. (2020b) using a mask labeled Large Region (LR) 71, with f sky = 0.71. We extrapolate this Gaussian realization map at 353 GHz to other frequencies using a standard MBB (eq. 1) with the β dust and T dust maps from Planck Collaboration et al. (2016b). Note that because we are using the LR71 mask including the full Galactic emission, power will be grossly overestimated in the cleanest regions, but because we are using the Fisher information matrix, the derivative will cancel the absolute amplitude of the foreground residual spectrum. We mix the 200 power-law spectra Gaussian realizations with the same weights calculated by CCA, as described above. This is labeled method B. These two methods are the extremes in realism when accounting for foreground residuals in a full pipeline for the posterior of r, but in real observations, the data analysis will be closer to method A than B, because we could use foreground models like the ones from pysm, which is significantly more realistic than a Gaussian random field. We calculate all of our power spectra using D = ( + 1)C /2π instead of C .

RESULTS
To represent the SO-SAT instrument specifications, we consider a survey over a disk-shaped region covering f sky = 0.1, which is a disk with a radius of 36.8 • . These masks are apodized by smoothing with a 1 • FWHM Gaussian kernel. We define 192 of such regions distributed over the full sky, with centers matching the pixels in a healpix pixelization with N side = 4. In each of these regions, we run our pipeline to get the uncertainties on r through component separation, power spectrum estimation, and the Fisher matrix analysis.
As mentioned in Sec. 3, the thermal dust is modeled in the mixing matrix with eq. 1, with the spectral index β dust and the temperature T dust potentially both free parameters. However, because the SO-SAT has a limited range of observing frequencies (ν ≤ 280 GHz) where thermal dust dominates, constraining the two thermal dust spectral parameters (β dust and T dust in eq. 1) at the same time is very difficult. A perturbation of ±1 K in the temperature of a fiducial MBB with β dust = 1.53, T dust = 23 K, and an anchor frequency of 27 GHz, results in a change of ∼ 1% in the SED at 280 GHz. To be able to constrain T dust precisely, more frequency channels are required in the range where T dust affects the behavior of the MBB, ν > 300 GHz. We therefore find that for the fiducial six frequency channels SO-SAT configuration we need to fix the value of T dust to 19 K, a representative value of the dust temperature away from the Galactic plane region, and leave only the dust spectral index β dust as a free parameter.
As a first example, Fig. 1 shows the beam-and maskcorrected reconstructed binned power spectra of the three sources for the SO-SAT mask in the fiducial six channels with baseline noise configuration. The dashed lines show the fiducial power spectra for the three components in the model, calculated for the two foregrounds by running the power spectrum estimation in the considered region over the simulated model with only synchrotron or only thermal dust emission, respectively. The individual points represent the reconstructed angular power spectra for each source for each of the 200 Monte Carlo simulations.
We note the discrepancies between the model spectra and the estimated spectra at high multipoles, clearly visible in the synchrotron and to a lesser extent in the CMB. This arises from the fact that the channels that constrain synchrotron are mostly the 27 and 39 GHz channels, which have low spatial resolution. Therefore, at small angular scales, the synchrotron reconstruction is very noisy. The difficulty in accurately measuring the polarized foregrounds in this region, because they are weaker, has an impact on the scatter of the reconstructed spectra for the thermal dust and synchrotron. This large scatter on the foregrounds spectra also somewhat impacts the CMB reconstruction.
As another example, we consider in Fig. 2 the case for the masked region centered at Galactic coordinates  Fig. 3) with six SO-SAT channels in the baseline noise level. The left is the EE power spectrum, and the right is the BB power spectrum. The dashed blue, green, and red lines are the fiducial spectra for the CMB, synchrotron, and thermal dust emission, respectively. The simulated CMB has r = 0. The points with the same color scheme correspond to 200 Monte Carlo simulations of the reconstructed spectra. The reference frequency for the foregrounds is 145 GHz. Note that there is no significant bias in the estimation of the CMB BB spectrum (blue points in the right-hand-side plot). l = 0.0 • , b = 9.6 • , very close to the Galactic Center (GC). In this case, there is a sizeable bias seen in the CMB BB estimated spectrum, the product of the large impact of polarized foreground residual in the GC region. The scatter on the spectra reconstruction of the foreground components, however, is smaller compared with the SO-SAT mask case. The comparison between Figs. 1 and 2 nicely illustrates the kind of trade-off between foreground bias and foreground scatter that this work aims to investigate.

Results for σ(r)
Our estimation of σ(r) (marginalized over A fore for method A, using the pure foregrounds to construct the foreground residual spectrum) over the 192 f sky = 0.1 circular regions is shown in Fig. 3, for the SO-SAT fiducial baseline noise configuration. Each dot is located at the center of the region and its color represents the estimated σ(r) value according to the color scale. In the background, we plot the polarization intensity P = Q 2 + U 2 for the Planck 353 GHz frequency map (Planck Collaboration et al. 2020c). The light-green contours show the SO-SAT mask footprint. We also show the histogram of σ(r) over the masked regions in the lower-left corner.
In general, the σ(r) values are of the order of a few 10 −2 for regions in the Galactic plane, due to fore-ground contamination. The regions with the lowest σ(r) are close to the north and south Galactic poles. Most of them coincide with the location of the SO-SAT mask (shown as green contours in Fig. 3), with some few exceptions of regions with high decl., not observable from the SO site in the Atacama desert in Chile, above the SO-SAT mask contours at the right and left edges. The estimation over the SO-SAT mask yields σ(r) = 4.9 × 10 −3 (2.7 × 10 −3 before marginalization over A fore ). Table 2 lists the σ(r) constraints for both the SO-SAT mask and the best result (for method A) of the 192 masked regions considered in this work, as well as the result for method B for the same masks, for the SO-SAT fiducial baseline noise configuration.
The σ(r) estimated in the forecast done by SO (Ade et al. 2019) finds σ(r) = 1.9 × 10 −3 in the Fisher matrix configuration. Our result is ∼ 2.5 times that value. We can attribute this difference to several factors. First, the SO Fisher forecast is a different method from what we do in this work. Their approach is a direct likelihood Monte Carlo Markov Chain (MCMC) of the simulated observations' auto-and cross-spectra, modeling the foregrounds and the CMB at the same time, skipping component separation, an approach followed in, e.g. BICEP2/Keck Collaboration et al. (2015). The foregrounds are modeled as power laws in multipole space, Figure 2. Analogous to Fig. 1. Reconstructed power spectra for the high-foreground masked region centered at Galactic coordinates l = 0.0 • , b = 9.6 • , close to the GC. The reference frequency is 145 GHz. The simulated CMB has r = 0. Note the substantial bias due to foreground contamination in the estimation of the CMB BB spectrum (blue points in the right-hand-side plot).  as well as an MBB for thermal dust and a power law for synchrotron in frequency space. Having the former constraint as a stronger prior would lead to improved agreement between the simulated and (likelihood) modeled foreground components, and in turn to an improved estimate of the CMB signal. This is a constraint we do not impose in our component separation modeling (we only impose the SEDs in frequency space). Also, the fact that we perform component separation in a previous and separate step does amplify the noise to some level when inverting matrices at the point of reconstructing the power spectra of the components, which would explain some of the increased scatter of the D BB estimates, which translate into a larger σ(r) in the Fisher forecasts. Another difference between our modeling and the one done in Ade et al. (2019) is that the Fisher matrix of the D BB is calculated as the second derivatives of the likelihood around the fiducial parameters, while in our case this is calculated empirically over 200 Monte Carlo iterations. Finally, the effect of the A fore parameter does increase σ(r) in our approach when marginalizing over it, as it does to other component separation methods in Ade et al. (2019) that do marginalize over foreground residuals, such as xForecast, BFoRe and ILC. These methods see an increment in the error bar to σ(r) ∼ 3.4 × 10 −3 . In Fig. 4 (top), we show some examples of the correlation between r and A fore for method A. A fore is better constrained in regions with high foreground contamination. At the same time, r and A fore are anticorrelated in these same regions and that anticorrelation diminishes toward zero as we move away from the GC towards the foreground-free Galactic pole regions. Fig. 4 (top) shows this for four representative regions: one close to the north Galactic pole, one close to the south Galactic pole, one close to the GC, and the SO-SAT mask. This correlation allows us to quantify the bias introduced by the foreground contamination.
In Fig. 4 (bottom) we compare the σ(r) constraint in all 192 circular masked regions with and without marginalization over the A fore parameter for method A 2 . ∆σ(r) is the difference between the marginalized and unmarginalized σ(r) error bar. We can see that the foreground contaminated regions along the Galactic plane have a much larger ∆σ(r). Weak polarized foreground regions in the Galactic poles have ∆σ(r) of a few 10 −3 , while strong polarized foreground regions in the Galactic plane have ∆σ(r) of a few 10 −2 . In practice, this shows the effect of marginalizing over A fore and how it leads to a σ(r) parameter that is reflective of both scatter as well as bias due to strong foreground contamination.

Different instrumental configurations
Besides the fiducial SO-SAT instrumental configuration with baseline noise level, we consider two extra instrumental configurations: the fiducial SO-SAT channels but with the lower noise goal performance; and the fiducial SO-SAT channels with the goal noise performance supplemented with the three CCAT-prime channels at 350, 410, and 850 GHz, labeled the extended configuration. The results for these configurations are listed in Table 2.
The improvement introduced by the lower noise level is evident. The σ(r) constraint over the SO-SAT mask improves from σ(r) = 4.9 × 10 −3 to σ(r) = 3.3 × 10 −3 when comparing baseline and goal levels. Fig. 5 (top) shows the confidence contours and posterior probabilities for three representative masked regions for this instrumental configuration, analogous to Fig. 4 (top).
The resulting forecast for the extended frequency configuration on the SO-SAT mask is σ(r) = 3.1 × 10 −3 , which is slightly better than the constraint achieved with only the six fiducial SO-SAT channels. For this configuration, the spectral parameter T dust was still kept fixed to 19 K. In Sec. 5.1, we discuss why our σ(r) forecast is about the same (or improves a little) even though in this case we are including additional information from three extra frequency channels to constrain the thermal  . Top: constraints on the r and A fore parameters for the baseline noise level using method A. 1σ and 2σ contours, along with the marginalized posterior probability for each parameter, for four representative masked regions: one close to the Galactic north pole, one close to the Galactic south pole, one close to the GC, and the SO-SAT mask. Bottom: ∆σ(r) over all 192 masked regions in Galactic coordinates calculated with method A. ∆σ(r) is the difference between the σ(r) marginalized over A fore and the σ(r) without marginalization.
dust. Fig. 5 (bottom) shows the confidence contours and posterior probabilities for the same three representative regions for this extended instrumental configuration.

Adding extra high-frequency channels in CCA
We have verified that the most important contributor to the polarized foreground contamination in the fiducial six-SO-SAT-frequency channel configuration is the thermal dust. We have estimated the dust and synchrotron residuals separately by mixing simulations that only contain each corresponding foreground with the CCA-estimated weights. In these, the synchrotron residual is significantly smaller than the thermal dust one, so from this point forward, we focus our discussion on the effect of thermal dust. However, a similar analy- sis can be repeated for synchrotron and how its residual can be mitigated. The high-frequency channels from CCAT-prime are highly dust dominated, and as such, they can measure the thermal dust precisely. They also have the potential to introduce strong foreground contamination to the estimated CMB signal if we do not model the SED perfectly. Whether adding such channels gives better or worse results in terms of CMB recovery critically depends on the complexity of the foreground sky in the considered region and the adequacy of the adopted component separation model to represent it. In our case, the balance is not quite so favorable, as the simulated thermal dust has a spatially variable T dust and β dust , while our model adopts constant values within a given . Example of the contribution per frequency channel for the reconstructed CMBxCMB D BB when we run the extended configuration with goal noise levels in the SO-SAT mask. We show the term in eq. 8, which is the weight times the dV vector, which corresponds to the power spectra of the full mix minus the power spectra of a noise realization. We only plot the 9 auto-BB power spectra for the 9 frequency channels, we omit the other 36 cross-power spectra for clarity, but the CCA uses all 45 auto-and cross-power spectra in its calculation.
patch. We note that the small improvement in σ(r) when adding the CCAT-prime high-frequency channels is a consequence of dust residuals contaminating the primordial BB, which is due to this mismatch between the simulated model and component separation modeling.
First, we have checked that if we consider simpler simulations, e.g. a spatially constant MBB SED, the agreement between the simulation and the component separation modeling is much closer for thermal dust, therefore the thermal dust residuals are smaller and there is a noticeable improvement in σ(r) when using the extra information from the high-frequency CCAT-prime channels.
To illustrate the issue in our simulation with spatially variable spectral parameters, we can look at the weights that are used to mix the observed data power spectra d V , described in eq. 8. Fig. 6 shows the product of the weights and the data vector d V =C x (ˆ ) −C n (ˆ ) for all multipole binsˆ , where d V is the noise-biassubtracted full-mix signal. This is for the CMB D BB within the SO-SAT mask with the extended configuration and goal noise levels. For clarity, this figure shows only the contributions for the auto-power spectra for each of the 9 frequency channels, even though the CCA method uses all 45 auto-and cross-spectra between the frequency channels. This figure illustrates the relative contribution from each channel to reconstruct the CMB angular power spectrum D BB . For example, the 93, 145, and 225 GHz channels contribute the most to the CMB, which makes sense because these are the frequencies where the CMB emission is strongest. Conversely, the low and high frequencies of 27, 39, 350, 410, and 850 GHz have the least relative contribution, which also makes sense. Ideally the CCAT-prime channels weights would be small, but just enough to subtract the dust out of the the channels where the CMB is strong. For the thermal dust, the mismatch in the simulated SEDs and the component separation modeled SEDs in CCA will result in errors in the weights, with the effect of high-frequency channels contaminating the CMB reconstruction with thermal dust residuals. If the weights for the high-frequency channels are small, but they are multiplied by a very strong thermal dust power spectra, such as the one present in the 350, 410, or 850 GHz channels, they will introduce some dust residual in the CMB reconstruction, which will contaminate the B-modes in a significant way.
As another example, we can look at the BB thermal dust residuals directly. In Fig. 7, we plot the reconstructed CMB BB spectrum as triangles in the SO-SAT mask for the fiducial (blue) versus extended (red) configurations, both with goal noise levels. Shown in the plot is the mean over the 200 Monte Carlo simulations, while the error bar corresponds to the standard deviation of that. We also plot the BB thermal dust residuals (squares) for the same two configurations. These are calculated with the maps of only the simulated thermal dust, mixed with the CCA-calculated weights used for the full simulation for the corresponding case. This is done in the exact same manner as how we calculate C BB,foregrounds in eq. 10. Both configurations have roughly the same error bar in their estimation of the CMB, as well as roughly the same dust residual. The nine-channels extended configuration has a slightly smaller nondiagonal term in the Fisher matrix, which means that r and A fore are slightly less correlated, and this configuration has a slight edge over the 6-channel fiducial configuration. We do not see the dust residuals diminishing substantially when adding the extra information from three high-frequency channels, which suggests that our thermal dust component separation modeling is not detailed enough.
We know that parametric methods might not have enough complexity to model the dust residuals that will bias r (Hensley & Bull 2018 CMB reconstructed, fiducial 6 channels Thermal dust residual, fiducial 6 channels CMB reconstructed, extended 9 channels Thermal dust residual, extended 9 channels in the literature that increase the modeling complexity for the thermal dust SED by, for example, using an improved pixel-based foreground parameterization component separation, which allows for the partition of the sky into disjoint patches that can be modeled independently (Errard & Stompor 2019). In recent years the moment expansion has shown promise when performing thermal dust component separation (Chluba et al. 2017b;Mangilli et al. 2021;Remazeilles et al. 2021;Azzoni et al. 2021). In these methods, the dust spectral parameters are modeled as some value with a small perturbation, so we can account for a higher-order expansion when describing the spatial variability of the dust SED. This is also handy when describing deviations from a perfect MBB spectrum, as we should expect from the measured dust frequency decorrelation (Planck Collaboration et al. 2017;Clark et al. 2021). We can always increase the complexity of the modeling by adding more parameters to the model, but this will increase the degrees of freedom and ultimately dilute the statistical certainty of the estimates. This analysis stresses the importance of careful modeling in the component separation. Dust-monitor channels can be included for improving the dust modeling (e.g. fitting T dust ), but our analysis suggests that maybe they should not be used naively in the CMB reconstruction itself. We are not saying that using the extra high-frequency information is not useful but rather that our component separation modeling is not good enough in dealing with the residuals that arise by modeling as spatially constant something that is actually spatially variable. Extra care is warranted when modeling components and minimizing the mismatch between observations and models, especially if one wishes to include such high-frequency channels in the component separation analysis.

How Correlated Is σ(r) with the foreground Contamination?
Given our analysis, we now attempt to answer the question that we posed at the beginning of this work: Is the region of sky where the foregrounds are at a minimum also the region that yields the smallest possible σ(r)?
In Fig. 8, we plot the σ(r) forecast (for method A) versus the polarized intensity of the thermal dust for each of the 192 masked regions we consider in this work. The polarized thermal dust intensity is measured by calculating the BB angular power spectrum of the 2018 data release Planck 353 GHz frequency map (Planck Collaboration et al. 2020c) as prescribed in Planck Collaboration et al. (2020b, and fitting a power law to it in the multipole range = 40 − 599, given by D BB = A BB dust ( /80) α BB +2 . In the figure, we plot the A BB dust amplitude. In blue, we show the fiducial configuration with baseline noise levels, the red points are the 10 2 10 3 10 4 A BB dust [ K 2 CMB ] 10 2 (r) SO SAT 6 channels baseline SO SAT 6 channels goal SO SAT extended goal Figure 8. Correlation between the σ(r) forecast for method A and the amplitude of polarized thermal dust BB power spectrum over the 192 masked regions for the σ(r) value after marginalization over A fore . The blue points correspond to the fiducial SO-SAT channels with the baseline noise level, the red points are for the same with the goal noise levels, and the green points are the extended configuration with the goal noise levels. same as the goal noise level, and the green points are the extended configuration with the goal noise levels.
This figure shows that when polarized thermal dust is weak (roughly when A BB dust < 10 2 µK 2 ), the forecasted σ(r) is roughly constant, meaning that in the cleanest regions of the sky, we can obtain a similar BB detection performance, regardless of the exact foreground contamination. Then, for regions with stronger foregrounds, σ(r) increases slightly with A BB dust , until we reach regions where σ(r) increases exponentially, roughly where A BB dust > 10 2 µK 3 . Regarding the difference in σ(r) for regions with similar dust contamination, we have checked that this is an effect due to the uncertainty on A fore , either through a higher synchrotron contamination or more spatial complexity of the dust spectral parameters in the region.

CONCLUSIONS
In our work, we have determined the best regions of the sky to survey for primordial B-modes, and we find that the least dust-contaminated regions of the sky, close to the south Galactic pole, have a similar performance. We have performed this forecast for the σ(r) error bar for different survey areas over the full sky, using an SO-SAT-like experiment as our example. We consider the fiducial configurations that are scheduled to begin observations in the next couple of years, as well as analyzing the effect of adding the thermal dust-monitor high-frequency channels that will be available with the operation of the CCAT-prime project. In our analysis, we use the CCA method, which works in harmonic space minimizing the observed data covariance, by finding the optimal spectral parameters for the modeled foreground components. We perform our analysis over 200 Monte Carlo simulations and estimate the σ(r) error bar by using the Fisher information matrix. We perform forecasts for 192 circular masked regions covering the entire celestial sphere, each equivalent to f sky = 0.1, and we include the more realistic SO-SAT mask as well.
In our results the SO-SAT survey could measure r with an error bar of σ(r) = 4.9 − 5.2 × 10 −3 using the SO-SAT mask with the fiducial six-channel configuration, a baseline noise level, and depending on the method how we account for foreground residuals. When we consider the more sensitive goal noise levels, our forecast improves to σ(r) = 3.3 − 3.4 × 10 −3 using the SO-SAT mask. These values are considering a model where we account for the foreground residuals, and marginalizing over this increases the value of σ(r) somewhat. However, this paper should not be interpreted as an official SO forecast, as the aim of this work is different: to compare performance for distinct survey regions over the full sky.
We also analyze the potential impact of including three additional high-frequency channels from the CCAT-prime experiment that monitor the strong thermal dust emission. In our analysis, adding these extra high-frequency channels improves the measurement of the CMB BB spectrum slightly. We have investigated the reason for this, and conclude that a naive analysis will leak small dust residuals from these extra channels that have a very strong dust component. Even if the polarized dust residuals are very small, they are still quite significant for B-modes. Any mismatch in the thermal dust SED when we model it in our component separation will create residuals that will affect the reconstruction of the CMB B-mode spectrum. This analysis emphasizes the need to accurately model the SED of the CMB foregrounds, especially for the thermal dust polarization.
We attempt to answer the question posed at the beginning of this work, that is, to measure r, is it better to use the region that is cleanest of foregrounds, or is there a sweet spot where the foreground amplitudes are optimal from the point of view of facilitating their removal? We know the most important factor for mitigating foreground residuals is to model the SED correctly. From our analysis, there appears to be a range of sky regions where the foreground contamination is weak and the σ(r) error bar is roughly the same. This seems to support the notion that having the smallest possible amount of foregrounds and having a slightly higher amount of foreground will yield a similar measurement of r, meaning that component separation is doing its job by mitigating BB foreground contamination in a range. One does not need the absolute best possible region for a Bmode survey to yield a good result. Nevertheless, there is a limit to this, and for regions with higher foreground contamination, the achievable σ(r) increases.
We conclude by stressing the importance of characterizing the SED for polarized foregrounds in the component separation analysis, in particular for polarized thermal dust residuals, which will significantly bias a measurement of B-modes if not accounted for. From our results, we can conclude that the cleanest regions of the sky will yield the overall best results in terms of measuring the primordial B-modes. There is some tolerance for regions with a slightly higher polarized foreground contamination, because the component separation is able to fulfill its purpose. Our goal is to mitigate the polarized dust residual as much as possible, and to do this, we need to improve the modeling of polarized foregrounds SEDs in the component separation. In parallel, the instrumentation that improves the sensitivity is another important factor. The need to better understand foregrounds and their SED is imperative, and new experiments that will monitor ancillary frequencies for synchrotron and thermal dust, such as QUIJOTE (Génova-Santos et al. 2015), C-BASS (Jones et al. 2018), and the already-mentioned CCAT-prime will be vital in the endeavor to detect primordial B-modes.

ACKNOWLEDGMENTS
C.H.C. acknowledges the funding from Becas Chile/CONICYT. C.H.C. and K.M.H. are supported by NASA through ATP award NNX17AF87G and by NSF through AAG awards 1815887 and 2009870. We thank Stefano Camera, David Alonso, Christian Reichardt, and Michael D. Niemack for useful suggestions that improved this work. We thank the anonymous referee for their useful comments. This work uses public data from Simons Observatory. Some of the computing for this project was performed on the HPC cluster at the Research Computing Center at the Florida State University (FSU).