© The Authors 2024

1. Introduction

Since the first results were released in 2013, the Planck satellite’s measurements of the cosmic microwave background (CMB) anisotropies have provided highly precise constraints on cosmological models. These measurements have tested the cosmological-constant-dominated cold dark matter (ΛCDM) model, given tight constraints on its parameters, and ruled out many plausible extensions. As a consequence, the best-fitting six-parameter ΛCDM model is now frequently used as the standard reference to be compared to new observational results as well as in combination with other data sets to provide further constraints.

Since the last Planck Collaboration cosmological analysis in 2018 (Planck Collaboration VI 2020), the very last version of the Planck data processing, called NPIPE, was released as the Planck Public Release 4 (PR4) and extensively detailed in Planck Collaboration Int. LVII (2020). In addition to drawing on previously neglected data from the repointing periods, NPIPE processed the entire set of Planck channels within the same framework, including the latest versions of corrections for systematics and data treatment.

In this paper, our objective is to enhance the precision on cosmological parameters through the utilization of PR4 data. Indeed, we expect to achieve a better sensitivity on almost all cosmological parameters owing to improved map sensitivity. Additionally, we look for better internal consistency for the lensing amplitude affecting the primordial CMB anisotropies. We thus derive constraints on cosmology using both low-ℓ and high-ℓ likelihoods based on Planck PR4. The only part that still relyies on PR3 (also known as Planck 2018) is the low-ℓ temperature likelihood, Commander, as we do not anticipate significant improvements at large scales in temperature between PR3 and PR4. On the other hand, our analysis includes the large scales in polarization from PR4 for which the NPIPE processing provides a significant improvement compared to PR3.

Since the foregrounds dominate polarization at large scales for the low-ℓ likelihood, namely, the LOw-ℓ LIkelihood on POlarized Power spectra (LoLLiPoP), we have made use of component-separated CMB maps processed by Commander using the whole range of Planck polarized frequencies from 30 to 353 GHz. This has been extensively discussed in Tristram et al. (2021, 2022), where it was combined with the BICEP2/Keck likelihood (Ade et al. 2021) to provide constraints on the tensor-to-scalar ratio, r.

For the high-ℓ power-spectrum analysis, HiLLiPoP, we use a multi-frequency Gaussian likelihood approximation using sky maps at three frequencies (100, 143, and 217 GHz), while the channel at 353 GHz is used to derive a template for the dust power spectrum contaminating the CMB signal at large scales. HiLLiPoP is one of the likelihoods developed within the Planck Collaboration and used to analyse previous Planck data sets (Planck Collaboration XV 2014; Planck Collaboration XI 2016). Here, we describe a new version adapted to PR4 and called “HiLLiPoP V4.2.” It differs from the previous one essentially by using a larger sky fraction (covering 75% of the sky) and a refined model for the foregrounds (particularly for point sources and dust emission). We specifically use the high-frequency instrument maps called “detsets,” which are made up of splits of the detectors at each frequency into specific subsets. We compute the cross-spectra for each of the CMB modes (TT, TE, EE) by cross-correlating the two detset maps at each of the three Planck channels dominated by the CMB (100, 143, and 217 GHz), together with their associated covariance. As illustrated in Fig. 1, the variance of the cross-spectra is close to the expected sample variance for 75% of the sky in temperature for TT, while the impact of the Planck noise in polarization is more visible in TE and EE. However, at those scales (ℓ < 2000), Planck PR4 is the most sensitive data set for CMB anisotropies as of today.

Fig. 1. Uncertainties on each angular cross-power spectrum (blue lines) and their combination (red line) for the PlanckTT (top), TE (middle), and EE (bottom) data, compared to sample variance for 75% of the sky (black dashed line).

The cross-spectra are then co-added into cross-frequency spectra and compared through a Gaussian likelihood to a model taking into account Galactic as well as extragalactic residual emission on the top of the CMB signal. As opposed to other Planck likelihoods, HiLLiPoP considers all cross-frequency power spectra. Even if the Planck PR4 data set is dominated by CMB anisotropies over the entire range of multipoles considered in the high-ℓ likelihood (30 < ℓ < 2500), the use of all cross-frequency spectra allows us to check the robustness of the results with respect to our knowledge of the astrophysical foregrounds. Indeed, even if the basic ΛCDM parameters are insignificantly affected by the details of the foreground modeling, the constraints on extensions to ΛCDM might depend more critically on the accuracy of the foreground description. Moreover, future ground-based experiments, measuring smaller scales than those accessible by Planck, will be even more sensitive to extragalactic foregrounds.

We begin this paper by summarizing the Planck PR4 pipeline (NPIPE), focusing on the improvements as compared to PR3 in Sect. 2. In Sect. 3, we explain how the angular power spectra were calculated and we describe the masks, multipole ranges, pseudo-C_ℓ algorithm, and covariance matrix we used. The LoLLiPoP likelihood is briefly described in Sect. 4, with reference to Tristram et al. (2021, 2022). The HiLLiPoP likelihood is described in Sect. 5, including details of foreground modeling and instrumental effects. Results for the parameters for the ΛCDM model are described and commented on in Sect. 6. Constraints on foreground parameters and instrumental parameters are discussed in Sects. 7 and 8, respectively. Section 9 is dedicated to consistency checks with respect to previous Planck results and Sect. 10 to the combination with other datasets. Finally we explore some extensions to ΛCDM in Sect. 11, specifically the lensing consistency parameter, A_L, the curvature, Ω_K, the effective number of neutrino species, N_eff, and the sum of neutrino masses, ∑m_ν.

2. The Planck PR4 data set

The Planck sky measurements used in this analysis are the PR4 maps available from the Planck Legacy Archive¹ (PLA) and from the National Energy Research Scientific Computing Center (NERSC)². They have been produced with the NPIPE processing pipeline, which creates calibrated frequency maps in temperature and polarization from both the Planck Low-Frequency Instrument (LFI) and the High-Frequency Instrument (HFI) data. As described in Planck Collaboration Int. LVII (2020), NPIPE processing includes data from the repointing periods that were neglected in previous data releases. There were additionally several improvements, resulting in lower levels of noise and systematics in both frequency and component-separated maps at essentially all angular scales, as well as notably improved internal consistency between the various frequencies. Moreover, PR4 also provides a set of “end-to-end” Monte Carlo simulations processed with NPIPE, which enables the characterization of potential biases and the uncertainties associated with the pipeline.

To compute unbiased estimates of the angular power spectra, we perform cross-correlations of two independent splits of the data. As shown in Planck Collaboration Int. LVII (2020), the most appropriate split for the Planck data is represented by the detset maps, comprising two subsets of maps with nearly independent noise characteristics, made by combining half of the detectors at each frequency. This was obtained by processing each split independently, in contrast to the split maps produced in the previous Planck releases. We note that time-split maps (made from, e.g., “odd-even rings” or “half-mission data”) share the same instrumental detectors and, therefore, exhibit noise correlations due to identical spectral bandpasses and optical responses. As a consequence, the use of time-split maps gives rise to systematic biases in the cross-power spectra (see Sect. 3.3.3 in Planck Collaboration V 2020), as well as underestimation of the noise levels in computing the half-differences (which needed to be compensated by a rescaling of the noise in PR3, as described in Appendix A.7 of Planck Collaboration III 2020). For this reason, we performed the cross-correlation using detset splits only.

Nevertheless, in order to verify the level of noise correlation between detsets, we computed the detset cross-power spectra from the half-ring difference maps, which we show in Fig. 2. The spectra are computed on 75% of the sky and are fully compatible with zero, ensuring that any correlated noise is much smaller than the uncorrelated noise over the range of multipoles from ℓ = 30 to 2500. As discussed above, this test is not sensitive to correlations at scales smaller than the half-ring period. Indeed, if both halves of a ring are affected by the same systematic effect, it will vanish in the half-ring difference map and, thus, it will not be tested in cross-correlation with another detset.

Fig. 2. Detset cross-spectra for half-ring differences computed on 75% of the sky, divided by their uncertainties. From top to bottom we show: TT, EE, TE, and ET. Spectra are binned with Δℓ = 40. The projections on the right show the distribution for each unbinned spectrum over the range ℓ = 30–2500.

3. Planck PR4 angular power spectra

3.1. Large-scale polarized power spectra

The foregrounds are stronger in polarization relative to the CMB than in temperature and cleaning the Planck frequencies using C_ℓ templates in the likelihood (as done at small scales) is not accurate enough, especially at large angular scales. In order to clean sky maps of polarized foregrounds, we used the Commander component-separation code (Eriksen et al. 2008), with a model that includes three polarized components, namely the CMB, synchrotron emission, and thermal dust emission. Commander was run on each detset map independently, as well as on each realization from the PR4 Monte Carlo simulations.

We then computed unbiased estimates of the angular power spectra by cross-correlating the two detset-cleaned maps. We computed the power spectra using an extension of the quadratic maximum-likelihood estimator (Tegmark & de Oliveira-Costa 2001) adapted for cross-spectra in Vanneste et al. (2018). At multipoles below 40, this has been shown to produce unbiased polarized power spectra with almost optimal errors. We used downgraded N_side = 16 maps (Górski et al. 2005) after convolution with a cosine apodizing kernel $b_{\ell }=\frac{1}{2}\{1+cos\pi (\ell -1)/(3N_{\mathrm{side}}-1)\}$. The signal is then corrected with the PR4 transfer function, to compensate for the filtering induced by the degeneracies between the signal and the templates for systematics used in the mapmaking procedure (see Planck Collaboration Int. LVII 2020).

The resulting power spectrum estimated on the cleanest 50% of the sky is plotted in Fig. 3 up to ℓ = 30 (for more details, see Tristram et al. 2021). We also performed the same estimation on each of the PR4 simulations and derived the ℓ-by-ℓ covariance matrix that was then used to propagate uncertainties in LoLLiPoP, the low-ℓ CMB likelihood described in Sect. 4.

Fig. 3. EE power spectrum of the CMB computed on 50% of the sky with the PR4 maps at low multipoles (Tristram et al. 2021). The Planck 2018 ΛCDM model is plotted in black. The grey band represents the associated sample variance. Error bars are deduced from the PR4 Monte Carlo simulations.

3.2. Small-scale power spectra

3.2.1. Sky fractions

For small scales (ℓ > 30), we used detset maps at frequencies of 100, 143, and 217 GHz, and we selected only a fraction of the sky to reduce the contamination from Galactic foregrounds. The main difference with respect to the masks used for the previous versions of HiLLiPoP (Couchot et al. 2017b) lies in two points: the new Galactic masks allow for a larger sky fraction and the point-source mask is common to all three frequencies. The resulting masks applied to each frequency are made of a combination of four main components, which we now describe.

Galactic mask. We applied a mask to remove the region of strongest Galactic emission, adapted to each frequency. We can keep a larger sky fraction at the lowest frequency (100 GHz) where the emission from the Galactic sources is low. Since Planck uncertainty is dominated by sample variance up to multipole ℓ ≃ 1800 in temperature (and ℓ ≃ 1100 in TE polarization), this allows us to reduce the sampling variance by ensuring a larger sky fraction. However, we removed a larger fraction of the sky for the highest frequency channel (217 GHz), since it is significantly more contaminated by Galactic dust emission.

