© The Authors 2025

1. Introduction

Precise measurements of cosmic microwave background (CMB) polarization are particularly essential for modern cosmology and for exploring new physics beyond the ΛCDM model (Planck Collaboration VI 2020). The detection of the primordial B-mode polarization would reveal the mystery of the cosmic inflation era. The amplitude of the primordial B modes is commonly described by the ratio of the tensor power spectrum to the scalar power spectrum at some pivot scale, called the tensor-to-scalar ratio. Successive ground-based CMB experiments over the past decade have made great efforts to increase their sensitivities to enable a significant detection of the tensor-to-scalar ratio, such as BICEP/Keck Array (BICEP/Keck Collaboration 2021), the Atacama Cosmology Telescope (ACT) (Louis et al. 2017), the Cosmology Large Angular Scale Surveyor (CLASS) (Essinger-Hileman et al. 2014; Harrington et al. 2016), and the South Pole Telescope (SPT) (Sayre et al. 2020). The current best constraint on the tensor-to-scalar ratio is r < 0.032 (95% C.L.) and was obtained from the combination of Planck, BAO, and BICEP/Keck data (Tristram et al. 2022). Future CMB experiments such as the Ali CMB Polarization Telescope (AliCPT) (Li et al. 2019; Li et al. 2018), the Simons Observatory (Simons Observatory Collaboration 2019), the Q&U Bolometric Interferometer for Cosmology (QUBIC) (Mennella et al. 2019), LiteBIRD (Hazumi et al. 2019), and CMB-S4 (Abazajian et al. 2016, 2022; CMB-S4 Collaboration 2022) will continue to pursue the target of higher sensitivities on constraining r.

Precise measurements of the CMB polarization are challenging not only because it is weaker than the instrumental noise, but also because of complicated contaminants from astrophysics and instruments. Diffuse galactic foregrounds such as synchrotron and thermal dust emit strong polarized radiation at microwave frequency bands with a significant amplitude relative to the CMB polarization signal. Component separation methods are implemented on the multifrequency sky observations to remove the foregrounds (mostly) based on the fact that CMB and foregrounds have different spectral energy distributions (SED). Our imperfect knowledge of the actual optical response of the telescopes results in instrumental systematic effects, such as the unknown difference between the response of two polarized detectors that might transform a fraction of the CMB temperature modes into polarization modes. This is called the T → P leakage. The leakage must be mitigated especially because the temperature modes and the polarization modes peak at similar angular scales so that the leakage of the former has a higher magnitude than the signal of the latter. Additionally, the gravitational lensing effect where the matter field gravitationally deflects the CMB photons can covert part of E modes into B modes at intermediate and small scales, giving rise to lensing B modes that dominate the primordial B modes (Blanchard & Schneider 1987; Zaldarriaga & Seljak 1998). Fortunately, the lensing convergence can be reconstructed using either internal CMB anisotropies or external large-scale structure (LSS) tracers such as galaxy surveys, 21 cm observations, and measurements of the cosmic infrared background (CIB) (Hu 2001; Hu & Okamoto 2002; Okamoto & Hu 2003; Carron & Lewis 2017; Maniyar et al. 2021; Planck Collaboration IV 2020a).

We consider a typical configuration for a small-aperture ground-based CMB telescope. Specifically, we assumed a total of 7000 polarized detectors, evenly split between the 95 GHz and 150 GHz frequency bands. Each detector had a noise equivalent temperature (NET) of $350\mathrm{\mu }\text{K}\sqrt{s}$. The full width at half maximum (FWHM) beam sizes were assumed to be 19 arcminutes and 11 arcminutes for the 95 GHz and 150 GHz bands, respectively. This setup reflects the typical design of a third-generation ground-based CMB experiment (S3) (Abazajian et al. 2016). We focus on exploring the impact of the beam mismatch, which is a difference in the beam shape or beam center between the two detectors of a detector pair that leads to the T → P leakage. An analysis technique called the deprojection was exploited to filter this contamination out (BICEP2 Collaboration 2015; Sheehy 2019). We used the deprojection algorithm to handle the time-ordered data (TOD) with the beam systematics, and we evaluated the effects of undeprojected residual systematics on the follow-up steps of our whole pipeline, including component separation, lensing reconstruction, and the estimation of the tensor-to-scalar ratio. The residual systematic error was found to be negligible given the sensitivity of S3.

The paper is organized as follows. In Sect. 2 we first review the formalism of beam mismatches modeled by six coefficients. Then, we introduce our deprojection pipeline to mitigate the systematics induced by beam mismatches. Finally, we summarize the simulations used to assess the effectiveness of deprojection. In Sect. 3 we introduce the foreground cleaning pipelines including Needlet ILC (NILC) and constrained ILC (cILC), which we used to clean the T or E modes and B modes, respectively. We present the results of the mock data after deprojection and foreground cleaning, carry out a lensing reconstruction process on the foreground-cleaned maps, and constrain the tensor-to-scalar ratio r using the cleaned B-mode power spectrum in Sect. 4. Finally, we conclude in Sect. 5.

2. Beam systematics and deprojection

The beam is the response of an antenna to sky radiation as a function of angle. Basically, each pixel on the focal plane includes a pair of two orthogonally polarized antennas that observe the same direction on the sky. The time streams from detectors a and b (a pair) are summed to obtain the total intensity (temperature) measurement, and they are differenced to cancel the unpolarized component. The radiation polarization is measured in this way. The typical beams of S3 are circular Gaussians with FWHM = 19 arcmin for 95 GHz, and 11 arcmin for 150 GHz, where there is no systematics in the ideal case. However, when the beams of the two antennas are mismatched, which means that their instrumental responses are different, the temperature signal is not completely canceled and could induce a spurious polarization signal in the pair difference data of the detector pair¹. For the ith detector, the TOD stream reads

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\tau }_i(t)& =g_i{\displaystyle \int d}\nu {\displaystyle \int d}\mathrm{\Omega }\prime A_i(\nu )B_i[\widehat{\mathit{n}}(t)-\widehat{\mathit{n}}\prime ]\hfill \\ \hfill & [I{\prime }_{\nu }(\widehat{\mathit{n}}\prime )+Q{\prime }_{\nu }(\widehat{\mathit{n}}\prime )cos2{\psi }_i(t)+U{\prime }_{\nu }(\widehat{\mathit{n}}\prime )sin2{\psi }_i(t)].\hfill \end{array}\end{array}$$(1)

