© J. R. Eskilt 2022

1. Introduction

Parity violation has so far only been observed in the weak interaction, but there are extensions of the standard model that introduce more parity-violating effects. One popular example is a pseudo-scalar axion-like field ϕ which couples to the electromagnetic tensor and can cause parity-violating physics for electromagnetic waves (Turner & Widrow 1988; Carroll et al. 1990; Harari & Sikivie 1992). The Chern-Simons term couples the field ϕ to the electromagnetic tensor F_μν of the form $\mathcal{L}\supset \frac{1}{4}{g}_{\varphi \gamma }\varphi F_{\mu \nu }{\stackrel{\sim }{F}}^{\mu \nu }$, where ${\stackrel{\sim }{F}}_{\mu \nu }$ is the dual tensor of F_μν and g_ϕγ is the coupling constant (see the review by Marsh 2016). It is possible to show that this term changes the phase velocity of the left- and right-handed circular polarization of photons, which effectively causes a rotation of the linear polarization of electromagnetic waves by an angle β. This can be detected in the polarization of the cosmic microwave background (CMB). The standard model of cosmology, ΛCDM, predicts no intrinsic EB correlation of the CMB.

A cosmic birefringence angle β would rotate the intrinsic EE power spectrum, $C_{\mathcal{l}}^{EE}$, and BB power spectrum, $C_{\mathcal{l}}^{BB}$, into an observed EB power spectrum, $C_{\ell }^{EB,o}=\frac{sin(4\beta )}{2}(C_{\ell }^{\mathit{EE}}-C_{\ell }^{\mathit{BB}})$, where “o” denotes the observed value and ℓ is the multipole moment. However, the problem of measuring β this way is two-fold: there is polarized Galactic emission that conceals the cosmic signal and no instruments are perfectly calibrated. The polarization miscalibration angle of an instrument, α, is perfectly degenerate with the cosmic birefringence angle, β, when solely the CMB is analyzed. A lack of strong priors on α for an instrument will necessarily give a large degree of uncertainty on β. This is seen in the original Planck 2015 analysis, where the authors found β = 0.29 ± 0.05(stat.)±0.28(syst.) (Planck Collaboration Int. XLIX 2016). Here, the latter systematic uncertainty is the miscalibration uncertainty expected from the calibration before launch. The other problem is that the Galactic foreground emission is polarized; without a proper foreground removal, the polarized foreground emission will inevitably contaminate the measurements.

A novel method to mitigate these two problems was developed by Minami et al. (2019) and Minami & Komatsu (2020a). The assumption is that the mechanism that can cause cosmic birefringence mostly affects the CMB photons. This could either be because the effect was the strongest shortly after the last-scattering surface or because the effect is stronger over longer distances. This means that cosmic birefringence has a negligible effect on the Galactic foreground emission and that the linearly polarized CMB emission is rotated by α + β, while the polarized foreground emission is only rotated by α. This method allows for the degeneracy between α and β to be broken, but the intrinsic foreground EB power spectrum needed to be specified. Assuming that the intrinsic EB correlation of the foreground is zero, Minami & Komatsu (2020b) used the high-frequency instrument (HFI) polarized maps from the Planck public release 3 (PR3) to find a measurement of β = 0.34° ±0.14° (68% C.L.), which corresponds to the statistical significance of 2.4σ. They emphasized the possibility that the intrinsic foreground EB can create an effective angle γ which skews the measurements β → β − γ and α → α + γ; this leaves the sum α + β invariant. The authors pointed out that the positive TB and TE inferred from the Planck data suggest EB > 0 (Huffenberger et al. 2020); hence, γ > 0, which would make the measured β a lower bound.

An astrophysical model for the intrinsic EB of dust was described by Clark et al. (2021). Their work showed evidence for the applied Galactic mask to greatly affect the magnitude and sign of the EB power spectra of the polarized foreground emission. With a Galactic mask of sky fraction f_sky ≈ 0.7, they found a robustly positive $C_{\mathcal{l}}^{EB,\text{dust}}$ for ℓ ≲ 500 using their filamentary dust model, which would also suggest a lower bound of β if one ignores the intrinsic dust EB. On the other hand, a smaller Galactic mask gave EB correlations with a fluctuating sign over multipoles ℓ. This mask dependence was confirmed by Diego-Palazuelos et al. (2022) where a declining measurement of β was found when ignoring the intrinsic EB power spectra of the foreground using polarized data from the Planck public release 4 (PR4; Planck Collaboration Int. LVII 2020), often called the “NPIPE” release.

In Minami & Komatsu (2020b) and Diego-Palazuelos et al. (2022), the cosmic birefringence angle β was assumed to be independent of the photon frequency ν. In this paper, we explore a possible frequency dependence. This is not the first time a frequency-dependent cosmic birefringence angle has been considered in the literature. Gubitosi et al. (2014) used polarization data of WMAP7, BOOMERanG, QUAD, and BICEP to look for any hints of a frequency-dependent β. Other constraints have come from looking at Lorentz violating origins, such as Kahniashvili et al. (2008) using WMAP data. Galaverni et al. (2015) also constrained several models of cosmic birefringence using CMB and other astrophysical data at a much wider frequency range. No detection of a frequency-dependent β has been found, but this work generalizes the methods of Diego-Palazuelos et al. (2022) and Minami & Komatsu (2020b) to search for such a frequency dependence. These authors assumed that β was constant for all frequencies, which gave them a hint of a non-zero β. Motivated by their result, the goal of this work is to learn more about the frequency dependence of that signal. Also, for the first time, we include the low-frequency instrument (LFI) polarized data to put constraints on the frequency dependence of the signal. All uncertainties quoted in this paper are at a 68% confidence level.

2. Theory

In this section, the focus is mainly on an axion-like field that couples to electromagnetism. We also go on to show how we measure the cosmic birefringence signal to be consistent with frequency independence, which we go on to show is consistent with the predictions of an axion-like field. Therefore, we do not go into the details of Faraday rotation or any other models that would cause a frequency-dependent signal.

In this paper, we assume that only the linearly polarized CMB emission has been rotated by the cosmic birefringence angle, β, while the polarized foreground emission from the Milky Way has not been rotated. We are, therefore, looking for a potential mechanism that is either substantial on large length scales so that its effect on the polarized foreground emission is negligible or it is, instead, an effect that was only present around the time of recombination or shortly after.

An axion-like field ϕ that couples to electromagnetism can be the cause of cosmic birefringence. We make the assumption that the field is homogeneous and only varies over time. The dispersion relation for electromagnetic waves becomes (Carroll et al. 1990; Harari & Sikivie 1992)

$$\begin{array}{c}\hfill {\omega }_{\pm }^2=k^2\pm k{g}_{\varphi \gamma }\dot{\varphi },\end{array}$$(1)

where ω is the angular frequency, k is the wavenumber, and g_ϕγ is the coupling constant; in addition, + is right-handed circular polarization and − is left-handed. Performing a series expansion of ω_± in small $g_{\varphi \gamma }\dot{\varphi }$ gives:

$$\begin{array}{c}\hfill {\omega }_{\pm }=k\pm \frac{g_{\varphi \gamma }\dot{\varphi }}{2}-\frac{{\left(g_{\varphi \gamma }\dot{\varphi }\right)}^2}{8k}\pm \frac{{\left(g_{\varphi \gamma }\dot{\varphi }\right)}^3}{16k^2}+\mathcal{O}({\left(g_{\varphi \gamma }\dot{\varphi }\right)}^4/k^3).\end{array}$$(2)

Our convention for the cosmic birefringence angle β is that it rotates the linear polarization clockwise on the sky. We then get

$$\begin{array}{cc}\hfill \beta & =-{\displaystyle \int }\mathrm{d}t\frac{{\omega }_{+}-{\omega }_{-}}{2}\hfill \\ \hfill & =-{\displaystyle \int }\mathrm{d}t[\frac{g_{\varphi \gamma }\dot{\varphi }}{2}+\frac{{\left(g_{\varphi \gamma }\dot{\varphi }\right)}^3}{16k^2}+\mathcal{O}({\left(g_{\varphi \gamma }\dot{\varphi }\right)}^5/k^4)].\hfill \end{array}$$(3)

The first term is independent of the wavenumber. Neglecting any electromagnetic contribution, the equation of motion for ϕ is given by the Klein-Gordon equation (Marsh 2016),

$$\begin{array}{c}\hfill \varphi +3H\dot{\varphi }+m^2\varphi =0,\end{array}$$(4)

