MGCAMB with massive neutrinos and dynamical dark energy

TL;DR

This paper introduces an upgraded version of MGCAMB, a tool for testing gravity theories with cosmological data, now compatible with latest datasets and including massive neutrinos and dynamical dark energy.

Contribution

The paper presents a major upgrade to MGCAMB, integrating massive neutrinos and dynamical dark energy, and improving code flexibility for cosmological tests.

Findings

Consistent constraints with the ΛCDM model from recent observations.

Enhanced code compatibility with latest CosmoMC and datasets.

Successful implementation of massive neutrinos and dynamical dark energy.

Abstract

We present a major upgrade of MGCAMB, a patch for the Einstein-Boltzmann solver CAMB used for phenomenological tests of general relativity against cosmological datasets. This new version is compatible with the latest CosmoMC code and includes a consistent implementation of massive neutrinos and dynamical dark energy. The code has been restructured to make it easier to modify to fit the custom needs of specific users. To illustrate the capabilities of the code, we present joint constraints on the modified growth, massive neutrinos and the dark energy equation of state from the latest cosmological observations, including the recent galaxy counts and weak lensing measurements from the Dark Energy Survey, and find a good consistency with the {\Lambda}CDM model.

Tables

Table 1. Table 1: The 68 % CL uncertainties and best fit values of parameters obtained using the original CosmoMC, compared to the results from the GR limit of MGCosmoMC for two different values of a trans subscript 𝑎 trans a_{\rm trans} which sets the scale factor beyond which the modified set of equations is evolved. The bound on the net mass of neutrinos is at the 95% CL.

Parameter	CosmoMC	MGCosmoMC ($a_{\mathrm{trans}}=0.01$)	MGCosmoMC ($a_{\mathrm{trans}}=0.001$)
${\omega }_b$	$0.02237\pm 0.00014$	$0.02237\pm 0.00014$	$0.02237\pm 0.00014$
${\omega }_c$	$0.1178\pm 0.0011$	$0.1178\pm 0.0011$	$0.1178\pm 0.0011$
$100{\theta }_{\mathrm{MC}}$	$1.04095\pm 0.00030$	$1.0410\pm 0.0003$	$0.104095\pm 0.00030$
$\tau$	$0.076\pm 0.015$	$0.075\pm 0.015$	$0.075\pm 0.015$
$\sum m_{\nu }$ (95 % CL)	eV	eV	eV
$n_s$	$0.9684\pm 0.0042$	$0.9684\pm 0.0042$	$0.9684\pm 0.0042$
$\ln \left({10}^{10}{A}_s\right)$	$3.080\pm 0.028$	$3.079\pm 0.027$	$3.080\pm 0.029$
Best fit: $-\log \left(\mathrm{Like}\right)$	$7023.371$	$7023.607$	$7023.964$

Table 2. Table 2: The 68 % CL uncertainties and best fit values of parameters constrained using MGCosmoMC. The bound on the net mass of neutrinos is at the 95% CL. Model 0 corresponds to Λ Λ \Lambda CDM with massive neutrinos. Model 1, in addition, includes modified growth on the Λ Λ \Lambda CDM background, while Model 2 adds a varying DE equation of state using the CPL parameterization.

Parameter	Model 0	Model 1	Model 2
${\omega }_b$	$0.02237\pm 0.00014$	$0.02239\pm 0.00014$	$0.02231\pm 0.00016$
${\omega }_c$	$0.1178\pm 0.0011$	$0.1175\pm 0.0011$	$0.1183\pm 0.0013$
$100{\theta }_{\mathrm{MC}}$	$1.0409\pm 0.0003$	$1.0410\pm 0.0003$	$1.0408\pm 0.0003$
$\tau$	$0.075\pm 0.015$	$0.067\pm 0.017$	$0.072\pm 0.018$
$n_s$	$0.969\pm 0.004$	$0.969\pm 0.004$	$0.967\pm 0.005$
$\ln \left({10}^{10}{A}_s\right)$	$3.08\pm 0.03$	$3.06\pm 0.03$	$3.07\pm 0.03$
$\sum m_{\nu }$ (95 % CL)	eV	eV	eV
${\mu }_0-1$	$0$	$-0.09\pm 0.30$	$-0.07\pm 0.29$
${\gamma }_0-1$	$0$	$0.46\pm 0.79$	$0.43\pm 0.77$
${\Sigma }_0-1$	$0$	$0.01\pm 0.06$	$0.02\pm 0.07$
$w_0$	$-1$	$-1$	$-0.84\pm 0.13$
$w_a$	$0$	$0$	$-0.48\pm 0.36$
${\chi }^2$	$7023.37$	$7023.04$	$7024.44$
$\Delta {\chi }^2$	-	$-0.33$	$+1.07$

Equations

$$d{s}^2=a(\tau {)}^2\left[-(1+2\Psi )d{\tau }^2+(1-2\Phi )d{x}^2\right],$$

$$T_0^0+\delta T_0^0$$

$$T_i^0+\delta T_i^0$$

$$T_j^i+\delta T_j^i$$

$$\dot{\delta }$$

$$\dot{\theta }$$

$$k^2\Phi =-4\pi G{a}^2\rho \Delta ,$$

$$\rho \Delta \equiv \rho \delta +\frac{3\mathcal{H}}{k^2}(\rho +P)\theta .$$

$$k^2(\Phi -\Psi )=12\pi G{a}^2(\rho +P)\sigma ,$$

$$k^2\Psi =-4\pi G{a}^2\left[\rho \Delta +3(\rho +P)\sigma \right].$$

