## ISiTGR : Testing deviations from GR at cosmological scales including dynamical dark energy, massive neutrinos, functional or binned parametrizations, and spatial curvature

Cristhian Garcia-Quintero, ∗ Mustapha Ishak, † Logan Fox, and Jason Dossett Department of Physics, The University of Texas at Dallas, Richardson, Texas 75080, USA (Dated: December 10, 2019)

We introduce a new version of the Integrated Software in Testing General Relativity ( ISiTGR ) which is a patch to the software CAMB and CosmoMC . ISiTGR is intended to test deviations from GR at cosmological scales using cosmological data sets. While doing so, it allows for various extensions to the standard flat ΛCDM model. In this new release, we have support for the following: 1) dynamical dark energy parametrizations with a constant or time-dependent equation of state; 2) a consistent implementation of anisotropic shear to model massive neutrinos throughout the full formalism; 3) multiple commonly-used parametrizations of modified growth (MG) parameters; 4) functional, binned and hybrid time- and scale-dependencies for all MG parameters; 5) spatially flat or curved backgrounds. ISiTGR is designed to allow cosmological analyses to take full advantage of ongoing and future surveys to test simultaneously or separately various extensions to the standard model. We describe here the formalism and its implementation in the CMB code, the Integrated SachsWolfe (ISW) effect, and the 3x2 point statistics. Next, we apply ISiTGR to current data sets from Planck-2018, Planck-2015, Dark Energy Survey YR1 release, Baryonic Acoustic Oscillations (BAO), Redshift Space Distortions (BAO/RSD) from the BOSS Data Release 12, the 6DF Galaxy Survey and the SDSS Data Release 7 Main Galaxy Sample, and Supernova from the Pantheon compilation, joint SNLS/SDSS data analysis and the Hubble Space Telescope. We derive constraints on MG parameters for various combinations of the five features above and find that GR is consistent with current data sets in all cases. The code is made publicly available at https://github.com/mishakb/ ISiTGR .

PACS numbers: 95.36.+x,98.80.Es,04.50.Kd

∗ Electronic address:

† Electronic address:

## Contents

| I. Introduction   | I. Introduction                                                                   |   3 |
|-------------------|-----------------------------------------------------------------------------------|-----|
| II.               | Growth equations including anisotropic shear and spatial curvature                |   4 |
|                   | A. Growth Equations                                                               |   4 |
|                   | B. Massive Neutrino contributions                                                 |   5 |
| III.              | Effective Dynamical Dark Energy parametrization and evolution                     |   5 |
| IV.               | Modified Growth Equations and MG parameters including time and scale dependencies |   6 |
|                   | A. The ( µ ( a, k ) , γ ( a, k )) parametrization                                 |   6 |
|                   | B. The ( µ ( a, k ) , Σ( a, k )) parametrization                                  |   6 |
|                   | C. The ( Q ( a, k ) , R ( a, k ) , D ( a, k )) parametrization                    |   6 |
|                   | D. MG parameter functional dependencies in time and scale                         |   7 |
|                   | E. MG parameter binned dependencies in time and scale                             |   9 |
| V.                | Implementation of modified gravity equations in ISiTGR patch to CAMB and CosmoMC  |  10 |
|                   | A. CMB implementation and the synchronous gauge                                   |  10 |
|                   | B. CAMB variables and potential evolution                                         |  10 |
| VI.               | Modified power spectra and correlation functions                                  |  12 |
|                   | A. 3 × 2 point power spectra and correlation functions                            |  12 |
|                   | B. Integrated Sachs-Wolfe (ISW) effect                                            |  13 |
| VII.              | Applications and results from current available data sets                         |  14 |
|                   | A. Data sets                                                                      |  14 |
|                   | B. Constraints and correlations                                                   |  14 |
|                   | Summary and concluding remarks                                                    |  21 |
|                   | Acknowledgments                                                                   |  21 |
| A.                | Spatially flat case                                                               |  21 |
|                   | 1. Formalism and MG equations                                                     |  21 |
|                   | 2. Implementation                                                                 |  22 |
| B.                | Parametrization summary table                                                     |  25 |
| C.                | Flowcharts for the ISiTGR code                                                    |  26 |
|                   | References                                                                        |  27 |

## I. INTRODUCTION

Continuous progress is being made toward precision cosmology with a number of ongoing and planned surveys and missions [1-10] 1 The resulting complementary and precise observations have opened the door to testing gravity physics (General Relativity (GR)) at cosmological scales, see, e.g., the reviews [11-15].

One of the chief motivations to test GR at cosmological scales is the pressing question of cosmic acceleration and the dark energy associated with it, see e.g. [16-22]. However, testing GR at large scales is a well-motivated and justified objective in its own right.

There are, at least, two primary routes to testing gravity at large scales. The first is to model departure from GR in a phenomenological way by adding parameters that would signal such a deviation at the level of the growth of large-scale structure in the Universe. Such parameters would take some expected values in GR, often one or zero, but will deviate from them otherwise. Interestingly, different models of gravity that have the same expansion history can still exhibit distinct growth rates of large-scale structures, which can be used as a discriminant between gravity theories, see e.g. [23-34]. The second approach is to develop analysis pipelines and simulations specific to some proposed modified gravity theories such as the well-known f ( R ), the DGP (Dvali-Gabadadze-Porrati) [35], or other modified gravity models. Understandably, this second approach has been progressing at a slower pace because it requires more involved development and resources. The two methods complement each other in the effort to test gravity at cosmic scales. For further discussion on both approaches, we refer the reader to the following partial list of reviews [11-15] and references therein. In view of recent interest and developments, it is also worth mentioning the method of exploring inconsistencies between data sets as a way to test departures from the standard model, see e.g. [36-40].

In this paper, we introduce and describe a new version of ( ISiTGR ) ( I ntegrated S oftware i n T esting G eneral R elativity ) [41, 42] which is a patch for the widely-used software packages CAMB (Code for Anisotropies in the Microwave Background [43]) and CosmoMC (Cosmological Monte Carlo [44]). ISiTGR follows the first approach described above to constrain departure from GR based on various pairs of modified growth (MG) parameters. These parametrize the strength of the coupling between the gravitational potentials and the spacetime sources, as well as the relationships between the two potentials, as we describe further below in the paper.

In this new release of ISiTGR , we included contributions from anisotropic shear stress throughout all of the formalism to consistently account for contributions from massive neutrinos and radiation. We also included dynamical dark energy for the background with constant or time-varying equations of state. Furthermore, we expanded the support for several pairs of existing MG parametrizations as needed by various types of cosmological probes and surveys. These features have been made to work consistently in a spatially flat or curved background. Additionally, we implement functional as well as binned methods for the time and scale dependencies of MG parameters. In this way, the new version of ISiTGR has been designed to suit the needs of analyses that intend to test various aspects of extended models using incoming and future data sets, and makes it possible to constrain such extensions separately or simultaneously . ISiTGR has been cited or used in over 50 papers and has been applied to CFHTLens, KidS-450, 2dF and Planck data, see for example [45-47].

There are over a dozen other codes that test deviation from GR or specific MG models at cosmological scales. These include for example MGCAMB [30, 48] that is similar to ISiTGR and is built on the top of CAMB ; hi class [49] that is built on CLASS [50, 51] and based on the Horndeski models (and beyond) [52]; EFTCAMB [53, 54] which follows an approach inspired by Effective Field Theory perturbations applied to dark energy; and EoS class that is based on the equation of state approach applied to Horndeski models [55]. We refer the reader for overviews of codes in Refs. [15, 56].

The paper is organized as follows. In section II, we describe the growth equations in a flat or curved background with anisotropic shear stress. Effective dynamical dark energy evolution with constant or time-varying equations of state is summarized in section III. In section IV, we describe modified growth (MG) equations and various MG parameters as used in various surveys and probes. We describe there as well the time and scale dependencies for MG parameters as implemented in functional and binned forms. Section V is a brief overall description of how the modifications have been implemented in the code, while section VI describes power spectra and correlation functions for the observables. In section VII, we apply ISiTGR to current available data sets to constrain MG parameters in different forms and backgrounds. We summarize in section VIII.

1 e.g. Kilo-Degree Survey (KiDS) [1], Dark Energy Survey (DES) [2], Hyper Suprime-Cam (HSC) [3], Large Synoptic Survey Telescope (LSST) [4], Dark Energy Spectroscopic Instrument (DESI) [5], Euclid [6], Square Kilometre Array (SKA) [7], Wide Field Infrared Spectroscopic Telescope (WFIRST) [8], Simons Observatory [9], CMB-S4 [10] and many others.

## II. GROWTH EQUATIONS INCLUDING ANISOTROPIC SHEAR AND SPATIAL CURVATURE

## A. Growth Equations

As discussed above, a modification to GR at cosmological scales can affect the growth rate of large scale structure. This can be phenomenologically modeled by changing the linearly perturbed Einstein equations. Specifically, one can change the resulting Poisson-like equations involving the gravitational scalar potentials. Let us first briefly review the key equations in the GR case where we allow for spatial hyper-surfaces to be flat or curved and we keep the anisotropic shear stress coming from, for example, the Neutrino sector.

We start with the perturbed FriedmannLemaˆtreRobertsonWalker (FLRW) metric written in the general conformal Newtonian gauge given by

$$d s ^ { 2 } = a ( \tau ) ^ { 2 } [ - ( 1 + 2 \Psi ) d \tau ^ { 2 } + ( 1 - 2 \Phi ) \gamma _ { i j } d x ^ { i } d x ^ { j } ], \quad \quad \quad ( 1 ) \\ \cdot \quad \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot$$

where Φ and Ψ are scalar gravitational potentials describing the scalar mode of the metric perturbations, x i 's are the comoving coordinates, τ is conformal time, and a τ ( ) is the scale factor. γ ij is the 3-dimensional metric, which can be written in the Cartesian coordinates ( x, y, z ) as [57]

$$\gamma _ { i j } = \delta _ { i j } \left [ 1 + \frac { K } { 4 } \left ( x ^ { 2 } + y ^ { 2 } + z ^ { 2 } \right ) \right ] ^ { - 2 },$$

where K = -Ω k H 2 0 is the spatial curvature and we use units where c = 1. H 0 is the Hubble parameter (constant) today.

In a non-flat FLRW Universe the Fourier modes can be generalized as eigen-functions, G , of the Laplacian operator such that ∇ 2 G glyph[vector] ( k, glyph[vector] x ) = -k G glyph[vector] 2 ( k, glyph[vector] x ). One can then expand perturbations in terms of the eigen-function G and its spatial covariant derivatives, see e.g. [58, 59].

The first order perturbed Einstein equation gives two equations. The first equation derives from the combination of the time-space and time-time perturbed equations and gives a Poisson equation relating the potential, Φ, and the gauge-invariant, rest-frame over-density, ∆ . i The second equation derives from the traceless, space-space component of the equations and gives a relation between the two metric potentials involving the shear stress, σ i (where i denotes a particular matter species). The two equations read:

$$\left ( k ^ { 2 } - 3 K \right ) \Phi \ = \ - 4 \pi G a ^ { 2 } \sum _ { i } \rho _ { i } \Delta _ { i },$$

i

$$k ^ { 2 } ( \Psi - \Phi ) \ = \ - 1 2 \pi G a ^ { 2 } \sum _ { i } \rho _ { i } ( 1 + w _ { i } ) \sigma _ { i },$$

where ρ i is the density for matter species i .

The gauge-invariant, rest-frame overdensity, ∆ , is a key quantity since its evolution describes the growth of inhoi mogeneities (structures) in the Universe. It is defined as

$$\Delta _ { i } & = \delta _ { i } + 3 \mathcal { H } \frac { q _ { i } } { k },$$

where H = ˙ a/a is the Hubble factor in conformal time; δ i = δρ /ρ i ¯ is the fractional overdensity; and q i is the heat flux related to the divergence of the peculiar velocity, θ i , by θ i = k q i 1+ w .

i From conservation of the energy momentum tensor of the perturbed matter fluids, the fractional overdensity and heat flux for uncoupled fluid species or the mass-averaged quantities for all the fluids evolve as [60]:

$$\dot { \delta } \ = \ - k q + 3 ( 1 + w ) \dot { \Phi } + 3 \mathcal { H } ( w - \frac { \delta P } { \delta \rho } ) \delta, \\.$$

$$\frac { \dot { q } } { k } \ = \ - \mathcal { H } ( 1 - 3 w ) \frac { q } { k } + \frac { \delta P } { \delta \rho } \delta + \left ( 1 + w \right ) \left ( \Psi - \sigma \right ),$$

where w = p/ρ is the equation of state of the fluid. Next, combining the two equations above, one can express the evolution of ∆ (or ∆ ) by i

$$\dot { \Delta } = 3 ( 1 + w ) \left ( \dot { \Phi } + \mathcal { H } \Psi \right ) + 3 \mathcal { H } w \Delta - \left [ k ^ { 2 } + 3 \left ( \mathcal { H } ^ { 2 } - \dot { \mathcal { H } } \right ) \right ] \frac { q } { k } - 3 \mathcal { H } ( 1 + w ) \sigma.$$

Combining the growth equations (3) and (4), along with the evolution equations (8) and a τ ( ), the growth history of large scale structures in the Universe can be fully described.

## B. Massive Neutrino contributions

When considering consistently the contributions from massive neutrinos to MG equations, one needs to include the anisotropic shear stress throughout the evolution equations. In CAMB it is defined as

$$\Pi _ { i } = \frac { 3 } { 2 } ( 1 + w _ { i } ) \sigma _ { i }.$$

In previous implementations of MG codes [30, 41], the anisotropic shear stress contributions related with Π for ˙ massive neutrinos were usually neglected. In view of recent interest and developments in the neutrino sector, we have consistently included the anisotropic shear stress contributions at all levels and calculations involving MG modifications. These are described in section V.

