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

HAL
open science

# Cosmological parameters derived from the final (PR4) Planck data release

M. Tristram, A.J Banday, M Douspis, X. Garrido, K.M Górski, S.

Henrot-Versille, L.T Hergt, S Ilić, R Keskitalo, G Lagache, et al.

► To cite this version:

M. Tristram, A.J Banday, M Douspis, X. Garrido, K.M Górski, et al. Cosmological parameters derived from the final (PR4) Planck data release. Astronomy & Astrophysics - A&A, 2024, 682, pp.A37. 10.1051/0004-6361/202348015. hal-04231374

HAL Id: hal-04231374

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

Submitted on 6 Feb 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Cosmological parameters derived from the final *Planck* data release (PR4)

M. Tristram<fnref n="1" />, A. J. Banday<fnref n="2" />, M. Douspis<fnref n="3" />, X. Garrido<fnref n="1" />, K. M. Górski<fnref n="4" />,<fnref n="5" />, S. Henrot-Versillé<fnref n="1" />,
L. T. Hergt<fnref n="6" />, S. Ilić<fnref n="1" />,<fnref n="7" />, R. Keskitalo<fnref n="8" />,<fnref n="9" />, G. Lagache<fnref n="10" />, C. R. Lawrence<fnref n="4" />, B. Partridge<fnref n="11" />, and D. Scott<fnref n="6" />

Received 15 September 2023 / Accepted 29 October 2023

## ABSTRACT

We present cosmological parameter constraints using maps from the last *Planck* data release (PR4). In particular, we detail an upgraded version of the cosmic microwave background likelihood, HiLLiPoP, that is based on angular power spectra and relies on a physical modeling of the foreground residuals in the spectral domain. This new version of the likelihood retains a larger sky fraction (up to 75%) and uses an extended multipole range. Using this likelihood, along with low-$\ell$ measurements from LoLLiPoP, we derived constraints on $\Lambda$CDM parameters that are in good agreement with previous *Planck* 2018 results, but with smaller uncertainties by 10% to 20%. We demonstrate that the foregrounds can be accurately described in the spectral domain, with a negligible impact on $\Lambda$CDM parameters. We also derived constraints on single-parameter extensions to $\Lambda$CDM, including $A_L$, $\Omega_K$, $N_{\text{eff}}$, and $\sum m_\nu$. Noteworthy results from this updated analysis include a lensing amplitude value of $A_L = 1.039 \pm 0.052$, which is more closely aligned with theoretical expectations within the $\Lambda$CDM framework. Additionally, our curvature measurement, $\Omega_K = -0.012 \pm 0.010$, is now fully consistent with a flat universe and our measurement of $S_8$ is closer to the measurements derived from large-scale structure surveys (at the $1.5\sigma$ level). We also added constraints from PR4 lensing, making this combination the most tightly constrained data set currently available from *Planck*. Additionally, we explored the addition of baryon acoustic oscillation data, which tightens the limits on some particular extensions to the standard cosmology.

**Key words.** cosmic background radiation – methods: data analysis – cosmological parameters – cosmology: observations

## 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 ($\Lambda$CDM) model, given tight constraints on its parameters, and ruled out many plausible extensions. As a consequence, the best-fitting six-parameter $\Lambda$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-$\ell$ and high-$\ell$ likelihoods based on *Planck* PR4. The only part that still relies on PR3 (also known as *Planck* 2018) is the low-$\ell$ 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-$\ell$ likelihood, namely, the $LOw-\ell$ Likelihood on

<footnote n="1">¹ Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France
e-mail: tristram@ijclab.in2p3.fr</footnote>

<footnote n="2">² IRAP, Université de Toulouse, CNRS, CNES, UPS, Toulouse, France</footnote>

<footnote n="3">³ Université Paris-Saclay, CNRS, Institut d'Astrophysique Spatiale, 91405 Orsay, France</footnote>

<footnote n="4">⁴ Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA, USA</footnote>

<footnote n="5">⁵ Nicolaus Copernicus Academy and Superior School, Ul. Nowogrodzka 47A, 00-695 Warszawa, Poland</footnote>

<footnote n="6">⁶ Department of Physics & Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver,
British Columbia V6T 1Z1, Canada</footnote>

<footnote n="7">⁷ Centre National d'Études Spatiales – Centre Spatial de Toulouse, 18 Avenue Édouard Belin, 31401 Toulouse Cedex 9, France</footnote>

<footnote n="8">⁸ Computational Cosmology Center, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA</footnote>

<footnote n="9">⁹ Space Sciences Laboratory, University of California, Berkeley, CA 94720, USA</footnote>

<footnote n="10">¹⁰ Aix-Marseille Université, CNRS, CNES, LAM, Marseille, France</footnote>

