Performance forecasts for the primordial gravitational wave detection pipelines for AliCPT-1

The key of SZ method is to construct two mutually orthogonal scalar (pseudo-scalar) fields $\mathcal{E}$ and $\mathcal{B}$. For an incomplete sky observation case, we can define the partial-sky harmonic coefficients of the $\mathcal{E}$- and $\mathcal{B}$-fields (with overhead tilde) as [66],

where $W(\hat{n})$ is the window function. These expressions can be expanded and simplified further for implementation and full form expressions can be found in [67].

The partial-sky pure harmonic coefficients $\tilde{\mathcal{E}}_{\ell m}$ and $\tilde{\mathcal{B}}_{\ell m}$ are related to the full-sky harmonic coefficients $E_{\ell m}$ and $B_{\ell m}$ as follows,

where $K^{ru}_{\ell^{\prime}m^{\prime}\ell m}$, with $r,u\in\{E,\ B\}$, represents the mixing kernels of pure fields. And then we can define pseudo estimators of $\mathcal{E}$- and $\mathcal{B}$- fields as,

They relate to the CMB power spectra $C_{\ell}^{EE}$ and $C_{\ell}^{BB}$ by the following relations

where the notation $\langle\cdot\cdot\cdot\rangle$ denotes the ensemble average. The exact expressions for the mixing kernels $K^{ru}_{\ell^{\prime}m^{\prime}\ell m}$ and mixing matrices $\mathcal{M}^{ru}_{\ell\ell^{\prime}}$ can be found in [68, 69]. Therefore, we can construct the unbiased estimators of the actual power spectrum as follows,

This estimator is adopted in the McMfL pipeline. A slight variation to the PCL-SZ estimator is to construct pure-$B$ modes but not pure-$E$ modes. Then we can get a mixing matrix relating the partial sky $\begin{pmatrix}\tilde{C}_{\ell}^{EE}&\tilde{C}^{\mathcal{BB}}_{\ell}\end{pmatrix}^{t}$ and the actual power spectrum $\begin{pmatrix}\hat{C}_{\ell}^{EE}&\hat{C}^{BB}_{\ell}\end{pmatrix}^{t}$. This variant of the PCL-SZ estimator has been used in the ABS and the TF pipelines. In this work we use the PCL-SZ method implemented in python package of the NaMaster [70].

In this work, we also adopt a novel estimator to reconstruct $B$-mode power spectrum form $Q$, $U$ maps, called the PCL-TC estimator. This method uses the template cleaning (TC) method proposed in [71, 62, 64] to obtain a leakage free $B$-mode map. The main idea of the TC method is to use the high-precision $E$-mode signal to construct the $E$-to-$B$ leakage template, and then subtract a linear fit of this leakage template from the original $B$-mode signal. This gives us the $B$-mode map in our observation window without the partial sky $E$-to-$B$ leakage.

Since the cleaned $B$ map is essentially the isolated $B$-mode pseudo-scalar field, we can apply the scalar PCL estimator to reconstruct its power spectrum. For the cleaned $B$ map, the pseudo-multipoles are defined as

From pseudo-multipole coefficients, we can define the following PCL power spectrum estimate

By inverting the mode-mode coupling matrix between pseudo power spectrum and CMB power spectrum, we can built an unbiased estimator of CMB $B$-mode power spectrum

where

Here $w_{\ell}$ is the power spectrum of the window function $W(\hat{n})$. The PCL-TC estimator is used in the GLS and cILC pipelines. The scalar PCL is implemented with NaMaster.

This filtered $B$-mode power can be written in terms of unfiltered $E$- and $B$-mode power as:

Here $F_{BB,\ell}^{2}$ and $F_{BE,\ell}^{2}$ represent the auto and cross-transfer functions. The auto transfer is responsible for power suppression and the cross transfer is responsible for leakage. The transfer functions are computed from 50 filtered CMB only simulations. The transfer functions are computed as:

To compute the transfer functions, the CMB simulations from ancillary data discussed in Sec. 2 are used. In this work, two different computation approaches are adopted. For the McMfL pipeline, all the power spectra calculated from the maps are using the PCL-SZ method. The leakage power ($D^{BB}_{\ell,L}$) and filtered $E$-mode power spectrum ($D^{EE}_{\ell,\text{ filt}}$), are computed from $T$- and $E$-only filtered CMB maps. The power suppressed $B$-mode spectrum $\bar{D}^{BB}_{\ell,S}$ is computed by subtracting $\bar{D}^{BB}_{\ell,L}$ from the $B$-mode spectrum of the filtered CMB maps. The unfiltered $B$-mode spectrum, $D^{BB}_{\ell}$ is computed from the corresponding unfiltered CMB realization. For all cases the BPSapo mask described in Sec. A is used. In Fig. 2 we show the mean auto- and cross-transfer function used in the McMfL pipeline for all frequency channel pairs. We find some difference only for $\ell<40$.

