Detecting Primordial Gravitational Waves: a forecast study on optimizing frequency distribution of next generation ground-based CMB telescope

Probing primordial gravitational waves is one of the core scientific objectives of the next generation CMB polarization experiment. Integrating more detector modules on the focal plane and performing high accurate observations are the main directions of the next generation CMB polarization telescope, like CMB S4. Also, multi-band observation is required by foreground analysis and reduction, as it is understood that foregrounds have become the main obstacles of CMB polarization measurements. However, ground observation is limited by the atmospheric window and can be usually carried out in one or two bands, like what BICEP or Keck array have done in the south pole. In this paper, we forecast the sensitivity of tensor-to-scalar ratio r that may be achieved by a multi-frequency CMB polarization experiment, basing on which to provide guidance for further expanding frequency bands and optimize the focal plane of a telescope. At the same time, the realization of having two frequency bands in one atmospheric window is discussed. With fixed number of detectors, the simulation results show that, in order to get a good limit, more frequency bands are needed. Better constraints can be obtained when it includes at least three bands, i.e., one CMB channel (95 GHz) + one dust channel (high frequency) and one synchrotron channel (low frequency). For example, 41 + 95 + 220 GHz, which is better than only focusing around the CMB band, like 85 + 105 + 150 GHz, and 95 + 135 + 155 GHz, and this frequency combination is even better than the combination of 41 + 95 + 150 + 220 GHz. As CMB S4 plans to consider two frequency bands in each atmospheric window, and along this way, we find that one CMB band and more bands in synchrotron and dust channels are helpful, for example, 2 bands in lower frequency, 30 + 41 GHz, 2 bands in higher frequency, 220 + 270 GHz, i.e. 30 + 41 + 95 + 220 + 270 GHz, can get better constraints, and in this case, more detectors are asked to be assigned in the CMB channel.


