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

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σ 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.


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][2][3][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][7][8][9][10][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 where F µν 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][14][15][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][19][20] able to disentangle the instrumental polarization angle from the CB effect1 , hint at a detection of a CB angle β ∼ 0.3 • at the level of 2.5 -3σ.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][24][25][26][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][30][31], these analyses hinting at an isotropic CB effect, suggest the need to extend the electromagnetic sector of the standard model L SM em with a parity-violating term L CS known as a Chern-Simons term [32] 2 .In Eq. (1.2), λ/f is a coupling with the dimension of the inverse of an energy, ϕ is a new scalar (or pseudo-scalar) field and F µν 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 ϕ is taken to be homogeneous [35,36].Fluctuations of ϕ around its homogeneous part, are instead able to produce anisotropic CB [37,38].Hence, the CB effect can be seen as a tracer of the existence of a new cosmological field ϕ (typically referred as an axion) acting as dark matter or dark energy [39][40][41][42] which might also play a role in alleviating the Hubble tension, see e.g.[43,44] and reference therein.See also [45][46][47][48][49].
Other works, e.g.[36,50], have constrained the axion's parameters through the CMB EB 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.[51].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 an harmonic estimator based on the approach presented in [52] and [25].The aim is to apply this estimator to the most recent Planck data, namely PR3 [5] and PR4 (also known as NPIPE) [53], 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).
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.

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 X ℓm , instead we observe their sum with the rotation-induced modes, δa X ℓm .Here X = E, B for E-and B-modes respectively: 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 4 .The expressions for the rotation-induced E-modes and B-modes, as derived in [38], are given as: ) where equation (2.3) is different from zero for ℓ + ℓ ′ + L even, equation (2.4) is different from zero for ℓ+ℓ ′ +L odd, α LM are the spherical harmonic coefficients of the CB field, and ξ LM ℓmℓ ′ m ′ and H L ℓℓ ′ are defined in terms of Wigner-3j symbols as follows: The rotation-induced modes generate correlations between ℓ − ℓ ′ pairs with ℓ ̸ = ℓ ′ , leading to, as anticipated in section 1, parity-violating cross-correlations.These non-standard crosscorrelations caused by CB can be characterized by the following general structure: where X = {T, E, B} and contains the information about the primordial spectra before the rotation (see table 1 of [52]).Another way to write equation (2.7) is in terms of the rotational invariants (i.e.quantities independent of m) [38], where, (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 ℓ , so that [52]: 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): DLM,XX ′ ,map and the other is the inverse variance weighting average estimator from [54]: where G ℓℓ ′ is defined as: The first expression of the estimates for the α LM coefficients for a fixed ℓ − ℓ ′ pair is obtained inverting equation (2.11) and substituting the expression of the estimator for the rotational invariants of equation (2.12): where we defined 38].Note that we are not using the over-hat symbol since, before ending up with the final estimates of the α 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 (α LM ) XX ′ ℓℓ ′ defined in equation (2.14), ending up with: where is the analytic variance associated to each (α LM ) ℓℓ ′ 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): where A represents the considered cross-correlation, i.e.XX ′ , and C ℓℓ ′ AA ′ is the entry AA ′ of the covariance matrix of the rotational invariants (see equation 33 of [52] 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 fullsky 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 α LM coefficients, we followed 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 α 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 α 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: where the α 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 α LM coefficients have been obtained as: Next, we will focus specifically on the EB cross-correlation.

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: where the last equivalence holds since the term δa E ℓm δa B, * ℓ ′ m ′ 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 ℓ + ℓ ′ + L is even (see equation (2.3)).Equation (2.7) for the specific case of the EB cross-correlation now reads: where C EE ℓ refers to the Z XX ′ ℓℓ ′ term (eq.(2.7)) when XX ′ = EB.Following the logical steps previously described, the first expression of the estimator for the α LM coefficients before the subtraction of the mean field bias (eq.(2.15)) is calculated as: and the expression for the inverse variance of the estimator for the EB cross-correlation is: Focusing only on the numerator in equation (2.21) (meaning the unnormalized estimator, UN ): and exploiting properties of Wigner-3j symbols, it is possible to re-write equation (2.23) as5 : where ±2 Y ℓm are spin-2 spherical harmonics and C XX,map ℓ indicates the analytical expression for the power spectrum recovered from CMB maps, evaluated as: with C XX ℓ the cosmological signal, W ℓ the window function, and N 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)).
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: with a E, * ℓm and a B, * ℓm defined as: By making use of equations (2.26) and (2.27), the unnormalized estimator for the α LM coefficients can be written as: At this point we define the quantity inside the square brackets as a "map-like" object: and, after performing the complex conjugate of equation (2.30), we obtain: where the second equivalence holds since maps are real objects, i.e. m ′ * (α) = m ′ (α).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 ′ (α) 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 α LM estimates, thus the estimates of the α LM coefficients before the de-biasing procedure are obtained as: To end up with the final estimates of the spherical harmonic coefficients of the CB field, i.e. αLM , we have to subtract the mean field bias (eq.(2.17)) from equation (2.33), as presented in equation (2.18).