Of these, g_i is the gain of the detector, A_i(ν) is the response of this detector at different frequencies, and $B_i[\widehat{\mathit{n}}(t)-\widehat{\mathit{n}}\prime ]$ is the beam function, representing the spatial response of the detector. $[I,Q,U](\widehat{\mathit{n}},\nu )$ are the Stokes parameters of the cosmic microwave background radiation at frequency ν and position $\widehat{\mathit{n}}$. The terms $\widehat{\mathit{n}}(t)$ and ψ_i(t) are determined by the scan strategy. In an ideal case, $B(\widehat{\mathit{n}})$ and A_i(ν) are delta functions. The beam systematic errors are included in the following way:

1. The pointing error refers to the deviation between the pointing direction read by the instrument ${\widehat{\mathit{n}}}_{\text{Ins}}$ and the true pointing direction of the telescope $\widehat{\mathit{n}}$.

2. The band pass mismatch for a pair of detectors a and b means A_a(ν)≠A_b(ν), which can introduce false signals when taking the difference between the readings of the detectors.

3. The beam mismatch for a pair of detectors a and b similarly refers to B_a(r)≠B_b(r), which can also introduce false signals in the detector difference.

4. The gain mismatch arises when g_a ≠ g_b. Like the previous two points, this can introduce false signals in the detector difference.

5. The detector polarization angle calibration is typically detected by two detectors set to an angle of 90°. However, due to the precision of the detector manufacturing, the preferred polarization angle can deviate, that is, ψ_a(t)−ψ_b(t)≠90°, which can introduce systematic errors during the subsequent map-making process.

The first four errors above are all equivalent to an overall miscalibration of the detectors and a differential response between a pair of detectors at the same frequency, where we refer to the latter as beam mismatch in the following. A detector polarization angle miscalibration would introduce the E → B leakage, which is beyond the scope of this work. Given the high amplitude of the CMB temperature anisotropies relative to the polarization signal, the leakage from T to P must be carefully taken into account. The T → P leakage due to beam mismatches can be written as the convolution of the unpolarized signal and the differential beam within a detector pair,

$$\begin{array}{c}\hfill d_{T\to P}=T(\widehat{\mathit{n}})\ast [B_a(\widehat{\mathit{n}})-B_b(\widehat{\mathit{n}})]=T(\widehat{\mathit{n}})\ast B_{\delta }(\widehat{\mathit{n}}),\end{array}$$(2)

where $T(\widehat{\mathit{n}})$ refers to the temperature anisotropies, $B_a(\widehat{\mathit{n}})$ and $B_b(\widehat{\mathit{n}})$ are the beams of two detectors, and $B_{\delta }(\widehat{\mathit{n}})$ is the differential beam of the pair.

As in BICEP2 Collaboration (2015) (hereafter BKIII), we modeled $B_{\delta }(\widehat{\mathit{n}})$ as the difference of two elliptical Gaussian beams for simplicity. This can be parameterized by six parameters (or differential modes): the gain difference δg, the differences of the pointing center coordinates (with respect to the nominal beam center) δx and δy, the Gaussian beam-width difference δσ, and the differences of the plus- and cross-ellipticity δp and δc. The plus-ellipticity describes an ellipse whose major axis is oriented horizontally or vertically, and the cross-ellipticity describes an ellipse whose major axis is oriented ±45° with respect to the x-axis (see more details in BKIII, Appendix A). For each differential mode k = {g, x, y, σ, p, c}, δk = k_a − k_b is the difference of the beam parameter between two detectors. $B_{\delta }(\widehat{\mathit{n}})$ is the linear combination of the beams of the six modes, where each mode independently produces a pattern of the spurious polarized signal, also called a leakage template. For instance, since the beam of differential gain $B_{\delta g}(\widehat{\mathit{n}})$ is simply the scaled version of the circular Gaussian beam, its leakage is the scaled version of the beam-smoothed temperature $\stackrel{\sim }{T}(\widehat{\mathit{n}})$ according to Eq. (2). Under the assumption that the mismatch is small, except the third and higher orders, all the templates are the beam-smoothed map of $T(\widehat{\mathit{n}})$ and its first and second spatial derivatives (see Table 1) (Hu et al. 2003). The map-making procedure is a linear operation, meaning that the total leakage is the linear combination of the leakage templates of these differential modes. Therefore, we filtered the leakage out by fitting these templates to our data and then subtracting them.

The leakage templates we used are from the Planck 100 and 143 GHz reobserved TOD, which we used to clean 95 and 150 GHz, respectively. The reobserved templates were produced as follows:

1. We first generated CMB, foreground, and noise simulations of the Planck HFI bands as introduced below.

2. We coadded all components and reconvolved the coadded sky maps with the S3 nominal beams. For Planck 100 GHz, we divided a_ℓm by the 100 GHz beam function and then multiplied them by the S3 95 GHz Gaussian beam. (This can be done on the input maps since the Planck and S3 nominal beams are circularly symmetric.)

3. We computed the first and second derivatives of the Planck T maps².

4. We rescanned the T maps and their derivatives for each HFI band using the same pointing matrix as in S3.

For the realistic observations, we reobserved the PR4 (NPIPE) maps (Planck Collaboration Int. LVII 2020) using the S3 scan strategy. Our templates contained the noise component of the Planck HFI bands, which is not involved in the true T → P leakage, but we assumed that the noise is much lower than the temperature signal. We produced simulated S3 TOD from the CMB-plus-foreground coadded maps using Eq. (1), and we computed the pair-difference time stream of the TOD by differencing the a and b detector time streams within a pair: τ_diff(t) = (τ_a − τ_b)(t)/2. We finally fit the templates to the pair difference TOD and subtracted the fitted templates. Noise is not involved in the data since we expect that the noise is uncorrelated to the templates from Planck HFI maps, and the separation of noise ensured clearer analyses of the different components in the following foreground-cleaning pipeline.

The assumed S3 experiment is located in the northern hemisphere. In its first observing season, it will scan 17% of the sky where the thermal dust emission is clean, centering at RA = 170° and Dec = 40° in the celestial coordinate system. The polarized maps have a map depth of about 18 μK-arcmin at 95 GHz and 24 μK-arcmin at 150 GHz.

