Searching for scalar gravitational interactions in current and future cosmological data

Modified gravity theories often contain a scalar field of gravitational strength which interacts with matter. We examine constraints on the range and the coupling strength of a scalar gravitational degree of freedom using a subset of current data that can be safely analyzed within the linear perturbation theory. Using a model-independent implementation of scalar-tensor theories in MGCAMB in terms of two functions of the scale factor describing the mass and the coupling of the scalar degree of freedom, we derive constraints on the $f(R)$, generalized chameleon, Symmetron and Dilaton models. Since most of the large scale structure data available today is from relatively low redshifts, only a limited range of observed scales is in the linear regime, leading to relatively weak constraints. We then perform a forecast for a future large scale structure survey, such as Large Synoptic Survey Telescope (LSST), which will map a significant volume at higher redshifts, and show that it will produce much stronger constraints on scalar interactions in specific models. We also perform a principal component analysis and find that future surveys should be able to provide tight constraints on several eigenmodes of the scalar mass evolution.


I. INTRODUCTION
A non-vanishing cosmological constant, Λ, is the simplest and the most common explanation of the observed cosmic acceleration [1,2]. Because, gravitationally, Λ is equivalent to the large vacuum energy predicted in particle physics, its value requires a technically unnatural fine-tuning [3,4] in order to be consistent with observations. The cosmological constant could be embedded in a larger class of dark energy models, where dynamics dictate the value of the vacuum energy. Because of the absence of apparent violation of Lorentz invariance in the Universe, dark energy is commonly described by the field theory of a scalar. Usually, some degree of fine-tuning of the parameters of the model must be introduced.
Another explanation could be provided by a modification of the laws of gravity on large scales. Such modifications generically involve a scalar degree of freedom which can lead to dynamical dark energy when the range of the scalar interaction is cosmological. As a result, scalartensor models with couplings to matter represent a wellmotivated and versatile class of dark energy. Theories describing the behaviour of the scalar field involve conformal [5] and disformal couplings to matter [6,7]. It turns out that the disformal coupling is severely constrained by local experiments and cosmological observations [8,9]. On the other hand, the conformal couplings, albeit large on cosmological scales, can be screened in the local environment where none of their effects, such as deviations from Newton's law, have been uncovered.
In this paper, we will focus on scalar-tensor models with screening mechanisms that are broadly classified to be of chameleon type [10,11], i.e. where either the mass of the scalar and/or its coupling to matter has a dependence of the local matter density. Specifically, we will consider three types of models with the chameleon property: the f (R) theories [12,13], the environmentally dependent Dilatons [14] and the Symmetron [15]. The latter two models use the Damour-Polyakov mechanism for screening [16]. We will take advantage of the fact that these three very different types of models can be described using the same formalism defined in terms of two dynamical functions m(a) and β(a), where a is the scale factor [17,18]. The first one represents the mass of the scalar in the cosmological background at the redshift 1 + z = a −1 and the second one is the coupling of the scalar to matter. The growth of cosmological perturbations in these models in the linear regime and on sub-horizon scales can be entirely described using a single function, (k, a) = 2β 2 (a)/[1 + m 2 (a)a 2 /k 2 ], which appears in the modification of Newton's constant and in the modified relation between the curvature and the gravitational potential. While in this paper we shall restrict ourselves to observables which are sensitive to the linear regime only, we note that, given m(a) and β(a), one can also reconstruct the full non-linear dynamics of the models. Namely, using a known evolution of the background matter density, ρ(a), one can express the mass and the coupling as functions of local matter density: m(ρ(x, t)) and β(ρ(x, t)) and use them to perform N-body simulations of these models or to analyse local gravitational tests.
Modified gravity (MG) and its comparison with dark energy has been investigated using various cosmological probes in the last ten years . Some models of modified gravity, such as the f (R) theories, have been strongly constrained by observations both cosmological and astrophysical. The strongest bound on the range of the scalar interaction, expressed in terms of the parameter f R0 , is at the level of 10 −7 and comes from astrophysical tests of modified gravity using the period of cepheids or the gas dynamics of dwarf galaxies [45][46][47]. The cosmological bounds are less effective, at the level of 10 −5 [43,48]. On the other hand, dilatons and symmetrons have not been constrained as systematically as f (R) on cosmological scales. Only a few tests have been performed using the [m(a), β(a)] parameterisation [42]. The strongest bounds on dilatons and symmetrons still spring from local gravitational tests [49]. Local tests of gravity for the chameleon-type models of modified gravity imply that the range of the scalar interaction cannot exceed 1 Mpc, implying that linear analyses are limited to probing only some of the features of the chameleon screening mechanisms. On the other hand, studying effects of modified gravity on shorter scales requires the use of either semianalytical methods suited to the quasi-linear regime of cosmological perturbations or N-body simulations, both of which are model-specific. Here, in order to keep our analysis as model-independent as possible, we shall restrict ourselves to observables which are sensitive to the linear regime only.
The range of scales that are safely in the linear regime at low redshifts is quite limited. Most of the large scale structure data available today is from relatively low redshifts and provides only weak constraints on scalar-tensor models unless one considers information from non-linear scales. The only way to do so is to run N-body simulations for specific models. On the other hand, future surveys, such as LSST [50] and Euclid [51], will provide a high volume of data from higher redshifts at which the range of linear scales is significantly larger, allowing one to deduce stronger constraints on scalar interactions not only for specific models but in a more general model-independent way. In this paper, we start by deriving constraints on f (R), Symmetron and Dilatons models from the subset of today's data that can be safely considered to be in the linear regime. Then we perform a Fisher forecast for the same models assuming data from a future LSST-like survey in combination with other types of data expected over the next 5-10 years to show that they will be significantly tighter. Finally, we perform a principal component analysis (PCA) forecast of m(a) for the same future data, assuming that β(a) is a slowly varying function that can be taken to be a O(1) constant over the range of redshifts relevant to LSST.