The GLS and cILC pipelines do not use $F^{2}_{BE,\ell}$ to correct for the $E$-to-$B$ leakage. The leakage correction used in these pipelines is discussed in the next section. For computing $F^{2}_{BB,\ell}$ for GLS and cILC, we compute $\bar{D}^{BB}_{\ell,S}$ from $B$-mode only filtered $QU$ map using PCL-TC method with the inverse noise weighted mask, UNPinv discussed in Sec. A. The unfiltered spectra $D^{BB}_{\ell}$ are computed from the unfiltered full-sky input map. The full-sky power spectrum is computed with anafast function from healpy and binned with the same binning used in PCL-TC pipeline. The transfer functions are computed for each pair of input full sky and filtered CMB maps and use the mean transfer function from 50 simulations in our power spectrum correction. In the left plot of Fig. 3 we show the mean and standard deviation of the auto transfer functions $F^{2}_{BB,\ell}$ for AliCPT 90 GHz and 150 GHz channels. In the multipole range of our interest $F^{2}_{BB,\ell}$ is identical for all channels. The small number of simulations used in these calculations can impact the estimates for the transfer functions. However, we expect the mean estimates of the transfer function to converge fairly soon. In future studies, we intend to repeat the calculations with significant increase in number of simulations.

Correcting the leakage using $F^{2}_{BE,\ell}$ is the standard process of correcting leakage. We consider a novel option for filtering leakage corrections using prior information from Planck $E$-maps to construct a template of the leakage. In the case of the GLS and the cILC pipelines where we obtain a foreground cleaned $B$-map, in harmonic space we can write:

where $B^{\rm GLS/cILC}_{\ell m,\text{ filt}}$ is the GLS/cILC cleaned filtered $B$-mode harmonic coefficient, $\bar{B}^{s}_{\ell m}$ is the suppressed $B$ modes of the CMB signal, $\bar{B}^{n}_{\ell m}$ is the suppressed residual noise in the cleaned $B$-mode, $\bar{B}^{f}_{\ell m}$ is the suppressed residual foreground in the cleaned $B$-mode, $L^{s}_{\ell m}$ is the leakage term of the CMB $E$ modes due to filtering and $L^{n\prime}_{\ell m}$ is the projected leakage of the $E$-mode noise. For our analysis the main contribution to the cleaned $B$-map comes from the AliCPT channel, implying $L^{n\prime}_{\ell m}<L^{s}_{\ell m}$, as AliCPT 90 and 150 GHz noise level is smaller than the $E$-mode CMB signal. Also note that $B^{\rm GLS/cILC}_{\ell m,\text{ filt}}$ should pick up a leakage contribution from the $E$-mode foregrounds, but this should be removed by foreground cleaning methods as it will have different frequency scaling compared with the CMB. So any residual leakage from foreground should be small.

The Planck experiment measures $E$-modes which are independent of filtering effects of AliCPT. The Planck $E$-mode maps can therefore be propagated through the AliCPT observation pipeline to be ‘reconstructed’ to get the $B$-mode leakage template $\bm{L}$. First we obtain the NILC cleaned $T$- and $E$-modes from Planck HFI and WMAP K band combination. Then we get the $IQU$ map from these $T$- and $E$-modes. This $IQU$ map is then ‘reconstructed’ by propagating through the AliCPT simulation pipeline. The $B$-map obtained from the filtered $T$ and $E$ mode only $IQU$ map is the leakage template. The leakage template used in the DC1 analysis is shown in the right plot of Fig. 3. Note that Planck $E$-modes have much larger noise contamination than AliCPT. The filtering effect is linear, so in harmonic space:

where the spherical harmonic coefficients $L_{\ell m}$ of the template is a sum of the contribution from the $E$-mode signal $L^{s}_{\ell m}$, and the $E$-mode noise $L^{n}_{\ell m}$. Since we have filtered the NILC cleaned maps to construct the template, the leakage of any residual $E$-mode foreground should be very small. The noise leakage term $L^{n}_{\ell m}$ is different from $L^{n\prime}_{\ell m}$ in Eq. (4.4), as the dominant contribution in the later is from the AliCPT channels which are not used in constructing the leakage template. Finally the cross spectrum of $\bm{B}^{\rm GLS/cILC}_{\text{filt}}$ and $\bm{L}$ will give us the leakage bias at power spectrum level:

