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



Article Number  A14  
Number of page(s)  31  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201525814  
Published online  20 September 2016 
Planck 2015 results
XIV. Dark energy and modified gravity
^{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, 68 Melrose Road, Muizenberg, Cape Town, South Africa
^{4} Agenzia Spaziale Italiana Science Data Center, via del Politecnico snc, 00133 Roma, Italy
^{5} Aix Marseille Université, CNRS, LAM (Laboratoire d’Astrophysique de Marseille) UMR 7326, 13388 Marseille, France
^{6} Aix Marseille Université, Centre de Physique Théorique, 163 Avenue de Luminy, 13288 Marseille, France
^{7} Astrophysics Group, Cavendish Laboratory, University of Cambridge, J J Thomson Avenue, Cambridge CB3 0HE, UK
^{8} Astrophysics & Cosmology Research Unit, School of Mathematics, Statistics & Computer Science, University of KwaZuluNatal, Westville Campus, Private Bag X54001, 4000 Durban, South Africa
^{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, Ireland
^{13} California Institute of Technology, Pasadena, CA 91195 California, 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 91195 California, USA
^{17} Consejo Superior de Investigaciones Científicas (CSIC), 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} Departamento de Astrofísica, Universidad de La Laguna (ULL), 38206 La Laguna, Tenerife, Spain
^{22} Departamento de Física, Universidad de Oviedo, Avda. Calvo Sotelo s/n, 33003 Oviedo, Spain
^{23} Department of Astronomy and Astrophysics, University of Toronto, 50 Saint George Street, Toronto, Ontario, Canada
^{24} Department of Astrophysics/IMAPP, Radboud University Nijmegen, PO Box 9010, 6500 GL Nijmegen, The Netherlands
^{25} Department of Physics & Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, British Columbia, Canada
^{26} Department of Physics and Astronomy, Dana and David Dornsife College of Letter, Arts and Sciences, University of Southern California, Los Angeles, CA 90089, USA
^{27} Department of Physics and Astronomy, University College London, London WC1E 6BT, UK
^{28} Department of Physics and Astronomy, University of Sussex, Brighton BN1 9QH, UK
^{29} Department of Physics, Florida State University, Keen Physics Building, 77 Chieftan Way, Tallahassee, Florida, USA
^{30} Department of Physics, Gustaf Hällströmin katu 2a, University of Helsinki, 00100 Helsinki, Finland
^{31} Department of Physics, Princeton University, NJ 08544 Princeton, New Jersey, USA
^{32} Department of Physics, University of California, Berkeley, CA 91195 California, USA
^{33} Department of Physics, University of California, Santa Barbara, CA 91195 California, USA
^{34} Department of Physics, University of Illinois at UrbanaChampaign, 1110 West Green Street, Urbana, Illinois, USA
^{35} Dipartimento di Fisica e Astronomia G. Galilei, Università degli Studi di Padova, via Marzolo 8, 35131 Padova, Italy
^{36} Dipartimento di Fisica e Scienze della Terra, Università di Ferrara, via Saragat 1, 44122 Ferrara, Italy
^{37} Dipartimento di Fisica, Università La Sapienza, P.le A. Moro 2, Roma, Italy
^{38} Dipartimento di Fisica, Università degli Studi di Milano, via Celoria 16, Milano, Italy
^{39} Dipartimento di Fisica, Università degli Studi di Trieste, via A. Valerio 2, Trieste, Italy
^{40} Dipartimento di Matematica, Università di Roma Tor Vergata, via della Ricerca Scientifica 1, Roma, Italy
^{41} Discovery Center, Niels Bohr Institute, Blegdamsvej 17, Copenhagen, Denmark
^{42} Discovery Center, Niels Bohr Institute, Copenhagen University, Blegdamsvej 17, Copenhagen, Denmark
^{43} European Space Agency, ESAC, Planck Science Office, Camino bajo del Castillo, s/n, Urbanización Villafranca del Castillo, Villanueva de la Cañada, 28692 Madrid, Spain
^{44} European Space Agency, ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands
^{45} Gran Sasso Science Institute, INFN, viale F. Crispi 7, 67100 L’ Aquila, Italy
^{46} HGSFP and University of Heidelberg, Theoretical Physics Department, Philosophenweg 16, 69120 Heidelberg, Germany
^{47} Helsinki Institute of Physics, Gustaf Hällströmin katu 2, University of Helsinki, 00100 Helsinki, Finland
^{48} INAF–Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, Padova, Italy
^{49} INAF–Osservatorio Astronomico di Roma, via di Frascati 33, Monte Porzio Catone, Italy
^{50} INAF–Osservatorio Astronomico di Trieste, via G.B. Tiepolo 11, Trieste, Italy
^{51} INAF/IASF Bologna, via Gobetti 101, Bologna, Italy
^{52} INAF/IASF Milano, via E. Bassini 15, Milano, Italy
^{53} INFN, Sezione di Bologna, via Irnerio 46, 40126 Bologna, Italy
^{54} INFN, Sezione di Roma 1, Università di Roma Sapienza, Piazzale Aldo Moro 2, 00185 Roma, Italy
^{55} INFN, Sezione di Roma 2, Università di Roma Tor Vergata, via della Ricerca Scientifica 1, Roma, Italy
^{56} INFN/National Institute for Nuclear Physics, via Valerio 2, 34127 Trieste, Italy
^{57} IPAG: Institut de Planétologie et d’Astrophysique de Grenoble, Université Grenoble Alpes, IPAG; CNRS, IPAG, 38000 Grenoble, France
^{58} IUCAA, Post Bag 4, Ganeshkhind, Pune University Campus, 411 007 Pune, India
^{59} Imperial College London, Astrophysics group, Blackett Laboratory, Prince Consort Road, London, SW7 2AZ, UK
^{60} Infrared Processing and Analysis Center, California Institute of Technology, Pasadena, CA 91125, USA
^{61} Institut Néel, CNRS, Université Joseph Fourier Grenoble I, 25 rue des Martyrs, Grenoble, France
^{62} Institut Universitaire de France, 103 Bd SaintMichel, 75005 Paris, France
^{63} Institut d’Astrophysique Spatiale, CNRS, Univ. ParisSud, Université ParisSaclay, Bât. 121, 91405 Orsay Cedex, France
^{64} Institut d’Astrophysique de Paris, CNRS (UMR7095), 98 bis Boulevard Arago, 75014 Paris, France
^{65} Institut für Theoretische Teilchenphysik und Kosmologie, RWTH Aachen University, 52056 Aachen, Germany
^{66} Institute for Space Sciences, 077125 BucharestMagurale, Romania
^{67} Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
^{68} Institute of Theoretical Astrophysics, University of Oslo, Blindern, 0316 Oslo, Norway
^{69} Instituto de Astrofísica de Canarias, C/Vía Láctea s/n, La Laguna, Tenerife, Spain
^{70} Instituto de Física de Cantabria (CSICUniversidad de Cantabria), Avda. de los Castros s/n, Santander, Spain
^{71} Istituto Nazionale di Fisica Nucleare, Sezione di Padova, via Marzolo 8, 35131 Padova, Italy
^{72} Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, California, USA
^{73} Jodrell Bank Centre for Astrophysics, Alan Turing Building, School of Physics and Astronomy, The University of Manchester, Oxford Road, Manchester, M13 9PL, UK
^{74} Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA
^{75} Kavli Institute for Cosmology Cambridge, Madingley Road, Cambridge, CB3 0HA, UK
^{76} Kazan Federal University, 18 Kremlyovskaya St., 420008 Kazan, Russia
^{77} LAL, Université ParisSud, CNRS/IN2P3, Orsay, France
^{78} LERMA, CNRS, Observatoire de Paris, 61 Avenue de l’Observatoire, 75000 Paris, France
^{79} Laboratoire AIM, IRFU/Service d’Astrophysique – CEA/DSM – CNRS – Université Paris Diderot, Bât. 709, CEASaclay, 91191 GifsurYvette Cedex, France
^{80} Laboratoire Traitement et Communication de l’Information, CNRS (UMR 5141) and Télécom ParisTech, 46 rue Barrault, 75634 Paris Cedex 13, France
^{81} Laboratoire de Physique Subatomique et Cosmologie, Université GrenobleAlpes, CNRS/IN2P3, 53 rue des Martyrs, 38026 Grenoble Cedex, France
^{82} Laboratoire de Physique Théorique, Université ParisSud 11 & CNRS, Bâtiment 210, 91405 Orsay, France
^{83} Lawrence Berkeley National Laboratory, Berkeley, CA 94720 California, USA
^{84} Lebedev Physical Institute of the Russian Academy of Sciences, Astro Space Centre, 84/32 Profsoyuznaya st., Moscow, GSP7, 117997, Russia
^{85} MaxPlanckInstitut für Astrophysik, KarlSchwarzschildStr. 1, 85741 Garching, Germany
^{86} McGill Physics, Ernest Rutherford Physics Building, McGill University, 3600 rue University, Montréal, QC, H3A 2T8, Canada
^{87} National University of Ireland, Department of Experimental Physics, Maynooth, Co. Kildare, Ireland
^{88} Nicolaus Copernicus Astronomical Center, Bartycka 18, 00716 Warsaw, Poland
^{89} Niels Bohr Institute, Blegdamsvej 17, Copenhagen, Denmark
^{90} Niels Bohr Institute, Copenhagen University, Blegdamsvej 17, 2100 Copenhagen, Denmark
^{91} Nordita (Nordic Institute for Theoretical Physics), Roslagstullsbacken 23, 106 91 Stockholm, Sweden
^{92} Optical Science Laboratory, University College London, Gower Street, WC1E 6 BT London, UK
^{93} SISSA, Astrophysics Sector, via Bonomea 265, 34136, Trieste, Italy
^{94} SMARTEST Research Centre, Università degli Studi eCampus, via Isimbardi 10, 22060 Novedrate (CO), Italy
^{95} School of Physics and Astronomy, Cardiff University, Queens Buildings, The Parade, Cardiff, CF24 3AA, UK
^{96} School of Physics and Astronomy, University of Nottingham, Nottingham NG7 2RD, UK
^{97} Sorbonne UniversitéUPMC, UMR7095, Institut d’Astrophysique de Paris, 98 bis Boulevard Arago, 75014 Paris, France
^{98} Space Research Institute (IKI), Russian Academy of Sciences, Profsoyuznaya Str, 84/32, 117997 Moscow, Russia
^{99} Space Sciences Laboratory, University of California, Berkeley, California, USA
^{100} Special Astrophysical Observatory, Russian Academy of Sciences, Nizhnij Arkhyz, Zelenchukskiy region, 369167 KarachaiCherkessian Republic, Russia
^{101} Stanford University, Dept of Physics, Varian Physics Bldg, 382 via Pueblo Mall, Stanford, California, USA
^{102} SubDepartment of Astrophysics, University of Oxford, Keble Road, Oxford OX1 3RH, UK
^{103} The Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, AlbaNova, 106 91 Stockholm, Sweden
^{104} Theory Division, PHTH, CERN, 1211, Geneva 23, Switzerland
^{105} UPMC Univ. Paris 06, UMR7095, 98 bis Boulevard Arago, 75014 Paris, France
^{106} Universität Heidelberg, Institut für Theoretische Astrophysik, Philosophenweg 12, 69120 Heidelberg, Germany
^{107} Université de Toulouse, UPSOMP, IRAP, 31028 Toulouse Cedex 4, France
^{108} Universities Space Research Association, Stratospheric Observatory for Infrared Astronomy, MS 23211, Moffett Field, CA 94035, USA
^{109} University of Granada, Departamento de Física Teórica y del Cosmos, Facultad de Ciencias, 18010 Granada, Spain
^{110} University of Granada, Instituto Carlos I de Física Teórica y Computacional, 18010 Granada, Spain
^{111} University of Heidelberg, Institute for Theoretical Physics, Philosophenweg 16, 69120 Heidelberg, Germany
^{112} Warsaw University Observatory, Aleje Ujazdowskie 4, 00478 Warszawa, Poland
^{⋆}
Corresponding author: Valeria Pettorino, email: v.pettorino@thphys.uniheidelberg.de
Received: 5 February 2015
Accepted: 19 February 2016
We study the implications of Planck data for models of dark energy (DE) and modified gravity (MG) beyond the standard cosmological constant scenario. We start with cases where the DE only directly affects the background evolution, considering Taylor expansions of the equation of state w(a), as well as principal component analysis and parameterizations related to the potential of a minimally coupled DE scalar field. When estimating the density of DE at early times, we significantly improve present constraints and find that it has to be below ~2% (at 95% confidence) of the critical density, even when forced to play a role for z < 50 only. We then move to general parameterizations of the DE or MG perturbations that encompass both effective field theories and the phenomenology of gravitational potentials in MG models. Lastly, we test a range of specific models, such as kessence, f(R) theories, and coupled DE. In addition to the latest Planck data, for our main analyses, we use background constraints from baryonic acoustic oscillations, typeIa supernovae, and local measurements of the Hubble constant. We further show the impact of measurements of the cosmological perturbations, such as redshiftspace distortions and weak gravitational lensing. These additional probes are important tools for testing MG models and for breaking degeneracies that are still present in the combination of Planck and background data sets. All results that include only background parameterizations (expansion of the equation of state, early DE, general potentials in minimallycoupled scalar fields or principal component analysis) are in agreement with ΛCDM. When testing models that also change perturbations (even when the background is fixed to ΛCDM), some tensions appear in a few scenarios: the maximum one found is ~2σ for Planck TT+lowP when parameterizing observables related to the gravitational potentials with a chosen time dependence; the tension increases to, at most, 3σ when external data sets are included. It however disappears when including CMB lensing.
Key words: dark energy / cosmic background radiation / cosmology: theory / gravitation
© ESO, 2016
1. Introduction
The cosmic microwave background (CMB) is a key probe of our cosmological model (Planck Collaboration XIII 2016), providing information on the primordial Universe and its physics, including inflationary models (Planck Collaboration XX 2016) and constraints on primordial nonGaussianities (Planck Collaboration XVII 2016). In this paper we use the 2015 data release from Planck^{1} (Planck Collaboration I 2016) to perform a systematic analysis of a large set of dark energy and modified gravity theories.
Observations have long shown that only a small fraction of the total energy density in the Universe (around 5%) is in the form of baryonic matter, with the dark matter needed for structure formation accounting for about another 26%. In one scenario, the dominant component, generically referred to as dark energy (DE), brings the total close to the critical density and is responsible for the recent phase of accelerated expansion. In another scenario, the accelerated expansion arises, partly or fully, owing to a modification of gravity on cosmological scales. Elucidating the nature of this DE and testing General Relativity (GR) on cosmological scales are major challenges for contemporary cosmology, both on the theoretical and experimental sides (e.g., LSST Science Collaboration et al. 2009; Amendola et al. 2013; Clifton et al. 2012; Joyce et al. 2015; Huterer et al. 2015).
In preparation for future experimental investigations of DE and MG, it is important to determine what we already know about these models at different epochs in redshift and different length scales. CMB anisotropies fix the cosmology at early times, while additional cosmological data sets further constrain on how DE or MG evolve at lower redshifts. The aim of this paper is to investigate models for dark energy and modified gravity using Planck data in combination with other data sets.
The simplest model for DE is a cosmological constant, Λ, first introduced by Einstein (1917) to keep the Universe static, but soon dismissed when the Universe was found to be expanding (Lemaître 1927; Hubble 1929). This constant has been reintroduced several times over the years in attempts to explain several astrophysical phenomena including, most recently, the flat spatial geometry implied by the CMB and supernova observations of a recent phase of accelerated expansion (Riess et al. 1998; Perlmutter et al. 1999). A cosmological constant is described by a single parameter, the inclusion of which brings the model (ΛCDM) into excellent agreement with the data. ΛCDM still represents a good fit to a wide range of observations, more than 20 years after it was introduced. Nonetheless, theoretical estimates for the vacuum density are many orders of magnitude larger than its observed value. In addition, Ω_{Λ} and Ω_{m} are of the same order of magnitude at present only, which marks our epoch as a special time in the evolution of the Universe (the “coincidence problem”). This lack of a clear theoretical understanding has motivated the development of a wide variety of alternative models. Those models which are close to ΛCDM are in broad agreement with current constraints on the background cosmology, but the perturbations may still evolve differently, and hence it is important to test their predictions against CMB data.
We had to face at least three difficulties within this paper. First, there appears to be a vast array of possibilities in the literature and no agreement yet in the scientific community on a comprehensive framework for discussing the landscape of models. A second complication is that robust constraints on DE come from a combination of different data sets working in concert. Hence we have to be careful in the choice of the data sets so that we do not find apparent hints for nonstandard models that are, in fact, due to systematic errors. A third area of concern is the fact that numerical codes available at present for DE and MG are not as well tested in these scenarios as for ΛCDM, especially given the accuracy reached by the data. Furthermore, in some cases, we need to rely on stability routines that deserve further investigation to assure that they are not excluding more models than required.
To navigate the range of modelling possibilities, we adopt the following threepart approach.

1.
Background parameterizations. Here we consider onlyparameterizations of backgroundlevel cosmologicalquantities. Perturbations are always included, but theirevolution depends only on the background. This set includes models involving expansions, parameterizations,or principal component analyses of the equation of statew ≡ p/ρ of a DE fluid with pressure p and energy density ρ. Early DE also belongs to this class.

2.
Perturbation parameterizations. Here the perturbations themselves are parameterized or modified explicitly, not only as a consequence of a change in background quantities. There are two main branches we consider: firstly, effective field theory for DE (EFT, e.g. Gubitosi et al. 2013), which has a clear theoretical motivation, since it includes all theories derived when accounting for all symmetry operators in the Lagrangian, written in unitary gauge, i.e. only in terms of metric perturbations. This is a very general classification that has the advantage of providing a broad overview of (at least) all universally coupled DE models. However, a clear disadvantage is that the number of free parameters is large and the constraints are consequently weak. Moreover, with the currently available numerical codes, one needs to rely on stability routines, which are not fully tested and may discard more models than necessary. As a complementary approach, we include a more phenomenological class of models, obtained by directly parameterizing two independent functions of the gravitational potentials. This approach can, in principle, probe all degrees of freedom at the background and perturbation level (e.g. Kunz 2012) and is easier to handle in numerical codes. While the connection to physical models is less obvious here than in EFT, this approach allows us to gain a more intuitive understanding of the general constraining power of the data.

3.
Examples of particular models. Here we focus on a selection of theories that have already been discussed in the literature and are better understood theoretically; these can partly be considered as applications of previous cases for which the CMB constraints are more informative, because there is less freedom in any particular theory than in a more general one.
The CMB is the cleanest probe of large scales, which are of particular interest for modifications to gravity. We will investigate the constraints coming from Planck data in combination with other data sets, addressing strengths and potential weaknesses of different analyses. Before describing in detail the models and data sets that correspond to our requirements, in Sect. 2 we first address the main question that motivates our paper, discussing why CMB is relevant for DE. We then present the specific model parameterizations in Sect. 3. The choice of data sets is discussed in detail in Sect. 4 before we present results in Sect. 5 and discuss conclusions in Sect. 6.
2. Why is the CMB relevant for dark energy?
The CMB anisotropies are largely generated at the lastscattering epoch, and hence can be used to pin down the theory at early times. In fact many forecasts of future DE or MG experiments are for new data plus constraints from Planck. However, there are also several effects that DE and MG models can have on the CMB, some of which are to:

1.
change the expansion history and hence distance to the lastscattering surface, with a shift in the peaks, sometimes referred toas a geometrical projection effect (Hu &White 1996);

2.
cause the decay of gravitational potentials at late times, affecting the lowmultipole CMB anisotropies through the integrated SachsWolfe (ISW) effect (Sachs & Wolfe 1967; Kofman & Starobinskii 1985);

3.
enhance the crosscorrelation between the CMB and largescale structure, through the ISW effect (Giannantonio et al. 2008);

4.
change the lensing potential, through additional DE perturbations or modifications of GR (Acquaviva & Baccigalupi 2006; Carbone et al. 2013);

5.
change the growth of structure (Peebles 1984; Barrow & Saich 1993) leading to a mismatch between the CMBinferred amplitude of the fluctuations A_{s} and latetime measurements of σ_{8} (Kunz et al. 2004; Baldi & Pettorino 2011);

6.
impact small scales, modifying the damping tail in ${\mathit{C}}_{\mathit{\ell}}^{\mathit{TT}}$, giving a measurement of the abundance of DE at different redshifts (Calabrese et al. 2011; Reichardt et al. 2012);

7.
affect the ratio between odd and even peaks if modifications of gravity treat baryons and cold dark matter differently (Amendola et al. 2012);

8.
modify the lensing Bmode contribution, through changes in the lensing potential (Amendola et al. 2014);

9.
modify the primordial Bmode amplitude and scale dependence, by changing the sound speed of gravitational waves (Antolini et al. 2013; Amendola et al. 2014; Raveri et al. 2015).
In this paper we restrict our analysis to scalar perturbations. The dominant effects on the temperature power spectrum are due to lensing and the ISW effect, as can be seen in Fig. 1, which shows typical power spectra of temperature anisotropies and lensing potential for modified gravity models. Different curves correspond to different choices of the μ and η functions, which change the relation between the metric potentials and the sources, as well as introducing a gravitational slip; we will define these functions in Sect. 3.2.2, Eqs. (4) and (6), respectively. Spectra are obtained using a scaleindependent evolution for both μ and η. The two parameters in the figure then determine the change in amplitude of μ and η with respect to the ΛCDM case, in which E_{11} = E_{22} = 0 and μ = η = 1.
Fig. 1 Typical effects of modified gravity on theoretical CMB temperature (top panel) and lensing potential (bottom panel) power spectra. An increase (or decrease) of E_{22} with respect to zero introduces a gravitational slip, higher at present, when Ω_{de} is higher (see Eqs. (4) and (6)); this in turns changes the Weyl potential and leads to a higher (or lower) lensing potential. On the other hand, whenever E_{11} and E_{22} are different from zero (quite independently of their sign) μ and η change in time: as the dynamics in the gravitational potential is increased, this leads to an enhancement in the ISW effect. Note also that even when the temperature spectrum is very close to ΛCDM (as for E_{11} = E_{22} = 0.5) the lensing potential is still different with respect to ΛCDM, shown in black. 
3. Models and parameterizations
We now provide an overview of the models addressed in this paper. Details on the specific parameterizations will be discussed in Sect. 5, where we also present the results for each specific method.
We start by noticing that one can generally follow two different approaches: (1) given a theoretical set up, one can specify the action (or Lagrangian) of the theory and derive background and perturbation equations in that framework; or (2) more phenomenologically, one can construct functions that map closely onto cosmological observables, probing the geometry of spacetime and the growth of perturbations. Assuming spatial flatness for simplicity, the geometry is given by the expansion rate H and perturbations to the metric. If we consider only scalartype components the metric perturbations can be written in terms of the gravitational potentials Φ and Ψ (or equivalently by any two independent combinations of these potentials). Cosmological observations thus constrain one “background” function of time H(a) and two “perturbation” functions of scale and time Φ(k,a) and Ψ(k,a) (e.g., Kunz 2012). These functions fix the metric, and thus the Einstein tensor G_{μν}. Einstein’s equations link this tensor to the energymomentum tensor T_{μν}, which in turn can be related to DE or MG properties.
Throughout this paper we will adopt the metric given by the line element $\mathrm{d}{\mathit{s}}^{\mathrm{2}}\mathrm{=}{\mathit{a}}^{\mathrm{2}}\mathrm{[}\mathrm{}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{2}\mathrm{\Psi}\mathrm{)}\mathrm{d}{\mathit{\tau}}^{\mathrm{2}}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{}\mathrm{2}\mathrm{\Phi}\mathrm{\left)}\mathrm{d}{\mathit{x}}^{\mathrm{2}}\mathrm{\right]}\mathit{.}$(1)The gauge invariant potentials Φ and Ψ are related to the Bardeen (1980) potentials Φ_{A} and Φ_{H} and to the Kodama & Sasaki (1984) potentials Ψ_{KS} and Φ_{KS} in the following way: Ψ = Φ_{A} = Ψ_{KS} and Φ = −Φ_{H} = −Φ_{KS}. Throughout the paper we use a metric signature (− , + + +) and follow the notation of Ma & Bertschinger (1995); the speed of light is set to c = 1, except where explicitly stated otherwise.
We define the equation of state , where and are the average pressure and energy density. The sound speed c_{s} is defined in the fluid rest frame in terms of pressure and density perturbations as $\mathit{\delta p}\mathrm{\left(}\mathit{k,a}\mathrm{\right)}\mathrm{=}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\left(}\mathit{k,a}\mathrm{\right)}\mathit{\delta \rho}\mathrm{\left(}\mathit{k,a}\mathrm{\right)}$. The anisotropic stress σ(k,a) (equivalent to π_{T} in the notation of Kodama & Sasaki 1984) is the scalar part of the offdiagonal spacespace stress energy tensor perturbation. The set of functions { H,Φ,Ψ } describing the metric is formally equivalent to the set of functions $\mathrm{\left\{}\mathit{w,}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{,\sigma}\mathrm{\right\}}$ (Ballesteros et al. 2012).
Specific theories typically cover only subsets of this function space and thus make specific predictions for their form. In the following sections we will discuss the particular theories that we consider in this paper.
3.1. Background parameterizations
The first main “category” of theories we describe includes parameterizations of background quantities. Even when we are only interested in constraints on background parameters, we are implicitly assuming a prescription for Dark Energy fluctuations. The conventional approach, that we adopt also here, is to choose a minimallycoupled scalar field model (Wetterich 1988; Ratra & Peebles 1988), also known as quintessence, which corresponds to the choice of a restframe sound speed ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}\mathrm{1}$ (i.e., equal to the speed of light) and σ = 0 (no scalar anisotropic stress). In this case the relativistic sound speed suppresses the dark energy perturbations on subhorizon scales, preventing it from contributing significantly to clustering.
Background parameterizations discussed in this paper include:

(w_{0}, w_{a}) Taylor expansion at first order (and potentially higher orders);

Principal Component Analysis of w(a) (Huterer & Starkman 2003), that allows to estimate constraints on w in independent redshift bins;

general parameterization of any minimally coupled scalar field in terms of three parameters ϵ_{s}, ζ_{s}, ϵ_{∞}. This is a novel way to describe minimally coupled scalar field models without explicitly specifying the form of the potential (Huang et al. 2011);

Dark Energy density as a function of z (including parameterizations such as early Dark Energy).
The specific implementation for each of them is discussed in Sect. 5.1 together with corresponding results. We will conclude the background investigation by describing, in Sect. 5.1.6, a compressed Gaussian likelihood that captures most of the constraining power of the Planck data applied to smooth Dark Energy or curved models (following Mukherjee et al. 2008). The compressed likelihood is useful for example to include more easily the Planck CMB data in Fisherforecasts for future largescale structure surveys.
3.2. Perturbation parameterizations
Modified gravity models (in which gravity is modified with respect to GR) in general affect both the background and the perturbation equations. In this subsection we go beyond background parameterizations and identify two different approaches to constrain MG models, one more theoretically motivated and a second more phenomenological one. We will not embark on a fullscale survey of DE and MG models here, but refer the reader to e.g. Amendola et al. (2013) for more details.
3.2.1. Modified gravity and effective field theory
The first approach starts from a Lagrangian, derived from an effective field theory (EFT) expansion (Cheung et al. 2008), discussed in Creminelli et al. (2009) and Gubitosi et al. (2013) in the context of DE. Specifically, EFT describes the space of (universally coupled) scalar field theories, with a Lagrangian written in unitary gauge that preserves isotropy and homogeneity at the background level, assumes the weak equivalence principle, and has only one extra dynamical field besides the matter fields conventionally considered in cosmology. The action reads: $\begin{array}{ccc}\mathit{S}& \mathrm{=}& \int {\mathrm{d}}^{\mathrm{4}}\mathit{x}\sqrt{\mathrm{}\mathit{g}}\{\frac{{\mathit{m}}_{\mathrm{0}}^{\mathrm{2}}}{\mathrm{2}}\left[\mathrm{1}\mathrm{+}\mathrm{\Omega}\mathrm{\left(}\mathit{\tau}\mathrm{\right)}\right]\mathit{R}\mathrm{+}\mathrm{\Lambda}\mathrm{\left(}\mathit{\tau}\mathrm{\right)}\mathrm{}{\mathit{a}}^{\mathrm{2}}\mathit{c}\mathrm{\left(}\mathit{\tau}\mathrm{\right)}\mathit{\delta}{\mathit{g}}^{\mathrm{00}}\\ & & \u2001\mathrm{+}\frac{{\mathit{M}}_{\mathrm{2}}^{\mathrm{4}}\mathrm{\left(}\mathit{\tau}\mathrm{\right)}}{\mathrm{2}}{\mathrm{(}{\mathit{a}}^{\mathrm{2}}\mathit{\delta}{\mathit{g}}^{\mathrm{00}}\mathrm{)}}^{\mathrm{2}}\mathrm{}\mathit{M\u0305}\begin{array}{c}\mathrm{3}\\ \mathrm{1}\end{array}\mathrm{\left(}\mathit{\tau}\mathrm{\right)}\mathrm{2}{\mathit{a}}^{\mathrm{2}}\mathit{\delta}{\mathit{g}}^{\mathrm{00}}\mathit{\delta}{{\mathit{K}}_{\mathit{\mu}}^{\mathit{\mu}}}^{}\\ & & \u2001\mathrm{}\frac{\mathit{M\u0305}\begin{array}{c}\mathrm{2}\\ \mathrm{2}\end{array}\mathrm{\left(}\mathit{\tau}\mathrm{\right)}}{\mathrm{2}}{\mathrm{(}\mathit{\delta}{{\mathit{K}}_{\mathit{\mu}}^{\mathit{\mu}}}^{\mathrm{)}}}^{\mathrm{2}}\mathrm{}\frac{\mathit{M\u0305}\begin{array}{c}\mathrm{2}\\ \mathrm{3}\end{array}\mathrm{\left(}\mathit{\tau}\mathrm{\right)}}{\mathrm{2}}\mathit{\delta}{\mathit{K}}_{\mathit{\nu}}^{\mathit{\mu}}\mathit{\delta}{\mathit{K}}_{\mathit{\mu}}^{\mathit{\nu}}\mathrm{+}\frac{{\mathit{a}}^{\mathrm{2}}\mathit{M\u0302}\mathrm{2}\mathrm{\left(}\mathit{\tau}\mathrm{\right)}}{\mathrm{2}}\mathit{\delta}{\mathit{g}}^{\mathrm{00}}\mathit{\delta}{\mathit{R}}^{\mathrm{\left(}\mathrm{3}\mathrm{\right)}}\\ & & \u2001\mathrm{+}{\mathit{m}}_{\mathrm{2}}^{\mathrm{2}}\mathrm{\left(}\mathit{\tau}\mathrm{\right)}\left({\mathit{g}}^{\mathit{\mu \nu}}\mathrm{+}{\mathit{n}}^{\mathit{\mu}}{\mathit{n}}^{\mathit{\nu}}\right){\mathit{\partial}}_{\mathit{\mu}}\mathrm{(}{\mathit{a}}^{\mathrm{2}}{\mathit{g}}^{\mathrm{00}}\mathrm{)}{\mathit{\partial}}_{\mathit{\nu}}\mathrm{(}{\mathit{a}}^{\mathrm{2}}{\mathit{g}}^{\mathrm{00}}\mathrm{)}\begin{array}{c}\mathrm{\u23ab}\\ \mathrm{\u23aa}\\ \mathrm{\u23aa}\\ \mathrm{\u23ac}\\ \mathrm{\u23aa}\\ \mathrm{\u23aa}\\ \mathrm{\u23ad}\end{array}\\ & & \u2001\mathrm{+}{\mathit{S}}_{\mathrm{m}}[{\mathit{\chi}}_{\mathit{i}}\mathit{,}{\mathit{g}}_{\mathit{\mu \nu}}]\mathit{.}\end{array}$(2)Here R is the Ricci scalar, δR^{(3)} is its spatial perturbation, ${\mathit{K}}_{\mathit{\nu}}^{\mathit{\mu}}$ is the extrinsic curvature, and m_{0} is the bare (reduced) Planck mass. The matter part of the action, S_{m}, includes all fluid components except dark energy, i.e., baryons, cold dark matter, radiation, and neutrinos. The action in Eq. (2) depends on nine timedependent functions (Bloomfield et al. 2013), here $\mathrm{\left\{}\mathrm{\Omega}\mathit{,c,}\mathrm{\Lambda}\mathit{,}\mathit{M\u0305}\begin{array}{c}\mathrm{3}\\ \mathrm{1}\end{array}\mathit{,}\mathit{M\u0305}\begin{array}{c}\mathrm{4}\\ \mathrm{2}\end{array}\mathit{,}\mathit{M\u0305}\begin{array}{c}\mathrm{2}\\ \mathrm{3}\end{array}\mathit{,}{\mathit{M}}_{\mathrm{2}}^{\mathrm{4}}\mathit{,}\mathit{M\u0302}{}^{\mathrm{2}}\mathit{,}{\mathit{m}}_{\mathrm{2}}^{\mathrm{2}}\mathrm{\right\}}$, whose choice specifies the theory. In this way, EFT provides a direct link to any scalar field theory. A particular subset of EFT theories are the Horndeski (1974) models, which include (almost) all stable scalartensor theories, universally coupled to gravity, with secondorder equations of motion in the fields and depend on five functions of time (Gleyzes et al. 2013; Bellini & Sawicki 2014; Piazza et al. 2014).
Although the EFT approach has the advantage of being very versatile, in practice it is necessary to choose suitable parameterizations for the free functions listed above, in order to compare the action with the data. We will describe our specific choices, together with results for each of them, in Sect. 5.2.
3.2.2. MG and phenomenological parameterizations
The second approach adopted in this paper to test MG is more phenomenological and starts from the consideration that cosmological observations probe quantities related to the metric perturbations, in addition to the expansion rate. Given the line element of Eq. (1), the metric perturbations are determined by the two potentials Φ and Ψ, so that we can model all observationally relevant degrees of freedom by parameterizing these two potentials (or, equivalently, two independent combinations of them) as functions of time and scale. Since a nonvanishing anisotropic stress (proportional to Φ−Ψ) is a generic signature of modifications of GR (Mukhanov et al. 1992; Saltas et al. 2014), the parameterized potentials will correspond to predictions of MG models.
Various parameterizations have been considered in the literature. Some of the more popular (in longitudinal gauge) are:
 1.
Q(a,k), which modifies the relativistic Poisson equation through extraDE clustering according to$\mathrm{}{\mathit{k}}^{\mathrm{2}}\mathrm{\Phi}\mathrm{\equiv}\mathrm{4}\mathit{\pi G}{\mathit{a}}^{\mathrm{2}}\mathit{Q}\mathrm{\left(}\mathit{a,}{k}\mathrm{\right)}\mathit{\rho}\mathrm{\Delta}\mathit{,}$(3)where Δ is the comoving density perturbation;
 2.
μ(a,k) (sometimes also called Y(a,k)), which modifies the equivalent equation for Ψ rather than Φ: $\mathrm{}{\mathit{k}}^{\mathrm{2}}\mathrm{\Psi}\mathrm{\equiv}\mathrm{4}\mathit{\pi G}{\mathit{a}}^{\mathrm{2}}\mathit{\mu}\mathrm{\left(}\mathit{a,}{k}\mathrm{\right)}\mathit{\rho}\mathrm{\Delta ;}$(4)
 3.
Σ(a,k), which modifies lensing (with the lensing/Weyl potential being Φ + Ψ), such that $\mathrm{}{\mathit{k}}^{\mathrm{2}}\mathrm{(}\mathrm{\Phi}\mathrm{+}\mathrm{\Psi}\mathrm{)}\mathrm{\equiv}\mathrm{8}\mathit{\pi G}{\mathit{a}}^{\mathrm{2}}\mathrm{\Sigma}\mathrm{\left(}\mathit{a,}{k}\mathrm{\right)}\mathit{\rho}\mathrm{\Delta ;}$(5)
 4.
η(a,k), which reflects the presence of a nonzero anisotropic stress, the difference between Φ and Ψ being equivalently written as a deviation of the ratio^{2}$\mathit{\eta}\mathrm{\left(}\mathit{a,}{k}\mathrm{\right)}\mathrm{\equiv}\mathrm{\Phi}\mathit{/}\mathrm{\Psi}\mathit{.}$(6)
The four functions above are certainly not independent. It is sufficient to choose two independent functions of time and scale to describe all modifications with respect to General Relativity (e.g. Zhang et al. 2007; Amendola et al. 2008b). Popular choices include: (μ,η), which have a simple functional form for many theories; (μ,Σ), which is more closely related to what we actually observe, given that CMB lensing, weak galaxy lensing and the ISW effect measure a projection or derivative of the Weyl potential Φ + Ψ. Furthermore, redshift space distortions constrain the velocity field, which is linked to Ψ through the Euler equation of motion.
All four quantities, Q, μ, Σ, and η, are free functions of time and scale. Their parameterization in terms of the scale factor a and momentum k will be specified in Sect. 5.2.2, together with results obtained by confronting this class of models with data.
3.3. Examples of particular models
The last approach is to consider particular models. Even though these are in principle included in the case described in Sect. 3.2.1, it is nevertheless still useful to highlight some well known examples of specific interest, which we list below.

Minimallycoupled models beyond simple quintessence.Specifically, we consider “kessence” models, which are definedby an arbitrary sound speed ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$ in addition to a free equation of state parameter w (ArmendarizPicon et al. 2000).

An example of a generalized scalar field model (Deffayet et al. 2010) and of Lorentzviolating massive gravity (Dubovsky 2004; Rubakov & Tinyakov 2008), both in the “equation of state” formalism of Battye & Pearson (2012).

Universal “fifth forces”. We will show results for f(R) theories (Wetterich 1995; Capozziello 2002; Amendola et al. 2007; De Felice & Tsujikawa 2010), which form a subset of all models contained in the EFT approach.

Nonuniversal fifth forces. We will illustrate results for coupled DE models (Amendola 2000), in which dark matter particles feel a force mediated by the DE scalar field.
All these particular models are based on specific actions, ensuring full internal consistency. The reviews by Amendola et al. (2013), Clifton et al. (2012), Joyce et al. (2015) and Huterer et al. (2015) contain detailed descriptions of a large number of models discussed in the literature.
4. Data
We now discuss the data sets we use, both from Planck and in combination with other experiments. As mentioned earlier, if we combine many different data sets (not all of which will be equally reliable) and take them all at face value, we risk attributing systematic problems between data sets to genuine physical effects in DE or MG models. On the other hand, we need to avoid bias in confirming ΛCDM, and remain open to the possibility that some tensions may be providing hints that point towards DE or MG models. While discussing results in Sect. 5, we will try to assess the impact of additional data sets, separating them from the Planck baseline choice, keeping in mind caveats that might appear when considering some of them. For a more detailed discussion of the data sets we refer to Planck Collaboration XIII (2016).
4.1. Planck data sets
4.1.1. Planck lowℓ data
The 2013 papers used WMAP polarization measurements (Bennett et al. 2013) at multipoles ℓ ≤ 23 to constrain the optical depth parameter τ. The corresponding likelihood was denoted “WP” in the 2013 papers.
For the present release, we use in its place a Planck polarization likelihood that is built through lowresolution maps of Stokes Q and U polarization measured by LFI at 70 GHz (excluding data from Surveys 2 and 4), foregroundcleaned with the LFI 30 GHz and HFI 353 GHz maps, used as polarized synchrotron and dust templates, respectively (see Planck Collaboration XI 2016).
The foregroundcleaned LFI 70 GHz polarization maps are processed, together with the temperature map from the Commander component separation algorithm over 94% of the sky (see Planck Collaboration IX 2016 for further details), using the lowℓPlanck temperaturepolarization likelihood. This likelihood is pixelbased, extends up to multipoles ℓ = 29 and masks the polarization maps with a specific polarization mask, which uses 46% of the sky. Use of this likelihood is denoted as “lowP” hereafter.
The Planck lowP likelihood, when combined with the highℓPlanck temperature one, provides a best fit value for the optical depth τ = 0.078 ± 0.019, which is about 1σ lower than the value inferred from the WP polarization likelihood, i.e., τ = 0.089 ± 0.013, in the Planck 2013 papers (see also Planck Collaboration XIII 2016). However, we find that the LFI 70 GHz and WMAP polarization maps are extremely consistent when both are cleaned with the HFI 353 GHz polarized dust template, as discussed in more detail in Planck Collaboration XI (2016).
4.1.2. Planck highℓ data
Following Planck Collaboration XV (2014), the highℓ part of the likelihood (30 <ℓ< 2500) uses a Gaussian approximation, $\mathrm{}\mathrm{log}\mathrm{\mathcal{L}}\mathrm{\left(}\mathit{C\u0302}\mathrm{\right}\mathit{C}\mathrm{\left(}\mathit{\theta}\mathrm{\right)}\mathrm{)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}\mathit{C\u0302}\mathrm{}\mathit{C}\mathrm{\left(}\mathit{\theta}\mathrm{\right)}{\mathrm{)}}^{\mathrm{T}}\mathrm{\xb7}{\mathrm{C}}^{1}\mathrm{\xb7}\mathrm{(}\mathit{C\u0302}\mathrm{}\mathit{C}\mathrm{(}\mathit{\theta}\mathrm{\left)}\mathrm{\right)}\mathrm{+}\mathrm{const}\mathit{.}\mathit{,}$(7)with Ĉ the data vector, C(θ) the model with parameters θ and C the covariance matrix. The data vector consists of the temperature power spectra of the best CMB frequencies of the HFI instrument. Specifically, as discussed in Planck Collaboration XI (2016), we use 100 GHz, 143 GHz and 217 GHz halfmission crossspectra, measured on the cleanest part of the sky, avoiding the Galactic plane, as well as the brightest point sources and regions where the CO emission is the strongest. The point source masks are specific to each frequency. We retain, 66% of the sky for the 100 GHz map, 57% for 143 GHz, and 47% for 217 GHz. All the spectra are corrected for beam and pixel window functions. Not all crossspectra and multipoles are included in the data vector; specifically, the TT100 × 143 and 100 × 217 crossspectra, which do not bring much extra information, are discarded. Similarly, we only use multipoles in the range 30 < ℓ < 1200 for 100 × 100 and 30 < ℓ < 2000 for 143 × 143, discarding modes where the S/N is too low. We do not coadd the different crossfrequency spectra, since, even after masking the highest dustcontaminated regions, each crossfrequency spectrum has a different, frequencydependent residual foreground contamination, which we deal with in the model part of the likelihood function.
The model, C(θ) can be rewritten as ${\mathit{C}}_{\mathit{\mu ,\nu}}\mathrm{\left(}\mathit{\theta}\mathrm{\right)}\mathrm{=}\frac{{\mathit{C}}^{\mathrm{cmb}}\mathrm{+}{\mathit{C}}_{\mathit{\mu ,\nu}}^{\mathrm{fg}}\mathrm{\left(}\mathit{\theta}\mathrm{\right)}}{\sqrt{{\mathit{A}}_{\mathit{\mu}}{\mathit{A}}_{\mathit{\nu}}}}\mathit{,}$(8)where C^{cmb} is the set of CMB C_{ℓ}s, which is independent of frequency, ${\mathit{C}}_{\mathit{\mu ,\nu}}^{\mathrm{f}\mathit{g}}\mathrm{\left(}\mathit{\theta}\mathrm{\right)}$ is the foreground model contribution to the crossfrequency spectrum μ × ν, and A_{μ} the calibration factor for the μ × μ spectrum. We retain the following contributions in our foreground modelling: dust; clustered cosmic infrared background (CIB); thermal SunyaevZeldovich (tSZ) effect; kinetic SunyaevZeldovich (kSZ) effect; tSZCIB crosscorrelations; and point sources. The dust, CIB and point source contributions are the dominant contamination. Specifically, dust is the dominant foreground at ℓ< 500, while the diffuse point source term (and CIB for the 217 × 217) dominates the small scales. All our foreground models are based upon smooth C_{ℓ} templates with free amplitudes. All templates but the dust are based on analytical models, as described in Planck Collaboration XI (2016). The dust is based on a mask difference of the 545 GHz map and is well described by a power law of index n = −2.63, with a wide bump around ℓ = 200. A prior for the dust amplitude is computed from the crossspectra with the 545 GHz map. We refer the reader to Planck Collaboration XI (2016) for a complete description of the foreground model. The overall calibration for the 100 × 100 and 217 × 217 power spectra free to vary within a prior measured on a small fraction of the sky near the Galactic pole.
The covariance matrix C accounts for the correlations due to the mask and is computed following the equations in Planck Collaboration XV (2014). The fiducial model used to compute the covariance is based on a joint fit of ΛCDM and nuisance parameters. The covariance includes the nonGaussianity of the noise, but assumes Gaussian statistics for the dust. The nonwhiteness of the noise is estimated from the difference between the cross and autohalf mission spectra and accounted for in an approximate manner in the covariance. Different Monte Carlo based corrections are applied to the covariance matrix calculation to account for inaccuracies in the analytic formulae at large scales (ℓ < 50) and when dealing with the point source mask. Beamshape uncertainties are folded into the covariance matrix. A complete description of the computation and its validation is discussed in Planck Collaboration XI (2016).
The TT unbinned covariance matrix is of size about 8000 × 8000. When adding the polarization, the matrix has size 23 000 × 23 000, which translates into a significant memory requirement and slows the likelihood computation considerably. We thus bin the data and covariance matrix, using a variable binsize scheme, to reduce the data vector dimension by about a factor of ten. We checked that for the ΛCDM model, including single parameter classical extensions, the cosmological and nuisance parameter fits are identical with or without binning.
4.1.3. Planck CMB lensing
Gravitational lensing by largescale structure introduces dependencies in CMB observables on the latetime geometry and clustering, which otherwise would be degenerate in the primary anisotropies (Hu 2002; Lewis & Challinor 2006). This provides some sensitivity to dark energy and latetime modifications of gravity from the CMB alone. The source plane for CMB lensing is the lastscattering surface, so the peak sensitivity is to lenses at z ≈ 2 (i.e., halfway to the lastscattering surface) with typical sizes of order 10^{2} Mpc. Although this peak lensing redshift is rather high for constraining simple latetime dark energy models, CMB lensing deflections at angular multipoles ℓ ≲ 60 have sources extending to low enough redshift that DE becomes dynamically important (e.g., Pan et al. 2014).
The main observable effects of CMB lensing are a smoothing of the acoustic peaks and troughs in the temperature and polarization power spectra, the generation of significant nonGaussianity in the form of a nonzero connected 4point function, and the conversion of Emode to Bmode polarization. The smoothing effect on the power spectra is included routinely in all results in this paper. We additionally include measurements of the power spectrum ${\mathit{C}}_{\mathit{\ell}}^{\mathit{\phi \phi}}$ of the CMB lensing potential φ, which are extracted from the Planck temperature and polarization 4point functions, as presented in Planck Collaboration XV (2016) and discussed further below. Lensing also produces 3point nonGaussianity, which peaks in squeezed configurations, due to the correlation between the lensing potential and the ISW effect in the largeangle temperature anisotropies. This effect has been measured at around 3σ with the fullmission Planck data (Planck Collaboration XV 2016; Planck Collaboration XXI 2016). Although in principle this is a further probe of DE (Verde & Spergel 2002) and MG (Acquaviva et al. 2004), we do not include these T–φ correlations in this paper as the likelihood was not readily available. We plan however to test this effect in future work.
The construction of the CMB lensing likelihood we use in this paper is described fully in Planck Collaboration XV (2016); see also Planck Collaboration XIII (2016). It is a simple Gaussian approximation in the estimated ${\mathit{C}}_{\mathit{\ell}}^{\mathit{\phi \phi}}$ bandpowers, covering the multipole range 40 ≤ ℓ ≤ 400. The ${\mathit{C}}_{\mathit{\ell}}^{\mathit{\phi \phi}}$ are estimated from the fullmission temperature and polarization 4point functions, using the SMICA componentseparated maps (Planck Collaboration IX 2016) over approximately 70% of the sky. A large number of tests of internal consistency of the estimated ${\mathit{C}}_{\mathit{\ell}}^{\mathit{\phi \phi}}$ to different data cuts (e.g., whether polarization is included, or whether individual frequency bands are used in place of the SMICA maps) are reported in Planck Collaboration XV (2016). All such tests are passed for the conservative multipole range 40 ≤ ℓ ≤ 400 that we adopt in this paper. For multipoles ℓ> 400, there is marginal evidence of systematic effects in reconstructions of the lensing deflections from temperature anisotropies alone, based on curlmode tests. Reconstructing the lensing deflections on large angular scales is very challenging because of the large “meanfield” due to survey anisotropies, which must be carefully subtracted with simulations. We conservatively adopt a minimum multipole of ℓ = 40 here, although the results of the null tests considered in Planck Collaboration XV (2016) suggest that this could be extended down to ℓ = 8. For Planck, the multipole range 40 ≤ ℓ ≤ 400 captures the majority of the S/N on ${\mathit{C}}_{\mathit{\ell}}^{\mathit{\phi \phi}}$ for ΛCDM models, although this restriction may be more lossy in extended models. The Planck 2014 lensing measurements are the most significant to date (the amplitude of ${\mathit{C}}_{\mathit{\ell}}^{\mathit{\phi \phi}}$ is measured at greater than 40σ), and we therefore choose not to include lensing results from other CMB experiments in this paper.
4.1.4. Planck CMB polarization
The TE and EE likelihood follows the same principle as the TT likelihood described in Sect. 4.1.2. The data vector is extended to contain the TE and EE crosshalfmission power spectra of the same 100 GHz, 143 GHz, and 217 GHz frequency maps. Following Planck Collaboration Int. XXX (2016), we mask the regions where the dust intensity is important, and retain 70%, 50%, and 41% of the sky for our three frequencies. We ignore any other polarized galactic emission and in particular synchrotron, which has been shown to be negligible, even at 100 GHz. We use all of the crossfrequency spectra, using the multipole range 30 <ℓ < 1000 for the 100 GHz crossspectra and 500 <ℓ < 2000 for the 217 GHz crossspectra. Only the 143 × 143 spectrum covers the full 30 < ℓ < 2000 range. We use the same beams as for the TT spectra and do not correct for leakage due to beam mismatch. A complete description of the beam mismatch effects and correction is described in Planck Collaboration XI (2016).
The model is similar to the TT one. We retain a single foreground component accounting for the polarized emission of the dust. Following Planck Collaboration Int. XXX (2016), the dust C_{ℓ} template is a power law with index n = −2.4. A prior for the dust amplitude is measured in the crosscorrelation with the 353 GHz maps. The calibration parameters are fixed to unity.
The covariance matrix is extended to polarization, as described in Planck Collaboration XI (2016), using the correlation between the TT, TE, and EE spectra. It is computed similarly to the TT covariance matrix, as described in Sect. 4.1.2.
In this paper we will only show results that include CMB highℓ polarization data where we find that it has a significant impact. DE and MG can in principle also affect the Bmode power spectrum through lensing of Bmodes (if the lensing Weyl potential is modified) or by changing the position and amplitude of the primordial peak (Antolini et al. 2013; Pettorino & Amendola 2015), including modifications of the sound speed of gravitational waves (Amendola et al. 2014; Raveri et al. 2015). Due to the unavailability of the likelihood, results from Bmode polarization are left to future work.
4.2. Background data combination
We identify a first basic combination of data sets that we mostly rely on, for which we have a high confidence that systematics are under control. Throughout this paper, we indicate for simplicity with “BSH” the combination BAO + SNIa + H_{0}, which we now discuss in detail.
4.2.1. Baryon acoustic oscillations
Baryon acoustic oscillations (BAO) are the imprint of oscillations in the baryonphoton plasma on the matter power spectrum and can be used as a standard ruler, calibrated to the CMBdetermined sound horizon at the end of the drag epoch. Since the acoustic scale is so large, BAO are largely unaffected by nonlinear evolution. As in the cosmological parameter paper, Planck Collaboration XIII (2016), BAO is considered as the primary data set to break parameter degeneracies from CMB measurements and offers constraints on the background evolution of MG and DE models. The BAO data can be used to measure both the angular diameter distance D_{A}(z), and the expansion rate of the Universe H(z) either separately or through the combination ${\mathit{D}}_{\mathrm{V}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}{\left[\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{\mathrm{2}}{\mathit{D}}_{\mathrm{A}}^{\mathrm{2}}\mathrm{(}\mathit{z}\mathrm{)}\frac{\mathit{cz}}{\mathit{H}\mathrm{\left(}\mathit{z}\mathrm{\right)}}\right]}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{\xb7}$(9)As in Planck Collaboration XIII (2016) we use: the SDSS Main Galaxy Sample at z_{eff} = 0.15 (Ross et al. 2015); the Baryon Oscillation Spectroscopic Survey (BOSS) “LOWZ” sample at z_{eff} = 0.32 (Anderson et al. 2014); the BOSS CMASS (i.e. “constant mass” sample) at z_{eff} = 0.57 of Anderson et al. (2014); and the sixdegreeField Galaxy survey (6dFGS) at z_{eff} = 0.106 (Beutler et al. 2011). The first two measurements are based on peculiar velocity field reconstructions to sharpen the BAO feature and reduce the errors on the quantity D_{V}/r_{s}; the analysis in Anderson et al. (2014) provides constraints on both D_{A}(z_{eff}) and H(z_{eff}). In all cases considered here the BAO observations are modelled as distance ratios, and therefore provide no direct measurement of H_{0}. However, they provide a link between the expansion rate at low redshift and the constraints placed by Planck at z ≈ 1100.
4.2.2. Supernovae
TypeIa supernovae (SNe) are among the most important probes of expansion and historically led to the general acceptance that a DE component is needed (Riess et al. 1998; Perlmutter et al. 1999). Supernovae are considered as “standardizable candles” and so provide a measurement of the luminosity distance as a function of redshift. However, the absolute luminosity of SNe is considered uncertain and is marginalized out, which also removes any constraints on H_{0}.
Consistently with Planck Collaboration XIII (2016), we use here the analysis by Betoule et al. (2013) of the “Joint Lightcurve Analysis” (JLA) sample. JLA is constructed from the SNLS and SDSS SNe data, together with several samples of low redshift SNe. Cosmological constraints from the JLA sample^{3} are discussed by Betoule et al. (2014), and as mentioned in Planck Collaboration XIII (2016) the constraints are consistent with the 2013 and 2104 Planck values for standard ΛCDM.
4.2.3. The Hubble constant
The CMB measures mostly physics at the epoch of recombination, and so provides only weak direct constraints about lowredshift quantities through the integrated SachsWolfe effect and CMB lensing. The CMBinferred constraints on the local expansion rate H_{0} are model dependent, and this makes the comparison to direct measurements interesting, since any mismatch could be evidence of new physics.
Here, we rely on the reanalysis of the Riess et al. (2011, hereafter R11) Cepheid data made by Efstathiou (2014, hereafter E14). By using a revised geometric maser distance to NGC 258 from Humphreys et al. (2013), E14 obtains the following value for the Hubble constant: ${\mathit{H}}_{\mathrm{0}}\mathrm{=}\mathrm{(}\mathrm{70.6}\mathrm{\pm}\mathrm{3.3}\mathrm{)}\hspace{0.17em}\mathrm{km}\hspace{0.17em}{\mathrm{s}}^{1}\hspace{0.17em}{\mathrm{Mpc}}^{1}\mathit{,}$(10)which is within 1σ of the Planck TT+lowP estimate. In this paper we use Eq. (10) as a conservative H_{0} prior. We note that the 2015 Planck TT+lowP value is perfectly consistent with the 2013 Planck value (Planck Collaboration XVI 2014) and so the tension with the R11 H_{0} determination is still present at about 2.4σ. We refer to the cosmological parameter paper Planck Collaboration XIII (2016) for a more comprehensive discussion of the different values of H_{0} present in the literature.
4.3. Perturbation data sets
The additional freedom present in MG models can be calibrated using external data that test perturbations in particular. In the following we describe other available data sets that we included in the grid of runs for this paper.
4.3.1. Redshift space distortions
Observations of the anisotropic clustering of galaxies in redshift space permit the measurement of their peculiar velocities, which are related to the Newtonian potential Ψ via the Euler equation. This, in turn, allows us to break a degeneracy with gravitational lensing that is sensitive to the combination Φ + Ψ. Galaxy redshift surveys now provide very precise constraints on redshiftspace clustering. The difficulty in using these data is that much of the signal currently comes from scales where nonlinear effects and galaxy bias are significant and must be accurately modelled (see, e.g., the discussions in Bianchi et al. 2012; GilMarín et al. 2012). Moreover, adopting the wrong fiducial cosmological model to convert angles and redshifts into distances can bias measurements of the rateofgrowth of structure (Reid et al. 2013; Howlett et al. 2015). Significant progress in the modelling has been achieved in the last few years, so we shall focus here on the most recent (and relatively conservative) studies. A compilation of earlier measurements can be found in the references above.
In linear theory, anisotropic clustering along the line of sight and in the transverse directions measures the combination f(z)σ_{8}(z), where the growth rate is defined by $\mathit{f}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\frac{\mathrm{d}\mathrm{ln}{\mathit{\sigma}}_{\mathrm{8}}}{\mathrm{d}\mathrm{ln}\mathrm{a}}\mathrm{\xb7}$(11)where σ_{8} is calculated including all matter and neutrino density perturbations. Anisotropic clustering also contains geometric information from the AlcockPaczynski (AP) effect (Alcock & Paczynski 1979), which is sensitive to ${\mathit{F}}_{\mathrm{AP}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}{\mathit{D}}_{\mathrm{A}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathit{H}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathit{.}$(12)In addition, fits which constrain RSD frequently also measure the BAO scale, D_{V}(z) /r_{s}, where r_{s} is the comoving sound horizon at the drag epoch, and D_{V} is given in Eq. (9). As in Planck Collaboration XIII (2016) we consider only analyses which solve simultaneously for the acoustic scale, F_{AP} and fσ_{8}.
The Baryon Oscillation Spectroscopic Survey (BOSS) collaboration have measured the power spectrum of their CMASS galaxy sample (Beutler et al. 2014) in the range k = 0.01−0.20 h Mpc^{1}. Samushia et al. (2014) have estimated the multipole moments of the redshiftspace correlation function of CMASS galaxies on scales > 25 h^{1} Mpc. Both papers provide tight constraints on the quantity fσ_{8}, and the constraints are consistent. The Samushia et al. (2014) result was shown to behave marginally better in terms of smallscale bias compared to mock simulations, so we choose to adopt this as our baseline result. Note that when we use the data of Samushia et al. (2014), we exclude the measurement of the BAO scale, D_{V}/r_{s}, from Anderson et al. (2014), to avoid double counting.
The Samushia et al. (2014) results are expressed as a 3 × 3 covariance matrix for the three parameters D_{V}/r_{s}, F_{AP} and fσ_{8}, evaluated at an effective redshift of z_{eff} = 0.57. Since Samushia et al. (2014) do not apply a density field reconstruction in their analysis, the BAO constraints are slightly weaker than, though consistent with, those of Anderson et al. (2014).
4.3.2. Galaxy weak lensing
The distortion of the shapes of distant galaxies by largescale structure along the line of sight (weak gravitational lensing or cosmic shear) is particularly important for constraining DE and MG, due to its dependence on the growth of fluctuations and the two scalar metric potentials.
Currently the largest weak lensing (WL) survey is the Canada France Hawaii Telescope Lensing Survey (CFHTLenS), and we make use of two data sets from this survey:

1.
2D CFHTLenS data (Kilbingeret al. 2013), whose shearcorrelation functions ξ^{±} are estimated in the angular range 0.9 to 296.5 arcmin;

2.
the tomographic CFHTLenS blue galaxy sample (Heymans et al. 2013), whose data have an intrinsic alignment signal consistent with zero, eliminating the need to marginalize over any additional nuisance parameters, and where the shear correlation functions are estimated in six redshift bins, each with an angular range 1.7 <θ< 37.9 arcmin.
Since these data are not independent we do not combine them, but rather check the consistency of our results with each. The galaxy lensing convergence power spectrum, ${\mathit{P}}_{\mathit{ij}}^{\mathit{\kappa}}\mathrm{\left(}\mathit{\ell}\mathrm{\right)}$, can be written in terms of the Weyl potential, P_{Φ + Ψ}, by $\begin{array}{ccc}{\mathit{P}}_{\mathit{ij}}^{\mathit{\kappa}}\mathrm{\left(}\mathit{\ell}\mathrm{\right)}\mathrm{\approx}\mathrm{2}{\mathit{\pi}}^{\mathrm{2}}\mathit{\ell}\mathrm{\int}\frac{\mathrm{d}\mathit{\chi}}{\mathit{\chi}}{\mathit{g}}_{\mathit{i}}\left(\mathit{\chi}\right){\mathit{g}}_{\mathit{j}}\left(\mathit{\chi}\right){\mathit{P}}_{\mathrm{\Phi}\mathrm{+}\mathrm{\Psi}}\left(\mathit{\ell}\mathit{/}\mathit{\chi ,\chi}\right)\mathit{,}& & \end{array}$(13)where we have made use of the Limber approximation in flat space, and χ is the comoving distance. The lensing efficiency is given by $\begin{array}{ccc}{\mathit{g}}_{\mathit{i}}\mathrm{\left(}\mathit{\chi}\mathrm{\right)}\mathrm{=}{\mathrm{\int}}_{\mathit{\chi}}^{\mathrm{\infty}}\mathrm{d}{\mathit{\chi}}^{\mathrm{\prime}}{\mathit{n}}_{\mathit{i}}\left({\mathit{\chi}}^{\mathrm{\prime}}\right)\frac{{\mathit{\chi}}^{\mathrm{\prime}}\mathrm{}\mathit{\chi}}{{\mathit{\chi}}^{\mathrm{\prime}}}\mathit{,}& & \end{array}$(14)where n_{i}(χ) is the radial distribution of source galaxies in bin i. In the case of no anisotropic stress and no additional clustering from the DE, the convergence power spectrum can be written in the usual form $\begin{array}{ccc}{\mathit{P}}_{\mathit{ij}}^{\mathit{\kappa}}\mathrm{\left(}\mathit{\ell}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{9}}{\mathrm{4}}{\mathrm{\Omega}}_{\mathrm{m}}^{\mathrm{2}}{\mathit{H}}_{\mathrm{0}}^{\mathrm{4}}{\mathrm{\int}}_{\mathrm{0}}^{\mathrm{\infty}}\frac{{\mathit{g}}_{\mathit{i}}\mathrm{\left(}\mathit{\chi}\mathrm{\right)}{\mathit{g}}_{\mathit{j}}\mathrm{\left(}\mathit{\chi}\mathrm{\right)}}{{\mathit{a}}^{\mathrm{2}}\mathrm{\left(}\mathit{\chi}\mathrm{\right)}}\mathit{P}\mathrm{(}\mathit{\ell}\mathit{/}\mathit{\chi ,\chi}\mathrm{)}\mathrm{d}\mathit{\chi}\mathit{.}\end{array}$(15)However, in this paper we always use the full Weyl potential to compute the theoretical WL predictions. The convergence can also be written in terms of the correlation functions ξ^{±} via $\begin{array}{ccc}{\mathit{\xi}}_{\mathit{i,j}}^{\mathrm{\pm}}\mathrm{\left(}\mathit{\theta}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi}}\mathrm{\int}\mathrm{d}\mathit{\ell}\hspace{0.17em}\mathit{\ell}\hspace{0.17em}{\mathit{P}}_{\mathit{ij}}^{\mathit{\kappa}}\mathrm{\left(}\mathit{\ell}\mathrm{\right)}{\mathit{J}}_{\mathrm{\pm}}\mathrm{\left(}\mathit{\ell \theta}\mathrm{\right)}\mathit{,}& & \end{array}$(16)where the Bessel functions are J_{+} = J_{0} and J_{−} = J_{4}.
In this paper we need to be particularly careful about the contribution of nonlinear scales to ξ^{±}, since the behaviour of MG models in the nonlinear regime is not known very precisely. The standard approach is to correct the power spectrum on nonlinear scales using the Halofit fitting function. Since its inception, there have been several revisions to improve the agreement with Nbody simulations. We use the following convention to label the particular Halofit model:

1.
the original model of Smithet al. (2003);

2.
an update from higher resolution Nbody simulations, to include the effect of massive neutrinos (Bird et al. 2012);

3.
an update to improve the accuracy on small scales^{4};

4.
an update from higher resolution Nbody simulations, including DE cosmologies with constant equation of state (Takahashi et al. 2012).
Given this correction, one can scale the Weyl potential transfer functions by the ratio of the nonlinear to linear matter power spectrum ${\mathit{T}}_{\mathrm{\Phi}\mathrm{+}\mathrm{\Psi}}\mathrm{\left(}\mathit{k,z}\mathrm{\right)}\mathrm{\to}{\mathit{T}}_{\mathrm{\Phi}\mathrm{+}\mathrm{\Psi}}\mathrm{\left(}\mathit{k,z}\mathrm{\right)}\sqrt{\frac{{\mathit{P}}_{\mathit{\delta}}^{\mathrm{nonlin}}\mathrm{\left(}\mathit{k,z}\mathrm{\right)}}{{\mathit{P}}_{\mathit{\delta}}^{\mathrm{lin}}\mathrm{\left(}\mathit{k,z}\mathrm{\right)}}}\mathrm{\xb7}$(17)Both Kilbinger et al. (2013) and Heymans et al. (2013) quote a “conservative” set of cuts to mitigate uncertainty over the nonlinear modelling scheme. For the 2D analysis of Kilbinger et al. (2013) angular scales θ< 17′ are excluded for ξ^{+}, and θ< 54′ for ξ^{−}. For the tomographic analysis of Heymans et al. (2013), angular scales θ< 3′ are excluded for ξ^{+} for any bin combination involving the two lowest redshift bins, and no cut is applied for the highest four redshift bins. For ξ^{−}, angular scales θ< 30′ are excluded for any bin combination involving the four lowest redshift bins, and θ< 16′ for the highest two bins.
Fig. 2 Ω_{m}–σ_{8} constraints for tomographic lensing from Heymans et al. (2013), using a very conservative angular cut, as described in the text (see Sect. 4.3.2). We show results using linear theory, nonlinear corrections from Halofit (HL) versions 1, 4, marginalization over baryonic AGN feedback (BF), and intrinsic alignment (IA) (the latter two using nonlinear corrections and Halofit 4). Coloured points indicate H_{0} values from WL+HL4. 
These cuts, however, may be insufficient for our purposes, since we are interested in extensions to ΛCDM. We therefore choose a very conservative set of cuts to mitigate the total contribution from nonlinear scales. In order to select these cuts we choose the baseline Planck TT+lowP ΛCDM cosmology as described in Planck Collaboration XIII (2016), for which one can use Eq. (15). The cuts are then chosen by considering $\mathrm{\Delta}{\mathit{\chi}}^{\mathrm{2}}\mathrm{=}\mathrm{}{\mathit{\chi}}_{\mathrm{lin}}^{\mathrm{2}}\mathrm{}{\mathit{\chi}}_{\mathrm{nonlin}}^{\mathrm{2}}\mathrm{}$ of the WL likelihood as a function of angular cut. In order for this to remain Δχ^{2}< 1 for each of the Halofit versions, we find it necessary to remove ξ^{−} entirely from each data set, and exclude θ< 17′ for ξ^{+} for both the 2D and tomographic bins. We note that a similar approach to Kitching et al. (2014) could also be followed using 3D CFHTLenS data, where the choice of cut is more well defined in kspace, however the likelihood for this was not available at the time of this paper.
On small scales, the effects of intrinsic alignments and baryonic feedback can also become significant. In order to check the robustness of our cuts to these effects we adopt the same methodology of MacCrann et al. (2015). Using the same baseline model and choosing Halofit version 4, we scale the matter power spectrum by an active galactic nuclei (AGN) component, derived from numerical simulations (van Daalen et al. 2011), marginalizing over an amplitude α_{AGN}. The AGN baryonic feedback model has been shown by HarnoisDéraps et al. (2015) to provide the best fit to smallscale CFHTLens data. For intrinsic alignment we adopt the model of Bridle & King (2007), including the additional nonlinear alignment contributions to ξ^{±}, and again marginalizing over an amplitude α_{IA}. For more details on this procedure, we refer the reader to MacCrann et al. (2015).
The robustness of our ultraconservative cuts to nonlinear modelling, baryonic feedback and intrinsic alignment marginalization, is illustrated in Fig. 2 for the tomographic data, with similar constraints obtained from 2D data. Assuming the same base ΛCDM cosmology, and applying priors of Ω_{b}h^{2} = 0.0223 ± 0.0009, n_{s} = 0.96 ± 0.02, and 40 km s^{1} Mpc^{1}<H_{0}< 100 km s^{1} Mpc^{1} to avoid overfitting the model, we find that the WL likelihood is insensitive to nonlinear physics. We therefore choose to adopt the tomographic data with the ultraconservative cuts as our baseline data set.
Table of models tested in this paper.
4.4. Combining data sets
We show for convenience in Table 1 the schematic summary of models. All models have been tested for the combinations: Planck, Planck+BSH, Planck+WL, Planck+BAO/RSD and Planck+WL+BAO/RSD. Throughout the text, unless otherwise specified, Planck refers to the baseline Planck TT+lowP combination. The effects of CMB lensing and Planck TT,TE,EE polarization have been tested on all runs above and are, in particular, used to constrain the amount of DE at early times. For each of them we indicate the section in which the model is described and the corresponding figures. In addition, all combinations in the table have been tested with and without CMB lensing. The impact of Planck highℓ polarization has been tested on all models for the combination Planck+BAO+SNe+H_{0}.
5. Results
We now proceed by illustrating in detail the models and parameterizations described in Sect. 3, through presenting results for each of them. The structure of this section is as follows. We start in Sect. 5.1 with smooth dark energy models that are effectively parameterized by the expansion history of the Universe alone. In Sect. 5.2 we study the constraints on the presence of nonnegligible dark energy perturbations, both in the context of general modified gravity models described through effective field theories and with phenomenological parameterizations of the gravitational potentials and their combinations, as illustrated in Sect. 3.2.2. The last part, Sect. 5.3, illustrates results for a range of particular examples often considered in the literature.
5.1. Background parameterizations
In this section, we consider models where DE is a generic quintessencelike component with equation of state w ≡ p/ρ, where p and ρ are the spatially averaged (background) DE pressure and density. Although it is important to include, as we do, DE perturbations, models in this section have a sound speed that is equal to the speed of light, which means that they are smooth on subhorizon scales (see Sect. 3.1 for more details). We start with Taylor expansions and a principal component analysis of w in a fluid formalism, then consider actual quintessence models parameterized through their potentials and finally study the limits that can be put on the abundance of DE density at early times. At the end of the subsection we provide the necessary information to compress the Planck CMB power spectrum into a 4parameter Gaussian likelihood for applications where the full likelihood is too unwieldy.
5.1.1. Taylor expansions of w and w_{0},w_{a} parameterization
Fig. 3 Parameterization { w_{0},w_{a} } (see Sect. 5.1.1). Marginalized posterior distributions for w_{0}, w_{a}, H_{0} and σ_{8} for various data combinations. The tightest constraints come from the Planck TT+lowP+BSH combination, which indeed tests background observations, and is compatible with ΛCDM. 
If the dark energy is not a cosmological constant with w = −1 then there is no reason why w should remain constant. In order to test a timevarying equation of state, we expand w(a) in a Taylor series. The first order corresponds to the { w_{0},w_{a} } case, also discussed in Planck Collaboration XIII (2016): $\mathit{w}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{=}{\mathit{w}}_{\mathrm{0}}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{}\mathit{a}\mathrm{)}{\mathit{w}}_{\mathit{a}}\mathit{.}$(18)We use the parameterized postFriedmann (PPF) model of Hu & Sawicki (2007) and Fang et al. (2008) to allow for values w< −1 (note that there is another PPF formalism discussed in Baker et al. 2014a). Marginalized posterior distributions for w_{0}, w_{a}, H_{0} and σ_{8} are shown in Fig. 3 and the corresponding 2D contours can be found in Fig. 4 for w_{a} vs. w_{0} and for σ_{8} vs. Ω_{m}. Results from Planck TT+lowP+BSH data are shown in blue and corresponds to the combination we consider the most secure, which in this case also gives the strongest constraints. This is expected, since the BAO and SNe data included in the BSH combination provide the best constraints on the background expansion rate.
Fig. 4 Marginalized posterior distributions of the (w_{0},w_{a}) parameterization (see Sect. 5.1.1) for various data combinations. The best constraints come from the priority combination and are compatible with ΛCDM. The dashed lines indicate the point in parameter space (− 1,0) corresponding to the ΛCDM model. CMB lensing and polarization do not significantly change the constraints. Here Planck indicates Planck TT+lowP. 
Results for weak lensing (WL) and redshift space distortions (RSD) are also shown, both separately and combined. The constraints from these probes are weaker, since we are considering a smooth dark energy model where the perturbations are suppressed on small scales. While the WL data appear to be in slight tension with ΛCDM, according to the green contours shown in Fig. 4, the difference in total χ^{2} between the bestfit in the { w_{0},w_{a} } model and in ΛCDM for Planck TT+lowP+WL is Δχ^{2} = −5.6, which is not very significant for 2 extra parameters (for normal errors a 2σ deviation corresponds to a χ^{2} absolute difference of 6.2). The WL contributes a Δχ^{2} of −2.0 and the $\mathrm{\Delta}{\mathit{\chi}}_{\mathrm{CMB}}^{\mathrm{2}}\mathrm{=}\mathrm{}\mathrm{3.3}$ (virtually the same as when using Planck TT+lowP alone, for which $\mathrm{\Delta}{\mathit{\chi}}_{\mathrm{CMB}}^{\mathrm{2}}\mathrm{=}\mathrm{}\mathrm{3.2}$, which seems to indicate that WL is not in tension with Planck TT+lowP within a (w_{0},w_{a}) cosmology). However, as also discussed in Planck Collaboration XIII (2016), these data combinations prefer very high values of H_{0}, which is visible also in the third panel of Fig. 3. The combination Planck TT+lowP+BSH, on the other hand, is closer to ΛCDM, with a total χ^{2} difference between (w_{0},w_{a}) and ΛCDM of only −0.8. We also show in Fig. 5 the equation of state reconstructed as a function of redshift from the linear expansion in the scale factor a for different combinations of data.
Fig. 5 Reconstructed equation of state w(z) as a function of redshift (see Sect. 5.1.1), when assuming a Taylor expansion of w(z) to firstorder (N = 1 in Eq. (19)), for different combinations of the data sets. The coloured areas show the regions which contain 95% of the models. The central blue line is the median line for Planck TT+lowP+BSH. Here Planck indicates Planck TT+lowP. 
One might wonder whether it is reasonable to stop at first order in w(a). We have therefore tested a generic expansion in powers of the scale factor up to order N: $\mathit{w}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{=}{\mathit{w}}_{\mathrm{0}}\mathrm{+}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}\mathrm{(}\mathrm{1}\mathrm{}\mathit{a}{\mathrm{)}}^{\mathit{i}}{\mathit{w}}_{\mathit{i}}\mathit{.}$(19)We find that all parameters are very stable when allowing higher order polynomials; the w_{i} parameters are weakly constrained and going from N = 1 (the linear case) to N = 2 (quadratic case) to N = 3 (cubic expansion) does not improve the goodness of fit and stays compatible with ΛCDM, which indicates that a linear parameterization is sufficient.
5.1.2. 1parameter varying w
A simple example of a varying w model that can be written in terms of one extra parameter only (instead of w_{0},w_{a}) was proposed in Gott & Slepian (2011), Slepian et al. (2014), motivated in connection to a DE minimallycoupled scalar field, slowly rolling down a potential $\frac{\mathrm{1}}{\mathrm{2}}{\mathit{m}}^{\mathrm{2}}{\mathit{\phi}}^{\mathrm{2}}$, analogous to the one predicted in chaotic inflation (Linde 1983). More generally, one can fully characterize the background by expanding a varying equation of state $\mathit{w}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{\equiv}\mathrm{}\mathrm{1}\mathrm{+}\mathit{\delta w}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{\approx}\mathrm{}\mathrm{1}\mathrm{+}\mathit{\delta}{\mathit{w}}_{\mathrm{0}}\mathrm{\times}{\mathit{H}}_{\mathrm{0}}^{\mathrm{2}}\mathit{/}{\mathit{H}}^{\mathrm{2}}\mathrm{\left(}\mathit{z}\mathrm{\right)}$, where: $\frac{{\mathit{H}}^{\mathrm{2}}\mathrm{\left(}\mathit{z}\mathrm{\right)}}{{\mathit{H}}_{\mathrm{0}}^{\mathrm{2}}}\mathrm{\approx}{\mathrm{\Omega}}_{\mathrm{m}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{\mathrm{3}}\mathrm{+}{\mathrm{\Omega}}_{\mathrm{de}}{\left[\frac{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{\mathrm{3}}}{{\mathrm{\Omega}}_{\mathrm{m}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{\mathrm{3}}\mathrm{+}{\mathrm{\Omega}}_{\mathrm{de}}}\right]}^{\mathit{\delta}{\mathit{w}}_{\mathrm{0}}\mathit{/}{\mathrm{\Omega}}_{\mathrm{de}}}\mathit{,}$(20)at first order in δw_{0}, which is then the only extra parameter. Marginalized posterior contours in the plane h–δw_{0} are shown in Fig. 6. The tightest constraints come from the combination Planck TT+lowP+lensing+BSH that gives δw_{0} = −0.008 ± 0.068 at 68% confidence level, which slightly improves constraints found by Aubourg et al. (2015).
Fig. 6 Marginalized posterior contours in the h–δw_{0} plane are shown for 1parameter varying w models (see Sect. 5.1.2) for different data combinations. Here Planck indicates Planck TT+lowP. 
5.1.3. Principal Component Analysis on w(z)
A complementary way to measure the evolution of the equation of state, which is better able to model rapid variations, proceeds by choosing w in N fixed bins in redshift and by performing a principal component analysis to uncorrelate the constraints. We consider N = 4 different bins in z and assume that w has a constant value p_{i} in each of them. We then smooth the transition from one bin to the other such that: $\begin{array}{ccc}\mathit{w}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}{\mathit{p}}_{\mathit{i}\mathrm{}\mathrm{1}}\mathrm{+}\mathrm{\Delta}\mathit{w}\left(\mathrm{tanh}\left[\frac{\mathit{z}\mathrm{}{\mathit{z}}_{\mathit{i}}}{\mathit{s}}\right]\mathrm{+}\mathrm{1}\right)\mathrm{for}\hspace{0.17em}\mathit{z}\mathit{<}{\mathit{z}}_{\mathit{i}}\mathit{,}\mathrm{i}\mathit{\u03f5}\mathrm{\left\{}\mathrm{1}\mathrm{,}\mathrm{4}\mathrm{\right\}}\mathit{,}& & \end{array}$(21)with Δw ≡ (p_{i}−p_{i−1}) / 2, a smoothing scale s = 0.02, and a binning z_{i} = (0,0.2,0.4,0.6,1.8). We have tested also a larger number of bins (up to N = 18) and have found no improvement in the goodness of fit.
The constraints on the vector p_{i}(i = 1,...,N) of values that w(z) can assume in each bin is difficult to interpret, due to the correlations between bins. To uncorrelate the bins, we perform a principal component analysis (Huterer & Starkman 2003; Huterer & Cooray 2005; Said et al. 2013). We first run COSMOMC (Lewis & Bridle 2002) on the original binning values p_{i}; then extract the covariance matrix that refers to the parameters we want to constrain: $\mathit{C}\mathrm{\equiv}\mathrm{\u27e8}\mathit{p}{\mathit{p}}^{\mathrm{T}}\mathrm{\u27e9}\mathrm{}\mathrm{\u27e8}\mathit{p}\mathrm{\u27e9}\mathrm{\u27e8}{\mathit{p}}^{\mathrm{T}}\mathrm{\u27e9}\mathit{,}$(22)where p is the vector of parameters p_{i} and p^{T} is its transpose. We calculate the Fisher matrix, F = C^{1}, and diagonalize it, F = O^{T}ΛO, where Λ is diagonal and O is the orthogonal matrix whose columns are the eigenvectors of the Fisher matrix. We then define $\begin{array}{c}{}_{\mathrm{\U0010ff65}}\\ \mathit{W}\end{array}\mathrm{=}{\mathit{O}}^{\mathrm{T}}{\mathrm{\Lambda}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{O}$ (e.g., Huterer & Cooray 2005) and normalize this such that its rows sum up to unity; this matrix can be used to find the new vector $\mathit{q}\mathrm{=}\begin{array}{c}{}_{\mathrm{\U0010ff65}}\\ \mathit{W}\end{array}\mathit{p}$ of uncorrelated parameters that describe w(z). This choice of $\begin{array}{c}{\mathrm{\u02dc}}_{}\\ \mathit{W}\end{array}$ has been shown to be convenient, since most of the weights (i.e., the rows of $\begin{array}{c}{\mathrm{\U0010ff65}}_{}\\ \mathit{W}\end{array}$) are found to be positive and fairly well localized in redshift. In Fig. 7 (lower panel) we show the weights for each bin as a function of redshift. Because they overlap only partially, we can assume the binning to be the same as the original one and attach to each of them error bars corresponding to the mean and standard deviations of the q values. The result is shown in Fig. 7, top panel. The equation of state is compatible with the ΛCDM value w = −1. Note however that this plot contains more information than a Taylor expansion to first order.
Fig. 7 PCA analysis constraints (described in Sect. 5.1.3). The top panel shows the reconstructed equation of state w(z) after the PCA analysis. Vertical error bars correspond to mean and standard deviations of the q vector parameters, while horizontal error bars are the amplitude of the original binning. The bins are not exactly independent but are rather smeared out as illustrated in the bottom panel. The bottom panel shows the PCA corresponding weights on w(z) as a function of redshift for the combination Planck TT+lowP+BSH. In other words, error bars in the top panel correspond therefore to the errors in the q parameters, which are linear combinations of the p parameters, i.e. a smeared out distribution with weights shown in the lower panel. 
5.1.4. Parameterization for a weaklycoupled canonical scalar field
We continue our investigation of background parameterizations by considering a slowly rolling scalar field. In this case, as in inflation, we can avoid writing down an explicit potential V(φ) and instead parameterize w(a) at late times, in the presence of matter, as (Huang et al. 2011) $\mathit{w}\mathrm{=}\mathrm{}\mathrm{1}\mathrm{+}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{\u03f5}}_{\mathrm{s}}{\mathit{F}}^{\mathrm{2}}\left(\frac{\mathit{a}}{{\mathit{a}}_{\mathrm{de}}}\right)\mathit{,}$(23)where the “slope parameter” ϵ_{s} is defined as: ${\mathit{\u03f5}}_{\mathrm{s}}\mathrm{\equiv}{{\mathit{\u03f5}}_{\mathrm{V}}}_{\mathit{a}\mathrm{=}{\mathit{a}}_{\mathrm{de}}}\mathit{,}$(24)with ${\mathit{\u03f5}}_{\mathrm{V}}\mathrm{\equiv}\mathrm{(}\frac{\mathrm{d}\mathrm{ln}\mathit{V}}{\mathrm{d}\mathit{\phi}}{\mathrm{)}}^{\mathrm{2}}{\mathit{M}}_{\mathrm{P}}^{\mathrm{2}}\mathit{/}\mathrm{2}$ being a function of the slope of the potential. Here ${\mathit{M}}_{\mathrm{P}}\mathrm{\equiv}\mathrm{1}\mathit{/}\sqrt{\mathrm{8}\mathit{\pi G}}$ is the reduced Planck mass and a_{de} is the scale factor where the total matter and DE densities are equal. The function F(x) in Eq. (23) is defined as: $\mathit{F}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{\equiv}\frac{\sqrt{\mathrm{1}\mathrm{+}{\mathit{x}}^{\mathrm{3}}}}{{\mathit{x}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\mathrm{}\frac{\mathrm{ln}\left({\mathit{x}}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{+}\sqrt{\mathrm{1}\mathrm{+}{\mathit{x}}^{\mathrm{3}}}\right)}{{\mathit{x}}^{\mathrm{3}}}\mathrm{\xb7}$(25)Equation (23) parameterizes w(a) with one parameter ϵ_{s}, while a_{de} depends on Ω_{m} and ϵ_{s} and can be derived using an approximated fitting formula that facilitates numerical computation (Huang et al. 2011). Positive (negative) values of ϵ_{s} correspond to quintessence (phantom) models.
Equation (23) is only valid for lateUniverse slowroll (ϵ_{V} ≲ 1 and ${\mathit{\eta}}_{\mathrm{V}}\mathrm{\equiv}{\mathit{M}}_{\mathrm{P}}^{\mathrm{2}}{\mathit{V}}^{\mathrm{\prime \prime}}\mathit{/}\mathit{V}\mathrm{\ll}\mathrm{1}$) or the moderateroll (ϵ_{V} ≲ 1 and η_{V} ≲ 1) regime. For quintessence models, where the scalar field rolls down from a very steep potential, at early times ϵ_{V}(a) ≫ 1, however the fractional density Ω_{φ}(a) → 0 and the combination ϵ_{V}(a)Ω_{φ}(a) aprroaches a constant, defined to be a second parameter ϵ_{∞} ≡ lim_{a → 0}ϵ_{V}(a)Ω_{φ}(a).
One could also add a third parameter ζ_{s} to capture the timedependence of ϵ_{V} via corrections to the functional dependence of w(a) at late time. This parameter is defined as the relative difference of $\mathrm{d}\sqrt{{\mathit{\u03f5}}_{\mathrm{V}}{\mathrm{\Omega}}_{\mathit{\phi}}}\mathit{/}\mathrm{d}\mathit{y}$ at a = a_{de} and at a → 0, where $\mathit{y}\mathrm{\equiv}\mathrm{(}\mathit{a}\mathit{/}{\mathit{a}}_{\mathrm{de}}{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathit{/}\sqrt{\mathrm{1}\mathrm{+}\mathrm{(}\mathit{a}\mathit{/}{\mathit{a}}_{\mathrm{de}}{\mathrm{)}}^{\mathrm{3}}}$. If ϵ_{∞} ≪ 1, ζ_{s} is proportional to the second derivative of lnV(φ), but for large ϵ_{∞}, the dependence is more complicated (Huang et al. 2011). In other words, while ϵ_{s} is sensitive to the late time evolution of 1 + w(a), ϵ_{∞} captures its early time behaviour. Quintessence/phantom models can be mapped into ϵ_{s}–ϵ_{∞} space and the classification can be further refined with ζ_{s}. For ΛCDM, all three parameters are zero.
In Fig. 8 we show the marginalized posterior distributions at 68.3% and 95.4% confidence levels in the parameter space ϵ_{s}–Ω_{m}, marginalizing over the other parameters. In Fig. 9 we show the current constraints on quintessence models projected in ϵ_{s}–ϵ_{∞} space. The constraints are obtained by marginalizing over all other cosmological parameters. The models here include exponentials V = V_{0}exp(−λφ/M_{P}) (Wetterich 1988), cosines from pseudoNambu Goldstone bosons (pnGB) V = V_{0} [1 + cos(λφ/M_{P})] (Frieman et al. 1995; Kaloper & Sorbo 2006), power laws V = V_{0}(φ/M_{P})^{− n} (Ratra & Peebles 1988), and models motivated by supergravity (SUGRA) V = V_{0}(φ/M_{P})^{− α}exp [(φ/M_{P})^{2}] (Brax & Martin 1999). The model projection is done with a fiducial Ω_{m} = 0.3 cosmology. We have verified that variations of 1% compared to the fiducial Ω_{m} lead to negligible changes in the constraints.
Mean values and uncertainties for a selection of cosmological parameters are shown in Table 2, for both the 1parameter case (i.e., ϵ_{s} only, with ϵ_{∞} = 0 and ζ_{s} = 0, describing “thawing” quintessence/phantom models, where in the early Universe) and the 3parameter case (general quintessence/phantom models where an earlyUniverse fastrolling phase is allowed). When we vary the data sets and theoretical prior (between the 1parameter and 3parameter cases), the results are all compatible with ΛCDM and mutually compatible with each other. Because ϵ_{s} and ϵ_{∞} are correlated, caution has to be taken when looking at the marginalized constraints in the table. For instance, the constraint on ϵ_{s} is tighter for the 3parameter case, because in this case flatter potentials are preferred in the late Universe in order to slowdown larger from the early Universe. A better view of the mutual consistency can be obtained from Fig. 9. We find that the addition of polarization data does not have a large impact on these DE parameters. Adding polarization data to Planck+BSH shifts the mean of ϵ_{s} by −1 / 6σ and reduces the uncertainty of ϵ_{s} by 20%, while the 95% upper bound on ϵ_{∞} remains unchanged.
Fig. 8 Marginalized posterior distributions showing 68% and 95% C.L. constraints on Ω_{m} and ϵ_{s} for scalar field models (see Sect. 5.1.4). The dashed line for ϵ_{s} = 0 is the ΛCDM model. Here Planck indicates Planck TT+lowP. 
Fig. 9 Marginalized posterior distributions at 68% C.L. and 95% C.L. in the parameter space of ϵ_{s} and ϵ_{∞} for scalar field modes (see Sect. 5.1.4). We have computed ϵ_{s} and ϵ_{∞} for various quintessence potentials V(φ), with the functional forms of V(φ) labelled on the figure. The field φ is in reduced Planck mass M_{P} units. The normalization of V(φ) is computed using Ω_{m} = 0.3. Here Planck indicates Planck TT+lowP. 
5.1.5. Dark energy density at early times
Quintessence models can be divided into two classes, namely cosmologies with or without DE at early times. Although the equation of state and the DE density are related to each other, it is often convenient to think directly in terms of DE density rather than the equation of state. In this section we provide a more direct estimate of how much DE is allowed by the data as a function of time. A key parameter for this purpose is Ω_{e}, which measures the amount of DE present at early times (“early dark energy,” EDE) (Wetterich 2004). Early DE parameterizations encompass features of a large class of dynamical DE models. The amount of early DE influences CMB peaks and can be strongly constrained when including smallscale measurements and CMB lensing. Assuming a constant fraction of Ω_{e} until recent times (Doran & Robbers 2006), the DE density is parameterized as: ${\mathrm{\Omega}}_{\mathrm{de}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{=}\frac{{\mathrm{\Omega}}_{\mathrm{de}}^{\mathrm{0}}\mathrm{}{\mathrm{\Omega}}_{\mathrm{e}}\mathrm{(}\mathrm{1}\mathrm{}{\mathit{a}}^{\mathrm{}\mathrm{3}{\mathit{w}}_{\mathrm{0}}}\mathrm{)}}{{\mathrm{\Omega}}_{\mathrm{de}}^{\mathrm{0}}\mathrm{+}{\mathrm{\Omega}}_{\mathrm{m}}^{\mathrm{0}}{\mathit{a}}^{\mathrm{3}{\mathit{w}}_{\mathrm{0}}}}\mathrm{+}{\mathrm{\Omega}}_{\mathrm{e}}\mathrm{(}\mathrm{1}\mathrm{}{\mathit{a}}^{\mathrm{}\mathrm{3}{\mathit{w}}_{\mathrm{0}}}\mathrm{)}\mathit{.}$(26)This expression requires two parameters in addition to those of ΛCDM, namely Ω_{e} and w_{0}, while ${\mathrm{\Omega}}_{\mathrm{m}}^{\mathrm{0}}\mathrm{=}\mathrm{1}\mathrm{}{\mathrm{\Omega}}_{\mathrm{de}}^{\mathrm{0}}$ is the present matter abundance. The strongest constraints to date were discussed in Planck Collaboration XVI (2014), finding Ω_{e}< 0.010 at 95% CL using Planck combined with WMAP polarization. Here we update the analysis using Planck 2015 data. In Fig. 10 we show marginalized posterior distributions for Ω_{e} for different combination of data sets; the corresponding marginalized limits are shown in Table 3, improving substantially current constraints, especially when the Planck TT,TE,EE+lowP polarization is included, leading to Ω_{e}< 0.0036 at 95% confidence level for Planck TT,TE,EE+lowP+BSH.
Fig. 10 Marginalized posterior distributions for Ω_{e} for the early DE parameterization in Eq. (26) and for different combinations of data (see Sect. 5.1.5). Here Planck indicates Planck TT+lowP. 
As first shown in Pettorino et al. (2013), bounds on Ω_{e} can be weaker if DE is present only over a limited range of redshifts. In particular, EDE reduces structure growth in the period after last scattering, implying a smaller number of clusters as compared to ΛCDM, and therefore a weaker lensing potential to influence the anisotropies at high ℓ. It is possible to isolate this effect by switching on EDE only after last scattering, at a scale factor a_{e} (or equivalently for redshifts smaller than z_{e}). Here we adopt the parameterization “EDE3” proposed in Pettorino et al. (2013) to which we refer for more details: ${\mathrm{\Omega}}_{\mathrm{de}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \mathrm{for}\mathit{a}\mathrm{\le}{\mathit{a}}_{\mathrm{e}}\mathrm{;}\\ & \mathrm{for}{\mathit{a}}_{\mathrm{e}}\mathit{<}\mathit{a}\mathit{<}{\mathit{a}}_{\mathrm{c}}\mathrm{:}\\ & \mathrm{for}\mathit{a}\mathit{>}{\mathit{a}}_{\mathrm{c}}\mathit{.}\end{array}$(27)In this case, early dark energy is present in the time interval a_{e}<a<a_{c}, while outside this interval it behaves as in ΛCDM, including the radiation contribution, unlike in Eq. (26). During that interval in time, there is a nonnegligible EDE contribution, parameterized by Ω_{e}. The constant a_{c} is fixed by continuity, so that the parameters {Ω_{e},a_{e}} fully determine how much EDE there was and how long its presence lasted. We choose four fixed values of a_{e} corresponding to z_{e} = 10, 50, 200 and 1000 and include Ω_{e} as a free parameter in MCMC runs for each value of a_{e}. Results are shown in Fig. 11 where we plot Ω_{e} as a function of the redshift z_{e} at which DE starts to be nonnegligible. The smaller the value of z_{e}, the weaker are the constraints, though still very tight, with Ω_{e} ≲ 2% (95% CL) for z_{e} ≈ 50.
Fig. 11 Amount of DE at early times Ω_{e} as a function of the redshift z_{e} after which early DE is nonnegligible (see Eq. (27), Sect. 5.1.5) for different combinations of data sets. The heights of the columns give the limit at 95% CL on Ω_{e}, as obtained from Monte Carlo runs for the values z_{e} = 10, 50, 200 and 1000. The width of the columns has no physical meaning and is just due to plotting purposes. Here Planck indicates Planck TT+lowP. 
5.1.6. Compressed likelihood
Before concluding the set of results on background parameterizations, we discuss here how to reduce the full likelihood information to few parameters. As discussed for example in Kosowsky et al. (2002) and Wang & Mukherjee (2007), it is possible to compress a large part of the information contained in the CMB power spectrum into just a few numbers^{5}: here we use specifically the CMB shift parameter R (Efstathiou & Bond 1999), the angular scale of the sound horizon at last scattering ℓ_{A} (or equivalently θ_{∗}), as well the baryon density ω_{b} and the scalar spectral index n_{s}. The first two quantities are defined as $\mathit{R}\mathrm{\equiv}\sqrt{{\mathrm{\Omega}}_{\mathrm{m}}{\mathit{H}}_{\mathrm{0}}^{\mathrm{2}}}\hspace{0.17em}{\mathit{D}}_{\mathrm{A}}\mathrm{\left(}{\mathit{z}}_{\mathrm{\ast}}\mathrm{\right)}\mathit{/}\mathit{c,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{\ell}}_{\mathrm{A}}\mathrm{\equiv}\mathit{\pi}{\mathit{D}}_{\mathrm{A}}\mathrm{\left(}{\mathit{z}}_{\mathrm{\ast}}\mathrm{\right)}\mathit{/}{\mathit{r}}_{\mathrm{s}}\mathrm{\left(}{\mathit{z}}_{\mathrm{\ast}}\mathrm{\right)}\mathrm{=}\mathit{\pi}\mathit{/}{\mathit{\theta}}_{\mathrm{\ast}}\mathit{,}$(28)where D_{A}(z) is the comoving angular diameter distance to redshift z, z_{∗} is the redshift for which the optical depth is unity and r_{s}(z_{∗}) = r_{∗} is the comoving size of the sound horizon at z_{∗}. These numbers are effectively observables and they apply to models with either nonzero curvature or a smooth DE component (Mukherjee et al. 2008). It should be noted, however, that the constraints on these quantities, especially on R, are sensitive to changes in the growth of perturbations. This can be seen easily with the help of the “dark degeneracy” (Kunz 2009), i.e., the possibility to absorb part of the dark matter into the dark energy, which changes Ω_{m} without affecting observables. For this reason the compressed likelihood presented here cannot be used for models with low sound speed or modifications of gravity (and is therefore located at the end of this “background” section).
The marginalized mean values and 68% confidence intervals for the compressed likelihood values are shown in Table 4 for Planck TT+lowP. The posterior distribution of { R,ℓ_{A},ω_{b},n_{s} } is approximately Gaussian, which allows us to specify the likelihood easily by giving the mean values and the covariance matrix, as derived from a Monte Carlo Markov chain (MCMC) approach, in this case from the grid chains for the wCDM model. Since these quantities are very close to observables directly derivable from the data, and since smoothly parameterized DE models are all compatible with the Planck observations to a comparable degree, they lead to very similar central values and essentially the same covariance matrix. The Gaussian likelihood in { R,ℓ_{A},ω_{b},n_{s} } given by Table 4 is thus useful for combining Planck temperature and lowℓ polarization data with other data sets and for inclusion in Fisher matrix forecasts for future surveys. This is especially useful when interested in parameters such as { w_{0},w_{a} }, for which the posterior is very nonGaussian and cannot be accurately represented by a direct covariance matrix (as can be seen in Fig. 4).
The quantities that make up the compressed likelihood are supposed to be “early Universe observables” that describe the observed power spectrum and are insensitive to late time physics. However, lensing by largescale structure has an important smoothing effect on the C_{ℓ} and is detected at over 10σ in the power spectrum (see Sect. 5.2 of Planck Collaboration XIII 2016). We checked by comparing different MCMC chains that the compressed likelihood is stable for ΛCDM, wCDM and the { w_{0},w_{a}} model. However, the “geometric degeneracy” in curved models is broken significantly by the impact of CMB lensing on the power spectrum (see Fig. 25 of Planck Collaboration XIII 2016) and for nonflat models one needs to be more careful. For this reason we also provide the ingredients for the compressed likelihood marginalized over the amplitude A_{L} of the lensing power spectrum in the lower part of Table 4. Marginalizing over A_{L} increases the errors in some variables by over 20% and slightly shifts the mean values, giving a more conservative choice for models where the impact of CMB lensing on the power spectrum is nonnegligible.
We notice that the constraints on R,ℓ_{A},ω_{b},n_{s} given in Table 4 for Planck TT+lowP data are significantly weaker than those predicted by Table II of Mukherjee et al. (2008), which were based on the “Planck Blue Book” specifications (Planck Collaboration 2005). This is because these forecasts also used highℓ polarization. If we derive the actual Planck covariance matrix for the Planck TT,TE,EE+lowP likelihood then we find constraints that are about 50% smaller than those given above, and are comparable and even somewhat stronger than those quoted in Mukherjee et al. (2008). The mean values have of course shifted to represent what Planck has actually measured.
5.2. Perturbation parameterizations
Up to now we have discussed in detail the ensemble of background parameterizations, in which DE is assumed to be a smooth fluid, minimally interacting with gravity. General modifications of gravity, however, change both the background and the perturbation equations, allowing for contribution to clustering (via a sound speed different than unity) and anisotropic stress different from zero. Here we illustrate results for perturbation degrees of freedom, approaching MG from two different perspectives, as discussed in Sect. 3. First we discuss results for EFT cosmologies, with a “topdown” approach that starts from the most general action allowed by symmetry and selects from there interesting classes belonging to socalled “Horndeski models”, which, as mentioned in Sect. 3.2.1, include almost all stable scalartensor theories, universally coupled, with secondorder equations of motion in the fields. We then proceed by parameterizing directly the gravitational potentials and their combinations, as illustrated in Sect. 3.2.2. In this way we can test more phenomenologically their effect on lensing and clustering, in a “bottomup” approach from observations to theoretical models.
5.2.1. Modified gravity: EFT and Horndeski models
The first of the two approaches described in Sect. 3.2.1 adopts effective field theory (EFT) to investigate DE (Gleyzes et al. 2013; Gubitosi et al. 2013), based on the action of Eq. (2). The parameters that appear in the action, when choosing the nine timedependent functions $\mathrm{\left\{}\mathrm{\Omega}\mathit{,c,}\mathrm{\Lambda}\mathit{,}\mathit{M\u0305}\begin{array}{c}\mathrm{3}\\ \mathrm{1}\end{array}\mathit{,}\mathit{M\u0305}\begin{array}{c}\mathrm{4}\\ \mathrm{2}\end{array}\mathit{,}\mathit{M\u0305}\begin{array}{c}\mathrm{2}\\ \mathrm{3}\end{array}\mathit{,}{\mathit{M}}_{\mathrm{2}}^{\mathrm{4}}\mathit{,}\mathit{M\u0302}{}^{\mathrm{2}}\mathit{,}{\mathit{m}}_{\mathrm{2}}^{\mathrm{2}}\mathrm{\right\}}$, describe the effective DE. The full background and perturbation equations for this action have been implemented in the publicly available Boltzmann code EFTCAMB (Hu et al. 2014a; Raveri et al. 2014)^{6}. Given an expansion history (which we fix to be ΛCDM, i.e., effectively w = −1) and an EFT function Ω(a), EFTCAMB computes c and Λ from the Friedmann equations and the assumption of spatial flatness (Hu et al. 2014b). As we have seen in Sect. 5.1, for smooth DE models the constraints on the DE equation of state are compatible with w = −1; hence this choice is not a limitation for the following analysis. In addition, EFTCAMB uses a set of stability criteria in order to specify whether a given model is stable and ghostfree, i.e. without negative energy density for the new degrees of freedom. This will automatically place a theoretical prior on the parameter space while performing the MCMC analysis.
The remaining six functions, $\mathit{M\u0305}\begin{array}{c}\mathrm{3}\\ \mathrm{1}\end{array}$, $\mathit{M\u0305}\begin{array}{c}\mathrm{4}\\ \mathrm{2}\end{array}$, $\mathit{M\u0305}\begin{array}{c}\mathrm{2}\\ \mathrm{3}\end{array}$, ${\mathit{M}}_{\mathrm{2}}^{\mathrm{4}}$, , ${\mathit{m}}_{\mathrm{2}}^{\mathrm{2}}$, are internally redefined in terms of the dimensionless parameters α_{i} with i running from 1 to 6: $\begin{array}{ccc}& & {\mathit{\alpha}}_{\mathrm{1}}^{\mathrm{4}}\mathrm{=}\frac{{\mathit{M}}_{\mathrm{2}}^{\mathrm{4}}}{{\mathit{m}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{H}}_{\mathrm{0}}^{\mathrm{2}}}\mathit{,}\u2001{\mathit{\alpha}}_{\mathrm{2}}^{\mathrm{3}}\mathrm{=}\frac{\mathit{M\u0305}\begin{array}{c}\mathrm{3}\\ \mathrm{1}\end{array}}{{\mathit{m}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{H}}_{\mathrm{0}}}\mathit{,}\u2001{\mathit{\alpha}}_{\mathrm{3}}^{\mathrm{2}}\mathrm{=}\frac{\mathit{M\u0305}\begin{array}{c}\mathrm{2}\\ \mathrm{2}\end{array}}{{\mathit{m}}_{\mathrm{0}}^{\mathrm{2}}}\mathit{,}\\ & & {\mathit{\alpha}}_{\mathrm{4}}^{\mathrm{2}}\mathrm{=}\frac{\mathit{M\u0305}\begin{array}{c}\mathrm{2}\\ \mathrm{3}\end{array}}{{\mathit{m}}_{\mathrm{0}}^{\mathrm{2}}}\mathit{,}\u2001{\mathit{\alpha}}_{\mathrm{5}}^{\mathrm{2}}\mathrm{=}\frac{\mathit{M\u0302}\mathrm{2}}{{\mathit{m}}_{\mathrm{0}}^{\mathrm{2}}}\mathit{,}\u2001{\mathit{\alpha}}_{\mathrm{6}}^{\mathrm{2}}\mathrm{=}\frac{{\mathit{m}}_{\mathrm{2}}^{\mathrm{2}}}{{\mathit{m}}_{\mathrm{0}}^{\mathrm{2}}}\mathrm{\xb7}\end{array}$We will always demand that $\begin{array}{ccc}{\mathit{m}}_{\mathrm{2}}^{\mathrm{2}}& \mathrm{=}& \mathrm{0}\mathrm{(}\mathrm{orequivalently}{\mathit{\alpha}}_{\mathrm{6}}^{\mathrm{2}}\mathrm{=}\mathrm{0}\mathrm{)}\mathit{,}\\ \mathit{M\u0305}\begin{array}{c}\mathrm{2}\\ \mathrm{3}\end{array}& \mathrm{=}& \end{array}$which eliminates models containing higherorder spatial derivatives (Gleyzes et al. 2013, 2015). In this case the nine functions of time discussed above reduce to a minimal set of five functions of time that can be labelled {α_{M},α_{K},α_{B},α_{T},α_{H}}, in addition to the Planck mass ${\mathit{M}}_{\mathrm{\ast}}^{\mathrm{2}}$ (the evolution of which is determined by H and α_{M}), and an additional function of time describing the background evolution, e.g., H(a). The former are related to the EFT functions via the following relations: $\begin{array}{ccc}& & \\ & & \\ & & \\ & & \\ & & {\mathit{M}}_{\mathrm{\ast}}^{\mathrm{2}}{\mathit{\alpha}}_{\mathrm{T}}\mathrm{=}\mathrm{}\mathit{M\u0305}\begin{array}{c}\mathrm{2}\\ \mathrm{2}\end{array}\mathrm{;}\\ & & {\mathit{M}}_{\mathrm{\ast}}^{\mathrm{2}}{\mathit{\alpha}}_{\mathrm{H}}\mathrm{=}\mathrm{2}\mathit{M\u0302}\mathrm{2}\mathrm{}\mathit{M\u0305}\begin{array}{c}\mathrm{2}\\ \mathrm{2}\end{array}\mathit{.}\end{array}$These five α functions are closer to a physical description of the theories under investigation (Bellini & Sawicki 2014). For example: α_{T} enters in the equation for gravitational waves, affecting their speed and the position of the primordial peak in Bmode polarization; α_{M} affects the lensing potential, but also the amplitude of the primordial polarization peak in Bmodes (Amendola et al. 2014; Raveri et al. 2015; Pettorino & Amendola 2015). It is then possible to relate the desired choice for the Horndeski variables to an appropriate choice of the EFT functions, $\begin{array}{ccc}& & \\ & & \\ & & \\ & & \mathrm{4}{\mathit{M}}_{\mathrm{2}}^{\mathrm{4}}\mathrm{=}{\mathit{M}}_{\mathrm{\ast}}^{\mathrm{2}}{\mathrm{\mathscr{H}}}^{\mathrm{2}}{\mathit{\alpha}}_{\mathrm{K}}\mathrm{}\mathrm{2}\mathit{c,}\\ & & \\ & & \end{array}$where ℋ is the conformal Hubble function, m_{0} the bare Planck mass and M_{∗} the effective Plank mass. Fixing α_{M} corresponds to fixing M_{∗} through Eq. (37). Once α_{T} has been chosen, Ω is obtained from Eq. (38). Finally, α_{B} determines $\mathit{M\u0305}\begin{array}{c}\mathrm{3}\\ \mathrm{1}\end{array}$ via Eq. (34), while the choice of α_{H} fixes via Eq. (42). In this way, our choice of the EFT functions can be guided by the selection of different “physical” scenarios, corresponding to turning on different Horndeski functions.
To avoid possible consistency issues with higher derivatives, we set^{7}$\mathit{M\u0305}\begin{array}{c}\mathrm{2}\\ \mathrm{3}\end{array}\mathrm{=}\mathit{M\u0305}\begin{array}{c}\mathrm{2}\\ \mathrm{2}\end{array}\mathrm{=}\mathrm{0}$ in order to satisfy Eq. (30). From Eqs. (39) and (31) this implies α_{T} = 0, so that tensor waves move with the speed of light. In addition, we set α_{H} = 0 so as to remain in the original class of Horndeski theories. As a consequence, from Eq. (42) and ${\mathit{M}}_{\mathrm{\ast}}^{\mathrm{2}}\mathrm{=}{\mathit{m}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{\Omega}\mathrm{)}$ from Eq. (31). For simplicity we also turn off all other higherorder EFT operators and set $\mathit{M\u0305}\begin{array}{c}\mathrm{3}\\ \mathrm{1}\end{array}\mathrm{=}{\mathit{M}}_{\mathrm{2}}^{\mathrm{4}}\mathrm{=}\mathrm{0}$. Comparing Eqs. (32) and (34), this implies α_{B} = −α_{M}.
In summary, in the following we consider Horndeski models in which α_{M} = −α_{B}, α_{K} is fixed by Eq. (33), with M_{2} = 0 as a function of c and α_{T} = α_{H} = 0. We are thus considering nonminimally coupled “Kessence” type models, similar to the ones discussed in Sawicki et al. (2013).
The only free function in this case is α_{M}, which is linked to Ω through: ${\mathit{\alpha}}_{\mathrm{M}}\mathrm{=}\frac{\mathit{a}}{\mathrm{\Omega}\mathrm{+}\mathrm{1}}\frac{\mathrm{d\Omega}}{\mathrm{da}}\mathrm{\xb7}$(43)By choosing a nonzero α_{M} (and therefore a time evolving Ω) we introduce a nonminimal coupling in the action (see Eq. (2)), which will lead to nonzero anisotropic stress and to modifications of the lensing potential, typical signatures of MG models. Here we will use a scaling ansatz, α_{M} = α_{M0}a^{β}, where α_{M0} is the value of α_{M} today, and β> 0 determines how quickly the modification of gravity decreases in the past.
Integrating Eq. (43) we obtain $\mathrm{\Omega}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{=}\mathrm{exp}\left\{\frac{{\mathit{\alpha}}_{\mathrm{M}\mathrm{0}}}{\mathit{\beta}}{\mathit{a}}^{\mathit{\beta}}\right\}\mathrm{}\mathrm{1}\mathit{,}$(44)which coincides with the builtin exponential model of EFTCAMB for Ω_{0} = α_{M0}/β. The marginalized posterior distributions for the two parameters Ω_{0} and β are plotted in Fig. 12 for different combinations of data. For α_{M0} = 0 we recover ΛCDM. For small values of Ω_{0} and for β = 1, the exponential reduces to the builtin linear evolution in EFTCAMB, $\mathrm{\Omega}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{=}{\mathrm{\Omega}}_{\mathrm{0}}\hspace{0.17em}\mathit{a}\mathit{.}$(45)The results of the MCMC analysis are shown in Table 5. For both the exponential and the linear model we use a flat prior Ω_{0} ∈ [0,1]. For the scaling exponent β of the exponential model we use a flat prior β ∈ (0,3]. For β → 0 the MG parameter α_{M} remains constant and does not go to zero in the early Universe, while for β = 3 the scaling would correspond to M functions in the action (2) which are of the same order as the relative energy density between DE and the dark matter background, similar to the suggestion in Bellini & Sawicki (2014). An important feature visible in Fig. 12 is the sharp cutoff at β ≈ 1.5. This cutoff is due to “viability conditions” that are enforced by EFTCAMB and that reject models due to a set of theoretical criteria (see Hu et al. 2014b for a full list of theoretical priors implemented in EFTCAMB). Disabling some of these conditions allows to extend the acceptable model space to larger β, and we find that the constraints on α_{M0} continue to weaken as β grows further, extending Fig. 12 in the obvious way. We prefer however to use here the current public EFTCAMB version without modifications. A better understanding of whether all stability conditions implemented in the code are really necessary or exclude a larger region than necessary in parameter space will have to be addressed in the future. The posterior distribution of the linear evolution for Ω is shown in Fig. 13 and is compatible with ΛCDM. Finally, it is interesting to note that in both the exponential and the linear expansion, the inclusion of WL data set weakens constraints with respect to Planck TT+lowP alone. This is due to the fact that in these EFT theories, WL and Planck TT+lowP are in tension with each other, WL preferring higher values of the expansion rate with respect to Planck.
Fig. 12 Marginalized posterior distributions at 68% and 95% C.L. for the two parameters α_{M0} and β of the exponential evolution, Ω(a) = exp(Ω_{0}a^{β})−1.0, see Sect. 5.2.1. Here α_{M0} is defined as Ω_{0}β and the background is fixed to ΛCDM. Ω_{M0} = 0 corresponds to the ΛCDM model also at perturbation level. Note that Planck means Planck TT+lowP. Adding WL to the data sets results in broader contours, as a consequence of the slight tension between the Planck and WL data sets. 
Fig. 13 Marginalized posterior distribution of the linear EFT model background parameter, Ω, with Ω parameterized as a linear function of the scale factor, i.e., Ω(a) = α_{M0}a, see Sect. 5.2.1. The equation of state parameter w_{de} is fixed to −1, and therefore, Ω_{0} = 0 will correspond to the ΛCDM model. Here Planck means Planck TT+lowP. Adding CMB lensing to the data sets does not change the results significantly; highℓ polarization tightens the constraints by a few percent, as shown in Table 5. 
Marginalized mean values and 68% CL intervals for the EFT parameters, both in the linear model, α_{M0}, and in the exponential one, { α_{M0},β } (see Sect. 5.2.1).
5.2.2. Modified gravity and the gravitational potentials
The second approach used in this paper to address MG is more phenomenological and, as described in Sect. 3.2.2, starts from directly parameterizing the functions of the gravitational potentials listed in Eqs. (3)−(6). Any choice of two of those functions will fully parameterize the deviations of the perturbations from a smooth DE model and describe the cosmological observables of an MG model.
In Simpson et al. (2013) the amplitude of the deviation with respect to ΛCDM was parameterized similarly to the DErelated case that we will define as case 1 below, but using μ(a) and Σ(a) instead of μ(a,k) and η(a,k)^{8}. They found the constraints μ_{0}−1 = 0.05 ± 0.25 and Σ_{0}−1 = 0.00 ± 0.14 using RSD data from the WiggleZ Dark Energy Survey (Blake et al. 2011) and 6dF Galaxy Survey (6dFGS; Beutler et al. 2012), together with CFHTLenS WL. Baker et al. (2014b) provided forecasts on μ_{0}−1 and Σ_{0}−1 for a future experiment that combines galaxy clustering and tomographic weak lensing measurements. The amplitude of departures from the standard values was parameterized as in Simpson et al. (2013), but a possible scale dependence was introduced. In Zhao et al. (2010), the authors constrained μ_{0} and η_{0} and derived from those the limits on Σ_{0}, using WMAP5 data along with CFHTLenS and ISW data. Together with a principal component analysis, they also constrained μ and η assuming a time evolution of the two functions, introducing a transition redshift z_{s} where the functions move smoothly from an early time value to a late time one; they obtained ${\mathit{\mu}}_{\mathrm{0}}\mathrm{=}\mathrm{1.}{\mathrm{1}}_{0.34}^{\mathrm{+}\mathrm{0.62}}\mathit{,}{\mathit{\eta}}_{\mathrm{0}}\mathrm{=}\mathrm{0.9}{\mathrm{8}}_{1.0}^{\mathrm{+}\mathrm{0.73}}$ for z_{s} = 1 and μ_{0} = 0.87 ± 0.12,η_{0} = 1.3 ± 0.35 for z_{s} = 2. A similar parametrization was also used in Daniel et al. (2010) in terms of μ_{0} and ϖ (equivalent to μ_{0}−1 and η_{0}−1 in our convention) using WMAP5, Union2, COSMOS and CFHTLenS data, both binning these functions in redshift and assuming a time evolution (different from the one we will assume in the following), obtaining −0.83 <μ_{0}< 2.1 and −1.6 <ϖ< 2.7 at 95% confidence level for their present values. In Macaulay et al. (2013) the authors instead parameterized Ψ / Φ (the inverse of η) as (1−ζ) and use RSD data from 6dFGS, BOSS, LRG, WiggleZ and VIPERS galaxy redshift surveys to constrain departures from ΛCDM; they did not assume a functional form for the time evolution of ζ, but rather constrained its value at two different redshifts (z = 0 and z = 1), finding a 2σ tension with the ΛCDM limit (ζ = 0) at z = 1.
In this paper, we choose the pair of functions μ(a,k) (related to the Poisson equation for Ψ) and η(a,k) (related to the gravitational slip), as defined in Eqs. (4) and (6), since these are the functions directly implemented in the publicly available code MGCAMB^{9} (Zhao et al. 2009; Hojjati et al. 2011) integrated in the latest version of CosmoMC.
Other functional choices can be easily derived from them (Baker et al. 2014b). We then parameterize μ and η as follows. Since the Planck CMB data span three orders of magnitude in ℓ, it seems sensible to allow for two scales to be present: $\begin{array}{ccc}\mathit{\mu}\mathrm{\left(}\mathit{a,k}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{+}{\mathit{f}}_{\mathrm{1}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\frac{\mathrm{1}\mathrm{+}{\mathit{c}}_{\mathrm{1}}\mathrm{(}\mathit{\lambda H}\mathit{/}\mathit{k}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{1}\mathrm{+}\mathrm{(}\mathit{\lambda H}\mathit{/}\mathit{k}{\mathrm{)}}^{\mathrm{2}}}\mathrm{;}& & \\ \mathit{\eta}\mathrm{\left(}\mathit{a,k}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{+}{\mathit{f}}_{\mathrm{2}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\frac{\mathrm{1}\mathrm{+}{\mathit{c}}_{\mathrm{2}}\mathrm{(}\mathit{\lambda H}\mathit{/}\mathit{k}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{1}\mathrm{+}\mathrm{(}\mathit{\lambda H}\mathit{/}\mathit{k}{\mathrm{)}}^{\mathrm{2}}}\mathrm{\xb7}& & \end{array}$For large length scales (small k), the two functions reduce to μ → 1 + f_{1}(a)c_{1} and η → 1 + f_{2}(a)c_{2}; for small length scales (large k), one has μ → 1 + f_{1}(a) and η → 1 + f_{2}(a). In other words, we implement scale dependence in a minimal way, allowing μ and η to go to two different limits for small and large scales. Here the f_{i} are functions of time only, while the c_{i} and λ parameters are constants. The c_{i} give us information on the scale dependence of μ and η, but the f_{i} measure the amplitude of the deviation from standard GR, corresponding to μ = η = 1.
We choose to parameterize the time dependence of the f_{i}(a) functions as

1.
coefficients related to the DE density,f_{i}(a) = E_{ii}Ω_{DE}(a);

2.
timerelated evolution, f_{i}(a) = E_{i1} + E_{i2}(1−a).
The first choice is motivated by the expectation that the contribution of MG to clustering and to the anisotropic stress is proportional to its effective energy density, as is the case for matter and relativistic particles. The second parameterization provides a complementary approach to the first: E_{i1} describes the MG contribution at late times, while E_{i2} is relevant at early times. Therefore the adoption of the timerelated evolution allows, in principle, for deviations from the standard behaviour also at high redshift, while the parameterization connected to the DE density leads by definition to (μ,η) → 1 at high redshift, since the redshift evolution is tied to that of Ω_{DE}(z).
Fig. 14 68% and 95% contour plots for the two parameters { μ_{0}−1,η_{0}−1 } obtained by evaluating Eqs. (46) and (47) at the present time when no scale dependence is considered (see Sect. 5.2.2). We consider both the DErelated (left panel) and timerelated evolution cases (right panel). Results are shown for the scaleindependent case (c_{1} = c_{2} = 1). In the labels, Planck stands for Planck TT+lowP. 
For case 1 (referred to as “DErelated” parameterization) we then have five free parameters, E_{11}, c_{1}, E_{22}, c_{2}, and λ, while for case 2 (the “timerelated” parameterization) we have two additional parameters, E_{12} and E_{21}. The choice above looks very similar to the BZ parameterization (Bertschinger & Zukin 2008) for the quasistatic limit of f(R) and scalartensor theories. However, we emphasise that Eqs. (46) and (47) should not be seen as a quasistatic limit of any specific theory, but rather as a (minimal) way to allow for (arbitrary) scale dependence, since the data cover a sufficiently wide range of scales. Analogously to the EFT approach discussed in the previous section, we set the background evolution to be the same as in ΛCDM, so that w = −1. In this way the additional parameters purely probe the perturbations.
The effect of the E_{ii} parameters on the CMB temperature and lensing potential power spectra has been shown in Fig. 1 for the “DErelated” choice. In the temperature spectrum the amplitude of the ISW effect is modified; the lensing potential changes more than the temperature spectrum for the same amplitude of the E_{ii} parameter.
We ran Monte Carlo simulations to compare the theoretical predictions with different combinations of the data for both cases 1 and 2. For both choices we tested whether scale dependence plays a role (via the parameters c and λ) with respect to the scaleindependent case in which we fix c_{1} = c_{2} = 1. Results show that a scale dependence of μ and η does not lead to a significantly smaller χ^{2} with respect to the scaleindependent case, both for the DErelated and timerelated parameterizations. Therefore there is no gain in adding c_{i} and λ as extra degrees of freedom. For this reason, in the following we will mainly show results obtained for the scaleindependent parameterization.
Table 6 shows results for the DErelated case for different combinations of the data. Adding the BSH data sets to the Planck TT+lowP data does not significantly increase the constraining power on MG parameters; Planck polarization also has little impact. On the contrary, the addition of RSD data tightens the constraints significantly. The WL contours, including the ultraconservative cut that removes dependence on nonlinear physics, result in weaker constraints. In the table, μ_{0}−1 and η_{0}−1 are obtained by reconstructing Eqs. (46) and (47) from E_{11} and E_{22} at the present time. In addition, the present value of the Σ parameter, defined in Eq. (5), can be obtained from μ and η as Σ = (μ/ 2)(1 + η) using Eqs. (4) and (6).
Some tension appears, in particular, when plotting the marginalized posterior distributions in the planes (μ_{0}−1, η_{0}−1) and (μ_{0}−1, Σ_{0}−1), as shown in Figs. 14 and 15. Here the constraints on the two parameters that describe the perturbations in MG are simultaneously taken into account. In Fig. 14, left and right panels refer to the DErelated and timerelated parameterizations defined in Sect. 5.2.2, respectively, while the dashed lines indicate the values predicted in ΛCDM. Interestingly, results appear similar in both parameterizations. In the DErelated case (left panel), the ΛCDM point lies at the border of the 2σ contour, already when considering Planck TT+lowP alone. More precisely, when looking at the goodness of fit, with respect to the standard ΛCDM assumption, the MG scenario (which includes two extra parameters E_{11} and E_{22}) leads to an improvement of Δχ^{2} = −6.3 when using Planck TT+lowP (similarly divided between lowP and TT) and of Δχ^{2} = −6.4 when including BSH (with a $\mathrm{\Delta}{\mathit{\chi}}_{\mathrm{CMB}}^{\mathrm{2}}\mathrm{~}\mathrm{}\mathrm{5.6}$ equally divided between TT and lowP). When Planck data (TT+lowP) are combined also with WL data, the tension increases to Δχ^{2} = −10.6 (with the CMB still contributing about the same amount, $\mathrm{\Delta}{\mathit{\chi}}_{\mathrm{CMB}}^{\mathrm{2}}\mathrm{=}\mathrm{}\mathrm{6.0}$). When considering Planck TT+lowP+BAO/RSD, Δχ^{2} = −8.1 with respect to ΛCDM while, when combining both WL and BAO/RSD, the tension is maximal, with Δχ^{2} = −10.8 and ${\mathit{\chi}}_{\mathrm{CMB}}^{\mathrm{2}}\mathrm{=}\mathrm{}\mathrm{6.9}$. There is instead less tension for the timerelated parameterization, as is visible in the right panel of Fig. 14.
Fig. 15 Marginalized posterior distributions for 68% and 95% C.L. for the two parameters {μ_{0}−1,Σ_{0}−1} obtained by evaluating Eqs. (46) and (47) at the present time in the DErelated parametrization when no scale dependence is considered (see Sect. 5.2.2). Σ is obtained as Σ = (μ/ 2)(1 + η). The timerelated evolution would give similar contours. In the labels, Planck stands for Planck TT+lowP. 
Once the behaviour of the coefficients in the two parameterizations is known, we can use Eq. (46) to reconstruct the evolution of μ and η with scale factor (or redshift, equivalently). In Fig. 16 we choose to show the linear combination 2 [μ(z,k)−1] + [η(z,k)−1], which corresponds approximately to the maximum degeneracy line in the 2 dimensional μ−1,η−1 parameter space, which allows us to better visualize the joint constraints on μ and η and their maximal allowed departure from ΛCDM. As expected, the DErelated dependence forces the combination to be compatible with ΛCDM in the past, when the DE density is negligible; the timerelated parameterization, instead allows for a larger variation in the past.
The tension can be understood by noticing that the best fit power spectrum corresponds to a value of μ and η (E_{11} = −0.3, E_{22} = 2.2 for Planck TT+lowP) close to the thick long dashed line shown in Fig. 1 for demonstration. This model leads to less power in the CMB at large scales and a higher lensing potential, which is slightly preferred by the data points with respect to ΛCDM. This explains also why the MG parameters are somewhat degenerate with the lensing amplitude A_{L} (which is an “unphysical” parameter redefining the lensing amplitude that affects the CMB power spectrum). As discussed in Planck Collaboration XIII (2016) (see for ex. Sect. 5.1.2), ΛCDM would lead to a value of A_{L} (Calabrese et al. 2008) somewhat larger than 1. When varying it in MG, we find a mean value of ${\mathit{A}}_{\mathrm{L}}\mathrm{=}\mathrm{1.11}{\mathrm{6}}_{0.13}^{\mathrm{+}\mathrm{0.095}}$ which is compatible with A_{L} = 1 at 1σ. The price to pay is the tension with ΛCDM in MG parameter space, which compensates the need for a higher A_{L} that one would have in ΛCDM. The CMB lensing likelihood extracted from the 4point function of the Planck maps Planck Collaboration XV (2016) on the other hand does not prefer a higher lensing potential and agrees well with ΛCDM. For this reason the tension is reduced when we add CMB lensing, as shown in Fig. 17. We also note that constraints for this class of model are sensitive to the estimation of the optical depth τ. Smaller values of τ tend to shift the results further away from ΛCDM.
In order to have a quick overall estimate of the tension for all cases discussed above, we then show in Table 7 the marginalized mean and 68% CL errors for the linear combination 2 [μ(z,k)−1] + [η(z,k)−1]. In the table, we indicate in brackets, for convenience, the “tension” with ΛCDM for each case. This is the maximum allowed tension, since it is calculated along the maximum degeneracy direction. The DErelated parameterization is more in tension with ΛCDM than the timerelated one. The maximum tension reaches 3σ when including WL and BAO/RSD, being therefore mainly driven by external data sets. The inclusion of CMB lensing shifts results towards ΛCDM, as discussed.
Marginalized mean values and 68% C.L. errors on the present day value of the function 2 [μ(z,k)−1] + [η(z,k)−1], which corresponds to the (approximate) maximum degeneracy line identified within the 2 dimensional posterior distributions.
Finally, in general, μ and η depend not only on redshift but also on scale, via the parameters (c_{i}, λ). When marginalizing over them, constraints become weaker, as expected. The comparison with the scaleindependent case is shown in Fig. 18 for Planck TT+lowP+BSH and different values of k. When allowing for scale dependence, the tension with ΛCDM is washed out by the weakening of the constraints and the goodness of fit does not improve with respect to the scale independent case.
Fig. 16 Redshift dependence of the function 2 [μ(z,k)−1] + [η(z,k)−1], defined in Eqs. (46) and (47), which corresponds to the maximum degeneracy line identified within the 2 dimensional posterior distributions. This combination shows the strongest allowed tension with ΛCDM. The left panel refers to the DErelated case while the right panel refers to the timerelated evolution (see Sect. 5.2.2). In both panels, no scale dependence is considered. The coloured areas show the regions containing 68% and 95% of the models. In the labels, Planck stands for Planck TT+lowP. 
Fig. 17 68% and 95% marginalised posterior distributions for the two parameters { μ_{0}−1,η_{0}−1 } obtained by evaluating Eqs. (46) and (47) at the present time when no scale dependence is considered (see Sect. 5.2.2). Here we show the effect of CMB lensing, which shifts the contours towards ΛCDM. In the labels, Planck stands for Planck TT+lowP. 
Fig. 18 68% and 95% contour plots for the two parameters {μ_{0}(k)−1, η_{0}(k)−1} obtained by evaluating Eqs. (46) and (47) at the present time for the DErelated parameterization (see Sect. 5.2.2). We consider both the scaleindependent and scaledependent cases, choosing k values of 10^{10} Mpc^{1} and 10^{2} Mpc^{1}. 
5.3. Further examples of particular models
Quite generally, DE and MG theories deal with at least one extra degree of freedom that can usually be associated with a scalar field. For “standard” DE theories the scalar field couples minimally to gravity, while in MG theories the field can be seen as the mediator of a fifth force in addition to standard interactions. This happens in scalartensor theories (including f(R) cosmologies), massive gravity, and all coupled DE models, both when matter is involved or when neutrino evolution is affected. Interactions and fifth forces are therefore a common characteristic of many proposed models, the difference being whether the interaction is universal (i.e., affecting all species with the same coupling, as in scalartensor theories) or is different for each species (as in coupled DE, Wetterich 1995; Amendola 2000 or growing neutrino models; Fardon et al. 2004; Amendola et al. 2008a). In the following we will test well known examples of particular models within all these classes.
5.3.1. Minimally coupled DE: sound speed and kessence
In minimally coupled quintessence models, the sound speed is ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}\mathrm{1}$ and DE does not contribute significantly to clustering. However, in socalled “kessence” models, the kinetic term in the action is generalised to an arbitrary function of (∇φ)^{2} (ArmendarizPicon et al. 2000): the sound speed can then be different from the speed of light and if c_{s} ≪ 1, the DE perturbations can become nonnegligible on subhorizon scales and impact structure formation. To test this scenario we have performed a series of analyses where we allow for a constant equation of state parameter w and a constant speed of sound ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$ (with a uniform prior in log c_{s}). We find that the limits on w do not change from the quintessence case and that there is no significant constraint on the DE speed of sound using current data. This can be understood as follows: on scales larger than the sound horizon and for w close to −1, DE perturbations are related to dark matter perturbations through Δ_{DE} ≃ (1 + w)Δ_{m}/ 4 and inside the sound horizon they stop growing because of pressure support (see e.g., Creminelli et al. 2009; Sapone & Kunz 2009). In addition, at early times the DE density is much smaller than the matter density, with ρ_{DE}/ρ_{m} = [(1−Ω_{m}) / Ω_{m}] a^{− 3w}. Since the relative DE contribution to the perturbation variable Q(a,k) defined in Eq. (3) scales like ρ_{DE}Δ_{DE}/ (ρ_{m}Δ_{m}), in kessence type models the impact of the DE perturbations on the total clustering is small when 1 + w ≈ 0. For the DE perturbations in kessence to be detectable, the sound speed would have had to be very small, and  1 + w  relatively large.
5.3.2. Massive gravity and generalized scalar field models
We now give two examples of subclasses of Horndeski models, written in terms of an alternative pair of DE perturbation functions (with respect to μ and η used before, for example), given by the anisotropic stress σ and the entropy perturbation Γ: $\mathit{w}\mathrm{\Gamma}\mathrm{=}\frac{\mathit{\delta p}}{\mathit{\rho}}\mathrm{}\frac{\mathrm{d}\mathit{p}}{\mathrm{d}\mathit{\rho}}\mathit{\delta}\mathit{.}$(48)When Γ = 0 the perturbations are adiabatic, that is $\mathit{\delta p}\mathrm{=}\frac{\mathrm{d}\mathit{p}}{\mathrm{d}\mathit{\rho}}\mathit{\delta \rho}$.
For this purpose, it is convenient to adopt the “equation of state” approach described in Battye et al. (2015), Soergel et al. (2015). The gaugeinvariant quantities Γ and σ can be specified in terms of the other perturbation variables, namely δρ, θ, h and η in the scalar sector, and their derivatives.
We then show results for two limiting cases in this formalism, corresponding to Lorentzviolating massive gravity (LVMG) for which (σ ≠ 0,Γ = 0) and generalized scalar field models (GSF) in which the anisotropic stress is zero (σ = 0,Γ ≠ 0).
Lorentzviolating massive gravity (LVMG).
If the Lagrangian is ℒ ≡ ℒ(g_{μν}) (i.e. only written in terms of metric perturbations, as in the EFT action) and one imposes time translation invariance (but not spatial translational invariance), one finds that this corresponds to an extra degree of freedom, ξ^{i}, that has a physical interpretation as an elastic medium, or as Lorentzviolating massive gravity (Dubovsky 2004; Rubakov & Tinyakov 2008; Battye & Pearson 2013). In this case, the scalar equations are characterized by Γ = 0 (the model is adiabatic) and a nonvanishing anisotropic stress: $\mathit{\sigma}\mathrm{=}\mathrm{(}\mathit{w}\mathrm{}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{)}\left[\frac{\mathit{\delta}}{\mathrm{1}\mathrm{+}\mathit{w}}\mathrm{}\mathrm{3}\mathit{\eta}\right]\mathit{,}$(49)including one degree of freedom, the sound speed ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$, which can be related equivalently to the shear modulus of the elastic medium or the Lorentz violating mass. Tensor (gravitational wave) equations will also include a mass term. The low sound speed may lead to clustering of the DE fluid, which allows the data to place constraints on ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$. But as w approaches −1, the DE perturbations are suppressed and the limits on the sound speed weaken. We can take this degeneracy between 1 + w and ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$ into account by using the combination ${\mathit{\lambda}}_{\mathrm{c}}\mathrm{=}\mathrm{}\mathrm{1}\mathrm{+}\mathit{w}{\mathrm{}}^{\mathit{\alpha}}{\mathrm{log}}_{\mathrm{10}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$ in the MCMC analysis, where α = 0.35 was chosen to decorrelate w and λ_{c}. With this, we find Planck TT+lowP+lensing gives lower limit of λ_{c}> −1.6 at 2σ and a tighter one when including BAO/RSD and WL, with λ_{c}> −1.3 at 2σ. For any w ≠ −1 these limits can be translated into limits on ${\mathrm{log}}_{\mathrm{10}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$ by computing λ_{c}/  1 + w  ^{α}. The ΛCDM limit is however fully compatible with the data, i.e. there is no detection of any deviation from w = −1 (and in this limit ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$ is unconstrained).
Generalized scalar field models (GSF).
One can allow for generalized scalar fields by considering a Lagrangian ℒ ≡ ℒ(φ,∂_{μ}φ,∂_{μ}∂_{ν}φ,g_{μν},∂_{α}g_{μν}), in which the dependence on the scalar fields is made explicit, imposing full reparameterization invariance (x^{μ} → x^{ν} + ξ^{μ}), allowing for only linear couplings in ∂_{α}g_{μν} and secondorder field equations. In this case the anisotropic stresses are zero and $\begin{array}{ccc}\mathit{w}\mathrm{\Gamma}& \mathrm{=}& \mathrm{(}\mathit{\alpha}\mathrm{}\mathit{w}\mathrm{)}\{\mathit{\delta}\mathrm{}\mathrm{3}{\mathit{\beta}}_{\mathrm{1}}\mathrm{\mathscr{H}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{w}\mathrm{)}\mathit{\theta}\mathrm{}\frac{\mathrm{3}{\mathit{\beta}}_{\mathrm{2}}\mathrm{\mathscr{H}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{w}\mathrm{)}}{\mathrm{2}{\mathit{k}}^{\mathrm{2}}\mathrm{}\mathrm{6}\mathrm{(}\mathrm{\mathscr{H}\u0307}\mathrm{}{\mathrm{\mathscr{H}}}^{\mathrm{\in}}\mathrm{)}}\mathit{h\u0307}\\ & & \end{array}$(50)This has three extra parameters (α,β_{1},β_{2}), in addition to w. If β_{1} = 1 and β_{2} = 0 this becomes the generalized kessence model. An example of this class of models is “kinetic gravity braiding” (Deffayet et al. 2010) and similar to the nonminimally coupled kessence discussed via EFT in Sect. 5.2.1. The α parameter in Eq. (50) can be now interpreted as a sound speed, unconstrained as in results above. There are however two additional degrees of freedom, β_{1} and β_{2}. RSD data are able to constrain them, with the addition of Planck lensing and WL making only a minor change to the joint constraints. As in the LVMG case, we use a new basis γ_{i} =  1 + w  ^{αi}β_{i} in the MCMC analysis, where α_{1} = 0.2 and α_{2} = 1 were chosen to decorrelate w and γ_{i}. The resulting 2σ upper limits are γ_{1}< 0.67 and γ_{2}< 0.61 (for w> −1), γ_{2}< 2.4 (for w< −1) for the combination of Planck TT+lowP+lensing+BAO/RSD+WL. As for the LVMG case, there is no detection of a deviation from ΛCDM and for w = −1 there are no constraints on β_{1} and β_{2}.
5.3.3. Universal couplings: f(R) cosmologies
A wellinvestigated class of MG models is constituted by the f(R) theories that modify the EinsteinHilbert action by substituting the Ricci scalar with a more general function of itself: $\mathit{S}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}{\mathit{\kappa}}^{\mathrm{2}}}\mathrm{\int}{\mathrm{d}}^{\mathrm{4}}\mathit{x}\sqrt{\mathrm{}\mathit{g}}\mathrm{(}\mathit{R}\mathrm{+}\mathit{f}\mathrm{(}\mathit{R}\mathrm{\left)}\mathrm{\right)}\mathrm{+}\mathrm{\int}{\mathrm{d}}^{\mathrm{4}}\mathit{x}{\mathrm{\mathcal{L}}}_{\mathrm{M}}\mathrm{\left(}{\mathit{\chi}}_{\mathit{i}}\mathit{,}{\mathit{g}}_{\mathit{\mu \nu}}\mathrm{\right)}\mathit{,}$(51)where κ^{2} = 8πG. f(R) cosmologies can be mapped to a subclass of scalartensor theories, where the coupling of the scalar field to the matter fields is universal.
For a fixed background, the Friedmann equation provides a secondorder differential equation for f(R(a)) (see e.g., Song et al. 2007a; Pogosian & Silvestri 2008). One of the initial conditions is usually set by requiring $\underset{\mathit{R}\mathrm{\to}\mathrm{\infty}}{\mathrm{lim}}\frac{\mathit{f}\mathrm{\left(}\mathit{R}\mathrm{\right)}}{\mathit{R}}\mathrm{=}\mathrm{0}\mathit{,}$(52)and the other initial (or boundary condition), usually called B_{0}, is the present value of $\mathit{B}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\frac{{\mathit{f}}_{\mathit{RR}}}{\mathrm{1}\mathrm{+}{\mathit{f}}_{\mathit{R}}}\frac{\mathit{H}\mathit{R\u0307}}{\mathit{H\u0307}\mathrm{}{\mathit{H}}^{\mathrm{2}}}\mathrm{\xb7}$(53)Here, f_{R} and f_{RR} are the first and second derivatives of f(R), and a dot means a derivative with respect to conformal time. Higher values of B_{0} suppress power at large scales in the CMB power spectrum, due to a change in the ISW effect. This also changes the CMB lensing potential, resulting in slightly smoother peaks at higher ℓs (Song et al. 2007b; Schmidt 2008; Bertschinger & Zukin 2008; Marchini et al. 2013).
It is possible to restrict EFTcamb to describe f(R)cosmologies. Given an evolution history for the scale factor and the value of B_{0}, EFTcamb effectively solves the Friedmann equation for f(R). It then uses this function at the perturbation level to evolve the metric potentials and matter fields. The merit of EFTcamb over the other available similar codes is that it checks the model against some stability criteria and does not assume the quasistatic regime, where the scales of interest are still linear but smaller than the horizon and the time derivatives are ignored.
As shown in Fig. 19, there is a degeneracy between the optical depth, τ, and the f(R) parameter, B_{0}. Adding any structure formation probe, such as WL, RSD or CMB lensing, breaks the degeneracy. Figure 20 shows the likelihood of the B_{0} parameter using EFTcamb, where a ΛCDM background evolution is assumed, i.e., w_{DE} = −1.
As the different data sets provide constraints on B_{0} that vary by more than four orders of magnitude, we show plots for log _{10}B_{0}; to make these figures we use a uniform prior in log _{10}B_{0} to avoid distorting the posterior due to prior effects. However, for the limits quoted in the tables we use B_{0} (without log ) as the fundamental quantity and quote 95% limits based on B_{0}. In this way the upper limit on B_{0} is effectively given by the location of the drop in probability visible in the figures, but not influenced by the choice of a lower limit of log _{10}B_{0}. Overall this appears to be the best compromise to present the constraints on the B_{0} parameter. In the plots, the GR value (B_{0} = 0) is reached by a plateau stretching towards minus infinity.
Fig. 19 68% and 95% contour plots for the two parameters, { Log_{10}(B_{0}),τ } (see Sect. 5.3.3). There is a degeneracy between the two parameters for Planck TT+lowP+BSH. Adding lensing will break the degeneracy between the two. Here Planck indicates Planck TT+lowP. 
Fig. 20 Likelihood plots of the f(R) theory parameter, B_{0} (see Sect. 5.3.3). CMB lensing breaks the degeneracy between B_{0} and the optical depth, τ, resulting in lower upper bounds. 
Finally, we note that f(R) models can be studied also with the MGcamb parametrization, assuming the quasistatic limit. We find that for the allowed range of the B_{0} parameter, the results with and without the quasistatic approximation are the same within the uncertainties. The 95% confidence intervals are reported in Table 8. These values show an improvement over the WMAP analysis made with MGcamb (B_{0}< 1 (95 % C.L.) in Song et al. 2007a) and are similar to the limits obtained in Marchini & Salvatelli (2013), Hu et al. (2013; see also Dossett et al. 2014, where data from WiggleZ were used; Cataneo et al. 2015 where they considered galaxy clusters).
5.3.4. Nonuniversal couplings: coupled dark energy
Universal couplings discussed in the previous subsection generally require screening mechanisms to protect baryonic interactions in high density environments, where local measurements are tightly constraining (see e.g. Khoury 2010 and Vikram et al. (2014) for astrophysical constraints). An alternative way to protect baryons is to allow for nonuniversal couplings, in which different species can interact with different strengths: baryons are assumed to be minimally coupled to gravity while other species (e.g., dark matter or neutrinos) may feel a “fifth force”, with a range at cosmological scales.
A fifth force between dark matter particles, mediated by the DE scalar field, is the key ingredient for the coupled DE scenario Amendola (2000). In the Einstein frame, the interaction is described by the Lagrangian $\mathrm{\mathcal{L}}\mathrm{=}\mathrm{}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\partial}}^{\mathit{\mu}}\mathit{\phi}{\mathit{\partial}}_{\mathit{\mu}}\mathit{\phi}\mathrm{}\mathit{V}\mathrm{\left(}\mathit{\phi}\mathrm{\right)}\mathrm{}\mathit{m}\mathrm{\left(}\mathit{\phi}\mathrm{\right)}\mathit{\psi \u0305}\mathit{\psi}\mathrm{+}{\mathrm{\mathcal{L}}}_{\mathrm{kin}}\mathrm{\left[}\mathit{\psi}\mathrm{\right]}\mathit{,}$(54)in which the mass of matter fields ψ is not a constant (as in the standard cosmological model), but rather a function of the DE scalar field φ. A coupling between matter and DE can be reformulated in terms of scalartensor theories or f(R) models (Wetterich 1995, 2015; Pettorino & Baccigalupi 2008) via a Weyl scaling from the Einstein frame (where matter is coupled and gravity is standard) to the Jordan frame (where the gravitational coupling to the Ricci scalar is modified and matter is uncoupled). This is exactly true when the contribution of baryons is neglected.
Dark matter (indicated with the subscript c) and DE densities are then not conserved separately, but coupled to each other: $\begin{array}{ccc}{\mathit{\rho}}_{\mathit{\phi}}^{\mathrm{\prime}}& \mathrm{=}& \mathrm{}\mathrm{3}\mathrm{\mathscr{H}}{\mathit{\rho}}_{\mathit{\phi}}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{w}}_{\mathit{\phi}}\mathrm{)}\mathrm{+}\mathit{\beta}{\mathit{\rho}}_{\mathrm{c}}{\mathit{\phi}}^{\mathrm{\prime}}\mathit{,}\\ {\mathit{\rho}}_{\mathrm{c}}^{\mathrm{\prime}}& \mathrm{=}& \mathrm{}\mathrm{3}\mathrm{\mathscr{H}}{\mathit{\rho}}_{\mathrm{c}}\mathrm{}\mathit{\beta}{\mathit{\rho}}_{\mathrm{c}}{\mathit{\phi}}^{\mathrm{\prime}}\mathit{.}\end{array}$(55)Here each component is treated as a fluid with stress energy tensor ${{\mathit{T}}^{\mathit{\nu}}}_{\mathrm{\left(}\mathit{\alpha}\mathrm{\right)}\mathit{\mu}}\mathrm{=}\mathrm{(}{\mathit{\rho}}_{\mathit{\alpha}}\mathrm{+}{\mathit{p}}_{\mathit{\alpha}}\mathrm{)}{\mathit{u}}_{\mathit{\mu}}{\mathit{u}}^{\mathit{\nu}}\mathrm{+}{\mathit{p}}_{\mathit{\alpha}}{\mathit{\delta}}_{\mathit{\mu}}^{\mathit{\nu}}$, where u_{μ} = (−a,0,0,0) is the fluid 4velocity and w_{α} ≡ p_{α}/ρ_{α} is the equation of state. Primes denote derivatives with respect to conformal time and β is assumed, for simplicity, to be a constant. This choice corresponds to a Lagrangian in which dark matter fields have an exponential mass dependence m(φ) = m_{0}exp^{− βφ} (originally motivated by Weyl scaling scalartensor theories), where m_{0} is a constant. The DE scalar field (expressed in units of the reduced Planck mass M = (8πG_{N})^{− 1 / 2}) evolves according to the KleinGordon equation, which now includes an extra term that depends on the density of cold dark matter: ${\mathit{\phi}}^{\mathrm{\prime \prime}}\mathrm{+}\mathrm{2}\mathrm{\mathscr{H}}{\mathit{\phi}}^{\mathrm{\prime}}\mathrm{+}{\mathit{a}}^{\mathrm{2}}\frac{\mathrm{d}\mathit{V}}{\mathrm{d}\mathit{\phi}}\mathrm{=}{\mathit{a}}^{\mathrm{2}}\mathit{\beta}{\mathit{\rho}}_{\mathrm{c}}\mathit{.}$(56)Following Pettorino & Baccigalupi (2008), we choose an inverse powerlaw potential defined as V = V_{0}φ^{− α}, with α and V_{0} being constants. The amplitude V_{0} is fixed thanks to an iterative routine, as implemented by Amendola et al. (2012), Pettorino et al. (2012). To a first approximation α only affects latetime cosmology. For numerical reasons, the iterative routine finds the initial value of the scalar field, in the range α ≥ 0.03, which is close to the ΛCDM value α = 0 and extends the range of validity with respect to past attempts; the equation of state w is approximately related to α via the expression (Amendola et al. 2012): w = −2 / (α + 2) so that a value of α = 0.03 corresponds approximately to w(α = 0.03) = −0.99. The equation of state w ≡ p/ρ is not an independent parameter within coupled DE theories, being degenerate with the flatness of the potential. Dark matter particles then feel a fifth force with an effective gravitational constant G_{eff} that is stronger than the Newtonian one by a factor of β^{2}, i.e. $\begin{array}{ccc}{\mathit{G}}_{\mathrm{eff}}\mathrm{=}\mathit{G}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{2}{\mathit{\beta}}^{\mathrm{2}}\mathrm{)}\mathit{,}& & \end{array}$(57)so that a value of β = 0 recovers the standard gravitational interaction. The coupling affects the dynamics of the gravitational potential (and therefore the late ISW effect), hence the shape and amplitude of perturbation growth, and shifts the position of the acoustic peaks to larger multipoles, due to an increase in the distance to the lastscattering surface; furthermore, it reduces the ratio of baryons to dark matter at early times with respect to its present value, since coupled dark matter dilutes faster than in an uncoupled model. The strength of the coupling is known to be degenerate with a combination of Ω_{c}, n_{s} and H(z) (Amendola & Quercellini 2003; Pettorino & Baccigalupi 2008; Bean et al. 2008; Amendola et al. 2012). Several analyses have previously been carried out, with hints of coupling different from zero, e.g., by Pettorino (2013), who found β = 0.036 ± 0.016 (using Planck 2013 + WMAP polarization + BAO) different from zero at 2.2σ (the significance increasing to 3.6σ when data from HST were included).
The marginalized posterior distribution, using Planck 2015 data, for the coupling parameter β is shown in Fig. 21, while the corresponding mean values are shown in Table 9. Planck TT data alone gives constraints compatible with zero coupling and the slope of the potential is consistent with a cosmological constant value of α = 0 at 1.3σ. When other data sets are added, however, both the value of the coupling and the slope of the potential are pushed to nonzero values, i.e., further from ΛCDM. In particular, Planck+BSH gives a value which is ~ 2.5σ in tension with ΛCDM, while, separately, Planck+WL+BAO/RSD gives a value of the coupling β compatible with the one from Planck+BSH and about 2.3σ away from ΛCDM.
Fig. 21 Marginalized posterior distribution for the coupling β (see Sect. 5.3.4). The value corresponding to standard gravity is zero. Results and goodness of fit are discussed in the text. 
When comparing with ΛCDM, however, the goodness of fit does not improve, despite the additional parameters. Only the ${\mathit{\chi}}_{\mathrm{BAO}\mathit{/}\mathrm{RSD}}^{\mathrm{2}}$ improves by ≈ 1 in CDE with respect to ΛCDM, the difference not being significant enough to justify the additional parameters. The fact that the marginalized likelihood does not improve, despite the apparent 2σ tension, may hint at some dependence on priors: for example, the first panel in Fig. 22 shows that there is some degeneracy between the coupling β and the potential slope α; while contours are almost compatible with ΛCDM in the 2 dimensional plot, the marginalization over α takes more contributions from higher values of β, due to the degeneracy, and seems to give a slight more significant peak in the one dimensional posterior distribution shown in Fig. 21. Whether other priors also contribute to the peak remains to be understood. In any case, the goodness of fit does not point towards a preference for nonzero coupling. Degeneracy between the coupling and other cosmological parameters is shown in the other panels of the same figure, with results compatible with those discussed in Amendola et al. (2012) and Pettorino (2013). Looking at the conservation equations (i.e., Eqs. (55) and (56)), larger positive values of β correspond to a larger transfer of energy from dark matter to DE (effectively adding more DE in the recent past, with roughly Ω_{φ} ∝ β^{2} for an inverse powerlaw potential) and therefore lead to a smaller Ω_{m} today; as a consequence, the distance to the lastscattering surface and the expansion rate are modified, with H′/H = −3 / 2(1 + w_{eff}), where w_{eff} is the effective equation of state given by the ratio of the total pressure over total (weighted) energy density of the coupled fluid; a larger coupling prefers larger H_{0} and higher σ_{8}.
The addition of polarization tightens the bounds on the coupling, increasing the tension with ΛCDM, reaching 2.8σ and 2.7σ for Planck+BSH and Planck+WL+BAO/RSD, respectively. Also in this case the overall χ^{2} does not improve between coupled DE and ΛCDM.
Fig. 22 Marginalized posterior distribution for coupled DE and different combinations of the data sets (see Sect. 5.3.4). Here Planck refers to Planck TT+lowP. We show the degeneracy of the coupling β with α, Ω_{m}, σ_{8} and H_{0}. 
6. Conclusions
The quest for Dark Energy and Modified Gravity is far from over. A variety of different theoretical scenarios have been proposed in literature and need to be carefully compared with the data. This effort is still in its early stages, given the variety of theories and parameterizations that have been suggested, together with a lack of well tested numerical codes that allow us to make detailed predictions for the desired range of parameters. In this paper, we have provided a systematic analysis covering a general survey of a variety of theoretical models, including the use of different numerical codes and observational data sets. Even though most of the weight in the Planck data lies at high redshift, Planck can still provide tight constraints on DE and MG, especially when used in combination with other probes. Our focus has been on the scales where linear theory is applicable, since these are the most theoretically robust. Overall, the constraints that we find are consistent with the simplest scenario, ΛCDM, with constraints on DE models (including minimallycoupled scalar field models or evolving equation of state models) and MG models (including effective field theory, phenomenological parameterizations, f(R) and coupled DE models) that are significantly improved with respect to past analyses. We discuss here our main results, drawing our conclusions for each of them and summarizing the storyline we have followed in this paper to discuss DE and MG.
Our journey started from distinguishing background and perturbation parameterizations. In the first case, the background is modified (which in turn affects the perturbations), leading to the following results.

1.
The equation of state w(z) as a function of redshift has been tested for a variety of parameterizations.

(a)
In (w_{0},w_{a}), Planck TT+lowP+BSH is compatible with ΛCDM, as well as BAO/RSD. When adding WL to Planck TT+lowP, both WL and CMB prefer the (w_{0},w_{a}) model with respect to ΛCDM at about 2σ, although with a preference for high values of H_{0} (third panel of Fig. 3) that are excluded when including BSH.

(b)
We have reconstructed the equation of state in redshift, testing a Taylor expansion up to the third order in the scale factor and by doing a Principal Component Analysis of w(z) in different redshifts bins. In addition, we have tested an alternative parametrization, that allows to have a varying w(z) that depends on one parameter only. All tests on time varying w(z) are compatible with ΛCDM for all data sets tested.

(a)
“Background” Dark Energy models are generally of quintessence type where a scalar field rolls down a potential. We have shown via the (ϵ_{s},ϵ_{∞}) parameterization, related respectively to late and early time evolution, that the quintessence/phantom potential at low redshift must be relatively flat: dlnV/ dφ < 0.35 /M_{P} for quintessence; and dlnV/ dφ < 0.68 /M_{P} for phantom models. A zero slope (ΛCDM) remains consistent with the data and compared to previous studies, the uncertainty has been reduced by about 10%. We have produced a new plot (Fig. 9) that helps to visualize minimally coupled scalar field models, similarly to analogous plots often used in inflationary theories.
Information on DE, complementary to (w_{0},w_{a}), comes from asking whether there can be any DE at early times. First, we have obtained constraints on early DE parametrizations, assuming a constant DE relative density at all epochs until it matches ΛCDM in recent times: we have improved previous constraints by a factor ~3−4, leading to Ω_{e} < 0.0071 at 95% C.L. from Planck TT+lowP+lensing+BSH and Ω_{e} < 0.0036 for Planck TT,TE,EE+lowP+BSH. In addition, we have also asked how much such tight constraints are weakened when the fraction of early DE is only present in a limited range in redshift and presented a plot of Ω_{e}(z) as a function of z_{e}, the redshift starting from which a fraction Ω_{e} is present. Also in this case constraints are very tight, with Ω_{e} ≲ 2% (95% C.L.) even for z_{e} as late as ≈50.
The background is then forced to be very close to ΛCDM, unless the tight constraints on early DE can somehow be evaded in a realistic model by counter balancing effects.
In the second part of the paper, we then moved on to understanding what Planck can say when the evolution of the perturbations is modified independently of the background, as is the case in most MG models. For that, we followed two complementary approaches: one (topdown) that starts from a very general EFT action for DE (Sect. 5.2.1); the other (bottomup) that starts from parameterizing directly observables (Sect. 5.2.2). In both cases we have assumed that the background is exactly ΛCDM, in order to disentangle the effect of perturbations. We summarize here our results.

1.
Starting from EFT theories for DE, which include (almost) alluniversally coupled models in MG via nine generic functions oftime, we have discussed how to restrict them to Horndeskitheories, described in terms of five free functions of time. Usingthe publicly available code EFTcamb, we have then varied three ofthese functions (in the limits allowed by the code) whichcorrespond to a nonminimally coupled Kessence model (i.e.α_{B}, α_{M}, α_{K} are varying functions of the scale factor). We have found limits on the present value α_{M}_{0}< 0.052 at 95% C.L. (in the linear EFT approximation), in agreement with ΛCDM. Constraints depend on the stability routines included in the code, which will need to be further tested in the future, together with allowing for a larger set of choices for the Horndeski functions, not available in the present version of the numerical code.

2.
When starting from observables, two functions of time and scale are required to describe perturbations completely, in any model. Among the choices available (summarized in Sect. 3.2.2), we choose μ(a,k) and η(a,k) (other observables can be derived from them). In principle, constraints on these functions are dependent on the chosen parameterization, which needs to be fixed. We have tested two different time dependent parameterizations (DErelated and timerelated) and both lead to similar results, although the first is slightly more in tension than the other with ΛCDM (Fig. 14). In this framework, ΛCDM lies at the 2σ limit when Planck TT+lowP+BSH is considered, the tension increasing to about 3σ when adding WL and BAO/RSD to Planck TT+lowP. As discussed in the text, the mild tension with Planck TT+lowP is related to lower power in the TT spectrum and a larger lensing potential in the MG model, with respect to ΛCDM. The inclusion of CMB lensing shifts all contours back to ΛCDM. We have reconstructed the two observables in redshift for both parameterizations, along the maximum degeneracy line (Fig. 16). When scale dependence is also included, constraints become much weaker and the goodness of fit does not improve, indicating that the data do not seem to need the addition of additional scale dependent parameters.
The last part of the paper discusses a selection of particular MG models of interest in literature.

1.
We first commented on the simple case of a minimally coupled scalar field in which not only the equation of state is allowed to vary but the sound speed of the DE fluid is not forced to be 1, as it would be in the case of quintessence. Such a scenario corresponds to kessence type models. As expected, given that the equation of state is very close to the ΛCDM value, the total impact of DE perturbations on the clustering is small.

2.
We adopt an alternative way to parameterize observables (the equation of state approach) in terms of gauge invariant quantities Γ and σ. We have used this approach to investigate Lorentzbreaking massive gravity and generalised scalar fields models, updating previous bounds.

3.
As a concrete example of universally coupled theories, we have considered f(R) models, written in terms of B(z), conventionally related to the first and second derivatives of f(R) with respect to R. Results are compatible with ΛCDM. Such theories assume that some screening mechanism is in place, in order to satisfy current bounds on baryonic physics at solar system scales.

4.
Alternatively to screening mechanisms, one can assume that the coupling is not universal, such that baryons are still feeling standard gravitational attraction. As an example of this scenario, we have considered the case in which the dark matter evolution is coupled to the DE scalar field, feeling an effective fifth force stronger than gravity by a factor β^{2}. Constraints on coupled dark energy show a tension with ΛCDM at the level of about 2.5σ, slightly increasing when including polarization. The apparent tension, however, seems to hint at a dependence on priors, partly related to the degeneracy between the coupling and the slope of the background potential (and possibly others not identified here). Future studies will need to identify the source of tension and possibly disentangle background from perturbation effects.
There are several ways in which the analysis can be extended. We have made an effort to (at least start) to put some order in the variety of theoretical frameworks discussed in literature. There are of course scenarios not included in this picture that deserve future attention, such as additional cosmologies within the EFT (and Horndeski) framework, Galileons (see for example Barreira et al. 2014), other massive gravity models (see de Rham 2014 for a recent review), general violations of Lorentz invariance as a way to modify GR (Audren et al. 2015), nonlocal gravity (which, for some choices of the action, appears to fit Planck 2013 data sets Dirian et al. 2015, as well as ΛCDM, although there is no connection to a fundamental theory available at present); models of bigravity (Hassan & Rosen 2012) appear to be affected by instabilities in the gravitational wave sector (Cusin et al. 2015 see also Akrami et al. 2015) and are not considered in this paper. In addition to extending the range of theories, which requires new numerical codes, future tests should verify whether all the assumptions (such as stability constraints, as pointed out in the text) in the currently available codes are justified. Further promising input may come from data sets such as WL and BAO/RSD, that allow to tighten considerably constraints on MG models in which perturbations are modified. With the data available at the time of this paper, there seems to be no significant trend (at more than 3 standard deviations) that compensates any possible tension in σ_{8} or H_{0} by favouring Modified Gravity; nevertheless, this issue will need to be further investigated in the future. We also anticipate that these constraints will strengthen with future releases of the Planck data, including improved likelihoods for polarization and new likelihoods, not available at the time of this paper, such as ISW, ISWlensing and Bmode polarization, all of which can be used to further test MG scenarios.
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).
A CosmoMC likelihood module for the JLA sample is available at http://supernovae.in2p3.fr/sdss_snls_jla/ReadMe.html
There are also alternative approaches that compress the power spectra directly, like e.g. PICO (Fendt & Wandelt 2006).
http://www.lorentz.leidenuniv.nl/~hu/codes/, version 1.1, Oct. 2014.
Because of the way EFTCAMB currently implements these equations internally, it is not possible to satisfy Eq. (30) otherwise.
The parameterization of μ and Σ in Simpson et al. (2013) uses Ω_{DE}(a) / Ω_{DE} instead of Ω_{DE}(a); their μ_{0} and Σ_{0} correspond to our μ_{0}−1 and Σ_{0}−1 respectively.
Available at http://www.sfu.ca/~aha25/MGCAMB.html (Feb. 2014 version), see appendix A of Zhao et al. (2009) for a detailed description of the implementation.
Acknowledgments
It is a pleasure to thank Luca Amendola, Emilio Bellini, Diego Blas, Sarah Bridle, Noemi Frusciante, Catherine Heymans, Alireza Hojjati, Bin Hu, Thomas Kitching, Niall MacCrann, Marco Raveri, Ignacy Sawicki, Alessandra Silvestri, Fergus Simpson, Christof Wetterich and Gongbo Zhao for interesting discussions on theories, external data sets and numerical codes. Part of the analysis for this paper was run on the Andromeda and Perseus clusters of the University of Geneva and on WestGrid computers in Canada. We deeply thank Andreas Malaspinas for invaluable help with the Andromeda and Perseus Clusters. This research was partly funded by the DFG TransRegio TRR33 grant on “The Dark Universe” and by the Swiss NSF. 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.
References
 Acquaviva, V., & Baccigalupi, C. 2006, Phys. Rev. D, 74, 103510 [NASA ADS] [CrossRef] [Google Scholar]
 Acquaviva, V., Baccigalupi, C., & Perrotta, F. 2004, Phys. Rev. D, 70, 023515 [NASA ADS] [CrossRef] [Google Scholar]
 Akrami, Y., Hassan, S. F., Könnig, F., SchmidtMay, A., & Solomon, A. R. 2015, Phys. Lett. B, 748, 37 [NASA ADS] [CrossRef] [Google Scholar]
 Alcock, C., & Paczynski, B. 1979, Nature, 281, 358 [NASA ADS] [CrossRef] [Google Scholar]
 Amendola, L. 2000, Phys. Rev. D, 62, 043511 [CrossRef] [Google Scholar]
 Amendola, L., & Quercellini, C. 2003, Phys. Rev. D, 68, 023514 [NASA ADS] [CrossRef] [Google Scholar]
 Amendola, L., Gannouji, R., Polarski, D., & Tsujikawa, S. 2007, Phys. Rev. D, 75, 083504 [NASA ADS] [CrossRef] [Google Scholar]
 Amendola, L., Baldi, M., & Wetterich, C. 2008a, Phys. Rev. D, 78, 023015 [NASA ADS] [CrossRef] [Google Scholar]
 Amendola, L., Kunz, M., & Sapone, D. 2008b, J. Cosmol. Astropart. Phys., 0804, 013 [Google Scholar]
 Amendola, L., Pettorino, V., Quercellini, C., & Vollmer, A. 2012, Phys. Rev. D, 85, 103008 [NASA ADS] [CrossRef] [Google Scholar]
 Amendola, L., Appleby, S., Bacon, D., et al. 2013, Liv. Rev. Rel., 16, 6 [Google Scholar]
 Amendola, L., Ballesteros, G., & Pettorino, V. 2014, Phys. Rev. D, 90, 043009 [NASA ADS] [CrossRef] [Google Scholar]
 Anderson, L., Aubourg, É., Bailey, S., et al. 2014, MNRAS, 441, 24 [NASA ADS] [CrossRef] [Google Scholar]
 Antolini, C., Martinelli, M., Fantaye, Y., & Baccigalupi, C. 2013, J. Cosmol. Astropart. Phys., 2, 024 [NASA ADS] [CrossRef] [Google Scholar]
 ArmendarizPicon, C., Mukhanov, V. F., & Steinhardt, P. J. 2000, Phys. Rev. Lett., 85, 4438 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Aubourg, É., Bailey, S., Bautista, J. E., et al. 2015, Phys. Rev. D, 92, 123516 [NASA ADS] [CrossRef] [Google Scholar]
 Audren, B., Blas, D., Ivanov, M., Lesgourgues, J., & Sibiryakov, S. 2015, JCAP, 03, 016 [NASA ADS] [CrossRef] [Google Scholar]
 Baker, T., Ferreira, P., & Skordis, C. 2014a, Phys. Rev. D, 89, 024026 [NASA ADS] [CrossRef] [Google Scholar]
 Baker, T., Ferreira, P. G., Leonard, C. D., & Motta, M. 2014b, Phys. Rev. D, 90, 124030 [NASA ADS] [CrossRef] [Google Scholar]
 Baldi, M., & Pettorino, V. 2011, MNRAS, 412, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Ballesteros, G., Hollenstein, L., Jain, R. K., & Kunz, M. 2012, J. Cosmol. Astropart. Phys., 1205, 038 [NASA ADS] [CrossRef] [Google Scholar]
 Bardeen, J. 1980, Phys. Rev. D, 22, 1882 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Barreira, A., Li, B., Baugh, C. M., & Pascoli, S. 2014, J. Cosmol. Astropart. Phys., 8, 059 [Google Scholar]
 Barrow, J. D., & Saich, P. 1993, MNRAS, 262, 717 [NASA ADS] [CrossRef] [Google Scholar]
 Battye, R. A., & Pearson, J. A. 2012, J. Cosmol. Astropart. Phys., 1207, 019 [NASA ADS] [CrossRef] [Google Scholar]
 Battye, R. A., & Pearson, J. A. 2013, Phys. Rev. D, 88, 084004 [NASA ADS] [CrossRef] [Google Scholar]
 Battye, R. A., Moss, A., & Pearson, J. A. 2015, J. Cosmol. Astropart. Phys., 4, 048 [NASA ADS] [CrossRef] [Google Scholar]
 Bean, R., Flanagan, E. E., Laszlo, I., & Trodden, M. 2008, Phys. Rev. D, 78, 123514 [NASA ADS] [CrossRef] [Google Scholar]
 Bellini, E., & Sawicki, I. 2014, J. Cosmol. Astropart. Phys., 7, 50 [Google Scholar]
 Bennett, C. L., Larson, D., Weiland, J. L., et al. 2013, ApJS, 208, 20 [Google Scholar]
 Bertschinger, E., & Zukin, P. 2008, Phys. Rev. D, 78, 024015 [NASA ADS] [CrossRef] [Google Scholar]
 Betoule, M., Marriner, J., Regnault, N., et al. 2013, A&A, 552, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Betoule, M., Kessler, R., Guy, J., et al. 2014, A&A, 568, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Beutler, F., Blake, C., Colless, M., et al. 2011, MNRAS, 416, 3017 [NASA ADS] [CrossRef] [Google Scholar]
 Beutler, F., Blake, C., Colless, M., et al. 2012, MNRAS, 423, 3430 [NASA ADS] [CrossRef] [Google Scholar]
 Beutler, F., Saito, S., Seo, H.J., et al. 2014, MNRAS, 443, 1065 [NASA ADS] [CrossRef] [Google Scholar]
 Bianchi, D., Guzzo, L., Branchini, E., et al. 2012, MNRAS, 427, 2420 [NASA ADS] [CrossRef] [Google Scholar]
 Bird, S., Viel, M., & Haehnelt, M. G. 2012, MNRAS, 420, 2551 [NASA ADS] [CrossRef] [Google Scholar]
 Blake, C., Kazin, E. A., Beutler, F., et al. 2011, MNRAS, 418, 1707 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Bloomfield, J., Flanagan, É. É., Park, M., & Watson, S. 2013, J. Cosmol. Astropart. Phys., 8, 10 [Google Scholar]
 Brax, P. H., & Martin, J. 1999, Phys. Lett. B, 468, 40 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Bridle, S., & King, L. 2007, New J. Phys., 9, 444 [NASA ADS] [CrossRef] [Google Scholar]
 Calabrese, E., Slosar, A., Melchiorri, A., Smoot, G. F., & Zahn, O. 2008, Phys. Rev. D, 77, 123531 [NASA ADS] [CrossRef] [Google Scholar]
 Calabrese, E., Huterer, D., Linder, E. V., Melchiorri, A., & Pagano, L. 2011, Phys. Rev. D, 83, 123504 [NASA ADS] [CrossRef] [Google Scholar]
 Capozziello, S. 2002, Int. J. Mod. Phys. D, 11, 483 [Google Scholar]
 Carbone, C., Baldi, M., Pettorino, V., & Baccigalupi, C. 2013, J. Cosmol. Astropart. Phys., 9, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Cataneo, M., Rapetti, D., Schmidt, F., et al. 2015, Phys. Rev. D, 92, 044009 [NASA ADS] [CrossRef] [Google Scholar]
 Cheung, C., Fitzpatrick, A. L., Kaplan, J., Senatore, L., & Creminelli, P. 2008, J. High Energy Phys., 3, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Clifton, T., Ferreira, P. G., Padilla, A., & Skordis, C. 2012, Phys. Rep., 513, 1 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Creminelli, P., D’Amico, G., Noreña, J., & Vernizzi, F. 2009, J. Cosmol. Astropart. Phys., 2, 018 [Google Scholar]
 Cusin, G., Durrer, R., Guarato, P., & Motta, M. 2015, J. Cosmol. Astropart. Phys., 05, 030 [NASA ADS] [CrossRef] [Google Scholar]
 Daniel, S. F., Linder, E. V., Smith, T. L., et al. 2010, Phys. Rev. D, 81, 123508 [NASA ADS] [CrossRef] [Google Scholar]
 De Felice, A., & Tsujikawa, S. 2010, Liv. Rev. Rel., 13, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Deffayet, C., Pujolas, O., Sawicki, I., & Vikman, A. 2010, J. Cosmol. Astropart. Phys., 1010, 026 [Google Scholar]
 de Rham, C. 2014, Liv. Rev. Rel., 17, 7 [Google Scholar]
 Dirian, Y., Foffa, S., Kunz, M., Maggiore, M., & Pettorino, V. 2015, J. Cosmol. Astropart. Phys., 04, 044 [NASA ADS] [CrossRef] [Google Scholar]
 Doran, M., & Robbers, G. 2006, J. Cosmol. Astropart. Phys., 06, 026 [Google Scholar]
 Dossett, J., Hu, B., & Parkinson, D. 2014, J. Cosmol. Astropart. Phys., 3, 046 [NASA ADS] [CrossRef] [Google Scholar]
 Dubovsky, S. L. 2004, J. High Energy Phys., 10, 76 [Google Scholar]
 Efstathiou, G. 2014, MNRAS, 440, 1138 [NASA ADS] [CrossRef] [Google Scholar]
 Efstathiou, G., & Bond, J. 1999, MNRAS, 304, 75 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Einstein, A. 1917, Sitzungsberichte der Königlich Preußischen Akademie der Wissenschaften (Berlin), 142 [Google Scholar]
 Fang, W., Hu, W., & Lewis, A. 2008, Phys. Rev. D, 78, 087303 [NASA ADS] [CrossRef] [Google Scholar]
 Fardon, R., Nelson, A. E., & Weiner, N. 2004, J. Cosmol. Astropart. Phys., 0410, 005 [Google Scholar]
 Fendt, W. A., & Wandelt, B. D. 2006, ApJ, 654, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Frieman, J. A., Hill, C. T., Stebbins, A., & Waga, I. 1995, Phys. Rev. Lett., 75, 2077 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Giannantonio, T., Scranton, R., Crittenden, R. G., et al. 2008, Phys. Rev., D, 77, 123520 [Google Scholar]
 GilMarín, H., Wagner, C., Verde, L., Porciani, C., & Jimenez, R. 2012, J. Cosmol. Astropart. Phys., 11, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Gleyzes, J., Langlois, D., Piazza, F., & Vernizzi, F. 2013, J. Cosmol. Astropart. Phys., 8, 025 [Google Scholar]
 Gleyzes, J., Langlois, D., Piazza, F., & Vernizzi, F. 2015, JCAP, 02, 018 [NASA ADS] [CrossRef] [Google Scholar]
 Gott, J. R., & Slepian, Z. 2011, MNRAS, 416, 907 [NASA ADS] [CrossRef] [Google Scholar]
 Gubitosi, G., Piazza, F., & Vernizzi, F. 2013, J. Cosmol. Astropart. Phys., 1302, 032 [Google Scholar]
 HarnoisDéraps, J., van Waerbeke, L., Viola, M., & Heymans, C. 2015, MNRAS, 450, 1212 [NASA ADS] [CrossRef] [Google Scholar]
 Hassan, S. F., & Rosen, R. A. 2012, J. High Energy Phys., 2, 126 [NASA ADS] [CrossRef] [Google Scholar]
 Heymans, C., Grocutt, E., Heavens, A., et al. 2013, MNRAS, 432, 2433 [NASA ADS] [CrossRef] [Google Scholar]
 Hojjati, A., Pogosian, L., & Zhao, G.B. 2011, J. Cosmol. Astropart. Phys., 1108, 005 [Google Scholar]
 Horndeski, G. W. 1974, Int. J. Theor. Phys., 10, 363 [CrossRef] [MathSciNet] [Google Scholar]
 Howlett, C., Ross, A., Samushia, L., Percival, W., & Manera, M. 2015, MNRAS, 449, 848 [NASA ADS] [CrossRef] [Google Scholar]
 Hu, B., Liguori, M., Bartolo, N., & Matarrese, S. 2013, Phys. Rev. D, 88, 123514 [NASA ADS] [CrossRef] [Google Scholar]
 Hu, B., Raveri, M., Frusciante, N., & Silvestri, A. 2014a, Phys. Rev. D, 89, 103530 [NASA ADS] [CrossRef] [Google Scholar]
 Hu, B., Raveri, M., Frusciante, N., & Silvestri, A. 2014b, ArXiv eprints [arXiv:1405.3590] [Google Scholar]
 Hu, W. 2002, Phys. Rev. D, 65, 023003 [NASA ADS] [CrossRef] [Google Scholar]
 Hu, W., & Sawicki, I. 2007, Phys. Rev. D, 76, 104043 [NASA ADS] [CrossRef] [Google Scholar]
 Hu, W., & White, M. 1996, ApJ, 471, 30 [NASA ADS] [CrossRef] [Google Scholar]
 Huang, Z., Bond, J. R., & Kofman, L. 2011, ApJ, 726, 64 [NASA ADS] [CrossRef] [Google Scholar]
 Hubble, E. 1929, Proc. Nat. Acad. Sci., 15, 168 [NASA ADS] [CrossRef] [Google Scholar]
 Humphreys, L., Reid, M., Moran, J., Greenhill, L., & Argon, A. 2013, ApJ, 775, 13 [Google Scholar]
 Huterer, D., & Cooray, A. 2005, Phys. Rev. D, 71, 023506 [NASA ADS] [CrossRef] [Google Scholar]
 Huterer, D., & Starkman, G. 2003, Phys. Rev. Lett., 90, 031301 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Huterer, D., Kirkby, D., Bean, R., et al. 2015, Astropart. Phys., 63, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Joyce, A., Jain, B., Khoury, J., & Trodden, M. 2015, Phys. Rep., 568, 1 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Kaloper, N., & Sorbo, L. 2006, J. Cosmol. AstroPart. Phys., 4, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Khoury, J. 2010, Conference Proceedings, 22nd Rencontres de Blois on Particle Physics and Cosmology [arXiv:1011.5909] [Google Scholar]
 Kilbinger, M., Fu, L., Heymans, C., et al. 2013, MNRAS, 430, 2200 [NASA ADS] [CrossRef] [Google Scholar]
 Kitching, T., Heavens, A. F., Alsing, J., et al. 2014, MNRAS, 442, 1326 [NASA ADS] [CrossRef] [Google Scholar]
 Kodama, H., & Sasaki, M. 1984, Prog. Theor. Phys. Suppl., 78, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Kofman, L. A., & Starobinskii, A. A. 1985, Sov. Astron. Lett., 11, 271 [NASA ADS] [Google Scholar]
 Kosowsky, A., Milosavljevic, M., & Jimenez, R. 2002, Phys. Rev. D, 66, 063007 [NASA ADS] [CrossRef] [Google Scholar]
 Kunz, M. 2009, Phys. Rev. D, 80, 123001 [NASA ADS] [CrossRef] [Google Scholar]
 Kunz, M. 2012, Comptes Rendus Physique, 13, 539 [NASA ADS] [CrossRef] [Google Scholar]
 Kunz, M., Corasaniti, P.S., Parkinson, D., & Copeland, E. J. 2004, Phys. Rev. D, 70, 041301 [NASA ADS] [CrossRef] [Google Scholar]
 Lemaître, G. 1927, Annales de la Société Scientifique de Bruxelles, 47, 49 [NASA ADS] [Google Scholar]
 Lewis, A., & Bridle, S. 2002, Phys. Rev. D, 66, 103511 [NASA ADS] [CrossRef] [Google Scholar]
 Lewis, A., & Challinor, A. 2006, Phys. Rep., 429, 1 [Google Scholar]
 Linde, A. D. 1983, Phys. Lett. B, 129, 177 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 LSST Science Collaboration, Abell, P. A., Allison, J., et al. 2009, ArXiv eprints [arXiv:0912.0201] [Google Scholar]
 Ma, C.P., & Bertschinger, E. 1995, ApJ, 455, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Macaulay, E., Wehus, I. K., & Eriksen, H. K. 2013, Phys. Rev. Lett., 111, 161301 [NASA ADS] [CrossRef] [Google Scholar]
 MacCrann, N., Zuntz, J., Bridle, S., Jain, B., & Becker, M. R. 2015, MNRAS, 451, 2877 [NASA ADS] [CrossRef] [Google Scholar]
 Marchini, A., & Salvatelli, V. 2013, Phys. Rev. D, 88, 027502 [NASA ADS] [CrossRef] [Google Scholar]
 Marchini, A., Melchiorri, A., Salvatelli, V., & Pagano, L. 2013, Phys. Rev. D, 87, 083527 [NASA ADS] [CrossRef] [Google Scholar]
 Mukhanov, V. F., Feldman, H., & Brandenberger, R. H. 1992, Phys. Rep., 215, 203 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Mukherjee, P., Kunz, M., Parkinson, D., & Wang, Y. 2008, Phys. Rev. D, 78, 083529 [NASA ADS] [CrossRef] [Google Scholar]
 Pan, Z., Knox, L., & White, M. 2014, MNRAS, 445, 2941 [NASA ADS] [CrossRef] [Google Scholar]
 Peebles, P. J. E. 1984, ApJ, 284, 439 [NASA ADS] [CrossRef] [Google Scholar]
 Perlmutter, S., Aldering, G., Goldhaber, G., et al. 1999, ApJ, 517, 565 [NASA ADS] [CrossRef] [Google Scholar]
 Pettorino, V. 2013, Phys. Rev. D, 88, 063519 [NASA ADS] [CrossRef] [Google Scholar]
 Pettorino, V., & Amendola, L. 2015, Phys. Lett. B, 742, 353 [NASA ADS] [CrossRef] [Google Scholar]
 Pettorino, V., & Baccigalupi, C. 2008, Phys. Rev. D, 77, 103003 [NASA ADS] [CrossRef] [Google Scholar]
 Pettorino, V., Amendola, L., Baccigalupi, C., & Quercellini, C. 2012, Phys. Rev. D, 86, 103507 [NASA ADS] [CrossRef] [Google Scholar]
 Pettorino, V., Amendola, L., & Wetterich, C. 2013, Phys. Rev. D, 87, 083009 [NASA ADS] [CrossRef] [Google Scholar]
 Piazza, F., Steigerwald, H., & Marinoni, C. 2014, J. Cosmol. Astropart. Phys., 5, 43 [Google Scholar]
 Planck Collaboration 2005, ESA publication ESASCI(2005)/01 [arXiv:astroph/0604069] [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 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. XXX. 2016, A&A, 586, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pogosian, L., & Silvestri, A. 2008, Phys. Rev. D, 77, 023503 [NASA ADS] [CrossRef] [Google Scholar]
 Ratra, B., & Peebles, P. J. E. 1988, Phys. Rev. D, 37, 3406 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Raveri, M., Hu, B., Frusciante, N., & Silvestri, A. 2014, Phys. Rev. D, 90, 043513 [NASA ADS] [CrossRef] [Google Scholar]
 Raveri, M., Baccigalupi, C., Silvestri, A., & Zhou, S.Y. 2015, Phys. Rev. D, 91, 061501 [NASA ADS] [CrossRef] [Google Scholar]
 Reichardt, C. L., de Putter, R., Zahn, O., & Hou, Z. 2012, ApJ, 749, L9 [NASA ADS] [CrossRef] [Google Scholar]
 Reid, M. J., Braatz, J. A., Condon, J. J., et al. 2013, ApJ, 767, 154 [NASA ADS] [CrossRef] [Google Scholar]
 Riess, A. G., Filippenko, A. V., Challis, P., et al. 1998, AJ, 116, 1009 [NASA ADS] [CrossRef] [Google Scholar]
 Riess, A. G., Macri, L., Casertano, S., et al. 2011, ApJ, 730, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Ross, A. J., Samushia, L., Howlett, C., et al. 2015, MNRAS, 449, 835 [NASA ADS] [CrossRef] [Google Scholar]
 Rubakov, V. A., & Tinyakov, P. G. 2008, Phys. Uspekhi, 51, 759 [Google Scholar]
 Sachs, R. K., & Wolfe, A. M. 1967, ApJ, 147, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Said, N., Baccigalupi, C., Martinelli, M., Melchiorri, A., & Silvestri, A. 2013, Phys. Rev. D, 88, 043515 [NASA ADS] [CrossRef] [Google Scholar]
 Saltas, I. D., Sawicki, I., Amendola, L., & Kunz, M. 2014, Phys. Rev. Lett., 113, 191101 [NASA ADS] [CrossRef] [Google Scholar]
 Samushia, L., Reid, B. A., White, M., et al. 2014, MNRAS, 439, 3504 [NASA ADS] [CrossRef] [Google Scholar]
 Sapone, D., & Kunz, M. 2009, Phys. Rev. D, 80, 083519 [NASA ADS] [CrossRef] [Google Scholar]
 Sawicki, I., Saltas, I. D., Amendola, L., & Kunz, M. 2013, J. Cosmol. Astropart. Phys., 1301, 004 [NASA ADS] [CrossRef] [Google Scholar]
 Schmidt, F. 2008, Phys. Rev. D, 78, 043002 [NASA ADS] [CrossRef] [Google Scholar]
 Simpson, F., Heymans, C., Parkinson, D., et al. 2013, MNRAS, 429, 2249 [NASA ADS] [CrossRef] [Google Scholar]
 Slepian, Z., Gott, J. R., & Zinn, J. 2014, MNRAS, 438, 1948 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, R., Peacock, J. A., Jenkins, A., et al. 2003, MNRAS, 341, 1311 [NASA ADS] [CrossRef] [Google Scholar]
 Soergel, B., Giannantonio, T., Weller, J., & Battye, R. A. 2015, J. Cosmol. Astropart. Phys., 2, 037 [NASA ADS] [CrossRef] [Google Scholar]
 Song, Y.S., Hu, W., & Sawicki, I. 2007a, Phys. Rev. D, 75, 044004 [NASA ADS] [CrossRef] [Google Scholar]
 Song, Y.S., Peiris, H., & Hu, W. 2007b, Phys. Rev. D, 76, 063517 [NASA ADS] [CrossRef] [Google Scholar]
 Takahashi, R., Sato, M., Nishimichi, T., Taruya, A., & Oguri, M. 2012, ApJ, 761, 152 [NASA ADS] [CrossRef] [Google Scholar]
 van Daalen, M. P., Schaye, J., Booth, C., & Vecchia, C. D. 2011, MNRAS, 415, 3649 [NASA ADS] [CrossRef] [Google Scholar]
 Verde, L., & Spergel, D. N. 2002, Phys. Rev. D, 65, 043007 [NASA ADS] [CrossRef] [Google Scholar]
 Vikram, V., Sakstein, J., Davis, C., & Neil, A. 2014, ArXiv eprints [arXiv:1407.6044] [Google Scholar]
 Wang, Y., & Mukherjee, P. 2007, Phys. Rev. D, 76, 103533 [NASA ADS] [CrossRef] [Google Scholar]
 Wetterich, C. 1988, Nucl. Phys. B, 302, 668 [NASA ADS] [CrossRef] [Google Scholar]
 Wetterich, C. 1995, A&A, 301, 321 [NASA ADS] [Google Scholar]
 Wetterich, C. 2004, Phys. Lett. B, 594, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Wetterich, C. 2015, Lect. Notes Phys., 892, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, P., Liguori, M., Bean, R., & Dodelson, S. 2007, Phys. Rev. Lett., 99, 141302 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Zhao, G.B., Pogosian, L., Silvestri, A., & Zylberberg, J. 2009, Phys. Rev. D, 79, 083513 [NASA ADS] [CrossRef] [Google Scholar]
 Zhao, G.B., Giannantonio, T., Pogosian, L., et al. 2010, Phys. Rev. D, 81, 103510 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Marginalized mean values and 68% CL intervals for the EFT parameters, both in the linear model, α_{M0}, and in the exponential one, { α_{M0},β } (see Sect. 5.2.1).
Marginalized mean values and 68% C.L. errors on the present day value of the function 2 [μ(z,k)−1] + [η(z,k)−1], which corresponds to the (approximate) maximum degeneracy line identified within the 2 dimensional posterior distributions.
All Figures
Fig. 1 Typical effects of modified gravity on theoretical CMB temperature (top panel) and lensing potential (bottom panel) power spectra. An increase (or decrease) of E_{22} with respect to zero introduces a gravitational slip, higher at present, when Ω_{de} is higher (see Eqs. (4) and (6)); this in turns changes the Weyl potential and leads to a higher (or lower) lensing potential. On the other hand, whenever E_{11} and E_{22} are different from zero (quite independently of their sign) μ and η change in time: as the dynamics in the gravitational potential is increased, this leads to an enhancement in the ISW effect. Note also that even when the temperature spectrum is very close to ΛCDM (as for E_{11} = E_{22} = 0.5) the lensing potential is still different with respect to ΛCDM, shown in black. 

In the text 
Fig. 2 Ω_{m}–σ_{8} constraints for tomographic lensing from Heymans et al. (2013), using a very conservative angular cut, as described in the text (see Sect. 4.3.2). We show results using linear theory, nonlinear corrections from Halofit (HL) versions 1, 4, marginalization over baryonic AGN feedback (BF), and intrinsic alignment (IA) (the latter two using nonlinear corrections and Halofit 4). Coloured points indicate H_{0} values from WL+HL4. 

In the text 
Fig. 3 Parameterization { w_{0},w_{a} } (see Sect. 5.1.1). Marginalized posterior distributions for w_{0}, w_{a}, H_{0} and σ_{8} for various data combinations. The tightest constraints come from the Planck TT+lowP+BSH combination, which indeed tests background observations, and is compatible with ΛCDM. 

In the text 
Fig. 4 Marginalized posterior distributions of the (w_{0},w_{a}) parameterization (see Sect. 5.1.1) for various data combinations. The best constraints come from the priority combination and are compatible with ΛCDM. The dashed lines indicate the point in parameter space (− 1,0) corresponding to the ΛCDM model. CMB lensing and polarization do not significantly change the constraints. Here Planck indicates Planck TT+lowP. 

In the text 
Fig. 5 Reconstructed equation of state w(z) as a function of redshift (see Sect. 5.1.1), when assuming a Taylor expansion of w(z) to firstorder (N = 1 in Eq. (19)), for different combinations of the data sets. The coloured areas show the regions which contain 95% of the models. The central blue line is the median line for Planck TT+lowP+BSH. Here Planck indicates Planck TT+lowP. 

In the text 
Fig. 6 Marginalized posterior contours in the h–δw_{0} plane are shown for 1parameter varying w models (see Sect. 5.1.2) for different data combinations. Here Planck indicates Planck TT+lowP. 

In the text 
Fig. 7 PCA analysis constraints (described in Sect. 5.1.3). The top panel shows the reconstructed equation of state w(z) after the PCA analysis. Vertical error bars correspond to mean and standard deviations of the q vector parameters, while horizontal error bars are the amplitude of the original binning. The bins are not exactly independent but are rather smeared out as illustrated in the bottom panel. The bottom panel shows the PCA corresponding weights on w(z) as a function of redshift for the combination Planck TT+lowP+BSH. In other words, error bars in the top panel correspond therefore to the errors in the q parameters, which are linear combinations of the p parameters, i.e. a smeared out distribution with weights shown in the lower panel. 

In the text 
Fig. 8 Marginalized posterior distributions showing 68% and 95% C.L. constraints on Ω_{m} and ϵ_{s} for scalar field models (see Sect. 5.1.4). The dashed line for ϵ_{s} = 0 is the ΛCDM model. Here Planck indicates Planck TT+lowP. 

In the text 
Fig. 9 Marginalized posterior distributions at 68% C.L. and 95% C.L. in the parameter space of ϵ_{s} and ϵ_{∞} for scalar field modes (see Sect. 5.1.4). We have computed ϵ_{s} and ϵ_{∞} for various quintessence potentials V(φ), with the functional forms of V(φ) labelled on the figure. The field φ is in reduced Planck mass M_{P} units. The normalization of V(φ) is computed using Ω_{m} = 0.3. Here Planck indicates Planck TT+lowP. 

In the text 
Fig. 10 Marginalized posterior distributions for Ω_{e} for the early DE parameterization in Eq. (26) and for different combinations of data (see Sect. 5.1.5). Here Planck indicates Planck TT+lowP. 

In the text 
Fig. 11 Amount of DE at early times Ω_{e} as a function of the redshift z_{e} after which early DE is nonnegligible (see Eq. (27), Sect. 5.1.5) for different combinations of data sets. The heights of the columns give the limit at 95% CL on Ω_{e}, as obtained from Monte Carlo runs for the values z_{e} = 10, 50, 200 and 1000. The width of the columns has no physical meaning and is just due to plotting purposes. Here Planck indicates Planck TT+lowP. 

In the text 
Fig. 12 Marginalized posterior distributions at 68% and 95% C.L. for the two parameters α_{M0} and β of the exponential evolution, Ω(a) = exp(Ω_{0}a^{β})−1.0, see Sect. 5.2.1. Here α_{M0} is defined as Ω_{0}β and the background is fixed to ΛCDM. Ω_{M0} = 0 corresponds to the ΛCDM model also at perturbation level. Note that Planck means Planck TT+lowP. Adding WL to the data sets results in broader contours, as a consequence of the slight tension between the Planck and WL data sets. 

In the text 
Fig. 13 Marginalized posterior distribution of the linear EFT model background parameter, Ω, with Ω parameterized as a linear function of the scale factor, i.e., Ω(a) = α_{M0}a, see Sect. 5.2.1. The equation of state parameter w_{de} is fixed to −1, and therefore, Ω_{0} = 0 will correspond to the ΛCDM model. Here Planck means Planck TT+lowP. Adding CMB lensing to the data sets does not change the results significantly; highℓ polarization tightens the constraints by a few percent, as shown in Table 5. 

In the text 
Fig. 14 68% and 95% contour plots for the two parameters { μ_{0}−1,η_{0}−1 } obtained by evaluating Eqs. (46) and (47) at the present time when no scale dependence is considered (see Sect. 5.2.2). We consider both the DErelated (left panel) and timerelated evolution cases (right panel). Results are shown for the scaleindependent case (c_{1} = c_{2} = 1). In the labels, Planck stands for Planck TT+lowP. 

In the text 
Fig. 15 Marginalized posterior distributions for 68% and 95% C.L. for the two parameters {μ_{0}−1,Σ_{0}−1} obtained by evaluating Eqs. (46) and (47) at the present time in the DErelated parametrization when no scale dependence is considered (see Sect. 5.2.2). Σ is obtained as Σ = (μ/ 2)(1 + η). The timerelated evolution would give similar contours. In the labels, Planck stands for Planck TT+lowP. 

In the text 
Fig. 16 Redshift dependence of the function 2 [μ(z,k)−1] + [η(z,k)−1], defined in Eqs. (46) and (47), which corresponds to the maximum degeneracy line identified within the 2 dimensional posterior distributions. This combination shows the strongest allowed tension with ΛCDM. The left panel refers to the DErelated case while the right panel refers to the timerelated evolution (see Sect. 5.2.2). In both panels, no scale dependence is considered. The coloured areas show the regions containing 68% and 95% of the models. In the labels, Planck stands for Planck TT+lowP. 

In the text 
Fig. 17 68% and 95% marginalised posterior distributions for the two parameters { μ_{0}−1,η_{0}−1 } obtained by evaluating Eqs. (46) and (47) at the present time when no scale dependence is considered (see Sect. 5.2.2). Here we show the effect of CMB lensing, which shifts the contours towards ΛCDM. In the labels, Planck stands for Planck TT+lowP. 

In the text 
Fig. 18 68% and 95% contour plots for the two parameters {μ_{0}(k)−1, η_{0}(k)−1} obtained by evaluating Eqs. (46) and (47) at the present time for the DErelated parameterization (see Sect. 5.2.2). We consider both the scaleindependent and scaledependent cases, choosing k values of 10^{10} Mpc^{1} and 10^{2} Mpc^{1}. 

In the text 
Fig. 19 68% and 95% contour plots for the two parameters, { Log_{10}(B_{0}),τ } (see Sect. 5.3.3). There is a degeneracy between the two parameters for Planck TT+lowP+BSH. Adding lensing will break the degeneracy between the two. Here Planck indicates Planck TT+lowP. 

In the text 
Fig. 20 Likelihood plots of the f(R) theory parameter, B_{0} (see Sect. 5.3.3). CMB lensing breaks the degeneracy between B_{0} and the optical depth, τ, resulting in lower upper bounds. 

In the text 
Fig. 21 Marginalized posterior distribution for the coupling β (see Sect. 5.3.4). The value corresponding to standard gravity is zero. Results and goodness of fit are discussed in the text. 

In the text 
Fig. 22 Marginalized posterior distribution for coupled DE and different combinations of the data sets (see Sect. 5.3.4). Here Planck refers to Planck TT+lowP. We show the degeneracy of the coupling β with α, Ω_{m}, σ_{8} and H_{0}. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.