II. THE MODEL
We consider scalar-tensor theories defined by the action where g µν is the Einstein frame metric, ψ are the matter fields that follow geodesics of A 2 (φ)g µν , and L φ is the scalar field Lagrangian given by The action in Eq. (1) is a Generalized Brans-Dicke (GBD) theory [5] that includes a potential for the scalar field. In all GBD, the scalar field mediates an additional gravitational interaction between massive particles. The net force on a test mass is given by where Ψ is the Newtonian potential. Since solar system and laboratory tests severely constrain the presence of the scalar force, GBD can only be viable if either the coupling of the scalar field to matter is always negligible, or if there is a dynamical screening mechanism that suppresses the force in dense environments. The latter can be accomplished with appropriately chosen functional forms of A(φ) and V (φ). Because of its coupling to matter, the scalar field dynamics are determined by an effective potential which takes into account the presence of the conserved matter density ρ of the environment For some forms of V (φ) and A(φ), the effective potential can have a density dependent minimum, φ(ρ). The scalar force will be screened if either the mass of the field happens to be extremely large or the coupling happens to be negligibly small at the minimum of V eff (φ). Such models can be broadly classified as "Generalized Chameleons" (GC), and include the original chameleon model [10], f (R), dilatons [14] and symmetrons [15]. We note that the GC scalar-tensor theories considered in this work are viable only if the field stays at the minimum of the effective potential V eff (φ) [18]. In this case, the effective dark energy equation of state is indistinguishable from −1 and the expansion history practically the same as in the ΛCDM model. Furthermore, as long as the scalar field is at its density dependent minimum, φ(ρ), the theory can be described parametrically from the sole knowledge of the mass function m(ρ) and the coupling β(ρ) at the minimum of the potential [17,18] where we have identified the mass as the second derivative and the coupling It is often simpler to characterize the functions m(ρ) and β(ρ) using the time evolution of the matter density of the Universe where a is the scale factor whose value now is a 0 = 1. This allows one to describe characteristic models in a simple way and the full dynamics can be recovered from the time evolution of the mass and coupling functions, m(a), β(a).

