Planck constraints on the tensor-to-scalar ratio

We present constraints on the tensor-to-scalar ratio r using Planck data. We use the latest release of Planck maps (PR4), processed with the NPIPE code, which produces calibrated frequency maps in temperature and polarization for all Planck channels from 30 GHz to 857 GHz using the same pipeline. We compute constraints on r using the BB angular power spectrum, and also discuss constraints coming from the TT spectrum. Given Planck's noise level, the TT spectrum gives constraints on r that are cosmic-variance limited (with $\sigma$(r)=0.093), but we show that the marginalized posterior peaks toward negative values of r at about the 1.2$\sigma$ level. We derive Planck constraints using the BB power spectrum at both large angular scales (the"reionization bump") and intermediate angular scales (the"recombination bump") from $\ell$=2 to 150, and find a stronger constraint than that from TT, with $\sigma$(r)=0.069. The Planck BB spectrum shows no systematic bias, and is compatible with zero, given both the statistical noise and the systematic uncertainties. The likelihood analysis using B modes yields the constraint r<0.158 at 95% confidence using more than 50% of the sky. This upper limit tightens to r<0.069 when Planck EE, BB, and EB power spectra are combined consistently, and tightens further to r<0.056 when the Planck TT power spectrum is included in the combination. Finally, combining Planck with BICEP2/Keck 2015 data yields an upper limit of r<0.044.


Introduction
Gravitational waves entering the horizon between the epoch of recombination and the present day generate a tensor contribution to the large-scale cosmic microwave background (CMB) anisotropy. Hence, primordial tensor fluctuations contribute to the CMB anisotropies, both in temperature (T ) and in polarization (E and B modes).
As described in Planck Collaboration VI (2020) and Planck Collaboration X (2020), the comoving wavenumbers of tensor modes probed by the CMB temperature anisotropy power spectrum have k < ∼ 0.008 Mpc −1 , with very little sensitivity to higher wavenumbers because gravitational waves decay on subhorizon scales. The corresponding multipoles in the harmonic domain are < ∼ 100, for which, in temperature, the scalar perturbations dominate with respect to tensor modes. The tensor component can be fitted together with the scalar one, and the precision of the Planck constraint is limited by the cosmic variance of the large-scale anisotropies.
In polarization, the EE and T E spectra also contain a tensor signal coming from the last-scattering and reionization epochs.
However, the addition of Planck polarization constraints at > ∼ 30 does not significantly change the results coming from temperature and low-polarization data (see Planck Collaboration XIII 2016). BB power spectra are treated differently in determining the tensor contribution, since the model does not predict any primordial scalar fluctuations in BB. As a consequence, a primordial B-mode signal would be a direct signature of tensor modes. However, depending on the amplitude of the tensor-to-scalar ratio, such a signal may be masked by E-mode power transformed to B-mode power through lensing by gravitational potentials along the line of sight (so-called "BB lensing"). BB lensing has been measured with high accuracy by Planck in both harmonic (Planck Collaboration VIII 2020) and map (Planck Collaboration Int. XLI 2016) domains. But a primordial BB tensor signal has not yet been detected.
The scalar and tensor CMB angular power spectra are plotted in Fig. 1 for the Planck 2018 cosmology and for two values of the tensor-to-scalar ratio, namely r = 0.1 and r = 0.01. For further discussion of the tensor-to-scalar ratio and its implications for inflationary models, see Planck Collaboration XXII (2014), Planck Collaboration XX (2016), and Planck Collaboration X (2020). Note that the signal from tensor modes in EE is similar to that in BB modes, which makes EE (in particular at low multipoles) an important data set for tensor constraints. In this paper, we will make use of a polarized E-B likelihood, which consistently includes the correlated polarization fields E and B, and covers the range of multipoles where tensor modes can be constrained using Planck data (i.e., from = 2 to = 150).  At present the tightest B-mode constraints on r come from the BICEP/Keck measurements (BK15; BICEP2 Collaboration 2018), which cover approximately 400 deg 2 centred on RA = 0 h , Dec = −57. • 5. These measurements probe the peak of the Bmode power spectrum at around = 100, corresponding to gravitational waves with k ≈ 0.01 Mpc −1 that enter the horizon during recombination (i.e., somewhat smaller than the scales contributing to the Planck temperature constraints on r). The results of BK15 give a limit of r < 0.07 at 95 % confidence, which tightens to r < 0.06 in combination with Planck temperature and other data sets.
Planck Collaboration V (2020) presented Planck B-mode constraints from the 100-and 143-GHz HFI channels with a 95 % upper limit of r < 0.41 (at a somewhat larger pivot scale, as described in the next section), using only a limited number of multipoles around the so-called "reionization bump" (2 ≤ ≤ 29). Using Planck NPIPE maps (Planck Collaboration Int. LVII 2020), called Planck Release 4 (PR4), we are now able to constrain the BB power spectrum for a much larger number of modes, including both the reionization bump at large angular scales ( < ∼ 30) and the so-called "recombination bump" at intermediate scales (50 < ∼ < ∼ 150). In this paper, we first describe, in Sect. 2, the cosmological model used throughout the analysis. We then detail the data and the likelihoods in Sect. 3. Section 4 focuses on constraints from T T and in particular the impact of the low-data in temperature. Section 5 gives constraints from the BB angular power spectrum using Planck data, while results from the full set of polarization power spectra are given in Sect. 6. Finally, in Sect. 7, we combine all data sets to provide the most robust constraints on BB coming from Planck and in combination with other CMB data sets, such as the results from the BICEP/Keck Collaboration.