We built Galactic masks using the Planck 353-GHz map as a tracer of the thermal dust emission in intensity. In practice, we smoothed the Planck 353-GHz map to increase the signal-to-noise ratio (S/N) before applying a threshold that depends on the frequency. Masks are then apodized using a 1° ​.0 Gaussian taper for power spectra estimation. For polarization, Planck dust maps show that the diffuse emission is strongly related to the Galactic magnetic field at large scales (Planck Collaboration Int. XIX 2015). However, at the smaller scales that matter here (ℓ > 30), the orientation of dust grains is driven by local turbulent magnetic fields that produce a polarization intensity approximately proportional to the total intensity dust map. We thus used the same Galactic mask for polarization as for temperature.

CO mask. We applied a mask for CO line emission. We considered the combination of maps of the two lines in the Planck frequency bands at 115 and 230 GHz. We smoothed the Planck reconstructed CO maps to 30 arcmin before applying a threshold at 2 K km s⁻¹. The resulting masks are then apodized at 15 arcmin. The CO masks remove 17% and 19% of the sky at 100 and 217 GHz, respectively, although the removed pixels largely fall within the Galactic masks.

Point-sources mask. We used a common mask for the three CMB frequencies to cover strong sources (both radio and infrared). In contrast to the masks used in Plik or CamSpec, the point-source mask used in our analysis relies on a more refined procedure that preserves Galactic compact structures and ensures the completeness level at each frequency, but with a higher flux cut on sources (approximately 340, 250, and 200 mJy at 100, 143, and 217 GHz, respectively). The consequence is that these masks leave a slightly greater number of unmasked extragalactic sources, but more accurately preserve the power spectra of dust emission (see Sect. 5.2). We apodized these masks with a Gaussian taper of 15 arcmin. We produce a single point-source mask as the combination of the three frequency masks; in total, this removes 8.3% of the sky.

Large objects. We masked a limited number of resolved objects in the sky, essentially nearby galaxies including the LMC, SMC, and M 31, as well as the Coma cluster. This removes less than 0.4% of the sky.

We used the same mask for temperature and polarization. Even though masking point sources in polarization is not mandatory (given the Planck noise in EE and TE); this makes the computation of the covariance matrix much simpler while not removing a significant part of the sky.

The Galactic masks ultimately used for HiLLiPoP V4.2 cover 20%, 30%, and 45% of the sky for the 100, 143, and 217 GHz channels, respectively. After combining with the other masks, the effective sky fraction used for computing cross-spectra are 75%, 66%, and 52%, respectively (see Fig. 4). The sky fractions retained for the likelihood analysis are about 5% larger than the ones used in the previous version of HiLLiPoP. Before extending the sky fraction used in the likelihood, we have checked the robustness of the results and the goodness-of-fit (through estimating χ²) using various combinations of Galactic masks (see Sect. 9).

Fig. 4. Sky masks used for HiLLiPoP V4.2 as a combination of a Galactic mask (blue, green, and red for the 100, 143, and 217 GHz channel, respectively), a CO mask, a point-source mask, and a mask removing nearby galaxies. The effective sky fractions remaining at 100, 143 and 217 GHz are 75%, 66%, and 52%, respectively.

3.2.2. PR4 small-scale spectra

We used Xpol (an extension to polarization of Xspect, described in Tristram et al. 2005) to compute the cross-power spectra in temperature and polarization (TT, EE, and TE). Overall, Xpol is a pseudo-C_ℓ method (see e.g., Hivon et al. 2002; Brown et al. 2005) that also computes an analytical approximation of the C_ℓ covariance matrix directly from data³. Using the six maps presented in Sect. 2, we derived the 15 cross-power spectra for each CMB mode, as outlined below: one each for 100 × 100, 143 × 143, and 217 × 217; and four each for 100 × 143, 100 × 217, and 143 × 217.

From the coefficients of the spherical harmonic decomposition of the (I,Q,U) masked maps ${\stackrel{\mathbf{\sim }}{\mathit{a}}}_{\ell m}^X=\{{\stackrel{\sim }{a}}_{\ell m}^T,{\stackrel{\sim }{a}}_{\ell m}^E,{\stackrel{\sim }{a}}_{\ell m}^B\}$, we form the pseudo cross-power spectra between map i and map j,

$$\begin{array}{c}\hfill {\stackrel{\sim }{\mathit{C}}}_{\ell }^{\mathit{ij}}=\frac{1}{2\ell +1}{\displaystyle \sum _m}{\stackrel{\mathbf{\sim }}{\mathit{a}}}_{\ell m}^{i\ast }{\stackrel{\mathbf{\sim }}{\mathit{a}}}_{\ell m}^j,\end{array}$$(1)

where the vector ${\stackrel{\mathbf{\sim }}{\mathit{C}}}_{\ell }$ includes the four modes $\{{\stackrel{\sim }{C}}_{\ell }^{TT},{\stackrel{\sim }{C}}_{\ell }^{EE},{\stackrel{\sim }{C}}_{\ell }^{\mathit{TE}},{\stackrel{\sim }{C}}_{\ell }^{ET}\}$. We note that the TE and ET cross-power spectra do not carry the same information, since computing T from map i and E from map j is different from computing E from map j and T from i. They were computed independently and averaged afterwards using their relative weights for each cross-frequency. The pseudo-spectra are then corrected for beam and sky fraction using a mode-mixing coupling matrix, M, which depends on the masks used for each set of maps (Peebles 1973; Hivon et al. 2002),

$$\begin{array}{c}\hfill {\stackrel{\mathbf{\sim }}{\mathit{C}}}_{\ell }^{\mathit{ij}}=(2{\ell }^{\prime }+1){\mathsf{M}}_{\ell {\ell }^{\prime }}^{\mathit{ij}}{\mathit{C}}_{{\ell }^{\prime }}^{\mathit{ij}}.\end{array}$$(2)

The Planck data set suffers from leakage of T to E and B, essentially due to beam mismatch between the detectors used to construct the (I, Q, U) maps. We debiased the beam leakage together with the beam transfer function using the beam window functions evaluated with QuickPol (Hivon et al. 2017). We used the QuickPol transfer functions specifically evaluated for PR4, since data cuts, glitch flagging, and detector noise weights all differ from earlier Planck releases. Once corrected, the cross-spectra are inverse-variance averaged for each frequency pair in order to form six unbiased (though correlated) estimates of the angular power spectrum.

The resulting cross-frequency spectra are plotted in Fig. 5 with respect to the C_ℓ average. For TT, the agreement between the different spectra is better than 20 μK², except (as expected) for the 100 × 100 and the 217 × 217 cases, which are affected by residuals from point sources and Galactic emission (for the latter). In EE, only the 217 × 217 case is affected by Galactic emission residuals at low multipoles, but the spectra are still consistent at the few μK² level. For TE and ET, we can see various features at the level of 10 μK² (especially for the 100T × 100E and 217E × 217T spectra). Even though the consistency between the cross-frequencies is very good, the likelihood presented in Sect. 5 will take into account those residuals from foreground emission.

Fig. 5. Frequency cross-power spectra with respect to the mean spectra for TT, EE, TE, and ET. Spectra are binned with Δℓ = 40 for this figure.

3.2.3. Multipole ranges

The HiLLiPoP likelihood covers the multipoles starting from ℓ_min = 30 up to ℓ_max = 2500 in temperature and ℓ_max = 2000 in polarization. The multipoles below ℓ < 30 are considered in the low-ℓ likelihoods (Commander and LoLLiPoP, see Sect. 4).

Table 1 gives the HiLLiPoP multipole ranges, [ℓ_min, ℓ_max], considered for each of the six cross-frequencies in TT, TE, and EE. The multipole ranges used in the likelihood analysis have been chosen to limit the contamination by Galactic dust emission at low ℓ and instrumental noise at high ℓ. In practice, we ignore the lowest multipoles for cross-spectra involving the 217 GHz map, where dust contamination is the highest, and cut out multipoles higher than ℓ = 1500 for cross-spectra involving the 100 GHz channel given its high noise level.

In total, the number of multipoles considered is now 29 758 for TT + TE + EE, to be compared to the number in the HiLLiPoP analysis of PR3, which was 25 597. The spectra are sample-variance limited up to ℓ ≃ 1800 in TT and ℓ ≃ 1100 in TE, while the EE mode is essentially limited by instrumental noise.

3.2.4. The covariance matrix

We use a semi-analytical estimate of the C_ℓ covariance matrix computed using Xpol. The matrix captures the ℓ-by-ℓ correlations between all the power spectra involved in the analysis. The computation relies directly on data for the estimates. It follows that contributions from noise (correlated and uncorrelated), sky emission (from astrophysical and cosmological origin), and the sample variance are implicitly taken into account in this computation without relying on any model or simulations.

The covariance matrix Σ of the cross-power spectra is directly related to the covariance $\stackrel{\sim }{\mathrm{\Sigma }}$ of the pseudo cross-power spectra through the coupling matrices:

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{{\ell }_1{\ell }_2}^{ab,cd}\equiv ⟨\mathrm{\Delta }{C}_{\ell }^{\mathit{ab}}\mathrm{\Delta }{C}_{{\ell }^{\prime }}^{cd\ast }⟩={\left(M_{\ell {\ell }_1}^{\mathit{ab}}\right)}^{-1}{\stackrel{\sim }{\mathrm{\Sigma }}}_{{\ell }_1{\ell }_2}^{ab,cd}{\left(M_{{\ell }^{\prime }{\ell }_2}^{cd\ast }\right)}^{-1},\end{array}$$(3)

with (a, b, c, d)∈{T, E} for each map.

The matrix $\stackrel{\sim }{\mathrm{\Sigma }}$, which gives the correlations between the pseudo cross-power spectra (ab) and (cd), is an N-by-N matrix (where N = 4ℓ_max) and reads:

$$\begin{array}{cc}& {\stackrel{\sim }{\mathrm{\Sigma }}}_{\ell {\ell }^{\prime }}^{ab,cd}\equiv ⟨\mathrm{\Delta }{\stackrel{\sim }{C}}_{\ell }^{\mathit{ab}}\mathrm{\Delta }{\stackrel{\sim }{C}}_{{\ell }^{\prime }}^{cd\ast }⟩=⟨{\stackrel{\sim }{C}}_{\ell }^{\mathit{ab}}{\stackrel{\sim }{C}}_{{\ell }^{\prime }}^{cd\ast }⟩-{\stackrel{\sim }{C}}_{\ell }^{\mathit{ab}}{\stackrel{\sim }{C}}_{{\ell }^{\prime }}^{cd\ast }\hfill \\ \hfill & ={\displaystyle \sum _{m{m}^{\prime }}}\frac{⟨{\stackrel{\sim }{a}}_{\ell m}^a{\stackrel{\sim }{a}}_{{\ell }^{\prime }{m}^{\prime }}^{c\ast }⟩⟨{\stackrel{\sim }{a}}_{\ell m}^{b\ast }{\stackrel{\sim }{a}}_{{\ell }^{\prime }{m}^{\prime }}^d⟩+⟨{\stackrel{\sim }{a}}_{\ell m}^a{\stackrel{\sim }{a}}_{{\ell }^{\prime }{m}^{\prime }}^{d\ast }⟩⟨{\stackrel{\sim }{a}}_{\ell m}^{b\ast }{\stackrel{\sim }{a}}_{{\ell }^{\prime }{m}^{\prime }}^c⟩}{(2\ell +1)(2{\ell }^{\prime }+1)},\hfill \end{array}$$

