Issue 
A&A
Volume 594, October 2016
Planck 2015 results



Article Number  A20  
Number of page(s)  65  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201525898  
Published online  20 September 2016 
Planck 2015 results
XX. Constraints on inflation
^{1} APC, AstroParticule et Cosmologie, Université Paris Diderot, CNRS/IN2P3, CEA/lrfu, Observatoire de Paris, Sorbonne Paris Cité, 10 rue Alice Domon et Léonie Duquet, 75205 Paris Cedex 13, France
^{2} Aalto University Metsähovi Radio Observatory and Dept of Radio Science and Engineering, PO Box 13000, 00076 Aalto, Finland
^{3} African Institute for Mathematical Sciences, 6−8 Melrose Road, Muizenberg, 7945 Cape Town, South Africa
^{4} Agenzia Spaziale Italiana Science Data Center, via del Politecnico snc, 00133 Roma, Italy
^{5} AixMarseille Université, CNRS, LAM (Laboratoire d’Astrophysique de Marseille) UMR 7326, 13388 Marseille, France
^{6} Astrophysics Group, Cavendish Laboratory, University of Cambridge, J J Thomson Avenue, Cambridge CB30 HE, UK
^{7} Astrophysics & Cosmology Research Unit, School of Mathematics, Statistics & Computer Science, University of KwaZuluNatal, Westville Campus, Private Bag X54001, 4000 Durban, South Africa
^{8} Atacama Large Millimeter/submillimeter Array, ALMA Santiago Central Offices, Alonso de Cordova 3107, Vitacura, 763 0355 Casilla, Santiago, Chile
^{9} CGEE, SCS Qd 9, Lote C, Torre C, 4° andar, Ed. Parque Cidade Corporate, CEP 70308200 Brasília, DF, Brazil
^{10} CITA, University of Toronto, 60 St. George St., Toronto, ON M5S 3H8, Canada
^{11} CNRS, IRAP, 9 Av. colonel Roche, BP 44346, 31028 Toulouse Cedex 4, France
^{12} CRANN, Trinity College, Dublin 2, Ireland
^{13} California Institute of Technology, Pasadena, CA 91125, USA
^{14} Centre for Theoretical Cosmology, DAMTP, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK
^{15} Centro de Estudios de Física del Cosmos de Aragón (CEFCA), Plaza San Juan, 1, planta 2, 44001 Teruel, Spain
^{16} Computational Cosmology Center, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA
^{17} Consejo Superior de Investigaciones Científicas (CSIC), 28049 Madrid, Spain
^{18} DSM/Irfu/SPP, CEASaclay, 91191 GifsurYvette Cedex, France
^{19} DTU Space, National Space Institute, Technical University of Denmark, Elektrovej 327, 2800 Kgs. Lyngby, Denmark
^{20} Département de Physique Théorique, Université de Genève, 24 quai E. Ansermet, 1211 Genève 4, Switzerland
^{21} Dark Cosmology Centre, Niels Bohr Institute, University of Copenhagen, Juliane Maries Vej 30, 2100 Copenhagen, Denmark
^{22} Departamento de Astrofísica, Universidad de La Laguna (ULL), 38206 La Laguna, Tenerife, Spain
^{23} Departamento de Física, Universidad de Oviedo, Avda. Calvo Sotelo s/n, 33007 Oviedo, Spain
^{24} Department of Astronomy and Astrophysics, University of Toronto, 50 Saint George Street, Toronto, Ontario, ON M5S 3H41, Canada
^{25} Department of Astrophysics/IMAPP, Radboud University Nijmegen, PO Box 9010, 6500 GL Nijmegen, The Netherlands
^{26} Department of Physics & Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, British Columbia, Canada
^{27} Department of Physics and Astronomy, Dana and David Dornsife College of Letter, Arts and Sciences, University of Southern California, Los Angeles, CA 90089, USA
^{28} Department of Physics and Astronomy, University College London, London WC1E 6BT, UK
^{29} Department of Physics and Astronomy, University of Sussex, Brighton BN1 9QH, UK
^{30} Department of Physics, Florida State University, Keen Physics Building, 77 Chieftan Way, Tallahassee, Florida, USA
^{31} Department of Physics, Gustaf Hällströmin katu 2a, University of Helsinki, 00560 Helsinki, Finland
^{32} Department of Physics, Princeton University, Princeton, NJ 08540, USA
^{33} Department of Physics, University of California, Berkeley, CA 94720, USA
^{34} Department of Physics, University of California, Santa Barbara, California, USA
^{35} Department of Physics, University of Illinois at UrbanaChampaign, 1110 West Green Street, Urbana, Illinois, USA
^{36} Dipartimento di Fisica e Astronomia G. Galilei, Università degli Studi di Padova, via Marzolo 8, 35131 Padova, Italy
^{37} Dipartimento di Fisica e Astronomia, ALMA MATER STUDIORUM, Università degli Studi di Bologna, viale Berti Pichat 6/2, 40127 Bologna, Italy
^{38} Dipartimento di Fisica e Scienze della Terra, Università di Ferrara, via Saragat 1, 44122 Ferrara, Italy
^{39} Dipartimento di Fisica, Università La Sapienza, P.le A. Moro 2, 00133 Roma, Italy
^{40} Dipartimento di Fisica, Università degli Studi di Milano, via Celoria, 16, 00133 Milano, Italy
^{41} Dipartimento di Fisica, Università degli Studi di Trieste, via A. Valerio 2, 00133 Trieste, Italy
^{42} Dipartimento di Matematica, Università di Roma Tor Vergata, via della Ricerca Scientifica, 1, 00133 Roma, Italy
^{43} Discovery Center, Niels Bohr Institute, Blegdamsvej 17, 1165 Copenhagen, Denmark
^{44} Discovery Center, Niels Bohr Institute, Copenhagen University, Blegdamsvej 17, 1165 Copenhagen, Denmark
^{45} European Southern Observatory, ESO Vitacura, Alonso de Cordova 3107, Vitacura, Casilla 19001, Santiago, Chile
^{46} European Space Agency, ESAC, Planck Science Office, Camino bajo del Castillo, s/n, Urbanización Villafranca del Castillo, 28691 Villanueva de la Cañada, Madrid, Spain
^{47} European Space Agency, ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands
^{48} Gran Sasso Science Institute, INFN, viale F. Crispi 7, 67100 L’ Aquila, Italy
^{49} HGSFP and University of Heidelberg, Theoretical Physics Department, Philosophenweg 16, 69120 Heidelberg, Germany
^{50} Helsinki Institute of Physics, Gustaf Hällströmin katu 2, University of Helsinki, 00560 Helsinki, Finland
^{51} INAF–Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, 35131 Padova, Italy
^{52} INAF–Osservatorio Astronomico di Roma, via di Frascati 33, 00040 Monte Porzio Catone, Italy
^{53} INAF–Osservatorio Astronomico di Trieste, via G.B. Tiepolo 11, 34127 Trieste, Italy
^{54} INAF/IASF Bologna, via Gobetti 101, 40127 Bologna, Italy
^{55} INAF/IASF Milano, via E. Bassini 15, 20133 Milano, Italy
^{56} INFN, Sezione di Bologna, via Irnerio 46, 40126 Bologna, Italy
^{57} INFN, Sezione di Roma 1, Università di Roma Sapienza, P.le Aldo Moro 2, 00185 Roma, Italy
^{58} INFN, Sezione di Roma 2, Università di Roma Tor Vergata, via della Ricerca Scientifica, 1, 00185 Roma, Italy
^{59} INFN/National Institute for Nuclear Physics, via Valerio 2, 34127 Trieste, Italy
^{60} IPAG: Institut de Planétologie et d’Astrophysique de Grenoble, Université Grenoble Alpes, IPAG, 38000 Grenoble, France; CNRS, IPAG, 38000 Grenoble, France
^{61} IUCAA, Post Bag 4, Ganeshkhind, Pune University Campus, 411 007 Pune, India
^{62} Imperial College London, Astrophysics group, Blackett Laboratory, Prince Consort Road, London, SW7 2AZ, UK
^{63} Infrared Processing and Analysis Center, California Institute of Technology, Pasadena, CA 91125, USA
^{64} Institut Néel, CNRS, Université Joseph Fourier Grenoble I, 25 rue des Martyrs, 38042 Grenoble, France
^{65} Institut Universitaire de France, 103 bd SaintMichel, 75005 Paris, France
^{66} Institut d’Astrophysique Spatiale, CNRS, Univ. ParisSud, Université ParisSaclay, Bât. 121, 91405 Orsay Cedex, France
^{67} Institut d’Astrophysique de Paris, CNRS (UMR 7095), 98bis boulevard Arago, 75014 Paris, France
^{68} Institut für Theoretische Teilchenphysik und Kosmologie, RWTH Aachen University, 52056 Aachen, Germany
^{69} Institute for Space Sciences, BucharestMagurale, Romania
^{70} Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
^{71} Institute of Theoretical Astrophysics, University of Oslo, Blindern, 0371 Oslo, Norway
^{72} Instituto de Astrofísica de Canarias, C/Vía Láctea s/n, 28205 La Laguna, Tenerife, Spain
^{73} Instituto de Física de Cantabria (CSICUniversidad de Cantabria), Avda. de los Castros s/n, 39005 Santander, Spain
^{74} Istituto Nazionale di Fisica Nucleare, Sezione di Padova, via Marzolo 8, 35131 Padova, Italy
^{75} Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, California, USA
^{76} Jodrell Bank Centre for Astrophysics, Alan Turing Building, School of Physics and Astronomy, The University of Manchester, Oxford Road, Manchester, M13 9PL, UK
^{77} Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA
^{78} Kavli Institute for Cosmology Cambridge, Madingley Road, Cambridge, CB3 0HA, UK
^{79} Kazan Federal University, 18 Kremlyovskaya St., 420008 Kazan, Russia
^{80} LAL, Université ParisSud, CNRS/IN2P3, Orsay, France
^{81} LERMA, CNRS, Observatoire de Paris, 61 avenue de l’Observatoire, 75000 Paris, France
^{82} Laboratoire AIM, IRFU/Service d’Astrophysique – CEA/DSM – CNRS – Université Paris Diderot, Bât. 709, CEASaclay, 91191 GifsurYvette Cedex, France
^{83} Laboratoire Traitement et Communication de l’Information, CNRS (UMR 5141) and Télécom ParisTech, 46 rue Barrault 75634 Paris Cedex 13, France
^{84} Laboratoire de Physique Subatomique et Cosmologie, Université GrenobleAlpes, CNRS/IN2P3, 53 rue des Martyrs, 38026 Grenoble Cedex, France
^{85} Laboratoire de Physique Théorique, Université ParisSud 11 & CNRS, Bâtiment 210, 91405 Orsay, France
^{86} Lawrence Berkeley National Laboratory, Berkeley, California, USA
^{87} Lebedev Physical Institute of the Russian Academy of Sciences, Astro Space Centre, 84/32 Profsoyuznaya st., 117997 Moscow, GSP7, Russia
^{88} Leung Center for Cosmology and Particle Astrophysics, National Taiwan University, 10617 Taipei, Taiwan
^{89} MaxPlanckInstitut für Astrophysik, KarlSchwarzschildStr. 1, 85741 Garching, Germany
^{90} McGill Physics, Ernest Rutherford Physics Building, McGill University, 3600 rue University, Montréal, QC, H3A 2T8, Canada
^{91} National University of Ireland, Department of Experimental Physics, Maynooth, Co. Kildare, Ireland
^{92} Nicolaus Copernicus Astronomical Center, Bartycka 18, 00716 Warsaw, Poland
^{93} Niels Bohr Institute, Blegdamsvej 17, 1165 Copenhagen, Denmark
^{94} Niels Bohr Institute, Copenhagen University, Blegdamsvej 17, 1165 Copenhagen, Denmark
^{95} Nordita (Nordic Institute for Theoretical Physics), Roslagstullsbacken 23, 106 91 Stockholm, Sweden
^{96} Optical Science Laboratory, University College London, Gower Street, London, UK
^{97} SISSA, Astrophysics Sector, via Bonomea 265, 34136 Trieste, Italy
^{98} SMARTEST Research Centre, Università degli Studi eCampus, via Isimbardi 10, 22060 Novedrate (CO), Italy
^{99} School of Physics and Astronomy, Cardiff University, Queens Buildings, The Parade, Cardiff, CF24 3AA, UK
^{100} School of Physics and Astronomy, University of Nottingham, Nottingham NG7 2RD, UK
^{101} Simon Fraser University, Department of Physics, 8888 University Drive, Burnaby BC, Canada
^{102} Sorbonne UniversitéUPMC, UMR 7095, Institut d’Astrophysique de Paris, 98bis boulevard Arago, 75014 Paris, France
^{103} Space Research Institute (IKI), Russian Academy of Sciences, Profsoyuznaya Str, 84/32, 117997 Moscow, Russia
^{104} Space Sciences Laboratory, University of California, Berkeley, CA 94720, USA
^{105} Special Astrophysical Observatory, Russian Academy of Sciences, Nizhnij Arkhyz, Zelenchukskiy region, 369167 KarachaiCherkessian Republic, Russia
^{106} Stanford University, Dept of Physics, Varian Physics Bldg, 382 via Pueblo Mall, Stanford, California, USA
^{107} SubDepartment of Astrophysics, University of Oxford, Keble Road, Oxford OX1 3RH, UK
^{108} Sydney Institute for Astronomy, School of Physics A28, University of Sydney, NSW 2006, Australia
^{109} The Oskar Klein Centre for Cosmoparticle Physics, Department of Physics,Stockholm University, AlbaNova, 106 91 Stockholm, Sweden
^{110} Theory Division, PHTH, CERN, 1211 Geneva 23, Switzerland
^{111} UPMC Univ Paris 06, UMR7095, 98bis boulevard Arago, 75014 Paris, France
^{112} Université de Toulouse, UPSOMP, IRAP, 31028 Toulouse Cedex 4, France
^{113} Universities Space Research Association, Stratospheric Observatory for Infrared Astronomy, MS 23211, Moffett Field, CA 94035, USA
^{114} University of Granada, Departamento de Física Teórica y del Cosmos, Facultad de Ciencias, 18071 Granada, Spain
^{115} University of Granada, Instituto Carlos I de Física Teórica y Computacional, 18071 Granada, Spain
^{116} Warsaw University Observatory, Aleje Ujazdowskie 4, 00478 Warszawa, Poland
^{⋆}
Corresponding author: Martin Bucher, email: bucher@apc.univparis7.fr; Fabio Finelli, email: finelli@iasfbo.inaf.it
Received: 14 February 2015
Accepted: 3 March 2016
We present the implications for cosmic inflation of the Planck measurements of the cosmic microwave background (CMB) anisotropies in both temperature and polarization based on the full Planck survey, which includes more than twice the integration time of the nominal survey used for the 2013 release papers. The Planck full mission temperature data and a first release of polarization data on large angular scales measure the spectral index of curvature perturbations to be n_{s} = 0.968 ± 0.006 and tightly constrain its scale dependence to dn_{s}/ dlnk = −0.003 ± 0.007 when combined with the Planck lensing likelihood. When the Planck highℓ polarization data are included, the results are consistent and uncertainties are further reduced. The upper bound on the tensortoscalar ratio is r_{0.002}< 0.11 (95% CL). This upper limit is consistent with the Bmode polarization constraint r< 0.12 (95% CL) obtained from a joint analysis of the BICEP2/Keck Array and Planck data. These results imply that V(φ) ∝ φ^{2} and natural inflation are now disfavoured compared to models predicting a smaller tensortoscalar ratio, such as R^{2} inflation. We search for several physically motivated deviations from a simple powerlaw spectrum of curvature perturbations, including those motivated by a reconstruction of the inflaton potential not relying on the slowroll approximation. We find that such models are not preferred, either according to a Bayesian model comparison or according to a frequentist simulationbased analysis. Three independent methods reconstructing the primordial power spectrum consistently recover a featureless and smooth over the range of scales 0.008 Mpc^{1} ≲ k ≲ 0.1 Mpc^{1}. At large scales, each method finds deviations from a power law, connected to a deficit at multipoles ℓ ≈ 20−40 in the temperature power spectrum, but at an uncompelling statistical significance owing to the large cosmic variance present at these multipoles. By combining power spectrum and nonGaussianity bounds, we constrain models with generalized Lagrangians, including Galileon models and axion monodromy models. The Planck data are consistent with adiabatic primordial perturbations, and the estimated values for the parameters of the base Λ cold dark matter (ΛCDM) model are not significantly altered when more general initial conditions are admitted. In correlated mixed adiabatic and isocurvature models, the 95% CL upper bound for the nonadiabatic contribution to the observed CMB temperature variance is  α_{non  adi}  < 1.9%, 4.0%, and 2.9% for CDM, neutrino density, and neutrino velocity isocurvature modes, respectively. We have tested inflationary models producing an anisotropic modulation of the primordial curvature power spectrum findingthat the dipolar modulation in the CMB temperature field induced by a CDM isocurvature perturbation is not preferred at a statistically significant level. We also establish tight constraints on a possible quadrupolar modulation of the curvature perturbation. These results are consistent with the Planck 2013 analysis based on the nominal mission data and further constrain slowroll singlefield inflationary models, as expected from the increased precision of Planck data using the full set of observations.
Key words: cosmic background radiation / cosmology: theory / early Universe / inflation
© ESO, 2016
1. Introduction
The precise measurements by Planck^{1} of the cosmic microwave background (CMB) anisotropies covering the entire sky and over a broad range of scales, from the largest visible down to a resolution of approximately 5′, provide a powerful probe of cosmic inflation, as detailed in the Planck 2013 inflation paper (Planck Collaboration XXII 2014, hereafter PCI13). In the 2013 results, the robust detection of the departure of the scalar spectral index from exact scale invariance, i.e., n_{s}< 1, at more than 5σ confidence, as well as the lack of the observation of any statistically significant running of the spectral index, were found to be consistent with simple slowroll models of inflation. Singlefield inflationary models with a standard kinetic term were also found to be compatible with the new tight upper bounds on the primordial nonGaussianity parameters f_{NL} reported in Planck Collaboration XXVI (2014). No evidence of isocurvature perturbations as generated in multifield inflationary models (PCI13) or by cosmic strings or topological defects was found (Planck Collaboration XXV 2014). The Planck 2013 results overall favoured the simplest inflationary models. However, we noted an amplitude deficit for multipoles ℓ ≲ 40 whose statistical significance relative to the sixparameter base Λ cold dark matter (ΛCDM) model is only about 2σ, as well as other anomalies on large angular scales but also without compelling statistical significance (Planck Collaboration XXIII 2014). The constraint on the tensortoscalar ratio, r< 0.12 at 95% CL, inferred from the temperature power spectrum alone, combined with the determination of n_{s}, suggested models with concave potentials.
This paper updates the implications for inflation in the light of the Planck full mission temperature and polarization data. The Planck 2013 cosmology results included only the nominal mission, comprising the first 14 months of the data taken, and used only the temperature data. However, the full mission includes the full 29 months of scientific data taken by the cryogenically cooled high frequency instrument (HFI; which ended when the He/He supply for the final stage of the cooling chain ran out) and the approximately four years of data taken by the low frequency instrument (LFI), which covered a longer period than the HFI because the LFI did not rely on cooling down to 100 mK for its operation. For a detailed discussion of the new likelihood and a comparison with the 2013 likelihood, we refer the reader to Planck Collaboration XI (2016) and Planck Collaboration XIII (2016), but we mention here some highlights of the differences between the 2013 and 2015 data processing and likelihoods: (1) improvements in the data processing such as beam characterization and absolute calibration at each frequency result in a better removal of systematic effects and (2) the 2015 temperature highℓ likelihood uses halfmission crosspower spectra over more of the sky, owing to less aggressive Galactic cuts. The use of polarization information in the 2015 likelihood release contributes to the constraining power of Planck in two principal ways: (1) the measurement of the Emode polarization at large angular scales (presently based on the 70 GHz channel) constrains the reionization optical depth, τ, independently of other estimates using ancillary data; and (2) the measurement of the TE and EE spectra at ℓ ≥ 30 at the same frequencies used for the TT spectra (100, 143, and 217 GHz) helps break parameter degeneracies, particularly for extended cosmological models (beyond the baseline sixparameter model). A full analysis of the Planck lowℓ polarization is still in progress and will be the subject of another forthcoming set of Planck publications.
The Planck 2013 results have sparked a revival of interest in several aspects of inflationary models. We mention here a few examples without the ambition to be exhaustive. A lively debate arose on the conceptual problems of some of the inflationary models favoured by the Planck 2013 data (Ijjas et al. 2013, 2014; Guth et al. 2014; Linde 2014). The interest in the R^{2} inflationary model originally proposed by Starobinsky (1980) increased, since its predictions for cosmological fluctuations (Mukhanov & Chibisov 1981; Starobinsky 1983) are compatible with the Planck 2013 results (PCI13). It has been shown that supergravity motivates a potential similar to the Einstein gravity conformal representation of the R^{2} inflationary model in different contexts (Ellis et al. 2013a,b; Buchmüller et al. 2013; Farakos et al. 2013; Ferrara et al. 2013b). A similar potential can also be generated by spontaneous breaking of conformal symmetry (Kallosh & Linde 2013b).
The constraining power of Planck also motivated a comparison between large numbers of inflationary models (Martin et al. 2014) and stimulated different perspectives on how best to compare theoretical inflationary predictions with observations based on the parameterized dependence of the Hubble parameter on the scale factor during inflation (Mukhanov 2013; Binétruy et al. 2015; GarciaBellido & Roest 2014). The interpretation of the asymmetries on large angular scales (Planck Collaboration XXIII 2014) also prompted a reanalysis of the primordial dipole modulation (Lyth 2013; Liddle & Cortês 2013; Kanno et al. 2013) of curvature perturbations during inflation.
Another recent development has been the renewed interest in possible tensor modes generated during inflation, sparked by the BICEP2 results (BICEP2 Collaboration 2014a,b). The BICEP2 team suggested that the Bmode polarization signal detected at 50 <ℓ< 150 at a single frequency (150 GHz) might be of primordial origin. However, a crucial step in this possible interpretation was excluding an explanation based on polarized thermal dust emission from our Galaxy. The BICEP2 team put forward a number of models to estimate the likely contribution from dust, but at the time relevant observational data were lacking, and this modelling involved a high degree of extrapolation. If dust polarization were negligible in the observed patch of 380 deg^{2}, this interpretation would lead to a tensortoscalar ratio of for a scaleinvariant spectrum. A value of r ≈ 0.2, as suggested by BICEP2 Collaboration (2014b), would have obviously changed the Planck 2013 perspective according to which slowroll inflationary models are favoured, and such a high value of r would also have required a strong running of the scalar spectral index, or some other modification from a simple powerlaw spectrum, to reconcile the contribution of gravitational waves to temperature anisotropies at low multipoles with the observed TT spectrum.
The interpretation of the Bmode signal in terms of gravitational waves alone presented in BICEP2 Collaboration (2014b) was later cast in doubt by Planck measurements of dust polarization at 353 GHz (Planck Collaboration Int. XIX 2015; Planck Collaboration Int. XX 2015; Planck Collaboration Int. XXI 2015; Planck Collaboration Int. XXII 2015). The Planck measurements characterized the frequency dependence of intensity and polarization of the Galactic dust emission, and moreover showed that the polarization fraction is higher than expected in regions of low dust emission. With the help of the Planck measurements of Galactic dust properties (Planck Collaboration Int. XIX 2015), it was shown that the interpretation of the Bmode polarization signal in terms of a primordial tensor signal plus a lensing contribution was not statistically preferred to an explanation based on the expected dust signal at 150 GHz plus a lensing contribution (see also Flauger et al. 2014a; Mortonson & Seljak 2014). Subsequently, Planck Collaboration Int. XXX (2016) extrapolated the PlanckBmode power spectrum of dust polarization at 353 GHz over the multipole range 40 <ℓ< 120 to 150 GHz, showing that the Bmode polarization signal detected by BICEP2 could be entirely due to dust.
More recently, a BICEP2/Keck ArrayPlanck (BKP) joint analysis (BICEP2/Keck Array and Planck Collaborations 2015, herafter BKP) combined the highsensitivity Bmode maps from BICEP2 and Keck Array with the Planck maps at higher frequencies where dust emission dominates. A study of the crosscorrelations of all these maps in the BICEP2 field found the absence of any statistically significant evidence for primordial gravitational waves, setting an upper limit of r< 0.12 at 95% CL (BKP). Although this upper limit is numerically almost identical to the Planck 2013 result obtained combining the nominal mission temperature data with WMAP polarization to remove parameter degeneracies (Planck Collaboration XVI 2014; Planck Collaboration XXII 2014), the BKP upper bound is much more robust against modifications of the inflationary model, since B modes are insensitive to the shape of the predicted scalar anisotropy pattern. In Sect. 13 we explore how the recent BKP analysis constrains inflationary models.
This paper is organized as follows. Section 2 briefly reviews the additional information on the primordial cosmological fluctuations encoded in the polarization angular power spectrum. Section 3 describes the statistical methodology as well as the Planck and other likelihoods used throughout the paper. Sections 4 and 5 discuss the Planck 2015 constraints on scalar and tensor fluctuations, respectively. Section 6 is dedicated to constraints on the slowroll parameters and provides a Bayesian comparison of selected slowroll inflationary models. In Sect. 7 we reconstruct the inflaton potential and the Hubble parameter as a Taylor expansion of the inflaton in the observable range without relying on the slowroll approximation. The reconstruction of the curvature perturbation power spectrum is presented in Sect. 8. The search for parameterized features is presented in Sect. 9, and combined constraints from the Planck 2015 power spectrum and primordial nonGaussianity derived in Planck Collaboration XVII (2016) are presented in Sect. 10. The analysis of isocurvature perturbations combined and correlated with curvature perturbations is presented in Sect. 11. In Sect. 12 we study the implications of relaxing the assumption of statistical isotropy of the primordial fluctuations. We discuss two examples of anisotropic inflation in light of the tests of isotropy performed in Planck Collaboration XVI (2016). Section 14 presents some concluding remarks.
2. What new information does polarization provide?
This section provides a short theoretical overview of the extra information provided by polarization data over that of temperature alone. (More details can be found in White et al. 1994; Ma & Bertschinger 1995; Bucher 2015, and references therein.) In Sect. 2 of the Planck 2013 inflation paper (PCI13), we gave an overview of the relation between the inflationary potential and the threedimensional primordial scalar and tensor power spectra, denoted as and respectively. (The scalar variable ℛ is defined precisely in Sect. 3.) We shall not repeat the discussion there, instead referring the reader to PCI13 and references therein.
Under the assumption of statistical isotropy, which is predicted in all simple models of inflation, the twopoint correlations of the CMB anisotropies are described by the angular power spectra and where ℓ is the multipole number. (See Kamionkowski et al. 1997; Zaldarriaga & Seljak 1997; Seljak & Zaldarriaga 1997; Hu & White 1997; Hu et al. 1998 and references therein for early discussions elucidating the role of polarization.) In principle, one could also envisage measuring and , but in theories where parity symmetry is not explicitly or spontaneously broken, the expectation values for these cross spectra (i.e., the theoretical cross spectra) vanish, although the observed realizations of the cross spectra are not exactly zero because of cosmic variance.
The CMB angular power spectra are related to the threedimensional scalar and tensor power spectra via the transfer functions and so that the contributions from scalar and tensor perturbations are (1)and (2)respectively, where The scalar and tensor primordial perturbations are uncorrelated in the simplest models, so the scalar and tensor power spectra add in quadrature, meaning that (3)Roughly speaking, the form of the linear transformations encapsulated in the transfer functions and probe the late time physics, whereas the primordial power spectra and are solely determined by the primordial Universe, perhaps not so far below the Planck scale if largefield inflation turns out to be correct.
Fig. 1 Comparison of transfer functions for the scalar and tensor modes. The CMB transfer functions and , where , define the linear transformations mapping the primordial scalar and tensor cosmological perturbations to the CMB anisotropies as seen by us on the sky today. These functions are plotted for two representative values of the multipole number: ℓ = 2 (in black) and ℓ = 65 (in red). 
To better understand this connection, it is useful to plot and compare the shapes of the transfer functions for representative values of ℓ and characterize their qualitative behavior. Referring to Fig. 1, we emphasize the following qualitative features:
 1.
For the scalar mode transfer functions, of which only and are nonvanishing (because to linear order, a threedimensional scalar mode cannot contribute to the B mode of the polarization), both transfer functions start to rise at more or less the same small values of k (due to the centrifugal barrier in the Bessel differential equation), but falls off much faster at large k and thus smooths sharp features in to a lesser extent than This means that polarization is more powerful than temperature for reconstructing possible sharp features in the scalar primordial power spectrum provided that the required signaltonoise is available.
 2.