Cosmological model
We use the base-ΛCDM model, which has been established over the last couple of decades to be the simplest viable cosmological model, in particular with the Planck results (e.g., Planck Collaboration VI 2020). In this model, we assume purely adiabatic, scale-invariant perturbations at very early times, with curvature-mode (scalar) and tensor-mode power spectra parameterized by where A s and A t are the initial super-horizon amplitudes for curvature and tensor perturbations, respectively. The primordial spectral indexes for scalar (n s ) and tensor (n t ) perturbations are taken to be constant. This means that we assume no "running," i.e., a pure power-law spectrum with dn s /d ln k = 0. We set the pivot scale at k 0 = 0.05 Mpc −1 , which roughly corresponds to approximately the middle of the logarithmic range of scales probed by Planck; with this choice, n s is not strongly degenerate with the amplitude parameter A s . Note that for historical reasons, the definitions of n s and n t differ, so that a scale-invariant scalar spectrum corresponds to n s = 1, while a scale-invariant tensor spectrum corresponds to n t = 0. The late-time parameters, on the other hand, determine the linear evolution of perturbations after they re-enter the Hubble radius. We use the basis (Ω b , Ω c , θ * , τ) following the approach in Planck cosmological studies (Planck Collaboration VI 2020), where Ω b is the density parameter of baryons, Ω c is the density parameter of cold dark matter, θ * is the observed angular size of the sound horizon at recombination, and τ is the reionization optical depth.
The amplitude of the small-scale linear CMB power spectrum is proportional to A s e −2τ . Because Planck measures this amplitude very accurately, there is a tight linear constraint between τ and ln A s . For this reason, we usually adopt ln A s as a base parameter with a flat prior; ln A s has a significantly more Gaussian posterior than A s . A linear parameter redefinition then allows the degeneracy between τ and A s to be explored efficiently. Note that the degeneracy between τ and A s is broken by the relative amplitudes of large-scale temperature and polarization CMB anisotropies and by the effect of CMB lensing.
We define r ≡ A t /A s , the primordial tensor-to-scalar ratio defined explicitly at the scale k 0 = 0.05, Mpc −1 . Our constraints are only weakly sensitive to the tensor spectral index, n t (which is assumed to be close to zero). Note that the Planck Collaboration also discussed r constraints for k 0 = 0.002, Mpc −1 (Planck Collaboration V 2020). Given the definitions in Eqs. (1) and (2), the relation scales like (0.05/0.002) −r/8 , which means that r 0.002 is lower by 4 % at r 0.1 compared to r 0.05 and less than 0.4 % lower for r < 0.01.
In this work, we will use an effective tensor-to-scalar ratio r eff , which we extend into the negative domain by modifying the Boltzmann-solver code CLASS 1 (Blas et al. 2011). While negative tensor amplitudes are unphysical, this approach will allow us to derive posteriors without boundaries, facilitating detection of potential biases, and enabling us to determine a more accurate statistical definition of the constraints on r. With r eff we will be able to independently discuss both the uncertainty of r (σ r ) and corresponding upper limits (depending on the maximum a posteriori probability). In the rest of this paper, we will simply write r as the effective tensor-to-scalar ratio, and report upperlimits for positive tensor amplitudes, for which r eff = r. We use 95 % confidence levels when reporting upper limits, and a 68 % confidence interval with the maximum a posteriori probability.