## III. EFFECTIVE DYNAMICAL DARK ENERGY PARAMETRIZATION AND EVOLUTION

We model further below the time evolution of MG parameters via the effective dark energy density parameter's dependence on the scale factor. This is given by:

$$\Omega _ { D E } ( a ) \equiv \rho _ { D E } ( a ) / \rho _ { c } ( a ).$$

Once the effective equation of state for the dark energy is specified, we can compute ρ DE ( a ), and thus obtain an expression for Ω DE ( a ):

$$\Omega _ { D E } ( a ) = \Omega _ { \Lambda } \left ( \frac { H _ { 0 } } { H } \right ) ^ { 2 } \exp \left [ - 3 \int _ { a _ { 0 } } ^ { a } [ 1 + w ( a ^ { \prime } ) ] \frac { d a ^ { \prime } } { a ^ { \prime } } \right ].$$

We solve for this dark energy parameter for four different models: standard ΛCDM model; w CDM model, the Chevallier-Polarski-Linder (CPL) parametrization (see [61, 62]); and the pivot equation of state for the CPL parametrization described in [63]. The corresponding solutions to Eq. (11) for each of these models along with the corresponding dark energy densities are given in Table (I).

We note that a pivot dark energy equation of state based on the CPL parametrization is implemented in ISiTGR . This is done through a linear transformation giving w = w p +( a p -a w ) a , where a p is the scale factor related to a pivot redshift z p = 1 /a p -1. The pivot scale factor is determined by [64]

$$a _ { p } = 1 + C ( w _ { 0 }, w _ { a } ) / C ( w _ { a }, w _ { a } ),$$

where C represents the covariance matrix. The choice of a p is such that w p and w a are decorrelated so it minimizes the uncertainty in w a ( ), but user should be aware of some limitations or misinterpretations expressed in, e.g. see [65].

TABLE I: Dark energy parametrizations given in ISiTGR : ΛCDM, w CDM, the ( w 0 , w a ) parametrization and a pivot dark energy equation of state ( w p , w a ). Once the user chooses one of these models, ISiTGR selects the corresponding dark energy density and the corresponding dark energy evolution.

| Dark Energy model                                          | Dark Energy density                                           | Dark Energy evolution                                                                                                                                                                      |
|------------------------------------------------------------|---------------------------------------------------------------|--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| w = - 1 w = w 0 ( a ) = w 0 +(1 - a ) a ) = w p +( a p - a | ρ v ( t ) = ρ (0) v (0) - 3(1+                                | Ω DE ( a ) = Ω v ( H 0 H ) 2 Ω DE ( a ) = Ω v ( H 0 H ) 2 a - 3(1+ w 0 ) ( a ) = Ω v ( H 0 H ) 2 a - 3(1+ w 0 + w a ) e 3 w a ( a ( a ) = Ω v ( H 0 ) 2 a - 3(1+ w p + a p w a ) e 3 w a ( |
|                                                            | ρ v ( t ) = ρ v a w 0 )                                       |                                                                                                                                                                                            |
| w w a                                                      | ρ v ( t ) = ρ (0) v a - 3(1+ w 0 + w a ) e 3 w a ( a - 1)     | Ω DE - 1)                                                                                                                                                                                  |
| w ( ) w a                                                  | ρ v ( t ) = ρ (0) v a - 3(1+ w p + a p w a ) e 3 w a ( a - 1) | Ω DE H a - 1)                                                                                                                                                                              |

## IV. MODIFIED GROWTH EQUATIONS AND MG PARAMETERS INCLUDING TIME AND SCALE DEPENDENCIES

As usual, we implement deviations from GR using modifications of the linearly perturbed Einstein equations and the resulting modified versions of the metric potential equations (3) and (4). Again, to implement ISiTGR for both spatially flat and curved backgrounds, we do not limit the formalism to the K = 0 -case (However, the explicit equations for the flat case are given in Appendix A). We also keep the shear terms non-zero throughout the equations. Further below we provide the functional dependencies on time (scale factor) and scale, and allow for a time-dependent equation of state for dark energy.

## A. The ( µ a, k ( ) , γ a, k ( )) parametrization

The first parameter ( µ a, k ( )) enters the two MG equations below, (13) and (14), and quantifies the strength of the gravitational coupling between the potentials and the sources. The second parameter, γ a, k ( ), (sometimes also noted as η a, k ( )), is called the slip parameter [26] and quantifies the difference between the two gravitational potentials. At late times when anisotropic shear can be assumed to be negligible then γ a, k ( ) ≡ Φ Ψ. The modified equations read: /

$$( k ^ { 2 } - 3 K ) \Psi = - 4 \pi G a ^ { 2 } \mu ( a, k ) \sum _ { i } \left [ \rho _ { i } \Delta _ { i } + 3 \left ( \frac { k ^ { 2 } - 3 K } { k ^ { 2 } } \right ) \rho _ { i } ( 1 + w _ { i } ) \sigma _ { i } \right ]$$

and

$$k ^ { 2 } ( \Phi - \gamma ( a, k ) \Psi ) = 1 2 \pi G a ^ { 2 } \mu ( a, k ) \sum _ { i } \rho _ { i } ( 1 + w _ { i } ) \sigma _ { i }.$$

GR is recovered when µ a, k ( ) and γ a, k ( ) are equal to unity. This parametrization has been used in for example [66, 67] with evolutions that we describe further below.

## B. The ( µ a, k ( ) , Σ( a, k )) parametrization

Some observables such as gravitational lensing or the Integrated Sachs-Wolfe (ISW) effect involve the Weyl potential,

$$\Phi _ { W } & = ( \Phi + \Psi ) / 2,$$

which governs the motion of light-like particles. Taking advantage of this, (3) and (4) can be combined, and another parameter, Σ( a, k ), can be defined to directly probe modifications to the Weyl potential (see e.g. [68]). Consequently, this parameter can be defined through

$$k ^ { 2 } ( \Phi + \Psi ) = - 4 \pi G a ^ { 2 } \Sigma ( a, k ) \sum _ { i } \left [ \frac { 2 \rho _ { i } \Delta _ { i } } { 1 - 3 K / k ^ { 2 } } + 3 \rho _ { i } ( 1 + w _ { i } ) \sigma _ { i } \right ].$$

Deviation from GR is now measured using Eqs. (13) and (16). Moreover, Σ( a, k ) also takes the value of unity in the GR case. In the case of negligible or zero shear, it follows that

$$\Sigma ( k, a ) = \frac { \mu ( k, a ) [ 1 + \gamma ( k, a ) ] } { 2 }.$$

This parametrization or a similar one have been used in, for example, [47, 68].

## C. The ( Q a, k , ( ) R a, k , ( ) D a,k ( )) parametrization

Another parametrization similar to the above was introduced in [32] with a parameter Q a, k ( ) to characterize the gravitational strength and R a, k ( ) as a gravitational slip parameter. The MG equations read:

$$( k ^ { 2 } - 3 K ) \Phi = - 4 \pi G a ^ { 2 } Q ( a, k ) \sum _ { i } \rho _ { i } \Delta _ { i }$$

and

$$k ^ { 2 } ( \Psi - R ( a, k ) \Phi ) = - 1 2 \pi G a ^ { 2 } Q ( a, k ) \sum _ { i } \rho _ { i } ( 1 + w _ { i } ) \sigma _ { i }.$$

Similarly, at late times when anisotropic stress is negligible, Ψ = R Φ.

Again, it is also possible to combine the two equations above and define the parameter D a,k ( ) = Q a, k ( )(1 + R a, k ( )) / 2 to be used instead of R a, k ( ). This not only avoids a strong degeneracy between Q a, k ( ) and R a, k ( ), but also gives a parameter which can be directly probed by lensing and ISW observations. The combined equation reads

$$k ^ { 2 } ( \Phi + \Psi ) = \left ( \frac { - 8 \pi G a ^ { 2 } } { 1 - 3 K / k ^ { 2 } } \right ) D ( a, k ) \sum _ { i } \rho _ { i } \Delta _ { i } - 1 2 \pi G a ^ { 2 } Q ( a, k ) \sum _ { i } \rho _ { i } ( 1 + w _ { i } ) \sigma _ { i }.$$

Then, for this parametrization, one can use (18) and (20) instead of (19). We implemented both pairs in ISiTGR . Here it is important to note that Σ( a, k ) = D a,k ( ) only in the zero anisotropic shear stress case, and that while D a,k ( ) ≡ Q a, k ( )(1 + R a, k ( )) / 2 is a definition, the relationship Σ( a, k ) = µ a, k ( )(1 + γ a, k ( )) / 2 also holds only in the zero anisotropic shear case.

## D. MG parameter functional dependencies in time and scale

We implement in ISiTGR the dependencies of each set of MG parameters on both time and scale. For the time evolution of MG parameters, we mainly use the dark energy density time evolution, described by the parameter Ω DE ( a ), as shown further below for various pairs of parameters. In this way, the contribution to clustering and anisotropic stress by MG effects is proportional to their effective dark energy density. One needs to ensure though that the time evolution of the effective dark energy density is governed by its corresponding effective equation of state, as we describe in section III. Note that some limitations and failures of this proportionality of MG parameters with the parameter Ω DE ( a ) have been reported and discussed in [69] and [70]. It was reported there that such parametrization can miss deviations from GR and the proportionality does not hold for all cosmic evolution including the epoch of interest of cosmic acceleration. Then, the user should use these with some caution and complement them with binning methods.

It was shown in [71] and other works that scale dependence can provide further insights to discriminate between gravity theories. For this dependence, we implemented the function used in, for example, [66] and [72], with the additional factor

$$S _ { i } ( a, k ) = \frac { 1 + c _ { i } \left ( \lambda H ( a ) / k \right ) ^ { 2 } } { 1 + \left ( \lambda H ( a ) / k \right ) ^ { 2 } },$$

where the index i = 1 2 stands for the first and the second MG parameters and , H a ( ) is the Hubble parameter. The form of S i ( a, k ) is such that at small scales (large k ) we have S i → 1, while at large scales (small k ) we get S i → c i . Therefore, at small scales (21) has no effect on the MG parameters, while for large scales it makes MG parameters proportional to c 1 and c 2 . Therefore S i provides information about how MG parameters evolve at large scales. Moreover, each MG parameter may have a different scale dependent evolution, since c 1 is independent of c 2 (recovering GR when c 1 = c 2 = 1, regardless of the value of λ ). However, even if c i ∼ 0, the λ parameter will still provide some scale dependence. In this case, we have that S i → 1 for large k , but for small k we observe that S i → 0. Therefore, if c i is negligible, then at large scales the MG parameters will go to their GR limit value, but at small scales the MG parameters still evolve in time. Furthermore, we can recover GR if λ = 0 regardless of the values that c 1 and c 2 can take.

Next, we apply these functional dependencies to each set of MG parameters.

## 1. Time and scale functional dependencies for ( µ a, k ( ) , γ a, k ( ))

For the ( µ , γ ) case, we implemented the parametrization as given in the recent Planck analyses [66, 67]. The MG parameters are implemented explicitly as

$$\mu ( a, k ) = 1 + E _ { 1 1 } \Omega _ { D E } ( a ) \left [ \frac { 1 + c _ { 1 } \left ( \lambda H ( a ) / k \right ) ^ { 2 } } { 1 + \left ( \lambda H ( a ) / k \right ) ^ { 2 } } \right ]$$

and

$$\gamma ( a, k ) = 1 + E _ { 2 2 } \Omega _ { D E } ( a ) \left [ \frac { 1 + c _ { 2 } \left ( \lambda H ( a ) / k \right ) ^ { 2 } } { 1 + \left ( \lambda H ( a ) / k \right ) ^ { 2 } } \right ].$$

As we mentioned before, the scale dependent implementation is such that, for small k we have that µ → 1+ c 1 E 11 Ω DE and γ → 1 + c 2 E 22 Ω DE . Furthermore, for large k we note that µ → 1 + E 11 Ω DE and γ → 1 + E 22 Ω DE . In all cases, Ω DE becomes negligible at early times (high redshift) so MG parameters go to the GR value of 1. So MG effects here are modeled to be negligible at early times.

## 2. Time and scale functional dependencies for ( µ a, k ( ) , Σ( a, k ))

We implemented here the following time and scale dependencies:

$$\mu ( a, k ) = 1 + \mu _ { 0 } \frac { \Omega _ { D E } ( a ) } { \Omega _ { \Lambda } } \left [ \frac { 1 + c _ { 1 } \left ( \lambda H ( a ) / k \right ) ^ { 2 } } { 1 + \left ( \lambda H ( a ) / k \right ) ^ { 2 } } \right ]$$

and

$$\Sigma ( a, k ) = 1 + \Sigma _ { 0 } \frac { \Omega _ { D E } ( a ) } { \Omega _ { \Lambda } } \left [ \frac { 1 + c _ { 2 } \left ( \lambda H ( a ) / k \right ) ^ { 2 } } { 1 + \left ( \lambda H ( a ) / k \right ) ^ { 2 } } \right ].$$

This time evolution has been used in several recent works including the Dark Energy Survey (DES) [73]. The effective dark energy density is used in the ratio such that Ω DE ( a / ) Ω Λ = 1 today, so the parameters take their GR values today. For the scale dependence, the same reasoning and limits apply as in the previous parametrization. Here and in [73], it is worth noting that GR is recovered when µ 0 = 0 and Σ 0 = 0. Also, in [73], the scale dependence is not modeled and their equations don't have the terms '1' as in the RHS of our equations (24) and (25) because they defined their MG equations (13) and (16) with 1 + µ a, k ( ) and 1 + Σ( a, k ) instead. In other words, the two parametrizations are the same except that we absorbed the terms '1' from their MG equations into our two parameters. This is handled by ISiTGR and leads to exactly the same constraints on MG parameters (see section VII).

## 3. Time and scale functional dependencies for ( Q a, k , ( ) R a, k , ( ) D a,k ( ))

We implemented the same time and scale dependencies for these 3 parameters with the same limiting cases and behaviors. The equations read:

$$Q ( a, k ) = 1 + Q _ { 0 } \frac { \Omega _ { D E } ( a ) } { \Omega _ { \Lambda } } \left [ \frac { 1 + c _ { 1 } \left ( \lambda H ( a ) / k \right ) ^ { 2 } } { 1 + \left ( \lambda H ( a ) / k \right ) ^ { 2 } } \right ],$$

$$R ( a, k ) = 1 + R _ { 0 } \frac { \Omega _ { D E } ( a ) } { \Omega _ { \Lambda } } \left [ \frac { 1 + c _ { 2 } \left ( \lambda H ( a ) / k \right ) ^ { 2 } } { 1 + \left ( \lambda H ( a ) / k \right ) ^ { 2 } } \right ]$$

and

$$D ( a, k ) = 1 + D _ { 0 } \frac { \Omega _ { D E } ( a ) } { \Omega _ { \Lambda } } \left [ \frac { 1 + c _ { 2 } \left ( \lambda H ( a ) / k \right ) ^ { 2 } } { 1 + \left ( \lambda H ( a ) / k \right ) ^ { 2 } } \right ].$$

Again, we can go from the ( Q a, k ( ), R a, k ( )) to the ( Q a, k ( ), D a,k ( )) parametrization by using D a,k ( ) = Q [1 + R a, k ( )] / 2, but we recall that we cannot use Σ( a, k ) = µ a, k ( )[1 + γ a, k ( )( a, k ))] / 2 unless the anisotropic shear stress is zero.