There will be a small correlation between $L^{n}_{\ell m}$ and $L^{n\prime}_{\ell m}$ but we found it to be small as $L^{n\prime}_{\ell m}$ mostly gets contribution from the AliCPT channels which are absent in $L^{n}_{\ell m}$. We subtract $C_{\ell}^{L\text{ bias}}$ to correct for the $E$-to-$B$ filtering leakage. The contribution of $L^{n\prime}_{\ell m}$ to the final power spectrum is smaller than that of $L^{s}_{\ell m}$ and is corrected in the mean from the noise debiasing performed with filtered noise simulations. We debias the power spectrum obtained from the foreground cleaned maps for the noise bias and the leakage bias. Finally the power suppression is corrected by dividing the binned debiased power spectrum by the auto transfer functions, $F^{2}_{BB,\ell}$. The Planck based correction method has been adopted in the GLS and cILC pipelines.

In general the power spectrum of a filtered observation map can be written as:

where $c$ is the index over various emission components and $N_{\ell,\text{ filt}}^{BB}$ is the filtered $B$-mode noise power spectrum. In case of the GLS and cILC method we assume that the map contains only the CMB signal, but for McMfL it is CMB and foreground. The noise bias $N_{\ell,\text{ filt}}^{BB}$ is computed from the 50 noise simulations mentioned in Sec. 2. The estimate of the corrected $B$-mode spectrum is then given by:

The leakage contribution $\sum_{c}D^{BB}_{\ell,L,c}$, for McMfL is computed as $F_{BE,\ell}^{2}\sum_{c}D^{EE}_{\ell,c}$, using the $E$ modes of the CMB and foregrounds. While, in the GLS or cILC pipeline we use $D_{\ell}^{L\text{ bias}}$ computed with the leakage template as $\sum_{c}D^{BB}_{\ell,L,c}$. For the McMfL pipeline this procedure is applied to each $\nu_{1}\times\nu_{2}$ spectra, but we do not show the frequency indices explicitly in Eq. (4.8).

The filtering corrections using transfer function and cross spectra from leakage template represents a ‘backward correction process’ to recover the corrected $B$-mode power spectrum. Alternatively, we can completely avoid the correction of the $B$-mode power spectrum in favor of the ‘forward modeling’ of the filtering effects, in fitting the CMB $B$-mode parameters like $r$ and $A_{\rm lens}$. We use this ‘forward modeling’ of filtering effects for the ABS and template-fitting parameter estimations. The method is discussed in detail in Paper III, and also is elaborated later in Sec. 5.1.

Now we will provide a brief summary of the forward modelling approach. The filtering effects are determined by a linear mapping of the filtering matrix. Thus, one can use a set of simulations to construct several templates for each independent component of the filtered CMB, including the tensor-only component, the $E$-to-$B$ leakage component, and the lensing one. The tensor-only component is built by the mean of the band powers by averaging over 50 realizations of ‘reconstructed’ maps, each of them being generated with a fiducial value of $r=0.01$ and entering into the observation pipeline. Thus, an arbitrary tensor-only band power can be achieved by simple rescaling such template by a factor of $r/0.01$, so that by fitting the parameter $r$, one can directly match the template to data without correcting for filtering effects. In the same way, we can build for the $E$-to-$B$ leakage template by using 50 realization $TE$-only maps, which are generated by the fiducial $\Lambda$CDM model. The $E$-to-$B$ filtering leakage contribution in the $BB$ powers does not depend on $r$, and we can compute the leakage contribution using $E$-mode only simulations. We correct the leakage contribution by directly subtracting the $E$-to-$B$ leakage band powers from the data without any fitting procedure.

In the likelihood analysis of the ABS and the TF pipelines, based on the constructed templates, the simulation tests have validated this template-fitting approach, showing the same performance in the $r$ estimate as that of the ‘backward correction’ approach. It should also be noted that during the TF procedures, the fitted foreground parameters in fact characterize the filtered foreground components.

The analytic blind separation (ABS) method is a novel, computationally efficient method for the blind separation of the CMB from foregrounds [72], which has been tested extensively for extracting the CMB temperature and polarization $E$- and $B$-mode power spectra from the simulated maps [73, 47]. The ABS estimator acts directly on measured multi-frequency cross band powers, which extracts CMB power spectra by projecting the data onto CMB subspace. In the presence of instrumental noise and a limited number of frequency channels, the ILC and ABS would yield slightly different results [74, 75, 47].

Assuming that the measured cross band power spectrum between the $i$- and $j$-th frequency channels in the multipole bin $\ell$ is given by $D^{\rm obs}_{ij}(\ell)$, the unique solution of the CMB in the noise-free case is achieved by the Sylvester’s determinant theorem and reads

