© 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 non-Gaussianities (Planck Collaboration XVII 2016). In this paper we use the 2015 data release from Planck¹ (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 non-standard 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 three-part approach.

• 1.

  Background parameterizations. Here we consider onlyparameterizations of background-level cosmologicalquantities. Perturbations are always included, but theirevolution depends only on the background. This set in-cludes 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 last-scattering 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 low-multipole CMB anisotropies through the integrated Sachs-Wolfe (ISW) effect (Sachs & Wolfe 1967; Kofman & Starobinskii 1985);

• 3.

  enhance the cross-correlation between the CMB and large-scale 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 CMB-inferred amplitude of the fluctuations A_s and late-time measurements of σ₈ (Kunz et al. 2004; Baldi & Pettorino 2011);

• 6.

  impact small scales, modifying the damping tail in , 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 B-mode contribution, through changes in the lensing potential (Amendola et al. 2014);

• 9.

  modify the primordial B-mode 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 scale-independent 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₁₁ = E₂₂ = 0 and μ = η = 1.

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 scalar-type 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 energy-momentum 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 (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 . The anisotropic stress σ(k,a) (equivalent to π_T in the notation of Kodama & Sasaki 1984) is the scalar part of the off-diagonal space-space stress energy tensor perturbation. The set of functions { H,Φ,Ψ } describing the metric is formally equivalent to the set of functions (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 minimally-coupled scalar field model (Wetterich 1988; Ratra & Peebles 1988), also known as quintessence, which corresponds to the choice of a rest-frame sound speed (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 sub-horizon scales, preventing it from contributing significantly to clustering.

Background parameterizations discussed in this paper include:

• (w₀, 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 Fisher-forecasts for future large-scale 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 full-scale 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: (2)Here R is the Ricci scalar, δR⁽³⁾ is its spatial perturbation, is the extrinsic curvature, and m₀ 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 time-dependent functions (Bloomfield et al. 2013), here , 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 scalar-tensor theories, universally coupled to gravity, with second-order 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 non-vanishing 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(3)where Δ is the comoving density perturbation;

• 2.

  μ(a,k) (sometimes also called Y(a,k)), which modifies the equivalent equation for Ψ rather than Φ: (4)

• 3.

  Σ(a,k), which modifies lensing (with the lensing/Weyl potential being Φ + Ψ), such that (5)

• 4.

  η(a,k), which reflects the presence of a non-zero anisotropic stress, the difference between Φ and Ψ being equivalently written as a deviation of the ratio²(6)

In the equations above, ρΔ = ρ_mΔ_m + ρ_rΔ_r so that the parameters Q, μ, or Σ quantify the deviation of the gravitational potentials from the value expected in GR due to perturbations of matter and relativistic particles. At low redshifts, where most DE models become relevant, we can neglect the relativistic contribution. The same is true for η, where we can neglect the contribution of relativistic particles to the anisotropic stress at late times.

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.

• Minimally-coupled models beyond simple quintessence.Specifically, we consider “k-essence” models, which are definedby an arbitrary sound speed in addition to a free equation of state parameter w (Armendariz-Picon et al. 2000).

• An example of a generalized scalar field model (Deffayet et al. 2010) and of Lorentz-violating 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.

• Non-universal 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 low-resolution maps of Stokes Q and U polarization measured by LFI at 70 GHz (excluding data from Surveys 2 and 4), foreground-cleaned with the LFI 30 GHz and HFI 353 GHz maps, used as polarized synchrotron and dust templates, respectively (see Planck Collaboration XI 2016).

The foreground-cleaned 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 temperature-polarization likelihood. This likelihood is pixel-based, 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, (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 half-mission cross-spectra, 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 cross-spectra and multipoles are included in the data vector; specifically, the TT100 × 143 and 100 × 217 cross-spectra, 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 co-add the different cross-frequency spectra, since, even after masking the highest dust-contaminated regions, each cross-frequency spectrum has a different, frequency-dependent residual foreground contamination, which we deal with in the model part of the likelihood function.

The model, C(θ) can be rewritten as (8)where C^cmb is the set of CMB C_ℓs, which is independent of frequency, is the foreground model contribution to the cross-frequency 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 Sunyaev-Zeldovich (tSZ) effect; kinetic Sunyaev-Zeldovich (kSZ) effect; tSZ-CIB cross-correlations; 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 cross-spectra 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 non-Gaussianity of the noise, but assumes Gaussian statistics for the dust. The non-whiteness of the noise is estimated from the difference between the cross- and auto-half 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. Beam-shape 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 bin-size 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 large-scale structure introduces dependencies in CMB observables on the late-time 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 late-time modifications of gravity from the CMB alone. The source plane for CMB lensing is the last-scattering surface, so the peak sensitivity is to lenses at z ≈ 2 (i.e., half-way to the last-scattering surface) with typical sizes of order 10² Mpc. Although this peak lensing redshift is rather high for constraining simple late-time 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 non-Gaussianity in the form of a non-zero connected 4-point function, and the conversion of E-mode to B-mode 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 of the CMB lensing potential φ, which are extracted from the Planck temperature and polarization 4-point functions, as presented in Planck Collaboration XV (2016) and discussed further below. Lensing also produces 3-point non-Gaussianity, which peaks in squeezed configurations, due to the correlation between the lensing potential and the ISW effect in the large-angle temperature anisotropies. This effect has been measured at around 3σ with the full-mission 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 bandpowers, covering the multipole range 40 ≤ ℓ ≤ 400. The are estimated from the full-mission temperature and polarization 4-point functions, using the SMICA component-separated maps (Planck Collaboration IX 2016) over approximately 70% of the sky. A large number of tests of internal consistency of the estimated 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 curl-mode tests. Reconstructing the lensing deflections on large angular scales is very challenging because of the large “mean-field” 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 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 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 cross-half-mission 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 cross-frequency spectra, using the multipole range 30 <ℓ < 1000 for the 100 GHz cross-spectra and 500 <ℓ < 2000 for the 217 GHz cross-spectra. 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 cross-correlation 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 B-mode power spectrum through lensing of B-modes (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 B-mode 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 + SN-Ia + H₀, which we now discuss in detail.

4.2.1. Baryon acoustic oscillations

Baryon acoustic oscillations (BAO) are the imprint of oscillations in the baryon-photon plasma on the matter power spectrum and can be used as a standard ruler, calibrated to the CMB-determined 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 (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 six-degree-Field 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₀. 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

Type-Ia 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₀.

Consistently with Planck Collaboration XIII (2016), we use here the analysis by Betoule et al. (2013) of the “Joint Light-curve 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³ 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 low-redshift quantities through the integrated Sachs-Wolfe effect and CMB lensing. The CMB-inferred constraints on the local expansion rate H₀ 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 re-analysis 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: (10)which is within 1σ of the Planck TT+lowP estimate. In this paper we use Eq. (10) as a conservative H₀ 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₀ 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₀ 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 redshift-space 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; Gil-Marín et al. 2012). Moreover, adopting the wrong fiducial cosmological model to convert angles and redshifts into distances can bias measurements of the rate-of-growth 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)σ₈(z), where the growth rate is defined by (11)where σ₈ is calculated including all matter and neutrino density perturbations. Anisotropic clustering also contains geometric information from the Alcock-Paczynski (AP) effect (Alcock & Paczynski 1979), which is sensitive to (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σ₈.

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 redshift-space correlation function of CMASS galaxies on scales > 25 h^-1 Mpc. Both papers provide tight constraints on the quantity fσ₈, and the constraints are consistent. The Samushia et al. (2014) result was shown to behave marginally better in terms of small-scale 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σ₈, 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 large-scale 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, , can be written in terms of the Weyl potential, P_{Φ + Ψ}, by (13)where we have made use of the Limber approximation in flat space, and χ is the comoving distance. The lensing efficiency is given by (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 (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 (16)where the Bessel functions are J₊ = J₀ and J₋ = J₄.

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 N-body 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 N-body simulations, to include the effect of massive neutrinos (Bird et al. 2012);

• 3.

  an update to improve the accuracy on small scales⁴;

• 4.

  an update from higher resolution N-body 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 (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.

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 of the WL likelihood as a function of angular cut. In order for this to remain Δχ²< 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 k-space, 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 Harnois-Déraps et al. (2015) to provide the best fit to small-scale 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 ultra-conservative 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 Ω_bh² = 0.0223 ± 0.0009, n_s = 0.96 ± 0.02, and 40 km s^-1 Mpc^-1<H₀< 100 km s^-1 Mpc^-1 to avoid over-fitting the model, we find that the WL likelihood is insensitive to nonlinear physics. We therefore choose to adopt the tomographic data with the ultra-conservative cuts as our baseline data set.

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₀.

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 non-negligible 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 quintessence-like 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 sub-horizon 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 sub-section we provide the necessary information to compress the Planck CMB power spectrum into a 4-parameter Gaussian likelihood for applications where the full likelihood is too unwieldy.

5.1.1. Taylor expansions of w and w₀,w_a parameterization

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 time-varying equation of state, we expand w(a) in a Taylor series. The first order corresponds to the { w₀,w_a } case, also discussed in Planck Collaboration XIII (2016): (18)We use the parameterized post-Friedmann (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₀, w_a, H₀ and σ₈ are shown in Fig. 3 and the corresponding 2D contours can be found in Fig. 4 for w_a vs. w₀ and for σ₈ 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.

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 χ² between the best-fit in the { w₀,w_a } model and in ΛCDM for Planck TT+lowP+WL is Δχ² = −5.6, which is not very significant for 2 extra parameters (for normal errors a 2σ deviation corresponds to a χ² absolute difference of 6.2). The WL contributes a Δχ² of −2.0 and the (virtually the same as when using Planck TT+lowP alone, for which , which seems to indicate that WL is not in tension with Planck TT+lowP within a (w₀,w_a) cosmology). However, as also discussed in Planck Collaboration XIII (2016), these data combinations prefer very high values of H₀, 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 χ² difference between (w₀,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.

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: (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. 1-parameter varying w

A simple example of a varying w model that can be written in terms of one extra parameter only (instead of w₀,w_a) was proposed in Gott & Slepian (2011), Slepian et al. (2014), motivated in connection to a DE minimally-coupled scalar field, slowly rolling down a potential , 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 , where: (20)at first order in δw₀, which is then the only extra parameter. Marginalized posterior contours in the plane h–δw₀ are shown in Fig. 6. The tightest constraints come from the combination Planck TT+lowP+lensing+BSH that gives δw₀ = −0.008 ± 0.068 at 68% confidence level, which slightly improves constraints found by Aubourg et al. (2015).

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: (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: (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 (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 of uncorrelated parameters that describe w(z). This choice of has been shown to be convenient, since most of the weights (i.e., the rows of ) 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.

5.1.4. Parameterization for a weakly-coupled 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) (23)where the “slope parameter” ϵ_s is defined as: (24)with being a function of the slope of the potential. Here 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: (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 late-Universe slow-roll (ϵ_V ≲ 1 and ) or the moderate-roll (ϵ_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 time-dependence of ϵ_V via corrections to the functional dependence of w(a) at late time. This parameter is defined as the relative difference of at a = a_de and at a → 0, where . 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₀exp(−λφ/M_P) (Wetterich 1988), cosines from pseudo-Nambu Goldstone bosons (pnGB) V = V₀ [1 + cos(λφ/M_P)] (Frieman et al. 1995; Kaloper & Sorbo 2006), power laws V = V₀(φ/M_P)^{− n} (Ratra & Peebles 1988), and models motivated by supergravity (SUGRA) V = V₀(φ/M_P)^{− α}exp [(φ/M_P)²] (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 1-parameter case (i.e., ϵ_s only, with ϵ_∞ = 0 and ζ_s = 0, describing “thawing” quintessence/phantom models, where in the early Universe) and the 3-parameter case (general quintessence/phantom models where an early-Universe fast-rolling phase is allowed). When we vary the data sets and theoretical prior (between the 1-parameter and 3-parameter 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 3-parameter case, because in this case flatter potentials are preferred in the late Universe in order to slow-down 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.

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 small-scale measurements and CMB lensing. Assuming a constant fraction of Ω_e until recent times (Doran & Robbers 2006), the DE density is parameterized as: (26)This expression requires two parameters in addition to those of ΛCDM, namely Ω_e and w₀, while 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.

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: (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 non-negligible 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 non-negligible. 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.

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⁵: 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 (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 non-zero 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₀,w_a }, for which the posterior is very non-Gaussian 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 large-scale 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₀,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 non-flat 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 non-negligible.

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 “top-down” approach that starts from the most general action allowed by symmetry and selects from there interesting classes belonging to so-called “Horndeski models”, which, as mentioned in Sect. 3.2.1, include almost all stable scalar-tensor theories, universally coupled, with second-order 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 “bottom-up” 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 time-dependent functions , 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)⁶. 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 ghost-free, 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, , , , , , , are internally redefined in terms of the dimensionless parameters α_i with i running from 1 to 6: We will always demand that which eliminates models containing higher-order 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 (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: 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 B-mode polarization; α_M affects the lensing potential, but also the amplitude of the primordial polarization peak in B-modes (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, where ℋ is the conformal Hubble function, m₀ 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 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⁷ 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 from Eq. (31). For simplicity we also turn off all other higher-order EFT operators and set . 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₂ = 0 as a function of c and α_T = α_H = 0. We are thus considering non-minimally coupled “K-essence” 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: (43)By choosing a non-zero α_M (and therefore a time evolving Ω) we introduce a non-minimal coupling in the action (see Eq. (2)), which will lead to non-zero anisotropic stress and to modifications of the lensing potential, typical signatures of MG models. Here we will use a scaling ansatz, α_M = α_M0a^β, 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 (44)which coincides with the built-in exponential model of EFTCAMB for Ω₀ = α_M0/β. The marginalized posterior distributions for the two parameters Ω₀ and β are plotted in Fig. 12 for different combinations of data. For α_M0 = 0 we recover ΛCDM. For small values of Ω₀ and for β = 1, the exponential reduces to the built-in linear evolution in EFTCAMB, (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,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.

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 DE-related case that we will define as case 1 below, but using μ(a) and Σ(a) instead of μ(a,k) and η(a,k)⁸. They found the constraints μ₀−1 = 0.05 ± 0.25 and Σ₀−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 μ₀−1 and Σ₀−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 μ₀ and η₀ and derived from those the limits on Σ₀, using WMAP-5 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 for z_s = 1 and μ₀ = 0.87 ± 0.12,η₀ = 1.3 ± 0.35 for z_s = 2. A similar parametrization was also used in Daniel et al. (2010) in terms of μ₀ and ϖ (equivalent to μ₀−1 and η₀−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 <μ₀< 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⁹ (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: For large length scales (small k), the two functions reduce to μ → 1 + f₁(a)c₁ and η → 1 + f₂(a)c₂; for small length scales (large k), one has μ → 1 + f₁(a) and η → 1 + f₂(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.

  time-related 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 time-related 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).

For case 1 (referred to as “DE-related” parameterization) we then have five free parameters, E₁₁, c₁, E₂₂, c₂, and λ, while for case 2 (the “time-related” parameterization) we have two additional parameters, E₁₂ and E₂₁. The choice above looks very similar to the BZ parameterization (Bertschinger & Zukin 2008) for the quasi-static limit of f(R) and scalar-tensor theories. However, we emphasise that Eqs. (46) and (47) should not be seen as a quasi-static 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 “DE-related” 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 scale-independent case in which we fix c₁ = c₂ = 1. Results show that a scale dependence of μ and η does not lead to a significantly smaller χ² with respect to the scale-independent case, both for the DE-related and time-related 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 scale-independent parameterization.

Table 6 shows results for the DE-related 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, μ₀−1 and η₀−1 are obtained by reconstructing Eqs. (46) and (47) from E₁₁ and E₂₂ 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 (μ₀−1, η₀−1) and (μ₀−1, Σ₀−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 DE-related and time-related 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 DE-related 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₁₁ and E₂₂) leads to an improvement of Δχ² = −6.3 when using Planck TT+lowP (similarly divided between lowP and TT) and of Δχ² = −6.4 when including BSH (with a equally divided between TT and lowP). When Planck data (TT+lowP) are combined also with WL data, the tension increases to Δχ² = −10.6 (with the CMB still contributing about the same amount, ). When considering Planck TT+lowP+BAO/RSD, Δχ² = −8.1 with respect to ΛCDM while, when combining both WL and BAO/RSD, the tension is maximal, with Δχ² = −10.8 and . There is instead less tension for the time-related parameterization, as is visible in the right panel of Fig. 14.

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 DE-related dependence forces the combination to be compatible with ΛCDM in the past, when the DE density is negligible; the time-related 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₁₁ = −0.3, E₂₂ = 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 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 4-point 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 DE-related parameterization is more in tension with ΛCDM than the time-related 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.

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 scale-independent 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.

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 scalar-tensor 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 scalar-tensor 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 k-essence

In minimally coupled quintessence models, the sound speed is and DE does not contribute significantly to clustering. However, in so-called “k-essence” models, the kinetic term in the action is generalised to an arbitrary function of (∇φ)² (Armendariz-Picon 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 non-negligible on sub-horizon 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 (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 k-essence type models the impact of the DE perturbations on the total clustering is small when 1 + w ≈ 0. For the DE perturbations in k-essence 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 Γ: (48)When Γ = 0 the perturbations are adiabatic, that is .

For this purpose, it is convenient to adopt the “equation of state” approach described in Battye et al. (2015), Soergel et al. (2015). The gauge-invariant 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 Lorentz-violating massive gravity (LVMG) for which (σ ≠ 0,Γ = 0) and generalized scalar field models (GSF) in which the anisotropic stress is zero (σ = 0,Γ ≠ 0).

Lorentz-violating 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, ξⁱ, that has a physical interpretation as an elastic medium, or as Lorentz-violating 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 non-vanishing anisotropic stress: (49)including one degree of freedom, the sound speed , 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 . 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 into account by using the combination 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 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 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 second-order field equations. In this case the anisotropic stresses are zero and (50)This has three extra parameters (α,β₁,β₂), in addition to w. If β₁ = 1 and β₂ = 0 this becomes the generalized k-essence model. An example of this class of models is “kinetic gravity braiding” (Deffayet et al. 2010) and similar to the non-minimally coupled k-essence 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, β₁ and β₂. 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 α₁ = 0.2 and α₂ = 1 were chosen to decorrelate w and γ_i. The resulting 2σ upper limits are γ₁< 0.67 and γ₂< 0.61 (for w> −1), γ₂< 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 β₁ and β₂.

5.3.3. Universal couplings: f(R) cosmologies

A well-investigated class of MG models is constituted by the f(R) theories that modify the Einstein-Hilbert action by substituting the Ricci scalar with a more general function of itself: (51)where κ² = 8πG. f(R) cosmologies can be mapped to a subclass of scalar-tensor theories, where the coupling of the scalar field to the matter fields is universal.

For a fixed background, the Friedmann equation provides a second-order 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 (52)and the other initial (or boundary condition), usually called B₀, is the present value of (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₀ 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₀, 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 quasi-static 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₀. Adding any structure formation probe, such as WL, RSD or CMB lensing, breaks the degeneracy. Figure 20 shows the likelihood of the B₀ parameter using EFTcamb, where a ΛCDM background evolution is assumed, i.e., w_DE = −1.

As the different data sets provide constraints on B₀ that vary by more than four orders of magnitude, we show plots for log ₁₀B₀; to make these figures we use a uniform prior in log ₁₀B₀ to avoid distorting the posterior due to prior effects. However, for the limits quoted in the tables we use B₀ (without log ) as the fundamental quantity and quote 95% limits based on B₀. In this way the upper limit on B₀ 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 ₁₀B₀. Overall this appears to be the best compromise to present the constraints on the B₀ parameter. In the plots, the GR value (B₀ = 0) is reached by a plateau stretching towards minus infinity.

Finally, we note that f(R) models can be studied also with the MGcamb parametrization, assuming the quasi-static limit. We find that for the allowed range of the B₀ parameter, the results with and without the quasi-static 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₀< 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. Non-universal 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 non-universal 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 (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 scalar-tensor 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: (55)Here each component is treated as a fluid with stress energy tensor , where u_μ = (−a,0,0,0) is the fluid 4-velocity 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₀exp^{− βφ} (originally motivated by Weyl scaling scalar-tensor theories), where m₀ 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 Klein-Gordon equation, which now includes an extra term that depends on the density of cold dark matter: (56)Following Pettorino & Baccigalupi (2008), we choose an inverse power-law potential defined as V = V₀φ^{− α}, with α and V₀ being constants. The amplitude V₀ is fixed thanks to an iterative routine, as implemented by Amendola et al. (2012), Pettorino et al. (2012). To a first approximation α only affects late-time 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 β², i.e. (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 last-scattering 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 non-zero 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.

When comparing with ΛCDM, however, the goodness of fit does not improve, despite the additional parameters. Only the 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 non-zero 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 Ω_φ ∝ β² for an inverse power-law potential) and therefore lead to a smaller Ω_m today; as a consequence, the distance to the last-scattering 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₀ and higher σ₈.

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 χ² does not improve between coupled DE and ΛCDM.

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 minimally-coupled 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 story-line 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₀,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₀,w_a) model with respect to ΛCDM at about 2σ, although with a preference for high values of H₀ (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.

“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₀,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 (top-down) that starts from a very general EFT action for DE (Sect. 5.2.1); the other (bottom-up) 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 non-minimally coupled K-essence model (i.e.α_B, α_M, α_K are varying functions of the scale factor). We have found limits on the present value α_M0< 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 (DE-related and time-related) 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 k-essence 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 Lorentz-breaking 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 β². 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), non-local 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 σ₈ or H₀ 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, ISW-lensing and B-mode polarization, all of which can be used to further test MG scenarios.

---

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/INSU-IN2P3-INP (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/planck-collaboration.

References

All Tables

All Figures