$$k^2\Psi =-4\pi G\mu (a,k)a^2\left[\rho \Delta +3(\rho +P)\sigma \right],$$

$$k^2[\Phi -\gamma (a,k)\Psi ]=12\pi G\mu (a,k)a^2(\rho +P)\sigma .$$

$$k^2(\Phi +\Psi )=-4\pi G\Sigma (a,k)a^2\left[2\rho \Delta +3(\rho +P)\sigma \right].$$

$$d{s}^2=a(\tau {)}^2\left[-d{\tau }^2+({\delta }_{ij}+h_{ij}(\mathbf{x},\tau ))d{x}^{i}d{x}^j\right],$$

$$h_{ij}=\int d^{3}k{e}^{i\mathbf{k}\cdot \mathbf{x}}[{\hat{\mathbf{k}}}_i{\hat{\mathbf{k}}}_{j}h(\mathbf{k},\tau )+({\hat{\mathbf{k}}}_i{\hat{\mathbf{k}}}_j-{\delta }_{ij}/3)6\eta (\mathbf{k},\tau )].$$

$$\Phi$$

$$\Psi$$

$$\dot{\delta }$$

$$\dot{\theta }$$

$$k^2\eta +k\mathcal{H}\mathcal{Z}=-4\pi G{a}^2\rho \delta .$$

$$\mathcal{Z}=k\alpha -\frac{3\dot{\eta }}{k}.$$

$$\alpha =\left\{\eta +\frac{\mu a^2}{2k^2}\left[\gamma \rho \Delta +3(\gamma -1)(\rho +P)\sigma \right]\right\}\frac{1}{\mathcal{H}},$$

$$\begin{array}{rl}\eta & =\mathcal{H}\alpha -\frac{\mu a^2}{2k^2}\left\{\gamma \rho \Delta +3(\gamma -1)\rho (1+w)\sigma \right\}\\ & =\mathcal{H}\alpha -\frac{\mu a^2}{2k^2}\Gamma ,\end{array}$$

$$\Gamma =\gamma \rho \Delta +3(\gamma -1)\rho (1+w)\sigma .$$

$$\begin{array}{rl}(\rho \Delta {)}^{\cdot }& =-3\mathcal{H}\rho \Delta -(1+w)\rho \theta \left[1+\frac{3}{k^2}({\mathcal{H}}^2-\dot{\mathcal{H}})\right]\\ & -3\mathcal{H}\rho (1+w)\sigma -(1+w)\rho k\mathcal{Z}.\end{array}$$

$$\begin{split}\dot{\eta}=&\frac{1}{2}\frac{a^{2}}{3\rho a^{2}\mu\gamma(1+w)/2+k^{2}}\\ &\times\biggl{\{}\rho(1+w)\mu\gamma\theta\left[1+3\frac{\mathcal{H}^{2}-\dot{\mathcal{H}}^{2}}{k^{2}}\right]\\ &+\rho\Delta\left[\mathcal{H}\mu(\gamma-1)-\dot{\mu}\gamma-\mu\dot{\gamma}\right]\\ &+3\mu(1-\gamma)\rho(1+w)\dot{\sigma}\\ &+k^{2}\alpha\left[\rho\mu\gamma(1+w)-2\left(\frac{\mathcal{H}^{2}-\dot{\mathcal{H}}}{a^{2}}\right)\right]\\ &+3\mathcal{H}\mu(\gamma-1)(1+w)\rho\sigma(3w+2)\\ &-3(1+w)\rho\sigma\left[\dot{\mu}(\gamma-1)-\dot{\gamma}\mu+\mu(1-\gamma)\frac{\dot{w}}{1+w}\right]\biggr{\}},\end{split}$$

$$\dot{\alpha }=\Phi +\Psi -\eta ,$$

$$(1+w)\theta =kq,\frac{3}{2}(1+w)\sigma =\Pi ,$$

$$\begin{split}\dot{\eta}&=\frac{1}{2}\frac{a^{2}}{3\rho a^{2}\mu\gamma(1+w)/2+k^{2}}\biggl{\{}\rho\mu\gamma kq\left[1+3\frac{\mathcal{H}^{2}-\dot{\mathcal{H}}}{k^{2}}\right]\\ &+\rho\Delta\left[\mathcal{H}\mu(\gamma-1)-\dot{\mu}\gamma-\mu\dot{\gamma}\right]\\ &+2\mu(1-\gamma)\rho\dot{\Pi}\\ &+k^{2}\alpha\left[\rho\mu\gamma(1+w)-2\left(\frac{\mathcal{H}^{2}-\dot{\mathcal{H}}}{a^{2}}\right)\right]\\ &+2\rho\Pi\left[\mathcal{H}(\gamma-1)(3w+2)\mu-\dot{\mu}(\gamma-1)-\dot{\gamma}\mu\right]\biggr{\}}.\end{split}$$

$$C_{\ell }^{\mathrm{TT}}=(4\pi {)}^2\int dk{k}^2|{\Delta }_{\ell }^{\mathrm{T}}(k){|}^2{P}_{\mathcal{R}}(k),$$

$${\Delta }_{\ell }^{\mathrm{T}}(k)={\int }_0^{{\tau }_0}d\tau S_{\mathrm{T}}(k,\tau )j_{\ell }(k\tau ),$$