Data and simulations
The sky measurements used in this analysis are the PR4 maps available from the Planck Legacy Archive 2 (PLA). They have been produced with the NPIPE processing pipeline, which creates calibrated frequency maps in temperature and polarization from the Planck Low Frequency Instrument (LFI) and High Frequency Instrument (HFI) data. As described in Planck Collaboration Int. LVII (2020), NPIPE processing includes 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.
NPIPE achieves an overall lower noise level in part by incorporating the data acquired during the 4-minute spacecraft repointing manoeuvres that take place between the 30-to-70-min stable science scans. Residual systematics are suppressed using a mitigation strategy that combines aspects of both LFI and HFI processing pipelines. Most importantly, gain fluctuations, bandpass mismatch, and other systematics are formulated into time-domain templates that are fitted and subtracted as a part of the mapmaking process. Degeneracies between sky polarization and systematic templates are broken by constructing a prior of the polarized foreground sky using the extreme polarizationsensitive frequencies (30,217,and 353 GHz).
Moreover, the PR4 release comes with 400 simulations of signal, noise, and systematics, component-separated into CMB maps, which allow for an accurate characterization of the noise and systematic residuals in the Planck maps. This is important because Planck polarization data are cosmic-variancedominated only for a few multipoles at very large scales in EE ( < 8, as shown in Fig. 2). These simulations, even though limited in number, represent a huge effort in terms of CPU time. They are essential in order to compute the following two additional quantities.
1. The end-to-end transfer function from the data reduction (including TOI processing, mapmaking and template fitting for mitigation of systematics, component separation, and power-spectrum estimation). The transfer function is defined as the ratio between the output and the input, averaged over all the simulations (see section 4.3 of Planck Collaboration Int. LVII 2020, for the details). 2. The covariance of the data (here we will use the cross-power spectra), which is the only way to propagate uncertainties when those are dominated by systematics (from the instrument or from foregrounds).
Note that these two quantities estimated from the simulations are directly related to two different characteristics of the final parameter posteriors: the bias of the mean (the transfer function); and the width of the posterior (as propagated into parameter constraints by the covariance matrix in the likelihood). They can be separated from each other, meaning that one systematic effect 2 https://pla.esac.esa.int can easily produce a significant bias without any strong impact on the variance, while another effect can produce a large increase of the variance with no associated bias.  The NPIPE simulations include the systematic effects relevant for polarization studies, specifically analogue-to-digitalconverter nonlinearities, gain fluctuations, bandpass mismatch between detectors, correlated noise (including 4-K line residuals), and full-beam convolutions for each detector.
The use of a polarization prior in NPIPE processing causes a suppression of large-scale ( < 20) CMB polarization, which needs to be corrected. As explained in Planck Collaboration Int. LVII (2020), allowing for a non-trivial transfer function is a compromise between measuring very noisy but unbiased large-scale polarization from all low-modes, and filtering out the modes that are most affected by the calibration uncertainties left in the data by the Planck scan strategy. As detailed in Planck Collaboration Int. LVII (2020), the transfer function to correct for this bias is determined from simulations. It is then used to correct the power spectrum estimates, just as instrumental beam and pixel effects must be deconvolved. Due to the fact that E modes dominate the CMB polarization, the simulations do not yield a definitive measurement of the Bmode transfer function. We have nevertheless chosen to conservatively deconvolve the E-mode transfer function from the Bmode spectrum in order to provide a robust upper limit on the true B-mode spectrum. Indeed, in the situation where primordial B-mode power is not detected, the transfer function correction essentially increases the variance estimate at low multipoles, which propagates the uncertainty induced by the degeneracy between the sky and the systematic templates used in NPIPE. Note that this uncertainty is small compared to the impact of systematics in the error budget (see Fig. 2).
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 detector-set (hereafter "detset") maps, comprising two subsets of maps at each frequency, with nearly independent noise characteristics, made by combining half of the detectors. Note that time-split maps (made from, e.g., "odd-even rings" or "halfmission data") share the same instrumental detectors, and therefore exhibit noise correlations due to identical spectral bandpasses and optical responses. Using time-split maps would result in systematic biases in the cross-power spectra, as well as underestimation of the noise properties when computing the halfdifferences.
Uncertainties at the power-spectrum level are dominated by noise and systematics, as illustrated in Fig. 2. Thanks to the NPIPE processing, we are now able to show the impact of the systematics at low . This is illustrated by comparing the PR4 end-to-end noise (based on the Monte Carlo simulations, including instrumental noise, systematics, and foreground uncertainties, and corrected for the transfer function both in EE and BB) with the propagation of the statistical noise coming from the analytic pixel-pixel covariance matrix. The systematic uncertainties dominate at < ∼ 15, then slowly decrease so that the effective uncertainties converge towards the analytic estimate at higher multipoles.