by expanding the four-point Gaussian correlation using Isserlis’ formula (or Wick’s theorem). We compute $\stackrel{\sim }{\mathrm{\Sigma }}$ for each pseudo cross-spectra block independently, which includes ℓ-by-ℓ correlation and four spectral mode correlations {TT, EE, TE, ET}.

Each two-point correlation of pseudo-a_ℓms can be expressed as the convolution of C_ℓ with a kernel that depends on the polarization mode considered:

$$\begin{array}{cc}& ⟨{\stackrel{\sim }{a}}_{\ell m}^{T_a\ast }{\stackrel{\sim }{a}}_{{\ell }^{\prime }{m}^{\prime }}^{T_b}⟩={\displaystyle \sum _{{\ell }_1{m}_1}}{C}_{{\ell }_1}^{T_a{T}_b}{W}_{\ell m{\ell }_1{m}_1}^{0,T_a}{W}_{{\ell }^{\prime }{m}^{\prime }{\ell }_1{m}_1}^{0,T_b\ast },\hfill \\ \hfill & ⟨{\stackrel{\sim }{a}}_{\ell m}^{E_a\ast }{\stackrel{\sim }{a}}_{{\ell }^{\prime }{m}^{\prime }}^{E_b}⟩=\frac{1}{4}{\displaystyle \sum _{{\ell }_1{m}_1}}\{C_{{\ell }_1}^{E_a{E}_b}{W}_{\ell m{\ell }_1{m}_1}^{+,E_a\ast }{W}_{{\ell }^{\prime }{m}^{\prime }{\ell }_1{m}_1}^{+,E_b}+C_{{\ell }_1}^{B_a{B}_b}{W}_{\ell m{\ell }_1{m}_1}^{-,E_a\ast }{W}_{{\ell }^{\prime }{m}^{\prime }{\ell }_1{m}_1}^{-,E_b}\},\hfill \\ \hfill & ⟨{\stackrel{\sim }{a}}_{\ell m}^{T_a\ast }{\stackrel{\sim }{a}}_{{\ell }^{\prime }{m}^{\prime }}^{E_b}⟩=\frac{1}{2}{\displaystyle \sum _{{\ell }_1{m}_1}}{C}_{{\ell }_1}^{T_a{E}_b}{W}_{\ell m{\ell }_1{m}_1}^{0,T_a\ast }{W}_{{\ell }^{\prime }{m}^{\prime }{\ell }_1{m}_1}^{+,E_b},\hfill \end{array}$$

where the kernels W⁰, W⁺, and W⁻ are defined as linear combinations of products of Y_ℓm of spin 0 and ±2, weighted by the spherical transform of the window function in the pixel domain (the apodized mask). As suggested in Efstathiou (2006), by neglecting the gradients of the window function and applying the completeness relation for spherical harmonics (Varshalovich et al. 1988), we can reduce the products of four Ws into kernels similar to the coupling matrix M defined in Eq. (2). In the end, the blocks of the Σ matrices are:

$$\begin{array}{cc}\hfill {\mathrm{\Sigma }}^{T_a{T}_b,T_c{T}_d}& \simeq C_{\ell {\ell }^{\prime }}^{T_a{T}_c}{C}_{\ell {\ell }^{\prime }}^{T_b{T}_d}{\mathsf{M}}_{TT,TT}+C_{\ell {\ell }^{\prime }}^{T_a{T}_d}{C}_{\ell {\ell }^{\prime }}^{T_b{T}_c}{\mathsf{M}}_{TT,TT},\hfill \\ \hfill {\mathrm{\Sigma }}^{E_a{E}_b,E_c{E}_d}& \simeq C_{\ell {\ell }^{\prime }}^{E_a{E}_c}{C}_{\ell {\ell }^{\prime }}^{E_b{E}_d}{\mathsf{M}}_{EE,EE}+C_{\ell {\ell }^{\prime }}^{E_a{E}_d}{C}_{\ell {\ell }^{\prime }}^{E_b{E}_c}{\mathsf{M}}_{EE,EE},\hfill \\ \hfill {\mathrm{\Sigma }}^{T_a{E}_b,T_c{E}_d}& \simeq C_{\ell {\ell }^{\prime }}^{T_a{T}_c}{C}_{\ell {\ell }^{\prime }}^{E_b{E}_d}{\mathsf{M}}_{TE,TE}+C_{\ell {\ell }^{\prime }}^{T_a{E}_d}{C}_{\ell {\ell }^{\prime }}^{E_b{T}_c}{\mathsf{M}}_{TT,TT},\hfill \\ \hfill {\mathrm{\Sigma }}^{T_a{T}_b,T_c{E}_d}& \simeq C_{\ell {\ell }^{\prime }}^{T_a{T}_c}{C}_{\ell {\ell }^{\prime }}^{T_b{E}_d}{\mathsf{M}}_{TT,TT}+C_{\ell {\ell }^{\prime }}^{T_a{E}_d}{C}_{\ell {\ell }^{\prime }}^{T_b{T}_c}{\mathsf{M}}_{TT,TT},\hfill \\ \hfill {\mathrm{\Sigma }}^{T_a{T}_b,E_c{E}_d}& \simeq C_{\ell {\ell }^{\prime }}^{T_a{E}_c}{C}_{\ell {\ell }^{\prime }}^{T_b{E}_d}{\mathsf{M}}_{TT,TT}+C_{\ell {\ell }^{\prime }}^{T_a{E}_d}{C}_{\ell {\ell }^{\prime }}^{T_b{E}_c}{\mathsf{M}}_{TT,TT},\hfill \\ \hfill {\mathrm{\Sigma }}^{E_a{E}_b,T_c{E}_d}& \simeq C_{\ell {\ell }^{\prime }}^{E_a{T}_c}{C}_{\ell {\ell }^{\prime }}^{E_b{E}_d}{\mathsf{M}}_{TE,TE}+C_{\ell {\ell }^{\prime }}^{E_a{E}_d}{C}_{\ell {\ell }^{\prime }}^{E_b{T}_c}{\mathsf{M}}_{TE,TE},\hfill \end{array}$$

which are thus directly related to the measured auto- and cross-power spectra (see the appendix in Couchot et al. 2017b). In practice, to avoid any correlation between C_ℓ estimates and their covariance, we used a smoothed version of each measured power spectrum (using a Gaussian filter with σ_ℓ = 5) to estimate the covariance matrix.

We finally average the cross-power spectra covariance matrix to form the full cross-frequency power-spectra matrices for the three modes {TT, TE, EE}. The resulting covariance matrix (Fig. 6) has 29 758 × 29 758 elements and is symmetric as well as positive definite.

Fig. 6. Full HiLLiPoP covariance matrix, including all correlations in multipoles between cross-frequencies and power spectra.

This semi-analytical estimation has been tested against the Monte Carlo simulations. In particular, we tested how accurate the approximations are in the case of a non-ideal Gaussian signal (due to the presence of small foregrounds residuals), Planck’s realistic (low) level of pixel–pixel correlated noise, and the apodization length used for the mask. We found no deviation to the sample covariance estimated from the 1000 realizations of the full focal plane Planck simulations that include anisotropic correlated noise and foreground residuals. To go further and check the detailed impact from the sky mask (including the choice of the apodization length), we simulated CMB maps from the Planck best-fit ΛCDM angular power spectrum, to which we added realistic anisotropic Gaussian noise (non-white, but without correlation) corresponding to each of the six data set maps. We then computed their cross-power spectra using the same foreground masks as for the data. A total of 15 000 sets of cross-power spectra were produced. When comparing the diagonal of the covariance matrix from the analytical estimation with the corresponding simulated variance, a precision better than a few percent is found (see Couchot et al. 2017b). Since we are using a Gaussian approximation of the likelihood, the uncertainty of the covariance matrix will not bias the estimation of the cosmological parameters. The percent-level precision obtained here will then only propagate into a sub-percent error on the variance of the recovered cosmological parameters.

4. Large-scale CMB likelihoods: LoLLiPoP and Commander

The Planck low-ℓ polarization likelihood, LoLLiPoP, is based on cross-spectra. It was first applied to Planck PR3 EE data for investigating the reionization history in Planck Collaboration Int. XLVII (2016). It was then upgraded to PR4 data and described in detail in Tristram et al. (2021, 2022), where it was used to derive constraints on the tensor-to-scalar ratio. LoLLiPoP can include EE, BB, and EB cross-power spectra calculated on component-separated CMB detset maps processed by Commander from the PR4 frequency maps. Here, we focus solely on the E-mode component.

Systematic effects are considerably reduced in cross-correlation compared to auto-correlation and LoLLiPoP is based on cross-power spectra for which the bias is zero when the noise is uncorrelated between maps. It uses the approximation presented in Hamimeche & Lewis (2008), modified as described in Mangilli et al. (2015) to apply to cross-power spectra. The idea is to apply a change of variable C_ℓ → X_ℓ, so that the new variable X_ℓ is nearly Gaussian-distributed. Similarly to Hamimeche & Lewis (2008), we define

$$\begin{array}{c}\hfill X_{\ell }=\sqrt{C_{\ell }^{\mathrm{f}}+O_{\ell }}g\left(\frac{{\stackrel{\sim }{C}}_{\ell }+O_{\ell }}{C_{\ell }+O_{\ell }}\right)\sqrt{C_{\ell }^{\mathrm{f}}+O_{\ell }},\end{array}$$(4)

where $g(x)=\sqrt{2(x-ln(x)-1)}$, ${\stackrel{\sim }{C}}_{\ell }$ are the measured cross-power spectra, C_ℓ are the power spectra of the model to be evaluated, $C_{\mathcal{l}}^{\text{f}}$ is a fiducial CMB model, and O_ℓ are the offsets needed in the case of cross-spectra. In the case of auto-power spectra, the offsets, O_ℓ, are given by the noise bias effectively present in the measured power spectra. For cross-power spectra, the noise bias is zero, and we use effective offsets defined from the C_ℓ noise variance:

$$\begin{array}{c}\hfill \mathrm{\Delta }{C}_{\ell }\equiv \sqrt{\frac{2}{2\ell +1}}{O}_{\ell }.\end{array}$$(5)

The distribution of the new variable X_ℓ can be approximated as Gaussian, with a covariance given by the covariance of the C_ℓs. The likelihood function of the C_ℓ given the data ${\stackrel{\sim }{C}}_{\ell }$ is then

$$\begin{array}{c}\hfill -2lnP(C_{\ell }|{\stackrel{\sim }{C}}_{\ell })={\displaystyle \sum _{\ell {\ell }^{\prime }}}{X}_{\ell }^{\mathsf{T}}{\mathsf{M}}_{\ell {\ell }^{\prime }}^{-1}{X}_{{\ell }^{\prime }}.\end{array}$$(6)

