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

We present constraints on cosmological parameters using maps from the last Planck data release (PR4). In particular, we detail an upgraded version of the cosmic microwave background likelihood, HiLLiPoP , based on angular power spectra and relying on a physical modelling 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-ℓ measurements from LoLLiPoP , we derive constraints on Λ CDM parameters that are in good agreement with previous Planck 2018 results, but with 10% to 20% smaller uncertainties. We demonstrate that the foregrounds can be accurately described in the spectral domain with only negligible impact on Λ CDM parameters. We also derive constraints on single-parameter extensions to Λ CDM including A L , Ω K , N e ff , and (cid:80) m ν . Noteworthy results from this updated analysis include a lensing amplitude value of A L = 1 . 039 ± 0 . 052, which aligns more closely with theoretical expectations within the Λ CDM framework. Additionally, our curvature measurement, Ω K = − 0 . 012 ± 0 . 010, now demonstrates complete consistency 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 σ level). 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.


Introduction
Since the first results were released in 2013, the Planck satellite's measurements of the cosmic microwave background (CMB) anisotropies have provided highly precise constraints on cosmological models.These measurements have tested the cosmological-constant-dominated cold dark matter (ΛCDM) model, given tight constraints on its parameters, and ruled out many plausible extensions.As a consequence, the best-fitting 6parameter ΛCDM model is now frequently used as the standard reference to be compared to new observational results, and when combining 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).As well as including 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 better sensitivity on almost all cosmological parameters owing to improved map sensitivity.Additionally, we look for better internal consistency for the lensing amplitude affecting the primordial CMB anisotropies.We thus derive constraints on cosmology using both low-ℓ and high-ℓ likelihoods based on Planck PR4.The only part still relying on PR3 (also known as Planck 2018) is the low-ℓ temperature likelihood, Commander, as we do not anticipate significant improvements at large scales in temperature between PR3 and PR4.On the other hand, our analysis includes the large scales in polarization from PR4 for which the NPIPE processing provides a significant improvement compared to PR3.
Since the foregrounds dominate polarization at large scales, for the low-ℓ likelihood, LoLLiPoP, we make use of componentseparated 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) and Tristram et al. (2022), where it was combined with the BICEP2/Keck likelihood (Ade et al. 2021) in order to provide constraints on the tensor-to-scalar ratio r.
For the high-ℓ power-spectrum analysis, HiLLiPoP, we use a multi-frequency Gaussian likelihood approximation using sky maps at three frequencies (100, 143 and 217 GHz), while the channel at 353 GHz is used to derive a template for the dust power spectrum contaminating the CMB signal at large scales.HiLLiPoP is one of the likelihoods developed within the Planck collaboration and used to analyse previous Planck data sets (Planck Collaboration XV 2014;Planck Collaboration XI 2016).Here, we describe a new version adapted to PR4 and called "HiLLiPoP V4.2."It differs from the previous one essentially by using a larger sky fraction (covering 75 % of the sky) and a refined model for the foregrounds (in particular for point sources and dust emission).We specifically use High-Frequency Instrument "detsets", which are splits of the detectors at each frequency into specific subsets.We compute cross-spectra for each of the CMB modes (T T , T E, EE), 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 T T , while the impact of the Planck noise in polarization is more visible in T E and EE.However, at those scales (ℓ < 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 top of the CMB signal.As opposed to other Planck likelihoods, HiLLiPoP considers all cross-frequency power spectra.Even if the Planck PR4 data set is dominated by CMB anisotropies over the entire range of multipoles considered in the high-ℓ likelihood (30 < ℓ < 2500), using all crossfrequency spectra allows us to check the robustness of the results with respect to our knowledge of the astrophysical foregrounds.Indeed, even if the basic ΛCDM parameters are insignificantly affected by the details of the foreground modelling, the constraints on extensions to ΛCDM might depend more critically on the accuracy of the foreground description.Moreover, future ground-based experiments, measuring smaller scales than those accessible by Planck, will be even more sensitive to extragalactic foregrounds.
We begin this paper by summarizing the Planck PR4 pipeline (NPIPE), focusing on the improvements as compared to PR3 (Sect.2).Then, in Sect.3, we explain how the angular power spectra are calculated, and describe the masks we use, the multipole ranges, the pseudo-C ℓ algorithm, and the covariance matrix.The LoLLiPoP likelihood is briefly described in Sect.4, with reference to Tristram et al. (2021Tristram et al. ( , 2022)).The HiLLiPoP likelihood is described in Sect.5, including details of foreground modelling and instrumental effects.Results on the parameters for the ΛCDM model are described and commented on in Sect.6. Constraints on foreground parameters and instrumental parameters are discussed in Sects.7 and 8, respectively.Section 9 is dedicated to consistency checks with respect to previous Planck results.Finally we explore some extensions to ΛCDM in Sect.11, specifically the lensing consistency parameter A L , the curvature Ω K , the effective number of neutrino species N eff , and the sum of neutrino masses m ν .

The Planck PR4 data set
The Planck sky measurements used in this analysis are the PR4 maps available from the Planck Legacy Archive 1 (PLA) and from the National Energy Research Scientific Computing Center (NERSC). 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 systematics in both frequency and component-separated maps at essentially all angular scales, as well as notably improved internal consistency between the various frequencies.Moreover, PR4 also provides a set of "End-to-End" Monte Carlo simulations processed with NPIPE, which enables the characterization of potential biases and the uncertainties associated with the pipeline.
To compute unbiased estimates of the angular power spectra, we perform cross-correlations of two independent splits of the data.As shown in Planck Collaboration Int.LVII (2020), the most appropriate split for the Planck data is represented by the detset maps, comprising two subsets of maps with nearly independent noise characteristics, made by combining half of the detectors at each frequency.This was obtained by processing each split independently, in contrast to the split maps produced in the previous Planck releases.We note that time-split maps (made from, e.g., "odd − even rings" or "half-mission data") share the same instrumental detectors, and therefore exhibit noise correlations due to identical spectral bandpasses and optical responses.As a consequence, the use of time-split maps gives rise to systematic biases in the cross-power spectra (see section 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 cross-correlate using detset splits only.
Nevertheless, in order to verify the level of noise correlation between detsets, we computed the detset cross-power spectra from the half-ring difference maps, which we show in Fig. 2. The spectra are computed on 75 % of the sky and are fully compatible with zero, ensuring that any correlated noise is much smaller than the uncorrelated noise over the range of multipoles from ℓ = 30 to 2500.As discussed above, this test is not sensitive to correlations at scales smaller than the half-ring period.Indeed, if both halves of a ring are affected by the same systematic effect, it will vanish in the half-ring difference map and thus will not be tested in cross-correlation with another detset.

Planck PR4 angular power spectra
3.1.Large-scale polarized power spectra Foregrounds are stronger in polarization relative to the CMB than in temperature, and cleaning the Planck frequencies using C ℓ templates in the likelihood (as done at small scales) is not accurate enough, especially at large angular scales.In order to clean sky maps of polarized foregrounds, we use 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 is run on each detset map independently, as well as on each realization from the PR4 Monte Carlo simulations.
We then compute unbiased estimates of the angular power spectra by cross-correlating the two detset-cleaned maps.We compute 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 use downgraded N side = 16 maps (Górski et al. 2005) after convolution with a cosine apodizing kernel b ℓ = 1 2 {1 + cos π(ℓ − 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 ℓ = 30 (for more details, see Tristram et al. 2021).We also performed the same estimation on each of the PR4 simulations and derive the ℓ-by-ℓ covariance matrix that is then used to propagate uncertainties in LoLLiPoP, the low-ℓ CMB likelihood described in Sect. 4.

Sky fractions
For small scales (ℓ > 30), we are using detset maps at frequencies 100, 143, and 217 GHz, and we select only a fraction of the sky in order 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 apply a mask to remove the region of strongest Galactic emission, adapted to each frequency.We can keep a larger sky fraction at the lowest frequency (100 GHz) where the emission from the Galactic sources is low.Since Planck uncertainty is dominated by sample variance up to multipole ℓ ≃ 1800 in temperature (and ℓ ≃ 1100 in T E polarization), this allows us to reduce the sampling variance by ensuring a larger sky fraction.However, we remove a larger fraction of the sky for the highest frequency channel (217 GHz) since it is significantly more contaminated by Galactic dust emission.
We build Galactic masks using the Planck 353-GHz map as a tracer of the thermal dust emission in intensity.In practice, we smooth the Planck 353-GHz map to increase the signal-to-noise ratio 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 use the same Galactic mask for polarization as for temperature.
CO mask.We apply a mask for CO line emission.We consider the combination of maps of the two lines in the Planck frequency bands at 115 and 230 GHz.We smooth the Planck reconstructed CO maps to 30 arcmin before applying a threshold at 2 K km s −1 .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 use a common mask for the three CMB frequencies to cover strong sources (both radio and infrared).In contrast to the masks used in Plik or CamSpec, the point-source mask used in our analysis relies on a more refined procedure that preserves Galactic compact structures and ensures the completeness level at each frequency, but with a higher flux cut on sources (approximately 340, 250, and 200 mJy at 100, 143, and 217 GHz, respectively).The consequence is that these masks leave a slightly greater number of unmasked extragalactic sources, but more accurately preserve the power spectra of dust emission (see Sect. 5.2).We apodize 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 mask a limited number of resolved objects in the sky, essentially nearby galaxies including the LMC, SMC, and M31, as well as the Coma cluster.This removes less than 0.4 % of the sky.
We use the same mask for temperature and polarization.Even though masking point sources in polarization is not mandatory (given the Planck noise in EE, and T E); 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 crossspectra 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 χ 2 ) using various combinations of Galactic masks (see Sect. 9).Fig. 4: Sky masks used for HiLLiPoP V4.2 as a combination of a Galactic mask (blue, green, and red for the 100, 143, and 217 GHz channel, respectively), a CO mask, a point-source mask, and a mask removing nearby galaxies.The effective sky fractions remaining at 100, 143 and 217 GHz are 75 %, 66 %, and 52 %, respectively.

PR4 small-scale spectra
We use Xpol (an extension to polarization of Xspect, described in Tristram et al. 2005) to compute the cross-power spectra in temperature and polarization (T T , EE, and T E).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. 3 Using the six maps presented in Sect.2, we derive 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 ãX ℓm = {ã T ℓm , ãE ℓm , ãB ℓm }, we form the pseudo cross-power spectra between map i and map j, where the vector C ℓ includes the four modes We note that the T E and ET crosspower 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 are 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), 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 debias the beam leakage together with the beam transfer function using the beam window functions evaluated with QuickPol (Hivon et al. 2017).We use 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 crossspectra 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 T T , the agreement between the different spectra is better than 20 µK 2 , except (as expected) for the 100×100 and the 217×217 cases, which are affected by residuals from point sources and Galactic emission (for the latter).In EE, only the 217×217 case is affected by Galactic emission residuals at low multipoles, but the spectra are still consistent at the few µK 2 level.For T E and ET , we can see various features at the level of 10 µK 2 (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.

Multipole ranges
The HiLLiPoP likelihood covers the multipoles starting from ℓ min = 30 up to ℓ max = 2500 in temperature and ℓ max = 2000 in polarization.The multipoles below ℓ < 30 are considered in the low-ℓ likelihoods (Commander and LoLLiPoP, see Sect. 4).
Table 1 gives the HiLLiPoP multipole ranges, [ℓ min , ℓ max ], considered for each of the six cross-frequencies in T T , T E, and EE.The multipole ranges used in the likelihood analysis have been chosen to limit the contamination by Galactic dust emission at low ℓ and instrumental noise at high ℓ.In practice, we ignore the lowest multipoles for cross-spectra involving the 217 GHz map, where dust contamination is the highest, and cut out multipoles higher than ℓ = 1500 for cross-spectra involving the 100 GHz channel given its high noise level.

Channels
Table 1: Multipole ranges used in the HiLLiPoP analysis and corresponding number of ℓs available (n ℓ = ℓ max − ℓ min + 1).The total number of ℓs across all spectra is 29 758.

The covariance matrix
We use a semi-analytical estimate of the C ℓ covariance matrix computed using Xpol.The matrix captures the ℓ-by-ℓ correlations between all the power spectra involved in the analysis.The computation relies directly on data for the estimates.It follows that contributions from noise (correlated and uncorrelated), sky emission (from astrophysical and cosmological origin), and the sample variance are implicitly taken into account in this computation without relying on any model or simulations.
The covariance matrix Σ of the cross-power spectra is directly related to the covariance Σ of the pseudo cross-power spectra through the coupling matrices: with (a, b, c, d) ∈ {T, E} for each map.The matrix Σ, which gives the correlations between the pseudo cross-power spectra (ab) and (cd), is an N-by-N matrix (where N = 4ℓ max ) and reads by expanding the four-point Gaussian correlation using Isserlis' formula (or Wick's theorem).We compute Σ for each pseudo cross-spectra block independently, which includes ℓ-by-ℓ correlation and four spectral mode correlations {T T, EE, T E, ET }.
Each two-point correlation of pseudo-a ℓm s can be expressed as the convolution of C ℓ with a kernel that depends on the polarization mode considered: where the kernels W 0 , W + , and W − are defined as linear combinations of products of Y ℓm of spin 0 and ±2, weighted by the spherical transform of the window function in the pixel domain (the apodized mask).As suggested in Efstathiou (2006), by neglecting the gradients of the window function and applying the completeness relation for spherical harmonics (Varshalovich et al. 1988), we can reduce the products of four Ws into kernels similar to the coupling matrix M defined in Eq. ( 2).In the end, the blocks of the Σ matrices are , which are thus directly related to the measured auto-and crosspower spectra (see the appendix in Couchot et al. 2017b).In practice, to avoid any correlation between C ℓ estimates and their covariance, we use a smoothed version of each measured power spectrum (using a Gaussian filter with σ ℓ = 5) to estimate the covariance matrix.
We finally average the cross-power spectra covariance matrix to form the full cross-frequency power-spectra matrices for the three modes {T T, T E, 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 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 (nonwhite, 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.
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 are focusing only on the E-mode component.Systematic effects are considerably reduced in crosscorrelation compared to auto-correlation, and LoLLiPoP is based on cross-power spectra for which the bias is zero when the noise is uncorrelated between maps.It uses the approximation presented in Hamimeche & Lewis (2008), modified as described in Mangilli et al. (2015) to apply to cross-power spectra.The idea is to apply a change of variable C ℓ → X ℓ so that the new variable X ℓ is nearly Gaussian-distributed. Similarly to Hamimeche & Lewis (2008), we define where g(x) = √ 2(x − ln(x) − 1), C ℓ are the measured crosspower spectra, C ℓ are the power spectra of the model to be evaluated, C f ℓ is a fiducial CMB model, and O ℓ are the offsets needed in the case of cross-spectra.In the case of auto-power spectra, the offsets O ℓ are given by the noise bias effectively present in the measured power spectra.For cross-power spectra, the noise bias is zero, and we use effective offsets defined from the C ℓ noise variance: (5) The distribution of the new variable X ℓ can be approximated as Gaussian, with a covariance given by the covariance of the C ℓ s.The likelihood function of the C ℓ given the data C ℓ is then Uncertainties are incorporated into the C ℓ covariance matrix M ℓℓ ′ , which is evaluated after applying the same pipeline (including Commander component-separation and cross-spectrum estimation on each simulation) to the Monte Carlo simulations provided in PR4.While foreground emission and the cleaning procedure are kept fixed in the simulations (so that we cannot include uncertainties arising from an imperfect foreground model), the resulting C ℓ covariance consistently includes CMB sample variance, statistical noise, and systematic residuals, as well as uncertainties from the foreground-cleaning procedure, together with the correlations induced by masking.We further marginalize 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.
LoLLiPoP is publicly available on GitHub. 4In this work, we consider only the information from E modes, and restrict the multipole range from ℓ = 2 to ℓ = 30.
To cover the low multipoles (ℓ < 30) in temperature, we make use of the Commander T T 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).

Small-scale CMB likelihood: HiLLiPoP
This section describes HiLLiPoP (High-ℓ Likelihood on Polarized Power spectra), including the models used for the foreground residuals and the instrumental systematic residuals.HiLLiPoP 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. 5 We make use of the 15 cross-spectra computed from the six detset maps at 100, 143, and 217 GHz (see Sect. 3).From those 15 crossspectra (one each for 100×100, 143×143, and 217×217; four each for 100×143, 100×217, and 143×217), we derive six crossfrequency spectra after recalibration and co-addition, and these are compared 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 use the full semi-analytic covariance matrix that includes ℓ-by-ℓ correlation, and {T T, T E, EE} mode correlation as described in Sect.3.2.4.

The likelihood approximation
On the full-sky, the distribution of auto-spectra is a scaled-χ 2 with 2ℓ + 1 degrees of freedom.The distribution of the crossspectra is slightly different (see appendix A in Mangilli et al. 2015); however, above ℓ = 30 the number of modes is large enough that we can safely assume that the C ℓ are Gaussiandistributed.Consequently, for high multipoles the resulting likelihood can be approximated by a multivariate Gaussian, including correlations between the values of C ℓ arising from the cutsky, and reads where ℓ denotes the residual of the estimated crosspower spectrum C ℓ with respect to the model C ℓ , which depends on the frequencies {i, j} and is described in the next section.The matrix Σ = RR T is the full covariance matrix that includes the instrumental variance from the data as well as the cosmic variance from the model.The latter is directly proportional to the model so that the matrix Σ should, in principle, depend on the model.In practice, given our current knowledge of the cosmological parameters, the theoretical power spectra typically differ from each other at each ℓ by less than they differ from the observed C ℓ , so that we can expand Σ around a reasonable fiducial model.As described in Planck Collaboration XV (2014), the additional terms in the expansion are small if the fiducial model is accurate and leaving it out entirely does not bias the likelihood.Using a fixed covariance matrix Σ, we can drop the constant term ln |Σ| and recover nearly optimal variance (see Carron 2013).Within the approximations discussed above, we expect the likelihood to be χ 2 -distributed with a mean equal to the number of degrees of freedom n dof = n ℓ − n p (n ℓ being the number of band powers in the power spectra and n p the number of fitted parameters) and a variance equal to 2n dof .

The model
We now present the model ( Ĉℓ ) 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 of 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 We include in the foregrounds, for the temperature likelihood, contributions from: -Galactic dust; cosmic infrared background (CIB); thermal (tSZ) and kinetic (kSZ) Sunyaev-Zeldovich components; -Poisson-distributed point sources from radio and infrared star-forming galaxies; 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 components, we take into account the bandpass response using effective frequencies as listed in table 4 of Planck Collaboration IX (2014).Galactic emission from freefree 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 were not able to detect any residuals from Galactic synchrotron or free-free emission.In the following, we therefore neglect these contributions.
Galactic dust emission.At frequencies above 100 GHz, Galactic emission is dominated by dust.The dust template is fitted on the Planck 353-GHz data using a power-law model.In practice, we compute the 353-GHz cross-spectra Ĉ353Ax353B ℓ for each pair of masks (M i , M j ) associated with the cross-spectra ν i × ν j (Fig. 7).We then subtract the Planck best-fit CMB power spectrum and fit a power-law model with a free constant Aℓ α d + B in the range ℓ = [30, 1500] for T T , to account for the unresolved point sources at 353 GHz.A simple power law is used to fit the EE and T E power spectra in the range ℓ = [30,1000].Thanks to the use of the point-source mask (described in Sect.3.2.1),our Galactic dust residual power spectrum is much simpler than in the case of other Planck likelihoods.Indeed, the point-source masks used in the Planck PR3 analysis removes some Galactic structures and bright cirrus, which induces an artificial knee in the residual dust power spectra around ℓ = 200 (see section 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 T E power spectra are directly comparable to those derived in Planck Collaboration Int.XXX (2016), with indices of α d = −2.3 and −2.4 for EE and T E, respectively, the indices for T T vary with the sky fraction considered, ranging from α d = −2.2 down to −2.6 for the largest sky fraction.Fig. 7: Dust power spectra, D ℓ = ℓ(ℓ + 1)C ℓ /2π, at 353 GHz for T T (top), EE (middle), and T E (bottom).The power spectra are computed from cross-correlation between detset maps at 353 GHz for different sets of masks, as defined in Sect.3.2.1, and further corrected for the CMB power spectrum (solid black line) and CIB power spectrum (dashed black line).The coloured dashed lines are simple fits, as described in the text.
For each polarization mode (T T , EE, T E), we then extrapolate the dust templates at 353 GHz for each cross-mask to the cross-frequency considered: where 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 powerlaw at high multipoles.The template is rescaled by A CIB , the amplitude of the contamination at our reference frequency (ν 0 = 143 GHz) and ℓ = 3000.The emission law is modelled by a modified blackbody a CIB ν = ν β CIB B ν (T ) with a fixed temperature (T = 25 K) and a variable index β CIB .We use a strong prior β 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), Thermal Sunyaev-Zeldovich (tSZ) effect.The template for the tSZ emission comes from the halo model fitted on Planck measurements in Planck Collaboration XXII (2016) and used more recently with PR4 data in Tanimura et al. (2022).The tSZ signal is parameterized by a single amplitude A tSZ , corresponding to the amplitude of the tSZ signal at our reference frequency (ν 0 = 143 GHz) at ℓ = 3000, where Kinetic Sunyaev-Zeldovich (kSZ) effect.The kSZ emission is parameterized by A kSZ , the amplitude at ℓ = 3000, scaling a fixed template including homogeneous and patchy reionization components from Shaw et al. (2012) and Battaglia et al. (2013), Thermal SZ×CIB correlation.The cross-correlation between the thermal SZ and the CIB is parameterized as with ξ the correlation coefficient rescaling the template D tSZ×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 rad ν ∝ ν −β s (Tucci et al. 2011) while we use a IR ν = ν β IR B ν (T ) (Béthermin et al. 2012) for infrared dusty star-forming galaxies.The residual cross-power spectra for point sources are finally Following Lagache et al. (2020), radio source emission is dominated at frequencies above about 100 GHz by radio quasars whose spectral indices can vary from −1.0 to 0.0 (Planck Collaboration XIII 2011; Planck Collaboration Int.VII 2013).We constrain the SED by fixing β s = − 0.8, following results from Reichardt et al. (2021).For infrared dusty starforming galaxies, we adopt β IR identical to β CIB and T = 25 K.The C ℓ s are then converted into D ℓ s such that the amplitudes A rad and A IR refer to the amplitude of D 3000 at 143 GHz.In 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 choose to set priors on the SED parameters so that the correlation between the amplitudes of residuals is 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.

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 sample five inter-calibration coefficients while fixing as the reference the calibration of the most sensitive map (the first detset at 143 GHz, 143A).In addition, we sample a Planck calibration parameter A Planck with a strong prior, A Planck = N(1.0000,0.0025), in order 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).
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.
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, 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.

Results on the 6-parameter ΛCDM model
In this section, we describe the constraints on cosmological parameters in the ΛCDM model using the Planck PR4 data.In addition to HiLLiPoP (hlp), we also make use of the Commander low-ℓ likelihood (lowT, see Planck Collaboration IV 2020) and the polarized low-ℓ EE likelihood LoLLiPoP (lolE, discussed in Sect.4).We define the following combination of likelihoods for the rest of the paper: -TT, lowT+hlpTT; -TE, lowT+lolE+hlpTE; -EE, lolE+hlpEE; -TTTEEE, lowT+lolE+hlpTTTEEE.Note that for "TT", we only use data from temperature and combine lowT+hlpTT; this is in contrast to Planck Collaboration VI (2020) and Rosenberg et al. (2022), in which low-ℓ data from EE are systematically added in order to constrain the reionization optical depth.
The model for the CMB is computed by numerically solving the background and perturbation equations for a specific cosmological model using CAMB (Lewis et al. 2000;Howlett et al. 2012). 6In this paper, we consider a ΛCDM model with six free parameters describing: the current physical densities of baryons (Ω b h 2 ) and cold dark matter (Ω c h 2 ); the angular acoustic scale (θ * ); the reionization optical depth (τ); and the amplitude and spectral index of the primordial scalar spectrum (A s and n s ).Here h is the dimensionless Hubble constant, h = H 0 /(100 km s −1 Mpc −1 ).
In addition, we fit for six inter-calibration parameters, seven foreground residual amplitudes in temperature (c T dust , A radio , A IR , A CIB , A tSZ , A kSZ , and ξ SZ×CIB ), plus one in polarization (c P dust ), and three foreground spectral indices (β T dust , β P dust , and β CIB ).Foreground and instrumental parameters are listed in Table A.1, together with their respective priors.
To quantify the agreement between the data and the model, we computed the χ 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 χ 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 χ 2 value for hlpTT (corresponding to about 2.7 σ), while the hlpTE and hlpEE χ 2 values are compatible with unity, at 1.8 σ and 0.1 σ, respectively.For the full combination hlpTTTEEE, we obtain a χ 2 = 30495 for a data size of 29768, corresponding to a 3.02 σ deviation.As described in Rosenberg et al. (2022), where the goodness of fit is also somewhat poor (4.07 σ for TT and 4.46 σ for the TTTEEE), this could be explained by a slight misestimation of the instrumental noise, rather than a bias that could be fit by an improved foreground model or a different cosmology.However, we emphasize that the level of this divergence is small, since the recovered reduced-χ 2 , χ 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 Co-added CMB power spectra are shown in Figs. 9 and 10, for T T , T E, and EE; they are compared to the best-fit obtained with the full TTTEEE combination.Planck spectra are binned with ∆ℓ = 30 for the plots, but considered ℓ-by-ℓ in the likelihood.The plots also show the residuals relative to the ΛCDM best-fit to TTTEEE, as well as the normalized residuals.We cannot identify any deviation from statistical noise or any bias from foreground residuals.
In Fig. 11, we compare the constraints on ΛCDM parameters obtained using TT, TE, and EE and their combination.We find very good consistency between TT and TE, while EE constraints are wider, with a deviation in the acoustic scale θ * toward lower values.This feature of the Planck PR4 data was previously reported in Rosenberg et al. (2022), in which the authors studied the correlation with other parameters and concluded that this is likely due to parameter degeneracies coupling to residual systematics in EE.However, the deviation of θ * between EE and TT is now reduced with the increase of the sky fraction enabled by HiLLiPoP V4.2, though still present at the 1.6 σ level.In addition, we have checked that this shift in θ * is not related to any super-sample lensing effect (as described in Manzotti et al. 2014), or to any aberration correction (see Jeong et al. 2014), both of which are negligible for the large sky fraction considered in the Planck data set.We note that, interestingly, θ * is the only parameter that deviates in EE; the others, including H 0 , are compatible with TT at much better than 1 σ.Given the weak sensitivity of the Planck EE spectra as compared to T T and T E, discrepancies in the EE parameter reconstruction will have little impact on overall cosmological parameter results.
The HiLLiPoP V4.2 constraints on ΛCDM cosmological parameters are summarized in Table  For the constraint on the Hubble constant, we obtain 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 compatible with PR3 results (σ 8 = 0.8120 ± 0.0073) but lower by 0.5 σ.The matter density, Ω m also shifts by roughly 1 σ, so that S 8 ≡ σ 8 (Ω m /0.3) 0.5 = 0.819 ± 0.014.Before discussing results on the foreground parameters (Sect.7) and instrumental parameters (Sect.8), we show in Fig. 8 the correlation matrix for the fitted parameters.We can see that foreground parameters are only weakly correlated with the cosmological parameters and the inter-calibrations.This strengthens the robustness of the results with respect to the foreground model and ensures very low impact on cosmology.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 ℓ.We report mean values and symmetrical 68 % confidence intervals.

Foreground parameters
All Planck cross-spectra are dominated by the CMB signal at all the scales we consider.This is illustrated for T T in Fig. B.1 of Appendix B, where we show each component of the model fitted in the likelihood with the best-fit parameters for the six crossfrequencies.It is also true for T E and EE.Thanks to the multifrequency 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 with a population spectral index for the radio power fixed to β s = −0.8,close to the value recovered by the SPT team (β s = −0.76±0.15, Reichardt et al. 2021).Allowing β s to vary in Planck data, gives β s = −0.54± 0.08, with a corresponding increase of the amplitude A radio .This also impacts the SZ-CIB cross-correlation amplitude with a significant increase of ξ.
We obtain a high-significance detection of CIB anisotropies, with amplitudes at 143 GHz and ℓ = 3000 given by for the clustered and Poisson parts, respectively.We note that these amplitudes cannot be directly compared to values in previous works because they strongly depend on the prior used for the β CIB index for the former and on the flux cut applied by the point-source mask for the latter.The thermal Sunyaev-Zeldovich effect is also significantly detected, with an amplitude at 143 GHz and ℓ = 3000 of A tSZ = (5.9± 1.7) µK 2 .
(21) This is close to (but somewhat higher than) what is reported in Reichardt et al. (2021), with A tSZ = (3.42± 0.54) µK 2 , even though the uncertainties are larger.However, it is more closely comparable with ACTpol results, A tSZ = (5.29 ± 0.66) µ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: We note that those last results are about 10 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 find rescaling coefficients c dust to be 1.08 ± 0.03 and 1.20 ± 0.03 for 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, β P dust = N(1.59,0.02), we recover an amplitude higher than expected.While changing the models as described above, the impact on ΛCDM parameters is 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 modelling, with variations limited to less than 2 % (4 % for n s ).

Consistency between Planck likelihoods
We now 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 obtain similar χ 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 ∆χ 2 values are lower than 100 for 29758 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 modelling of foreground residuals in cross-frequency spectra, we find consistent constraints at better than the 1 σ level.This gives confidence in the robustness of our cosmological constraints.
Figure 15 shows the 1-d posterior distributions for the ΛCDM parameters using different sky fractions.We also compare to the posteriors obtained from Planck PR3 and with CamSpec PR4 (where we use LoLLiPoP instead of the polarized low-ℓ constraint from PR3 used in Rosenberg et al. 2022).We find good consistency between the different likelihoods and between the two data sets (PR3 and PR4).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 Ω c h 2 , for which HiLLiPoP on PR4 finds a value 1.0 σ lower.Associated with Commander and LoLLiPoP, CamSpec on PR4 also gives lower Ω c h 2 by −0.45 σ.The spectral index n s is found to be a bit higher with HiLLiPoP by 0.7 σ.
As discussed in Sect.6, we obtain a slightly higher value for the Hubble constant (+0.6 σ) with h = 0.6766 ± 0.0053, compared to h = 0.6727 ± 0.0060 for PR3.The amplitude of density fluctuations, σ 8 , and the matter density, Ω 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 (θ * ).

Combination with other data sets
We now present some results of our new likelihood in combination with CMB lensing measurements using the Planck PR4 data (Carron et al. 2022).We specifically use the conservative range recommended in Carron et al. (2022), consisting of nine power bins between multipoles of 8 and 400.The addition of the C ϕϕ ℓ 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.
Table 5 presents the constraints on the 6-parameter ΛCDM model when adding lensing and BAO data.Figure 16 shows the posterior distribution for the particular subset Ω b h 2 , Ω m , σ 8 , and H 0 .

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

Gravitational lensing, A L
We sample the phenomenological extension A L in order to check the consistency of the Planck PR4 data set with the smoothing of the power spectra by weak gravitational lensing as predicted by the ΛCDM model.A mild preference for A L > 1 was seen in the Planck PR1 data (Planck Collaboration XVI 2014) and  In Rosenberg et al. (2022), the CamSpec likelihood associated with low-ℓ likelihoods from Planck 2018 also showed a decrease in the A L parameter in Planck PR4 data compared to PR3 data, reducing the difference from unity from 2.4 σ to 1.7 σ.When LoLLiPoP is adopted as the low-ℓ polarized likelihood, instead of the low-ℓ likelihoods from Planck 2018, the constraint on A L from CamSpec changed from A L = 1.095 ± 0.056 to A L = 1.075 ± 0.058, still a 1.3 σ difference from unity.We compare the posteriors for Plik (PR3), CamSpec (PR4), and HiLLiPoP (PR4) in Fig. 18.Previously, when there was a preference for A L > 1, adding A L as a seventh parameter could lead to shifts in other cosmological parameters (e.g., Planck Collaboration Int.LI 2017).However, we confirm that with HiLLiPoP on PR4, the ΛCDM parameters are only affected through a very slight increase of the error bars, without significantly affecting the mean posterior values.
With the PR4 lensing reconstruction described in Carron et al. (2022), the amplitude of the lensing power spectrum is 1.004 ± 0.024 relative to the Planck 2018 best-fit model.When combining CMB lensing with TTTEEE we then recover a tighter constraint on A L , with A L = 1.037 ± 0.037 (TTTEEE+lensing). (36) 11.2.Curvature, Ω K For the spatial curvature parameter, we report a significant difference with respect to Planck Collaboration VI (2020), which used PR3 and reported a mild preference for closed models (i.e., Ω K < 0).Indeed, with HiLLiPoP V4.2, the measurements are consistent with a flat universe (Ω K = 0) for all spectra.
As noticed in Rosenberg et al. (2022), with Planck PR4, the constraint on Ω K is more precise and shifts toward zero, along the so-called geometrical degeneracy with H 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 which is only 1.2 σ discrepant from zero.This is to be compared to Ω K = −0.044+0.018 −0.015 obtained for Plik on PR3 (Planck Collaboration VI 2020) and Ω 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 -Ω K plane at low H 0 and negative Ω 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±3.60)km s −1 Mpc −1 , only 1.3 σ away from the constraint with fixed Ω K = 0.The combination of TTTEEE with lensing yields the improved constraint This is now compatible with the baryon acoustic oscillation measurements from SDSS, which are consistent with a flat Universe and give Ω K = −0.0022± 0.0022 (Alam et al. 2021).Finally, the mean posterior for the combination of Planck PR4 TTTEEE with lensing and BAO is This is consistent with our Universe being spatially flat to within a 1 σ accuracy of 0.16 % (see Fig. 20). Figure 21 shows the posteriors for TT, TE, EE, and their combination when we consider the N eff extension.Both TT and TE are compatible with similar uncertainties, while EE is not sensitive to N eff .The mean posterior for TTTEEE is The uncertainties are comparable to Planck 2018 results (N eff = 2.92 ± 0.19, Planck Collaboration VI 2020) with a slight shift toward higher values, closer to the theoretical expectation N eff = 3.044 (Akita & Yamaguchi 2020;Froustey et al. 2020;Bennett et al. 2021), which was also reported with CamSpec analysis based on PR4 data (N eff = 3.00 ± 0.21, Rosenberg et al. 2022).11.4.Sum of the neutrino masses, m ν Figure 22 shows the posterior distribution for the sum of the neutrino masses, m ν .There is no detection of the effects of neutrino mass and we report an upper limit of m ν < 0.39 eV (95 % CL, TTTEEE).( 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: m ν < 0.26 eV at 95 % CL.Our constraint is comparable to CamSpec, which gives m ν < 0.36 eV at 95 % CL.As explained in Couchot et al. (2017a) and Planck Collaboration VI (2020), this is directly related to the value of A L .Indeed, the correlation between A L and m ν pushes the peak posterior of m ν toward negative values when A L is fixed to unity; the data, however, prefer values of A L larger than 1.With HiLLiPoP V4.2, the value of A L reported in this work is more compatible with unity (A L = 1.039 ± 0.052, see Sect.11.1), thus, the posterior for m ν is shifted to higher values, with a peak closer to zero, increasing the upper limit accordingly.Figure 23 shows constraints in the m ν -τ plane when combining our new likelihood with with CMB lensing and BAO data.This combination further strengthens the limits to m ν < 0.26 eV (95 % CL, TTTEEE+lensing), (42) m ν < 0.11 eV (95 % CL, TTTEEE+lensing+BAO). ( 43) This is slightly tighter than the upper limit from Planck 2018 ( m ν < 0.12 eV) and getting close to the lower-limit for the inverted mass hierarchy ( m ν > ∼ 0.1 eV, see e.g., Jimenez et al. 2022).

Conclusions
In this paper, we have derived cosmological constraints using CMB anisotropies from the final Planck data release (PR4).We detailed a new version of a CMB high-ℓ likelihood based on cross-power spectra computed from the PR4 maps.This version of HiLLiPoP, labelled V4.2, uses more sky (75 %) and a wider range of multipoles.Our likelihood makes use of physicallymotivated 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 modelling on cosmological ΛCDM parameters is negligible.
Combined with the low-ℓ EE likelihood LoLLiPoP, we derived constraints on ΛCDM and find good consistency with Planck 2018 results (based on PR3) with better goodness-of-fit and higher sensitivity (from 10 % to 20 %, depending on the parameters).In particular, we now constrain the reionization optical depth at the 10 % level.We found a value for the Hubble constant consistent with previous CMB measurements and thus still in tension with distance-ladder results.We also obtained a lower value for S 8 , alleviating the CMB versus large-scale structure tension to 1.5 σ.
We found good consistency with the other published CMB likelihood analysis based on PR4, CamSpec (Rosenberg et al. 2022), which relies on a procedure to clean power spectra prior to constructing the likelihood.The consistency of the results using two different approaches reinforces the robustness of the results obtained with Planck data.
We also add constraints from PR4 lensing, making the combination the most constraining data set that is currently available from Planck.Additionally we explore adding baryon acoustic oscillation data, which tightens limits on some particular extensions to the standard cosmology.
We provided constraints on some extensions to ΛCDM, including the lensing amplitude A L , the curvature Ω K , the effective number of relativistic species N eff , and the sum of the neutrino masses m ν .For both A L and Ω K , our results show a significant reduction of the so-called "tensions" with standard ΛCDM, together with a reduction of the uncertainties.The final constraints indeed are fully compatible with ΛCDM predictions.In particular, with the new version of the likelihood presented in this work, we report A L = 1.039 ± 0.052, entirely compatible with the ΛCDM prediction.The better agreement is explained both by the improvement of the Planck maps thanks to the NPIPE processing (with less noise and better systematic control in polarization) and the use of the LoLLiPoP and HiLLiPoP likelihoods.

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

Fig. 2 :
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 T T , EE, T E, and ET .Spectra are binned with ∆ℓ = 40.The projections on the right show the distribution for each unbinned spectrum over the range ℓ = 30-2500.

Fig. 3 :
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.

Fig. 5 :
Fig. 5: Frequency cross-power spectra with respect to the mean spectra for T T , EE, T E, and ET .Spectra are binned with ∆ℓ = 40 for this figure.

Fig. 6 :
Fig.6: Full HiLLiPoP covariance matrix, including all correlations in multipoles between cross-frequencies and power spectra.
is a modified blackbody with T d fixed to 19.6 K, while c dust and β d are sampled independently for temperature and polarization.We use Gaussian priors for the spectral indices β d from Planck Collaboration Int.XXII (2015), which gives β T d = N(1.51,0.01) and β T d = N(1.59,0.02) for temperature and polarization, respectively.The coefficient c dust allows us to propagate the uncertainty from fitting the 353-GHz dust spectrum with a power law.We sample c dust with a Gaussian prior, c dust = N(1.0,0.1). .
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 to 20 %, depending on the parameter.The reionization optical depth is now constrained at close to the 10 % level: τ = 0.058 ± 0.006.(14) This is the result of the NPIPE treatment of the PR4 data associated with the low-ℓ likelihood LoLLiPoP (see Planck Collaboration Int.LVII 2020).
(17)    Compared to PR3 (S 8 = 0.834±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 ± 0.019, for ΛCDM with fixed m ν ,Abbott et al. 2022), decreasing the CMB versus largescale structure tension on S 8 from 2.1 σ to 1.5 σ.

Fig. 8 :Fig. 9 :Fig. 10 :
Fig. 8: 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.

Fig. 13 :Fig. 14 :
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.

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

Table 2 :
.1.χ 2 values compared to the size of the data vector (n d ) for each of the Planck HiLLiPoP likelihoods.Here δσ(χ

Table 4 :
Relative variation and improvement in the error bars between Planck 2018 and this work for each cosmological parameter.

Table 5 :
(Couchot et al. 2017c 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 ℓ, with the addition of CMB lensing and BAO constraints.We report mean values and symmetrical 68 % confidence intervals.sincetheanalysis 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).Posterior distributions for A L .With Planck PR4, we find results even more compatible with unity compared to previous releases.Indeed for TTTEEE, we now obtainA L = 1.039 ± 0.052,(35)which is compatible with the ΛCDM expectation (at the 0.7 σ level).As shown in Table6, while the results for EE and TE are compatible with unity, the A L value for TT is still high by 0.8 σ.Figure17shows posterior distributions of A L for each of the mode-spectra and for the TTTEEE combination using Planck PR4.

Table 6 :
Mean values and 68 % confidence intervals for A L .The significance of the deviation from unity is given in the last column.

Table B .
1: χ 2 values for each cross-spectrum compared to the size of the data vector (n d ).