## E. MG parameter binned dependencies in time and scale

Besides the functional forms described above, ISiTGR implements a binning method where the time (redshift) and scale dependencies of the MG parameters are modeled by different parameters in each bin. The advantages of using binning methods have been highlighted in [47, 74, 75].

Moreover, ISiTGR uses two different ways to evolve the MG parameters using binning methods. The first one is the traditional binning method in which one evolves the MG parameters in two different predefined redshift and scale bins. Unlike some treatments, we provide some extra control in the transition between the scale bins by the addition of hyperbolic functions, see e.g. [47, 76, 77]. The second approach incorporated into ISiTGR is the hybrid method, in which we still have two predefined bins for redshift and scale, but we allow an independent monotonic functional evolution for the MG parameters in each bin, this produces a smoother transition between the bins. In the current version of ISiTGR , we implement binning methods for the ( µ γ , ), ( µ ,Σ) and ( Q D , ) parametrizations.

## 1. Traditional binning

Here MG parameters are binned in both redshift, z, and wavenumber (scale), k. A total of four bins are created by using two redshift bins and two scale bins. The scale bins are k ≤ k c and k &gt; k c , while the redshift bins are 0 &lt; z ≤ z div and z div &lt; z ≤ z TGR . For redshifts z &gt; z TGR the MG parameters take their GR value of 1 at all scales. If X a,k ( ) represents any MG parameter in ISiTGR (e.g. µ , γ , Σ, Q , or D ), then the binned form of the parameter is written as

$$X ( a, k ) = \frac { 1 + X _ { z _ { 1 } } ( k ) } { 2 } + \frac { X _ { z _ { 2 } } ( k ) - X _ { z _ { 1 } } ( k ) } { 2 } \tanh \left ( \frac { z - z _ { d i v } } { z _ { t w } } \right ) + \frac { 1 - X _ { z _ { 2 } } ( k ) } { 2 } \tanh \left ( \frac { z - z _ { T G R } } { z _ { t w } } \right ),$$

where z div is the specific redshift at which the transition between the two bins occurs and z TGR is the redshift below which GR is to be tested. In the ISiTGR code, z TGR = 2 z div is hard-coded giving equally sized bins. Also, z tw acts as a transition width for the hyperbolic tangent function which is used in order to make the transition between the bins smooth. It is hard-coded as z tw = 0 05. The binning in scale is implemented via the two parameters . X z 1 ( k ) and X z 2 ( k ) as follows:

$$X _ { z _ { 1 } } ( k ) = \frac { X _ { 2 } + X _ { 1 } } { 2 } + \frac { X _ { 2 } - X _ { 1 } } { 2 } \tanh \left ( \frac { k - k _ { c } } { k _ { t w } } \right )$$

and

$$X _ { z _ { 2 } } ( k ) = \frac { X _ { 4 } + X _ { 3 } } { 2 } + \frac { X _ { 4 } - X _ { 3 } } { 2 } \tanh \left ( \frac { k - k _ { c } } { k _ { t w } } \right ),$$

where k tw is the transition width between k bins and is hard-coded as k tw = k / c 10. Therefore, since k c quantifies the scale dependence, such a value for k tw ensures that the transition between bins will occur early, before the scale at which the transition occurs. However, the user can change this feature easily.

The method defines 4 parameters X i for each parameter in the pairs ( µ a, k ( ), γ a, k ( )), ( µ a, k ( ),Σ( a, k )), ( R a, k ( ), Q a, k ( )) or ( R a, k ( ), D a,k ( )), so there are a total of 8 parameters. As we can observe in Table V, most of these can be constrained at almost the same level of significance as for the functional forms.

## 2. Hybrid binning method

This method keeps the same binning in the redshift as the traditional one described above, but it replaces the binning in scale by a functional form within each redshift bin. In other words, the functions X z 1 ( k ) and X z 2 ( k ) follow a monotonic evolution inside each redshift bin. In this case, the functions X z 1 ( k ) and X z 2 ( k ) are implemented as

$$X _ { z _ { 1 } } ( k ) = X _ { 1 } e ^ { - k / k _ { c } } + X _ { 2 } ( 1 - e ^ { - k / k _ { c } } )$$

and

$$X _ { z _ { 2 } } ( k ) = X _ { 3 } e ^ { - k / k _ { c } } + X _ { 4 } ( 1 - e ^ { - k / k _ { c } } ).$$

Of course, the exponential functional form can be replaced by other functions. As shown in the original ISiTGR paper, using hybrid parametrization produces a smoother matter power spectrum.

A summary table of the parametrizations covered by ISiTGR is given in Table VI in Appendix B. An illustrating diagram for the various parametrizations (see Fig. 9), as well as a flowchart that shows the parameter files in ISiTGR (see Fig. 10) are given in Appendix C.

## V. IMPLEMENTATION OF MODIFIED GRAVITY EQUATIONS IN ISITGR PATCH TO CAMB AND COSMOMC

## A. CMB implementation and the synchronous gauge

We first describe here the implementation of the MG formalism in the CMB software code CAMB [43] to calculate various CMB temperature anisotropy and polarization auto and cross spectra ( C TT glyph[lscript] , C TE glyph[lscript] , C EE glyph[lscript] , C BB glyph[lscript] ) as well as the three-dimensional matter power spectrum P δ ( k, z ). These are powerful probes to constrain both the growth history of structure in the Universe and the expansion history of the Universe. We describe the overall formalism and some key changes to the code in a self-contained way here but refer the reader for more details in the technical documentation of ISiTGR provided in the github repository https://github.com/mishakb/ISiTGR .

We recall that the package CAMB is written in the synchronous gauge where the perturbed FLRW metric is written as:

$$d s ^ { 2 } & = a ( \tau ) ^ { 2 } [ - d \tau ^ { 2 } + ( \gamma _ { i j } + h _ { i j } ) d x ^ { i } d x ^ { j } ], & \cdot \cdot$$

where h ij represents the metric perturbation in this gauge. The metric potentials are defined from the trace ( h ) and traceless ( h +6 ) part of the metric perturbation, following the notation of [60]. η h ij is explicitly expanded in terms of G , described further above, giving the following form for a single mode [60]

$$h _ { i j } = \frac { h } { 3 } \gamma _ { i j } G + ( h + 6 \eta ) ( k ^ { - 2 } G _ { | i j } + \frac { 1 } { 3 } \gamma _ { i j } G ).$$

Next, combining the perturbed Einstein's equations as discussed to get Eqs. (3) and (4), we have in the synchronous gauge [60]

$$\left ( k ^ { 2 } - 3 K \right ) \left ( \eta - \mathcal { H } \alpha \right ) \, = \, - 4 \pi G a ^ { 2 } \sum _ { i } \rho _ { i } \Delta _ { i },$$

$$k ^ { 2 } ( \dot { \alpha } + 2 \mathcal { H } \alpha - \eta ) \ = \ - 1 2 \pi G a ^ { 2 } \sum _ { i } \rho _ { i } ( 1 + w _ { i } ) \sigma _ { i },$$

where α = ( ˙ h +6˙) 2 η / k 2 . Next, using the gauge invariance of ∆ i and σ i and comparing equations (3) and (4) versus (36) and (37), we can see that the potentials in the two gauges are related by

$$\Phi = \eta - \mathcal { H } \alpha$$

and and it follows (39) that

$$\Psi = \dot { \alpha } + \mathcal { H } \alpha.$$

## B. CAMB variables and potential evolution

Additionally, CAMB defines two quantities that are used throughout the code when evolving the perturbations. These are

$$\sigma _ { \text{CAMB} } \equiv k \alpha = \frac { k ( \eta - \Phi ) } { \mathcal { H } },$$

$$\mathcal { Z } _ { \text{CAMB} } \equiv \dot { h } / 2 k = \sigma _ { \text{CAMB} } - 3 \frac { \dot { \eta } } { k }$$

$$\dot { \sigma } _ { C A M B } = k \Psi - \mathcal { H } \sigma _ { C A M B }.$$

CAMB evolves the metric potential η ( kη exactly) as well as the matter perturbations, δ i , the heat flux, q i , and the shear stress σ i (via Π i = 3 2 (1 + w i ) σ i ) for each matter species in the synchronous gauge according to the evolution equations as described in [60] and the documentation of CAMB . Additionally, the CAMB variables, σ CAMB , ˙ σ CAMB and Z CAMB are evaluated at each time step. Further, we give below the expressions for the key quantities ˙ η and ˙ α for each pair of MG parameters and taking into account massive neutrino, spatial curvature, and dynamical dark energy equation parameters.

In this new version of ISiTGR we account consistently for contributions from massive neutrinos and radiation as they enter the terms of the form 8 πGa 2 ∑ i ρ i Π (3 i w i +1) found in the expressions for the time evolution of the potentials, where the anisotropic shear Π i is defined as Π i = 3 2 (1 + w i ) σ i .

## 1. ˙ η for various MG parametrizations

Taking into account the specific modifications to the growth equations in each parametrization and not neglecting the anisotropic shear contribution leads to different expressions for ˙. η We note again that some relationships like (17) only hold when the anisotropic shear contribution is neglected. We provide the full expressions below.

First, taking the time derivative of (38) gives

$$\dot { \eta } = \dot { \mathcal { H } } \alpha + \mathcal { H } \dot { \alpha } + \dot { \Phi }. \\ \dots \dots \dots \dots \dots$$

Next, we need to analyze each parametrization separately.

For the ( µ, γ ) pair, we take the derivative of (14) and put it into (43) to obtain