Uncertainties are incorporated into the C_ℓ covariance matrix M_ℓℓ′, which is evaluated after applying the same pipeline (including Commander component-separation and cross-spectrum estimation on each simulation) to the Monte Carlo simulations provided in PR4. While the foreground emission and the cleaning procedure are kept fixed in the simulations (so that we cannot include uncertainties arising from an imperfect foreground model), the resulting C_ℓ covariance consistently includes CMB sample variance, statistical noise, and systematic residuals, as well as uncertainties from the foreground-cleaning procedure, together with the correlations induced by masking. We further marginalized the likelihood over the unknown true covariance matrix (as proposed in Sellentin & Heavens 2016) in order to propagate the uncertainty in the estimation of the covariance matrix caused by a limited number of simulations. We note that LoLLiPoP is publicly available on GitHub⁴. In this work, we only considered the information from E modes and restricted the multipole range from ℓ = 2 to ℓ = 30.

To cover the low multipoles (ℓ < 30) in the temperature, we made use of the CommanderTT likelihood. It is based on a Bayesian posterior sampling that combines astrophysical component separation and likelihood estimation and employs Gibbs sampling to map out the full joint posterior (Eriksen et al. 2008). It was extensively used in previous Planck analyses (Planck Collaboration XV 2014; Planck Collaboration XI 2016). For the 2018 analysis, the version which is used in this work, Commander makes use of all Planck frequency channels, with a simplified foreground model including CMB, a unique low-frequency power-law component, thermal dust, and CO line emission (see Planck Collaboration V 2020).

5. Small-scale CMB likelihood: HiLLiPoP

This section describes the High-ℓ Likelihood on Polarized Power spectra (HiLLiPoP), including the models used for the foreground residuals and the instrumental systematic residuals. It was developed for the Planck 2013 results and then applied to PR3 and PR4 (e.g., Planck Collaboration XI 2016; Couchot et al. 2017c; Tristram et al. 2021). Here, we focus on the latest version of HiLLiPoP, released as V4.2⁵. We made use of the 15 cross-spectra computed from the six detset maps at 100, 143, and 217 GHz (see Sect. 3). From those 15 cross-spectra (one each for 100 × 100, 143 × 143, and 217 × 217; four each for 100 × 143, 100 × 217, and 143 × 217), we derived six cross-frequency spectra after recalibration and co-addition and compared them to the model. Using all cross-frequencies allows us to break some degeneracies in the foreground domain. However, because Planck spectra are dominated by sample variance, the six cross-frequency spectra are highly correlated. We used the full semi-analytic covariance matrix that includes the ℓ-by-ℓ correlation and {TT, TE, EE} mode correlation, as described in Sect. 3.2.4.

5.1. The likelihood approximation

On the full-sky, the distribution of auto-spectra is a scaled-χ² with 2ℓ+1 degrees of freedom. The distribution of the cross-spectra is slightly different (see Appendix A in Mangilli et al. 2015); however, above ℓ = 30, the number of modes is large enough that we can safely assume that the ${\stackrel{\sim }{C}}_{\ell }$ are Gaussian-distributed. Consequently, for high multipoles the resulting likelihood can be approximated by a multivariate Gaussian, including correlations between the values of C_ℓ arising from the cut-sky, and is expressed as:

$$\begin{array}{c}\hfill -2ln\mathcal{L}={\displaystyle \sum _{\begin{array}{c}i⩽j\hfill \\ i\prime ⩽j\prime \hfill \end{array}}\sum _{\ell {\ell }^{\prime }}}{\mathit{R}}_{\ell }^{\mathit{ij}}{\left[{\mathrm{\Sigma }}^{-1}\right]}_{\ell {\ell }^{\prime }}^{ij,i\prime j\prime }{\mathit{R}}_{{\ell }^{\prime }}^{i\prime j\prime }+ln|\mathrm{\Sigma }|,\end{array}$$(7)

where ${\mathit{R}}_{\ell }^{\mathit{ij}}={\stackrel{\mathbf{\sim }}{\mathit{C}}}_{\ell }^{\mathit{ij}}-{\mathit{C}}_{\ell }^{\mathit{ij}}$ denotes the residual of the estimated cross-power spectrum ${\stackrel{\mathbf{\sim }}{\mathit{C}}}_{\ell }$ with respect to the model C_ℓ, which depends on the frequencies {i, j} and is described in the next section. The matrix Σ = ⟨RR^T⟩ is the full covariance matrix that includes the instrumental variance from the data as well as the cosmic variance from the model. The latter is directly proportional to the model so that the matrix Σ should, in principle, depend on the model. In practice, given our current knowledge of the cosmological parameters, the theoretical power spectra typically differ from each other at each ℓ by less than they differ from the observed ${\stackrel{\sim }{C}}_{\ell }$, so that we can expand Σ around a reasonable fiducial model. As described in Planck Collaboration XV (2014), the additional terms in the expansion are small if the fiducial model is accurate and leaving it out entirely does not bias the likelihood. Using a fixed covariance matrix Σ, we can drop the constant term ln|Σ| and recover nearly optimal variance (see Carron 2013). Within the approximations discussed above, we expect the likelihood to be χ²-distributed with a mean equal to the number of degrees of freedom n_d.o.f. = n_ℓ − n_p (n_ℓ being the number of band powers in the power spectra and n_p the number of fitted parameters) and a variance equal to 2n_d.o.f..

5.2. The model

We now present the model (${\widehat{\mathit{C}}}_{\ell }$) used in the likelihood of Eq. (7). The foreground emission is mitigated by masking the part of the sky with high foreground signal (Sect. 3.2.1) and using an appropriate choice for the multipole range (Sect. 3.2.3). However, our likelihood function explicitly takes into account residuals of foreground emission in the power spectra, together with the CMB model and instrumental systematic effects. In practice, we consider the model and the data in the form D_ℓ = ℓ(ℓ + 1)C_ℓ/2π. In the foregrounds, for the temperature likelihood, we include the contributions from the: (1) Galactic dust; (2) cosmic infrared background (CIB); (3) thermal (tSZ) and kinetic (kSZ) Sunyaev–Zeldovich components; (4) Poisson-distributed point sources from radio and infrared star-forming galaxies; (5) the correlation between CIB and the tSZ effect (tSZ × CIB).

We highlight that this new version of HiLLiPoP, labelled V4.2, now includes a model for two point-source components, namely, dusty star-forming galaxies and radio sources. Consequently, the term “CIB” hereafter refers to the clustered part only. For all the components, we take into account the bandpass response using effective frequencies as listed in Table 4 of Planck Collaboration IX (2014). Galactic emission from free-free or synchrotron radiation is supposed to be weak at the frequencies considered here (above 100 GHz). Nevertheless, we implemented a model for such emission and we were not able to detect any residuals from Galactic synchrotron or free-free emission. Therefore, in the following, we neglect these contributions.

Galactic dust emission. At frequencies above 100 GHz, Galactic emission is dominated by dust. The dust template is fitted on the Planck 353-GHz data using a power-law model. In practice, we compute the 353-GHz cross-spectra, ${\widehat{\mathit{C}}}_{\ell }^{353A\times 353B}$, for each pair of masks (M_i, M_j) associated with the cross-spectra ν_i × ν_j (Fig. 7). We then subtract the Planck best-fit CMB power spectrum and fit a power-law model with a free constant Aℓ^α_d + B, in the range of ℓ = [30, 1500] for TT, to account for the unresolved point sources at 353 GHz. A simple power law is used to fit the EE and TE power spectra in the range ℓ = [30, 1000]. Thanks to the use of the point-source mask (described in Sect. 3.2.1), our Galactic dust residual power spectrum is much simpler than in the case of other Planck likelihoods. Indeed, the point-source masks used in the Planck PR3 analysis removes some Galactic structures and bright cirrus, which induces an artificial knee in the residual dust power spectra around ℓ = 200 (see Sect. 3.3.1 in Planck Collaboration XI 2016). In contrast with our point-source mask, the Galactic dust power spectra are fully compatible with power laws (Fig. 7). While the EE and TE power spectra are directly comparable to those derived in Planck Collaboration Int. XXX (2016), with indices of α_d = −2.3 and −2.4 for EE and TE, respectively, the indices for TT vary with the sky fraction considered, ranging from α_d = −2.2 down to −2.6 for the largest sky fraction.

Fig. 7. Dust power spectra, Dℓ = ℓ(ℓ + 1)Cℓ/2π, at 353 GHz for TT (top), EE (middle), and TE (bottom). The power spectra are computed from a cross-correlation between the detset maps at 353 GHz for different sets of masks, as defined in Sect. 3.2.1, and further corrected for the CMB power spectrum (solid black line) and CIB power spectrum (dashed black line). The coloured dashed lines are simple fits, as described in the text.

For each polarization mode (TT, EE, TE), we then extrapolated the dust templates at 353 GHz for each cross-mask to the cross-frequency considered:

$$\begin{array}{c}\hfill D_{\ell }^{\mathrm{dust}}(\nu \times \nu \prime )=c_{\mathrm{dust}}\frac{a_{\nu }^{\mathrm{dust}}}{a_{353}^{\mathrm{dust}}}\frac{a_{\nu \prime }^{\mathrm{dust}}}{a_{353}^{\mathrm{dust}}}{\mathcal{D}}_{\ell }^{\mathrm{dust}}(M_{\nu },M_{\nu \prime }),\end{array}$$(8)

where $a_{\nu }^{\text{dust}}={\nu }^{{\beta }_{\text{d}}}{B}_{\nu }(T_{\text{d}})$ is a modified blackbody with T_d fixed to 19.6 K, while c_dust and β_d are sampled independently for temperature and polarization. We use Gaussian priors for the spectral indices β_d from Planck Collaboration Int. XXII (2015), which gives ${\beta }_{\text{d}}^T=\mathcal{N}(1.51,0.01)$ and ${\beta }_{\text{d}}^T=\mathcal{N}(1.59,0.02)$ for the temperature and polarization, respectively. The coefficient c_dust allows us to propagate the uncertainty from fitting the 353-GHz dust spectrum with a power law. We sample c_dust with a Gaussian prior, c_dust = 𝒩(1.0, 0.1).

Cosmic infrared background (CIB). We use a template based on the halo model fitted on Planck and Herschel data (Planck Collaboration XXX 2014), extrapolated with a power-law at high multipoles. The template is rescaled by A^CIB, the amplitude of the contamination at our reference frequency (ν₀ = 143 GHz) and ℓ = 3000. The emission law is modelled by a modified blackbody $a_{\nu }^{\text{CIB}}={\nu }^{{\beta }_{\text{CIB}}}{B}_{\nu }(T)$ with a fixed temperature (T = 25 K) and a variable index β_CIB. We use a strong prior β_CIB = 𝒩(1.75, 0.06) (Planck Collaboration XXX 2014) and assume perfect correlation between the emission in the frequency range considered (from 100 to 217 GHz),

$$\begin{array}{c}\hfill D_{\ell }^{\mathrm{CIB}}(\nu \times \nu \prime )=A^{\mathrm{CIB}}\frac{a_{\nu }^{\mathrm{CIB}}}{a_{{\nu }_0}^{\mathrm{CIB}}}\frac{a_{\nu \prime }^{\mathrm{CIB}}}{a_{{\nu }_0}^{\mathrm{CIB}}}{\mathcal{D}}_{\ell }^{\mathrm{CIB}}.\end{array}$$(9)

