# Planck constraints on Cosmic Birefringence and its cross-correlation with the CMB

G Zagatti, M Bortolami, A Gruppuso, P Natoli, L Pagano, G Fabbian

► To cite this version:

G Zagatti, M Bortolami, A Gruppuso, P Natoli, L Pagano, et al. Planck constraints on Cosmic Birefringence and its cross-correlation with the CMB. Journal of Cosmology and Astroparticle Physics, 2024, 05, pp.034. 10.1088/1475-7516/2024/05/034. [hal-04446865]

HAL Id: hal-04446865

https://hal.science/hal-04446865v1

Submitted on 17 Jul 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

Planck constraints on cosmic birefringence and its cross-correlation with the CMB

To cite this article: G. Zagatti et al JCAP05(2024)034

View the [article online](#) for updates and enhancements.

You may also like

- [Planck PR4 (NPIPE) map-space cosmic birefringence](#)
Raelyn M. Sullivan, Arefe Abghari, Patricia Diego-Palazuelos et al.

- [Cosmic birefringence: cross-spectra and cross-bispectra with CMB anisotropies](#)
Alessandro Greco, Nicola Bartolo and Alessandro Gruppuso

- [Constraints on cosmic birefringence using E-mode polarisation](#)
Arefe Abghari, Raelyn M. Sullivan, Lukas T. Hergt et al.
RECEIVED: January 29, 2024

REVISED: March 26, 2024

ACCEPTED: March 26, 2024

PUBLISHED: May 3, 2024

# Planck constraints on cosmic birefringence and its cross-correlation with the CMB

G. Zagatti$^{a,}$b M. Bortolami,$^{a,b}$ A. Gruppuso,$^{c,d,a}$ P. Natoli,$^{a,b}$ L. Pagano$^{a,b,e}$
and G. Fabbian$^{f,g}$

E-mail: giorgia.zagatti@unife.it, marco.bortolami@unife.it,
alessandro.gruppuso@inaf.it, paolo.natoli@unife.it,
luca.pagano@unife.it, fabbiang@cardiff.ac.uk

**ABSTRACT:** Cosmic birefringence is the in-vacuo, frequency independent rotation of the polarization plane of linearly polarized radiation, induced by a parity-violating term in the electromagnetic Lagrangian. We implement a harmonic estimator for the birefringence field that only relies on the CMB E to B mode cross-correlation, thus suppressing the effect of cosmic variance from the temperature field. We derive constraints from *Planck* public releases 3 and 4, revealing a cosmic birefringence power spectrum consistent with zero at about $2\sigma$ up to multipole $L = 1500$. Moreover, we find that the cross-correlations of cosmic birefringence with the CMB T-, E- and B-fields are also well compatible with null. The latter two cross-correlations are provided here for the first time up to $L = 1500$.

**KEYWORDS:** CMBR polarisation, CMBR theory

ARXIV ePRINT: 2401.11973

<figure><img src="image_2.png" /></figure>

<footnote>$^a$Dipartimento di Fisica e Scienze della Terra, Università degli Studi di Ferrara,
via Saragat 1, I-44122 Ferrara, Italy</footnote>

<footnote>$^b$Istituto Nazionale di Fisica Nucleare, Sezione di Ferrara,
via Saragat 1, I-44122 Ferrara, Italy</footnote>

<footnote>$^c$Istituto Nazionale di Astrofisica — Osservatorio di Astrofisica e Scienza dello Spazio di Bologna,
via Gobetti 101, I-40129 Bologna, Italy</footnote>

<footnote>$^d$Istituto Nazionale di Fisica Nucleare, Sezione di Bologna,
viale Berti Pichat 6/2, I-40127 Bologna, Italy</footnote>

<footnote>$^e$Institut d'Astrophysique Spatiale, CNRS, Univ. Paris-Sud, Université Paris-Saclay,
Bât. 121, 91405 Orsay cedex, France</footnote>

<footnote>$^f$School of Physics and Astronomy, Cardiff University,
The Parade, Cardiff, Wales CF24 3AA, U.K.</footnote>

<footnote>$^g$Center for Computational Astrophysics, Flatiron Institute,
162 Fifth Avenue, New York, NY 10010, U.S.A.</footnote>
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>Harmonic estimator</td><td>3</td></tr><tr><td></td><td>2.1 EB-estimator</td><td>5</td></tr><tr><td>3</td><td>Angular power spectrum and de-biasing procedure</td><td>8</td></tr><tr><td>4</td><td>Data set and simulations</td><td>9</td></tr><tr><td>5</td><td>Results</td><td>10</td></tr><tr><td></td><td>5.1 Full mission</td><td>11</td></tr><tr><td></td><td>5.2 Data splits</td><td>15</td></tr><tr><td></td><td>5.3 Consistency checks</td><td>17</td></tr><tr><td></td><td>5.4 CB cross-correlation</td><td>18</td></tr><tr><td>6</td><td>Sensitivities of future experiments</td><td>19</td></tr><tr><td>7</td><td>Conclusions</td><td>21</td></tr><tr><td>A</td><td>Construction of the estimator</td><td>23</td></tr><tr><td>B</td><td>Validation tests</td><td>25</td></tr><tr><td></td><td>B.1 Validation with a rotation signal</td><td>25</td></tr><tr><td></td><td>B.2 Validation without a rotation signal</td><td>26</td></tr><tr><td>C</td><td>Jack-knife tests</td><td>26</td></tr></table>

## 1 Introduction

The Cosmic Microwave Background (CMB), a radiation that marks the transition from an opaque to a transparent Universe, is a key observable for investigating cosmological physics. For decades, CMB experiments [1–4] mainly focused on the temperature field of the CMB radiation, whose information was extracted almost completely by the *Planck* satellite [5] up to few arc-minute scale. Much of the current and future experimental effort is devoted to measuring the polarization part of the CMB radiation [6–11] which is linearly polarized at the 1–10% level due to Thomson scattering. The polarization field is usually decomposed into two linear polarization modes: the E-mode, which is parity-even and couples to both scalar and tensor perturbations, and the B-mode, which is parity-odd and exclusively couples to tensor perturbations [12]. In the standard scenario, the polarization pattern of the CMB is described by the Maxwell’s electromagnetism, which preserves parity symmetry. In such a case the electromagnetic Lagrangian is described by

$$ \mathcal{L}_{em}^{SM} = -\frac{1}{4} F_{\mu\nu} F^{\mu\nu}, \qquad (1.1) $$
where $F_{\mu\nu}$ is the electromagnetic tensor that contains the electric and magnetic fields. Since eq. (1.1) satisfies parity symmetry, it is possible to show that the CMB cross-correlations TB and EB are expected to be zero.

However, there are recent claims of deviations from null of the latter cross-correlations [13–16], which are consistent with the Cosmic Birefringence (CB) effect [17], i.e. the rotation of linear polarisation plane of photons during propagation. Specifically these papers, which are based on Planck data and make use of a new technique [18–20] able to disentangle the instrumental polarization angle from the CB effect,<fnref n="1" /> hint at a detection of a CB angle $\beta \sim 0.3^\circ$ at the level of 2.5–3$\sigma$. The latter CB angle is also called isotropic, meaning that it does not depend on the direction of observation. There is also an anisotropic CB effect, which instead depends on the direction of observations, that is currently found to be well compatible with null [23–27]. See also [28] for a review of the CB effect from CMB observations.

If not due to systematic effects of instrumental or astrophysical origin [29–31], these analyses hinting at an isotropic CB effect, suggest the need to extend the electromagnetic sector of the standard model $\mathcal{L}_{em}^{SM}$ with a parity-violating term $\mathcal{L}_{CS}$

$$ \mathcal{L}_{em} = \mathcal{L}_{em}^{SM} + \mathcal{L}_{CS} = -\frac{1}{4} F_{\mu\nu} F^{\mu\nu} - \frac{\lambda}{4f} \phi F_{\mu\nu} \tilde{F}^{\mu\nu}, \qquad (1.2) $$

known as a Chern-Simons term [32].<fnref n="2" /> In eq. (1.2), $\lambda/f$ is a coupling with the dimension of the inverse of an energy, $\phi$ is a new scalar (or pseudo-scalar) field and $\tilde{F}^{\mu\nu}$ is the dual electromagnetic tensor. With such an extension, it is possible to describe an isotropic CB effect, and consequently a CMB EB cross-correlation compatible with observations, when $\phi$ is taken to be homogeneous [35, 36]. Fluctuations of $\phi$ around its homogeneous part, are instead able to produce anisotropic CB [37–39]. Hence, the CB effect can be seen as a tracer of the existence of a new cosmological field $\phi$ (typically referred as an axion) acting as dark matter or dark energy [40–43] which might also play a role in alleviating the Hubble tension, see e.g. [44, 45] and reference therein. See also [46–50].

Other works, e.g. [36, 51], have constrained the axion's parameters through the CMB EB<fnref n="3" /> power spectrum or the isotropic CB effect. However, these models can be put on additional tests considering also the anisotropic CB and the cross-correlations between the anisotropic CB and the CMB field, see e.g. [52]. Therefore it is essential to provide updated constraints on the latter observables. For this reason in this work, we focus on the anisotropic CB effect implementing a harmonic estimator based on the approach presented in [53] and [25]. The aim is to apply this estimator to the most recent Planck data, namely PR3 [5] and PR4 (also known as NPIPE) [54], and to provide new constraints on the CB power spectrum and the CB cross-correlation with the CMB fields. In particular the cross-correlation between anisotropic CB and the CMB E- and B-fields are given here for the first time up to $L = 1500$ (previously in [27], considering a different technique, those were provided only at very low multipoles).

<footnote n="1">¹Otherwise the uncertainty of the instrumental polarization angle has to be assessed independently in the total error budget, see e.g. [21, 22].</footnote>

<footnote n="2">²For other extensions, see e.g. [33, 34].</footnote>

<footnote n="3">³Current detections of the isotropic CB effect are based on the CMB cross-correlation EB and not TB, since the latter has a lower signal to noise ratio, at least a factor of 2 for Planck [22].</footnote>
This paper is organized as follows. In sections 2 and 3, we present the methodology employed to estimate the CB power spectrum, describing the structure of the estimator for the spherical harmonic coefficients and the de-biasing procedure necessary to obtain the final estimate. In section 4, we describe the CMB data and simulation sets utilized in this study. In section 5, we present the results of applying our pipeline to Planck CMB polarization maps and the cross-correlations of the CB and the CMB temperature and polarization fields. Furthermore, in section 6, we forecast the sensitivity of forthcoming CMB experiments, such as the LiteBIRD satellite, the Simons Observatory, and CMB-S4 to the EB cross-correlation. We conclude in section 7. The full calculations leading to the final expression of the estimator can be found in appendix A, and validation tests for our pipeline are presented in appendix B.

## 2 Harmonic estimator

This section introduces the impact of CB on CMB observations, in order to provide the structure of the harmonic estimator used in this work.

The primary effect is that the observed CMB polarization field carries the information of the rotation field [38]. Consequently, we do not observe the primordial E-modes and B-modes, denoted as $a_{\ell m}^X$, instead we observe their sum with the rotation-induced modes, $\delta a_{\ell m}^X$. Here $X = E, B$ for E- and B-modes respectively:

$$a_{\ell m}^{E,tot} = a_{\ell m}^{E} + \delta a_{\ell m}^{E}, \qquad (2.1)$$

$$a_{\ell m}^{B,tot} = a_{\ell m}^{B} + \delta a_{\ell m}^{B} \approx \delta a_{\ell m}^{B}, \qquad (2.2)$$

where the second equivalence in equation (2.2) holds since we are assuming that the B-modes generated on the last scattering surface are null.<fnref n="4" /> The expressions for the rotation-induced E-modes and B-modes, as derived in [38], are given as:

$$\delta a_{\ell m}^{B} = 2 \sum_{LM} \sum_{\ell' m'} \alpha_{LM} a_{\ell' m'}^{E} \xi_{\ell m \ell' m'}^{LM} H_{\ell \ell'}^{L}, \qquad (2.3)$$

$$\delta a_{\ell m}^{E} = 2i \sum_{LM} \sum_{\ell' m'} \alpha_{LM} a_{\ell' m'}^{E} \xi_{\ell m \ell' m'}^{LM} H_{\ell \ell'}^{L}, \qquad (2.4)$$

where equation (2.3) is different from zero for $\ell + \ell' + L$ even, equation (2.4) is different from zero for $\ell + \ell' + L$ odd, $\alpha_{LM}$ are the spherical harmonic coefficients of the CB field, and $\xi_{\ell m \ell' m'}^{LM}$ and $H_{\ell \ell'}^{L}$ are defined in terms of Wigner-3j symbols as follows:

$$\xi_{\ell m \ell' m'}^{LM} = (-1)^m \sqrt{\frac{(2\ell + 1)(2L + 1)(2\ell' + 1)}{4\pi}} \begin{pmatrix} \ell & L & \ell' \\ -m & M & m' \end{pmatrix}, \qquad (2.5)$$

$$H_{\ell \ell'}^{L} = \begin{pmatrix} \ell & L & \ell' \\ 2 & 0 & -2 \end{pmatrix}. \qquad (2.6)$$

The rotation-induced modes generate correlations between $\ell - \ell'$ pairs with $\ell \neq \ell'$, leading to, as anticipated in section 1, parity-violating cross-correlations. These non-standard cross-correlations caused by CB can be characterized by the following general structure:

$$\langle a_{\ell m}^{X,tot} a_{\ell' m'}^{X',tot,*} \rangle = 2 \sum_{LM} \alpha_{LM} Z_{\ell \ell'}^{XX'} \xi_{\ell m \ell' m'}^{LM} H_{\ell \ell'}^{L}, \qquad (2.7)$$

<footnote n="4">⁴In this study we present the calculations and the results based on this assumption. However, the entire analysis can be generalized to the case of non-zero B-modes on the last scattering surface.</footnote>
where $X = \{T, E, B\}$ and $Z_{\ell\ell'}^{XX'}$ contains the information about the primordial spectra before the rotation (see table 1 of [53]). Another way to write equation (2.7) is in terms of the rotational invariants (i.e. quantities independent of $m$) [38], $D_{\ell\ell'}^{LM,XX'}$:

$$ \langle a_{\ell m}^{X, tot} a_{\ell' m'}^{X', tot, *} \rangle = \sum_{LM} D_{\ell\ell'}^{LM,XX'} \xi_{\ell m \ell' m'}^{LM}, \qquad (2.8) $$

where,

$$ D_{\ell\ell'}^{LM,XX'} = 2\alpha_{LM}Z_{\ell\ell'}^{XX'}H_{\ell\ell'}^{L}. \qquad (2.9) $$

The starting point to get an expression for the estimator are the rotational invariants, indeed. Their definition in equation (2.9) refers to the primordial signal, thus when moving to observations we have to account for the window function, approximated with a Gaussian symmetric beam, $W_\ell$, so that [53]:

$$ D_{\ell\ell'}^{LM,XX',map} = D_{\ell\ell'}^{LM,XX'} W_{\ell} W_{\ell'} = 2\alpha_{LM} Z_{\ell\ell'}^{XX'} H_{\ell\ell'}^{L} W_{\ell} W_{\ell'}, \qquad (2.10) $$

where we use the superscript “map” to denote quantities recovered from CMB maps.

Having the analytical definition of the observed rotational invariants, we can now provide for two different expressions of the associated estimators. In this work we use an over-hat symbol to indicate estimated quantities. The first expression is directly derived from equation (2.10):

$$ \hat{D}_{\ell\ell'}^{LM,XX',map} = 2\hat{\alpha}_{LM} Z_{\ell\ell'}^{XX'} H_{\ell\ell'}^{L} W_{\ell} W_{\ell'}, \qquad (2.11) $$

and the other is the inverse variance weighting average estimator from [55]:

$$ \hat{D}_{\ell\ell'}^{LM,XX',map} = (G_{\ell\ell'})^{-1} \sum_{mm'} a_{\ell m}^{X, map} a_{\ell' m'}^{X', map, *} \xi_{\ell m \ell' m'}^{LM}, \qquad (2.12) $$

where $G_{\ell\ell'}$ is defined as:

$$ G_{\ell\ell'} = \sum_{mm'} (\xi_{\ell m \ell' m'}^{LM})^2 = \frac{(2\ell + 1)(2\ell' + 1)}{4\pi}. \qquad (2.13) $$

The first expression of the estimates for the $\alpha_{LM}$ coefficients for a fixed $\ell-\ell'$ pair is obtained inverting equation (2.11) and substituting the expression of the estimator for the rotational invariants of equation (2.12):

$$ (\bar{\alpha}_{LM})_{\ell\ell'}^{XX'} = \frac{\hat{D}_{\ell\ell'}^{LM,XX',map}}{F_{\ell\ell'}^{L,XX'}} = \frac{(G_{\ell\ell'})^{-1} \sum_{mm'} a_{\ell m}^{X, map} a_{\ell' m'}^{X', map, *} \xi_{\ell m \ell' m'}^{LM}}{F_{\ell\ell'}^{L,XX'}}, \qquad (2.14) $$

where we defined $F_{\ell\ell'}^{L,XX'} = 2Z_{\ell\ell'}^{XX'} H_{\ell\ell'}^{L} W_{\ell} W_{\ell'}$ [38]. Note that we are not using the over-hat symbol since, before ending up with the final estimates of the $\alpha_{LM}$ coefficients, we have to encode for a de-biasing procedure.

The final expression of the harmonic estimator before the de-biasing procedure at the level of the spherical harmonic coefficients, has been obtained applying the definition of the inverse variance weighting average to the $(\bar{\alpha}_{LM})_{\ell\ell'}^{XX'}$ defined in equation (2.14), ending up with:

$$ \bar{\alpha}_{LM}^{XX'} = \frac{\sum_{\ell\ell'} (\bar{\alpha}_{LM})_{\ell\ell'}^{XX'} / (\sigma_L^2)_{\ell\ell'}^{XX'}}{\sum_{\ell\ell'} 1 / (\sigma_L^2)_{\ell\ell'}^{XX'}}, \qquad (2.15) $$
where $(\sigma_L^2)_{\ell\ell'}^{XX'}$ is the analytic variance associated to each $(\bar{\alpha}_{LM})_{\ell\ell'}$ estimate and can be obtained after a proper re-scaling of the variance of the rotational invariants (see eq. (2.11) for the re-scaling):

$$ (\sigma_L^2)_{\ell\ell'}^A = \frac{C_{AA'}^{\ell\ell'}}{G_{\ell\ell'} F_{\ell\ell'}^{L,A} F_{\ell\ell'}^{L,A'}}, \qquad (2.16) $$

where $A$ represents the considered cross-correlation, i.e. $XX'$, and $C_{AA'}^{\ell\ell'}$ is the entry $AA'$ of the covariance matrix of the rotational invariants (see equation 33 of [53] for the full expression of the covariance matrix).

In this study, we employ the analytic variance to normalize the estimator. While it is important to note that the analytic variance is rigorously justified under the conditions of full-sky observations with homogeneous noise, our investigation has demonstrated that even for cut-sky observations, the analytic expression approximates the true variance of the estimator effectively.

After obtaining the initial estimates for the $\alpha_{LM}$ coefficients, we follow the approach outlined in [24], which encodes for a de-biasing procedure to derive the final estimates of the spherical harmonic coefficients of the CB field. Having the de-biased $\alpha_{LM}$ coefficients allows us to evaluate the CB power spectrum (section 3) and the map of the CB field (section 5.4).

To end up with the un-biased estimates of the $\alpha_{LM}$ we have to subtract the mean field bias. The mean field is a contribution, at the level of maps, coming from mask effects, not homogeneous noise, and other signals of the map different from the CB field.

The evaluation of the mean field bias entirely relies on simulations:

$$ \alpha_{LM}^{\text{bias,MF}} = \langle \bar{\alpha}_{LM} \rangle_{\text{sim}}, \qquad (2.17) $$

where the $\bar{\alpha}_{LM}$ are the spherical harmonic coefficients of the CB field evaluated over the simulations of the CMB maps from equation (2.14), and averaged over the entire simulation set.

Thus, the final estimates of the $\alpha_{LM}$ coefficients have been obtained as:

$$ \hat{\alpha}_{LM} = \bar{\alpha}_{LM} - \alpha_{LM}^{\text{bias,MF}}. \qquad (2.18) $$

Next, we will focus specifically on the EB cross-correlation.

## 2.1 EB-estimator

In the following, we are going to focus on the EB cross-correlation only, showing how the information contained in the CMB polarization field can be used to develop an estimator for the spherical harmonic coefficients of the CB field.

The EB cross-correlation is induced by a rotation of the primordial EE power spectrum, while the TB cross-correlation is generated from the rotation of the primordial TE power spectrum. Observations involving the CMB temperature field are affected by the cosmic variance. For the Planck satellite, as well as for the forthcoming CMB experiments, the signal-to-noise ratio for the EB signal is larger than the one for TB. For this reason, in this study we implement the estimator based on the information coming from the EB signal only.

Furthermore, in section 3, starting from the estimates of the spherical harmonic coefficients of the rotation field, we present the procedure to estimate the CB power spectrum.
The EB cross-correlation, assuming that B-modes on the last scattering surface are null, is:

$$ \begin{aligned} \langle a_{\ell m}^{E,tot} a_{\ell' m'}^{B,tot,*} \rangle &= \langle (a_{\ell m}^E + \delta a_{\ell m}^E) \delta a_{\ell' m'}^{B,*} \rangle = \\ &= \langle a_{\ell m}^E \delta a_{\ell' m'}^{B,*} \rangle + \langle \delta a_{\ell m}^E \delta a_{\ell' m'}^{B,*} \rangle \approx \langle a_{\ell m}^E \delta a_{\ell' m'}^{B,*} \rangle, \end{aligned} \quad (2.19) $$

where the last equivalence holds since the term $\langle \delta a_{\ell m}^E \delta a_{\ell' m'}^{B,*} \rangle$ is second order and we neglect it. From equation (2.19) we see that the EB cross-correlation is determined by the rotated B-modes, thus it is different from zero if $\ell + \ell' + L$ is even (see equation (2.3)).

Equation (2.7) for the specific case of the EB cross-correlation now reads:

$$ \langle a_{\ell m}^{E,tot} a_{\ell' m'}^{B,tot,*} \rangle = 2 \sum_{LM} \alpha_{LM} C_{\ell}^{EE} \xi_{\ell m \ell' m'}^{LM} H_{\ell \ell'}^{L}, \quad (2.20) $$

where $C_{\ell}^{EE}$ refers to the $Z_{\ell\ell'}^{XX'}$ term (eq. (2.7)) when $XX' = EB$.

Following the logical steps previously described, the first expression of the estimator for the $\alpha_{LM}$ coefficients before the subtraction of the mean field bias (eq. (2.15)) is calculated as:

$$ \bar{\alpha}_{LM}^{EB} = \frac{\sum_{\ell\ell'} (\bar{\alpha}_{LM})_{\ell\ell'}^{EB} / (\sigma_L^2)_{\ell\ell'}^{EB}}{\sum_{\ell\ell'} 1 / (\sigma_L^2)_{\ell\ell'}^{EB}}, \quad (2.21) $$

and the expression for the inverse variance of the estimator for the EB cross-correlation is:

$$ \sigma_L^{-2} = \sum_{\ell \ge \ell'} (1 + \delta_{\ell\ell'})^{-1} G_{\ell\ell'} \left\{ \frac{(F_{\ell\ell'}^{L,EB})^2}{C_{\ell}^{EE, map} C_{\ell'}^{BB, map}} + \frac{(F_{\ell\ell'}^{L,BE})^2}{C_{\ell}^{BB, map} C_{\ell'}^{EE, map}} \right\}. \quad (2.22) $$

Focusing only on the numerator in equation (2.21) (meaning the unnormalized estimator, UN):

$$ \bar{\alpha}_{LM}^{UN} = \sum_{\ell\ell'} \left\{ \frac{F_{\ell\ell'}^{L,EB} \sum_{mm'} a_{\ell m}^{E, map} a_{\ell' m'}^{B, map,*} \xi_{\ell m \ell' m'}^{LM}}{C_{\ell}^{EE, map} C_{\ell'}^{BB, map}} \right\}, \quad (2.23) $$

and exploiting properties of Wigner-3j symbols, it is possible to re-write equation (2.23) as:<fnref n="5" />

$$ \bar{\alpha}_{LM}^{UN} = \int d\hat{n} Y_{LM} \left[ \sum_{\ell m} \frac{C_{\ell}^{EE} a_{\ell m}^{E, map,*} W_{\ell-2} Y_{\ell m}}{C_{\ell}^{EE, map}} \sum_{\ell' m'} \frac{a_{\ell' m'}^{B, map,*} W_{\ell'+2} Y_{\ell' m'}}{C_{\ell'}^{BB, map}} + \sum_{\ell m} \frac{C_{\ell}^{EE} a_{\ell m}^{E, map,*} W_{\ell+2} Y_{\ell m}}{C_{\ell}^{EE, map}} \sum_{\ell' m'} \frac{a_{\ell' m'}^{B, map,*} W_{\ell'-2} Y_{\ell' m'}}{C_{\ell'}^{BB, map}} \right], \quad (2.24) $$

where $\pm_2 Y_{\ell m}$ are spin-2 spherical harmonics and $C_{\ell}^{XX, map}$ indicates the analytical expression for the power spectrum recovered from CMB maps, evaluated as:

$$ C_{\ell}^{XX, map} = C_{\ell}^{XX} W_{\ell}^2 + N_{\ell}^{XX}, \quad (2.25) $$

with $C_{\ell}^{XX}$ the cosmological signal, $W_{\ell}$ the window function, and $N_{\ell}^{XX}$ the noise power spectrum. Note that we are using the superscript UN to indicate that we refer to the first estimate of the estimator (eq. (2.21)) without its normalization (eq. (2.22)).

<footnote n="5">⁵Details about the calculations are provided in appendix A.</footnote>
Following the approach of [25], it is possible to re-write the expression of the estimator in a way so that the computational time is remarkably reduced, defining the two following new objects:

$$Q^E \pm iU^E = \sum_{\ell m} (C_{\ell}^{EE} \bar{a}_{\ell m}^{E,*})_{\pm 2} Y_{\ell m}, \qquad (2.26)$$

$$Q^B \pm iU^B = \sum_{\ell m} (\pm i \bar{a}_{\ell m}^{B,*})_{\pm 2} Y_{\ell m}, \qquad (2.27)$$

with $\bar{a}_{\ell m}^{E,*}$ and $\bar{a}_{\ell m}^{B,*}$ defined as:

$$\bar{a}_{\ell m}^{E,*} = \frac{a_{\ell m}^{E, map,*}}{C_{\ell}^{EE, map}} W_{\ell}, \qquad (2.28)$$

$$\bar{a}_{\ell m}^{B,*} = \frac{a_{\ell m}^{B, map,*}}{C_{\ell}^{BB, map}} W_{\ell}. \qquad (2.29)$$

By making use of equations (2.26) and (2.27), the unnormalized estimator for the $\alpha_{LM}$ coefficients can be written as:

$$\bar{\alpha}_{LM}^{UN} = \int d\hat{n} Y_{LM} [2(Q^E U^B - U^E Q^B)]. \qquad (2.30)$$

At this point we define the quantity inside the square brackets as a “map-like” object:

$$m'(\alpha) = 2(Q^E U^B - U^E Q^B), \qquad (2.31)$$

and, after performing the complex conjugate of equation (2.30), we obtain:

$$\bar{\alpha}_{LM}^{UN,*} = \int d\hat{n} Y_{LM}^* m'^*(\alpha) = \int d\hat{n} Y_{LM}^* m'(\alpha), \qquad (2.32)$$

where the second equivalence holds since maps are real objects, i.e. $m'^*(\alpha) = m'(\alpha)$. We need to perform the complex conjugate of equation (2.30) since otherwise we do not have the correct relation that allows to move from the map to the spherical harmonic coefficients.

The above equation is what allows the reduction of the computational time since, having defined the map-like object $m'(\alpha)$ in eq. (2.31) the computation of the associated spherical harmonic coefficients is straightforward.

A word of caution before proceeding. Equation (2.32) provides for the *complex conjugate* of the unnormalized $\alpha_{LM}$ estimates, thus the estimates of the $\alpha_{LM}$ coefficients before the de-biasing procedure are obtained as:

$$\bar{\alpha}_{LM} = \frac{(\bar{\alpha}_{LM}^{UN,*})^*}{\sigma_L^{-2}}. \qquad (2.33)$$

To end up with the final estimates of the spherical harmonic coefficients of the CB field, i.e. $\hat{\alpha}_{LM}$, we have to subtract the mean field bias (eq. (2.17)) from equation (2.33), as presented in equation (2.18).
### 3 Angular power spectrum and de-biasing procedure

Having the estimates of the $\alpha_{LM}$ coefficients, the $\alpha\alpha$ power spectrum can be written as:

$$C_L^{\hat{\alpha}\hat{\alpha}} = \frac{1}{f_{\text{sky}}} \frac{1}{2L+1} \sum_M \hat{\alpha}_{LM} \hat{\alpha}_{LM}^*, \qquad (3.1)$$

where $f_{\text{sky}}$ is the sky fraction of the mask used for the analysis.

The estimate of the CB power spectrum in equation (3.1) is intrinsically biased, that is, $C_L^{\hat{\alpha}\hat{\alpha}} \neq \hat{C}_L^{\alpha\alpha}$. The bias comes from the diagonal contribution of $\ell - \ell'$ pairs with $\ell = \ell'$ when combining together the two estimates of the $\alpha_{LM}$ coefficients and from the off-diagonal contribution from sources different from the rotation induced by CB. We thus encode for a de-biasing procedure [56], now at the level of the power spectrum, to, first, select the diagonal contribution (i.e. the contribution from $\ell - \ell'$ pairs with $\ell = \ell'$) and, second, among the off-diagonal contributions, select the contribution coming only from the rotation induced by CB:

$$\hat{C}_L^{\alpha\alpha} = C_L^{\hat{\alpha}\hat{\alpha}} - C_L^{\text{bias}}. \qquad (3.2)$$

The bias term accounts for two different contributions:

* the isotropic bias term, $C_L^{\text{bias,iso}}$, which is an analytic bias term calculated on data, so that it selects the diagonal contribution coming from the $\ell - \ell'$ pairs with $\ell = \ell'$;

* the Monte Carlo bias term, $C_L^{\text{bias,MC}}$. This term is based on Monte Carlo simulations, in order to describe the off-diagonal signal coming from other contributions different from CB, such as not homogeneous noise, cut-sky effects and the contribution coming from lensing. At the level of the spherical harmonics coefficients of the CB field, lensing has no contribution since CB and lensing are orthogonal effects ([57]). However, at the power spectrum level, this assertion does not hold true. Therefore, we include a bias term to account for and subtract the contribution of lensing.

For this specific case, meaning for the application of the pipeline to Planck data, the de-biasing procedure to obtain the final estimate of the power spectrum can entirely rely on simulations, without the need of the analytic computation of the bias. Despite that, in this work we present the general de-biasing procedure that encodes for the analytic bias term too.

The isotropic bias term is defined as:

$$C_L^{\text{bias,iso}} = \langle \hat{\alpha}_{LM} \hat{\alpha}_{LM}^* \rangle. \qquad (3.3)$$

From a general point of view, the estimator (eq. (2.23)) probes three disconnected Wick contractions [56]: EE-BB, EB-EB, and EB-BE. The cross-correlation terms (i.e. EB-EB and EB-BE) are negligible with respect to the auto-correlation terms (i.e. EE-BB), both with and without a rotation signal induced by CB. Under this assumption, the equation for the isotropic bias term reduces to:

$$C_L^{\text{bias,iso}} = \frac{1}{\sigma_L^{-2} \sigma_L^{-2}} \sum_{\ell' \ge \ell} (1 + \delta_{\ell\ell'})^{-1} G_{\ell\ell'} \left\{ \frac{(F_{\ell\ell'}^{L,EB})^2 \hat{C}_{\ell}^{EE,\text{map}} \hat{C}_{\ell'}^{BB,\text{map}}}{C_{\ell'}^{BB,\text{map}} C_{\ell}^{EE,\text{map}} C_{\ell'}^{BB,\text{map}} C_{\ell}^{EE,\text{map}}} + \frac{(F_{\ell\ell}^{L,BE})^2 \hat{C}_{\ell'}^{EE,\text{map}} \hat{C}_{\ell}^{BB,\text{map}}}{C_{\ell}^{BB,\text{map}} C_{\ell'}^{EE,\text{map}} C_{\ell}^{BB,\text{map}} C_{\ell'}^{EE,\text{map}}} \right\}, \qquad (3.4)$$
where $\hat{C}_{\ell}^{XX, map}$ with $XX = \{EE, BB\}$ are the power spectra estimated from the CMB polarization maps, and corrected by the $1/f_{sky}$ of the applied mask (see section 4 for details about the masks employed in this analysis). Instead, we use the $C_{\ell}^{XX, map}$ (without the over-hat symbol) for the analytic expression of the power spectrum including both the cosmological signal and the noise contributions (see equation (2.25)).

The Monte Carlo bias term is based on simulations with the aim of evaluating the off-diagonal signal induced by not homogeneous noise, cut-sky effects and lensing. For this reason, the $C_{L}^{bias,MC}$ is computed over a set of simulations that resemble the data, masked with the fiducial analysis mask used for data themselves, and which do not contain a rotation signal induced by CB:

$$C_{L}^{bias,MC} = \langle C_{L}^{\hat{\alpha}\hat{\alpha}} - C_{L}^{bias,iso} \rangle_{\text{sims}}. \qquad (3.5)$$

Here the brackets indicate the average computed over the simulations. Since, by construction, all the simulations do not have off-diagonal contributions coming from CB and the diagonal contribution is erased by the isotropic bias term for each simulation, the only off-diagonal signal has to come from the correlations induced by not homogeneous noise, cut-sky effects and lensing.

The full expression for the bias term in equation (3.2) has to encode for both the isotropic and Monte Carlo bias terms (eqs. (3.4) and (3.5)):

$$C_{L}^{bias} = C_{L}^{bias,iso} + C_{L}^{bias,MC}. \qquad (3.6)$$

This methodology has been validated and the results of this validation are reported in appendix B. Furthermore, in appendix C we perform jack-knife tests on *Planck* data to test for the robustness against the presence of systematic errors. In the following sections we present the dataset used for the analysis and the application of the aforementioned pipeline to *Planck* data.

# 4 Data set and simulations

In this section we describe the data products employed during this work. The results that will be presented in the following have been obtained on the Public Release 3 (PR3) [5] and the Public Release 4 (NPIPE) [54] data products of the *Planck* satellite, which contain CMB data and simulation maps at the HEALPix<fnref n="6" /> [58] resolution of $N_{side}= 2048$. The CMB maps employed for the main results have been cleaned using the official *Planck* component separation method Commander [59]. However, we also present a comparison between the different component separation methods as SEVEM, SMICA and NILC.

*Planck* NPIPE has 400 CMB polarization+noise simulations and 100 CMB temperature+noise simulations for the component separation method Commander, and 600 CMB+noise simulations for the component separation method SEVEM; all for *full mission* data and for the two data splits, A and B.

The available simulations in PR3 are CMB-only and noise simulations. The former accounts for 1000 CMB Monte Carlo simulations obtained using the *Planck* ΛCDM best-fit

<footnote n="6">⁶http://healpix.sourceforge.net.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Different masks used in this work. From left to right: *Planck* NPIPE standard mask which retains 78% of the sky; *Planck* PR3 standard mask which retains 75% of the sky; NPIPE standard plus galactic mask with 20% of sky coverage; NPIPE standard plus galactic mask with 40% sky coverage; NPIPE standard plus galactic mask with 60% sky coverage.</figcaption></figure>

model. For what concerns noise, there are 300 noise simulations for the *half-mission 1*, 300 for the *half-mission 2* and 300 for the *full mission*. The simulation set used in this work is divided in three subsets, obtained adding the CMB Monte Carlo simulations and the 300 noise simulations for the two half missions and for the full mission.

The CMB maps employed for the main analysis of *Planck* NPIPE data products are masked with the *Planck* fiducial analysis mask corresponding to a sky fraction of $f_{\text{sky}} = 78\%$ (first mask in the first row of figure 1), while the mask used for *Planck* PR3 data products is characterized by a sky fraction of $f_{\text{sky}} = 75\%$ (second mask in the first row of figure 1). As presented in section 5.3, we also test for different sky coverages, meaning for different masks applied to CMB data and simulations. Figure 1 shows the masks used in this work.

# 5 Results

We now present the application of our pipeline to *Planck full mission* data and to different choices of data splits. Our primary findings have been obtained from *Planck* NPIPE and PR3 data products, cleaned with the component separation method *Commander*. Furthermore, we present a consistency check comparing the CB power spectra obtained from *full mission* data products cleaned with the other component separation methods, i.e. *SEVEM*, *SMICA* and *NILC* for *Planck* PR3 and only *SEVEM* for *Planck* NPIPE. Subsequently, with the obtained $\alpha_{LM}$ estimates, we discuss the procedure and present the results to end up with the map of the CB field and its cross-correlations with the CMB temperature and polarization fields.

From now on, we use [*A*], [*B*] superscripts to indicate whether a quantity has been evaluated from split A or B of *Planck* NPIPE data splits, respectively, and [*1*] and [*2*] superscripts to distinguish among the *half-mission 1* or *half-mission 2* of *Planck* PR3. If no superscript is present, it means that the quantity has been evaluated from *full mission* data.
All the following results that are presented have been obtained using CMB polarization maps at the full resolution of $N_{side}= 2048$. Furthermore, unless otherwise stated, we exclude from our analysis the first 50 CMB multipoles ($\ell_{min}^{CMB} = 50$) and the maximum multipole included in the analysis is $\ell_{max}^{CMB} = 2000$. The CB power spectrum has been evaluated from $L_{min}^{CB} = 0$ up to $L_{max}^{CB} = 1500$.

The main steps of the analysis are the following:

* evaluate the first estimates of the spherical harmonic coefficients of the CB field (eq. (2.33)) from all the CMB maps, i.e. both data and simulations;

* compute the mean field bias (eq. (2.17)) averaging the $\bar{\alpha}_{LM}$ coefficients estimated from simulations only;

* subtract the mean field bias from all the $\bar{\alpha}_{LM}$ estimates, both the ones from data and that from simulations;

* calculate the $C_L^{\hat{\alpha}\hat{\alpha}}$ (eq. (3.1)) and the $C_L^{\text{bias,iso}}$ (eq. (3.4)) for both data and simulations of Planck data products;

* split the $C_L^{\hat{\alpha}\hat{\alpha}}$ and the $C_L^{\text{bias,iso}}$ evaluated from simulations only, into two equally sized sets:

- Set A, used for the evaluation of the Monte Carlo bias term (eq. (3.5));

- Set B, used to obtain a final set of fully de-biased simulations. In particular, for each simulation of this set, we evaluate the unbiased CB power spectrum, $\hat{C}_L^{\alpha\alpha}$, subtracting the isotropic bias term, $C_L^{\text{bias,iso}}$ evaluated from the same set, and the Monte Carlo bias term, $C_L^{\text{bias,MC}}$, calculated from set A;

* subtract the isotropic bias term evaluated from data, $C_L^{\text{bias,iso}}$, and the Monte Carlo bias term evaluated from set A, $C_L^{\text{bias,MC}}$, from the biased $\alpha\alpha$ power spectrum calculated on Planck data, in order to obtain the estimated CB power spectrum, $\hat{C}_L^{\alpha\alpha}$ (eq. (3.2)).

## 5.1 Full mission

We present the application of our pipeline in case of a CB power spectrum obtained combining together $\alpha_{LM}$ coefficients estimated from Planck NPIPE full mission data.

The equation for the $\alpha\alpha$ power spectrum before the de-bias is the same of equation (3.1), where $f_{\text{sky}}$ is the sky fraction of the mask (first mask in the first row of figure 1) applied to CMB polarization maps, both data and simulations, and corresponding to $f_{\text{sky}} = 0.78$.

In figure 2 we plot the power spectra of the different terms. The black and green curves represent the biased $\alpha\alpha$ power spectrum and the isotropic bias term, respectively, evaluated using the $\alpha_{LM}$ estimates from Planck NPIPE data, while the orange curve represents the Monte Carlo bias term, computed by averaging the de-biased $\alpha\alpha$ power spectra evaluated over the first 200 CMB+noise simulations, that is, evaluated from the simulation set A. Even though we show the estimate of the $\alpha\alpha$ power spectrum in figure 3, we display it with the gray curve also in figure 2.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Power spectra evaluated from *full mission* NPIPE data and simulations: $\alpha\alpha$ power spectrum before the de-biasing procedure, evaluated from *Planck* data (black curve); isotropic bias term evaluated from *Planck* data (green curve); Monte Carlo bias term evaluated from the simulation set A (orange curve); $\alpha\alpha$ power spectrum after the de-biasing procedure (gray curve).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Results evaluated from *full mission* NPIPE data. *Upper panel:* de-biased $\alpha\alpha$ power spectrum (black line), with the $1\sigma$ (dark green area), $2\sigma$ (light green area) and $3\sigma$ (gray area) confidence levels. *Lower panel:* de-biased $\alpha\alpha$ power spectrum after binning with 100 multipoles per bin and excluding the first 8 multipoles.</figcaption></figure>
The black curve in the upper panel of figure 3 is the $\alpha\alpha$ power spectrum evaluated for each multipole after the de-biasing procedure (eq. (3.2)). The shaded areas represent the $1\sigma$, $2\sigma$ and $3\sigma$ confidence intervals, obtained after the computation of the variance of the fully de-biased $\alpha\alpha$ power spectra evaluated over the simulation set B. In order to better visualize the results, in the lower panel of figure 3 we also plot the de-biased CB power spectrum after binning with 100 multipoles per bin, excluding the first 8 multipoles; the error bars are at $1\sigma$.

In figure 2 the rotation signal induced by CB before the de-biasing procedure, i.e. $C_L^{\hat{\alpha}\hat{\alpha}}$, is different from zero. We can understand the reason of this strong deviation from zero since the single estimate for the $\alpha_{LM}$ coefficients goes as:

$$ \hat{\alpha}_{LM} \propto a_{\ell m}^{E, map} a_{\ell' m'}^{B, map,*}, \qquad (5.1) $$

and when combining two estimates to obtain the power spectrum before the de-bias, we end up with:

$$ C_L^{\hat{\alpha}\hat{\alpha}} \propto (a_{\ell_1 m_1}^{E, map} a_{\ell_2 m_2}^{B, map,*})(a_{\ell_3 m_3}^{E, map,*} a_{\ell_4 m_4}^{B, map}). \qquad (5.2) $$

The only non negligible contributions come from:

$$ a_{\ell_1 m_1}^{E, map} a_{\ell_3 m_3}^{E, map,*} \sim C_{\ell_1}^{EE, map}, \qquad (5.3) $$

$$ a_{\ell_2 m_2}^{B, map} a_{\ell_4 m_4}^{B, map,*} \sim C_{\ell_2}^{BB, map}. \qquad (5.4) $$

The observed EE and BB power spectra encode for two contributions; the cosmological signal, i.e. $C_{\ell}^{EE}$ and $C_{\ell}^{BB}$, and the noise contribution, i.e $N_{\ell}^{EE}$ and $N_{\ell}^{BB}$, so that:

$$ C_{\ell}^{EE, map} = C_{\ell}^{EE} + N_{\ell}^{EE}, \qquad (5.5) $$

$$ C_{\ell}^{BB, map} = C_{\ell}^{BB} + N_{\ell}^{BB}. \qquad (5.6) $$

In case of $a_{\ell m}^{X, map}$ observed from the same data set, the noise auto-correlates. For the Planck sensitivity the dominant contribution to the observed EE and BB power spectra comes from the noise itself, and the bias in figure 2 is dominated by the auto-correlation of the noise. Nevertheless, it is worth stressing that even in a signal dominated regime the de-bias procedure is necessary for removing the signal auto-correlation.

The CB power spectrum after the de-biasing procedure obtained applying our pipeline to full mission Planck NPIPE data is compatible with zero with a Probability To Exceed (PTE) of 84.35%. This was determined from the $\chi^2$ of the null hypothesis computed in harmonic space, considering a covariance matrix evaluated from a set of fully debiased binned simulations.

To further quantify the implications of our findings, we estimated an upper limit on the amplitude, $A_{CB}$, of a scale-invariant CB power spectrum in $D_L \equiv L(L+1)C_L/2\pi$. This was done by minimizing the gaussian likelihood evaluated using the same covariance matrix utilized for the computation of the $\chi^2$ described above. This analysis yielded an upper bound of $A_{CB} < 0.09 \text{ deg}^2$, which is compatible with [27].

Our results, indicating a null detection of anisotropic CB, are in concordance with those obtained from ACT-Pol ([25]) and SPT-Pol ([24]) analyses. Furthermore, [25] and [24] also
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Left panel: de-biased $\alpha\alpha$ power spectrum after binning with 100 multipoles per bin and excluding the first 8 multipoles using *Planck* NPIPE *full mission* data. We compare the different component separation methods; *Commander* in red and *SEVEM* in green. Right panel: same as the left panel, but for *Planck* PR3 data and for the four component separation methods available for PR3; in blue we show the power spectrum obtained using *SMICA* component separation method, in orange *NILC*, in green *SEVEM* and in red *Commander*.</figcaption></figure>

reported constraints on the amplitude of a scale-invariant CB power spectrum, setting upper limits of $A_{CB} < 0.033 \text{ deg}^2$. Our upper bound is less stringent than the ones evaluated from ACT and SPT, remaining however robust against foregrounds thanks to the component separation analysis. Furthermore, it is important to notice that, even though their estimates are favoured by the lower noise, especially at small scales, *Planck* results are competitive with ACT-Pol and SPT-Pol due to the wider sky fraction covered.

In the left and right panels of figure 4 we show the de-biased $\alpha\alpha$ power spectrum evaluated from *full mission Planck* NPIPE and PR3 data products respectively, for the different component separation methods. In the left panel of figure 4 we show the de-biased power spectrum after binning with 100 multipoles per bin and excluding the first 8 multipoles, evaluated from *Planck* NPIPE data products, for the two component separation methods available, i.e., for *Commander*, in red, and for *SEVEM*, in green. The estimated CB power spectrum is consistent among the different component separation methods. In the right panel of figure 4 we show the same de-biased $\alpha\alpha$ power spectrum evaluated from *Planck* PR3 data products for the four component separation methods, meaning *SMICA* (blue), *NILC* (orange), *SEVEM* (green) and *Commander* (red). The CB power spectrum estimated from *Planck* PR3 *Commander* is compatible with zero with a PTE of 6.61%. This low PTE reflects the behaviour observed in the data at large multipoles (see right panel of figure 4), where there appears to be an excess of power. We speculate that this could be addressed to a mismatch between the noise in the CMB simulations and the one of data of *Planck* PR3 data products.

In the left and right panels of figure 5 we also show the correlation matrix evaluated from the $\alpha\alpha$ power spectra estimated from the 400 CMB+noise simulations of *Planck* NPIPE and from the 300 CMB+noise simulations of *Planck* PR3, respectively, after binning with 100 multipoles per bin. The multipole $L_{cent}$ reported in the axes of the matrix is the center of the multipole bin. The correlation coefficients reported in figure 5 are expressed in terms of percentage.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Left panel: correlation matrix evaluated from the 400 CB power spectra estimated from the simulations of *Planck* NPIPE data products. Right panel: correlation matrix evaluated from the 300 CB power spectra estimated from the simulations of *Planck* PR3 data products. All simulations are binned with 100 multipoles per bin. The correlation coefficients are expressed in terms of percentage.</figcaption></figure>

## 5.2 Data splits

We summarize the results of the application of our pipeline to different combinations of data splits. In particular, we show the CB power spectrum obtained:

* auto-correlating $\alpha_{LM}$ coefficients estimated from data split A (half mission 1) of *Planck* NPIPE (PR3) data products;

* auto-correlating $\alpha_{LM}$ coefficients estimated from data split B (half mission 2) of *Planck* NPIPE (PR3) data products;

* cross-correlating $\alpha_{LM}$ coefficients evaluated from data split A (half mission 1) and data split B (half mission 2) of *Planck* NPIPE (PR3) data products.

A word of caution concerning the case of the cross-correlation (since both auto-correlations follow exactly the same pipeline of the *full mission* case). The main differences with respect to the previously described case are the following:

* the calculation of the $\alpha\alpha$ power spectrum before the debias (eq. (3.1)), as well as the one of the isotropic bias term (eq. (3.4)), involves estimates of the $\alpha_{LM}$ coefficients coming from the two data splits (or from the two half-missions, if we work with *Planck* PR3 data products);

* all quantities involving window functions and noise curves, i.e. $F_{\ell\ell'}^{L,EB}$, $F_{\ell\ell'}^{L,BE}$ and the analytic power spectra, $C_{\ell}^{XX,\text{map}}$, must be evaluated separately for each data split (or each half-mission).
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Left panel: de-biased $\alpha\alpha$ power spectrum after binning with 100 multipoles per bin, excluding the first 8 multipoles from the analysis, obtained auto-correlating $\alpha_{LM}$ estimates from data split A (blue), data split B (orange) and cross-correlating estimates from the two data splits (green) of Planck NPIPE data products. Right panel: same as the left panel, for the auto-correlation of $\alpha_{LM}$ estimates from half mission 1 (blue), half mission 2 (orange) and for the cross-correlation from the two half missions (green) of Planck PR3 data products.</figcaption></figure>

Thus, the biased power spectrum of equation (3.1) is now calculated as:

$$C_L^{\hat{\alpha}\hat{\alpha}} = \frac{1}{f_{\text{sky}}} \frac{1}{2L+1} \sum_M \hat{\alpha}_{LM}^{[A]} \hat{\alpha}_{LM}^{[B],*} \quad (5.7)$$

And, since we are combining together estimates coming from two different data sets, the normalization of the $\hat{\alpha}_{LM}$ coefficients is different depending on whether we are dealing with data split A (half-mission 1) or data split B (half-mission 2) estimates. The distinction follows the same notation, i.e. we indicate the inverse variance from the data split A (half-mission 1) as $(\sigma_L^{-2})^{[A]} ((\sigma_L^{-2})^{[1]})$ and as $(\sigma_L^{-2})^{[B]} ((\sigma_L^{-2})^{[2]})$ the one from the data split B (half-mission 2). The same modification applies to the isotropic bias term (3.4), which now reads as:

$$C_L^{\text{bias,iso}} = \langle \hat{\alpha}_{LM}^{[A]} \hat{\alpha}_{LM}^{[B],*} \rangle \quad (5.8)$$

In the left panel of figure 6 we show the de-biased CB power spectrum after binning with 100 multipoles per bin, excluding the first 8 multipoles, for the different combinations of data splits of Planck NPIPE; the error bars are at $1\sigma$. More precisely, the presented results show the $\alpha\alpha$ power spectra obtained auto-correlating $\alpha_{LM}$ estimates from split A (blue), auto-correlating estimates from split B (orange) and cross-correlating $\alpha_{LM}$ estimates from both split A and split B (green).

In the right panel of figure 6 we show the application of the pipeline to the different choices of data splits of Planck PR3 data products; the auto-correlation from half mission 1 (blue), the auto-correlation from half mission 2 (orange) and the cross-correlation between half mission 1 and half mission 2 (green).

The results presented in figure 6 have been obtained using the Planck NPIPE and PR3 data products cleaned with the component separation method Commander. In this section we present the results for Commander only since our analysis shows consistency among the
<table><caption>**Table 1.** Probability To Exceed for the CB power spectra estimated for the different data splits presented in figure 6. For completeness, in the last row we also report the PTEs for the CB power spectra estimated from *Planck* NPIPE and *Planck* PR3 *full mission* data products.</caption><thead><tr><th colspan="2">NPIPE Commander</th><th colspan="2">PR3 Commander</th></tr></thead><tbody><tr><td>α<sup>[A]</sup>α<sup>[A]</sup></td><td>71.01%</td><td>α<sup>[1]</sup>α<sup>[1]</sup></td><td>75.95%</td></tr><tr><td>α<sup>[B]</sup>α<sup>[B]</sup></td><td>7.77%</td><td>α<sup>[2]</sup>α<sup>[2]</sup></td><td>50.48%</td></tr><tr><td>α<sup>[A]</sup>α<sup>[B]</sup></td><td>83.08%</td><td>α<sup>[1]</sup>α<sup>[2]</sup></td><td>31.88%</td></tr><tr><td>αα</td><td>84.35%</td><td>αα</td><td>6.61%</td></tr></tbody></table>

different component separation methods (see left and right panels of figure 4). The CB power spectra estimated for the different data splits are compatible with zero with the PTE reported in table 1.

Regarding the low PTEs associated with $α^{[2]}α^{[2]}$ and $α^{[B]}α^{[B]}$, it is crucial to recognize that they correspond to distinct types of data splits. In the case of *Planck* PR3, the splits are time-based, whereas for *Planck* NPIPE we are working with detector-based data splits.

Regarding the lowest PTEs listed in table 1, specifically those associated with $α^{[2]}α^{[2]}$ and $α^{[B]}α^{[B]}$, it is crucial to recognize that they correspond to distinct types of data splits. In the case of *Planck* PR3, the splits are time-based, whereas for *Planck* NPIPE we are working with detector-based data splits.

## 5.3 Consistency checks

In the following, we consider the specific case of *full mission Planck* NPIPE data products and we go through three consistency checks. In particular, we compare the de-biased $αα$ power spectrum:

* for four different choices of the minimum CMB multipole included in the analysis, considering the cases where $ℓ_{min}^{CMB} = 10, ℓ_{min}^{CMB} = 30, ℓ_{min}^{CMB} = 50$ and $ℓ_{min}^{CMB} = 100$;

* for three different choices of the maximum CMB multipole included in the analysis, encoding for $ℓ_{max}^{CMB} = 1500, ℓ_{max}^{CMB} = 2000$ and $ℓ_{max}^{CMB} = 2500$;

* for different masks applied to *Planck* NPIPE CMB polarization maps, corresponding to the sky-fractions of $f_{sky} = 78.0\%, 59.2\%, 39.7\%, 20.3\%$.

With reference to figure 7, it is possible to conclude that the obtained results exhibit consistency across various choices of minimum and maximum CMB multipole included in the analysis.

The CB power spectrum results consistent also for the different masks applied to CMB data and simulations, as displayed in figure 8. It is worth noting that, despite the results showing strong consistency across different masks, the entire analysis normalized the estimator using its analytic variance instead of the true variance. As we observe a smaller portion of
<figure><img src="image_9.png" /><figcaption>**Figure 7.** CB power spectrum evaluated from *Planck* NPIPE *full mission* data products, binning with 100 multipoles per bin and excluding the first 8 multipoles. **Left panel:** power spectra for different values of $\ell_{\min}^{\text{CMB}}$ included in the analysis and $\ell_{\max}^{\text{CMB}} = 2000$ for all cases. We indicate in blue the power spectrum obtained with $\ell_{\min}^{\text{CMB}} = 10$, in orange the one with $\ell_{\min}^{\text{CMB}} = 30$, in green the one with $\ell_{\min}^{\text{CMB}} = 50$ and in red the one with $\ell_{\min}^{\text{CMB}} = 100$. **Right panel:** power spectra for different values of $\ell_{\max}^{\text{CMB}}$ included in the analysis and $\ell_{\min}^{\text{CMB}} = 50$ for all cases. We indicate in blue, orange and green the CB power spectra obtained including CMB multipoles up to $\ell_{\max}^{\text{CMB}} = 1500$, $\ell_{\max}^{\text{CMB}} = 2000$ and $\ell_{\max}^{\text{CMB}} = 2500$, respectively.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** CB power spectrum evaluated from *Planck* NPIPE *full mission* data products, binning over 100 multipoles per bin and excluding the first 8 multipoles, for different masks applied to *Planck* maps. In blue we show the power spectrum obtained with the standard mask, corresponding to a $f_{\text{sky}} = 78.0\%$; in orange, green and red we show the power spectrum obtained applying masks corresponding to $f_{\text{sky}} = 59.2\%, 39.7\%, 20.3\%$ respectively.</figcaption></figure>

the sky, the difference between the analytic and the true variance becomes more pronounced. Specifically, we start observing a deviation from the analytic variance when applying the galactic mask of $f_{\text{sky}} = 20\%$.

## 5.4 CB cross-correlation

In this section, we present the CB field map, derived using the de-biased estimates of its spherical harmonic coefficients (eq. (2.18)). This map is illustrated in figure 9. We
<figure><img src="image_11.png" /><figcaption>**Figure 9.** CB map with a 1deg smoothing evaluated from *Planck* NPIPE data products on the left, and the one evaluated from *Planck* PR3 data products on the right.</figcaption></figure>

then delve into the cross-correlation analysis between this field and the CMB temperature and polarization fields. Figures 10 detail these results, showcasing data from both the *Planck* NPIPE (in green) and PR3 (in red) datasets.

The CB maps evaluated from *Planck* NPIPE (left panel in figure 9) and PR3 (right panel in figure 9) data products are processed using a 1° FWHM Gaussian beam smoothing. This processing follows the subtraction of the mean field term from the $\alpha_{LM}$ estimates, based on both *Planck* NPIPE (map on the left) and PR3 (map on the right) data. The final maps have been masked using the standard masks of *Planck* NPIPE (first mask in figure 1) and *Planck* PR3 (second mask in figure 1). A notable feature in both CB maps is the correlation of smaller angle multipoles (yellow regions) with the *Planck* scanning strategy, as documented in [60].

For the analysis of the cross-correlation between CB and CMB temperature and polarization fields, we employed the Pymaster Python package [61], which facilitated the calculation of the cross-spectra. In figure 10, the upper panel illustrates the cross-correlation between CB and the CMB temperature field ($\alpha T$) displayed in band-powers, for both *Planck* NPIPE (green points) and PR3 (red points) datasets. The lower panels display the cross-correlation of the CB field with the CMB polarization fields ($\alpha E$ and $\alpha B$). The power spectra displayed in figure 10 are binned with 100 multipoles per bin, excluding the first 8 multipoles and the error bars have been evaluated from simulations. Note that the $\alpha T$ power spectrum is not represented in the units used for $\alpha E$ and $\alpha B$, rather it is showed in band-powers to allow for a more direct comparison with the $\alpha T$ power spectrum presented in [24].

Table 2 presents the Probability To Exceed (PTE) values for the cross-correlation power spectra between the CB and CMB fields. Notably, the $\alpha B$ power spectrum from NPIPE, illustrated in the lower right panel of figure 10, exhibits minimal scatter and this is reflected in its high PTE value (98.75%).

# 6 Sensitivities of future experiments

In this section we present forecasts for forthcoming CMB experiments. The numerical computation of the variance of the estimator is crucial. Not only because it enables the calculation of the $\alpha_{LM}$ coefficients, but also serves as a rapid tool for forecasting future CMB experi-
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Cross-spectra between the CMB fields (temperature and polarization) and the CB field for NPIPE, in green, and PR3, in red, data products. *Upper panel:* $\alpha T$ cross-spectrum in band-powers; *Lower left panel:* $\alpha E$ cross-spectrum; *Lower right panel:* $\alpha B$ cross-spectrum.</figcaption></figure>

<table><caption>**Table 2.** Probability To Exceed for the cross-spectra of the CB field with the temperature and polarization CMB fields. For completeness, in the last row we also report the PTEs for the $\alpha\alpha$ auto-spectrum evaluated in the previous sections, for Planck NPIPE and PR3 data products.</caption><thead><tr><th></th><th>NPIPE Commander</th><th>PR3 Commander</th></tr></thead><tbody><tr><td>$\alpha T$</td><td>8.27%</td><td>18.71%</td></tr><tr><td>$\alpha E$</td><td>79.37%</td><td>16.31%</td></tr><tr><td>$\alpha B$</td><td>98.75%</td><td>56.15%</td></tr><tr><td>$\alpha\alpha$</td><td>84.85%</td><td>6.61%</td></tr></tbody></table>
<table><caption>**Table 3.** Instrumental specifications of the considered CMB experiments.</caption><thead><tr><th colspan="4">CMB experiments specifications</th></tr><tr><th>Experiment</th><th>σ<sub>T</sub><br>(μKarcmin)</th><th>σ<sub>P</sub><br>(μKarcmin)</th><th>θ<sub>fwhm</sub><br>(arcmin)</th></tr></thead><tbody><tr><td>Planck</td><td>40</td><td>56.57</td><td>5</td></tr><tr><td>LiteBIRD</td><td>2.2</td><td>3.26</td><td>30</td></tr><tr><td>Simons Observatory LAT</td><td>6</td><td>8.49</td><td>1.4</td></tr><tr><td>CMB-S4</td><td>3</td><td>4.24</td><td>1</td></tr></tbody></table>

ments. By assessing $\sigma_L$ (eq. (2.22)), which reflects the sensitivity of an experiment to specific cross-correlations, we can identify the most promising avenues for detecting CB signatures.

We present forecasts for the LiteBIRD satellite [9], the Simons Observatory [10] and CMB-S4 [11]. Together with the sensitivities of these forthcoming CMB experiments, we also plot the one associated to the Planck satellite in order to provide a straightforward comparison. In the following, we briefly list the instrumental specifications used for the evaluation of the $\sigma_L$ for each experiment.

The LiteBIRD satellite, whose launch is predicted for the late 2020s, is composed of three telescopes which cover a total frequency range of 34–448 GHz. Its expected total sensitivity is ~2.2μK·arcmin with an angular resolution of 30 arcmin. The ground-based Simons Observatory will cover the frequency range 27-280 GHz and is composed by the Small Aperture Telescopes (SATs) and the Large Aperture Telescope (LAT). In this work we focus on the Large Angular Telescope [62]. Simons Observatory LAT is expected to have a total sensitivity of 6μK·arcmin and an angular resolution of 1.4 arcmin. CMB-S4 is composed of 21 telescopes, 3 are large aperture telescopes and 18 are small aperture telescopes, and it will cover the frequency range 30-300 GHz. Its predicted total sensitivity is of 3μK·arcmin with an angular resolution of 1 arcmin. All the instrumental specifications, relevant for the analysis presented here, are summarised in table 3.

Figure 11 shows the improvement of the LiteBIRD satellite with respect to Planck. Further improvements will be granted by ground based experiments. It is worth also to notice that the sensitivity that characterizes each experiment is not the only factor entering in the computation of $\sigma_L$, but also the angular resolution plays an important role. This can be seen comparing LiteBIRD and Simons Observatory LAT (as well as LiteBIRD and CMB-S4). Of course, experiments with higher resolution will grant access to smaller angular scales.

## 7 Conclusions

In this study we have built an estimator for the spherical harmonic coefficients of the CB field which implements the method described in [53] and exploits the information contained in the CMB EB cross-correlation. Based on the latter, we have built a pipeline aimed at extracting the angular power spectrum of CB from CMB polarization maps. We applied it to both Planck PR3 and NPIPE data products, considering full-mission data and different data
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Variance of the estimator evaluated with the instrumental specifications of the four CMB experiments reported in table 3.</figcaption></figure>

splits. In all cases, our analysis consistently found that the CB power spectrum is compatible with zero at a significance level of approximately $2\sigma$. As expected, *Planck* NPIPE full mission data are slightly more constraining than the corresponding PR3 products.

Our findings agree with the CB power spectrum estimated from other analysis of *Planck* data [23, 26, 27] and various other experiments, including ACT [25], POLARBEAR [63], BICEP2/Keck Array [8] and SPT [24]. Moreover, we carried out a series of consistency checks reinforcing the reliability of our analysis, which we showed to be robust against: 1) the different component separation methods considered; 2) the different choices of minimum and maximum CMB multipoles included in the analysis; 3) the different masks applied to CMB maps.

Additionally, we employed the spherical harmonic coefficients of the CB field estimated from our pipeline to cross-correlate with the CMB temperature and polarization fields, producing $\alpha T$, $\alpha E$, and $\alpha B$ power spectra up to $L = 1500$. Of these, the first two are predicted to be non-null in several models providing a further mean to constrain axion parameters, see e.g. [50, 64].

We have also presented forecast for future CMB experiments showing that they will achieve sensitivities to anisotropic CB order of magnitudes better than what is currently available. In particular, the *LiteBIRD* satellite will reach a factor of $\sim 25$ improvement with respect to the *Planck* at power spectrum level, while *CMB-S4* [11] will be able to reach an improvement of a factor of $\sim 1000$.

The code and the pipeline developed for this work are publicly available,<fnref n="7" /> along with products employed.<fnref n="8" /> Upon requests we can provide additional products such as birefringence maps or $\alpha_{LM}$ coefficients.

<footnote n="7">⁷The code is available on GitHub at https://github.com/paganol/alpha_lm.</footnote>

<footnote n="8">⁸CB spectra from *Planck* PR3 and NPIPE, as well as its cross-correlations with the CMB T-, E- and B-fields are available on GitHub at https://github.com/giorgiazagatti/CB_Planck Maps spectra.git.</footnote>
## Acknowledgments

We thank Margherita Lembo for useful discussions on the implementation of the estimator. We acknowledge the financial support from the INFN InDark initiative and from the COSMOS network (www.cosmosnet.it) through the ASI (Italian Space Agency) Grants 2016-24-H.0 and 2016-24-H.1-2018, as well as 2020-9-HH.0 (participation in LiteBIRD phase A). GZ acknowledges support from ASI/INFN grant no. 2021-43-HH.0. GZ acknowledges financial support by “Bando Giovani anno 2023 per progetti di ricerca finanziati con il contributo 5x1000 anno 2021”. GF acknowledges the support of the European Research Council under the Marie Skłodowska Curie actions through the Individual Global Fellowship No. 892401 PiCOGAMBAS. This work is supported in part by the MUR PRIN2022 Project “BROWSEPOL: Beyond standaRd mOdel With coSmic microwavE background POLarization”-2022EJNZ53 financed by the European Union — Next Generation EU. We acknowledge the use of numpy [65] and matplotlib [66] software packages, and the use of computing facilities at CINECA. Some of the results in this paper have been derived using the NaMaster [61], the healpy [67] and HEALPix [68] packages.

## A Construction of the estimator

In section 2.1 we re-write the part of equation (2.33) without its normalization (i.e., without $\sigma_L^{-2}$), that we indicate as $\bar{\alpha}_{LM}^{UN}$, in order to reduce the computation time.

$$ \bar{\alpha}_{LM}^{UN} = \sum_{\ell \ge \ell'} (1 + \delta_{\ell \ell'})^{-1} \left\{ \frac{F_{\ell \ell'}^{L,EB} \sum_{mm'} a_{\ell m}^{E, map} a_{\ell' m'}^{B, map,*} \xi_{\ell m \ell' m'}^{LM}}{C_{\ell}^{EE, map} C_{\ell'}^{BB, map}} + \frac{F_{\ell \ell'}^{L,BE} \sum_{mm'} a_{\ell m}^{B, map} a_{\ell' m'}^{E, map,*} \xi_{\ell m \ell' m'}^{LM}}{C_{\ell}^{BB, map} C_{\ell'}^{EE, map}} \right\}. \quad (A.1) $$

In this appendix we show the calculations to obtain the implemented expression of the harmonic estimator for the CB field.

Since the second term inside the summation is equal to the first one with $\ell - \ell'$ inverted, an equivalent expression of (A.1), making the expression of the $F_{\ell \ell'}^{L,EB}$ term explicit, is:

$$ \bar{\alpha}_{LM}^{UN} = \sum_{\ell \ell'} \frac{2H_{\ell \ell'}^{L} C_{\ell}^{EE} W_{\ell} W_{\ell'} \sum_{mm'} a_{\ell m}^{E, map} a_{\ell' m'}^{B, map,*} \xi_{\ell m \ell' m'}^{LM}}{C_{\ell}^{EE, map} C_{\ell'}^{BB, map}}. \quad (A.2) $$

In the following, we re-write the $2H_{\ell \ell'}^{L} \xi_{\ell m \ell' m'}^{LM}$ term. Specifying the expression of $\xi_{\ell m \ell' m'}^{LM}$ we end up with:

$$ 2H_{\ell \ell'}^{L} \sqrt{\frac{(2\ell + 1)(2L + 1)(2\ell' + 1)}{4\pi}} \begin{pmatrix} \ell & L & \ell' \\ -m & M & m' \end{pmatrix}, \quad (A.3) $$

and, exploiting the definition of $H_{\ell \ell'}^{L}$ and the properties of the Wigner-3j symbols, recalling that $\ell + L + \ell'$ must be even in the case of the EB cross-correlation induced by CB, it is
also possible to write:

$$2H_{\ell\ell'}^{L} = \begin{pmatrix} \ell & L & \ell' \\ +2 & 0 & -2 \end{pmatrix} + \begin{pmatrix} \ell & L & \ell' \\ -2 & 0 & +2 \end{pmatrix},$$

so that we can use the definition of the triple integral to re-write the product (A.3) as:

$$\sqrt{\frac{(2\ell + 1)(2L + 1)(2\ell' + 1)}{4\pi}} \left[ \begin{pmatrix} \ell & L & \ell' \\ +2 & 0 & -2 \end{pmatrix} + \begin{pmatrix} \ell & L & \ell' \\ -2 & 0 & +2 \end{pmatrix} \right] \begin{pmatrix} \ell & L & \ell' \\ -m & M & m' \end{pmatrix} \\ = \int d\hat{n} Y_{LM}(\hat{n}) [-2Y_{\ell-m}(\hat{n})_{+2}Y_{\ell'm'}(\hat{n}) + {}_{+2}Y_{\ell-m}(\hat{n})_{-2}Y_{\ell'm'}(\hat{n})]. \quad (A.4)$$

Therefore, the estimator can be written as:

$$\bar{\alpha}_{LM}^{UN} = \int d\hat{n} Y_{LM} \left[ \sum_{\ell m} \frac{C_{\ell}^{EE} (-1)^m a_{\ell m}^{E, map} W_{\ell-2} Y_{\ell-m}}{C_{\ell}^{EE, map}} \sum_{\ell' m'} \frac{a_{\ell' m'}^{B, map, *} W_{\ell'+2} Y_{\ell' m'}}{C_{\ell'}^{BB, map}} \right. \\ \left. + \sum_{\ell m} \frac{C_{\ell}^{EE} (-1)^m a_{\ell m}^{E, map} W_{\ell+2} Y_{\ell-m}}{C_{\ell}^{EE, map}} \sum_{\ell' m'} \frac{a_{\ell' m'}^{B, map, *} W_{\ell'-2} Y_{\ell' m'}}{C_{\ell'}^{BB, map}} \right] \quad (A.5)$$

$$= \int d\hat{n} Y_{LM} \left[ \sum_{\ell m} \frac{C_{\ell}^{EE} a_{\ell m}^{E, map, *} W_{\ell-2} Y_{\ell m}}{C_{\ell}^{EE, map}} \sum_{\ell' m'} \frac{a_{\ell' m'}^{B, map, *} W_{\ell'+2} Y_{\ell' m'}}{C_{\ell'}^{BB, map}} \right] \quad (A.6)$$

$$+ \sum_{\ell m} \frac{C_{\ell}^{EE} a_{\ell m}^{E, map, *} W_{\ell+2} Y_{\ell m}}{C_{\ell}^{EE, map}} \sum_{\ell' m'} \frac{a_{\ell' m'}^{B, map, *} W_{\ell'-2} Y_{\ell' m'}}{C_{\ell'}^{BB, map}}, \quad (A.7)$$

where the second equivalence follows from $(a_{\ell m}^X)^* = (-1)^m a_{\ell-m}^X$.

At this point we define two new objects:

$$Q^E \pm iU^E = \sum_{\ell m} (C_{\ell}^{EE} \bar{a}_{\ell m}^{E,*})_{\pm 2} Y_{\ell m}, \quad (A.8)$$

$$Q^B \pm iU^B = \sum_{\ell m} (\pm i \bar{a}_{\ell m}^{B,*})_{\pm 2} Y_{\ell m}, \quad (A.9)$$

with $\bar{a}_{\ell m}^{E,*}$ and $\bar{a}_{\ell m}^{B,*}$ defined as:

$$\bar{a}_{\ell m}^{E,*} = \frac{a_{\ell m}^{E, map, *}}{C_{\ell}^{EE, map}} W_{\ell}, \quad (A.10)$$

$$\bar{a}_{\ell m}^{B,*} = \frac{a_{\ell m}^{B, map, *}}{C_{\ell}^{BB, map}} W_{\ell}. \quad (A.11)$$

Re-writing equation (A.7) in terms of (A.8) and (A.9), we obtain:

$$\bar{\alpha}_{LM}^{UN} = \int d\hat{n} Y_{LM} [2(Q^E U^B - U^E Q^B)]. \quad (A.12)$$

Performing the complex conjugate of the above equation and recognizing a “map-like” object in the term inside the square brackets, we obtain:

$$\bar{\alpha}_{LM}^{UN,*} = \int d\hat{n} Y_{LM}^* m'(\alpha). \quad (A.13)$$
This final expression is particularly useful since, having defined $m'(\alpha)$ in terms of the real Q-like and U-like objects from (A.8) and (A.9), we can directly obtain the associated spherical harmonic coefficients.

A word of caution before proceeding. This procedure allows us to obtain an estimate for the complex conjugate and unnormalized $\alpha_{LM}$ coefficients. We recover the final estimate for the $\alpha_{LM}$ coefficients as:

$$ \bar{\alpha}_{LM} = \frac{(\bar{\alpha}_{LM}^{UN,*})^*}{\sigma_L^{-2}}. \qquad (A.14) $$

## B Validation tests

In this appendix we show some validation tests for our pipeline. In particular, we show that the application of the implemented estimator to a rotated CMB polarization map recovers the same input rotation. On the contrary, we also test the de-biasing procedure, showing that the pipeline produces a CB power spectrum compatible with zero in case of un-rotated input CMB polarization maps.

All the presented results have been obtained using a set of 100 CMB full-sky simulations at a resolution of $N_{side}= 512$, with $\ell_{max}^{CMB} = 512$, $\theta_{fwhm} = 30 arcmin$, $\sigma_{noise}^T = 1\mu Karcmin$ and $\sigma_{noise}^P = \sqrt{2}\sigma_{noise}^T$.

### B.1 Validation with a rotation signal

In order to validate the pipeline in case of a non-zero rotation signal in the input CMB polarization maps, we use a CB power spectrum that is scale invariant in band powers:

$$ D_L^{\alpha\alpha} = \mathcal{A} = \frac{L(L+1)}{2\pi} C_L^{\alpha\alpha}. $$

After a proper choice of the amplitude, $\mathcal{A} = 0.006 \text{ deg}^2$, we can end up with the fiducial CB power spectrum used for the test of the estimator, by inverting the previous relation:

$$ C_L^{\alpha\alpha, \text{fid}} = \frac{2\pi\mathcal{A}}{L(L+1)}. \qquad (B.1) $$

At this point we generate 100 realizations of CMB maps from the fiducial power spectrum and we rotate each CMB realization accordingly to the associated CB field realization. We are generating different CB realizations for each simulation, each one obtained from the same input $\alpha\alpha$ power spectrum. Once we have the rotated CMB maps, we convolve for the beam and add the noise to each realization.

Having the set of rotated CMB+noise simulations, we apply the pipeline described in this work to the simulated maps, evaluating the $\alpha\alpha$ power spectrum before the de-biasing procedure, $C_L^{\hat{\alpha}\hat{\alpha}}$, and the isotropic bias term, $C_L^{\text{bias,iso}}$. Since for the validation part we are working in the ideal case of full-sky and white noise, only the computation of the isotropic bias term is needed.

In the upper panel of figure 12 we compare the input signal and the average of the estimated CB power spectra computed over the simulations. In the lower panel we show the difference between the estimated and the input power spectra. The difference is compatible with zero at $2\sigma$.
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Validation with rotation signal: *Upper panel:* input power spectrum (red curve), average of the estimated $\alpha\alpha$ power spectra (black curve), the $1\sigma$, $2\sigma$ and $3\sigma$ confidence intervals are the dark green, light green and gray areas, respectively. *Lower panel:* difference between the input and the recovered power spectrum.</figcaption></figure>

## B.2 Validation without a rotation signal

For this second part of the validation we apply our pipeline to the 100 full-sky CMB polarization maps without rotating them.

In figure 13 we compare the average of the $\alpha\alpha$ power spectra computed over the simulations (black curve) with the expected zero rotation signal (red curve). The shaded areas are the $1\sigma$, $2\sigma$ and $3\sigma$ confidence intervals. The average de-biased power spectrum is compatible with zero at $2\sigma$.

## C Jack-knife tests

In this appendix, we detail the outcomes of two different jack-knife tests conducted on Planck NPIPE data.

For the first test, we consider the maps obtained from the difference between split A and B, for both data and simulations:

$$m_{JK} = m_A - m_B, \qquad (C.1)$$

and we extract the CB spectrum applying the pipeline presented in this study on this new set of maps. Ideally, $m_{JK}$ should contain only noise and therefore its CB spectrum should be compatible with zero. However, in case of systematics, deviations from null might show up.
<figure><img src="image_15.png" /><figcaption>Figure 13. Average of the estimated $\alpha\alpha$ power spectra (black curve) from simulated CMB realizations without rotation. The $1\sigma$, $2\sigma$ and $3\sigma$ confidence intervals are the dark green, light green and gray areas, respectively.</figcaption></figure>

In the left panel of figure 14, we show the CB power spectrum estimated from the set of $m_{JK}$, which is found to be consistent with a null hypothesis, exhibiting a PTE of 33.13%. From this analysis, we do not find a significant impact of systematic effects on the presented results.

The second jack-knife test consists in the estimation of the CB power spectrum from a set of $\alpha_{LM}$ obtained as:

$$ \hat{\alpha}_{LM}^{[JK]} = \hat{\alpha}_{LM}^{[A]} - \hat{\alpha}_{LM}^{[B]}, \qquad (C.2) $$

where $\alpha_{LM}^{[A]}$ and $\alpha_{LM}^{[B]}$ are the estimates of the spherical harmonic coefficients evaluated from Planck NPIPE data splits A and B, respectively.

The pipeline followed to end up with the final estimate of the CB power spectrum is the same described in this work, with the only exception on how the isotropic bias term, $C_L^{bias,iso}$ (see eq. (3.3)), has been evaluated. Specifically, the starting point for its evaluation is:

$$ C_L^{bias,iso} = \langle \hat{\alpha}_{LM}^{[JK]} \hat{\alpha}_{LM}^{[JK],*} \rangle, \qquad (C.3) $$

and, since the $\hat{\alpha}_{LM}^{[JK]}$ estimates have been obtained in terms of the estimates from each NPIPE data split (eq. (C.2)), we can express the isotropic bias term needed for the debiasing procedure for this test, in terms of the one evaluated from the $slpitA \times splitA$, $splitB \times splitB$ and $splitA \times splitB$ cases, as:

$$ C_L^{bias,iso} = C_L^{bias,iso,A \times A} - 2C_L^{bias,iso,A \times B} + C_L^{bias,iso,B \times B}. \qquad (C.4) $$

In the right panel of figure 14 we show the binned CB power spectrum obtained from the new set of spherical harmonic coefficients, $\hat{\alpha}_{LM}^{[JK]}$. Also for this case, the CB power spectrum is compatible with zero with a PTE of 29.83%.
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Jack-knife tests. Left panel: CB power spectrum estimated from a new set of maps, $m_{JK}$, obtained subtracting the CMB maps of splitB from the ones of splitA. Right panel: CB power spectrum estimated from a new set of spherical harmonic coefficients estimates, $\hat{\alpha}_{LM}^{[JK]}$, obtained as the difference of the estimates from splitA and splitB.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Left panel: correlation matrix for the first jack-knife test. Right panel: correlation matrix for the second jack-knife test. The correlation coefficients are expressed in terms of percentage.</figcaption></figure>

In figure 15 we report the correlation matrices associated to the first and the second jack-knife tests presented in this study, in the left and right panel, respectively. Even though by eye we do not expect a better PTE associated to the first jack-knife test (power spectrum showed in the left panel of figure 14), we can explain it observing its correlation matrix. In fact, it is possible to appreciate a higher degree of correlation.
References

[1] J.C. Mather et al., *A preliminary measurement of the Cosmic Microwave Background spectrum by the Cosmic Background Explorer (COBE) satellite*, Astrophys. J. Lett. 354 (1990) L37 [INSPIRE].

[2] W.C. Jones et al., *A measurement of the angular power spectrum of the CMB temperature anisotropy from the 2003 flight of BOOMERANG*, Astrophys. J. 647 (2006) 823 [astro-ph/0507494] [INSPIRE].

[3] S. Hanany et al., MAXIMA-1: a measurement of the cosmic microwave background anisotropy on angular scales of 10 arcminutes to 5 degrees, Astrophys. J. Lett. 545 (2000) L5 [astro-ph/0005123] [INSPIRE].

[4] WMAP collaboration, *Nine-year Wilkinson Microwave Anisotropy Probe (WMAP) observations: cosmological parameter results*, Astrophys. J. Suppl. 208 (2013) 19 [arXiv:1212.5226] [INSPIRE].

[5] PLANCK collaboration, Planck 2018 results. I. Overview and the cosmological legacy of Planck, Astron. Astrophys. 641 (2020) A1 [arXiv:1807.06205] [INSPIRE].

[6] ACT collaboration, *The Atacama Cosmology Telescope: DR4 maps and cosmological parameters*, JCAP 12 (2020) 047 [arXiv:2007.07288] [INSPIRE].

[7] SPT-3G collaboration, *Measurement of the CMB temperature power spectrum and constraints on cosmology from the SPT-3G 2018 TT, TE, and EE dataset*, Phys. Rev. D 108 (2023) 023510 [arXiv:2212.05642] [INSPIRE].

[8] BICEP and KECK collaborations, *Improved constraints on primordial gravitational waves using Planck, WMAP, and BICEP/Keck observations through the 2018 observing season*, Phys. Rev. Lett. 127 (2021) 151301 [arXiv:2110.00483] [INSPIRE].

[9] LITEBIRD collaboration, *Probing cosmic inflation with the LiteBIRD cosmic microwave background polarization survey*, PTEP 2023 (2023) 042F01 [arXiv:2202.02773] [INSPIRE].

[10] SIMONS OBSERVATORY collaboration, *The Simons observatory: science goals and forecasts*, JCAP 02 (2019) 056 [arXiv:1808.07445] [INSPIRE].

[11] K. Abazajian et al., CMB-S4 science case, reference design, and project plan, arXiv:1907.04473 [INSPIRE].

[12] S. Dodelson, *Modern cosmology*, Academic Press, Amsterdam, The Netherlands (2003) [INSPIRE].

[13] Y. Minami and E. Komatsu, *New extraction of the cosmic birefringence from the Planck 2018 polarization data*, Phys. Rev. Lett. 125 (2020) 221301 [arXiv:2011.11254] [INSPIRE].

[14] P. Diego-Palazuelos et al., *Cosmic birefringence from Planck public release 4*, in the proceedings of the 56th Rencontres de Moriond on cosmology, (2022) [arXiv:2203.04830] [INSPIRE].

[15] J.R. Eskilt, *Frequency-dependent constraints on cosmic birefringence from the LFI and HFI Planck data release 4*, Astron. Astrophys. 662 (2022) A10 [arXiv:2201.13347] [INSPIRE].

[16] J.R. Eskilt and E. Komatsu, *Improved constraints on cosmic birefringence from the WMAP and Planck cosmic microwave background polarization data*, Phys. Rev. D 106 (2022) 063503 [arXiv:2205.13962] [INSPIRE].

[17] S.M. Carroll, G.B. Field and R. Jackiw, *Limits on a Lorentz and parity violating modification of electrodynamics*, Phys. Rev. D 41 (1990) 1231 [INSPIRE].
[18] Y. Minami et al., *Simultaneous determination of the cosmic birefringence and miscalibrated polarization angles from CMB experiments*, PTEP 2019 (2019) 083E02 [arXiv:1904.12440] [INSPIRE].

[19] Y. Minami, *Determination of miscalibrated polarization angles from observed cosmic microwave background and foreground EB power spectra: application to partial-sky observation*, PTEP 2020 (2020) 063E01 [arXiv:2002.03572] [INSPIRE].

[20] Y. Minami and E. Komatsu, *Simultaneous determination of the cosmic birefringence and miscalibrated polarization angles II: including cross frequency spectra*, PTEP 2020 (2020) 103E02 [arXiv:2006.15982] [INSPIRE].

[21] L. Pagano et al., *CMB polarization systematics, cosmological birefringence and the gravitational waves background*, Phys. Rev. D 80 (2009) 043522 [arXiv:0905.1651] [INSPIRE].

[22] PLANCK collaboration, *Planck intermediate results. XLIX. Parity-violation constraints from polarization data*, Astron. Astrophys. 596 (2016) A110 [arXiv:1605.08633] [INSPIRE].

[23] D. Contreras, P. Boubel and D. Scott, *Constraints on direction-dependent cosmic birefringence from Planck polarization data*, JCAP 12 (2017) 046 [arXiv:1705.06387] [INSPIRE].

[24] SPT collaboration, *Searching for anisotropic cosmic birefringence with polarization data from SPTpol*, Phys. Rev. D 102 (2020) 083504 [arXiv:2006.08061] [INSPIRE].

[25] T. Namikawa et al., *Atacama Cosmology Telescope: constraints on cosmic birefringence*, Phys. Rev. D 101 (2020) 083527 [arXiv:2001.10465] [INSPIRE].

[26] A. Gruppuso, D. Molinari, P. Natoli and L. Pagano, *Planck 2018 constraints on anisotropic birefringence and its cross-correlation with CMB anisotropy*, JCAP 11 (2020) 066 [arXiv:2008.10334] [INSPIRE].

[27] M. Bortolami et al., *Planck constraints on cross-correlations between anisotropic cosmic birefringence and CMB polarization*, JCAP 09 (2022) 075 [arXiv:2206.01635] [INSPIRE].

[28] E. Komatsu, *New physics from the polarized light of the cosmic microwave background*, Nature Rev. Phys. 4 (2022) 452 [arXiv:2202.13919] [INSPIRE].

[29] S.E. Clark, C.-G. Kim, J.C. Hill and B.S. Hensley, *The origin of parity violation in polarized dust emission and implications for cosmic birefringence*, Astrophys. J. 919 (2021) 53 [arXiv:2105.00120] [INSPIRE].

[30] A.J. Cukierman, S.E. Clark and G. Halal, *Magnetic misalignment of interstellar dust filaments*, Astrophys. J. 946 (2023) 106 [arXiv:2208.07382] [INSPIRE].

[31] P. Diego-Palazuelos et al., *Robustness of cosmic birefringence measurement against galactic foreground emission and instrumental systematics*, JCAP 01 (2023) 044 [arXiv:2210.07655] [INSPIRE].

[32] S.-S. Chern and J. Simons, *Characteristic forms and geometric invariants*, Annals Math. 99 (1974) 48 [INSPIRE].

[33] G. Gubitosi et al., *A constraint on Planck-scale modifications to electrodynamics with CMB polarization data*, JCAP 08 (2009) 021 [arXiv:0904.3201] [INSPIRE].

[34] M. Lembo et al., *Cosmic microwave background polarization as a tool to constrain the optical properties of the universe*, Phys. Rev. Lett. 127 (2021) 011301 [arXiv:2007.08486] [INSPIRE].

[35] K. Murai, F. Naokawa, T. Namikawa and E. Komatsu, *Isotropic cosmic birefringence from early dark energy*, Phys. Rev. D 107 (2023) L041302 [arXiv:2209.07804] [INSPIRE].
[36] J.R. Eskilt et al., *Constraints on early dark energy from isotropic cosmic birefringence*, Phys. Rev. Lett. 131 (2023) 121001 [arXiv:2303.15369] [INSPIRE].

[37] M. Li and X. Zhang, *Cosmological CPT violating effect on CMB polarization*, Phys. Rev. D 78 (2008) 103516 [arXiv:0810.0403] [INSPIRE].

[38] M. Kamionkowski, *How to de-rotate the cosmic microwave background polarization*, Phys. Rev. Lett. 102 (2009) 111302 [arXiv:0810.1286] [INSPIRE].

[39] G. Gubitosi et al., *Using CMB data to constrain non-isotropic Planck-scale modifications to electrodynamics*, JCAP 11 (2011) 003 [arXiv:1106.6049] [INSPIRE].

[40] S.M. Carroll, *Quintessence and the rest of the world*, Phys. Rev. Lett. 81 (1998) 3067 [astro-ph/9806099] [INSPIRE].

[41] S. Panda, Y. Sumitomo and S.P. Trivedi, *Axions as quintessence in string theory*, Phys. Rev. D 83 (2011) 083506 [arXiv:1011.5877] [INSPIRE].

[42] D.J.E. Marsh, *Axion cosmology*, Phys. Rept. 643 (2016) 1 [arXiv:1510.07633] [INSPIRE].

[43] E.G.M. Ferreira, *Ultra-light dark matter*, Astron. Astrophys. Rev. 29 (2021) 7 [arXiv:2005.03254] [INSPIRE].

[44] L.M. Capparelli, R.R. Caldwell and A. Melchiorri, *Cosmic birefringence test of the Hubble tension*, Phys. Rev. D 101 (2020) 123529 [arXiv:1909.04621] [INSPIRE].

[45] E. Di Valentino et al., *In the realm of the Hubble tension — a review of solutions*, Class. Quant. Grav. 38 (2021) 153001 [arXiv:2103.01183] [INSPIRE].

[46] F. Finelli and M. Galaverni, *Rotation of linear polarization plane and circular polarization from cosmological pseudo-scalar fields*, Phys. Rev. D 79 (2009) 063002 [arXiv:0802.4210] [INSPIRE].

[47] M.A. Fedderke, P.W. Graham and S. Rajendran, *Axion dark matter detection with CMB polarization*, Phys. Rev. D 100 (2019) 015040 [arXiv:1903.02666] [INSPIRE].

[48] H. Nakatsuka, T. Namikawa and E. Komatsu, *Is cosmic birefringence due to dark energy or dark matter? A tomographic approach*, Phys. Rev. D 105 (2022) 123509 [arXiv:2203.08560] [INSPIRE].

[49] M. Galaverni, F. Finelli and D. Paoletti, *Redshift evolution of cosmic birefringence in CMB anisotropies*, Phys. Rev. D 107 (2023) 083529 [arXiv:2301.07971] [INSPIRE].

[50] A. Greco, N. Bartolo and A. Gruppuso, *Probing axions through tomography of anisotropic cosmic birefringence*, JCAP 05 (2023) 026 [arXiv:2211.06380] [INSPIRE].

[51] A. Greco, N. Bartolo and A. Gruppuso, *A new solution for the observed isotropic cosmic birefringence angle and its implications for the anisotropic counterpart through a Boltzmann approach*, arXiv:2401.07079 [INSPIRE].

[52] R.R. Caldwell, V. Gluscevic and M. Kamionkowski, *Cross-correlation of cosmological birefringence with CMB temperature*, Phys. Rev. D 84 (2011) 043504 [arXiv:1104.1634] [INSPIRE].

[53] V. Gluscevic, M. Kamionkowski and A. Cooray, *De-rotation of the cosmic microwave background polarization: full-sky formalism*, Phys. Rev. D 80 (2009) 023510 [arXiv:0905.1687] [INSPIRE].

[54] PLANCK collaboration, *Planck intermediate results. LVII. Joint Planck LFI and HFI data processing*, Astron. Astrophys. 643 (2020) A42 [arXiv:2007.04997] [INSPIRE].
[55] A.R. Pullen and M. Kamionkowski, Cosmic microwave background statistics for a direction-dependent primordial power spectrum, Phys. Rev. D 76 (2007) 103529 [arXiv:0709.1144] [INSPIRE].

[56] V. Gluscevic, D. Hanson, M. Kamionkowski and C.M. Hirata, First CMB constraints on direction-dependent cosmological birefringence from WMAP-7, Phys. Rev. D 86 (2012) 103529 [arXiv:1206.5546] [INSPIRE].

[57] V. Gluscevic, M. Kamionkowski and A. Cooray, De-rotation of the cosmic microwave background polarization: full-sky formalism, Phys. Rev. D 80 (2009) 023510 [arXiv:0905.1687] [INSPIRE].

[58] K.M. Górski et al., HEALPix — a framework for high resolution discretization, and fast analysis of data distributed on the sphere, Astrophys. J. 622 (2005) 759 [astro-ph/0409513] [INSPIRE].

[59] PLANCK collaboration, Planck 2018 results. IV. Diffuse component separation, Astron. Astrophys. 641 (2020) A4 [arXiv:1807.06208] [INSPIRE].

[60] PLANCK collaboration, Planck 2013 results. VIII. HFI photometric calibration and mapmaking, Astron. Astrophys. 571 (2014) A8 [arXiv:1303.5069] [INSPIRE].

[61] LSST DARK ENERGY SCIENCE collaboration, A unified pseudo-$C_\ell$ framework, Mon. Not. Roy. Astron. Soc. 484 (2019) 4127 [arXiv:1809.09603] [INSPIRE].

[62] N. Zhu et al., The Simons observatory large aperture telescope receiver, Astrophys. J. Supp. 256 (2021) 23 [arXiv:2103.02747] [INSPIRE].

[63] POLARBEAR collaboration, POLARBEAR constraints on cosmic birefringence and primordial magnetic fields, Phys. Rev. D 92 (2015) 123509 [arXiv:1509.02461] [INSPIRE].

[64] A. Greco, N. Bartolo and A. Gruppuso, Cosmic birefringence: cross-spectra and cross-bispectra with CMB anisotropies, JCAP 03 (2022) 050 [arXiv:2202.04584] [INSPIRE].

[65] C.R. Harris et al., Array programming with NumPy, Nature 585 (2020) 357 [arXiv:2006.10256] [INSPIRE].

[66] J.D. Hunter, Matplotlib: a 2D graphics environment, Comput. Sci. Eng. 9 (2007) 90 [INSPIRE].

[67] A. Zonca et al., healpy: equal area pixelization and spherical harmonics transforms for data on the sphere in python, J. Open Source Softw. 4 (2019) 1298 [INSPIRE].

[68] K.M. Górski et al., HEALPix — a framework for high resolution discretization, and fast analysis of data distributed on the sphere, Astrophys. J. 622 (2005) 759 [astro-ph/0409513] [INSPIRE].