$$\lambda \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colots \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \col: \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colot \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \col. \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \col.: \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \col. \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \color \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \col on \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \col:: \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colons \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colo \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colON \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colan \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colun \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colen \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colcon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colone \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \col { \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \$$

where we have defined β k ≡ (1 -3 K k 2 ) -1 with β k = 1 for the spatially flat case. We have also added the subscript label ( µ γ , ) to specify that we are deriving the expression for ˙ η in this particular parametrization. Using (39) and (13), we get

$$\dot { \alpha } = - \frac { \mu } { 2 k ^ { 2 } } \sum _ { i } \left [ \beta _ { k } \Delta _ { i } \hat { \rho } _ { i } + 2 \Pi _ { i } \hat { \rho } _ { i } \right ] - \mathcal { H } \alpha.$$

In order to work with the quantities used by CAMB , we rewrite (8) as

$$\dot { \Delta _ { i } } & = 3 ( 1 + w _ { i } )$$

where f 1 ≡ 1 + 3( H -H 2 ˙ ) k 2 and, for now, we use the superscript ( N ) to denote that in this relation the heat flux q i is still given in the Newtonian gauge; this will be dealt with below. Note that we now use ∆ ˙ i for the anisotropic stress used in CAMB .

Next, putting equations (45) and (46) into (44), and using the continuity equation

$$\begin{array} {$$

as well as Eqs. (38) and (39), and the relationship

$$q _ { i } ^ { ( N ) } = q _ { i } ^ { ( S )$$

for q i between the two gauges, we obtain the final expression for ˙ η ( µ,γ ) as

$$\text{between the two gauges, we obtain the final expression for } \dot { \eta }$$

where we defined

$$f _ { \mu, \gamma } \equiv k ^ { 2 } + \frac {$$

Following similar steps, using the corresponding equations for the ( µ, Σ) parametrization, we derive

$$\text{rnowing similar steps, using one corresponding equations for use } \mu, \, \exp \text{parameterization, we receive} \\ \hat { \eta } _ { ( \mu, \Sigma ) } & = \frac { 1 } { 2 f _ { \mu, \Sigma } } \left \{ k \beta _ { k } ( 2 \Sigma - \mu ) f _ { 1 } \sum _ { i } q _ { i } ^ { ( s ) } \hat { \rho } _ { i } + \beta _ { k } [ ( \mu - 2 \dot { \Sigma } ) + 2 \mathcal { H } ( \Sigma - \mu ) ] \sum _ { i } \Delta _ { i } \hat { \rho } _ { i } \\ & + 2 ( \mu - \Sigma ) \sum _ { i } \mathbb { I } _ { i } \hat { \rho } _ { i } + 2 [ ( \mu - \dot { \Sigma } ) + \mathcal { H } \beta _ { k } ( 2 \Sigma - \mu ) - \mathcal { H } \mu ] \sum _ { i } \Pi _ { i } \hat { \rho } _ { i } \\ & + 2 \mathcal { H } ( \Sigma - \mu ) \sum _ { i } \Pi _ { i } ( 1 + 3 w _ { i } ) \hat { \rho } _ { i } + k ^ { 2 } \alpha \left [ \beta _ { k } ( 2 \Sigma - \mu ) \sum _ { i } ( 1 + w _ { i } ) \hat { \rho } _ { i } - 2 ( \mathcal { H } ^ { 2 } - \mathcal { H } ) \right ] \right \}, \\ \text{where}$$

where

$$f _ { \mu, \Sigma } \equiv k ^ { 2 } + \frac { 3 } { 2 } \beta _ { k } ( 2 \Sigma - \mu ) \sum _ { i } \hat { \rho } _ { i } ( 1 + w _ { i } ).$$

Again, it is important to recall that we cannot go from (51) to (49) just by using the usual relation Σ = µ 2 (1 + γ ), because this is only valid when contributions to the anisotropic stress are zero for all species. These are more general equations including non-zero shear terms, and they are valid for dynamical dark energy and for flat or curved spaces. Finally, we derive the expression for the ( Q a, k , D a, k ( ) ( )) parametrization as

$$\dot { \eta } _ { ( Q, D ) } & = - \frac { 1 } { 2 f _ { Q } } \left \{ 2 k ^ { 2 } \alpha ( \mathcal { H } ^ { 2 } - \dot { \mathcal { H } } ) + [ 2 \mathcal { H } ( D - Q ) + \dot { Q } ] \beta _ { k } \sum _ { i } \Delta _ { i } \hat { \rho } _ { i } \\ & - k ^ { 2 } \alpha \sum _ { i } \beta _ { k } Q ( 1 + w _ { i } ) \hat { \rho } _ { i } - k \beta _ { k } Q f _ { 1 } \sum _ { i } q _ { i } ^ { ( S ) } \hat { \rho } _ { i } - 2 \mathcal { H } Q ( \beta _ { k } - 1 ) \sum _ { i } \hat { \rho } _ { i } \Pi _ { i } \right \}, \\.$$

where

$$f _ { Q } \equiv k ^ { 2 } + \frac { 3 } { 2 } \beta _ { k } Q \sum _ { i } ( 1 + w _ { i } ) \hat { \rho } _ { i }.$$

In this case, this expression is also valid for the pair ( Q a, k , R a, k ( ) ( )) upon using that R a.k ( ) = 2 D a,k ( ) Q a,k ( ) -1, which is a definition rather than a relationship.

Finally, we made some comparisons of the angular and matter power spectra ( C TT glyph[lscript] and P k ( )) between ISiTGR and MGCAMB [48] and found them to be in an agreement to 10 -4 or better for some of the scale-independent parametrizations listed in section IV.

## VI. MODIFIED POWER SPECTRA AND CORRELATION FUNCTIONS

## A. 3 × 2 point power spectra and correlation functions

The current version of ISiTGR consistently implements the often-used 3 × 2 point statistics (lensing - lensing, galaxy - galaxy lensing and galaxy - galaxy) by modifying the respective transfer functions via modification to the Weyl potential (for the first two) satisfying, for example, the MG equation (16). The angular power spectra are then calculated as

$$P _ { \kappa \kappa } ^ { i j } ( \ell ) = \int _ { 0 } ^ { \chi _ { H } } d \chi \frac { q _ { \kappa } ^ { i } ( \chi ) q _ { \kappa } ^ { j } ( \chi ) } { \chi ^ { 2 } } P _ { \Phi _ { w }, \Phi _ { w } } \left ( \frac { \ell + 1 / 2 } { f _ { K } ( \chi ) }, \chi \right ),$$

for lensing - lensing and

$$P _ { \delta \kappa } ^ { i j } ( \ell ) = \int _ { 0 } ^ { \chi _ { H } } d \chi \frac { q _ { \delta } ^ { i } \left ( \frac { \ell + 1 / 2 } { f _ { K } ( \chi ) }, \chi \right ) q _ { \kappa } ^ { j } ( \chi ) } { \chi ^ { 2 } } P _ { \delta, \phi w } \left ( \frac { \ell + 1 / 2 } { f _ { K } ( \chi ) }, \chi \right ),$$

for the galaxy-galaxy lensing, where P Φ W , Φ W is the Weyl potential power spectrum, and f K ( χ ) here is the comoving angular diameter distance. The lensing efficiency function and the radial weight function are given by

$$q _ { \kappa } ^ { i } ( \chi ) = \chi \int _ { \chi } ^ { \chi _ { H } } d \chi ^ { \prime } n ^ { i } ( \chi ) \frac { f _ { K } ( \chi ^ { \prime } - \chi ) } { f _ { K } ( \chi ) }$$

and

$$q _ { \delta } ^ { i } \left ( \frac { \ell + 1 / 2 } { f _ { K } ( \chi ) }, \chi \right ) = b ^ { i } ( k, z ( \chi ) ) n ^ { i } ( \chi ),$$

respectively. For the radial weight function, b i ( k, z ( χ )) corresponds to the galaxy bias. It is worth mentioning that we are using the two Newtonian Gauge potentials directly to compute the power spectrum, so we do not use the equality k 2 Φ W = -3 2 (1 + z ) ( H 0 c ) Ω m δ , which is only valid when Φ = Ψ.

The corresponding lensing-lensing (or shear-shear) 2-point correlation functions are calculated as

$$\xi _ { \pm } ^ { i j } ( \theta ) = \int \frac { d \ell \, \ell } { 2 \pi } P _ { \kappa \kappa } ^ { i j } ( \ell ) J _ { 2 \pm 2 } ( \ell \theta ),$$

where J n is the n th order Bessel function of the first kind. Similarly, the galaxy - galaxy lensing 2-point correlation function is given by

$$\gamma _ { t } ^ { i j } ( \theta ) = \int \frac { d \ell \ \ell } { 2 \pi } P _ { \delta \kappa } ^ { i j } ( \ell ) J _ { 2 } ( \ell \theta )$$

and the clustering 2D correlation function is given by

$$w ^ { i } ( \theta ) = \int \frac { d \ell \ \ell } { 2 \pi } P _ { \delta \delta } ^ { i i } ( \ell ) J _ { 0 } ( \ell \theta ),$$

where P ii δδ ( glyph[lscript] ) is the matter power spectrum. We refer the reader to, for example, [78, 79] for details on the 3 × 2 point formalism.

## B. Integrated Sachs-Wolfe (ISW) effect

In order to propagate the MG equation changes to the ISW effect, we follow CAMB 's implementation to calculate the derivatives of MG potentials. This is done for Φ by recalling Eq. ˙ (43) and using ˙ η and α = σ CAMB k as

$$\dot { \Phi } = \dot { \eta } - \dot { \mathcal { H } } \alpha - \mathcal { H } \dot { \alpha }. \\.. \ \cdot \. \ \cdot \. \ \cdot \. \ \cdot \. \ \cdot \. \ \cdot \. \ \cdot \. \$$

Then, one obtains Ψ by directly taking the derivative of (13) which depends only on ˙ µ , so it can be used for both the ( µ, γ ) and ( µ, Σ) parametrizations. This gives

$$\dot { \Psi } = - \frac { \dot { \mu } } { 2 k ^ { 2 } } \sum _ { i } \left [ \beta _ { k } \Delta _ { i } \hat { \rho } _ { i } + 2 \Pi _ { i } \hat { \rho } _ { i } \right ] + \frac { \mu } { 2 k ^ { 2 } } \sum _ { i } \left [ ( \beta _ { k } \Delta _ { i } + 2 \Pi _ { i } ) \dot { \hat { \rho } } _ { i } + ( \beta _ { k } \dot { \Delta } _ { i } + 2 \dot { \Pi } _ { i } ) \hat { \rho } _ { i } \right ].$$

Next, we can substitute (46), (47) and (48) into the above equation to obtain

$$\cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdot \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \colon \\ \dot { \Psi } } & = - \frac { \dot { \mu } } { 2 k ^ { 2 } } \sum _ { i } [ \beta _ { k } \Delta _ { i } \hat { \rho } _ { i } + 2 \Pi _ { i } \hat { \rho } _ { i } ] + \frac { \mu } { 2 k ^ { 2 } } \sum _ { i } \left \{ \mathcal { H } \beta _ { k } \Delta _ { i } \hat { \rho } _ { i } - 2 \dot { \Pi } _ { i } \hat { \rho } _ { i } + k \beta _ { k } f _ { 1 } q _ { i } ^ { ( S ) } \hat { \rho } _ { i } } \dot { \rho } _ { i } \\ & \quad + \beta _ { k } ( 1 + w _ { i } ) \hat { \rho } _ { i } \dot { [ k } ^ { 2 } \alpha f _ { 1 } - 3 ( \dot { \Phi } + \mathcal { H } \Psi ) ] + 2 \mathcal { H } \Pi _ { i } \hat { \rho } _ { i } ( 1 + 3 w _ { i } ) + 2 \beta _ { k } \mathcal { H } \hat { \rho } _ { i } \Pi _ { i } \hat { \rho } _ { i } ( 1 + 3 w _ { i } ) + 2 \beta _ { k } \mathcal { H } \hat { \rho } _ { i } \Pi _ { i } \right \}. \\ \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \ddots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \\ \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots$$

Thus, the MG changes to the ISW are propagated through (62) and (64).

Finally, for the ISW effect for the ( Q,D ) and ( Q,R ) parametrizations, we derive and use another expression for Ψ. ˙ Following similar steps we can obtain

$$\dot { \Psi } & = - \dot { \Phi } - \frac { 1 } { k ^ { 2 } } \sum _ { i } \left \{ \beta _ { k } \dot { D } \Delta _ { i } \hat { \rho } _ { i } + \dot { Q } \Pi _ { i } \hat { \rho } _ { i } - \beta _ { k } D \mathcal { H } \Delta _ { i } \hat { \rho } _ { i } - Q \mathcal { H } \Pi _ { i } ( 3 w _ { i } + 1 ) \hat { \rho } _ { i } + Q \dot { \Pi } _ { i } \hat { \rho } _ { i } \\ & - 2 \mathcal { H }. D \mathcal { H } \Pi _ { i } \hat { \rho } _ { i } - k \mathcal { H }. D f. a ^ { ( S ) } \hat { \rho } _ { i } + \mathcal { H }. ( 1 + n _ { i }. \lambda \hat { \rho } _ { i }! 3 D ( \dot { \Phi } + \mathcal { H } \Psi ) - k ^ { 2 } \alpha f. D l }$$

-

2

β D

(

S

)

k

H

Π ˆ

i

ρ

i

-

kβ Df q

k

1

2

i

ˆ +

ρ

i

β

k

(1 +

w

i

) ˆ [3

ρ

i

D

( ˙

Φ+

H

Ψ)

-

k αf D

1

]

.

Here, it is important to mention that firstly, the code computes Φ and afterwards, it computes Ψ in terms of Φ. ˙ ˙ ˙

## VII. APPLICATIONS AND RESULTS FROM CURRENT AVAILABLE DATA SETS

In this section we apply ISiTGR to current cosmological data sets using different MG parametrizations. As a start, we reproduce and compare to some of the results from the Planck-2018 [67], Planck-2015 [66] and DES-2018 [73] papers about extended models. We then derive new results involving the different features of ISiTGR such as spatial curvature, dynamical dark energy, and massive neutrinos along with MG parameters. Moreover, we derive constraints and correlations for the new binning methods implemented in ISiTGR for MG parameters.

For all results, in addition to MG parameters, we also vary the six core cosmological parameters: Ω b h 2 and Ω c h 2 , the baryon and cold dark matter physical density parameters, respectively; θ , the ratio of the sound horizon to the angular diameter distance of the surface of last scattering; τ , the reionization optical depth; n s , the spectral index; and ln(10 10 A s ), the amplitude of the primordial power spectrum.

## A. Data sets

We combine several current data sets including measurements from the Cosmic Microwave Background (CMB) Planck mission. Specifically, the Planck-2018 TT, TE and EE spectra at glyph[lscript] ≥ 30 likelihood, the lowglyph[lscript] temperature Commander likelihood, the lowglyph[lscript] SimAll EE likelihood and CMB lensing measurements, used in [67]. We label this combination as TTTEEE + lowE + CMBlens (2018). Moreover, we use the Planck-2015 likelihood presented in [80], where the Planck TT likelihood for multipoles 30 ≤ glyph[lscript] ≤ 2508 and the joint TT, EE, BB, and TE likelihood for 2 ≤ glyph[lscript] ≤ 30 were used. We refer to this likelihood combination as TT + lowP (2015) . 2 We also use in some cases the Planck-2015 CMB lensing measurements for temperature from [81]. We also use the likelihood included into CosmoMC for clustering and lensing data from the Dark Energy Survey (DES) Year 1 [79], but with the non-linear data points removed since these are not yet accurately modeled for MG theories, see e.g. [73] for a discussion. We also add measurements coming from Type Ia supernovae. Namely, the Pantheon sample data presented in [82], which combines 279 SNe Ia (0 03 . &lt; z &lt; 0 68) with useful distance estimates of SNe Ia from SDSS, SNLS, various low-z and . HST samples, giving a total of 1048 SNe Ia ranging from (0 01 . &lt; z &lt; 2 3). . In this work, we refer to these data sets as DES and Pantheon, respectively. Finally, we also consider measurements from Baryon Acoustic Oscillations (BAO) and Redshift Space Distortions (BAO/RSD). Specifically, we use BAO and BAO/RSD measurements coming from the BOSS Data Release 12 [83]. Moreover, we use another two BAO data sets from the 6DF Galaxy Survey [84] and the SDSS Data Release 7 Main Galaxy Sample [85].

## B. Constraints and correlations

In the case of the ( µ, γ ) parametrization, we reproduce some results for the constraints and scale dependence as shown in Planck's 2015 dark energy and modified gravity paper [66], while for the ( µ, Σ) parametrization we reproduce some of the results of the DES 2018 constraints on extended cosmological models. Then, we add the new features of ISiTGR in combination with the MG parameters. Finally, we obtain constraints and correlations for the binning methods for ( µ, γ ) and ( µ, Σ).

## 1. Results for ( µ, γ ) parametrization

We evolve these parameters in time and scale according to Eqs. (22) and (23) as in, for example, [66]. Our results are shown in Fig. 1 and Table II, and are in good agreement for the time dependence with [67] and [66] but we find some differences when adding scale dependence, see Fig. 10. Specifically, the contour plots agree for the large scales with k = 10 -7 Mpc -1 but not for small scales with large k = 10 Mpc 2 -1 as shown in their figure 18 [66]. With a closer look, one can see that for such large k the MG parameters should reduce to µ → 1 + E 11 Ω DE and γ → 1 + E 22 Ω DE , so they becomes similar to the scale independent case, as we find in our Fig. 10. For large scales (small k ) the MG parameters should reduce to µ → 1 + c 1 E 11 Ω DE and γ → 1 + c 2 E 22 Ω DE . Therefore, after marginalizing over c 1 and c 2 constraints become weaker and our results are in agreement for this regime. This difference and further study of scale-dependence will be explored further in a separate study.

2 We use the public Planck likelihood package files plik dx11dr2 HM v18 TT.clik and lowl SMW 70 dx11d 2014 10 03 v5c Ap.clik, unless otherwise specified.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00290v2-figure-25ceed6d-2ac4-4302-838c-f587efc9ba1e.jpg)

-

TTTEEE + lowE + CMBlens (2018)

.

1 5

1

-

1

= 0)

(

µ z

0 5

.

0

-

0 5

.

-

1

+ BAO + SNe

+ BAO/RSD + DES

0

1

2

3

1

0 00

.

TT + lowP (2015)

+ BSH

+ BAO/RSD

.

0 25

0 50

.

0 75

.

1 00

Σ(

z

= 0)

-

1

TTTEEE + lowE + CMBlens (2018)

+ BAO + SNe

+ BAO/RSD + DES

.

0 00

0 25

.

0 50

.

0 75

.

1 00

.

.

2

.

1 5

1

-

1

= 0)

(

µ z

0 5

.

0

-

0 5

.

-

1

2

.

1 5

1

-

1

= 0)

(

µ z

0 5

.

0

-

0 5

.

-

-

1

η z

(

= 0)

-

1

Σ(

z

= 0)

-

1

FIG. 1: Left: 68% and 95% confidence contour plots for µ 0 and η 0 in the ( µ, γ ) parametrization with DE-time evolution and no scale dependence. Figures at the top show constraints with Planck-2015 data while figures at the bottom show constraints using Planck-2018 data. This is in good agreement with results from Planck-2015 [66] (see left panel of Fig. 14 there) and Planck-2018 [67]. Right: Contour plots for µ 0 and Σ 0 in the ( µ, Σ) parametrization with similar evolution. This is again in good agreement with Fig. 15 of Planck-2015 [66] and Planck-2018 [67] (see our Table II).

| Parameter      | Planck-2015                | Planck-2015 + BSH          | Planck-2015 + BAO/RSD      | Planck-2018                   | Planck-2018 + BAO + SNe    | Planck-2018 + BAO/RSD + DES   |
|----------------|----------------------------|----------------------------|----------------------------|-------------------------------|----------------------------|-------------------------------|
| E 11           | 0 . 07 +0 . 31 - 0 . 69    | 0 . 03 +0 . 30 - 0 . 63    | - 0 . 25 +0 . 19 - 0 . 32  | 0 . 13 +0 . 46 - 0 . 66       | 0 . 12 +0 . 42 - 0 . 71    | - 0 . 10 +0 . 29 - 0 . 53     |
| E 22           | 1 . 1 ± 1 . 4              | 1 . 1 ± 1 . 3              | 1 . 80 ± 0 . 88            | 0 . 30 +0 . 77 - 1 . 6        | 0 . 4 ± 1 . 1              | 0 . 5 ± 1 . 0                 |
| µ ( z = 0) - 1 | 0 . 05 +0 . 22 - 0 . 48    | 0 . 02 +0 . 21 - 0 . 44    | - 0 . 18 +0 . 14 - 0 . 22  | 0 . 09 +0 . 31 - 0 . 46       | 0 . 09 +0 . 29 - 0 . 50    | - 0 . 07 +0 . 20 - 0 . 36     |
| η ( z = 0) - 1 | 0 . 77 +0 . 94 - 1 . 1     | 0 . 78 ± 0 . 88            | 1 . 25 ± 0 . 62            | 0 . 20 +0 . 50 - 1 . 1        | 0 . 25 ± 0 . 80            | 0 . 35 ± 0 . 70               |
| Σ( z = 0) - 1  | 0 . 29 ± 0 . 15            | 0 . 28 ± 0 . 14            | 0 . 29 ± 0 . 13            | 0 . 07 +0 . 12 - 0 . 10       | 0 . 082 +0 . 12 - 0 . 096  | - 0 . 007 +0 . 076 - 0 . 060  |
| Ω b h 2        | 0 . 02253 ± 0 . 00027      | 0 . 02253 ± 0 . 00021      | 0 . 02246 ± 0 . 00022      | 0 . 02244 ± 0 . 00015         | 0 . 02249 ± 0 . 00014      | 0 . 02247 ± 0 . 00014         |
| Ω c h 2        | 0 . 1173 ± 0 . 0025        | 0 . 1176 ± 0 . 0012        | 0 . 1179 ± 0 . 0014        | 0 . 1190 ± 0 . 0015           | 0 . 1185 ± 0 . 0010        | 0 . 1187 ± 0 . 0010           |
| 100 θ          | 1 . 04126 ± 0 . 00051      | 1 . 04121 ± 0 . 00044      | 1 . 04115 ± 0 . 00044      | 1 . 04104 ± 0 . 00031         | 1 . 04110 ± 0 . 00029      | 1 . 04103 ± 0 . 00029         |
| ln(10 10 A s ) | 3 . 057 ± 0 . 040          | 3 . 058 ± 0 . 041          | 3 . 052 ± 0 . 038          | 3 . 031 +0 . 018 - 0 . 016    | 3 . 032 ± 0 . 017          | 3 . 042 +0 . 015 - 0 . 017    |
| n s            | 0 . 9719 ± 0 . 0072        | 0 . 9711 ± 0 . 0046        | 0 . 9698 ± 0 . 0047        | 0 . 9665 ± 0 . 0047           | 0 . 9679 ± 0 . 0040        | 0 . 9669 ± 0 . 0039           |
| τ              | 0 . 064 ± 0 . 020          | 0 . 064 ± 0 . 021          | 0 . 061 ± 0 . 019          | 0 . 0494 +0 . 0083 - 0 . 0072 | 0 . 0504 ± 0 . 0079        | 0 . 0548 +0 . 0073 - 0 . 0082 |
| H 0            | 68 . 6 ± 1 . 2             | 68 . 44 ± 0 . 57           | 68 . 26 ± 0 . 66           | 67 . 80 ± 0 . 66              | 68 . 04 ± 0 . 46           | 67 . 93 ± 0 . 46              |
| σ 8            | 0 . 816 +0 . 033 - 0 . 052 | 0 . 814 +0 . 029 - 0 . 048 | 0 . 792 +0 . 020 - 0 . 024 | 0 . 813 +0 . 031 - 0 . 046    | 0 . 811 +0 . 028 - 0 . 050 | 0 . 800 +0 . 023 - 0 . 036    |

TABLE II: Marginalized mean values and 1σ errors for cosmological parameters from various data set combinations. We refer to the combination TTTEEE + lowE + CMBlens (2018) as Planck-2018, while we label TT + lowP (2015) as Planck-2015.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00290v2-figure-fe488014-1886-42f6-9ca9-65e27df030a4.jpg)

-

FIG. 2: Results on scale- dependent evolution for µ and γ (called η in ISiTGR ). We find that for large scales the constraints become weaker than in the scale-independent case, as expected. Here, we plot the cases in which k = 10 -7 Mpc -1 and k = 10 -10 Mpc -1 , as well as the scale-independent case , all for z = 0. As we mentioned in section IV D, for large scales µ → 1+ c 1 E 11 Ω DE and γ → 1+ c 2 E 22 Ω DE . Therefore, after marginalizing over c 1 and c 2 constraints become weaker. However, for small scales we have that µ → 1+ E 11 Ω DE and γ → 1+ E 22 Ω DE , then practically this reduce to the scale independent case as we found. The plot in this case of k = 10 2 Mpc -1 is not in agreement with that of figure 18 of Planck-2015 [66]. See text for discussion.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00290v2-figure-ca45faef-1b67-480c-815e-b84608a7e933.jpg)