A. Evolution of linear perturbations
While the scalar-tensor theories considered in this work predict the same expansion history as ΛCDM, the existence of the additional scalar interaction gives them distinguishing features in the evolution of linear matter and metric perturbations. More specifically, the attractive force mediated by the scalar enhances the overall growth of inhomogeneities. In addition, the relation between the curvature perturbation Φ and the Newtonian potential Ψ is modified [52]. Both of these effects can be captured in terms of two phenomenological functions employed in MGCAMB 1 [25,34,54], parametrizing effective modifications to the Poisson and the anisotropy Einstein equations in Fourier space. Namely, one defines µ(a, k) and γ(a, k), such that where ∆ is the comoving matter density contrast 2 . In the quasi-static approximation, whose validity is discussed below, functions µ(k, a) and γ(k, a) can be expressed in terms of m(a) and β(a) as [18].
The conformal factor A 2 (φ) that appears in Eq. (11) is indistinguishable from unity for viable models within the class of scalar-tensor theories considered in this paper, and can be safely ignored. ΛCDM is recovered when → 0 and µ = γ = 1.
1 MGCAMB is a publicly available patch to CAMB [53]. 2 These equations are valid at late times, when the contribution of relativistic species can be neglected. Equations used in MG-CAMB are more general and are valid at all times [35].
In the quasi-static approximation, the equation governing the evolution of matter density contrast δ reads δ + Hδ − 3 2 Ω m H 2 µ(k, a)δ = 0 (14) where is the derivative with respect to conformal time and H = a /a. Two regimes can be distinguished. When the mode k is outside the Compton wavelength of the scalar field, i.e. k am(a), 1 and the growth is not modified. Inside the Compton wavelength, k am(a), gravity is enhanced by 1+2β 2 (a), implying more growth. In addition, in the Symmetron and Dilatons models, the coupling β(a) depends on the matter density and controls the transition to the enhanced growth.
Since is a manifestly non-negative number, the growth is generically enhanced. Also, generically, γ < 1 in these models. At the same time, the relation between the lensing potential Φ + Ψ and the matter density is effectively unchanged. Namely, if one defines Σ(k, a) as then Σ = A 2 (φ) (16) and is effectively unity for all viable modes. Thus, a clear detection of Σ = 1 would not only signal a breakdown of ΛCDM but would rule out the entire class of GBD models. We note that, even though Σ is constrained to be very close to unity in viable GBD models, its time derivative, Σ, can, in principle, be non-negligible and affect the observables via the Integrated Sachs-Wolfe (ISW) effect. When functions m(a) and β(a) are regular, which is the case for chameleon models such as f (R), and for dilatons, the error introduced by working in the quasi-static approximation scales as H/k [55]. For models such as the Symmetron, in which the functions m(a) and β(a) vanish with a power n < 1 for a > a and are zero for a < a (and thus have a diverging derivative at a ), the accuracy is reduced to (H/k) n [55,56]. In what follows, we briefly motivate specific functional forms of (m(a) and β(a) adopted for the analysis in Sections III and IV. Given the forms of m(a) and β(a), the predictions for the observables can be calculated using MGCAMB [54]. Plots of the CMB temperature anisotropy and the matter power spectra for a few representative models are shown in Fig. 1.
Among theories exhibiting chameleon screening are the f (R) class of models [57,58] described by the action where the function f (R) is designed to depart from the Einstein-Hibert form at smaller values of the curvature R. As a specific example, we take the form proposed by Hu and Sawicki (HS) [12], where Λ is the cosmological constant term, R 0 is the value of the curvature today and f R0 ≡ (1 − df /dR) R=R0 . As argued in [12,59,60], all viable f (R) models should be of such "disappearing cosmological constant" type [60], and models similar to HS were proposed in [59,60]. For all f (R) models, β(a) = 1/ √ 6, while the mass function is model dependent. In the HS model, we have where Ω Λ and Ω m are the dark energy and matter density fractions today, and m 0 is a mass scale that can be expressed in terms of f R0 as [18] Local tests of gravity require f R0 < ∼ 10 −6 [49], while astrophysical constraints from dwarf galaxies imply that f R0 < ∼ 10 −7 [47]. These bounds depend on accurate modelling of non-linear physics. In what follows, we will derive the constraint on f R0 from current cosmological data using only information from linear scales, and also forecast constraints expected from future surveys like LSST.
Representative CMB and matter power spectra for f (R) are shown in Fig. 1. A notable effect on the CMB spectrum is the suppression of power at small multipoles, which is due to the reduced Integrated Sachs-Wolfe (ISW) effect. The magnitude of the ISW effect is proportional to the net change in the gravitational potential along the line of sight. In ΛCDM, the change in the potential is a reduction caused by the onset of cosmic acceleration. In f (R), the additional scalar force enhances the potential which, combined with the decay due to acceleration, leads to a smaller net change and, thus, a smaller ISW effect. The other notable impact of f (R) on the CMB spectrum is the enhanced lensing, which has the effect of slightly dumping the peaks. The enhanced growth is more evident in the plot of P (k). Qualitatively, these features are common to all GBD models. Another relevant example is the environmentally dependent Dilaton [14], where the screening mechanism is of the Damour-Polyakov type [16]. This model, inspired by string theory in the large string coupling limit, has an exponentially runaway potential with the value of V 0 set to generate the current acceleration of the Universe, while the coupling function is In dense environments, the minimum of the effective potential approaches φ = φ , and the coupling function β(a) vanishes. The coefficient A 2 has to be large to satisfy local tests of gravity; typically A 2 > ∼ 10 6 . These models can be described by a mass function given by and, assuming matter domination, a coupling function where β 0 = Ω Λ /Ω m ∼ 2.7 is related to V 0 , and is determined by requiring that φ plays the role of dark energy. We will present our constraints on the mass in terms of a scalar-force range parameter ξ 0 , defined as where m 0 = m(a = 1). We show representative CMB and matter power spectra for the Dilaton model in Fig. 1, with parameter values being large on purpose to exaggerate the qualitative features of the model.

D. Functions m(a) and β(a) for symmetrons
Another example of a GBD model with the Damour-Polyakov screening mechanims is the Symmetron [15], where the scalar field has a quartic potential, and a coupling function, When matter density is large, the effective potential has a minimum at φ = 0 and A(φ) → 1, thus decoupling the scalar from matter. At lower densities, the effective potential acquires a non-zero minimum, activating the scalar force. For cosmological densities, the transition occurs at where ρ m is the matter density today. Thus, one can work with a , along with m and β , as the three free parameters of the theory. At a > a , the model can be described by and while β(a) = 0 for a < a . As in the case of dilatons, we represent our bounds in terms of a range parameter ξ , defined as Representative CMB and matter power spectra for this model are shown in Fig. 1.

E. Generalized Chameleon models
In our forecasts, we will also consider generalized models of chameleon type [61] defined by In practically all viable chameleon models, the coupling function is expected to vary extremely slowly at redshifts probed by large scale structure surveys. Thus, for all practical purposes, it can be taken to be a constant of order unity.

F. Binned Model
As discussed so far, for any of the aforementioned models, each with its own theoretical motivation, one can determine the functional forms of m(a) and β(a). This effectively reduces the two free functions m(a) and β(a) to a handful of parameters. However, one might be interested in knowing how well the two functions are constrained in general, without regard for any specific model. One can then proceed by discretizing either of the two functions in bins of redshift space and treating the amplitude in each bin as a free parameter to be constrained.
Varying both, the coupling and the mass functions, simultaneously would be redundant, since their effect is largely degenerate. Since it is the mass parameter that affects the shape of the matter power spectrum, we fix β(a) to a constant value of order unity and bin m(a) in redshift. If a non-zero m −1 (a) were detected, it would signal the presence of a scalar interaction and further investigation would be required to determine if the variation occurs in β(a), m(a) or both.
While a binning scheme gives a model independent (rather a far less model dependent) treatment of m(a), the larger number of parameters (values of m in each bin) results in weaker constraints on the individual parameters. To extract useful information, we apply the Principal Components Analysis (PCA) technique (reviewed in Section IV G). The resulting Principal Components (PCs) are linear combinations of the original bin values and the propagated uncertainty (from original errors on the bins) in their values can inform us about those PCs that are best constrained by data and the number of degrees of freedom the can potentially be constrained.

III. CONSTRAINTS FROM CURRENT DATA
In this Section, we use a combination of currently available CMB, lensing and Baryonic Acoustic Oscillation (BAO) data, as well as measurements of the matter power spectrum, to derive constraints on the GBD parameters. To compute the observables, we implemented the parametrizations described in the previous Section in MGCAMB. We then use it with an appropriately modified version of CosmoMC [62] to obtain the posterior distributions for the model parameters. Since current data is unable to simultaneously constrain multiple GBD parameters, we will only consider models from the previous Section for which meaningful constraints are possible.

A. The datasets used in the analysis
We use the measurements of CMB temperature anisotropy from the second data release of the Planck survey [63] in the form of the full Planck TT high-likelihood (30 < < 2500) along with the low-polarization ( < 30). We refer to the above datasets as PLC. We also consider the Planck 2015 lensing potential spectrum [64] extracted from mode-coupling correlations, and refer to this dataset as CMBLens.
In addition to inducing higher order correlations, lensing by large scale structures affects the TT spectrum at higher , slightly damping the oscillatory features. In [65], and subsequently in [63], the lensing contribution to TT was quantified via an amplitude A L multiplying the lensing power spectrum in the calculation of the theoretical prediction for TT. The parameter A L was used to quantify the significance of detection of the lensing contribution to TT. However, instead of measuring the expected value of A L = 1, since the lensing contribution to TT is calculated from the same model as the rest of the spectrum, the best fit value obtained for LCDM from the PLC dataset in [63] was A L = 1.22 ± 0.10, or two stan- dard deviations away from the expectation. As discussed in [65] (see also [66]) this is due to an apparent tension between the higher-and lower-data when trying to fit LCDM to Planck TT data. To negate the effect of this tension, the parameter A L was sometimes co-varied with other parameters when deriving constraints on LCDM in [63]. In what follows, we take the view that A L is not a physical parameter and should be held fixed to 1 when deriving constraints on cosmological models. However, we also investigate and discuss the effect of co-varying A L in the case of f (R).
For BAO measurements, we used data from the 6dF survey [67] and from SDSS, specifically the MGS [68] and BOSS data releases (LOWZ and CMASS) [69].
We also use the matter power spectrum (referred to as MPK) from SDSS LRG DR4 [70], but only on linear scales, k ≤ 0.1 h/Mpc. We are aware of the fact that nonlinear corrections can play a role even at k < ∼ 0.1Mpc and that a proper treatment of the bias and the redshift space distortions (RSD) must take them into account. This was studied at length in [71] for the SDSS DR9 power spectrum and it was found that the differences in the upper bounds on neutrino masses obtained using four different RSD models were under 20%. Based on this, we expect that bounds on the GBD parameters (such as f R0 ) obtained from MPK are accurate to within 30%, which is sufficient given that constraints form current data are relatively weak.
Finally, we consider the weak lensing data from the Canada France Hawaii Telescope Lensing Survey 2DCFHTLenS [72], referred to as WL. To avoid dealing with non-linear scales, we adopt a conservative cut and exclude θ < 30 from the measurements of the correlation function ξ ± , which corresponds to k < 0.1 h/Mpc scales.

B. Constraints on f (R)
The Hu-Sawicki f (R) model has two parameters, f R0 and n. In what follows, we fix n = 1 because that is a common choice in the literature, and also because the two parameters are highly correlated and the current data cannot simultaneously constrain both. We chose a flat prior on log 10 f R0 within the [−7, 0] range. We have checked that changing the range of the flat prior does not mν , we see that the constraint on fR 0 become weaker when the neutrino mass is varied. The datasets are labeled according to the notation introduced in Sec. III A. The symbol + means that we add data on top of the PLC+BAO dataset. For example, +lensing means PLC+BAO+lensing. affect our results. Fig. 2 shows constraints on f R0 for different combinations of datasets described in Sec. III A, after marginalizing over all the other cosmological parameters. We considered the case in which the total neutrino mass is fixed at m ν = 0.06 eV (solid lines), and the case where it can vary within 0 ≤ m ν ≤ 1 eV (dashed lines). The results from Fig. 2 are summarized in Table I.
We can see that the combination of PLC and BAO datasets (blue lines) only weakly constrains the model. Modified gravity affects the CMB temperature anisotropy spectrum in two ways: it affects the lowpower spectrum through the ISW effect and enhances the damping at high-due to the enhancement in clustering and, as a consequence, the lensing potential. Thus, the observed lack of power at low-multipoles and the apparent preference of enhanced lensing in CMB TT, when compared to the ΛCDM prediction, can be reconciled by a non-zero f R0 . This is the reason for the peak in  Fig. 2 the PLC+BAO likelihood. Adding the CMBLens data (red lines) tightens the constraint substantially. The enhancement of growth due to the extra scalar interaction affects the lensing potential measured by Planck, which is known to be in excellent agreement with the LCDM prediction [64]. Thus, the weak preference for larger f R0 coming from PLC+BAO is overwhelmed by the stronger CMBLens data that is consistent with f R0 = 0. The constraint becomes even tighter after adding the MPK and WL datasets (green lines).
The dashed lines in Fig. 2 show the impact of covarying the combined mass of neutrinos, m ν , along with f R0 . Massive neutrinos suppress the growth and can partially compensate for the enhanced clustering in f (R), slightly weakening the bounds on f R0 . The extent of the degeneracy can be inferred from Fig. 3 which shows the joint confidence contours for the two parameters. We see that, although the constraint on f R0 becomes tighter as we add the LSS data, the constraint on m ν remains roughly the same. This is because we are restricting our analysis to linear scales, while the effect of massive neutrinos becomes more relevant on smaller scales and, hence, causes only a small degradation of f R0 constraints.
Up to this point, we kept the unphysical lensing amplitude parameter A L fixed at its expected value of 1. However, one may wonder if the discrepancy in A L observed in the ΛCDM model also persists in f (R), and what effect co-varying A L has on the bounds on f R0 . The results for two different combinations of data are shown in Fig. 4. Although it seems that, in the case of PLC+BAO, the lensing amplitude tension has been reconciled, we argue that this is not due to a genuine signal of modified gravity. As discussed previously, the PLC+BAO data yields a peak in the likelihood of f R0 because the preference for enhanced lensing and the lack of power at low in C T T can be reconciled with a nonzero f R0 . The enhanced lensing appears to cure the A L problem and this is depicted in Fig. 4, where we see that there is a strong degeneracy between A L and f R0 for large values of the latter (blue contours). However such large values of f R0 are ruled out once we add the datasets that probe clustrering (green contours). Still, the value of A L when co-fit with f R0 is in better disagreement with the prediction. For the combination of all data we find A L = 1.11 +0.20 −0.14 68 % C.L., all datasets.
The results of the analysis with varying A L are summarized in Table II.

C. Constraints on the Symmetron model
In this Subsection we derive constraints on the inverse mass parameter, ξ , defined in Eq. (31), which represents the Compton wavelength of the scalar interaction. We fix the other two Symmetron parameters, taking a = 0.25 and β = 1, since current data is unable to constrain them simultaneously with ξ . solid line) as well as after marginalizing over a varying m ν (blue dashed line). We find an upper bound of ξ < 1.5 × 10 −3 at 95 % C.L, which corresponds to a Compton wavelength of ∼ a few Mpc. Our bounds are summarized in Table III. As mentioned above, current data is unable to simultaneously constrain all the model parameters because they are highly correlated. We also note that one cannot derive meaningful constraints for smaller values of coupling constant β as the modification of growth is relatively small for the scales and redshifts currently probed. Further, since a sets the onset of modified growth, we would see tighter constraints on ξ for smaller a values. Nevertheless, as we will show in Sec. IV, future surveys with larger sky and deeper redshift coverage will be able to constrain ξ along with the other two parameters. the inverse mass parameter ξ 0 defined by Eq. (25), and fix β 0 to a constant. Fig. 6 shows the posterior distribution for ξ 0 with the current value of the coupling parameter fixed at β 0 = 5. We find an upper bound of ξ 0 < 3 × 10 −3 (95 % C.L.). As for symmetrons, the sensitivity to the coupling is weak due to the lack of data on linear scales. However, as we will see in Sec. IV, constraints will improve significantly with future surveys. Our results for the Dilaton model are summarized in Table III IV. FORECASTS Constraints on scalar gravitational interactions derived in the previous Section, using current information available on linear scales, are relatively weak when compared to bounds available from astrophysical tests. With improved redshift resolution, depth and sky coverage that future surveys will provide, the number of modes in the linear regime will dramatically increase. Thus, it is interesting to know if future constraints from linear scales can become compatible with astrophysical bounds.