$$\begin{array}{rl}{S}_{\mathrm{T}}(k,\tau )=& g\left({\Delta }_{\mathrm{T}0}+2\dot{\alpha }+\frac{{\dot{v}}_b}{k}+\frac{{\Pi }^{\mathrm{pol}}}{4}+\frac{3{\ddot{\Pi }}^{\mathrm{pol}}}{4k^2}\right)\\ & +e^{-\kappa }(\dot{\eta }+\ddot{\alpha })+\dot{g}\left(\alpha +\frac{v_b}{k}+\frac{3{\dot{\Pi }}^{\mathrm{pol}}}{2k^2}\right)\\ & +\frac{3\ddot{g}{\Pi }^{\mathrm{pol}}}{4k^2},\end{array}$$

$$\ddot{\alpha }+\dot{\eta }=\dot{\Phi }+\dot{\Psi },$$

Full text

1. I Introduction

2. II Modified Growth framework

3. III The new MGCAMB patch

4. III.1 Massive Neutrinos

5. III.2 The CMB source function and the weak lensing transfer function

6. III.3 The GR limit of MGCAMB

7. IV Joint constraints on modified growth, massive neutrinos and the dark energy equation of state

8. IV.1 The datasets

9. IV.2 The GR limit consistency check

10. IV.3 Results

11. V Discussion

12. A Other parameterizations

13. A.1 The $\mu,\Sigma$ parameterization

14. A.2 The $Q,R$ parameterization

15. B Linear DES data

16. C Galaxy Clustering - Weak lensing theory with Weyl potential.

I Introduction

Since the discovery of cosmic acceleration two decades ago (Perlmutter et al., 1999; Riess et al., 1998), explaining it has been one of the primary goals of cosmology. The broadly accepted working model, in which the accelerated expansion is driven by the cosmological constant $\Lambda$ and most of the remaining density is in cold dark matter (CDM), fits remarkably well a plethora of observations, such as the cosmic microwave background (CMB) anisotropies (Bennett et al., 2013; Ade et al., 2016a), baryon acoustic oscillations (BAO) (Percival et al., 2002; Beutler et al., 2011a; Alam et al., 2017), type Ia supernovae (Conley et al., 2011; Suzuki et al., 2012), galaxy clustering (Zehavi et al., 2011) and galaxy lensing (Heymans et al., 2012; Hildebrandt et al., 2017). However, the nature of CDM remains unknown, and the observed value of $\Lambda$ requires technically unnatural fine-tuning to reconcile it with the large vacuum energy density expected from quantum field theory (Burgess, 2017). These questions, along with new opportunities for testing gravity on cosmological scales afforded by the current and upcoming surveys Abbott et al. (2005); Abell et al. (2009); Laureijs et al. (2011), have motivated extensive studies of extensions of General Relativity (GR) Silvestri and Trodden (2009); Clifton et al. (2012); Ishak (2018).

Gravitational potentials evolve differently in alternative gravity theories compared to $\Lambda$CDM, leading to different predictions for the growth of cosmic structures. At linear order in cosmological perturbations, one can search for modified growth patterns phenomenologically, by introducing two functions that parameterize the altered relations between the Newtonian potential and the curvature perturbation, and between the matter density contrast and the Newtonian potential Amendola et al. (2008); Zhang et al. (2007); Hu and Sawicki (2007); Bertschinger and Zukin (2008); Bean and Tangmatitham (2010); Pogosian et al. (2010). Modified Growth with CAMB (MGCAMB) (Zhao et al., 2009; Hojjati et al., 2011) is a patch for the popular Einstein-Boltzmann solver CAMB Lewis et al. (2000), enabling the calculation of cosmological observables for a given form of such phenomenological functions in manner suitable for constraining them with CosmoMC Lewis and Bridle (2002); Lewis (2013) or other similar Monte-Carlo Markov Chain algorithms.

The original version of MGCAMB (Zhao et al., 2009) and its update in Hojjati et al. (2011) were largely based on the assumption that modifications of gravity appear well after the radiation-matter equality, and that the role played by the anisotropic stress in relativistic particle species is negligible. This limited the accuracy of modelling the effects of massive neutrinos. On the other hand, the neutrino mass can no longer be neglected in cosmological predictions, as the upcoming surveys are expected to probe masses close to the measured difference of $0.05$ eV between the masses of different neutrino flavours. Massive neutrinos contribute to the expansion rate as matter, but stream out of smaller gravitational potentials, suppressing the growth on small scales. This effect can be partially degenerate with those of modifications of gravity and, therefore, must be accounted for.

Another limitation of the prior versions of MGCAMB was the assumption of constant dark energy density. Generally, modified gravity (MG) theories predict modifications of the expansion history along with the modified growth of structures. The ability to study the covariance of the two can be important for ruling out broad classes of alternative gravity theories Espejo et al. (2018).

This paper presents a major update of MGCAMB111This new patch is available at https://github.com/sfu-cosmo/MGCAMB allowing for dynamical dark energy and accurate modelling of massive neutrinos. This version is compatible with the latest CosmoMC and is restructured to make it easier to customize to work with different parameterizations of phenomenological functions. The modifications of gravity can now be introduced at arbitrarily high redshifts, as long as the phenomenological functions continuously approach their $\Lambda$CDM values in the past. Users wishing to study models that do not approach the GR limit in the past will need to introduce the corresponding changes in the initial conditions.

To demonstrate the capabilities of this new version, we present joint constraints on the modified growth, massive neutrinos and the dark energy equation of state from the latest cosmological observations, including the Planck 2015 CMB data Aghanim et al. (2016), the Joint Light Curve analysis (JLA) supernovae (Betoule et al., 2014), the BAO measurements from 6dF (Beutler et al., 2011b) and SDSS DR7 (Ross et al., 2015), the measurements of the Hubble parameter, the angular-diameter distance and the redshift space distortions from BOSS DR12 (Alam et al., 2017), and the recent galaxy clustering and weak lensing measurements from the Dark Energy Survey (DES) (Abbott et al., 2018a).