Ω

-

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00290v2-figure-9aa65079-1d38-4156-9b10-286d1c129d8a.jpg)

-

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00290v2-figure-3299558b-cf5b-4c33-b7f3-a02d8527bb27.jpg)

-

FIG. 3: 68% and 95% confidence contours for the parameters µ 0 -1, η 0 -1 and Σ 0 -1 with the addition of the extra feature of spatial curvature implemented in ISiTGR .

For the ( µ, γ ) parametrization without scale dependence, we set λ = 0 and vary the parameters E 11 and E 22 from which µ and γ are constructed. We do not include the weak lensing plots here since the weak lensing likelihood has been replaced by the DES likelihood implemented in the current version of CosmoMC . Also, Planck-2018 reported a bug in the Weyl potential which may have affected the previous weak lensing results in [66]. It is worth mentioning that, for this specific case, we use the combination of BAO from [86], Supernovae Type Ia data from the Joint Light-Curve Analysis [87], and the Hubble constant measurements obtained from the Hubble Space Telescope in [88]. We refer to this combination as BSH, as in previous analysis [66].

When adding the scale dependence to MG parameters, we not only vary the parameters E 11 and E 22 for the time dependence, but also c 1 , c 2 and λ for the scale dependence. In Fig. 10 we show the corresponding results. As we can see, for large scales the constraints become weaker because of the marginalization over the parameters c 1 and c 2 .

µ

FIG. 4: 68% and 95% confidence contour plots for µ 0 and Σ 0 in the ( µ, Σ) parametrization with DE-time evolution and no scale dependence. This is in agreement with results from DES 2018 [73] (Figure 3 there). Some small differences may be due to different priors, the likelihood/sampling methods or the process related to removing the non-linear data. We can see that adding the DES data to the EXT data set removes any tension with GR.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00290v2-figure-4985c066-c0cd-4bd3-95c5-282cd372cab6.jpg)

TABLE III: Marginalized mean values of and 1σ errors on MG and cosmological parameters using DES+EXT data sets. The results are well within 1σ agreement with [73].

| Parameter      | EXT                           | EXT + DES                     |
|----------------|-------------------------------|-------------------------------|
| µ 0            | 0 . 14 ± 0 . 27               | 0 . 03 +0 . 30 - 0 . 63       |
| Σ 0            | 0 . 29 +0 . 15 - 0 . 13       | 0 . 053 +0 . 059 - 0 . 066    |
| Σ m ν          | 0 . 181 +0 . 088 - 0 . 12     | 0 . 146 +0 . 062 - 0 . 11     |
| Ω b h 2        | 0 . 02244 ± 0 . 00021         | 0 . 02242 ± 0 . 00021         |
| Ω c h 2        | 0 . 1163 +0 . 0018 - 0 . 0016 | 0 . 1165 +0 . 0015 - 0 . 0013 |
| 100 θ          | 1 . 04127 ± 0 . 00042         | 1 . 04123 ± 0 . 00042         |
| ln(10 10 A s ) | 3 . 041 ± 0 . 042             | 3 . 082 +0 . 033 - 0 . 038    |
| n s            | 0 . 9723 +0 . 0050 - 0 . 0057 | 0 . 9718 ± 0 . 0048           |
| τ              | 0 . 058 +0 . 021 - 0 . 023    | 0 . 078 +0 . 017 - 0 . 020    |

Finally, we extend the results on this parametrization by adding spatial curvature to the ( µ, γ ) case. We show the corresponding constraints of the MG parameters in terms of Ω k in Fig. 3. Since we are not considering contributions from massive neutrinos in this case, we take into account the derive parameter Σ -1 at z = 0. When the stress shear is negligible, the relationship Σ = µ 2 (1 + η ) is satisfied.

wa

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00290v2-figure-bc5b5b20-a669-4f5e-b33d-bdc0eb2d9763.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00290v2-figure-aa5f105d-f8f7-4428-8f13-4ceeeab3915e.jpg)

wa

FIG. 5: Left: 68% and 95% confidence contour plots for w 0 and w a for various data sets. We label the combination TT + lowP (2015) + CMBlens + BAO + BAO/RSD + Pantheon as EXT (for external data set) as done in DES 2018 so a direct comparison can be made. Right: Confidence contour plots for the pivot equation of state that uses w p and w a . We can see that in the right-hand side plot, the parameters w p and w a are decorrelated. We used a pivot redshift of z p = 0 21 for EXT . data, and the pivot redshift given in DES 2018, z p = 0 20, for EXT+DES data. .