Thermal Sunyaev–Zeldovich (tSZ) effect. The template for the tSZ emission comes from the halo model fitted on Planck measurements in Planck Collaboration XXII (2016) and used more recently with PR4 data in Tanimura et al. (2022). The tSZ signal is parameterized by a single amplitude A^tSZ, corresponding to the amplitude of the tSZ signal at our reference frequency (ν₀ = 143 GHz) at ℓ = 3000,

$$\begin{array}{c}\hfill D_{\ell }^{\mathrm{tSZ}}(\nu \times \nu \prime )=A^{\mathrm{tSZ}}\frac{a_{\nu }^{\mathrm{tSZ}}}{a_{{\nu }_0}^{\mathrm{tSZ}}}\frac{a_{\nu \prime }^{\mathrm{tSZ}}}{a_{{\nu }_0}^{\mathrm{tSZ}}}{\mathcal{D}}_{\ell }^{\mathrm{tSZ}},\end{array}$$(10)

where $a_{\nu }^{\mathrm{tSZ}}=x[e^x+1]/[e^x-1]-4$ (with x = hν/k_BT_CMB).

Kinetic Sunyaev–Zeldovich (kSZ) effect. The kSZ emission is parameterized by A^kSZ, the amplitude at ℓ = 3000, scaling a fixed template that includes homogeneous and patchy reionization components from Shaw et al. (2012) and Battaglia et al. (2013),

$$\begin{array}{c}\hfill D_{\ell }^{\mathrm{kSZ}}(\nu \times \nu \prime )=A^{\mathrm{kSZ}}{\mathcal{D}}_{\ell }^{\mathrm{kSZ}}.\end{array}$$(11)

Thermal SZ × CIB correlation. The cross-correlation between the thermal SZ and the CIB is parameterized as:

$$\begin{array}{cc}\hfill D_{\ell }^{\mathrm{tSZ}\times \mathrm{CIB}}(\nu \times \nu \prime )& =-\xi \sqrt{A^{\mathrm{tSZ}}{A}^{\mathrm{CIB}}}\hfill \\ \hfill & \times (\frac{a_{\nu }^{\mathrm{tSZ}}{a}_{\nu \prime }^{\mathrm{CIB}}+a_{\nu }^{\mathrm{CIB}}{a}_{\nu \prime }^{\mathrm{tSZ}}}{a_{{\nu }_0}^{\mathrm{tSZ}}{a}_{{\nu }_0}^{\mathrm{CIB}}}){\mathcal{D}}_{\ell }^{\mathrm{tSZ}\times \mathrm{CIB}},\hfill \end{array}$$(12)

with ξ the correlation coefficient rescaling the template ${\mathcal{D}}_{\mathcal{l}}^{\text{tSZ}\times \text{CIB}}$ from Addison et al. (2012).

Point sources. Point-source residuals in CMB data sets consist of a combination of the emission coming from radio and infrared sources. For earlier Planck data releases, HiLLiPoP used different point-source masks adapted to each frequency. This would require the estimation of the flux cut for each mask in order to use a physical model for the two point-source components. Since the flux-cut estimates are subject to large uncertainties, we used to fit one amplitude for the Poisson term at each cross-frequency in previous HiLLiPoP versions. In this new version of HiLLiPoP, we adopt a common mask for point sources (see Sect. 3.2.1). We then consider a flat Poisson-like power spectrum for each component and use a power law to describe the spectral energy distribution (SED) for the radio sources as $a_{\nu }^{\mathrm{rad}}\propto {\nu }^{-{\beta }_{\mathrm{s}}}$ (Tucci et al. 2011), while we use $a_{\nu }^{\text{IR}}={\nu }^{{\beta }_{\text{IR}}}{B}_{\nu }(T)$ (Béthermin et al. 2012) for infrared dusty star-forming galaxies. The residual cross-power spectra for point sources are:

$$\begin{array}{c}\hfill C_{\ell }^{\mathrm{PS}}(\nu \times \nu \prime )=A^{\mathrm{rad}}\frac{a_{\nu }^{\mathrm{rad}}}{a_{{\nu }_0}^{\mathrm{rad}}}\frac{a_{\nu \prime }^{\mathrm{rad}}}{a_{{\nu }_0}^{\mathrm{rad}}}+A^{\mathrm{IR}}\frac{a_{\nu }^{\mathrm{IR}}}{a_{{\nu }_0}^{\mathrm{IR}}}\frac{a_{\nu \prime }^{\mathrm{IR}}}{a_{{\nu }_0}^{\mathrm{IR}}}.\end{array}$$(13)

Following Lagache et al. (2020), the radio source emission is dominated at frequencies above about 100 GHz by radio quasars whose spectral indices can vary from −1.0 to 0.0 (Planck Collaboration XIII 2011; Planck Collaboration Int. VII 2013). We constrain the SED by fixing β_s = −0.8, following results from Reichardt et al. (2021). For infrared dusty star-forming galaxies, we adopt β_IR identical to β_CIB and T = 25 K. The C_ℓs are then converted into D_ℓs such that the amplitudes A^rad and A^IR refer to the amplitude of D₃₀₀₀ at 143 GHz. In the polarization, we do not include any contribution from point sources, since it is negligible compared to Planck noise for both components (Tucci et al. 2004; Lagache et al. 2020).

With the frequencies and the range of multipoles used in the HiLLiPoP likelihood, the foreground residuals are small in amplitude and mostly degenerate in the SED domain. As a result, we chose to set priors on the SED parameters, so that the correlation between the amplitudes of residuals would be significantly reduced. The optimization of the foreground model and, in particular, the determination of the priors adopted for the baseline analysis have been driven by astrophysical knowledge and results from the literature. We have extensively tested the impact of the priors using the ΛCDM model as a baseline (without any of its extensions). The results of these tests are discussed in Sect. 8.

5.3. Instrumental effects

The main instrumental effects that we propagate to the likelihood are the calibration uncertainties of each of the frequency maps in temperature and polarization (through the polarization efficiency). As a consequence, we sampled five inter-calibration coefficients, while fixing the calibration of the most sensitive map (the first detset at 143 GHz, 143A) as the reference. In addition, we sampled a Planck calibration parameter A_Planck with a strong prior, A_Planck = 𝒩(1.0000, 0.0025) to propagate the uncertainty coming from the absolute calibration based on the Planck orbital dipole.

We also allow for a recalibration of the polarized maps using polar efficiencies for each of the six maps considered. Those coefficients have been re-estimated in the NPIPE processing and we expect them to now be closer to unity and consistent within a frequency channel (Planck Collaboration Int. LVII 2020). By default, we fixed the polarization efficiencies to their best-fit values (unity at 100 and 143 GHz and 0.975 at 217 GHz; see Sect. 8 for details).

The angular power spectra have been corrected for beam effects using the beam window functions, including the beam leakage, estimated with QuickPol (see Sect. 3.2.2). With the improvement of the beam-estimation pipeline in Planck Collaboration XI (2016), the associated uncertainties have been shown to be negligible in Planck data and are ignored in this analysis.

A discrete sampling of the sky can lead to a small additive (rather than multiplicative) noise contribution known as the “subpixel” effect. Its amplitude depends on the temperature gradient within each pixel. With a limited number of detectors per frequency (and even more so per detset), the Planck maps are affected by the subpixel effect. However, the estimation of the size of the effect using QuickPol (Hivon et al. 2017), assuming fiducial spectra including CMB and foreground contributions, has shown it to be small (Planck Collaboration V 2020) and it is therefore neglected in this work.

6. Results on the six-parameter ΛCDM model

In this section, we describe the constraints on cosmological parameters in the ΛCDM model using the Planck PR4 data. In addition to HiLLiPoP (hlp), we also make use of the Commander low-ℓ likelihood (lowT, see Planck Collaboration IV 2020) and the polarized low-ℓ EE likelihood LoLLiPoP (lolE, discussed in Sect. 4). We define the following combination of likelihoods for the rest of the paper:

• TT, lowT+hlpTT;

• TE, lowT+lolE+hlpTE;

• EE, lolE+hlpEE;

• TTTEEE, lowT+lolE+hlpTTTEEE.

We note that for “TT”, we only used the temperature data and combined lowT+hlpTT; this is in contrast to Planck Collaboration VI (2020) and Rosenberg et al. (2022), in which low-ℓ data from EE are systematically added in order to constrain the reionization optical depth.

The model for the CMB was computed by numerically solving the background and perturbation equations for a specific cosmological model using CAMB (Lewis et al. 2000; Howlett et al. 2012)⁶. In this paper, we consider a ΛCDM model with six free parameters describing: the current physical densities of baryons (Ω_bh²) and cold dark matter (Ω_ch²); the angular acoustic scale (θ_*); the reionization optical depth (τ); and the amplitude and spectral index of the primordial scalar spectrum (A_s and n_s). Here, h is the dimensionless Hubble constant, h = H₀/(100 km s⁻¹ Mpc⁻¹).

In addition, we fit six inter-calibration parameters, seven foreground residual amplitudes in the temperature ($c_{\text{dust}}^{\text{T}}$, A_radio, A_IR, A_CIB, A_tSZ, A_kSZ, and ξ_{SZ × CIB}), plus one in polarization ($c_{\text{dust}}^{\text{P}}$), as well as three foreground spectral indices (${\beta }_{\text{dust}}^{\text{T}}$, ${\beta }_{\text{dust}}^{\text{P}}$, and β_CIB). Foreground and instrumental parameters are listed in Table A.1, together with their respective priors.

To quantify the agreement between the data and the model, we computed the χ² values with respect to the best-fit model for each of the data sets using Cobaya (Torrado & Lewis 2021) with its adaptive, speed-hierarchy-aware MCMC sampler (Lewis & Bridle 2002; Lewis 2013). The χ² values and the number of standard deviation from unity are given in Table 2. The goodness-of-fit is better than for previous Planck releases, but we still found a relatively large χ² value for hlpTT (corresponding to about 2.7σ), while the hlpTE and hlpEE χ² values are compatible with unity, at 1.8σ and 0.1σ, respectively. For the full combination hlpTTTEEE, we obtained χ² = 30 495 for a data size of 29 768, corresponding to a 3.02σ deviation. As described in Rosenberg et al. (2022), where the goodness of fit is also somewhat poor (4.07σ for TT and 4.46σ for the TTTEEE), this could be explained by a slight misestimation of the instrumental noise, rather than a bias that could be fit by an improved foreground model or a different cosmology. However, we emphasize that the level of this divergence is small, since the recovered reduced-χ², χ²/n_d = 1.02, shows that the semi-analytical estimation of the covariance of the data is accurate at the percent level. The goodness-of-fit values for individual cross-spectra are given in Table B.1.

Co-added CMB power spectra are shown in Figs. 8 and 9, for TT, TE, and EE; they are compared to the best-fit obtained with the full TTTEEE combination. Planck spectra are binned with Δℓ = 30 for the plots, but considered ℓ-by-ℓ in the likelihood. The plots also show the residuals relative to the ΛCDM best-fit to TTTEEE, as well as the normalized residuals. We cannot identify any deviation from statistical noise or any bias from foreground residuals.

Fig. 8. Maximum-likelihood frequency-co-added temperature power spectrum for HiLLiPoP V4.2. For the purposes of this figure, the power spectrum is binned with Δℓ = 30. The middle panel shows the residuals with respect to the fiducial base-ΛCDM cosmology and the bottom panel shows the residuals normalized by the uncertainties.