Our mock dataset consisted of 300 simulated sky maps at N_side = 1024 in seven frequency bands, including the 95 and 150 GHz dual bands, the Planck HFI 100, 143, 217, 353 GHz bands, and the WMAP K band. The additional Planck and WMAP channels were for the purpose of cleaning the foregrounds more thoroughly. Each of the sky maps contained CMB, foregrounds, and white noise. First, we drew the six ΛCDM cosmological parameters within one standard deviation of the Planck 2018 best-fit cosmological constraints (Planck Collaboration VI 2020), from which we generated the lensed CMB power spectrum by CAMB (Lewis et al. 2000). The CMB realizations were generated from the lensed power spectrum with the tensor-to-scalar ratio r = 0, using the synfast facility from lenspyx (Reinecke et al. 2023). The Galactic foregrounds containing thermal dust, synchrotron, free-free, spinning dust, and CO emissions were produced using the Planck Sky Model (PSM) package (Delabrouille et al. 2013), which is the same as the foregrounds used in previous papers (Dou et al. 2024; Han et al. 2023b). The thermal dust polarization maps were generated by scaling the Planck 2018 GNILC polarized-dust template (Planck Collaboration IV 2020) to different frequencies using a modified blackbody SED with the dust temperature and spectral indices adopted from the best fit of the Planck 2015 GNILC dust maps (Planck Collaboration Int. XLVIII 2016). The synchrotron polarization template was based on the Planck 2018 SMICA map and was scaled by a power-law SED with a fixed β_s of −3.08. We did not use the current best-fit values to avoid confirmation biases during foreground removal. For the 95 and 150 GHz bands alone, we generated the TOD simulations of CMB plus foregrounds following Eq. (1) and the S3 scan strategy, where we assumed the stochastic beam mismatch parameters whose mean and uncertainty are listed in Table 1. The values were chosen to be at the same level as the measured beam parameters of the BICEP2 instruments (BICEP2 Collaboration 2014; BICEP2/Keck Collaboration 2015). They measured the beam shape parameters for each detector using a chopped thermal source mounted on a mast, and calibrated the differential gains using the cross-correlation of temperature maps for individual detectors with Planck. We created white-noise simulations assuming that the realistic 1/f noise can be significantly removed by pair differencing and TOD polynomial filtering (BICEP2 Collaboration 2014). The white-noise simulations unrelated to systematics were generated from the noise covariance matrices of the S3 and WMAP-K bands (Bennett et al. 2013). For the Planck HFI bands, we adopted 300 FFP10 noise simulations from the Planck Legacy Archive³. The deprojection process was performed on the CMB-plus-foreground data assuming that the effects on the noise component are negligible.

3. Foreground-cleaning pipeline

We implemented the same foreground-cleaning methods as Han et al. (2023b) to qualify the deprojection performance on the foreground-cleaned maps. The T and E modes were cleaned by the NILC method, and the B modes were cleaned by the cILC method because the foreground dominates the B-mode signal.

The maps after foreground cleaning (NILC or cILC) consisted of four ingredients: the CMB signal, the residual foregrounds, the residual noise, and the systematic residual. For the purpose of examining the efficacy of the foreground cleaning methods, each component in the output map can theoretically be obtained by projecting ILC weights on the maps of the input component, given the linearity of ILC methods. However, the input CMB and foregrounds are mixed up during deprojection so that the systematic contamination on the CMB and foregrounds can never be separated afterward. Instead, we first subtracted the residual noise from the total residual (the foreground-cleaned map minus the input CMB map), which was the sum of the residual foregrounds and the systematic residual. We then projected the ILC weights on the deprojection residual maps at the 95 and 150 GHz bands to obtain the deprojection residual in the foreground-cleaned maps. Finally, we subtracted it from the sum to compute the residual foregrounds.

3.1. The NILC pipeline

The needlet internal linear combination (NILC) method (Delabrouille et al. 2009; Basak & Delabrouille 2012, 2013) is a widely used blind component-separation technique that linearly combines the raw data from different frequencies while minimizing the variance in the needlet space. Needlets are a special form of wavelets that permit a localization in harmonic and real (pixel) space. They are therefore suitable for component separation when the properties of the sky components vary significantly with the sky positions and scales. We adopted the NILC method for cleaning foregrounds in T- and E-mode maps, the details of which are described below.

The sky observations d_ν(p) at a frequency band ν and sky pixel p were first transformed into spherical harmonic coefficients using HEALPix⁴ (Gorski et al. 2005): d_ℓm^ν = ∫d_ν(p)Y_ℓm^*(p)dΩ. The raw maps with different beams were then reconvolved into the common resolution by $d_{\ell m}^{\nu ,\mathrm{out}}=d_{\ell m}^{\nu ,\mathrm{in}}{b}_{\ell }^{\mathrm{out}}/b_{\ell }^{\mathrm{in}}$, where b_ℓ is the beam transfer function. The common beam resolution (b_ℓ^out) was chosen as a Gaussian beam with FWHM = 11 arcmin, which is the resolution of the 150 GHz band. For the WMAP K band and 95 GHz, which have larger beam sizes than the common beam (b_ℓⁱⁿ < b_ℓ^out), we cut off the harmonic coefficients of these two bands with ℓ greater than 350 and 1200, respectively, to avoid an overamplification of the instrumental noise. The harmonic coefficients can be decomposed into a set of filtered maps,

$$\begin{array}{c}\hfill d_{\ell m}^{\nu ,j}=h_{\ell }^j{d}_{\ell m}^{\nu },\end{array}$$(3)

where the needlet bands h_ℓ^j satisfy ∑_j(h_ℓ^j)² = 1, defining a localization scheme in harmonic space. We adopted the cosine needlet bands as shown in Fig. 1 of Han et al. (2023b). The eight needlet bands peak at 15, 30, 60, 120, 300, 700, 1200, and 2000 with ℓ_max = 2000. For the WMAP K band, whose coefficients with ℓ > 350 were cut off, the needlets whose peaks (ℓ_peak) are greater than 350 would not be applied to the K-band data (and so is for 95 GHz band). For each frequency band, the filtered maps were then transformed back to real space, forming a set of needlet maps. The needlet coefficient of the jth needlet and the kth pixel after transformation is given by

$$\begin{array}{c}\hfill b_j^{\nu }({\widehat{p}}_{\mathit{jk}})=\sqrt{\frac{4\pi }{N_j}}{\displaystyle \sum _{\ell m}}{h}_{\ell }^j{d}_{\ell m}^{\nu }{Y}_{\ell m}({\widehat{p}}_{\mathit{jk}}),\end{array}$$(4)

where N_j denotes the number of pixels of the jth needlet map, and ${\widehat{p}}_{\mathit{jk}}$ denotes the pixel center of the kth pixel of the jth needlet map.