D. Constraints on the Dilaton model
In what follows, we perform a series of Fisher forecasts for the model parameters described in the previous Section, using, where possible, the current bounds on model parameters as fiducial values in the forecast. Where there was no upper bound, we use fiducial values motivated by a combination of theoretical considerations and existing constraints from non-linear scales. We also perform a principal component analysis (PCA) of m(a) for a fixed order unity coupling β, to see how well future datasets can constrain an evolving mass parameter.

A. The data assumed in the forecast
The data we consider in our forecast include CMB temperature anisotropy (T) and polarization (E) power spectra with characteristics of the Planck survey, weak lensing shear (WL) and galaxy number count (GC) from an LSST-like survey [50], with the survey parameters adopted from [73], and their cross-correlations. In some cases, we compare this to constraints expected from the Dark Energy Survey (DES) [74].
Theoretical power spectra are calculated assuming the LSST (DES) GC data is partitioned into 10 (4) tomographic redshift bins, while the WL shear field is split into 6 (4) tomographic redshift bins. In addition, we assume a flat FRW geometry and vary h, Ω c h 2 , Ω b h 2 , τ , n s , w and A s , together with the modified gravity parameters. The fiducial values of the cosmological parameters are taken to be the Planck 2015 best fit results. To calculate the WL and GC auto-and cross-correlation spectra in our scalar-tensor models, we have applied the MG-CAMB patch to CAMBSources [75]. The details of the implementation are described in [25,35].