Angular power spectrum and de-biasing procedure
Having the estimates of the α LM coefficients, the αα power spectrum can be written as: where f sky is the sky fraction of the mask used for the analysis.The estimate of the CB power spectrum in equation The bias comes from the diagonal contribution of ℓ − ℓ ′ pairs with ℓ = ℓ ′ when combining together the two estimates of the α LM coefficients and from the off-diagonal contribution from sources different from the rotation induced by CB.We thus encode for a debiasing procedure [55], now at the level of the power spectrum, to, first, select the diagonal contribution (i.e. the contribution from ℓ − ℓ ′ pairs with ℓ = ℓ ′ ) and, second, among the off-diagonal contributions, select the contribution coming only from the rotation induced by CB: The bias term accounts for two different contributions: • the isotropic bias term, C bias,iso L , which is an analytic bias term calculated on data, so that it selects the diagonal contribution coming from the ℓ − ℓ ′ pairs with ℓ = ℓ ′ ; • the Monte Carlo bias term, C bias,M C L .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 ( [56]).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 debiasing 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: From a general point of view, the estimator (eq.(2.23)) probes three disconnected Wick contractions [55]: 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: where Ĉ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 XX,map ℓ (without the overhat 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 bias,M C L 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: 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)): This methodology has been validated and the results of this validation are reported in appendix B. In the following sections we present the data set used for the analysis and the application of the aforementioned pipeline to Planck data.

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) [53] data products of the Planck satellite, which contain CMB data and simulation maps at the HEALPix6 [57] resolution of N side = 2048.The CMB maps employed for the main results have been cleaned using the official Planck component separation method Commander [58].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 tempera-ture+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 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 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 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.

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

Mask standard
Galactic mask: 20% Galactic mask: 40% Galactic mask: 60% 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 (ℓ CM B min = 50) and the maximum multipole included in the analysis is ℓ CM B max = 2000.The CB power spectrum has been evaluated from L CB min = 0 up to L CB max = 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 α LM coefficients estimated from simulations only; • subtract the mean field bias from all the α LM estimates, both the ones from data and that from simulations; • calculate the C α α L (eq.(3.1)) and the C bias,iso L (eq.(3.4)) for both data and simulations of Planck data products; • split the C α α L and the C bias,iso L 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, Ĉαα L , subtracting the isotropic bias term, C bias,iso L evaluated from the same set, and the Monte Carlo bias term, C bias,M C L , calculated from set A; • subtract the isotropic bias term evaluated from data, C bias,iso L , and the Monte Carlo bias term evaluated from set A, C bias,M C L , from the biased αα power spectrum calculated on Planck data, in order to obtain the estimated CB power spectrum, Ĉαα L (eq. (3.2)).