When the data were transformed to needlet space, we applied the internal linear combination method to the multichannel data. First we estimated the data covariance matrix across the frequencies by averaging the needlet coefficient product $b_j^{{\nu }_1}({\widehat{p}}_{\mathit{jk}})b_j^{{\nu }_2}({\widehat{p}}_{\mathit{jk}})$ over a disk of pixels centered at pixel k. The empirical data covariance between ν₁ and ν₂ is written as

$$\begin{array}{c}\hfill {\widehat{C}}_{\mathit{jk}}^{{\nu }_1\times {\nu }_2}=\frac{1}{n_k}{\displaystyle \sum _{k^{\prime }}}{w}_j(k,k^{\prime })b_j^{{\nu }_1}({\widehat{p}}_{j{k}^{\prime }})b_j^{{\nu }_2}({\widehat{p}}_{j{k}^{\prime }}),\end{array}$$(5)

where w_j(k, k′) are the weights that select the region around a pixel k to average out, and n_k is the number of selected pixels. The NILC solution is a linear combination of the needlet coefficients of all frequencies,

$$\begin{array}{c}\hfill b_j^{\mathrm{NILC}}({\widehat{p}}_{\mathit{jk}})={\displaystyle \sum _{\nu }}{w}_{\nu ,j}^{\mathrm{NILC}}({\widehat{p}}_{\mathit{jk}})b_j^{\nu }({\widehat{p}}_{\mathit{jk}}),\end{array}$$(6)

where the NILC weights are given by

$$\begin{array}{c}\hfill w_{\nu ,j}^{\mathrm{NILC}}({\widehat{p}}_{\mathit{jk}})={\left[\frac{{\widehat{\mathit{C}}}_{\mathit{jk}}^{-1}\mathit{a}}{{\mathit{a}}^t{\widehat{\mathit{C}}}_{\mathit{jk}}^{-1}\mathit{a}}\right]}_{\nu }.\end{array}$$(7)

Here, ${\widehat{\mathit{C}}}_{\mathit{jk}}$ is the n_ν × n_ν empirical data covariance matrix, where n_ν is the number of frequency channels, and a is the spectral response vector that characterizes the variation in the signal with frequencies. Because the CMB signal is frequency independent in thermodynamic temperature units, we set a as a unit column vector with a length of n_ν. The relation ${\displaystyle \sum _{\nu }}{a}_{\nu }{w}_{\nu ,j}^{\mathrm{NILC}}=1$ ensures the precise recovery of the CMB component.

Finally, we performed an inverse needlet transformation to obtain the NILC-cleaned map,

$$\begin{array}{c}\hfill {\widehat{s}}_{\ell m}^{\mathrm{NILC}}={\displaystyle \sum _{\mathit{jk}}}{b}_j^{\mathrm{NILC}}({\widehat{p}}_{\mathit{jk}})\sqrt{\frac{4\pi }{N_j}}{h}_{\ell }^j{Y}_{\ell m}({\widehat{p}}_{\mathit{jk}}).\end{array}$$(8)

The recovered CMB map in real space is given by ${\widehat{s}}_{\mathrm{NILC}}(p)={\displaystyle \sum _{\ell m}}{\widehat{s}}_{\ell m}^{\mathrm{NILC}}{Y}_{\ell m}(p)$. We used a mask that preserved the pixels with a noise variance of the 150 GHz band smaller than 20 μK pixels, whose sky fraction was about 10% (see Fig. 2).

3.2. The cILC pipeline

Following the ILC method, the constrained ILC (cILC) method (Remazeilles et al. 2011, 2021) was developed to further mitigate the residual foreground contaminants by adding constraints on the ILC weights. We used the cILC method in harmonic space to produce foreground-cleaned B-mode maps because the foregrounds are several orders of magnitude higher than the CMB signal for B modes, and the standard ILC method is therefore insufficient to reduce the residual foregrounds to a level of r = 0.01 primordial tensor modes under the sensitivities of S3.

We modeled the B-mode observation at a frequency ν and pixel p as the sum of n_c sky components and the thermal noise,

$$\begin{array}{c}\hfill d_{\nu }(p)={\displaystyle \sum _c}{A}_{\nu c}{s}_c(p)+n_{\nu }(p),\end{array}$$(9)

where s_c(p) is the astrophysical emission template of component c, and n_ν(p) is the noise. The n_ν × n_c mixing matrix A stores the spectral response information that describes the SED of component c at frequency ν. In our implementation, we only considered three dominant sky components: the lensed CMB, the polarized thermal dust, and the polarized synchrotron. The first column of the mixing matrix refers to the mixing vector of CMB, a. The second column corresponding to dust emissions was calculated assuming a modified blackbody SED with fixed T_dust = 19.6 K and β_dust = 1.59. The third column, that is, the mixing vector of the synchrotron, follows the power law SED with fixed β_sync = −3 in Rayleigh-Jeans brightness temperature units.

The cILC weights are given by

$$\begin{array}{c}\hfill {\mathit{w}}^{\mathrm{cILC}}=\mathit{e}{\left({\mathit{A}}^T{\widehat{\mathit{C}}}^{-1}\mathit{A}\right)}^{-1}{\mathit{A}}^T{\widehat{\mathit{C}}}^{-1},\end{array}$$(10)

which satisfy the constraints of nulling the foreground components (dust and synchrotron),

$$\begin{array}{c}\hfill {\displaystyle \sum _{\nu }}{w}_{\nu ,\ell }^{\mathrm{cILC}}{A}_{\nu ,c}=e_c,\end{array}$$(11)

where e is a row vector with length n_c, whose element e_c is one for the CMB and zero for the foreground components. The empirical data covariance matrix $\widehat{\mathit{C}}$ is computed by

$$\begin{array}{c}\hfill {\widehat{C}}_{\ell }^{{\nu }_1\times {\nu }_2}=\frac{1}{{\sum }_{{\ell }^{\prime }={\ell }_{min}}^{{\ell }_{max}}(2\ell +1)}{\displaystyle \sum _{{\ell }^{\prime }={\ell }_{min}}^{{\ell }_{max}}}{\displaystyle \sum _{m=-{\ell }^{\prime }}^{{\ell }^{\prime }}}{d}_{{\ell }^{\prime }m}^{{\nu }_1\ast }{d}_{{\ell }^{\prime }m}^{{\nu }_2},\end{array}$$(12)