B. Fisher analysis
For a given model, one can calculate the Fisher matrix [76] to determine how well future surveys can constrain its parameters. The inverse of the Fisher matrix provides a lower bound on the covariance matrix of the model parameters via the Cramér-Rao inequality, C ≥ F −1 . For zero-mean Gaussian-distributed observables, such as the angular correlations C XY , the Fisher matrix is given by where p a is the a th parameter of our model andC is the "observed" covariance matrix with elementsC XY that include contributions from noise: Eq. (34) assumes that all fields X(n) are measured over contiguous regions covering a fraction f sky of the sky. The value of the lowest multipole can be approximately inferred from min ≈ π/(2f sky ). The noise matrix N XY includes the statistical noise as well as the expected systematic errors. We refer the reader to [25,35] for the details of the Fisher matrix calculations for the individual experiments considered in our analysis. C. The f(R) forecast In Fig. 7 we show 1σ constraints on parameters of the Hu-Sawicki f (R) model, as expected from LSST+ (LSST WL + LSST GC + Planck CMB). Recall that current data is unable to constrain f R0 unless one assumes a fixed value for n, since the two parameters are highly degenerate. Thus, the forecast in Fig. 7 depends strongly on the assumed fiducial value, indicated with a on the plot. What we see is that for n ∼ 1 or smaller, future data will be able to constrain both parameters simultaneously. Fig. 7 also shows the importance of including the crosscorrelation between WL and GC. The information from GC alone is largely diluted by the unknown galaxy bias. Weak lensing, while not sensitive to the bias, is plagued by degeneracies coming from projection effects. Combining them helps determine the bias and break the degeneracies coming from projections. Fig. 8 compares joint 1σ constraints on f R0 and the combined mass of neutrinos, m ν , as expected from LSST+ vs those expected from DES+. We see that LSST+ can reduce uncertainties in both parameters by a factor of 3. The plot shows the effect of marginalizing over n, however the outcome depends on the assumed fiducial value of n (which is n = 1). Fig. 9 shows the bounds on the parameters of the Symmetron parameters expected from LSST+. As a fiducial model, we assume β = 1 and a mass scale of ξ = 10 −3 , which corresponds to a range of a few Mpc. Current data is unable to constrain ξ if a = 0.5 or larger. For this reason, the bound on ξ in Sec. III was derived for a  Table IV for a quantitative comparison.