Polarized sky masks
Foreground residuals in the foreground-cleaned maps dominate the polarized CMB signal near the Galactic plane. To avoid contamination from these residuals in the cosmological analysis, we mask the Galactic plane. We use a series of different retained sky fractions (from 30 % to 70 %) to check the consistency of our results with respect to foreground residuals (Fig. 3). The masks used in this analysis are a combination of a mask for polarization intensity (to avoid polarized foreground residuals), a mask for total intensity (to avoid potential temperatureto-polarization leakage residuals), and the confidence mask for component separation provided by the Planck Colaboration. The intensity mask is obtained by thresholding the 353-GHz intensity map (which traces dust) scaled to 143 GHz, and the 30-GHz intensity map (which traces synchrotron) scaled to 100 GHz. The polarization map is constructed similarly. Both foreground trac-ers are smoothed beforehand with a 10 • Gaussian window function.
The impact of the emission of extragalactic polarized sources on the power spectra is negligible, given the Planck resolution and noise level. The confidence mask for component separation ensures the masking of the strongest sources, which could also produce residuals through temperature-to-polarization leakage. Table 1 summarizes the likelihoods used in this analysis, which are described below.

Low-temperature likelihood
We use the Planck public low-temperature-only likelihood based on the PR3 CMB map recovered from the componentseparation procedure (specifically Commander) described in detail in Planck Collaboration V (2020). At large angular scales, Planck temperature maps are strongly signal-dominated, and there is no expected gain in updating this likelihood with the PR4 data.
As discussed in Planck Collaboration XX (2016), the lowtemperature data from Planck have a strong impact on the r posterior and the derivation of the corresponding constraints. This is because the deficit of power in the measured C s at low-in temperature (see the discussions in Planck Collaboration XVI 2014 and Planck Collaboration Int. LI 2017) lowers the probability of tensor models, which add power at low multipoles. This shifts the maximum in the posterior of r toward low values (or even negative values when using r eff , as we will show in Sect. 4).

High-likelihood
At small angular scales ( > 30), we use the HiLLiPoP likelihood, which can include the T T , T E, and/or EE power spectra. HiLLiPoP has been used as an alternative to the public Planck likelihood in the 2013 and 2015 Planck releases (Planck Collaboration XV 2014; Planck Collaboration XI 2016), and is described in detail in Couchot et al. (2017a). In this paper, the HiLLiPoP likelihood is applied to the PR4 detset maps at 100, 143, and 217 GHz. We focus on the T T spectra, since there is no additional information at small scales in T E or EE for tensor modes. We only make use of T E in Sect. 7 in order to help constrain the spectral index n s . The likelihood is a spectrum-based Gaussian approximation, with semi-analytic estimates of the C covariance matrix based on the data. The cross-spectra are debiased from the effects of the mask and the beam leakage using Xpol (a generalization to polarization of the algorithm presented in Tristram et al. 2005 3 ) before being compared to the model, which includes CMB and foreground residuals. The beam window functions are evaluated using QuickPol (Hivon et al. 2017), adapted to the PR4 data. These adaptations include an evaluation of the beam-leakage effect, which couples temperature and polarization modes due to the beam mismatch between individual detectors.
The model consists of a linear combination of the CMB power spectrum and several foregrounds residuals. These are: -Galactic dust (estimated directly from the 353-GHz channel);  (2013); a tSZ-CIB correlation consistent with both models above; and unresolved point sources as a Poisson-like power spectrum with two components (extragalactic radio galaxies and infrared dusty galaxies).
On top of the cosmological parameters associated with the computation of the CMB spectrum, with HiLLiPoP we sample seven foreground amplitudes (one per emission source, the spectral energy density rescaling the amplitude for each crossfrequency being fixed) and six nuisance parameters (one overall calibration factor plus intercalibrations for each map). See Appendix A for more details.