<footnote n="11">¹¹ Department of Astronomy, Haverford College, Haverford, PA 19041, USA</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Uncertainties on each angular cross-power spectrum (blue lines) and their combination (red line) for the *Planck* *TT* (top), *TE* (middle), and *EE* (bottom) data, compared to sample variance for 75% of the sky (black dashed line).</figcaption></figure>

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-$\ell$ 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 ($\ell$ < 2000), *Planck* PR4 is the most sensitive data set for CMB anisotropies as of today.

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-$\ell$ likelihood (30 < $\ell$ < 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 $\Lambda$CDM parameters are insignificantly affected by the details of the foreground modeling, the constraints on extensions to $\Lambda$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_\ell$ 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 $\Lambda$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 $\Lambda$CDM in Sect. 11, specifically the lensing consistency parameter, $A_L$, the curvature, $\Omega_K$, the effective number of neutrino species, $N_{\text{eff}}$, and the sum of neutrino masses, $\sum m_\nu$.

## 2. The Planck PR4 data set

The *Planck* sky measurements used in this analysis are the PR4 maps available from the *Planck* Legacy Archive<fnref n="1" /> (PLA) and from the National Energy Research Scientific Computing Center (NERSC)<fnref n="2" />. 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 systems 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

<footnote n="1">¹ pla.esac.esa.int</footnote>

<footnote n="2">² portal.nersc.gov/project/cmb/planck2020</footnote>
<figure><img src="image_4.png" /><figcaption>**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 $\Delta\ell = 40$. The projections on the right show the distribution for each unbinned spectrum over the range $\ell = 30–2500$.</figcaption></figure>

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 $\ell = 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.

# 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_\ell$ 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.

<figure><img src="image_5.png" /><figcaption>**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.</figcaption></figure>

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_{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 $\ell = 30$ (for more details, see Tristram et al. 2021). We also performed the same estimation on each of the PR4 simulations and derived the $\ell$-by-$\ell$ covariance matrix that was then used to propagate uncertainties in LoLLiPoP, the low-$\ell$ CMB likelihood described in Sect. 4.

## 3.2. Small-scale power spectra

### 3.2.1. Sky fractions

For small scales ($\ell > 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 $\ell \approx 1800$ in temperature (and $\ell \approx 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 extra-galactic 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).

### 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

<figure><img src="image_6.png" /><figcaption>**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.</figcaption></figure>

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<fnref n="3" />. 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 $\tilde{a}_{\ell m}^X = \{\tilde{a}_{\ell m}^T, \tilde{a}_{\ell m}^E, \tilde{a}_{\ell m}^B\}$, we form the pseudo cross-power spectra between map i and map j,

$$ \tilde{C}_{\ell}^{ij} = \frac{1}{2\ell + 1} \sum_{m} \tilde{a}_{\ell m}^{i*} \tilde{a}_{\ell m}^{j}, \qquad (1) $$

where the vector $\tilde{C}_\ell$ includes the four modes $\{\tilde{C}_\ell^{TT}, \tilde{C}_\ell^{EE}, \tilde{C}_\ell^{TE}, \tilde{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),

$$ \tilde{C}_{\ell}^{ij} = (2\ell' + 1)M_{\ell\ell'}^{ij}C_{\ell'}^{ij}. \qquad (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

<footnote n="3">³ gitlab.in2p3.fr/tristram/Xpol</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Frequency cross-power spectra with respect to the mean spectra for *TT*, *EE*, *TE*, and *ET*. Spectra are binned with $\Delta\ell = 40$ for this figure.</figcaption></figure>

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.

### 3.2.3. Multipole ranges

The HiLLiPoP likelihood covers the multipoles starting from $\ell_{min} = 30$ up to $\ell_{max} = 2500$ in temperature and $\ell_{max} = 2000$ in polarization. The multipoles below $\ell < 30$ are considered in the low-$\ell$ likelihoods (Commander and LoLLiPoP, see Sect. 4).

Table 1 gives the HiLLiPoP multipole ranges, [$\ell_{min}$, $\ell_{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 $\ell$ and instrumental noise at high $\ell$. 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 $\ell = 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 $\ell \approx 1800$ in *TT* and $\ell \approx 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_\ell$ covariance matrix computed using Xpol. The matrix captures the $\ell$-by-$\ell$ 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.

<table><caption>**Table 1.** Multipole ranges used in the HiLLiPoP analysis and corresponding number of $\ell$s available ($n_\ell = \ell_{max} - \ell_{min} + 1$).</caption><thead><tr><th>Channels</th><th>TT</th><th>TE</th><th>EE</th></tr></thead><tbody><tr><td>100 × 100</td><td>[30,1500]</td><td>[30,1500]</td><td>[100,1200]</td></tr><tr><td>100 × 143</td><td>[30,1500]</td><td>[30,1500]</td><td>[30,1500]</td></tr><tr><td>100 × 217</td><td>[250,1500]</td><td>[100,1500]</td><td>[250,1500]</td></tr><tr><td>143 × 143</td><td>[50,2000]</td><td>[30,2000]</td><td>[30,2000]</td></tr><tr><td>143 × 217</td><td>[250,2500]</td><td>[200,2000]</td><td>[250,2000]</td></tr><tr><td>217 × 217</td><td>[250,2500]</td><td>[300,2000]</td><td>[250,2000]</td></tr><tr><td></td><td>10 646</td><td>9816</td><td>9296</td></tr></tbody></table>

**Notes.** The total number of $\ell$s across all spectra is 29 758.

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

$$ \Sigma_{\ell_1 \ell_2}^{ab,cd} = \langle \Delta C_{\ell}^{ab} \Delta C_{\ell'}^{cd*} \rangle = (M_{\ell \ell_1}^{ab})^{-1} \tilde{\Sigma}_{\ell_1 \ell_2}^{ab,cd} (M_{\ell' \ell_2}^{cd*})^{-1}, \quad (3) $$

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

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

$$ \tilde{\Sigma}_{\ell\ell'}^{ab,cd} = \langle \Delta \tilde{C}_{\ell}^{ab} \Delta \tilde{C}_{\ell'}^{cd*} \rangle = \langle \tilde{C}_{\ell}^{ab} \tilde{C}_{\ell'}^{cd*} \rangle - \tilde{C}_{\ell}^{ab} \tilde{C}_{\ell'}^{cd*} \\ = \sum_{mm'} \frac{\langle \tilde{a}_{\ell m}^{a} \tilde{a}_{\ell' m'}^{c*} \rangle \langle \tilde{a}_{\ell m}^{b*} \tilde{a}_{\ell' m'}^{d} \rangle + \langle \tilde{a}_{\ell m}^{a} \tilde{a}_{\ell' m'}^{d*} \rangle \langle \tilde{a}_{\ell m}^{b*} \tilde{a}_{\ell' m'}^{c} \rangle}{(2\ell + 1)(2\ell' + 1)}, $$

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

Each two-point correlation of pseudo-$a_{\ell m}$s can be expressed as the convolution of $C_\ell$ with a kernel that depends on the
polarization mode considered:

$$ \langle \tilde{a}_{\ell m}^{T_a} \tilde{a}_{\ell' m'}^{T_b} \rangle = \sum_{\ell_1 m_1} C_{\ell_1}^{T_a T_b} W_{\ell m \ell_1 m_1}^{0, T_a} W_{\ell' m' \ell_1 m_1}^{0, T_b^*}, $$

$$ \langle \tilde{a}_{\ell m}^{E_a^*} \tilde{a}_{\ell' m'}^{E_b} \rangle = \frac{1}{4} \sum_{\ell_1 m_1} \left\{ C_{\ell_1}^{E_a E_b} W_{\ell m \ell_1 m_1}^{+ , E_a^*} W_{\ell' m' \ell_1 m_1}^{+ , E_b} + C_{\ell_1}^{B_a B_b} W_{\ell m \ell_1 m_1}^{- , E_a^*} W_{\ell' m' \ell_1 m_1}^{- , E_b} \right\}, $$

$$ \langle \tilde{a}_{\ell m}^{T_a^*} \tilde{a}_{\ell' m'}^{E_b} \rangle = \frac{1}{2} \sum_{\ell_1 m_1} C_{\ell_1}^{T_a E_b} W_{\ell m \ell_1 m_1}^{0, T_a^*} W_{\ell' m' \ell_1 m_1}^{+ , E_b}, $$

where the kernels $W^0$, $W^+$, and $W^-$ are defined as linear combinations of products of $Y_{\ell m}$ of spin 0 and $\pm 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 $W$s into kernels similar to the coupling matrix $M$ defined in Eq. (2). In the end, the blocks of the $\Sigma$ matrices are:

$$ \begin{aligned} \Sigma_{T_a T_b, T_c T_d} &\approx C_{\ell \ell'}^{T_a T_c} C_{\ell \ell'}^{T_b T_d} M_{TT,TT} + C_{\ell \ell'}^{T_a T_d} C_{\ell \ell'}^{T_b T_c} M_{TT,TT}, \\ \Sigma_{E_a E_b, E_c E_d} &\approx C_{\ell \ell'}^{E_a E_c} C_{\ell \ell'}^{E_b E_d} M_{EE,EE} + C_{\ell \ell'}^{E_a E_d} C_{\ell \ell'}^{E_b E_c} M_{EE,EE}, \\ \Sigma_{T_a E_b, T_c E_d} &\approx C_{\ell \ell'}^{T_a T_c} C_{\ell \ell'}^{E_b E_d} M_{TE,TE} + C_{\ell \ell'}^{T_a E_d} C_{\ell \ell'}^{E_b T_c} M_{TT,TT}, \\ \Sigma_{T_a T_b, T_c E_d} &\approx C_{\ell \ell'}^{T_a T_c} C_{\ell \ell'}^{T_b E_d} M_{TT,TT} + C_{\ell \ell'}^{T_a E_d} C_{\ell \ell'}^{T_b T_c} M_{TT,TT}, \\ \Sigma_{T_a T_b, E_c E_d} &\approx C_{\ell \ell'}^{T_a E_c} C_{\ell \ell'}^{T_b E_d} M_{TT,TT} + C_{\ell \ell'}^{T_a E_d} C_{\ell \ell'}^{T_b E_c} M_{TT,TT}, \\ \Sigma_{E_a E_b, T_c E_d} &\approx C_{\ell \ell'}^{E_a T_c} C_{\ell \ell'}^{E_b E_d} M_{TE,TE} + C_{\ell \ell'}^{E_a E_d} C_{\ell \ell'}^{E_b T_c} M_{TE,TE}, \end{aligned} $$

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_\ell$ estimates and their covariance, we used a smoothed version of each measured power spectrum (using a Gaussian filter with $\sigma_\ell = 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.

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.

<figure><img src="image_8.png" /><figcaption>Fig. 6. Full HiLLiPoP covariance matrix, including all correlations in multipoles between cross-frequencies and power spectra.</figcaption></figure>

## 4. Large-scale CMB likelihoods: LoLLiPoP and Commander

The Planck low-$\ell$ 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_\ell \rightarrow X_\ell$, so that the new variable $X_\ell$ is nearly Gaussian-distributed. Similarly to Hamimeche & Lewis (2008), we define

$$ X_{\ell} = \sqrt{C_{\ell}^{T} + O_{\ell}} g\left(\frac{\tilde{C}_{\ell} + O_{\ell}}{C_{\ell} + O_{\ell}}\right) \sqrt{C_{\ell}^{T} + O_{\ell}}, \quad (4) $$

where $g(x) = \sqrt{2(x - \ln(x) - 1)}$, $\tilde{C}_\ell$ are the measured cross-power spectra, $C_\ell$ are the power spectra of the model to be evaluated, $C_\ell^T$ is a fiducial CMB model, and $O_\ell$ are the offsets needed in the case of cross-spectra. In the case of auto-power spectra, the offsets, $O_\ell$, 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_\ell$ noise variance:

$$ \Delta C_{\ell} \equiv \sqrt{\frac{2}{2\ell + 1} O_{\ell}}. \quad (5) $$

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

$$ -2 \ln P(C_{\ell} | \tilde{C}_{\ell}) = \sum_{\ell \ell'} X_{\ell}^{T} M_{\ell \ell'}^{-1} X_{\ell'}. \quad (6) $$
Uncertainties are incorporated into the $C_{\ell}$ covariance matrix $M_{\ell\ell'}$, 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_{\ell}$ 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<fnref n="4" />. In this work, we only considered the information from E modes and restricted the multipole range from $\ell = 2$ to $\ell = 30$.

To cover the low multipoles ($\ell < 30$) in the temperature, we made use of the Commander TT 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-$\ell$ 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 $\ell$-by-$\ell$ 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-$\chi^2$ with $2\ell + 1$ degrees of freedom. The distribution of the cross-spectra is slightly different (see Appendix A in Mangilli et al. 2015); however, above $\ell = 30$, the number of modes is large enough that we can safely assume that the $\tilde{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_{\ell}$ arising from the cut-sky, and is expressed as:

$$-2 \ln \mathcal{L} = \sum_{\substack{i \le j \\ i' \le j'}} \sum_{\ell \ell'} \mathbf{R}_{\ell}^{ij} [\Sigma^{-1}]_{\ell \ell'}^{ij i' j'} \mathbf{R}_{\ell'}^{i' j'} + \ln |\Sigma|, \quad (7)$$

where $\mathbf{R}_{\ell}^{ij} = \tilde{C}_{\ell}^{ij} - C_{\ell}^{ij}$ denotes the residual of the estimated cross-power spectrum $\tilde{C}_{\ell}$ with respect to the model $C_{\ell}$, which depends on the frequencies {$i, j$} and is described in the next section. The matrix $\Sigma = \langle \mathbf{RR}^T \rangle$ 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 $\Sigma$ 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 $\ell$ by less than they differ from the observed $\tilde{C}_{\ell}$, so that we can expand $\Sigma$ 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 $\Sigma$, we can drop the constant term $\ln|\Sigma|$ and recover nearly optimal variance (see Carron 2013). Within the approximations discussed above, we expect the likelihood to be $\chi^2$-distributed with a mean equal to the number of degrees of freedom $n_{\text{d.o.f.}} = n_{\ell} - n_{\text{p}}$ ($n_{\ell}$ being the number of band powers in the power spectra and $n_{\text{p}}$ the number of fitted parameters) and a variance equal to $2n_{\text{d.o.f.}}$.

## 5.2. The model

We now present the model ($\hat{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_{\ell} = \ell(\ell + 1)C_{\ell}/2\pi$. In the foreground, 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 band-pass 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

<footnote n="4">⁴ github.com/planck-npipe/lollipop</footnote>

<footnote n="5">⁵ github.com/planck-npipe/hillipop</footnote>
<figure><img src="image_9.png" /><figcaption>Fig. 7. Dust power spectra, $D_{\ell} = \ell(\ell + 1)C_{\ell}/2\pi$, 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.</figcaption></figure>

practice, we compute the 353-GHz cross-spectra, $\hat{C}_{\ell}^{353A\times353B}$, for each pair of masks ($M_i, M_j$) associated with the cross-spectra $v_i \times v_j$ (Fig. 7). We then subtract the Planck best-fit CMB power spectrum and fit a power-law model with a free constant $A\ell^{\alpha_d} + B$, in the range of $\ell$ = [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 $\ell$ = [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 $\ell$ = 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 $\alpha_d$ = -2.3 and -2.4 for EE and TE, respectively, the indices for TT vary with the sky fraction considered, ranging from $\alpha_d$ = -2.2 down to -2.6 for the largest sky fraction.

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:

$$D_{\ell}^{\text{dust}}(v \times v') = c_{\text{dust}} \frac{a_{v}^{\text{dust}}}{a_{353}^{\text{dust}}} \frac{a_{v'}^{\text{dust}}}{a_{353}^{\text{dust}}} D_{\ell}^{\text{dust}}(M_{v}, M_{v'}), \quad (8)$$

where $a_{v}^{\text{dust}} = v^{\beta_d} B_{v}(T_d)$ is a modified blackbody with $T_d$ fixed to 19.6 K, while $c_{\text{dust}}$ and $\beta_d$ are sampled independently for temperature and polarization. We use Gaussian priors for the spectral indices $\beta_d$ from Planck Collaboration Int. XXII (2015), which

gives $\beta_d^T = N(1.51, 0.01)$ and $\beta_d^T = N(1.59, 0.02)$ for the temperature and polarization, respectively. The coefficient $c_{\text{dust}}$ allows us to propagate the uncertainty from fitting the 353-GHz dust spectrum with a power law. We sample $c_{\text{dust}}$ with a Gaussian prior, $c_{\text{dust}} = N(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^{\text{CIB}}$, the amplitude of the contamination at our reference frequency ($v_0$ = 143 GHz) and $\ell$ = 3000. The emission law is modelled by a modified blackbody $a_{v}^{\text{CIB}} = v^{\beta_{\text{CIB}}} B_{v}(T)$ with a fixed temperature ($T$ = 25 K) and a variable index $\beta_{\text{CIB}}$. We use a strong prior $\beta_{\text{CIB}} = N(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),

$$D_{\ell}^{\text{CIB}}(v \times v') = A^{\text{CIB}} \frac{a_{v}^{\text{CIB}}}{a_{v_0}^{\text{CIB}}} \frac{a_{v'}^{\text{CIB}}}{a_{v_0}^{\text{CIB}}} D_{\ell}^{\text{CIB}}. \quad (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^{\text{tSZ}}$, corresponding to the amplitude of the tSZ signal at our reference frequency ($v_0$ = 143 GHz) at $\ell$ = 3000,

$$D_{\ell}^{\text{tSZ}}(v \times v') = A^{\text{tSZ}} \frac{a_{v}^{\text{tSZ}}}{a_{v_0}^{\text{tSZ}}} \frac{a_{v'}^{\text{tSZ}}}{a_{v_0}^{\text{tSZ}}} D_{\ell}^{\text{tSZ}}, \quad (10)$$

where $a_{v}^{\text{tSZ}} = x[e^x + 1]/[e^x - 1] - 4$ (with $x = h\nu/k_B T_{\text{CMB}}$).

**Kinetic Sunyaev-Zeldovich (kSZ) effect.** The kSZ emission is parameterized by $A^{\text{kSZ}}$, the amplitude at $\ell$ = 3000, scaling a fixed template that includes homogeneous and patchy reionization components from Shaw et al. (2012) and Battaglia et al. (2013),

$$D_{\ell}^{\text{kSZ}}(v \times v') = A^{\text{kSZ}} D_{\ell}^{\text{kSZ}}. \quad (11)$$

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

$$D_{\ell}^{\text{tSZ} \times \text{CIB}}(v \times v') = -\xi \sqrt{A^{\text{tSZ}} A^{\text{CIB}}} \times \left( \frac{a_{v}^{\text{tSZ}} a_{v'}^{\text{CIB}} + a_{v}^{\text{CIB}} a_{v'}^{\text{tSZ}}}{a_{v_0}^{\text{tSZ}} a_{v_0}^{\text{CIB}}} \right) D_{\ell}^{\text{tSZ} \times \text{CIB}}, \quad (12)$$

with $\xi$ the correlation coefficient rescaling the template $D_{\ell}^{\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_v^{\text{rad}} \propto v^{-\beta_s}$ (Tucci et al. 2011), while we use $a_v^{\text{IR}} = v^{\beta_{\text{IR}}} B_v(T)$ (Béthermin et al. 2012) for infrared dusty star-forming galaxies. The residual cross-power spectra for point sources are:

$$C_{\ell}^{\text{PS}}(\nu \times \nu') = A^{\text{rad}} \frac{a_{\nu}^{\text{rad}} a_{\nu'}^{\text{rad}}}{a_{\nu 0}^{\text{rad}} a_{\nu 0}^{\text{rad}}} + A^{\text{IR}} \frac{a_{\nu}^{\text{IR}} a_{\nu'}^{\text{IR}}}{a_{\nu 0}^{\text{IR}} a_{\nu 0}^{\text{IR}}}. \quad (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 $\beta_s = -0.8$, following results from Reichardt et al. (2021). For infrared dusty star-forming galaxies, we adopt $\beta_{\text{IR}}$ identical to $\beta_{\text{CIB}}$ and $T = 25$ K. The $C_{\ell s}$ are then converted into $D_{\ell s}$ such that the amplitudes $A^{\text{rad}}$ and $A^{\text{IR}}$ refer to the amplitude of $D_{3000}$ 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_{\text{Planck}}$ with a strong prior, $A_{\text{Planck}} = N(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-$\ell$ likelihood (lowT, see Planck Collaboration IV 2020) and the polarized low-$\ell$ EE likelihood LoLLiPoP (loIe, discussed in Sect. 4). We define the following combination of likelihoods for the rest of the paper:

- TT, lowT+hlpTT;

- TE, lowT+loIe+hlpTE;

- EE, loIe+hlpEE;

- TTTEE, lowT+loIe+hlpTTTEE.

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-$\ell$ 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)<fnref n="6" />. In this paper, we consider a ΛCDM model with six free parameters describing: the current physical densities of baryons ($\Omega_b h^2$) and cold dark matter ($\Omega_c h^2$); the angular acoustic scale ($\theta_*$); the reionization optical depth ($\tau$); 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_0/(100 \text{ km s}^{-1} \text{ Mpc}^{-1})$.

In addition, we fit six inter-calibration parameters, seven foreground residual amplitudes in the temperature ($c_{\text{dust}}^T$, $A_{\text{radio}}$, $A_{\text{IR}}$, $A_{\text{CIB}}$, $A_{\text{SZ}}$, $A_{\text{KSZ}}$, and $\xi_{\text{SZ}\times\text{CIB}}$), plus one in polarization ($c_{\text{dust}}^P$, as well as three foreground spectral indices ($\beta_{\text{dust}}^T$, $\beta_{\text{dust}}^P$, and $\beta_{\text{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 $\chi^2$ 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 $\chi^2$ 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 $\chi^2$ value for hlpTT (corresponding to about 2.7$\sigma$), while the hlpTE and hlpEE $\chi^2$ values are compatible with unity, at 1.8$\sigma$ and 0.1$\sigma$, respectively. For the full combination hlpTTTEE, we obtained $\chi^2 = 30495$ for a data size of 29768, corresponding to a 3.02$\sigma$ deviation. As described in Rosenberg et al. (2022), where the goodness of fit is also somewhat poor (4.07$\sigma$ for TT and 4.46$\sigma$ for the TTTEE), 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

<footnote n="6">⁶ One can equally well use CLASS (Blas et al. 2011) instead, except that the definition of $\theta_*$ differs slightly between the two codes.</footnote>
<table><caption>**Table 2.** $\chi^2$ values compared to the size of the data vector ($n_d$) for each of the Planck HiLLiPoP likelihoods.</caption><thead><tr><th>Likelihood</th><th>χ<sup>2</sup></th><th>n<sub>d</sub></th><th>χ<sup>2</sup>/n<sub>d</sub></th><th>δσ(χ<sup>2</sup>)</th></tr></thead><tbody><tr><td>hlpEE</td><td>9289</td><td>9296</td><td>1.00</td><td>0.05</td></tr><tr><td>hlpTE</td><td>10 071</td><td>9816</td><td>1.03</td><td>1.82</td></tr><tr><td>hlpTT</td><td>11 044</td><td>10 646</td><td>1.04</td><td>2.73</td></tr><tr><td>hlpTTTEEE</td><td>30 495</td><td>29 758</td><td>1.02</td><td>3.02</td></tr></tbody></table>

Notes. Here $\delta\sigma(\chi^2) = (\chi^2/n_d - 1)/\sqrt{2/n_d}$.

<figure><img src="image_10.png" /><figcaption>**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 $\Delta\ell = 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.</figcaption></figure>

reduced-$\chi^2$, $\chi^2/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 $\Delta\ell = 30$ for the plots, but considered $\ell$-by-$\ell$ 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.

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 $\theta_*$ 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 systems in EE. However, the deviation of $\theta_*$ 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\sigma$ level. In addition, we have checked that this shift in $\theta_*$ is not related to any super-sample lensing effect (as described in Manzotti et al.

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** As in Fig. 8, but for TE (top) and EE (bottom) power spectra.</figcaption></figure>

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, $\theta_*$ is the only parameter that deviates in EE; the others, including $H_0$, are compatible with TT at much better than $1\sigma$. Given the weak sensitivity of the Planck EE spectra as compared to TT and TE, discrepancies in the EE parameter reconstruction will have little impact on the overall cosmological parameter results.

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:

$$ \tau = 0.058 \pm 0.006. \qquad (14) $$

This is the result of the NPIPE treatment of the PR4 data associated with the low-$\ell$ likelihood LoLLiPoP (see Planck Collaboration Int. LVII 2020).
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Posterior distributions for the cosmological parameters using power spectra from *Planck* PR4 with TT (lowT+hlpTT), TE (lowT+loIe+hlpTE), EE (loIe+hlpEE), and TTTEE (lowT+loIe+hlpTTTEE).</figcaption></figure>

<table><caption>**Table 3.** Parameter constraints in the 6-parameter ΛCDM model for each data set and their combination, using HiLLiPoP V4.2 in addition to Commander and LoLLiPoP at low $\ell$.</caption><thead><tr><th>Parameter</th><th>TT</th><th>TE</th><th>EE</th><th>TTTEE</th></tr></thead><tbody><tr><td>Ω<sub>b</sub>h<sup>2</sup></td><td>0.02224 ± 0.00025</td><td>0.02238 ± 0.00020</td><td>0.02283 ± 0.00081</td><td>0.02226 ± 0.00013</td></tr><tr><td>Ω<sub>c</sub>h<sup>2</sup></td><td>0.1183 ± 0.0024</td><td>0.1172 ± 0.0017</td><td>0.1169 ± 0.0035</td><td>0.1188 ± 0.0012</td></tr><tr><td>100θ<sub>*</sub></td><td>1.04123 ± 0.00046</td><td>1.04151 ± 0.00041</td><td>1.04001 ± 0.00059</td><td>1.04108 ± 0.00026</td></tr><tr><td>log(10<sup>10</sup>A<sub>s</sub>)</td><td>3.073 ± 0.061</td><td>3.030 ± 0.028</td><td>3.077 ± 0.039</td><td>3.040 ± 0.014</td></tr><tr><td>n<sub>s</sub></td><td>0.9678 ± 0.0072</td><td>0.9689 ± 0.0091</td><td>0.9797 ± 0.0110</td><td>0.9681 ± 0.0039</td></tr><tr><td>τ</td><td>0.0753 ± 0.0322</td><td>0.0572 ± 0.0065</td><td>0.0582 ± 0.0066</td><td>0.0580 ± 0.0062</td></tr><tr><td>H<sub>0</sub></td><td>67.89 ± 1.11</td><td>68.49 ± 0.76</td><td>68.49 ± 1.91</td><td>67.64 ± 0.52</td></tr><tr><td>σ<sub>8</sub></td><td>0.8186 ± 0.0221</td><td>0.7973 ± 0.0129</td><td>0.8149 ± 0.0189</td><td>0.8070 ± 0.0065</td></tr><tr><td>S<sub>8</sub></td><td>0.826 ± 0.024</td><td>0.795 ± 0.021</td><td>0.814 ± 0.044</td><td>0.819 ± 0.014</td></tr><tr><td>Ω<sub>m</sub></td><td>0.3059 ± 0.0147</td><td>0.2983 ± 0.0099</td><td>0.2995 ± 0.0226</td><td>0.3092 ± 0.0070</td></tr></tbody></table>

**Notes.** We report mean values and symmetrical 68% confidence intervals.

For the constraint on the Hubble constant, we obtain:

$$H_0 = (67.64 \pm 0.52) \text{ km s}^{-1} \text{ Mpc}^{-1}, \qquad (15)$$

which is consistent with previous *Planck* results and still significantly lower than the local distance-ladder measurements, which typically range from $H_0 = 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:

$$\sigma_8 = 0.8070 \pm 0.0065, \qquad (16)$$

compatible with PR3 results ($\sigma_8 = 0.8120 \pm 0.0073$) but lower by $0.5\sigma$. The matter density, $\Omega_m$, also shifts by roughly $1\sigma$, such

that:

$$S_8 \equiv \sigma_8(\Omega_m/0.3)^{0.5} = 0.819 \pm 0.014. \qquad (17)$$

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

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
<figure><img src="image_13.png" /><figcaption>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.</figcaption></figure>

strengthens the robustness of the results with respect to the foreground model and ensures very low impact on cosmology.

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

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

$$A_{\text{radio}} = (63.3 \pm 4.7) \mu\text{K}^2, \qquad (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 \pm 0.15$, Reichardt et al. 2021). Allowing $β_s$ to vary in Planck data, gives $β_s = -0.54 \pm 0.08$, with a corresponding increase in the amplitude $A_{\text{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:

$$A_{\text{CIB}} = (1.03 \pm 0.34) \mu\text{K}^2, \qquad (19)$$

$$A_{\text{IR}} = (6.07 \pm 0.63) \mu\text{K}^2, \qquad (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

<figure><img src="image_14.png" /><figcaption>Fig. 12. Posteriors for foreground amplitudes. Units are μK² normalized at $ℓ = 3000$ and $ν = 143$ GHz.</figcaption></figure>

the $β_{\text{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:

$$A_{\text{tSZ}} = (5.9 \pm 1.7) \mu\text{K}^2. \qquad (21)$$

This is close to (but somewhat higher than) what is reported in Reichardt et al. (2021), with $A_{\text{tSZ}} = (3.42 \pm 0.54) \mu\text{K}^2$, even though the uncertainties are larger. However, it is more closely comparable with ACTpol results, $A_{\text{tSZ}} = (5.29 \pm 0.66) \mu\text{K}^2$ (Choi et al. 2020).

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

$$A_{\text{kSZ}} < 7.6 \mu\text{K}^2 \quad (\text{at 95\% CL}); \qquad (22)$$

$$\xi_{\text{SZ}\times\text{CIB}} = 0.46 \pm 0.30. \qquad (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_{\text{dust}}$, to be $1.08 \pm 0.03$ and $1.20 \pm 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, $β_{\text{dust}}^P = N(1.59, 0.02)$, we ended up recovering an amplitude that was higher than expected.
<figure><img src="image_15.png" /><figcaption>**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.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**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.</figcaption></figure>

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.

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:

$$c_{100A} = 1.003 \pm 0.007, \qquad (24)$$

$$c_{100B} = 1.004 \pm 0.007, \qquad (25)$$

$$c_{143B} = 1.004 \pm 0.006, \qquad (26)$$

$$c_{217A} = 1.001 \pm 0.008, \qquad (27)$$

$$c_{217B} = 1.001 \pm 0.008. \qquad (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\sigma$, except for the two maps at 217 GHz, which differ from unity by about $2\sigma$:

$$\eta_{100A} = 0.994 \pm 0.013; \qquad (29)$$

$$\eta_{100B} = 0.987 \pm 0.013; \qquad (30)$$

$$\eta_{143A} = 1.016 \pm 0.013; \qquad (31)$$

$$\eta_{143B} = 1.001 \pm 0.010; \qquad (32)$$

$$\eta_{217A} = 0.978 \pm 0.013; \qquad (33)$$

$$\eta_{217B} = 0.972 \pm 0.014. \qquad (34)$$

Fixing polarization efficiencies to 1.00, 1.00, and 0.975 (at 100, 143, and 217 GHz, respectively) increases the $\chi^2$ by $\Delta\chi^2 = 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 $\chi^2$ for the fits, demonstrating that the model used in HiLLiPoP V4.2 is valid for the considered sky fraction. For the TTTEEE likelihood, the $\Delta\chi^2$ 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\sigma$ 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-$\ell$ 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).

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 $\Omega_{c}h^2$, for which HiLLiPoP on PR4 finds a value $1.0\sigma$ lower. Associated with Commander and LoLLiPoP, CamSpec on PR4 also gives lower $\Omega_{c}h^2$ by $-0.45\sigma$. The spectral index $n_s$ is found to be a bit higher with HiLLiPoP by $0.7\sigma$.
<figure><img src="image_17.png" /><figcaption>**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 loIe used at low $\ell$.</figcaption></figure>

<table><caption>**Table 4.** Relative variation and improvement in the error bars between Planck 2018 and this work for each cosmological parameter.</caption><thead><tr><th>Parameter</th><th>Δ/σ</th><th>Δσ</th></tr></thead><tbody><tr><td>Ω<sub>b</sub>h<sup>2</sup></td><td>-0.67</td><td>-13.7%</td></tr><tr><td>Ω<sub>c</sub>h<sup>2</sup></td><td>-0.99</td><td>-15.2%</td></tr><tr><td>100θ<sub>*</sub></td><td>-0.01</td><td>-16.1%</td></tr><tr><td>log(10<sup>10</sup>A<sub>s</sub>)</td><td>-0.30</td><td>-12.0%</td></tr><tr><td>n<sub>s</sub></td><td>+0.75</td><td>-11.0%</td></tr><tr><td>τ</td><td>+0.44</td><td>-21.4%</td></tr><tr><td>H<sub>0</sub></td><td>+0.61</td><td>-13.7%</td></tr><tr><td>σ<sub>8</sub></td><td>-0.70</td><td>-11.5%</td></tr><tr><td>S<sub>8</sub></td><td>-0.89</td><td>-14.2%</td></tr><tr><td>Ω<sub>m</sub></td><td>-0.79</td><td>-16.1%</td></tr></tbody></table>

As discussed in Sect. 6, we obtain a slightly higher value for the Hubble constant (+0.6σ) with $h = 0.6766 \pm 0.0053$, compared to $h = 0.6727 \pm 0.0060$ for PR3. The amplitude of density fluctuations, $\sigma_8$, and the matter density, $\Omega_m$, are lower by 0.7σ and 0.8σ, respectively, so that $S_8$ 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 ($\theta_*$).

## 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

<figure><img src="image_18.png" /><figcaption>**Fig. 16.** Posterior distributions for some parameters using TTTEEE in combination with lensing and BAO.</figcaption></figure>

$C_{\ell}^{\phi\phi}$ 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-$\alpha$ auto, and Ly-$\alpha$QSO, 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 $\Lambda$CDM model when adding lensing and BAO data. Figure 16 shows the posterior distribution for the particular subset $\Omega_b h^2$, $\Omega_m$, $\sigma_8$, and $H_0$.

## 11. Extensions

We now discuss constraints on some extensions to the base-$\Lambda$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 $\Lambda$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:

$$A_L = 1.039 \pm 0.052, \qquad (35)$$

which is compatible with the $\Lambda$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
<table><caption>**Table 5.** Parameter constraints in the 6-parameter ΛCDM model for each data set and their combination, using HiLLiPoP V4.2 in addition to Commander and LoLLiPoP at low $\ell$, with the addition of CMB lensing and BAO constraints.</caption><thead><tr><th>Parameter</th><th>TTTEEE</th><th>TTTEEE +lensing</th><th>TTTEEE +lensing+BAO</th></tr></thead><tbody><tr><td>Ω<sub>b</sub>h<sup>2</sup></td><td>0.02226 ± 0.00013</td><td>0.02226 ± 0.00013</td><td>0.02229 ± 0.00012</td></tr><tr><td>Ω<sub>c</sub>h<sup>2</sup></td><td>0.1188 ± 0.0012</td><td>0.1190 ± 0.0011</td><td>0.1186 ± 0.0009</td></tr><tr><td>100θ<sub>*</sub></td><td>1.04108 ± 0.00026</td><td>1.04107 ± 0.00025</td><td>1.04111 ± 0.00024</td></tr><tr><td>log(10<sup>10</sup>A<sub>s</sub>)</td><td>3.040 ± 0.014</td><td>3.045 ± 0.012</td><td>3.048 ± 0.012</td></tr><tr><td>n<sub>s</sub></td><td>0.9681 ± 0.0039</td><td>0.9679 ± 0.0038</td><td>0.9690 ± 0.0035</td></tr><tr><td>τ</td><td>0.0580 ± 0.0062</td><td>0.0590 ± 0.0061</td><td>0.0605 ± 0.0059</td></tr><tr><td>H<sub>0</sub></td><td>67.64 ± 0.52</td><td>67.66 ± 0.49</td><td>67.81 ± 0.38</td></tr><tr><td>σ<sub>8</sub></td><td>0.8070 ± 0.0065</td><td>0.8113 ± 0.0050</td><td>0.8118 ± 0.0050</td></tr><tr><td>S<sub>8</sub></td><td>0.819 ± 0.014</td><td>0.824 ± 0.011</td><td>0.821 ± 0.009</td></tr><tr><td>Ω<sub>m</sub></td><td>0.3092 ± 0.0070</td><td>0.3092 ± 0.0066</td><td>0.3071 ± 0.0051</td></tr></tbody></table>

**Notes.** We report mean values and symmetrical 68% confidence intervals.

<table><caption>**Table 6.** Mean values and 68% confidence intervals for $A_L$.</caption><thead><tr><th>Likelihood</th><th>A<sub>L</sub></th><th>ΔA<sub>L</sub></th></tr></thead><tbody><tr><td>TT</td><td>1.075 ± 0.102</td><td>0.73σ</td></tr><tr><td>TE</td><td>0.937 ± 0.158</td><td>-0.40σ</td></tr><tr><td>EE</td><td>0.912 ± 0.150</td><td>-0.59σ</td></tr><tr><td>TTTEEE</td><td>1.039 ± 0.052</td><td>0.75σ</td></tr></tbody></table>

**Notes.** The significance of the deviation from unity is given in the last column.

<figure><img src="image_19.png" /><figcaption>**Fig. 17.** Posterior distributions for $A_L$.</figcaption></figure>

0.8σ. Figure 17 shows posterior distributions of $A_L$ for each of the mode-spectra and for the TTTEEE combination using Planck PR4.

In Rosenberg et al. (2022), the CamSpec likelihood associated with low-$\ell$ 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-$\ell$ polarized likelihood, instead of the low-$\ell$ likelihoods from Planck 2018, the constraint on $A_L$ from CamSpec changed from $A_L = 1.095 \pm 0.056$ to $A_L = 1.075 \pm 0.058$, still a 1.3σ difference from unity. We compare the posteriors for Plik (PR3), CamSpec (PR4), and HiLLiPoP (PR4) in Fig. 18.

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

<figure><img src="image_20.png" /><figcaption>**Fig. 18.** Posterior distributions for $A_L$ from HiLLiPoP PR4, compared to CamSpec (PR4) and Plik (PR3).</figcaption></figure>

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

$$A_L = 1.037 \pm 0.037 \quad (\text{TTTEEE}+\text{lensing}). \qquad (36)$$

## 11.2. Curvature, $\Omega_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., $\Omega_K < 0$). Indeed, with HiLLiPoP V4.2, the measurements are consistent with a flat universe ($\Omega_K = 0$) for all spectra.

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

$$\Omega_K = -0.012 \pm 0.010, \qquad (37)$$

which is only 1.2σ discrepant from zero.

This is to be compared to $\Omega_K = -0.044^{+0.018}_{-0.015}$ obtained for Plik on PR3 (Planck Collaboration VI 2020) and
<figure><img src="image_21.png" /><figcaption>Fig. 19. Posterior distributions in the $\Omega_K-H_0$ plane using HiLLiPoP PR4, compared to CamSpec (PR4) and Plik (PR3).</figcaption></figure>

$\Omega_K = -0.025^{+0.013}_{-0.010}$ obtained with CamSpec on PR4 (Rosenberg et al. 2022).

As a consequence, the tail of the 2-d posterior in the $H_0-\Omega_K$ plane at low $H_0$ and negative $\Omega_K$ is no longer favoured. Indeed, when fitting for a non-flat Universe, the recovered value for the Hubble constant is $H_0 = (63.03 \pm 3.60) \text{ km s}^{-1} \text{ Mpc}^{-1}$, only $1.3\sigma$ away from the constraint with fixed $\Omega_K = 0$.

The combination of TTTEEE with lensing yields the improved constraint

$$ \Omega_K = -0.0078 \pm 0.0058 \quad (\text{TTTEEE}+\text{lensing}). \quad (38) $$

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

$$ \Omega_K = 0.0000 \pm 0.0016 \quad (\text{TTTEEE}+\text{lensing}+\text{BAO}). \quad (39) $$

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

### 11.3. Effective number of relativistic species, $N_{\text{eff}}$

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

$$ N_{\text{eff}} = 3.08 \pm 0.17. \quad (40) $$

The uncertainties are comparable to Planck 2018 results ($N_{\text{eff}} = 2.92 \pm 0.19$, Planck Collaboration VI 2020) with a slight shift toward higher values, closer to the theoretical expectation $N_{\text{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_{\text{eff}} = 3.00 \pm 0.21$, Rosenberg et al. 2022).

<figure><img src="image_22.png" /><figcaption>Fig. 20. Posterior distributions in the $\Omega_K-H_0$ plane using Planck PR4 TTTEEE (i.e., lowT+loIe+hlpTTTEEE) in combination with lensing and BAO.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 21. Posterior distributions for $N_{\text{eff}}$. The vertical dashed line shows the theoretical expectation ($N_{\text{eff}} = 3.044$).</figcaption></figure>

### 11.4. Sum of the neutrino masses, $\sum m_\nu$

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

$$ \sum m_\nu < 0.39 \text{ eV} \quad (95\% \text{ CL, TTTEEE}). \quad (41) $$

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: $\sum m_\nu < 0.26 \text{ eV}$ at 95% CL. Our constraint is comparable to CamSpec, which gives $\sum m_\nu < 0.36 \text{ 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 $\sum m_\nu$ pushes the peak posterior of $\sum m_\nu$ 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 \pm 0.052$, see Sect. 11.1), thus, the posterior
<figure><img src="image_24.png" /><figcaption>**Fig. 22.** Posterior distributions for $\sum m_v$. Units are electronvolts.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Fig. 23.** Posterior distributions in the $\sum m_v-\tau$ plane using Planck PR4 TTTEE (i.e., lowT+loI+E+hlpTTTEE) in combination with lensing and BAO.</figcaption></figure>

for $\sum m_v$ is shifted to higher values, with a peak closer to zero, increasing the upper limit accordingly.

Figure 23 shows constraints in the $\sum m_v-\tau$ plane when combining our new likelihood with with CMB lensing and BAO data. This combination further strengthens the limits to:

$$ \sum m_v < 0.26 \text{ eV} \quad (95\%\text{CL, TTTEE}+\text{lensing}), \quad (42) $$

$$ \sum m_v < 0.11 \text{ eV} \quad (95\%\text{CL, TTTEE}+\text{lensing}+\text{BAO}). \quad (43) $$

This is slightly tighter than the upper limit from Planck 2018 ($\sum m_v < 0.12 \text{ eV}$) and getting close to the lower limit for the inverted mass hierarchy ($\sum m_v \ge 0.1 \text{ 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-$\ell$ 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 $\Lambda$CDM parameters is negligible.

Combined with the low-$\ell$ EE likelihood LoLLiPoP, we derived constraints on $\Lambda$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_8$, alleviating the CMB versus large-scale structure tension to $1.5\sigma$.

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 $\Lambda$CDM, including the lensing amplitude, $A_L$, the curvature, $\Omega_K$, the effective number of relativistic species, $N_{\text{eff}}$, and the sum of the neutrino masses, $\sum m_v$. For both $A_L$ and $\Omega_K$, our results show a significant reduction of the so-called "tensions" within the standard $\Lambda$CDM, along with a reduction of the uncertainties. Indeed, the final constraints are fully compatible with $\Lambda$CDM predictions. In particular, with the new version of the likelihood presented in this work, we report $A_L = 1.039 \pm 0.052$, which is entirely compatible with the $\Lambda$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.

**Acknowledgements.** 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 (gitlab.in2p3.fr/xQML); Xpol, for the computation of large-scale power spectra (gitlab.in2p3.fr/tristram/Xpol); Cobaya, for the sampling of the likelihoods (github.com/Cobaya Sampler); and CLASS (github.com/lesgourg/class_public) and CAMB (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 (alliancecan.ca). 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 (pla.esac.esa.int). Both likelihoods LoLLiPoP and HiLLiPoP based on PR4 are publicly available on GitHub (github.com/planck-npipe) as external likelihoods for Cobaya.

## References

* Abbott, T. M. C., Aguena, M., Alarcon, A., et al. 2022, Phys. Rev. D, 105, 023520

* Abdalla, E., Abellán, G. F., Aboubrahim, A., et al. 2022, J. High Energy Astrophys., 34, 49

* Addison, G. E., Dunkley, J., & Spergel, D. N. 2012, MNRAS, 427, 1741
Ade, P. A. R., Ahmed, Z., Amiri, M., et al. 2021, Phys. Rev. Lett., 127, 151301

Akita, K., & Yamaguchi, M. 2020, J. Cosmol. Astropart. Phys., 2020, 012

Alam, S., Aubert, M., Avila, S., et al. 2021, Phys. Rev. D, 103, 083533

Battaglia, N., Natarajan, A., Trac, H., Cen, R., & Loeb, A. 2013, ApJ, 776, 83

Bennett, J. J., Buldgen, G., de Salas, P. F., et al. 2021, J. Cosmol. Astropart. Phys., 2021, 073

Béthermin, M., Daddi, E., Magdis, G., et al. 2012, ApJ, 757, L23

Beutler, F., Blake, C., Colless, M., et al. 2011, MNRAS, 416, 3017

Blas, D., Lesgourgues, J., & Tram, T. 2011, J. Cosmol. Astropart. Phys., 2011, 034

Brown, M. L., Castro, P. G., & Taylor, A. N. 2005, MNRAS, 360, 1262

Carron, J. 2013, A&A, 551, A88

Carron, J., Mirmelstein, M., & Lewis, A. 2022, J. Cosmol. Astropart. Phys., 2022, 039

Choi, S. K., Hasselfield, M., Ho, S.-P. P., et al. 2020, J. Cosmol. Astropart. Phys., 2020, 045

Couchot, F., Henrot-Versillé, S., Perdereau, O., et al. 2017a, A&A, 606, A104

Couchot, F., Henrot-Versillé, S., Perdereau, O., et al. 2017b, A&A, 602, A41

Couchot, F., Henrot-Versillé, S., Perdereau, O., et al. 2017c, A&A, 597, A126

Efstathiou, G. P. 2006, MNRAS, 370, 343

Eriksen, H. K., Jewell, J. B., Dickinson, C., et al. 2008, ApJ, 676, 10

Froustey, J., Pitrou, C., & Volpe, M. C. 2020, J. Cosmol. Astropart. Phys., 2020, 015

Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759

Hamimeche, S., & Lewis, A. 2008, Phys. Rev. D, 77, 103013

Hivon, E., Górski, K. M., Netterfield, C. B., et al. 2002, ApJ, 567, 2

Hivon, E., Mottet, S., & Ponthieu, N. 2017, A&A, 598, A25

Howlett, C., Lewis, A., Hall, A., & Challinor, A. 2012, J. Cosmol. Astropart. Phys., 2012, 027

Jeong, D., Chluba, J., Dai, L., Kamionkowski, M., & Wang, X. 2014, Phys. Rev. D, 89, 023003

Jimenez, R., Pena-Garay, C., Short, K., Simpson, F., & Verde, L. 2022, J. Cosmol. Astropart. Phys., 2022, 006

Lagache, G., Béthermin, M., Montier, L., Serra, P., & Tucci, M. 2020, A&A, 642, A232

Lewis, A. 2013, Phys. Rev. D, 87, 103529

Lewis, A., & Bridle, S. 2002, Phys. Rev. D, 66, 103511

Lewis, A., Challinor, A., & Lasenby, A. 2000, ApJ, 538, 473

Mangilli, A., Plaszczyński, S., & Tristram, M. 2015, MNRAS, 453, 3174

Manzotti, A., Hu, W., & Benoit-Lévy, A. 2014, Phys. Rev. D, 90, 023003

Peebles, P. J. E. 1973, ApJ, 185, 413

Planck Collaboration XIII. 2011, A&A, 536, A13

Planck Collaboration IX. 2014, A&A, 571, A9

Planck Collaboration XV. 2014, A&A, 571, A15

Planck Collaboration XVI. 2014, A&A, 571, A16

Planck Collaboration XXX. 2014, A&A, 571, A30

Planck Collaboration XI. 2016, A&A, 594, A11

Planck Collaboration XIII. 2016, A&A, 594, A13

Planck Collaboration XXII. 2016, A&A, 594, A22

Planck Collaboration III. 2020, A&A, 641, A3

Planck Collaboration IV. 2020, A&A, 641, A4

Planck Collaboration V. 2020, A&A, 641, A5

Planck Collaboration VI. 2020, A&A, 641, A6

Planck Collaboration Int. VII. 2013, A&A, 550, A133

Planck Collaboration Int. XIX. 2015, A&A, 576, A104

Planck Collaboration Int. XXII. 2015, A&A, 576, A107

Planck Collaboration Int. XXX. 2016, A&A, 586, A133

Planck Collaboration Int. XLVII. 2016, A&A, 596, A108

Planck Collaboration Int. LI. 2017, A&A, 607, A95

Planck Collaboration Int. LVII. 2020, A&A, 643, A42

Reichardt, C. L., Patil, S., Ade, P. A. R., et al. 2021, ApJ, 908, 199

Rosenberg, E., Gratton, S., & Efstathiou, G. 2022, MNRAS, 517, 4620

Ross, A. J., Samushia, L., Howlett, C., et al. 2015, MNRAS, 449, 835

Sellentin, E., & Heavens, A. F. 2016, MNRAS, 456, L132

Shaw, L. D., Rudd, D. H., & Nagai, D. 2012, ApJ, 756, 15

Tanimura, H., Douspis, M., Aghanim, N., & Salvati, L. 2022, MNRAS, 509, 300

Tegmark, M., & de Oliveira-Costa, A. 2001, Phys. Rev. D, 64, 063001

Torrado, J., & Lewis, A. 2021, J. Cosmol. Astropart. Phys., 2021, 057

Tristram, M., Macías-Pérez, J. F., Renault, C., & Santos, D. 2005, MNRAS, 358, 833

Tristram, M., Banday, A. J., Górski, K. M., et al. 2021, A&A, 647, A128

Tristram, M., Banday, A. J., Górski, K. M., et al. 2022, Phys. Rev. D, 105, 083524

Tucci, M., Martínez-González, E., Toffolatti, L., González-Nuevo, J., & De Zotti, G. 2004, MNRAS, 349, 1267

Tucci, M., Toffolatti, L., de Zotti, G., & Martínez-González, E. 2011, A&A, 533, 57

Vanneste, S., Henrot-Versillé, S., Louis, T., & Tristram, M. 2018, Phys. Rev. D, 98, 103526

Varshalovich, D. A., Moskalev, A. N., & Khersonskii, V. K. 1988, Quantum Theory of Angular Momentum (Singapore: World Scientific)
## 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.

<table><caption>**Table A.1.** Instrumental and foreground parameters for the HiLLiPoP likelihood with their respective priors. Amplitudes refer to $D_{\ell} = \ell(\ell + 1)C_{\ell}/2\pi$ for $\ell = 3000$ at 143 GHz, except for dust coefficients, $c_{\text{dust}}$, for which the priors are found by rescaling the dust power spectrum at 353 GHz.</caption><thead><tr><th>Name</th><th>Definition</th><th>Prior</th><th>Mean</th></tr></thead><tbody><tr><td>A<sub>planck</sub></td><td>Absolute calibration</td><td>1.0000 ± 0.0025</td><td>0.9997 ± 0.0024</td></tr><tr><td>c<sub>100A</sub></td><td>Map recalibration (100A)</td><td>[0.9,1.1]</td><td>1.003 ± 0.007</td></tr><tr><td>c<sub>100B</sub></td><td>Map recalibration (100B)</td><td>[0.9,1.1]</td><td>1.004 ± 0.007</td></tr><tr><td>c<sub>143A</sub></td><td>Map recalibration (143A)</td><td>1.0 (fixed)</td><td></td></tr><tr><td>c<sub>143B</sub></td><td>Map recalibration (143B)</td><td>[0.9,1.1]</td><td>1.004 ± 0.006</td></tr><tr><td>c<sub>217A</sub></td><td>Map recalibration (217A)</td><td>[0.9,1.1]</td><td>1.001 ± 0.008</td></tr><tr><td>c<sub>217B</sub></td><td>Map recalibration (217B)</td><td>[0.9,1.1]</td><td>1.001 ± 0.008</td></tr><tr><td>η<sub>100-A</sub></td><td>Cross-polarization (100-A)</td><td>1.000 (fixed)</td><td></td></tr><tr><td>η<sub>100-B</sub></td><td>Cross-polarization (100-B)</td><td>1.000 (fixed)</td><td></td></tr><tr><td>η<sub>143-A</sub></td><td>Cross-polarization (143-A)</td><td>1.000 (fixed)</td><td></td></tr><tr><td>η<sub>143-B</sub></td><td>Cross-polarization (143-B)</td><td>1.000 (fixed)</td><td></td></tr><tr><td>η<sub>217-A</sub></td><td>Cross-polarization (217-A)</td><td>0.975 (fixed)</td><td></td></tr><tr><td>η<sub>217-B</sub></td><td>Cross-polarization (217-B)</td><td>0.975 (fixed)</td><td></td></tr><tr><td>c<sup>T</sup><sub>dust</sub></td><td>Rescaling for Galactic dust in temperature</td><td>1.0 ± 0.1</td><td>1.08 ± 0.03</td></tr><tr><td>c<sub>dust</sub></td><td>Rescaling for Galactic dust in polarization</td><td>1.0 ± 0.1</td><td>1.20 ± 0.03</td></tr><tr><td>A<sub>radio</sub></td><td>Amplitude for radio sources</td><td>[0,150]</td><td>63.3 ± 4.7</td></tr><tr><td>A<sub>IR</sub></td><td>Amplitude for IR sources</td><td>[0,150]</td><td>6.07 ± 0.63</td></tr><tr><td>A<sub>CIB</sub></td><td>Amplitude for the CIB</td><td>[0,20]</td><td>1.03 ± 0.34</td></tr><tr><td>A<sub>tSZ</sub></td><td>Amplitude for the tSZ effect</td><td>[0,50]</td><td>5.87 ± 1.66</td></tr><tr><td>A<sub>kSZ</sub></td><td>Amplitude for the kSZ effect</td><td>[0,50]</td><td>&lt; 7.6</td></tr><tr><td>ξ<sub>SZ×CIB</sub></td><td>Cross-correlation tSZ×CIB</td><td>[-1,1]</td><td>0.46 ± 0.30</td></tr><tr><td>β<sup>T</sup><sub>dust</sub></td><td>Spectral index for dust in temperature</td><td>1.51 ± 0.01</td><td>1.51 ± 0.01</td></tr><tr><td>β<sub>dust</sub></td><td>Spectral index for dust in polarization</td><td>1.59 ± 0.02</td><td>1.59 ± 0.02</td></tr><tr><td>β<sub>CIB</sub></td><td>Spectral index for CIB</td><td>1.75 ± 0.06</td><td>1.85 ± 0.06</td></tr><tr><td>β<sub>radio</sub></td><td>Spectral index for radio sources</td><td>-0.8</td><td></td></tr></tbody></table>
## 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 $\chi^2$ values are given in Table B.1.

<figure><img src="image_26.png" /><figcaption>**Fig. B.1.** Best-fit model for each cross-frequency power spectrum in temperature, including emission from CMB, dust, tSZ, kSZ, CIB, SZxCIB, 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 $\Delta\ell = 20$ for this plot.</figcaption></figure>

<table><caption>**Table B.1.** $\chi^2$ values for each cross-spectrum compared to the size of the data vector ($n_d$).</caption><thead><tr><th rowspan="2">Cross-spectrum</th><th colspan="2">TT</th><th colspan="2">EE</th><th colspan="2">TE</th><th colspan="2">ET</th></tr><tr><th>χ<sup>2</sup>/n<sub>d</sub></th><th>δσ(χ<sup>2</sup>)</th><th>χ<sup>2</sup>/n<sub>d</sub></th><th>δσ(χ<sup>2</sup>)</th><th>χ<sup>2</sup>/n<sub>d</sub></th><th>δσ(χ<sup>2</sup>)</th><th>χ<sup>2</sup>/n<sub>d</sub></th><th>δσ(χ<sup>2</sup>)</th></tr></thead><tbody><tr><td>100A×100B</td><td>1590.0 / 1471</td><td>2.19</td><td>1079.1 / 1101</td><td>-0.47</td><td>1597.4 / 1471</td><td>2.33</td><td>1450.1 / 1471</td><td>-0.39</td></tr><tr><td>100A×143A</td><td>1616.5 / 1471</td><td>2.68</td><td>1551.5 / 1471</td><td>1.48</td><td>1564.5 / 1471</td><td>1.72</td><td>1490.8 / 1471</td><td>0.37</td></tr><tr><td>100A×143B</td><td>1605.1 / 1471</td><td>2.47</td><td>1431.3 / 1471</td><td>-0.73</td><td>1396.4 / 1471</td><td>-1.38</td><td>1520.6 / 1471</td><td>0.92</td></tr><tr><td>100B×143A</td><td>1596.3 / 1471</td><td>2.31</td><td>1485.7 / 1471</td><td>0.27</td><td>1535.2 / 1471</td><td>1.18</td><td>1615.8 / 1471</td><td>2.67</td></tr><tr><td>100B×143B</td><td>1576.5 / 1471</td><td>1.94</td><td>1495.5 / 1471</td><td>0.45</td><td>1466.9 / 1471</td><td>-0.08</td><td>1614.6 / 1471</td><td>2.65</td></tr><tr><td>100A×217A</td><td>1379.1 / 1251</td><td>2.56</td><td>1331.5 / 1251</td><td>1.61</td><td>1478.0 / 1401</td><td>1.45</td><td>1432.3 / 1401</td><td>0.59</td></tr><tr><td>100A×217B</td><td>1364.5 / 1251</td><td>2.27</td><td>1278.4 / 1251</td><td>0.55</td><td>1481.3 / 1401</td><td>1.52</td><td>1445.3 / 1401</td><td>0.84</td></tr><tr><td>100B×217A</td><td>1336.8 / 1251</td><td>1.71</td><td>1283.0 / 1251</td><td>0.64</td><td>1507.3 / 1401</td><td>2.01</td><td>1545.9 / 1401</td><td>2.74</td></tr><tr><td>100B×217B</td><td>1335.0 / 1251</td><td>1.68</td><td>1218.3 / 1251</td><td>-0.65</td><td>1466.8 / 1401</td><td>1.24</td><td>1505.7 / 1401</td><td>1.98</td></tr><tr><td>143A×143B</td><td>2108.5 / 1951</td><td>2.52</td><td>1995.3 / 1971</td><td>0.39</td><td>2014.7 / 1971</td><td>0.70</td><td>1972.1 / 1971</td><td>0.02</td></tr><tr><td>143A×217A</td><td>2324.1 / 2251</td><td>1.09</td><td>1647.4 / 1751</td><td>-1.75</td><td>1847.9 / 1801</td><td>0.78</td><td>1868.7 / 1801</td><td>1.13</td></tr><tr><td>143A×217B</td><td>2327.3 / 2251</td><td>1.14</td><td>1853.6 / 1751</td><td>1.73</td><td>1746.9 / 1801</td><td>-0.90</td><td>1898.1 / 1801</td><td>1.62</td></tr><tr><td>143B×217A</td><td>2351.2 / 2251</td><td>1.49</td><td>1725.3 / 1751</td><td>-0.43</td><td>1812.0 / 1801</td><td>0.18</td><td>1835.4 / 1801</td><td>0.57</td></tr><tr><td>143B×217B</td><td>2321.0 / 2251</td><td>1.04</td><td>1799.5 / 1751</td><td>0.82</td><td>1696.4 / 1801</td><td>-1.74</td><td>1862.6 / 1801</td><td>1.03</td></tr><tr><td>217A×217B</td><td>2283.6 / 2251</td><td>0.49</td><td>1732.8 / 1751</td><td>-0.31</td><td>1625.4 / 1701</td><td>-1.30</td><td>1725.8 / 1701</td><td>0.43</td></tr></tbody></table>