where m is the mass of the axion-like particle, while H is the Hubble parameter. There is only a finite mass interval of m where we can detect an axion-like particle with our method. If the mass is smaller than the Hubble expansion today, H₀ ≈ 10⁻³³ eV, then the field is generally too slowly varying by the Hubble friction to be detected. We note that a large coupling constant, g_ϕγ, can make fields with smaller masses detectable (Fujita et al. 2021). On the other hand, to get a significant rotation from β, we need the field to start oscillating after decoupling which yields m ≲ 10⁻²⁸ eV (Marsh 2016; Arvanitaki et al. 2010). Fields that oscillate before the decoupling quickly get suppressed since the era of decoupling was not instantaneous, meaning that not all polarized waves were rotated equally. The method behind this analysis is therefore able to probe ultra-light axion-like fields with a mass range of 10⁻³³ eV ≲ m ≲ 10⁻²⁸ eV and possibly even lower, depending on the amplitude of g_ϕγ.

If the mass is much larger than H₀, the axion-like field has decayed, ϕ(t₀)≈0, leading the cosmic birefringence angle to be β = g_ϕγϕ_rec/2, where ϕ_rec is the value of the field at recombination that is expected to have a negligible difference to its primordial value. A detection of cosmic birefringence in polarized CMB data can therefore constrain an axion-like field in the (g_ϕγ, ϕ_rec)-plane (Marsh 2016).

We go on to derive a crude upper bound on the first frequency-dependent term of β in Eq. (3) to show that the frequency dependence of β cannot be detected by the sensitivity level of Planck. We assume that the Hubble parameter is constant H = H_rec ≈ 10⁻²⁸ eV, and ϕ is critically damped, namely, m = 3H/2. The field is then ϕ = ϕ_rece^−3H t/2, where we take t = 0 at the time of recombination and approximate today to be at infinite time t → ∞. The contribution from the second term in Eq. (3) is:

$$\begin{array}{c}\hfill {\displaystyle |\int \mathrm{d}t\frac{{\left(g_{\varphi \gamma }\dot{\varphi }\right)}^3}{16k^2}|=\frac{3|g_{\varphi \gamma }{\varphi }_{\mathrm{rec}}{|}^3}{256{\pi }^2}{\left(\frac{H_{\mathrm{rec}}}{(1+z)\nu }\right)}^2={(1.7\times {10}^{-60})}^{°},}\end{array}$$(5)

where the frequency is ν = k/2π, and we used the lowest Planck frequency band that we consider in this work, namely, ν = 30 GHz; also, we assume that the redshift at recombination is z ≈ 1100, so that the initial frequency at recombination is ν_i = (1 + z)ν. We also use the result of Minami & Komatsu (2020b) that $\beta =-\int \mathrm{d}t\frac{g_{\varphi \gamma }\dot{\varphi }}{2}=\frac{g_{\varphi \gamma }{\varphi }_{\mathrm{rec}}}{2}\approx 0.3°$. The Hubble parameter, H, is certainly not constant and it will decrease over time, causing the field ϕ to become underdamped. In fact, Eq. (4) can be solved exactly if we assume the expansion parameter a ∝ tⁿ for some value of n. The resulting equation takes the form of Bessel functions (Marsh 2016). The point is however that changing the mass m and allowing H to vary with time will not magnify the first frequency-dependent term by around 60 orders of magnitude. It is, therefore, safe to assume that Planck’s sensitivity level will not pick up any potential frequency dependence of β caused by an ultra-light axion-like field coupled to electromagnetism.

There are, however, other mechanisms that could cause cosmic birefringence that we mention here in brief (Gubitosi & Paci 2013). If the observed isotropic cosmic birefringence was caused by Faraday rotation originating from primordial magnetic fields (Subramanian 2016), we would expect β_ν ∝ ν⁻². In the case where the birefringence angle is linear in frequency, β_ν ∝ ν, we could be observing Lorentz violating electrodynamics (Shore 2005). Such a measurement could be translated into a test of Lorentz violation (Kahniashvili et al. 2008). We also consider the possibility of β_ν ∝ ν², which has been predicted from quantum gravity theories to potentially modify the dispersion relation of photons (Myers & Pospelov 2003). We are therefore able to learn a lot about the origin of the measured cosmic birefringence result of Diego-Palazuelos et al. (2022) by quantifying the frequency dependence of the signal.

3. Method

Our analysis method follows that of Minami et al. (2019) and Minami & Komatsu (2020a). We start with the main assumption that the effect of cosmic birefringence, β, only substantially affects the CMB photons. However, an instrumental polarization miscalibration angle alters all linearly polarized plane waves, which we characterize by an instrument-dependent angle α. For a single detector, we can therefore write the spherical harmonics coefficients for the E- and B-modes (Seljak & Zaldarriaga 1997; Kamionkowski et al. 1997) as:

$$\begin{array}{cc}\hfill E_{\ell m}^o=& cos(2\alpha )E_{\ell m}^{\mathrm{fg}}-sin(2\alpha )B_{\ell m}^{\mathrm{fg}}\hfill \\ \hfill & +cos(2\alpha +2\beta )E_{\ell m}^{\mathrm{CMB}}-sin(2\alpha +2\beta )B_{\ell m}^{\mathrm{CMB}}+E_{\ell m}^{\mathrm{noise}},\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill B_{\ell m}^o=& cos(2\alpha )B_{\ell m}^{\mathrm{fg}}+sin(2\alpha )E_{\ell m}^{\mathrm{fg}}\hfill \\ \hfill & +cos(2\alpha +2\beta )B_{\ell m}^{\mathrm{CMB}}+sin(2\alpha +2\beta )E_{\ell m}^{\mathrm{CMB}}+B_{\ell m}^{\mathrm{noise}}.\hfill \end{array}$$(7)

Here, “fg” denotes the foreground component. We define the power spectra as $C_{\ell }^{\mathit{XY}}=\frac{1}{2\ell +1}\sum _m{X}_{\ell m}{Y}_{\ell m}^{\ast }$. To account for the effect of noise, we also define the ensemble average of the power spectra as $\langle C_{\mathcal{l}}^{XY}\rangle ={\delta }_{mm}{\delta }_{\mathcal{l}{\mathcal{l}}^{\prime }}\langle X_{\mathcal{l}m}{Y}_{{\mathcal{l}}^{\prime }{m}^{\prime }}^{*}\rangle$. Taking the equations above, we can derive the following equation (Minami et al. 2019):

$$\begin{array}{cc}\hfill ⟨C_{\ell }^{EB,o}⟩=& \frac{tan(4\alpha )}{2}(⟨C_{\ell }^{EE,o}⟩-⟨C_{\ell }^{BB,o}⟩)\hfill \\ \hfill & +\frac{sin(4\beta )}{2cos(4\alpha )}(⟨C_{\ell }^{EE,\mathrm{CMB}}⟩-⟨C_{\ell }^{BB,\mathrm{CMB}}⟩)\hfill \\ \hfill & +\frac{1}{cos(4\alpha )}⟨C_{\ell }^{EB,\mathrm{fg}}⟩+\frac{cos(4\beta )}{cos(4\alpha )}⟨C_{\ell }^{EB,\mathrm{CMB}}⟩,\hfill \end{array}$$(8)

where “o” denotes the observed value. We have thus been able to derive an equation that does not explicitly depend on the E- and B-modes of the polarized foregrounds, $\langle C_{\mathcal{l}}^{EE,\text{fg}}\rangle$ and $\langle C_{\mathcal{l}}^{BB,\text{fg}}\rangle$. And we can use the foreground to our advantage to break the degeneracy between the miscalibration angle α and the cosmic birefringence angle β. Unfortunately, we are left with the intrinsic EB-mode of the polarized foreground emission in our equation. For lack of an astrophysical model, this was set to zero in the initial analysis of Planck Data Release 3 in Minami & Komatsu (2020b), but they gave a discussion of how a non-zero EB would influence the analysis (Minami et al. 2019; Minami & Komatsu 2020a). We later show how the filamentary dust model of Clark et al. (2021) can create non-zero EB correlations of dust and how to implement an ansatz motivated by their model into our equations.

The last term in Eq. (8) can be assumed to be zero since ΛCDM does not predict any parity-violating correlations at the last scattering surface. However, chiral gravitational waves (Thorne et al. 2018) would predict intrinsic EB correlations of the CMB emissions. Thus, we could, in principle, use Eq. (8) to probe this term, but for the purpose of this work, we set $\langle C_{\mathcal{l}}^{EB,\text{CMB}}\rangle =0$.