Large-scale polarized likelihood
We construct a polarized E-B likelihood based on power spectra, focusing on the large scales where the tensor signal is dominant.
Because it carries very little information about the tensor modes, we do not include the T E spectrum in this analysis.
In polarization, especially at large angular scales, foregrounds are stronger relative to the CMB than in temperature, and cleaning the Planck frequencies using C templates in the likelihood (as done in temperature) is not accurate enough. 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, and thermal dust emission. Commander was run on each detset map independently, as well as on each realization from the PR4 Monte Carlo simulations. Maps are available on the PLA in HEALPix 4 format (Górski et al. 2005) at a resolution N side = 2048.
To compute unbiased estimates of the angular power spectra, we calculate the cross-correlation of the two detset maps. We make use of two different angular cross-power spectra estimators (described below), which are then concatenated to produce a full-multipole-range power spectrum. There is no information loss in this process, since the covariances are deduced using Monte Carlo simulations including the correlations over the entire multipole range.
-For multipoles 2 ≤ ≤ 35, 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). 5 At multipoles below 40, it has been shown to produce unbiased polarized power spectra with almost optimal errors. We use downgraded N side = 16 maps after convolution with a cosine apodizing kernel b = 1 2 1 + cos π l−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 in the mapmaking procedure (see Sect. 3.1).
-For multipoles 35 < < 300, we compute power spectra with a classical pseudo-C estimator Xpol (Sect. 3.3.2). We used N side = 1024 maps and the native beam of Commander maps (i.e., 5 ). In this case, we apodize the mask (see Sect. 3.2) with a 1 • Gaussian taper. Given the low signal-tonoise ratio in polarization, we bin the spectra with ∆ = 10.
The EE, BB, and EB power spectra estimates are presented in Fig. 4 for 50 % of the sky, which provides the best combination of sensitivity and freedom from foreground residuals. Power spectra computed on different sky fractions (using masks from Sect. 3.2) are compared in Fig. B.1. A simple χ 2 test on the first 34 multipoles shows no significant departure from the Planck 2018 ΛCDM model for any of these spectra. The most extreme multipole in the BB spectrum is = 6, which, for a Gaussian distribution, would correspond conditionally to a 3.4 σ outlier (reducing to 2.3 σ after taking into account the look-elsewhere effect, including the first 34 multipoles). However, at such low multipoles, the distribution is not Gaussian and the "probability to exceeed" values are certainly higher than the numbers of σ would suggest. In EE, the largest deviation from the model is for = 17 at 3.1 σ and in EB it is = 19 at 2.7 σ.
The C covariance matrix is computed from the PR4 Monte Carlos. For each simulation, we compute the power spectra using both estimators. The statistical distribution of the recovered C then naturally includes the effect of the components included in the Monte Carlo, namely the CMB signal, instrumental noise, Planck systematic effects incorporated in the PR4 simulations (see Sect. 3.1), component-separation uncertainties, and foreground residuals. The residual power spectra (both for the simulations and the data) are shown in Fig. B.2. Given the Planck noise level in polarization, we focus on multipoles below = 150, which contain essentially all the information on tensor modes in the Planck CMB angular power spectra. At those scales, and given Planck noise levels, the likelihood function needs to consistently take into account the two polarization fields E and B, as well as all correlations between multipoles and modes (EE, BB, and EB).
LoLLiPoP (LOw-LIkelihood on POlarized Power-spectra) is a Planck low-polarization likelihood based on cross-spectra, and was previously applied to Planck EE data for investigating the reionization history in Planck Collaboration Int. XLVII  Figure 4: EE, BB, and EB power spectra of the CMB computed on 50 % of the sky with the PR4 maps at low (left panels) and intermediate multipoles (right panels). Grey bands represent the associated cosmic variance. Error bars are deduced from the PR4 Monte Carlo simulations. Correlations between data points are given in Appendix C. A simple χ 2 test shows no significant departure from the model for any of these spectra.
(2016). The version used here is updated to use cross-spectra calculated on component-separated CMB detset maps processed by Commander from the PR4 frequency maps. Systematic effects are considerably reduced in cross-correlation compared to auto-correlation, and LoLLiPoP is based on crosspower 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 model, and O are the offsets needed in the case of cross-spectra. For multi-dimensional CMB modes (here we restrict ourselves to E and B fields only), the C generalise to C , a 2 × 2 matrix of power spectra, and the g function is applied to the eigenvalues of C −1/2 C C −1/2 (with C −1/2 the square root of the positive-definite matrix C).
In the case of auto-spectra, the offsets O are given by the noise bias effectively present in the measured power spectra. For crosspower spectra, the noise bias is zero, and we use effective offsets defined from the C noise variance: The distribution of the new variable X ≡ vecp(X ), the vector of distinct elements of 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. The resulting C covariance thus consistently includes CMB sample variance, statistical noise, and systematic residuals, as well as foreground-cleaning uncertainties, together with the correlations induced by masking. These uncertainties are then propagated through the likelihood up to the level of cosmological parameters. Figures of the correlation matrices are given in Appendix C.
Using this approach, we are able to derive three different likelihoods, one using only information from E modes (lowlE), one using only information from B modes (lowlB), and one using EE+BB+EB spectra (lowlEB). We have used these likelihoods from = 2 up to = 300 with a nominal range of = [2, 150], since multipoles above 150 do not contribute to the result due to the Planck noise (see Sect. 5).
The approach used in this paper is different from the one used for the Planck 2018 results. Indeed, in Planck Collaboration V (2020), the probability density of the polarized spectra at low multipoles was modelled with a polynomial function adjusted on simulations in which only τ is varied, with all other cosmological parameters in a ΛCDM model fixed to the Planck 2018 best-fit values. As a consequence, the probability density is not proportional to the likelihood L(Ω model |C data ) when the model is not ΛCDM (and in particular for our case ΛCDM+r), and even in the ΛCDM case it neglects correlations with other parameters that affect the posterior on τ. In addition, the simulations used in Planck Collaboration V (2020) were generated with the same CMB realization for the mapmaking solution. Cosmic variance was included afterwards by adding CMB realizations on top of noise-only maps, neglecting correlations between foregrounds or systematic templates and the CMB. The information in polarization at low-was then extracted using a polynomial function fitted to the distribution from simulations. While this is supposed to empirically take into account the effects of systematics on the likelihood shape, it does not include -by-correlations, and is limited in the C power that one can test (for example imposing a strong prior on the EE power at = 3). As a consequence, the combination of those two effects reduces the covariance, especially at low multipoles, leading to error bars (especially on τ) that are underestimated.