| Parameter   | EXT                        | EXT + DES                    | Planck-2018 + BAO +SNe    | Planck-2018 + BAO/RSD + DES   |
|-------------|----------------------------|------------------------------|---------------------------|-------------------------------|
| w 0         | - 0 . 965 +0 . 07 - 0 . 08 | - 0 . 972 +0 . 068 - 0 . 081 | - 0 . 956 ± 0 . 082       | - 0 . 74 ± 0 . 20             |
| w a         | - 0 . 24 +0 . 45 - 0 . 27  | - 0 . 27 +0 . 44 - 0 . 29    | - 0 . 32 +0 . 34 - 0 . 28 | - 0 . 83 ± 0 . 59             |
| w p         | - 1 . 005 ± 0 . 041        | - 0 . 999 ± 0 . 043          | -                         | -                             |
| H 0         | 67 . 84 ± 0 . 75           | 67 . 82 ± 0 . 78             | 68 . 28 ± 0 . 83          | 66 . 2 +1 . 7 - 1 . 9         |
| σ 8         | 0 . 812 +0 . 016 - 0 . 013 | 0 . 801 +0 . 016 - 0 . 014   | 0 . 825 ± 0 . 011         | 0 . 801 +0 . 015 - 0 . 017    |
| S 8         | 0 . 823 ± 0 . 014          | 0 . 811 ± 0 . 014            | 0 . 835 ± 0 . 011         | 0 . 834 ± 0 . 013             |

TABLE IV: 68% confidence limits for the DE parameters and cosmological parameters. The constraints are well within the 1σ level agreement with the constraints presented from [67, 73].

## 2. Results for ( µ, Σ) parametrization

We use here Eqs. (24) and (25) to evolve µ and Σ, but set λ = 0 to eliminate the scale dependence. Furthermore, we removed the non-linear data from DES measurements as suggested in [66] and [73], until we get ∆ χ 2 &lt; 1 where ∆ χ 2 ≡ ( d NL -d L ) T C -1 ( d NL -d L ) and d L represents the vector non-linear data predictions and d L represents the vector linear data predictions.

We let the neutrino mass vary freely and just for this case we consider the Planck likelihood file plik lite v18 TT.CLIK instead of plik dx11dr2 HM v18 TT.clik. This makes it possible to compare some of our results to those of DES 2018 analysis of extended models [73]. We label the combination TT + lowP (2015) + CMBlens + BAO + BAO/RSD + Pantheon as EXT (for external data set). The priors on µ 0 and Σ 0 were set as in [73]. In Fig. 4 and Table III, we show our results from ISiTGR and find them in overall agreement with the DES 2018 analysis using CosmoSIS [89] and MGCAMB [90]. Some small differences between the constraints may be due to other parameter priors (for example the neutrino mass), the likelihood/sampling method differences or the process related with removing the non-linear data.

FIG. 6: 68% and 95% confidence contours for the constraints obtained using the MG parameters µ 0 and Σ 0 for the ( µ, Σ) parametrization, combined with the dark energy equation of state ( w , w 0 a ), and allowing the neutrino mass to vary. GR is found to be consistent with current data sets.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00290v2-figure-6f50a507-fd5b-4c93-95c1-31c39c89c62d.jpg)

## 3. Results for different dark energy parametrizations

In addition to the usual MG parameter approach, we also show our results for the new dark energy parametrizations implemented in ISiTGR . We present our results for dark energy parameter constraints in Table IV. The constraints shown for EXT and EXT + DES can be compared with the results obtained in [73], while the constraints using Planck-2018 data should be compared with results in [67]. The contraints are found within 1σ agreement with these works. See Fig. 5 for the corresponding contour plots.

Furthermore, ISiTGR is able to combine these new dark energy parametrizations with MG parameters and massive neutrinos. We show in Fig. 6 the results obtained when these features are combine.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00290v2-figure-76a93276-a173-45a3-8d7f-b549b033ac01.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00290v2-figure-49fc8a0d-3f73-4a05-8e5c-582c895a5b63.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00290v2-figure-fa9169c7-4dff-4a12-8ff9-9d25167a3c8e.jpg)

µ

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00290v2-figure-b91a593e-9e3a-4b4b-8f15-2d6212b34a8f.jpg)

FIG. 7: 68% and 95% confidence contours for both the traditional binning method and the hybrid binning method. To obtain the constraints we use the combination TT + lowP (2015) + Pantheon + BAO + BAO/RSD in a joint analysis. The z-bins are 0 &lt; z ≤ 1 and 1 &lt; z ≤ 2, with GR assumed for z &gt; 2. The k-bins are k ≤ 0 01 and . k &gt; 0 01. . We plot the binning parameters for the ( µ, γ ) parametrization (in ISiTGR , γ is referred to as η ).

FIG. 8: 68% and 95% confidence contours for both the traditional binning method and the hybrid binning method, using the binning parameters for the ( µ, Σ) parametrization. We use the combination TT + lowP (2015)+Pantheon+BAO+BAO/RSD in a joint analysis. The z-bins are 0 &lt; z ≤ 1 and 1 &lt; z ≤ 2 with GR assumed for z &gt; 2. Also, the k-bins are k ≤ 0 01 and . k &gt; 0 01.) .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00290v2-figure-59d5e2e4-7222-4c76-9cab-2128abac2d2b.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00290v2-figure-d4df5962-8dc3-4f36-8584-5a4b49603598.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00290v2-figure-866b36ad-e5ad-4576-8591-cc2cacf2b899.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00290v2-figure-9060ee3f-faa7-473b-8065-8ab940f402a9.jpg)

TABLE V: 95% confidence limits for the MG parameters in both the traditional binning method and the hybrid binning method. We present the results for the new parametrizations included in this new version of ISiTGR , the ( µ, γ ) and ( µ, Σ) parametrizations. We found some small tension between the MG parameters and GR when looking at the 68% confident limits. However, this disagreement disappears when considering the 95% confidence limits.

| Constraints for the MG parameters using the binning methods in the ( µ, γ ) parametrization   | Constraints for the MG parameters using the binning methods in the ( µ, γ ) parametrization   | Constraints for the MG parameters using the binning methods in the ( µ, γ ) parametrization   | Constraints for the MG parameters using the binning methods in the ( µ, γ ) parametrization   | Constraints for the MG parameters using the binning methods in the ( µ, γ ) parametrization   | Constraints for the MG parameters using the binning methods in the ( µ, γ ) parametrization   | Constraints for the MG parameters using the binning methods in the ( µ, γ ) parametrization   | Constraints for the MG parameters using the binning methods in the ( µ, γ ) parametrization   | Constraints for the MG parameters using the binning methods in the ( µ, γ ) parametrization   |
|-----------------------------------------------------------------------------------------------|-----------------------------------------------------------------------------------------------|-----------------------------------------------------------------------------------------------|-----------------------------------------------------------------------------------------------|-----------------------------------------------------------------------------------------------|-----------------------------------------------------------------------------------------------|-----------------------------------------------------------------------------------------------|-----------------------------------------------------------------------------------------------|-----------------------------------------------------------------------------------------------|
| Traditional binning Hybrid binning                                                            | µ 1 1 . 06 +0 . 52 - 0 . 45 1 . 02 +0 . 55 - 0 . 51                                           | µ 2 1 . 08 +0 . 37 - 0 . 35 1 . 06 +0 . 37 - 0 . 37                                           | µ 3 1 . 04 +0 . 42 - 0 . 38 1 . 05 +0 . 52 - 0 . 47                                           | µ 4 0 . 91 +0 . 28 - 0 . 25 0 . 92 +0 . 36 - 0 . 36                                           | η 1 1 . 07 +1 . 0 - 0 . 90 1 . 2 +1 . 1 - 1 . 1                                               | η 2 1 . 03 +0 . 72 - 0 . 63 1 . 03 +0 . 81 - 0 . 70                                           | η 3 1 . 04 +0 . 87 - 0 . 73 < 1 . 96                                                          | η 4 1 . 30 +0 . 69 - 0 . 65 1 . 30 +0 . 92 - 0 . 87                                           |
| Constraints for the MG parameters using the binning methods in the ( µ, Σ) parametrization    | Constraints for the MG parameters using the binning methods in the ( µ, Σ) parametrization    | Constraints for the MG parameters using the binning methods in the ( µ, Σ) parametrization    | Constraints for the MG parameters using the binning methods in the ( µ, Σ) parametrization    | Constraints for the MG parameters using the binning methods in the ( µ, Σ) parametrization    | Constraints for the MG parameters using the binning methods in the ( µ, Σ) parametrization    | Constraints for the MG parameters using the binning methods in the ( µ, Σ) parametrization    | Constraints for the MG parameters using the binning methods in the ( µ, Σ) parametrization    | Constraints for the MG parameters using the binning methods in the ( µ, Σ) parametrization    |
| Traditional binning Hybrid binning                                                            | µ 1 1 . 11 +0 . 84 - 0 . 93 1 . 05 +0 . 82 - 1 . 0                                            | µ 2 1 . 63 +0 . 77 - 0 . 83 1 . 58 +0 . 85 - 0 . 95                                           | µ 3 1 . 01 +0 . 75 - 0 . 83 0 . 98 +0 . 76 - 0 . 89                                           | µ 4 < 1 . 07 < 1 . 21                                                                         | Σ 1 1 . 043 +0 . 096 - 0 . 099 1 . 07 +0 . 19 - 0 . 21                                        | Σ 2 1 . 11 +0 . 15 - 0 . 16 1 . 08 +0 . 18 - 0 . 18                                           | Σ 3 1 . 017 +0 . 056 - 0 . 055 1 . 01 +0 . 11 - 0 . 10                                        | Σ 4 1 . 046 +0 . 094 - 0 . 10 1 . 05 +0 . 12 - 0 . 13                                         |

## 4. Results for binning method for time and scale dependencies of MG parameters

We derive results for the binning methods for the parametrizations ( µ, γ ) and ( µ, Σ). Our results are shown in Figs. 7 and 8, as well as Table V. Again, we find that binning methods complement very well functional methods. All results show that GR is consistent with current datasets.

and

## VIII. SUMMARY AND CONCLUDING REMARKS

We have described in this paper the new version of the ISiTGR patch to test deviations from GR at cosmological scales using various data sets. This version now assembles the following capabilities:

- 1. Dynamical dark energy parametrizations with a constant or time-dependent equation of state (also allowing the use of a pivot redshift to optimize constraints on its parameters).
- 2. A consistent implementation through all formalism to account for anistotropic shear to model, for example, contributions from massive neutrinos.
- 3. Spatially flat or curved backgrounds.
- 4. Multiple commonly used parametrizations of modified growth (MG) parameters to accommodate various types of data sets.
- 5. Functional and binned time- and scale-dependencies of the MG parameters.

As we are moving rapidly toward an era of precision cosmology with a plethora of incoming and future highly constraining data sets. It is important to have a software like ISiTGR that can constrain separately or simultaneously various extensions to the standard cosmological model.

Finally, We have also provided a number of results from applying ISiTGR to various available data sets and using various model extensions listed above. We find that GR is consistent with current data at cosmological scales.

The code is made publicly available and ready to use at https://github.com/mishakb/ISiTGR .

## Acknowledgments

We thank Shahab Joudaki, Weikang Lin and Eric Linder for useful comments on the manuscript. We thank Antony Lewis for clarifying some points about CosmoMC and Matteo Martinelli for clarifying some points about some of the MG results in [66]. MI acknowledges that this material is based upon work supported in part by the U.S. Department of Energy, Office of Science, under Award Number DE-SC0019206 and the National Science Foundation under grant AST-1517768. CGQ gratefully acknowledges a PhD scholarship from the Mexican National Council for Science and Technology (CONACYT).

The authors acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources that have contributed to the research results reported within this paper. URL: http://www.tacc.utexas.edu

## Appendix A: Spatially flat case

For convenience and clarity we write explicit forms of various equations for the spatially flat case ( K = 0).

## 1. Formalism and MG equations

For K = 0 the various Poisson-like equations for the metric potentials, and their corresponding modifications with MG parameters, are

## ( µ, γ ) parametrization:

$$k ^ { 2 } \Psi = - 4 \pi G a ^ { 2 } \mu ( a, k ) \sum _ { i } \left [ \rho _ { i } \Delta _ { i } + 3 \rho _ { i } ( 1 + w _ { i } ) \sigma _ { i } \right ]$$

$$k ^ { 2 } ( \Psi - \gamma ( a, k ) \Phi ) & = - 1 2 \pi G a ^ { 2 } \mu ( a, k ) \sum _ { i } \rho _ { i } ( 1 + w _ { i } ) \sigma _ { i }. & ( A 2 )$$

( µ, Σ) parametrization:

and

( Q,R ) parametrization:

and

( Q,D ) parametrization:

$$k ^ { 2 } \Phi = - 4 \pi G a ^ { 2 } Q ( a, k ) \sum _ { i } \rho _ { i } \Delta _ { i }$$

and

$$k ^ { 2 } ( \Phi + \Psi ) = - 8 \pi G a ^ { 2 } D ( a, k ) \sum _ { i } \rho _ { i } \Delta _ { i } - 1 2 \pi G a ^ { 2 } Q ( a, k ) \sum _ { i } \rho _ { i } ( 1 + w _ { i } ) \sigma _ { i }.$$

Note that Σ is defined directly as a modification to the unmodified version of (A4), while D is defined in terms of R and Q as D ≡ Q 2 ( R + 1). Σ can similarly be related to µ and γ , but the relation only takes the simple form Σ = µ 2 ( γ +1) in the limit of zero anisotropic shear ( σ i = 0) . It should also be noted that it is only in the case of zero shear that we have Σ = D .