Introduction
Detection and measurement of primordial gravitational waves (PGWs) provide a key test of inflation which has been regarded as a cornerstone of modern cosmology. Some of its predictions, like the flat universe, the Gaussian, adiabatic and nearly scale-invariant initial condition, etc, are in good agreement with current observational data. CMB polarization observation has the potential to reveal signatures of PGWs from inflation. There are two CMB polarization modes, polarized E-mode and B-mode, and for the primordial B-mode, which is mainly sourced from tensor perturbations of inflation on degree-scale, can serve as a unique observational window for the detection of primordial gravitational waves [1][2][3][4].
CMB study has gained lots of achievements in the past fifty years. With three generations of space mission, COBE, WMAP [5,6] and Planck [7,8], the accuracy of CMB measurement has been improved step by step, especially the temperature anisotropy, which has made a great contribution to pushing the study of modern cosmology into precision epoch. Now, CMB research is entering a new era, the socalled stage 3 [9][10][11][12], in which more detectors are planned for improving the sensitivity, and the key scientific goal is searching for primordial gravitational waves. Also, with the collaborative efforts of CMB-S4 [13,14], which are focused on ground-based CMB observations, the detector sensitivity on CMB will be logically improved ambitiously with at least ten times more detectors. It will develop larger detector arrays and carry out joint observation at Antarctica, Atacama Desert and maybe also Northern Hemisphere stations, so as to further enhance statistical samples and improve observation accuracy.
Despite its significance, in fact, detecting the PGWs is extremely challenging, not only the signal of primordial Bmode is very weak, theoretically, it is about at least two to four orders smaller than CMB temperature anisotropies, but also the measurement of primordial B-modes is heavily contaminated by the polarized foreground emission [15,16], i.e. synchrotron and thermal dust [17][18][19] from the Milky Way. In addition, lensing [20][21][22] B-modes generated from E-modes and intervening large scale structures will dominate on the sub-degree scale, which also makes obstacles for the detection. CMB signal itself is frequency independent, but the foreground radiation is frequency-dependent, in order to do foreground component subtraction, it is necessary to carry out observation in several bands to eliminate frequency-dependent foreground signals.
However, for ground-based CMB experiments [23][24][25][26][27][28], the frequency bands that could be observed on earth are limited to the microwave atmospheric windows. There are strong absorptions from oxygen and precipitable water vapor (PWV), also these latter contribute a lot to the brightness temperature of the detectors loading in the microwave band. Therefore, ground-based CMB experiments are usually carried out at high altitude with very dry air [29]. Even on high ridges, subject to the presence of the atmosphere, observations are usually performed only in a few band windows: below 50 GHz, around 95, 150 GHz and above 200 GHz [29,30]. For example, BICEP1 has observed in two frequencies of 100 GHz and 150 GHz, BICEP2 in 150GHz, even for South Pole telescope, a large aperture CMB camera, hold only three channels of 95, 150 and 220 GHz, centered around atmospheric windows, far fewer than space projects. The modularity of the focal plane allowed it to be broken into many different frequency configurations. Therefore, how to increase the observation frequency band reasonably in the limited atmospheric window is very important to carry out ground-based multi-frequency CMB observation in the future. Considering that setting two frequency bands in a wider atmospheric window can effectively increase the number of observation frequency bands, this can be a way to realize a ground-based CMB observation in multi-frequencies and could be important to characterize systematic effects.
In addition, based on current technology, the noise level of a single quantum sensor has almost reached the limit of photon noise [31][32][33]. In order to further improve the detection accuracy, the next generation of telescopes will continue to increase the number of integrated detectors on the focal plane to more than 10,000 [34][35][36][37][38]. Such a large array of detectors provides a precondition for multi-band observation, thus how to optimize the number of detectors distributed in each frequency channel becomes an unavoidable question.
In this paper, we plan to do such kind of estimation. As the crucial parameter for characterizing PGWs is tensor-toscalar ratio r , its uncertainty σ r from measurement should quantify the sensitivity of a CMB polarization experiment [14,39,40] to the detection of PGWs. We do forecast on σ r by adopting Fisher information matrix approach [41][42][43][44][45][46] for different design of the focal plane so as to make the comparison. The frequencies we consider include 30 GHz, 41 GHz, 95 GHz, 150 GHz, 220 GHz and 270 GHz. Additionally, in order to add in more frequency bands to get better constraints, we also consider the bands of 85 GHz, 105 GHz, 135 GHz and so on. These bands can be realized by antenna-coupled filters. Finally, by comparing the achieved σ r we estimate the detection sensitivity for each focal plane model and obtain the optimal solution for the frequency channels distribution. In the literature, there are lots of studies in this aspect and they show useful and interesting results [39,47,48].
This paper is organized as follows. In Sect. 2 we quantify polarized foreground components, synchrotron and thermal dust, and instrumental noise, they are frequency-dependent and contribute to the sensitivity of r . In fact, it is the balance between these components as well as weak gravitational lensing B-mode that will determine the final sensitivity of r detection. In Sect. 3, we briefly introduce atmospheric windows that limit ground-based CMB observation experiments. In Sect. 4, we introduce the Fisher matrix method for the prediction of σ r on the tensor-to-ration r in detail. We summarize the results and give a detailed analysis of the optimization of the focal plane in Sects. 5 and 6. Finally, we briefly conclude in Sect. 7.

Foreground models and instrumental noise
In fact, all polarization signals observed by CMB telescope include what we usually call foreground, i.e. radiation from the Milky Way itself, the atmospheric emission, telescope instrument radiation and so on, will be included in the observation data. Foreground radiation, such as Galactic synchrotron radiation, thermal dust radiation can dominate at a certain range of frequency channels and will interfere with CMB polarization observation. The effects can be mitigated by foregrounds reduction treatment and the efficiency will greatly affect CMB signal detection. In the study of this paper, we mainly consider the thermal dust and synchrotron as these two components are very strong in the process of polarization observation [41,49,50]. Foreground simulation and reduction have been extensively studied in the literature (e.g., [44,51]).
The noise, generated by atmosphere and instrument, is also the item that must be controlled to improve the constraint on r . This section mainly gives a detailed description of the foreground and noise.