Constraints from TT
To derive constraints on the tensor-to-scalar ratio from the temperature power spectrum, we use the high-HiLLiPoP likelihood for 30 ≤ ≤ 2500, and the Commander likelihood (lowT) in temperature for < 30, with a prior on the reionization optical depth to break the degeneracy with the scalar amplitude A s . We use a Gaussian prior τ = 0.055 ± 0.009. For the base-ΛCDM model, using PR4 data, we obtain the same results as presented in Planck Collaboration VI (2020).
We now describe the results obtained when fitting the tensorto-scalar ratio r in addition to the six ΛCDM parameters (Ω b h 2 , Ω c h 2 , θ * , A s , n s , τ). In Planck Collaboration X (2020), the constraint from T T is reported as r 0.002 < 0.10 (95 % CL) using PR3 data. This is much lower than the expected 2 σ upper bound on r. Indeed, when we calculate r eff as proposed in Sect. 2, we find that the maximum of the posterior is in the negative region by about 1.7 σ. That the maximum happens to fall at negative values is the major reason for the apparently strong constraint on r.
With PR4 data, we find that the maximum of the posterior is negative by less than 1.2 σ when using HiLLiPoP in temperature (hlpTT) along with lowT. As discussed in Planck Collaboration X (2020), this result is related to the lowdeficit in the temperature power spectrum. Indeed, removing lowT from the likelihood moves the maximum of the posterior closer to zero, as illustrated in Fig. 5. The corresponding posterior maximum and 68 % confidence interval are r 0.05 = +0.031 ± 0.120 (hlpTT+τ-prior), (7a) r 0.05 = −0.131 ± 0.093 (hlpTT+lowT), (7b) r 0.05 = −0.101 ± 0.094 (hlpTT+lowT+τ-prior). (7c) Using the temperature power spectrum from PR4, we recover the same constraints on other parameters, in particular the scalar spectral tilt n s , as found using PR3 data (see Appendix D). With the full posterior distribution on r, we get a better idea of the intrinsic sensitivity of the Planck temperature data to r. Using only high-data, with a prior on the reionization optical depth τ, we find σ r = 0.12 for T T . Note that we find σ r = 0.43 for T E, indicating that T E is much less constraining for r than T T . When adding information from low multipoles in temperature, σ r reduces to 0.094, but at the price of pushing the maximum distribution towards negative values. The fact that the distribution peaks in the non-physical domain can be considered as a statistical fluctuation (with a significance between 1 and 2 σ, depending on the data set used), which on its own is not a serious problem. However, the fact that this behaviour is strongly related to the deficit of power at low-in temperature is worth noting.