## 2. Implementation

Following section V, we can write equations for the synchronous gauge potentials for K = 0 as

$$k ^ { 2 } ( \eta - \mathcal { H } \alpha ) = - 4 \pi G a ^ { 2 } \sum _ { i } \rho _ { i } \Delta _ { i }$$

and

$$k ^ { 2 } ( \dot { \alpha } + 2 \mathcal { H } \alpha - \eta ) = - 1 2 \pi G a ^ { 2 } \sum _ { i } \rho _ { i } ( 1 + w _ { i } ) \sigma _ { i }, & & ( A 1 0 )$$

where α = ( ˙ h +6˙) 2 η / k 2 , and again these potentials are related to the the Newtownian gauge potentials by

$$\Phi & = \eta - \mathcal { H } \alpha & ( \mathbf A 1 1 )$$

and

The expressions for ˙ η in the flat case become

$$k ^ { 2 } \Psi = - 4 \pi G a ^ { 2 } \mu ( a, k ) \sum _ { i } \left [ \rho _ { i } \Delta _ { i } + 3 \rho _ { i } ( 1 + w _ { i } ) \sigma _ { i } \right ]$$

$$k ^ { 2 } ( \Phi + \Psi ) = - 4 \pi G a ^ { 2 } \Sigma ( a, k ) \sum _ { i } \left [ 2 \rho _ { i } \Delta _ { i } + 3 \rho _ { i } ( 1 + w _ { i } ) \sigma _ { i } \right ].$$

$$k ^ { 2 } \Phi = - 4 \pi G a ^ { 2 } Q ( a, k ) \sum _ { i } \rho _ { i } \Delta _ { i }$$

$$k ^ { 2 } ( \Psi - R ( a, k ) \Phi ) & = - 1 2 \pi G a ^ { 2 } Q ( a, k ) \sum _ { i } \rho _ { i } ( 1 + w _ { i } ) \sigma _ { i }. & ( A 6 )$$

$$\Psi & = \dot { \alpha } + \mathcal { H } \alpha. & & ( \mathbf A ^ { } { 1 } 2 )$$

( µ, γ ) parametrization:

where

$$\iota _ { \mu, \gamma } & = \frac { 1 } { 2 f _ { \mu, \gamma } } \Big \{ k \mu \gamma f _ { 1 } \sum _ { i } q _ { i } ^ { ( S ) } \hat { \rho } _ { i } + \sum _ { i } \Delta _ { i } \hat { \rho } _ { i } [ \mathcal { H } \mu ( \gamma - 1 ) - \dot { \mu } \gamma - \mu \gamma ] \\ & + 2 \mu ( 1 - \gamma ) \sum _ { i } \mathbb { \Pi } _ { i } \hat { \rho } _ { i } + k ^ { 2 } \alpha [ - 2 ( \mathcal { H } ^ { 2 } - \dot { \mathcal { H } } ) + \mu \gamma \sum _ { i } \hat { \rho } _ { i } ( 1 + w _ { i } ) ] \\ & - 2 [ \mu \dot { \gamma } + \dot { \mu } ( \gamma - 1 ) ] \sum _ { i } \Pi _ { i } \hat { \rho } _ { i } + 2 \mathcal { H } \mu \sum _ { i } \Pi _ { i } \hat { \rho } _ { i } ( \gamma - 1 ) ( 3 w _ { i } + 2 ) \Big \},$$

$$f _ { \mu, \gamma } \equiv k ^ { 2 } + \frac { 3 } { 2 } \mu \gamma \sum _ { i } ( 1 + w _ { i } ) \hat { \rho } _ { i }.$$

( µ, Σ) parametrization:

$$\{ \mu, \mathcal { Z } \} \text{ parametricization} \colon \\ \hat { \eta } _ { ( \mu, \Sigma ) } & = \frac { 1 } { 2 f _ { \mu, \Sigma } } \left \{ k ( 2 \Sigma - \mu ) f _ { 1 } \sum _ { i } q _ { i } ^ { ( s ) } \hat { \rho } _ { i } + [ ( \mu - 2 \dot { \Sigma } ) + 2 \mathcal { H } ( \Sigma - \mu ) ] \sum _ { i } \Delta _ { i } \hat { \rho } _ { i } \\ & + 2 ( \mu - \Sigma ) \sum _ { i } \hat { \Pi } _ { i } \hat { \rho } _ { i } + 2 [ ( \mu - \dot { \Sigma } ) + \mathcal { H } ( 2 \Sigma - \mu ) - \mathcal { H } \mu ] \sum _ { i } \Pi _ { i } \hat { \rho } _ { i } \\ & + 2 \mathcal { H } ( \Sigma - \mu ) \sum _ { i } \Pi _ { i } ( 1 + 3 w _ { i } ) \hat { \rho } _ { i } + k ^ { 2 } \alpha \left [ ( 2 \Sigma - \mu ) \sum _ { i } ( 1 + w _ { i } ) \hat { \rho } _ { i } - 2 ( \mathcal { H } ^ { 2 } - \dot { \mathcal { H } } ) \right ] \right \}, \\ \text{where}$$

where

$$f _ { \mu, \Sigma } & \equiv k ^ { 2 } + \frac { 3 } { 2 } ( 2 \Sigma - \mu ) \sum _ { i } \hat { \rho } _ { i } ( 1 + w _ { i } ).$$

## ( Q,D ) parametrization:

where

$$f _ { Q } \equiv k ^ { 2 } + \frac { 3 } { 2 } Q \sum _ { i } ( 1 + w _ { i } ) \hat { \rho } _ { i }.$$

Note that ˙ η ( Q,R ) can be obtained via the equation R = 2 D Q -1.

Finally we need the derivative of the Weyl potential, Φ + Ψ, to compute ISW effect contributions. ˙ ˙ During the derivation of the expressions for ˙, η ˙ Φ was derived, so we just need expressions for Ψ. ˙ Explicitly, in each case, these are given by

( µ, γ ) or ( µ, Σ) parametrizations:

$$\dot { \Psi } & = - \frac { \dot { \mu } } { 2 k ^ { 2 } } \sum _ { i } [ \Delta _ { i } \hat { \rho } _ { i } + 2 \Pi _ { i } \hat { \rho } _ { i } ] + \frac { \mu } { 2 k ^ { 2 } } \sum _ { i } \left \{ \mathcal { H } \Delta _ { i } \hat { \rho } _ { i } - 2 \Pi _ { i } \hat { \rho } _ { i } + k f _ { 1 } q _ { i } ^ { ( S ) } \hat { \rho } _ { i } \\ & \quad + ( 1 + w _ { i } ) \hat { \rho } _ { i } [ k ^ { 2 } \alpha f _ { 1 } - 3 ( \dot { \Phi } + \mathcal { H } \Psi ) ] + 2 \mathcal { H } \Pi _ { i } \hat { \rho } _ { i } ( 1 + 3 w _ { i } ) + 2 \mathcal { H } \hat { \rho } _ { i } \Pi _ { i } \right \}. \\ \cdots \quad \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots.$$

Note that this expression works for both of these parametrizations since it only involves µ , however the explicit form of Φ will depend on the particular parametrization used. ˙

$$\dot { \eta } _ { ( Q, D ) } & = - \frac { 1 } { 2 f _ { Q } } \left \{ 2 k ^ { 2 } \alpha ( \mathcal { H } ^ { 2 } - \dot { \mathcal { H } } ) + [ 2 \mathcal { H } ( D - Q ) + \dot { Q } ] \sum _ { i } \Delta _ { i } \hat { \rho } _ { i } \\ & - k ^ { 2 } \alpha \sum _ { i } Q ( 1 + w _ { i } ) \hat { \rho } _ { i } - k Q f _ { 1 } \sum _ { i } q _ { i } ^ { ( S ) } \hat { \rho } _ { i } \right \},$$

## ( Q,D ) or ( Q,R ) parametrizations:

$$\dot { \Psi } = - \dot { \Phi } & - \frac { 1 } { k ^ { 2 } } \sum _ { i } \left \{ \dot { D } \Delta _ { i } \hat { \rho } _ { i } + \dot { Q } \Pi _ { i } \hat { \rho } _ { i } - D \mathcal { H } \Delta _ { i } \hat { \rho } _ { i } - Q \mathcal { H } \Pi _ { i } ( 3 w _ { i } + 1 ) \hat { \rho } _ { i } + Q \dot { \Pi } _ { i } \hat { \rho } _ { i } \\ & - 2 D \mathcal { H } \Pi _ { i } \hat { \rho } _ { i } - k D f _ { 1 } q _ { i } ^ { ( S ) } \hat { \rho } _ { i } + ( 1 + w _ { i } ) \hat { \rho } _ { i } [ 3 D ( \dot { \Phi } + \mathcal { H } \Psi ) - k ^ { 2 } \alpha f _ { 1 } D ] \right \}.$$

Note that this expression works for both of these parametrizations because the relation R = 2 D Q -1 may be used.

Appendix B: Parametrization summary table

| set of MG parameters   | MG first order perturbed equations                                                                                                              | functional form of MG parameters                                                                                                                                                                                                          |
|------------------------|-------------------------------------------------------------------------------------------------------------------------------------------------|-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| ( µ, γ )               | ( k 2 - 3 K )Ψ = - 4 πGa 2 µ ( a, k ) [ ρ ∆+3 k 2 - 3 K k 2 ( ρ + p ) σ ] k 2 [Φ - γ ( a, k )Ψ] = 12 πGa 2 µ ( a, k )( ρ + p ) σ                | µ = 1 + E 11 Ω DE ( a ) [ 1+ c 1 ( λH/k ) 2 1+( λH/k ) 2 ] γ = 1 + E 22 Ω DE ( a ) [ 1+ c 2 ( λH/k ) 2 1+( λH/k ) 2 ]                                                                                                                     |
| ( µ, Σ)                | ( k 2 - 3 K )Ψ = - 4 πGa 2 µ ( a, k ) [ ρ ∆+3 k 2 - 3 K k 2 ( ρ + p ) σ ] . k 2 [Φ +Ψ] = - 4 πGa 2 Σ( a, k ) [ 2 1 - 3 K/k 2 ρ ∆+3( ρ + p ) σ ] | µ = 1 + µ 0 Ω DE ( a ) Ω Λ [ 1+ c 1 ( λH/k ) 2 1+( λH/k ) 2 ] Σ = 1 +Σ 0 Ω DE ( a ) Ω Λ [ 1+ c 2 ( λH/k ) 2 1+( λH/k ) 2 ]                                                                                                                |
| ( Q,D )                | ( k 2 - 3 K )Φ = - 4 πGa 2 Q ( a, k ) ρ ∆ k 2 [Φ +Ψ] = - 8 πGa 2 D ( a, k ) ρ ∆ - 12 πGa 2 Q ( a, k )( ρ + p )                                  | Q = 1 + Q 0 Ω DE ( a ) Ω Λ [ 1+ c 1 ( λH/k ) 2 1+( λH/k ) 2 ] D = 1 + D 0 Ω DE ( a ) Ω Λ [ 1+ c 2 ( λH/k ) 2 1+( λH/k ) 2 ] σ Q = [ Q 0 e - k/k c + Q ∞ (1 - e - k/k c ) - 1] a s +1                                                      |
| ( Q,R )                | ( k 2 - 3 K )Φ = - 4 πGa 2 Q ( a, k ) ρ ∆ k 2 [Ψ - R ( a, k )Φ] = - 12 πGa 2 Q ( a, k )( ρ + p ) σ                                              | Q = 1 + Q 0 Ω DE ( a ) Ω Λ [ 1+ c 1 ( λH/k ) 2 1+( λH/k ) 2 ] R = 1 + R 0 Ω DE ( a ) Ω Λ [ 1+ c 2 ( λH/k ) 2 1+( λH/k ) 2 ] Q = [ Q 0 e - k/k c + Q ∞ (1 - e - k/k c ) - 1] a s +1 R = [ R 0 e - k/k c + R ∞ (1 - e - k/k c ) - 1] a s +1 |

TABLE VI: Summary table of the parametrizations used in ISiTGR.

FIG. 9: Diagram that illustrates the various parametrizations of the ISiTGR code.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00290v2-figure-a915191c-32ee-4582-9e39-3b2fd838d8cb.jpg)

FIG. 10: Flowchart that shows the different files that need to be modified in order to use the ISiTGR patch.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00290v2-figure-b929f72a-6abf-4415-bf42-a6a7b26483df.jpg)