The paper is organized as follows. In Section II, we briefly review the framework for phenomenological tests of gravity on cosmological scales. In Section III, we detail the implementation of the MGCAMB patch and the accuracy tests. Then, in Section IV, we demonstrate the use of MGCAMB by deriving joint constraints on modified gravity, massive neutrinos and the dark energy equation of state from the latest cosmological datasets. We conclude with a discussion in Section V.

II Modified Growth framework

When describing deviations from GR on large scales, we assume that the universe is well described by a FRW metric with small perturbations. We adopt the conformal Newtonian gauge and write the line element as

[TABLE]

where $\Psi$ and $\Phi$ are the scalar gravitational potential and the curvature perturbation, respectively, and both depend on the conformal time $\tau$ and the comoving coordinate ${\mathbf{x}}$. We consider perturbations of the total energy-momentum tensor $T^{\mu}_{\nu}$ denoted as

[TABLE]

where $\delta$ is the density contrast, $v$ is the velocity field, $\delta P$ the pressure perturbation and $\pi^{i}_{j}$ denotes the traceless anisotropic stress tensor. The energy momentum tensor components evolve according to the conservation equation, $T^{\mu}_{\nu;\mu}=0$; working in Fourier space, the perturbations obey Bertschinger (2006)

[TABLE]

where $\theta$ is the divergence of the velocity field. These equations hold for the combination of all cosmological fluids or for any decoupled subset of fluids such as CDM, the photon-baryon fluid and massive neutrinos after decoupling.

To close the system of equations for cosmological perturbations, one needs two additional equations relating the perturbed energy momentum tensor and the metric potentials $\Phi$ and $\Psi$. In $\Lambda$CDM, one can combine the ${}^{0}_{0}$ and the divergence of the ${}^{0}_{i}$ Einstein equations to obtain

[TABLE]

where $\Delta$ is the gauge-invariant density contrast,

[TABLE]

One can also use the traceless part of the ${}^{i}_{j}$ Einstein equation, given by

[TABLE]

and use this to write an equation relating $\Psi$ and $\Delta$:

[TABLE]

Eqs. (9) and (10), combined with the conservation equations (5) and (6), can be used in Einstein-Boltzmann solvers, such as CAMB and CLASS, to compute the cosmological observables.

Einstein’s equation would be modified in alternative gravity theories. For example, in scalar-tensor theories of gravity, the scalar field alters the relation between $\Phi$ and $\Psi$, and mediates an additional force between massive particles, which enhances the growth of structure. This can be phenomenologically modelled by introducing two functions of time and scale, $\mu(a,k)$ and $\gamma(a,k)$, that encode possible modification of (9) and (10), defined via

[TABLE]

Given the form of $\mu(a,k)$ and $\gamma(a,k)$, Eqs (11) and (12), along with (5) and (6), can be used to evolve the system of equations and compute the cosmological observables of interest. Note that $\mu(a,k)$ is introduced as a modification of Eq. (10) that relates $\Psi$ and $\Delta$, instead of (7), because matter perturbations respond to the gradients of the Newtonian potential $\Psi$. This makes $\mu$ a parameter that directly controls the strength of the gravitational interaction.

The gravitational slip function $\gamma$ is not directly related to cosmological observables (Amendola et al., 2008, 2013) and, therefore, is generally difficult to constrain without fixing $\mu$ or making additional assumptions. Instead, it is often more informative to work with $\mu$ paired with function $\Sigma$ that modifies the relation between the lensing potential $(\Phi+\Psi)$ and $\Delta$ via

[TABLE]

In the limit of negligible anisotropic stress, $\Sigma$ is simply related to $\mu$ and $\gamma$ through $\Sigma=\mu(1+\gamma)/2$. One can break the degeneracy between $\mu$ and $\Sigma$ and constrain both of them independently by combining data from clustering surveys with measurements of weak lensing Song et al. (2011); Simpson et al. (2013). While it might be obvious to many readers, we would still like to note that constraints on $\mu$ depend on whether one marginalizes over $\gamma$ or $\Sigma$.

The above-mentioned framework is implemented in MGCAMB, which can be used for two purposes. One can adopt functional forms of $\mu$ and $\gamma$, or $\mu$ and $\Sigma$, and fit the function parameters to data to search for any departure from the $\Lambda$CDM values of $\mu=\gamma=\Sigma=1$. One can go even further and reconstruct the functions from the data using the correlated prior approach Crittenden et al. (2009, 2012); Zhao et al. (2017); Wang et al. (2018), or the Gaussian Process Reduction Holsclaw et al. (2010a, b). Another way to use MGCAMB is to study predictions of specific theories for certain cosmological observables. This application is limited by the fact that deriving the analytical forms for $\mu$, $\gamma$ and $\Sigma$ in a given theory requires adopting the quasi-static approximation (QSA) Silvestri et al. (2013). The validity of the QSA depends on the observable and the strength of modifications introduced by the theory on near-horizon scales. As a rule, QSA tends to work in most viable models Silvestri et al. (2013); Sawicki and Bellini (2015) and thus MGCAMB can be a good starting point for looking at the characteristic observational signatures of a theory.

III The new MGCAMB patch

The set of equations used in the new MGCAMB patch is based on and has a large overlap with the previous versions (Zhao et al., 2009; Hojjati et al., 2011). However, there are several important differences and, for the sake of completeness, we present the entire formalism.

