Planck 2018 results. X. Constraints on inflation

We report on the implications for cosmic inflation of the 2018 Release of the Planck CMB anisotropy measurements. The results are fully consistent with the two previous Planck cosmological releases, but have smaller uncertainties thanks to improvements in the characterization of polarization at low and high multipoles. Planck temperature, polarization, and lensing data determine the spectral index of scalar perturbations to be $n_\mathrm{s}=0.9649\pm 0.0042$ at 68% CL and show no evidence for a scale dependence of $n_\mathrm{s}.$ Spatial flatness is confirmed at a precision of 0.4% at 95% CL with the combination with BAO data. The Planck 95% CL upper limit on the tensor-to-scalar ratio, $r_{0.002}<0.10$, is further tightened by combining with the BICEP2/Keck Array BK14 data to obtain $r_{0.002}<0.064$. In the framework of single-field inflationary models with Einstein gravity, these results imply that: (a) slow-roll models with a concave potential, $V"(\phi)<0,$ are increasingly favoured by the data; and (b) two different methods for reconstructing the inflaton potential find no evidence for dynamics beyond slow roll. Non-parametric reconstructions of the primordial power spectrum consistently confirm a pure power law. A complementary analysis also finds no evidence for theoretically motivated parameterized features in the Planck power spectrum, a result further strengthened for certain oscillatory models by a new combined analysis that includes Planck bispectrum data. The new Planck polarization data provide a stringent test of the adiabaticity of the initial conditions. The polarization data also provide improved constraints on inflationary models that predict a small statistically anisotropic quadrupolar modulation of the primordial fluctuations. However, the polarization data do not confirm physical models for a scale-dependent dipolar modulation.