Synchrotron and thermal dust
Synchrotron is a significant component in the polarized foreground emission and mainly dominates at low frequency [5,52,53]. Synchrotron emission is sourced by relativistic cosmic ray electrons spiraling in the Galactic magnetic field, and it contributes to both temperature and polarization measurements. The intensity and power spectrum of synchrotron rely on cosmic ray density and magnetic field strength [52]. Therefore, the spectral index of synchrotron changes with the observed position and frequency, e.g., β s ≈ −2.7 around 1 GHz while it is close to −3 at higher frequencies given by the recent Planck results [54,55]. Fortunately, both theory and experiment indicate that the synchrotron spectrum is approximate to the power-law spectrum at the frequencies above 20 GHz [56,57]. Therefore the synchrotron spectrum is usually parameterized as with reference frequency ν s = 30 GHz and reference multipole s = 80. Where superscript 'p' stands for polarization. The reference frequency ν s can also be taken as 0.408 GHz, see the study in reference [40,58]. As the power spectrum is a region-dependent quantity, the result from different region statistics will also be different. We will consider two cases of a large sky coverage with f sky = 0.7 and a small sky coverage with f sky = 0.05, the parameters of synchrotron and thermal dust are summarized in Table 1. Thermal dust is another non-negligible foreground polarization emission, which mainly contributes at high frequencies [16]. We adopt a Modified Black Body for dust polarization spectra following Planck Collaboration at T d = 19.6 K [59] where d = 80, ν d = 353 GHz and B ν (T ) represents Planck function, which is given by Note that, it is more common to use antenna temperature for describing foregrounds and hence one needs to convert the unit from antenna to thermodynamic temperature by g ν = x 2 e x /(e x − 1) 2 , where x = hν/(kT CMB ), T CMB = 2.73K . Therefore, substituting g ν into Eqs. (1) and (2), finally we get and where B ν (T ) = e hν/kT − 1.
Since our goal is to estimate the uncertainty on r by Fisher information matrix, the temperature spectra are not within our consideration.

Instrumental noise
In addition to foregrounds, instrumental noise is another nonnegligible contamination in the measurement of primordial polarized B-modes. Instrumental noise is usually assumed to be white and Gaussian, for a multi-frequency experiment with m frequencies, its total noise power spectra read as [39,60]: where N XX,inst ,m denotes the instrumental noise power spectrum with m th frequency of X = E, B given by: The subscript m in w and θ 2 indicate that they are frequencydependent. In addition, the expression of w is where Y refers to detector efficiency. In our study, we assume an ideal case where yield Y = 1. Finally, we calculate the value of the Noise Equivalent Temperature (NET) which is frequency-dependent as we will see below.

Noise equivalent temperature
In general, the sensitivity of the detector is normally quoted as NET. In order to calculate the spectrum of instrumental noise in Eq. (7), we need to calculate the NET first (with units μK √ s). For CMB experiments, using Planck blackbody equation: NET which is related to the optical transfer function τ (ν) and the atmospheric opacity (ν) can be finally expressed as [33,61]: where NEP (with units W Hz −1/2 ) is defined at the level of absorbed power, and where NEP photon and NEP phonon represent photon noise and phonon noise, respectively. Since noise is mainly dominated by the former, only photon noise will be considered below and it can be written approximately as [62] where P load is the power in a band of width ν, and it is directly related to the optical efficiency η, the effective bandwidth ν, as well as the internal power of instrument Finally, for our cases, the NET for each frequency are shown in Table 2.

Angular power spectra
In this section, we will visually display the spectra through the figure. We have shown the angular power spectra ( + 1) C BB /2π of foregrounds, primordial B-modes, lensing and instrumental noise in Fig. 1. For comparison, we plot the spectra of 41 GHz, 95 GHz and 150 GHz in Fig. 1, where the latest results of Planck are used for the foreground spectra [59]. Picture (a), (c) and (e) represent the angular power spectra for sky coverage of 5% and picture (b), (d) and (f) represent the angular power spectra for sky coverage of 70%. Comparing the synchrotron and thermal dust, the figure shows that the synchrotron dominates at 41 GHz while the dust dominates at 95 GHz and 150 GHz. Moreover, the instrumental noise increases exponentially at high-multipole as shown in Eq. (7).