as long as $M<N_{f}$, where $M\equiv rank\left(\mathcal{D}_{ij}^{\rm fore}(\ell)\right)$, and we choose $f_{i}=1$ for all channels so as to satisfy the emission law of the CMB in the units of thermodynamic temperature. ${\bf E}^{\mu}$ and $\lambda_{\mu}$ stand for the $\mu$-th eigenvector and associated eigenvalue of $\mathcal{D}^{\rm obs}_{ij}(\ell)$. We adopt the normalization condition for eigenvectors, ${\bf E}^{\mu}\cdot{\bf E}^{\nu}=\delta_{\mu\nu}$. This is similar to the idea of the harmonic-space ILC, in that the data are projected onto certain basis vectors to uniquely extract the CMB signal.

In the presence of instrumental noise, which is assumed to be an uncorrelated Gaussian distribution with zero mean and rms levels of $\sigma_{\mathcal{D},i}^{\rm noise}$ for the $i$-th frequency channel, $i=1,2\cdots N_{f}$, the noise-free solution of Eq. 5.1 then can be recast into

where

Here ${\bf\tilde{E}}^{\mu}$ and $\tilde{\lambda}_{\mu}$ are the $\mu$-th eigenvector and corresponding eigenvalue of $\tilde{D}^{\rm obs}_{ij}$, respectively. Note that in the above expressions we have chosen to not write the $\ell$ index explicitly. As seen, the ABS method thresholds the eigenvalues to keep only signal dominated modes, and we choose $\tilde{\lambda}_{\rm cut}=1$ here, as suggested in  [73]. The one important essence of the ABS is to introduce a positive shift parameter, $\mathcal{S}$, in Eq. (5.2), corresponding to that we artificially add a positive CMB signal to the data and then subtract it out of the solution. This is particularly important for low signal-to-noise regime and responsible for stabilizing the computation as long as it is large enough. From tests, we choose $\mathcal{S}=10\mu{\rm K}^{2}$ in our analysis.

The seven frequency channel DC1 simulations are analyzed on a sky patch with $f_{\rm sky}\sim 5\%$, by adopting the ABSapo mask. Details of the mask is further discussed in Sec. A and shown in Fig. 10. The plot in Fig. 4 presents the ABS derived CMB $BB$ band powers in the range of $\ell\in[20,200]$ with 9 $\ell$-bins, where the filtering effects on both band powers are not corrected since we use forwarding modelling of filtering effects in the ABS pipeline as mentioned in Sec. 4. Note that, the chi-squared value when fitting a model to data does not change if the model-predicted band powers and associated covariance have already incorporated filtering effects (refer to Paper III for details).

The estimation on tensor-to-scalar ratio $r$ through the ABS pipeline can be decomposed into two sequential steps. In short, we first recover the CMB band powers in each multipole bins through ABS estimator, and then we fit the recovered powers to the fiducial cosmological model to obtain the final estimate on $r$. The overall $BB$ band powers are modeled as a superposition of the tensor-only and the $TE$-only components, which reads

Here the tensor-only template $\bar{D}^{BB}_{r=0.01}$ consists merely of the filtered primordial $B$-mode signal, which is built from the mean of the band powers computed from 50 realization filtered CMB polarization maps with a fixed $r=0.01$, where the ensemble average of $BB$ spectrum from the $E$-to-$B$ leakage has been subtracted by using the leakage template, $\bar{D}^{BB}_{\rm TE~{}only}$. Such leakage template is constructed by the mean of $BB$ power spectra of 50 $TE$-only filtered CMB simulation maps, where the input CMB spectra are generated using the fiduical cosmological parameters but we explicitly null the $BB$ spectrum with setting $r=0$. All these maps used for building the templates are based on the ancillary data sets mentioned in Sec. 2.

In the ABS analysis, we employ a Gaussian likelihood (defined in Eq. (B.1) and discussed in Sec. B) to estimate the tensor-to-scale ratio, with using 9 $\ell$-bins in the multipole range $20\leq\ell\leq 200$. We stimate the covariance by using the ABS-derived band powers from 50 realization maps with the same generation process as for DC1 data, keeping the same foreground sky and the noise level as in DC1 data but setting the fiducial value of $r=0.01$ for $\Lambda$CDM model. These realizations are produced by combining the different component filtered maps from the ancillary data listed in Sec. 2. Since the CMB tensor contribution to the covariance is extremely minor compared with the noise and other contributions, the choice of fiducial value of $r$ has almost no changes on the estimated covariance, as long as it does not deviate greatly from the true value. As a result, various uncertainties contributed from the noise, CMB signal, $E$-to-$B$ leakage and TOD filtering effects as well as foreground residuals after the ABS cleaning have been appropriately included in the covariance.