Introduction
This paper, one of a set associated with the 2018 release of data from the Planck 1 mission, presents the implications for cosmic inflation of the 2018 Planck measurements of the cosmic microwave background (CMB) anisotropies. In terms of data, this paper updates Planck Collaboration XXII (2014) (henceforth PCI13), which was based on the temperature data of the nominal Planck mission ("PR1"), including the first 14 months of observations, and Planck Collaboration XX (2016) (henceforth PCI15), which used temperature data and a first set of polarization data from the full Planck mission ("PR2"), comprising 29 and 52 months of observations for the high frequency instrument (HFI) and low frequency instrument (LFI), respectively.
The ideas underlying cosmic inflation were developed during the late 1970s and early 1980s in order to remedy a number of defects of the hot big-bang cosmological model (e.g., the horizon, smoothness, flatness, and monopole problems) (Brout et al. 1978;Starobinsky 1980;Kazanas 1980;Sato 1981;Guth 1981;Linde 1982;Albrecht & Steinhardt 1982;Linde 1983). Subsequently, it was realized that, on account of quantum vacuum fluctuations, cosmic inflation also provides a means to generate the primordial cosmological perturbations (Mukhanov & Chibisov 1981Hawking 1982;Guth & Pi 1982;Starobinsky 1982;Bardeen et al. 1983;Mukhanov 1985). The development of cosmic inflation is one of the major success stories of modern cosmology, and in this paper we explore how the latest 2018 release of the Planck data constrains inflationary models.
Planck data currently provide the best constraints on the CMB anisotropies, except on very small angular scales beyond the resolution limit of Planck. The Planck data set has enabled a precision characterization of the primordial cosmological perturbations and has allowed cosmological parameters to be constrained at the sub-percent level. One of the main data products, described in more detail in the following section, is the Planck T T , T E, EE, and lensing power spectra, which are shown in Fig. 1, together with the residuals from the six-parameter concordance Λ cold dark matter (ΛCDM) model using the best-fit parameter values.
In order to provide a quantitative estimate of the improvement achieved by Planck, as well as to show where Planck stands compared to an ultimate cosmic-variance-limited survey, we consider an idealistic estimator for the number of modes [i.e., the effective number of a m 's measured (Planck Collaboration I 2016)]: where C XY (∆C XY ) is the (error on the) angular power spectrum of the XY channel (Planck Collaboration I 2016; Scott et al. 2016). The number of modes measured by Planck is 1 430 000 and 109 000 for temperature (XY = T T up to = 2500) and polarization (XY = EE up to = 2000), respectively (Planck Collaboration I 2018). The number of modes measured is increased by approximately a factor of 7 (570) for temperature (polarization) with respect to the WMAP 9-year measurement, but there is still a factor of 3 (40) to gain for a cosmic-variancelimited experiment up to = 2500 accessing 70 % of the sky. The additional modes measured by Planck play a key role in improving the constraints on the initial conditions for the cosmological perturbations and on models of inflation with respect to previous measurements of CMB anisotropies.
Planck data have also greatly improved the constraints on bispectral non-Gaussianity, both for the "local" pattern, as predicted by many inflationary models, and for other templates such as the equilateral one, as analysed and reported in detail in the dedicated Planck non-Gaussianity papers (Planck Collaboration XXIV 2014;Planck Collaboration XVII 2016;Planck Collaboration IX 2018). The constraints on the non-Gaussianity parameter f NL are limited by a combination of cosmic variance and instrumental noise. An order-of-magnitude estimate for the signal-to-noise ratio for the local pattern (with f loc NL = 1) is given by S N For the local shape, the logarithm enters because most of the signal derives from detecting the modulation of the smallscale power by the large-scale CMB anisotropy, highlighting the importance of full-sky maps for this kind of analysis. For other shapes such as equilateral, one instead has (S /N) 2 ∼ Ω sky 2 max . Planck has significantly sharpened the constraints on f NL , largely on account of its measurement of high multipoles with higher signal-to-noise ratio compared to past surveys. Some improvement has also been obtained from including polarization.
The Planck measurements have significantly constrained the physics of inflation. The hypothesis of adiabatic Gaussian scalar fluctuations with a power spectrum described by a simple power law, which is the key prediction of the standard single-field slow-roll inflationary models, has been tested to unprecedent accuracy (PCI13; PCI15;Planck Collaboration XXIV 2014; Planck Collaboration XVII 2016). Planck has set tight constraints on the amount of inflationary gravitational waves by exploiting the shape of the CMB temperature spectrum (PCI13). These results have inspired a resurgence of activity in inflationary model building. For more details, see, for example, the following review articles and references therein: Linde (2015); Martin et al. (2014a); Guth et al. (2014); and Burgess et al. (2013). Planck analysis and interpretation have also sparked a debate on the likelihood of initial conditions for some inflationary models (Ijjas et al. 2013;Ijjas & Steinhardt 2016;Linde 2017), which is primarily of theoretical interest and is not addressed in this paper. In combination with more sensitive B-mode ground-based polarization measurements, as from BICEP-Keck Array (BKP), Planck has convincingly ruled out the slow-roll inflationary model with a quadratic potential (PCI15). In terms of physics beyond the simplest slow-roll inflationary models, the pre-Planck hints of a running spectral index (Hou et al. 2014) or of large non-Gaussianities (Bennett et al. 2012) have disappeared as a result of the Planck measurements. How to interpret anomalies on the largest angular scales and at high multipoles is a question motivating the search for new nonstandard inflationary models. We discuss how the Planck 2018 release data further test these ideas.
This paper is organized as follows. In Sect. 2 we describe the statistical methodology, the Planck likelihoods, and the complementary data sets used in this paper. In Sect. 3 we discuss the updated constraints on the spectral index of the scalar perturbations, on spatial curvature, and on the tensor-to-scalar ratio.  For each panel we also show the residuals with respect to this baseline best fit. Plotted are D = ( + 1)C /(2π) for T T and T E, C for EE, and L 2 (L + 1) 2 C φφ L /(2π) for lensing. For T T , T E, and EE, the multipole range 2 ≤ ≤ 29 shows the power spectra from Commander (T T ) and SimAll (T E, EE), while at ≥ 30 we display the co-added frequency spectra computed from the Plik cross-half-mission likelihood, with foreground and other nuisance parameters fixed to their best-fit values in the base-ΛCDM cosmology. For the Planck lensing potential angular power spectrum, we show the conservative (orange dots; used in the likelihood) and aggressive (grey dots) cases. Note some of the different horizontal and vertical scales on either side of = 30 for the temperature and polarization spectra and residuals. Section 4 is devoted to constraining slow-roll parameters and to a Bayesian model comparison of inflationary models, taking into account the uncertainties in connecting the inflationary expansion to the subsequent big-bang thermalized era. In Sect. 5 the potential for standard single-field inflation is reconstructed using two different methodologies. Section 6 describes the primordial power spectrum reconstruction using three different approaches. In Sect. 7, the parametric search for features in the primordial scalar power spectrum is described, including a dedicated study of the axion monodromy model. In Sect. 8, the Planck power spectrum data are combined with information from the Planck bispectrum in a search for oscillations in the primordial spectra. The constraints on isocurvature modes are summarized in Sect. 9. Section 10 updates and extends the constraints on anisotropic inflationary models of inflation. We summarize our conclusions in Sect. 11, highlighting the key results and the legacy of Planck for inflation.

Methodology and data
The general theoretical background and analysis methods applied in this paper closely match those of the previous Planck inflation papers (PCI13; PCI15). Consequently, in this section we provide only a brief summary of the methodology and focus on changes in the Planck likelihood relative to previous releases.

Cosmological models and inference
For well over a decade, the base-ΛCDM model has been established as the simplest viable cosmological model. Its six free parameters can be divided into primordial and late-time parameters. The former describe the state of perturbations on observable scales (corresponding to a wavenumber range of 10 −4 Mpc −1 k 10 −1 Mpc −1 today) prior to re-entering the Hubble radius around recombination. In base ΛCDM, the initial state of perturbations is assumed to be purely adiabatic and scalar, with the spectrum of curvature perturbations given by the power law ln P R (k) = ln A s + (n s − 1) ln(k/k * ) ≡ ln P 0 (k), where k * denotes an arbitrary pivot scale. The late-time parameters, on the other hand, determine the linear evolution of perturbations after re-entering the Hubble radius. By default we use the basis (ω b , ω c , θ MC , τ) 2 for the late-time parameters, but occasionally also consider non-minimal late-time cosmologies. Because of the inflationary perspective of this paper, we are mainly interested in exploring modifications of the primordial sector and their interpretation in terms of the physics of the inflationary epoch.
Perturbations produced by generic single-field slow-roll models of inflation are typically well approximated by the following form of the adiabatic scalar and tensor components: ln P R (k) = ln P 0 (k) + 1 2 d ln n s d ln k ln(k/k * ) 2 + 1 6 d 2 ln n s d ln k 2 ln(k/k * ) 3 + . . . , ln P t (k) = ln(rA s ) + n t ln(k/k * ) + . . . , which allows for a weak scale dependence of the scalar spectral index, n s , modelled by a running, d ln n s /d ln k, or a running of the running, d 2 ln n s /d(ln k) 2 . 3 The power spectrum parameterization in Eq. 4 can be extended to address wider classes of inflation-related questions, (e.g., the search for isocurvature perturbations, specific primordial features in the spectra, etc.), as described in subsequent sections. We also go beyond simple functions to parameterize the primordial power spectrum.
In the spirit of reconstructing the primordial spectrum from the data, we consider some general parameterizations (e.g., taking the power spectrum as an interpolation between knots of freely varied amplitudes at fixed or varying wavenumbers). One could argue that the primordial power spectra are merely intermediate quantities and assess theories directly from more fundamental parameters. By using the slow-roll approximation, or by evolving the mode equations to obtain exact numerical predictions for the spectra without resorting to the slow-roll approximation, we can relate the primordial perturbations to the dynamics of the Hubble parameter during inflation or to the inflaton potential and its derivatives, thus constraining these quantities directly.
For any given model, theoretical predictions of CMBrelated and other cosmological observables are calculated using appropriately modified versions of the Boltzmann codes CAMB (Lewis et al. 2000) or CLASS (Blas et al. 2011). As in (PCI13; PCI15), we compare models M 1 and M 2 by the difference in the logarithm of the likelihood of their best fits, or effective ∆χ 2 ≡ 2 [ln L max (M 1 ) − ln L max (M 2 )]. We apply Bayesian statistical methods to infer the posterior probability distributions of the model parameters and select between competing models (Trotta 2008), using either the Metropolis-Hastings Markov-chain Monte Carlo (MCMC) sampling algorithm, as implemented in CosmoMC (Lewis & Bridle 2002) and Monte Python (Audren et al. 2013), or software based on nested sampling (Skilling 2004), such as MultiNest (Feroz et al. 2009(Feroz et al. , 2013 or PolyChord (Handley et al. 2015). The latter can simultaneously estimate the Bayesian evidence E i of a model M i , allowing the comparison between different models via the Bayes factor, B = E 2 /E 1 , where |ln B| > 5 is commonly considered 2 Refer to Table 1 for definitions. 3 Unless explicitly stated otherwise, we adopt a default pivot scale k * = 0.05 Mpc −1 in this work. As in previous Planck releases, we will also quote the tensor-to-scalar ratio r 0.002 at k * = 0.002 Mpc −1 in order to facilitate comparison with earlier primordial tensor-mode constraints. "strong" evidence for or against the respective model (Jeffreys 1998;Trotta 2007a).

Planck data
The Planck data processing has improved in a number of key aspects with respect to the previous 2015 cosmological release. We briefly summarize the main points here, referring the interested reader to Planck Collaboration II (2018) and Planck Collaboration III (2018) for details.
The flagging procedure in the LFI 2018 pipeline has been made more aggressive, in particular for the first 200 operational days. However, the most important improvement in the LFI pipeline is in the calibration approach. Whereas in the 2015 release, the main calibration source for LFI was the Planck orbital dipole (i.e., the amplitude modulation of the CMB dipole induced by the satellite orbit) of each single radiometer model, the 2018 procedure also includes the Galactic emission along with the orbital dipole in the calibration model and becomes iterative (Planck Collaboration II 2018).
The HFI data for the 2018 release have also been made more conservative, cutting approximately 22 days of observations under non-stationary conditions with respect to 2015. The main change in the HFI data processing is the use of a new map making and calibration algorithm called SRoll, whose first version was introduced in Planck Collaboration Int. XLVI (2016) for the initial analysis of HFI polarization on large angular scales. This algorithm employs a generalized polarization destriper which uses the redundancy in the data to extract several instrumental systematic-effect parameters directly from the data (Planck Collaboration III 2018).
These improvements have a minor impact on Planck temperature maps, but are much more important for polarization, particularly on large angular scales, allowing, for instance, the removal of the high-pass filtering adopted in the 2015 study of isotropy and statistics Planck Collaboration XVI (2016).
In the following, we summarize the essentials of the Planck inputs used in this paper (i.e., the Planck likelihood approach to the information contained in the 2-point statistics of the temperature and polarization maps and the Planck CMB lensing likelihood). As for previous cosmological releases, the Planck likelihood approach is hybridized between low-and high-multipole regions, which therefore are summarized separately below. We refer the interested reader to the relevant papers Planck Collaboration V (2018) (henceforth PPL18) and Planck Collaboration VIII (2018) (henceforth PPLe18) for a more complete description of these inputs.

Planck low-likelihood
As in the Planck 2015 release, several options are available for evaluating the temperature likelihood on large angular scales, each with its own computational complexity and approximations. One option is based on the Commander framework and implements full Bayesian sampling of an explicit parametric model that includes both the cosmological CMB signal and non-cosmological astrophysical signals, such as thermal dust, CO, and low-frequency foregrounds. This framework is described in earlier papers [see Planck Collaboration XI (2016) and Planck Collaboration Int. XLVI (2016) and references therein for details]. The only changes since the 2015 implementation concern the data and model selection. As described Table 1. Baseline and optional late-time parameters, primordial power spectrum parameters, and slow-roll parameters. All primordial quantities are evaluated at a pivot scale of k * = 0.05 Mpc −1 , unless otherwise stated. Parameter Definition Fourth potential slow-roll parameter in Planck Collaboration IV (2018), we only use the Planck 2018 data in the current data release, whereas the previous 2015 version additionally included WMAP (Bennett et al. 2013) and Haslam (Haslam et al. 1982) observations. With fewer frequencies available, this requires a simpler model, and in particular we now fit for only a single low-frequency foreground component, rather than individual synchrotron, free-free, and spinning dust emission components, and we only fit a single CO component, rather than for individual CO line components at 100, 217, and 353 GHz. On the one hand, this results in fewer internal foreground degeneracies compared to the 2015 version, and a likelihood that only depends on Planck data, but at the same time the simpler foreground modelling also requires a slightly larger Galactic mask. Overall, the two versions are very compatible in terms of the recovered CMB power spectra, as discussed in PPL18. For additional details on the Commander temperature analysis, see Planck Collaboration IV (2018).
The HFI low-polarization likelihood is based on the fullmission HFI 100-GHz and 143-GHz Q and U low-resolution maps cleaned through a template-fitting procedure with LFI 30 GHz and HFI 353 GHz information 4 used as tracers of polarized synchrotron and thermal dust, respectively (see PPL18 for details about the cleaning procedure). The likelihood method, called SimAll, represents a follow-up of the SimBaL algorithm presented in Planck Collaboration Int. XLVI (2016) and uses the FFP10 simulations to construct an empirical probability for the EE and BB spectra. The method is based on the quadratic maximum likelihood estimation of the cross-spectrum between 100 and 143 GHz, and its multipole range spans from = 2 to = 29. We do not build a likelihood for T E, given the poor statistical 4 The polarized synchrotron component is fitted only at 100 GHz, being negligible at 143 GHz. For the polarized dust component, following the prescription in Planck Collaboration III (2018), the low-HFI polarization likelihood used the 353-GHz map constructed only from polarization-sensitive bolometers. consistency of the T E spectrum, as discussed in PPL18. Further details about the method and validation tests are presented in PPL18. When combined with the low-temperature likelihood (based on the Commander CMB solution), the low-polarization likelihood implies τ = 0.050 ± 0.009 and r < 0.36 at 95 % CL.
As an alternative to the Commander and SimAll low-likelihood, an update of the joint temperature and polarization pixelbased low-LFI likelihood used in 2015 is part of this Planck data release. Its methodology (see PPL18 for details) is similar to that of 2015 (Planck Collaboration XI 2016), i.e., a pixel-based approach in T QU at N side = 16, and employs the Commander solution in temperature along with the LFI 70-GHz linear polarization maps, foreground cleaned using the Planck 30-GHz and 353-GHz channels as tracer templates for synchrotron and dust, respectively. This 2018 version allows for a larger sky fraction in polarization (66.4 %, compared to the previous 46 %) and retains the sky surveys 2 and 4 that were excluded in 2015. By performing a two-parameter estimate for A s and τ restricted to < 30, we find using this likelihood that τ = 0.063 ± 0.020 and ln(10 10 A s ) = 2.975±0.056 at 68 % CL. The latter value has been derived by varying the T T , EE, and T E CMB spectra.

Planck high-likelihood
The 2018 baseline high-likelihood (Plik) is an update of the 2015 baseline version. The CamSpec likelihood is also used to explore alternative data cuts and modelling of the data and is described below. Both approaches implement a Gaussian likelihood approximation using cross-spectra between the 100-, 143-, and 217-GHz maps. Plik covers the multipoles 30 ≤ ≤ 2509 in temperature and 30 ≤ ≤ 1997 in polarization (i.e., for T E and EE). In order to avoid noise bias, the high-likelihood relies only on half-mission map cross-spectra, which have been demonstrated to be largely free of correlated noise. The spectra are computed on masked maps in order to reduce the anisotropic 6 Planck Collaboration: Constraints on Inflation Galactic contamination (dominated by dust emission), and in the case of T T also strong point sources and CO emission. The Plik masks, identical to the 2015 masks, are tailored to each frequency channel and differ in temperature and polarization to take into account differing foreground behaviour and channel beams. The Plik intensity (polarization) masks effectively retain 66, 57, and 47 % (70, 50, and 41 %) of the sky after apodization for the 100, 143, and 217 GHz channels, respectively. Unlike in 2015, when the map beams were computed for an average sky fraction, they are now computed for the exact sky fraction used at each frequency. The data vector used in the likelihood approximation discards multipoles that are highly contaminated by foregrounds or have low signal-to-noise ratios.
The Plik power spectra are binned using the same scheme as in 2015. Unbinned likelihoods are also available. When forming the data vector, individual cross-frequency spectra are not co-added. This allows for independent exploration of the calibration, nuisance, and foreground parameter space for each crossspectrum using dedicated templates in the theory vector.
The Plik (and CamSpec) covariance matrices are computed for a fixed fiducial CMB including the latest estimate of the foreground and systematics, which are all assumed to be Gaussian. As verified in 2015, after the masks have been applied this is a reasonable assumption. The covariance matrix computation uses an approximation to account for mask-induced correlations. Plik uses only the large Galactic mask in the analytic computation and then takes extra correlations due to the point-source mask into account using a Monte Carlo estimate of the extra variance induced. Missing pixels are ignored in the covariance. In 2015 its was shown that this approach (i.e., Gaussian approximation and approximate covariance) induced only a less than 0.1σ bias on n s (from the 30 < < 100 modes).
The Plik noise model has been re-estimated on the latest HFI maps using the same methodology as in 2015, based on a comparison between noise-biased auto-spectra and crossspectra. This procedure avoids correlated glitch residuals, which had biased previous noise estimates (Planck Collaboration XI 2016), particularly in polarization at 500. The 2018 HFI data processing pipeline has refined the maps used in the likelihood relative to 2015. For example, an improved destriping procedure reduced the residual scatter in the polarization maps, in particular at 143 GHz (yielding about 12 % lower noise on the half-mission cross-spectrum). More stringent selection cuts resulted in the discarding of the last 1000 rings of data, increasing the noise in the temperature half-mission spectrum by about 3 %. Also, a higher threshold was imposed on the conditioning of the T QU intra-pixel noise covariance matrix for a pixel to be considered well-measured, resulting in more missing pixels relative to 2015.
The data modelling has also significantly improved, in particular for polarization, making cosmological constraints from polarization more reliable. In 2015, the polarization spectra (T E and EE) displayed relatively large inter-frequency disagreements. A plausible explanation (at least for T E) was the temperature-to-polarization leakage induced by beam and calibration differences (so-called "beam leakage"). The beam-leakage modelling has improved substantially in 2018 (Hivon et al. 2017) so that we can now propagate the beams, gain differences, polarization angles, etc. to compute a reliable template for the beam leakage and thus remove these leakage effects. These improvements substantially reduce the T E interfrequency disagreements.
We also reassessed the estimates of the polarization efficiency for the polarized channels. Comparing different data-based estimates demonstrates that the ground-based polarization efficiency uncertainty estimates (of the order of a fraction of a percent) were too optimistic by a factor of 5 to 10. Correcting for the observed polarization efficiency errors (at the percent level) very significantly reduces the EE inter-frequency disagreements. This calibration correction relies on cosmological priors (using the T T best-fit cosmology). Calibrating using either the T E or the EE spectra yields generally consistent results, except at 143 Ghz where there is disagreement at more than 2σ. At this level, this discrepancy can be caused either by a statistical fluctuation, or by an unknown residual.
The Plik baseline likelihood implements a map-based calibration. The T E calibration is deduced from the T T and EE calibrations, including at 143 GHz. Other improvements over the 2015 version are the following. The dust model has been improved in temperature and polarization, using also the latest version of the 353-GHz maps. The level of synchrotron contamination in the 100-GHz and 143-GHz maps has been estimated and shown to be negligible. Sub-pixel noise has been included in T T and EE (and demonstrated to have a negligible effect on the cosmological parameters). Finally, a correlated component of the noise has been observed in the end-to-end HFI simulation, affecting the large scales and very small scales of the EE auto-frequency spectra. The large-scale contribution affects the dust correction and the n s constraints. We constructed an empirical model of this correlated noise from our simulations, which is included in the Plik likelihood.
CamSpec was the baseline for the 2013 release and was described in detail in Planck Collaboration XV (2014), and used cross-spectra formed from detector-set temperature maps using data from the nominal mission period. It was extended for the 2015 release to include both polarization and temperaturepolarization cross-spectra and to use the data from the full mission period. Similarly to Plik, CamSpec switched from detectorset cross-spectra to cross-spectra formed from frequency maps constructed from separate halves of the full mission data in order to mitigate the effects of noise correlated between detectors. In 2015, the foreground modelling was also modified and the sky fraction retained at each frequency was increased, using common masks with Plik in temperature. CamSpec used a more conservative mask in polarization than Plik.
Differently from Plik, CamSpec corrects each T E and EE cross-frequency spectrum with a fixed dust and temperature-topolarization leakage template before co-adding them to form the EE and T E components of its data vector and bases its noise estimate on differences between maps constructed using alternating pointing periods. Note also that CamSpec uses an individualspectrum-based calibration scheme, where the T E calibrations are not fixed to be those inferred from the T T and EE ones.
In the 2018 release further improvements in the CamSpec foreground modelling have been implemented. The dust model in temperature has been updated in a way similar to Plik. CamSpec now uses a richer model of the cosmic infrared background, allowing for the exploration of any impact on the cosmological parameters. As explained above, the noise modelling was also modified. Further modifications of the masking have been made for polarization, still using the same masks for each frequency channel, but different from the Plik mask. As we discussed above, beam-leakage and polarization-efficiency corrections are applied to the individual polarization spectra before their addition for inclusion in the likelihood. More details on the Plik and CamSpec likelihoods can be found in PPL18.
As in 2015, the high-Plik and CamSpec likelihoods are in excellent agreement for temperature. The different assumptions for the polarization-efficiency parameters and the masks (Planck Collaboration XI 2016) propagate into differences in cosmological parameter estimates. For the baseline ΛCDM model, the difference in cosmological parameters between the Plik likelihood and the CamSpec likelihood (using joint TT,TE,EE in combination with Commander, SimAll, and lensing) is at most 0.5σ (for Ω b h 2 ) (Planck Collaboration XIII 2016, henceforth PCP15). Similar differences in cosmological parameters occur in extended cosmological models. The differences between the Plik and CamSpec parameters are dominated by calibration-model differences for the joint TT,TE,EE and TEonly cases, and by mask differences for the EE-only case. To a large extent, the CamSpec results can be reproduced within the Plik framework simply by changing in Plik the calibration model (for T E) and the polarization mask (for EE). Below we use Plik as the Planck baseline high-likelihood. CamSpec results are used to assess the residual uncertainty from modelling and mask choices. We quote values obtained with CamSpec only for a few cases.
We use the following conventions for naming the Planck likelihoods: (i) Planck TT+lowE denotes the combination of the high-TT likelihood at multipoles ≥ 30 and the lowtemperature-only Commander likelihood, plus the SimAll low-EE-only likelihood in the range 2 ≤ ≤ 29; (ii) Planck TE and Planck EE denote the TE and EE likelihood at ≥ 30, respectively; (iii) Planck TT,TE,EE+lowE denotes the combination of the combined likelihood using T T , T E, and EE spectra at ≥ 30, the low-temperature Commander likelihood, and the low-SimAll EE likelihood; and (iv) Planck TT,TE,EE+lowP denotes the combination of the likelihood using T T , T E, and EE spectra at ≥ 30 and the alternative joint temperaturepolarization likelihood at 2 ≤ ≤ 29, based on the temperature Commander map and the 70-GHz polarization map. Unless otherwise stated, high-results are based on the Plik likelihood and low-polarization information is based on SimAll.

Planck CMB lensing likelihood
The Planck 2018 lensing likelihood, presented in PPLe18, uses the lensing trispectrum to estimate the power spectrum of the lensing potential C φφ L . This signal is extracted using a minimum-variance combination of a full set of temperature-and polarization-based quadratic lensing estimators (Okamoto & Hu 2003) applied to the SMICA CMB map over approximately 70 % of the sky using CMB multipoles 100 ≤ ≤ 2048, as described in PPLe18. We use the lensing bandpower likelihood, with bins spanning lensing multipoles 8 ≤ L ≤ 400, which has been validated with numerous consistency tests. Because its multipole range has been extended down to L = 8 (compared to L = 40 for the Planck 2015 analysis), the statistical power of the lensing likelihood used here is slightly greater.

Non-Planck data
While the data derived exclusively from Planck observations are by themselves already extremely powerful at constraining cosmology, external data sets can provide helpful additional information. The question of consistency between Planck and external data sets is discussed in detail in Planck Collaboration VI (2018) (henceforth PCP18). Here we focus on two data sets that are particularly useful for breaking degeneracies and whose errors can be assessed reliably. We consider the measurement of the CMB B-mode polarization angular power spectrum by the BICEP2/Keck Array collaboration and measurements of the baryon acoustic oscillation (BAO) scale. The supplementary Bmode data provide independent constraints on the tensor sector, which are better than those that can be derived from the Planck data alone (based on the shape of the scalar power spectrum). The BAO data, on the other hand, do not directly constrain the primordial perturbations. These data, however, provide invaluable low-redshift information that better constrain the late-time cosmology, especially in extensions of ΛCDM, and thus allow degeneracies to be broken.

BICEP2/Keck Array 2014 B-mode polarization data
Although Planck measured the CMB polarization over the full sky, its polarization sensitivity in the cosmological frequency channels is not sufficient to compete with current suborbital experiments surveying small, particularly low-foreground patches of the sky very deeply using many detectors. In PCI15, constraints on r using the joint BICEP2/Keck Array and Planck (BKP) analysis (BICEP2/Keck Array and Planck Collaborations 2015) were reported. Here we make use of the more recent B-mode polarization data available from the analysis of the BICEP2/Keck field (Ade et al. 2016, henceforth BK14). The BK14 likelihood draws on data from the new Keck array at 95 GHz in addition to the 150-GHz channel, as well as from Planck and WMAP data to remove foreground contamination. The BK14 observations measure B-mode polarization using 11 auto-and 55 cross-spectra between the BICEP2/Keck maps at 95 and 150 GHz, the WMAP maps at 23 and 33 GHz, and the Planck maps at 30, 44, 70, 100, 143, 217, and 353 GHz, using nine bins in multipole number. By using B-mode information only within the BK14 likelihood, a 95 % upper limit of r < 0.09 is found (BK14). An improved likelihood is imminent (Ade et al. 2018, BK15) including more data and most notably the addition of a new 220-GHz channel; however, this version of the paper does not include results with the BK15 likelihood.

Baryon acoustic oscillations
Acoustic oscillations of the baryon-photon fluid prior to recombination are responsible for the acoustic peak structure of the CMB angular power spectra. The counterpart to the CMB acoustic peaks in the baryon distribution are the BAOs, which remain imprinted into the matter distribution to this day. In the position-space picture, the BAOs of the power spectrum correspond to a peak in the correlation function, defining a characteristic, cosmology-dependent length scale that serves as a standard ruler and can be extracted (e.g., from galaxy redshift surveys). The transverse information of a survey constrains the ratio of the comoving angular diameter distance and the sound horizon at the drag epoch (i.e., when the baryon evolution becomes unaffected by coupling to the photons), D M /r d , whereas the line-of-sight information yields a measurement of H(z)r d . Sometimes, these two observables are combined to form the direction-averaged For our BAO data compilation, we use the measurements of D V /r d from the 6dF survey at an effective redshift z eff = 0.106 (Beutler et al. 2011), and the SDSS Main Galaxy Sample at z eff = 0.15 (Ross et al. 2015), plus the final interpreta-8 Planck Collaboration: Constraints on Inflation tion of the SDSS III DR12 data (Alam et al. 2016), with separate constraints on H(z)r d and D M /r d in three correlated redshift bins at z eff = 0.38, 0.51, and 0.61. In Addison et al. (2018), the same set of BAO data combined with either non-Planck CMB data or measurements of the primordial deuterium fraction was shown to favour a cosmology fully consistent with, but independent of, Planck data.

Planck 2018 results for the main inflationary observables
As in PCI13 and PCI15, we start by describing Planck measurements of the key inflationary parameters. Some of the results reported in this section can be found in the Planck Legacy Archive. 5

Results for the scalar spectral index
Planck temperature data in combination with the EE measurement at low multipoles determine the scalar spectral tilt in the ΛCDM model as n s = 0.9626 ± 0.0057 (68 % CL, Planck TT+lowE) . (6) This result for n s is compatible with the Planck 2015 68 % CL value n s = 0.9655 ± 0.0062 for Planck TT+lowP (PCP15). The slightly lower value for n s is mainly driven by a corresponding shift in the average optical depth τ, now determined as which is to be compared with the Planck 2015 value τ = 0.078 ± 0.022 (PCP15). This more precise determination of τ is due to better noise sensitivity of the HFI 100-and 143-GHz channels employed in the low-SimAll polarization likelihood, compared to the joint temperature-polarization likelihood based on the LFI 70-GHz channel in 2015. Because of the degeneracy between the average optical depth and the amplitude of the primordial power spectrum, A s and σ 8 are also lower than in the Planck 2015 release. These shifts from the Planck 2015 values for the cosmological parameters have been anticipated with the first results from the HFI largeangular polarization pattern (Planck Collaboration Int. XLVI 2016; Planck Collaboration Int. XLVII 2016). The trend toward smaller values for (n s ,τ) with respect to the Planck 2015 release also occurs for different choices for the lowlikelihood. By substituting Commander and SimAll with the updated joint temperature-polarization pixel likelihood coming from the LFI 70-GHz channel, we obtain in combination with high-temperature data: Although with larger errors, these latter results are consistent with the shifts induced by a determination of a lower optical depth than in the Planck 2015 release, 6 as found in Eqs. (6) and (7). Given this broad agreement and the consistency in the values of τ derived with SimAll separately from the three cross-spectra 70 × 100, 70 × 143, and 100 × 143 (PPL18), we will mainly use the baseline low-likelihood in the rest of the paper. As anticipated in 2015, the information in the high-polarization Planck data is powerful for breaking degeneracies in the parameters and to further decrease parameter uncertainties compared to temperature data alone. The addition of high-polarization leads to a tighter constraint on n s : The shift in n s (and, more generally, in the cosmological parameters of the base-ΛCDM model) obtained when Planck lensing is combined with TT,TE,EE+lowE is smaller than in 2015 because of the improved polarization likelihoods. The combination with lensing is, however, powerful for breaking parameter degeneracies in extended cosmological models, and, therefore, for this 2018 release we will consider the full information contained in temperature, polarization, and lensing, i.e., TT,TE,EE+lowE+lensing, as the baseline Planck data set. Figure 3 shows a comparison of the Planck 2018 baseline results with those from alternative likelihoods and from the 2015 baseline for the ΛCDM cosmological parameters. As in 2013 and 2015, BAO measurements from galaxy surveys are consistent with Planck. When BAO data are combined, we obtain for the base-ΛCDM cosmology: n s = 0.9665 ± 0.0038 (14) (68 % CL, Planck TT,TE,EE+lowE+lensing+BAO) .
The combination with BAO data decreases (increases) the marginalized value of Ω c h 2 (Ω b h 2 ) obtained by Planck, and this effect is compensated for by a shift in n s towards slightly larger values.  Table 2. Confidence limits for the cosmological parameters in the base-ΛCDM model from Planck temperature, polarization, and temperature-polarization cross-correlation separately and combined, in combination with the EE measurement at low multipoles.
3.2. Ruling out n s = 1 One of the main findings drawn from previous Planck releases was that the scale-independent Harrison-Zeldovich (HZ) spectrum (Harrison 1970;Zeldovich 1972;) is decisively ruled out. This conclusion is reinforced in this release: in standard ΛCDM late-time cosmology, the scalar spectral index from Table 2  Simple one-parameter modifications of the cosmological model are not sufficient to reconcile a scale-invariant power spectrum with Planck data. For instance, when the effective number of neutrino species N eff is allowed to float for a cosmology with a scale-invariant spectrum, the effective ∆χ 2 with respect to the power-law spectrum are ∆χ 2 = 12.9, 27.5, and 30.2, respectively.
When instead the assumption of flat spatial sections is relaxed, 7 we obtain ∆χ 2 = 11. 8, 28.8, and 40.9, respectively, for the same data sets. Therefore, the corresponding closed cosmological models fitting Planck TT+lowE The Planck 2018 data are consistent with a vanishing running of the scalar spectral index. Using Planck 7 For non-flat models the power spectra encode the eigenvalues of the corresponding Laplacian operator of the spatial sections and scale invariance holds at scales much smaller than the curvature radius. 8 This is not a new result based on the Planck 2018 release, but just an update of a similar conclusion also reached with the Planck 2015 data. Compared to the flat ΛCDM tilted model, we obtain ∆χ 2 = 12. 3, 34.8, and 45 with Planck 2015 TT+lowP, Planck 2015 TT,TE,EE+ lowP, and Planck 2015 TT,TE,EE+lowP+lensing, respectively. Therefore, even with Planck 2015 data, a closed model with n s = 1 provides a worse fit than tilted ΛCDM and is not compelling as claimed in Ooba et al. (2017).
TT,TE,EE+lowE+lensing we obtain These results are consistent with, and improve on, the Planck 2015 result, dn s /d ln k = −0.008 ± 0.008 (PCP15). As discussed in PCI13 and PCI15, a better fit to the temperature low-deficit was found in 2015, thanks to a combination of non-negative values for the running and the running of the running. The Planck 2018 release has significantly reduced the parameter volume of this extension of the base-ΛCDM model. The Planck 2018 TT(TT,TE,EE)+lowE+lensing constraints for the model including running of running are all at 68 % CL. It is interesting to note that the high-temperature data still allow a sizable value for the running of the running, although slightly decreased with respect to the Planck 2015 results (PCI15). However, when high-Planck 2018 polarization data are also included, dn s /d ln k and d 2 n s /d ln k 2 are tightly constrained.
The model including a scale-dependent running can produce a better fit to the low-deficit at the cost of an increase of power at small scales; this latter effect is constrained in this release. As an example of a model with suppression only on large scales, we also reconsider the phenomenological model with an exponential cutoff: which can be motivated by a short stage of inflation (Contaldi et al. 2003;Cline et al. 2003) (see also Kuhnel & Schwarz (2010), Hazra et al. (2014), and Gruppuso et al. (2016) for other types of large-scale suppressions). We do not find any statistically significant detection of k c using either logarithmic or linear priors and for different values of λ c , with any combination of Planck baseline likelihoods. We have also checked that these results depend only weakly on the exclusion of the EE quadrupole in SimAll and are stable to the substitution of Commander and SimAll with the joint temperature-polarization likelihood based on the 70-GHz channel.

Constraints on spatial curvature
Since the vast majority of inflation models predict that the Universe has been driven towards spatial flatness, constraints on the spatial curvature provide an important test of the standard scenario. Therefore in this subsection we extend the base-ΛCDM model with the addition of the spatial curvature parameter, Ω K . For the case of Planck TT,TE,EE+lowE+lensing, we find a constraint of Ω K = −0.011 +0.013 −0.012 (95 % CL) .
The inclusion of Planck lensing information only weakly breaks the geometrical degeneracy (Efstathiou & Bond 1999) which results in the same primary fluctuations while varying the total matter density parameter, Λ, and H 0 . The degeneracy can be effectively broken with the addition of BAO data, in which case Planck TT,TE,EE+lowE+lensing+BAO gives Although Ω K is one of the cosmological parameters exhibiting some differences between Plik and Camspec, the constraints in Eqs. (20) and (21) are quite robust due to the inclusion of lensing (and BAO) information.
A constraint on the curvature parameter can be translated into a constraint on the radius of curvature, R K , via both at 95 % confidence. These lengths are considerably greater than our current (post-inflation) particle horizon, at 13.9 Gpc.
Our tightest constraint, Eq. (21), tells us that our observations are consistent with spatial flatness, with a precision of about 0.4 %. However, even if inflation has driven the background curvature extremely close to zero, the presence of fluctuations implies a fundamental "cosmic variance" for measurements of curvature confined to our observable volume. In particular, the known amplitude of fluctuations implies a standard deviation for Ω K of roughly 2 × 10 −5 (Waterhouse & Zibin 2008). Therefore our best constraint is still a factor of roughly 10 2 above the cosmic variance limit for a flat universe. A future measurement of negative curvature above the cosmic variance floor would point to open inflation (Gott 1982;Gott & Statler 1984;Bucher et al. 1995;Yamamoto et al. 1995;Ratra & Peebles 1995;Lyth & Stewart 1990), while a measurement of positive curvature could pose a problem for the inflationary paradigm due to the difficulty of producing closed inflationary models (Kleban & Schillo 2012).
Alternatively, excess spatial curvature might be evidence for the intriguing possibility that there was "just enough" inflation to produce structure on the largest observable scales. Indeed an upper limit on spatial curvature implies a lower limit on the total number of e-folds of inflation (see, e.g., Komatsu et al. 2009). We can relate these limits to the number of e-folds of inflation, N * = N(k * ), after scale k * left the Hubble radius during inflation, to be given explicitly in Eq. (47). We define the (constant) curvature scale, k K , as the inverse of the comoving radius of curvature, i.e., In the absence of special initial conditions, inflation will begin with a curvature parameter of order unity. Equation (25) then implies that k K ∼ aH at the start of inflation, i.e., the curvature scale is "exiting the horizon" at that time. Then the lower limit on the number of e-folds of inflation will simply be N K ≡ N(k K ), i.e., the number of e-folds after scale k K left the Hubble radius during inflation. With Eq. (47) this gives 9 With the pivot scale of k * = 0.002 Mpc −1 (for comparison with the values in Sect. 4.2) and our tightest upper limit on Ω K from Eq. (21), this becomes That is, our constraint on spatial curvature implies that inflation must have lasted at least about 5 e-folds longer than required to produce the pivot scale k * . Equation (27) provides a modelindependent comparison between the e-folds required to solve the flatness problem (to current precision) and to produce largescale fluctuations (at scale k * ). We stress that this limit assumes a unity curvature parameter at the start of inflation (although the dependence on this assumption, being logarithmic, is weak).
For comparison with the result of Komatsu et al. (2009), we can simplify to the case of instantaneous thermalization and constant energy density during inflation. Then we find where T end is the reheating temperature.

Constraints on the tensor-to-scalar ratio
This subsection updates constraints on the tensor-to-scalar ratio r assuming that the tensor tilt satisfies the consistency relation, n t = −r/8, which is the case for slow-roll inflation driven by a single scalar field with a canonical kinetic term. By combining Planck temperature, low-polarization, and lensing we obtain r 0.002 < 0.10 (95 % CL, Planck TT+lowE+lensing) . (29) This constraint slightly improves on the corresponding Planck 2015 95 % CL bound, i.e., r 0.002 < 0.11 (PCI15), and is unchanged when high-polarization data are also combined. Note that by using CAMspec instead of Plik as the high-joint temperature-polarization likelihood, we obtain a slightly looser bound, i.e., r 0.002 < 0.14 at 95 % CL. By including the Planck B-mode information at 2 < < 30 in the low-polarization likelihood, the 95 % CL constraint is essentially unchanged.
Since inflationary gravitational waves contribute to CMB temperature anisotropies mostly at 100, the low-temperature deficit contributes in a nontrivial way to the Planck bound on r. By excising the 2 ≤ ≤ 29 temperature data, the constraint on r with Planck TT,TE,EE+lensing+lowEB relaxes to r 0.002 < 0.16 (95 % CL) .
This result improves on the 2015 95 % CL result, i.e., r 0.24 (PCI15), because of the inclusion of high-polarization and of the improved determination of τ.
Since this Planck constraint on r relies on temperature and E-mode polarization, the Planck-only limit depends somewhat on the underlying cosmological model. Table 3 shows the constraints on n s and r for a few important extensions of ΛCDM plus tensors, which include a non-zero running, a non-zero spatial curvature, and a non-minimal neutrino sector. We observe that the bound on r is relaxed by at most 30 % when the scale dependence of the scalar tilt is allowed to vary. In all the other extensions the Planck r bound is modified at most by 10 %, demonstrating the constraining power of the Planck 2018 release in reducing the degeneracy of the tensor-to-scalar ratio with other cosmological parameters. As far as the scalar tilt is concerned, we find the largest shift (by roughly 1σ higher) when the assumption of spatial flatness is relaxed.
A B-mode polarization measurement can further tighten the constraint on r and help in reducing its degeneracies with other cosmological parameters that may appear when using only temperature and E-mode polarization data. After the release of the first BICEP-Keck Array-Planck (BKP) joint cross-correlation, constraints on r from B-mode polarization data alone have become tighter than those based on Planck data alone, thanks to the inclusion of the 95-GHz channel (BK14). By combining the Planck 2018 and BK14 data we obtain r 0.002 < 0.064 (95 % CL, Planck TT,TE,EE +lowE+lensing+BK14) .
Note that by using CAMspec instead of Plik as high-TT,TE,EE likelihood, we obtain a slightly looser bound, i.e., r 0.002 < 0.075 at 95 % CL. The effectiveness of the combination with the BKP likelihood in constraining r is also remarkable in the extensions of ΛCDM plus tensors, as can be seen from 3.6. Beyond the tensor-to-scalar ratio consistency condition The increasing constraining power of B-mode polarization data allows us to set upper bounds on r without imposing the consistency condition for the tensor tilt, n t = −r/8, which is motivated by standard slow-roll single-scalar-field inflation. Deviations can occur in multifield inflation (Bartolo et al. 2001;Wands et al. 2002;Byrnes & Wands 2006), in the models with generalized Lagrangians (Garriga & Mukhanov 1999;Kobayashi et al. 2010), in gauge inflation (Maleknejad et al. 2013), or in a more radical way in alternative models to inflation (Gasperini & Veneziano 1993;Boyle et al. 2004;Brandenberger et al. 2007).
As the current data do not lead to a detection of a non-zero tensor amplitude, virtually any value of n t would give a good fit as long as r is close enough to zero. Therefore, as in PCI15, we characterize the tensor perturbations by two well-constrained parameters that we choose to be r at two different scales, (r k 1 , r k 2 ), with k 1 = 0.002 Mpc −1 and k 2 = 0.02 Mpc −1 , and assume a power-law power spectrum. We call this two-parameter extension of the ΛCDM model the "ΛCDM+r 0.002 +r 0.02 " model. We also quote our results in terms of (r˜k, n t ), calculated from the primary parameters as n t = [ln(r k 2 /r k 1 ) / ln(k 2 /k 1 )] + n s − 1 and r˜k = A t1 (k/k 1 ) n t , where A t1 = r k 1 A s (k 1 /k * ) n s −1 is the tensor perturbation amplitude at k 1 and k * = 0.05 Mpc −1 is the pivot scale used for scalar perturbations. Fork we choose 0.01 Mpc −1 , which corresponds roughly to the decorrelation scale of r and n t when using the Planck and BK14 data.
The constraints on the derived tensor parameters are r 0.01 < 0.091 and −0.34 < n t < 2.63 at 95 % CL. The left and right panels of Fig. 5 show the two-dimensional contours for the primary parameters (r 0.002 , r 0.02 ) and the derived ones (r 0.01 , n t ), respectively. The consistency condition, n t = −r/8, denoted by the dashed lines, is fully compatible with Planck+BK14 data. However, a very blue tensor tilt with n t 2 and r 0.01 0.05 is still within the 68 % CL region. Indeed, despite the larger amplitude of the primordial tensor power spectrum at small scales for blue n t , the tensor modes are suppressed when re-entering the Hubble radius, which leads to damping of the observational signal at high k. This explains why the 68 % CL CMB constraint on r 0.02 is about an order of magnitude weaker than the one on r 0.002 .
A stochastic background of gravitational waves (GWs) with a blue tensor tilt could be further constrained at much shorter wavelengths, such as those probed by ground-based interferometers dedicated to the direct detection of GWs. For example, assuming a scale-invariant tensor spectrum and using the frequency range (20-85.8) Hz, which corresponds to the wavenumbers k = 2π f = (1.3-5.5) × 10 16 Mpc −1 , LIGO and Virgo set an upper bound on the GW density parameter of Ω GW ( f ) ≤ 1.7 × 10 −7 at 95 % CL (Abbott et al. 2017). While these scales are likely to be dominated by astrophysical sources, such as GWs from binary mergers, we next examine what constraints the LIGO&Virgo upper bound sets on primordial tensor perturbations, if we assume that they had a power-law spectrum all the way from CMB scales to ultra-short scales. We refer the interested reader to  and Cabass et al. (2016) for the use of alternative data on short scales or for additional assumptions on the energy-momentum tensor of the stochastic background of GWs averaged over wavelengths. We obtain a conservative upper limit on the primordial contribution by demanding that the GW density from our scaledependent primordial tensor perturbations Abbott et al. 2017;Cabass et al. 2016), stays below the above-quoted limit at least at k = 1.3×10 16 Mpc −1 ( f = 20 Hz). The posterior probability densities when this constraint is included in the analysis as a half-Gaussian prior are compared with those obtained by Planck+BK14 alone in Figs. 4 and 5. LIGO&Virgo sets a very high upper bound 10 on r at ultrahigh k, separated from CMB scales by a factor of 10 18 in k. Due . 68 % and 95 % CL constraints on tensor perturbations in the ΛCDM+r 0.002 +r 0.02 model, i.e., when the inflationary consistency relation is relaxed. Filled contours in the left panel show the results for our independent primary parameters r 0.002 and r 0.02 , which have uniform priors, and in the right panel for the derived parameters n t and r 0.01 , which have non-uniform priors. The dotted lines assume uniform priors on r 0.01 and n t , calculated as in Fig. 4. The scale k = 0.01 Mpc −1 is near the decorrelation scale of (n t , r) for the Planck+BK14 data. In both panels the dashed black line indicates the inflationary consistency condition, n t = −r 0.01 /8.
to the long arm length, this effectively provides a cutoff for n t and excludes the bluest spectra that were allowed by the CMB alone, leading to r 0.002 < 0.069 r 0.02 < 0.099 or r 0.01 < 0.080 and −0.62 < n t < 0.53. The consistency condition n t = −r/8 is also compatible with these tighter constraints, as can be seen by comparing the red contours and dashed black lines in Fig. 5. As LIGO&Virgo pushes r 0.02 down (and we assume a power-law tensor spectrum), the upper bound on r 0.002 becomes weaker than with the CMB alone. This is not surprising, since the system is analogous to a see-saw with a pivot point at k ∼ 0.01 Mpc −1 , where the data are the most sensitive to the tensor perturbations (taking into account also the transfer function from primordial tensor perturbations to the observable Bmode signal). Once one end of the see-saw is pushed down the other end can go up without disturbing the spectrum too much at the middle scales. We will observe analogous behaviour with isocurvature perturbations, for which we also assume a powerlaw spectrum and have only an upper bound (not a detection); see Sect. 9.3.

Implications for single-field slow-roll inflationary models
In this section we discuss the implications of the Planck 2018 likelihood for standard single-field slow-roll inflation. We first update the results for the Hubble flow functions (HFFs) i and the potential slow-roll parameters obtained by the analytic perturbative expansion in terms of the HFFs for the primordial spectra of fluctuations. For definitions of the HFF hierarchy and the potential slow-roll parameters see Table 1. We then present a Bayesian comparison for a representative selection of standard slow-roll inflationary models.

Constraints on slow-roll parameters
Exploiting the approximate analytic expressions for the primordial power spectrum of scalar and tensor fluctuations obtained by the Green's function method (Stewart & Lyth 1993;Gong & Stewart 2001;Leach et al. 2002), we can construct constraints on the slow-roll parameters. When restricting to parameters first order in the HFFs, we obtain with Planck TT,TE,EE+lowE+lensing(+BK14)  As can be seen from Fig. 6, the 95 % CL allowed contours are mostly in the region of concave potentials when BK14 is combined with Planck 2018 data. When contributions to the primordial power spectra that are second-order in the HFFs are included, for Planck TT,TE,EE+lowE+lensing(+BK14) we obtain the following constraints on the slow-roll HFFs: and on the slow-roll potential parameters we obtain: The marginalized 68 % and 95 % CLs for the slow-roll HFF and potential parameters, allowing 3 0, with Planck data alone or in combination with BK14, are displayed in Fig. 7.
Parameter range Prior type Table 4. Priors for cosmological parameters used in the Bayesian comparison of inflationary models.

Implications for selected slow-roll inflationary models
The predictions for (n s , r) to first order in the slow-roll approximation for a few inflationary models are shown in Fig. 8, which updates figure 12 of PCI15 and figure 1 of PCI13 with the same notation. These predictions are calculated for scale k = 0.002 Mpc −1 and include an uncertainty in the number of e-folds of 50 < N * < 60.
In the following we discuss the implications of the Planck 2018 data release for selected slow-roll inflationary models by taking into account the uncertainties in the entropy generation stage. As in PCI15, we use the primordial power spectra of cosmological fluctuations generated during slow-roll inflation parameterized by the HFFs, i , to second order, which can be expressed in terms of the parameters of the inflationary model and the number of e-folds to the end of inflation, N * (Liddle & Leach 2003;Martin & Ringeval 2010), given by (PCI13) where ρ end is the energy density at the end of inflation, a 0 H 0 is the present Hubble scale, V * is the potential energy when k * left the Hubble radius during inflation, w int characterizes the effective equation of state between the end of inflation and the thermalization energy scale ρ th , and g th is the number of effective bosonic degrees of freedom at the energy scale ρ th . We fix g th = 10 3 and end = 1, and we use modified routines of the public code ASPIC 11 (Martin et al. 2014b). In order to make contact with Fig. 8, we consider the pivot scale k * = 0.002 Mpc −1 in this subsection. We assume the uniform priors for the cosmological parameters listed in Table 4, and logarithmic priors on 10 10 A s (over the interval [e 2.5 , e 3.7 ]) and ρ th (over the interval [(1 TeV) 4 , ρ end ]). Prior ranges for additional parameters in the inflationary models considered are listed in Table 5. In this paper we restrict ourselves to the case w int = (p − 2)/(p + 2), when the potential can be approximated as V(φ) ∝ φ p during the coherent oscillation regime after inflation, or simply w int = 0 when the potential considered describes only the inflationary stage. For data we use the full constraining power of Planck, i.e., Planck TT,TE,EE+lowE+lensing, in combination with BK14. The ∆χ 2 and the Bayesian evidence values for a selection of inflationary models with respect to the R 2 model (Starobinsky 1980;Mukhanov & Chibisov 1981;Starobinsky 1983) are shown in Table 5. Figure 9 shows the resulting marginalized probability densities of n s and r at k = 0.002 Mpc −1 for a few inflationary models with the above specified priors, compared to the corresponding 68 % and 95 % CL  Marginalized joint two-dimensional 68 % and 95 % CL regions for combinations of ( 1 , 2 , 3 ) (upper panels) and ( V , η V , ξ 2 V ) (lower panels) for Planck TT,TE,EE+lowE+lensing (red contours), compared with Planck TT,TE,EE+lowE+lensing+BK14 (blue contours). Table 5. Bayesian comparison for a selection of slow-roll inflationary models with w int fixed (see text for more details). We quote 0.3 as the error on the Bayes factor. Models are strongly disfavoured when ln B < −5.

Power-law potential
Marginalized joint 68 % and 95 % CL regions for n s and r at k = 0.002 Mpc −1 from Planck alone and in combination with BK14 or BK14 plus BAO data, compared to the theoretical predictions of selected inflationary models. Note that the marginalized joint 68 % and 95 % CL regions assume dn s /d ln k = 0.
limits obtained from a ΛCDM-plus-tensor fit. We refer the interested reader to PCI15 for a concise description of the inflationary models studied here and we limit ourselves here to a summary of the main results of this analysis.
-The inflationary predictions (Mukhanov & Chibisov 1981;Starobinsky 1983) originally computed for the R 2 model (Starobinsky 1980) to lowest order, are in good agreement with Planck 2018 data, confirming the previous 2013 and 2015 results. The 95 % CL allowed range 49 < N * < 58 is compatible with the R 2 basic predictions N * = 54, corresponding to T reh ∼ 10 9 GeV (Bezrukov & Gorbunov 2012). A higher reheating temperature T reh ∼ 10 13 GeV, as predicted in Higgs inflation (Bezrukov & Shaposhnikov 2008), is also compatible with the Planck data. -Monomial potentials (Linde 1983 Pl (φ/M Pl ) p with p ≥ 2 are strongly disfavoured with respect to the R 2 model. For these values the Bayesian evidence is worse than in 2015 because of the smaller level of tensor modes allowed by BK14. Models with p = 1 or p = 2/3 (Silverstein & Westphal 2008;McAllister et al. 2010McAllister et al. , 2014 are more compatible with the data. -There are several mechanisms which could lower the predictions for the tensor-to-scalar ratio for a given potential V(φ) in single-field inflationary models. Important examples are a subluminal inflaton speed of sound due to a nonstandard kinetic term (Garriga & Mukhanov 1999), a nonminimal coupling to gravity (Spokoiny 1984;Lucchin et al. 1986;Salopek et al. 1989;Fakir & Unruh 1990), or an additional damping term for the inflaton due to dissipation in other degrees of freedom, as in warm inflation (Berera 1995;Bastero-Gil et al. 2016). In the following we report on the constraints for a non-minimal coupling to gravity of the type To be more specific, a quartic potential, which would be excluded at high statistical significance for a minimally-coupled scalar inflaton as seen from Table 5, can be reconciled with Planck and BK14 data for ξ > 0: we obtain a 95 % CL lower limit log 10 ξ > −1.6 with ln B = −1.6.
-Natural inflation (Freese et al. 1990;Adams et al. 1993) is disfavoured by the Planck 2018 plus BK14 data with a Bayes factor ln B = −4.2.
-Within the class of hilltop inflationary models (Boubekeur & Lyth 2005) we find that a quartic potential provides a better fit than a quadratic one. In the quartic case we find the 95 % CL lower limit log 10 (µ 2 /M Pl ) > 1.1. -D-brane inflationary models (Kachru et al. 2003;Dvali et al. 2001;García-Bellido et al. 2002) provide a good fit to Planck and BK14 data for a large portion of their parameter space. -For the simple one parameter class of inflationary potentials with exponential tails (Goncharov & Linde 1984;Stewart 1995;Dvali & Tye 1999;Burgess et al. 2002;Cicoli et al. 2009) we find ln B = −1.0. -Planck 2018 data strongly disfavour the hybrid model driven by logarithmic quantum corrections in spontaneously broken supersymmetric (SUSY) theories (Dvali et al. 1994), with ln B = −5.0.
. Marginalized probability densities of the scalar tilt n s (top panel) and r (bottom panel) at k = 0.002 Mpc −1 for natural, R 2 , hilltop quartic, and V(φ) ∝ φ 2/3 inflation, obtained by marginalizing over the uncertainties in the entropy generation stage, compared to the corresponding 68 % and 95 % CL limits obtained from a ΛCDM-plus-tensor fit.
-Planck and BK14 data set tight constraints on α attractors Ferrara et al. 2013). We obtain log 10 α E 1 < 1.5 and log 10 α E 2 < 1.2 at 95 % CL for the Emodel. We obtain slightly tighter 95 % CL bounds for the T-model, i.e., log 10 α T 1 < 1.1 and log 10 α T 2 < 1.0. Given the relation |R K | = 2/(3α) between the curvature of the Kähler geometry R K and α in some of the T-models motivated by supergravity, Planck and BK14 data imply a lower bound on |R K |, which is still in the low-curvature regime. The discrete set of values α = i/3 with an integer i in the range [1,7] motivated by maximal supersymmetry (Ferrara & Kallosh 2016;Kallosh et al. 2017) is compatible with the current data.

Taylor expansion of V(φ) in the observable region
In this section, as in section 6 of PCI13 and section 7.1 of PCI15, we try to reconstruct the inflaton potential only in its observable window, making no assumptions about the end of inflation. The motivation for being so conservative is that what happens after the inflaton rolls down beyond this range might not be captured by the simplest descriptions. More elaborate treatments would be required, for instance, in the case of a non-trivial potential shape before the end of inflation, a waterfall transition involving extra scalar fields, or several short inflationary stages between the time at which CMB scales exit the Hubble radius and the nucleosynthesis epoch. The analysis of this section relies, however, on the assumption that the potential is smooth enough inside the observable window to be described by a Taylor expansion up to order four. Note that this assumption is much weaker than assuming that a Taylor expansion is valid up to the end of inflation. However, it excludes from the analysis potentials with sharp features in the observable window, such as those studied in the next sections.
We perform the Taylor expansion around the value φ * of the inflaton field evaluated precisely at the time t * when the pivot scale k * = 0.05 Mpc −1 fulfills the relation k * = a(t * )H(t * ). We separately study the cases where the expansion is performed at order n = 2, n = 3, or n = 4. We compute the primordial spectrum with a full integration of the Fourier mode evolution, using the inflationary module of the CLASS code. Although this method assumes no slow-roll approximation at any point, we speed up the convergence of the Markov Chain by taking flat priors not directly on the five Taylor coefficients {V, V φ , . . . , V φφφφ }, but on combinations of them matching the definitions of the potential slow-roll parameters { V , η V , ξ 2 V , 3 V } presented in table 2 of PCI13. Even beyond the slow-roll approximation, these combinations provide nearly linear contributions to the tilt, running, running of the running, etc., of the scalar and tensor spectrum. Hence, they are directly related to observable quantities and well constrained by the data. Instead, if we ran with flat priors on {V, V φ , . . . , V φφφφ }, the convergence would be plagued by complicated parameter degeneracies.
The results of this analysis are shown in the panels of Fig. 10 and Table 6 for n = 2, 3, and 4, using two data sets for each: Planck TT,TE,EE+lowE alone; or Planck TT,TE,EE+lowE+BAO+BK14. The plot in Fig. 10 deliberately has a lot of white space because, for the sake of comparison, we plotted it over the same parameter ranges as the same plot in PCI15. We notice some significant improvement. Comparing Planck TT,TE,EE+lowE results from 2015 and 2018, we find that error bars on individual parameters have typically been reduced by 30 % thanks to improved polarization data. Including BK data provides further constraining power. Comparing Planck 2015 TT+lowP+BAO and Planck 2018 TT,TE,EE+lowE+BAO+BK14, we find that error bars shrink by a factor of 2, or even a factor of 3 in the case n = 4. The new data tend to resolve degeneracies which previously appeared in the n = 4 case and could be understood as a compensation mechanism between potentials with large running of the tilt, running of the running, tensor contribution, etc. The param-  Table 6. Numerical reconstruction of the potential slowroll parameters beyond any slow-roll approximation, when the potential is Taylor-expanded to nth order, using Planck TT,TE,EE+lowE+BAO+BK14. We also show the corresponding bounds on some related parameters (here n s , dn s /d ln k, and r 0.002 are derived from the numerically computed primordial spectra). All error bars are at the 68 % CL and all upper bounds at the 95 % CL. The effective χ 2 value of model n is given relative to model n − 1.
eters ξ 2 V and 3 V are perfectly compatible with zero (see Fig. 10 and Table 6), and so are V φφφ and V φφφφ (see the contours on the parameters {V, V φ , . . . , V φφφφ } in Fig. 11). This is consistent with the fact that the new data set brings no evidence for running or running of the running. It also explains why the results of this section are close to those of Sect. 4.1, obtained under the slow-roll approximation. Similar to 2015, the best-fit value of running for n = 3 is negative, but has moved down from −0.013 to −0.008, and remains compatible with zero at the 1.2σ level. For n = 4, the trend observed in 2015 to fit the data slightly better with a non-zero tensor contribution has disappeared. The decrease of the minimum effective χ 2 when moving from n = 2 to n = 3 is insignificant and even smaller than in 2015, showing that the data do not require anything more complicated than an approximately parabolic shape for the inflaton potential within the observable window. This can be checked by considering the random sample of well-fitting potentials presented in Fig. 12. Actually, for n = 4, a few of the plotted potentials have a non-parabolic "spoon-like" shape (with a kink and a plateau), because non-negligible values of |V φφφφ | are still allowed. However, this sub-class of models is by no means preferred over simpler parabolic-like potentials with a negligible |V φφφφ |; otherwise, we would have obtained a better χ 2 eff for n = 4. Hence one should not take from Fig. 12 the message that special potentials with a kink and a plateau are favoured by the Planck data. Comparing this plot to figure 15 of PCI15, we see that the models with the largest V(φ) amplitude are excluded by stronger bounds on the tensor modes.
Finally, it is interesting to notice that the predictions for the parameters of the minimal ΛCDM model, such as n s or τ, remain extremely stable when increasing the freedom in the inflaton potential.

Taylor expansion of H(φ) in the observable region
To assess the robustness of our method, in this section we repeat the analysis with a Taylor expansion of the Hubble function H(φ) in the observable window, as we did in 2015. We refer the reader to section 7.2 of PCI15 for a precise description of this analysis, and we recall that the difference with respect to the V(φ) reconstruction is more than a mere change of priors. For each value of n, the new parameterization covers a slightly different range of potentials, and, more importantly, it naturally includes a marginalization over the uncertainty in the initial value of the derivativeφ when the inflaton enters the observable window. Instead, in the previous analysis,φ was assumed to have reached the inflaton attractor solution, i.e., there was an implicit assump- Taylor-expanded at order n = 2, 3, and 4 in the observable region, making no assumption about the end of inflation, and using Planck TT,TE,EE+lowE+BAO+BK14. In natural units (where √ 8πM Pl = 1). We use the same scales as in PCI15. Note that there is another branch of solutions that is symmetric under tion that inflation started well before that time. In the analysis based on H(φ), inflation models with a minimal duration are not excluded by the priors.
The improvement with respect to the 2015 results is even more impressive in this case. Bounds on the n = 4 parameters are typically 3 to 4 times stronger compared with 2015, as can be checked from Table 7 and Fig. 13. We found that a factor of 2 improvement comes from switching to the new set of lowlikelihoods, and another factor of 2 from adding the BK likelihood. On the other hand, the use of more recent high-and BAO likelihoods has a modest impact.
A consequence of these improved constraints can be seen in For n = 4, the previously bestfitting models included many scenarios starting with a fast-roll stage, producing a tail with large V(φ) before pivot-scale crossing. These models are now excluded by better polarization data and tensor constraints.
Going beyond the parabolic approximation for H(φ) does not improve the goodness-of-fit: as in the potential-based analysis of Sect. 5.1, the ∆χ 2 s between n = 2, n = 3, and n = 4 are negligible, and the parameters ξ 2 H and 3 H related to H φφφ and H φφφφ are compatible with zero.

Taylor expansion of full V(φ)
We now present a new analysis with less conservative assumptions than in the previous subsections. We switch to the hypothesis that the inflaton potential is very smooth not only within its observable window, but also until the end of inflation, such that its whole shape can be captured by a Taylor expansion. We further assume that inflation ends when the first slow-roll condition is violated ( V = 1), without invoking any other field. Finally, we fix the number of e-folds between Hubble crossing of the pivot scale and the end of inflation to N * = 55, which implicitly relies  Table 7. Numerical reconstruction of the Hubble slow-roll parameters beyond any slow-roll approximation, using Planck TT,TE,EE+lowE+BAO+BK14. We also show the corresponding bounds on some related parameters (here n s , dn s /d ln k, and r 0.002 are derived from the numerically computed primordial spectra). All error bars are at the 68 % CL and all upper bounds at the 95 % CL. The effective χ 2 value of model n is given relative to model n − 1. on the hypothesis that no further inflationary stage took place at a later epoch. Technically, the analysis pipeline for this case is similar to that of Sect. 5.1, except for an extra step in which the CLASS inflationary module integrates the background equations until the end of inflation, goes backwards in time by 55 e-folds, and imposes that the Hubble crossing for the pivot scale k * = 0.05 Mpc −1 matches that time.
These models are much more constrained than those of Sects. 5.1 and 5.2, since the e-fold condition is imposed in addition to having a potential with a good shape within the observable window. The constraining power is then sufficient 21 Planck Collaboration: Constraints on Inflation for running the MCMC chains directly with flat priors on Our results are presented in Figs. 15 and 16 and in Table 8. Models with a purely quadratic potential always predict a large tensor-to-scalar ratio (of roughly 0.15) in tension with the Planck data, and even more so with the Planck+BK data. Thus the effective χ 2 is poor in the n = 2 case and improves considerably when adding some freedom in going to n = 3. Indeed, the presence of an additional cubic term allows us to reach smaller values of the tensor-to-scalar ratio for roughly the same scalar tilt, and lowers χ 2 eff by more than 10. Instead, when also adding a quartic term, we find no significant improvement in the goodness of fit, and the coefficient of the φ 4 term is consistent with zero. These findings are consistent with the global picture that Planck data prefer potentials which are concave in the observable window. The blue and green curves in the lower left panel of Fig. 16 illustrate the preference of the Planck+BAO+BK14 data for potentials with an inflection point, appearing qualitatively similar to scalar field potentials associated with spontaneous symmetry breaking models, hilltop models, new inflation, natural inflation, etc.
In these runs, the value of the scalar tilt running is always very precisely constrained around a value of dn s /d ln k −6 × 10 −4 . This does not come as a surprise if we keep in mind that these bounds are not imposed directly by the data, but rather by the class of inflationary potentials considered here, with potential parameters fixed by observational bounds on the amplitude and tilt of the scalar and tensor spectra. In other words, the running is not directly measured, but rather predicted as a function of the scalar/tensor amplitudes and scalar tilt. Interestingly, if combinations of future CMB and large-scale structure data with a wide lever arm in wavenumber space could become directly sensitive to such tiny values (which would require a factor of around 10 improvement in sensitivity compared to current CMB+BAO data), a very large class of currently successful inflationary models could be either confirmed or ruled out.   Table 8. Numerical reconstruction of the potential parameters beyond any slow-roll approximation, when the potential is Taylor-expanded to nth order, trusted until the end of inflation, and using Planck high-TT,TE,EE+lowE+BAO+BK14. We also show the corresponding bounds on some related parameters (here n s , dn s /d ln k, and r 0.002 are derived from the numerically computed primordial spectra). All error bars are at the 68 % CL and all upper bounds at the 95 % CL. The effective χ 2 value of model n is given relative to model n − 1.

Free-form potential reconstruction
As a complementary analysis to the previous three subsections, we next perform a free-form reconstruction of the inflationary potential with cubic splines, in a manner akin to the reconstructions of PCI15 and Sect. 6.2.1. Further plots and theoretical detail can be found in Handley et al. (2018). A free-form reconstruction usually proceeds by parameterizing the function of interest via a spline and taking the locations of the interpolation knots as free parameters in a posterior distribution. These are then varied along with any other model parameters, and then marginalized out to yield a model-independent reconstruction of the function of interest. The analysis is run 22 Planck Collaboration: Constraints on Inflation  for differing numbers of knots, N, and the Bayesian evidence is computed to allow for model comparison to determine how many knots are appropriate from the perspective of the data.
To reconstruct the inflationary potential V(φ), one cannot take a linear interpolating spline (as in Sect. 6.2.1), since the equations of motion in general depend on first (and sometimes second) derivatives of V. We therefore choose to parameterize the second derivative of the log-potential as a linear spline. The log-potential is computed by integrating this function twice, yielding a function with two additional free parametersa global offset and a gradient. Our reconstruction function is therefore Here LS(φ; θ) is a standard linear spline dependent on N knots, ln V * is the potential at the pivot scale, and dln V * /dφ is the gradient of the log-potential at the pivot scale. In general, any reconstruction of the potential will be sensitive only to the observable window of inflation in φ ∈ [φ min , φ max ], where φ min and φ max are defined as the field values when the largest and smallest observable scales k min and k max exit the Hubble radius during inflation. Regions of the potential outside these φ values are unconstrained by current CMB data. In our analysis, we take k min = 10 −4 Mpc −1 and k max = 10 −0.3 Mpc −1 , which encompasses the multipole range constrained by Planck (see Sect. 6.2.1). The locations φ 1 , . . . , φ N of the reconstruction knots should be distributed throughout this observable window. Whilst the locations φ 1 , . . . φ N and heights d 2 ln V 1 /dφ 2 , . . . , d 2 ln V N /dφ 2 themselves influence the size of the observable window, a reasonable approach is to first estimate it using the unperturbed potential (i.e., setting N = 0). This gives an alternative range φ ∈ [φ min ,φ max ]. The priors on all our variables are indicated in Table 9.
Our results are detailed in Fig. 17. The Bayesian evidence shows that the reconstruction preferred by the data is that using N = 1, corresponding to a constant non-zero d 2 ln V/dφ 2 . This indicates that the Planck data do not significantly constrain the  The top-right panel shows Bayes factors for the free-form potential reconstruction. The preferred reconstruction has N = 1, corresponding to a constant nonzero d 2 ln V 1 /dφ 2 . The remaining panels show reconstructions for the N = 8 knot case, focusing on the scalar primordial power spectrum, and the inflationary slow-roll parameters ε V and η V . Red lines indicate sample trajectories from the prior, whilst black lines are from the posterior. Technically the slow-roll parameters are defined as functions of φ, but we instead substitute this for the Hubble-radius-exit value to make for clearer comparison between posterior samples. In all plots, the approximate link between and k is via the Limber approximation, k/D A , where D A = r * /θ * is the comoving angular distance to recombination, which is at comoving distance r * .
inflationary potential within the window any further than up to a quadratic term in a Taylor expansion. It is illuminating, however, to consider adding further structure to the potential, and Fig. 17 shows reconstructions for N = 8.
Considering the predictive posterior of the primordial power spectrum, we see that our parameterization is sufficient to exhibit the deficit at 30, cosmic variance at low , and the loss of resolution at high , as seen in Sect. 6.2.1. Consistent with the rest of the analyses, ε V is unconstrained, whilst the Planck data provide relatively powerful constraints on η V within the observable window of inflation.
Indirect constraint [2,4] Table 9. Parameters of the free-form potential reconstruction analysis and details of the priors. There is a further prior constraint in that we require that the inflaton should evolve in an inflating phase throughout the observable window, that the inflaton should be rolling downhill from negative to positive φ throughout, and that any primordial power spectra generated sit in the range 2 < ln 10 10 P R (k) < 4.

Primordial power spectrum reconstruction
This section reports results for the non-parametric reconstruction of the primordial scalar power spectrum using the new Planck 2018 likelihoods, as well as comparisons with the previously reported results for the Planck 2013 and 2015 releases. The objective here is to search for deviations from a simple power-law primordial power spectrum (i.e., P R (k) = A s (k/k * ) n s −1 ) in a manner that does not presuppose any particular theoretical model giving rise to such deviations. This work is complementary to the searches considered in Sect. 7, where particular functional forms for such deviations motivated by theory are investigated.
Here we apply three distinct nonparametric methods. In 2013 only the first method was used to reconstruct the primordial power spectrum, the so-called "penalized likelihood" method, for which the 2018 results are presented in Sect. 6.1. In 2015 two additional methods were also used: a linear spline method (discussed in Sect. 6.2) for which both the number of knots and their positions were allowed to vary, and ideas from Bayesian model selection were applied to determine the appropriate number of knots; and a method using cubic splines (discussed in Sect. 6.3). Although the discussion below includes some description of each method in order to make the paper self-contained, for details the reader is referred to the 2013 and 2015 papers.
Here we specify only those details specific to the 2018 analysis or different from the choices in the 2013 and 2015 analyses.

Penalized likelihood
The underlying idea of the penalized likelihood approach is to add a term to the log-likelihood that penalizes deviations from a perfect power-law spectrum. We parameterize the power spectrum as where P 0 (k) = A s (k/k * ) n s −1 , and add the following term to −2 ln L: where κ = ln k. The interval [κ min , κ max ] is chosen to cover the range over which the likelihood is able to constrain the data. The two α terms serve to pin the reconstruction to the simple power law where the data have almost no constraining power. One may imagine that α > 0 should be infinite, but for numerical reasons a large but finite value is used to simplify the numerics. Numerically, for each λ the dimension of f is chosen to be so large that the continuum version of the penalty given in Eq. (53) has been accurately approximated. For more details see Gauthier & Bucher (2012) and the extensive references therein to prior literature, as well as PCI13 and PCI15.
In Fig. 18 we show the results using Planck TT+lowE and in Fig. 19 we show the results for Planck TT,TE,EE+lowE. In both cases we have assumed the usual base-ΛCDM model specified in PCP18, except that the power spectrum is now parameterized by a set of spline points. In addition to these spline points, we also maximize the likelihood with respect to the dimensionless Hubble parameter, h, and the baryon, Ω b h 2 , and CDM, Ω c h 2 , densities. All other cosmological and nuisance parameters are the same as those quoted in PCP18.
For the TT-only case, the maximum deviations are 1.55σ, 2.10σ, 1.80σ, and 1.65σ for λ = 10 3 , 10 4 , 10 5 , and 10 6 , respectively, for which the probabilities to exceed are 13 %, 28 %, 31 %, and 23 % (where we have taken into account the lookelsewhere effect). Similarly, for the TT,TE,EE case, the maximum deviations are 2.07σ, 1.77σ, 1.77σ, and 1.08σ for λ = 10 3 , 10 4 , 10 5 , and 10 6 , respectively, for which the probabilities to exceed are 29 %, 23 %, 32 %, and 25 %. We consequently find no statistically significant evidence for a deviation from the simple power-law hypothesis. This result is consistent with the results previously reported for the Planck 2013 and 2015 releases using essentially the same method. It is likewise consistent with the results below in Sects. 6.2 and 6.3, which use different methods.

Bayesian reconstruction
To reconstruct the primordial power spectrum of curvature perturbations, we follow the methodology of section 8.2 of PCI15, using an N-point interpolating logarithmic spline with the positions of the knots considered as free parameters in the full posterior distribution. The positions of the points in the (k, P) plane are treated as likelihood parameters with loguniform priors. Further, the k-positions are sorted a priori such that k 1 < k 2 < · · · < k N , with k 1 and k N fixed. We compute posteriors and evidence values (conditioned on N) using PolyChord (Handley et al. 2015a,b), also varying all cosmological and nuisance parameters. We then use evidence values for each model to correctly marginalize out the number of knots N.
To plot our reconstructions of P(k), we compute the marginalized posterior distribution of ln P conditioned on k. The iso-probability confidence intervals are then plotted in the (k, P) plane (see, e.g., Fig. 20), using code recorded in Handley (2018). To quantify the constraining power of a given experiment, we use the conditional Kullback-Leibler (KL) divergence as exemplified by Hee et al. (2016). For two distributions P(θ) and Q(θ), the KL divergence is defined as and may be interpreted as the information gain in moving from a prior Q to a posterior P (Raveri et al. 2016). For our reconstructions, we compute the KL divergence of each distribution conditioned on k and N, and then marginalize over N using evidence values to produce a k-dependent number which quantifies the compression or information that each data set provides at each value of k. Further plots and theoretical detail can be found in Handley et al. (2018).  In PCI15, our analysis focussed predominantly on the TT+lowTEB data set. Here we present results for TT,TE,EE+lowE+lensing. First, in updating to the lowE likelihood, we find that there is a marked tightening in the constraint on the amplitude of the reconstructed spectrum at all values of k. The improvement in the constraint can be seen directly in the predictive posterior plots (Fig. 20, top-left panel, and Fig. 21), and is quantified in Fig. 20 (bottom-right) via the KL divergence. The reason for the high-constraint provided by a low-likelihood change is due to the reduced uncertainty on τ that SimAll EE provides. This can be seen by examining the shifts in the underlying cosmological parameters in Fig. 2.
Upon adding TE and EE data, we find that the hint of a feature at 30 is still present, in spite of the additional constraining power provided by polarization. Using polarization data, the N = 3 case is now the most strongly favoured model by the evidence criterion. This indicates that there is some scope for models which account for low-cosmic variance to be preferred in a Bayesian sense. The other underlying cosmological parameters are unaffected by the additional degrees of freedom in the primordial power spectrum provided by the reconstruction.
In order to combine Planck polarization data with BK14, we also allow the tensor-to-scalar ratio r to vary, and fix the tensor tilt n t via the inflationary consistency condition. As can be seen in the bottom-left panel of Fig. 20, upon adding BK14, the effect of the low-deficit is softened, but with otherwise little change to the reconstruction.

Free-form search for features
Next we examine the effect that sharp features in the primordial power spectrum can have on cosmological parameters. We model sharp features in the spectrum as a variable number of tophat functions with varying widths, heights, and locations. On top of the traditional A s , n s parameterization of the power spectrum, we place N sharp top-hat features into the spectrum at locations k i with widths d i and heights h i (i = 1, . . . , N). That is, we set where the square brackets in the summation denote a logical truth function as introduced by Graham et al. (1994). For values of N = 0, . . . , 8, we treat the variables in parameterization (55) as parameters in a posterior distribution along with the traditional cosmological and Planck nuisance parameters, with priors as detailed in Table 10. We run with both linear and logarithmic priors on the k-locations of the features, as this alters the sensitivity to the type of features uncovered. We sample the posteriors using PolyChord (Handley et al. 2015a,b). Figures 22 and 23 show our results. With the linear priors case, there are statistically insignificant features corresponding to the peaks of the T T spectrum, which arise due to the enhanced cosmic variance at these locations. With the logarithmic priors case, a stronger but still statistically insignificant feature is detected at 30, with a small deficit and surrounding enhancement of power. This case reproduces the results found in Sect. 6.2.1. In both cases, the Bayesian evidence shows preference for a no-features spectrum, and steadily declines as more Parameter Prior Range Sorted log-uniform −4 < log 10 k 1 · · · log 10 k N < −0.3 d i Log-uniform 0 < log 10 d 1 , · · · , log 10 d N < 1 h i Uniform −1 < h 1 , . . . , h N < 1 ln(10 10 A s ) Uniform 2 < ln(10 10 A s ) < 4 n s Uniform 0.8 < n s < 1. 2   Table 10. Priors for the search for sharp features in the primordial power spectrum. Units for k and d are Mpc −1 .
features are added. The cosmological parameters remain unperturbed despite the introduction of features.

Cubic spline reconstruction
In this subsection we update the third method of reconstruction used in PCI15, in which ln P R (ln(k)) was expanded in cubic splines localized in ln(k) about uniformly spaced "knots," {ln(k b ), b = 1, . . . , N}, whose range was chosen to cover all relevant cosmological scales, from 10 −4 Mpc −1 to O(1) Mpc −1 . We single out the standard scalar power spectrum pivot scale as a "pivot knot" b = p, with k p = k * = 0.05 Mpc −1 . Its associated power ln A s = ln P R (k * ) is assigned a uniformly distributed prior. A tilted primordial power spectrum P R,fid ≡ A s (k/k * ) n s,fid −1 , with fixed spectral index n s,fid is used as the fiducial baseline from which deviations are measured, expressed in terms of N − 1 relative spectral shape parameters: For the results presented here, n s,fid = 0.967 was chosen. As in PCI15 we continue to use cubic splines for the k-space modes we expand in, with natural boundary conditions (i.e., vanishing second derivatives at the first and last knots). The treatment here is therefore quite analogous to that in Sect. 5.4, where the inflaton potential rather than the curvature power spectrum is expanded in cubic splines. Knot numbers up to 18 were reported in PCI15, and it was shown that 12 were sufficient to capture the variations desired by the Planck CMB data. The mode functions were also varied. For example, linear interpolation leads to similar reconstructions as long as enough knots are used. A weak uniform prior (−1 ≤ q b ≤ 1) was imposed on q b . Outside of the spline coverage region [k 1 , k n ] we set ln P R (k)/P R,fid (k) to be q 1 for k < k 1 and q n for k > k n . The prior on q b and boundary condition choices have little impact on the reconstructions over most of the k-range. The current Planck TT,TE,EE+lowE+lensing+BK14 data give only upper limits to the allowed values of the tensor amplitude, r < 0.07. Consequently, adding shape degrees of freedom to the tensor power spectrum would yield a completely priordriven result. Instead we adopt the standard model power-law parameterization for tensors P t (k) = rA s (k/k * ) n t with the tensor spectral index constrained by the consistency relation n t = −r/8. Without B-mode constraints and with enough knots one could deform the primordial scalar spectrum to mimic a tensor contribution to the CMB power. However, this near-degeneracy is broken with direct B-mode observations, even effectively so if there are only upper limits as for the BK14 data. Our reconstructions here focus on letting r float over a prior range 0 ≤ r ≤ 1, but the posterior is strongly constrained by the BK14 data.
The joint probability distributions of {q b , b p}, ln A s , and the other cosmic and nuisance parameters are determined by

PR)
Planck TT,TE,EE+lowE+lensing As more knots are added, the 30 feature in the C temperature spectrum is visible as a dip to lower power.    CosmoMC modified to incorporate the n-knot parameterization for fixed knot number N. Figure 24 shows the reconstruction. Apart from the mean and 1σ and 2σ limits on the ensemble of trajectories allowed by the posterior probability, we also show a set of individual trajectories with parameters taken from 1σ samples to illustrate the knot-to-knot coherence (dashed curves). The tensor trajectories are straight lines, as required by the adopted tensor power model. In spite of the extra scalar shape freedom in the k-space region over which the tensor modes affect the CMB, the 12 knot reconstruction still leads to a strong constraint of r < 0.084, rather close to the r < 0.07 limit obtained if the only shape parameter is n s . In fact we find that the current limits on r are such that the scalar-power reconstructions are not sensitive to the details of the r distribution. To illustrate this, the lower panel of Fig. 24 shows the spectrum when r is fixed at the tiny value of r = 0.001. One could regard this as a theoretical prior for lowenergy inflation models or a forecast for a future in which r is measured or tightly constrained by B-mode experiments.
In PCI15, the main cubic spline reconstruction included non-CMB data to help pin down cosmological parameters such as H 0 , τ, and the late-time expansion history. The improvements in the data from 2015 to 2018, especially the decreased errors on τ, result in no non-CMB data being needed. Although τ and ln A s are about 90 % correlated, as they are in the standard power-law model, neither are very correlated with the q b . The strongest is about 40 % for q 3 at k 0.0006 Mpc −1 corresponding to ∼ 10 where reionization is kicking in. The second strongest is about 30 % for k 0.02 Mpc −1 , similar to the correlation of τ and n s in the standard power-law model. (The correlations among the q b are also relatively small, except for the high k bands b = 10 and 11, where the data are not constraining.) The middle panel in Fig. 24 shows the effect of adding the BAO constraint. Although apparently visually identical, there are slight differences. For example, the 1σ error on q b at k 0.02 Mpc −1 decreases by about 6 %, from 0.0086 to 0.0083, while at k 0.1 Mpc −1 the decrease is about 5 %. The restriction to r = 0.001 does not change the error bars over the floating r case. At intermediate k for modes 5 to 8 the errors on q b are so close to zero that the reconstruction is quite compatible with a simple power law, corresponding to a straight line in Fig. 24. This was also a main result of the 2015 Planck reconstructions.
The errors on the q b grow above ±0.1 for b = 4 as a consequence of increased cosmic variance, giving more freedom in the allowed trajectories. Unfortunately this is also the region of relevance to the T T power spectrum deficit in the 20-30 range. The most significant deviation from zero occurs for q 4 at k 0.0014 Mpc −1 : −0.259 ± 0.125, −0.259 ± 0.123, and −0.235 ± 0.122 for the three cases. Thus the anomaly in terms of deviation from the power law of the standard model hovers at around the 2σ level. More precisely, the 95 % upper confidence limits on q 4 are −0.015, −0.020, and −0.00023, for the respective cases. This 2σ level of the anomaly was also the conclusion of the 2015 Planck reconstructions. Therefore, even though the low-k deficit is robust against the various choices for the reconstruction, we conclude that it is not statistically significant. The associated T T , T E, EE, and BB power spectra responses to the allowed primordial power variations are derived from the mode expansion, and match the D XY data well, in particular following the dip in T T at = 20-30.
In Fig. 25 we show the reconstructed power spectra using only the TT, TE, and EE data in conjunction with BK14. The fixed r = 0.001 cases look very similar. Except at high k, the polarization data using either EE or TE also enforce a nearly uniform n s (k) over a broad range in k, with values in excellent agreement with those obtained from TT alone, from TE and EE used in combination, and from the combined TT,TE,EE results. For example, the ±0.0083 and ±0.0059 1σ errors at k 0.02 Mpc −1 32 Planck Collaboration: Constraints on Inflation 10 1 10 2 10 3 10 4 ℓ k ≡ kD rec 10 0 10 1 10 2  Reconstructed primordial scalar power spectrum derived using Planck TT,TE,EE+lowE+lensing+BK14 data and 12 knots for the cubic spline interpolation (with positions marked as ∆ at the bottom of each panel). Mean (ensemble-averaged) spectra are heavy lines, allowed ±1σ and ±2σ regions for trajectories are the shaded regions, and the dashed lines denote selected trajectories with parameters sampled within the ±1σ posterior. Below the scalar power is the tensor power reconstruction. The addition of the BAO likelihood shown in the middle panel makes almost no visual difference to the reconstructions. In the bottom panel, fixing the tensor-to-scalar ratio to r = 0.001 also shows only small differences in reconstruction. Knot positions in k roughly translate to multipoles through kD rec , where D rec is the comoving distance to recombination.
10 1 10 2 10 3 10 4 ℓ k ≡ kD rec 10 0 10 1 10 2     figure), with cubic-spline interpolation and the Planck TT,TE,EE+lowE+lensing+BK14+BAO data for the two cases of floating r and r fixed at 0.001. Sample 1σ trajectories for the floating r case allow wide variability, which is naturally greatly diminished if r is fixed to r = 0.001. and k 0.1 Mpc −1 , respectively, increase only slightly for TT only, to ±0.012 and ±0.0067, but, more significantly, to ±0.017 and ±0.069 with EE alone. The deficit region remains about the same, with the TT,TE,EE result for q 4 of −0.259 ± 0.125 quoted above changing slightly for TT alone, to −0.253 ± 0.137, but with no hint of an anomaly for EE alone, at −0.132 ± 0.445. If just the TE cross data are used, the values are closer to the TT case, namely, −0.232 ± 0.162, now under a 2σ excursion from the tilted fiducial model.
As in PCI15, we can use the P R (k) ∝ H 2 / and P t (k)/P R (k) 16 reconstructions to get an idea of the history of the acceleration of the Universe as a function of time over the significant number of e-folds of the cosmic expansion that the CMB data probe, codified by the dynamical slow-roll parameter (k) = −Ḣ/H 2 | k=aH , considered as a function of aH, the value of the wavenumber at Hubble crossing. Results with floating r and r fixed to 0.001 are shown in Fig. 26. For the dynamical time variable we use k = aH for the horizontal axis for ease in comparing with the P R (k) curves of Fig. 24. The wide spread in the trajectories for the floating r case is a consequence of being able to fit the n s (k) shape by a combination of (k) and Top: Reconstructed shape of the single-field inflaton potential from the cubic-spline power spectra mode-expansion using 12 knots and the Planck TT,TE,EE+lowE+lensing+BK14+BAO data. Bottom: Result when r is fixed at 0.001. Instead of plotting as a function of wavenumber k we plot ln V(φ)/V pivot about a pivot field value φ pivot . Note that the range on the φ axis is quite different for the small r case than the floating case. d ln (k)/d ln k. When r is small, n s (k) is almost entirely determined by d ln /d ln k, and the (k) values cluster near r/16.
The Hamilton-Jacobi energy constraint equation relates the potential to and H via V = 3M 2 P H 2 (1 − /3). Figure 27 shows the reconstructed inflationary potential shapes in the region over which the allowed inflationary potentials are constrained by the data for the floating r and fixed r cases. Instead of using k for the horizontal axis, we translate into inflaton-field φ-space using the relation between φ and √ . It is referenced to the pivot position φ pivot . For the vertical axis we plot ln V/V pivot , with the overall normalization V pivot removed. Its value is set by r, hence there is a distribution of constant V pivot amplitudes to superimpose if we want the total V. The radically different visual appearance for the floating r and fixed r cases is due to the observable k range being compressed through the smallness of into a small precisely determined field range, whereas this range has a distribution in the floating r case. One can monitor whether the shapes of the individual realizations of the potential trajectories bend upwards or downwards or do both, an indication of convexity. The sample trajectories shown are not exclusively convex or concave, and a measure of the probability that they are convex can be made from the ensemble. As indicated in Fig. 27 for the 12 knot case, the ensemble-averaged potentials are roughly exponential, with individual trajectories bending away from the mean, but with no strong tendency for convexity or concavity. (The rougly 50 % probability changes somewhat depending upon the combination of data used, whether TT,TE,EE or the individual data sets.) The standard cosmological parameter determinations are highly robust to the addition of these spline shape degrees of freedom. The mean values change little and the error bars grow slightly, by around 10 % for ln A s , τ, and H 0 . The largest error increase is for σ 8 , with σ 8 = 0.811 ± 0.0059 becoming 0.813 ± 0.010. The main conclusions of this section on and V, and P R (k), remain as in PCI15, but the results have been noticeably sharpened by the improvements in the Planck 2018 data sets.

Search for primordial features in the Planck power spectrum
The "bottom-up" power spectrum reconstruction methods of the previous section are an excellent way to search for coarse features in the spectrum, but lack the resolution to detect the higher-frequency features generically predicted by various physical mechanisms [see, e.g., Chluba et al. (2015) for a review]. It is therefore useful to complement power spectrum reconstruction with a "top-down" approach by fitting specific feature models to the data. In this section we will analyse a representative range of power spectrum templates which parameterize features in terms of a handful of new parameters. With Planck's temperature and polarization data, we have two essentially independent probes of features at our disposal and will pay particular attention to examining the consistency between the two .

Global oscillation models
Periodic or quasi-periodic modulations of the power spectrum which extend over the entire observable range of wavenumbers can occur in a variety of models [cf., e.g., Danielsson (2002); Martin & Brandenberger (2003); Bozza et al. (2003); Chen (2012); Jackson & Shiu (2013)]. A general parameterization of models with a sinusoidal modulation of the primordial power spectrum reads where X ∈ {log, lin, rf}. Defining κ ≡ k/k * 12 , we consider the logarithmic oscillation model, given by Ξ log ≡ ln κ, and the linear oscillation model, Ξ lin ≡ κ. In addition, we investigate a logarithmic model with running frequency, Ξ rf ≡ ln κ (1 + α rf ln κ). For 0 ≤ α rf < ∼ 0.01, this is a good approximation for the scalar power spectrum in the axion monodromy model (Flauger et al. 2017a), which will be analysed in more detail below, but here we allow for a wider range of the running parameter α rf , including negative values (i.e., decreasing frequency with increasing k).

Localized oscillatory features: inflation with a step
A sudden transient event in the evolution of the inflation field, triggered by a sharp feature in the inflaton potential, or a sharp 12 Throughout this section, we take k * = 0.05 Mpc −1 . turn in field space, generically leads to a localized oscillatory feature in the power spectrum (Adams et al. 2001;Chen et al. 2007;Achúcarro et al. 2011;Miranda et al. 2012;Bartolo et al. 2013). As an example of this class of feature models, we consider here the case of a tanh-step in an otherwise smooth inflaton potential (Adams et al. 2001), whose power spectrum can be parameterized as (Miranda & Hu 2014) where the first-and second-order terms are given by with window functions W 0 (x) = 1 2x 4 18x − 6x 3 cos 2x + 15x 2 − 9 sin 2x , (60) and damping function In this model, the parameter A s determines the amplitude of the oscillatory feature, the step scale k s sets the position of the step in k-space, and the damping parameter x s determines the width of the envelope function.

Models with suppressed power at large scales
The apparent lack of power at the largest scales in the temperature power spectrum with respect to the expectation of base ΛCDM serves as a motivation for models with a suppression of primordial perturbations below a cutoff scale k c . Physically, this effect may be due to fluctuations at the largest observable scales being generated at the onset of the inflationary phase after a prior era of, e.g., kinetic or radiation domination (Vilenkin & Ford 1982;Contaldi et al. 2003), or due to an isolated event such as kink the inflaton potential (Starobinsky 1992). In these scenarios, the primordial spectrum can generally be analytically approximated by an expression of the form ln P Y R (k) = ln P 0 with Y ∈ {kin, rad, kink}, where Υ Y is a function with ln Υ Y → 0 in the limit k k c Y that describes the shape of the cutoff and the transition to a power-law spectrum at smaller scales.

Initial kinetic domination
If inflation is preceded by an era dominated by the kinetic energy of the inflaton field (i.e., fast roll), we have with D c (y) = e i y H (2) where H (2) n denotes the Hankel function of the second kind (Contaldi et al. 2003 Table 11. Prior ranges for the parameters of the feature model templates of Sect. 7.1.

Kink in the inflaton potential (Starobinsky model)
A kink in the inflaton potential, first discussed by Starobinsky (1992), leads to a spectrum approximately given by with the parameter R c expressing the ratio of the slopes of the inflaton potential before and after the kink (Sinha & Souradeep 2006).

Data analysis
We employ a modified version of CAMB with suitably increased numerical precision settings to calculate the CMB angular power spectra for the feature models. Since variations of the primordial spectrum may be degenerate with late-time cosmology parameters (Obied et al. 2017), we explore a parameter space consisting of the base-ΛCDM parameters and the respective additional free parameters of the feature models (see Table 11 for the prior ranges). Note that we take primordial tensor perturbations to be absent in our analysis. In the results presented in Sect. 7.3, nuisance parameters are assumed to be uncorrelated with the feature parameters and kept fixed to their base-ΛCDM best-fit values. In order to maximize sensitivity to narrow features, we use only the unbinned versions of the Planck high-likelihoods in the following combinations: (i) temperature data, Planck TT(unbinned)+lowE; (ii) E-polarization data only, Planck EE(unbinned)+lowE; and (iii) temperature plus polarization data, Planck TT,TE,EE(unbinned)+lowE.
For all combinations of feature models and data, the parameter space is sampled with the nested sampling algorithm as implemented in MultiNest. The improvement in the fit due to the introduction of a feature is quantified by the effective ∆χ 2 ≡ −2(ln L best fit ΛCDM − ln L best fit feature ). Being more complex than a power-law spectrum, feature models will in general have a negative ∆χ 2 . However, determining whether the improvement in fit is due to overfitting scatter in the data or due to an actual feature is not straightforward and requires model-dependent simulations (PCI15) or analytic estimates  to determine the expected ∆χ 2 under the null-hypothesis of an underlying power-law spectrum. In the Bayesian approach, a feature model's general performance relative to base ΛCDM can be expressed in terms of the Bayesian evidence E (Trotta 2007a), which is also evaluated by MultiNest.

Feature candidates and their evidence
We list the best-fit effective ∆χ 2 and Bayes factors with respect to a power-law spectrum in Tables 12 and 13. Examining the effective ∆χ 2 for the feature models previously considered in PCI15 reveals only minor differences, with a general trend towards smaller improvements due to features. The ∆χ 2 of the oscillation and step models fall well within the expected range of ∆χ 2 ∼ 10 found in PCI15. Of note are the relatively high values of the radiation and kink cutoff models for polarization-only data, partially driven by the high quadrupole of the EE data. However, the best-fit parameters and spectra (see Fig. 28) do not match their counterparts in the temperature data at all, which strongly suggests that this is not a physical effect. The same observation can also be made for the step model: the best fit to the EE data is clearly out of phase with the temperature best fit.
We find a similar conclusion for the three oscillation models. As can be seen from the profile likelihood of the frequency parameters in Fig. 29, the likelihood peaks in the modulation frequencies do not match up between the TT and EE data sets. Furthermore, the preferred modulation amplitude for the EE data is in all cases much larger than that for the TT or TT,TE,EE data-given that the polarization data are noisier, this behaviour would be expected for a procedure that is overfitting the data.
Consequently, the Bayesian evidence for all combinations of models and data lies between barely worth mentioning and substantial evidence against the feature model on the Jeffreys scale. This implies that, currently, the Planck data do not show a preference for the feature models considered here.
It may also be worth pointing out that in models with oscillations linear in k, the wavelength of the corresponding modulation of the angular power spectra matches that of the CMB's acoustic oscillations, ∆ 300, if log 10 ω lin 1.15. One might therefore suspect that features with frequencies around this value and carefully tuned amplitudes and phases could in principle mimic the (unphysical) effect of a lensing parameter, A L 1. However, for a model with k-independent modulation amplitude A lin , we find no correlation between A L and A lin for frequencies in the vicinity of the BAO frequency, due to different dependence of the 36 Planck Collaboration: Constraints on Inflation Step Kin cutoff Rad cutoff  Kink cutoff  TT  EE  TT,TE,EE  TT  EE  TT,TE,EE  TT  EE  TT,TE,EE  TT  EE  TT,TE, Table 13. Same as Table 12, but for the oscillatory feature models.
∆D 's. Explaining the lensing discrepancy would thus require a model with a tuned scale-dependent linear modulation of the primordial spectrum. Additionally, while the phenomenology of A L and linear modulation models is similar for temperature and polarization spectra individually, the two scenarios are in principle distinguishable by a combination of temperature and polarization data. This is due to the phase difference of the acoustic peaks in T T , T E, and EE, which leads to similar phase differences for the residuals when varying A L -unlike modifications of the primordial spectrum. However, for features with an amplitude chosen to resemble the apparent lensing excess in the Planck TT data, the Planck TE and EE data are not sensitive enough to make this distinction.

Axion monodromy
As in section 10.3 of PCI15, we next derive constraints on the underlying parameters in axion monodromy inflation (Silverstein & Westphal 2008;McAllister et al. 2010;Kaloper et al. 2011;Flauger et al. 2017a), which within string theory motivates a broad class of inflationary potentials of the form where µ, Λ 0 , f, and φ 0 are constants which have dimensions of mass, while C 0 , p, p Λ , p f , and γ 0 are dimensionless. In the literature, one can find theoretically motivated models with p = 3, 2, 4/3, 1, and 2/3 (Silverstein & Westphal 2008;McAllister et al. 2010McAllister et al. , 2014. In the following, we neglect a possible amplitude drift in the modulation amplitude by fixing C 0 = p Λ = 0, focussing instead on a possible frequency drift p f , as was done in previous analyses (Peiris et  Due to its oscillating nature, a numerical study of this model is restrictive (Peiris et al. 2013). As such, we employ the semianalytic template (Flauger et al. 2017a) used in previous analyses, namely (70) We neglect the effect of small oscillations in the tensor primordial spectrum, and approximate it as a power law with a very small spectral index n t (fixed by the single-field slow-roll self-consistency condition). The most well studied case to date is for p = 4/3, but given the high tensor-to-scalar ratio predicted by this model and the current upper bounds on r given in Sect. 3.5, we extend our study to the cases of p = 1 and p = 2/3. Furthermore, to completely specify this template, we assume instantaneous reheating, which, for a pivot scale of k * = 0.05 Mpc −1 , corresponds to N * ≈ 57.5, and φ 0 = 12.38M Pl with φ end = 0.59M Pl . This leads to definite predictions for (r, n s ); namely, (0.0922, 0.971) for p = 4/3, (0.0692, 0.974) for p = 1, and (0.0462, 0.977) for p = 2/3. To constrain this model, we carry out a Bayesian analysis using a modified version of CLASS (Lesgourgues 2011;Blas et al. 2011), which has been adapted to allow for a full parameter exploration, using the aforementioned template. As part of these modifications, special care needs to be taken to ensure that a correct sampling ∆k in wavenumber space is chosen, at two different levels in the Boltzmann code: when computing an interpolation table for the primordial spectrum of scalars and tensors; and when performing the integral over the squared photon transfer functions multiplied by the primordial spectra to get the multipoles C . This sampling needs to be fine enough to guarantee 37 Planck Collaboration: Constraints on Inflation Step 10 9 Fig. 28. Best-fit and central 95 % CL regions for the primordial power spectrum in the three cutoff and the step models for TT data (red curves), EE data (green), and TT,TE,EE data (blue). Lin osc log 10 ω X Fig. 29. Profile likelihood of the frequency parameter in the three oscillatory feature models for TT (red curves), EE (green), and TT,TE,EE data (blue). The dotted grey line in the bottom panels marks the frequency for which the linear oscillation model leads to a modulation of the angular power spectra whose wavelength roughly matches that of the CMB's acoustic oscillations. Note the lack of alignment between the temperature and polarization likelihood peaks in the vicinity of this frequency.  that no features are smoothed out or lost in this convolution, and we checked carefully that this is the case in our runs. The grid of values at which the C 's are actually computed and not just interpolated also needs to be refined. We fit to the data the five cosmological parameters {ω b , ω c , θ, A s , τ} plus the frequency f of the underlying axion decay constant, the frequency drift p f , and the oscillation amplitude δn s . We adopt the same priors used in previous analyses: −4 ≤ log 10 ( f /M Pl ) ≤ −1 for the frequency; −0.75 < p f < 1 for the frequency drift; and an upper bound on the amplitude of δn s < 0.5. Furthermore, for the phase parameter ∆φ we take a uniform prior of −π < ∆φ < π.
In all three cases, we find two expected asymptotic behaviours. First, when the frequency is very high (which means that f is small in our parameterization), the oscillations in the primordial spectrum are smoothed out in the angular power spectrum, and the oscillation amplitude parameter δn s becomes irrelevant and unconstrained. Second, in the limit of a very small amplitude parameter δn s , the oscillations become undetectable and the parameter f is also unconstrained. In all cases, no preferred frequency drift is found, which is compatible with previous analyses.
We recover the complex structures (highlighted in red in Figs. 30, 31, and 32) found in previous analyses in the frequencyfrequency drift parameter space, which, as was discussed in PCI15, arise due to underlying modulations in the data and the model (Easther et al. 2005). These structures become more apparent as we reduce the index p.
We perform a χ 2 comparison with the minimal 6-parameter ΛCDM model, and find ∆χ 2 (p=4/3)/ΛCDM = 0.4, ∆χ 2 (p=1)/ΛCDM = 0.6 and ∆χ 2 (p=2/3)/ΛCDM = 1.2. The reason for higher χ 2 in the axion monodromy models, despite the addition of extra parameters, is that the predicted r values are in tension with the CMB data. This shows that, overall, axion monodromy models are disfavoured due to their high tensor-mode amplitude. In order to check specifically whether the data give any hint of oscillatory patterns in the primordial spectrum matching the axion monodromy template, as well as to compare with the results discussed in the previous subsection, we fitted the data with ΛCDM+r models in which r and n s were fixed to the same values as in the axion monodromy model with p = 2/3, 1, and 4/3. In each case, the comparison between axion monodromy and ΛCDM+r with the same (r, n s ) gives ∆χ 2 (p=4/3)/ΛCDM+r = −7.8, ∆χ 2 (p=1)/ΛCDM+r = −7.6 and ∆χ 2 (p=2/3)/ΛCDM+r = −8. That is, in all cases we find ∆χ 2 ∼ 10, which is compatible with the general results shown in Table 13. With three more free parameters, these improvements are statistically insignificant, and we conclude that the data show no preference for axion monodromy models.
39 Planck Collaboration: Constraints on Inflation

Approach
This section establishes constraints on oscillatory models using the power spectrum and the bispectrum simultaneously. Oscillatory features can appear in multiple correlation functions (Chen et al. 2008;Meerburg et al. 2009;Flauger et al. 2010;Flauger & Pajer 2011;Achúcarro et al. 2011;Achúcarro et al. 2014;Flauger et al. 2017a,b) (see, e.g., Chluba et al. (2015) for a recent review). More powerful constraints result when spectra of various orders are combined (Palma 2015;Mooij et al. 2016;Gong & Yamaguchi 2017). Past work has suggested that the statistical weight of the oscillations in the bispectrum (or higher-order correlation functions) is less than that in the power spectrum ; however, counterexamples exist as well (see, e.g., Behbahani & Green 2012). The analysis in Sect. 7 used the Planck data to establish stringent constraints on the presence of features in the power spectrum. The 2015 Planck data were analyzed to constrain non-Gaussianties containing features (Planck Collaboration XVII 2016), where, as in the power spectrum analysis, several candidate features were identified at low statistical significance. The analysis here focuses on the location, or frequency, of the feature. Joint analyses of the power spectrum and bispectrum were discussed in several studies Meerburg et al. 2016). We apply some of the tools developed there to the Planck temperature and polarization data. The analysis here is incomplete and limited in several respects. First, we analyse the bispectrum keeping all cosmological parameters fixed. Second, the parameters varied in the bispectrum are not varied in the Bayesian sense. The bispectrum is analyzed using a best-fit analysis based on how well a template shape fits the data. Third, the data suggest the primary bispectrum is close to zero and its covariance dominated by the scalar contributions in the power spectrum. We find that an ideal Bayesian analysis is not computationally feasible (see, e.g., Verde et al. 2013).
The output from the bispectrum analysis for features provides us with a map that specifies the significance of a feature in units of σ, given the location (frequency) and the phase of the feature. We can turn this map into a likelihood, which we can simply add to that of the power spectrum; in other words, we take where c represents the standard cosmological parameters, f the foregrounds, ω P,B the frequency, A P,B the amplitude, and φ P,B the phase of the modulation in, respectively, the power spectrum and the bispectrum. We assume vanishing covariance between the power spectrum and the bispectrum, which has been shown to be a good approximation Meerburg et al. 2016). Furthermore, the likelihood ln L BS is not normalized (more precisely ln L tot is not normalized in a universe with a non-zero bispectrum). Strictly speaking, we do not have a likelihood that measures A B with a certain probability. Furthermore, several studies have shown that the frequency parameter, in combination with the amplitude and the phase, does not obey a χ 2 fitting to the data Meerburg et al. 2014b,c,a;Meerburg 2014;Easther & Flauger 2014). Removing the frequency from the search results in a χ 2 distribution with two degrees of freedom. Since ln L tot is rather large (of order 10 4 when combining all data), we can change the equation above by limiting ourselves only to improvements that are driven by ω B ∼ ω P ≡ ω, that is, Assuming that φ B and A B are well described by a twoparameter χ 2 distribution, we can now convert our σ map into a χ 2 improvement via or, in terms of the likelihoods, −2 ln L tot = −2 ln L PS (c, f, ω, A P , φ P |dat) We will use the above expression to derive the posterior of the joint fit.

Models
We will focus on two models: the local or linear feature model and the log feature model. For the log model we set {A P , ω, φ P , A B , φ B } = {A log , ω log , φ log , B log ,φ log } and use for the power spectrum with P 0 (k) = A s (k/k * ) n s −1 . For the bispectrum we use cos ω log log k ĩ k 0 +φ log . For the linear model we follow  with and In both models we choosek 0 = 1 Mpc −1 , which is different from the choice in Sect. 7 for the linear model. As a result the linear frequencies can be related using ω lin,Sect8 = 10ω lin,Sect7 . The pivot scale is set to the usual value k * = 0.05 Mpc −1 .

Power spectrum
Our analysis uses a modified version of CAMB (Lewis et al. 2000), which is capable of adaptively changing the sampling in both k and depending on the frequency of the feature, allowing us to scan a wide range of frequencies. As in the previous section, we use the unbinned versions of the Planck high-likelihoods for temperature plus polarization, in combination with lensing and large-scale temperature and polarization (i.e., lowE). We compared the power spectrum results for the limited frequency range considered in Sect. 7 for the log and linear model, and find excellent agreement (sampled with Multinest). We developed a bispectrum likelihood module based on Eq. (74) using the 2015 data analysis (Planck Collaboration XVII 2016) for both the log and linear feature models, obtained using optimal estimators following Münchmeyer et al. (2014Münchmeyer et al. ( , 2015 and . For the log model, the frequency range is set to 10 ≤ ω log ≤ 1000, while for the linear model we consider the frequency range 10 ≤ ω lin ≤ 3000. This joint analysis excludes the very low frequencies known to (weakly) correlate with cosmological parameters. The cosmological parameters are held fixed in the bispectrum analysis. We consider amplitudes 0 ≤ A log,lin ≤ 0.9; the highest amplitudes will only be allowed for high frequencies where projection suppresses the power of the oscillating part in the power spectrum significantly. The phase is varied and marginalized over in the joint analysis. We use the PolyChord sampler (Handley et al. 2015;Handley et al. 2015b), which is powerful enough to include foregrounds (with n live = 512).

Bispectrum
The bispectrum likelihood is derived from the posterior distributions generated in Planck Collaboration IX (2018). Although the linear bispectrum of Eq. (78) can easily be factorized, the log bispectrum of Eq. (76) is not of the factorized form. Using modal techniques developed by Fergusson & Shellard (2009) and Fergusson et al. (2010Fergusson et al. ( , 2012, any shape can be factorized, with a close-to-optimal estimator. The modal method converts the angular-average bispectrum into a set of factorizable orthogonal mode functions. These functions can be directly constrained using foreground-cleaned CMB maps. From these measurements, a large number of bispectra can be reconstructed and constrained by appropriately weighting the mode functions. The convergence of this method, in terms of how many mode functions are required to accurately reconstruct the shape of interest, depends on the choice of the mode functions. In the 2015 analysis, two different mode functions were used: a polynomial-based reconstruction; and a trigonometric-based reconstruction. The latter was developed by Münchmeyer et al. (2014Münchmeyer et al. ( , 2015 and relies on expanding around linear oscillations. The polynomialbased reconstruction is extremely powerful for most bispectra, but is non-optimal for oscillatory bispectra, which require a large number of modes (e.g., more than 2000 for ω log = 50). Trigonometric modes allow for faster convergence and provide good reconstruction for much higher frequencies, both for linear-and log-type modulations. For low frequencies, both methods can be compared and results show excellent agreement (Planck Collaboration XVII 2016). In addition, both methods were developed independently, which provides further confidence in the results. In the analysis presented here, we use the results obtained using the trigonometric mode functions. Further details can be found in Planck Collaboration IX (2018).

Estimating significance
Next we will estimate the significance of the improvements driven by the joint analysis. For this purpose we generate 100 mock spectra as in Meerburg et al. (2016) without features and perform an analysis jointly with true CMB power spectrum data, i.e., we use the same power spectrum likelihood (real data) in combination with the simulated bispectrum likelihoods (mock data). We will do this for both the linear and the log models, with 100 simulations in total. Each analysis requires a similar amount of time as does the real data analysis, using about 12 000 CPU hours for the linear feature and about 40 000 CPU hours for the log feature per simulation. More details on the simulated spectra can be found in (Planck Collaboration IX 2018). These simulations help us assess the statistical significance of our results. Improvement in fit is given in units of χ 2 compared to a no-feature model as defined the previous section [i.e., ∆χ 2 −2(ln L best fit ΛCDM − ln L best fit feature )]. The left panel of Fig. 33 shows the typical best-fit improvement from a set of simulations for the log feature model. This first analysis shows that the best fit in the data is perfectly consistent with 41 Planck Collaboration: Constraints on Inflation a standard ΛCDM universe, without features, with P(∆χ 2 ≥ ∆χ 2 data ) = 28 %. This outcome is not unexpected, given earlier analyses for the power spectrum (see, e.g., Meerburg et al. 2014c;Easther & Flauger 2014;Benetti 2013;Miranda & Hu 2014;Planck Collaboration XX 2016;Hazra et al. 2016) and the significance of features in the bispectrum alone (Planck Collaboration XIII 2016). The lookelsewhere effect lowers the significance of features and by jointly constraining features in the power spectrum and the bispectrum it is possible to alleviate some of this suppression. To quantify this, we consider the following two questions: 1) considering the various frequencies with ∆χ 2 improvement over no features in the joint analysis, how many of those were present in the power spectrum analysis only; and 2) what is the mean improvement, in units of ∆χ 2 , of these fits? We compare the results of the simulations, which do not contain any real features, to the data.
Before we answer these questions we need a criterion to decide if two frequencies will be considered the same or not, i.e., we need a frequency correlation measure. We will consider a simple ansatz, which will have an analytical solution and will serve to estimate the correlation between frequencies, by defining Next we marginalize over phase, defining x max cos ∆ω 12 log (x max ) − cos ∆ω 12 log (x min ) x min + ∆ω 12 sin ∆ω 12 log (x max ) x max − sin ∆ω 12 log (x min ) x min , where ∆ω 12 = ω 1 − ω 2 . For linear oscillations we can derive a similar measure, with g lin [ω 1 , ω 2 ] = π 2∆ω 12 [sin (2∆ω 12 x max ) − sin (2∆ω 12 x min )] .
The correlation is given by The parameters x min and x max play a role in determining the correlation length. Although strictly speaking they correspond to the minimum and maximum scales observable in the CMB, they can be used to model the correlator to allow for shifts in the frequency coming from a non-optimal analysis. We argue that this is reasonable given the low number of peaks in the analysis. We tested the above on various nearby peaks in the data and found that demanding Cor(ω 1 , ω 2 ) ≤ 0.1 is generally sufficient to effectively identify independent peaks. We explored the sensitivity of the results to the correlation criterion of 0.1. First we increased it to 0.3 and found that in this case many peaks were missed when counting the number of aligned peaks (with little effect on determining the peaks). When we lowered the criterion to 0.01, we obtained many aligned peaks that should not be aligned. Small changes in the correlation criterion have minimal effect on the results presented here. Ideally, more simulations should be generated, which would help to establish the best choice for the correlation criterion. We found that the choice of x min does not affect the correlator as long as x min 1. We set x max = 0.05 for linear oscillations, which roughly correlates peaks with ∆ω lin ∼ 10, which is within the tails of the observed widths of the peaks in the power and bispectrum analysis. For log oscillations we choose x max = 1, which has the advantage that the correlator has no zero-crossings near the peak (but hardly effects the correlation length). We find ∆ω log ∼ 1 at ω log = 100, which seems reasonable in light of the power and bispectrum peak widths (see Fig. 29).
In Fig. 34 we show the number of peaks in the joint analysis that have improved (left panel) as well as their mean improvement (right) over a no-feature analysis. We find P(#peaks ≥ #peaks data ) = 16 % and those peaks do not lead to significant improvements in the joint χ 2 , with P(∆χ 2 peaks ≥ ∆χ 2 peaks,data ) = 27 %. Assuming that these 100 simulations provide a fair sample of the noisy data, we conclude that there are no significant features present in both the power spectrum and the bispectrum for the models considered within the chosen range of feature parameters.
We carry out the same analysis for linear features and show the results in Fig. 33 (right panel), deriving a typical best fit from 100 simulated noisy spectra. The true best fit, as derived from the joint analysis of the 2015 bispectrum and the 2018 power spectrum, shows a relatively small improvement with P(∆χ 2 ≥ ∆χ 2 data ) = 82 %. Further correlated-features searches find 5 features within the frequency window which may be considered aligned, with a mean ∆χ 2 of 12.3, as illustrated in Fig. 35. Compared to 100 simulated noisy spectra, we obtain P(#peaks ≥ #peaks data ) = 28 % and P(∆χ 2 peaks ≥ ∆χ 2 peaks,data ) = 53 %. Similar to the log feature, the analysis shows that although the number of aligned peaks in the real data is relatively high, the mean improvement is not improbable. Since the overall improvement from fitting these aligned peaks does not exceed the 3σ threshold, we conclude that there is no statistically significant evidence for any of these features.
We conclude that the simple parameterization considered in this analysis does not provide any evidence for features. The two models analyzed are representive of a broad class and have well-studied phenomenological spectra; however, other classes of models exist. Features, for example, can have scale dependence [i.e., an "envelope" (Chen 2010;Achucarro et al. 2014;Torrado et al. 2017)]. Likewise, more realistic modelling of axion models shows that the frequency could depend on scale [i.e., "running" (Flauger et al. 2017a,b)]. Both these possibilities could substantially change the spectra and likely the joint analysis and significance.
Here we have not imposed an explicit relation between the amplitude of the bispectrum and the frequency of the power spectrum. In the simplest form of axion monodromy, one has f NL = A log ω 2 log /8. On account of the quadratic scaling with frequency and the fact that the constraint on the amplitude tends to become poorer as the frequency increases due to projection (see, e.g., Fig. 30), it was already pointed out in Planck Collaboration XX (2016) that there is no evidence for this relation in the data, and with the current data set this situation remains unchanged.

Background and modeling
Single-field inflation with a canonical kinetic term gives rise to primordial super-Hubble comoving curvature perturbations, R. In this case the relative number densities of the various particle species are spatially constant, i.e., the perturbations are adiabatic. Typically photons are chosen as a reference species. Then adiabaticity implies that for every particle species with number density n i the quantity δ(n i /n γ ) vanishes. However, in addition to R, multi-field inflation can stimulate isocurvature modes, I i , where at primordial times n i /n γ varies spatially (Linde 1985;Polarski & Starobinsky 1994;Linde & Mukhanov 1997;García-Bellido & Wands 1996). In this section we consider all possible non-decaying modes of this type (Bucher et al. 2000): cold dark matter density isocurvature (CDI); baryon density isocurvature (BDI); and neutrino density isocurvature (NDI) modes. For completeness, we also constrain the fourth nondecaying mode, neutrino velocity isocurvature (NVI), although there are no known mechanisms to excite it. Finally, we consider compensated isocurvature perturbations (CIP) between baryons and CDM (Grin et al. 2011a,b). In this case, opposite BDI and CDI perturbations cancel in such a way that the total matter isocurvature perturbation vanishes and there is no first-order isocurvature signal in the CMB. However, we utilize a higherorder lensing-like effect from this mode to obtain constraints on CIP from Planck temperature and polarization power spectra. We find the most powerful power-spectra-based constraints on this mode by employing the very low-L lensing potential reconstruction in Sect. 9.5, but leave the use of Planck trispectra in constraining CIP for future work.
As the positions of the peaks and dips of the CMB angular power spectra in the density isocurvature models are roughly in opposite phase compared to the pure adiabatic (ADI) spectrum, the primordial CDI, BDI, and NDI modes leave a very distinctive observational imprint on the CMB, whereas the imprint of the NVI mode more closely resembles the pure ADI mode; see, e.g., figure 43 in PCI15. Prior to the detection of CMB anisotropies, studies such as  and Efstathiou & Bond (1986, 1987 discussed the possibility that isocurvature perturbations were the sole source of cosmological fluctuations. However, at least after the detection of the first acoustic peak in T T , it became clear that the density isocurvature mode(s) had to be subdominant (Enqvist et al. 2000, while the adiabatic mode led to a good agreement with observations. Several pre-Planck isocurvature constraints were obtained (Stompor et al. 1996;Pierpaoli et al. 1999;Langlois & Riazuelo 2000;Amendola et al. 2002;Peiris et al. 2003;Valiviita & Muhonen 2003;Bucher et al. 2004;Moodley et al. 2004;Beltran et al. 2004;Kurki-Suonio et al. 2005;Dunkley et al. 2005;Bean et al. 2006;Trotta 2007b;Keskitalo et al. 2007;Komatsu et al. 2009;Valiviita & Giannantonio 2009).
The mixture of curvature and isocurvature perturbations can be uncorrelated, but typically an arbitrary amount of correlation arises between them if the trajectory in field space is curved between Hubble radius exit and the end of multi-field inflation (Gordon et al. 2001). In extreme cases, such as the simplest curvaton models, there is full correlation or full anticorrelation between R and I. In the following subsections, we start with the generic case of generally correlated adiabatic and CDI, NDI, or NVI perturbations. Then we deal with various special CDI (or BDI) cases with no correlation or full (anti)correlation.
We parameterize the primordial perturbations as in PCI15, following the notation described there. The primary perturba-  Table 14. Constraints on mixed adiabatic and isocurvature models. We report 95 % CL intervals or upper bounds on the isocurvature fraction β iso at three scales (k low = 0.002 Mpc −1 , k mid = 0.050 Mpc −1 , and k high = 0.100 Mpc −1 ), the scale-independent correlation fraction, cos ∆, and the non-adiabatic contribution to the CMB temperature variance, α non-adi . Here ∆χ 2 is the difference between the χ 2 of the best-fit mixed and pure adiabatic models. In the last column we give the difference between the log of Bayesian evidences.  Fig. 37. Constraints on the primordial isocurvature fraction, β iso , at k low = 0.002 Mpc −1 and k high = 0.100 Mpc −1 ; the primordial correlation fraction, cos ∆; the isocurvature spectral index, n II ; and the correlation spectral index, n RI = (n RR + n II )/2, for the generally correlated mixed ADI+CDI model (a), for the ADI+NDI model (b), and for the ADI+NVI model (c). All these parameters are derived, and the distributions shown here result from a uniform prior on the primary parameters shown in Fig. 36. However, the effect of the non-flat derived-parameter priors is negligible for all parameters except for n II (and n RI ) where the prior biases the distribution toward unity. Note that these spectral indices are not well constrained, since we do not have a detection of non-zero isocurvature or correlation amplitude. With a sufficiently small isocurvature or correlation amplitude, an arbitrarily small or large spectral index leads to a very good fit to the data, since the model is then practically adiabatic over the range covered by the Planck data.
tion parameters scanned by MultiNest (in addition to the four standard ΛCDM background cosmological parameters and the Planck nuisance parameters) are the primordial abiabatic perturbation power and isocurvature perturbation power at two scales, corresponding to k 1 = k low = 0.002 Mpc −1 and k 2 = k high = 0.1 Mpc −1 , namely, P (1) RR , P (2) RR , P (1) II , P (2) II , and the correlation power between R and I at k 1 , i.e., P (1) RI . We assume a powerlaw form for the adiabatic and isocurvature power spectra and denote the spectral indices that can be calculated from the primary parameters by n RR and n II . The correlation spectrum is also assumed to obey a power law, with spectral index n RI = (n RR + n II )/2. Thus P (2) RI is not an independent parameter. This ensures that the correlation fraction cos ∆ = P RI /(P RR P II ) 1/2 stays inside the interval (−1, 1) at every k, as long as we reject any P (1) RI which does not obey this requirement. While the correlation fraction is k-independent in our modelling, the primordial isocurvature fraction β iso (k) = P II (k)/ [P RR (k) + P II (k)] depends on k, unless n II = n RR . We also report β iso at an intermediate scale, k mid = 0.05 Mpc −1 .
We do not separately quote constraints on BDI or the total matter density isocurvature (MDI), since these modes are observationally indistinguishable from the CDI case. 13 13 If we assume no NVI or NDI perturbations, then the MDI perturbation (i.e., the spatial perturbation in the relative number densities of Numerical results for various isocurvature models and selected derived parameters are reported in Table 14, utilizing various data combinations. The table is divided into three main sections: generally correlated models (discussed in Sect. 9.2); oneisocurvature-parameter CDI models (discussed in Sects. 9.4.1 and 9.4.2); and, finally, two-isocurvature-parameter CDI models (discussed in Sects. 9.4.3, 9.4.4, and 9.4.5). For generally correlated CDI we study the stability of constraints (see Sect. 9.3) by using several different subsets of the Planck data: (1) only high-TT; (2) high-TT+lensing; (3) TT,TE,EE; and (4) TT,TE,EE+lensing. For comparison, some Planck 2015 and WMAP results are also cited. Table 14 also includes comparmatter particles and photons) is As we will see, the posteriors for Ω c h 2 and Ω b h 2 are insensitive to the assumed initial conditions. Thus it is a good approximation to use the mean values obtained in the generally correlated mixed adiabatic and CDI model with TT,TE,EE+lowE+lensing data, namely Ω c /Ω m 0.842, Ω b /Ω m 0.158, Ω c /Ω b 5.33, and (Ω c /Ω b ) 2 28.4. For example, to convert our CDI upper bound on P II to a BDI bound, we should multiply the constraint by (Ω c /Ω b ) 2 = 28.4, and to convert the CDI P RI to BDI, we should multiply the constraint by Ω c /Ω b 5.33. If β iso 1, then this also can be converted to a BDI constraint by multiplying the CDI constraint by 28.4. The constraint on cos ∆ will be the same for the CDI and BDI cases, since the conversion factor cancels out.  Fig. 38. Posterior probability density of the observable nonadiabatic fraction of the CMB temperature variance, assuming a generally correlated mixed adiabatic and isocurvature model. These results used Planck TT+lowE+lensing data (dashed lines) and TT,TE,EE+lowE+lensing data (solid lines).
isons to the pure adiabatic model in terms of the difference in the best-fit χ 2 and the natural logarithm of the Bayesian evidence ("model probability") ratios ln B, negative ln B being evidence against the mixed models. 14

Results for generally correlated adiabatic and isocurvature modes
This subsection explores mixed adiabatic and isocurvature models where only one isocurvature mode at a time is considered. We consider the CDI, NDI, and NVI modes using the Planck 2018 TT(,TE,EE)+lowE(+lensing) data. All five primordial perturbation power amplitudes (of which three describe the isocurvature perturbations) are free parameters. It follows that n II and n RR are independent and cos ∆ varies between −1 and +1. The constraints for the primary perturbation parameters and the derived parameter P (2) RI are shown in Fig. 36. In all three cases the Planck TT+lowE+lensing results are very similar to the previous results from the Planck 2015 TT+lowP+lensing likelihood. As expected, the lower value of τ preferred by the 2018 (lowE) data is reflected in the adiabatic amplitudes P (1) RR and P (2) RR . For CDI and NDI there is no significant shift in the constraints on isocurvature parameters, but we find slightly tighter constraints than in 2015. For NVI, a minor shift towards more negative correlations is observed (see the last two panels of Fig. 36c). As in 2015, adding the high-TE,EE data significantly tightens the constraints in all three cases.
When fitting the generally correlated three-isocurvatureparameter models, the Planck data are consistent with null detection, i.e., with the pure adiabatic model, P (1) II , P (2) II = (0, 0) and P (1) RI = 0 (and P (2) RI = 0). The natural logarithm of the ratio of model probabilities [i.e., the Bayes factor ln B = ln(P ISO /P ADI )] is below −10.9, corresponding to odds of less than 1 : 54 000 for all three (CDI, NDI, NVI) models. If there were an undetected subdominant isocurvature contribution to the primordial perturbations, a negative correlation between R and I would be favoured, in particular for NDI and NVI (see the last two panels of Figs. 36a,b,c). With our sign convention, this leads to a neg- 14 The values of ln B depend on the priors. We adopt uniform priors in the range (15, 40) × 10 −10 for the adiabatic, (0, 100) × 10 −10 for the isocurvature, and (−100, 100) × 10 −10 for the primordial correlation power parameter. ative contribution to the Sachs-Wolfe effect and hence reduces the amplitude of the temperature angular power spectrum at low multipoles. Figure 37 updates the 2015 Planck constraints on the derived primordial fractions and spectral indices. At large scales we find with Planck TT,TE,EE+lowE+lensing that β iso (k low ) < 2.5 % for the CDI, 7.4 % for the NDI, and 6.8 % for the NVI model, all at 95 % CL. Figure 38 shows the non-adiabatic fraction in the observed CMB temperature variance, defined as where The non-adiabatic fraction |α non−adi | is below 1.7 % with Planck TT,TE,EE+lowE+lensing data for all three cases at 95 % CL.
Since the Planck data do not allow a significant isocurvature contribution, the determination of standard cosmological parameters depends only very weakly on the assumed initial conditions, as seen in Fig. 39. We place this result in historical perspective in Fig. 40 (and Table 14) where the parameter determinations of the mixed CDI model and the pure adiabatic model are compared to the pre-Planck constraints set by the WMAP 9year data. 15 Planck has dramatically tightened the constraint on the adiabatic spectral index. Its value is now 8.4σ below unity (scale-invariance) in the pure ADI case. Allowing for generally correlated CDI reduces the significance of this detection only slightly, to 7σ, whereas the WMAP 9-year data were consistent with a blue tilt as large as n RR = 1.06 at 95 % CL. The non-adiabatic contribution to the CMB temperature variance is constrained (about zero) 5 times more tightly than by WMAP. Finally, the allowed range for the sound horizon angle has shrunk by a factor of 10 in the CDI case, thanks to the Planck data covering more acoustic peaks beyond the first three peaks detected by WMAP.

Role of lensing parameter A L and likelihood choices
The small-scale primordial CDI amplitude is extremely sensitive to the details of the high-temperature and polarization power spectra and to choices made in constructing the likelihoods. Therefore the general CDI model serves as a robustness test of the Planck data and likelihoods. We now discuss a few curious aspects related to CMB lensing and likelihoods.
Lensing smooths the peaks of the CMB power spectra. This effect is taken into account in our theoretical predictions for the mixed adiabatic and isocurvature models by first calculating the total unlensed CMB spectra as a sum of adiabatic, isocurvature, and correlation C 's, and performing a similar summation for the lensing potential power spectrum (Seljak 1996;Lewis & Challinor 2006). The total lensing potential is then used to lens the total CMB spectra. Starting with the WMAP data, accounting for CMB lensing became necessary for calculating constraints on isocurvature models, as Valiviita et al. (2012) showed that there is a strong degeneracy between the lensing effect and  the CDI contribution in the generally correlated mixed models. Fixing n II = 1 or n II = n RR (as is done in the next subsection) makes this degeneracy disappear. This is because in these models the CDI contribution modifies only the low-part of the angular power spectra. The transfer function mapping the primordial CDI mode to the T T (and EE) angular power is suppressed by a factor (k/k eq ) −2 ∼ ( / eq ) −2 relative to the adiabatic mode. Therefore, to be observable at high , the CDI mode must be blue tilted (n II > 1). A blue-tilted CDI mode affects the total angular power spectra in a manner somewhat similar to lensing. Since the acoustic peaks of the CDI mode have the opposite phase compared to the adiabatic mode, a CDI admixture can "smooth" the peaks and dips of adiabatic acoustic oscillations. The NDI mode does not have precisely the opposite phase and is not damped relative to the adiabatic mode (see figure 43 in PCI15). Thus we expect a weaker impact of lensing on the primordial NDI amplitude than in the CDI case. Therefore in this subsection we explore the general CDI model as an example.
Starting with the Planck 2013 release, the consistency of the smoothing effect with the adiabatic ΛCDM model has been rou-tinely tested by multiplying the lensing power spectrum by a phenomenological lensing consistency parameter, A L , prior to lensing the unlensed CMB spectra (Calabrese et al. 2008). The expectation is that A L = 1. However, the Planck temperature and polarization data prefer a higher level of lensing-like smoothing (A L > 1) than expected in the adiabatic ΛCDM model. In the 2018 release (PCP18) we have Adding the lensing reconstruction data pulls these constraints towards A L = 1 (see also Table 15). The measurement of A L when TT,TE,EE data are included depends on the calibration of the polarization channels. This procedure and the details of the sky masks differ between the Planck baseline Plik TT,TE,EE and the alternative Planck CamSpec TT,TE,EE likelihood, as discussed in PPL18 and PCP18. The Planck CamSpec TT,TE,EE likelihood prefers a smaller value of A L than Plik, but still lying about 2σ above unity.  Fig. 37a (red curves), except now without the lensing likelihood, which to some extent hides the differences between other likelihoods and the effect of A L . Black solid contours show the results using the same likelihood as in the reference case, but now for the mixed adiabatic and CDI model when simultaneously allowing A L to vary. The remaining two curves are for the mixed adiabatic and CDI model (with A L = 1), but now changing the low-likelihood from Commander TT+SimAll EE to LFI 70 GHz T,E,B (grey), or high-likelihood from Plik to CamSpec (blue).
Given the above motivation, we check the response of the generally correlated CDI model to the various possible choices of likelihoods available in the Planck 2018 release, and, on the other hand, we gauge how the baseline Plik likelihood reacts when allowing A L to vary. For clarity, in Fig. 41 we restrict the analysis to high-TT,TE,EE and low-TT,EE(,BB) data without the lensing reconstruction data, but in Table 14 we also report TT+lowE and TT+lowP results, and include lensing in some cases.
We notice a considerable variation in the constraints on the isocurvature power at high k, P (2) II , which corresponds to the high-region in the observed power spectra. In Fig. 41 compare the red dashed reference contours (obtained with the Planck baseline Plik TT,TE,EE+lowE likelihood) for the CDI model (where A L = 1) with the solid black contours (obtained with the same data) for the CDI+A L model, where A L is allowed to vary. In many cases, adding an extra free parameter is expected to weaken the constraints on the other parameters, but in this case adding A L tightens the 95 % CL constraint on P (2) II by a factor of 2.5 from 28.6 × 10 −10 to 11.4 × 10 −10 . This is reflected in the derived primordial isocurvature fraction β iso (k high ), whose upper bound changes from 0.58 to 0.36, according to Table 14. Therefore, we conclude that, when A L = 1, the CDI mode par-tially accounts for the extra lensing-like smoothing effect required by the Planck TT(,TE,EE) data. Once we allow the lensing amplitude to vary, there is not much need for the CDI contribution at high , which should be kept in mind when interpreting the results. In Table 14 we report four cases where A L is allowed to vary. In all the other cases we have fixed A L = 1. In these cases the constraints at high k are "conservative", i.e., weaker than the Planck data were expected to be capable of (Finelli et al. 2018), due to CDI (or NDI) partially fitting the lensing anomaly. Furthermore, again comparing the red dashed and black solid contours, we observe a slight weakening of the constraint for P (1) II in the CDI+A L model. This is due to the rigidity of the assumed power-law spectrum. When the high-data allow much less CDI, the low-(low-k) CDI amplitude can be larger without much affecting the middle-range of the CMB power spectra between the first and third acoustic peaks, which is the most sensitive region to departures from adiabaticity. (This is the same see-saw effect discussed in the end of Sect. 3.6 in the case of tensor perturbations.) Finally, in the P (1) RR , P (2) RR panel we see a minor shift toward smaller amplitudes, which is an indication of the well-known degeneracy between A L and the overall primordial perturbation amplitude. From Table 14 it is obvious that adding the lensing data reduces the differences discussed above. This is again as expected, since the lensing data favour values of A L only mildly above unity. For example, with Planck TT,TE,EE+lowE+lensing we obtain β iso (k high ) < 0.49 for CDI+A L and 0.47 for the CDI model.
We now proceed to a comparison of the likelihoods. The grey shaded contours in Fig. 41 indicate the results for the same CDI model as the red dashed contours, but changing the low-likelihood from the combination Commander TT+SimAll EE to the LFI 70 GHz pixel-based low T,E,B, which is by its methodology and construction very similar to the 2015 baseline lowlikelihood (dotted contours). Indeed, this can be seen in the results: most of the isocurvature parameters follow more closely the 2015 results with this likelihood combination than with the 2018 baseline. This implies that when it comes to isocurvature, not much has changed in high-TT. 2018 lowP favours slightly smaller values of the optical depth τ than the 2015 version, hence the small shift towards smaller values of the adiabatic amplitude in the P (1) RR , P (2) RR panel. With respect to the red dashed contours, the grey contours prefer higher adiabatic amplitudes and have a long degeneracy line in the P (1) RR , P (2) RR plane. This is due to lowP having a higher central value and larger uncertainty on τ.
Finally, the blue shaded contours in Fig. 41 represent the results when using the CamSpec likelihood, to be compared to the red dashed contours obtained by the baseline Plik likelihood. All the other parameters shown are relatively stable against the high-likelihood, but P (2) II stands out. CamSpec leads to an upper bound of 12.6 × 10 −10 , whereas the baseline Plik result was 28.6 × 10 −10 , or for β iso (k high ) 0.38 versus 0.58 at 95 % CL, ac-cording to Table 14. This difference is not surprising, given the different responses of these likelihoods to A L in the adiabatic ΛCDM+A L model, and keeping in mind the A L -CDI degeneracy in the CDI model. However, this difference is not as concerning as it might appear at first sight: all the cases shown in Fig. 41 are fully consistent with zero isocurvature. It is only the upper bound that varies, with the baseline Plik likelihood and CDI model with A L = 1 leading to the most conservative (i.e., weakest or safest) upper bounds.
For CDI the 2015 release Planck high-TT data favoured a negative correlation fraction but the preliminary high-TT,TE,EE data favoured a slightly positive correlation. This was confirmed using only the high-Plik likelihood (and a prior on τ), as shown by the red curves in the top panel of Fig. 42. Including the low-data (black curves) did not significantly alter this tension between TT and TT,TE,EE results. In the present 2018 Planck release this tension has disappeared. Both high-TT and TT,TE,EE data lead to a correlation fraction posterior peaking at zero, as demonstrated in the bottom panel of Fig. 42. Including the low-TT data (black dashed curve) still shifts the posterior slightly towards negative values, due to the low T T power at low multipoles in the data.

Specific CDI models
In this subsection we constrain CDI models with only one or two isocurvature parameters. The two-parameter cases were not studied in the 2013 and 2015 Planck releases.
First we fix n II to unity and assume no correlation between the CDI and adiabatic modes ("axion"), or we fix n II = n RR and assume full (anti)correlation between the CDI and adiabatic modes ("curvaton I/II"). These models are less sensitive to any residual systematic effects in the high-data (such as the determination of polarization efficiencies or foreground modeling) than the generally correlated models, since CDI now modifies the angular power spectra insignificantly at 200 (see figure 43 in PCI15). As seen in the middle section of Table 14, the Bayesian evidence values for the one-parameter extensions of the adiabatic ΛCDM model are higher than for the threeparameter extensions, but all Bayes factors fall below −5. None of the one-parameter extensions improve χ 2 over the adiabatic ΛCDM model. The two-parameter extensions in the bottom section of Table 14 are even more strongly disfavoured, except for the uncorrelated case with free n II ("axion II"), which is actually the only model that improves the best-fit χ 2 by slightly more than the number of extra parameters.

Uncorrelated ADI+CDI ("axion I")
Particularly insensitive to any 30 data is the "axion I" case, since the CDI transfer function has a (k/k eq ) −2 suppression and there is no correlation component whose amplitude would be higher than that of the isocurvature alone and hence would modify the adiabatic spectrum. The "axion I" case is achieved in our parameterization by setting P RI = 0 and P (2) II = P (1) II . Thus the only varied isocurvature parameter is P (1) II . This uncorrelated case with n II = 1 is a good approximation for many multi-field inflationary models where the slow-roll parameter (in the isocurvature field perturbation direction) η ss is negligible and the background trajectory in field space is straight between Hubble radius exit and the end of inflation. The predictions for the spectral indices (to first order in the slow-roll parameters) are n RR = 1 − 6 + 2η σσ and n II = 1 − 2 + 2η ss , where 50 ≥ 0 and η σσ is the second slow-roll parameter in the "adiabatic" direction (i.e., along the trajectory) in the field space. (An exact match with our model would require η ss = .) The axion model [see, e.g., a recent review by Marsh (2016) and references therein], which was originally proposed to solve the strong CP problem and provides a dark matter candidate, can produce this type of isocurvature modes with n II 1 under the following assumptions (PCI13; PCI15): the Peccei-Quinn symmetry should be broken before inflation; it should not be restored by quantum fluctuations of the inflaton nor by thermal fluctuations when the Universe reheats; and axions produced through the misalignment angle should form a significant fraction of the dark matter. Table 14 indicates a slight tightening of the "axion I" constraints using TT+lowE+lensing with respect to 2015 TT+lowP+lensing. This is due to the change of the baseline lowdata from the 2015 LFI 70-GHz pixel-based T,E,B to the 2018 combination of Commander TT and SimAll EE, which in the generally correlated cases also gave tighter constraints at low k. As expected, the addition of high-TE,EE data only marginally improves the constraints, since the standard (non-isocurvature) parameters are better constrained now. For n II = 1 uncorrelated CDI, we obtain β iso (k mid ) < 0.038 0 ≤ α non−adi < 1.55 % (95 % CL, Planck TT,TE,EE +lowE+lensing).
Using equation (73) of PCI13, we convert the constraint on the primordial isocurvature fraction to a bound on the inflationary energy scale. If all the dark matter is in axions, the above β iso (k mid ) constraint corresponds to the same limit we quoted in 2015, that is, H inf < 0.86 × 10 7 GeV f a 10 11 GeV 0.408 (95 % CL) , (89) where H inf is the expansion rate at Hubble radius exit of the scale corresponding to k mid and f a is the Peccei-Quinn symmetrybreaking energy scale.

Fully (anti)correlated ADI+CDI ("curvatonI/II")
If n II = n RR , the low-data are maximally sensitive to the fully correlated isocurvature perturbations. In this case the correlation component is a geometric average of the adiabatic and isocurvature components, and hence much larger than the isocurvature component alone. We achieve this case in our parameterization by setting P (2) II = P (2) RR /P (1) RR P (1) II and P (1) RI = ± P (1) RR P (1) II 1/2 , i.e., cos ∆ = ±1. The only isocurvature parameter to be varied is again P (1) II . Since n II = n RR , the derived isocurvature fraction β iso is independent of k. A physically motivated example of this type of model is the simplest curvaton model, where a light scalar field χ that is subdominant (and hence irrelevant for the inflationary dynamics) starts to oscillate at the bottom of its potential after the end of inflation, causing its average energy density to evolve like non-relativistic matter. Once fully (or almost fully) dominating the energy density of the Universe, this curvaton field decays either to CDM or to other species (Mollerach 1990;Linde & Mukhanov 1997;Enqvist & Sloth 2002;Moroi & Takahashi 2001;Lyth & Wands 2002;Bartolo & Liddle 2002;Lyth et al. 2003). The amount of isocurvature and non-Gaussianity present after curvaton decay depends on the "curvaton decay fraction," r D = 3ρ χ /(3ρ χ + 4ρ radiation ), evaluated at curvaton decay time. Under a number of (very) restrictive assumptions discussed in PCI15, the curvaton model can lead to fully (anti)correlated CDI (or BDI) and adiabatic perturbations.
Not surprisingly, both in the fully correlated and anticorrelated cases, the constraint on β iso is much (about 40 times) stronger than in the uncorrelated case. At 95 % CL, Planck TT,TE,EE+lowE+lensing leads to β iso < 0.00095 (Fully correlated, "curvaton I"), (90) β iso < 0.00107 (Fully anticorrelated, "curvaton II"), both rounded to 0.001 in Table 14. As in 2015, the TT data favour anticorrelation, due to the low power in the low-temperature compared to the expectation of the adiabatic ΛCDM model. But when the TE,EE data (which do not particularly favour negative correlation) are added, a very tight (one part per thousand) constraint on the primordial isocurvature fraction results.
Fully correlated perturbations are obtained, e.g., in case 4 described in Gordon & Lewis (2003). Many models giving anticorrelation produce too large an isocurvature fraction to be consistent with the above limit, but case 9 of Gordon & Lewis (2003) survives. After the curvaton decay, the primordial isocurvature fraction in these models will be β iso 9(1 −r) 2 /[r 2 D + 9(1 −r) 2 ], wherer = r D for the fully correlated CDI case andr = r D /R c ≥ 1 for the fully anticorrelated CDI case, and R c =ρ c /(ρ c +ρ b ) is the CDM fraction of the total non-relativistic matter.
On the other hand, the nonlinearity parameter describing non-Gaussianity is (Sasaki et al. 2006) where ∆ 2 s = δχ 2 s /χ 2 is the small-scale variance of the curvaton perturbations, or the ratio of the energy density carried by the curvaton particles to the energy density of the curvaton field (if there is significant production of curvaton particles). The parameter f local NL cannot be smaller than −5/4, which is obtained when r D = 1 and ∆ 2 s = 0, as implicitly assumed, for example, in Bartolo et al. (2004a,b). The above β iso limits correspond to the following r D and f local NL constraints (assuming ∆ 2 s = 0): 16 0.98982 < r D ≤ 1 ("curvaton I") 0.84347 ≤ r D < 0.85129 ("curvaton II") Even with the maximal allowed isocurvature fraction, the local non-Gaussianity in the curvaton model is well within the observational Planck limits presented in Planck Collaboration IX (2018). The residual isocurvature peturbations in the two studied curvaton models set much tighter constraints on the curvaton decay fraction than do constraints on the observed (consistent with zero) non-Gaussianity.
9.4.3. Uncorrelated ADI+CDI with free n II ("axion II") Axion models do not necessarily produce nearly scale-invariant isocurvature perturbations. In particular, even highly bluetilted spectra (in the observable CMB range) are possible. For example, Kasuya & Kawasaki (2009) construct a model with n II = 2-4. This motivates studying a two-isocurvatureparameter model, where adiabatic and isocurvature modes are uncorrelated, but the isocurvature fraction and spectral index are free to vary. In our parameterization this is achieved by setting P (1) RI = 0, and varying P (1) II and P (2) II independently. The results for this model are presented in the first two rows of the third section of Table 14. The low-temperature data do not favour any extra contribution beyond the (already too high) abiabatic contribution, whereas the fit to the high-temperature and polarization data can be improved slightly by the "smoothing" caused by the CDI mode. This leads to a very blue isocurvature spectrum. Planck TT,TE,EE+lowE+lensing gives at 95 % CL 1.55 < n II < 3.67, consistent with the recent findings of Chung & Upadhye (2017). Even the very large upper bound β iso (k high ) < 77 % corresponds to a contribution of less than order 1% to the observable CMB T T (or EE) power spectra at 1400. The uncertainty in the Planck T T spectrum at these high multipoles is ∆D T T ∼ 10 µK 2 and the actual spectrum is D T T ∼ 1000 µK 2 . Thus the allowed CDI contribution is only of the same 1 % order as the observable uncertainty. Consequently the non-adiabatic contribution to the observed CMB temperature variance, α non−adi , is also vanishingly small, between 7×10 −4 and 7 × 10 −3 .

("curvaton III")
Apart from the extremes of ±100 % correlation, some curvaton models predict an arbitrary degree of correlation. The generic feature of most curvaton models is that the isocurvature and adiabatic spectral indices are equal. This is because both perturbations typically arise from the same source. In the nextto-simplest models, the correlation fraction can be written as cos ∆ = √ λ/(1 + λ), where λ = (8/9)r 2 D * (M Pl /χ * ) 2 . Therefore, the model is fully correlated only if λ 1, in which case the results of "curvaton I" apply. If the slow-roll parameter * is very close to zero or the curvaton field valueχ * is large compared to the Planck mass, this model leads to almost uncorrelated perturbations and the constraints are well approximated by "axion I." Any other case leads to an arbitrary degree of positive correlation between the CDI and adiabatic modes.
Modulated reheating with thermal or non-thermal production of gravitinos can lead to positive or negative correlation, respectively (Takahashi et al. 2009). While the correlation could in principle be arbitrarily large, the observational constrains on β iso favour only small correlations.
Arbitrarily correlated ADI+CDI with n II = n RR is also a good approximation for those two-field (or multi-field) slowroll models [e.g., double quadratic inflation (Langlois 1999;Beltrán et al. 2005)] where the trajectory in field space is curved between the Hubble radius exit of perturbations during inflation and the end of inflation. The fraction of isocurvature perturbations converted to adiabatic depends on how the trajectory is curved and this part of the adiabatic perturbations will be fully (anti)correlated with the isocurvature modes, whereas the adiabatic perturbations already present at Hubble radius exit are uncorrelated with isocurvature modes to first order in the slow-roll parameters, and only slightly correlated to second or-der. [See, e.g., Gordon et al. (2001), Amendola et al. (2002), van Tent (2004, and Byrnes & Wands (2006).] The result is a non-zero correlation between isocurvature and total adiabatic perturbations. The spectral indices of both components are typically 1 − O(slow-roll parameters), which is well approximated by n II = n RR since the data indicate n RR 0.965.
As expected, the Planck data favour negative correlations, since these n II = n RR models modify only the low-part of the CMB spectra, where T T power is lower than predicted by the adiabatic ΛCDM model. With TT,TE,EE+lowE+lensing we find, at 95 % CL, β iso < 0.039 and −0.41 < cos ∆ < 0.31. 9.4.5. Fully (anti)correlated ADI+CDI with free n II The remaining two-parameter CDI extensions of the adiabatic ΛCDM model are those where the perturbations are fully (anti)correlated, as in the simplest curvaton models, but the isocurvature spectral index is not fixed to the adiabatic one. In this case the free isocurvature parameters are P (1) II and P (2) II , while P (1) RI = ± P (1) RR P (1) . These models are somewhat difficult to motivate, since full (anti)correlation typically implies that the curvature and isocurvature perturbations have their origin in (the decay products of) the same field. Then one would expect equal spectral indices, as in the curvaton model. The conversion of isocurvature perturbations to adiabatic ones (e.g., between Hubble radius exit and the end of inflation, or by curvatontype decay, or by reheating/thermalization) should be scale dependent in order to obtain n II n RR . Slow-roll two-field inflation leads to an exact match, n II = n RR , in the case where cos 2 ∆ = 1. [See, e.g., Byrnes & Wands (2006)]. Nevertheless, for completeness we report constraints on these phenomenological models in the last four rows of Table 14. Since the low-TT data favour negative correlation, a larger isocurvature fraction is allowed in the fully anticorrelated case at low k. This leads to scale-invariant isocurvature perturbations being in the favoured region of parameter space, namely −0.28 < n II < 1.86 with TT,TE,EE+lowE+lensing at 95 % CL. In contrast, in the fully correlated case the low-TT data disfavour any isocurvature contribution, and hence prefer a blue spectrum, with 1.37 < n II < 3.65.

Compensated BDI-CDI mode
This subsection presents constraints on uncorrelated adiabatic and scale-invariant CIP modes and discusses the strong degeneracy between the phenomenological lensing parameter A L and the CIP amplitude (Valiviita 2017). Assuming that there are no NVI or NDI perturbations, the total matter density isocurvature perturbation I MDI , given by Eq. (83), vanishes if This mode, where the anticorrelated CDI and BDI perturbations cancel even though their individual amplitudes can be large, is called a compensated baryon and cold dark matter isocurvature mode. The CIP mode does not leave a linear-order isocurvature signal in the CMB or matter power spectra (Gordon & Lewis 2003), although it modifies the trispectrum (Grin et al. 2011a(Grin et al. , 2014. However, at the next order there is a smoothing effect on the high-T T , T E, and EE spectra. A formal derivation can be found, for example, in Smith et al. (2017). Here we summarize the heuristic arguments of Muñoz et al. (2016).   , and conservative Planck 2018 lensing data (black points in grey boxes), along with the best-fit models to the Planck data: the best-fit adiabatic ΛCDM model to 2018 TT+lowE (green dashed line); the best-fit ΛCDM+A L model to 2018 TT+lowE (magenta solid line, A L = 1.26); and the best-fit ΛCDM+CIP model to 2018 TT,TE,EE+lowE and conservative lensing data (black solid line, ∆ 2 rms = 0.0036) and to 2015 TT,TE,EE+lowP and conservative lensing data (black dotted line, ∆ 2 rms = 0.0071). As CIP modifies only the very low-L part of the lensing power spectrum, the conservative 2015 lensing data (40 ≤ L ≤ 400) are insensitive to CIP even when ∆ 2 rms = 0.0071. On the other hand, the first two data points of the 2015 aggressive lensing data disfavour the large CIP amplitude (Smith et al. 2017), which gives a very good fit to all the other data. Planck 2018 conservative lensing data cover the range 8 ≤ L ≤ 400 and consequently disfavour CIP variances ∆ 2 rms 0.004.
On scales larger than the sound horizon, condition (95) is preserved until last scattering and can be written as Consequently, CIP can be described as a large scale modulation of the baryon and CDM density (Muñoz et (97) Here the overbar denotes an average over the whole sky and ∆(n) a small perturbation about this average in the directionn, as illustrated in Fig. 43. In patches of sky where the CMB photons originate from baryon-overdense regions, the odd acoustic peaks at high-are more pronounced relative to the even peaks compared to the patches where the photons originate from baryonunderdense regions. Averaging over the sky leads to a lensinglike smoothing of the high-peaks. A convenient measure of CIP is the variance ∆ 2 rms ≡ |∆(n)| 2 P I BDI I BDI . If ∆ is a Gaussian random variable, the observed angular power of T T , T E, or EE will be where Ω b (∆) = (1 + ∆)Ω b and Ω c (∆) =Ω c −Ω b ∆. For brevity, we will denote the power spectrum in the integrand by C | ∆=δ . For each δ it can be calculated by assuming adiabatic initial conditions. Approximating the integrand by the first three terms of 53 Planck Collaboration: Constraints on Inflation its Taylor series about ∆ = 0, we end up with In the following we describe parameter scans where we vary the six standard (adiabatic) ΛCDM parameters, the Planck nuisance parameters, and the CIP variance ∆ 2 rms , calling this oneparameter extension of the ΛCDM model the "ΛCDM+CIP" model. We evaluate the right-hand side of Eq. (99) at each point in parameter space using a finite-difference approximation for the second derivative: where δ should be "sufficiently small." In practice, good numerical accuracy is achieved if δ is of order ∆ 2 rms . So at each point in our MultiNest scan we set δ = ∆ 2 rms for the point currently under evaluation, and thus the result of Eq. (99) simplifies to With this method each angular power spectra evaluation takes twice as long as for the pure adiabatic case since the spectra are now an average of two spectra, resulting from different values of Ω b and Ω c . 17 Unlike the high-T T , T E, and EE spectra, the high-L lensing potential power spectrum is virtually unaffected by CIP. Instead, CIP modifies the low multipoles of [L(L + 1)] 2 C φφ L /(2π) by, approximately, adding a term ∆ 2 rms × (L/0.053) −2 . For details, see table II in Smith et al. (2017). As illustrated in Fig. 44, when using the Planck 2015 conservative lensing data (40 ≤ L ≤ 400) this term does not affect the results. In contrast, the Planck 2018 conservative lensing data also contain the range 8 ≤ L < 40 and thus CIP variances ∆ 2 rms 0.004 fit the first data point of the 2018 lensing power spectrum (8 ≤ L ≤ 400) worse than in ΛCDM. However, even in this case the joint fit of the ΛCDM+CIP model to the TT, TE, EE, and lensing data is better than that of the ΛCDM model, the improvement being of the same order as for the ΛCDM+A L model.
The top panel of Fig. 45 shows the A L -∆ 2 rms degeneracy in the ΛCDM+A L +CIP model and how it can be broken by the lensing data. The value A L = 1 provides a good fit to the TT+lowE data, if ∆ 2 rms 0.016, and to the TT,TE,EE+lowE data, if ∆ 2 rms 0.010. The ΛCDM+CIP model (where A L = 1) with ∆ 2 rms 0.008 provides a better simultaneous fit to the Planck 2015 TT,TE,EE and conservative lensing data (40 ≤ L ≤ 400) than does the ΛCDM+A L model. When using the Planck 2018 conservative lensing data (8 ≤ L ≤ 400), the best-fit value of ∆ 2 rms decreases to 0.0036. This is due to the extra term ∝ L −2 brought by CIP to the lensing power estimator, as discussed above and shown in Fig. 44.
Since Planck TT,TE,EE+lowE and lensing data can be fit well by A L = 1 in the CIP model, we show in the bottom panel of Fig. 45 the one-dimensional posterior of ∆ 2 rms in the ΛCDM+CIP model. A non-zero value of ∆ 2 rms is preferred at the 2.9σ (2.5σ) level by Planck 2018 (2015) TT+lowE(lowP) data and at the 2.6σ (1.8σ) level by the TT,TE,EE+lowE(lowP) data. Without 17 Actually, we call CAMB three times (i.e., also withΩ b andΩ c , or ∆ = 0), in order to obtain some auxiliary parameters such as σ 8 correctly, although we do not report them here.

Constraints
Best-fit Data and model 1000∆ Table 15. Comparison of ΛCDM+CIP, ΛCDM+A L , and ΛCDM+A L +CIP models with various Planck datasets, when using the baseline Plik likelihood at high . The first two columns ("Constraints") are the 68% CL ranges or the 95% CL upper bounds on 1000∆ 2 rms (highlighted in bold for ΛCDM+CIP) and A L . The remaining columns give the best-fit pameter values, and the difference of the best-fit χ 2 and the difference of the log of the Bayesian evidence with respect to the pure adiabatic ΛCDM model. A negative ∆χ 2 means that the quoted model fits the data better than ΛCDM, while a positive ln B means that the Bayesian model comparison favours the quoted model, when adopting the uniform priors 0 ≤ 1000∆ 2 rms < 75 and 0.3 ≤ A L ≤ 1.7.
lensing the 2018 data thus more strongly favour the non-zero CIP amplitude than the 2015 data, which is as we would expect, since the favoured A L value in the ΛCDM+A L model has also increased. When using 2018 TT,TE,EE+lowE and the 2018 conservative lensing data the significance decreases to 2.0σ, while switching to the aggressive lensing data (8 ≤ L ≤ 2048) leads to 2.1σ. The 68 % CL ranges of ∆ 2 rms in the ΛCDM+CIP model, obtained with the baseline high-Plik likelihood in combination of other Planck data, are highlighted in Table 15. Replacing Plik with CamSpec (in particular CamSpec TT,TE,EE) leads to somewhat lower values, and reduced significance above zero: 2.7σ, 1.9σ, 1.7σ, and 1.9σ, respectively.
In order to check that the preference for ∆ 2 rms > 0 or A L > 1 is not just a parameter-space volume effect upon marginalization over other parameters, we also report in Table 15 the difference of χ 2 between the best fit in extended models and the base adiabatic ΛCDM model. With all data sets, all three extended models lead to an improvement of χ 2 which clearly exceeds the number of extra parameters of the model (1 for ΛCDM+CIP and ΛCDM+A L , and 2 for ΛCDM+A L +CIP). Although the inclusion of lensing data reduces this improvement of fit, the ΛCDM+CIP model gives a rather impressive ∆χ 2 = −4 with Planck TT,TE,EE+lowE and aggressive lensing data.
Since we observe a moderate preference for a non-zero CIP amplitude, it might be tempting to "solve" the Planck lensing anomaly by using CIP. However, this explanation seems quite unlikely, since in our treatment the CIP and adiabatic perturbations should be uncorrelated with each other, whereas CDI and BDI should be fully anticorrelated (and have a few orders of magnitude larger amplitude than the adiabatic modes while keeping the perturbations nearly Gaussian). It is difficult to imagine a physical model that could lead to this situation. For example, some variants of curvaton model would naturally lead to anticorrelated CDI and BDI, but in these models there would be a correlation with the adiabatic mode too (Gordon & Lewis 2003;He et al. 2015). The above-studied compensated BDI-CDI mode falls into a similar category to NVI: it is an interesting theoretical setup, but a compelling early-Universe model for stimulating this mode has still to be discovered.
Nevertheless, the baseline Planck Plik TT,TE,EE+lowE plus conservative lensing result, ∆ 2 rms = 0.0037 +0.0016 −0.0021 , is fully compatible with current complementary observations, in particular, the WMAP 95 % CL trispectrum constraint, ∆ 2 rms 0.012 (Grin et al. 2014), and the upper bound, ∆ 2 rms 0.006, following from the direct measurements of the variation of the baryon fraction in galaxy clusters (Holder et al. 2010;Grin et al. 2014). It will be interesting to learn what other future CMB anisotropy (Abazajian et al. 2016;Valiviita 2017;Finelli et al. 2018) and complementary measurements, such as observations of the distribution of neutral hydrogen using 21 cm absorption lines (Gordon & Pritchard 2009), BAO (Soumagnac et al. 2016(Soumagnac et al. , 2018, or CMB spectral distortion anisotropies (Haga et al. 2018), will tell us about the possible contribution of CIP to the primordial perturbations.

Constraints on anisotropic models of inflation
In this section we will test specific physical models for statistical anisotropy in the primordial fluctuations. More phenomenological multipole-or map-space tests are performed in the companion paper, Planck Collaboration VII (2018). Here we update the results of the 2015 release (PCI15) with polarization and new temperature analyses. Incorporating polarization into these tests is particularly important, due to the mild statistical significance of temperature anomalies such as the dipolar asymmetry. Polarization offers the potential to confirm or refute a physical origin for such anomalies via the measurement of independent fluctuation modes. We perform such a new test with k-space dipolar modulation models. In cases such as quadrupolar asymmetry, where no detection has been claimed with temperature, polarization offers the prospect of tightening existing constraints.
Some asymmetry models predict a modification to the isotropic power spectra, in addition to a dipolar or quadrupolar asymmetry. In other words, for these models, as well as non-zero off-diagonal multipole covariance elements, we expect departures in the diagonal elements relative to the standard ΛCDM prediction. Therefore the isotropic spectra can provide independent tests of such models even using temperature data alone (Contreras et al. 2018). The curvaton dipole modulation model we examine in Sect. 10.1.1 exhibits this property, and can be constrained via its predictions for isotropic isocurvature power. Similarly, some versions of the quadrupolar modulation model we study in Sect. 10.2 modify the isotropic spectra via a monopole term. In both cases these isotropic constraints will be important in narrowing the viable parameter space.

Dipolar asymmetry
A dipolar temperature power asymmetry has long been observed at the largest scales in the CMB (Eriksen et al. 2004), although its statistical significance is not high and is subject to a posteriori (look-elsewhere) corrections (Bennett et al. 2011;Planck Collaboration XXIII 2014;Planck Collaboration XVI 2016). Nevertheless, its large-scale character suggests potential links with inflationary physics and various models have been proposed to explain it. In this subsection we examine several physical models for a dipolar modulation. Some models where a generic CDM density isocurvature (CDI) or tensor component is dipole modulated have already been ruled out due to their isotropic predictions (Contreras et al. 2018), so we do not consider these further here.

Curvaton model
First we update our 2015 study (PCI15) of a specific inflationary model for the dipolar asymmetry: namely, the modulated curvaton model of Erickcek et al. (2009). In that study we showed that that model could not explain the observed asymmetry. Here, we generalize the curvaton model to allow for a non-scale-invariant uncorrelated CDI component. In addition, we treat the power spectrum (isotropic) constraints in a fully unified way with the asymmetry likelihood. Finally, we incorporate polarization. The modulated curvaton model employs a gradient in a background curvaton field to explain the observed large-scale power asymmetry. The curvaton, via coupling κ, produces nearly scaleinvariant CDI fluctuations, as well as a fraction, ξ, of the adiabatic fluctuations. Both of these components will be modulated. Up to a sign, ξ is equal to the correlation parameter, and is also a measure of the amplitude of dipolar modulation. The isocurvature fraction can be written in terms of these two parameters as Full details of this model and our treatment of it can be found in Erickcek et al. (2009) and PCI15. Using the dipolar asymmetry estimator from PCI15 we find the posteriors for the dipolar modulation parameters κ and ξ; the results are presented in Fig. 46 (red contours). We see that a substantial amount of asymmetry (as measured by amplitude ξ) can be captured by the model. This preference for asymmetry simply means that the curvaton model can explain the well-known dipolar asymmetry in temperature. However, isocurvature constraints from the power spectra via Eq. (103), which we refer to as the isotropic constraints, can provide independent information (Contreras et al. 2018). This is also shown in Fig. 46, with the blue contours. Here we see that the asymmetry and isotropic posteriors only weakly overlap, and the independent isotropic data do not support the presence of asymmetry for this model. No evidence for asymmetry (i.e., no preference for ξ > 0) is present in the joint constraints, which treat the isotropic and asymmetry data as independent. In other words, we have no reason to prefer this model over base ΛCDM. The model can explain the well-known dipolar asymmetry: note the preference for ξ > 0 in the asymmetry constraint (red contours and curves). However, the modulation preferred by the asymmetry constraint is reduced substantially when the isotropic constraint (blue) is added (black). The asymmetry constraint here uses SMICA, while the isotropic constraint uses Planck TT,TE,EE+lowE+lensing. Resolution is reduced at very small κ due to the sampling in β iso .

Adiabatic models
In the presence of a sufficiently large bispectrum it is possible that a long-wavelength mode can induce a dipolar asymmetry in the two-point function across our observable volume, although such scenarios appear to require fine tuning (Byrnes et al. 2016b). Nevertheless, examples have been constructed which satisfy the Planck f NL constraints (Byrnes et al. 2016a). In this subsection we consider adiabatic models of this type, in which the isotropic power spectra agree with standard ΛCDM, while a scale-dependent dipolar asymmetry is present in the off-diagonal multipole covariance (Contreras et al. 2017(Contreras et al. , 2018. As proposed in Contreras et al. (2017), we fit the asymmetry model parameters to the temperature data and then use those parameters to predict the asymmetry in polarization. We then compare those predictions with the Planck polarization data as a test for a physical modulation. Importantly, a position-(or k-) space model for the modulation is needed for reliable polarization predictions-it is not enough to restrict considerations to multipole space (Contreras et al. 2017) 18 .
As discussed in detail in Contreras et al. (2017), we take a portion R lo (x) of the adiabatic primordial fluctuations to be spa-18 E.g., if we observe a dipolar modulation in T at, say, 5 % to = 65, there is no reason to expect a modulation of the same amplitude and to the same scale in E, due to the different T and E transfer functions (Contreras et al. 2017). tially linearly modulated according to where R lo (x) is statistically isotropic with power spectrum P lo R (k), A ≤ 1 andd are the amplitude and direction of modulation, respectively, and r LS is the comoving radius to last scattering. This leads, to a good approximation, to the total temperature or polarization multipole covariance C m m ≡ a m a * m (105) to first order in A. Here C is the usual ΛCDM anisotropy power spectrum; δC ≡ 2(C lo + C lo ), where C lo is the power spectrum calculated in the usual way from P lo R (k); ∆X M is the multipole decomposition of An ·d; and the ξ M m m coefficients couple to ± 1 via In principle the scale dependence of the asymmetry spectrum P lo R (k) is completely free, but here we take three phenomenological forms which are capable of producing a large-scale asymmetry with a small number of parameters. First, we consider a simple power-law modulation, where P 0 R (k) is the usual ΛCDM spectrum, and n lo s and k lo 0 are the tilt and pivot scale of the modulation. We consider only red asymmetry tilts with n lo s ≤ n s , and choose k lo 0 = 1.5×10 −4 Mpc −1 . We also consider a tanh model, defined according to This spectrum approaches that of ΛCDM on scales larger than k c , with a width determined by ∆ ln k. That is, scales well above the cutoff k c will be modulated with amplitude A, and scales below will be unmodulated. Finally, we consider a model with a linear gradient in the scalar tilt, n s , across our volume. In this case the asymmetry spectrum can be written as with modulation amplitude ∆n s . There will be an implicit dependence on the pivot scale k * for this model. Given the multipole covariance, Eq. (106), we can construct a maximum likelihood estimator for the modulation, ∆X M . In the noise-free, full-sky case this takes the form (Moss et al. 2011;Planck Collaboration XVI 2016) where the cosmic variance of the estimator is given by The modifications we use to deal with realistic skies are described in detail in Planck Collaboration XVI (2016) and Contreras et al. (2017).
To decide whether the polarization data support the modulation model or not, we consider the quantityÔ j0 , which is the ratio of the maximum likelihood for modulation model j to that of ΛCDM (Contreras et al. 2017). In Fig. 47 we plot for the three adiabatic models histograms ofÔ j0 calculated for 300 statistically isotropic polarization simulations (sharing the required T E correlation with the real T data) added to the Planck temperature data (red outlines). This indicates our expectation forÔ j0 for the scenario that the temperature asymmetry is due to a statistical fluctuation and not to a physical modulation. We also plot in Fig. 47 histograms for 300 polarization simulations modulated with the best-fit parameters from the Planck temperature data (black outlines), to represent the scenario that the asymmetry is due to a physical modulation. In both cases the polarization simulations contain realistic levels of noise for Planck. By comparing the isotropic and modulated histograms, we can see that the quantityÔ j0 can serve to distinguish the two scenarios, but only relatively weakly for Planck noise (Contreras et al. 2017). The blue lines indicate the values using the actual SMICA polarization data (the results for the other component-separation methods are similar). We see that for these models the data do not help to decide whether we have a physical modulation or not, with pvalues of 43 %, 30 %, and 57 % for the power-law, tanh, and n s gradient models, respectively, relative to the isotropic simulations.

Quadrupolar asymmetry
We will next explore models that predict a quadrupolar direction dependence in the primordial power spectrum. In PCI15 we found no evidence for such a modulation, but several inflationary models have been constructed which predict this effect. Therefore it is important to extend those results with the improved polarization data. We now attempt to reduce the effect of unresolved point sources using the bias-hardened estimator approach of Planck Collaboration XV (2016). In PCI15 we pointed out that some models of quadrupolar asymmetry predict a modification to the angular power spectra as well. Here we will account for such modifications in our analysis, increasing the constraining power of temperature data, in particular for tilted models with non-scale-invariant modulation spectral index.
We assume a modulation of the primordial comoving curvature power spectrum of the form which can be rewritten as Here g 2m (k) ≡ 8π 15 g(k) Y * 2m (d), with g 2m (k) satisfying g 2,−m (k) = (−1) m g * 2m (k). We parameterize the scale dependence of the modulation as g(k) = g * (k/k * ) q , with pivot scale k * = 0.05 Mpc −1 . For q 0, in addition to producing a quadrupolar modulation of the anisotropies, this model affects the CMB isotropic power spectra via the term g(k)/3 in Eq. (114). We therefore consider a joint constraint with the isotropic power spectra likelihood to improve constraints over the modulation alone.
As in PCI15, we obtain constraints on the modulation parameters by forming quadratic maximum-likelihood estimates,ĝ 2m , for the data and simulations. For this we use the componentseparated data and 300 simulations provided by NILC, SEVEM, SMICA, and COMMANDER. For brevity we only show the SMICA results. We can then compute a covariance G and likelihood as We then evaluate the marginalized (over the angles) posterior for g * . For the isotropic constraints we simply include the modulation parameters in a CosmoMC run using Planck TT,TE,EE+lowE data, and then evaluate the marginalized (over all other ΛCDM parameters) posterior for g * . For q > 0 the effect of g * on the isotropic spectra occurs mainly at high , and is highly degenerate with n s . This degeneracy leads to slightly less stringent constraints than what one would achieve with a fixed n s . We show the marginalized posteriors for this case in the top panels of Fig. 48, where we see that the isotropic constraints are roughly comparable in strength to (and fully consistent with) the constraints from the asymmetry data.
For q < 0 the isotropic constraints are much more constraining than the modulation constraints, as seen in the bottom panels of Fig. 48. This is because for large scales the factor k * /k can become large and a negative g * will decrease isotropic power on those scales, which is compensated for by increasing A s and τ. Strongly negative g * values are disallowed by predicting unphysical negative power spectra at low . Note that even the parameter ranges in which the power spectra are reduced to close to zero are likely beyond the perturbative regime for the models in question, and so should be approached with caution. The isotropic constraints still prefer a slightly negative g * , likely due to being able to fit the power deficit at large scales. The joint constraint in this case is then greatly improved by the isotropic data.
Minimum-χ 2 and p values (relative to isotropic simulations) for g * are presented in Table 16. The addition of polarization does not affect the temperature results greatly.
Finally, when allowing the completely general form of quadrupolar modulation, i.e., with no restriction on the g 2m , we present results for the quantity g 2 ≡ m |g 2m | 2 /5 in  47. Histograms of the quantityÔ j0 for the tanh, power-law, and n s gradient modulation models using Planck temperature data combined with 300 statistically isotropic polarization simulations (red outlines) or 300 polarization simulations modulated according to the best-fit parameters from the temperature data (black). The blue lines indicate the values for the actual SMICA polarization data. A large value relative to the isotropic (red) simulations would indicate that the modulation model is preferred over ΛCDM. Marginalized posteriors for quadrupolar modulation parameter g * , using SMICA data for the T T +EE asymmetry constraints (orange curves) and Planck TT,TE,EE+lowE for the isotropic constraints (blue curves), which probe the modification to the power spectrum via Eq. (114). Top: constraints for q = 2 and 1 (left and right, respectively). Bottom: constraints for q = −1 and −2 (left and right, respectively). Strongly negative g * is suppressed for q < 0, due to the unphysical prediction of negative power. Table 16. Minimum-χ 2 g * values for quadrupolar modulation, determined from the SMICA foreground-cleaned maps. Also given are p values, defined as the fraction of isotropic simulations with larger |g * | than the data. The T T results use min = 2 and max = 1200, while EE uses min = 2 and max = 850. These results indicate that the data are consistent with cosmic variance in statistically isotropic skies.

Conclusions
This paper summarizes the status of cosmic inflation in light of the Planck 2018 release. The main improvements are in the Planck polarization likelihoods. The 2018 release now includes a low-HFI polarization likelihood based on the 100-and 143-GHz channels. This likelihood is now the baseline, whereas the Planck 2015 likelihood was based only on the LFI 70-GHz channel data, which also have been updated in this release. Corrections for beam-leakage effects, which had been flagged in the 2015 release as the main limitation of the T E and EE data at that time, have improved the accuracy of the high-polarization likelihoods. Our analyses focus on the results obtained using the Planck baseline likelihoods alone, but results supplemented by the BK14 likelihood (when tensors are included) and a compilation of BAO likelihoods are also given in order to help break cosmological parameter degeneracies. We summarize the main results of this paper in the form of responses to a number of key questions.
1. What is the value of the scalar tilt? Using a characterization of polarization anisotropy better at all multipoles in this release, we find that n s = 0.9649 ± 0.0042 at 68 % CL, including the full information provided by Planck (TT,TE,EE+lowE+lensing). The 2018 uncertainty is approximately 2/3 of that obtained with the Planck 2015 baseline likelihood. Importantly, this determination rules out perfect scale invariance (i.e., n s = 1) at 8.4σ. From an inflationary perspective, this result is consistent with slow-roll inflation evolving towards a natural exit. 2. Does n s depend on the wavelength?
We investigated the possibility of a running spectral index, as well as a running of the running [i.e., the next two (subleading) terms in a power series expansion of ln(P R ) in ln(k)], corresponding to non-negligible third-and fourthorder derivatives of the inflationary potential. Starting with its first 2013 cosmological release, Planck has removed any hint of a running spectral index, which had been suggested by pre-Planck data and would have pointed to inflationary models beyond the slow-roll approximation. Planck 2018 sets dn s /d ln k = −0.005±0.013 as the tightest 95 % CL constraint, when d 2 n s /d ln k 2 = 0. No hints of further extensions, such as running of the running, are found with Planck 2018 data. These results are consistent with the simplest slow-roll dynamics for the inflaton. A detection of running at the level predicted by slow-roll models will require a combination of future ambitious CMB anisotropy experiments and galaxy surveys. 3. Is the Universe spatially flat?
Most simple models of inflation predict a spatially flat universe, although inflationary models with a minimum degree of fine tuning producing a hyperbolic universe have been constructed. Planck has been the first experiment to constrain the spatial curvature at the percent level without any external information, thanks to the CMB lensing likelihood. Although negative values of Ω K ∼ −0.01 provide a non-statistically significant improvement to the fit of Planck temperature and polarization data (compared to the minimal ΛCDM model), Planck 2018 data including lensing constrain Ω K = −0.011 +0.013 −0.012 at 95 % CL. Combining with BAO data further tightens the uncertainty, constraining Ω K to lie within 0.4 % of a flat spatial geometry (at 95 % CL). 4. Are tensor modes required?
Inflationary models predict that tensor modes were also excited during the nearly exponential expansion, with a power spectrum amplitude proportional to the energy scale of inflation. Using the measurement of CMB temperature and E-mode polarization anisotropies from the quadrupole into the acoustic peak region, Planck has reduced the degeneracy between the tensor-to-scalar ratio r and n s , establishing the bound r 0.002 < 0.10 at 95 % CL, assuming n t = −r/8 as predicted by the simplest inflationary models. When the Planck likelihood is combined with the B-mode polariza-59 Planck Collaboration: Constraints on Inflation tion likelihood of the BICEP2-Keck experiment, a tight 95 % CL upper limit of r < 0.064 is obtained, corresponding to a 95 % CL bound on the energy scale of inflation of V 1/4 * < 1.7 × 10 16 GeV. Planck 2018 and BK14 data also set tight bounds on gravitational waves generated in the early Universe when r and n t are varied independently, complementary to the results obtained by the direct detection interferometers LIGO and VIRGO at much higher frequencies. 5. Which inflationary models are best able to account for the data?
Starting with the 2013 release using only a part of the data, Planck has substantially tightened the constraints on slowroll inflationary models, ruling out hybrid models with n s > 1 and power-law inflation (PCI13). In combination with the BK14 data, Planck 2018 now strongly disfavours monomial models with V(φ) ∝ φ p and p ≥ 2 and low-scale SUSY models. Within the representative cases studied in this paper, inflationary models such as R 2 , T and E α-attractor models, D-brane inflation, and those having a potential with exponential tails provide good fits to Planck and BK14 data. We used two methods to reconstruct the inflaton potential beyond the slow-roll approximation: by Taylor expanding the inflaton potential or Hubble parameter in the observable region; and through a free-form reconstruction of the potential with cubic splines. No statistically significant detection beyond the second derivative of the potential was found, suggesting that the slow-roll approximation is adequate for the Planck 2018 likelihood in combination with the BK14 data. 6. What model-independent constraints can be placed on the primordial power spectrum?
We reported on three different methods for the nonparametric reconstruction of the primordial power spectrum (penalized likelihood, a Bayesian spline reconstruction. and a method based on cubic splines). All three methods give broadly consistent results. In no case is any statistically significant evidence for a deviation from a pure power law found. The constraints on the deviations are at the fewpercent level for wavenumbers in the range 0.005 Mpc −1 k 0.2 Mpc −1 probed by the CMB, the precise constraint depending on the level of smoothing allowed. 7. Is there evidence for features in the primordial power spectrum?
We explored several classes of theoretically motivated parametric models with strong departures from a power law for the primordial power spectra and tested their predictions using combinations of Planck temperature and polarization power spectra. We also carried out an analysis using bispectrum data as well. No statistically significant evidence for features was found. 8. Were the primordial cosmological perturbations solely adiabatic?
A key question is whether the primordial cosmological fluctuations consisted exclusively of adiabatic growing-mode perturbations or whether isocurvature perturbations, possibly correlated with the adiabatic mode and with each other, were also excited. The new polarization data has helped to sharpen constraints on the allowed isocurvature fraction compared to the Planck 2015 results. In correlated mixed adiabatic and isocurvature models, the 95 % CL upper bound for the non-adiabatic contribution to the observed CMB temperature variance is |α non-adi | < 1.3 %, 1.7 %, and 1.7 % for CDM, neutrino density, and neutrino velocity isocurvature, respectively. For this release we also report constraints on a scale-invariant compensated baryon-CDM isocurvature mode, which is uncorrelated with the adiabatic mode. This mode would cause an additional lensing-like smoothing at high and modify the lensing potential at 40. By using the temperature, polarization, and lensing data, we obtain the constraint ∆ 2 rms = 0.0037 +0.0016 −0.0021 at 68 % CL for the variance of the baryon isocurvature density perturbation. A detection of isocurvature modes would suggest the need for a theory beyond single-field inflation, which is able to excite only one mode. 9. Were the primordial fluctuations statistically isotropic?
The Planck analysis has confirmed evidence at low statistical significance of anomalies in the CMB temperature anisotropies on large angular scales that are not alleviated in models with nontrivial topology or an anisotropic expansion (Planck Collaboration XVIII 2016). This motivates an exploration of inflation-based models giving such violation of statistical isotropy. We have found no statistically significant evidence in favour of a curvaton model for dipolar asymmetry (compared to the base-ΛCDM model), nor any evidence for a quadrupolar asymmetry in the temperature or polarization anistropies. Theoretical models producing the observed temperature dipolar asymmetry make a prediction for the polarization dipolar asymmetry. We tested whether the fit to the temperature dipolar asymmetry gives a prediction for the polarization asymmetry consistent with the data. We found no statistically significant evidence that the pattern seen in temperature is repeated in polarization. However, the discriminating power of this test is weak, due to the low polarization signal-to-noise ratio on large angular scales.
The Planck 2013, 2015, and 2018 releases have substantially improved the constraints on the space of inflationary models, as described above. Future CMB polarization data will be crucial for further constraining those inflationary models that currently provide an adequate fit to Planck and other data. Forthcoming Emode polarization data will be decisive for determining whether the intriguing features in the temperature power spectrum, such as the deficit at 20-30, the smaller average amplitude at 40, and other anomalies at higher multipoles require new physics or whether these features are simply the result of statistical fluctuations plus instrumental noise. Improved measurements of the B modes promise to constrain inflation even more tightly and it will be interesting to see how the search for B modes evolves. One possibility would be a convincing detection of inflationary gravitational waves, but a tighter upper limit of r 10 −3 is also an achievable outcome. Either case would substantially advance our understanding of inflation and the constraints on the physics of the very early Universe.