In CAMB, cosmological perturbations are evolved in synchronous gauge, with the line element given by (Bertschinger, 2006)

[TABLE]

where

[TABLE]

The Newtonian gauge potentials $\Phi$ and $\Psi$ are related to the synchronous gauge potentials $\eta$ and $h$ through

[TABLE]

where $\alpha=(\dot{h}+6\dot{\eta})/2k^{2}$. The synchronous gauge perturbations of the energy-momentum tensor evolve according to

[TABLE]

In order to evolve the perturbations, one needs to compute the quantity $\dot{h}$, or the quantity $\mathcal{Z}$ defined in CAMB as $\mathcal{Z}\equiv\dot{h}/2k$. In CAMB, this is done using the ${}^{0}_{0}$ Einstein equation in synchronous gauge

[TABLE]

In MGCAMB, the Einstein equations are modified and, hence, one needs an alternative way to compute $\mathcal{Z}$. From the definition of $\alpha$, we have

[TABLE]

To find $\alpha$ and $\dot{\eta}$, we start by substituting Eqs. (16) and (17) into the modified Einstein equations (11) and (12), and combine the resulting equations to write the following expression for $\alpha$:

[TABLE]

where we now include the factor $8\pi G$ in the definition of density and pressure, e.g. $8\pi G\rho\to\rho$. To derive an equation for $\dot{\eta}$, we first rearrange the equation above to solve for $\eta$, obtaining

[TABLE]

where we defined

[TABLE]

Next, we would like to differentiate Eq. (23) with respect to $\tau$. To compute $(\rho\Delta)^{\cdot}$, we combine the conservation equations (18) and (19) to obtain

[TABLE]

Finally, taking the derivative of Eq. (23) and using Eqs. (21) and (25), leaves us with the equation for $\dot{\eta}$:

[TABLE]

where we replaced $\dot{\alpha}$ with

[TABLE]

and used the modified Einstein equations to express $\Phi$ and $\Psi$ in terms of the energy-momentum perturbations.

Following the notation in CAMB, we introduce fluxes $q$ and the anisotropic stress $\Pi$ related to the velocity divergence $\theta$ and the anisotropic stress $\sigma$ through

[TABLE]

to rewrite the equation for $\dot{\eta}$ as

[TABLE]

The background and perturbation variables of energy-momentum appearing in the above equations are sums over the uncoupled fluid components, e.g. $\rho=\rho_{\rm b+\gamma}+\rho_{\rm CDM}+\rho_{\nu}$, etc.

There are two notable differences between Eq. (29) and its counterpart in the previous MGCAMB patch. First, the factor $(3w+2)$ in the last line corrects a typo present in the previous version222This factor is $(3w)$ in the last term of Eq. (36) of (Hojjati et al., 2011).. As this correction is proportional to $\Pi$, it has a negligible effect at late times when the anisotropic stress is small. More importantly, the pre-factors of $\alpha$ and $q$ are now generalized to allow for an arbitrary expansion history. The expression for $\dot{\eta}$ in the previous version of MGCAMB assumed dark energy with the equation of state $w_{\rm DE}=-1$.

As in previous versions, the present MGCAMB patch assumes that GR is recovered deep in the radiation era. The code starts with the same initial conditions as CAMB and evolves the original CAMB system of equations up to a certain value of the scale factor, $a_{\rm trans}$, set by the parameter GRtrans. After that, the code evolves the alternative equations described above. Unlike the previous version, the present patch has no restriction on how early the switch can happen, as long as it happens after the time at which the initial conditions are set for the smallest values of $k$ and the phenomenological functions are such that the GR limit is approached continuously.

Computing $\mathcal{Z}$ requires knowing the quantities $\delta$, $q$, $\Pi$ and $\dot{\Pi}$, which can be a challenging problem depending on the epoch at which $\dot{\eta}$ is evaluated. For example, at late times, CAMB stops evolving the full set of Boltzmann equations for photons and neutrinos and uses the radiation streaming approximation (RSA) instead (Blas et al., 2011). In the RSA, $\delta_{\gamma}$ and $\delta_{\nu}$ are computed by using approximated versions of $\mathcal{Z}$ and $\dot{\mathcal{Z}}$ that do not include radiation. However, the current CAMB implementation of RSA uses the ${}^{i}_{i}$ Einstein equation. Since we do not have all the modified Einstein equations, we opt to use the RSA implementation from an older version of CAMB, which did not depend on the ${}^{i}_{i}$ Einstein equation. For small values of GRtrans, in order to preserve the accuracy we had to increase the time at which the RSA is switched on, which slows down MGCAMB with respect to the default CAMB by a factor of two.

Before last scattering, CAMB uses a second order tight coupling expansion in which the computation of $q_{\gamma}$ and $q_{b}\equiv v_{b}$ requires the knowledge of $\mathcal{Z}$ and $\sigma^{*}\equiv k\alpha$. We resolve this by using the values of these quantities computed at the previous time-step.

In addition to the ($\mu,\gamma$) parameterization, MGCAMB offers options to work with ($\mu,\Sigma$) introduced in Sec. II and the ($Q,R$) functions of (Bean and Tangmatitham, 2010). More details on their implementation are given in Appendix A.

III.1 Massive Neutrinos

The default CAMB code calculates the quantities $\dot{\Pi}$ after computing $\mathcal{Z}$. In our case, $\dot{\Pi}$ is required to compute $\mathcal{Z}$. The computation of $\dot{\Pi}_{\gamma}$ and $\dot{\Pi}_{r}$ in the previous versions of MGCAMB ignored the contribution from massive neutrinos. In the current version, we compute $\dot{\Pi}_{\nu}$ by integrating the neutrinos equations before the computation of $\mathcal{Z}$. This is safe, since the equations for $\dot{\Pi}_{\nu}$ do not depend on $\mathcal{Z}$.