where ℓ_min = min[0.6ℓ,ℓ−7], and ℓ_max = max[1.4ℓ,ℓ+7], and we excluded the center multipole (ℓ′≠ℓ) to mitigate the ILC bias. The input maps were reconvolved to a common beam of 11 arcmin. The analysis was performed with ℓ_max = 2000, while the harmonic coefficients of WMAP-K and 95 GHz bands were not involved as ℓ is greater than 350 and 1200, respectively. The cILC-cleaned CMB map was finally obtained by ${\widehat{s}}_{\ell m}^{\mathrm{cILC}}={\displaystyle \sum _{\nu }}{w}_{\nu ,\ell }^{\mathrm{cILC}}{d}_{\nu ,\ell m}$. To reduce the foreground contaminants, we used a smaller mask (f_sky ≈ 7%) to clean the B maps, which was produced by removing the pixels with a noise variance in the 150 GHz band larger than 10 μK-pixel and a declination above 65° (see the third row of Fig. 4).

4. Results

In this section, we detail the results of the fitted beam parameters from deprojection and the deprojected maps of the S3 dual bands. We then combine the S3 channels with four Planck HFI channels and WMAP K band in the foreground-cleaning pipeline to extract the CMB signal, and we reconstruct the lensing potential power spectrum from the foreground-cleaned maps. We finally evaluate the impacts of the systematic residual on the CMB power spectra, lensing potential, and the tensor-to-scalar ratio r.

4.1. Deprojection results

The assumed S3 telescope contains 1728 detector pairs on the focal plane for the 95 and 150 GHz bands, with a field-of-view width ∼30°. We considered one observing season (lasting about half a year), which was divided into 3183 scansets (each lasting roughly 30 minutes). A scanset covers tens of constant elevation scans (CES) that the telescope operates from one end to the other at a speed of ∼4° /s. Between the scansets, the telescope makes elevation nodding (elnod) scans to calibrate the relative gains of the detectors. We generated time streams from the simulated CMB plus foregrounds following the S3 scan trajectory. We derived the six leakage templates from the Planck simulations as described in Sect. 2. Although the input beam parameters are constant with time, we fit the time-ordered templates to the S3 TOD every 30 scansets for each detector pair to account for the time-varying beam systematics in realistic instruments. The histogram of the fitted beam mismatch parameters of 106 time chunks for one detector pair of a realization at 150 GHz is shown in Fig. 1. The scatter across time slices would be averaged in the map-making process, and the deviation of the mean from the true values (input parameters) due to noise fluctuations in templates leads to the undeprojected residual.

Fig. 1. Distribution of the recovered beam mismatch parameters of 106 time trunks for one detector pair of a noise-free simulation at 150 GHz. The dashed red lines represent the input beam parameter, and the dashed black lines indicate the mean over all time series. There is a bias on the differential plus-ellipticity of Δ ≈ 7.3 × 10−3 due to the cosmological TE correlation.

We note a systematic bias on the recovered differential plus-ellipticity δp at a level of Δ = δp_recov − δp_input ≈ 7.3 × 10⁻³. Due to the CMB-inherent TE correlation, the deprojection template of δp correlates with the polarized signal in the pair-difference time stream, resulting in a bias on the deprojection coefficient even when it is free of beam systematics (BICEP2 Collaboration 2015). This bias does not impair the filtering of T → P leakage, but leads to an additional filtering of cosmological E modes (the effects on the B modes are negligible). We corrected for the bias using the averaged deprojection coefficient estimated from systematics-free simulations, which is ⟨δp_recov⟩ = 7.2 × 10⁻³⁵. We also checked the strong correlation between the time-averaged recovered coefficients and the input coefficients for all detector pairs in Appendix A.

Finally, we subtracted the fitted templates from the pair difference data for all time trunks of detector pairs. We simplified the map-making process by averaging the time streams of detectors pointed at each map pixel and adding white-noise maps to obtain the observed sky maps.

4.2. Frequency maps

We performed the deprojection procedure on CMB-plus-foregrounds TOD with beam systematics of the S3 95 and 150 GHz bands and implemented the map-making to obtain the frequency maps. The 150 GHz B-mode maps before and after the deprojection pipeline are plotted in the first row of Fig. 2, where the left plot shows the input CMB plus foregrounds without beam systematics, the middle plot shows the CMB plus foregrounds with beam systematics, and the right plot shows that after deprojection. The beam systematics with the level considered in this work only substantially affect the B modes, which induce a spurious signal mainly at small scales. The deprojection residual maps shown in the second row of Fig. 2 illustrate that the deprojection can recover the input T-, E-, or B-mode signal. The deprojection has no effect on the T modes, but the noise fluctuation in the deprojection template may lead to a residual T-to-P leakage on the E and B modes. However, our results show that the residual leakage almost vanishes in deprojected maps.

Fig. 2. Maps of CMB plus foregrounds at S3 150 GHz. First row: Input B-mode CMB and foregrounds without beam systematics (left), CMB plus foregrounds with beam systematics (middle) and after deprojection (right). Second row: Difference between the deprojected map and the input map without beam systematics for the T (left), E (middle), and B (right) modes. The results of 95 GHz are similar to the 150 GHz results.

We used the PCL-TC estimator (as detailed in Appendix B) to estimate the power spectra of 150 GHz maps, as shown in Fig. 3. The band powers have a multipole bin size of 30 and ℓ_max = 700. Only in the BB power spectrum can the beam systematics bias the power spectrum of input CMB plus foregrounds significantly, while the power spectra of deprojected maps for all modes fit the input spectra well. The residual after deprojection (plotted as the dashed cyan curve) for BB is about 2∼3 orders of magnitude lower than the input CMB plus foreground power spectrum (solid black curves).

Fig. 3. TT, EE, TE, and BB power spectra of the 150 GHz CMB plus foregrounds map with beam systematics before (blue curves) and after deprojection (red curves) for one realization. The input CMB+FG. power spectrum without beam systematics is shown as a solid black curve for comparison with the deprojected one, and their difference, the residual after deprojection, is shown as a dashed cyan curve. The results of 95 GHz are similar to the 150 GHz results.

4.3. Foreground-cleaned maps

After propagating the maps with beam systematics through the deprojection pipeline, we used the NILC method to clean the foregrounds of T- and E-mode maps and used the cILC approach to clean the B-mode maps. The cleaned maps are plotted in Fig. 4. The left column shows the input CMB maps without adding beam systematics, the middle column shows the NILC-cleaned T and E maps from seven bands, and the right column shows the difference between the cleaned map and the input CMB. The residual maps are dominated by the noise fluctuation for the B modes, while the NILC maps for the T and E modes are dominated by signal.