Fig. 9. As in Fig. 8, but for TE (top) and EE (bottom) power spectra.

In Fig. 10, we compare the constraints on ΛCDM parameters obtained using TT, TE, and EE and their combination. We find very good consistency between TT and TE, while EE constraints are wider, with a deviation in the acoustic scale θ_* toward lower values. This feature of the Planck PR4 data was previously reported in Rosenberg et al. (2022), in which the authors studied the correlation with other parameters and concluded that this is likely due to parameter degeneracies coupling to residual systematics in EE. However, the deviation of θ_* between EE and TT is now reduced with the increase of the sky fraction enabled by HiLLiPoP V4.2, though still present at the 1.6σ level. In addition, we have checked that this shift in θ_* is not related to any super-sample lensing effect (as described in Manzotti et al. 2014), or to any aberration correction (see Jeong et al. 2014), both of which are negligible for the large sky fraction considered in the Planck data set. We note that, interestingly, θ_* is the only parameter that deviates in EE; the others, including H₀, are compatible with TT at much better than 1σ. Given the weak sensitivity of the PlanckEE spectra as compared to TT and TE, discrepancies in the EE parameter reconstruction will have little impact on the overall cosmological parameter results.

Fig. 10. Posterior distributions for the cosmological parameters using power spectra from Planck PR4 with TT (lowT+hlpTT), TE (lowT+lolE+hlpTE), EE (lolE+hlpEE), and TTTEEE (lowT+lolE+hlpTTTEEE).

The HiLLiPoP V4.2 constraints on ΛCDM cosmological parameters are summarized in Table 3. As compared to the last Planck cosmological results in Planck Collaboration VI (2020), the constraints are tighter, with no major shifts. The error bars are reduced by 10–20%, depending on the parameter. The reionization optical depth is now constrained at close to the 10% level:

$$\begin{array}{c}\hfill \tau =0.058\pm 0.006.\end{array}$$(14)

This is the result of the NPIPE treatment of the PR4 data associated with the low-ℓ likelihood LoLLiPoP (see Planck Collaboration Int. LVII 2020).

For the constraint on the Hubble constant, we obtain:

$$\begin{array}{c}\hfill H_0=(67.64\pm 0.52)\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1},\end{array}$$(15)

which is consistent with previous Planck results and still significantly lower than the local distance-ladder measurements, which typically range from H₀ = 70 to 76, depending on the data set and the calibration used for the first step of the distance ladder (see for instance Abdalla et al. 2022).

The amplitude of density fluctuations is:

$$\begin{array}{c}\hfill {\sigma }_8=0.8070\pm 0.0065,\end{array}$$(16)

compatible with PR3 results (σ₈ = 0.8120 ± 0.0073) but lower by 0.5σ. The matter density, Ω_m, also shifts by roughly 1σ, such that:

$$\begin{array}{c}\hfill S_8\equiv {\sigma }_8{({\mathrm{\Omega }}_{\mathrm{m}}/0.3)}^{0.5}=0.819\pm 0.014.\end{array}$$(17)

Compared to PR3 (S₈ = 0.834 ± 0.016), this shift to a lower value of S₈ brings it closer to the measurements derived from galaxy clustering and weak lensing from the Dark Energy Survey Year 3 analysis (S₈ = 0.782 ± 0.019, for ΛCDM with fixed ∑m_ν, Abbott et al. 2022), decreasing the CMB versus large-scale structure tension on S₈ from 2.1σ to 1.5σ.

Before discussing results on the foreground parameters (Sect. 7) and instrumental parameters (Sect. 8), we show in Fig. 11 the correlation matrix for the fitted parameters. We can see that foreground parameters are only weakly correlated with the cosmological parameters and the inter-calibrations. This strengthens the robustness of the results with respect to the foreground model and ensures very low impact on cosmology.

Fig. 11. Correlation matrix for the fitted parameters of the combined HiLLiPoP likelihood TTTEEE. The first block corresponds to cosmological parameters from the ΛCDM model, the second block gathers the foreground parameters, and the last block shows the instrumental parameters.

7. Foreground parameters

All Planck cross-spectra are dominated by the CMB signal at all the scales we consider. This is illustrated for TT in Fig. B.1, where we show each component of the model fitted in the likelihood with the best-fit parameters for the six cross-frequencies. It is also true for TE and EE. Thanks to the multi-frequency analysis, we are able to break degeneracies related to the fact that some foreground-component power spectra are very similar. The resulting marginalized posteriors are plotted in Fig. 12. With the choice made for the multipole range and sky fraction, the Planck PR4 data set is sensitive to the CIB, the tSZ, and residual point sources (radio at 100 GHz and infrared at 217 GHz). Very low multipoles are sensitive to residuals from Galactic dust emission, especially at 217 GHz.

Fig. 12. Posteriors for foreground amplitudes. Units are μK2 normalized at ℓ = 3000 and ν = 143 GHz.

We detect the emission of radio point sources at better than 16σ. The preferred radio power in D_ℓ at ℓ = 3000 for 143 GHz is:

$$\begin{array}{c}\hfill A_{\mathrm{radio}}=(63.3\pm 4.7)\mathrm{\mu }{\mathrm{K}}^2,\end{array}$$(18)

with a population spectral index for the radio power fixed to β_s = −0.8, close to the value recovered by the SPT team (β_s = −0.76 ± 0.15, Reichardt et al. 2021). Allowing β_s to vary in Planck data, gives β_s = −0.54 ± 0.08, with a corresponding increase in the amplitude A_radio. This also impacts the SZ-CIB cross-correlation amplitude with a significant increase of ξ.

We obtain a high-significance detection of CIB anisotropies, with amplitudes at 143 GHz and ℓ = 3000, given by:

$$\begin{array}{cc}\hfill A_{\mathrm{CIB}}& =(1.03\pm 0.34)\mathrm{\mu }{\mathrm{K}}^2,\hfill \end{array}$$(19)

$$\begin{array}{cc}\hfill A_{\mathrm{IR}}& =(6.07\pm 0.63)\mathrm{\mu }{\mathrm{K}}^2,\hfill \end{array}$$(20)

for the clustered and Poisson parts, respectively. We note that these amplitudes cannot be directly compared to values in previous works because they strongly depend on the prior used for the β_CIB index for the former and on the flux cut applied by the point-source mask for the latter.

The thermal Sunyaev–Zeldovich effect is also significantly detected, with an amplitude at 143 GHz and ℓ = 3000 of:

$$\begin{array}{c}\hfill A_{\mathrm{tSZ}}=(5.9\pm 1.7)\mathrm{\mu }{\mathrm{K}}^2.\end{array}$$(21)

This is close to (but somewhat higher than) what is reported in Reichardt et al. (2021), with A_tSZ = (3.42 ± 0.54) μK², even though the uncertainties are larger. However, it is more closely comparable with ACTpol results, A_tSZ = (5.29 ± 0.66) μK² (Choi et al. 2020).

We find an upper-limit for the kSZ effect, while the correlation between tSZ and CIB is compatible with zero:

$$\begin{array}{cc}& A_{\mathrm{kSZ}}<7.6\mathrm{\mu }{\mathrm{K}}^2(\text{at}95\%\text{CL});\hfill \end{array}$$(22)

$$\begin{array}{cc}\hfill & {\xi }_{\mathrm{SZ}\times \mathrm{CIB}}=0.46\pm 0.30.\hfill \end{array}$$(23)

We note that those last results are about ten times less sensitive than the constraints from ground-based CMB measurements, such as those from SPT or ACTpol.

For the residuals of Galactic dust emission, with priors on the spectral indices driven by Planck Collaboration Int. XXII (2015), we found the rescaling coefficients, c_dust, to be 1.08  ±  0.03 and 1.20 ± 0.03 for the temperature and polarization, respectively. This indicates that we recover slightly more dust contamination than our expectations derived from the measurements at 353 GHz, especially in polarization. To estimate the impact on the reconstructed parameters (both cosmological and from foregrounds), we sampled the dust amplitudes at each frequency. The constraints are shown in Fig. 13 for temperature (top) and polarization (bottom). The figure illustrates that we have a good fit of the dust emission in temperature, while we are marginally sensitive to dust residuals in polarization. This explains why, given our prior on the SED for the polarized dust emission, ${\beta }_{\text{dust}}^P=\mathcal{N}(1.59,0.02)$, we ended up recovering an amplitude that was higher than expected.

Fig. 13. Amplitude of the dust emission relative to 353 GHz for a modified-blackbody dust model (blue line) as a function of the effective frequency (computed as the geometric mean of the two frequencies involved), compared to a fit using one amplitude per frequency (black dots). The top panel is for temperature and the bottom panel for polarization.

As discussed in Sect. 5.2, HiLLiPoP V4.2 also includes a two-component model for point sources. Figure 14 shows how the model, as the sum of the two point-source components, matches with the fit with one amplitude for each cross-frequency.

Fig. 14. Point-source model as a function of the effective frequency (computed as the geometric mean of the two frequencies involved), compared to the fit of one amplitude per cross-spectrum.

While adjusting the models as described above, the impact on ΛCDM parameters was shown to be very limited. We experienced variations of less than 0.11σ for all ΛCDM parameters, with the exception of n_s, can vary by 0.18σ when changing the model for point sources. Error bars on ΛCDM parameters are also stable with respect to foreground modeling, with variations limited to less than 2% (4% for n_s).

8. Instrumental parameters

Inter-calibration parameters are fitted in HiLLiPoP with respect to the first detset at 143 GHz (see Sect. 5.3). The inter-calibrations are recovered at better than the percent level and are compatible with unity. Using the full TTTEEE likelihood, we find:

$$\begin{array}{cc}& c_{100\mathrm{A}}=1.003\pm 0.007,\hfill \end{array}$$(24)

$$\begin{array}{cc}\hfill & c_{100\mathrm{B}}=1.004\pm 0.007,\hfill \end{array}$$(25)

$$\begin{array}{cc}\hfill & c_{143\mathrm{B}}=1.004\pm 0.006,\hfill \end{array}$$(26)

$$\begin{array}{cc}\hfill & c_{217\mathrm{A}}=1.001\pm 0.008,\hfill \end{array}$$(27)

$$\begin{array}{cc}\hfill & c_{217\mathrm{B}}=1.001\pm 0.008.\hfill \end{array}$$(28)

HiLLiPoP also allows us to fit for the polarization efficiency even though, by default, those are fixed. Using the full TTTEEE likelihood, we constrain the polarization efficiencies for each map at the percent level. The mean posteriors show polarization efficiencies compatible with unity at better than 1σ, except for the two maps at 217 GHz, which differ from unity by about 2σ:

$$\begin{array}{cc}& {\eta }_{100\mathrm{A}}=0.994\pm 0.013;\hfill \end{array}$$(29)

$$\begin{array}{cc}\hfill & {\eta }_{100\mathrm{B}}=0.987\pm 0.013;\hfill \end{array}$$(30)

$$\begin{array}{cc}\hfill & {\eta }_{143\mathrm{A}}=1.016\pm 0.013;\hfill \end{array}$$(31)

$$\begin{array}{cc}\hfill & {\eta }_{143\mathrm{B}}=1.001\pm 0.010;\hfill \end{array}$$(32)

