Determination of miscalibrated polarization angles from observed CMB and foreground $EB$ power spectra: Application to partial-sky observation

We study a strategy to determine miscalibrated polarization angles of cosmic microwave background (CMB) experiments using the observed $EB$ polarization power spectra of CMB and Galactic foreground emission. We apply the methodology of Minami et al.(2019) developed for full-sky observations to ground-based experiments such as Simons Observatory. We take into account the $E$-to-$B$ leakage and $\ell$-to-$\ell$ covariance due to partial sky coverage using the public code NaMaster. We show that our method yields an unbiased estimate of miscalibrated angles. Our method also enables simultaneous determination of miscalibrated angles and the intrinsic $EB$ power spectrum of polarized dust emission when the latter is proportional to $\sqrt{C_\ell^{EE}C_\ell^{BB}}$ and $C_\ell^{BB}$ is proportional to $C_\ell^{EE}$.


Introduction
B-mode polarization of the cosmic microwave background (CMB) is a probe of primordial gravitational waves generated during inflation [2]. To detect this, we need to remove spurious B-mode signals such as those from gravitational lensing [3], patchy reionization [4,5], Eto-B leakage due to partial sky coverage [6,7], and miscalibration of polarization angles of detectors [8]. In this paper, we focus on estimation of the miscalibrated polarization angles from CMB experiments.
The miscalibrated polarization angles, α, yield not only a spurious power spectrum of the B-mode polarization, but also a cross power spectrum of E-and B-mode polarization [8].
The traditional method to estimate α is to relate the observed EB power spectrum to the theoretical model of the difference between the E-and B-mode power spectra of CMB [9]. This approach is thus not applicable in foreground-dominated frequency bands. In Ref. [1], we solved this problem by relating the observed EB power spectrum to the observed difference between the E-and B-mode power spectra in the sky. However, the analysis given there was limited to full-sky observations by a satellite experiment such as LiteBIRD [10]. In this paper, we apply the methodology developed in Ref. [1] to ground-based observations with a partial sky coverage. While we use Simons Observatory [11] as an example, our method can be applied to any other experiments including Simons Array [12], BICEP Array [13], and CMB-S4 [14].
Partial sky breaks orthogonality of spin-2 spherical harmonics, resulting in the leakage of E-mode to B-mode and the correlation between different multipoles. In this paper we take these effects into account by using the public code NaMaster [15].
The rest of the paper is organized as follows. In Sect. 2 we review the methodology of Ref. [1] and describe the extension to partial-sky data. In Sect. 3 we describe the way to generate sky simulations for validation of our method. In Sect. 4 we validate our method using sky simulations, and present the main results. In Sect. 5 we extend our method to determine the miscalibration angle and the intrinsic EB correlation of thermal dust emission simultaneously. We conclude in Sect. 6.

Review of the methodology
We review the methodology of Ref. [1] developed for full-sky observations. In this section, all spherical harmonics coefficients and power spectra are calculated from full-sky data without a mask.
In Ref. [1], we took into account both the miscalibration angle, α, and the "cosmic birefringence" [16][17][18][19][20] angle, β, which only rotates CMB polarization. When we consider noise ("N"), foreground ("fg"), and CMB ("CMB") components, spherical harmonics coefficients of the observed ('o') E-and B-mode polarization are related to the intrinsic ones by In this paper, we use the notations that all spherical harmonics coefficients and power spectra have been multiplied by the appropriate beam smoothing function.
When we define the power spectra with ensemble average as The current data show no evidence for non-zero EB correlation from the foreground emission [21,22] and from the CMB [23]. Ignoring C EB,fg and C EB,CMB , we can determine α and β simultaneously.
It is however possible that future experiments may find a non-zero EB correlation. We revisit the term C EB,fg in Sect. 5 for the case that we do not ignore C EB,fg .
To determine α and β simultaneously, we used a log-likelihood function given by where C EE,CMB,th b 2 and C BB,CMB,th b 2 are the best-fitting ΛCDM theoretical power spectra multiplied by the beam transfer functions, b 2 . As for the variance in the denominator, we used an approximated variance given by [1] Var To determine α and β, we minimized Eq. 4 with respect to α and β, given C EB,o , C EE,o − C BB,o , C EE,CMB,th , and C BB,CMB,th .