Before applying this to real data, we need to generalize the method to incorporate multiple frequency bands (Minami & Komatsu 2020a). This work is based on the NPIPE pipeline which splits the detectors of each of the HFI frequency bands and the 70 GHz band into two groups, namely the A and B split. Therefore, for each of these frequency bands, we have two maps of the sky and two miscalibration angles per band, which act as the average over the detectors for a given split. We also analyze the frequency dependence of a possible cosmic birefringence angle β_ν. Thus, we denote i as a specific frequency band and data split for the miscalibration angle α_i, but for β_i, it only denotes the frequency of a given band. Following Minami & Komatsu (2020a) with the exception of having a frequency-dependent β_i, we now write observed spherical harmonics coefficients for polarization. For completeness, we derive the equations based on the assumption that we have access to the EB power spectrum of the polarized foreground emission. Assuming that the CMB has no intrinsic EB correlations, we get

$$\begin{array}{cc}\hfill \left[\begin{array}{c}⟨C_{\ell }^{E_i{E}_j,o}⟩\\ ⟨C_{\ell }^{B_i{B}_j,o}⟩\end{array}\right]=& \mathbf{R}({\alpha }_i,{\alpha }_j)\left[\begin{array}{c}⟨C_{\ell }^{E_i{E}_j,\mathrm{fg}}⟩\\ ⟨C_{\ell }^{B_i{B}_j,\mathrm{fg}}⟩\end{array}\right]+\mathbf{D}({\alpha }_i,{\alpha }_j)\left[\begin{array}{c}⟨C_{\ell }^{E_i{B}_j,\mathrm{fg}}⟩\\ ⟨C_{\ell }^{B_i{E}_j,\mathrm{fg}}⟩\end{array}\right]\hfill \\ \hfill & +\mathbf{R}({\alpha }_i+{\beta }_i,{\alpha }_j+{\beta }_j)\left[\begin{array}{c}⟨C_{\ell }^{E_i{E}_j,\mathrm{CMB}}⟩\\ ⟨C_{\ell }^{B_i{B}_j,\mathrm{CMB}}⟩\end{array}\right]+{\delta }_{i,j}\left[\begin{array}{c}⟨C_{\ell }^{E_i{E}_j,\mathrm{noise}}⟩\\ ⟨C_{\ell }^{B_i{B}_j,\mathrm{noise}}⟩\end{array}\right],\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill ⟨C_{\ell }^{E_i{B}_j,o}⟩=& R^T({\alpha }_i,{\alpha }_j)\left[\begin{array}{c}⟨C_{\ell }^{E_i{E}_j,\mathrm{fg}}⟩\\ ⟨C_{\ell }^{B_i{B}_j,\mathrm{fg}}⟩\end{array}\right]+D^T({\alpha }_i,{\alpha }_j)\left[\begin{array}{c}⟨C_{\ell }^{E_i{B}_j,\mathrm{fg}}⟩\\ ⟨C_{\ell }^{B_i{E}_j,\mathrm{fg}}⟩\end{array}\right]\hfill \\ & +R^T({\alpha }_i+{\beta }_i,{\alpha }_j+{\beta }_j)\left[\begin{array}{c}⟨C_{\ell }^{E_i{E}_j,\mathrm{CMB}}⟩\\ ⟨C_{\ell }^{B_i{B}_j,\mathrm{CMB}}⟩\end{array}\right].\hfill \end{array}$$(10)

Here,

$$\begin{array}{cc}\hfill \mathbf{R}({\theta }_i,{\theta }_j)& =\left[\begin{array}{cc}cos(2{\theta }_i)cos(2{\theta }_j)& sin(2{\theta }_i)sin(2{\theta }_j)\\ sin(2{\theta }_i)sin(2{\theta }_j)& cos(2{\theta }_i)cos(2{\theta }_j)\end{array}\right],\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill \mathbf{D}({\theta }_i,{\theta }_j)& =\left[\begin{array}{cc}-cos(2{\theta }_i)sin(2{\theta }_j)& -sin(2{\theta }_i)cos(2{\theta }_j)\\ sin(2{\theta }_i)cos(2{\theta }_j)& cos(2{\theta }_i)sin(2{\theta }_j)\end{array}\right],\hfill \end{array}$$(12)

$$\begin{array}{cc}\hfill R({\theta }_i,{\theta }_j)& =\left[\begin{array}{c}cos(2{\theta }_i)sin(2{\theta }_j)\\ -sin(2{\theta }_i)cos(2{\theta }_j)\end{array}\right],\hfill \end{array}$$(13)

$$\begin{array}{cc}\hfill D({\theta }_i,{\theta }_j)& =\left[\begin{array}{c}cos(2{\theta }_i)cos(2{\theta }_j)\\ -sin(2{\theta }_i)sin(2{\theta }_j)\end{array}\right].\hfill \end{array}$$(14)

We are not working with auto-spectra, so we can discard the noise term in Eq. (9). Combining Eqs. (9) and (10), we can eliminate $\left[\begin{array}{c}⟨C_{\ell }^{E_i{E}_j,\mathrm{fg}}⟩\\ ⟨C_{\ell }^{B_i{B}_j,\mathrm{fg}}⟩\end{array}\right]$ explicitly to end up with one equation,

$$\begin{array}{cc}\hfill ⟨C_{\ell }^{E_i{B}_j,o}⟩=& R^T({\alpha }_i,{\alpha }_j){\mathbf{R}}^{-1}({\alpha }_i,{\alpha }_j)\left[\begin{array}{c}⟨C_{\ell }^{E_i{E}_j,o}⟩\\ ⟨C_{\ell }^{B_i{B}_j,o}⟩\end{array}\right]+[R^T({\alpha }_i+{\beta }_i,{\alpha }_j+{\beta }_j)\hfill \\ & -R^T({\alpha }_i,{\alpha }_j){\mathbf{R}}^{-1}({\alpha }_i,{\alpha }_j)\mathbf{R}({\alpha }_i+{\beta }_i,{\alpha }_j+{\beta }_j)]\hfill \\ & \times \left[\begin{array}{c}⟨C_{\ell }^{E_i{E}_j,\mathrm{CMB}}⟩\\ ⟨C_{\ell }^{B_i{B}_j,\mathrm{CMB}}⟩\end{array}\right]+Z^T({\alpha }_i,{\alpha }_j)\left[\begin{array}{c}⟨C_{\ell }^{E_i{B}_j,\mathrm{fg}}⟩\\ ⟨C_{\ell }^{B_i{E}_j,\mathrm{fg}}⟩\end{array}\right],\hfill \end{array}$$(15)

where

$$\begin{array}{cc}\hfill Z^T({\alpha }_i,{\alpha }_j)& =D^T({\alpha }_i,{\alpha }_j)-R^T({\alpha }_i,{\alpha }_j){\mathbf{R}}^{-1}({\alpha }_i,{\alpha }_j)\mathbf{D}({\alpha }_i,{\alpha }_j)\hfill \\ \hfill & =\frac{2[cos(2{\alpha }_i)cos(2{\alpha }_j),sin(2{\alpha }_i)sin(2{\alpha }_j)]}{cos(4{\alpha }_i)+cos(4{\alpha }_j)}·\hfill \end{array}$$(16)

We now switch from the ensemble average spectra, ⟨C_ℓ⟩, to the estimated power spectra, C_ℓ. Excluding auto-spectra, we can group all combinations of (i, j) into a vector of observed spectra ${\mathit{C}}_{\mathcal{l}}^o={[C_{\mathcal{l}}^{E_i{E}_j,o},C_{\mathcal{l}}^{B_i{B}_j,o},C_{\mathcal{l}}^{E_i{B}_j,o}]}^T$ and vectors of beam-smoothed theoretical ΛCDM spectra ${\mathit{C}}_{\mathcal{l}}^{\mathrm{\Lambda }\text{CDM}}={[C_{\mathcal{l}}^{E_i{E}_j,\mathrm{\Lambda }\text{CDM}},C_{\mathcal{l}}^{B_i{B}_j,\mathrm{\Lambda }\text{CDM}}]}^T$. If we wish to include the intrinsic EB-modes of the foreground, we can make vectors of the polarized foreground EB-spectra ${\mathit{C}}_{\mathcal{l}}^{\text{fg}}={[C_{\mathcal{l}}^{E_i{B}_j,\text{fg}},C_{\mathcal{l}}^{B_i{E}_j,\text{fg}}]}^T$. As an example, we later use ten maps which gives us that the length of ${\mathit{C}}_{\mathcal{l}}^o$ is 3 ⋅ (10 ⋅ 10 − 10) = 270, where we excluded the auto spectra, while ${\mathit{C}}_{\mathcal{l}}^{\mathrm{\Lambda }\text{CDM}}$ and ${\mathit{C}}_{\mathcal{l}}^{\text{fg}}$ have 2 ⋅ (10 ⋅ 10 − 10) = 180 elements for each multipole ℓ. To solve this equation for all combinations of input maps, we group them into block diagonal matrices:

$$\begin{array}{cc}\hfill {\mathbf{A}}_{\mathit{ij}}=& [-R^T({\alpha }_i,{\alpha }_j){\mathbf{R}}^{-1}({\alpha }_i,{\alpha }_j),1],\hfill \end{array}$$(17)

$$\begin{array}{cc}\hfill {\mathbf{B}}_{\mathit{ij}}=& [R^T({\alpha }_i+{\beta }_i,{\alpha }_j+{\beta }_j)-R^T({\alpha }_i,{\alpha }_j){\mathbf{R}}^{-1}({\alpha }_i,{\alpha }_j)\hfill \\ \hfill & \times \mathbf{R}({\alpha }_i+{\beta }_i,{\alpha }_j+{\beta }_j)],\hfill \end{array}$$(18)

$$\begin{array}{cc}\hfill {\mathbf{Z}}_{\mathit{ij}}=& Z^T({\alpha }_i,{\alpha }_j),\hfill \end{array}$$(19)

so that we are solving the equation

$$\begin{array}{c}\hfill {\mathit{v}}_{\ell }^T\equiv \mathbf{A}{\mathit{C}}_{\ell }^o-\mathbf{B}{\mathit{C}}_{\ell }^{\mathrm{\Lambda }\mathrm{CDM}}-\mathbf{Z}{\mathit{C}}_{\ell }^{\mathrm{fg}}=0.\end{array}$$(20)

The covariance matrix for this equation is

$$\begin{array}{cc}\hfill {\mathbf{M}}_{\ell }=& \mathbf{A}\mathrm{Cov}({\mathit{C}}_{\ell }^o,{\mathit{C}}_{\ell }^o){\mathbf{A}}^T+\mathbf{B}\mathrm{Cov}({\mathit{C}}_{\ell }^{\mathrm{\Lambda }\mathrm{CDM}},{\mathit{C}}_{\ell }^{\mathrm{\Lambda }\mathrm{CDM}}){\mathbf{B}}^T\hfill \\ \hfill & +\mathbf{Z}\mathrm{Cov}({\mathit{C}}_{\ell }^{\mathrm{fg}},{\mathit{C}}_{\ell }^{\mathrm{fg}}){\mathbf{Z}}^T\hfill \\ \hfill & -\mathbf{A}\mathrm{Cov}({\mathit{C}}_{\ell }^o,{\mathit{C}}_{\ell }^{\mathrm{\Lambda }\mathrm{CDM}}){\mathbf{B}}^T-\mathbf{B}\mathrm{Cov}({\mathit{C}}_{\ell }^{\mathrm{\Lambda }\mathrm{CDM}},{\mathit{C}}_{\ell }^o){\mathbf{A}}^T\hfill \\ \hfill & -\mathbf{A}\mathrm{Cov}({\mathit{C}}_{\ell }^o,{\mathit{C}}_{\ell }^{\mathrm{fg}}){\mathbf{Z}}^T-\mathbf{Z}\mathrm{Cov}({\mathit{C}}_{\ell }^{\mathrm{fg}},{\mathit{C}}_{\ell }^o){\mathbf{A}}^T,\hfill \end{array}$$(21)

where we have assumed that there are no correlations between ΛCDM power spectra and the polarized foreground power spectra. Unfortunately, we do not have access to $C_{\mathcal{l}}^{EB,\text{fg}}$ explicitly. In the next section, however, we show how we can use a filamentary dust model to relate the EB spectra of the foreground to the E-modes and B-modes. In this work, we thus set ${\mathit{C}}_{\mathcal{l}}^{\text{fg}}=0$, and in the next section incorporate the effect of the foreground into the matrices A and B. We are also not including the ΛCDM power spectrum ${\mathit{C}}_{\mathcal{l}}^{\mathrm{\Lambda }\text{CDM}}$ in the covariance matrix in Eq. (21) since we find a negligible difference in the posterior distribution when we include it. Thus, we only use ${\mathbf{M}}_{\mathcal{l}}=\mathbf{A}\text{Cov}({\mathit{C}}_{\mathcal{l}}^o,{\mathit{C}}_{\mathcal{l}}^o){\mathbf{A}}^T$.

As in previous works (Diego-Palazuelos et al. 2022; Minami & Komatsu 2020b), we bin the spectra into groups of Δℓ = 20 multipoles to minimize noise and cosmic variance. We focus only on high-ℓ multipoles by using ℓ_min = 51 and ℓ_max = 1490, which gives the number of bins N_bins = 72. Explicitly, we get

$$\begin{array}{cc}& C_{\mathrm{b}}^{\mathit{XY}}=\frac{1}{\mathrm{\Delta }\ell }{\displaystyle \sum _{\ell \in b}{C}_{\ell }^{\mathit{XY}}}\hfill \end{array}$$(22)

$$\begin{array}{cc}& \mathrm{Cov}(C_{\mathrm{b}}^{\mathit{XY}},C_{\mathrm{b}}^{\mathit{ZW}})=\frac{1}{\mathrm{\Delta }{\ell }^2}{\displaystyle \sum _{\ell \in b}\mathrm{Cov}(C_{\ell }^{\mathit{XY}},C_{\ell }^{\mathit{ZW}}),}\hfill \end{array}$$(23)

where b is the bin number. To account for masks, we divide the covariance matrix, M_b, by the sky fraction f_sky which we define as

$$\begin{array}{c}\hfill f_{\mathrm{sky}}=\frac{1}{N_{\mathrm{pix}}}\frac{{\left({\sum }_{i=1}^{N_{\mathrm{pix}}}{w}_i^2\right)}^2}{{\sum }_{i=1}^{N_{\mathrm{pix}}}{w}_i^4}·\end{array}$$(24)

Here, w_i is the weight of the apodized mask at pixel i and N_pix is the number of pixels (Hivon et al. 2002; Challinor & Chon 2005).

We solve the equation by sampling over β_i and α_i using a Markov chain Monte Carlo to evaluate:

$$\begin{array}{c}\hfill lnL=-\frac{1}{2}{\displaystyle \sum _{b=1}^{N_{\mathrm{bin}}}({\mathit{v}}_{\mathrm{b}}^T{\mathbf{M}}_{\mathrm{b}}^{-1}{\mathit{v}}_{\mathrm{b}}+ln|{\mathbf{M}}_{\mathrm{b}}|),}\end{array}$$(25)

where M_b is the binned covariance matrix. The latter term ln|M_b| in the log-likelihood was used in Diego-Palazuelos et al. (2022) but not in Minami & Komatsu (2020a). In simulations, we found that by adding the term, we obtained posteriors for the angles that were better aligned with the input values, especially for smaller values of f_sky.

To get the theoretical ΛCDM power spectra of the CMB, we use CAMB¹ (Lewis et al. 2000) using the cosmological parameters from Planck 2018 (Planck Collaboration VI 2020). We then beam smooth the theoretical ΛCDM spectra so that

$$\begin{array}{c}\hfill {\mathit{C}}_{\ell }^{\mathrm{\Lambda }\mathrm{CDM}}=\left[\begin{array}{c}{C}_{\ell }^{EE,\mathrm{\Lambda }\mathrm{CDM}}{W}_{\ell }^{E_i{E}_j,E_i{E}_j}{w}_{\mathrm{pix},\ell }^2,C_{\ell }^{BB,\mathrm{\Lambda }\mathrm{CDM}}{W}_{\ell }^{B_i{B}_j,B_i{B}_j}{w}_{\mathrm{pix},\ell }^2\end{array}\right],\end{array}$$(26)

where $W_{\mathcal{l}}^{XY,XY}$ are the beam window matrices and w_pix, ℓ are the pixel window functions. Both i and j denote the frequency and data split.