$$\begin{array}{cc}\hfill & {\eta }_{217\mathrm{A}}=0.978\pm 0.013;\hfill \end{array}$$(33)

$$\begin{array}{cc}\hfill & {\eta }_{217\mathrm{B}}=0.972\pm 0.014.\hfill \end{array}$$(34)

Fixing polarization efficiencies to 1.00, 1.00, and 0.975 (at 100, 143, and 217 GHz, respectively) increases the χ² by Δχ² = 36 for 29 758 data points. However, this choice has no effect on either the ΛCDM parameters or the foreground parameters.

9. Consistency between Planck likelihoods

We go on to investigate the impact of the increased sky fraction used in this new version of HiLLiPoP. We repeat the analysis using more conservative Galactic masks reducing the sky fraction at each frequency by 5% (labelled “XL”) or 10% (labelled “L”) with respect to our baseline (“XXL”, which masks, 20%, 30%, and 45% at 100, 143, and 217 GHz, respectively; see Sect. 3.2.1 for more details). Within ΛCDM, we obtained similar χ² for the fits, demonstrating that the model used in HiLLiPoP V4.2 is valid for the considered sky fraction. For the TTTEEE likelihood, the Δχ² values are lower than 100 for 29 758 data points.

The other Planck likelihood using PR4 data is CamSpec and is described in detail in Rosenberg et al. (2022). Although CamSpec is focused on cleaning procedures to build co-added polarization spectra rather than modeling of foreground residuals in cross-frequency spectra, we find consistent constraints at better than the 1σ level. This gives confidence in the robustness of our cosmological constraints.

Figure 15 shows the 1D posterior distributions for the ΛCDM parameters using different sky fractions. We also make a comparison with the posteriors obtained from Planck PR3 and those of CamSpec PR4 (where we used LoLLiPoP instead of the polarized low-ℓ constraint from PR3 used in Rosenberg et al. 2022). We find good consistency between the different likelihoods and between the two data sets (PR3 and PR4).

Fig. 15. Posterior distributions for the cosmological parameters from PR4 for HiLLiPoP (using different sky fractions labelled L, XL, and XXL) and CamSpec, as compared to Planck 2018 (Plik PR3). Likelihoods are considered for the combination of TT+TE+EE, with lowT and lolE used at low ℓ.

Table 4 shows the relative difference in the cosmological parameters between Planck 2018 (Planck Collaboration VI 2020) and this work, together with the gain in accuracy. The largest difference with respect to Planck 2018 appears for Ω_ch², for which HiLLiPoP on PR4 finds a value 1.0σ lower. Associated with Commander and LoLLiPoP, CamSpec on PR4 also gives lower Ω_ch² by −0.45σ. The spectral index n_s is found to be a bit higher with HiLLiPoP by 0.7σ.

As discussed in Sect. 6, we obtain a slightly higher value for the Hubble constant (+0.6σ) with h = 0.6766 ± 0.0053, compared to h = 0.6727 ± 0.0060 for PR3. The amplitude of density fluctuations, σ₈, and the matter density, Ω_m, are lower by 0.7σ and 0.8σ, respectively, so that S₈ is also lower by about 0.9σ. The error bars shrink by more than 10%, with a noticeable gain of 20% for the acoustic scale (θ_*).

10. Combination with other data sets

We now present some results of our new likelihood in combination with CMB lensing measurements using the Planck PR4 data (Carron et al. 2022). We specifically use the conservative range recommended in Carron et al. (2022), consisting of nine power bins between multipoles of 8 and 400. The addition of the $C_{\mathcal{l}}^{\varphi \varphi }$ information means that we are using all the power spectra available from PR4; hence TTTEEE+lensing provides the best Planck-only cosmological constraints currently available.

We supplement this with measurements of the baryon acoustic oscillations (BAOs). This includes data from 6dF (Beutler et al. 2011), SDSS DR7 (specifically MGS, Ross et al. 2015), and SDSS DR16 (LRG, ELG, QSO, Ly-α auto, and Ly-αxQSO, Alam et al. 2021), which also incorporates some constraints on the growth of structures through redshift-space distortions.

Table 5 presents the constraints on the 6-parameter ΛCDM model when adding lensing and BAO data. Figure 16 shows the posterior distribution for the particular subset Ω_bh², Ω_m, σ₈, and H₀.

Fig. 16. Posterior distributions for some parameters using TTTEEE in combination with lensing and BAO.

11. Extensions

We now discuss constraints on some extensions to the base-ΛCDM model.

11.1. Gravitational lensing, A_L

We sample the phenomenological extension A_L in order to check the consistency of the Planck PR4 data set with the smoothing of the power spectra by weak gravitational lensing as predicted by the ΛCDM model. A mild preference for A_L > 1 was seen in the Planck PR1 data (Planck Collaboration XVI 2014) and since the analysis of Planck PR2 data (Planck Collaboration XI 2016; Planck Collaboration XIII 2016), HiLLiPoP has provided a significantly lower A_L value than the public Planck likelihood Plik, but still slightly higher than unity. The tension was at the 2.2σ level for PR3 (Couchot et al. 2017c).

With Planck PR4, we find the results to be even more compatible with unity compared to previous releases. Indeed for TTTEEE, we obtain:

$$\begin{array}{c}\hfill A_{\mathrm{L}}=1.039\pm 0.052,\end{array}$$(35)

which is compatible with the ΛCDM expectation (at the 0.7σ level). As shown in Table 6, while the results for EE and TE are compatible with unity, the A_L value for TT is still high by 0.8σ. Figure 17 shows posterior distributions of A_L for each of the mode-spectra and for the TTTEEE combination using Planck PR4.

Fig. 17. Posterior distributions for AL.

In Rosenberg et al. (2022), the CamSpec likelihood associated with low-ℓ likelihoods from Planck 2018 also showed a decrease in the A_L parameter in Planck PR4 data compared to PR3 data, reducing the difference from unity from 2.4σ to 1.7σ. When LoLLiPoP is adopted as the low-ℓ polarized likelihood, instead of the low-ℓ likelihoods from Planck 2018, the constraint on A_L from CamSpec changed from A_L = 1.095 ± 0.056 to A_L = 1.075 ± 0.058, still a 1.3σ difference from unity. We compare the posteriors for Plik (PR3), CamSpec (PR4), and HiLLiPoP (PR4) in Fig. 18.

Fig. 18. Posterior distributions for AL from HiLLiPoP PR4, compared to CamSpec (PR4) and Plik (PR3).

Previously, when there was a preference for A_L > 1, adding A_L as a seventh parameter could lead to shifts in other cosmological parameters (e.g., Planck Collaboration Int. LI 2017). However, we confirm that with HiLLiPoP on PR4, the ΛCDM parameters are only affected through a very slight increase of the error bars, without significantly affecting the mean posterior values.

With the PR4 lensing reconstruction described in Carron et al. (2022), the amplitude of the lensing power spectrum is 1.004 ± 0.024 relative to the Planck 2018 best-fit model. When combining CMB lensing with TTTEEE we could then recover a tighter constraint on A_L, with

$$\begin{array}{c}\hfill A_{\mathrm{L}}=1.037\pm 0.037\text{(TTTEEE+lensing)}.\end{array}$$(36)

11.2. Curvature, Ω_K

For the spatial curvature parameter, we report a significant difference with respect to Planck Collaboration VI (2020), which used PR3 and reported a mild preference for closed models (i.e., Ω_K < 0). Indeed, with HiLLiPoP V4.2, the measurements are consistent with a flat universe (Ω_K = 0) for all spectra.

As noticed in Rosenberg et al. (2022), with Planck PR4, the constraint on Ω_K is more precise and shifts toward zero, along the so-called geometrical degeneracy with H₀ (Fig. 19). Indeed, with HiLLiPoP V4.2 on PR4, the posterior is more symmetrical and the mean value of the posterior for TTTEEE is:

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_K=-0.012\pm 0.010,\end{array}$$(37)

Fig. 19. Posterior distributions in the ΩK–H0 plane using HiLLiPoP PR4, compared to CamSpec (PR4) and Plik (PR3).

which is only 1.2σ discrepant from zero.

This is to be compared to ${\mathrm{\Omega }}_K=-0.{044}_{-0.015}^{+0.018}$ obtained for Plik on PR3 (Planck Collaboration VI 2020) and ${\mathrm{\Omega }}_K=-0.{025}_{-0.010}^{+0.013}$ obtained with CamSpec on PR4 (Rosenberg et al. 2022).

As a consequence, the tail of the 2-d posterior in the H₀–Ω_K plane at low H₀ and negative Ω_K is no longer favoured. Indeed, when fitting for a non-flat Universe, the recovered value for the Hubble constant is H₀ = (63.03 ± 3.60) km s⁻¹ Mpc⁻¹, only 1.3σ away from the constraint with fixed Ω_K = 0.

The combination of TTTEEE with lensing yields the improved constraint

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_K=-0.0078\pm 0.0058\text{(TTTEEE+lensing)}.\end{array}$$(38)

This is now compatible with the baryon acoustic oscillation measurements from SDSS, which are consistent with a flat Universe and give Ω_K = −0.0022 ± 0.0022 (Alam et al. 2021). Finally, the mean posterior for the combination of Planck PR4 TTTEEE with lensing and BAO is:

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_K=0.0000\pm 0.0016\text{(TTTEEE+lensing+BAO)}.\end{array}$$(39)

This is consistent with our Universe being spatially flat to within a 1σ accuracy of 0.16% (see Fig. 20).

Fig. 20. Posterior distributions in the ΩK–H0 plane using Planck PR4 TTTEEE (i.e., lowT+lolE+hlpTTTEEE) in combination with lensing and BAO.

11.3. Effective number of relativistic species, N_eff

Figure 21 shows the posteriors for TT, TE, EE, and their combination when we consider the N_eff extension. Both TT and TE are compatible with similar uncertainties, while EE is not sensitive to N_eff. The mean posterior for TTTEEE is:

$$\begin{array}{c}\hfill N_{\mathrm{eff}}=3.08\pm 0.17.\end{array}$$(40)

Fig. 21. Posterior distributions for Neff. The vertical dashed line shows the theoretical expectation (Neff = 3.044).

The uncertainties are comparable to Planck 2018 results (N_eff = 2.92 ± 0.19, Planck Collaboration VI 2020) with a slight shift toward higher values, closer to the theoretical expectation N_eff = 3.044 (Akita & Yamaguchi 2020; Froustey et al. 2020; Bennett et al. 2021), which was also reported with the CamSpec analysis based on PR4 data (N_eff = 3.00 ± 0.21, Rosenberg et al. 2022).

11.4. Sum of the neutrino masses, ∑m_ν

Figure 22 shows the posterior distribution for the sum of the neutrino masses, ∑m_ν. There is no detection of the effects of neutrino mass and we report an upper limit of:

$$\begin{array}{c}\hfill {\displaystyle \sum }{m}_{\nu }<0.39\text{eV}(95\%\text{CL, TTTEEE}).\end{array}$$(41)

Fig. 22. Posterior distributions for ∑mν. Units are electronvolts.