Constraints from BB
To derive constraints on the tensor-to-scalar ratio from BB using the PR4 maps, we sample the likelihood with a fixed ΛCDM model based on the Planck 2018 best fit, to which we add tensor fluctuations with a free amplitude parametrized by the tensorto-scalar ratio r. We use the LoLLiPoP likelihood described in Sect. 3.3.3, restricted to BB only (referred as "lowlB"). As discussed in Sect. 3.3.3, we construct the C covariance matrix using the PR4 Monte Carlo simulations, which include CMB signal, foreground emission, realistic noise, and systematic effects. Before giving the final constraints coming from the Planck BB spectra, we should distinguish between the two different regimes, corresponding to large scales (the reionization bump) and intermediate scales (the recombination bump). Across the reionization bump, uncertainties are dominated by systematic residuals, as discussed in Sect. 3.1, while foreground residuals may bias the results. Across the recombination bump, uncertainties are dominated by statistical noise; however, systematic effects, as well as foreground residuals, can still bias constraints on r. In order to test the effects of potential foreground residuals, we calculate the posterior distributions of r using various Galactic masks, as described in Sect. 3.2. While large sky fractions ( f sky > 60 %) show deviations from r = 0, the posteriors for 40, 50, and 60 % of the sky are consistent with zero ( Fig. E.1). As a robustness test, we also calculate the posterior distribution when changing the range of multipoles (Fig. E.2) and find consistent results, with posteriors compatible with r = 0. Multipoles above 150 do not contribute to the result, since the noise in BB is too high. For the rest of this paper, unless otherwise noted, we use a sky fraction of 50%, and compute the likelihood over the range of multipoles from = 2 to = 150.
Both results are obtained over 50 % of the sky, with multipoles in the range = [2, 35] for the former and = [50, 150] for the latter. With these ranges of multipoles, and given the statistics of the PR4 maps, we can see that the reionization bump (σ r = 0.110) and the recombination bump (σ r = 0.113) contribute equally to the overall Planck sensitivity to the tensor-toscalar ratio. We can combine the results from the two bumps in order to give the overall constraints on the tensor-to-scalar ratio from the Planck BB spectrum (Fig. 6). The full constraint on r from the PR4 BB spectrum over 50 % of the sky, including correlations between all multipoles between = 2 and = 150, is r 0.05 = 0.033 ± 0.069 (lowlB).
This is fully compatible with no tensor signal, and we can derive an upper limit by integrating the posterior distribution out to 95 %, after applying the physical prior r > 0, which yields r 0.05 < 0.158 (95 % CL, lowlB).
This result can be compared with the BICEP2/Keck Array constraints (BICEP2 Collaboration 2018) of r 0.05 < 0.072 (95 % CL, BK15), with σ r = 0.02 compared to σ r = 0.069 for the Planck result presented in this analysis

Additional constraints from polarization
As shown in Fig. 1, the EE tensor spectrum is similar in amplitude to the BB tensor spectrum, even though the scalar mode in EE is stronger. Given that noise dominates the tensor signal at all multipoles in both EE and BB, we expect the likelihood for EE to give useful constraints on r. We thus present the constraints from polarized low-data ( < 150) using different combinations of the LoLLiPoP likelihood (specifically EE, BB, and EE+BB+EB) in Fig. 7. We emphasize that EE+BB+EB is a likelihood of the correlated polarization fields E and B and not the combination of individual likelihoods (see Sect. 3.3.3).
The first thing to notice is that the posterior distribution for EE peaks at r = 0.098±0.097, while the other modes give results compatible with zero within 1 σ. Given the lower sensitivity of lowlE to r (σ r 0.10) compared to that of lowlB (σ r 0.07), this is mitigated when adding the information from other modes. The posterior distributions for r give r 0.05 = 0.033 ± 0.069 (lowlB), (14) r 0.05 = −0.031 ± 0.046 (lowlEB).
As a consistency check, Fig. 7 also shows the constraints when fitting the BB tensor model on the EB data power spectrum, which is compatible with zero (r = −0.012 ± 0.068) as expected.  : Posterior distributions for r from Planck polarized low-data ( < 150) using LoLLiPoP and the EE, BB, and EE+BB+EB spectra. The dashed black line is obtained from EB data by fitting a BB tensor model. The sky fraction used here is f sky = 50 % Using polarization data, Planck's sensitivity to the tensor-toscalar ratio reaches σ r = 0.046. Combining all Planck polarization modes (EE, BB, and EB) out to = 150 leads to the following upper limit: r 0.05 < 0.069 (95 % CL, lowlEB).
Note that this constraint is almost independent of the other ΛCDM parameters, and in particular the reionization optical depth τ. To demonstrate this, using the same data set (lowlB and lowlEB), we derive 2-dimensional constraints for τ and r and plot them in Fig. 8. The constraint is stable when sampling for τ. Indeed, in this case, we obtain r 0.05 = 0.025 ± 0.064 (lowlB), (17) r 0.05 = −0.015 ± 0.045 (lowlEB), and for the reionization optical depth compatible with lowlE results, while lowlB shows no detection of τ, since BB is dominated by noise.