We first use five polarization-sensitive frequency bands from the Planck satellite. From the HFI, we adopt all the polarization-sensitive bands, namely the 100, 143, 217, and 353 GHz channels, along with the 70 GHz band from the LFI. The NPIPE pipeline splits the data for each frequency band into two groups of maps, A and B. For the 70 GHz and higher frequency channels, these data splits are detector splits. This gives one average miscalibration angle for each set of detectors and, hence, we get ten miscalibration angles α_i where i = 70A,  …,353B. To get the observed power spectra ${\mathit{C}}_{\mathcal{l}}^o$ from these maps, we use PolSpice² (Chon et al. 2004).

We use the same masks as in Diego-Palazuelos et al. (2022). They all include a carbon monoxide (CO) mask that excludes pixels where there are bright CO emission lines. Unpolarized CO emission can cause intensity to polarization leakage, and so pixels with CO emission brighter than 45 K_RJ km s⁻¹ are excluded. These masks also exclude point sources that are the union of the polarized HFI maps. Four Galactic masks that cover 5%, 10%, 20%, and 30% are combined with the CO and point-source masks, leaving one mask with only a combined CO and point-source mask. This gives five masks with a sky coverage of f_sky = 0.93, 0.90, 0.85, 0.75, and 0.63.

We look at the frequency dependence of a possible cosmic birefringence angle β_ν where ν is the frequency of the photons. To test the frequency dependence of the birefringence angle, we will sample β_ν in two different ways. First, we sample β_ν for each frequency band individually for ν ∈ {70, 100, 143, 217, 353} GHz. To minimize the loss of the signal-to-noise ratio, we also chose to fit the cosmic birefringence angle to a power law β_ν = β₀(ν/ν₀)ⁿ. Here, ν₀ is a reference frequency that we set to 150 GHz, since we find it to give the smallest error bars on β₀.

In the last section, we analyze the inclusion of the 30 and 44 GHz frequency bands separately. This is due to their polarized foreground emission being dominated by synchrotron emission rather than dust emission. For the 30 and 44 GHz bands, the data splits A and B are, in fact, time splits, unlike the other detectors which are detector splits. We, therefore, sample one miscalibration angle α_i for each of the two frequency bands. The other bands, 70 GHz and higher frequency bands, still have two miscalibration angles each. Therefore, using all seven polarized Planck maps gives us a total of 12 miscalibration angles.

Even though the 70 GHz maps contain more polarized synchrotron emission than dust emission, we made the choice of including the 70 GHz channel in the analysis with the HFI bands. This is partly to see if the measured birefringence angle β_ν at 70 GHz is also heavily mask-dependent, which is hypothesized to be due to EB correlations caused by polarized dust emission which is sub-dominate at 70 GHz. The 70 GHz channel is also a band where there is a less total foreground contribution, and so, we might hope to get a stronger signal of β_ν at this frequency.

4. Filament model for EB correlations in dust

The largest uncertainty of Eq. (15) is the intrinsic EB correlations of the polarized foreground emission. Unlike the other power spectra, we do not have a reliable and robust theory of what this could be. However, as in Diego-Palazuelos et al. (2022) we make an ansatz based on the filamentary dust model of EB from Clark et al. (2021) as our best estimate of the intrinsic dust EB correlations.

First, we start by relating the intrinsic effect of the foreground’s EB correlations to the miscalibration angle α. By defining the correlation ratio as:

$$\begin{array}{c}\hfill r_{\ell }^{\mathit{XY}}\equiv \frac{C_{\ell }^{\mathit{XY}}}{\sqrt{C_{\ell }^{\mathit{XX}}{C}_{\ell }^{\mathit{YY}}}},\end{array}$$(27)

we can write $C_{\ell }^{EB,\mathrm{fg}}=r_{\ell }^{EB,\mathrm{fg}}\sqrt{C_{\ell }^{EE,\mathrm{fg}}{C}_{\ell }^{EB,\mathrm{fg}}}$. The intrinsic EE and BB-correlations of dust emission are approximately proportional to each other, $C_{\mathcal{l}}^{EB,\text{fg}}=\xi C_{\mathcal{l}}^{EE,\text{fg}}$ where ξ ≈ 0.5 (Planck Collaboration XI 2020). This gives (Minami et al. 2019)

$$\begin{array}{c}\hfill C_{\ell }^{EB,\mathrm{fg}}=\frac{r_{\ell }^{EB,\mathrm{fg}}\sqrt{\xi }}{1-\xi }(C_{\ell }^{EE,\mathrm{fg}}-C_{\ell }^{EB,\mathrm{fg}}).\end{array}$$(28)

For the sake of argument, let us assume that the correlation ratio is independent of multipoles, $r_{\mathcal{l}}^{EB,\text{fg}}=r^{EB,\text{fg}}$. Then we can define an effective angle γ, $sin(4\gamma )/2=r^{EB,\mathrm{fg}}\sqrt{\xi }/(1-\xi )$. This angle is completely degenerate with the miscalibration angle α. By ignoring the intrinsic EB correlations of the foreground emission, instead of measuring α and β, we are measuring α + γ and β − γ, respectively, which leaves the sum of α + β as invariant.

Interstellar dust dominates the Galactic polarized emission at frequencies ⪆100 GHz (Planck Collaboration XI 2020). The polarized dust emission was hypothesized in Huffenberger et al. (2020) to cause EB correlations when the filaments have their long axis partially misaligned with interstellar magnetic fields. The misalignment angle can be characterized by a multipole-dependent effective angle ${\psi }_{\mathcal{l}}^{\text{dust}}$ (Clark et al. 2021). This angle works as an averaged misalignment angle for filaments on the given angular scale. The estimate for this angle was given as:

$$\begin{array}{c}\hfill {\psi }_{\ell }^{\mathrm{dust}}=\frac{1}{2}arctan\left(\frac{C_{\ell }^{\mathit{TB}}}{C_{\ell }^{\mathit{TE}}}\right)·\end{array}$$(29)

With this angle, Clark et al. (2021) estimated the EB correlations to be

$$\begin{array}{c}\hfill C_{\ell }^{EB,\mathrm{dust}}=|r_{\ell }^{TB,\mathrm{dust}}|\sqrt{C_{\ell }^{EE,\mathrm{dust}}{C}_{\ell }^{BB,\mathrm{dust}}}sin(4{\psi }_{\ell }^{\mathrm{dust}}).\end{array}$$(30)

Due to the approximate proportionality between $C_{\mathcal{l}}^{EE,\text{dust}}$ and $C_{\mathcal{l}}^{BB,\text{dust}}$, and to avoid taking the square root of noisy data, we set explicitly $C_{\mathcal{l}}^{BB,\text{dust}}\to C_{\mathcal{l}}^{EE,\text{dust}}$. We therefore use the ansatz from Diego-Palazuelos et al. (2022):

$$\begin{array}{c}\hfill C_{\ell }^{EB,\mathrm{dust}}=A_{\ell }{C}_{\ell }^{EE,\mathrm{dust}}sin\left(4{\psi }_{\ell }^{\mathrm{dust}}\right),\end{array}$$(31)

where A_ℓ is equivalent to $|r_{\ell }^{\mathit{TB}}|$ which we then sample over different multipole ranges along with α_i and β_i. For small angles, we can relate this to the effective angle γ_ℓ:

$$\begin{array}{c}\hfill {\gamma }_{\ell }\approx A_{\ell }\frac{C_{\ell }^{EE,\mathrm{dust}}}{C_{\ell }^{EE,\mathrm{dust}}-C_{\ell }^{BB,\mathrm{dust}}}\frac{C_{\ell }^{TB,\mathrm{dust}}}{C_{\ell }^{TE,\mathrm{dust}}}·\end{array}$$(32)

To get the spectra needed to compute γ_ℓ, we take the cross-power spectrum of the A and B splits at 353 GHz with the respective mask. The polarization of the CMB emission is much less intense at this frequency, so it is expected that this frequency band describes polarized dust emission. Taking out A_ℓ, we calculate Eq. (32) with the 353 GHz spectra and apply it to all frequency bands. Since we are taking ratios of spectra, the spectral energy density of the power spectra cancels out for other frequency bands. Thus, we expect γ_ℓ to be independent of the frequency. This has been confirmed numerically by allowing A_ℓ to vary with frequency, which gives similar posterior distributions of the sampled parameters.

For the 70 GHz band, polarized thermal dust no longer dominates, and we expect synchrotron emission to be of similar or slightly higher amplitude than polarized dust emissions (Planck Collaboration Int. XXII 2015; Planck Collaboration IV 2020; BeyondPlanck Collaboration 2020). Nevertheless, we still apply the filamentary dust model to the 70 GHz band because we still expect a dust contribution. As we comment on more later, synchrotron emission has not been found to contain any intrinsic EB correlations (Martire et al. 2022), and so, we assume that polarized dust emission is the only foreground component that is expected to contain any intrinsic EB for the 70 GHz channel.