An uniform prior on $r$ is assumed in the ABS analysis for $r\in[0,1]$. We sample the posterior distribution of the parameters using the Monte-Carlo Markov chain (MCMC) ensemble sampler dynesty [76, 77] with dynamic nested sampling scheme. As seen from Fig. 4, the recovered amplitudes of CMB $BB$ band power over all $\ell$-bins are almost consistent with the theoretical expectations within $1\sigma$ level, which can provide an accurate measurement of $r$ statistically. Note that, the band powers in the range of $80\leq\ell\leq 140$ become negative, which is caused by large noise fluctuations in these bins as the mean noise subtraction has been applied to the data. We also find an increased statistical uncertainty as $\ell$ increases, which is mainly due to the noise band power goes like $\ell^{2}$. By combining all 9 $\ell$-bin measurements from DC1 data, we find the likelihood peaks at near the input and the best-fit value is $r=0.019$, which is listed in Table 4. The ABS likelihood yields the constraint, $r<0.082$ at 95% confidence with $\sigma(r)=0.024$, fully consistent with the input $r=0.023$. The results for the null test will be shown in Table 5, where in the absence of the primordial $B$-modes, bounds on $r$ can be derived, $r<0.048$ at 68% CL and $r<0.069$ at 95% CL (with a best-fit in zero), respectively. The normalized posterior for $r$ is shown in Fig. 9 for both DC1 and null test cases.

The generalised least squares (GLS) pipeline is a foreground cleaning pipeline where we assume our multi-frequency maps to be a linear mixture of three emission components, CMB, synchrotron and thermal dust. The cleaned $B$-mode map is obtained by solving the linear system. This method has been discussed in detail in Paper III. We only summarise the pipeline here. For the GLS method, the observations, $\bm{d}(p)$, at pixel $p$, are modelled as:

Here $\bm{A}$ denotes the mixing matrix of the three emission components, $\bm{s}$ represents the templates for the three components. The mixing matrix for the foreground components is computed by assuming a power-law emission for the synchrotron, and the thermal dust is modeled as a modified black body. We adopt Planck priors to compute the mixing matrix. The GLS pipeline assumes this linear system as a model for the observed sky and solves the linear system of Eq. (5.5) with inverse noise variance ($N$) weights. The GLS weights are given by [78]:

and the template for component $c$ given by:

All DC1 bands are used in the analysis. The maps are pre-processed by masking for compact point sources and then interpolating in the masked spots. The UNPapo mask shown in Fig. 10 is used for further analysis of these maps. Refer to Sec. A and Paper III for point source pre-processing and mask choices. The GLS weights are computed and applied in harmonic space. The power spectrum is computed from the GLS cleaned CMB map with the PCL-TC estimator and debiased for noise and TOD filtering leakage. Finally we correct for power suppression due to the filtering by correcting with the transfer function. We show the power spectrum results for the GLS in Fig. 5. For the parameter estimation we will use the power spectra estimates in multipole range 40-200.

We will use a Gaussian likelihood to fit for the value of $r$ and other model parameters $\theta$. The Gaussian likelihood has the form given in Eq. (B.1). For the GLS pipeline, the theoretical power spectrum $C_{\ell}^{\rm th}$ is modeled as:

where we have introduced $A_{d}^{\rm res}$ as a nuisance parameter for the residual foreground contamination. For DC1, we use $A_{\rm lens}=0$ as these simulations are lensing free. The primary foreground contaminant for our experiment is thermal dust, so our foreground residuals only model the polarized dust emission. We adopt the model for residual dust power spectrum based on [79]:

Here $W^{\rm GLS}_{{\rm CMB},\nu,\ell}$ denotes the GLS weights for the CMB in frequency channel $\nu$. The dust frequency scaling in Eq. (5.9) is in brightness temperature units, while our observations are in thermodynamic temperature units. To convert the dust model to thermodynamic units we use a $C_{\nu}/U_{\nu}$ factor, where $C_{\nu}$ is the color correction factor to convert from dust spectrum to IRAS reference spectrum, $I_{\nu}=\nu^{-1}$, and $U_{\nu}$ is the unit conversion factor to convert thermodynamic temperature to brightness temperature with IRAS reference spectrum. Note that the $C_{\nu}/U_{\nu}$ factor is independent of the choice of reference spectrum. The color correction and unit conversion factors are calculated with bandpass integration following the prescription in [80]. The GLS weights are multiplied to compute the residual dust power in the GLS CMB map. We set $\alpha=-2.54$, $\beta_{d}=1.59$ and $T_{d}=19.6$ K, from [81, 79].