III.2 The CMB source function and the weak lensing transfer function

The CMB temperature angular power spectrum is given by (Zaldarriaga and Seljak, 1997)

[TABLE]

where $P_{\cal R}(k)$ is the primordial curvature perturbation power spectrum and the CMB temperature transfer function is given by

[TABLE]

where $S_{\rm T}(k,\tau)$ is the source term and $j_{\ell}(x)$ are the spherical Bessel functions. In synchronous gauge, the source is

[TABLE]

where $\kappa$ is the optical depth, $g$ is the visibility function and $\Pi^{\rm pol}$ is the polarization term. In CAMB, the calculation of the ISW term, $(\dot{\eta}+\ddot{\alpha})$, assumes GR and has to be replaced in MGCAMB. This was already done in the previous versions, but the contribution of massive neutrinos was neglected. We introduce massive neutrinos properly in the current version, using the following prescription. The ISW term can be written as

[TABLE]

where terms on the right hand side can be computed separately. The first term is determined by taking a derivative of Eq. (16), giving

[TABLE]

Then, $\dot{\Psi}$ can be obtained by differentiating the modified Poisson equation (11),

[TABLE]

where $(\rho a^{2}\Delta)^{\cdot}$ is determined from Eq. (25) and

[TABLE]

completing the set of required equations in a form suitable for implementation in CAMB.

The last modification concerns the implementation of the weak lensing transfer function, since CAMB assumes GR to calculate the Weyl potential $\Phi+\Psi$. Instead, in MGCAMB, we compute $\Psi$ and $\Phi$ from Eqs. (11) and (12). The Weyl potential is used to compute the lensing correlation functions in the analysis presented in Section IV.

III.3 The GR limit of MGCAMB

We have checked the output of MGCAMB in the GR limit, when $\mu=\gamma=1$, for a wide range of neutrino masses and values of the scale factor at which modifications are switched on. Fig. 1 shows the maximum relative difference in $C_{\ell}^{\rm TT}$ and $P(k)$ between CAMB and MGCAMB for $0.05$ eV $\leq\sum m_{\nu}\leq 0.5$ eV and $0.0001\leq a_{\rm trans}\leq 0.1$. In all cases the deviations are below 0.1%. In order to achieve this accuracy for $a_{\rm trans}<0.005$, we had to delay the time at which the RSA in MGCAMB is switched on by a factor of 20, which doubles the running time of the code.

IV Joint constraints on modified growth, massive neutrinos and the dark energy equation of state

The MGCAMB patch has been implemented in the Markov Chain Monte Carlo (MCMC) engine CosmoMC Lewis and Bridle (2002); Lewis (2013), and is called MGCosmoMC333Available at https://github.com/sfu-cosmo/MGCosmoMC. To demonstrate its use, we derive joint constraints of massive neutrinos, modified growth, and the DE equation of state from the datasets included in the current version of CosmoMC. We consider three models: $\Lambda$CDM with the sum of neutrino masses $\sum m_{\nu}$ as an additional parameter (hereafter referred to as Model 0), a model with $\mu$ and $\gamma$ varied along with $\sum m_{\nu}$ and $\Lambda$ playing the role of dark energy (Model 1), and a model with the DE equation of state varying in addition to $\sum m_{\nu}$, $\mu$ and $\gamma$ (Model 2).

In Models 1 and 2, we adopt the ($\mu,\gamma$) parameterization used by the Planck collaboration (Ade et al., 2016b; Aghanim et al., 2018), i.e.

[TABLE]

where $\Omega_{\rm DE}(a)=\rho_{\rm DE}/\rho_{\rm tot}$. We present the results in terms of the derived quantities $\mu_{0}\equiv\mu(a=1)$, $\gamma_{0}\equiv\gamma(a=1)$ and $\Sigma_{0}\equiv\mu_{0}(1+\gamma_{0})/2$. In Model 2, we adopt the CPL parameterization (Chevallier and Polarski, 2001; Linder, 2003) of the DE equation of state,

[TABLE]

In all cases, we also varied the six “vanilla” $\Lambda$CDM parameters, running four parallel chains until the Gelman/Rubin convergence statistics reached $R<0.01$.

IV.1 The datasets

Our analysis made use of the Planck CMB temperature and polarization anisotropy spectra in combination with other datasets. Although the 2018 results were recently released Aghanim et al. (2018), the latest Planck likelihood code was not available at the time of writing, so we used the 2015 version Aghanim et al. (2016). Specifically, we used the Planck 2015 TT, TE and EE likelihood for multipoles in the range $30\leq\ell\leq 2508$ along with the lowTEB polarization likelihood for multipoles in the range $2\leq\ell\leq 29$. We also used the Planck CMB lensing measurements from the minimum variance combination of temperature and polarization with the conservative cut444We use the files smicadx12_Dec5_ftl_mv2_ndclpp_p_teb_consext8 of $40\leq\ell<400$.

We have combined CMB data with the Type Ia supernovae luminosity distance measurements from JLA (Betoule et al., 2014), the BAO measurements from the 6dF galaxy survey (Beutler et al., 2011b) (at $z=0.106$), the SDSS DR7 Main Galaxy Sample (MGS) (Ross et al., 2015) (at $z=0.15$), the measurements of the Hubble parameter $H(z_{i})$, the angular-diameter distance $d_{A}(z_{i})$ and the redshift space distortion measurements of $f(z_{i})\sigma_{8}(z_{i})$ at $z_{i}=\{0.38,0.51,0.61\}$ provided by the BOSS DR12 (Alam et al., 2017). As usually done in the literature, we assumed that the 6DF and MGS measurements are independent from BOSS DR12.