Combined results
Up to this point, the constraints on r have been derived relative to a fixed fiducial ΛCDM spectrum based on the Planck 2018 results. Including the Planck temperature likelihoods (both lowT and hlpTT) in a combined analysis of the Planck CMB spectra allows us to properly propagate uncertainties from other cosmological parameters to r, as well as to self-consistently derive constraints in the n s -r plane. In this section, we combine the lowT and hlpTT with the low-polarized likelihood lowlEB to sample the parameter space of the ΛCDM+r model. The comparison of contours at 68 % and 95 % confidence levels between PR3 and PR4 data is presented in Fig. F.1 of Appendix F.
We also include the BK15 constraints from BICEP2 Collaboration (2018). When combining Planck and BK15, we neglect the correlation between the two data sets and simply multiply the likelihood distributions. This is justified because the BK15 spectra are estimated on 1 % of the sky, while the Planck analysis is derived from 50 % of the sky.  Figure 10 shows the constraints in the r-n s plane for Planck data in combination with BK15. The constraints from the full combination of Planck data are comparable to those from BK15. The addition of the high-T E likelihood produces tighter constraints on the spectral index n s (as already reported in Planck Collaboration VI 2020).

Conclusions
In this paper, we have derived constraints on the amplitude of tensor perturbations using Planck PR4 data. We investigated the intrinsic sensitivity of the T T spectrum, which is cosmicvariance limited, and found σ r = 0.094 using the full range of multipoles. We noted the impact of the low-anomaly, which pushes the maximum posterior distribution toward negative values of r eff at roughly the 1 σ level.
For the first time, we analyzed the Planck BB spectrum for r and obtained σ r = 0.069, which is lower than in temperature. The Planck B-mode spectrum, being dominated by noise, gives a constraint on r that is fully compatible with zero from both low and intermediate multipoles, in other words from both the reionization and recombination peaks. Multipoles above 150 do not contribute to the result, since the noise in BB is too high.
Using an appropriate likelihood in polarization, we showed that the Planck EE spectrum is also sensitive to the amplitude of the tensor-to-scalar ratio r. The combined constraints from Planck EE and BB, including EB correlations, lead to a sensitivity on r of σ r = 0.046, two times better than in temperature. We also investigated the impact of foreground residuals using different Galactic cuts and by varying the range of multipoles used in the polarized likelihood. Finally, by combining temperature and polarization constraints, we derived the posterior distribution on r marginalized over the ΛCDM cosmological parameters as well as the nuisance parameters from the T T highlikelihood. The result gives an upper-limit of r < 0.056 at the 95 % confidence level using Planck data only. In combination with the BICEP/Keck measurements from 2015, this constraint is further reduced to r < 0.044 (95 % CL), the tightest limit on r to date. the full sky) in grey. The spectra show a remarkable consistency with the model, and are largely insensitive to sky fraction. The BB and EB spectra are dominated by noise. At low multipoles, the BB spectrum computed on the largest sky fraction (70 %) exhibits an excess at very low multipoles ( ≤ 5), attributable to Galactic residuals. For both BB and EB, a reduction of the sky fraction corresponds to a larger dispersion.

Appendix C: Cross-spectrum correlation matrix
Uncertainties are propagated to the likelihood function through the C covariance matrices (see Sect. 3.3.3). We use Monte Carlo simulations to estimate the covariance over the multipoles considered in this analysis (2 ≤ ≤ 150). Four hundred simulations are available, so the covariance matrices are accurate at the 5 % level. This is enough to propagate the uncertainties of the C s to r as well as to τ. Figure C.1 shows the correlations of the C BB computed from the covariance matrix. Correlations from bin to bin are below 15 %, except for the next-to-neighbour bins at < 40. Figure C.2 shows the correlations computed from the full covariance matrix, including EE, BB, and EB spectra for all multipole bins from 2 to 150.

Appendix D: ΛCDM+r parameters for PR3 and PR4
In this paper, we use a new data set, PR4, and an alternative likelihood, HiLLiPoP. Figure D.1 shows that we obtain essentially the same parameter values for the ΛCDM+r model based on the temperature power spectrum as are obtained from PR3 and plik (available on the Planck Legacy Archive). The differences between PR3 and PR4 are primarily seen in polarization.