To compute the covariance matrix, $\bm{M}_{\rm fid}$, defined in Sec. B, we use 300 simulations of CMB and noise, with fixed foreground simulations that include thermal dust, synchrotron and point sources. We set $r=0.03$ in the fiducial $C_{\ell}$s. The noise simulations use WMAP 9 year noise variance maps and for HFI we use FFP10 simulations. The noise maps for AliCPT bands are simulated from the noise variance maps. Note that the foreground simulations used in these set of simulations are different from those used in DC1. These simulations are propagated through the GLS pipeline to obtain foreground cleaned $B$-mode maps. These maps are used to obtain the power spectrum and the covariance matrix. In this way we incorporate the variance contribution of the cosmic variance, noise and foreground residuals. In the covariance matrix, we only keep the the elements where the two bins are separated from each other by one or less.

The simulations used to compute the covariance matrix are not TOD filtered as TOD simulation and filtering are computationally expensive. This would imply we would have an additional contribution to the covariance matrix from the filtering. We have verified from CMB and noise only filtered simulations that additional variance due to the filtering and filtering corrections is an subdominant to the leading contributions to the covariance matrix computed here. The filtering effect can only lead to a small correction in the covariance of the first bin. By avoiding the filtering step we can perform 300 simulations which is sufficient for us to compute the covariance matrix with the primary contributions. However, in future analysis we would like to obtain the covariance matrix by propagating sufficient number of filtered maps through the GLS pipeline.

We fit the $B$-mode power spectrum from the GLS cleaned $B$-map in the multipole range 40-200. For the model used here we adopt an uniform prior on $r$ in the range $0<r<1$ and uniform prior on nuisance parameter $A_{d}^{\rm res}$ in the range $0<A_{d}^{\rm res}<10$. We use emcee [82] to estimate the posterior distribution of $r$. The posterior distribution for $r$ for the GLS pipeline for DC1 case is shown on the left in Fig. 6, while the null test is shown on the right of Fig. 6. The maximum a posteriori (MAP) best-fit value of $r=0.023^{+0.026}_{-0.013}$ while the minimum mean squared error (MMSE) estimate of $r=0.030_{-0.020}^{+0.019}$ for the GLS pipeline in the DC1 case. In the null test case the posterior distribution maximizes at $r=0$, the 95% upper limit on $r$ in this case is $<0.043$.

Internal linear combination (ILC) is a time tested component separation method that has been used since COBE to obtain foreground cleaned CMB maps from multi-frequency observations [83, 84, 85, 86, 87, 88, 89, 90, 91]. The ILC method calculates weights for each frequency channel by minimizing the total variance. The only assumption in the ILC method is the frequency dependence of the CMB (mixing vector), which is unity in all frequency bands if the maps are in thermodynamic temperature units. The traditional ILC methods are blind foreground cleaning methods. They make no assumptions about the frequency scaling of the foreground. If the mixing vector for the CMB is $A_{\nu,\text{ CMB}}$, the traditional ILC weights $W_{\nu}^{\rm ILC}$ only requires the constraint $\sum_{\nu}W_{\nu}^{\rm ILC}A_{\nu,\text{ CMB}}=1$, to preserve the CMB signal. The ILC weights minimize the contribution of both the foreground and noise to reduce the total variance. However, the variance due to foreground and noise in the $B$-mode case for DC1 is too large for traditional ILC methods to minimize them both, leading to large residual foreground signal.

The constrained ILC (cILC) method [92] is a semi-blind foreground cleaning technique. It uses priors on foreground information to set additional constraints that remove the average dust and synchrotron emission captured by the modeling of the data. For cILC we model the observed data as in Eq. (5.5), where we include models for the average thermal dust and synchrotron emissions. Then we place additional constraints: $\sum_{\nu}W_{\nu}^{\rm cILC}A_{\nu,\text{ dust/sync}}=0$, to remove the average dust and synchrotron emissions from the map. The cILC computes weights by minimizing the variance of the cleaned maps with total three constraints. The constraints on the cILC weights, $W^{\rm cILC}_{\nu}$ can be written together as:

where $A_{\nu,c}$ is the same mixing matrix from Eq. (5.5) and $e_{c}$ is 1 for CMB like for usual ILC method and 0 for the two foreground components modeled. This ensures that for very noisy data, the average dust and synchrotron captured by the model is removed with priority. The cILC weights obtained with these constraints are given by [92, 48]:

with $\bm{C}$ as the covariance matrix of the data. The cILC cleaned CMB map is then given by:

We pre-process and mask the data identically to that for the GLS pipeline discussed in Sec. 5.2. The cILC is implemented in harmonic space. The detailed implementation of the cILC for the AliCPT DC1 is discussed in Paper III. We obtain $B$-mode power spectrum from the cILC cleaned map. We estimate the noise bias from 50 filtered noise simulations by projecting the noise with the cILC weights. The leakage bias is computed from the cross-spectrum between the Planck-WMAP leakage template and the cILC cleaned $B$-map as discussed in Sec. 4. The power spectrum is debiased from the noise and TOD filtering leakage bias. Finally we correct the power suppression with transfer functions. The final power spectrum for cILC is shown in Fig. 5 on the right.

Similar to the case of GLS pipeline in Eq. (5.8), we maximize the Gaussian likelihood function of Eq. (B.1) with the theoretical model given by:

again with $A_{\rm lens}=0$. We fit the power spectra in the cILC case multipole range 40-200. The covariance matrix is computed for fiducial power spectrum with $r=0.03$ using 300 simulations prepared in the same way as described in Sec. 5.2. These maps are propagated through the cILC pipeline to obtain foreground cleaned $B$-mode maps. We then calculate the power spectra and their covariance matrix. As with the GLS pipeline, only elements for bins separated by one or less is retained in the covariance matrix for the likelihood estimation. We obtain the posterior distribution for $r$ by performing an MCMC in the parameter space using emcee. The normalized posterior distribution of $r$ is shown in Fig. 9 where the DC1 case is plotted on the left and the null test case is plotted on the right. The MAP best-fit value of $r=0.024^{+0.017}_{-0.015}$ while the MMSE estimate for $r=0.025\pm 0.016$ in the DC1 case. The posterior distribution for the null test case peaks at $r=0$, and the 95% upper limit is $<0.050$.

In the previous subsections, we perform foreground removal with multi-frequency observations to get the cleaned map (or directly the spectra in case of the ABS), then we estimate the $BB$ power spectra and finally fit the spectra to give the constraint on tensor scalar ratio ($r$). We can also start with our filtering corrected, multi-frequency, auto and cross power spectra and parameterize the foreground emission at the power spectra level, and then constrain the foreground and cosmological parameters simultaneously. This method is conceptually close to the multi-detector multi-component SMICA method [93, 94] widely used for CMB power spectrum estimation and component separation on various data sets [95, 96, 97], but explicitly uses parametric models of both the CMB and the foreground emissions. It has been used for BICEP/Keck and Planck joint analysis [46]. We implemented a multi-component multi-frequency spectra based likelihood (McMfL) with the parameterized angular power spectra of CMB, uncorrelated galactic dust, uncorrelated galactic synchrotron, and spatial correlation between dust and synchrotron. The total $BB$ auto- and cross-spectra between maps at frequencies $\nu_{1}$ and $\nu_{2}$ can be modeled as:

where $A_{d}$ ($A_{s}$) is the power amplitude for dust (synchrotron) at pivot frequency ($\nu_{\rm pivot}$) at 353 GHz (23 GHz) and pivot angular scale $\ell_{*}=80$. Here, $\Delta_{d}^{\prime}$ ($\Delta_{s}^{\prime}$) describes the decorrelation of the dust (synchrotron) pattern, which we set to 1. The spectral index in $\ell$ domain for dust (synchrotron), $\alpha_{d}$ ($\alpha_{s}$) can be written as,

where $\alpha_{d/s,0}$ is the scale independent dust/synchrotron spectral index, and $\alpha_{d/s,r}$ is the corresponding running index. The level of spatial correlation between dust and synchrotron is given by $\epsilon$. The factors $f_{d}^{\nu}$ (resp. $f_{s}^{\nu}$) scales the dust (resp. synchrotron) power from pivot frequency 353 GHz (23 GHz) to observed frequency band with bandpass function, $G(\nu)$, and are defined as follows,

with

Here, $\beta_{d}$ ($\beta_{s}$) is the dust (synchrotron) spectral index, $T_{d}$ is the dust temperature, which we set 19.6 K. For the CMB spectra $D_{\ell,{\rm CMB}}^{BB}$, we use camb [98] to generate the angular power spectra with tensor-to-scalar ratio ($r$). The parameter space for the likelihood involves {$r$, $A_{d}$, $A_{s}$, $\beta_{d}$, $\beta_{s}$, $\alpha_{d,0}$, $\alpha_{d,0}$, $\alpha_{d,r}$, $\alpha_{s,r}$, $\epsilon$}, while the other cosmological parameters fixed to the best-fit values from Planck 2018 [99].