Despite the increase in sensitivity associated with PR4, the constraint is slightly weaker (the upper limit is larger) than the one reported for Planck 2018: ∑m_ν < 0.26 eV at 95% CL. Our constraint is comparable to CamSpec, which gives ∑m_ν < 0.36 eV at 95% CL.

As explained in Couchot (2017a) and Planck Collaboration VI (2020), this is directly related to the value of A_L. Indeed, the correlation between A_L and ∑m_ν pushes the peak posterior of ∑m_ν toward negative values when A_L is fixed to unity; the data, however, prefer values of A_L larger than 1. With HiLLiPoP V4.2, the value of A_L reported in this work is more compatible with unity (A_L = 1.039 ± 0.052, see Sect. 11.1), thus, the posterior for ∑m_ν is shifted to higher values, with a peak closer to zero, increasing the upper limit accordingly.

Figure 23 shows constraints in the ∑m_ν–τ plane when combining our new likelihood with with CMB lensing and BAO data. This combination further strengthens the limits to:

$$\begin{array}{cc}& {\displaystyle \sum }{m}_{\nu }<0.26\text{eV}(95\%\text{CL, TTTEEE+lensing}),\hfill \end{array}$$(42)

$$\begin{array}{cc}\hfill & {\displaystyle \sum }{m}_{\nu }<0.11\text{eV}(95\%\text{CL, TTTEEE+lensing+BAO}).\hfill \end{array}$$(43)

Fig. 23. Posterior distributions in the ∑mν–τ plane using Planck PR4 TTTEEE (i.e., lowT+lolE+hlpTTTEEE) in combination with lensing and BAO.

This is slightly tighter than the upper limit from Planck 2018 (∑m_ν < 0.12 eV) and getting close to the lower limit for the inverted mass hierarchy (∑m_ν ≳ 0.1 eV, see e.g., Jimenez et al. 2022).

12. Conclusions

In this paper, we have derived cosmological constraints using CMB anisotropies from the final Planck data release (PR4). We detailed a new version of a CMB high-ℓ likelihood based on cross-power spectra computed from the PR4 maps. This version of HiLLiPoP, labelled V4.2, uses more sky (75%) and a wider range of multipoles. Our likelihood makes use of physically-motivated models for foreground-emission residuals. Using only priors on the foreground spectral energy distributions, we found amplitudes for residuals consistent with expectations. Moreover, we have shown that the impact of this modeling on cosmological ΛCDM parameters is negligible.

Combined with the low-ℓ EE likelihood LoLLiPoP, we derived constraints on ΛCDM and find good consistency with Planck 2018 results (based on PR3) with better goodness-of-fit and higher sensitivity (from 10% to 20%, depending on the parameters). In particular, we now constrain the reionization optical depth at the 10% level. We found a value for the Hubble constant consistent with previous CMB measurements and thus still in tension with distance-ladder results. We also obtained a lower value for S₈, alleviating the CMB versus large-scale structure tension to 1.5σ.

We found good consistency with the other published CMB likelihood analysis based on PR4, CamSpec (Rosenberg et al. 2022), which relies on a procedure to clean power spectra prior to constructing the likelihood. The consistency of the results using two different approaches reinforces the robustness of the results obtained with Planck data.

We also add constraints from PR4 lensing, making the combination the most constraining data set that is currently available from Planck. Additionally we explore adding baryon acoustic oscillation data, which tightens limits on some particular extensions to the standard cosmology.

We provide constraints on a number of extensions to ΛCDM, including the lensing amplitude, A_L, the curvature, Ω_K, the effective number of relativistic species, N_eff, and the sum of the neutrino masses, ∑m_ν. For both A_L and Ω_K, our results show a significant reduction of the so-called "tensions" within the standard ΛCDM, along with a reduction of the uncertainties. Indeed, the final constraints are fully compatible with ΛCDM predictions. In particular, with the new version of the likelihood presented in this work, we report A_L = 1.039 ± 0.052, which is entirely compatible with the ΛCDM prediction. The better agreement is explained both by the improvement of the Planck maps thanks to the NPIPE processing (with less noise and better systematic control) and the use of the LoLLiPoP and HiLLiPoP likelihoods.

Acknowledgments

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). Some of the results in this paper have been derived using the HEALPix package. We acknowledge use of the following packages: xQML, for the computation of large-scale power spectra (https://gitlab.in2p3.fr/xQML); Xpol, for the computation of large-scale power spectra (https://gitlab.in2p3.fr/tristram/Xpol); Cobaya, for the sampling of the likelihoods (https://github.com/CobayaSampler); and CLASS (https://github.com/lesgourg/class_public) and CAMB (https://github.com/cmbant/CAMB) for calculating power spectra. We gratefully acknowledge support from the CNRS/IN2P3 Computing Center for providing computing and data-processing resources needed for this work. This research was enabled in part by support provided by the Digital Research Alliance of Canada (https://alliancecan.ca/en. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 788212). The Planck PR4 data are publicly available on the Planck Legacy Archive (https://pla.esac.esa.int). Both likelihoods LoLLiPoP and HiLLiPoP based on PR4 are publicly available on GitHub (https://github.com/planck-npipe) as external likelihoods for Cobaya.

References

Appendix A: Foregrounds and instrumental parameters

Here, we describe the “nuisance” parameters relating to foreground emission components and the instrument. They are listed in Table A.1 together with their prior and the recovered best-fit value for the TTTEEE combination.

Appendix B: Best-fit model components

Here, we present our results for the best-fitting model components for each cross-power spectrum. These are shown in Fig. B.1 and the corresponding χ² values are given in Table B.1.

Fig. B.1. Best-fit model for each cross-frequency power spectrum in temperature, including emission from CMB, dust, tSZ, kSZ, CIB, SZ×CIB, and Poisson-noise from radio sources and dusty galaxies. Negative components are shown as dashed lines. Vertical black dashed lines show the range of multipoles considered in HiLLiPoP V4.2. The bottom panels show the residuals normalized by the error bars. Data are binned with Δℓ = 20 for this plot.

All Tables

All Figures

	Fig. 1. Uncertainties on each angular cross-power spectrum (blue lines) and their combination (red line) for the PlanckTT (top), TE (middle), and EE (bottom) data, compared to sample variance for 75% of the sky (black dashed line).
In the text

	Fig. 2. Detset cross-spectra for half-ring differences computed on 75% of the sky, divided by their uncertainties. From top to bottom we show: TT, EE, TE, and ET. Spectra are binned with Δℓ = 40. The projections on the right show the distribution for each unbinned spectrum over the range ℓ = 30–2500.
In the text

	Fig. 3. EE power spectrum of the CMB computed on 50% of the sky with the PR4 maps at low multipoles (Tristram et al. 2021). The Planck 2018 ΛCDM model is plotted in black. The grey band represents the associated sample variance. Error bars are deduced from the PR4 Monte Carlo simulations.
In the text

	Fig. 4. Sky masks used for HiLLiPoP V4.2 as a combination of a Galactic mask (blue, green, and red for the 100, 143, and 217 GHz channel, respectively), a CO mask, a point-source mask, and a mask removing nearby galaxies. The effective sky fractions remaining at 100, 143 and 217 GHz are 75%, 66%, and 52%, respectively.
In the text

	Fig. 5. Frequency cross-power spectra with respect to the mean spectra for TT, EE, TE, and ET. Spectra are binned with Δℓ = 40 for this figure.
In the text

	Fig. 6. Full HiLLiPoP covariance matrix, including all correlations in multipoles between cross-frequencies and power spectra.
In the text

Fig. 7. Dust power spectra, Dℓ = ℓ(ℓ + 1)Cℓ/2π, at 353 GHz for TT (top), EE (middle), and TE (bottom). The power spectra are computed from a cross-correlation between the detset maps at 353 GHz for different sets of masks, as defined in Sect. 3.2.1, and further corrected for the CMB power spectrum (solid black line) and CIB power spectrum (dashed black line). The coloured dashed lines are simple fits, as described in the text.

In the text

	Fig. 8. Maximum-likelihood frequency-co-added temperature power spectrum for HiLLiPoP V4.2. For the purposes of this figure, the power spectrum is binned with Δℓ = 30. The middle panel shows the residuals with respect to the fiducial base-ΛCDM cosmology and the bottom panel shows the residuals normalized by the uncertainties.
In the text

	Fig. 9. As in Fig. 8, but for TE (top) and EE (bottom) power spectra.
In the text

	Fig. 10. Posterior distributions for the cosmological parameters using power spectra from Planck PR4 with TT (lowT+hlpTT), TE (lowT+lolE+hlpTE), EE (lolE+hlpEE), and TTTEEE (lowT+lolE+hlpTTTEEE).
In the text

	Fig. 11. Correlation matrix for the fitted parameters of the combined HiLLiPoP likelihood TTTEEE. The first block corresponds to cosmological parameters from the ΛCDM model, the second block gathers the foreground parameters, and the last block shows the instrumental parameters.
In the text

	Fig. 12. Posteriors for foreground amplitudes. Units are μK2 normalized at ℓ = 3000 and ν = 143 GHz.
In the text

	Fig. 13. Amplitude of the dust emission relative to 353 GHz for a modified-blackbody dust model (blue line) as a function of the effective frequency (computed as the geometric mean of the two frequencies involved), compared to a fit using one amplitude per frequency (black dots). The top panel is for temperature and the bottom panel for polarization.
In the text

	Fig. 14. Point-source model as a function of the effective frequency (computed as the geometric mean of the two frequencies involved), compared to the fit of one amplitude per cross-spectrum.
In the text

	Fig. 15. Posterior distributions for the cosmological parameters from PR4 for HiLLiPoP (using different sky fractions labelled L, XL, and XXL) and CamSpec, as compared to Planck 2018 (Plik PR3). Likelihoods are considered for the combination of TT+TE+EE, with lowT and lolE used at low ℓ.
In the text

	Fig. 16. Posterior distributions for some parameters using TTTEEE in combination with lensing and BAO.
In the text

	Fig. 17. Posterior distributions for AL.
In the text

	Fig. 18. Posterior distributions for AL from HiLLiPoP PR4, compared to CamSpec (PR4) and Plik (PR3).
In the text

	Fig. 19. Posterior distributions in the ΩK–H0 plane using HiLLiPoP PR4, compared to CamSpec (PR4) and Plik (PR3).
In the text

	Fig. 20. Posterior distributions in the ΩK–H0 plane using Planck PR4 TTTEEE (i.e., lowT+lolE+hlpTTTEEE) in combination with lensing and BAO.
In the text

	Fig. 21. Posterior distributions for Neff. The vertical dashed line shows the theoretical expectation (Neff = 3.044).
In the text

	Fig. 22. Posterior distributions for ∑mν. Units are electronvolts.
In the text

	Fig. 23. Posterior distributions in the ∑mν–τ plane using Planck PR4 TTTEEE (i.e., lowT+lolE+hlpTTTEEE) in combination with lensing and BAO.
In the text

Fig. B.1. Best-fit model for each cross-frequency power spectrum in temperature, including emission from CMB, dust, tSZ, kSZ, CIB, SZ×CIB, and Poisson-noise from radio sources and dusty galaxies. Negative components are shown as dashed lines. Vertical black dashed lines show the range of multipoles considered in HiLLiPoP V4.2. The bottom panels show the residuals normalized by the error bars. Data are binned with Δℓ = 20 for this plot.

In the text