We have also used the recent galaxy clustering and weak lensing measurements from the DES Year 1 results (Abbott et al., 2018a). This dataset consists of the measurements of the angular two-point correlation functions of galaxy clustering, cosmic shear and galaxy-galaxy lensing in a set of 20 logarithmic bins of angular separation in the range $2.5^{\prime}-250^{\prime}$. Since the MG formalism has no nonlinear prescription for structure formation, the angular separations probing the nonlinear scales were properly removed. To do so, we adopted the same method as in (Abbott et al., 2018b) and used the “standard” data DES cutoff as described in Appendix B. Moreover, we modified the DES likelihood code in order to compute the theoretical predictions of the cosmic shear and the galaxy-galaxy lensing correlation functions using the Weyl potential $k^{2}(\Phi+\Psi)/2$ instead of using the GR approximation $k^{2}(\Phi+\Psi)/2=k^{2}\Phi=-(3/2)\Omega_{m}H_{0}^{2}a^{-1}\delta_{m}$, as described in Appendix C. Finally, the covariance between the DES data and the 6DF, MGS and BOSS measurements is ignored, as the observations are carried on different sky patches (Abbott et al., 2018b).

IV.2 The GR limit consistency check

To assess the impact of the small systematic errors introduced by the approximations used in MGCAMB, we performed a consistency check of the GR limit of MGCosmoMC by comparing the results of three MCMC runs: 1) using the original CosmoMC code, 2) using MGCosmoMC with $\mu=\gamma=1$ and $a_{\rm trans}=0.01$ and 3) using MGCosmoMC with $\mu=\gamma=1$ and $a_{\rm trans}=0.001$. In all runs, we varied the six vanilla $\Lambda$CDM parameters and the mass of neutrinos. The results are summarized in Table 1. We can see that the best fit values and the confidence intervals for cosmological parameters are practically the same in all cases and, hence, the results are consistent.

IV.3 Results

The joint constraints derived on massive neutrinos, modified growth and the DE equation of state are summarized in Table 2 for Models 0, 1 and 2 defined at the beginning of this Section. Also, Fig. 2 shows the marginalized distributions of the relevant parameters, with their $\Lambda$CDM limits shown with dashed grey lines. We find the 95% CL bound on massive neutrinos to be

[TABLE]

Our bound on the neutrino mass for Model 0 is comparable to the DES Year 1 result of $0.29$ eV at 95% CL (Abbott et al., 2018a). The use of the CMB polarization data at high-$\ell$ in our analysis is the reason for the stronger constraint.

In Model 1, the effective (cosmological) Newton’s constant can vary at late times. Such variation can happen, for example, in scalar-tensor theories of gravity, where it would generally be scale-dependent. In such theories, the extra Yukawa force mediated by the scalar gravitational degree of freedom enhances the structure formation at scales below the Compton wavelength of the scalar field, which could negate the free streaming suppression of structure formation due to the non-zero neutrino mass. Since in our analysis we considered a scale-independent parameterization of $\mu,\gamma$, this degeneracy between $\sum m_{\nu}$ and $\mu_{0}-1$ is not present and the constraint on $\sum m_{\nu}$ is comparable to the one in Model 0.

In Model 2, $\mu$ is also scale-independent, however the DE density is time-dependent and the degeneracy between the dynamics of DE and the neutrino mass weakens the constraints on $\sum m_{\nu}$. The degeneracy between the neutrino mass and the CPL parameters $w_{0}$ and $w_{a}$ is evident from Fig. 2. The 95% C.L. bounds on the modified growth parameters are consistent with the $\Lambda$CDM limit and with the results obtained by DES (Abbott et al., 2018b). Note that, as expected, the bounds on $\gamma_{0}$ are generally weaker than those on $\mu_{0}$ and $\Sigma_{0}$, because there is no observable that can cleanly separate its effect from the latter two. Our constraints on the variation in the DE equation of state also indicate a good agreement with the $\Lambda$CDM model.

The analysis presented in this Section illustrates how the new MGCAMB patch allows one to derive simultaneous constraints on the neutrino mass, MG and the DE equation of state. The new patch also makes it easy to implement alternative parameterizations of the MG functions and the DE equation of state.

V Discussion

This paper presents a significant update of MGCAMB that features a consistent implementation of massive neutrinos and dynamical dark energy, as well as a new structure that renders the implementation of custom models easier. The new version also has no restriction on the value of the transition time at which the modifications to the linearized Einstein equations are switched on.

MGCAMB was the first publicly released modified Boltzmann solver for cosmological tests of gravity. Since its introduction in 2008 Zhao et al. (2009), MGCAMB has been used in over 100 works. A number of other codes have been introduced since, most notably ISiTGR Dossett et al. (2011), EFTCAMB Hu et al. (2014); Raveri et al. (2014) and hi_class Zumalacrregui et al. (2016). Of them, ISiTGR is close to MGCAMB in its spirit, also introducing phenomenological modifications of equations of motion using two functions $Q$ and $R$ defined in our Appendix A.2. EFTCAMB is based on the effective description of the background and perturbations solutions in general scalar-tensor theories Bloomfield et al. (2013); Gleyzes et al. (2013); Bloomfield (2013), while hi_class uses an alternative effective description of perturbations in scalar-tensor theories on a fixed background Bellini and Sawicki (2014).