Fig. 4. Maps of the input CMB (left column), the foreground-cleaned maps (middle column), and their difference, i.e., the residual maps (right column) for the T, E, and B modes. The first row shows T maps, the second row shows E maps, and the third row shows B maps.

For the B modes, we computed the systematic residual map by projecting the cILC weights on the 95 and 150 GHz deprojection residual maps plotted in Fig. 2. In this way, we separated the three components in the B-mode cILC residual map from each other, as shown in Fig. 5. The small-scale artifacts in the systematic contamination map are clearly visible. The residual on the boundaries would be removed through mask apodization in the power spectrum estimation.

Fig. 5. cILC residual B maps of the systematics (left), foreground (middle), and noise (right) components for one realization.

The power spectra of the NILC- and cILC-reconstructed maps are shown in Fig. 6. The orange curves refer to the residual noise power spectrum obtained from projecting the NILC or cILC weights on the input noise maps. The green curves, as mentioned in Sect. 3, represent the sum of the foreground and beam systematic residuals. For the B modes, the two contaminant components are separately plotted as red and purple curves. The foreground residual dominates the systematic residual at large scales (ℓ ≲ 300). The deprojection residual is negligible compared to the CMB B-mode power spectrum at all scales. Our foreground cleaning methods are sufficient to recover the CMB signal, and the residual contaminants remain much lower than the reconstructed CMB power spectrum for the T and E modes and lower than the noise fluctuation for the B modes.

Fig. 6. Power spectra including NILC-cleaned TT (left top), NILC-cleaned EE (right top), NILC-cleaned TE (left bottom), and cILC-cleaned BB (right bottom) power spectra with ℓmax = 1500. The input CMB (solid black curves), the reconstructed CMB (blue bars), and the residual noise (orange curves) power spectra are plotted in each panel. The green curves represent the sum of the systematic and foreground residuals. For the BB spectra, the deprojection and foreground residuals are separately shown as red and purple curves, respectively. The power spectra are binned over every 30 multipoles except for the BB spectra, where a logarithm binning at ℓ > 100 was adopted for clarity.

4.4. Lensing reconstruction

In addition to the primordial B-mode signal, S3 can also detect the lensing signal with a moderate significance because the noise performance is excellent. In what follows, we evaluate the systematic effects on the lensing reconstruction by using the simulated 4 modules*yr data. We fed the 301 simulation sets mentioned above through the lensing pipelines (Carron & Lewis 2017; Liu et al. 2022a; Han et al. 2023a) and tested their performance.

In the presence of lensing-induced correlations among different multipole moments, the off-diagonal terms of the covariance matrix of the CMB fields become nonzero. This allowed us to calculate the quadratic estimator (Hu & Okamoto 2002; Okamoto & Hu 2003) using pairs of filtered maps in their quadratic form (Carron & Lewis 2017; Maniyar et al. 2021), as described in Eqs. (3.4) to (3.8) of Liu et al. (2022b). In this approach, one map in the pair is Wiener filtered, while the other is inverse-variance filtered. To mitigate biases introduced by random disconnected noise, masking effects, and foreground residuals, we employed a Monte Carlo simulation method to estimate these contributions and performed a mean-field subtraction by averaging over the quadratic estimators, using 60 simulation sets.

Normalization was conducted in two stages: We first assessed the averaged noise level for the entire mock dataset and computed the normalization factor analytically, assuming an effective isotropic noise level. Then we corrected for normalization bias arising from noise inhomogeneities through numerical simulations. The raw power spectrum of the lensing is simply

$$\begin{array}{c}\hfill {\widehat{C}}_L^{\widehat{\varphi }\widehat{\varphi }}=\frac{1}{(2L+1)f_{\mathrm{sky}}}{\displaystyle \sum _{M=-L}^L}{\widehat{\varphi }}_{\mathit{LM}}{\widehat{\varphi }}_{\mathit{LM}}^{\ast },\end{array}$$(13)

where ${\widehat{\varphi }}_{\mathit{LM}}$ denotes the harmonic transformation of the reconstructed lensing potential. However, the outcome of the quadratic estimator contains not only the lensing potential signal, but also the Gaussian reconstruction noise sourced by the CMB and instrumental noise (N₀ bias) and the nonprimary couplings of the connected four-point function (Kesden et al. 2003) (N₁ bias). We calculated the realization-dependent N₀ bias (RDN0) using 240 simulation sets, while the N₁ bias was calculated analytically. After subtracting these biases, we obtained the final estimated power spectrum,

$$\begin{array}{c}\hfill {\widehat{C}}_L^{\varphi \varphi }={\widehat{C}}_L^{\widehat{\varphi }\widehat{\varphi }}-\mathrm{\Delta }{C}_L^{\widehat{\varphi }\widehat{\varphi }}{{|}_{\mathrm{RDN}0}-\mathrm{\Delta }{C}_L^{\widehat{\varphi }\widehat{\varphi }}|}_{\mathrm{N}1}.\end{array}$$(14)

In Fig. 7 we show the reconstructed lensing potential power spectrum from polarization data, which combines the EE, EB, and BB estimators. The discrepancy in the figure between the reconstructed data point and the theoretical prediction (black curve) originates in the Gaussian uncertainty of the single realization that represents the “real” data. Thus, it is expected that some of the data points may deviate by up to 1σ from the theoretical curve considering that each ℓ bin is independent.

Fig. 7. Reconstructed lensing potential power spectrum from the polarization estimator (green points). The colored boxes denote the 1σ error regions calculated from 240 simulation sets.

To evaluate the outcome of our lensing reconstruction, we estimated the signal-to-noise ratio (S/N) via the Fisher-matrix method,

$$\begin{array}{c}\hfill \mathrm{S}/\mathrm{N}=\sqrt{{\displaystyle \sum _{\ell ,{\ell }^{\prime }}}{C}_{\ell }{\mathbb{C}}_{\ell {\ell }^{\prime }}^{-1}{C}_{{\ell }^{\prime }}},\end{array}$$(15)

and we used the theoretical instead of the reconstruced C_ℓ in the numerator to ensure a stable prediction. ℂ_ℓℓ′ is the covariance matrix obtained from 240 simulation sets, which reads