Atmospheric windows
CMB ground-based experiments are sensitive to the atmosphere, which is a significant factor when determining the site of such experiments. On the one hand, the atmosphere absorbs CMB photons, which will reduce the CMB signal. On the other hand, atmospheric molecules can radiate photons at millimeter and sub-millimeter wavelengths to increase photon noise. In order to give a more intuitive description of the atmospheric window, we use the am code [63] to plot the curve of the atmospheric transmittance versus frequency shown in Fig. 2. The atmospheric window is divided into 4 parts by strong oxygen absorption lines around 60 GHz and 120 GHz and water vapor lines at 183 GHz as shown in the figure: 41 GHz band, 95 GHz band, 150 GHz band and a high-frequency band 220 GHz. In addition, the figure shows that the atmospheric transmittance decreases with the increase of precipitable water vapor (PWV). From these windows, in the following section, we will discuss the effect of adding more bands in one atmospheric window to constrain r . Moreover, we also plot the curve of the brightness temperature versus frequency at various values of PWV shown in Fig. 3. As we can see, the atmospheric brightness temperature is sensitive to PWV, especially the frequency above 150 GHz, and it increases with the increase of PWV.

The Fisher matrix formalism
The Fisher matrix formalism is generally considered as a relatively simple method to forecast the error in constraining parameters for future experiments, which is widely used in constraining cosmological parameters from data sets. In this section, the Fisher matrix approach is adopted for forecasting the uncertainty on the tensor-to-scalar ratio r . Assuming that the probability distribution of parameter x to be Gaussian, so the likelihood function for the n-dimensional random variable with mean vector μ is where C is the covariance and which can be expressed as [64].
Fisher matrix is the second-order term of likelihood function L which peaks at the fiducial value − → θ = − → θ fid after Taylor where θ i denotes ith parameter of parameter set θ. For CMB spectra, F i j can be expressed as [41,43] where f sky represents the sky coverage fraction and 'tr' denotes the trace of matrix. By definition, the conditional errors of cosmological parameters are given by the inverse of the Fisher matrix and the marginalized error on parameter θ i following from the Cramer-Rao inequality is given by [65] Now, let us try to illustrate the structure of covariance matrix C . As usual, we assume the covariance among theoretical CMB, foregrounds and instrumental noise are uncorrelated. Therefore, the covariance matrix C can be expressed as a m a t m , where a m is a vector, a m = a T , a E , a B m . Since foregrounds are considered, for each TT, EE, BB modes, there are three components in a m . For example, the harmonic coefficients of a CMB TT map can be written as (see for example [66] where a T,th m refers to the sum of primary CMB signal and lensing, while a T, f g m and a T,n m refer to the foregrounds and instrumental noise, respectively. For multi-frequency experiments, it is worthwhile to note that, in general, the angular power spectrum of instrumental noise is considered to be uncorrelated among different frequency channels, i.e., a X,n,ν m a X ,n,ν m = N X X,νν δ δ mm δ νν δ X X , where N X X,νν represents the frequency-dependent angular power spectrum of instrumental noise for X = T, E, B. Since the CMB signals, foregrounds and instrumental noise are assumed to be uncorrelated, we can get the expression for any frequency channels i and j Finally, based on the above assumptions, the covariance matrix C = a m a T m can be written as [43] If only BB power spectra are considered, for two frequency channels ν i and ν j Specifically, as only two foreground polarization components (i.e., synchrotron and thermal dust) are concerned in our study. Therefore, So far, we can obtain the error of r and foreground parameters by using Eqs. (17), (19) and (23). Parameters in Fisher calculation are summarized in the following: Foreground parameters: we take them to be free parameters when performing error estimation through Fisher matrix approach.
CMB parameters: {r } and standard cosmological model. Priors: we adopt the following priors on parameters of foreground components given by Planck publication [59] (a) Small sky coverage (41 GHz With the priors, we then can represent them with a diagonal covariance matrix and invert this matrix into a Fisher matrix as follows [67,68] and finally we add each diagonal element of this matrix into the appropriate diagonal element of the original Fisher matrix in Eq. (17).

Multipoles
Let's estimate the approximate range of multipoles first in this section before the numerical calculation of σ r . The primordial CMB temperature anisotropy has been well measured to multipoles ≈ 3000 at current stage [13], but the polarization anisotropy measurements are nowhere near as good as that of the temperature. However, due to high angular resolution and sensitivity of detectors, the scientific goals have set to P ≈ 5000 in the CMB-S4 for CMB polarization anisotropy. For , there's a minimum limit for a ground-based CMB experiment due to the presence of atmosphere, and also as the experiments on the ground cannot realize a full sky observation like satellite experiments, and therefore we take min = 50 in our study. On the other hand, we can calculate max by = max =0 (2 + 1) = N pix . For example, consider a full sky experiment with an angular resolution of 10 arcmin, which will give us about 1.5 million independent pixels, and we finally get max 1200. However, for constraining r we do not have to set max to be a quite large number. For a full sky coverage experiment if only lensing contamination is considered, then the multipoles greater than 150 will not contribute to the constraints on r [69]. When we take foregrounds and instrumental noise into account, as shown in Fig. 4 that means we only need to expand multipoles up to 600.

Forecast r and optimizing focal plane
The tensor-to-scalar ratio is a key parameter to characterize the primordial gravitational waves, and its uncertainty, σ r , given by data can be a sensitivity benchmark of a CMB polarization experiment. Therefore, we can optimize the detector distribution of the focal plane in the frequency-band by comparing the σ r that can be given by different focal plane detector integration schemes. Here, we consider a telescope with 10,000 detectors on the focal plane as a representative of the next generation ground-based CMB polarization telescope to study its band combination scheme and optimization of the number of detectors in each band distribution. It is worth noting that for space missions, the weight/size of detectors may be an important parameter because it varies for different frequency bands, whereas for ground-based experiments, it is not an important factor. For simplicity, we do not consider this factor in this paper.
The method of optimizing focal plane is based on the estimation of σ r for each model, since the smaller the σ r , the higher the detection sensitivity, thus the model which can provide the minimum σ r will be the best scheme. We consider two, three, four and five frequency bands on the focal plane models respectively, and we loop the detector numbers in each frequency channel while keeping the total detector number equal to N total = 10,000 when estimating the performance of the model.
In evaluating σ r , it is actually the process of converting the uncertainty of foreground, CMB lensing and instrumental noise into σ r . Foreground components are region dependent quantity, so they are related to the choice of f sky , like the description of dust power spectra listed in table 1 of [59], the amplitude of the power spectra will change for different f sky . Here, we consider two cases of f sky selection: a large sky coverage with f sky = 0.7 and a small sky coverage with f sky = 0.05.
Two different fiducial values of r (0.01 and 0) are considered, with which we simulate the power spectra for the Fisher matrix. We will start from a two-frequency model which contains 95 GHz and 150 GHz as our basic model, Case 0. Further, we will consider the combination with more frequency bands, such as three-frequency models: Case 1, Case 1*, Case 1**, Case 3 and Case 3*, four-frequency models: Case 2, Case 4 and Case 4* and five-frequency models: Case II and Case IV.

Results for large sky observation
In this section, we present the results, which include the following cases of focal plane schemes: Case 0: two frequency bands of 95 GHz and 150 GHz. Case 1, three frequency bands of 41 + 95 + 150 GHz, in which we include 41 GHz band in the synchrotron channel (low frequency) based on Case 0. Case 1*: three frequency bands of 95 + 150 + 220 GHz, in which we include 220 GHz band in the dust channel (high frequency).
Case 1**: three frequency bands of 41 + 95 + 220 GHz, in which we include a low band 41 GHz in the synchrotron channel and a high band 220 GHz in the dust channel simultaneously.
Case 2: four frequency bands of 41 + 95 + 150 + 220 GHz. While optimizing these focal plane models, we cycle the detector numbers N f i for each frequency band in the range of (0, 10,000), with the condition that f i i=1 N f i = 10,000 and step size of 200 to get the minimum σ r for each case.
We summarize the numerical results in Table 3. The first column shows the fiducial value of r while the second column shows the minimum σ r and the rest of the columns refer to the number distribution of detectors for each channel corresponding to the minimum σ r .
In Case 0, taking an example of r = 0.01, the optimal configuration of detector distribution is N 95 = 7400 and N 150 = 2600, which means that we can obtain a minimum σ r = 0.01431 if N 95 = 7400 and N 150 = 2600. This conclusion is also shown in (a) of Fig. 5, illustrating that as the increase of N 95 , σ r gradually decreases and until N 95 is close to 7400, σ r reaches the minimum value of 0.01431. Case 0 also shows that more detectors are needed to be assigned to 95 GHz channel in order to get a smaller σ r . The main reason is that the distribution of foreground contamination (dust plus synchrotron) and NET are relatively low around 95 GHz [48]. If we add a lower channel on the basis of Case 0, taking 41 GHz as an example in Case 1, the optimal configuration is N 41 = 600, N 95 = 6800, N 150 = 2600. It is clear that σ r has been significantly improved after adding 41 GHz. Taking r = 0.01 as an example, the limitation of σ r is decreased from 0.01431 to 0.00893 with an decrease rate of 38%. The reason is that when a low frequency is added, the parameters of synchrotron can be better constrained so that we can get a better result. As can be seen from Fig. 1, thermal dust has a large amount of contamination over a wide range of frequencies, especially at high frequencies. Therefore, adding a high-frequency channel to constrain the parameters of thermal dust should also be very helpful for improving σ r . In Case 1*, we add a high band, 220 GHz, which indeed decreases the σ r from 0.01431 to 0.00778 for the case of r = 0.01, with an decrease rate of 46% and the optimal configuration changes to N 95 = 4200, N 150 = 3000, N 220 = 2800 at the same time. Comparing Case 0 (95 + 150), Case 1 (41 + 95 + 150) and Case 1* (95 + 150 + 220), the results show that it is better to add a high-frequency channel like 220 GHz GHz. There is a great improvement compared to Case 0 as well while no improvement compared to Case 1**. Besides, since the step size is set to 200, this means the minimum number of detectors is 200 instead of 0 for each band. As Case 2 shows that N 150 is equal to 200, it indicates that the results may be improved if we remove 150 GHz and the results are shown in Case 1** which do have a little improvement. This means that we indeed should assign more detectors at 95 GHz and remove the band of 150 GHz to better constrain r when the total number of detectors is fixed. Finally, in order to compare these cases more intuitively, we show the results in bar chart (a) of Fig. 6. The chart shows clearly that the combination {41 + 95 + 220} GHz can give us the best limitation on r .

Adding more frequency channels within one atmospheric window
The design of the antenna-coupled filter can effectively realize the observation of two or more frequency bands in the same atmospheric window. The number of bands depends on the bandwidth design of the filter. This will help ground-based CMB telescope to achieve multi-band observa- Here, we try to increase the frequency bands within one atmospheric window. For example, we replace the 95 GHz band with a slightly lower frequency, 85 GHz, and a slightly higher frequency, 105 GHz, which belongs to the same atmospheric window shown in Fig. 2, to increase the frequency bands.
Similarly, we can also use 135 GHz and 155 GHz to replace 150 GHz. Based on these considerations, we further simulate the observations with more frequency bands, as listed in the following of Case 3, 3*, 4, 4*, and analyze their observational sensitivity. Now, let us calculate whether it is helpful to improve the detection sensitivity σ r after considering two bands in the same atmospheric window. The results are summarized in Table 4 and visualized in bar chart (b) of Fig. 6. Comparing the cases in Table 4 with Case 0 (95 + 150) in Table 3, it can be seen that there is no improvement for Case 3 while for Case 3*,4 and 4*, σ r are improved by 4.6%, 2.4% and 40.3% for fiducial r = 0.01, respectively. The results also indicate that it is advantageous to increase the frequency bands both at low and high frequencies. For example, comparing Case 4* (41 + 85 + 105 + 150) with Case 3 (85 + 105 + 150), σ r is decreased from 0.01492 to 0.00855 with an decrease around 43%. Meanwhile, the chart (b) also indicates that setting two bands within 95 GHz window or 150 GHz window is not very helpful (i.e., 85 and 105 GHz in 95 GHz window and 135, 155 GHz in 150 GHz window), which is much less than adding bands at both low and high frequencies, (e.g., {41 + 95 + 220} GHz). Again, the reason is that we can better constrain foreground parameters by observing at both lower and higher frequency bands.

Results for small sky observation
In this section, we will focus on a small sky coverage with sky fraction f sky = 0.05. As shown in Table 1, compared to the large sky coverage, small sky area has a smaller amplitude of foreground polarization contamination. In addition, the instrumental noise will also vary with the size of the sky coverage since it is the function of f sky (see Eq. 8). As we will show below, the conclusions in this section are basically the same as the previous section.
We summarize the results in Table 5 for Case 0, 1, 1*, 1** and Case 2 and list the minimum σ r as well as the number distribution of detectors. Moreover, the visualization of the table is available in Fig. 7. As shown in the table, for example, for r = 0.01, in order to obtain a minimum σ r in small sky area, the optimal combination in Case 0 is N 95 = 7400 and N 150 = 2600, which also requires more detectors to be assigned to 95 GHz channel like large sky coverage we discussed before. In addition, we can see that compared to Case 0 if a lower band like 41 GHz (Case 1) or a higher band like 220 GHz (Case 1*) is included, σ r will have an improvement of 10% (Case 1) or 23% (Case 1*) on r = 0.01. Similarly, if both high and low bands are considered, there will be a higher improvement such as Case 1** (41 + 95 + 220) and Case 2 (41 + 95 + 150 + 220). The reason is also that the foregrounds can be better limited after observing both high and low bands simultaneously. In addition, the bar chart (a) in Fig. 7 shows that adding 220 GHz is better than adding 41 GHz when expanding 2 bands (Case 0) to 3 bands (Case 1, 1*). Finally, we plot the curve of σ r with N 95 in (b) of Fig. 5 for Case 0, which indicates that σ r will decrease with the increase of N 95 and it will reach a minimum value of 0.00467 when N 95 is equal to 7400. Fig. 6 The comparison of σ r for different frequency combination models: two-frequency model, three-frequency model and four-frequency model for fiducial r = 0 and r = 0.01, respectively. The bar chart a is the visualization of Table 3, which does not include the case where more than one bands are contained in one atmospheric window shown in Fig. 2 while b considers this case (a) (b)

Adding more frequency channels within one atmospheric window
As large sky coverage, we try to add more bands within one atmospheric window in this section for small sky coverage with f sky = 0.05. As we will show below, the 95 GHz band and 150 GHz band will be replaced by {85+105} GHz and {135+155} GHz, respectively. As before, the cases in this section are: Case 3: 85+105+150 GHz Case 3*: 95+135+155 GHz Case 4: 85+105+135+155 GHz Case 4*: 41+85+105+150 GHz The results are summarized in Table 6. Comparing Case 0 (95 + 150) to Case 3 (85 + 105 + 150) and Case 3* (95 + 135 + 155), there is not much improvement on σ r by replacing 95 GHz with {85,105} GHz while a slight improvement to replace 150 GHz with {135,155} GHz. In addition, comparing Case 0 (95 + 150) to Case 4 (85 + 105 + 135 + 155), σ r is decreased from 0.00467 to 0.00456 with only 2.4% improvement, which again indicates that adding more frequency bands around 95 GHz and 150 GHz is not efficient for improving the constraint on r . The main reason is that in Case 4 there is not enough foreground information from the lower and higher frequency bands. As shown in (b) of Fig. 7, the combination {41+95+220} is still the best choice to improve σ r for small sky region.

Further improvement
In the CMB S4 experiments, four atmospheric windows will contain eight frequency bands, each window contains two frequency bands. We have found that adding a low band (41 GHz) or a high band (220 GHz) can obviously improve σ r . Since it is very helpful for limiting the contamination from synchrotron and thermal dust. Therefore, in this section, we will discuss the effect of further adding channels in lower frequency bands and higher frequency bands. As the calculation above shows that the channel with {41 + 95 + 220} GHz can give us an optimal σ r . We therefore further consider the case where the windows of 41  Tables 3 and 5, respectively, and they have been renamed here for convenience. Comparing Case I with Case II, σ r is decreased from 0.00422 to 0.00361 for r = 0.01 with an improvement of 14.5%, while comparing Case III with Case IV, σ r is decreased from 0.00319 to 0.00308 with an improvement of 3.4%. To be more intuitive, we visualize the results in Fig. 8 as well. In the bar chart of Fig. 8, we compare Case 0 (95 + 150) of Tables 3 and 5 with the case of (41 + 95 + 220) and (30 + 41 + 95 + 220 + 270), it clearly shows that it is crucial to contain both lower (41 GHz) and higher (220 GHz) bands and it is helpful to split them into 2 bands, {30, 41} GHz and {220, 270} GHz. Finally, it is worth noting that our conclusions are obtained by only considering two main Galactic components, thermal dust and synchrotron.

Constraints on foreground parameters
In this section, we discuss the constraints on foreground parameters and their correlation with r . The results of 68% limits with a fiducial model of r = 0.01 are summarized in Table 7. It shows that the parameters of foreground can get better constraints when the lower and higher frequency channels are taken into account, in detail, synchrotron can be further limited when adding a low frequency, 41 GHz, Fig. 8 The visualization of Table 8. The bar chart a and b represent the case of f sky = 0.7 and f sky = 0.05, respectively. Comparing the three Cases in the figure, the combination {30+41+95+220+270} can give us the optimal σ r (a) (b) Moreover, our results show that the correlation between r and the foreground parameters can not be neglected, as shown in Fig. 9, we plot the contours of foreground parameters for the frequency combination of {95 + 150} GHz and {41 + 95 + 220} GHz. We have noticed that the correlation is related to f sky and the number of frequency bands.

Summary
Galactic polarization emission from the Milky Way, such as synchrotron and thermal dust, is one of the main sources of pollution in CMB polarization measurements. Obtaining more information about foreground polarization radiation is very important for foreground reduction, so as to detect primordial B-modes. In this paper, we mainly study how to optimize the frequency band distribution and the number of detectors of each band for CMB polarization observation in the future.
The Fisher matrix method is adopted to simulate the sensitivity of a focal plane system with ten-thousand level detectors, and to judge the detection efficiency of the focal plane. First of all, the focal plane models with different frequency bands are calculated where the total number of detectors is fixed, such as the two-frequency model (95 + 150 GHz), the three-frequency model (41 + 95 + 220 GHz), the fourfrequency model (41 + 95 + 150 + 220 GHz) and fivefrequency model (30 + 41 + 95 + 220 + 270 GHz), the corresponding sensitivity are calculated.
The numerical results show that the combination {41 + 95 + 220} GHz is always the best to limit r when only considering three-frequency model, even better than {40 + 95 + 150 + 220} GHz. When we add more frequency bands in synchrotron channel (low frequency) and dust channel (high frequency), σ r can be further improved, e.g., {30 + 41 + 95 + 220 + 270} GHz, which is better than only focusing more on the CMB channel, like {85+105+150} GHz, {95 + 135 + 155} GHz and {85 + 105 + 135 + 155} GHz. This is because synchrotron can be better constrained by the observation at low frequency and thermal dust at high frequency. Therefore, the detection accuracy of r can be improved after better limiting foregrounds. Moreover, we found that in order to obtain better constraints on r , more detectors are needed to be assigned around CMB channel, as shown in the case of (a) 95+150 (b) 41+95+220 Fig. 9 The contours and the likelihood distributions for the foreground parameters and r at 68%, 95% and 99% CL of {95+150} GHz (a) and {41+95+220} (b) with f sky = 0.7 and fiducial r = 0.01 {41 + 95 + 220} GHz, more than 7000 detectors are asked to observe at 95 GHz.