Similarly to Diego-Palazuelos et al. (2022), when we use the filamentary dust model, we sample four multipole ranges for A_ℓ in combination with sampling over β_ν and α_i. For A_ℓ, we choose the ℓ-ranges 51 ≤ ℓ ≤ 130, 131 ≤ ℓ ≤ 210, 211 ≤ ℓ ≤ 510, 511 ≤ ℓ ≤ 1490.

We can then relate the intrinsic EB of the foreground to the EE and BB: $C_{\mathcal{l}}^{E_i{B}_j,\text{fg}}=2{\gamma }_{\mathcal{l}}(C_{\mathcal{l}}^{E_i{E}_j,\text{fg}}-C_{\mathcal{l}}^{B_i{B}_j,\text{fg}})$. In matrix notation, this becomes:

$$\begin{array}{c}\hfill \left[\begin{array}{c}{C}_{\ell }^{E_i{B}_j,\mathrm{fg}}\\ C_{\ell }^{B_i{E}_j,\mathrm{fg}}\end{array}\right]={\mathbf{F}}_{\ell }\left[\begin{array}{c}{C}_{\ell }^{E_i{E}_j,\mathrm{fg}}\\ C_{\ell }^{B_i{B}_j,\mathrm{fg}}\end{array}\right],{\mathbf{F}}_{\ell }=2{\gamma }_{\ell }\left[\begin{array}{cc}1& -1\\ 1& -1\end{array}\right].\end{array}$$(33)

We can now put this expression into Eqs. (9) and (10), and we end up with modified A and B matrices (Diego-Palazuelos et al. 2022):

$$\begin{array}{cc}\hfill {\mathbf{A}}_{\ell ,ij}^{\prime }=& [-{\mathrm{\Lambda }}_{\ell }^T({\alpha }_i,{\alpha }_j){\mathrm{\Lambda }}_{\ell }^{-1}({\alpha }_i,{\alpha }_j),1],\hfill \end{array}$$(34)

$$\begin{array}{cc}\hfill {\mathbf{B}}_{\ell ,ij}^{\prime }=& [R^T({\alpha }_i+{\beta }_i,{\alpha }_j+{\beta }_j)-{\mathrm{\Lambda }}_{\ell }^T({\alpha }_i,{\alpha }_j){\mathrm{\Lambda }}_{\ell }^{-1}({\alpha }_i,{\alpha }_j)\hfill \\ \hfill & \times \mathbf{R}({\alpha }_i+{\beta }_i,{\alpha }_j+{\beta }_j)],\hfill \end{array}$$(35)

where we have defined

$$\begin{array}{cc}& {\mathrm{\Lambda }}_{\ell }({\alpha }_i,{\alpha }_j)=\mathbf{R}({\alpha }_i,{\alpha }_j)+\mathbf{D}({\alpha }_i,{\alpha }_j){\mathbf{F}}_{\ell },\hfill \end{array}$$(36)

$$\begin{array}{cc}& {\mathrm{\Lambda }}_{\ell }^T({\alpha }_i,{\alpha }_j)=R^T({\alpha }_i,{\alpha }_j)+D^T({\alpha }_i,{\alpha }_j){\mathbf{F}}_{\ell }.\hfill \end{array}$$(37)

The equation we are solving for now is

$$\begin{array}{c}\hfill {\mathbf{A}}_{\ell }^{\prime }{\mathit{C}}_{\ell }^o-{\mathbf{B}}_{\ell }^{\prime }{\mathit{C}}_{\ell }^{\mathrm{\Lambda }\mathrm{CDM}}=0,\end{array}$$(38)

with the corresponding covariance matrix ${\mathbf{M}}_{\mathcal{l}}^{\prime }$ using ${\mathbf{A}}_{\mathcal{l}}^{\prime }$ and ${\mathbf{B}}_{\mathcal{l}}^{\prime }$ in Eq. (21) with Z = 0. As mentioned, we find no difference in the posterior distribution when including the ΛCDM power spectra in the covariance, and therefore, we set ${\mathbf{M}}_{\mathcal{l}}^{\prime }={\mathbf{A}}^{\prime }\text{Cov}({\mathit{C}}_{\mathcal{l}}^o,{\mathit{C}}_{\mathcal{l}}^o){\mathbf{A}}^{\prime T}$.

5. Results

In this section, we report the results of the analysis on the HFI bands and the 70 GHz band from the LFI. We first report the nearly full-sky, f_sky = 0.93, measurement of a frequency-independent β. When neglecting the intrinsic EB of the polarized dust emission, we find β = 0.31° ±0.10°. When using the filamentary dust model, we find β = 0.37° ±0.11°. We note that these are slightly different from those reported in Diego-Palazuelos et al. (2022) since we are including the 70 GHz channel.

For a frequency-dependent value of β_ν, we show our nearly full-sky, f_sky = 0.93, measurement in Fig. 1. This shows both the individual measurements of β_ν for each frequency band and the power law model, β_ν = β₀(ν/ν₀)ⁿ, where we use ν₀ = 150 GHz, since we find that this reference frequency gives the smallest uncertainty in β₀. The power law models are shown as 1σ-bands, and we show the results both when using the filamentary model for the polarized dust emission and when we ignore the intrinsic EB of the foreground emission.

Fig. 1. Constraints on βν for each frequency band considered at nearly full-sky, fsky = 0.93. The constraints of the power law measurement, βν = β0(ν/ν0)n, are shown as colored 1σ-bands. The individual measurements of βν for each frequency channel are shown as error bars.

For the individual measurements, we quote the nearly full-sky values in Table 1, while for the power law model, we find β₀ = 0.26° ±0.11° and $n=-0.{45}_{-0.82}^{+0.61}$ when ignoring $C_{\mathcal{l}}^{EB,\text{fg}}$, and, otherwise, β₀ = 0.33° ±0.12° and $n=-0.{37}_{-0.64}^{+0.49}$ when we use the filamentary dust model for $C_{\mathcal{l}}^{EB,\text{fg}}$.

It was reported in Diego-Palazuelos et al. (2022) that removing the 353 GHz maps from the analysis gave a drop in the measurement of the frequency-independent birefringence angle, which suggests that the 353 GHz band increases the measured value of β. On the other hand, we note that the individual measurement of β_ν for the 353 GHz channel yields a negative value as seen in Fig. 1. This leaves the role of the 353 GHz channel unclear. Therefore, we investigate the effect of the 353 GHz frequency band on the β_ν measurements in Fig. 2 where we ignore the intrinsic EB power spectra of the foreground emission. The 1σ-bands are in this case the frequency-independent measurements, where the gray band is the constraint taken from Diego-Palazuelos et al. (2022) when excluding both the 70 and 353 GHz channels (β = 0.08° ±0.21°) and the light green 1σ-band excludes only the 353 GHz channel (β = 0.19° ±0.19°). We also show the measurement of β_ν when we sample it individually for each frequency band. The figure indicates that the 353 GHz band does not bias the measurement of the other bands significantly. Instead, it decreases the individual measurement uncertainties for all frequency bands. This can be explained by the domination of polarized dust emission at 353 GHz which helps break the degeneracy between α and β through the cross-power spectra with lower frequency channels.

Fig. 2. Constraints on βν for each frequency band considered at nearly full-sky, fsky = 0.93, when we ignore the EB correlations of the foreground. The colored 1σ-bands are frequency-independent measurements, unlike the power law model shown in Fig. 1.

The measured angles of a frequency-dependent β_ν for ν ∈ {70, 100, 143, 217, 353} GHz are shown in Fig. 3 as a function of the fractional sky coverage f_sky. The upper panel shows the probability distribution for the cosmic birefringence angle β_ν for each frequency band when we ignore the intrinsic EB correlations of the foreground. We observe a decreasing β_ν as we mask more of the Galactic plane for 100, 143, and 217 GHz. This is consistent with Diego-Palazuelos et al. (2022) where a similar drop was shown for a single frequency-independent β. We do not see this decline for the 70 GHz band. This band has a higher contribution from polarized synchrotron emission than the other bands, therefore, the mask-dependent contribution to $C_{\mathcal{l}}^{EB,\text{dust}}$ is less significant. In the lower panel of Fig. 3, we show the probability distribution where we instead use the filament model for $C_{\mathcal{l}}^{EB,\text{dust}}$. The drop in β_ν when we ignore the foreground EB is mitigated when using this model. This reinforces the hypothesis that the intrinsic EB-correlation of dust is mostly positive for large Galactic masks as reported by Clark et al. (2021).