D. The Symmetron forecast
fixed a = 0.25. We perform a forecast using two different fiducial values: a = 0.25 and 0.5. In the former case, LSST+ clearly improves on the bound in Sec. III, even after marginalizing over a and β . It will also be able to provide a non-trivial bound on ξ for a = 0.5, which is the value assumed in much of the previous literature. The current and expected bounds are summarized in Table IV.
It is interesting to examine the possible degeneracy be- tween the Symmetron parameters and the total mass of neutrinos. Fig. 10 shows the joint uncertainties in ξ and m ν expected from LSST+ assuming a fiducial model with β = 1, ξ = 10 −3 , a = 0.5 and m ν = 0.06eV. It is clear from the figure that there is practically no degeneracy between ξ and m ν which is because they affect the growth on different scales. Fixing the other MG parameters in this model, as opposed to marginalizing over them, does not change the degree of degeneracy, neither it improves the constraints. Fig. 11 shows expected bounds on the Dilaton model parameters, with β 0 = 1 and ξ 0 = 10 −3 as the fiducial values. Similar to the Symmetron case, we find that an LSST-like survey can constrain the inverse mass parameter ξ 0 to a percent level accuracy which is a significant improvement over current constraints. Constraints on the coupling constant β 0 , however, are not as tight as those on β in the Symmetron case. This is due to a lesser impact of the Dilaton on the linear matter power spectrum. One can see from Fig. 1 that for the chosen fiducial values, P (k) would deviate from the LCDM prediction far less in the Dilaton case compared to the Symmetron. The bottom panel in Fig. 11 shows the expected joint constraints on the neutrino masses, which are tighter than those for the Symmetron. Again, this is because Dilatons have a much lesser impact on the growth on linear scales. scalar gravitational interactions with a next generation WL survey such as LSST. In Fig. 12 we show forecasted uncertainties on the parameters of the generalized Chameleon model for two fiducial values of r, assuming that the coupling is constant (s = 0). For a slower evolution with time (r = 1), the scalaron mass decreases slower and modification to growth extends back to larger redshifts, leading to significantly tighter constraints. Thus, while LSST+ can constrain the coupling, the mass and the time-variation of the scalaron mass simultaneously, the strength of the bounds depends strongly on the assumed fiducial model.