$$\begin{array}{c}\hfill {\mathbb{C}}_{\ell {\ell }^{\prime }}=\frac{1}{N-1}{\displaystyle \sum _{n=1}^{N=240}}[({\widehat{C}}_{\ell }^{\varphi \varphi }-{\overline{C}}_{\ell }^{\varphi \varphi })\times ({\widehat{C}}_{{\ell }^{\prime }}^{\varphi \varphi }-{\overline{C}}_{{\ell }^{\prime }}^{\varphi \varphi })],\end{array}$$(16)

where ${\overline{C}}_{\ell }^{\varphi \varphi }$ is the averaged lensing-potential power spectrum based on the simulation sets. The final S/N obtained from the multipole range of ℓ ∈ (20, 340) is S/N ≈4.2, and this number is consistent with Han et al. (2023b), where we reported an S/N ≈4.5 with the same foreground-removing and lensing- reconstruction pipeline using a database without systematics. We therefore conclude that the systematic residual does not contaminate the lensing reconstruction strongly.

4.5. Tensor-to-scalar ratio bias

In this section, we estimate the bias on the tensor-to-scalar ratio r due to the beam systematic residual by comparing the estimated r bias of the simulations with (after deprojection) and without the beam systematics.

We sampled the posterior distribution of the tensor-to-scalar ratio r from the averaged cILC-cleaned BB power spectra by an MCMC analysis. We assumed the Gaussian likelihood of r given the cleaned band powers ${\widehat{C}}_{{\ell }_b}$ as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill -2ln\mathcal{L}(r)={\displaystyle \sum _{{\ell }_b{\ell }_{b^{\prime }}}}& [{\widehat{C}}_{{\ell }_b}-r{C}_{{\ell }_b}^{r=1}-C_{{\ell }_b}^{\mathrm{lens}}]{\left[M_{\mathrm{fid}}^{-1}\right]}_{{\ell }_b{\ell }_{b^{\prime }}}\hfill \\ \hfill & [{\widehat{C}}_{{\ell }_{b^{\prime }}}-r{C}_{{\ell }_{b^{\prime }}}^{r=1}-C_{{\ell }_b}^{\mathrm{lens}}],\hfill \end{array}\end{array}$$(17)

where ℓ_b, ℓ_b′ are indices for multipole bins, rC_ℓb^r = 1 + C_ℓb^lens refers to the theoretical binned BB power spectrum with a multipole range of ℓ ∈ [40, 200] including 5 ℓ_b, $C_{\ell }^{\mathrm{lens}}$ is the B-mode spectrum due to gravitational lensing, $C_{\ell }^{r=1}$ is the B-mode spectrum sourced by tensor perturbations with r = 1, and ${[{\mathit{M}}_{\mathrm{fid}}]}_{{\ell }_b{\ell }_{b\prime }}=⟨({\widehat{C}}_{\mathrm{fid},{\ell }_b}-⟨{\widehat{C}}_{\mathrm{fid},{\ell }_b}⟩)({\widehat{C}}_{\mathrm{fid},{\ell }_{b\prime }}-⟨{\widehat{C}}_{\mathrm{fid},{\ell }_{b\prime }}⟩)⟩$ is the fiducial covariance matrix computed from 300 systematics-free simulations. We derived the posterior distribution function assuming a uniform prior of r ∈ [0, 1]. We constructed an MCMC chain consisting of 10 000 samples that satisfied the posterior distribution using the emcee (Foreman-Mackey et al. 2013) Python package. We obtained the 1σ confidence interval (C.I.) and the 95% upper limit of r for the deprojection- and systematics-free cases.

The distribution of r is shown in Fig. 8. The resulting 1σ C.I. for the systematics-free case is r = 0.019 ± 0.013 and the 95% C.I. is r < 0.043, while for the deprojection case, the results are nearly the same. We note that in the selected scales ℓ ∈ [40, 200] where the B-mode spectrum is sensitive to tensor perturbations, the deprojection residual is negligible compared to the tensor power spectrum $C_{\ell }^{r=1}$ and foreground contaminants (see Fig. 6). Therefore, the beam systematics imprint negligible effects on the measurements of r.

Fig. 8. Distribution of the tensor-to-scalar ratio for the cases with and without beam systematics. The dashed vertical lines represent the 95% upper bound, which is 0.043 for the two cases.

5. Conclusions

We forecast the systematic effects of beam mismatches on the forthcoming S3 observations. We adopted the deprojection technique in the map-making procedure to mitigate the T → P leakage induced by the random beam mismatches added to the S3 mock TOD, and we then assessed the performance of the deprojection on the foreground cleaning, lensing reconstruction, and r-estimation. First, we compared the CMB-plus-foregrounds after deprojection and the input CMB-plus-foregrounds without instrumental effects in terms of the TEB maps and the angular power spectra. The CMB and foregrounds are mixed up in deprojection and thus cannot be separated from each other. We found that the residual leakage can be ignored compared to the input CMB-plus-foregrounds.

Then, we imposed the foreground-cleaning methods, NILC for T and E modes and cILC for B modes, on the 300 simulations of seven frequency bands including the S3 bands, the Planck HFI bands, and the WMAP K band, where the S3 mock data are deprojected. The power spectrum of the biases consisting of the residual T → P leakage and foregrounds was shown to be negligible compared to the noise uncertainties for the TT, EE, TE, or BB power. We evaluated the effect of beam systematics on the lensing reconstruction. Our lensing-potential power spectrum reconstructed from the polarization estimator has a similar signal-to-noise ratio as that obtained in previous analyses without systematics, which demonstrates that the lensing reconstruction is almost unaffected by residual systematics. Finally, we set constraints on the tensor-to-scalar ratio r using the cleaned BB spectrum, where we obtained similar results with and without beam systematics because the deprojection technique mitigates the systematics to a negligible level. Our results justify the use of a deprojection to filter the T → P leakage due to beam mismatches in S3 experiments, and the residual leakage can be ignored in the following data analysis pipeline. The effects of far sidelobes are not involved in this work since they are related to the specific shielding system of the telescope (BICEP2/Keck Collaboration 2015) and the calibration or modeling of the realistic beam maps (Gallardo et al. 2018). The undeprojected residual due to the differential beams of the sidelobes will be calibrated in future works.

Acknowledgments

We thank Zirui Zhang, Siyu Li and Jing Jin for helpful discussions. This work is supported by the National Key R&D Program of China Grant No. 2021YFC2203102, 2020YFC2201603, NSFC No. 12325301 and 12273035. Jiakang Han thanks Stefano Camera for various support for this project. Some of the results in this paper have been derived using the pymaster, emcee, healpy, and HEALPix packages.