MGCAMB is best suited for model-independent constraints on $\mu$ and $\Sigma$, sometimes referred to as $G_{\rm matter}$ and $G_{\rm light}$, which are closely related to observables. The choice of the parameterizated forms of $\mu$ and $\Sigma$ can be informed by the QSA limit of particular types of modified gravity theories Pogosian and Silvestri (2016). One can also perform non-parametric reconstructions of $\mu$ and $\Sigma$ aided by a prior covariance derived from ensembles of modified gravity theories. Such priors can be obtained with the help of EFTCAMB as in Ref. Espejo et al. (2018).

With the present update, MGCAMB should remain a useful tool for cosmological tests of gravity, offering accuracy appropriate for data expected from the next generation surveys such as Euclid Laureijs et al. (2011) and LSST Abell et al. (2009).

Acknowledgements.

We thank Alireza Hojjati for co-developing and maintaining the earlier versions of MGCAMB. We also thank Marco Raveri, Matteo Martinelli and Meng-Xiang Lin for useful discussions. The work of AZ and LP is supported in part by the National Sciences and Engineering Research Council (NSERC) of Canada. AS acknowledges support from the NWO and the Dutch Ministry of Education, Culture and Science (OCW), and also from the D-ITP consortium, a program of the NWO that is funded by the OCW. GBZ is supported by NSFC Grants 1171001024 and 11673025, and the National Key Basic Research and Development Program of China (No. 2018YFA0404503).

Appendix A Other parameterizations

A.1 The $\mu,\Sigma$ parameterization

As mentioned in Sec. II, rather than working with $\mu$ and $\gamma$, it can be beneficial to constrain $\mu$ and $\Sigma$, with the latter defined in Eq. (13). In this version of MGCAMB, we implement ($\mu,\Sigma$) by mapping it onto ($\mu,\gamma$) using $\gamma=2\Sigma-\mu$, which agrees with Eq. (13) in the limit of negligible anisotropic stress ($\sigma\rightarrow 0$). In other words, we define $\Sigma$ as

[TABLE]

The circumstances in which the difference between this definition and the one in Eq. (13) can be important are not entirely clear to us. If necessary, it is relatively easy to add to MGCAMB a separate set of equations for the ($\mu,\Sigma$) parameterization based on Eq. (13).

A.2 The $Q,R$ parameterization

Another phenomenological parameterization was introduced in (Bean and Tangmatitham, 2010) in which modifications of gravity are encoded in functions $Q$ and $R$ defined through

[TABLE]

This parameterization is consistently implemented in MGCAMB. The corresponding equation for $\dot{\eta}$ is

[TABLE]

where the factor of $8\pi G$ is absorbed into $\rho$ and $P$. For the ISW effect, we replace Eq. (35) with

[TABLE]

Appendix B Linear DES data

Since the phenomenological parameterization implemented in MGCAMB has no prescription for nonlinear structure formation, in order to use the DES data, we remove the nonlinear data in the same way as was done in (Ade et al., 2016b; Abbott et al., 2018b). We define

[TABLE]

where ${\bf t}_{\rm NL}$ and ${\bf t}_{\rm L}$ represent the data vector containing the nonlinear and linear theory predictions, respectively, in the $\Lambda$CDM best-fit model. The nonlinear predictions are obtained using the Halofit model present in the default CAMB. We then find the data point that contributes the most to $\Delta\chi^{2}$ and remove it. We then repeat the procedure until $\Delta\chi^{2}$ is less than a threshold. We arbitrarily define three set of cuts on the data: a “soft” cut where the $\Delta\chi^{2}$ threshold is 10, a “standard” cut with the threshold set to 5 and finally an “aggressive” cut with $\Delta\chi^{2}<1$. The number of data points removed are 88, 118 and 178, respectively. As an example, in Fig. 3 we show the “standard” cut applied to the DES galaxy-galaxy angular correlation function $w_{ij}(\theta)$, where $i,j$ label the redshift bins. The blue and orange lines represent the nonlinear and linear theoretical predictions and the grey shaded lines show the data which is excluded by the above method.

Appendix C Galaxy Clustering - Weak lensing theory with Weyl potential.

Here we describe the modifications to the DES likelihood required for evaluating the weak lensing observables $\xi_{+}$, $\xi_{-}$ and $\gamma_{t}$. The standard DES likelihood in CosmoMC assumes $\Phi+\Psi=2\Phi$ and then relates the potential $\Phi$ to the density perturbation $\delta$ using the Poisson equation,

[TABLE]

In MGCAMB, the Poisson equation is modified and $\Phi=\Psi$ does not hold. Hence, we compute the Weyl potential power spectrum directly,

[TABLE]

where the tilde quantities are the transfer functions at redshift $z$ and $P_{\cal R}(k)$ is the primordial power spectrum. The cosmic shear correlations $\xi_{\pm}$ are then given by

[TABLE]

where $J_{0/4}(x)$ is the spherical Bessel function of order zero (fourth), $P_{\kappa}$ (in the Limber approximation) is given by

[TABLE]

where $q^{i}(\chi)$ is the lensing efficiency function,

[TABLE]

$n^{i}$ denotes the effective number density of galaxies normalized to one and the Weyl power spectrum is evaluated using the linear theory only. The nonlinear data is removed according to the procedure explained in App. B.

Similarly, to calculate the tangential shear of background galaxies around foreground galaxies, we define the Weyl-matter power spectrum as

[TABLE]

The tangential shear is then given by

[TABLE]