G. Principal Component Analysis of m(a)
In addition to considering particular functional forms of β(a) and m(a) as motivated by the scalar-tensor models mentioned in preceding sections, it is also interesting to treat the coupling and the mass as two general functions and ask what features of these two functions can be constrained by the future data. In principle, one could discretize the functions β(a) and m(a) into N bins in a and treat the bins values as free parameters. However, we find that even future data will not be able to simultaneously constrain m(a) and β(a) in a completely modelindependent way, since the two parameters are largely degenerate in their effect on the observables on linear scales, as they appear together in (a, k) (see Eq. (13)). For this reason, we fix β at a constant value of order unity and discretize m(a) into bins with m(a i ), i = 1, ..., N .
As with earlier forecasts, we can calculate the Fisher matrix, and invert it to find the covariance matrix, wherep i are the "fiducial" values, and parameters include the bins m(a i ), as well as the rest of cosmological parameters. We then isolate the N ×N block of the matrix, C m corresponding to the covariance of m(a i ) after marginalization over other parameters. Since the individual bins of m(a i ) bins are highly correlated, the covariance matrix for these parameters will be non-diagonal, and the value of m in any particular bin will be practically unconstrained. The Principal Component Analysis (PCA) [27,35,36,77,78] is a way to decorrelate the parameters and find their linear combinations that are best constrained by data. Namely, we solve an eigenvalue problem to find a matrix W m that diagonalizes C m : where W m ij ≡ê i (a j ) are the eigenvectors (or eigenmodes) and λ i 's are the eigenvalues. In the limit of large N , one can write an arbitrary m(a) as an expansion intoê i (a): in which case λ i can be interpreted as the variance of α i , It is customary to order the eigenmodes from the best constrained to the worst. Then the i th eigenmode is referred to as the i th principal component (PC). Typically, one finds that only the first few modes are well constrained by the data, while most of them are practically unconstrained. For our forecast, we partition m(a) into 11 bins, with 10 of them evenly spaced in redshift within z ∈ [0, 3], and the 11th bin ranging from z = 3 to z = 30. The last bin can be taken to be wide because the observables we work with are weakly sensitive to modifications at high redshifts. In what follows, we marginalize over the 11th bin, since it is largely degenerate with some of the cosmological parameters, most prominently with Ω m . We take the fiducial model to be β = 0.4 and m(a i ) = m 0 for all i, with m 0 = H 0 /ξc and ξ = 10 −3 , corresponding to m 0 = 0.2 h/Mpc.
The left panel in Fig. 13 shows the the forecasted uncertainties in the measurement of the eigenvectors for two cases: when β is fixed, and when β is marginalized over. In both cases, we marginalize over all cosmological parameters and the 11th m-bin. The right panel in Fig. 13 presents the first four best constrained eigenvectors of m(a) after marginalizing over β. One can interpret the best constrained mode (PC1) as that corresponding to a weighted average value of m(a). The second best constrained mode (PC2) has a single node and corresponds to the difference between the high-z and low-z values of m(a). The third best mode (PC3) has 3 nodes, PC4 has 4 nodes, and so on.
The eigenvalue plot demonstrates that marginalizing over β affects the first eigenmode of m(a), but not the others. This is because the main effect of a constant β is an overall rescaling of the strength of the 5th force. It is largely degenerate with the average value of m(a), but has no impact on the detectability of time-variation of m(a). After marginalizing over β, LSST+ can measure one mass parameters (the average m(a)) with an accuracy that is better than 0.01 h/Mpc, or about 5% of the fiducial m 0 , and another 3 parameters, describing more rapid evolution of the mass with time, with accuracy better than 0.02 h/Mpc, or 10% of the fiducial value.
The extrema of the eigenmodes indicate the "sweet spots" in redshift, or epochs at which variations in m(a) are best constrained with LSST+. It is evident from the right panel in Fig. 13, for instance the shape of PC2, that LSST+ is more sensitive to time-variations at z > 1.5. This is because at higher redshifts there is a larger number of Fourier modes that are still in the linear regime.