References

Appendix A: Recovered beam parameters versus the input ones

Figure A.1 shows the comparison between the fitted beam mismatch parameters in deprojection and the input values for 1728 detector pairs of a realization at 150 GHz. The red solid curves refer to y = x, and the red dashed curve was vertically offset by the estimated mean bias on δp (7.2 × 10⁻³). The Pearson correlation coefficient R approaching 1 for all parameters indicates a strong correlation between the deprojection coefficients and the input parameters, demonstrating the precision of deprojection. The standard deviation of the bias (Δ) on each coefficient is significantly lower than the median measurement uncertainty for individual detectors by BICEP2/Keck Array (see BICEP2/Keck Collaboration 2015, Table 3), but is too optimistic since the noise was not added to the data.

Fig. A.1. Scatter plot of the recovered beam mismatch parameters averaged over time trunks (x-axis) vs. the input values (y-axis) for all 1728 detector pairs of a simulation at 150 GHz. The red solid curves have a slope of 1 and a y-intercept of 0. The red dashed curve has been offset vertically by the mean bias on δp estimated from systematics-free simulations. The correlation coefficient R and the standard deviation of the bias on the deprojection coefficient σ(Δ) are exhibited.

Appendix B: PCL-TC estimator

The E → B mixing is introduced when the partial-sky Q and U maps are converted to the E and B modes in harmonic space (Zhao & Baskaran 2010; Ghosh et al. 2021). We used the template cleaning method (Liu et al. 2019a,b) to get rid of the resulting E → B leakage. To estimate the full-sky B-mode power spectrum from foreground-cleaned maps, we adopted the pseudo-C_ℓ method to account for mode-coupling induced by masking and multipole binning using the NaMASTER (Alonso et al. 2019)⁶ package. The whole process is called the “PCL-TC” estimation.

The template cleaning method is a pixel-based technique to remove E → B leakage due to partial-sky effects. The steps are as follows:

1. First convert the masked (Q, U) maps into (a_ℓm^E, a_ℓm^B) using the map2alm function of healpy (Zonca et al. 2019). Transform a_ℓm^B into the B map being contaminated by the E → B leakage, B_c.

2. Convert the E-only modes (a_ℓm^E, 0) back into maps (Q_E, U_E) using the alm2map function of healpy. Multiply them by the mask and convert them into $(a_{\ell m}^{E\prime },a_{\ell m}^{B\prime })$ using map2alm. The B map derived from $a_{\ell m}^{B\prime }$ is taken as our E → B leakage template, B_t.

3. Linearly fit the masked leakage template B_t to the masked contaminated map B_c in the pixel space. Finally subtract the fitted leakage to obtain the pure B-mode map: $B_{\mathrm{pure}}=B_c-\alpha B_t$ where α is the linear coefficient.

Deprojection, as any other filtering operation on Q and U time streams, can produce a small amount of E → B mixing. The E → B leakage due to deprojection can also be corrected for by taking advantage of the deprojected E-mode only data to get a B-mode leakage template (see Ghosh et al. 2022, Sect. 4.2).

After the foreground cleaning, we corrected for the common beam (11 arcmin) and pixel window functions in the a_ℓm’s. Finally we computed the full-sky bandpowers with a bin size of 30 using the decouple_cell module of pymaster.

All Tables

All Figures

	Fig. 1. Distribution of the recovered beam mismatch parameters of 106 time trunks for one detector pair of a noise-free simulation at 150 GHz. The dashed red lines represent the input beam parameter, and the dashed black lines indicate the mean over all time series. There is a bias on the differential plus-ellipticity of Δ ≈ 7.3 × 10−3 due to the cosmological TE correlation.
In the text

Fig. 2. Maps of CMB plus foregrounds at S3 150 GHz. First row: Input B-mode CMB and foregrounds without beam systematics (left), CMB plus foregrounds with beam systematics (middle) and after deprojection (right). Second row: Difference between the deprojected map and the input map without beam systematics for the T (left), E (middle), and B (right) modes. The results of 95 GHz are similar to the 150 GHz results.

In the text

Fig. 3. TT, EE, TE, and BB power spectra of the 150 GHz CMB plus foregrounds map with beam systematics before (blue curves) and after deprojection (red curves) for one realization. The input CMB+FG. power spectrum without beam systematics is shown as a solid black curve for comparison with the deprojected one, and their difference, the residual after deprojection, is shown as a dashed cyan curve. The results of 95 GHz are similar to the 150 GHz results.

In the text

	Fig. 4. Maps of the input CMB (left column), the foreground-cleaned maps (middle column), and their difference, i.e., the residual maps (right column) for the T, E, and B modes. The first row shows T maps, the second row shows E maps, and the third row shows B maps.
In the text

	Fig. 5. cILC residual B maps of the systematics (left), foreground (middle), and noise (right) components for one realization.
In the text

Fig. 6. Power spectra including NILC-cleaned TT (left top), NILC-cleaned EE (right top), NILC-cleaned TE (left bottom), and cILC-cleaned BB (right bottom) power spectra with ℓmax = 1500. The input CMB (solid black curves), the reconstructed CMB (blue bars), and the residual noise (orange curves) power spectra are plotted in each panel. The green curves represent the sum of the systematic and foreground residuals. For the BB spectra, the deprojection and foreground residuals are separately shown as red and purple curves, respectively. The power spectra are binned over every 30 multipoles except for the BB spectra, where a logarithm binning at ℓ > 100 was adopted for clarity.

In the text

	Fig. 7. Reconstructed lensing potential power spectrum from the polarization estimator (green points). The colored boxes denote the 1σ error regions calculated from 240 simulation sets.
In the text

	Fig. 8. Distribution of the tensor-to-scalar ratio for the cases with and without beam systematics. The dashed vertical lines represent the 95% upper bound, which is 0.043 for the two cases.
In the text

Fig. A.1. Scatter plot of the recovered beam mismatch parameters averaged over time trunks (x-axis) vs. the input values (y-axis) for all 1728 detector pairs of a simulation at 150 GHz. The red solid curves have a slope of 1 and a y-intercept of 0. The red dashed curve has been offset vertically by the mean bias on δp estimated from systematics-free simulations. The correlation coefficient R and the standard deviation of the bias on the deprojection coefficient σ(Δ) are exhibited.

In the text