We first calculate the bandpowers, ${D}^{BB,\text{ filt}}_{b,\nu_{1}\times\nu_{2}}$, using the PCL-SZ for each $\ell$-bin, for each frequency pair $\nu_{1}\times\nu_{2}$ from the observed maps. We set a bin width of $\Delta\ell=25$, and the bin window function is chosen as top-hat function. The BPSapo mask shown in Fig. 10 is used for all computations. These filtered bandpowers are then corrected for filtering by Eq. (4.8).

We adopt the Hamimeche & Lewis (HL) likelihood approximation [100] for the McMfL pipeline. The HL likelihood choice is justified as the degrees of freedom per bandpower bin is small due to the filtering and the small sky coverage. The HL likelihood is summarized in Sec. B, and the likelihood is given in Eq. (B.2). The HL likelihood is insensitive to the chosen fiducial model. The bandpass covariance matrix (BPCM), $\bm{M}_{\rm fid}$ only depends on the fiducial model, so we can pre-compute the BPCM. We use the 50 filtered CMB maps generated from our fiducial model ($\Lambda$CDM with tensor ($r=0.01$), which means the fiducial model is free from any foreground emissions) from the Ancillary Data to compute the fiducial power spectra, and the noise simulations in the Ancillary Data to compute the noise only power spectrum. These fiducial and noise bandpowers are then used to construct the BPCM following a semi-analytic method [101]. We also set the covariance between bandpowers separated by more than one $\ell$-bin to zero. The posterior distribution of model parameters is sampled using the modified CosmoMC³³3https://cosmologist.info/cosmomc/ [102]. Uniform prior is adopted for all parameters. The marginalized 2-D posterior distribution of $r$ as well as other foreground parameters is shown in Fig. 7 and Table 2. We obtained a best-fit value of $r=0.026\pm 0.019$. For the null test case, the constraint on $r$ is $r<0.042$ at $95\%$ confidence level.

The template fitting (TF) pipeline is an alternate pipeline to perform joint parametric fit of foreground, and filtering parameters along with the tensor-to-scalar ratio, similar to the McMfL pipeline. The TF and McMfL pipelines share the parametric philosophy of modeling the foreground and the CMB, and the availability of both provides an option to cross-check results obtained with the two pipelines. Despite their similarities they do have some key differences: 1) multiple variants of the foreground parameterization models are adopted in TF (see Paper III), 2) the TF pipeline uses Gaussian likelihood, whereas McMfL uses the HL likelihood, which leads to completely different estimates of the covariance, and 3) the TF uses the ‘forward modeling’ approach for filtering effects (see Sec. 4.4), and so no corrections are applied to the auto- and cross-spectra from multi-frequency maps during analysis. In the TF pipeline the leakage and suppression effects are already included in the CMB tensor-only template and the estimate of $r$ is derived from directly likelihood fitting the template to the DC1 data.

Here we report only the results of the TF pipeline using a more complicated foreground model (dubbed p16), as a representative example, different from that of McMfL in Eq. (5.14). The TF p16 model including 12 free parameters for the foreground, 3 parameters for the filtering effects on the foreground, together with one CMB parameter $r$ is designed to capture detailed features that we observe from simulations. The added parameters allow us to model the filtering-induced low-$\ell$ suppression for the foreground, the scale-dependence of the synchrotron-dust cross correlation, and the scale-dependence of $\beta_{d}$ and $\beta_{s}$. The scale-dependence of the synchrotron-dust cross correlation is modeled by replacing the constant $\epsilon$ in Eq. (5.14) with $\epsilon(\ell)=\epsilon_{2}\left({2}/{\ell}\right)^{\alpha_{\epsilon}}$, where the uniform priors are adopted on $\epsilon_{2}\in[-1,1]$ and $\alpha_{\epsilon}\in[0,2]$. Analogously, the scale-dependence of the spectral indices of $\beta_{d}$ and $\beta_{s}$ is modeled by replacing constant $\beta_{d}$ and $\beta_{s}$ with

respectively. The filtering-induced suppression factor in the foreground band-power amplitude is modeled as

with the uniform priors, $f_{0}\in[0,1]$, $\alpha_{F}\in[1,8]$ and $\ell_{F}\in[50,100]$. Here, this suppression factor acts only on the foreground model and not on the CMB, which is modeled through the CMB template defined in Eq. (5.4). Moreover, the filtering-induced suppression factor here is not determined from simulations. Instead we simultaneously estimate the foreground filtering parameters by sampling the total parameter space of the TF posterior. In this respect our approach to the filtering corrections is very different from that of the McMfL pipeline. Due to this difference in foreground modeling, the foreground parameters estimated from the TF, especially with respect to the dust and synchrotron amplitudes, are not guaranteed to be exactly the same as those from the McMfL.