V. DISCUSSION AND CONCLUSIONS
Modifications of gravity on cosmological scales can potentially explain the origin of cosmic acceleration. The Generalized Brans-Dicke theory, in which there is an additional scalar degree of freedom that mediates a fifth force, is one of the viable MG models that are able to fit observations after the required tuning of model parameters.
In this work, we have investigated the observational constraints on three MG models within the general framework of the GBD theory, namely, the f (R), the Symmetron and the Dilaton models, using latest observations of CMB, BAO, weak lensing and galaxy clustering. In all cases, we used observables on linear scales to avoid the complexities of the modelling of nonlinearities and redshift-space distortions.
We find that the ΛCDM model is consistent with all observations. Specifically, we find the constraint on f R0 , the model parameter in the Hu-Sawicki f (R) model, to be f R0 < 8 × 10 −5 (95% CL) when the sum of neutrino masses is fixed to be 0.06 eV. Since both massive neutrinos and MG models studied in this paper can alter the structure growth in a scale-dependent way, a degeneracy is expected. Therefore we perform another analysis with the neutrino mass varying, and we find that the constraint is diluted to f R0 < 1.0 × 10 −4 (95% CL). For the Symmetron model, the 95% CL upper limit is ξ < 1.8 × 10 −3 with β and a fixed at 1 and 0.25, respectively. For the Dilaton model, we find ξ 0 < 3 × 10 −3 at 95% CL when β 0 = 5. Tables II and III summarize the  current bounds. We have also performed a forecast for ongoing and upcoming imaging surveys including DES and LSST, and present the results in Sec. IV. A comparison between the current and future constraints on model parameters is shown in Table IV. As one can see, the improvement is significant and, despite the high level of degeneracy, more than one parameter can be constrained simultaneously. In the Hu-Sawicki model, the upper limit of f R0 is reduced by a factor of 6.7 and n can be constrained with ≈ 25% accuracy for n = 1. For the Dilaton model, current data is unable to constrain ξ 0 if β 0 = 1. However, we find that LSST+ can simultaneously constrain ξ 0 at ∼ few × 10 −5 and measure β 0 ∼ 1 with ≈ 20% accuracy. In the Symmetron model, the constraint on ξ is improved by a factor of 3, while simultaneously constraining a and β within a few percent. This is compatible with current bounds derived from astrophysical tests, such as the cluster profile [79], galactic dynamics and so on, which requires high-resolution N-body or hydrodynamical simulations [80,81] of the MG models. Additionally, to demonstrate the capabilities of an LSST-like survey, we have presented constraints on the Generalized Chameleon model in Fig. 12.
Given the power of future surveys, a modelindependent analysis will become possible. In this work, we performed a PCA study of m(a), to forecast the maximum number of parameters of the scalaron mass function that can be well determined. We find that an LSST-like survey will be able to measure the average mass parameter with an accuracy of 0.01 h/Mpc and another 3 parameters quantifying the time-variation of m(a) with an accuracy that is better than 0.02 h/Mpc. Finally, we note that future spectroscopic and HI surveys, such as eBOSS and SKA [82,83], will also provide powerful constraints on MG parameters that will be highly complementary to those from a photometric survey like LSST [84].