- [1] Kilo-Degree Survey , URL http://kids.strw.leidenuniv.nl/ .
- [2] Dark Energy Survey , URL https://www.darkenergysurvey.org .
- [3] Hyper Suprime-Cam , URL https://hsc.mtk.nao.ac.jp/ssp/ .
- [4] Large Synoptic Survey Telescope , URL https://www.lsst.org/ .
- [5] Dark Energy Spectroscopic Instrument , URL https://www.desi.lbl.gov/ .
- [6] Euclid , URL http://www.euclid-ec.org/ .
- [7] Square Kilometre Array , URL https://www.skatelescope.org/ .
- [8] Wide Field Infrared Spectroscopic Telescope , URL https://wfirst.gsfc.nasa.gov/ .
- [9] Simons Observatory , URL https://simonsobservatory.org/ .
- [10] CMB-S4 , URL https://cmb-s4.org/ .
- [11] T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, Phys. Rep. 513 , 1 (2012), 1106.2476.
- [12] K. Koyama, Rep. Progr. Phys. 79 , 046902 (2016), 1504.04623.
- [13] A. Joyce, B. Jain, J. Khoury, and M. Trodden, Phys. Rep. 568 , 1 (2015), 1407.0059.
- [14] A. Joyce, L. Lombriser, and F. Schmidt, Annu. Rev. Nucl. Part. Sci. 66 , 95 (2016), 1601.06133.
- [15] M. Ishak, Living Reviews in Relativity 22 , 1 (2019), 1806.10122.
- [16] S. Weinberg, Rev. Mod. Phys. 61 , 1 (1989).
- [17] S. M. Carroll, Living Rev. Relativ. 4 , 1 (2001), astro-ph/0004075.
- [18] V. Sahni and A. Starobinsky, Int. J. Mod. Phys. D 9 , 373 (2000), astro-ph/9904398.
- [19] P. J. E. Peebles and B. Ratra, Rev. Mod. Phys. 75 , 559 (2003), astro-ph/0207347.
- [20] M. Ishak, Found. Phys. 37 , 1470 (2007), astro-ph/0504416.
- [21] D. H. Weinberg, M. J. Mortonson, D. J. Eisenstein, C. Hirata, A. G. Riess, and E. Rozo, Phys. Rep. 530 , 87 (2013), 1201.2434.
- [22] D. Huterer and D. L. Shafer, Rep. Progr. Phys. 81 , 016901 (2018), 1709.01091.
- [23] E. V. Linder, Phys. Rev. D 72 , 043529 (2005), astro-ph/0507263.
- [24] K. Koyama, J. Cosmol. Astropart. Phys. 3 , 017 (2006), astro-ph/0601220.
- [25] P. Zhang, M. Liguori, R. Bean, and S. Dodelson, Phys. Rev. Lett. 99 , 141302 (2007), 0704.1932.
- [26] R. Caldwell, A. Cooray, and A. Melchiorri, Phys. Rev. D 76 , 023507 (2007), astro-ph/0703375.
- [27] E. V. Linder and R. N. Cahn, Astroparticle Physics 28 , 481 (2007), astro-ph/0701317.
- [28] D. Polarski and R. Gannouji, Phys. Lett. B 660 , 439 (2008), 0710.1510.
- [29] Y. Gong, M. Ishak, and A. Wang, Phys. Rev. D 80 , 023002 (2009), 0903.0001.
- [30] G.-B. Zhao, L. Pogosian, A. Silvestri, and J. Zylberberg, Phys. Rev. D 79 , 083513 (2009), 0809.3791.
- [31] V. Acquaviva and E. Gawiser, Phys. Rev. D 82 , 082001 (2010), 1008.3392.
- [32] R. Bean and M. Tangmatitham, Phys. Rev. D 81 , 083534 (2010), 1002.4197.
- [33] L. Lombriser, Phys. Rev. D 83 , 063519 (2011), 1101.0594.
- [34] J. N. Dossett, J. Moldenhauer, and M. Ishak, Phys. Rev. D 84 , 023012 (2011), 1103.1195.
- [35] G. Dvali, G. Gabadadze, and M. Porrati, Phys. Lett. B 485 , 208 (2000), hep-th/0005016.
- [36] M. Ishak, A. Upadhye, and D. N. Spergel, Phys. Rev. D 74 , 043513 (2006), astro-ph/0507184.
- [37] E. J. Ruiz and D. Huterer, Phys. Rev. D 91 , 063009 (2015), 1410.5832.
- [38] J. L. Bernal, L. Verde, and A. G. Riess, J. Cosmol. Astropart. Phys. 10 , 019 (2016), 1607.05617.
- [39] W. Lin and M. Ishak, Phys. Rev. D 96 , 023532 (2017), 1705.05303.
- [40] C. Garcia-Quintero, M. Ishak, L. Fox, and W. Lin, arXiv e-prints arXiv:1910.01608 (2019), 1910.01608.
- [41] J. N. Dossett, M. Ishak, and J. Moldenhauer, Phys. Rev. D 84 , 123001 (2011), 1109.4583.
- [42] J. N. Dossett and M. Ishak, Phys. Rev. D 86 , 103008 (2012), 1205.2422.
- [43] A. Lewis, A. Challinor, and A. Lasenby, Astrophys. J. 538 , 473 (2000), astro-ph/9911177.
- [44] A. Lewis and S. Bridle, Phys. Rev. D 66 , 103511 (2002), astro-ph/0205436.
- [45] S. Joudaki, A. Mead, C. Blake, et al., Mon. Not. R. Astron. Soc. 471 , 1259 (2017), 1610.04606.
- [46] S. Joudaki, C. Blake, A. Johnson, et al., Mon. Not. R. Astron. Soc. 474 , 4894 (2018), 1707.06627.

[47]

J. N. Dossett, M. Ishak, D. Parkinson, and T. M. Davis, Phys. Rev. D

92

,

023003 (2015), 1501.03119.

- [48] A. Zucca, L. Pogosian, A. Silvestri, and G. B. Zhao, J. Cosmol. Astropart. Phys. 2019 , 001 (2019), 1901.05956.
- [49] M. Zumalac´ arregui, E. Bellini, I. Sawicki, J. Lesgourgues, and P. G. Ferreira, J. Cosmol. Astropart. Phys. 8 , 019 (2017), 1605.06102.
- [50] J. Lesgourgues, ArXiv e-prints (2011), 1104.2932.
- [51] D. Blas, J. Lesgourgues, and T. Tram, J. Cosmol. Astropart. Phys. 7 , 034 (2011), 1104.2933.
- [52] G. W. Horndeski, Int. J. Theor. Phys. 10 , 363 (1974).
- [53] B. Hu, M. Raveri, N. Frusciante, and A. Silvestri, ArXiv e-prints (2014), 1405.3590.
- [54] B. Hu., M. Raveri, N. Frusciante, and A. Silvestri, Phys. Rev. D 89 , 103530 (2014), 1312.5742.
- [55] F. Pace, R. A. Battye, B. Bolliet, and D. Trinh, J. Cosmol. Astropart. Phys. 2019 , 018 (2019), 1905.06795.
- [56] E. Bellini, A. Barreira, N. Frusciante, B. Hu, S. Peirone, M. Raveri, M. Zumalac´rregui, A. Avilez-Lopez, M. Ballardini, a R. A. Battye, et al., Phys. Rev. D 97 , 023520 (2018), 1709.09135.
- [57] T. Mueller and F. Grave (2009), 0904.4184.
- [58] L. F. Abbott and R. K. Schaefer, Astrophys. J. 308 , 546 (1986).

- [59] M. Zaldarriaga, U. Seljak, and E. Bertschinger, Astrophys. J. 494 , 491 (1998), astro-ph/9704265.
- [60] C.-P. Ma and E. Bertschinger, Astrophys. J. 455 , 7 (1995), astro-ph/9506072.
- [61] M. Chevallier and D. Polarski, Int. J. Mod. Phys. D 10 , 213 (2001), gr-qc/0009008.
- [62] E. V. Linder, Phys. Rev. Lett. 90 , 091301 (2003), astro-ph/0208512.
- [63] A. Albrecht, G. Bernstein, R. Cahn, W. L. Freedman, J. Hewitt, W. Hu, J. Huth, M. Kamionkowski, E. W. Kolb, L. Knox, et al., ArXiv e-prints (2006), astro-ph/0609591.
- [64] A. Albrecht, L. Amendola, G. Bernstein, D. Clowe, D. Eisenstein, L. Guzzo, C. Hirata, D. Huterer, R. Kirshner, E. Kolb, et al., ArXiv e-prints (2009), 0901.0721.
- [65] E. V. Linder, Astropart. Phys. 26 , 102 (2006), 0604280.
- [66] P. A. R. Ade, others, and (Planck Collaboration), Astron. Astrophys. 594 , A14 (2016), 1502.01590.
- [67] Planck Collaboration, N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. B. Barreiro, and N. Bartolo, arXiv e-prints arXiv:1807.06209 (2018), 1807.06209.
- [68] F. Simpson, C. Heymans, D. Parkinson, C. Blake, M. Kilbinger, J. Benjamin, T. Erben, H. Hildebrandt, H. Hoekstra, T. D. Kitching, et al., Mon. Not. R. Astron. Soc. 429 , 2249 (2013), 1212.3339.
- [69] E. V. Linder, Phys. Rev. D 95 , 023518 (2017), 1607.03113.
- [70] E. V. Linder, J. Cosmol. Astropart. Phys. 2019 , 034 (2019), 1903.02010.
- [71] J. N. Dossett and M. Ishak, Phys. Rev. D 88 , 103008 (2013), 1311.0726.
- [72] E. Bertschinger and P. Zukin, Phys. Rev. D 78 , 024015 (2008), 0801.2431.
- [73] DES Collaboration, T. M. C. Abbott, F. B. Abdalla, S. Avila, M. Banerji, E. Baxter, K. Bechtol, M. R. Becker, E. Bertin, and J. Blazek, arXiv e-prints arXiv:1810.02499 (2018), 1810.02499.
- [74] G.-B. Zhao, T. Giannantonio, L. Pogosian, A. Silvestri, D. J. Bacon, K. Koyama, R. C. Nichol, and Y.-S. Song, Phys. Rev. D 81 , 103510 (2010), 1003.0001.
- [75] Y.-S. Song, G.-B. Zhao, D. Bacon, K. Koyama, R. C. Nichol, and L. Pogosian, Phys. Rev. D 84 , 083523 (2011), 1011.2106.
- [76] S. F. Daniel, E. V. Linder, T. L. Smith, R. R. Caldwell, A. Cooray, A. Leauthaud, and L. Lombriser, Phys. Rev. D 81 , 123508 (2010), 1002.1962.
- [77] S. F. Daniel and E. V. Linder, Phys. Rev. D 82 , 103523 (2010), 1008.0397.
- [78] E. Krause, T. Eifler, and J. Blazek, Mon. Not. R. Astron. Soc. 456 , 207 (2016), 1506.08730.
- [79] T. M. C. Abbott, others, and (DES Collaboration), ArXiv e-prints (2017), 1708.01530.
- [80] P. A. R. Ade, others, and (Planck Collaboration), Astron. Astrophys. 594 , A13 (2016), 1502.01589.
- [81] P. A. R. Ade, others, and (Planck Collaboration), Astron. Astrophys. 594 , A15 (2016), 1502.01591.
- [82] D. M. Scolnic, D. O. Jones, A. Rest, et al., ArXiv e-prints (2017), 1710.00845.
- [83] S. Alam, M. Ata, S. Bailey, et al., Mon. Not. R. Astron. Soc. 470 , 2617 (2017), 1607.03155.
- [84] F. Beutler, C. Blake, M. Colless, D. H. Jones, L. Staveley-Smith, L. Campbell, Q. Parker, W. Saunders, and F. Watson, Mon. Not. R. Astron. Soc. 416 , 3017 (2011), 1106.3366.
- [85] A. J. Ross, L. Samushia, C. Howlett, W. J. Percival, A. Burden, and M. Manera, Mon. Not. R. Astron. Soc. 449 , 835 (2015), 1409.3242.
- [86] L. Anderson, E. Aubourg, S. Bailey, et al., Mon. Not. R. Astron. Soc. ´ 441 , 24 (2014), 1312.4877.
- [87] M. Betoule, J. Marriner, N. Regnault, J.-C. Cuillandre, P. Astier, J. Guy, C. Balland, P. El Hage, D. Hardin, R. Kessler, et al., Astron. Astrophys. 552 , A124 (2013), 1212.4864.
- [88] A. G. Riess, L. Macri, S. Casertano, H. Lampeitl, H. C. Ferguson, A. V. Filippenko, S. W. Jha, W. Li, and R. Chornock, Astrophys. J. 730 , 119 (2011), 1103.2976.
- [89] J. Zuntz, M. Paterno, E. Jennings, D. Rudd, A. Manzotti, S. Dodelson, S. Bridle, S. Sehrish, and J. Kowalkowski, Astron.Comput. 12 , 45 (2015), 1409.3409.
- [90] A. Hojjati, L. Pogosian, and G.-B. Zhao, J. Cosmol. Astropart. Phys. 8 , 005 (2011), 1106.4543.
Cristhian Garcia-Quintero Mustapha Ishak Logan Fox Jason Dossett
ISiTGR: Testing deviations from GR at cosmological scales including dynamical dark energy, massive neutrinos, functional or binned parametrizations, and spatial curvature
We introduce a new version of the Integrated Software in Testing General
Relativity (ISiTGR) which is a patch to the software CAMB and CosmoMC. ISiTGR
is intended to test deviations from GR at cosmological scales using
cosmological data sets. While doing so, it allows for various extensions to the
standard flat $\Lambda$CDM model. In this new release, we have support for the
following: 1) dynamical dark energy parametrizations with a constant or
time-dependent equation of state; 2) a consistent implementation of anisotropic
shear to model massive neutrinos throughout the full formalism; 3) multiple
commonly-used parametrizations of modified growth (MG) parameters; 4)
functional, binned and hybrid time- and scale-dependencies for all MG
parameters; 5) spatially flat or curved backgrounds. ISiTGR is designed to
allow cosmological analyses to take full advantage of ongoing and future
surveys to test simultaneously or separately various extensions to the standard
model. We describe here the formalism and its implementation in the CMB code,
the Integrated Sachs-Wolfe (ISW) effect, and the 3x2 point statistics. Next, we
apply ISiTGR to current data sets from Planck-2018, Planck-2015, Dark Energy
Survey YR1 release, Baryonic Acoustic Oscillations (BAO), Redshift Space
Distortions (BAO/RSD) from the BOSS Data Release 12, the 6DF Galaxy Survey and
the SDSS Data Release 7 Main Galaxy Sample, and Supernova from the Pantheon
compilation, joint SNLS/SDSS data analysis and the Hubble Space Telescope. We
derive constraints on MG parameters for various combinations of the five
features above and find that GR is consistent with current data sets in all
cases. The code is made publicly available at
\url{https://github.com/mishakb/ISiTGR}.