Fig. 3. Constraints on a frequency-dependent β for various values of fsky. The upper panel assumes no intrinsic EB from the polarized foreground, while the lower panel shows the modeling of the foreground EB using the filament model.

For the power law measurement, β_ν = β₀(ν/ν₀)ⁿ, we show our result in Fig. 4 as a function of the sky fraction, f_sky. When β₀ = 0, n is fully degenerate. Since we find no theoretical models for large or small spectral indices, n, and no support from the individual measurements of β_ν that there should be a large spectral index n dependence, we set a flat prior |n|≤3. This prevents n from diverging around β₀ ≈ 0. When we remove the prior on n, we find no support for other solutions at higher |n|.

Fig. 4. Measurement of the power law, βν = β0(ν/ν0)n, for various values of fsky. The upper panel shows the measurement of β0, while the lower panel models the observed n. We set a flat prior, |n|≤3.

The decline in the cosmic birefringence angle parameter β₀ is also apparent at large Galactic masks when the intrinsic EB correlations of the polarized foreground emission are ignored. This drop is mostly mitigated when the filamentary dust model for EB is taken into account.

The spectral index n is shown in the lower panel of Fig. 4. When we ignore $C_{\mathcal{l}}^{EB,\text{fg}}$, we get results that are consistent with β₀ = 0° for larger sky cuts, and this generates broad uncertainties in n. The filamentary dust model yields a positive β₀, which allows us to probe the frequency dependence of the signal more accurately. For both nearly full-sky and large-sky cuts, we get a measurement that is consistent with a frequency-independent birefringence angle. The 68% C.L. lower limits for f_sky = 0.9 and f_sky = 0.85, touch n = −2, which would be expected from Faraday rotation caused by local magnetic fields.

The theoretical models considered in this work are based on integer values of n. Therefore, for nearly full-sky, f_sky = 0.93, we fix n to be the integer values, n ∈ { − 2, −1, 0, 1, 2}, and only measure β₀ in Table 2. We also give the corresponding Δχ² ≡ −2lnL with respect to n = 0, which shows that the frequency-invariant solution n = 0 is the favored integer solution.

In Fig. 5, we have included the posterior distribution of our power law assumption of β in the (n, β₀)-plane at nearly full-sky, f_sky = 0.93. The blue contour shows the probability distribution when we use the filamentary dust model for the EB foreground, while the green contour shows our measurement when we ignore the dust’s intrinsic EB. We see that when β₀ is close to zero, n is degenerate and seems to give n < 0. However, for larger values of β₀, we see support for a frequency-independent birefringence angle, n = 0.

Fig. 5. Constraints on the power law model βν = β0(ν/ν0)n for nearly full-sky fsky = 0.93. The green contour shows the probability distribution when we ignore the intrinsic EB of the foreground, while the blue contour shows that with the filament model for the foreground. The quoted values are for ignoring the EB correlations of the foreground. The equivalent values using the filament model for dust EB are β0 = 0.33° ±0.12° and $n=-0.{37}_{-0.64}^{+0.49}$.

This work is focused on probing the frequency dependence of β_ν assuming that the signal exists. Therefore, we check what constraints we get on n by having a flat prior of β₀ ≥ 0.1°. In that case, we find $n=-0.{43}_{-0.73}^{+0.58}$ when ignoring EB of the foreground emission and $n=-0.{34}_{-0.56}^{+0.48}$ when using the filament model.

6. Including all the Planck polarized frequency maps

In this section, we include the 30 and 44 GHz frequency channels from the LFI into our analysis. Thus, we are using all the polarized frequency maps from both the LFI and HFI of Planck. We treat the inclusion of these two frequency bands separately since they have a negligible contribution of polarized dust emission compared to the polarized synchrotron emission.

The EB power spectrum of synchrotron emission has not been detected yet. Martire et al. (2022) found that the EB of synchrotron emission is compatible with zero at 1σ using the 30 GHz band of Planck and the WMAP K-band for different sky fractions. The low signal-to-noise of these instruments puts weak constraints on it, but QUIJOTE and C-BASS aim to get a better understanding of the polarized synchrotron emission (Poidevin et al. 2018; Jones et al. 2018). Hence, we treat the full analysis of all Planck polarized maps with the inclusion of the 30 and 44 GHz channels separately, and we do not include any modeling of the intrinsic EB of either dust or synchrotron emission.

When including the two additional bands, we create a new point-source mask that is the union of all the LFI and HFI point-source masks. We combine it with the same CO and Galactic masks, and the new sky coverages become f_sky = 0.92, 0.88, 0.84, 0.74, and 0.62.

We present the nearly full-sky, f_sky = 0.92, measurements in Fig. 6. Here, we show the individual samples of β_ν for each frequency band, in addition to the power law measurement, β_ν = β₀(ν/ν₀)ⁿ, with ν₀ = 150 GHz shown as a 1σ-band. The values of the individual measurement are quoted in Table 3. For the power law model, we find ${\beta }_0={0.29}^{°}{}_{-{0.11}^{°}}^{+{0.10}^{°}}$ and $n=-0.{35}_{-0.47}^{+0.48}$. The posterior distribution for this measurement is shown in Fig. 7. We see that the individual measurements of β_ν at 30 and 44 GHz bands are positive, and they slightly increase the measured value of β₀ for the power law measurement as seen in Fig. 7 compared to Fig. 5. The inclusion of the two lowest frequency polarized LFI bands also help constrain better the spectral index n. The posterior distribution is slightly more centered around n = 0 with smaller uncertainties.

Fig. 6. Constraints on βν at nearly full-sky fsky = 0.92, where we have included the 30 and 44 GHz maps from the LFI. We ignore the intrinsic EB of the foreground emission, and the 1σ-band shows the constraints from the power law measurement, βν = β0(ν/ν0)n, where ν0 = 150 GHz. The dashed green line with the 1σ-band shows the frequency-independent β.

The correlation coefficient matrix for the sampled angles β_ν at nearly full-sky is shown in Fig. 8. We see little correlations between each LFI band and other frequencies, but the HFI channels give correlated measurements of β_ν. The strong correlations between β_ν for ν ∈ {100, 143, 217} GHz could be due to the high signal-to-noise of $C_{\mathcal{l}}^{\text{CMB},EE}$ for these bands. The CMB EE power spectrum is the same for all channels, and it strongly constrains the degeneracy α + β. Through cross-correlations with the 353 GHz band, we also get correlations with β₃₅₃. A combination of the low signal-to-noise ratio of $C_{\mathcal{l}}^{EE}$ and the foreground being synchrotron rather than dust at LFI could explain the lack of correlation between the LFI and HFI channels. To get a better understanding of the correlation coefficients, we could analyze simulations to see if the lower signal-to-noise level in the LFI bands are the cause of the smaller correlation coefficients as compared to the HFI bands. However, we leave such an analysis aside for future work.

We also try to fit the parameters of the power law model by using the covariance matrix of β_ν and the average measurements found in Table 3. We retrieve the values β₀ = 0.29° ±0.11° and $n=-0.{36}_{-0.47}^{+0.48}$, and we find that the shape of the posterior distribution is very consistent with the full (β₀, n) analysis seen in Fig. 7.

Fig. 7. Posterior distributions of n and β0 at nearly full-sky fsky = 0.92 for the power law measurement, βν = β0(ν/ν0)n, where ν0 = 150 GHz. We included the 30 and 44 GHz maps from the LFI, and we ignore the intrinsic EB of the foreground emission.

Fig. 8. Correlation coefficients between the individually sampled birefringence angles βν at nearly full-sky, fsky = 0.92.

We show β_ν as a function of the sky coverage, f_sky, in Fig. 9. This plot also includes the measurement when we assume no frequency dependence of β shown as the black line. At nearly full-sky, we measure this frequency-independent angle to be β = 0.33° ±0.10°. This figure can be compared to Fig. 1 in Diego-Palazuelos et al. (2022). The drop in the frequency-independent β as we mask more of the Galactic plane is partially mitigated by the inclusion of the LFI frequency bands. We also note that the individual measurements of β_ν for the HFI bands are generally lifted by the inclusion of 30 and 44 GHz bands as can be seen when comparing with Fig. 3. This is explained by the inclusion of a large number of new combinations of cross-power spectra with the 30 and 44 GHz channels. This supports the hypothesis that the EB power spectra of dust emission create a positive effective angle γ at large Galactic masks that gives a lower measured value of β. The 30 and 44 GHz should have a negligible contribution of dust compared to synchrotron emission.