Full mission
We present the application of our pipeline in case of a CB power spectrum obtained combining together α LM coefficients estimated from Planck NPIPE full mission data.
The equation for the αα power spectrum before the de-bias is the same of equation (3.1),where f 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 sky = 0.78.] Power spectra evaluated from full mission NPIPE data and simulations: αα 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); αα power spectrum after the de-biasing procedure (gray curve).
In figure 2 we plot the power spectra of the different terms.The black and green curves represent the biased αα power spectrum and the isotropic bias term, respectively, evaluated using the α LM estimates from Planck NPIPE data, while the orange curve represents the Monte Carlo bias term, computed by averaging the de-biased αα 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 αα power spectrum in figure 3, we display it with the gray curve also in figure 2. The black curve in the upper panel of figure 3 is the αα power spectrum evaluated for each multipole after the de-biasing procedure (eq.(3.2)).The shaded areas represent the 1σ, 2σ and 3σ confidence intervals, obtained after the computation of the variance of the fully de-biased αα 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σ.
In figure 2 the rotation signal induced by CB before the de-biasing procedure, i.e.C α α L , is different from zero.We can understand the reason of this strong deviation from zero since the single estimate for the α LM coefficients goes as: and when combining two estimates to obtain the power spectrum before the de-bias, we end up with: The only non negligible contributions come from: (5.4) The observed EE and BB power spectra encode for two contributions; the cosmological signal, i.e.C EE ℓ and C BB ℓ , and the noise contribution, i.e N EE ℓ and N BB ℓ , so that: (5.6) In case of a X,map ℓm 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 Fig. 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%.
In the left and right panels of figure 4 we show the de-biased αα 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 αα 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 αα 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.

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 α LM coefficients estimated from data split A (half mission 1 ) of Planck NPIPE (PR3) data products; • auto-correlating α LM coefficients estimated from data split B (half mission 2 ) of Planck NPIPE (PR3) data products;   • cross-correlating α 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 αα 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 α 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 L,EB ℓℓ ′ , F L,BE ℓℓ ′ and the analytic power spectra, C XX,map ℓ , must be evaluated separately for each data split (or each half-mission).Thus, the biased power spectrum of equation (3.1) is now calculated as: (5.7) And, since we are combining together estimates coming from two different data sets, the normalization of the α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 (halfmission 1 ) as (σ −2 L ) [A] ((σ −2 L ) [1] ) and as (σ −2 L ) [B] ((σ −2 L ) [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: 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σ.More precisely, the presented results show the αα power spectra obtained auto-correlating α LM estimates from split A (blue), autocorrelating estimates from split B (orange) and cross-correlating α 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 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. ] 71.01%α [1] α [1] 75.95% 7.77% α [2] α [2] 50.48%B] 83.08% α [1] α [2] 31.88% αα 84.35% αα 6.61% 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.
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.

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 ℓ CM B min = 10, ℓ CM B min = 30, ℓ CM B min = 50 and ℓ CM B min = 100; • for three different choices of the maximum CMB multipole included in the analysis, encoding for ℓ CM B max = 1500, ℓ CM B max = 2000 and ℓ CM B max = 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%.In blue we show the power spectrum obtained with the standard mask, corresponding to a f sky = 78.0%; in orange, green and red we show the power spectrum obtained applying masks corresponding to f sky = 59.2%, 39.7%, 20.3% respectively.
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 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 sky = 20%.

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 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 α 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 [59].For the analysis of the cross-correlation between CB and CMB temperature and polarization fields, we employed the Pymaster Python package, 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 (α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 (αE and α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 αT power spectrum is not represented in the units used for αE and αB, rather it is showed in band-powers to allow for a more direct comparison with the αT power spectrum presented in [24].

map -Planck NPIPE commander
Table 2 presents the Probability To Exceed (PTE) values for the cross-correlation power spectra between the CB and CMB fields.Notably, the α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%).

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 α LM coefficients, but also serves as a rapid tool for forecasting future CMB experiments.By assessing σ 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 σ 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 [60].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 σ 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.

Conclusions
In this study we have built an estimator for the spherical harmonic coefficients of the CB field which implements the method described in [52] 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 splits.In all cases, our analysis consistently found that the CB power spectrum is compatible with zero at a significance level of approximately 2σ.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 [61],  3. 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 αT , αE, and α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.[49,62].
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 ∼ 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 ∼ 1000.
The code and the pipeline developed for this work are publicly available7 , along with products employed 8 .Upon requests we can provide additional products such as birefringence maps or α LM coefficients.

A Construction of the estimator
In section 2.1 we re-write the part of equation (2.33) without its normalization (i.e., without σ −2 L ), that we indicate as α U N LM , in order to reduce the computation time.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 ℓ − ℓ ′ inverted, an equivalent expression of (A.1), making the expression of the F L,EB ℓℓ ′ term explicit, is: and, exploiting the definition of H L ℓℓ ′ and the properties of the Wigner-3j symbols, recalling that ℓ + L + ℓ ′ must be even in the case of the EB cross-correlation induced by CB, it is also possible to write: , so that we can use the definition of the triple integral to re-write the product (A.3)as: (2ℓ + 1)(2L + 1)(2ℓ ′ + 1) 4π Therefore, the estimator can be written as: where the second equivalence follows from (a X ℓm ) * = (−1) m a X ℓ−m .At this point we define two new objects: 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 αα 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 αα power spectrum before the de-biasing procedure, C α α L , and the isotropic bias term, C bias,iso L .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.average of the estimated αα power spectra (black curve), the 1σ, 2σ and 3σ confidence intervals are the dark green, light green and gray areas, respectively.Lower panel Difference between the input and the recovered power spectrum.
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σ.

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 αα power spectra computed over the simulations (black curve) with the expected zero rotation signal (red curve).The shaded areas are

Figure 1 :
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.

L [deg 2 ]Figure 3 :
Figure 3: Results evaluated from full mission NPIPE data.Upper panel De-biased αα power spectrum (black line), with the 1σ (dark green area), 2σ (light green area) and 3σ (gray area) confidence levels.Lower panel De-biased αα power spectrum after binning with 100 multipoles per bin and excluding the first 8 multipoles.

Figure 4 :
Figure4: Left panel De-biased αα 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.

Figure 5 :
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.

Figure 6 :
Figure6: Left panel De-biased αα power spectrum after binning with 100 multipoles per bin, excluding the first 8 multipoles from the analysis, obtained auto-correlating α 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 α 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.

Figure 7 :Figure 8 :
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 ℓ CM B min included in the analysis and ℓ CM B max = 2000 for all cases.We indicate in blue the power spectrum obtained with ℓ CM B min = 10, in orange the one with ℓ CM B min = 30, in green the one with ℓ CM B min = 50 and in red the one with ℓ CM B min = 100.Right panel Power spectra for different values of ℓ CM B max included in the analysis and ℓ CM B min = 50 for all cases.We indicate in blue, orange and green the CB power spectra obtained including CMB multipoles up to ℓ CM B max = 1500, ℓ CM B max = 2000 and ℓ CM B max = 2500, respectively.

Figure 9 :
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.

Figure 10 :
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 αT cross-spectrum in band-powers; Lower left panel αE cross-spectrum; Lower right panel αB cross-spectrum.

Figure 11 :
Figure 11: Variance of the estimator evaluated with the instrumental specifications of the four CMB experiments reported in table3.

Figure 12 :
Figure12: Validation with rotation signal: Upper panel Input power spectrum (red curve), average of the estimated αα power spectra (black curve), the 1σ, 2σ and 3σ confidence intervals are the dark green, light green and gray areas, respectively.Lower panel Difference between the input and the recovered power spectrum.
biased αα 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.

Table 2 :
Probability To Exceed for the cross-spectra of the CB field with the temperature and polarization CMB fields.

Table 3 :
. Instrumental specifications of the considered CMB experiments.