For the tensor modes, starts rising at about the same small k as and but falls off faster with increasing k than On the other hand, the polarization components, and have a shape completely different from any of the other transfer functions. The shape of and is much wider in ln(k) than the scalar polarization transfer function, with a variance ranging from 0.5 to 1.0 decades. These functions exhibit several oscillations with a period smaller than that for scalar transfer functions, due to the difference between the sound velocity for scalar fluctuations and the light velocity for gravitational waves (Polarski & Starobinsky 1996; Lesgourgues et al. 2000).
The inability of scalar modes to generate Bmode polarization (apart from the effects of lensing) has an important consequence. For the primordial tensor modes, polarization information, especially information concerning the Bmode polarization, offers powerful potential for discovery or for establishing upper bounds. Planck 2013 and WMAP established upper bounds on a possible tensor mode contribution using alone, but these bounds crucially relied on assuming a simple form for the scalar primordial power spectrum. For example, as reported in PCI13, when a simple power law was generalized to allow for running, the bound on the tensor contribution degraded by approximately a factor of two. The new joint BICEP2/Keck ArrayPlanck upper bound (see Sect. 13), however, is much more robust and cannot be avoided by postulating baroque models that alter the scale dependence of the scalar power spectrum.
3. Methodology
This section describes updates to the formalism used to describe cosmological models and the likelihoods used with respect to the Planck 2013 inflation paper (PCI13).
3.1. Cosmological model
The cosmological models that predict observables such as the CMB anisotropies rely on inputs specifying the conditions and physics at play during different epochs of the history of the Universe. The primordial inputs describe the power spectrum of the cosmological perturbations at a time when all the observable modes were situated outside the Hubble radius. The inputs from this epoch consist of the primordial power spectra, which may include scalar curvature perturbations, tensor perturbations, and possibly also isocurvature modes and their correlations. The late time (i.e., z ≲ 10^{4}) cosmological inputs include parameters such as ω_{b},ω_{c},Ω_{Λ}, and τ, which determine the conditions when the primordial perturbations become imprinted on the CMB and also the evolution of the Universe between last scattering and today, affecting primarily the angular diameter distance. Finally, there is a socalled “nuisance” component, consisting of parameters that determine how the measured CMB spectra are contaminated by unsubtracted Galactic and extragalactic foreground contamination. The focus of this paper is on the primordial inputs and how they are constrained by the observed CMB anisotropy, but we cannot completely ignore the other nonprimordial parameters because their presence and uncertainties must be dealt with in order to correctly extract the primordial information of interest here.
As in PCI13, we adopt the minimal sixparameter spatially flat base ΛCDM cosmological model as our baseline for the late time cosmology, mainly altering the primordial inputs, i.e., the simple powerlaw spectrum parameterized by the scalar amplitude and spectral index for the adiabatic growing mode, which in this minimal model is the only late time mode excited. This model has four free nonprimordial cosmological parameters (ω_{b},ω_{c},θ_{MC},τ; for a more detailed account of this model, we refer the reader to Planck Collaboration XIII 2016). On occasion, this assumption will be relaxed in order to consider the impact of more complex alternative late time cosmologies on our conclusions about inflation. Some of the commonly used cosmological parameters are defined in Table 1.
Primordial, baseline, and optional latetime cosmological parameters.
3.2. Primordial spectra of cosmological fluctuations
In inflationary models, comoving curvature (ℛ) and tensor (h) fluctuations are amplified by the nearly exponential expansion from quantum vacuum fluctuations to become highly squeezed states resembling classical states. Formally, this quantum mechanical phenomenon is most simply described by the evolution in conformal time, η, of the mode functions for the gaugeinvariant inflaton fluctuation, δφ, and for the tensor fluctuation, h: (4)with for scalars and (x,y) = (a,h) for tensors. Here a is the scale factor, primes indicate derivatives with respect to η, and and H = ȧ/a are the proper time derivative of the inflaton and the Hubble parameter, respectively. The curvature fluctuation, ℛ, and the inflaton fluctuation, δφ, are related via . Analytic and numerical calculations of the predictions for the primordial spectra of cosmological fluctuations generated during inflation have reached high standards of precision, which are more than adequate for our purposes, and the largest uncertainty in testing specific inflationary models arises from our lack of knowledge of the history of the Universe between the end of inflation and the present time, during the socalled “epoch of entropy generation”.
This paper uses three different methods to compare inflationary predictions with Planck data. The first method consists of a phenomenological parameterization of the primordial spectra of scalar and tensor perturbations according to: where A_{s} (A_{t}) is the scalar (tensor) amplitude and n_{s} (n_{t}), dn_{s}/ dlnk (dn_{t}/ dlnk), and d^{2}n_{s}/ dlnk^{2} are the scalar (tensor) spectral index, the running of the scalar (tensor) spectral index, and the running of the running of the scalar spectral index, respectively. h^{+ , ×} denotes the amplitude of the two polarization states (+ , ×) of gravitational waves and k_{∗} is the pivot scale. Unless otherwise stated, the tensortoscalar ratio, (7)is fixed to −8n_{t}, which is the relation that holds when inflation is driven by a single slowrolling scalar field with a standard kinetic term^{2}. We will use a parameterization analogous to Eq. (5) with no running for the power spectra of isocurvature modes and their correlations in Sect. 11.
Conventions and definitions for inflation physics.
The second method exploits the analytic dependence of the slowroll power spectra of primordial perturbations in Eqs. (5) and (6) on the values of the Hubble parameter and the hierarchy of its time derivatives, known as the Hubble flow functions (HFF): ϵ_{1} = −Ḣ/H^{2}, , with i ≥ 1. We will use the analytic power spectra calculated up to second order using the Green’s function method (Gong & Stewart 2001; Leach et al. 2002; see Habib et al. 2002; Martin & Schwarz 2003; and Casadio et al. 2006 for alternative derivations). The spectral indices and the relative scale dependence in Eqs. (5) and (6) are given in terms of the HFFs by: where C ≡ ln2 + γ_{E}−2 ≈ −0.7296 (γ_{E} is the EulerMascheroni constant). See the Appendix of PCI13 for more details. Primordial spectra as functions of the ϵ_{i} will be employed in Sect. 6, and the expressions generalizing Eqs. (8) to (11) for a general Lagrangian p(φ,X), where X ≡ −g^{μν}∂_{μ}φ∂_{ν}φ/ 2, will be used in Sect. 10. The good agreement between the first and second method as well as with alternative approximations of slowroll spectra is illustrated in the Appendix of PCI13.
The third method is fully numerical, suitable for models where the slowroll conditions are not well satisfied and analytical approximations for the primordial fluctuations are not available. Two different numerical codes, the inflation module of Lesgourgues & Valkenburg (2007) as implemented in CLASS (Lesgourgues 2011; Blas et al. 2011) and ModeCode (Adams et al. 2001; Peiris et al. 2003; Mortonson et al. 2009; Easther & Peiris 2012), are used in Sects. 7 and 10, respectively.^{3}
Conventions for the functions and symbols used to describe inflationary physics are defined in Table 2.
3.3. Planck data
The Planck data processing proceeding from timeordered data to maps has been improved for this 2015 release in various aspects (Planck Collaboration II 2016; Planck Collaboration VII 2016). We refer the interested reader to Planck Collaboration II (2016) and Planck Collaboration VII (2016) for details, and we describe here two of these improvements. The absolute calibration has been improved using the orbital dipole and more accurate characterization of the Planck beams. The calibration discrepancy between Planck and WMAP described in Planck Collaboration XXXI (2014) for the 2013 release has now been greatly reduced. At the time of that release, a blind analysis for primordial power spectrum reconstruction described a broad feature at ℓ ≈ 1800 in the temperature power spectrum, which was most prominent in the 217 × 217 GHz autospectra (PCI13). In work done after the Planck 2013 data release, this feature was shown to be associated with imperfectly subtracted systematic effects associated with the 4 K cooler lines, which were particularly strong in the first survey. This systematic effect was shown to potentially lead to 0.5σ shifts in the cosmological parameters, slightly increasing n_{s} and H_{0}, similarly to the case in which the 217 × 217 channel was excised from the likelihood (Planck Collaboration XV 2014; Planck Collaboration XVI 2014). The Planck likelihood (Planck Collaboration XI 2016) is based on the full mission data and comprises temperature and polarization data (see Fig. 2).
Fig. 2 PlanckTT (top), highℓTE (centre), and highℓEE (bottom) angular power spectra. Here . 
Planck lowℓ likelihood
The Planck lowℓ temperaturepolarization likelihood uses foregroundcleaned LFI 70 GHz polarization maps together with the temperature map obtained from the Planck 30 to 353 GHz channels by the Commander component separation algorithm over 94% of the sky (see Planck Collaboration IX 2016 for further details). The Planck polarization map uses the LFI 70 GHz (excluding Surveys 2 and 4) lowresolution maps of Q and U polarization from which polarized synchrotron and thermal dust emission components have been removed using the LFI 30 GHz and HFI 353 GHz maps as templates, respectively. (See Planck Collaboration XI 2016 for more details.) The polarization map covers the 46% of the sky outside the lowP polarization mask.
The lowℓ likelihood is pixelbased and treats the temperature and polarization at the same resolution of , or HEALpix (Górski et al. 2005) N_{side} = 16. Its multipole range extends from ℓ = 2 to ℓ = 29 in TT, TE, EE, and BB. In the 2015 Planck papers the polarization part of this likelihood is denoted as “lowP”.^{4} This Planck lowℓ likelihood replaces the Planck temperature lowℓ Gibbs module combined with the WMAP 9yr lowℓ polarization module used in the Planck 2013 cosmology papers (denoted by WP), which used lower resolution polarization maps at N_{side} = 8 (about ). With this Planckonly lowℓ likelihood module, the basic Planck results presented in this release are completely independent of external information.
The Planck lowmultipole likelihood alone implies τ = 0.067 ± 0.022 (Planck Collaboration XI 2016), a value smaller than the value inferred using the WP polarization likelihood, τ = 0.089 ± 0.013, used in the Planck 2013 papers (Planck Collaboration XV 2014). See Planck Collaboration XIII (2016) for the important implications of this decrease in τ for reionization. However, the LFI 70 GHz and WMAP polarization maps are in very good agreement when both are foregroundcleaned using the HFI 353 GHz map as a polarized dust template (see Planck Collaboration XI 2016 for further details). Therefore, it is useful to construct a noiseweighted combination to obtain a joint Planck/WMAP low resolution polarization data set, also described in Planck Collaboration XI (2016), using as a polarization mask the union of the WMAP P06 and Planck lowP polarization masks and keeping 74% of the sky. The polarization part of the combined low multipole likelihood is called lowP+WP. This combined low multipole likelihood gives (Planck Collaboration XI 2016).
Planck highℓ likelihood
Following Planck Collaboration XV (2014), and Planck Collaboration XI 2016 for polarization, we use a Gaussian approximation for the highℓ part of the likelihood (30 <ℓ< 2500), so that (12)where a constant offset has been discarded. Here Ĉ is the data vector, C(θ) is the model prediction for the parameter value vector θ, and ℳ is the covariance matrix. For the data vector, we use 100 GHz, 143 GHz, and 217 GHz halfmission crosspower spectra, avoiding the Galactic plane as well as the brightest point sources and the regions where the CO emission is the strongest. We retain 66% of the sky for 100 GHz, 57% for 143 GHz, and 47% for 217 GHz for the T masks, and respectively 70%, 50%, and 41% for the Q, U masks. Following Planck Collaboration XXX (2014), we do not mask for any other Galactic polarized emission. All the spectra are corrected for the beam and pixel window functions using the same beam for temperature and polarization. (For details see Planck Collaboration XI 2016.)
The model for the crossspectra can be written as (13)where C^{cmb}(θ) is the CMB power spectrum, which is independent of the frequency, is the foreground model contribution for the crossfrequency spectrum μ × ν, and c_{μ} is the calibration factor for the μ × μ spectrum. The model for the foreground residuals includes the following components: Galactic dust, clustered cosmic infrared background (CIB), thermal and kinetic SunyaevZeldovich (tSZ and kSZ) effect, tSZ correlations with CIB, and point sources, for the TT foreground modeling; and for polarization, only dust is included. All the components are modelled by smooth C_{ℓ} templates with free amplitudes, which are determined along with the cosmological parameters as the likelihood is explored. The tSZ and kSZ models are the same as in 2013 (see Planck Collaboration XV 2014), although with different priors (Planck Collaboration XI 2016; Planck Collaboration XIII 2016), while the CIB and tSZCIB correlation models use the updated CIB models described in Planck Collaboration XXX (2014). The point source contamination is modelled as Poisson noise with an independent amplitude for each frequency pair. Finally, the dust contribution uses an effective smooth model measured from high frequency maps. Details of our dust and noise modelling can be found in Planck Collaboration XI (2016). The dust is the dominant foreground component for TT at ℓ< 500, while the point source component, and for 217 × 217 also the CIB component, dominate at high ℓ. The other foreground components are poorly determined by Planck. Finally, our treatment of the calibration factors and beam uncertainties and mismatch are described in Planck Collaboration XI (2016).
The covariance matrix accounts for the correlation due to the mask and is computed following the equations in Planck Collaboration XV (2014), extended to polarization in Planck Collaboration XI (2016) and references therein. The fiducial model used to compute the covariance is based on a joint fit of base ΛCDM and nuisance parameters obtained with a previous version of the matrix. We iterate the process until the parameters stop changing. For more details, see Planck Collaboration XI (2016).
The joint unbinned covariance matrix is approximately of size 23 000×23 000. The memory and speed requirements for dealing with such a huge matrix are significant, so to reduce its size, we bin the data and the covariance matrix to compress the data vector size by a factor of 10. The binning uses varying bin width with Δℓ = 5 for 29 <ℓ< 100, Δℓ = 9 for 99 <ℓ< 1504, Δℓ = 17 for 1503 <ℓ< 2014, and Δℓ = 33 for 2013 <ℓ< 2509, and a weighting in ℓ(ℓ + 1) to flatten the spectrum. Where a higher resolution is desirable, we also use a more finely binned version (“bin3”, unbinned up to ℓ = 80 and Δℓ = 3 beyond that) as well as a completely unbinned version (“bin1”). We use odd bin sizes, since for an azimuthally symmetric mask, the correlation between a multipole and its neighbours is symmetric, oscillating between positive and negative values. Using the base ΛCDM model and singleparameter classical extensions, we confirmed that the cosmological and nuisance parameter fits with or without binning are indistinguishable.
As discussed in Planck Collaboration XI (2016) and Planck Collaboration XIII (2016), the TE and EE highℓ data are not free of small systematic effects, such as leakage from temperature to polarization. Although the propagated effects of these residual systematics on cosmological parameters are small and do not alter the conclusions of this paper, we mainly refer to Planck TT+lowP in combination with the Planck lensing or additional data sets as the most reliable results for this release.
Planck CMB bispectrum
We use measurements of the nonGaussianity amplitude f_{NL} from the CMB bispectrum presented in Planck Collaboration XVII (2016). NonGaussianity constraints have been obtained using three optimal bispectrum estimators: separable template fitting (also known as “KSW”), binned, and modal. The maps analysed are the Planck 2015 full mission sky maps, both in temperature and in E polarization, as cleaned with the four component separation methods SMICA, SEVEM, NILC, and Commander. The map is masked to remove the brightest parts of the Galaxy as well as the brightest point sources and covers approximately 70% of the sky. In this paper we mainly exploit the joint constraints on equilateral and orthogonal nonGaussianity (after removing the integrated SachsWolfe effectlensing bias), , from T only, and , from T and E (68% CL). For reference, the constraints on local nonGaussianity are from T only, and from T and E (68% CL). Starting from a Gaussian f_{NL}likelihood, which is an accurate assumption in the regime of small primordial nonGaussianity, we use these constraints to derive limits on the sound speed of the inflaton fluctuations (or other microscopic parameters of inflationary models; Planck Collaboration XXIV 2014). The bounds on the sound speed for various models are then used in combination with Planck power spectrum data.
Planck CMB lensing data
Some of our analysis includes the Planck 2015 lensing likelihood, presented in Planck Collaboration XV (2016), which utilizes the nonGaussian trispectrum induced by lensing to estimate the power spectrum of the lensing potential, . This signal is extracted using a full set of temperature and polarizationbased quadratic lensing estimators (Okamoto & Hu 2003) applied to the SMICA CMB map over approximately 70% of the sky, as described in Planck Collaboration IX (2016). We have used the conservative bandpower likelihood, covering multipoles 40 ≤ ℓ ≤ 400. This provides a measurement of the lensing potential power at the 40σ level, giving a 2.5%accurate constraint on the overall lensing power in this multipole range. The measurement of the lensing power spectrum used here is approximately twice as powerful as the measurement used in our previous 2013 analysis (Planck Collaboration XXII 2014; Planck Collaboration XVII 2014), which used temperatureonly data from the Planck nominal mission data set.
3.4. NonPlanck data
BAO data
Baryon acoustic oscillations (BAO) are the counterpart in the late time matter power spectrum of the acoustic oscillations seen in the CMB multipole spectrum (Eisenstein et al. 2005). Both originate from coherent oscillations of the photonbaryon plasma before these two components become decoupled at recombination. Measuring the position of these oscillations in the matter power spectra at different redshifts constrains the expansion history of the universe after decoupling, thus removing degeneracies in the interpretation of the CMB anisotropies.
In this paper, we combine constraints on (the ratio between the sphericallyaveraged distance scale D_{V} to the effective survey redshift, , and the sound horizon, r_{s}) inferred from 6dFGRS data (Beutler et al. 2011) at , the SDSSMGS data (Ross et al. 2015) at , and the SDSSDR11 CMASS and LOWZ data (Anderson et al. 2014) at redshifts and 0.32. For details see Planck Collaboration XIII (2016).
Joint BICEP2/Keck Array and Planck constraint on r
Since the Planck temperature constraints on the tensortoscalar ratio are close to the cosmic variance limit, the inclusion of data sets sensitive to the expected Bmode signal of primordial gravitational waves is particularly useful. In this paper, we provide results including the joint analysis crosscorrelating BICEP2/Keck Array observations and Planck (BKP). Combining the more sensitive BICEP2/Keck Array Bmode polarization maps in the approximately 400 deg^{2} BICEP2 field with the Planck maps at higher frequencies where dust dominates allows a statistical analysis taking into account foreground contamination. Using BB auto and crossfrequency spectra between BICEP2/Keck Array (150 GHz) and Planck (217 and 353 GHz), BKP find a 95% upper limit of r_{0.05}< 0.12.
3.5. Parameter estimation and model comparison
Much of this paper uses a Bayesian approach to parameter estimation, and unless otherwise specified, we assign broad tophat prior probability distributions to the cosmological parameters listed in Table 1. We generate posterior probability distributions for the parameters using either the MetropolisHastings algorithm implemented in CosmoMC (Lewis & Bridle 2002) or MontePython (Audren et al. 2013), the nested sampling algorithm MultiNest (Feroz & Hobson 2008; Feroz et al. 2009, 2013), or PolyChord, which combines nested sampling with slice sampling (Handley et al. 2015). The latter two also compute the Bayesian evidence needed for model comparison. Nevertheless, χ^{2} values are often provided as well (using CosmoMC’s implementation of the BOBYQA algorithm (Powell 2009) for maximizing the likelihood), and other parts of the paper employ frequentist methods when appropriate.
Fig. 3 Comparison of the marginalized joint 68% and 95% CL constraints on (n_{s},τ) (left panel), (n_{s},Ω_{b}h^{2}) (middle panel), and (n_{s},θ_{MC}) (right panel), for Planck 2013 (grey contours), Planck TT+lowP (red contours), Planck TT+lowP+lensing (green contours), and Planck TT, TE, EE+lowP (blue contours). 
Confidence limits on the parameters of the base ΛCDM model, for various combinations of Planck 2015 data, at the 68% confidence level.
4. Constraints on the primordial spectrum of curvature perturbations
One of the most important results of the Planck nominal mission was the determination of the departure from scale invariance for the spectrum of scalar perturbations at high statistical significance (Planck Collaboration XVI 2014; Planck Collaboration XXII 2014). We now update these measurements with the Planck full mission data in temperature and polarization.
4.1. Tilt of the curvature power spectrum
For the base ΛCDM model with a powerlaw power spectrum of curvature perturbations, the constraint on the scalar spectral index, n_{s}, with the Planck full mission temperature data is (14)This result is compatible with the Planck 2013 constraint, n_{s} = 0.9603 ± 0.0073 (Planck Collaboration XV 2014; Planck Collaboration XVI 2014). See Fig. 3 for the accompanying changes in τ, Ω_{b}h^{2}, and θ_{MC}. The shift towards higher values for n_{s} with respect to the nominal mission results is due to several improvements in the data processing and likelihood which are discussed in Sect. 3, including the removal of the 4 K cooler systematics. For the values of other cosmological parameters in the base ΛCDM model, see Table 3. We also provide the results for the base ΛCDM model and extended models online.^{5}
When the Planck highℓ polarization is combined with temperature, we obtain (15)together with τ = 0.079 ± 0.017 (68% CL), which is consistent with the TT+lowP results. The Planck highℓ polarization pulls τ up to a slightly higher value. When the Planck lensing measurement is added to the temperature data, we obtain (16)with τ = 0.066 ± 0.016 (68% CL). The shift towards slightly smaller values of the optical depth is driven by a marginal preference for a smaller primordial amplitude, A_{s}, in the Planck lensing data (Planck Collaboration XV 2016). Given that the temperature data provide a sharp constraint on the combination e^{− 2τ}A_{s}, a slightly lower A_{s} requires a smaller optical depth to reionization.
4.2. Viability of the HarrisonZeldovich spectrum
Even though the estimated scalar spectral index has risen slightly with respect to the Planck 2013 release, the assumption of a HarrisonZeldovich (HZ) scaleinvariant spectrum (Harrison 1970; Peebles & Yu 1970; Zeldovich 1972) continues to be disfavoured (with a modest increase in significance, from 5.1σ in 2013 to 5.6σ today), because the error bar on n_{s} has decreased. The value of n_{s} inferred from the Planck 2015 temperature plus largescale polarization data lies 5.6 standard deviations away from unity (with a corresponding Δχ^{2} = 29.9), if one assumes the base ΛCDM latetime cosmological model. If we consider more general reionization models, parameterized by a principal component analysis (Mortonson & Hu 2008) instead of τ (where reionization is assumed to have occurred instantaneously), we find Δχ^{2} = 14.9 for n_{s} = 1. Previously, simple oneparameter extensions of the base model, such as ΛCDM+N_{eff} (where N_{eff} is the effective number of neutrino flavours) or ΛCDM+Y_{P} (where Y_{P} is the primordial value of the helium mass fraction), could nearly reconcile the Planck temperature data with n_{s} = 1. They now lead to Δχ^{2} = 7.6 and 9.3, respectively. For any of the cosmological models that we have considered, the Δχ^{2} by which the HZ model is penalized with respect to the tilted model has increased since the 2013 analysis (PCI13) thanks to the constraining power of the full mission temperature data. Adding Planck highℓ polarization data further disfavours the HZ model: in ΛCDM, the χ^{2} increases by 57.8, for general reionization we obtain Δχ^{2} = 41.3, and for ΛCDM+N_{eff} and ΛCDM+Y_{P} we find Δχ^{2} = 22.5 and 24.0, respectively.
4.3. Running of the spectral index
The running of the scalar spectral index is constrained by the Planck 2015 full mission temperature data to (17)The combined constraint including highℓ polarization is (18)Adding the Planck CMB lensing data to the temperature data further reduces the central value for the running, i.e., dn_{s}/ dlnk = −0.0033 ± 0.0074 (68% CL, Planck TT+lowP+lensing).
The central value for the running has decreased in magnitude with respect to the Planck 2013 nominal mission (Planck Collaboration XVI 2014 found dn_{s}/ dlnk = −0.013 ± 0.009; see Fig. 4), and the improvement of the maximum likelihood with respect to a powerlaw spectrum is smaller, Δχ^{2} ≈ −0.8. Among the different effects contributing to the decrease in the central value of the running with respect to the Planck 2013 result, we mention a change in HFI beams at ℓ ≲ 200 (Planck Collaboration XIII 2016). Nevertheless, the deficit of power at low multipoles in the Planck 2015 temperature power spectrum contributes to a preference for slightly negative values of the running, but with low statistical significance.
Fig. 4 Marginalized joint 68% and 95% CL for (n_{s},dn_{s}/ dlnk) using Planck TT+lowP and Planck TT, TE, EE+lowP. Constraints from the Planck 2013 data release are also shown for comparison. For comparison, the thin black stripe shows the prediction for singlefield monomial chaotic inflationary models with 50 <N_{∗}< 60. 
Constraints on the primordial perturbation parameters for ΛCDM+r and ΛCDM+r+dn_{s}/ dlnk models from Planck.
The Planck constraints on n_{s} and dn_{s}/ dlnk are remarkably stable against the addition of the BAO likelihood. The combination with BAO shifts n_{s} to slighly higher values and shrinks its uncertainty by about 30% when only highℓ temperature is considered, and by only about 15% when highℓ temperature and polarization are combined. In slowroll inflation, the running of the scalar spectral index is connected to the third derivative of the potential (Kosowsky & Turner 1995). As was the case for the nominal mission results, values of the running compatible with the Planck 2015 constraints can be obtained in viable inflationary models (Kobayashi & Takahashi 2011).
When the running of the running is allowed to float, the Planck TT+lowP (Planck TT, TE, EE+lowP) data give: (19)at the pivot scale k_{∗} = 0.05 Mpc^{1}. Allowing for running of the running provides a better fit to the temperature spectrum at low multipoles, such that Δχ^{2} ≈ −4.8 (−4.9) for TT+lowP (TT, TE, EE+lowP), but is not statistically preferred over the simplest ΛCDM model.
Note that the inclusion of smallscale data such as Lyα might further constrain the running of the spectral index and its derivative. The recent analysis of the BOSS onedimensional Lyα flux power spectrum presented in PalanqueDelabrouille et al. (2015) and Rossi et al. (2015) was optimized for measuring the neutrino mass. It does not include constraints on the spectral index running, which would require new dedicated Nbody simulations. Hence we do not include Lyα constraints here.
In Sect. 7 on inflaton potential reconstruction we will show that the data cannot accomodate a significant running but are compatible with a larger running of the running.
4.4. Suppression of power on the largest scales
Although not statistically significant, the trend for a negative running or positive running of the running observed in the last subsection was driven by the lack of power in the Planck temperature power spectrum at low multipoles, already mentioned in the Planck 2013 release. This deficit could potentially be explained by a primordial spectrum featuring a depletion of power only at large wavelengths. Here we investigate two examples of such models.
We first update the analysis (already presented in PCI13) of a powerlaw spectrum multiplied by an exponential cutoff: (20)This simple parameterization is motivated by models with a short inflationary stage in which the onset of the slowroll phase coincides with the time when the largest observable scales exited the Hubble radius during inflation. The curvature spectrum is then strongly suppressed on those scales. We apply tophat priors on the parameter λ_{c}, controlling the steepness of the cutoff, and on the logarithm of the cutoff scale, k_{c}. We choose prior ranges λ_{c} ∈ [ 0,10 ] and ln(k_{c}/ Mpc^{1}) ∈ [ −12,−3 ]. For Planck TT+lowP (Planck TT, TE, EE+lowP), the bestfit model has λ_{c} = 0.50 (0.53), ln(k_{c}/ Mpc^{1}) = −7.98 (−7.98), n_{s} = 0.9647 (0.9649), and improves the effective χ^{2} by a modest amount, Δχ^{2} ≈ −3.4 (−3.4).
As a second model, we consider a broken powerlaw spectrum for curvature perturbations: (21)with A_{low} = A_{s}(k_{b}/k_{∗})^{− δ} to ensure continuity at k = k_{b}. Hence this model, like the previous one, has two parameters, and also suppresses power at large wavelengths when δ> 0. We assume tophat priors δ ∈ [ 0,2 ] and ln(k_{b}/ Mpc^{1}) ∈ [ −12,−3 ], and standard uniform priors for ln(10^{10}A_{s}) and n_{s}. The best fit to Planck TT+lowP (Planck TT, TE, EE+lowP) is found for n_{s} = 0.9658 (0.9647), δ = 1.14 (1.14), and ln(k_{b}/ Mpc^{1}) = −7.55 (−7.57), with a very small χ^{2} improvement of Δχ^{2} ≈ −1.9 (−1.6).
We conclude that neither of these two models with two extra parameters is preferred over the base ΛCDM model. (See also the discussion of a step inflationary potential in Sect. 9.1.1.)
5. Constraints on tensor modes
In this section, we focus on the Planck 2015 constraints on tensor perturbations. Unless otherwise stated, we consider that the tensor spectral index satisfies the standard inflationary consistency condition to lowest order in slow roll, n_{t} = −r/ 8. We recall that r is defined at the pivot scale k_{∗} = 0.05 Mpc^{1}. However, for comparison with other studies, we also report our bounds in terms of the tensortoscalar ratio r_{0.002} at k_{∗} = 0.002 Mpc^{1}.
5.1. Planck 2015 upper bound on r
The constraints on the tensortoscalar ratio inferred from the Planck full mission data for the ΛCDM+r model are: Table 4 also shows the bounds on n_{s} in each of these cases.
These results slightly improve over the constraint r_{0.002}< 0.12 (95% CL) derived from the Planck 2013 temperature data in combination with WMAP largescale polarization data (Planck Collaboration XVI 2014; Planck Collaboration XXII 2014). The constraint obtained by Planck temperature and polarization on large scales is tighter than the PlanckBmode 95% CL upper limit from the 100 and 143 GHz HFI channels, r< 0.27 (Planck Collaboration XI 2016). The constraints on r reported in Table 4 can be translated into upper bounds on the energy scale of inflation at the time when the pivot scale exits the Hubble radius using (26)This gives an upper bound on the Hubble parameter during inflation of H_{∗}/M_{pl}< 3.6 × 10^{5} (95% CL) for Planck TT+lowP.
These bounds are relaxed when allowing for a scale dependence of the scalar and tensor spectral indices. In that case, we assume that the tensor spectral index and its running are fixed by the standard inflationary consistency condition at second order in slow roll. We obtain with n_{s} = 0.9667 ± 0.0066 (68% CL). At the standard pivot scale, k_{∗} = 0.05 Mpc^{1}, the bound is stronger (r< 0.17 at 95% CL), because k_{∗} is closer to the scale at which n_{s} and r decorrelate. The constraint on r_{0.002} in Eq. (27) is 21% tighter than the corresponding Planck 2013 constraint. The mean value of the running in Eq. (28) is higher (lower in absolute value) than with Planck 2013 by 45%. Figures 5 and 6 clearly illustrate this significant improvement with respect to the previous Planck data release. Table 4 shows how bounds on (r,n_{s}, dn_{s}/ dlnk) are affected by the lensing reconstruction, BAO, or highℓ polarization data. The tightest bounds are obtained in combination with polarization: with n_{s} = 0.9644 ± 0.0049 (68% CL).
Fig. 5 Marginalized joint confidence contours for (n_{s},dn_{s}/ dlnk), at the 68% and 95% CL, in the presence of a nonzero tensor contribution, and using Planck TT+lowP or Planck TT, TE, EE+lowP. Constraints from the Planck 2013 data release are also shown for comparison. The thin black stripe shows the prediction of singlefield monomial inflation models with 50 <N_{∗}< 60. 
Fig. 6 Marginalized joint confidence contours for (n_{s},r), at the 68% and 95% CL, in the presence of running of the spectral indices, and for the same data combinations as in the previous figure. 
Neither the Planck full mission constraints in Eqs. (22)−(25) nor those including a running in Eqs. (27) and (29) are compatible with the interpretation of the BICEP2 Bmode polarization data in terms of primordial gravitational waves (BICEP2 Collaboration 2014b). Instead they are in excellent agreement with the results of the BICEP2/Keck ArrayPlanck crosscorrelation analysis, as discussed in Sect. 13.
5.2. Dependence of the r constraints on the lowℓ likelihood
The constraints on r discussed above are further tightened by adding WMAP polarization information on large angular scales. The Planck measurement of CMB polarization on large angular scales at 70 GHz is consistent with the WMAP 9year one, based on the K, Q, and Vbands (at 30, 40, and 60 GHz, respectively), once the Planck 353 GHz channel is used to remove the dust contamination, instead of the theoretical dust model used by the WMAP team (Page et al. 2007). (For a detailed discussion, see Planck Collaboration XI 2016.) By combining Planck TT data with LFI 70 GHz and WMAP polarization data on large angular scales, we obtain a 35% reduction of uncertainty, giving τ = 0.074 ± 0.012 (68% CL) and n_{s} = 0.9660 ± 0.060 (68% CL) for the base ΛCDM model. When tensors are added, the bounds become When tensors and running are both varied, we obtain r_{0.002}< 0.14 (95% CL) and dn_{s}/ dlnk = −0.010 ± 0.008 (68% CL) for Planck TT+lowP+WP. These constraints are all tighter than those based on Planck TT+lowP only.
5.3. The tensortoscalar ratio and the lowℓ deficit in temperature
As noted previously (Planck Collaboration XV 2014; Planck Collaboration XVI 2014; Planck Collaboration XXII 2014), the lowℓ temperature data display a slight lack of power compared to the expectation of the bestfit tensorfree base ΛCDM model. Since tensor fluctuations add power on small scales, the effect will be exacerbated in models allowing r> 0.
In order to quantify this tension, we compare the observed constraint on r to that inferred from simulated Planck data. In the simulations, we assume the underlying fiducial model to be tensorfree, with parameters close to the base ΛCDM bestfit values. We limit the simulations to mock temperature power spectra only and fit these spectra with an exact lowℓ likelihood for 2 ≤ ℓ ≤ 29 (see Perotto et al. 2006), and a highℓ Gaussian likelihood for 30 ≤ ℓ ≤ 2508 based on the frequencycombined, foregroundmarginalized, unbinned Planck temperature power spectrum covariance matrix. Additionally, we impose a Gaussian prior of τ = 0.07 ± 0.02.
Based on 100 simulated data sets, we find a 95% CL upper limit on the tensortoscalar ratio of . The corresponding constraint from real data (using lowℓCommander temperature data, the frequencycombined, foregroundmarginalized, unbinned Planck highℓ TT power spectrum, and the same prior on τ as above) reads r< 0.123, confirming that the actual constraint is tighter than what one would have expected. However, the actual constraint is not excessively unusual: out of the 100 simulations, 4 lead to an even tighter bound, corresponding to a significance of about 2σ. Thus, under the hypothesis of the base ΛCDM cosmology, the upper limit on r that we get from the data is not implausible as a chance fluctuation of the low multipole power.
To illustrate the contribution of the lowℓ temperature power deficit to the estimates of cosmological parameters, we show as an example in Fig. 7 how n_{s} shifts towards lower values when the ℓ< 30 temperature information is discarded (we will refer to this case as “Planck TT−lowT”). The shift in n_{s} is approximately −0.005 (or −0.003 when the lowP likelihood is replaced by a Gaussian prior τ = 0.07 ± 0.02). These shifts exceed those found in Sect. 4.4, where a primordial power spectrum suppressed on large scales was fitted to the data.
Fig. 7 Onedimensional posterior probabilities for n_{s} for the base ΛCDM model obtained by excluding temperature multipoles for ℓ< 30 (“TT−lowT”), while either keeping lowℓ polarization data, or in addition replacing them with a Gaussian prior on τ. 
Figure 8 displays the posterior probability for r for various combinations of data sets, some of which exclude the ℓ< 30TT data. This leads to the very conservative bounds r ≲ 0.24 and r ≲ 0.23 at 95% CL when combined with the lowP likelihood or with the Gaussian prior τ = 0.07 ± 0.020, respectively.
Fig. 8 Onedimensional posterior probabilities for r for various data combinations, either including or not including temperature multipoles for ℓ< 30, and compared with the baseline choice (Planck TT+lowP, black curve). 
5.4. Relaxing assumptions on the latetime cosmological evolution
Fig. 9 Marginalized joint 68% and 95% CL for (n_{s},r_{0.002}) using Planck TT+lowP+BAO (upper panel) and Planck TT, TE, EE+lowP (lower panel). 
Constraints on extensions of the ΛCDM+r cosmological model for Planck TT+lowP+lensing, Planck TT+lowP+BAO, and Planck TT, TE, EE+lowP.
As in the Planck 2013 release (PCI13), we now ask how robust the Planck results on the tensortoscalar ratio are against assumptions on the latetime cosmological evolution. The results are summarized in Table 5, and some particular cases are illustrated in Fig. 9. Constraints on r turn out to be remarkably stable for oneparameter extensions of the ΛCDM+r model, with the only exception the ΛCDM+r+Ω_{K} case in the absence of the late time information from Planck lensing or BAO data. The weak trend towards Ω_{K}< 0, i.e., towards a positively curved (closed) universe from the temperature and polarization data alone, and the wellknown degeneracy between Ω_{K} and H_{0}/Ω_{m} lead to a slight suppression of the SachsWolfe plateau in the scalar temperature spectrum. This leaves more room for a tensor component.
This further degeneracy when r is added builds on the negative values for the curvature allowed by Planck TT+lowP, at 95% CL (Planck Collaboration XIII 2016). The exploitation of the information contained in the Planck lensing likelihood leads to a tighter constraint, at 95% CL, which improves on the Planck 2013 results ( at 95% CL). However, due to the remaining degeneracies left by the uncertainties in polarization on large angular scales, a full appreciation of the improvement due to the full mission temperature and lensing data can be obtained by using lowP+WP, which leads to at 95% CL. Note that the negative values allowed for the curvature are decreased in magnitude when the running is allowed, suggesting that the lowℓ temperature deficit is contributing to the estimate of the spatial curvature.
The trend found for ΛCDM+r+Ω_{K} is even clearer when spatial curvature and the running of the spectral index are varied at the same time. In this case, the Planck temperature plus polarization data are compatible with r values as large as 0.19 (95% CL), at the cost of an almost 4σ deviation from spatial flatness (which, however, disappears as soon as lensing or BAO data are considered).
6. Implications for singlefield slowroll inflation
In this section we study the implications of Planck 2015 constraints on standard slowroll singlefield inflationary models.
6.1. Constraints on slowroll parameters
We first present the Planck 2015 constraints on slowroll parameters obtained through the analytic perturbative expansion in terms of the HFFs ϵ_{i} for the primordial spectra of cosmological fluctuations during slowroll inflation (Stewart & Lyth 1993; Gong & Stewart 2001; Leach et al. 2002). When restricting to first order in ϵ_{i}, we obtain When highℓ polarization is included we obtain ϵ_{1}< 0.0066 at 95% CL and at 68% CL. When secondorder contributions in the HFFs are included, we obtain When highℓ polarization is included we obtain ϵ_{1}< 0.011 at 95% CL, at 68% CL, and −0.32 <ϵ_{3}< 0.89 at 95% CL.
The potential slowroll parameters are obtained as derived parameters by using their exact expressions as function of ϵ_{i} (Leach et al. 2002; Finelli et al. 2010): where V(φ) is the inflaton potential, the subscript φ denotes the derivative with respect to φ, and M_{pl} = (8πG)^{− 1 / 2} is the reduced Planck mass (see also Table 2).
By using Eqs. (39) and (40) with ϵ_{3} = 0 and the primordial power spectra to lowest order in the HFFs, the derived constraints for the first two slowroll potential parameters are: When highℓ polarization is included we obtain ϵ_{V}< 0.0067 at 95% CL and at 68% CL. By using Eqs. (39)−(41) with ϵ_{4} = 0 and the primordial power spectra to second order in the HFFs, the derived constraints for the slowroll potential parameters are: When highℓ polarization is included we obtain ϵ_{V}< 0.011 at 95% CL, and and , both at 68% CL.
Fig. 10 Marginalized joint 68% and 95% CL regions for (ϵ_{1},ϵ_{2}) (top panel) and (ϵ_{V},η_{V}) (bottom panel) for Planck TT+lowP (red contours), Planck TT, TE, EE+lowP (blue contours), and compared with the Planck 2013 results (grey contours). 
In Figs. 10 and 11 we show the 68% CL and 95% CL of the HFFs and the derived potential slowroll parameters with and without the highℓ polarization and compare these values with the Planck 2013 results.
6.2. Implications for selected inflationary models
The predictions to lowest order in the slowroll approximation for (n_{s},r) at k = 0.002 Mpc^{1} of a few inflationary models with a representative uncertainty for the entropy generation stage (50 <N_{∗}< 60) are shown in Fig. 12. Figure 12 updates Fig. 1 of PCI13 with the same notation.
In the following we discuss the implications of Planck TT+lowP+BAO data for selected slowroll inflationary models by taking into account the uncertainties in the entropy generation stage. We model these uncertainties by two parameters, as in PCI13: the energy scale ρ_{th} by which the Universe has thermalized, and the parameter w_{int} which characterizes the effective equation of state between the end of inflation and the energy scale specified by ρ_{th}. We use the primordial power spectra of cosmological fluctuations generated during slowroll inflation parameterized by the HFFs, ϵ_{i}, to second order, which can be expressed in terms of the number of efolds to the end of inflation, N_{∗}, and the parameters of the considered inflationary model, using modified routines of the public code ASPIC^{6} (Martin et al. 2014). For the number of efolds to the end of inflation (Liddle & Leach 2003; Martin & Ringeval 2010) we use the expression (PCI13) (47)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 consider the pivot scale k_{∗} = 0.002 Mpc^{1}, g_{th} = 10^{3}, and ϵ_{end} = 1. We consider the uniform priors for the cosmological parameters listed in Table 6. We also consider a logarithmic prior on 10^{10}A_{s} (over the interval [ (e^{2.5},e^{3.7} ]) and ρ_{th} (over the interval [ (1 TeV)^{4},ρ_{end} ]). We consider both the case in which w_{int} is kept fixed at zero and the case in which it is allowed to vary with a uniform prior in the range −1 / 3 <w_{int}< 1 / 3.
Priors for cosmological parameters used in the Bayesian comparison of inflationary models.
We have validated the slowroll approach by crosschecking the Bayes factor computations against the fully numerical inflationary mode equation solver ModeCode coupled to the PolyChord sampler. For each inflationary model we provide in Table 7 and in the main text the Δχ^{2} value with respect to the base ΛCDM model, computed with the CosmoMC implementation of the BOBYQA algorithm for maximizing the likelihood, and the Bayesian evidence with respect to the R^{2} inflationary model (Starobinsky 1980), computed by CosmoMC connected to CAMB, using MultiNest as the sampler.
Fig. 11 Marginalized joint 68% and 95% CL regions for (ϵ_{1},ϵ_{2},ϵ_{3}) (top panels) and (bottom panels) for Planck TT+lowP (red contours), Planck TT, TE, EE+lowP (blue contours), and compared with the Planck 2013 results (grey contours). 
Power law potentials
We first investigate the class of inflationary models with a single monomial potential (Linde 1983): (48)in which inflation occurs for large values of the inflaton, φ>M_{pl}. The predictions for the scalar spectral index and the tensortoscalar ratio at first order in the slowroll approximation are n_{s}−1 ≈ −2(n + 2) / (4N_{∗} + n) and r ≈ 16n/ (4N_{∗} + n), respectively. By assuming a dust equation of state (i.e., w_{int} = 0) prior to thermalization, the cubic and quartic potentials are strongly disfavoured by lnB = −11.6 and lnB = −23.3, respectively. The quadratic potential is moderately disfavoured by lnB = −4.7. Other values, such as n = 4 / 3, 1, and 2 / 3, motivated by axion monodromy (Silverstein & Westphal 2008; McAllister et al. 2010), are compatible with Planck data with w_{int} = 0.
Results of the inflationary model comparison.
Small modifications occur when considering the effective equation of state parameter, w_{int} = (n−2) / (n + 2), defined by averaging over the coherent oscillation regime which follows inflation (Turner 1983). The Bayes factors are slightly modified when w_{int} is allowed to float, as shown in Table 7.
Hilltop models
In hilltop models (Boubekeur & Lyth 2005), with potential (49)the inflaton rolls away from an unstable equilibrium. The predictions to first order in the slowroll approximation are r ≈ 8p^{2}(M_{pl}/μ)^{2}x^{2p−2}/ (1−x^{p})^{2} and n_{s}−1 ≈ −2p(p−1)(M_{pl}/μ)^{2}x^{p−2}/ (1−x^{p})−3r/ 8, where x = φ_{∗}/μ. As in PCI13, the ellipsis in Eq. (49) and in what follows indicates higherorder terms that are negligible during inflation but ensure positiveness of the potential.
By sampling log _{10}(μ/M_{pl}) within the prior [ 0.30,4.85 ] for p = 2, we obtain log _{10}(μ/M_{pl}) > 1.02 (1.05) at 95% CL and lnB = −2.6 (−2.4) for w_{int} = 0 (allowing w_{int} to float).
An exact potential which could also apply after inflation, instead of the approximated one in Eq. (49), might be needed for a better comparison among different models. For μ/M_{pl} ≫ 1, hilltop models as defined in Eq. (49) by neglecting the additional terms denoted by the ellipsis lead to n_{s}−1 ≈ −3r/ 8, the same prediction as for the previously discussed linear potential, V(φ) ∝ φ. By considering a double well potential, V(φ) = Λ^{4} [ 1−φ^{2}/ (2μ^{2}) ] ^{2}, instead, we obtain a slightly worse Bayes factor than the hilltop p = 2 model, lnB = −3.1 (−2.3) for w_{int} = 0 (w_{int} allowed to vary). This different result can be easily understood. Although the double well potential is equal to the hilltop model for φ ≪ μ, it approximates V(φ) ∝ φ^{2} for μ/M_{pl} ≫ 1. Since a linear potential is a better fit to Planck than φ^{2}, the fit of the double well potential is therefore worse than the hilltop p = 2 case for μ/M_{pl} ≫ 1, and this partially explains the slightly different Bayes factors obtained.
In the p = 4 case, we obtain log _{10}(μ/M_{pl}) > 1.05 (1.02) at 95% CL and lnB = −2.8 (−2.6) for w_{int} = 0 (allowing w_{int} to float), assuming a prior range [−2,2 ] for log _{10}(μ/M_{pl}).
Natural inflation
In natural inflation (Freese et al. 1990; Adams et al. 1993) a nonperturbative shift symmetry is invoked to suppress radiative corrections leading to the periodic potential (50)where f is the scale which determines the curvature of the potential. We sample log _{10}(f/M_{pl}) within the prior [ 0.3,2.5 ] as in PCI13. We obtain log _{10}(f/M_{pl}) > 0.84 (> 0.83) at 95% CL and lnB = −2.4 (−2.3) for w_{int} = 0 (allowing w_{int} to vary).
Note that the superPlanckian value for f required by observations is not necessarily a problem for this class of models. When several fields φ_{i} with a cosine potential as in Eq. (50) and scales f_{i} appear in the Lagrangian, an effective singlefield inflationary trajectory can be found for a suitable choice of parameters (Kim et al. 2005). In such a setting, the superPlanckian value of the effective scale f required by observations can be obtained even if the original scales satisfy f_{i} ≪ M_{pl} (Kim et al. 2005).
Dbrane inflation
Inflation can arise from physics involving extra dimensions. If the standard model of particle physics is confined to our 3dimensional brane, the distance between our brane and antibrane can drive inflation. We consider the following parameterization for the effective potential driving inflation: (51)Sampling log _{10}(μ/M_{pl}) using a uniform prior over [−6,0.3 ] , we consider p = 4 (Kachru et al. 2003; Dvali et al. 2001) and p = 2 (GarciaBellido et al. 2002). The predictions for r and n_{s} can be obtained from the hilltop case with the substitution p → −p. These models agree with the Planck data with a Bayes factor of lnB = −0.4 (−0.6) and lnB = −0.7 (−0.9) for p = 4 and p = 2, respectively, for w_{int} = 0 (allowing w_{int} to vary).
Potentials with exponential tails
Exponential potentials are ubiquitous in inflationary models motivated by supergravity and string theory (Goncharov & Linde 1984; Stewart 1995; Dvali & Tye 1999; Burgess et al. 2002; Cicoli et al. 2009). We restrict ourselves to analysing the following class of potentials: (52)As for the hilltop models described earlier, the ellipsis indicates possible higherorder terms that are negligible during inflation but ensure positiveness of the potential. These models predict r ≈ 8q^{2}e^{− 2qφ/Mpl}/ (1−e^{− qφ/Mpl})^{2} and n_{s}−1 ≈ −q^{2}e^{− qφ/Mpl}(2 + e^{− qφ/Mpl}) / (1−e^{− qφ/Mpl})^{2} with a slowroll trajectory characterized by N ≈ f(φ/M_{pl})−f(φ_{end}/M_{pl}), with f(x) = (e^{qx}−qx) /q^{2}. By sampling log _{10}(q/M_{pl}) with a uniform prior over [−3,3 ], we obtain a Bayes factor of −0.6 for w_{int} = 0 (−0.9 when w_{int} is allowed to vary).
Fig. 12 Marginalized joint 68% and 95% CL regions for n_{s} and r at k = 0.002 Mpc^{1} from Planck compared to the theoretical predictions of selected inflationary models. Note that the marginalized joint 68% and 95% CL regions have been obtained by assuming dn_{s}/ dlnk = 0. 
Spontaneously broken SUSY
Hybrid models (Copeland et al. 1994; Linde 1994) predicting n_{s}> 1 are strongly disfavoured by the Planck data, as for the first cosmological release (PCI13). An example of a hybrid model predicting n_{s}< 1 is the case in which slowroll inflation is driven by loop corrections in spontaneously broken supersymmetric (SUSY) grand unified theories (Dvali et al. 1994) described by the potential (53)where α_{h}> 0 is a dimensionless parameter. Note that for α_{h} ≪ 1, this model leads to the same predictions as the powerlaw potential for p ≪ 1 to lowest order in the slowroll approximation. By sampling log _{10}(α_{h}) on a flat prior [−2.5,1 ], we obtain a Bayes factor of −2.2 for w_{int} = 0 (−1.7 when w_{int} is allowed to vary).
R^{2} inflation
The first inflationary model proposed (Starobinsky 1980), with action (54)still lies within the Planck 68% CL constraints, as for the Planck 2013 release (PCI13). This model corresponds to the potential (55)in the Einstein frame, which leads to the slowroll predictions n_{s}−1 ≈ −2 /N and r ≈ 12 /N^{2} (Mukhanov & Chibisov 1981; Starobinsky 1983).
After the Planck 2013 release, several theoretical developments supported the model in Eq. (54) beyond the original motivation of including quantum effects at oneloop (Starobinsky 1980). Noscale supergravity (Ellis et al. 2013a), the largefield regime of superconformal Dterm inflation (Buchmüller et al. 2013), or recent developments in minimal supergravity (Farakos et al. 2013; Ferrara et al. 2013b) can lead to a generalization of the potential in Eq. (55) (see Ketov & Starobinsky 2011 for a previous embedding of R^{2} inflation in F(ℛ) supergravity). The potential in Eq. (55) can also be generated by spontaneous breaking of conformal symmetry (Kallosh & Linde 2013b). This inflationary model has Δχ^{2} ≈ 0.8 (0.3) larger than the base ΛCDM model with no tensors for w_{int} = 0 (for w_{int} allowed to vary). We obtain 54 <N_{∗}< 62 (53 <N_{∗}< 64) at 95% CL for w_{int} = 0 (for w_{int} allowed to vary), compatible with the theoretical prediction, N_{∗} = 54 (Starobinsky 1980; Vilenkin 1985; Gorbunov & Panin 2011).
α attractors
We now study two classes of inflationary models motivated by recent developments in conformal symmetry and supergravity (Kallosh et al. 2013). The first class has been motivated by considering a vector rather than a chiral multiplet for the inflaton in supergravity (Ferrara et al. 2013a) and corresponds to the potential (Kallosh et al. 2013) (56)To lowest order in the slowroll approximation, these models predict and based on an inflationary trajectory characterized by N ≈ g(φ/M_{pl})−g(φ_{end}/M_{pl}) with . The relation between N and φ can be inverted through the use of the Lambert functions, as carried out for other potentials (Martin et al. 2014). By sampling log _{10}(α^{2}) with a flat prior over [ 0,4 ], we obtain log _{10}(α^{2}) < 1.7 (2.0) at 95% CL and a Bayes factor of −1.8 (−2) for w_{int} = 0 (for w_{int} allowed to vary).
The second class of models has been called superconformal α attractors (Kallosh et al. 2013) and can be understood as originating from a different generating function with respect to the first class. This second class is described by the following potential (Kallosh et al. 2013): (57)This is the simplest class of models with spontaneous breaking of conformal symmetry, and for α = m = 1 reduces to the original model introduced by Kallosh & Linde (2013b). The potential in Eq. (57) leads to the following slowroll predictions (Kallosh et al. 2013): where . The predictions of this second class of models interpolate between those of a largefield chaotic model, V(φ) ∝ φ^{2m}, for α ≫ 1 and the R^{2} model for α ≪ 1.
For α we adopt the same priors as for the previous class in Eq. (56). By fixing m = 1, we obtain log _{10}(α^{2}) < 2.3 (2.5) at 95% CL and a Bayes factor of −2.3 (−2.2) for w_{int} = 0 (when w_{int} is allowed to vary). When m is allowed to vary as well with a flat prior in the range [ 0,2 ], we obtain 0.02 <m< 1 (m< 1) at 95% CL for w_{int} = 0 (when w_{int} is allowed to vary).
Nonminimally coupled inflaton
Inflationary predictions are quite sensitive to a nonminimal coupling ξRφ^{2} of the inflaton to the Ricci scalar. One of the most interesting effects due to ξ ≠ 0 is to reconcile the quartic potential V(φ) = λφ^{4}/ 4 with Planck observations, even for ξ ≪ 1. Nonminimal coupling leads as well to attractor behaviour towards predictions similar to those in R^{2} inflation (Kaiser & Sfakianakis 2014; Kallosh & Linde 2013a).
The Higgs inflation model (Bezrukov & Shaposhnikov 2008), in which inflation occurs with and ξ ≫ 1 for φ ≫ φ_{0}, leads to the same predictions as the R^{2} model to lowest order in the slowroll approximation at tree level (see Barvinsky et al. 2008; and Bezrukov & Shaposhnikov 2009 for the inclusion of loop corrections). It is therefore in agreement with the Planck constraints, as for the first cosmological data release (PCI13).
We summarize below our findings for Planck lowP+BAO.

Monomial potentials with integral n> 2 are strongly disfavoured withrespect to R^{2}.

The Bayes factor prefers R^{2} over chaotic inflation with monomial quadratic potential by odds of 110:1 under the assumption of a dust equation of state during the entropy generation stage.

R^{2}inflation has the strongest evidence (i.e., the greatest Bayes factor) among the models considered here. However, care must be taken not to overinterpret small differences in likelihood lacking statistical significance.

The models closest to R^{2} in terms of evidence are brane inflation and exponential inflation, which have one more parameter than R^{2}. Both brane inflation considered in Eq. (51) and exponential inflation in Eq. (52) approximate the linear potential for a large portion of parameter space (for μ/M_{pl} ≫ 1 and q ≫ 1, respectively). For this reason these models have a higher Bayesian evidence (although not at a statistically significant level) than those that approximate a quadratic potential, as do αattractors, for instance.

In the models considered here, the Δχ^{2} obtained by allowing w to vary is modest (i.e., less than approximately 1.6 with respect to w_{int} = 0). The gain in the logarithm of the Bayesian evidence is even smaller, since an extra parameter is added.
7. Reconstruction of the potential and analysis beyond slowroll approximation
7.1. Introdution
In the previous section, we derived constraints on several types of inflationary potentials assumed to account for the inflaton dynamics between the time at which the largest observable scales crossed the Hubble radius during inflation and the end of inflation. The full shape of the potential was used in order to identify when inflation ends, and thus the field value φ_{∗} when the pivot scale crosses the Hubble radius.
In Sect. 6 of PCI13, we explored another approach, consisting of reconstructing the inflationary potential within its observable range without making any assumptions concerning the inflationary dynamics outside that range. Indeed, given that the number of efolds between the observable range and the end of inflation can always be adjusted to take a realistic value, any potential shape giving a primordial spectrum of scalar and tensor perturbations in agreement with observations is a valid candidate. Inflation can end abruptly by a phase transition, or can last a long time if the potential becomes very flat after the observable region has been crossed. Moreover, there could be a short inflationary stage responsible for the origin of observable cosmological perturbations, and another inflationary stage later on (but before nucleosynthesis), thus contributing to the total N_{∗}.
In Sect. 6 of PCI13, we performed this analysis with a full integration of the inflaton and metric perturbation modes, so that no slowroll approximation was made. The only assumption was that primordial scalar perturbations are generated by the fluctuations of a single inflaton field with a canonical kinetic term. Since in this approach one is only interested in the potential over a narrow range of observable scales (centred around the field value φ_{∗} when the pivot scale crosses the Hubble radius), it is reasonable to test relatively simple potential shapes described by a small number of free parameters.
This approach gave very similar results to calculations based on the standard slowroll analysis. This agreement can be explained by the fact that the Planck 2013 data already preferred a primordial spectrum very close to a power law, at least over most of the observable range. Hence the 2013 data excluded strong deviations from slowroll inflation, which would either produce a large running of the spectral index or imprint more complicated features on the primordial spectrum. However, this conclusion did not apply to the largest scales observable by Planck, for which cosmic variance and slightly anomalous data points remained compatible with significant deviations from a simple power law spectrum. The most striking result in Sect. 6 of PCI13 was that a less restricted functional form for the inflaton potential gave results compatible with a rather steep potential at the beginning of the observable window, leading to a “notsoslow” roll stage during the first few observable efolds. This explains the shape of the potential in Fig. 14 of PCI13 for a Taylor expansion at order n = 4 and in the region where φ−φ_{∗} ≤ −0.2. However, such features were only partially explored because the method used for potential reconstruction did not allow for an arbitrary value of the inflation velocity at the beginning of the observable window. Instead, our code imposed that the inflaton already tracked the inflationary attractor solution when the largest observable modes crossed the Hubble scale.
Given that the Planck 2015 data establish even stronger constraints on the primordial power spectrum than the 2013 results, it is of interest to revisit the reconstruction of the potential V(φ). Section 7.2 presents some new results following the same approach as in PCI13 (explained previously in Lesgourgues & Valkenburg 2007; and Mortonson et al. 2011). But in the present work, we also present some more general results, independent of any assumption concerning the initial field velocity when the inflaton enters the observable window. Following previous studies (Kinney 2002; Kinney et al. 2006; Peiris & Easther 2006a,b, 2008; Easther & Peiris 2006; Lesgourgues et al. 2008; Powell & Kinney 2007; Hamann et al. 2008; Norena et al. 2012), we reconstruct the Hubble function H(φ), which determines both the potential V(φ) through (60)and the solution φ(t) through (61)with H′(φ) = ∂H/∂φ. Note that these two relations are exact. In Sect. 7.3, we fit H(φ) directly to the data, implicitly including all canonical singlefield models in which the inflaton is rolling not very slowly (ϵ not much smaller than unity) just before entering the observable window, and the issue of having to start sufficiently early in order to allow the initial transient to decay is avoided. The only drawback in reconstructing H(φ) is that one cannot systematically test the simplest analytic forms for V(φ) in the observable range (for instance, polynomials of order n = 1,3,5,... in (φ−φ_{∗})). But our goal in this section is to explore how much one can deviate from slowroll inflation in general, independently of the shape of the underlying inflaton potential.
Numerical reconstruction of the potential slowroll parameters beyond any slowroll approximation when the potential is Taylor expanded to nth order, using Planck TT+lowP+BAO.
7.2. Reconstruction of a smooth inflaton potential
Following the approach of PCI13, we Taylor expand the inflaton potential around φ = φ_{∗} to order n = 2,3,4. To obtain fasterconverging Markov chains, instead of imposing flat priors on the Taylor coefficients { V,V_{φ},...,V_{φφφφ} }, we sample the potential slowroll (PSR) parameters related to the former as indicated in Table 2. We stress that this is just a choice of prior and does not imply any kind of slowroll approximation in the calculation of the primordial spectra.
The results are given in Table 8 (for Planck TT+lowP+BAO) and Fig. 13 (for the same data set and also for Planck TT, TE, EE+lowP). The second part of Table 8 shows the corresponding values of the spectral parameters n_{s}, dn_{s}/ dlnk, and r_{0.002} as measured for each numerical primordial spectrum (at the pivot scales k = 0.05 Mpc^{1} for the scalar and 0.002 Mpc^{1} for the tensor spectra), as well as the reionization optical depth. We also show in Fig. 14 the derived distribution of each coefficient V_{i} (with a nonflat prior) and in Fig. 15 the reconstructed shape of the bestfit inflation potentials in the observable window. Finally, the posterior distribution of the derived parameters r_{0.002} and dn_{s}/ dlnk is displayed in Fig. 16.
Fig. 13 Posterior distributions for the first four potential slowroll parameters when the potential is Taylor expanded to nth order using Planck TT+lowP+BAO (filled contours) or TT, TE, EE+lowP (dashed contours). The primordial spectra are computed beyond the slowroll approximation. 
Fig. 14 Posterior distributions for the coefficients of the inflation potential Taylor expanded to nth order (in natural units where ) reconstructed beyond the slowroll approximation using Planck TT+lowP+BAO (filled contours) or TT, TE, EE+lowP (dashed contours). The plot shows only half of the results; the other half is symmetric, with opposite signs for V_{φ} and V_{φφφ}. Note that, unlike Fig. 13, the parameters shown here do not have flat priors, since they are mapped from the slowroll parameters. 
Fig. 15 Observable range of the bestfit inflaton potentials, when V(φ) is Taylor expanded to the nth order around the pivot value φ_{∗} in natural units (where assuming a flat prior on ϵ_{V}, η_{V}, , and and using Planck TT+lowP+BAO. Potentials obtained under the transformation (φ−φ_{∗}) → (φ_{∗}−φ) leave the same observable signature and are also allowed. The sparsity of potentials with a small V_{0} = V(φ_{∗}) is due to our use of a flat prior on ϵ_{V} rather than on ln(V_{0}); in fact, V_{0} is unbounded from below in the n = 2 and 3 results. The axis ranges are identical to those in Fig. 20, to make the comparison easier. 
Fig. 16 Posterior distribution for the tensortoscalar ratio (at k = 0.002 Mpc^{1}) and for the running parameter dn_{s}/ dlnk (at k = 0.05 Mpc^{1}), for the potential reconstructions in Sects. 7.2 and 7.3. The V(φ) reconstruction gives the solid curves for Planck TT+lowP+BAO, or dashed for TT, TE, EE+lowP. The H(φ) reconstruction gives the dotted curves for Planck TT+lowP+BAO, or dasheddotted for TT, TE, EE+lowP. The tensortoscalar ratio appears as a derived parameter, but by taking a flat prior on either ϵ_{V} or ϵ_{H}, we implicitly also take a nearly flat prior on r. The same applies to dn_{s}/ dlnk. 
Figure 13 shows that bounds are very similar when temperature data are combined with either highℓ polarization data or BAO data. This gives a hint of the robustness of these results. For both data sets, the error bars on the PSR parameters are typically smaller by a factor of 1.5 than in PCI13.
Since potentials with n = 2 cannot generate a significant running, the bounds on the scalar spectral index and the tensortoscalar ratio and the bestfit models are very similar to those obtained with the ΛCDM+r model in Sect. 5 and Table 4. On the other hand, in the n = 3 model, results follow the trend of the previous ΛCDM+r+dn_{s}/ dlnk analysis. The data prefer potentials with V_{φ} and V_{φφφ} of the same sign, generating a significant negative running (as can be seen in Fig. 16). This trend for V_{φφφ} occurs because a scalar spectrum with negative running reduces the power on large scales, and provides a better fit to lowℓ temperature multipoles. However, such a running also suppresses power on small scales, so cannot be too large.
Fig. 17 Primordial spectra (scalar and tensor) of the bestfit V(φ) model with n = 4, for the Planck TT, TE, EE+lowP data set, compared to the primordial spectrum (scalar only) of the bestfit base ΛCDM model. The bestfit potential is initially very steep, as can be seen in Fig. 15 (note the typical shape of the green curves). The transition from “marginal slow roll” (ϵ_{V}(φ) between 0.01 and 1) to “full slow roll” (ϵ_{V}(φ) of order 0.01 or smaller) is responsible for the suppression of the largescale scalar spectrum. 
Fig. 18 Temperature and polarization spectra (total, scalar contribution, tensor contribution) of the bestfit V(φ) model with n = 4, for the Planck TT, TE, EE+lowP data set, compared to the spectra (scalar contribution only) of the bestfit base model. We also show the Planck lowℓ temperature data, which is driving the small differences between the two bestfit models. 
The n = 4 case possesses a new feature. The potential has more freedom to generate complicated shapes which would roughly correspond to a running of the running of the tilt (as studied in Sect. 4). The bestfit models now have V_{φ} and V_{φφφ} of opposite sign, and a large positive V_{φφφφ}. The preferred combination of these parameters allows for even more suppression of power on large scales, while leaving small scales nearly unchanged. This can be seen clearly from the shape of the scalar primordial spectrum corresponding to the bestfit models, for both data sets Planck TT+lowP+BAO and Planck TT, TE, EE+lowP. These two bestfit models are very similar, but in Fig. 17 we show the one for Planck TT, TE, EE+lowP, for which the trend is even more pronounced. Interestingly, the preferred models are such that power on large scales is suppressed in the scalar spectrum and balanced by a small tensor contribution, of roughly r_{0.002} ~ 0.05. This particular combination gives the best fit to the lowℓ data, shown in Fig. 18, while leaving the highℓ temperature spectrum identical to the best fit base ΛCDM model. Inflation produces such primordial perturbations with the family of green potentials displayed in Fig. 15. At the beginning of the observable range, the potential is very steep [ϵ_{V}(φ) decreases from O(1) to O(10^{2})], and produces a low amplitude of curvature perturbations (allowing a rather large tensor contribution, up to r_{0.002} ~ 0.3). Then there is a transition towards a second region with a much smaller slope, leading to a nearly powerlaw curvature spectrum with the usual tilt value n_{s} ≈ 0.96. In Fig. 15, one can check that the height of the n = 4 potentials varies in a definite range, while the n = 2 and 3 potentials can have arbitrarily small amplitude at the pivot scale, reflecting the posteriors on ϵ_{V} or r.
However, the improvement in χ^{2} between the base ΛCDM and n = 4 models is only 2.2 (for Planck TT+lowP+BAO) or 4.3 (for Planck TT, TE, EE+lowP). This is marginal and offers no statistically significant evidence in favour of these complicated models. This conclusion is also supported by the calculation of the Bayesian evidence ratios, shown in the last line of Table 8 (under the assumption of flat priors in the range [−1, 1] for and ): the evidence decreases each time that a new free parameter is added to the potential. At the 95% CL, r_{0.002} is still compatible with zero, and so are the higher order PSR parameters and . More freedom in the inflaton potential allows fitting the data better, but under the assumption of a smooth potential in the observable range, a simple quadratic form provides the best explanation of the Planck observations.
With the Planck TT+lowP+BAO and TT, TE, EE+lowP data sets, models with a large running or running of the running can be compatible with an unusually large value of the optical depth, as can be seen in Table 8. Including lensing information helps to break the degeneracy between the optical depth and the primordial amplitude of scalar perturbations. Hence the Planck lensing data could be used to strengthen the conclusions of this section.
Since in the n = 4 model, slow roll is marginally satisfied at the beginning of observable inflation, the reconstruction is very sensitive to the condition that there is an attractor solution at that time. Hence this case can in principle be investigated in a more conservative way using the H(φ) reconstruction method of the next section.
7.3. Reconstruction of a smooth Hubble function
In this section, we assume that the shape of the function H(φ) is well captured within the observable window by a polynomial of order n (corresponding to a polynomial inflaton potential of order 2n): (62)We vary n between 2 and 4. To avoid parameter degeneracies, as in the previous section we assume flat priors not on the Taylor coefficient H_{i}, but on the Hubble slowroll (HSR) parameters, which are related according to This is just a choice of prior. This analysis does not rely on the slowroll approximation.
Fig. 19 Posterior distributions for the first four Hubble slowroll parameters, when H(φ) is Taylor expanded to nth order, using Planck TT+lowP+BAO (filled contours) or TT, TE, EE+lowP (dashed contours). The primordial spectra are computed beyond the slowroll approximation. 
Numerical reconstruction of the Hubble slowroll parameters beyond the slowroll approximation, using Planck TT+lowP+BAO.
Fig. 20 Same as Fig. 15, when the Taylor expansion to nth order is performed on H(φ) instead of V(φ), and the potential is inferred from Eq. (60). 
Table 9 and Fig. 19 show our results for the reconstructed HSR parameters. Figure 20 shows a representative sample of potential shapes V(φ−φ_{∗}) derived using Eq. (60), for a sample of models drawn randomly from the chains, for the three cases n = 2,3,4.
Most of the discussion of Sect. 7.2 also applies to this section, and so will not be repeated. Results for Planck TT+lowP+BAO and TT, TE, EE+lowP are still very similar. The n = 2 case still gives results close to ΛCDM+r, and the n = 3 case to ΛCDM+r+dn_{s}/ dlnk. The type of potential preferred in the n = 4 case is very similar to the n = 4 analysis of the previous section, for the reasons explained in Sect. 7.2. There are, however, small differences, because the range of parametric forms for the potential explored by the two analyses differ. In the H(φ) reconstruction, the underlying potentials V(φ) are not polynomials. In the first approximation, they are close to polynomials of order 2n, but with constraints between the various coefficients. The main two differences with respect to the results of Sect. 7.2 are as follows:

The reconstructed potential shapes for n = 4 at the beginning of the observable window differ. Figure 20 shows that even steeper potentials are allowed than for the V(φ) method, with an even greater excursion of the inflaton field between Hubble crossing for the largest observable wavelengths and the pivot scale. This is because the H(φ) reconstruction does not rely on attractor solutions and automatically explores all valid potentials regardless of their initial field velocity.

The bestfit models are different, since they do not explore the same parametric families of potentials. In particular, for n = 4, the bestfit models have a negligible tensor contribution, but the distributions still have thick tails towards large tensortoscalar ratios, so that the upper bound on r_{0.002} is as high as in the previous n = 4 models, r_{0.002}< 0.32.
Note that can be significantly larger than unity for n = 4. This does not imply violation of slow roll within the observable range. By assumption, for all accepted models, ϵ_{H} must remain smaller than unity over that range. In fact, for most of the green potentials visible in Fig. 20, we checked that ϵ_{H} either has a maximum very close to unity near the beginning of the observable range or starts from unity. So the bestfit models (maximizing the power suppression at low multipoles) correspond either to inflation of short duration, or to models nearly violating slow roll just before the observable window. However, such peculiar models are not necessary for a good fit. Table 9 shows that the improvement in χ^{2} as n increases is negligible.
In summary, this section further establishes the robustness of our potential reconstruction and two main conclusions. Firstly, under the assumption that the inflaton potential is smooth over the observable range, we showed that the simplest parametric forms (involving only three free parameters including the amplitude V(φ_{∗}), no deviation from slow roll, and nearly power law primordial spectra) are sufficient to explain the data. No highorder derivatives or deviations from slow roll are required. Secondly, if one allows more freedom in the potential – typically, two more parameters – it is easy to decrease the largescale primordial spectrum amplitude with an initial stage of “marginal slow roll” along a steep branch of the potential followed by a transition to a less steep branch. This type of model can accommodate a large tensortoscalar ratio, as high as r_{0.002} ≈ 0.3.
8. reconstruction
In PCI13 (Sect. 7) we presented the results of a penalized likelihood reconstruction, seeking to detect any possible deviations from a homogeneous powerlaw form (i.e., ) for the primordial power spectrum (PPS) for various values of a smoothing parameter, λ. (For an extensive set of references to the prior literature concerning the methodology for reconstructing the power spectrum, see PCI13.) In the initial March 2013 preprint version of that paper, we reported evidence for a feature at moderate statistical significance around k ≈ 0.15 Mpc^{1}. However, in the November 2013 revision we retracted this finding, because subsequent tests indicated that the feature was no longer statistically significant when more aggressive cuts were made to exclude sky survey rings where contamination from electromagnetic interference from the 4 K cooler was largest, as indicated in the November 2013 “Note Added.”
In this section we report on results using the 2015 likelihood (Sect. 8.1) using essentially the same methodology as described in PCI13. (See Gauthier & Bucher 2012, and references therein for more technical details.) This method is also extended to include the EE and TE likelihoods in Sect. 8.1.2. As part of this 2015 release, we include the results of two other methods (see Sects. 8.2 and 8.3) to search for features. We find that all three methods yield broadly consistent reconstructions and reach the following main conclusion: there is no statistically significant evidence for any features departing from a simple powerlaw (i.e., ) PPS. Given the substantial differences between these methods, it is satisfying to observe this convergence.
8.1. Method I: penalized likelihood
8.1.1. Update with 2015 temperature likelihood
We repeated the same maximum likelihood analysis used to reconstruct the PPS in PCI13 using the updated Planck TT+lowP likelihood. Since we are interested in deviations from the nearly scaleinvariant model currently favoured by the parametric approach, we replaced the true PPS by a fiducial powerlaw spectrum , modulated by a small deviation function f(k): (65)The deviation function f(k)^{7} was represented by Bspline basis functions parameterized by n_{knot} control points , which are the values of f(k) along a grid of knot points κ_{i} = lnk_{i}.
Naively maximizing the Planck TT+lowP likelihood with respect to f results in overfitting to cosmic variance and noise in the data. Furthermore, due to the limited range of scales over which Planck measures the anisotropy power spectrum, the likelihood is very weakly dependent on f(k) at extremely small and large scales. To address these issues, the following two penalty functions were added to the Planck likelihood: (66)The first term on the righthand side of Eq. (66) is a roughness penalty, which disfavours f(κ) that “wiggle” too much. The last two terms drive f(κ) to zero for scales below κ_{min} and above κ_{max}. The values of λ and α represent the strengths of the respective penalties. The exact value of α is unimportant as long as it is large enough to drive f(κ) close to zero on scales outside [ κ_{min},κ_{max} ]. However, the magnitude of the roughness penalty λ controls the smoothness of the reconstruction.
Since the anisotropy spectrum depends linearly on the PPS, the NewtonRaphson method is well suited to optimizing with respect to f. However, a maximum likelihood analysis also has to take into account the cosmological parameters, Θ ≡ { H_{0},Ω_{b}h^{2},Ω_{c}h^{2} }^{8}. These additional parameters are not easy to include in the NewtonRaphson method since it is difficult to evaluate the derivatives ∂C_{ℓ}/∂Θ, ∂^{2}C_{ℓ}/∂Θ^{2}, etc., to the accuracy required by the method. Therefore a nonderivative method, such as the downhill simplex algorithm, is best suited to optimization over these parameters. Unfortunately the downhill simplex method is inefficient given the large number of control points in our parameter space. Since each method has its drawbacks, we combined the two methods to draw on their respective strengths. We define the function ℳ as (67)Given a set of nonPPS cosmological parameters Θ, ℳ is the value of the penalized log likelihood, minimized with respect to f using the NewtonRaphson method. The function ℳ is in turn minimized with respect to Θ using the downhill simplex method. In contrast to the analysis done in PCI13, the Planck lowℓ likelihood has been modified so that it can be included in the NewtonRaphson minimization. Thus the reconstructions presented here extend to larger scales than were considered in 2013.
Fig. 21 Planck TT+lowP likelihood primordial power spectrum (PPS) reconstruction results. Top four panels: reconstruction of the deviation f(k) using four different roughness penalties. The red curves represent the bestfit deviation f(k) using the Planck TT+lowP likelihood. f(k) = 0 would represent a perfectly featureless spectrum with respect to the fiducial PPS model, which is obtained from the bestfit base ΛCDM model with a powerlaw PPS. The vertical extent of the dark and light green error bars indicates the ± 1σ and ± 2σ errors, respectively. The width of the error bars represents the minimum reconstructible width (the minimum width for a Gaussian feature so that the mean square deviation of the reconstruction is less than 10%). The grey regions indicate where the minimum reconstructible width is undefined, indicating that the reconstruction in these regions is untrustworthy. The hatched region in the λ = 10^{6} plot shows where the fixing penalty has been applied. These hashed regions are not visible in the other three reconstructions, for which κ_{min} lies outside the range shown in the plots. For all values of the roughness penalty, all data points are within 2σ of the fiducial PPS except for the deviations around k ≈ 0.002 Mpc^{1} in the λ = 10^{3} and λ = 10^{4} reconstructions. Lower three panels: ± 1σ error bars of the three nonPPS cosmological parameters included in the maximum likelihood reconstruction. All values are consistent with their respective bestfit fiducial model values indicated by the dashed lines. 
Fig. 22 Planck TT, TE, EE+lowP likelihood primordial power spectrum reconstruction results. Top four panels: reconstruction of the deviation f(k) using four different roughness penalties. As in Fig. 21, the red curves represent the bestfit deviation f(k) and the height and width of the green error bars represent the error and minimum reconstructible width, respectively. For all values of the roughness penalty, the deviations are consistent with a featureless spectrum. Lower three panels: ± 1σ error bars of the three nonPPS cosmological parameters included in the maximum likelihood reconstruction. All values are consistent with their respective bestfit fiducial model values (indicated by the dashed lines). 
Figure 21 shows the bestfit PPS reconstruction using the Planck TT+lowP likelihood. The penalties in Eq. (66) introduce a bias in the reconstruction by smoothing and otherwise deforming potential features in the power spectrum. To assess this bias, we define the “minimum reconstructible width” (MRW) to be the minimum width of a Gaussian feature needed so that the integrated squared difference between the feature and its reconstruction is less than 1% of the integrated square of the input Gaussian, which is equivalent to 10% rms. Due to the combination of the roughness and fixing penalties, it is impossible to satisfy the MRW criterion too close to κ_{min} and κ_{max}. Wherever the MRW is undefined, the reconstruction is substantially biased and therefore suspect. An MRW cannot be defined too close to the endpoints κ_{min} and κ_{max} for two reasons: (1) lack of data; and (2) if a feature is too close to where the fixing penalty has been applied, the fixing penalty distorts the reconstruction. Consequently a larger roughness penalty decreases the range over which an MRW is well defined. The grey shaded areas in Fig. 21 show where the MRW is undefined and thus the reconstruction cannot be trusted. The cutoffs κ_{min} and κ_{max} have been chosen to maximize the range over which an MRW is defined for a given value of λ. The 1σ and 2σ error bars in Fig. 21 are estimated using the Hessian of the loglikelihood evaluated at the bestfit PPS reconstruction. More detailed MC investigations suggest that the nonlinear corrections to these error bars are small.
For the λ = 10^{5} and 10^{6} cases of the TT reconstruction, no deviation exceeds 2σ, so we do not comment on the probability of obtaining a worse fit. For the other cases, we use the maximum of the deviation, expressed in σ, of the plotted points as a metric of the quality of fit. Then using Monte Carlo simulations, we compute the pvalue, or the probability to obtain a worse fit, according to this metric. For λ = 10^{3} and 10^{4}, we obtain pvalues of 0.304 and 0.142, respectively, corresponding to 1.03 and 1.47σ. We thus conclude that the observed deviations are not statistically significant.
8.1.2. Penalized likelihood results with polarization
In Fig. 22 the bestfit reconstruction of the PPS from the Planck TT, TE, EE+lowP likelihood is shown. We observe that the reconstruction including polarization broadly agrees with the reconstruction obtained using temperature only. For the Planck TT, TE, EE+lowP likelihood, we obtain for λ = 10^{3}, 10^{4}, and 10^{5} the pvalues 0.166,0.107, and 0.045, respectively, corresponding to 1.38,1.61, and 2.00σ, and likewise conclude the absence of any statistically significant evidence for deviations from a simple powerlaw scalar primordial power spectrum.
8.2. Method II: Bayesian model comparison
In this section we model the PPS using a nested family of models where is piecewise linear in the ln(k) plane between a number of knots, N_{knots}, that is allowed to vary. The question arises as to how many knots one should use, and we address this question using Bayesian model comparison. A family of priors is chosen where both the horizontal and vertical positions of the knots are allowed to vary. We examine the “Bayes factor” or “Bayesian evidence” as a function of N_{knots} to decide how many knots are statistically justified. A similar analysis has been performed by Vázquez et al. (2012) and Aslanyan et al. (2014). In addition, we marginalize over all possible numbers of knots to obtain an averaged reconstruction weighted according to the Bayesian evidence.
The generic prescription is illustrated in Fig. 23. N_{knots} knots : i = 1,...,N_{knots} } are placed in the plane and the function is constructed by logarithmic interpolation (a linear interpolation in log log space) between adjacent points. Outside the horizontal range [ k_{1},k_{N} ] the function is extrapolated using the outermost interval.
Within this framework, base ΛCDM arises when N_{knots} = 2 – in other words, when there are two boundary knots and no internal knots, and the parameters and (in place of A_{s} and n_{s}) parameterize the simple powerlaw PPS. There are also, of course, the four standard cosmological parameters (Ω_{b}h^{2}, Ω_{c}h^{2}, 100θ_{MC}, and τ), as well as the numerous foreground parameters associated with the Planck highℓ likelihood, all of which are unrelated to the PPS. This simplest model can be extended iteratively by successively inserting an additional internal knot, thus requiring with each iteration two more variables to parameterize the new knot position.
Fig. 23 Linear spline reconstruction. The primordial power spectrum is reconstructed using N_{knots} interpolation points : i = 1,2,...,N_{knots} }. The end knots are fixed in k but allowed to vary in , whereas the internal knots can vary subject to the constraint that k_{1}<k_{2}< ··· <k_{Nknots}. The function is constructed within the range [ k_{1},k_{Nknots} ] by interpolating logarithmically between adjacent knots (i.e., linearly in log log space). Outside this range the function is extrapolated logarithmically. The function thus has 2N_{knots}−2 parameters. 
We run models for a variety of numbers of internal knots, N_{int} = N_{knots}−2, evaluating the evidence for N_{int}. Under the assumption that the prior is justified, the most likely, or preferred, model is the one with the highest evidence. Evidences are evaluated using the PolyChord sampler (Handley et al. 2015) in CAMB and CosmoMC. The use of PolyChord is essential, as the posteriors in this parameterization are often multimodal. Also, the ordered loguniform priors on the k_{i} are easy to implement within the PolyChord framework. All runs were performed with 1000 live points, oversampling the semislow and fast parameters by a factor of 5 and 100, respectively.
Prior for moveable knot positions.
Priors for the reconstruction parameters are detailed in Table 10. We report evidence ratios with respect to the base ΛCDM case. The cosmological priors remain the same for all models, and this part of the prior has almost no impact on the evidence ratios. The choice of prior on the reconstruction parameters does affect the Bayes factor. CosmoMC, however, puts an implicit prior on all models by excluding parameter choices that render the internal computational approximations in CAMB invalid. The baseline prior for the vertical position of the knots includes all of the range allowed by CosmoMC, so slighly increasing this prior range will not affect the evidence ratios. If one were to reduce the prior widths significantly, the evidence ratios would be increased. The allowed horizontal range includes all kscales accessible to Planck. Thus, altering this width would be unphysical.
After completion of an evidence calculation, PolyChord generates a representative set of samples of the posterior for each model, P(Θ) ≡ P(Θ  data,N_{int}). We may use this to calculate a marginalized probability distribution for the PPS: (68)This expression encapsulates our knowledge of at each value of k for a given number of knots. Plots of this PPS posterior are shown in Fig. 24 using Planck TT data.
If one considers the Bayesian evidence of each model, Fig. 25 shows that although no model is preferred over base ΛCDM, the case N_{int} = 1 is competitive. This model is analogous to the brokenpowerlaw spectrum of Sect. 4.4, although the models differ significantly in terms of the priors used. In this case, the additional freedom of one knot allows a reconstruction of the suppression of power at low ℓ. Adding polarization data does not alter the evidences significantly, although N_{int} = 1 is strengthened. We also plot a Planck TT run, but with the reduced vertical priors . As expected, this increases the evidence ratios, but does not alter the above conclusion.
For increasing numbers of internal knots, the Bayesian evidence monotonically decreases. Occam’s razor dictates, therefore, that these models should not be preferred, due to their higher complexity. However, there is an intriguing stable oscillatory feature, at 20 ≲ ℓ ≲ 50, that appears once there are enough knots to reconstruct it. This is a qualitative feature predicted by several inflationary models (discussed in Sect. 9), and a possible hint of new physics, although its statistical significance is not compelling.
A full Bayesian analysis marginalizes over all models weighted according to the normalized evidence Z_{Nint}, so that (69)as indicated in Fig. 26. This reconstruction is sensitive to how model complexity is penalized in the prior distribution.
8.3. Method III: cubic spline reconstruction
In this section we investigate another reconstruction algorithm based on cubic splines in the ln(k) plane, where (unlike for the approach of the previous subsection) the horizontal positions of the knots are uniformly spaced in ln(k) and fixed. A prior on the vertical positions (described in detail below) is chosen and the reconstructed power spectrum is calculated using CosmoMC for various numbers of knots. This method differs from the method in Sect. 8.1 in that the smoothness is controlled by the number of discrete knots rather than by a continuous parameter of a statistical model having a welldefined continuum limit. With respect to the Bayesian model comparison of Sect. 8.2, the assessment of model complexity differs because here the knots are not movable.
Fig. 24 Bayesian movable knot reconstructions of the primordial power spectrum using Planck TT data. The plots indicate our knowledge of the PPS for a given number of knots. The number of internal knots N_{int} increases (left to right and top to bottom) from 0 to 8. For each kslice, equal colours have equal probabilities. The colour scale is chosen so that darker regions correspond to lowerσ confidence intervals. 1σ and 2σ confidence intervals are also indicated (black curves). The upper horizontal axes give the approximate corresponding multipoles via ℓ ≈ kD_{rec}, where D_{rec} is the comoving distance to recombination. 
Let the horizontal positions of the n knots be given by k_{b}, where b = 1,...,n, spaced so that k_{b + 1}/k_{b} is independent of b. We single out a “pivot knot” b = p, so that k_{p} = k_{∗} = 0.05 Mpc^{1}, which is the standard scalar power spectrum pivot scale. For a given number of knots n we choose k_{1} and k_{n} so that the interval of relevant cosmological scales, taken to extend from 10^{4} Mpc^{1} to O(1) Mpc^{1}, is included. We now define the prior on the vertical knot coordinates. For the pivot point, we define , where lnA_{s} has a uniformly distributed prior, and for the other points with b ≠ p, we define the derived variable (70)where Here the spectral index n_{s,fid} is fixed. A uniform prior is imposed on each variable q_{b}(b ≠ p) and the constraint −1 ≤ q_{b} ≤ 1 is also imposed to force the reconstruction to behave reasonably near the endpoints, where it is hardly constrained by the data. The quantity is interpolated between the knots using cubic splines with natural boundary conditions (i.e., the second derivatives vanish at the first and the last knots). Outside [ k_{1},k_{n} ] we set (for k<k_{1}) and (for k>k_{n}). For most knots, we found that the upper and lower bounds of the q_{b} prior hardly affect the reconstruction, since the data sharpen the allowed range significantly. However, for superHubble scales (i.e., k ≲ 10^{4} Mpc^{1}) and very small scales (i.e., k ≳ 0.2 Mpc^{1}), which are only weakly constrained by the cosmological data, the prior dominates the reconstruction. For the results here, a fiducial spectral index n_{s,fid} = 0.967 for was chosen, which is close to the estimate from Planck TT+lowP+BAO. A different choice of n_{s,fid} leads to a trivial linear shift in the q_{b}.
Fig. 25 Bayes factor (relative to the base ΛCDM model) as a function of the number of knots for three separate runs. Solid line: Planck TT. Dashed line: Planck TT, TE, EE. Dotted line: Planck TT, with priors on the parameters reduced in width by a factor of 2 (). 
The possible presence of tensor modes (see Sect. 5) has the potential to bias and introduce additional uncertainty in the reconstruction of the primordial scalar power spectrum as parameterized above. Obviously, in the absence of a detection of tensors at high statistical significance, it is not sensible to model a possible tensor contribution with more than a few degrees of freedom. A complicated model would lead to prior dominated results. We therefore use the power law parameterization, , where the consistency relation n_{t} = −r/ 8 is enforced as a constraint.
Fig. 26 Bayesian reconstruction of the primordial power spectrum averaged over different values of N_{int} (as shown in Fig. 24), weighted according to the Bayesian evidence. The region 30 <ℓ< 2300 is highly constrained, but the resolution is lacking to say anything precise about higher ℓ. At lower ℓ, cosmic variance reduces our knowledge of The weights assigned to the lower N_{int} models outweigh those of the higher models, so no oscillatory features are visible here. 
Fig. 27 Reconstructed power spectra applied to the Planck 2015 data using 12 knots (with positions marked as Δ at the bottom of each panel) with cubic spline interpolation. Mean spectra as well as sample trajectories are shown for scalars and tensors, and ± 1σ and ± 2σ limits are shown for the scalars. The fiducial tensor spectrum corresponds, arbitrarily, to r = 0.13. Top: uniform prior, 0 ≤ r ≤ 1. Middle: fixed, r = 0.1. Bottom: fixed, r = 0.01. Data sets: Planck TT+lowP+BAO+SN+HST+z_{re} > 6 prior. D_{rec} is the comoving distance to recombination. 
Primordial tensor fluctuations contribute to CMB temperature and polarization angular power spectra, in particular at spatial scales larger than the recombination Hubble length, k ≲ (aH)_{rec} ≈ 0.005 Mpc^{1}. If a large number of knots in is included over that range, then a modified can mimic a tensor contribution, leading to a neardegeneracy. This can lead to large uncertainty in the tensor amplitude, r. Once r is measured or tightly constrained in Bmode experiments, this near degeneracy will be broken. As examples here, we do allow r to float, but also show what happens when r is constrained to take the values r = 0.1 and r = 0.01 in the reconstruction.
Figure 27 shows the reconstruction obtained using the 2015 Planck TT+lowP likelihood, BAO, SNIa, HST, and a z_{re} > 6 prior. Including these ancillary likelihoods improves the constraint on the PPS by helping to fix the cosmological parameters (e.g., H_{0}, τ, and the latetime expansion history), which in this context may be regarded as nuisance parameters. These results were obtained by modifying CosmoMC to incorporate the nknot parameterization of the PPS. Here 12 knots were used and the mean reconstruction as well as the 1σ and 2σ limits are shown. Some 1σ sample trajectories (dashed curves) are also shown to illustrate the degree of correlation or smoothing of the reconstruction. The tensor trajectories are also shown, but, as explained above, have been constrained to be straight lines. In the top panel r is allowed to freely float, and a wide range of r is allowed because of the neardegeneracy with the lowk scalar power. Two illustrative values of fixed r (i.e., r = 0.1 and r = 0.01) are also shown to give an idea of how much the reconstruction is sensitive to variations in r within the range of presently plausible values.
Fig. 28 Reconstructed 12knot power spectra with polarization included. Data sets in common: lowT+lowP+BAO+SN+HST+z_{re} > 6 prior. 
The reconstructions using the 2013 Planck likelihood in place of the 2015 likelihood are broadly consistent with the reconstruction shown in Fig. 27. To demonstrate robustness with respect to the interpolation scheme we tried using linear interpolation instead of cubic splines and found that the reconstruction was consistent provided enough knots (i.e., n_{knot} ≈ 14) were used. At intermediate k the reconstruction is consistent with a simple power law, corresponding to a straight line in Fig. 27. We observe that once k drops, so that the effective multipole being probed is below about 60, deviations from a power law appear, but the dispersion in allowed trajectories also rises as a consequence of cosmic variance. The power deficit at k ≈ 0.002 Mpc^{1} (i.e., ℓ_{k} ≡ kD_{rec} ≈ 30, where D_{rec} is the comoving distance to recombination) is largely driven by the power spectrum anomaly in the ℓ ≈ 20−30 range that has been evident since the early spectra from WMAP (Bennett et al. 2011), and verified by Planck.
Fig. 29 Reconstructed power spectra with the base ΛCDM best fit subtracted. The mean spectra shown are for the floating r and the two fixed r cases with 12 cubic spline knots. These should be contrasted with the running bestfit mean (green) and the similar looking uniform n_{s} case in which τ has been lowered from its bestfit base ΛCDM value to 0.04. Data points are the Planck 2015 Commander (ℓ< 30) and Plik (ℓ ≥ 30) temperature power spectrum. 
We also explore the impact of including the Planck polarization likelihood in the reconstruction. Figure 28 shows the reconstructed power spectra using various combinations of the polarization and temperature data. The ℓ < 30 treatments are the same in all cases, so this is mainly a test of the higher k region. What is seen is that, except at high k, the EE polarization data also enforce a nearly uniform n_{s}, consistent with that from TT, over a broad krange. When TE is used alone, or TE and EE are used in combination, the result is also very similar. The upper right panel shows the constraints from all three spectra together, and the errors on the reconstruction are now better than those from TT alone.
It is interesting to examine how the TT power spectrum obtained using the above reconstructions compares to the CMB data, in particular around the range ℓ ≈ 20–30, corresponding roughly to k_{4} ≈ 1.5 × 10^{3} Mpc^{1}. In Fig. 29 the differences in from the bestfit simple powerlaw model are plotted for various assumptions concerning r. We see that a better fit than the powerlaw model can apparently be obtained around ℓ ≈ 20−30. We quantify this improvement below.
Due to the degeneracy of scalar and tensor contributions to , the significance of the lowℓ anomaly depends on the tensor prior and whether polarization data are used. For k < 10^{3} Mpc^{1}, once more degeneracy appears: the shape of also depends on the reionization optical depth, τ. In Fig. 29 we also show the effect of replacing the bestfit τ for tilted base ΛCDM with a low value, while keeping A_{s}e^{− 2τ} unchanged. A low τ bends downward at ℓ ≲ 10. For the 12knot (or similar) runs, if τ is allowed to run into the (nonphysically) small values τ ≲ 0.04, a slight rise in at k ≈ 3 × 10^{4} Mpc^{1} is preferred to compensate the lowτ effect. This degeneracy can be broken to a certain extent using lowredshift data: z_{re} > 6 from quasar observations (Becker et al. 2001), BAO (SDSS), Supernova (JLA), and HST.
It is evident that allowing n_{s} to run is not what the data prefer. The bestfit running is also shown in Fig. 29. The kspace response in Fig. 27 shows that running does not capture the shape of the lowℓ residuals.
Reduced χ^{2} and pvalues for lowk knots (5 knots) and highk knots (6 knots, pivot knot excluded), with the null hypothesis being the bestfit powerlaw spectrum.
We have shown that the cubic spline reconstruction studied in this section consistently produces a dip in q_{4}, corresponding to k ≈ 1.5 × 10^{3} Mpc^{1}. We now turn to the question of whether this result is real or simply the result of cosmic variance. To assess the statistical significance of the departures of the mean reconstruction from a simple power law, we calculate the lowk and highk reduced χ^{2} for the five q_{b} values for scales below and six q_{b} values (b ≠ p) for scales above 50 /D_{rec}, respectively, indicating the corresponding pvalues (i.e., probability to exceed), for various data combinations, in Table 11. The highk fit is better than expected for reasons that we do not understand, but we attribute this situation to chance. The lowk region shows a poor fit, but in no case does the pvalue fall below 10%. Therefore, even though the lowk dip is robust against the various choices made for the reconstruction, we conclude that it is not statistically significant. The plot for the knot position of the dip (corresponding to q_{4}) in Fig. 30 does not contradict this conclusion.
Fig. 30 The degeneracy between τ and the knot variables q_{3} and q_{4} in the 12knot case shown in Fig. 27. 
Fig. 31 Slowroll parameter ϵ for reconstructed trajectories using 12 knots (marked as Δ at the bottom of the figure) with cubic spline interpolation. The mean values are shown for floating r and r fixed to be 0.1 and 0.01. Sample 1σ trajectories shown for the floating r case show wide variability, which is significantly diminished if r is fixed to r = 0.1, as shown. 
Fig. 32 Reconstructed singlefield inflaton potentials from the cubic spline power spectra mode expansion using 12 knots. 
Because of the r degeneracy associated with the scalar power, it is best when quoting statistics to use the fixed r cases, although for completeness we show the floating r case as well. There is also a smaller effect associated with the τ degeneracy, and the values quoted have restricted the redshift of reionization to exceed 6. The value z_{re} = 6.5 was used in Planck Collaboration XIII (2016). The significance of the lowk anomaly is meaningful only if an explicit r prior and lowredshift constraint on τ have been applied.
Finally, we relate the reconstructed calculated above to the trajectories of the slowroll parameter ϵ = −Ḣ/H^{2}  _{k = aH} plotted as a function of k (see Fig. 31). We also plot in Fig. 32 the reconstructed inflationary potential in the region over which the inflationary potential is constrained by the data. Here canonical singlefield inflation is assumed, and the value of r enters solely to fix the height of the potential at the pivot scale. This is not entirely selfconsistent, but justified by the lack of constraining power on the tensors at present.
8.4. Power spectrum reconstruction summary
The three nonparametric methods for reconstructing the primordial power spectrum explored here support the following two conclusions:

1.
Except possibly at low k, over the range of k where the CMB data best constrain the form of the primordial power spectrum, none of the three methods finds any statistically significant evidence for deviations from a simple powerlaw form. The fluctuations seen in this regime are entirely consistent with the expectations from cosmic variance and noise.

2.
At low k, all three methods reconstruct a power deficit at k ≈ 1.5–2.0 × 10^{3} Mpc^{1}, which can be linked to the dip in the TT angular power spectrum at ℓ ≈ 20−30. This agreement suggests that the reconstruction of this “anomaly” is not an artefact of any of the methods, but rather inherent in the CMB data themselves. However, the evidence for this feature is marginal since it is in a region of the spectrum where the fluctuations from cosmic variance are large.

3.
We have verified that the power deficit at ℓ = 20–30 is not substantially modified (a) by removing from the CMB pattern the hottest and coldest peaks selected by the KolmogorovSmirnov test studied in Sects. 4.5.3 and 4.5.4 of Planck Collaboration XVI (2016) or (b) by substituting the anomalously cold region around the Cold Spot with Gaussian constrained realizations.
9. Search for parameterized features
In this section, we explore the possibility of a radical departure from the nearscaleinvariant powerlaw spectrum of the standard slowroll scenario for a selection of theoretically motivated parameterizations of the spectrum (see Chluba et al. 2015 for a recent review).
9.1. Models
9.1.1. Step in the inflaton potential
A sudden, steplike feature in the inflaton potential (Adams et al. 2001) or the sound speed (Achúcarro et al. 2011) leads to a localized oscillatory burst in the scalar primordial power spectrum. A general parameterization describing both a tanhstep in the potential and in the warp term of a DBI model was proposed in Miranda & Hu (2014): (71)where the first and secondorder terms are given by with window functions and damping function (80)Due to the high complexity of this model, we focus on the limiting case of a step in the potential ().
9.1.2. Logarithmic oscillations
Logarithmic modulations of the primordial power spectrum generically appear, for example, in models with nonBunchDavies initial conditions (Martin & Brandenberger 2001; Danielsson 2002; Bozza et al. 2003), or, approximately, in the axion monodromy model, explored in more detail in Sect. 10. We assume a constant modulation amplitude and use (81)
9.1.3. Linear oscillations
A modulation linear in k can be obtained, for example, in boundary effective field theory models (Jackson & Shiu 2013), and is typically accompanied by a scaledependent modulation amplitude. We adopt the parameterization used in Meerburg & Spergel (2014), which allows for a strong scale dependence of the modulation amplitude: (82)
9.1.4. Cutoff model
If today’s largest observable scales exited the Hubble radius before the inflaton field reached the slowroll attractor, the amplitude of the primordial power spectrum is typically strongly suppressed at low k. As an example of such a model, we consider a scenario in which slow roll is preceded by a stage of kinetic energy domination. The resulting power spectrum was derived by Contaldi et al. (2003) and can be expressed as (83)with where denotes the Hankel function of the second kind. The power spectrum in this model is exponentially suppressed for wavenumbers smaller than the cutoff scale k_{c} and converges to a standard powerlaw spectrum for k ≫ k_{c}, with an oscillatory transition region for k ≳ k_{c}.
Parameters and prior ranges.
9.2. Analysis and results
Improvement in fit and Bayes factors with respect to powerlaw base ΛCDM for Planck TT+lowP and Planck TT, TE, EE+lowP data, as well as approximate probability to exceed the observed Δχ^{2} (pvalue), constructed from simulated Planck TT+lowP data.
Bestfit features parameters and parameter constraints on the remaining cosmological parameters for the four features models for Planck TT+lowP data.
We use MultiNest to evaluate the Bayesian evidence for the models, establish parameter constraints, and roughly identify the global maximum likelihood region of parameter space. The features model bestfit parameters and lnℒ are then obtained with the help of the CosmoMC minimization algorithm taking narrow priors around the MultiNest best fit. We assign flat prior probabilities to the parameters of the features models with prior ranges listed in Table 12. Note that throughout this section for the sake of maximizing sensitivity to very sharp features, the unbinned (“bin1”) versions of the highℓ TT and TT, TE, EE likelihoods are used instead of the standard binned versions.
Since the features considered here can lead to broad distortions of the CMB angular power spectrum degenerate with the late time cosmological parameters (Miranda & Hu 2014), in all cases we simultaneously vary primordial parameters and all the ΛCDM parameters, but keep the foreground parameters fixed to their bestfit values for the powerlaw base ΛCDM model.
We present the Bayes factors with respect to the powerlaw base ΛCDM model and the improvement in the effective χ^{2} over the powerlaw model in Table 13. For our choice of priors, none of the features models is preferred over a powerlaw spectrum. The bestfit power spectra are plotted in Fig. 33. While the cutoff and step model best fits reproduce the largescale suppression at ℓ ≈ 20−30 also obtained by direct power spectrum reconstruction in Sect. 8, the oscillation models prefer relatively high frequencies beyond the resolution of the reconstruction methods.
In addition to the four features models we also show in Fig. 33 the best fit of a model allowing for steps in both inflaton potential and warp (brown line). Note the strong resemblance to the reconstructed features of the previous section. The effective Δχ^{2} for this model is −12.1 (−11.5) for Planck TT+lowP (Planck TT, TE, EE+lowP) data at the cost of adding five new parameters, resulting in a lnBayes factor of −0.8 (−0.4). A similar phenomenology can be also be found for a model with a sudden change in the slope of the inflaton potential (Starobinsky 1992; Choe et al. 2004), which yields a bestfit Δχ^{2} = −4.5 (−4.9) for two extra parameters.
Fig. 33 Bestfit power spectra for the powerlaw (black curve), step (green), logarithmic oscillation (blue), linear oscillation (orange), and cutoff (red) models using Planck TT+lowP data. The brown curve is the best fit for a model with a step in the warp and potential (Eqs. (71)−(80)). 
As shown in Table 14, constraints on the remaining cosmological parameters are not significantly affected when allowing for the presence of features.
For the cutoff and step models, the inclusion of Planck smallscale polarization data does not add much in terms of direct sensitivity. The best fits lie in the same parameter region as for Planck TT+lowP data, and the Δχ^{2} and Bayes factors are not subject to major changes. The two oscillation models’ Planck TT+lowP best fits, on the other hand, also predict a nonnegligible signature in the polarization spectra at high ℓ. Therefore, if the features were real, one would expect an additional improvement in Δχ^{2} for Planck TT, TE, EE+lowP. This is not the case here. Though the linear oscillation model’s maximum Δχ^{2} does increase, the local Δχ^{2} in the Planck TT+lowP bestfit regions is in fact reduced for both models, and the global likelihood maxima occur at different frequencies (log _{10}ω_{log } = 1.25 and log _{10}ω_{lin} = 1.02) compared to their Planck TT+lowP counterparts.
In addition to the Bayesian model comparison analysis, we also approach the matter of the statistical relevance of the features models from a frequentist statistics perspective in order to give the Δχ^{2} numbers a quantitative interpretation. Assuming that the underlying was actually a featureless power law, we can ask how large an improvement to lnℒ the different features models would yield on average just by overfitting scatter from cosmic variance and noise. For this purpose, we simulate Planck power spectrum data sets consisting of temperature and polarization up to ℓ = 29 and unbinned temperature for 30 ≤ ℓ ≤ 2508, taking as input fiducial spectra the powerlaw base ΛCDM model’s bestfit spectra.
For each of these simulated Planck data sets, we perform the following procedure: (i) find the powerlaw ΛCDM model’s bestfit parameters with CosmoMC’s minimization algorithm; (ii) fix the nonprimordial parameters (ω_{b},ω_{c},θ_{MC},τ) to their respective bestfit values; (iii) using MultiNest, find the best fit of the features models;^{9} and (iv) extract the effective Δχ^{2} between powerlaw and features models.
The resulting distributions are shown in Fig. 34. Compared to the real data Δχ^{2} values from Table 13, they are biased towards lower values, since we do not vary the latetime cosmological parameters in the analysis of the simulated data. Nonetheless, the observed improvements in the fit do not appear to be extraordinarily large, with the respective (conservative) pvalues ranging between 0.09 and 0.50.
These observations lead to the conclusion that even though some of the peculiarities seen in the residuals of the Planck data with respect to a powerlaw primordial spectrum may be explained in terms of primordial features, none of the simple model templates considered here is required by Planck data. The simplicity of the powerlaw spectrum continues to give it an edge over more complicated initial spectra and the most plausible explanation for the apparent features in the data remains that we are just observing fluctuations due to cosmic variance at large scales and noise at small scales.
Fig. 34 Distribution of Δχ^{2} from 400 simulated Planck TT+lowP data sets. 
10. Implications of Planck bispectral constraints on inflationary models
The combination of power spectrum constraints and primordial nonGaussianity (NG) constraints, such as the Planck upper bound on the NG amplitude f_{NL} (Planck Collaboration XVII 2016), can be exploited to limit extensions to the simplest standard singlefield models of slowroll inflation. The next subsection considers inflationary models with a nonstandard kinetic term (Garriga & Mukhanov 1999), where the inflaton Lagrangian is a general function of the scalar inflaton field and its first derivative, i.e., L= P(φ,X), where X = −g^{μν}∂_{μ}φ∂_{ν}φ/ 2 (Garriga & Mukhanov 1999; Chen et al. 2007). Section 10.2 focuses on a specific example of a singlefield model of inflation with more general higherderivative operators, the socalled “Galileon inflation”. Section 10.3 presents constraints on axion monodromy inflation. See Planck Collaboration XVII (2016) for the analysis of other interesting nonstandard inflationary models, including warm inflation (Berera 1995), whose f_{NL} predictions can be constrained by Planck.
10.1. Inflation with a nonstandard kinetic term
This class of models includes kinflation (ArmendárizPicón 1999; Garriga & Mukhanov 1999) and DiracBornInfield (DBI) models introduced in the context of brane inflation (Silverstein & Tong 2004; Alishahiha et al. 2004; Chen 2005b,a). In these models inflation can take place despite a steep potential or may be driven by the kinetic term.
Moreover, one of the main predictions of inflationary models with a nonstandard kinetic term is that the inflaton perturbations can propagate with a sound speed c_{s}< 1. We show how the Planck combined measurement of the power spectrum and the nonlinearity parameter f_{NL} (Planck Collaboration XVII 2016) improves constraints on this class of models by breaking degeneracies between the parameters determining the observable power spectra. Such degeneracies (see, e.g., Peiris et al. 2007; Powell et al. 2009; Lorenz et al. 2008; Agarwal & Bean 2009; Baumann et al. 2015) are evident from the expressions for the power spectra. We adopt the same notation as Planck Collaboration XXIV (2014). At leading order in the slowroll parameters the scalar power spectrum depends additionally on the sound speed c_{s} via (Garriga & Mukhanov 1999) (86)which is evaluated at kc_{s} = aH. Correspondingly, the scalar spectral index (87)depends on an additional slowroll parameter , which describes the running of the sound speed. The usual consistency relation holding for the standard singlefield models of slowroll inflation (r = −8n_{t}) is modified to r ≈ −8n_{t}c_{s}, with n_{t} = −2ϵ_{1} as usual (Garriga & Mukhanov 1999), potentially allowing models which otherwise would predict a large tensortoscalar ratio for the KleinGordon case (Unnikrishnan et al. 2012).^{10}
At lowest order in the slowroll parameters, there are strong degeneracies between the parameters (A_{s},c_{s},ϵ_{1},ϵ_{2},s). This makes the constraints on these parameters from the power spectrum alone not very stringent, and for parameters like ϵ_{1} and ϵ_{2} less stringent compared with the standard case. However, combining the constraints on the power spectra observables with those on f_{NL} can also result in a stringent test for this class of inflationary models. Models where the inflaton field has a nonstandard kinetic term predict a high level of primordial NG of the scalar perturbations for c_{s} ≪ 1, (see, e.g., Chen et al. 2007). Primordial NG is generated by the higherderivative interaction terms arising from the expansion of the kinetic part of the Lagrangian, P(φ,X). There are two main contributions to the amplitude of the NG (i.e., to the nonlinearity parameter f_{NL}), coming from the inflaton field interaction terms and (Chen et al. 2007; Senatore et al. 2010). The NG from the first term scales as , while the NG arising from the other term is determined by a second parameter, (following the notation of Senatore et al. 2010). Each of these two interactions produces bispectrum shapes similar to the socalled equilateral shape (Babich et al. 2004) for which the signal peaks for equilateral triangles with k_{1} = k_{2} = k_{3}. (These two shapes are called, respectively, “EFT1” and “EFT2” in Planck Collaboration XVII 2016). However, the difference between the two shapes is such that the total signal is a linear combination of the two, leading to an “orthogonal” bispectral template (Senatore et al. 2010).
Fig. 35 (ϵ_{1},ϵ_{2}) 68% and 95% CL constraints for Planck data comparing the canonical Lagrangian case with c_{s} = 1 to the case of varying c_{s} with a uniform prior 0.024 <c_{s}< 1 derived from the Planck NG measurements. 
The equilateral and orthogonal NG amplitudes can be expressed in terms of the two “microscopic” parameters, c_{s} and (for more details see Planck Collaboration XVII 2016), according to Thus the measurements of and obtained in the companion paper (Planck Collaboration XVII 2016) provide a constraint on the sound speed, c_{s}, of the inflaton field. Such constraints allow us to combine the NG information with the analyses of the power spectra, since the sound speed is the NG parameter also affecting the power spectra.
In this subsection we consider three cases. In the first case we perform a general analysis as described above (focusing on the simplest case of a constant sound speed, s = 0), improving on PCI13 and Planck Collaboration XXIV (2014) by exploiting the full mission temperature and polarization data. The Planck constraints on primordial NG in general singlefield models of inflation provide the most stringent bound on the inflaton sound speed (Planck Collaboration XVII 2016):^{11}(91)We then use this information on c_{s} as a uniform prior 0.024 ≤ c_{s} ≤ 1 in Eq. (88) within the HFF formalism, as in PCI13. Figure 35 shows the joint constraints on ϵ_{1} and ϵ_{2}.Planck TT+lowP yields ϵ_{1}< 0.031 at 95% CL. No improvement in the upper bound on ϵ_{1} results when using Planck TT, TE, EE+lowP. This constraint improves the previous analysis in PCI13 and can be compared with the restricted case of c_{s} = 1, also shown in Fig. 35, with ϵ_{1}< 0.0068 at 95% CL. The limits on the sound speed from the constraints on primordial NG are crucial for deriving an upper limit on ϵ_{1}, because the relation between the tensortoscalar ratio and ϵ_{1} also involves the sound speed (see, e.g., Eq. (88)). This breaks the degeneracy in the scalar spectral index.
The other two cases analysed involve DBI models. The degeneracy between the different slowroll parameters can be broken for s = 0 or in the case where s ∝ ϵ_{2}. We first consider models defined by an action of the DBI form (92)where V(φ) is the potential and f(φ) describes the warp factor determined by the geometry of the extra dimensions. We follow an analogous procedure to exploit the NG limits derived in Planck Collaboration XVII (2016) on c_{s} in the case of DBI models: c_{s} ≥ 0.087 (at 95% CL). Assuming a uniform prior, 0.087 ≤ c_{s} ≤ 1, and s = 0, Planck TT+lowP gives ϵ_{1}< 0.024 at 95% CL, a 43% improvement with respect to PCI13. The addition of highℓ TE and EE does not improve the upper bound on ϵ_{1} for this DBI case.
Next we update the constraints on the particularly interesting case of infrared DBI models (Chen 2005b,a), where f(φ) ≈ λ/φ^{4}. (For details, see Silverstein & Tong 2004; Alishahiha et al. 2004; Chen et al. 2007, and references therein.) In these models the inflaton field moves from the IR to the UV side with an inflaton potential (93)From a theoretical point of view a wide range of values for β is allowed: 0.1 <β< 10^{9} (Bean et al. 2008). PCI13 dramatically restricted the allowed parameter space of these models in the limit where stringy effects can be neglected and the usual field theory computation of the primordial curvature perturbation holds (see Chen 2005a,c; Bean et al. 2008 for more details). In this limit of the IR DBI model, one finds (Chen 2005c; Chen et al. 2007) c_{s} ≈ (βN_{∗}/ 3)^{1}, n_{s}−1 = −4 /N_{∗}, and (In this model one can verify that s ≈ 1 /N_{∗} ≈ ϵ_{2}/ 3.) Combining the uniform prior on c_{s} with Planck TT+lowP, we obtain (94)and a preference for a high number of efolds: 78 <N_{∗}< 157 at 95% CL.
We now constrain the general case of the IR DBI model, including the “stringy” regime, which occurs when the inflaton extends back in time towards the IR side (Bean et al. 2008). The stringy phase transition is characterized by an interesting phenomenology altering the predictions for cosmological perturbations. A parameterization of the power spectrum of curvature perturbations interpolating between the two regimes is (Bean et al. 2008; see also Ma et al. 2013) (95)where A_{s} = 324π^{2}/ (n_{B}β^{4}) is the amplitude of the perturbations which depends on various microscopic parameters (n_{B} is the number of branes at the Bthroat; see Bean et al. 2008 for more details), while sets the stringy phase transition taking place at the critical efold N_{c}. (Here is the number of efolds to the end of IR DBI inflation.) The spectral index and its running are A prediction for the primordial NG in the stringy regime is not available. We assume the standard fieldtheoretic result for a primordial bispectrum of the equilateral type with an amplitude . By considering the same uniform prior on c_{s}, we obtain β< 0.77, , and x< 0.41 at 95% CL, which severely limits the general IR DBI model and strongly restricts the allowed parameter space.
10.2. Galileon inflation
As a further example of the implications of the NG constraints on (nonstandard) inflationary models we consider Galileon inflation Burrage et al. (2011; see also Kobayashi et al. 2010; Mizuno & Koyama 2010; Ohashi & Tsujikawa 2012). This represents a welldefined and wellmotivated model of inflation with more general higher derivatives of the inflaton field compared to the nonstandard kinetic term case analysed above. The Galileon models of inflation are based on the socalled “Galilean symmetry” (Nicolis et al. 2009), and enjoy some well understood stability properties (absence of ghost instabilities and protection from large quantum corrections). This makes the theory also very predictive, since observable quantities (scalar and tensor power spectra and higherorder correlators) depend on a finite number of parameters. From this point of view this class of models shares some of the same properties as the DBI inflationary models (Silverstein & Tong 2004; Alishahiha et al. 2004). The Galileon field arises naturally within fundamental physics constructions (e.g., de Rham & Gabadadze 2010b,a). These models also offer an interesting example of largescale modifications to Einstein gravity.
The Galileon model is based on the action (Deffayet et al. 2009a,b) (98)where Here X = −∇_{μ}φ∇^{μ}φ/ 2, (∇_{μ}∇_{ν}φ)^{2} = ∇_{μ}∇_{ν}φ∇^{μ}∇^{ν}φ, and (∇_{μ}∇_{ν}φ)^{3} = ∇_{μ}∇_{ν}φ∇^{μ}∇^{ρ}φ∇^{ν}∇_{ρ}φ. The coupling coefficients c_{i} are dimensionless and Λ is the cutoff of the theory. The case of interest includes a potential term V(φ) = V_{0} + λφ + (1 / 2)m^{2}φ^{2} + ... to drive inflation.
The predicted scalar power spectrum at leading order is (Ohashi & Tsujikawa 2012; Burrage et al. 2011; Tsujikawa et al. 2013; see also Kobayashi et al. 2011a; Gao & Steer 2011)^{12}(104)where and is the sound speed of the Galileon field. ϵ_{s} is different from the usual slowroll parameter ϵ_{1} and at leading order related according to The scalar spectral index (105)depends on the slowroll parameters ϵ_{1}, , and s = ċ_{s}/ (Hc_{s}). As usual the slowroll parameter s describes the running of the sound speed. In the following we restrict ourselves to the case of a constant sound speed with s = 0. The tensortoscalar ratio is (106)where we have introduced the parameter , which is related to the Galileon sound speed. The parameter can be either positive or negative. In the negative branch a blue spectral tilt for the primordial gravitational waves is allowed, contrary to the situation for standard slowroll models of inflation. We introduce such a quantity so that the consistency relation takes the form , with n_{t} = −2ϵ_{1}, analogous to Eq. (88). The measurements of primordial NG constrain which in turn constrains ϵ_{1} and η_{s} in Eq. (105). This is analogous to the constraints on ϵ_{1} and η of Eq. (87) in the previous subsection.
Galileon models of inflation predict interesting NG signatures (Burrage et al. 2011; Tsujikawa et al. 2013).^{13} We have verified (see also Creminelli et al. 2011) that bispectra can be generated with the same shapes as the “EFT1” and “EFT2” (Senatore et al. 2010; Chen et al. 2007) constrained in the companion paper (Planck Collaboration XVII 2016), which usually arise in models of inflation with nonstandard kinetic terms, with As explained in the previous subsection, the linear combinations of these two bispectra produce both equilateral and orthogonal bispectrum templates. Given Eqs. (104)–(108), we can proceed as in the previous section to exploit the limits on primordial NG in a combined analysis with the power spectra analysis. In Planck Collaboration XVII (2016) the constraint c_{s} ≥ 0.23 (95% CL) is obtained based on the constraints on and . One can proceed as described in Planck Collaboration XVII (2016) to constrain the parameter modifying the consistency relation, Eq. (106). Adopting a loguniform prior on A in the interval 10^{4} ≤ A ≤ 10^{4} and a uniform prior 10^{4} ≤ c_{s} ≤ 1, the Planck measurements on and constrain to be (95% CL) (Planck Collaboration XVII 2016). We also explore the possibility of the negative branch (corresponding to a blue tensor spectral index), finding (95% CL) (Planck Collaboration XVII 2016). By allowing a logarithmic prior on based on the f_{NL} measurements, Fig. 36 shows the joint constraints on ϵ_{1} and η_{s} for the n_{t}< 0 branch and for the n_{t}> 0 branch. Planck TT+lowP+BAO and the NG bounds on constrain ϵ_{1}< 0.036 at 95% CL for n_{t}< 0 (and  ϵ_{1}  < 0.041 for n_{t}> 0).
Fig. 36 Marginalized joint 68% and 95% CL for the Galileon parameters (ϵ_{1},η_{s}) for n_{t}< 0 (left panel) and n_{t}> 0 (right panel). 
10.3. Axion monodromy inflation
10.3.1. Introduction
The mechanism of monodromy inflation (Silverstein & Westphal 2008; McAllister et al. 2010; Kaloper et al. 2011; Flauger et al. 2014b) in string theory motivates a broad class of inflationary potentials of the form (109)Here μ, Λ_{0}, f_{0}, and φ_{0} are constants with the dimension of mass and C_{0}, p, p_{Λ}, p_{f}, and γ_{0} are dimensionless.
In simpler parameterizations used in prior analyses of oscillations from axion monodromy inflation (Peiris et al. 2013; Planck Collaboration XXII 2014; Easther & Flauger 2014; Jackson et al. 2014; Meerburg et al. 2014b,a; Meerburg & Spergel 2014; Meerburg 2014), one assumes p_{Λ} = p_{f} = 0, corresponding to a sinusoidal term with constant amplitude throughout inflation taken to be a periodic function of the canonicallynormalized inflaton φ. Taking p_{Λ} ≠ 0 and p_{f} ≠ 0 allows the magnitude and frequency, respectively, of the modulation to depend on φ. For example, the frequency is always a periodic function of an underlying angular axion field, but its relation to the canonically normalized inflaton field is modeldependent.
The microphysical motivation for p_{Λ} ≠ 0 and p_{f} ≠ 0 is that in string theory additional scalar fields, known as “moduli,” evolve during inflation. The inflationary potential depends on a subset of these fields. Because the magnitude and frequency of modulations are determined by the vacuum expectation values of moduli, both quantities are then naturally functions of φ. The case p_{Λ} = p_{f} = 0 corresponds to when these fields are approximately fixed, stabilized strongly by additional terms in the scalar potential. But in other cases, the axion potential that drives inflation also provides a leading term stabilizing the moduli. The exponential dependence of the magnitude in the potential of Eq. (109) arises because the modulations are generated nonperturbatively, e.g., by instantons. For this reason, the modulations can be undetectably small in this framework, although there are interesting regimes where they could be visible.
Specific examples studied thus far yield exponents p, p_{Λ}, and p_{f} that are rational numbers of modest size. For example, models with p = 3, 2, 4 / 3, 1, and 2 / 3 have been constructed (Silverstein & Westphal 2008; McAllister et al. 2010, 2014), or in another case p = 4 / 3, p_{Λ} = −1 / 3, and p_{f} = −1 / 3. Following Flauger et al. (2014b), we investigate the effect of a drift in frequency arising from p_{f}, neglecting a possible drift in the modulation amplitude by setting p_{Λ} = C_{0} = 0. Even in this restricted model, a parameter exploration using a fully numerical computation of the primordial power spectrum following the methodology of Peiris et al. (2013) is prohibitive, so we follow Flauger et al. (2014b) to study two templates capturing the features of the primordial spectra generated by this potential.
The first template, which we call the “semianalytic” template, is given by (110)The parameter f is higher than the underlying axion decay constant f_{0} of the potential by a few percent, but this difference will be neglected in this analysis. The quantity φ_{0} is some fiducial value for the scalar field, and φ_{k} is the value of the scalar field at the time when the mode with comoving momentum k exits the Hubble radius. At leading order in the slowroll expansion, in units where the reduced Planck mass M_{Pl} = 1, , where , and φ_{end} is the value of the scalar field at the end of inflation.
The second “analytic” template was derived by Flauger et al. (2014b) by expanding the argument of the trigonometric function in Eq. (110) in ln(k/k_{∗}), leading to (111)The relation between the empirical parameters in the templates and the potential parameters are approximated by , where (112)and b is the monotonicity parameter defined in Flauger et al. (2014b), providing relations converting bounds on c_{n} into bounds on the microphysical parameters of the potential. However, the analytic template can describe more general shapes of primordial spectra than just axion monodromy.
As discussed by Flauger et al. (2014b), there is a degeneracy between p (or alternatively n_{s}) and f. For both templates we fix p = 4 / 3 and also fix the tensor power spectrum to its form in the absence of oscillations. This is an excellent approximation because tensor oscillations are suppressed relative to the scalar oscillations by a factor α(f/M_{Pl})^{2} ≪ 1. A uniform prior −π< Δφ<π is adopted for the phase parameter of both templates as well as a prior 0 <δn_{s}< 0.7 for the modulation amplitude parameter.
Fig. 37 Constraints on the parameters of the analytic template, showing joint 68% and 95% CL. The dotted lines correspond to the frequencies showing the highestlikelihood improvements (see text). 
Fig. 38 Constraints on the parameters of the semianalytic template showing joint 68% and 95% CL. The solid lines on the lefthand panel mark the frequencies showing the highestlikelihood improvements (see text). 
In order to specify the semianalytic template, we assume instantaneous reheating, which for p = 4 / 3 corresponds to N_{∗} ≈ 57.5 for k_{∗} = 0.05 Mpc^{1}. We set φ_{0} = 12.38M_{Pl} with φ_{end} = 0.59M_{Pl}. We adopt uniform priors −4 < log _{10}(f/M_{Pl}) < −1 and −0.75 <p_{f}< 1 for the remaining parameters. The priors 0 < ln(α) < 6.9 and −2 <c_{1,2}< 2 specify the analytic template. The singlefield effective field theory becomes strongly coupled for α> 200. However, in principle the string construction remains valid in this regime.
10.3.2. Power spectrum constraints on monodromy inflation
We carry out a Bayesian analysis of axion monodromy inflation using a highresolution version of CAMB coupled to the PolyChord sampler (see Sect. 8.2). For our baseline analysis we conservatively adopt Planck TT+lowP, using the “bin1” highℓTT likelihood. In addition to the primordial template priors specified above, we marginalize over the standard priors for the cosmological parameters, the primordial amplitude, and foreground parameters.
Fig. 39 Frequency residuals for the ln(α) ≈ 3.5 likelihood peak, binned at Δℓ = 30. The ± 1σ errors are given by the square root of the diagonal elements of the covariance matrix. 
The marginalized joint posterior constraints on pairs of primordial parameters for the analytic and semianalytic templates are shown in Figs. 37 and 38, respectively.
The complex structures seen in these plots arise due to degeneracies in the likelihood frequency “beating” between underlying modulations in the data and the model (Easther et al. 2005). Parameter combinations where “beating” occurs over the largest k ranges lead to discrete local maxima in the likelihood. Fortuitous correlations in the observed realization of the C_{ℓ} can give the same effect.
The four frequencies picked out by these structures, ln(α) ≈ {3.5,5.4,6.0,6.8}, show improvements of Δχ^{2} ≈ {−9.7,−7.1,−12.2,−12.5} relative to ΛCDM, respectively. These frequencies are marked by dotted lines in Fig. 37, and by solid lines in Fig. 38 using Eq. (112). The semianalytic and analytic templates lead to selfconsistent results as expected, with analogous structures being picked out by the likelihood in each template. There is no evidence for a drifting frequency, p_{f} ≠ 0 or c_{n} ≠ 0. Thus, these parameters serve to smooth out structures in the marginalized posterior.
The improvement in χ^{2} is not compelling enough to suggest a primordial origin. Fitting a modulated model to simulations with a smooth spectrum can give rise to Δχ^{2} ~ −10 improvements (Flauger et al. 2014b). Furthermore, as the monodromy model contains only a single frequency, at least three of these structures must correspond to spurious fits to the noise. Considering the two models defined by the two templates and the parameter priors specified above, the Bayes factors calculated using PolyChord favours base ΛCDM over both templates by odds of roughly 8:1.
Compared to previous analyses of the linear (p = 1) axion monodromy model for WMAP9 (Peiris et al. 2013) and the 2013 Planck data (Planck Collaboration XXII 2014; Easther & Flauger 2014) the common frequencies are shifted slightly upward. The lower frequency in common appears shifted by a factor of order from α ≈ 28.9 to 31.8 and the higher frequency in common from α ≈ 210 to 223. Flauger et al. (2014b) suggest that the lower frequency (which had Δχ^{2} = −9 in PCI13) was associated with the 4 K cooler line systematic effects in the 2013 Planck likelihood. However, its presence at similar significance in the 2015 likelihood with improved handling of the cooler line systematics suggests that this explanation is not correct. The second frequency, which appeared with Δχ^{2} ≈ −20 in WMAP9 (Peiris et al. 2013) is still present but with much reduced significance, suggesting that the high multipoles do not give evidence for this frequency. Additionally, two higher frequencies are present, which if interpreted as being of primordial origin, correspond to a regime well beyond the validity of the singlefield effective field theory. If one of these frequencies were to be confirmed as primordial, a significantly improved understanding of the underlying string construction would need to be undertaken.
In order to check whether the improvement in fit at these four modulation frequencies is responding to residual foregrounds or other systematics, we examine the frequency residuals. Figure 39 shows the residuals of the data minus the model (including the bestfit foreground model) for the four PLIK frequency combinations binned at Δℓ = 30 for the lowest modulation frequency, ln(α) ≈ 3.5. This plot shows no significant frequency dependence, and thus there is no indication that the fit is responding to frequency dependent systematics. Furthermore, the plot does not show evidence that the improvement for this modulation frequency comes from the feature at ℓ ≈ 800, as suggested by Easther & Flauger (2014). This feature and another at ℓ ≈ 1500 are apparent at all frequency combinations. Similar plots for the three other modulation frequencies also do not show indications of frequency dependence.
In order to confirm whether any of the frequencies picked out here is of primordial origin, one can exploit independent information in the polarization data to perform a crosscheck of the temperature prediction, thus minimizing the “lookelsewhere” effect (Mortonson et al. 2009). Leaving a complete analysis of the independent information in the polarization for future work, we now check whether the temperatureonly result remains stable when highℓ polarization is added in the likelihood. In Fig. 40 we show a preliminary analysis using the PLIK temperature and polarization (TT, TE, and EE) “bin1” likelihood plus lowℓ polarization data. A comparison with the lefthand panels of Figs. 37 and 38 indicates slight differences from the Tonly analysis. However, all the four frequencies identified in the temperature are present when highℓ polarization is added. There is a maximum Δχ^{2} ≈ −8.0 improvement over ΛCDM. We also repeat the analysis using only the EE polarization “bin1” likelihood plus lowℓ temperature and polarization data. These results are presented in Fig. 41. The EEonly frequencies are offset with respect to the temperatureonly frequencies: the bestfit EEonly frequencies are at ln(α) ≈ {3.8,5.0,5.4,5.8,6.2}. The maximum improvement over ΛCDM for this case is Δχ^{2} ≈ −12.5.
Fig. 40 Constraints on the parameters of the analytic (top) and semianalytic (bottom) templates with the addition of highℓ polarization data in the likelihood, showing joint 68% and 95% CL. The lines mark the frequencies showing the highestlikelihood improvements identified in the baseline temperatureonly analysis. 
Fig. 41 Constraints on the parameters of the analytic (top) and semianalytic (bottom) templates with EEonly highℓ polarization data plus lowℓ temperature and polarization data, showing joint 68% and 95% CL. The lines mark the frequencies showing the highestlikelihood improvements identified in the baseline temperatureonly analysis. 
10.3.3. Predictions for resonant nonGaussianity
The lefthand panel of Fig. 42 presents derived constraints on the parameters of the potential in Eq. (109) calculated using the analytic template. Another crosscheck of primordial origin is available since the monodromy model predicts resonant NG, generating a bispectrum whose properties would be strongly correlated with that of the power spectrum (Chen et al. 2008; Flauger & Pajer 2011). Using the mapping (113)we use the analytic template to derive the posterior probability for the resonant NG signal predicted by constraints from the power spectrum, presented in the middle and right panels of Fig. 42.
Planck Collaboration XVII (2016) use an improved modal estimator to scan for resonant NG. The resolution of this scan is currently limited to ln(α) < 3.9, which potentially can probe the lowest frequency picked out in the power spectrum search. However, the modal estimator’s sensitivity (imposed by cosmic variance) of is significantly greater than the predicted value for this frequency from fits to the power spectrum, . Efforts to increase the resolution of the modal estimator are ongoing and may allow consistency tests of the significantly higher levels of resonant NG predicted by the higher frequencies in the future.
Fig. 42 Derived constraints on the parameters of the potential, Eq. (109), as well as the predicted resonant NG, , using the analytic template, showing joint 68% and 95% CL. The dotted lines mark the frequencies showing the highestlikelihood improvements (see text). 
10.3.4. Power spectrum and bispectrum constraints on axion inflation with a gauge field coupling
We now consider the case where the axion field is coupled to a gauge field. Such a scenario is physically well motivated. From an effective field theory point of view the derivative coupling is natural and must be included since it respects the same shift symmetry that leads to axion models of inflation (Anber & Sorbo 2010; Barnaby & Peloso 2011; Pajer & Peloso 2013). This type of coupling is also ubiquitous in string theory (see, e.g., Barnaby et al. 2012; Linde et al. 2013). The coupling term in the action is (Anber & Sorbo 2010; Barnaby & Peloso 2011; Barnaby et al. 2011) (114)where F_{μν} = ∂_{μ}A_{ν}−∂_{ν}A_{μ}, its dual is , and α is a dimensionless constant which, from an effective field theory perspective, is expected to be of order one. For the potential of the axion field, we will not investigate further the consequences of the oscillatory part of the potential, focusing on the coupling of the axion field to the U(1) gauge field (effectively setting Λ_{0} = 0).
The coupling of a pseudoscalar axion with the gauge field has interesting phenomenological consequences, both for density perturbations and primordial gravitational waves (Barnaby & Peloso 2011; Sorbo 2011; Barnaby et al. 2011, 2012; Meerburg & Pajer 2013; Ferreira & Sloth 2014). Gauge field quanta source the axion field via an inverse decay process δA + δA → δϕ, modifying the usual predictions already at the power spectrum level. Additionally, the inverse decay can generate a high level of primordial NG.
The parameter (115)characterizes the strength of theinverse decay effects. If ξ< 1 the coupling is too small to produce any modifications to the usual predictions of the uncoupled model. For previous constraints on ξ see Barnaby et al. (2011, 2012) and Meerburg & Pajer (2013). Using the slowroll approximation and neglecting the small oscillatory part of the potential, one can express (116)where N is, as usual, the number of efolds to the end of inflation. The scalar power spectrum of curvature perturbations is given by (117)where (Meerburg & Pajer 2013) (118)Here an asterisk indicates evaluation at the pivot scale, k_{∗} = 0.05 Mpc^{1}, and and n_{s}−1 = −2ϵ_{1}−ϵ_{2} are the amplitude and spectral index, respectively, of the standard slowroll power spectrum of vacuummode curvature perturbations (the usual power spectrum in the absence of the gaugecoupling). By numerically evaluating the function f_{2}(ξ) (defined in Eq. (3.27) of Barnaby et al. 2011), we created an analytical fit to this function accurate to better than 2% for 0.1 <ξ_{∗}< 7.^{14} In the following, unless stated otherwise, we fix p = 4 / 3 as in the previous subsection and assume instantaneous reheating so that N_{∗} ≈ 57.5 and the slowroll parameters ϵ_{1} and ϵ_{2} are fixed. For the tensor power spectrum we adopt the approximation (Barnaby et al. 2011) (119)where (120)Here and n_{t} = −2ϵ_{1} are the “usual” expressions for the tensor amplitude and tensor tilt in standard slowroll inflation.
The total bispectrum is (Barnaby et al. 2012) where the explicit expression for F_{inv.dec.}(k_{i}) (Barnaby et al. 2011; see also Meerburg & Pajer 2013) is reported in Planck Collaboration XVII (2016). This shows that the inverse decay effects and the resonant effects (which arise from the oscillatory part of the potential) simply “add up” in the bispectrum. The nonlinearity parameter is (122)The function f_{3}(ξ_{∗}) corresponds to the quantity f_{3}(ξ_{∗};1,1) defined in Eq. (3.29) of Barnaby et al. (2011). We have computed f_{3}(ξ_{∗}) numerically and used a fit with an accuracy of better than 2%.^{15} We use the observational constraint (68% CL) obtained in Planck Collaboration XVII (2016) from an analysis where only the inverse decay type NG is assumed present. We omit the explicit expression for the resonant bispectrum B_{res}, since it will not be used here.
We carried out an MCMC analysis of constraints on the (scalar and tensor) power spectra predicted by this model with the Planck TT+lowP likelihood, marginalizing over standard priors for the cosmological parameters and foreground parameters with the uniform priors and 0.1 ≤ ξ_{∗} ≤ 7.0.
The power spectrum constraint gives (123)Given that is exponentially sensitive to ξ, this translates into the prediction (using Eq. (122)) , which is significantly tighter than the current bispectrum constraint from Planck Collaboration XVII (2016). Indeed, importance sampling with the likelihood for , taken to be a Gaussian centred on the NG estimate (68% CL) (Planck Collaboration XVII 2016), changes the limit on ξ_{∗} only at the second decimal place.
We now derive constraints on model parameters using only the observational constraint on The constraints thus derived are applicable for generic p and also to the axion monodromy model discussed in Sect. 10.3, even in the case Λ_{0} ≠ 0. We follow the procedure described in Sect. 11 of Planck Collaboration XVII (2016). The likelihood for is taken to be a Gaussian centred on the NG estimate (68% CL) (Planck Collaboration XVII 2016). We use the expression of Eq. (122), where f_{3}(ξ_{∗}) is numerically evaluated. To find the posterior distribution for the parameter ξ_{∗} we choose uniform priors in the intervals and 0.1 ≤ ξ_{∗} ≤ 7.0. This yields 95% CL constraints for ξ_{∗} (for any value of p) of (124)If we choose a logconstant prior on ξ_{∗} we find (125)For both cases the results are insensitive to the upper limit chosen for the prior on ξ_{∗} since the likelihood quickly goes to zero for ξ_{∗}> 3. As the likelihood for ξ_{∗} is fairly flat, the tighter constraint seen for the logconstant case is mildly prior driven. The constraints from the bispectrum are consistent with, and slightly worse than, the result from the power spectrum alone.
Using a similar procedure and Eq. (116) one can also obtain a constraint on α/f. Adopting a logconstant prior^{16}2 ≤ α/f ≤ 100 and uniform priors 50 ≤ N_{∗} ≤ 70 and we obtain the 95% CL constraints (126)and (127)For example, for a linear potential, p = 1, if α ~ 1 as suggested by effective field theory, then the axion decay constant f is constrained to be (128)while for a potential with p = 4 / 3 we find (129)These limits are complementary to those derived in Sect. 10.3.
11. Constraints on isocurvature modes
Fig. 43 Angular power spectra for the scaleinvariant (i.e., n_{ℛℛ} = 1) pure adiabatic mode (ADI, green dashed curves) and for the scale invariant (n_{ℐℐ} = 1) pure isocurvature (CDI, NDI, or NVI) modes, with equal primordial perturbation amplitudes. The thick lines represent the temperature autocorrelation (TT) and the thin lines the Emode polarization autocorrelation (EE). 
In PCI13, we presented constraints on a number of simple models featuring a mixture of the adiabatic (ADI) mode and one type of isocurvature mode. We covered the cases of CDM density isocurvature (CDI), neutrino density isocurvature (NDI), and neutrino velocity isocurvature (NVI) modes (Bucher et al. 2000) with different assumptions concerning the correlation (Langlois 1999; Amendola et al. 2002) between the primordial adiabatic and isocurvature perturbations. Isocurvature modes, possibly correlated among themselves and with the adiabatic mode, can be generated in multifield models of inflation; however, at present a mechanism for exciting the neutrino velocity isocurvature mode is lacking. Section 11.2 shows how these constraints have evolved with the new Planck TT+lowP likelihoods, how much including the Planck lensing likelihood changes the results, and what extra information the Planck highℓ polarization contributes. A pure isocurvature mode as a sole source of perturbations has been ruled out (Enqvist et al. 2002), since, as can be seen from Fig. 43, any of the isocurvature modes leads to an acoustic peak structure for the temperature angular power very different from the adiabatic mode, which fits the data very well. The different phases and tilts of the various modes also occur in the polarization spectra, as shown in Fig. 43 for the E mode.^{17}
In Sect. 11.4 we add one extra degree of freedom to the generallycorrelated ADI+CDI model by allowing primordial tensor perturbations (assuming the inflationary consistency relation for the tilt of the tensor power spectrum and its running). Our main goal is to explore a possible degeneracy between tensor modes and negativelycorrelated CDI modes, tending to tilt the largescale temperature spectrum in opposite directions. In Sect. 11.5, we update the constraints on three special cases motivated by axion or curvaton scenarios.
The goal of this analysis is to test the hypothesis of adiabaticity and establish the robustness of the base ΛCDM model against different assumptions concerning initial conditions (Sect. 11.3). Adiabaticity is also an important probe of the inflationary paradigm, since any significant detection of isocurvature modes would exclude the possibility that all perturbations in the Universe emerged from quantum fluctuations of a single inflaton field, which can excite only one degree of freedom, the curvature (i.e., adiabatic) perturbation.^{18}
In this section, theoretical predictions were obtained with a modified version of the CAMB code (version Jul14) while parameter exploration was performed with the MultiNest nested sampling algorithm.
11.1. Parameterization and notation
A general mixture of the adiabatic mode and one isocurvature mode is described by the three functions and describing the curvature, isocurvature, and crosscorrelation power spectra, respectively. Our sign conventions are such that positive values for correspond to a positive contribution of the crosscorrelation term to the SachsWolfe component of the total temperature spectrum.
As in PCI13, we specify the amplitudes at two scales k_{1}<k_{2} and assume powerlaw behaviour, so that (130)where a,b = ℐ,ℛ and ℐ = ℐ_{CDI}, ℐ_{NDI}, or ℐ_{NVI}. We set k_{1} = 0.002 Mpc^{1} and k_{2} = 0.100 Mpc^{1}, so that [ k_{1},k_{2} ] spans most of the range constrained by the Planck data. The positive definiteness of the initial condition matrix imposes a constraint on its elements at any value of k: (131)We take uniform priors on the positive amplitudes, The correlation spectrum can be positive or negative. For a ≠ b we apply a uniform prior at large scales (at k_{1}): (134)but reject all parameter combinations violating the constraint in Eq. (131). To ensure that Eq. (131) holds for all k, we restrict ourselves to a scaleindependent correlation fraction: (135)Thus is a derived parameter^{19} given by (136)which in terms of spectral indices is equivalent to (137)The conservative baseline likelihood is Planck TT+lowP. The results obtained with Planck TT, TE, EE+lowP should be interpreted with caution because the data used in the 2015 release are known to contain some low level systematics, in particular arising from TtoE leakage, and it is possible that such systematics may be fit by the isocurvature autocorrelation and crosscorrelation templates. (See Planck Collaboration XIII 2016 for a detailed discussion.)
In what follows, we quote our results in terms of derived parameters identical to those in PCI13. We define the primordial isocurvature fraction as (138)Unlike the primordial correlation fraction cosΔ defined in Eq. (135), β_{iso} is scaledependent in the general case. We present bounds on this quantity at k_{low} = k_{1}, k_{mid} = 0.050 Mpc^{1}, and k_{high} = k_{2}.
We report constraints on the relative adiabatic (ab = ℛℛ), isocurvature (ab = ℐℐ), and correlation (ab = ℛℐ) components according to their contribution to the observed CMB temperature variance in various multipole ranges: (139)where (140)The ranges considered are (ℓ_{min},ℓ_{max}) = (2,20), (21,200),(201,2500), and (2,2500), where the last range describes the total contribution to the observed CMB temperature variance. Here α_{ℛℛ} measures the adiabaticity of the temperature angular power spectrum, a value of unity meaning “fully adiabatic initial conditions”. Values less than unity mean that some of the observed power comes from the isocurvature or correlation spectrum, while values larger than unity mean that some of the power is “cancelled” by a negativelycorrelated isocurvature contribution. The relative nonadiabatic contribution can be expressed as α_{non  adi} ≡ 1−α_{ℛℛ} = α_{ℐℐ} + α_{ℛℐ}.
Fig. 44 68% and 95% CL constraints on the primordial perturbation power in general mixed ADI+CDI a); ADI+NDI b); and ADI+NVI c) models at two scales, k_{1} = 0.002 Mpc^{1} (1) and k_{2} = 0.100 Mpc^{1} (2), for Planck TT+lowP (grey regions highlighted by dotted contours), Planck TT+lowP+lensing (blue), and Planck TT, TE, EE+lowP (red). In the first panels, we also show contours for the pure adiabatic base ΛCDM model with the corresponding colours of solid lines. 
Fig. 45 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 adiabatic spectral index, n_{ℛℛ}, the isocurvature spectral index, n_{ℐℐ}, and the correlation spectral index, n_{ℛℐ} = (n_{ℛℛ} + n_{ℐℐ}) / 2, with Planck TT+lowP data (dashed curves) and TT, TE, EE+lowP data (solid curves), for the generallycorrelated mixed ADI+CDI (black), ADI+NDI (red), and ADI+NVI (blue) models. All these parameters are derived, and the distributions shown here result from a uniform prior on the primary parameters, as detailed in Eqs. (132)–(134). However, the effect of the nonflat derivedparameter priors is negligible for all parameters except for n_{ℐℐ} (and n_{ℛℐ}) where the prior biases the distribution toward one. With TT+lowP, the flatness of β_{iso}(k_{high}) in the CDI case up to a “threshold” value of about 0.5 is a consequence of the (k/k_{eq})^{2} damping of its transfer function as explained in Footnote 17. 
11.2. Results for generallycorrelated adiabatic and one isocurvature mode (CDI, NDI, or NVI)
Results are reported as 2D and 1D marginalized posterior probability distributions. Numerical 95% CL intervals or upper bounds are tabulated in Table 16.
Figure 44 shows the Planck 68% and 95% CL contours for various 2D combinations of the primordial adiabatic and isocurvature amplitude parameters at large scales (k_{1} = 0.002 Mpc^{1}) and small scales (k_{2} = 0.100 Mpc^{1}) for (a) the generallycorrelated ADI+CDI; (b) ADI+NDI; and (c) ADI+NVI models. Overall, the results using Planck TT+lowP are consistent with the nominal mission results in PCI13, but slightly tighter. In the first panels of Figs. 44a−c we also show the constraints on the curvature perturbation power in the pure adiabatic case. Comparing the generallycorrelated isocurvature case to the pure adiabatic case with the same data combination summarizes neatly what the data tell us about the initial conditions. If the contours in the  plane were shifted significantly relative to the pure adiabatic case, the missing power could come either from the isocurvature and postive correlation contributions, or the extra adiabatic power could be cancelled by a negative correlation contribution. We can see that these shifts are small. The lowℓ temperature data continue to mildly favour a negative correlation (see in particular the bottom middle panel for each of the three models), since compared to the prediction of the bestfit adiabatic base ΛCDM model, the TT angular power at multipoles ℓ ≲ 40 is somewhat low. But the dotted grey shaded contours in the three middle top panels show that for Planck TT+lowP, the posterior peaks at values (,) entirely consistent with (0, 0), i.e., the pure adiabatic case is preferred. The bestfit values of (,) are (1.4 × 10^{11}, 4.7 × 10^{13}) for CDI, (1.2 × 10^{12}, 4.6 × 10^{10}) for NDI, and (1.6 × 10^{12}, 2.3 × 10^{10}) for NVI, while (,) ≈ (2.4 × 10^{9}, 2.1 × 10^{9}). It may appear from the bottomcentre panels of Fig. 44 that there is nonzero posterior probability for when , which would violate the positivity constraint, Eq. (131). However, the leftmost pixels of the plots are actually evaluated at values of large enough that the constraint is satisfied.
Including the Planck lensing likelihood does not significantly affect the nonadiabatic primordial powers, except for tightening the constraints on the adiabatic power (see the blue versus black contours in the first panels of Figs. 44a−c). Including the lensing () likelihood constrains the optical depth τ more tightly than the highℓ temperature and lowℓ polarization alone (Planck Collaboration XIII 2016). As there is a strong degeneracy between τ and the primordial (adiabatic) perturbation power (denoted in the other sections of this paper by A_{s}), it is natural that adding the lensing data leads to stronger constraints on . Moreover, replacing the lowℓ likelihood Planck lowP by Planck lowP+WP constrains τ better (Planck Collaboration XIII 2016). In the ADI+CDI case the effect of this replacement was very similar to adding the Planck lensing data (see also Table 16). Although the Planck lensing data do not directly constrain the isocurvature contribution,^{20} they can shift and tighten the constraints on some derived isocurvature parameters by affecting the favoured values of the standard parameters (present even in the pure adiabatic model). However this effect is small as confirmed in Table 16. Therefore, in the figures we do not show 1D posteriors of the derived isocurvature parameters for Planck TT+lowP+lensing, since they would be (almost) indistinguishable from Planck TT+lowP, as we see in Fig. 44 for the primary nonadiabatic parameters.
In contrast, the highℓpolarization data significantly tighten the bounds on isocurvature and crosscorrelation parameters, as seen by comparing the dotted grey and red contours in Fig. 44. The significant negative correlation previously allowed by the temperature data in the ADI+CDI and ADI+NDI models is now disfavoured. This is also clearly visible in the 1D posteriors of primordial and observable isocurvature and crosscorrelation fractions shown, respectively, in Figs. 45 and 46. Note how the cosΔ and α_{ℛℐ} parameters are driven towards zero by the inclusion of the highℓ TE, EE data (from the dashed to the solid lines) in the ADI+CDI and ADI+NDI cases. We also observed that when the lowP data are replaced by a simple Gaussian prior on the reionization optical depth (τ = 0.078 ± 0.019), the trend is similar. The highℓ (ℓ ≥ 30) Planck TT data allow a large negative correlation, while the highℓPlanck TE, EE data prefer positive correlation. This is clearly seen in Fig. 47 for the ADI+CDI case. The bestfit values show an even more dramatic effect. We find cosΔ = −0.55 with TT+lowP, and + 0.15 with TT, TE, EE+lowP.
Hence there is a competition between the temperature and polarization data that balances out and yields almost symmetric results about zero correlation (except in the ADI+NVI case). The isocurvature autocorrelation amplitude is also strongly reduced, especially in the ADI+CDI case. The bestfit values are slightly offset from , but the pure adiabatic model still lies inside the 68% CL (for ADI+CDI and ADI+NDI) or 95% CL (for ADI+NVI) regions. In summary, the highℓ polarization data exhibit a strong preference for adiabaticity, although one should keep in mind the possibility of unaccounted systematic effects in the polarization data, possibly leading to artificially strong constraints. For example, the tendency for polarization to shift the constraints towards positive correlation may be due to particular systematic effects that mimic modified acoustic peak structure, as we discussed in Sect. 11.1.
Fig. 46 Constraints on the fractional contribution of the adiabatic (ℛℛ), isocurvature (ℐℐ), and correlation (ℛℐ) components to the CMB temperature variance in various multipole ranges, as defined in Eq. (139), with Planck TT+lowP data (dashed curves) and with Planck TT, TE, EE+lowP data (solid curves). These are shown for the generallycorrelated mixed ADI+CDI (black), ADI+NDI (red), or ADI+NVI (blue) models. 
We also performed a parameter extraction with the Planck TT, TE, EE+lowP+lensing data, but this combination did not provide interesting new constraints. We found only a tightening of bounds on the standard adiabatic parameters as in the Planck TT+lowP+lensing case.
Fig. 47 Constraints on the primordial correlation fraction, cosΔ, in the mixed ADI+CDI model with Planck TT+lowP data (dashed black curve) compared to the case where Planck lowP data are not used, but replaced by a Gaussian prior τ = 0.078 ± 0.019 (dashed red curve). The same exercise is repeated with Planck TT, TE, EE data (solid curves) demonstrating that to a great extent the preferred value of cosΔ is driven by the highℓ data. 
We provide 95% CL upper limits or ranges for β_{iso}, cosΔ, and α_{ℛℛ} in Table 16. With Planck TT+lowP, the constraints on the nonadiabatic contribution to the temperature variance, 1−α_{ℛℛ}(2,2500), are (−1.5%, 1.9%), (−4.0%, 1.4%), and (−2.3%, 2.4%) in the ADI+CDI, ADI+NDI, and ADI+NVI cases, respectively.^{21} With Planck TT, TE, EE+lowP these tighten to (0.1%, 1.5%), (−0.1%, 2.2%), and (−2.0%, 0.8%). In the ADI+CDI case, zero is not in the 95% CL interval, but this should not be considered a detection of nonadiabaticity. For example, as mentioned above, is in the 68% CL region, and the bestfit values are (,) = (1.0 × 10^{13}, 3.5 × 10^{9}). Moreover, the improvement in χ^{2} with respect to the adiabatic model is only 5.3 with 3 extra parameters, so this is not a significant improvement of fit. Indeed, for all generallycorrelated mixed models the improvement in χ^{2} is very small. In particular, with Planck TT+lowP it does not even exceed the number of extra degrees of freedom, which is three (see Table 16).
Finally, we checked whether there is any Bayesian evidence for the presence of generallycorrelated adiabatic and isocurvature modes. In all cases and with all data combinations studied, the Bayesian model comparison supports the null hypothesis, i.e., adiabaticity. Indeed, the logarithm of the evidence ratio is lnB = ln(P_{ISO}/P_{ADI}) < −5 (i.e., odds of much greater than 150:1 in favour of pure adiabaticity within Planck’s accuracy and given the parameterization and prior ranges used in our analysis), except for ADI+NDI with Planck TT+lowP+lensing, for which the evidence ratio is slightly larger, −4.6, corresponding to odds of 1:100 for the ADI+NDI model compared to the pure adiabatic model.
Fig. 48 Constraints on selected “standard” cosmological parameters with Planck TT+lowP data for the generallycorrelated ADI+CDI (black), ADI+NDI (red), and ADI+NVI (blue) models compared to the pure adiabatic case (ADI, green dashed curves). 
11.3. Robustness of the determination of standard cosmological parameters
Another outcome of our analysis is the robustness of the determination of the standard cosmological parameters against assumptions on initial conditions. Figure 48 shows the 1D marginalized posteriors for several cosmological parameters (not all independent of each other) with the Planck TT+lowP data alone. For the first time, we observe that in the presence of one generallycorrelated isocurvature mode (CDI, NDI, or NVI), predictions for these parameters remain very stable with respect to the pure adiabatic case. Except for the ADI+NDI case, the posteriors neither broaden nor shift significantly. A small broadening is only observed in the sound horizon angle θ_{MC}, which is naturally the most sensitive parameter to tiny disturbances of the acoustic peak structure. In the ADI+NDI case, the peak of the posterior distribution for some parameters shifts slightly, but the largest shift (for Ω_{c}h^{2}) is less than 1σ.
It is striking that a scaleinvariant adiabatic spectrum (n_{ℛℛ} = 1) is excluded at many σ even when isocurvature modes are allowed: at 4.7σ (ADI+CDI), 5.0σ (ADI+NDI), and 5.4σ (ADI+NVI). This illustrates how much the constraining power of the CMB has improved. With WMAP data, there was still a strong degeneracy between, for example, the primordial isocurvature fraction and the adiabatic spectral index (Valiviita & Giannantonio 2009; Savelainen et al. 2013). This degeneracy nearly disappears with Planck TT+lowP, and even more so with Planck TT, TE, EE+lowP, as shown in the upper panel of Fig. 49. Contours in the (n_{ℛℛ}, cosΔ) space also shrink considerably, with some correlation remaining between these parameters in the ADI+CDI and ADI+NVI cases (Fig. 49, lower panel).
11.4. CDI and primordial tensor perturbations
A primordial tensor contribution adds extra temperature angular power at low multipoles, where the adiabatic base ΛCDM model predicts slightly more power than seen in the data. Hence allowing for a nonzero tensortoscalar ratio r might tighten the constraints on positivelycorrelated isocurvature, but degrade them in negativelycorrelated models. We test how treating r as a free parameter affects the constraints on isocurvature and how allowing for the generallycorrelated CDI mode affects the constraints on r. These cases are denoted as “CDI+r”. For comparison, we examine the pure adiabatic case in the same parameterization, and call it “ADI+r”. We also consider another approach where we fix r = 0.1. These cases are named “CDI+r = 0.1” and “ADI+r = 0.1”.
Fig. 49 Dependence of the determination of the adiabatic spectral index n_{ℛℛ} (called n_{s} in the other sections of this paper) on the primordial isocurvature fraction β_{iso} and correlation fraction cosΔ, with Planck TT+lowP data (dashed contours) and with Planck TT, TE, EE+lowP data (shaded regions). 
In the pure adiabatic case (where the curvature and tensor perturbations stay constant on superHubble scales), the primordial r is the same as the tensortoscalar ratio at the Hubble radius exit of perturbations during inflation, which we call . However, in the presence of an isocurvature component, is not constant in time even on superHubble scales (GarcíaBellido & Wands 1996). Instead, the isocurvature component may source , for example if the background trajectory in the field space is curved between Hubble exit and the end of inflation (Langlois 1999; Langlois & Riazuelo 2000; Gordon et al. 2001; Amendola et al. 2002). As a result, we will have at the primordial time , where is the curvature power at Hubble exit. That is, by the primordial time the curvature perturbation power is larger than at the Hubble radius exit time (Bartolo et al. 2001; Wands et al. 2002; Byrnes & Wands 2006). Thus we find a relation (Savelainen et al. 2013; Valiviita et al. 2012; Kawasaki & Sekiguchi 2008): (141)i.e., the tensortoscalar ratio at the primordial time (r) is smaller than the ratio at the Hubble radius exit time ().
The derivation of Eq. (141) assumes that the adiabatic and isocurvature perturbations are uncorrelated at Hubble radius exit (), and that all the possible primordial correlation (cosΔ ≠ 0) appears from the evolution of superHubble perturbations between Hubble exit and the primordial time. This is true to leading order in the slowroll parameters, but inflationary models that break slow roll might produce perturbations that are strongly correlated already at the Hubble radius exit time. In these cases the correlation would depend on the details of the particular model, such as the detailed shape of the potential and the interactions of the fields. However, a generic prediction of slowroll inflation is that, at Hubble radius exit, the crosscorrelation is very weak, and indeed is of the order of the slowroll parameters compared to the autocorrelations and (see, e.g., Byrnes & Wands 2006). Thus, for slowroll models, .
In our analysis, we fix the tensor spectral index by the leadingorder inflationary consistency relation, which now reads (Wands et al. 2002) (142)Assuming a uniform prior for r would lead to huge negative n_{t} whenever cos^{2}Δ was close to one. Therefore, when studying the CDI+r case we assume a uniform prior on at k = 0.05 Mpc^{1} (for details, see Savelainen et al. 2013).
Surprisingly, allowing for a generallycorrelated CDI mode (i.e., three extra parameters) hardly changes the constraints on r from those obtained in the pure adiabatic model. In Fig. 50 we demonstrate this in a “standard” plot of r_{0.002} versus adiabatic spectral index.
From Table 16 we notice that, with Planck TT+lowP and TT, TE, EE+lowP, fixing r to 0.1 tightens constraints on the primordial isocurvature fraction at large scales. This is as we expected, since both tensor and isocurvature perturbations add power at low ℓ, and the data do not prefer this. However, the shapes of the tensor spectrum and correlation spectrum are such that negative correlation cannot efficiently cancel the unwanted extra power over all scales produced by tensor perturbations (at ℓ ≲ 70). Therefore, the correlation fraction cosΔ is almost unaffected. However, when we allow r to vary, the cancelation mechanism works to some degree when using Planck TT+lowP data, leading to more negative cosΔ than without r: with varying r we have cosΔ in the range (−0.43, 0.20), while without r it is in (−0.30, 0.20), at 95% CL. As there is now some cancellation of power at large scales, the constraint on β_{iso}(k_{low}) weakens slightly from 0.041 without r to 0.043 with r. On the other hand, the highℓ polarization data constrain the correlation to be so close to zero that with Planck TT, TE, EE+lowP the results for cosΔ with and without r are almost identical.
The mean value of cosΔ in the CDI+r cases is −0.071 (TT+lowP) and −0.076 (TT, TE, EE+lowP). Therefore, 1−cos^{2}Δ ≈ 0.99, and so we do not expect a large difference between the primordial r and the Hubble radius exit value . The smallness of the difference is evident in Table 15. To summarize, CDI hardly affects the determination of r from the Planck data, and allowing for tensor perturbations hardly affects the determination of the nonadiabaticity parameters.
Fig. 50 68% and 95% CL constraints on the primordial adiabatic spectral index n_{ℛℛ} and the primordial tensortoscalar ratio r (more accurately, in the CDI+r model, the primordial tensortocurvature power ratio) at k = 0.002 Mpc^{1}. Filled contours are for generallycorrelated ADI+CDI and solid contours for the pure adiabatic model. 
95% CL upper bounds on the tensortoscalar ratio (actually the tensortocurvature power ratio) at the primordial time, r, and earlier, at the Hubble radius exit time during inflation, , at k = 0.05 Mpc^{1}.
Constraints on mixed adiabatic and isocurvature models.
11.5. Special CDI cases
Fig. 51 Uncorrelated ADI+CDI with n_{ℐℐ} = 1 (“axion”). 
Next we study three oneparameter CDI extensions to the adiabatic model. In all these extensions the isocurvature mode modifies only the largest angular scales, since we either fix n_{ℐℐ} to unity (“axion”) or to the adiabatic spectral index (“curvaton I/II”). As can be seen from Fig. 43, the polarization E mode at multipoles ℓ ≳ 200 will not be significantly affected by this type of CDI mode. Therefore, these models are much less sensitive to residual systematic effects in the highℓ polarization data than the generallycorrelated models.
11.5.1. Uncorrelated ADI+CDI (“Axion”)
We start with an uncorrelated mixture of adiabatic and CDI modes () and make the additional assumption that , i.e., we assume unit isocurvature spectral index, n_{ℐℐ} = 1. Constraints in the (n_{ℛℛ},β_{iso}) plane are presented in Fig. 51. This model is the only case for which our new results do not improve over bounds from PCI13. At k_{mid} = 0.050 Mpc^{1}, we find β_{iso}< 0.038 (95% CL, TT, TE, EE+lowP; see Table 16), compared with β_{iso}< 0.039 using Planck 2013 and lowℓ WMAP data. This is not surprising, since fixing n_{ℐℐ} to unity implies that bounds are dominated by measurements on very large angular scales, ℓ ≲ 30, as can easily be understood from Fig. 43. Hence the results are insensitive to the addition of better highℓ temperature data, or new highℓ polarization data.
We summarized in PCI13 why an uncorrelated CDI mode with n_{ℐℐ} ≈ 1 can be produced in axion models under a number of restrictive assumptions: the PecceiQuinn symmetry should be broken before inflation; it should not be restored by quantum fluctuations of the inflaton or by thermal fluctuations when the Universe reheats; and axions produced through the misalignment angle should contribute to a sizable fraction (or all) of the dark matter. Under all of these assumptions, limits on β_{iso} can be used to infer a bound on the energy scale of inflation, using Eq. (73) of PCI13. This bound is strongest when all the dark matter is assumed to be in the form of axions. In that case, the limit on β_{iso}(k_{mid}) for Planck TT, TE, EE+lowP gives (143)where H_{inf} is the expansion rate at Hubble radius exit of the scale corresponding to k_{mid} = 0.050 Mpc^{1} and f_{a} is the PecceiQuinn symmetrybreaking energy scale.
11.5.2. Fully correlated ADI+CDI (“Curvaton I”)
Another interesting special case of mixed adiabatic and CDI (or BDI) perturbations is a model where these perturbations are primordially fully correlated and their power spectra have the same shape. These cases are obtained by setting , which, by condition (136), implies that the corresponding statement holds at scale k_{2} and indeed at any scale. In addition, we set , i.e., n_{ℐℐ} = n_{ℛℛ}. From this it follows that β_{iso} is scaleindependent. Therefore, this model has only one primary nonadiabaticity parameter, .
A physically motivated example of this type of model is the curvaton model (Mollerach 1990; Linde & Mukhanov 1997; Enqvist & Sloth 2002; Moroi & Takahashi 2001; Lyth & Wands 2002; Lyth et al. 2003) with the following assumptions. (1) The average curvaton field value is sufficiently below the Planck mass when cosmologically interesting scales exit the Hubble radius during inflation. (2) At Hubble radius exit, the curvature perturbation from the inflaton is negligible compared to the perturbation caused by the curvaton. (3) The same is true for any inflaton decay products after reheating. This means that, after reheating, the Universe is homogeneous, except for the spatially varying entropy (i.e., isocurvature perturbation) due to the curvaton field perturbations. (4) Later, CDM is created from the curvaton decay and baryon number after curvaton decay. This corresponds to case 4 presented in Gordon & Lewis (2003). (5) The curvaton contributes a significant amount to the energy density of the Universe at the time of the curvaton’s decay to CDM, i.e., the curvaton decays late enough. (6) The energy density of curvaton particles possibly produced during reheating should be sufficiently low (Bartolo & Liddle 2002; Linde & Mukhanov 2006). (7) The smallscale variance of curvaton perturbations, , is negligible, so that it does not significantly contribute to the average energy density on CMB scales; see Eq. (102) in Sasaki et al. (2006). The last two conditions are necessary in order to have an almostGaussian curvature perturbation, as required by the Planck observations. Indeed, if they are not valid, a large follows, as discussed below. The conditions (6) and (7) are related, since curvaton particles would add a homogeneous component to the average energy density on large scales, and hence we can describe their effect by , where is the average energy density of the classical curvaton field on large scales; see Eq. (98) in Sasaki et al. (2006). Then the total energy density carried by the curvaton will be .
The amount of isocurvature and nonGaussianity present after curvaton decay depends on the “curvaton decay fraction” (144)evaluated at curvaton decay time. If conditions (6) and (7) do not hold, then the isocurvature perturbation disappears.^{22}
Fig. 52 Fully correlated ADI+CDI with n_{ℐℐ} = n_{ℛℛ} (“curvaton I”). Since the spectral indices are equal, the primordial isocurvature fraction β_{iso} is scaleindependent. 
The curvaton scenario presented here is one of the simplest to test against observations. It should be noted that at least the conditions (1)−(5) listed at the beginning of this subsection should be satisfied simultaneously. Indeed, if we relax some of these conditions, almost any type of correlation can be produced. For example, the relative correlation fraction can be written as , where . Therefore, the model is fully correlated only if λ ≫ 1. If the slowroll 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.
As seen in Fig. 52 and Table 16, the upper bound on the primordial isocurvature fraction in the fullycorrelated ADI+CDI model weakens slightly when we add the Planck lensing data to Planck TT+lowP, whereas adding highℓ TE, EE tightens the upper bound moderately. With all of these three data combinations, the pure adiabatic model gives an equally good bestfit χ^{2} as the fullycorrelated ADI+CDI model. Bayesian model comparison strengthens the conclusion that the data disfavour this model with respect to the pure adiabatic model.
The isocurvature fraction is connected to the curvaton decay fraction in Eq. (144) by (145)(see case 4 in Gordon & Lewis 2003). We can convert the constraints on β_{iso} from Table 16 into constraints on r_{D} and further into the nonGaussianity parameter assuming a quadratic potential for the curvaton and instantaneous decay^{23}(Sasaki et al. 2006): (146)If conditions (6) and (7) hold, i.e., , as implicitly assumed, e.g., in Bartolo et al. (2004a,b), then the smallest possible value of is −5 / 4, which is obtained when r_{D} = 1, and Eqs. (145) and (146) yield for the various Planck data sets (at 95% CL):^{24}Thus the results for the simplest curvaton model remain unchanged from those presented in PCI13. In in order to produce almost purely adiabatic perturbations, the curvaton should decay when it dominates the energy density of the Universe (r_{D}> 0.98), and the nonGaussianity parameter is constrained to close to its smallest possible value (), which is consistent with the result (68% CL, from T only) found in Planck Collaboration XVII (2016).
Fig. 53 Fully anticorrelated ADI+CDI with n_{ℐℐ} = n_{ℛℛ} (“curvaton II”). 
11.5.3. Fully anticorrelated ADI+CDI (“Curvaton II”)
The curvaton scenario or some other mechanism could also produce 100% anticorrelated perturbations, with n_{ℐℐ} = n_{ℛℛ}. The constraints in the (n_{ℛℛ},β_{iso}) plane are presented in Fig. 53. Examples of this kind of model are provided by cases 2, 3, and 6 in Gordon & Lewis (2003). These lead to a fixed, large amount of isocurvature, e.g., in case 2 to β_{iso} = 9 / 10, and are hence excluded by the data at very high significance. However, case 9 in Gordon & Lewis (2003), with a suitable r_{D} (i.e., r_{D}>R_{c}, where R_{c} = ρ_{c}/ (ρ_{c} + ρ_{b})), leads to fullyanticorrelated perturbations and might provide a good fit to the data. In this case CDM is produced by curvaton decay while baryons are created earlier from inflaton decay products and do not carry a curvature perturbation. We obtain a very similar expression to Eq. (145), namely (150)We convert this to an approximate constraint on r_{D} by fixing R_{c} to its bestfit value, R_{c} = 0.8437 (Planck TT+lowP), within this model. The results for the various Planck data sets are: After all the tests conducted in this section, both for the generallycorrelated CDI, NDI, and NVI cases as well as for the special CDI cases, we conclude that within the spatially flat base ΛCDM model, the initial conditions of perturbations are consistent with the hypothesis of pure adiabaticity, a conclusion that is also supported by the Bayesian model comparison. Moreover, Planck Collaboration XVII (2016) reports a null detection of isocurvature nonGaussianity, with polarization improving constraints significantly.
12. Statistical anisotropy and inflation
A key prediction of standard inflation, which in the present context includes all single field models of inflation as well as many multifield models, is that the stochastic process generating the primordial cosmological perturbations is completely characterized by its power spectrum, constrained by statistical isotropy to depend only on the multipole number ℓ. This statement applies at least to the accuracy that can be probed using the CMB given the limitations imposed by cosmic variance, since all models exhibit some level of nonGaussianity. Nevertheless, more general Gaussian stochastic processes can be envisaged for which one or more special directions on the sky are singled out, so that the expectation values for the temperature multipoles take the form (154)rather than the very special form (155)which is the only possibility consistent with statistical isotropy.
The most general form for a Gaussian stochastic process on the sphere violating the hypothesis of statistical isotropy in Eq. (154) is too broad to be useful, given that we have only one sky to analyse. For ℓ<ℓ_{max}, there are multipole expansion coefficients, compared with model parameters. Therefore, in order to make some progress on testing the hypothesis of statistical isotropy, we must restrict ourselves to examining only the simplest models violating statistical isotropy, for which the available data can establish meaningful constraints and for which one can hope to find a simple theoretical motivation.
12.1. Asymmetry: observations versus model building
In one simple class of statistically anisotropic models, we start with a map produced by a process respecting statistical isotropy, which becomes modulated by another field in the following manner to produce the observed sky map: (156)where denotes a position on the celestial sphere and is the outcome of the underlying statistically isotropic process before modulation. Roughly speaking, where the modulating field is positive, power on scales smaller than the scale of variation of is enhanced, whereas where is negative, power is suppressed. We refer to this as a “power asymmetry.” If , we have a model of dipolar modulation with amplitude A and direction , but higherorder or mixed modulation may also be considered, such as a quadrupole modulation or modulation by a scaleinvariant field , to name just a few special cases. Alternatively, and more closely tied to physical models, we can consider modulations of the position or kspace fluctuations.
In Planck Collaboration XXIII (2014) and Planck Collaboration XVI (2016), the details of constructing efficient estimators for statistical anisotropy, in particular in the presence of realistic data involving sky cuts and possibly incompletely removed foreground contamination, are considered in depth. In addition, the question of the statistical significance of any detected “anomalies” from the expectations of base ΛCDM is examined in detail. Importantly, in the absence of a particular inflationary model for such an observed anomaly, the significance should be corrected for the “multiplicity of tests” that could have resulted in similarlysignificant detections (i.e., for the “look elsewhere effect”), although applying such corrections can be ambiguous. In this paper, however, we consider only forms of statistical anisotropy that are predicted by specific inflationary models, and hence such corrections will not be necessary.
Several important questions can be posed regarding the link between statistical isotropy and inflation. In particular, we can ask the following questions. (1) Does a statistically significant finding of a violation of statistical isotropy falsify inflation? (2) If not, what sort of nonstandard inflation could produce the required departure from statistical isotropy? (3) What other perhaps noninflationary models could also account for the violation of statistical isotropy? In this section, we begin to address these questions by assessing the viability of an inflationary model for dipolar asymmetry, as well as by placing new limits on the presence of quadrupolar power asymmetry.
For the case of the observed dipolar asymmetry examined in detail in Planck Collaboration XVI (2016), there are two aspects that make inflationary model building difficult. First is the problem of obtaining a significant amplitude of dipole modulation. In Planck Collaboration XVI (2016) the asymmetry was found to have amplitude A ≈ 6−7% on scales 2 ≤ ℓ ≤ 64. This compares with the expected value of A = 2.9% on these scales due to cosmic variance in statistically isotropic skies. One basic strategy for incorporating the violation of statistical isotropy into inflation is to consider some form of multifield inflation and use one of the directions orthogonal to the direction of slow roll as the field responsible for the modulation. Obtaining the required modulation is problematic because most extra fields in multifield inflation become disordered in a nearly scaleinvariant way, just like the fluctuations in the field parallel to the direction of slow roll. What is needed resembles a pure gradient with no fluctuations of shorter wavelength. In Liddle & Cortês (2013) it was proposed that such a field could be produced using the supercurvature mode of open inflation. (See however the discussion in Kanno et al. 2013.) Also, in order to respect the f_{NL} constraints, one must avoid that the modulating field leave a direct imprint on the temperature anisotropy.
The second aspect which makes model building difficult for dipolar asymmetry is that the measured amplitude is strongly scale dependent, and on scales ℓ ≳ 100 no significant detection of a dipolar modulation amplitude is made (Planck Collaboration XVI 2016), once our proper motion has been taken into account (Planck Collaboration XXVII 2014). On the other hand, the simplest models are scalefree and produce statistical anisotropy of the type described by the ansatz in Eq. (156), for which the bulk of the statistical weight should be detected at the resolution of the survey. To resolve this difficulty, Erickcek et al. (2009) proposed modulating CDI fluctuations generated within the framework of a curvaton scenario, because, unlike adiabatic perturbations, CDI perturbations entering the Hubble radius before last scattering contribute negligibly to the CMB fluctuations (recall Fig. 43).
The situation for the quadrupolar power asymmetry is different from the dipolar case in that no detection is currently claimed. Model building is easier than the dipolar case since no pure gradient modes are required, but also more difficult in that anisotropy during inflation is needed. While the isotropy of the recent expansion of the Universe (i.e., since the CMB fluctuations were first imprinted) is tightly constrained, bounds on a possible anisotropic expansion at early times are much weaker. Ackerman et al. (2007) proposed using constraints on the quadrupolar statistical anisotropy of the CMB to probe the isotropy of the expansion during inflation – that is, during the epoch when the perturbations now seen in the CMB first exited the Hubble radius. Assuming an anisotropic expansion during inflation, Ackerman et al. (2007) computed its impact on the threedimensional power spectrum on superHubble scales by integrating the mode functions for the perturbations during inflation and beyond. Several sources of such anisotropy have been proposed, such as vector fields during inflation (Dimastrogiovanni et al. 2010; Soda 2012; Maleknejad et al. 2013; Schmidt & Hui 2013; Bartolo et al. 2013; Naruko et al. 2015), or an inflating solid or elastic medium (Bartolo et al. 2013).
12.2. Scaledependent modulation and idealized estimators
The ansatz in Eq. (156) expressed in angular space may be rewritten in terms of the multipole expansion and generalized to include scaledependent modulation by means of Wigner 3j symbols: (157)Because of the symmetry of the lefthand side, the coefficients acquire a phase (− 1)^{ℓ + ℓ′ + L} under interchange of ℓ and ℓ′. This is the most general form consistent with the hypothesis of Gaussianity. The usual isotropic power spectrum, which is the generic prediction of simple models of inflation, includes only the L = 0 term, where and the Wigner 3j symbol provides the δ_{ℓ,ℓ′}δ_{m,m′} factor. The coefficients with L> 0 introduce statistical anisotropy.
If we assume that there is a common vector (corresponding to L = 1 on the celestial sphere) that defines the direction of the anisotropy of the power spectrum for all the terms of L = 1, we may adopt a more restricted ansatz for the bipolar modulation, so that (158)where we assume that X_{M} is normalized (i.e., ∑ _{M}X_{M}X_{M}^{∗} = 1). In such a model, supposing that is theoretically determined, but the orientation of the unit vector X_{M} is random and isotropically distributed on the celestial sphere, we may construct the following quadratic estimator for the direction: (159)where the weights for the unbiased minimum variance estimator are given by (160)This construction, which for the L = 1 case may be found in Moss et al. (2011) and Planck Collaboration XVI (2016), may be readily generalized to L> 1 in the above way.
12.3. Constraining inflationary models for dipolar asymmetry
In this section, we confront with Planck data the modulated curvaton model of Erickcek et al. (2009), which attempts to explain the observed largescale power asymmetry via a gradient in the background curvaton field. In this model, the curvaton decays after CDM freeze out, which results in a nearlyscaleinvariant isocurvature component between CDM and radiation. In the viable version of this scenario, the curvaton contributes negligibly to the CDM density. A longwavelength fluctuation in the curvaton field initial value σ_{∗} is assumed, with amplitude Δσ_{∗} across our observable volume. This modulates the curvaton isocurvature fluctuations according to S_{σγ} ≈ 2δσ_{∗}/σ_{∗}. The curvaton produces all of the final CDI fluctuations, which are nearly scale invariant, as well as a component of the final adiabatic fluctuations. Hence both of these components will be modulated, and the parameter space of the model will be constrained by observations of the power asymmetry on large and small scales, as well as the fullsky CDI fraction. In practice, the very tight constraints on small scale power asymmetry obtained in Planck Collaboration XVI (2016) imply a small curvaton adiabatic component, which implies that the CDI and adiabatic fluctuations are only weakly correlated. This model easily satisfies constraints due to the CMB dipole, quadrupole, and nonGaussianity (Erickcek et al. 2009).
There are two main parameters that we constrain for this model. First, the fraction of adiabatic fluctuations due to the curvaton ξ is defined as (161)Here, and are the inflaton and curvaton primordial power spectra, respectively, and Σ_{σ} is the coupling from curvaton isocurvature to adiabatic fluctuations. (Up to a sign, ξ is equal to the correlation parameter.) Next, the coupling of curvaton to CDI, M_{CDIσ}, is determined by the constant κ ≡ M_{CDIσ}/R ≳ −1, where (162)and all density parameters are evaluated just prior to curvaton decay. The isocurvature fraction can be written in terms of these two parameters by (163)These parameters determine the modulation of the CMB temperature fluctuations via ΔC_{ℓ}/C_{ℓ} = 2K_{ℓ}Δσ_{∗}/σ_{∗}, where (Erickcek et al. 2009) (164)Here , , and are the adiabatic, CDI, and correlated power spectra calculated for unity primordial spectra.
Note that this modulated curvaton model contains some simple special cases. For κ = 0, we have a purely adiabatic (i.e., scaleinvariant) modulation. This is equivalent to a modulation of the scalar amplitude, A_{s}. On the other hand, if we take the limit κ → ∞, with fixed κ^{2}ξ (i.e., with fixed isocurvature fraction β_{iso}), we obtain a pure CDI modulation. For κ = ξ = 0 we have no modulation, i.e., we recover base ΛCDM. Therefore this model is particularly useful for examining a range of possible modulations within the context of a concrete framework.
In order to constrain this model, we use a formalism which was developed to determine the signatures of potential gradients in physical parameters in the CMB (Moss et al. 2011), and which is used to examine dipolar modulation and described in detail in Planck Collaboration XVI (2016). This approach is wellsuited to testing the modulated curvaton model since it can accommodate scaledependent modulations. Briefly, we write the temperature anisotropy covariance given a gradient ΔX_{M} in a parameter X as (165)where δC_{ℓℓ′} ≡ dC_{ℓ}/ dX + dC_{ℓ′}/ dX. Note that this covariance takes the form of Eqs. (157) and (158), with (166)We then construct a maximum likelihood estimator for the gradient components. We use C^{1} filtered data (Planck Collaboration XV 2016) and perform a meanfield subtraction, giving (167)Here f_{1M} is a normalization correction due to the applied mask, M(Ω), and is given by (168)The T_{ℓm} are the filtered data and . Note that the lack of aberration in the Planck Full Focal Plane simulations (Planck Collaboration XII 2016) is expected to have negligible effect on this analysis and on that of the quadrupolar modulation in the next subsection, since the CDI modulation is heavily suppressed for ℓ ≳ 500, whereas the effect of aberration has a very different ℓ dependence and will bias the modulation signal for ℓ ≲ 1000 at an insignificant level.
In practice, exploring the parameter space of the model is sped up dramatically by binning the estimator defined in Eq. (167) into bins of width Δℓ = 1, which means that the estimators only need to be calculated once (Planck Collaboration XVI 2016). Finally, for the modulated curvaton model we identify (169)Note that for our constraints we fix the curvaton gradient to its maximum value, Δσ_{∗}/σ_{∗} = 1. Therefore, our constraints are conservative, since smaller Δσ_{∗}/σ_{∗} would only reduce the modulation that this model could produce.
The temperature anisotropies measured by Planck constrain the modulated curvaton parameters κ and ξ via Eqs. (164) and (167). Figure 54 shows the constraints in this parameter space evaluating the estimator to ℓ_{max} = 1000. The maximum likelihood region corresponds to a band at κ ≳ 3. For parameters in this region, the model produces a largescale asymmetry via a mainlyCDI modulation. However, the amplitude of this largescale asymmetry is lower than the 6–7% actually observed (Planck Collaboration XVI 2016). The reason is that, had a CDI modulation produced all of the largescale asymmetry, the consequent smallscale asymmetry (due to the shape of the scaleinvariant CDI spectrum) would be larger than the Planck observations allow. The allowed CDI modulation is further reduced by the Planck 95% upper limit on an uncorrelated, scaleinvariant (“axion”type) isocurvature component, β_{iso}< 0.033, from Sect. 11. Imposing this constraint reduces the available parameter space in the κ–ξ plane via Eq. (163), as illustrated in Fig. 54.
Fig. 54 68% and 95% CL regions in the modulated isocurvature model parameter space using the Planck temperature data up to ℓ_{max} = 1000 (contours). The region above the dashed curve is ruled out by the Planck constraint on an uncorrelated, scaleinvariant isocurvature component. 
The best fit in Fig. 54 corresponds to Δχ^{2} = −6.8 relative to base ΛCDM, for two extra parameters. In order to assess how likely such an improvement would be in statistically isotropic skies, we note that the bestfit CDI modulation amplitude is very close to the mean amplitude expected due to cosmic variance, as calculated directly from Eq. (167). More precisely, since the amplitude is χ^{2} distributed with three degrees of freedom, i.e., MaxwellBoltzmann distributed, we conclude that about 44% of statistically isotropic skies will exhibit a measured (via Eq. (167)) isocurvature modulation larger than that of the actual sky.
To summarize, the modulated curvaton model can only produce a small part of the observed largescale asymmetry, and what it can produce is entirely consistent with cosmic variance in a statistically isotropic sky. Hence we must favour the base ΛCDM model over this model. Finally, note that further generalizing the model (e.g., to allow nonscaleinvariant CDI spectra) may allow more largescale asymmetry to be produced and hence result in an improved Δχ^{2}, at the expense of more parameters. On the other hand, the neutrino isocurvature modes are not expected to fit the observed asymmetry well due to their approximate scale invariance (see Fig. 43).
12.4. Constraints on quadrupolar asymmetry generated during inflation
In this section we assume a quadrupolar directional dependence of the primordial scalar power spectrum about some axis and having a scaledependent amplitude g(k). More specifically, we assume (170)which can be rewritten as (171)where (172)with g_{2M}(k) satisfying . In this analysis, we will treat the modulation scale dependence as a power law, g(k) = g_{∗}(k/k_{∗})^{q}, and consider five values of the spectral index, namely q = −2, −1, 0, 1, and 2. Importantly, for q ≠ 0 our constraints on g_{∗} will depend on the pivot scale, chosen as k_{∗} = 0.05 Mpc^{1} as elsewhere in this paper. Models have been proposed predicting both positive and negative g_{∗} (see, e.g. Tsujikawa 2014), so we keep the sign of g_{∗} free.
Often in the literature the term −g(k) / 3 is not included in the modulated power spectrum, Eq. (170). Our form sets the modulation monopole to zero, so that there is no correction to the isotropic power spectrum dependent on g(k). We do this because for large  q  the correction would require a joint analysis with the isotropic power spectrum likelihood. Inflationary models have been proposed which predict both forms. For example, the model in Ohashi et al. (2013) includes the modulation monopole, while the model in Libanov & Rubakov (2010) does not. For q = 0 our results apply to both forms due to the degeneracy of a scaleindependent correction to with the scalar amplitude, A_{s}. However, for nonzero tilt a joint analysis would yield tighter constraints on g_{∗} when the monopole correction is present, in which case our results will be conservative.
Given the anisotropic power spectrum of Eq. (171), the statistically anisotropic part of the CMB temperature fluctuations has the following expectation value (Ma et al. 2011): (173)where and denotes the temperature radiation transfer function.
The analysis is carried out using the foregroundcleaned CMB temperature maps Commander, NILC, SEVEM, and SMICA, where we apply the extended common mask UT76. The implementation details of the optimal estimator can be found in Sect. 5.3 of Planck Collaboration XVI (2016). However, here we apply an inverse variance weighted filter that assumes a simple white noise component, but optimally accounts for the mask in the same manner as Planck Collaboration XVII (2014) and Sects. 6.3 and 6.6 of Planck Collaboration XVI (2016). We estimate g_{2M} from the data at multipoles 2 ≤ ℓ ≤ 1200. The range of multipoles is chosen such that the impact of foreground residuals on the conclusions is insignificant. Neglecting very small scales, however, sacrifices little constraining power because those scales are noise dominated. This conclusion was based on realistic simulations containing residual foregrounds. Moreover, we estimate the statistical uncertainty in g_{2M} with various ℓ_{max} values using simulations.
Once we have obtained estimates for the five g_{2M} coefficients, we must determine values for the model parameters of interest, namely g_{∗} and . We assume the g_{2M} coefficients to be Gaussian distributed due to cosmic variance. We have explicitly checked this hypothesis with simulations. Hence the likelihood function is (174)where G is the g_{2M} covariance matrix, which is estimated using isotropic simulations. One approach to determining the model parameters would be to use this likelihood to calculate marginalized posterior distributions for g_{∗}, from which mean values and errors could be determined. However, we find that g_{∗} is so poorly constrained that the means and widths thus calculated strongly depend on the prior for g_{∗}. Two sensible priors are uniform in g_{∗} or proportional to [i.e., uniform in the Cartesian components of ]. In addition, we find that the posterior means are much closer to zero then the widths, which is due to the approximate degeneracy between a modulation and modulation ′), where ′ is orthogonal to . In such a situation the degree of consistency between the measured value of g_{∗} and the expectations of cosmic variance in statistically isotropic skies is unclear.
Instead we determine bestfit values for g_{∗} and by maximizing the likelihood over the three parameters using a grid approach. To characterize how unexpected our bestfit values are in statistically isotropic skies, we repeat the procedure replacing our estimates for g_{2M} from the data with estimates from 1000 isotropic simulations. We finally calculate pvalues, which give the fraction of simulations with a larger value of  g_{∗}  than the actual data. Note that from the Bayesian perspective the maximumlikelihood values amount to maximizing the posterior for g_{∗} given a uniform prior on g_{∗}, so that these values will change with a different prior. However, we have checked that the pvalues depend only very weakly on the choice of prior.
Minimumχ^{2}g_{∗} values for quadrupolar modulation, determined from the Commander, NILC, SEVEM, and SMICA foregroundcleaned maps.
Table 17 shows the g_{∗} values obtained by minimizing χ^{2} as well as the pvalues for the data compared to statistically isotropic simulations. Note that the constraints on g_{∗} are strongest for the most negative values of the exponent q. This is because for fixed g_{∗} the largest asymmetry over the range of observable scales occurs for q = −2, at the largest scales, due to the location of the pivot scale, k_{∗}. Our limits provide a stringent test of rotational symmetry during inflation. We find no sign of deviation from statistical isotropy.^{25}
13. Combination with BICEP2/Keck ArrayPlanck crosscorrelation
In this section we discuss the implications of the recent constraints on the primordial Bmode polarization from the crosscorrelation of the BICEP2 and Keck Array data at 150 GHz with the Planck maps at higher frequencies to characterize and remove the contribution from polarized thermal dust emission from our Galaxy (BKP). On its own, the BKP likelihood leads to a 95% CL upper limit of r< 0.12, compatible with and independent of the constraints obtained using the 2015 Planck temperature and large angular scale polarization in Sect. 5. (Note, however, that the BKP likelihood uses the HamimecheLewis approximation (Hamimeche & Lewis 2008), which requires the assumption of a fiducial model.) The BKP results are also compatible with the Planck 2013 Results (Planck Collaboration XVI 2014; Planck Collaboration XXII 2014). The posterior probability distribution for r obtained by BKP peaks away from zero at r ≈ 0.05, but the region of large posterior probability includes r = 0.
Fig. 55 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 its crosscorrelation with BICEP2/Keck Array and/or BAO data compared with the theoretical predictions of selected inflationary models. Note that the marginalized joint 68% and 95% CL regions have been obtained by assuming dn_{s}/ dlnk = 0. 
Here we combine the baseline twoparameter BKP likelihood using the lowest five Bmode bandpowers with the Planck 2015 likelihoods. The two BKP nuisance parameters are the Bmode amplitude and frequency spectral index of the polarized thermal dust emission. The combined analysis yields the following constraint on the tensortoscalar ratio: (175)further improving on the upper limits obtained from the different data combinations presented in Sect. 5.
By directly constraining the tensor mode, the BKP likelihood removes degeneracies between the tensortoscalar ratio and other parameters. Adding tensors and running, we obtain (176)which constitutes almost a 50% improvement over the Planck TT+lowP constraint quoted in Eq. (27). These limits on tensor modes are more robust than the limits using the shape of the spectrum alone because scalar perturbations cannot generate B modes irrespective of the shape of the scalar spectrum.
13.1. Implications of BKP on selected inflationary models
Using the BKP likelihood further strengthens the constraints on the inflationary parameters and models discussed in Sect. 6, as seen in Fig. 55. If we set ϵ_{3} = 0, the first slowroll parameter is constrained to ϵ_{1}< 0.0055 at 95% CL by Planck TT+lowP+BKP. With the same data combination, concave potentials are preferred over convex potentials with lnB = 3.8, which improves on the lnB = 2 result obtained from the Planck data alone.
Results of inflationary model comparison using the crosscorrelation between BICEP2/Keck Array and Planck.
Combining with the BKP likelihood strengthens the constraints on the selected inflationary models studied in Sect. 6. Using the same methodology as in Sect. 6 and adding the BKP likelihood gives a Bayes factor preferring R^{2} over chaotic inflation with monomial quadratic potential and natural inflation by odds of 403:1 and 270:1, respectively, under the assumption of a dust equation of state during the entropy generation stage. The combination with the BKP likelihood further penalizes the doublewell model compared to R^{2} inflation. However, adding BKP reduces the Bayes factor of the hilltop models compared to R^{2}, because these models can predict a value of the tensortoscalar ratio that better fits the statistically insignificant peak at r ≈ 0.05. See Table 18 for the Δχ^{2} and the Bayes factors of inflationary models with the same two cases of postinflationary evolution studied in Sect. 6. Note, however, that the Δχ^{2} are computed with respect to the best fit of baseline + tensors, unlike in Table 7.
13.2. Implications of BKP on scalar power spectrum
Fig. 56 Impact of BKP likelihood on scalar primordial power spectrum reconstruction. We show how including the BKP likelihood affects the reconstruction in Sect. 8.3. The top panel is to be compared with the reconstructions in Fig. 27, and we observe that including BKP has a minimal impact given the uncertainty in the reconstruction. The middle panel is to be compared with Fig. 31, and here we notice that including BKP excludes the trajectories with large values of ϵ. The bottom panel shows how the inflationary potential reconstructions are modified by BKP (to be compared with Fig. 32). 
The presence of tensors would, at least to some degree, require an enhanced suppression of the scalar power spectrum on large scales to account for the lowℓ deficit in the spectrum. We therefore repeat the analysis of an exponential cutoff studied in Sect. 4.4 with tensor perturbations included and the standard tensor tilt (i.e., n_{t} = −r/ 8). Allowing tensors does not significantly degrade the Δχ^{2} improvement found in Sect. 4.4 for Planck TT+lowP with a best fit at r ≈ 0. When the BKP likelihood is combined, we obtain Δχ^{2} = −4 with respect to the base ΛCDM model with a best fit at r ≈ 0.04. However, since this model contains 3 additional parameters, it is not preferred over base ΛCDM.
In Fig. 56 we show how the scalar primordial power spectrum reconstruction discussed in Sect. 8.3 is modified when the BKP likelihood is also included. While the power spectrum reconstruction hardly varies given the uncertainties in the method, the trajectories of the slowroll parameters are significantly closer to slow roll. When the 12knot reconstruction is carried out, the upper bound on the tensortoscalar ratio is r< 0.11 at 95% CL. The χ^{2} per degree of freedom for the 5 lowk and 6 highk knots are 1.14 and 0.22, respectively, corresponding to pvalues of 0.33 and 0.97.
13.3. Relaxing the standard singlefield consistency condition
We now relax the consistency condition (i.e., n_{t} = −r/ 8) and allow the tensor tilt to be independent of the tensortoscalar ratio. This fully phenomenological analysis with the BKP likelihood is complementary to the study of inflationary models with generalized Lagrangians in Sect. 10, which also predict modifications to the consistency condition n_{t} = −r/ 8 for a nearly scaleinvariant spectrum of tensor modes. In this subsection we adopt a phenomenological approach, thereby including radical departures from n_{t} ≲ 0, including values which are predicted in alternative models to inflation (Gasperini & Veneziano 1993; Boyle et al. 2004; Brandenberger et al. 2007). In Sect. 10 we folded in the Planckf_{NL} constraints (Planck Collaboration XVII 2016), whereas here we consider Planck and BKP likelihoods only. Complementary probes such as pulsar timing, direct detection of gravitational waves, and nucleosynthesis bounds could be used to constrain blue values for the tensor spectral index (Stewart & Brandenberger 2008), but here we are primarily interested in what CMB data can tell us.
We caution the reader that in the absence of a clear detection of a tensor component, joint constraints on r and n_{t} depend strongly on priors, or equivalently on the choice of parameterization. Nevertheless, the BKP likelihood has some constraining power over a range of scales more than a decade wide around k ≈ 0.01 Mpc^{1}, so the results are not entirely prior driven.
Fig. 57 Posterior probability density of the tensortoscalar ratio at two different scales. The inflationary consistency relation is relaxed and r_{0.002} and r_{0.020} are used as sampling parameters, assuming a powerlaw spectrum for primordial tensor perturbations. When the BKP likelihood is included in the analysis, the results with Planck TT+lowP+BAO and Planck TT, TE, EE+lowP coincide (dashed and solid red curves, respectively). 
Fig. 58 68% and 95% CL constraints on tensors when the inflationary consistency relation is relaxed, with Planck TT+lowP+BAO (blue dashed contours) and TT, TE, EE+lowP (blue shaded regions). The red colours are for the same data plus the BKP joint likelihood. The upper panel shows our independent primary parameters r_{0.002} and r_{0.020}. The lower panel shows the derived parameters n_{t} and r_{0.01}. The scale k = 0.01 Mpc^{1} is near the decorrelation scale of (n_{t}, r) for the Planck+BKP data. 
The commonly used (r, n_{t}) parameterization suffers from pathological behaviour around r = 0, which could be problematic for statistical sampling. We therefore use a parameterization specifying r at two different scales, (r_{k1},r_{k2}) (analogous to the treatment of primordial isocurvature in Sect. 11) as well as the more familiar (r, n_{t}) parameterization. We present results based on k_{1} = 0.002 Mpc^{1} and k_{2} = 0.02 Mpc^{1}, also quoting the amplitude at k = 0.01 Mpc^{1} for both parameterizations. This scale is close to the decorrelation scale for (r, n_{t}) for the Planck+BKP joint constraints. We obtain r_{0.002}< 0.07 (0.06) and r_{0.02}< 0.29 (0.31) at 95% CL from the twoscale parameterization with Planck TT+lowP+BAO+BKP (TT, TE, EE+lowP+BKP). Figure 57 illustrates the impact of the BKP likelihood on the onedimensional posterior probabilities for these two parameters. The derived constraint at k = 0.01 Mpc^{1} is r_{0.01}< 0.12 (0.12) at 95% CL with Planck TT+lowP+BAO+BKP (TT, TE, EE+lowP+BKP). The upper panel of Fig. 58 shows the relevant 2D contours for the tensortoscalar ratios at the two scales and the improvement due to the combination with the BKP likelihood. The lower panel shows the 2D contours in (r_{0.01},n_{t}) obtained by sampling with the twoscale parameterization. Figure 59 shows the 2D contours in (r_{0.01},n_{t}) obtained by the (r_{0.002},n_{t}) parameterization.
We conclude that positive values of the tensor tilt, n_{t}, are not statistically significantly preferred by the BKP joint measurement of Bmode polarization in combination with Planck data, a conclusion at variance with results reported using the BICEP2 data (Gerbino et al. 2014). However, the now firmly established contamination by polarized dust emission easily explains the discrepancy. Values of tensor tilt consistent with the standard singlefield inflationary consistency relation are compatible with the Planck+BKP constraints.
14. Conclusions
The Planck full mission temperature and polarization data are consistent with the spatially flat base ΛCDM model whose perturbations are Gaussian and adiabatic with a spectrum described by a simple power law, as predicted by the simplest inflationary models. For this release, the basic Planck results do not rely on external data. The first Planck polarization release at large angular scales from the LFI 70 GHz channel determines an optical depth of τ = 0.067 ± 0.022 (68% CL, Planck low multipole likelihood), a value smaller than the previous Planck 2013 result based on the WMAP9 polarization likelihood as delivered by the WMAP team. This Planck value of τ is consistent with an analysis of WMAP9 polarization data cleaned for polarized dust emission using the Planck 353 GHz data (Planck Collaboration XV 2014; Planck Collaboration XI 2016). The estimates of cosmological parameters from the full mission temperature data and polarization on large angular scales are consistent with those of the Planck 2013 release. The TE and EE spectra at ℓ ≥ 30 together with the lensing power spectrum lead to cosmological constraints in agreement with those obtained from temperature.
The Planck full mission temperature and largeangularscale polarization data rule out an exactly scaleinvariant spectrum of curvature perturbations at 5.6 σ. For the base ΛCDM model, the spectral index is measured to be n_{s} = 0.965 ± 0.006 (68% CL, Planck TT+lowP). No evidence for a running of the spectral index is found, with dn_{s}/ dlnk = −0.008 ± 0.008 (68% CL, Planck TT+lowP). By considering Planck TT+lowP+lensing we obtain n_{s} = 0.968 ± 0.006 and dn_{s}/ dlnk = −0.003 ± 0.007, both at 68% CL.
The Planck full mission data improve the upper bound on the tensortoscalar ratio to r_{0.002}< 0.10 (95% CL, Planck TT+lowP), a bound that changes only slightly when including the Planck lensing likelihood, the highℓ polarization likelihood, or the likelihood from the WMAP largeangularscale polarization map (dustcleaned with the Planck 353 GHz map). We showed how the lowℓ deficit in temperature contributes to the Planck upper bound on r_{0.002}, but this deficit is not a statistically significant anomaly within the base ΛCDM cosmology. Using the full mission Planck data, we find the upper bound on r_{0.002} stable, even when extended cosmological models or models with CDM isocurvature are considered. The Planck bound on r_{0.002} is consistent with the recent result r_{0.002}< 0.12 at 95% CL obtained by the BICEP2/Keck ArrayPlanck crosscorrelation analysis (BKP) which provides an estimate for the contamination from polarized dust emission (Planck Collaboration XXX 2014). By combining Planck TT+lowP with the BKP crosscorrelation likelihood, we obtain r_{0.002}< 0.08 at 95% CL.
The increased precision of the Planck full mission data reduces the area enclosed by the 95% confidence contour in the (n_{s},r) plane by 29%. We performed a Bayesian model comparison with the same methodology as in PCI13, taking into account reheating uncertainties by marginalizing over two extra parameters: the energy scale at thermalization, ρ_{th}, and the parameter w_{int} characterizing the average equation of state between the end of inflation and thermalization. Among the models considered using this approach, the R^{2} inflationary model proposed by Starobinsky (1980) is the most favoured. Due to its high tensortoscalar ratio, the quadratic model is now strongly disfavoured with respect to R^{2} inflation for Planck TT+lowP in combination with BAO data. By further including the BKP likelihood, this conclusion is confirmed, and natural inflation is also disfavoured.
We reconstructed the inflaton potential and the Hubble parameter evolution during the observable part of inflation using a Taylor expansion of the inflaton potential or H(φ). This analysis did not rely on the slowroll approximation, nor on any assumption about the end of inflation. When higherorder terms were allowed, both reconstructions led to a change in the slope of the potential at the beginning of the observable range, thus better fitting the lowℓ temperature deficit by turning on a nonzero running of running and accommodating r_{0.002} ≈ 0.2. These models, however, are not significantly favoured compared to lowerorder parameterizations that lead to slowroll evolution at all times.
Three distinct methods were used to reconstruct the primordial power spectrum. All three methods strongly constrain deviations from a featureless power spectrum over the range of scales 0.008 Mpc^{1} ≲ k ≲ 0.1 Mpc^{1}. More interestingly, they also independently find common patterns in the primordial power spectrum of curvature perturbations at k ≲ 0.008 Mpc^{1}. These patterns are related to the dip at ℓ ≈ 20–40 in the temperature power spectrum. This deviation from a simple powerlaw spectrum has weak statistical significance due to the large cosmic variance at low ℓ.
This direct reconstruction of the power spectrum is complemented by a search for parameterized features in physically motivated models. The models considered range from the minimal case of a kineticenergydominated phase preceding a short inflationary stage (with just one extra parameter), to a model with a steplike feature in the potential and in the sound speed (with five extra parameters). As with the Planck 2013 nominal mission data, these templates lead to an improved fit, up to Δχ^{2} ≈ 12. However, neither Bayesian model comparison nor a frequentistsimulationbased analysis shows any statistically significant preference over a simple power law.
We have updated the analysis that combines power spectrum constraints with those derived from the f_{NL} parameters (Planck Collaboration XVII 2016). New limits on the sound speed inferred from the full mission temperature and polarization data further constrain the slowroll parameters for generalized models, including DBI inflation. For the first time, we derived combined constraints on Galileon inflation, including the region of parameter space in which the predicted spectrum of gravitational waves has a blue spectral index.
Several models motivated by the axion monodromy mechanism in string theory predict oscillatory modulations and corresponding nonGaussianities, potentially detectable by Planck. A TTonly analysis picks up four possible modulation frequencies, which remain present when the highℓ polarization likelihood is included. An inspection of frequency residuals in the highℓTT likelihood does not reveal evidence of foregroundrelated systematics at similar frequencies. However, a Bayesian model comparison analysis prefers the smooth base ΛCDM model over modulated models, suggesting that the latter could simply be fitting the noise in the data. The monodromy model predicts resonant nonGaussian features correlated to power spectrum features. A partial analysis beyond the power spectrum was presented. We also constrained a possible pseudocoupling of the axion to gauge fields by requiring that nonGaussianities induced by inverse decay satisfy the Planck bounds on f_{NL}.
Section 11 reports on a search for possible deviations from purely adiabatic initial conditions by studying a range of models including isocurvature modes as well as possible correlations with the adiabatic mode. The Planck full mission temperature data are consistent with adiabaticity. The PlanckTT data place tight constraints on threeparameter extensions to the flat adiabatic base ΛCDM model, allowing arbitrarilycorrelated mixtures of the adiabatic mode with one isocurvature mode (of either the CDM, baryon, neutrino density, or neutrino velocity type). Adding the highℓTE and EE polarization data further squeezes the constraints, since polarization spectra contain additional shape and phase information on acoustic oscillations. The likelihood with polarization included is in agreement with adiabatic initial conditions. However, the tightening of the constraints after including polarization must be interpreted with caution because of possible systematic effects. For this reason we emphasize the more conservative Planck TT+lowP bounds in Table 16. The constraints on the six baseΛCDM cosmological model parameters remain stable when correlated isocurvature modes are allowed. The largest shifts occur for the neutrino density mode, but these shifts are not significant (i.e., are below 1σ). The constraints on the tensortoscalar ratio also remain stable when isocurvature modes are allowed.
Finally we examined the connection between inflation and statistical isotropy, a key prediction of the simplest inflationary models. We tested separately the two lowest moments of an anisotropic modulation of the primordial curvature power spectrum. We found that a modulated curvaton model proposed to explain the observed largescale dipolar power asymmetry cannot account for all of the asymmetry, and hence is not preferred over statistically isotropic base ΛCDM. The full mission temperature data place the tightest constraints to date on a quadrupolar modulation of curvature perturbations.
Planck (http://www.esa.int/Planck) is a project of the European Space Agency (ESA) with instruments provided by two scientific consortia funded by ESA member states and led by Principal Investigators from France and Italy, telescope reflectors provided through a collaboration between ESA and a scientific consortium led and funded by Denmark, and additional contributions from NASA (USA).
In this paper we use the conventions introduced in Planck Collaboration XIII (2016). We adopt the following labels for likelihoods: (i) Planck TT denotes the combination of the TT likelihood at multipoles ℓ ≥ 30 and a lowℓ temperatureonly likelihood based on the CMB map recovered with Commander; (ii) Planck TT–lowT denotes the TT likelihood at multipoles ℓ ≥ 30; (iii) Planck TT+lowP further includes the Planck polarization data in the lowℓ likelihood, as described in the main text; (iv) PlanckTE denotes the likelihood at ℓ ≥ 30 using the TE spectrum; and (v) Planck TT, TE, EE+lowP denotes the combination of the likelihood at ℓ ≥ 30 using TT, TE, and EE spectra and the lowℓ multipole likelihood. The label “τ prior” denotes the use of a Gaussian prior τ = 0.07 ± 0.02. The labels “lowT, P” and “lowEB” denote the lowℓ multipole likelihood and the Q,U pixel likelihood only, respectively.
The definition of f(k) used here differs from that of PCI13 in that exp(f) is used in place of 1 + f to ensure that the reconstructed primordial power spectrum is always nonnegative.
We use the more accurate relation (88) accounting for different epochs of freezeout for the scalar fluctuations (at sound horizon crossing, kc_{s} = aH) and tensor perturbations (at Hubble radius crossing, k = aH; Peiris et al. 2007; Powell et al. 2009; Lorenz et al. 2008; Agarwal & Bean 2009; Baumann et al. 2015).
This section uses results based on f_{NL} constraints from T and E (Planck Collaboration XVII 2016). In Planck Collaboration XVII (2016) it is shown that, although conservatively considered as preliminary, the f_{NL} constraints from T and E are robust, since they pass an extensive battery of validation tests and are in full agreement with Tonly constraints.
See also Mizuno & Koyama (2010), Gao & Steer (2011), Kobayashi et al. (2011b), De Felice & Tsujikawa (2013), and Regan et al. (2015).
The transfer function mapping the primordial CDI mode to is suppressed by a factor (k/k_{eq})^{2} ~ (ℓ/ℓ_{eq})^{2} relative to the ADI mode, where k_{eq} is the wavenumber of matterradiation equality. As seen in Fig. 43, there is a similar damping for the E mode in the CDI versus the ADI case. Therefore, to be observable at high ℓ, a CDI mode should be (highly) blue tilted. So, if the data favoured as small as possible a disturbance by CDI over all scales, then the CDI should have a spectral index, n_{ℐℐ}, of roughly three. In practice, the lowestℓ part of the data has very little weight due to cosmic variance, and thus we expect that the data should favour n_{ℐℐ} less than three, but significantly larger than one. This should be kept in mind when interpreting the results in the CDI case, i.e., one cannot expect strong constraints on the primordial CDI fraction at small scales, even if the data are purely adiabatic. The imprint of the baryon density isocurvature (BDI) mode in the CMB, at least at linear order, is indistinguishable from the CDI case, and hence we do not consider it separately as it can be described by . The trispectrum, however, can in principle be used to distinguish the BDI and CDI modes (Grin et al. 2014).
However, conversely, if no isocurvature was detected, the fluctuations could have been seeded either by single or multifield inflation, since later processes easily wash out inflationary isocurvature perturbations (Mollerach 1990; Weinberg 2004; Beltrán et al. 2005). An example is the curvaton model, in which perturbations can be purely isocurvature at Hubble exit during inflation, but are later converted to ADI if the curvaton or curvaton particles (Linde & Mukhanov 2006) dominate the energy density at the curvaton’s decay. For a summary of various curvaton scenarios, see, e.g., Gordon & Lewis (2003).
Given our ansatz of powerlaw primordial spectra, if we treated as an independent parameter as we do with , Eq. (131) would always be violated somewhere outside [ k_{1},k_{2} ]. In PCI13, we dealt with this by assuming that when maximal (anti)correlation is reached at some scale, the correlation remains at (−)100% beyond this scale. This introduced a kink in the crosscorrelation spectrum, located at a different wavenumber for each model. Even though the range [ k_{1},k_{2} ] was chosen to span most of the observable scales, this kink tended to impact the smallest (or largest) multipole values used in the analysis. In particular, the kink helped fit the dip in the temperature angular power in the multipole range ℓ ≈ 10–40.
Indeed, if curvaton particles are produced during reheating, they can be expected to survive and outweigh other particles at the moment of curvaton decay, but by how much depends on the details of the model. As the curvaton field (during its oscillations) and the curvaton particles have the same equation of state and they decay simultaneously, no isocurvature perturbations are produced.
The formula is often quoted or utilized, particularly in the older curvaton literature. This result, which follows from considering only squares of first order perturbations, is valid when r_{D} is close to zero (i.e., when is very large). However, when r_{D} is close to unity or , which is the case with the Planck measurements, the second and third terms of Eq. (146) are crucial. These follow from second order perturbation theory calculations. Coincidentally, if one erroneously uses the expression 5 / (4r_{D}) in the limit r_{D} → 1, one obtains the result + 5 / 4, whereas the correct formula (146) with leads to −5 / 4 when r_{D} → 1.
However, if was nonnegligible, then all the constraints on would shift upward. For example, with , our constraints on β_{iso} would translate to . On the other hand, the Planck constraint of can be converted to an upper bound (95% CL from T only) as shown in Planck Collaboration XVII (2016).
The constraints from the Planck 2013 data by Kim & Komatsu (2013) should be multiplied by a factor of in our normalization.
Acknowledgments
The Planck Collaboration acknowledges the support of: ESA; CNES and CNRS/INSUIN2P3INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MINECO, JA, and RES (Spain); Tekes, AoF, and CSC (Finland); DLR and MPG (Germany); CSA (Canada); DTU Space (Denmark); SER/SSO (Switzerland); RCN (Norway); SFI (Ireland); FCT/MCTES (Portugal); ERC and PRACE (EU). A description of the Planck Collaboration and a list of its members, indicating which technical or scientific activities they have been involved in, can be found at http://www.cosmos.esa.int/web/planck/planckcollaboration. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the US Department of Energy under Contract No. DEAC0205CH11231. Part of this work was undertaken at the STFC DiRAC HPC Facilities at the University of Cambridge, funded by UK BIS National Einfrastructure capital grants. We gratefully acknowledge the IN2P3 Computer Center (http://cc.in2p3.fr) for providing a significant amount of the computing resources and services needed for this work.
References
 Achúcarro, A., Gong, J.O., Hardeman, S., Palma, G. A., & Patil, S. P. 2011, JCAP, 1101, 030 [NASA ADS] [CrossRef] [Google Scholar]
 Ackerman, L., Carroll, S. M., & Wise, M. B. 2007, Phys. Rev. D, 75, 083502 [NASA ADS] [CrossRef] [Google Scholar]
 Adams, F. C., Bond, J. R., Freese, K., Frieman, J. A., & Olinto, A. V. 1993, Phys. Rev. D, 47, 426 [NASA ADS] [CrossRef] [Google Scholar]
 Adams, J. A., Cresswell, B., & Easther, R. 2001, Phys. Rev. D, 64, 123514 [NASA ADS] [CrossRef] [Google Scholar]
 Agarwal, N., & Bean, R. 2009, Phys. Rev. D, 79, 023503 [NASA ADS] [CrossRef] [Google Scholar]
 Alishahiha, M., Silverstein, E., & Tong, D. 2004, Phys. Rev. D, 70, 123505 [NASA ADS] [CrossRef] [Google Scholar]
 Amendola, L., Gordon, C., Wands, D., & Sasaki, M. 2002, Phys. Rev. Lett., 88, 211302 [NASA ADS] [CrossRef] [Google Scholar]
 Anber, M. M., & Sorbo, L. 2010, Phys. Rev. D, 81, 043534 [NASA ADS] [CrossRef] [Google Scholar]
 Anderson, L., Aubourg, E, Bailey, S., et al. 2014, MNRAS, 441, 24 [NASA ADS] [CrossRef] [Google Scholar]
 ArmendárizPicón, C., Damour, T., & Mukhanov, V. 1999, Phys. Lett. B, 458, 209 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Aslanyan, G., Price, L. C., Abazajian, K. N., & Easther, R. 2014, JCAP, 08, 052 [NASA ADS] [CrossRef] [Google Scholar]
 Audren, B., Lesgourgues, J., Benabed, K., & Prunet, S. 2013, JCAP, 1302, 001 [NASA ADS] [CrossRef] [Google Scholar]
 Barnaby, N., & Peloso, M. 2011, Phys. Rev. Lett., 106, 181301 [NASA ADS] [CrossRef] [Google Scholar]
 Barnaby, N., Pajer, E., & Peloso, M. 2012, Phys. Rev. D, 85, 023525 [NASA ADS] [CrossRef] [Google Scholar]
 Bartolo, N., Matarrese, S., Peloso, M., & Ricciardone, A. 2013, JCAP, 1308, 022 [NASA ADS] [CrossRef] [Google Scholar]
 Bartolo, N., Matarrese, S., & Riotto, A. 2001, Phys. Rev. D, 64, 123504 [NASA ADS] [CrossRef] [Google Scholar]
 Bartolo, N., Matarrese, S., & Riotto, A. 2004a, Phys. Rev. Lett., 93, 231301 [NASA ADS] [CrossRef] [Google Scholar]
 Babich, D., Creminelli, P., & Zaldarriaga, M. 2004, JCAP, 8, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Barnaby, N., Namba, R., & Peloso, M. 2011, JCAP, 4, 9 [Google Scholar]
 Bartolo, N., & Liddle, A. R. 2002, Phys. Rev. D, 65, 121301 [NASA ADS] [CrossRef] [Google Scholar]
 Bartolo, N., Matarrese, S., & Riotto, A. 2004b, Phys. Rev. D, 69, 043503 [NASA ADS] [CrossRef] [Google Scholar]
 Barvinsky, A. O., Kamenshchik, A. Y., & Starobinsky, A. A. 2008, JCAP, 0811, 021 [NASA ADS] [CrossRef] [Google Scholar]
 Baumann, D., Green, D., & Porto, R. A. 2015, JCAP, 1, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Bean, R., Chen, X., Peiris, H., & Xu, J. 2008, Phys. Rev. D, 77, 023527 [NASA ADS] [CrossRef] [Google Scholar]
 Becker, R. H., Fan, X., White, R. L., et al. 2001, AJ, 122, 2850 [NASA ADS] [CrossRef] [Google Scholar]
 Beltrán, M., GarcíaBellido, J., Lesgourgues, J., & Viel, M. 2005, Phys. Rev. D, 72, 103515 [NASA ADS] [CrossRef] [Google Scholar]
 Bennett, C. L., Hill, R. S., Hinshaw, G., et al. 2011, ApJS, 192, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Berera, A. 1995, Phys. Rev. Lett., 75, 3218 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Beutler, F., Blake, C., Colless, M., et al. 2011, MNRAS, 416, 3017 [NASA ADS] [CrossRef] [Google Scholar]
 Bezrukov, F., & Shaposhnikov, M. 2008, Phys. Lett. B, 659, 703 [NASA ADS] [CrossRef] [Google Scholar]
 Bezrukov, F., & Shaposhnikov, M. 2009, JHEP, 07, 089 [NASA ADS] [CrossRef] [Google Scholar]
 BICEP2 Collaboration. 2014a, ApJ, 792, 62 [NASA ADS] [CrossRef] [Google Scholar]
 BICEP2 Collaboration. 2014b, Phys. Rev. Lett., 112, 241101 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 BICEP2/Keck Array and Planck Collaborations. 2015, Phys. Rev. Lett., 114, 101301 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Binétruy, P., Kiritsis, E., Mabillard, J., Pieroni, M., & Rosset, C. 2015, J. Cosmol. Astropart. Phys., 2015, 33 [CrossRef] [Google Scholar]
 Blas, D., Lesgourgues, J., & Tram, T. 2011, JCAP, 1107, 034 [NASA ADS] [CrossRef] [Google Scholar]
 Boubekeur, L., & Lyth, D. 2005, JCAP, 0507, 010 [NASA ADS] [CrossRef] [Google Scholar]
 Boyle, L. A., Steinhardt, P. J., & Turok, N. 2004, Phys. Rev. D, 69, 127302 [NASA ADS] [CrossRef] [Google Scholar]
 Bozza, V., Giovannini, M., & Veneziano, G. 2003, JCAP, 0305, 001 [NASA ADS] [CrossRef] [Google Scholar]
 Brandenberger, R. H., Nayeri, A., Patil, S. P., & Vafa, C. 2007, Phys. Rev. Lett., 98, 231302 [NASA ADS] [CrossRef] [Google Scholar]
 Bucher, M. 2015, Int. J. Mod. Phys. D, 24, 1530004 [NASA ADS] [CrossRef] [Google Scholar]
 Bucher, M., Moodley, K., & Turok, N. 2000, Phys. Rev. D, 62, 083508 [NASA ADS] [CrossRef] [Google Scholar]
 Bucher, M., Moodley, K., & Turok, N. 2001, Phys. Rev. Lett., 87, 191301 [NASA ADS] [CrossRef] [Google Scholar]
 Buchmüller, W., Domcke, V., & Kamada, K. 2013, Phys. Lett., B726, 467 [Google Scholar]
 Burgess, C., Martineau, P., Quevedo, F., Rajesh, G., & Zhang, R. 2002, JHEP, 0203, 052 [Google Scholar]
 Burrage, C., de Rham, C., Seery, D., & Tolley, A. J. 2011, JCAP, 1, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Byrnes, C. T., & Wands, D. 2006, Phys. Rev. D, 74, 043529 [NASA ADS] [CrossRef] [Google Scholar]
 Casadio, R., Finelli, F., Kamenshchik, A., Luzzi, M., & Venturi, G. 2006, JCAP, 0604, 011 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, X. 2005a, JHEP, 8, 45 [Google Scholar]
 Chen, X. 2005b, Phys. Rev. D, 71, 063506 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, X. 2005c, Phys. Rev. D, 72, 123518 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, X., Huang, M.X., Kachru, S., & Shiu, G. 2007, JCAP, 1, 2 [Google Scholar]
 Chen, X., Easther, R., & Lim, E. A. 2008, JCAP, 0804, 010 [NASA ADS] [CrossRef] [Google Scholar]
 Chluba, J., Hamann, J., & Patil, S. P. 2015, Int. J. Mod. Phys. D, 24, 1530023 [NASA ADS] [CrossRef] [Google Scholar]
 Choe, J., Gong, J.O., & Stewart, E. D. 2004, JCAP, 0407, 012 [NASA ADS] [CrossRef] [Google Scholar]
 Cicoli, M., Burgess, C., & Quevedo, F. 2009, JCAP, 0903, 013 [NASA ADS] [CrossRef] [Google Scholar]
 Contaldi, C. R., Peloso, M., Kofman, L., & Linde, A. D. 2003, JCAP, 0307, 002 [Google Scholar]
 Copeland, E. J., Liddle, A. R., Lyth, D. H., Stewart, E. D., & Wands, D. 1994, Phys. Rev. D, 49, 6410 [NASA ADS] [CrossRef] [Google Scholar]
 Creminelli, P., D’Amico, G., Musso, M., Noreña, J., & Trincherini, E. 2011, JCAP, 2, 6 [Google Scholar]
 Danielsson, U. H. 2002, Phys. Rev. D, 66, 023511 [NASA ADS] [CrossRef] [Google Scholar]
 De Felice, A., & Tsujikawa, S. 2013, JCAP, 3, 30 [NASA ADS] [CrossRef] [Google Scholar]
 de Rham, C., & Gabadadze, G. 2010a, Phys. Rev. D, 82, 044020 [Google Scholar]
 de Rham, C., & Gabadadze, G. 2010b, Phys. Lett. B, 693, 334 [NASA ADS] [CrossRef] [Google Scholar]
 Deffayet, C., Deser, S., & EspositoFarèse, G. 2009a, Phys. Rev. D, 80, 064015 [NASA ADS] [CrossRef] [Google Scholar]
 Deffayet, C., EspositoFarèse, G., & Vikman, A. 2009b, Phys. Rev. D, 79, 084003 [NASA ADS] [CrossRef] [Google Scholar]
 Dimastrogiovanni, E., Bartolo, N., Matarrese, S., & Riotto, A. 2010, Adv. Astron., 752670 [Google Scholar]
 Dvali, G., & Tye, S. H. 1999, Phys. Lett. B, 450, 72 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Dvali, G. R., Shafi, Q., & Schaefer, R. K. 1994, Phys. Rev. Lett., 73, 1886 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Dvali, G., Shafi, Q., & Solganik, S. 2001, ArXiv eprints [arXiv:hepth/0105203] [Google Scholar]
 Easther, R., & Flauger, R. 2014, JCAP, 1402, 037 [NASA ADS] [CrossRef] [Google Scholar]
 Easther, R., & Peiris, H. 2006, JCAP, 0609, 010 [NASA ADS] [Google Scholar]
 Easther, R., & Peiris, H. V. 2012, Phys. Rev. D, 85, 103533 [NASA ADS] [CrossRef] [Google Scholar]
 Easther, R., Kinney, W. H., & Peiris, H. 2005, JCAP, 0505, 009 [NASA ADS] [Google Scholar]
 Eisenstein, D. J., Zehavi, I., Hogg, D. W., et al. 2005, ApJ, 633, 560 [NASA ADS] [CrossRef] [Google Scholar]
 Ellis, J., Nanopoulos, D. V., & Olive, K. A. 2013a, Phys. Rev. Lett., 111, 111301 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Ellis, J., Nanopoulos, D. V., & Olive, K. A. 2013b, JCAP, 1310, 009 [NASA ADS] [CrossRef] [Google Scholar]
 Enqvist, K., & Sloth, M. S. 2002, Nucl. Phys. B, 626, 395 [NASA ADS] [CrossRef] [Google Scholar]
 Enqvist, K., KurkiSuonio, H., & Valiviita, J. 2002, Phys. Rev. D, 65, 043002 [NASA ADS] [CrossRef] [Google Scholar]
 Erickcek, A. L., Hirata, C. M., & Kamionkowski, M. 2009, Phys. Rev. D, 80, 083507 [NASA ADS] [CrossRef] [Google Scholar]
 Farakos, F., Kehagias, A., & Riotto, A. 2013, Nucl. Phys. B, 876, 187 [NASA ADS] [CrossRef] [Google Scholar]
 Feroz, F., & Hobson, M. 2008, MNRAS, 384, 449 [NASA ADS] [CrossRef] [Google Scholar]
 Feroz, F., Hobson, M., & Bridges, M. 2009, MNRAS, 398, 1601 [NASA ADS] [CrossRef] [Google Scholar]
 Feroz, F., Hobson, M., Cameron, E., & Pettitt, A. 2013, ArXiv eprints [arXiv:1306.2144] [Google Scholar]
 Ferrara, S., Kallosh, R., Linde, A., & Porrati, M. 2013a, Phys. Rev. D, 88, 085038 [NASA ADS] [CrossRef] [Google Scholar]
 Ferrara, S., Kallosh, R., & Van Proeyen, A. 2013b, JHEP, 1311, 134 [NASA ADS] [CrossRef] [Google Scholar]
 Ferreira, R. Z., & Sloth, M. S. 2014, JHEP, 12, 139 [Google Scholar]
 Finelli, F., Hamann, J., Leach, S. M., & Lesgourgues, J. 2010, JCAP, 1004, 011 [NASA ADS] [CrossRef] [Google Scholar]
 Flauger, R., Hill, J. C., & Spergel, D. N. 2014a, JCAP, 1408, 039 [NASA ADS] [CrossRef] [Google Scholar]
 Flauger, R., & Pajer, E. 2011, JCAP, 1101, 017 [NASA ADS] [CrossRef] [Google Scholar]
 Flauger, R., McAllister, L., Silverstein, E., & Westphal, A. 2014b, ArXiv eprints [arXiv:1412.1814] [Google Scholar]
 Freese, K., Frieman, J. A., & Olinto, A. V. 1990, Phys. Rev. Lett., 65, 3233 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Gao, X., & Steer, D. A. 2011, JCAP, 12, 19 [Google Scholar]
 GarciaBellido, J., & Roest, D. 2014, Phys. Rev. D, 89, 103527 [NASA ADS] [CrossRef] [Google Scholar]
 GarcíaBellido, J., & Wands, D. 1996, Phys. Rev. D, 53, 5437 [NASA ADS] [CrossRef] [Google Scholar]
 GarciaBellido, J., Rabadan, R., & Zamora, F. 2002, JHEP, 0201, 036 [Google Scholar]
 Garriga, J., & Mukhanov, V. F. 1999, Phys. Lett. B, 458, 219 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Gasperini, M., & Veneziano, G. 1993, Astropart. Phys., 1, 317 [NASA ADS] [CrossRef] [Google Scholar]
 Gauthier, C., & Bucher, M. 2012, JCAP, 1210, 050 [NASA ADS] [CrossRef] [Google Scholar]
 Gerbino, M., Marchini, A., Pagano, L., et al. 2014, Phys. Rev. D, 90, 047301 [NASA ADS] [CrossRef] [Google Scholar]
 Goncharov, A., & Linde, A. D. 1984, Sov. Phys. JETP, 59, 930 [Google Scholar]
 Gong, J.O., & Stewart, E. D. 2001, Phys. Lett. B, 510, 1 [NASA ADS] [Google Scholar]
 Gorbunov, D., & Panin, A. 2011, Phys. Lett. B, 700, 157 [NASA ADS] [CrossRef] [Google Scholar]
 Gordon, C., & Lewis, A. 2003, Phys. Rev. D, 67, 123513 [NASA ADS] [CrossRef] [Google Scholar]
 Gordon, C., Wands, D., Bassett, B. A., & Maartens, R. 2001, Phys. Rev. D, 63, 023506 [NASA ADS] [CrossRef] [Google Scholar]
 Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Grin, D., Hanson, D., Holder, G. P., Doré, O., & Kamionkowski, M. 2014, Phys. Rev. D, 89, 023006 [NASA ADS] [CrossRef] [Google Scholar]
 Guth, A. H., Kaiser, D. I., & Nomura, Y. 2014, Phys. Lett. B, 733, 112 [NASA ADS] [CrossRef] [Google Scholar]
 Habib, S., Heitmann, K., Jungman, G., & MolinaParis, C. 2002, Phys. Rev. Lett., 89, 281301 [NASA ADS] [CrossRef] [Google Scholar]
 Hamann, J., Lesgourgues, J., & Valkenburg, W. 2008, JCAP, 0804, 016 [NASA ADS] [CrossRef] [Google Scholar]
 Hamimeche, S., & Lewis, A. 2008, Phys. Rev. D, 77, 103013 [NASA ADS] [CrossRef] [Google Scholar]
 Handley, W. J., Hobson, M. P., & Lasenby, A. N. 2015, MNRAS, 450, L61 [NASA ADS] [CrossRef] [Google Scholar]
 Harrison, E. R. 1970, Phys. Rev. D, 1, 2726 [NASA ADS] [CrossRef] [Google Scholar]
 Hu, W., Seljak, U., White, M. J., & Zaldarriaga, M. 1998, Phys. Rev. D, 57, 3290 [NASA ADS] [CrossRef] [Google Scholar]
 Hu, W., & White, M. J. 1997, Phys. Rev. D, 56, 596 [NASA ADS] [CrossRef] [Google Scholar]
 Ijjas, A., Steinhardt, P. J., & Loeb, A. 2013, Phys. Lett. B, 723, 261 [NASA ADS] [CrossRef] [Google Scholar]
 Ijjas, A., Steinhardt, P. J., & Loeb, A. 2014, Phys. Lett. B, 736, 142 [NASA ADS] [CrossRef] [Google Scholar]
 Jackson, M. G., & Shiu, G. 2013, Phys. Rev. D, 88, 123511 [NASA ADS] [CrossRef] [Google Scholar]
 Jackson, M. G., Wandelt, B., & Bouchet, F. 2014, Phys. Rev. D, 89, 023510 [NASA ADS] [CrossRef] [Google Scholar]
 Kachru, S., Kallosh, R., Linde, A. D., et al. 2003, JCAP, 0310, 013 [NASA ADS] [CrossRef] [Google Scholar]
 Kaiser, D. I., & Sfakianakis, E. I. 2014, Phys. Rev. Lett., 112, 011302 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Kallosh, R., & Linde, A. 2013a, JCAP, 1310, 033 [NASA ADS] [CrossRef] [Google Scholar]
 Kallosh, R., & Linde, A. 2013b, JCAP, 1306, 028 [NASA ADS] [CrossRef] [Google Scholar]
 Kallosh, R., Linde, A., & Roest, D. 2013, JHEP, 1311, 198 [Google Scholar]
 Kaloper, N., Lawrence, A., & Sorbo, L. 2011, JCAP, 1103, 023 [NASA ADS] [CrossRef] [Google Scholar]
 Kamionkowski, M., Kosowsky, A., & Stebbins, A. 1997, Phys. Rev. D, 55, 7368 [NASA ADS] [CrossRef] [Google Scholar]
 Kanno, S., Sasaki, M., & Tanaka, T. 2013, Progr. Theor. Exp. Phys., 2013, 111E01 [CrossRef] [Google Scholar]
 Kawasaki, M., & Sekiguchi, T. 2008, Prog. Theor. Phys., 120, 995 [NASA ADS] [CrossRef] [Google Scholar]
 Ketov, S. V., & Starobinsky, A. A. 2011, Phys. Rev. D, 83, 063512 [NASA ADS] [CrossRef] [Google Scholar]
 Kim, J., & Komatsu, E. 2013, Phys. Rev. D, 88, 101301 [NASA ADS] [CrossRef] [Google Scholar]
 Kim, J. E., Nilles, H. P., & Peloso, M. 2005, JCAP, 0501, 005 [NASA ADS] [CrossRef] [Google Scholar]
 Kinney, W. H. 2002, Phys. Rev. D, 66, 083508 [NASA ADS] [CrossRef] [Google Scholar]
 Kinney, W. H., Kolb, E. W., Melchiorri, A., & Riotto, A. 2006, Phys. Rev. D, 74, 023502 [NASA ADS] [CrossRef] [Google Scholar]
 Kobayashi, T., & Takahashi, F. 2011, JCAP, 1101, 026 [NASA ADS] [CrossRef] [Google Scholar]
 Kobayashi, T., Yamaguchi, M., & Yokoyama, J. 2010, Phys. Rev. Lett., 105, 231302 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Kobayashi, T., Yamaguchi, M., & Yokoyama, J. 2011a, Progr. Theor. Phys., 126, 511 [Google Scholar]
 Kobayashi, T., Yamaguchi, M., & Yokoyama, J. 2011b, Phys. Rev. D, 83, 103524 [NASA ADS] [CrossRef] [Google Scholar]
 Kosowsky, A., & Turner, M. S. 1995, Phys. Rev. D, 52, 1739 [NASA ADS] [CrossRef] [Google Scholar]
 Langlois, D. 1999, Phys. Rev. D, 59, 123512 [NASA ADS] [CrossRef] [Google Scholar]
 Langlois, D., & Riazuelo, A. 2000, Phys. Rev. D, 62, 043504 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Leach, S. M., Liddle, A. R., Martin, J., & Schwarz, D. J. 2002, Phys. Rev. D, 66, 023515 [NASA ADS] [CrossRef] [Google Scholar]
 Lesgourgues, J. 2011, ArXiv eprints [arXiv:1104.2932] [Google Scholar]
 Lesgourgues, J., & Valkenburg, W. 2007, Phys. Rev. D, 75, 123519 [NASA ADS] [CrossRef] [Google Scholar]
 Lesgourgues, J., Polarski, D., Prunet, S., & Starobinsky, A. A. 2000, Astron. Astrophys., 359, 414 [Google Scholar]
 Lesgourgues, J., Starobinsky, A. A., & Valkenburg, W. 2008, JCAP, 0801, 010 [NASA ADS] [CrossRef] [Google Scholar]
 Lewis, A., & Bridle, S. 2002, Phys. Rev. D, 66, 103511 [NASA ADS] [CrossRef] [Google Scholar]
 Libanov, M., & Rubakov, V. 2010, JCAP, 1011, 045 [NASA ADS] [CrossRef] [Google Scholar]
 Liddle, A. R., & Cortês, M. 2013, Phys. Rev. Lett., 111, 111302 [NASA ADS] [CrossRef] [Google Scholar]
 Linde, A. 2014, ArXiv eprints [arXiv:1402.0526] [Google Scholar]
 Linde, A. D. 1983, Phys. Lett. B, 129, 177 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Linde, A. D. 1994, Phys. Rev. D, 49, 748 [NASA ADS] [CrossRef] [Google Scholar]
 Linde, A. D., & Mukhanov, V. 2006, JCAP, 0604, 009 [NASA ADS] [CrossRef] [Google Scholar]
 Linde, A. D., & Mukhanov, V. F. 1997, Phys. Rev. D, 56, 535 [NASA ADS] [CrossRef] [Google Scholar]
 Liddle, A. R., & Leach, S. M. 2003, Phys. Rev. D, 68, 103503 [NASA ADS] [CrossRef] [Google Scholar]
 Linde, A., Mooij, S., & Pajer, E. 2013, Phys. Rev. D, 87, 103506 [NASA ADS] [CrossRef] [Google Scholar]
 Lorenz, L., Martin, J., & Ringeval, C. 2008, Phys. Rev. D, 78, 083513 [NASA ADS] [CrossRef] [Google Scholar]
 Lyth, D. H. 2013, JCAP, 1308, 007 [NASA ADS] [CrossRef] [Google Scholar]
 Lyth, D. H., & Wands, D. 2002, Phys. Lett. B, 524, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Lyth, D. H., Ungarelli, C., & Wands, D. 2003, Phys. Rev. D, 67, 023503 [NASA ADS] [CrossRef] [Google Scholar]
 Ma, C.P., & Bertschinger, E. 1995, ApJ, 455, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Ma, Y.Z., Efstathiou, G., & Challinor, A. 2011, Phys. Rev. D, 83, 083005 [NASA ADS] [CrossRef] [Google Scholar]
 Ma, Y.Z., Huang, Q.G., & Zhang, X. 2013, Phys. Rev. D, 87, 103516 [NASA ADS] [CrossRef] [Google Scholar]
 Maleknejad, A., SheikhJabbari, M., & Soda, J. 2013, Phys. Rept., 528, 161 [Google Scholar]
 Martin, J., & Brandenberger, R. H. 2001, Phys. Rev. D, 63, 123501 [NASA ADS] [CrossRef] [Google Scholar]
 Martin, J., & Ringeval, C. 2010, Phys. Rev. D, 82, 023511 [NASA ADS] [CrossRef] [Google Scholar]
 Martin, J., & Schwarz, D. J. 2003, Phys. Rev. D, 67, 083512 [NASA ADS] [CrossRef] [Google Scholar]
 Martin, J., Ringeval, C., & Vennin, V. 2014, Phys. Dark Univ., 56, 75 [CrossRef] [Google Scholar]
 McAllister, L., Silverstein, E., & Westphal, A. 2010, Phys. Rev. D, 82, 046003 [NASA ADS] [CrossRef] [Google Scholar]
 McAllister, L., Silverstein, E., Westphal, A., & Wrase, T. 2014, JHEP, 1409, 123 [Google Scholar]
 Meerburg, P. D. 2014, Phys. Rev. D, 90, 063529 [NASA ADS] [CrossRef] [Google Scholar]
 Meerburg, P. D., & Pajer, E. 2013, JCAP, 2, 17 [Google Scholar]
 Meerburg, P. D., & Spergel, D. N. 2014, Phys. Rev. D, 89, 063537 [NASA ADS] [CrossRef] [Google Scholar]
 Meerburg, P. D., Spergel, D. N., & Wandelt, B. D. 2014a, ArXiv eprints [arXiv:1406.0548] [Google Scholar]
 Meerburg, P. D., Spergel, D. N., & Wandelt, B. D. 2014b, Phys. Rev. D, 89, 063536 [NASA ADS] [CrossRef] [Google Scholar]
 Miranda, V., & Hu, W. 2014, Phys. Rev. D, 89, 083529 [NASA ADS] [CrossRef] [Google Scholar]
 Mizuno, S., & Koyama, K. 2010, Phys. Rev. D, 82, 103518 [NASA ADS] [CrossRef] [Google Scholar]
 Mollerach, S. 1990, Phys. Rev. D, 42, 313 [NASA ADS] [CrossRef] [Google Scholar]
 Moroi, T., & Takahashi, T. 2001, Phys. Lett. B, 522, 215 [NASA ADS] [CrossRef] [Google Scholar]
 Mortonson, M. J., & Hu, W. 2008, ApJ, 672, 737 [NASA ADS] [CrossRef] [Google Scholar]
 Mortonson, M. J., Dvorkin, C., Peiris, H. V., & Hu, W. 2009, Phys. Rev. D, 79, 103519 [NASA ADS] [CrossRef] [Google Scholar]
 Mortonson, M. J., Peiris, H. V., & Easther, R. 2011, Phys. Rev. D, 83, 043505 [NASA ADS] [CrossRef] [Google Scholar]
 Mortonson, M. J., & Seljak, U. 2014, JCAP, 1410, 035 [NASA ADS] [CrossRef] [Google Scholar]
 Moss, A., Scott, D., Zibin, J. P., & Battye, R. 2011, Phys. Rev. D, 84, 023014 [NASA ADS] [CrossRef] [Google Scholar]
 Mukhanov, V. 2013, Eur. Phys. J. C, 73, 2486 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mukhanov, V. F., & Chibisov, G. 1981, JETP Lett., 33, 532 [NASA ADS] [Google Scholar]
 Naruko, A., Komatsu, E., & Yamaguchi, M. 2015, J. Cosmol. Astropart. Phys., 2015, 45 [CrossRef] [Google Scholar]
 Nicolis, A., Rattazzi, R., & Trincherini, E. 2009, Phys. Rev. D, 79, 064036 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Norena, J., Wagner, C., Verde, L., Peiris, H. V., & Easther, R. 2012, Phys. Rev. D, 86, 023505 [NASA ADS] [CrossRef] [Google Scholar]
 Ohashi, J., & Tsujikawa, S. 2012, JCAP, 10, 35 [Google Scholar]
 Ohashi, J., Soda, J., & Tsujikawa, S. 2013, Phys. Rev. D, 87, 083520 [NASA ADS] [CrossRef] [Google Scholar]
 Okamoto, T., & Hu, W. 2003, Phys. Rev. D, 67, 083002 [NASA ADS] [CrossRef] [Google Scholar]
 Page, L., Hinshaw, G., Komatsu, E., et al. 2007, ApJS, 170, 335 [NASA ADS] [CrossRef] [Google Scholar]
 Pajer, E., & Peloso, M. 2013, Class. Quant. Grav., 30, 214002 [NASA ADS] [CrossRef] [Google Scholar]
 PalanqueDelabrouille, N., Yèche, C., Lesgourgues, J., et al. 2015, J. Cosmol. Astropart. Phys., 2015, 45 [NASA ADS] [CrossRef] [Google Scholar]
 Peebles, P., & Yu, J. 1970, ApJ, 162, 815 [NASA ADS] [CrossRef] [Google Scholar]
 Peiris, H., & Easther, R. 2006a, JCAP, 0607, 002 [NASA ADS] [CrossRef] [Google Scholar]
 Peiris, H., & Easther, R. 2006b, JCAP, 0610, 017 [NASA ADS] [CrossRef] [Google Scholar]
 Peiris, H. V., & Easther, R. 2008, JCAP, 0807, 024 [NASA ADS] [CrossRef] [Google Scholar]
 Peiris, H. V., Komatsu, E., Verde, L., et al. 2003, ApJS, 148, 213 [NASA ADS] [CrossRef] [Google Scholar]
 Peiris, H., Baumann, D., Friedman, B., & Cooray, A. 2007, Phys. Rev. D, 76, 103517 [NASA ADS] [CrossRef] [Google Scholar]
 Peiris, H., Easther, R., & Flauger, R. 2013, JCAP, 1309, 018 [NASA ADS] [CrossRef] [Google Scholar]
 Perotto, L., Lesgourgues, J., Hannestad, S., Tu, H., & Wong, Y. Y. 2006, JCAP, 0610, 013 [Google Scholar]
 Planck Collaboration XV. 2014, A&A, 571, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVI. 2014, A&A, 571, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVII. 2014, A&A, 571, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXII. 2014, A&A, 571, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXIII. 2014, A&A, 571, A23 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Planck Collaboration XXIV. 2014, A&A, 571, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXV. 2014, A&A, 571, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVI. 2014, A&A, 571, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVII. 2014, A&A, 571, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXX. 2014, A&A, 571, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXXI. 2014, A&A, 571, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration I. 2016, A&A, 594, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration II. 2016, A&A, 594, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration III. 2016, A&A, 594, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IV. 2016, A&A, 594, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration V. 2016, A&A, 594, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VI. 2016, A&A, 594, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VII. 2016, A&A, 594, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VIII. 2016, A&A, 594, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IX. 2016, A&A, 594, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration X. 2016, A&A, 594, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XI. 2016, A&A, 594, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XII. 2016, A&A, 594, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIV. 2016, A&A, 594, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XV. 2016, A&A, 594, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVI. 2016, A&A, 594, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVII. 2016, A&A, 594, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVIII. 2016, A&A, 594, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIX. 2016, A&A, 594, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XX. 2016, A&A, 594, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXI. 2016, A&A, 594, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXII. 2016, A&A, 594, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXIII. 2016, A&A, 594, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXIV. 2016, A&A, 594, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXV. 2016, A&A, 594, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVI. 2016, A&A, 594, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVII. 2016, A&A, 594, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVIII. 2016, A&A, 594, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XIX. 2015, A&A, 576, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XX. 2015, A&A, 576, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXI. 2015, A&A, 576, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXII. 2015, A&A, 576, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXX. 2016, A&A, 586, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Polarski, D., & Starobinsky, A. A. 1996, Class. Quant. Grav., 13, 377 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Powell, B. A., & Kinney, W. H. 2007, JCAP, 0708, 006 [NASA ADS] [CrossRef] [Google Scholar]
 Powell, B. A., Tzirakis, K., & Kinney, W. H. 2009, JCAP, 4, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Powell, M. J. D. 2009, The BoBYQA algorithm for bound constrained optimization without derivatives [Google Scholar]
 Regan, D., Anderson, G. J., Hull, M., & Seery, D. 2015, J. Cosmol. Astropart. Phys., 2015, 15 [CrossRef] [Google Scholar]
 Ross, A. J., Samushia, L., Howlett, C., et al. 2015, MNRAS., 449, 835 [NASA ADS] [CrossRef] [Google Scholar]
 Rossi, G., Yeche, C., PalanqueDelabrouille, N., & Lesgourgues, J. 2015, Phys. Rev. D, 92, 063505 [NASA ADS] [CrossRef] [Google Scholar]
 Sasaki, M., Valiviita, J., & Wands, D. 2006, Phys. Rev. D, 74, 103003 [NASA ADS] [CrossRef] [Google Scholar]
 Savelainen, M., Valiviita, J., Walia, P., Rusak, S., & KurkiSuonio, H. 2013, Phys. Rev. D, 88, 063010 [NASA ADS] [CrossRef] [Google Scholar]
 Schmidt, F., & Hui, L. 2013, Phys. Rev. Lett., 110, 011301 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Seljak, U., & Zaldarriaga, M. 1997, Phys. Rev. Lett., 78, 2054 [NASA ADS] [CrossRef] [Google Scholar]
 Senatore, L., Smith, K. M., & Zaldarriaga, M. 2010, JCAP, 1, 28 [NASA ADS] [CrossRef] [Google Scholar]
 Silverstein, E., & Tong, D. 2004, Phys. Rev. D, 70, 103505 [NASA ADS] [CrossRef] [Google Scholar]