Fig. 9. Constraints on βν when we sample them individually for each frequency band. We have included the 30 and 44 GHz maps from the LFI, and we ignore the intrinsic EB of the foreground emission. The black line shows β when we instead assume it is independent of frequency.

7. Conclusions

In this work, we employ one of the pipelines of Diego-Palazuelos et al. (2022) to explore the possibility of a frequency-dependent cosmic birefringence angle β_ν. First, we analyzed the HFI polarized frequency bands with the 70 GHz band from the LFI. At nearly full-sky we get a positive measurement of β_ν for all individual frequency bands except 353 GHz. As expected from earlier investigations (Diego-Palazuelos et al. 2022), the measured cosmic birefringence angle declines as a function of larger Galactic masks when we set the intrinsic EB correlations of the polarized foreground emission to zero. By using the EB dust ansatz of Diego-Palazuelos et al. (2022), motivated by the filamentary dust model of Clark et al. (2021), we can mitigate this effect.

We also sample a power law function of the birefringence angle, β_ν = β₀(ν/ν₀)ⁿ. Using both the filamentary dust model and neglecting the intrinsic dust EB, we get measurements that are inconsistent with a positive integer solution to the spectral index n, especially at large sky fractions. Some measurements at large sky coverages allow negative integer n solutions. However, we get the smallest uncertainty on n at f_sky = 0.93, which favors the frequency-independent solution, n = 0, over other integer solutions.

We then included the 30 and 44 GHz channels, thereby running the analysis on all the polarized Planck maps from both the LFI and HFI. These lowest frequency channels are dominated by polarized synchrotron emission. Measurements of the polarization of synchrotron emission indicate that the EB correlations must be small and is so far consistent with zero (Martire et al. 2022). So we did not try to model the intrinsic EB of the foreground emission.

The inclusion of all the LFI bands increased the mean value and tightened the constraints on the frequency-independent measurement of β at all sky fractions considered as compared to Diego-Palazuelos et al. (2022) which only used the HFI channels. Especially at large Galactic cuts, we find that including the 30 and 44 GHz bands raise the individual measurements of β_ν for the HFI bands. This supports the hypothesis of Diego-Palazuelos et al. (2022) and Clark et al. (2021) that large Galactic masks give a generally positive intrinsic EB of the polarized dust emission that biases the measurement of β to a lower value. The contribution of polarized dust emission at 30 and 44 GHz are negligible compared to synchrotron emission.

With all the polarized Planck maps, we measured β = 0.33° ±0.10° assuming β has no frequency dependence. We also sampled the parameters of the power law model, β_ν = β₀(ν/ν₀)ⁿ, at nearly full-sky for all the polarized LFI and HFI bands, and we found that the spectral index n is more constrained by adding the 30 and 44 GHz channels. The data favors n = 0 as the integer solution and does not favor Faraday rotation as the cause of the potential cosmic birefringence signal.

Diego-Palazuelos et al. (2022) did a thorough investigation of the instrumental systematics of the HFI frequency bands by looking at realistic simulations of the NPIPE processing. The impact of these systematics was found to be negligible for the measurement of β. However, such an analysis was not done for the LFI channels. This work has not analyzed the instrumental systematics of the LFI bands, and we leave this for future studies.

If the intrinsic EB of polarized dust emission was the cause of our positive measurement of β_ν, we would also have found a frequency-independent signal, as explained in Diego-Palazuelos et al. (2022). However, dust emission is not the dominant foreground contribution for the LFI channels; synchrotron is. Dust EB would, therefore, not directly explain a measurement of β_ν > 0 for the lowest frequency LFI bands. There have been measured correlations between polarized synchrotron and dust emission (Choi & Page 2015) which could create intrinsic foreground EB cross power spectra between synchrotron-dominated and dust-dominated channels. This could potentially bias our measurement of β_ν for frequencies where synchrotron emission dominates.

We used a filamentary dust model for the intrinsic EB of polarized dust emission, but to what extent can we trust this model to describe the EB of dust emission, especially at nearly full-sky? And does synchrotron emission have an intrinsic EB power spectrum? If the filamentary dust ansatz describes dust sufficiently and synchrotron has a negligible EB, the measurements presented in this work provide evidence for a mostly frequency-independent cosmic birefringence angle, but we cannot claim that we have a good enough understanding of the polarized foreground emission yet. Hence, we do not claim any statistical significance with regard to our cosmic birefringence measurements.

More work needs to be done to understand the physics of the polarized foreground emission, but strong priors on the miscalibration angles of future satellite missions can also help us ascertain whether we are measuring the foreground or a cosmic signal.

Acknowledgments

We are immensely grateful to Eiichiro Komatsu for guidance, suggestions, and discussions. We thank Hans Kristian Eriksen, Ingunn Wehus, and Duncan Watts for their comments on earlier drafts. We also thank Patricia Diego-Palazuelos and Yuto Minami for useful discussions and feedback on an earlier draft. We are also grateful to Matthieu Tristram and Raelyn M. Sullivan for useful discussions. Planck is a project of the European Space Agency (ESA) with instruments provided by two scientific consortia funded by ESA member states and led by Principal Investigators from France and Italy, telescope reflectors provided through a collaboration between ESA and a scientific consortium led and funded by Denmark, and additional contributions from NASA (USA). We acknowledge funding from the European Research Council (ERC) under the Horizon 2020 Research and Innovation Programme (Grant agreement No. 819478). Softwares: PolSpice (Chon et al. 2004), healpy (Gorski et al. 2005; Zonca et al. 2019), Matplotlib (Hunter 2007), NumPy (Harris et al. 2020), CAMB (Lewis et al. 2000), emcee (Foreman-Mackey et al. 2013).

References

All Tables

All Figures

	Fig. 1. Constraints on βν for each frequency band considered at nearly full-sky, fsky = 0.93. The constraints of the power law measurement, βν = β0(ν/ν0)n, are shown as colored 1σ-bands. The individual measurements of βν for each frequency channel are shown as error bars.
In the text

	Fig. 2. Constraints on βν for each frequency band considered at nearly full-sky, fsky = 0.93, when we ignore the EB correlations of the foreground. The colored 1σ-bands are frequency-independent measurements, unlike the power law model shown in Fig. 1.
In the text

	Fig. 3. Constraints on a frequency-dependent β for various values of fsky. The upper panel assumes no intrinsic EB from the polarized foreground, while the lower panel shows the modeling of the foreground EB using the filament model.
In the text

	Fig. 4. Measurement of the power law, βν = β0(ν/ν0)n, for various values of fsky. The upper panel shows the measurement of β0, while the lower panel models the observed n. We set a flat prior, |n|≤3.
In the text

Fig. 5. Constraints on the power law model βν = β0(ν/ν0)n for nearly full-sky fsky = 0.93. The green contour shows the probability distribution when we ignore the intrinsic EB of the foreground, while the blue contour shows that with the filament model for the foreground. The quoted values are for ignoring the EB correlations of the foreground. The equivalent values using the filament model for dust EB are β0 = 0.33° ±0.12° and $n=-0.{37}_{-0.64}^{+0.49}$.

In the text

	Fig. 6. Constraints on βν at nearly full-sky fsky = 0.92, where we have included the 30 and 44 GHz maps from the LFI. We ignore the intrinsic EB of the foreground emission, and the 1σ-band shows the constraints from the power law measurement, βν = β0(ν/ν0)n, where ν0 = 150 GHz. The dashed green line with the 1σ-band shows the frequency-independent β.
In the text

	Fig. 7. Posterior distributions of n and β0 at nearly full-sky fsky = 0.92 for the power law measurement, βν = β0(ν/ν0)n, where ν0 = 150 GHz. We included the 30 and 44 GHz maps from the LFI, and we ignore the intrinsic EB of the foreground emission.
In the text

	Fig. 8. Correlation coefficients between the individually sampled birefringence angles βν at nearly full-sky, fsky = 0.92.
In the text

	Fig. 9. Constraints on βν when we sample them individually for each frequency band. We have included the 30 and 44 GHz maps from the LFI, and we ignore the intrinsic EB of the foreground emission. The black line shows β when we instead assume it is independent of frequency.
In the text