Extension to partial-sky data
In this paper, we shall extend the methodology of Ref. [1] to partial-sky data. We focus on determination of α and ignore the cosmic birefringence angle, β, in Eq. 4. We obtain However, Eq. (6) ignores the -to-bin covariance due to partial sky coverage. We thus write where C(α) is a one-dimensional array of C EB,o − tan(4α) 2 C EE,o − C BB,o with from min to max , and C(α) is the corresponding covariance matrix. The explicit form of C(α) is The size of C(α) is ( max − min + 1)-by-( max − min + 1). To determine α, we minimize Eq. (7) with respect to α, given C EB,o and C EE,o − C BB,o .

Estimation of unbiased power spectra and their covariance matrices
Throughout this paper, we use a public code "NaMaster" [15] to estimate unbiased C from partial-sky data. Since we can describe the E-to-B leakage and -to-covariance due to a partial sky coverage by the spherical harmonics transform of the sky mask, we can correct them and estimate the unbiased C from partial-sky data [15,24,25].
We also use NaMaster to estimate the covariance matrices of all pairs of C EE,o , C BB,o , and C EB,o . NaMaster estimates the covariance matrices under certain approximations [26,27]. We validate the estimated covariance matrices against Monte-Carlo (MC) simulations in Appendix A.

Sky simulations
To validate our methodology, we use the "PySM" package [28] to produce realistic simulations of the microwave sky, with an experimental specification similar to the Large Aperture Telescope (LAT) of Simons Observatory (SO) [11] (see Table 1). We include two polarized Galactic foreground emission models: "s1" for synchrotron model and "d1" for dust emission model, as described in Ref. [28]. The noise is assumed to be white with standard deviation given by σ N = (π/10800)(w given in the "Polarization Sensitivity" row of Table 1. A CMB map is generated from the power spectra calculated by CAMB [30] using the latest Planck 2018 cosmological parameters for "TT,TE,EE+lowE+lensing" [31]: Ω b h 2 = 0.02237, Ω c h 2 = 0.1200, h = 0.6736, τ = 0.0544, A s = 2.100 × 10 −9 , and n s = 0.9649. The CMB map includes the lensed B-mode but does not include the primordial B-mode, i.e., the tensor-to-scalar ratio, r, is zero.
We choose a frequency band of 280 GHz for this study because it is the highest frequency of the SO LAT and is most affected by the polarized dust emission. The generated map is masked by LR42 galactic mask of Ref.
[22], because its sky coverage, f sky = 0.42, is similar to that of the SO LAT, f sky = 0.4 [11]. The edges of the mask are apodized with a 5 deg FWHM Gaussian to reduce the mask induced E-to-B leakage.
In this paper we set the input miscalibration angle α (α in ) to 0.33 deg, which corresponds to the calibration uncertainty of the Crab nebula (Tau A) [32].
In the following sections, we will show results both with a single realization and with MC simulations. With a single realization, we will test whether we can recover the input angles within the estimated uncertainties. With MC simulations, we will test whether the estimator is unbiased and whether the estimated uncertainty is valid.
For the MC simulations, we have generated many different realizations assuming Gaussian fluctuations with zero mean. The Gaussian foreground realizations have been generated from the power spectra of the PySM map, while the CMB and noise realizations have been generated with the same way described above. We have used the synfast function of HEALPix to generate Gaussian realizations from the input power spectra.

Results
First, we report the results of α determination from a simpler (but sub-optimal) loglikelihood given by Eq. (6). As the off-diagonal elements of the covariance matrix have negative values (see Appendix A.2), neglecting off-diagonal elements gives a conservative  We find that the sub-optimal method recovers correctly α in , whereas it overestimates uncertainty. Thus, we can use our simple (but sub-optimal) method as a conservative method.
As a reference, we also show the results from the traditional method [9] in the right two columns of Table 2 using the log-likelihood function of Here the traditional method relates the observed EB power spectrum to the theoretical model of the difference between the E-and B-mode power spectra of CMB. As for Var C EB,o , we use the diagonal elements of the covariance matrix, Cov C EB,o , C EB,o , estimated by NaMaster. Since this log-likelihood also ignores the -to-bin covariance, it is also a sub-optimal method.
There are two notable points: (1) With a single realization, the estimated uncertainty from our sub-optimal method is smaller than that from the traditional sub-optimal method because our sub-optimal method uses correlation between the observed power spectra; (2) The estimator from our sub-optimal method is unbiased while that from the sub-optimal traditional method is biased. This is because our method automatically includes foreground components, while the traditional method ignores them. Table 2: Recovered α and their 1σ uncertainties against α in = 0.33 deg with two methods: our simpler (but sub-optimal) method (6) and the traditional sub-optimal method (9). In both methods, we show results of both a single realization and MC simulations.
α (deg) from Eq. (7) Single realization MC simulations 0.302 ± 0.085 0.330 ± 0.086 Next, we report the results from the optimal log-likelihood as given in Eq. (7). We set the bin size to ∆ = 8, which leads to 600-bins from min = 200 to max = 5000. We show the recovered values of α and their 1σ uncertainty in Table 3. In a similar way with the sub-optimal methods, we derive results of both a single realization and MC simulations with 1000 realizations.
We find that our optimal method recovers correctly α in with a single realization and find that the estimator is unbiased. We also find that the the covariance matrix in Eq. (7) is valid because the uncertainty estimated by a single realization is consistent with the uncertainty derived from MC simulations.
The uncertainty on α is smaller than the uncertainty determined only with the diagonal elements of covariance matrix (6). This is because the off-diagonal covariance-matrix elements have negative values as described in Appendix A.2.

Intrinsic EB correlation of Galactic foreground emission
In this section, we show that we can determine the miscalibration angle, α, and the intrinsic EB power spectrum of polarized dust emission (if any) simultaneously. We suggested this idea in Ref. [1]; The key is to put thermal dust EB correlation parameters into a new rotation angle, γ. In a simple model [33], we can relate C EB,f g to C EE,f g and C BB,f g by where f c is a correlation coefficient. Since the dust foreground EE and BB power spectra have a similar shape in [22], we may write the BB spectrum as where ξ is, e.g., equal to A BB /A EE ≈ 0.5 [22]. When f c satisfies 0 ≤ 2f c √ ξ/(1 − ξ) ≤ 1, we can parametrize the foreground EB correlation with a rotation angle, γ, as sin(4γ)/2 = Thus we can put thermal dust EB correlation parameters into a rotation angle, γ.
Since dust foreground is rotated by α + γ and CMB is rotated by α, we can use Eq. (4) to determine α and γ simultaneously, if we replace α → α + γ and β → −γ. For partial-sky data, we have where C(α, γ) is a one dimensional array of C EB,o − tan(4α+4γ) 2 cos(4α+4γ) C EE,CMB,th b 2 − C BB,CMB,th b 2 from min to max . To test our idea, we prepare a dust foreground which follows Eq. (10) and Eq. (11) and ignore synchrotron foreground. As for the input dust foreground power spectrum, we use C EE,f g calculated from a full-sky map of PySM. Then we set C BB,f g = ξC EE,f g and C EB,f g = f c C EE,f g C BB,f g . Using these power spectra, we realize a dust foreground map assuming Gaussian fluctuations. As for the input parameters for the polarized dust emission, we set ξ = 0.5 as suggested in Ref. [22], and f c = 0.01 which is consistent with observed data [33]. From analytical calculation, we can interpret the values of foreground parameters as γ = 0.41 deg.
CMB and noise maps are generated in the same way as in Sect. 3. With these inputs, we determine α and γ simultaneously, using the log-likelihood function of Eq. (13); We use the same binning criteria as in Sect. 4. We show the recovered α, γ, and their 1σ uncertainties in Table 4. As in Sect. 4, we show results of both a single realization and MC simulations with 100 realizations.
We find that input values of α and γ are correctly recovered with a single realization and the estimators are unbiased. We also find that covariance matrix in Eq. (13) is Table 4: Recovered α and γ and their 1σ uncertainty against α in = 0.33 deg and γ in = 0.41 (f c = 0.01 and ξ = 0.5) from Eq. (13). We show results of both a single realization and MC simulations.
Single realization MC simulations Because of the simultaneous determination of two parameters, the uncertainties become larger compared to α only estimation (Table 3).

Discussion and conclusion
In this paper, we have studied a strategy to determine miscalibrated polarization angles of CMB experiments using the observed EB power spectra of CMB and Galactic foreground emission. We have extended the methodology of Ref. [1] developed for full-sky observations to partial-sky observations. We have corrected the E-to-B leakage and -to-bin covariance due to partial sky coverage using the framework of NaMaster [15].
Applying our method to simulated maps of CMB, realistic foreground emission [28], and instrumental noise with beam smearing similar to the SO LAT [11], we have found that the method correctly recovers the input values of α and estimates its reasonable uncertainty.
We have also developed a method to estimate the intrinsic EB power spectrum of polarized dust emission. This is possible when C EB,fg is proportional to C EE,fg C BB,fg and C BB,fg is proportional to C EE,fg . Then we can interpret the EB power spectrum of dust emission as if it were generated by a polarization rotation angle, γ. We can determine α and γ simultaneously, as α affects CMB, while α + γ affects polarized dust emission. We have found that this method correctly recovers the input values of α and γ and estimates their reasonable uncertainties.
Though we have assumed homogeneous foreground EB cross correlation parameters, f c and ξ, we can simply extend our method to deal with spatially varying foreground by replacing them with multipole-dependent ones, i.e., f c ( ) and ξ( ). These can be interpreted as γ( ) and simultaneously determined with α.
While we have applied our method to the SO, we can apply the method to other groundbased observations with a partial sky coverage. This new framework allows us to determine the miscalibration angle in foreground-dominant frequency bands, and to detect a non-zero EB power spectrum of polarized dust emission.
We have ignored some detailed experimental parameters. One of them is 1/f noise which contaminates low multipoles. Because the miscalibration angles are mainly determined by the information at high multipoles where cosmic variance is small, the effect of 1/f would be small. We leave studies with more detailed experimental parameters to future work.
A Validation of the covariance matrix of NaMaster In this section, we show two tests to validate the covariance matrix calculated by NaMaster [15]: (1) Validation of diagonal elements; (2) Validation of off-diagonal elements with binning.
We validate the covariance matrix calculated analytically with the NaMaster framework by comparing with that from Monte-Carlo (MC) simulations. Because it needs enormous data volume to simulate off-diagonal covariance matrix elements with MC simulation, we validate diagonal elements of the covariance matrix with N = 100 samples of unbinned power spectrum, and validate off-diagonal ones with N = 1000 MC samples of binned power spectrum.
In the MC simulations, we generated N full-sky maps assuming Gaussian distributions with zero mean as we have described in Sect. 3.

A.1 Validation of diagonal elements
We compare diagonal elements of the covariance matrix C(α) (8) calculated analytically with the NaMaster framework and those calculated with N = 100 MC-simulated samples. We show the comparison in Figure A1. We find that analytically-calculated covariancematrix elements are consistent with MC-simulated covariance-matrix elements at all multipoles.

A.2 Validation of off-diagonal elements
We compare off-diagonal elements of the covariance matrix C(α) (8) calculated analytically with the NaMaster framework with those from 1000 MC-simulated samples. In the calculation of both analytical covariance-matrix and MC-simulated power spectra, the covariance is binned with a bin size of ∆ = 8 with range of ( min , max ) = (200, 5000).
We show the comparison of the off-diagonal elements of covariance matrices in Figure A2. We find that analytically-calculated covariance-matrix elements are consistent with MCsimulated covariance-matrix elements at all distances between bins, ∆b. And, we can see that |∆b| ∼ 1 has a negative covariance. Therefore, ignoring the off-diagonal covariance-matrix elements (as in Eq. (6) and Eq. (9)) overestimates uncertainties on the angles. analytical 50th bin 1000 sample 50th bin analytical 100th bin 1000 sample 100th bin analytical 150th bin 1000 sample 150th bin Fig. A2: The covariance-matrix elements against the distance between bins, ∆b. The black markers are derived from 1,000 MC simulations, while the red markers are analytically calculated by NaMaster.