KiDS+GAMA: Constraints on Horndeski gravity from combined large-scale structure probes

TL;DR

This study uses combined large-scale structure data from KiDS and GAMA to constrain Horndeski gravity models, finding results compatible with Lambda-CDM and showing improved agreement with Planck measurements.

Contribution

First comprehensive joint analysis constraining Horndeski gravity using multiple probes, including cosmic shear, galaxy-galaxy lensing, and galaxy clustering, with publicly available analysis code.

Findings

Constraints on Horndeski functions consistent with Lambda-CDM

Improved S8 agreement with Planck in Horndeski models

Robustness of results to non-linear matter power spectrum modeling

Abstract

We present constraints on Horndeski gravity from a combined analysis of cosmic shear, galaxy-galaxy lensing and galaxy clustering from $450\,\mathrm{deg}^2$ of the Kilo-Degree Survey (KiDS) and the Galaxy And Mass Assembly (GAMA) survey. The Horndeski class of dark energy/modified gravity models includes the majority of universally coupled extensions to $\Lambda$CDM with one scalar field in addition to the metric. We study the functions of time that fully describe the evolution of linear perturbations in Horndeski gravity. Our results are compatible throughout with a $\Lambda$CDM model. By imposing gravitational wave constraints, we fix the tensor speed excess to zero and consider a subset of models including e.g. quintessence and $f(R)$ theories. Assuming proportionality of the Horndeski functions $\alpha_B$ and $\alpha_M$ (kinetic braiding and the Planck mass run rate, respectively) to the dark energy density fraction $\Omega_{\mathrm{DE}}(a) = 1 - \Omega_{\mathrm{m}}(a)$, we find for the proportionality coefficients $\hat{\alpha}_B = 0.20_{-0.33}^{+0.20} \,$ and $\, \hat{\alpha}_M = 0.25_{-0.29}^{+0.19}$. Our value of $S_8 \equiv \sigma_8 \sqrt{\Omega_{\mathrm{m}}/0.3}$ is in better agreement with the $Planck$ estimate when measured in the enlarged Horndeski parameter space than in a pure $\Lambda$CDM scenario. In our joint three-probe analysis we report a downward shift of the $S_8$ best fit value from the $Planck$ measurement of $\Delta S_8 = 0.016_{-0.046}^{+0.048}$ in Horndeski gravity, compared to $\Delta S_8 = 0.059_{-0.039}^{+0.040}$ in $\Lambda$CDM. Our constraints are robust to the modelling uncertainty of the non-linear matter power spectrum in Horndeski gravity. Our likelihood code for multi-probe analysis in both $\Lambda$CDM and Horndeski gravity is publicly available at http://github.com/alessiospuriomancini/KiDSHorndeski .

Figures

Tables

Table 1. Table 1 : Priors on the fit parameters. Rows 1–5 contain the priors on cosmological parameters, rows 6–9 the priors on astrophysical ‘nuisance’ parameters, rows 10–13 the priors on modified gravity parameters. All prior distributions are uniform within their ranges.

Parameter	Description	Prior range
${\omega }_{\mathrm{cdm}}$	Cold dark matter density	$[0.01,0.99]$
${\omega }_b$	Baryon density	$[0.019,0.026]$
$\ln \left({10}^{10}{A}_s\right)$	Scalar spectrum amplitude	$[1.7,5.0]$
$n_s$	Scalar spectral index	$[0.7,1.3]$
$h$	Dimensionless Hubble parameter	$[0.64,0.82]$
$A_{\mathrm{IA}}$	Intrinsic alignment amplitude	$[-6,6]$
$c_{\min }$	Baryonic feedback amplitude	$[2,4]$
${\mathrm{bias}}_{\mathrm{z1}}$	Galaxy bias of low-$z$ lens sample	$[0.1,5]$
${\mathrm{bias}}_{\mathrm{z2}}$	Galaxy bias of high-$z$ lens sample	$[0.1,5]$
${\hat{\alpha }}_B$	Prop. coeff. ${\alpha }_B\left(\tau \right)={\hat{\alpha }}_B{\Omega }_{\mathrm{DE}}\left(\tau \right)$	$[-2,2]$
${\hat{\alpha }}_M$	Prop. coeff. ${\alpha }_M\left(\tau \right)={\hat{\alpha }}_M{\Omega }_{\mathrm{DE}}\left(\tau \right)$	$[-2,2]$
${\Omega }_0$	Prop. coeff. $\Omega \left(\tau \right)={\Omega }_{0}a\left(\tau \right)$	$[0,1]$
log $k_s$	screening scale	$[-1,1]$

Table 2. Table 2 : Mean and marginalised 68 per cent credibility interval on the parameters listed obtained in Horndeski gravity with all different combinations of probes. The parameterization for the α 𝛼 \alpha functions is given by proportionality to Ω DE subscript Ω DE \Omega_{\mathrm{DE}} so that the free parameters become the proportionality coefficients (in this case α ^ B subscript ^ 𝛼 B \hat{\alpha}_{\mathrm{B}} and α ^ M subscript ^ 𝛼 M \hat{\alpha}_{\mathrm{M}} , since we fix α T = α K = 0 subscript 𝛼 T subscript 𝛼 K 0 \alpha_{\mathrm{T}}=\alpha_{\mathrm{K}}=0 ). The nuisance parameters (discussed in Secs.  4.3 ,  4.4 and 4.5 ) are the HMcode parameter c min subscript 𝑐 min c_{\mathrm{min}} , the effective bias values bias z ​ 1 subscript bias 𝑧 1 \mathrm{bias}_{z1} and bias z ​ 2 subscript bias 𝑧 2 \mathrm{bias}_{z2} for the low- z 𝑧 z and high- z 𝑧 z foreground samples, respectively, the intrinsic alignment amplitude A IA subscript 𝐴 IA A_{\mathrm{IA}} and the screening scale k s subscript 𝑘 𝑠 k_{s} . See Sec.  5.2 for a discussion of the χ 2 superscript 𝜒 2 \chi^{2} and number of degrees of freedom (d.o.f.) when Planck data are considered in the analysis.

Parameter	$P^E$	$P^E+P^{\mathrm{gm}}+P^{\mathrm{gg}}$	$P^E+P^{\mathrm{gm}}$	$P^E+P^{\mathrm{gg}}$	$P^{\mathrm{gm}}+P^{\mathrm{gg}}$	Planck	$P^E+P^{\mathrm{gm}}+P^{\mathrm{gg}}+$ Planck
${\omega }_{\mathrm{cdm}}$	${0.119}_{-0.025}^{+0.015}$	${0.121}_{-0.018}^{+0.016}$	${0.118}_{-0.020}^{+0.015}$	${0.119}_{-0.017}^{+0.013}$	${0.131}_{-0.022}^{+0.017}$	${0.117}_{-0.002}^{0.002}$	${0.117}_{-0.002}^{+0.002}$
${\mathrm{ln10}}^{10}{A}_s$	${3.01}_{-0.25}^{+0.30}$	${3.02}_{-0.19}^{+0.24}$	${3.10}_{-0.26}^{+0.24}$	${3.01}_{-0.18}^{+0.29}$	${3.06}_{-0.24}^{+0.25}$	${3.12}_{-0.01}^{0.01}$	${3.12}_{-0.01}^{+0.01}$
${\omega }_b$	${0.0225}_{-0.0037}^{+0.0035}$	${0.0224}_{-0.0037}^{+0.0038}$	${0.0221}_{-0.0033}^{+0.0034}$	${0.0223}_{-0.0036}^{+0.0038}$	${0.0226}_{-0.0038}^{+0.0037}$	${0.0225}_{0.0003}^{0.0003}$	${0.0225}_{-0.0002}^{+0.0002}$
$n_s$	${1.07}_{-0.10}^{+0.23}$	${1.13}_{-0.05}^{+0.17}$	${1.13}_{-0.07}^{+0.17}$	${1.14}_{-0.05}^{+0.16}$	${1.08}_{-0.09}^{+0.22}$	${0.97}_{-0.01}^{+0.01}$	${0.97}_{-0.01}^{+0.01}$
$h$	${0.73}_{-0.09}^{+0.09}$	${0.75}_{-0.03}^{+0.07}$	${0.77}_{-0.03}^{+0.05}$	${0.77}_{-0.02}^{+0.05}$	${0.75}_{-0.03}^{+0.07}$	${0.69}_{-0.01}^{+0.01}$	${0.69}_{-0.01}^{+0.01}$
$c_{\min }$	${3.20}_{-0.46}^{+0.80}$	${3.34}_{-0.37}^{+0.66}$	${3.28}_{-0.31}^{+0.72}$	${3.27}_{-0.29}^{+0.73}$	${3.16}_{-0.48}^{+0.84}$	-	${3.69}_{-0.08}^{+0.31}$
$A_{\mathrm{IA}}$	${0.44}_{-0.68}^{+1.16}$	${1.37}_{-0.44}^{+0.48}$	${1.39}_{-0.46}^{+0.51}$	${0.84}_{-0.61}^{+0.88}$	${1.51}_{-0.54}^{+0.51}$	-	${1.40}_{-0.35}^{+0.38}$
${\mathrm{bias}}_{z1}$	-	${1.01}_{-0.08}^{+0.06}$	${0.99}_{-0.09}^{+0.09}$	${1.01}_{-0.08}^{+0.05}$	${1.01}_{-0.09}^{+0.06}$	-	${1.02}_{-0.06}^{+0.04}$
${\mathrm{bias}}_{z2}$	-	${1.05}_{-0.09}^{+0.05}$	${1.06}_{-0.09}^{+0.06}$	${1.06}_{-0.10}^{+0.05}$	${1.05}_{-0.09}^{+0.06}$	-	${1.05}_{-0.06}^{+0.03}$
${\sigma }_8$	${0.86}_{-0.08}^{+0.07}$	${0.90}_{-0.05}^{+0.06}$	${0.93}_{-0.06}^{+0.07}$	${0.90}_{-0.05}^{+0.06}$	${0.95}_{-0.07}^{+0.07}$	${0.87}_{-0.02}^{+0.01}$	${0.86}_{-0.01}^{+0.01}$
${\Omega }_m$	${0.27}_{-0.04}^{+0.03}$	${0.26}_{-0.02}^{+0.02}$	${0.24}_{-0.03}^{+0.02}$	${0.24}_{-0.03}^{+0.02}$	${0.28}_{-0.03}^{+0.02}$	${0.29}_{-0.01}^{+0.01}$	${0.29}_{-0.01}^{+0.01}$
$S_8$	${0.803}_{-0.051}^{+0.054}$	${0.835}_{-0.039}^{+0.041}$	${0.830}_{-0.042}^{+0.053}$	${0.806}_{-0.044}^{+0.047}$	${0.905}_{-0.060}^{+0.057}$	${0.859}_{-0.032}^{+0.029}$	${0.843}_{-0.024}^{+0.021}$
${\hat{\alpha }}_B$	${0.05}_{-0.36}^{+0.30}$	${0.20}_{-0.33}^{+0.20}$	${0.26}_{-0.44}^{+0.33}$	${0.23}_{-0.33}^{+0.20}$	${0.28}_{-0.45}^{+0.22}$	${0.79}_{-0.71}^{+0.42}$	${0.36}_{-0.22}^{+0.18}$
${\hat{\alpha }}_M$	${0.46}_{-0.61}^{+0.25}$	${0.25}_{-0.29}^{+0.19}$	${0.31}_{-0.44}^{+0.26}$	${0.23}_{-0.30}^{+0.20}$	${0.36}_{-0.44}^{+0.24}$	${0.18}_{-0.52}^{+0.35}$	${0.15}_{-0.31}^{+0.13}$
$\log k_s$	${0.00}_{-1.00}^{+1.00}$	$-{0.02}_{-0.98}^{+0.74}$	$-{0.01}_{-0.99}^{+1.01}$	${0.01}_{-1.01}^{+0.99}$	${0.07}_{-1.06}^{+0.87}$	-	${0.11}_{-1.11}^{+0.89}$
${\chi }^2$	63.2	115.7	100.1	67.9	42.8	10697	10821
d.o.f.	40	88	78	48	38

Table 3. Table 3 : Mean and marginalised 68 credibility interval on the parameters listed obtained in Horndeski gravity with all different combinations of probes. The parameterization for the α 𝛼 \alpha functions is the planck_linear parameterization implemented in HiClass , where the α 𝛼 \alpha functions depend on only one function of time Ω ​ ( τ ) Ω 𝜏 \Omega(\tau) (see equations  5 ), which is assumed to be proportional to the scale factor, Ω ​ ( τ ) = Ω 0 ​ a ​ ( τ ) Ω 𝜏 subscript Ω 0 𝑎 𝜏 \Omega(\tau)=\Omega_{0}a(\tau) . The only free parameter as far as the α 𝛼 \alpha functions are concerned becomes then Ω 0 subscript Ω 0 \Omega_{0} . The nuisance parameters are the same as in Table 2 .

Parameter	$P^E$	$P^E+P^{\mathrm{gm}}+P^{\mathrm{gg}}$	$P^E+P^{\mathrm{gm}}$	$P^E+P^{\mathrm{gg}}$	$P^{\mathrm{gm}}+P^{\mathrm{gg}}$
${\omega }_{\mathrm{cdm}}$	${0.163}_{-0.052}^{+0.037}$	${0.138}_{-0.026}^{+0.017}$	${0.135}_{-0.030}^{+0.021}$	${0.136}_{-0.026}^{+0.021}$	${0.151}_{-0.033}^{+0.023}$
${\mathrm{ln10}}^{10}{A}_s$	${3.00}_{-0.60}^{+0.56}$	${3.06}_{-0.33}^{+0.34}$	${3.21}_{-0.43}^{+0.36}$	${3.06}_{-0.41}^{+0.35}$	${2.93}_{-0.32}^{+0.42}$
${\omega }_b$	${0.0225}_{-0.0037}^{+0.0038}$	${0.0225}_{-0.0035}^{+0.0038}$	${0.0224}_{-0.0036}^{+0.0038}$	${0.0223}_{-0.0036}^{+0.0035}$	${0.0224}_{-0.0037}^{+0.0038}$
$n_s$	${1.03}_{-0.08}^{+0.27}$	${1.10}_{-0.05}^{+0.20}$	${1.05}_{-0.10}^{+0.25}$	${1.09}_{-0.06}^{+0.21}$	${1.05}_{-0.10}^{+0.25}$
$h$	${0.73}_{-0.09}^{+0.08}$	${0.76}_{-0.02}^{+0.07}$	${0.76}_{-0.02}^{+0.06}$	${0.77}_{-0.02}^{+0.05}$	${0.76}_{-0.02}^{+0.06}$
$c_{\min }$	${3.07}_{-0.35}^{+0.93}$	${3.10}_{-0.38}^{+0.89}$	${3.21}_{-0.26}^{+0.79}$	${3.12}_{-0.47}^{+0.88}$	${3.21}_{-0.26}^{+0.79}$
$A_{\mathrm{IA}}$	${0.63}_{-0.60}^{1.03}$	${1.31}_{-0.39}^{+0.43}$	${1.35}_{-0.41}^{+0.45}$	${0.70}_{-0.61}^{+0.90}$	${1.35}_{-0.49}^{+0.43}$
${\mathrm{bias}}_{z1}$	-	${1.03}_{-0.11}^{+0.07}$	${0.97}_{-0.15}^{+0.14}$	${1.04}_{-0.11}^{+0.07}$	${1.07}_{-0.16}^{+0.09}$
${\mathrm{bias}}_{z2}$	-	${1.10}_{-0.13}^{+0.07}$	${1.15}_{-0.17}^{+0.09}$	${1.13}_{-0.15}^{+0.08}$	${1.10}_{-0.16}^{+0.10}$
${\sigma }_8$	${0.87}_{-0.14}^{+0.11}$	${0.93}_{-0.09}^{+0.08}$	${0.93}_{-0.10}^{+0.08}$	${0.92}_{-0.10}^{+0.08}$	${0.85}_{-0.07}^{+0.09}$
${\Omega }_m$	${0.35}_{-0.09}^{+0.07}$	${0.28}_{-0.03}^{+0.02}$	${0.27}_{-0.04}^{+0.03}$	${0.27}_{-0.04}^{+0.03}$	${0.31}_{-0.05}^{+0.03}$
$S_8$	${0.917}_{-0.106}^{+0.096}$	${0.895}_{-0.074}^{+0.048}$	${0.879}_{-0.066}^{+0.047}$	${0.866}_{-0.077}^{+0.051}$	${0.863}_{-0.043}^{+0.047}$
${\Omega }_0$	${0.54}_{-0.29}^{+0.46}$	${0.21}_{-0.21}^{+0.04}$	${0.31}_{-0.31}^{+0.07}$	${0.23}_{-0.23}^{+0.04}$	${0.44}_{-0.44}^{+0.56}$
log $k_s$	${0.18}_{-0.37}^{+0.48}$	${0.05}_{-1.05}^{+0.57}$	$-{0.24}_{-0.76}^{+0.28}$	${0.01}_{-1.01}^{+0.52}$	${0.06}_{-1.07}^{+0.86}$
${\chi }^2$	62.1	116.2	99.3	68.8	43.5
d.o.f.	41	89	79	49	39

Table 4. Table 4 : Mean and marginalised 68 per cent credibility interval on the parameters listed, obtained with our new likelihood in a Λ Λ \Lambda CDM scenario with the priors specified in Tab.  1 .

Param	$P^E$		$P^E+P^{\mathrm{gm}}+P^{\mathrm{gg}}$		$P^E+P^{\mathrm{gm}}$		$P^E+P^{\mathrm{gg}}$		$P^{\mathrm{gm}}+P^{\mathrm{gg}}$
	This work	vU18	This work	vU18	This work	vU18	This work	vU18	This work	vU18
${\omega }_{\mathrm{cdm}}$	${0.139}_{-0.041}^{+0.051}$	${0.134}_{-0.046}^{+0.044}$	${0.165}_{-0.045}^{+0.033}$	${0.153}_{-0.040}^{+0.030}$	${0.142}_{-0.043}^{+0.032}$	${0.132}_{-0.044}^{+0.032}$	${0.177}_{-0.038}^{+0.035}$	${0.169}_{-0.039}^{+0.034}$	${0.175}_{-0.051}^{+0.034}$	${0.159}_{-0.041}^{+0.035}$
${\mathrm{ln10}}^{10}{A}_s$	${2.65}_{-0.92}^{+0.21}$	${2.76}_{-1.06}^{+0.30}$	${2.47}_{-0.48}^{+0.44}$	${2.63}_{-0.44}^{+0.36}$	${2.63}_{-0.72}^{+0.43}$	${2.81}_{-0.78}^{+0.52}$	${2.22}_{-0.49}^{+0.14}$	${2.30}_{-0.53}^{+0.21}$	${2.53}_{-0.48}^{+0.41}$	${2.67}_{-0.45}^{+0.37}$
${\omega }_b$	${0.0225}_{-0.0033}^{+0.0030}$	${0.0224}_{-0.0034}^{+0.0036}$	${0.0225}_{-0.0033}^{+0.0030}$	${0.0225}_{-0.0035}^{+0.0035}$	${0.0225}_{-0.0031}^{+0.0036}$	${0.0224}_{-0.0034}^{+0.0036}$	${0.0225}_{-0.0029}^{+0.0034}$	${0.0225}_{-0.0035}^{+0.0034}$	${0.0226}_{-0.0029}^{+0.0035}$	${0.0225}_{-0.0033}^{+0.0034}$
$n_s$	${1.09}_{-0.08}^{+0.21}$	${1.11}_{-0.05}^{+0.19}$	${0.95}_{-0.21}^{+0.10}$	${0.97}_{-0.19}^{+0.15}$	${1.03}_{-0.15}^{+0.20}$	${1.08}_{-0.07}^{+0.22}$	${0.93}_{-0.19}^{+0.11}$	${0.97}_{-0.18}^{+0.14}$	${0.89}_{-0.18}^{+0.09}$	${0.93}_{-0.22}^{+0.08}$
$h$	${0.74}_{-0.06}^{+0.07}$	${0.74}_{-0.04}^{+0.08}$	${0.73}_{-0.07}^{+0.08}$	${0.73}_{-0.06}^{+0.09}$	${0.74}_{-0.04}^{+0.08}$	${0.74}_{-0.03}^{+0.08}$	${0.73}_{-0.07}^{+0.07}$	${0.73}_{-0.08}^{+0.04}$	${0.74}_{-0.05}^{+0.08}$	${0.73}_{-0.08}^{+0.09}$
$c_{\min }$	${3.24}_{-0.54}^{+0.74}$	${3.27}_{-0.23}^{+0.73}$	${2.87}_{-0.79}^{+0.34}$	${2.97}_{-0.71}^{+0.56}$	${3.19}_{-0.36}^{+0.79}$	${3.28}_{-0.22}^{+0.72}$	${2.98}_{-0.66}^{+0.69}$	${3.08}_{-0.39}^{+0.81}$	${2.66}_{-0.65}^{+0.25}$	${2.86}_{-0.84}^{+0.30}$
$A_{IA}$	${0.91}_{-0.57}^{+0.75}$	${0.92}_{-0.59}^{+0.78}$	${1.24}_{-0.36}^{+0.37}$	${1.27}_{-0.40}^{+0.39}$	${1.36}_{-0.37}^{+0.36}$	${1.46}_{-0.42}^{+0.41}$	${0.88}_{-0.47}^{+0.65}$	${0.88}_{-0.50}^{+0.70}$	${1.38}_{-0.43}^{+0.42}$	${1.38}_{-0.49}^{+0.46}$
${\mathrm{bias}}_{z1}$	-	-	${1.17}_{-0.18}^{+0.16}$	${1.12}_{-0.15}^{+0.14}$	${0.84}_{-0.21}^{+0.15}$	${0.78}_{-0.19}^{+0.14}$	${1.25}_{-0.17}^{+0.16}$	${1.21}_{-0.15}^{+0.14}$	${1.18}_{-0.19}^{+0.18}$	${1.13}_{-0.16}^{+0.15}$
${\mathrm{bias}}_{z2}$	-	-	${1.29}_{-0.19}^{+0.18}$	${1.25}_{-0.17}^{+0.16}$	${1.55}_{-0.33}^{+0.28}$	${1.45}_{-0.33}^{+0.27}$	${1.49}_{-0.19}^{+0.20}$	${1.45}_{-0.18}^{+0.18}$	${1.27}_{-0.20}^{+0.19}$	${1.23}_{-0.17}^{+0.16}$
${\Omega }_m$	${0.30}_{-0.08}^{+0.08}$	${0.29}_{-0.10}^{+0.07}$	${0.35}_{-0.06}^{+0.06}$	${0.33}_{-0.06}^{+0.05}$	${0.30}_{-0.07}^{+0.06}$	${0.29}_{-0.08}^{+0.06}$	${0.38}_{-0.07}^{+0.06}$	${0.36}_{-0.06}^{+0.06}$	${0.36}_{-0.08}^{+0.06}$	${0.34}_{-0.06}^{+0.05}$
${\sigma }_8$	${0.77}_{-0.16}^{+0.07}$	${0.80}_{-0.18}^{+0.09}$	${0.74}_{-0.09}^{+0.06}$	${0.78}_{-0.08}^{+0.06}$	${0.77}_{-0.12}^{+0.08}$	${0.81}_{-0.14}^{+0.09}$	${0.68}_{-0.08}^{+0.05}$	${0.70}_{-0.08}^{+0.05}$	${0.77}_{-0.10}^{+0.07}$	${0.80}_{-0.09}^{+0.07}$
$S_8$	${0.760}_{-0.038}^{+0.039}$	${0.761}_{-0.038}^{+0.040}$	${0.792}_{-0.031}^{+0.032}$	${0.800}_{-0.026}^{+0.030}$	${0.756}_{-0.035}^{+0.039}$	${0.769}_{-0.032}^{+0.037}$	${0.752}_{-0.036}^{+0.036}$	${0.759}_{-0.032}^{+0.036}$	${0.840}_{-0.040}^{+0.045}$	${0.835}_{-0.037}^{+0.038}$

Equations

$$S[g_{\mu \nu },\varphi ]$$

$${\mathcal{L}}_2$$

$${\mathcal{L}}_3$$

$${\mathcal{L}}_4$$

$${\mathcal{L}}_5$$

$$-\frac{1}{6}{G}_{5X}(\varphi ,X)\left[(\square \varphi {)}^3+2{\varphi }_{;\mu }{}^{\nu }{\varphi }_{;\nu }{}^{\alpha }{\varphi }_{;\alpha }{}^{\mu }-3{\varphi }_{;\mu \nu }{\varphi }^{;\mu \nu }\square \varphi \right].$$

$$\mathrm{d}{s}^2=-\left(1+2\frac{\Phi }{c^2}\right)c^2\mathrm{d}{t}^2+a^2\left(t\right)\left(1-2\frac{\Psi }{c^2}\right)\mathrm{d}{\mathbf{x}}^2,$$

$${\alpha }_M\equiv \frac{\mathrm{d}\ln M_{\ast }^2}{\mathrm{d}\ln a}.$$

$${\alpha }_i(\tau )={\hat{\alpha }}_i{\Omega }_{\mathrm{DE}}(\tau ).$$

$$\{\begin{array}{ll}{\alpha }_{\text{K}}=& \frac{3\left({\rho }_{\text{DE}}+p_{\text{DE}}\right)}{H^2}+\frac{3\Omega \left({\rho }_{\text{m}}+p_{\text{m}}\right)}{H^2\left(1+\Omega \right)}-\frac{{\Omega }^{\prime \prime }-2aH{\Omega }^{\prime }}{a^2{H}^2\left(1+\Omega \right)}\\ {\alpha }_{\text{M}}=& \frac{{\Omega }^{\prime }}{aH\left(1+\Omega \right)}\\ {\alpha }_{\text{B}}=& -{\alpha }_{\text{M}}\\ {\alpha }_{\text{T}}=& 0,\end{array}$$

$$\Phi$$

$$\frac{\Psi }{\Phi }$$

$$P_{ij}^E(\ell )$$

$$\times \mu {\left(\frac{\ell +1/2}{\chi },\chi \right)}^2{\left(\frac{1+\eta \left(\frac{\ell +1/2}{\chi },\chi \right)}{2}\right)}^2,$$

$$g_i(\chi )$$

$$P_{ij}^{\mathrm{gm}}(\ell )$$

$$\times \mu {\left(\frac{\ell +1/2}{\chi },\chi \right)}^2\left(\frac{1+\eta \left(\frac{\ell +1/2}{\chi },\chi \right)}{2}\right),$$

$$P_{ii}^{\mathrm{gg}}(\ell )=b_i{b}_i{\int }_0^{{\chi }_H}\mathrm{d}\chi \frac{n_{Fi}(\chi )n_{Fi}(\chi )}{{\chi }^2}{P}_{\delta }\left(\frac{\ell +1/2}{\chi };\chi \right)\mu {\left(\frac{\ell +1/2}{\chi },\chi \right)}^2,$$

$$-2\ln \mathcal{L}(\mathbf{\theta })=\sum _{\alpha ,\beta }{d}_{\alpha }(\mathbf{\theta })({\mathbf{C}}^{-1}{)}_{\alpha \beta }{d}_{\beta }(\mathbf{\theta }),$$

$$d_{\alpha }(\mathbf{\theta })=P_{\alpha }-⟨P_{\alpha }(\mathbf{\theta }){⟩}^{\mathrm{model}},$$

$$c(M,z)=c_{\min }\frac{1+z_{\mathrm{f}}}{1+z},$$

$$\begin{array}{rl}\mu (k,a)\to & 1+{\mu }_{\mathrm{MG}}(k,a)\exp (-(k/k_{\mathrm{s}}{)}^2)\\ \eta (k,a)\to & 1+{\eta }_{\mathrm{MG}}(k,a)\exp (-(k/k_{\mathrm{s}}{)}^2)\\ D_{+}(k,a)\to & D_{+\mathrm{GR}}(a)+D_{+\mathrm{MG}}(k,a)\exp (-(k/k_{\mathrm{s}}{)}^2).\end{array}$$

$$P_{\delta \mathrm{I}}(k,z)=-A_{\mathrm{IA}}{C}_1{\rho }_{\mathrm{crit}}\frac{{\Omega }_{\mathrm{m}}}{D_{+}(k,z)}{P}_{\delta }(k,z)\equiv F(k,z)P_{\delta }(k,z),$$

$$P_{ij}^{GI}(\ell )$$

$$\times \mu {\left(\frac{\ell +1/2}{\chi },\chi \right)}^2\left(\frac{1+\eta \left(\frac{\ell +1/2}{\chi },\chi \right)}{2}\right)+i\leftrightarrow j$$

$$P_{ij}^{II}(\ell )$$

$$\times \mu {\left(\frac{\ell +1/2}{\chi },\chi \right)}^2,$$

$$P_{ij}^{gI}(\ell )$$

$$\times \mu {\left(\frac{\ell +1/2}{\chi },\chi \right)}^2.$$

$${\varphi }_i=\frac{1}{c^2}\int \mathrm{d}\chi \left(\Phi +\Psi \right)g_i(\chi ),$$

$$⟨{\varphi }_i{\varphi }_j⟩(\ell )=\frac{1}{c^4}\int \frac{\mathrm{d}\chi }{{\chi }^2}{g}_i(\chi )g_j(\chi )P_{\Phi +\Psi }\left(\frac{\ell +1/2}{\chi };\chi \right),$$

$$\Delta \varphi =2\kappa ,$$

$$⟨{\kappa }_i{\kappa }_j⟩(\ell )=\frac{1}{4c^4}\int \frac{\mathrm{d}\chi }{{\chi }^2}{g}_i(\chi )g_j(\chi )\frac{P_{\Phi +\Psi }}{k^4}$$

$$\frac{P_{\Phi +\Psi }\left(\frac{\ell +1/2}{\chi }\right)}{k^4}={\left(\frac{3}{2}\frac{{\Omega }_{\mathrm{m}}{H}_0^2}{c^2}\right)}^2{P}_{\delta }\left(\frac{\ell +1/2}{\chi }\right)\left[1+2\eta +{\eta }^2\right]{\mu }^2,$$

$$⟨{\kappa }_i{\kappa }_j⟩(\ell )=$$

$$\displaystyle\beta(k,z)=\log\big{[}1+k/k_{\sigma}(z)\big{]}\times\zeta_{NL},$$

Full text

††pagerange: KiDS+GAMA: Constraints on Horndeski gravity from combined large-scale structure probes–11††pubyear: 2018

1 Introduction

The investigation of the accelerated expansion of the Universe is one of the main areas of active research in cosmology. A cosmological constant $\Lambda$ is an excellent fit to most observations and can be considered responsible for cosmic acceleration. However, while representing a key component of the concordance $\Lambda$CDM model, the cosmological constant still lacks of a solid theoretical understanding in terms of naturalness and interpretation as the energy density of the vacuum (see e.g. Martin, 2012, for a review). As an alternative to the cosmological constant, a dark energy component has been invoked, as a fluid with sufficiently negative pressure to drive the accelerated expansion of the Universe. Another possible interpretation of the observed acceleration may be as the consequence of a modification of the laws of gravity on large, cosmological scales, marking deviations from General Relativity (see e.g. Clifton et al., 2012, for a review).

A large fraction of models for dark energy and modified gravity is characterized by the presence of a scalar field in the Lagrange density, as an additional gravitational degree of freedom parallel to the metric. The mathematical formulation of the most general expression for the Lagrange density of such a scalar-tensor theory, with derivatives in the equations of motion not higher than second order, was first discussed in Horndeski (1974) and subsequently rediscovered in Nicolis et al. (2009) and Deffayet et al. (2011). This so-called Horndeski Lagrangian encompasses a remarkably large number of dark energy and modified gravity models (see Section 2 for a list of some of them); the condition on the highest order of the derivatives guarantees the stability of all these theories against ghost-like degrees of freedom (Ostrogradsky’s instabilities; see Woodard, 2015).

Current and future surveys aim at using different cosmological probes to investigate the true nature of the accelerated expansion (see Weinberg et al., 2013, for a review), with weak gravitational lensing and galaxy clustering at the forefront of these studies.

The weak gravitational lensing effect is particularly sensitive to the growth of cosmic structure, encoding therefore precious information on the acceleration as a dynamical effect in redshift, or equivalently time. Particularly rich in information on dark energy and modified gravity is the weak lensing effect produced by the large-scale structure (LSS) of the Universe, or cosmic shear. This is observed in the form of correlations of galaxy shapes due to the distortion of the cross-sectional shape of light bundles coming from background sources, caused by the gravitational tidal fields of the LSS in the foreground (see e.g. Bartelmann & Schneider, 2001; Hoekstra & Jain, 2008; Kilbinger, 2015, for reviews on the topic).

Besides considering correlations of galaxy shapes, information on the weak lensing effect can be extracted correlating the positions of foreground galaxies, which trace the large-scale structure, with the shapes of background galaxies. This galaxy-matter cross-correlation is often referred to as galaxy-galaxy lensing.

The shape distortions produced by weak gravitational lensing are difficult to measure, since the induced source galaxy ellipticities are at the percent level, and a number of systematic effects can mimic this signal (Mandelbaum, 2017). Despite these difficulties, cosmological results have already been derived by numerous cosmic shear surveys (Heymans et al., 2013; Jee et al., 2013; Hildebrandt et al., 2017; Hikage et al., 2018; Troxel et al., 2018a). Interestingly, the fiducial cosmic shear analyses of the Canada-France-Hawaii Lensing Survey 111http://www.cfhtlens.org/(CFHTLenS; Heymans et al., 2013; Joudaki et al., 2017) and the Kilo-Degree Survey 222http://kids.strw.leidenuniv.nl/(KiDS; Hildebrandt et al., 2017; Köhlinger et al., 2017), both prefer a cosmological model that is in mild tension with the best-fitting parameters obtained from Cosmic Microwave Background (CMB) measurements of the Planck 333http://sci.esa.int/planck/ satellite (Planck Collaboration et al., 2016a, b, 2018).

Galaxy-galaxy lensing measurements have also matured in recent years and can help constrain dark energy (see e.g. Kwan et al., 2017), in particular when combined with studies of galaxy clustering (see e.g. the recent analysis of Singh et al., 2018): the statistic describing at lowest order the galaxy spatial distribution is the two-point correlation function, which in the past provided early evidence for the $\Lambda$CDM model (Baugh, 1996; Eisenstein & Zaldarriaga, 2001; Saunders et al., 2000; Huterer et al., 2001; Hamilton & Tegmark, 2002; Cole et al., 2005; Tegmark et al., 2006). Today, studying the spatial distribution of galaxies and its evolution in time is crucial to analyse possible extensions to the cosmological concordance model. The evolution of the clustering of galaxies with redshift can put direct constraints on models for the evolution of density perturbations, which is key to discriminate between different dark energy/modified gravity theories. However, the interpretation of galaxy clustering is complicated by galaxy bias (Kaiser, 1984), the relation between the galaxy spatial distribution and the theoretically predicted matter distribution.

Modern optical imaging surveys measure the positions and ellipticities of millions of galaxies; from them, the galaxy overdensity field as well as the gravitational lensing shear field can be derived. The two-point auto- and cross-correlations of these two fields are the two-point correlation functions of cosmic shear, galaxy-galaxy lensing and galaxy clustering. A joint analysis of these correlation functions can break degeneracies between cosmological and nuisance parameters, leading to tighter cosmological constraints (Joachimi & Bridle, 2010). Several earlier studies have indeed considered such joint analyses (Cacciato et al., 2013; Mandelbaum et al., 2013; More et al., 2015; Kwan et al., 2017; Nicola et al., 2017), albeit very few in a modified gravity context. Among the latter, Joudaki et al. (2018) recently performed a combined analysis of cosmic shear tomography, galaxy-galaxy lensing tomography, and redshift-space multipole power spectra using data from KiDS-450 ($\sim 450~{}\mathrm{deg}^{2}$ of cosmic shear data from the KiDS survey) and two overlapping spectroscopic surveys, the 2-degree Field Lensing Survey 444http://2dflens.swin.edu.au/ (2dFLenS; Blake et al., 2016) and the Baryon Oscillation Spectroscopic Survey 555http://www.sdss3.org/surveys/boss.php (BOSS; Dawson et al., 2013). They found that none of the extended cosmologies considered were simultaneously favoured in a model selection sense and able to resolve the discordance with Planck, except for an evolving dark energy component with a time-dependent $w_{0}-w_{a}$ equation of state. Amon et al. (2018) presented a measurement of $E_{G}$, a statistic combining measurements of weak gravitational lensing, galaxy clustering and redshift space distortions, proposed as a consistency test of General Relativity (Zhang et al., 2007). They determined the value of $E_{G}$ using data from the KiDS, 2dFLenS, BOSS and Galaxy And Mass Assembly 666http://www.gama-survey.org/ (GAMA; Driver et al., 2009, 2011; Liske et al., 2015) surveys; their results show that measurements of the $E_{\mathrm{G}}$ statistic cannot be conducted as consistency checks of General Relativity until the aforementioned tension in cosmological parameters is resolved, and their $E_{\mathrm{G}}$ measurements favour a lower matter density cosmology than the CMB. Recently, DES Collaboration et al. (2018) presented a combined analysis of galaxy clustering and weak gravitational lensing from the first-year data of the Dark Energy Survey, targeting modifications of the metric potentials that would be a signal of modified gravity. They found that their constraints are compatible with a cosmological constant scenario.

For a flat $\Lambda$CDM model, van Uitert et al. (2018, hereafter vU18) exploited the overlap of the KiDS and GAMA surveys to produce constraints on cosmological parameters by cross-correlating cosmic shear measurements from KiDS-450, the galaxy-matter cross-correlation from KiDS-450 around two foreground samples of GAMA galaxies, and the angular correlation function of the same foreground galaxies. Their results are consistent with both the KiDS-450 (Hildebrandt et al., 2017; Köhlinger et al., 2017) and Planck analyses (Planck Collaboration et al., 2016a, b, 2018) when considering the cross-correlation of all three probes; with cosmic shear alone, their results are fully consistent with the KiDS analysis, while showing a similar tension with the Planck results.

In this work we extend the analysis of vU18 to the Horndeski class of dark energy and modified gravity models. The evolution of linear cosmological perturbations in Horndeski gravity can be fully described by four functions of time only (Gleyzes et al., 2013; Bellini & Sawicki, 2014). Assuming a time parameterization for these functions, it is possible to set constraints on the parameters describing their time evolution: this is the main goal of our study. We consider two parameterizations and concentrate in particular on a parameterization that sets these functions proportional to the dark energy density fraction $\Omega_{\mathrm{DE}}(a)=1-\Omega_{\mathrm{m}}(a)$; this choice allows the Horndeski functions to reproduce the late-time cosmic acceleration. One of these functions is largely unconstrained by large-scale structure probes, and uncorrelated with all other cosmological parameters, therefore we fix its value. We fix also the Horndeski function describing deviations of the gravitational wave speed from that of light, as it has been recently constrained by gravitational wave detections. We set constraints on the remaining two Horndeski functions, describing how the scalar and metric kinetic terms mix and how the effective gravitational constant evolves over time, respectively.

This paper is structured as follows. In Section 2 we review the Horndeski class of dark energy/modified gravity models and introduce the description of linear perturbations underlying our constraints. In Section 3 we present expressions for the power spectra of the three probes considered, assuming a generic modified gravity scenario. In Section 4, after briefly describing the data used in our analysis and the surveys from which they were obtained, we detail our methodology for cosmological inference, in particular how we modelled the theoretical power spectra for comparison with our data vector. In Section 5 we present our Horndeski gravity constraints, while the agreement of our $\Lambda$CDM results with the ones obtained by vU18 is shown in Appendix D. We also describe how our results compare with the Planck results, and present constraints obtained from a combined analysis with CMB measurements. Finally, we draw our conclusions in Section 6.

2 Horndeski theories

The most general scalar-tensor theory of gravity that is four-dimensional, Lorentz-invariant, local, energy-momentum conserving and has equations of motion with derivatives not higher than second order, can be written as (Horndeski, 1974):

[TABLE]

where $g=\mathrm{det}\,g_{\mu\nu}$, $R$ is the Ricci scalar and the integration in $\mathrm{d}^{4}x$ is carried out over the whole four-dimensional spacetime. The $G_{i}(\phi,X)$ are arbitrary functions of the additional scalar field $\phi$ and its kinetic term $X=-\frac{1}{2}\nabla_{\mu}\phi\nabla^{\mu}\phi$. The choice of the $G_{i}(\phi,X)$ functions specifies the single modified gravity model considered within the Horndeski class. The four contributions $\mathcal{L}_{i}$ to the gravitational sector depend on the $G_{i}(\phi,X)$ functions and on their partial derivatives, denoted with subscripts $\phi,X$, e.g. $G_{iX}=\partial G_{i}/\partial X$. The Horndeski Lagrangian only considers universal coupling between the metric and the matter fields (collectively described by $\psi_{m}$ and contained in the matter Lagrangian $\mathcal{L}_{m}$), which are therefore uncoupled to the scalar field. Most of the universally coupled models with one additional scalar degree of freedom belong to the Horndeski class. These include for example quintessence (Ratra & Peebles, 1988; Wetterich, 1988), Brans-Dicke models (Brans & Dicke, 1961), $k$-essence (Armendáriz-Picón et al., 1999; Armendariz-Picon et al., 2001), kinetic gravity braiding (Deffayet et al., 2010; Kobayashi et al., 2010; Pujolàs et al., 2011), covariant galileons (Nicolis et al., 2009; Deffayet et al., 2009), disformal and Dirac-Born-Infeld gravity (de Rham & Gabadadze, 2010; Zumalacárregui et al., 2013; Bettoni & Liberati, 2013), chameleons (Khoury & Weltman, 2004a, b), symmetrons (Hinterbichler & Khoury, 2010; Hinterbichler et al., 2011), Gauss-Bonnet couplings (Ezquiaga et al., 2016) and models screening the cosmological constant (Charmousis et al., 2012; Martín-Moruno et al., 2015), all variants of $f(R)$ (Carroll et al., 2004) and $f(G)$ (Carroll et al., 2005) theories.

We now turn to the description of linear cosmological perturbations in Horndeski gravity. Considering linear scalar perturbations on a Friedmann-Robertson-Walker metric, assuming spatial flatness on large scales one can write the line element in Newtonian gauge, using cosmic time $t$ and comoving coordinates $\mathbf{x}$, as

[TABLE]

with the Bardeen potentials $\Phi$ and $\Psi$ satisfying the condition $\Phi=\Psi$ in General Relativity in absence of anisotropic stress; for a modified gravity theory, this equality in general does not hold.

We fix the background expansion to $\Lambda$CDM and study the evolution of linear cosmological perturbations in Horndeski gravity. It has been shown that linear perturbations in Horndeski theories can be parameterized by means of four functions of (conformal) time $\tau$ only (Gleyzes et al., 2013; Bellini & Sawicki, 2014), which purely affect structure formation leaving the background expansion unchanged. A possible choice for these four functions is particularly informative on the physical meaning associated to each of them. We will collectively denote these four specific functions as $\alpha(\tau)$ and briefly comment on their physical interpretation here, referring to Bellini & Sawicki (2014) and references therein for a more complete description:

• •

$\alpha_{K}$ is the kineticity term, i.e. the kinetic energy of the scalar perturbations arising directly from the action. In the quasi-static approximation where time derivatives are negligible compared to space derivatives, $\alpha_{K}$ does not enter the equations of motion and is therefore largely unconstrained by cosmic shear and other large-scale structure probes (Bellini et al., 2016; Alonso et al., 2017; Spurio Mancini et al., 2018; Reischke et al., 2019). In our analysis we will fix $\alpha_{K}$ to zero, i.e. its General Relativity value. This does not affect the constraints on the other $\alpha$ functions, since $\alpha_{K}$ is largely uncorrelated with all other parameters (Spurio Mancini et al., 2018; Reischke et al., 2019) 777In practice, when running our MCMC chains we fix $\alpha_{K}$ to a very small value; for example, when using the parameterization $\alpha_{K}(\tau)=\hat{\alpha}_{K}\Omega_{\mathrm{DE}}(\tau)$ (see the end of Section 2), we set $\hat{\alpha}_{K}=0.01$. Observationally, the difference between predicted matter power spectra (which enter the expression for the power spectra of the cosmological probes in our analysis) for $\hat{\alpha}_{K}=0.01$ and $\hat{\alpha}_{K}=0$ would not be distinguishable, while setting $\hat{\alpha}_{K}=0$ produces numerical instabilities in HiClass. For this reason we fix $\hat{\alpha}_{K}$ to 0.01; we also notice that any other value for $\hat{\alpha}_{K}$ would be equally acceptable, as $\alpha_{K}$ is in any case unconstrained by large-scale structure probes, and largely uncorrelated with all other parameters.;

• •

$\alpha_{B}$ is the braiding term, which describes mixing of the scalar field with the metric kinetic term, leading to what is typically interpreted as a fifth force between massive particles influencing the structure formation rate. As shown by Bellini & Sawicki (2014), $\alpha_{B}$ introduces a new scale dependence of the theory beyond the Jeans length. The dark energy component clusters at small scales only if there is braiding, $\alpha_{B}\neq 0$ (Bellini et al., 2016);

• •

$\alpha_{M}$ is the Planck-mass run rate, defined by

[TABLE]

$\alpha_{M}$ describes the rate of evolution of the effective Planck mass $M_{*}^{2}$, defined as the dimensionless product of the normalization of the kinetic term for gravitons and $8\pi G_{N}$ measured on Earth, i.e. it describes a time evolution of the effective gravitational constant which modifies the growth of perturbations when there is non-zero braiding. A mere fixed redefinition of the Planck mass would not produce detectable effects on structure formation, as it could always be reabsorbed by an appropriate rescaling of the action. Theories with $\alpha_{M}\neq 0$ are non-minimally coupled. As shown in Saltas et al. (2014), $\alpha_{M}$ generates anisotropic stress and modifies the propagation of gravitational waves;

• •

$\alpha_{T}$ is the tensor speed excess, indicating deviations of the propagation speed of gravitational waves from the speed of light. Recently, very strong constraints have been placed on $\alpha_{T}$ by the measurement of the gravitational wave speed derived by the detection of the binary neutron star merger GW170817 and the associated gamma ray burst GRB170817A (Abbott et al., 2017a, b; Baker et al., 2017; Creminelli & Vernizzi, 2017; Ezquiaga & Zumalacárregui, 2017; Sakstein & Jain, 2017; Lombriser & Lima, 2017; Bettoni et al., 2017). Since the tensor speed has been found to be very close to that of light, the value of $\alpha_{T}$ at present time has been constrained to be very close to zero. While in principle this constraint may not apply to $\alpha_{T}$’s value at all times (see e.g. Amendola et al., 2018; Battye et al., 2018; de Rham & Melville, 2018), setting $\alpha_{T}$ to vanish across all cosmic history is the simplest choice and in this work we will stick to this option, following a common choice in recent literature (Kreisch & Komatsu, 2017; Spurio Mancini et al., 2018; Reischke et al., 2019; Denissenya & Linder, 2018). We also notice that, since we will be mainly interested in a parameterization of the $\alpha$ functions that make them vary linearly with the dark energy density fraction (as discussed later in this Section), for our purposes having constraints that make $\alpha_{T}$ vanish in our local Universe also imply that $\alpha_{T}$ vanishes across all cosmic history.

We remark that the other three $\alpha$ functions are still free to vary, although again, we fix $\alpha_{K}$ in our analysis because it is unconstrained by large-scale structure probes and uncorrelated with all other parameters. This leaves us with two Horndeski functions to constrain, $\alpha_{B}$ and $\alpha_{M}$. For details on the status of studies of Horndeski theories after the constraints on $\alpha_{T}$ from gravitational waves, we refer the reader to Ezquiaga & Zumalacárregui (2017), who identify the models within the Horndeski classes that are still viable after GW170817 (see also Kase & Tsujikawa, 2018), and Peirone et al. (2018), who show that even the strict bound on the present-day gravitational wave speed does not exclude at all nontrivial signatures of modified gravity that can be measured in linear cosmic perturbations. Excluding modifications to $\alpha_{T}$ (as we assume throughout our analysis, as well as fixing $\alpha_{K}$ and the expansion history) effectively reduces the large class of Horndeski models to a more limited set of theories; this includes prototypical examples of dark energy and modified gravity theories, such as quintessence, $f(R)$, Brans-Dicke, kinetic gravity braiding models, while more complicated theories such as Galileon models (Ezquiaga & Zumalacárregui, 2017; Ezquiaga & Zumalacárregui, 2018) are excluded.

These $\alpha(\tau)$ functions are all identically vanishing in $\Lambda$CDM: any clear detection of values different from zero for these functions would mark a deviation from the concordance cosmological model. In this paper, we aim at setting constraints on the $\alpha(\tau)$ functions. To achieve this goal, we first need to choose a parameterization for their time evolution. Our first choice is to consider their time evolution to trace that of the dark energy component, to which they are proportional. We set then:

[TABLE]

The first constraints we present in Section 5 will be on the proportionality coefficients $\hat{\alpha}_{i}$. This parameterization is the simplest and most common in the literature (as used e.g. in Planck Collaboration et al., 2016b) and can already provide comprehensive information on Horndeski gravity, as remarked by Gleyzes (2017).

The second parameterization we consider, also studied in Planck Collaboration et al. (2016b), represents a subclass of models, specifically k-essence conformally coupled to gravity. In this parameterization the four $\alpha(\tau)$ functions can be specified through a single function of time $\Omega\left(\tau\right)$ affecting three of the $\alpha(\tau)$ functions:

[TABLE]

where a prime denotes derivative with respect to conformal time. The parameterization for $\Omega(\tau)$ is given by $\Omega_{0}\,a(\tau)$: the only free modified gravity parameter becomes then $\Omega_{0}$.

3 Theoretical power spectra

In this section we summarise the mathematical expressions for the projected power spectra of the three cosmological probes considered in our analysis. Our data vector, which we share with vU18, is composed of estimates of these power spectra, as functions of the tomographic bin and the angular multipole $\ell$. There are multiple advantages in using estimators of power spectra, acting in Fourier space, over real space correlation functions, although the latter may be easier to measure (Köhlinger et al., 2016): we mention here in particular the fact that the covariance matrix of the power spectra is more diagonal than its real-space counterpart, which contributes to a cleaner separation of scales; furthermore, power spectrum estimators can be used to extract the $B$-mode part of the signal, which in cosmic shear studies serves as a systematic check because, in absence of systematics, it should be consistent with zero at lowest order.

For estimation of the power spectra the authors of vU18 followed the formalism originally developed for cosmic shear by Schneider et al. (2002) (and extended by vU18 to galaxy-galaxy lensing and galaxy clustering) based on simple integrals with appropriate weighting functions over the real-space correlation functions, which can be measured with existing public code such as TreeCorr888https://github.com/rmjarvis/TreeCorr (Jarvis et al., 2004), employed by vU18; we refer the reader to vU18 for details.

We conform to vU18 in referring to the power spectra associated to cosmic shear, galaxy-galaxy lensing and galaxy clustering as $P^{\mathrm{E}},P^{\mathrm{gm}},P^{\mathrm{gg}}$, respectively. In the cosmic shear case, the label stands for $E$-mode. vU18 also measured the $B$-mode power spectra and found a tentative signal; however, based on the assumption that a potential systematic affects $E$- and $B$-modes in the same way, they created a new data vector where they substracted the B-modes from the E-modes, and found that this correction shifted their main cosmological inference result by less than $0.5\sigma$. Since $P^{E}$ does not vary rapidly with $\ell$, vU18 only needed a few $\ell$-bins to capture most of the cosmological information. They used five logarithmically-spaced bins, whose logarithmic means range from $\ell=200$ to $\ell=1500$. $P^{\mathrm{gm}}$ and $P^{\mathrm{gg}}$ were estimated adopting the same $\ell$ ranges. Our data vector is unchanged with respect to the one used by vU18 and as such shares all the aforementioned properties. On the theoretical side, however, differently from vU18 we are working on constraints on Horndeski gravity: therefore, we need to present theoretical predictions for the power spectra in a modified gravity context, to be compared with the estimated power spectra in the cosmological inference process.

3.1 Modifications of Poisson equation and ratio of Bardeen potentials

In order to provide equations for the power spectra of the three probes in a general modified gravity scenario, we need to introduce the modifications to the Poisson equation and the ratio of the Bardeen potentials that distinguish a modified gravity theory from General Relativity in absence of anisotropic stress. In Newtonian gauge the linear perturbation equations in Fourier space in a general modified gravity scenario are given by (Saltas et al., 2014; Zhang et al., 2007; Amendola et al., 2008)

[TABLE]

where $\mu$ and $\eta$ are in general functions of both time and scale, and equal to unity in General Relativity in absence of anisotropic stress. In the quasi-static limit, where one neglects time derivatives in the Einstein equations for the perturbations (Sawicki & Bellini, 2015), an explicit formula can be found, connecting the $\alpha(\tau)$ functions directly to the $\mu$ and $\eta$ functions (see e.g. Appendix A in Alonso et al., 2017). However, the validity of the quasi-static approximation depends on the specific modified gravity model considered and the range of scales probed by the cosmological survey considered in the analysis (Baker & Bull, 2015). We stress here that we do not make use of the quasi-static approximation: as described in Section 4, we source the potential and density statistics directly from the Boltzmann code HiClass, which does not assume the quasi-static regime. The values of the $\alpha(\tau)$ functions at different times is a direct output of the code, while the values of $\mu$ and $\eta$ can be easily derived from the evolved potential and density statistics (as done e.g. in Spurio Mancini et al., 2018; Reischke et al., 2019).

The equations presented below are for this general modified gravity scenario, and the ones used in vU18 for $\Lambda$CDM can be deduced from ours by setting $\mu=\eta=1$ identically, which is equivalent to setting all $\alpha(\tau)$ functions to zero. For all our power spectrum estimators we assume the extended Limber approximation, in the form given by Loverde & Afshordi (2008), which uses $\ell+1/2$ in the argument of the matter power spectrum but no additional prefactors. Many recent papers have demonstrated that for cosmic shear these approximations are excellent on the scales that we consider (Kitching et al., 2017; Lemos et al., 2017; Kilbinger et al., 2017), and this is also true for galaxy-galaxy lensing and galaxy clustering power spectra (see vU18). A detailed derivation of the expression for the cosmic shear power spectrum is reported in Appendix A; expressions for galaxy-galaxy lensing and galaxy clustering power spectra can be derived following the same approach.

3.2 Cosmic shear

The weak lensing convergence power spectrum can be obtained from the 3D matter power spectrum $P_{\delta}$ via

[TABLE]

where indices $i,j$ label the tomographic bins, $H_{0}$ is the Hubble constant, $\Omega_{m}$ the present-day matter density parameter, $c$ the speed of light, $\chi$ the comoving distance (we assume spatial flatness here and throughout the paper), $a(\chi)$ the scale-factor, $\chi_{\rm H}$ the Hubble radius today, and $g_{i}(\chi)$ a geometric weight factor, which depends on the source redshift distribution $n_{i}(z)$ for tomographic bin $i$, $n_{i}(z)\,{\rm d}z=n_{i}(\chi)\,{\rm d}\chi$:

[TABLE]

The $n_{i}(\chi)$ distributions are normalised such that $\int\,\mathrm{d}\chi\,n_{i}(\chi)=1\,\forall i$.

3.3 Galaxy-galaxy lensing

The cross-correlation spectrum between lensing convergence and matter density has the following form

[TABLE]

with $n_{\rm F_{i}}(\chi)$ the redshift distribution of the foreground sample labelled by index $i$ (here indices $i$ and $j$ run over foreground samples and cosmic shear tomographic bins, respectively). Similarly to vU18, we use an effective linear bias $b_{i}$ for each foreground sample (see Section 4.5 for details).

3.4 Galaxy clustering

The angular power spectrum can be determined from the matter power spectrum via

[TABLE]

where, as above, $b_{i}$ corresponds to the effective bias of the sample and, as in vU18, we consider only auto-correlations of the tomographic bin power spectra. Our dataset consists of two GAMA foreground samples (see Section 4.2), one at lower and the other at higher redshift. Following Loverde & Afshordi (2008), for our two foreground samples we obtain threshold scales for the validity of the Limber approximation of $\ell\gtrsim 15$ and $\ell\gtrsim 25$, respectively. Since the minimum $\ell$ scale entering the analysis is 150, the Limber approximation is valid for both our samples at subpercent level (see also Kilbinger et al., 2017).

4 Data analysis

In order to obtain constraints on both standard cosmological parameters and those that describe the evolution of linear perturbations in Horndeski gravity (introduced in Section 2), we carry out cosmological inference in a Bayesian framework, sampling the likelihood

[TABLE]

where the indices $\alpha$, $\beta$ run over the probes considered, as well as the tomographic bins and the angular multipole $\ell$. The analytical covariance matrix $\mathbf{C}$ is calculated as outlined below in Section 4.1. Equation (12) assumes that the estimated power spectra are Gaussian distributed around their mean; it must be noted that this is an approximation, rigorously valid only at high $\ell$ as a consequence of the Central Limit theorem. For the multipole range considered here, this approximation is safe. The components of the vector multiplying the covariance matrix are calculated as

[TABLE]

where the dependence on cosmological parameters enters only in the calculation of the predicted power spectra, $\langle P_{\alpha}(\bm{\theta})\rangle^{\rm model}$, which is carried out according to Eqs. 8, 10, 11, while $P_{\alpha}$ denotes the measured power spectrum, with index $\alpha$ running over the cosmological probe and pair of tomographic bins considered.

To constrain the cosmological parameters, we use the sampler MontePython (Audren et al., 2013) 999Version 2.2.2 from https://github.com/baudren/montepython_public and build a likelihood module to analyse the three probes considered. The likelihood module is based on the one developed in Köhlinger et al. (2017) for the KiDS-450 quadratic estimator analysis. Our module has been developed considering as high priority the flexibility to choose different combinations of cosmological probes as well as the possibility to perform the analysis in either a standard $\Lambda$CDM cosmological scenario or an extended modified gravity scenario within the Horndeski framework. These two options are possible as a result of interfacing MontePython with HiClass (Zumalacárregui et al., 2017)101010https://github.com/miguelzuma/hi_class_public, a Boltzmann solver extending Class (Blas et al., 2011)111111Version 2.4.5 from https://github.com/lesgourg/class_public to Horndeski models, and thus allowing the user to source the matter power spectrum in either a $\Lambda$CDM or a Horndeski cosmological model, respectively (see Section 5.1 for details). As a sanity check, we verified that we obtain the same $\Lambda$CDM results when running with HiCLASS (with Horndeski corrections switched off) as when running directly with Class. Our code implementation is publicly accessible121212Our likelihood module will be made publicly available after acceptance of this paper by the Journal..

In contrast, the analysis of vU18 was carried out only in $\Lambda$CDM with a likelihood module built for the sampler CosmoMC  (Lewis & Bridle, 2002) for cosmological parameter estimation. The likelihood code of vU18 was based on the one used in Joudaki et al. (2017) for the fiducial KiDS cosmic shear analysis.

For an efficient evaluation of the likelihood we use two sampling methods, namely MultiNest (Feroz et al., 2009), a multimodal Nested Sampling (Skilling, 2006) algorithm included with a Python wrapper (PyMultiNest; Buchner et al., 2014) in MontePython, and highly parallelised affine invariant sampling through the CosmoHammer suite (Akeret et al., 2012), embedded in MontePython. CosmoHammer embeds in turn emcee, an implementation by Foreman-Mackey et al. (2013) of the affine invariant ensemble sampler by Goodman & Weare (2010). We verified that both MultiNest and CosmoHammer sampling methods yield inference results consistent between each other.

4.1 Covariance

As can be seen from equation (12), errors on cosmological parameters are determined by the parameter covariance matrix $\mathbf{C}$, or more precisely its inverse, the precision matrix $\mathbf{C}^{-1}$. Developing methods for efficient computation of precision matrices is an area of active research (see e.g. Taylor & Joachimi, 2014; Sellentin & Heavens, 2016; Schäfer & Reischke, 2016; Reischke et al., 2017; Heavens et al., 2017; Jeffrey & Abdalla, 2018). Our covariance matrix is unchanged with respect to the one used by vU18, i.e. is determined analytically, following a similar formalism as in Hildebrandt et al. (2017). Here we briefly review the main terms of the calculation, referring to vU18 for further details.

The analytical covariance matrix consists of three terms: (i) a Gaussian term accounting for the Gaussian contribution to sample variance, shape noise, and a mixed noise-sample variance term, estimated following Joachimi et al. (2008), (ii) an in-survey non-Gaussian term from the connected matter trispectrum, and (iii) a super-sample covariance term. The latter accounts for cosmic variance modes larger than the survey window and coupling to smaller modes within. One of the main advantages of the analytical approach to covariance estimation lies indeed in its ability to better account for super-sample covariance. To compute (ii) and (iii), vU18 followed Takada & Hu (2013) and the extension to galaxy-galaxy lensing and galaxy clustering presented e.g. in Krause & Eifler (2017).

The publicly available covariance that we employ makes idealistic assumptions about the shot noise contribution which leads to an under-prediction of variance on large scales where the survey geometry becomes relevant. This was recently shown to be the reason behind the relatively large $\chi^{2}$ values found for previous KiDS and DES analyses (Troxel et al., 2018b). As the posterior was largely unaffected, we can proceed with the available covariance matrix but have to keep the above limitations in mind when interpreting the goodness of fit.

We also notice that, similarly to vU18, we do not consider any parameter dependence for the covariance matrix. vU18 checked that this was choice was safe for their analysis. Here, the covariance matrix is calculated at the same fiducial model chosen by vU18. This is motivated by the fact that we choose a fiducial model for our priors that is given by the same $\Lambda$CDM model chosen by vU18.

4.2 Data

The Kilo-Degree Survey (de Jong et al., 2013) is an ongoing optical imaging survey that, once completed, will have spanned 1350 deg2 of the sky in four optical bands, $u$, $g$, $r$ and $i$. Our cosmic shear power spectra 131313Measured power spectra and covariance matrices used in this analysis and in vU18 are publicly available at http://kids.strw.leidenuniv.nl/sciencedata.php, which we have in common with vU18, are estimated from the KiDS-450 shape measurement catalogues (Hildebrandt et al., 2017; de Jong et al., 2017), containing shape measurements and photometric redshifts of 450 deg2 of data, with galaxies situated in regions of the sky overlapping with the observational window of GAMA. We use the samples associated with the same four KiDS tomographic source redshift bins adopted in Hildebrandt et al. (2017) and vU18, spanning in redshift $z_{\rm B}$ the intervals 0.1 $<z_{\rm B}\leq$ 0.3, 0.3 $<z_{\rm B}\leq$ 0.5, 0.5 $<z_{\rm B}\leq$ 0.7 and 0.7 $<z_{\rm B}\leq$ 0.9. The main properties of the source samples, such as their average redshift, number density and ellipticity dispersion, can be found in Table 1 of Hildebrandt et al. (2017). Importantly, we use the same redshift distributions used in vU18 and Hildebrandt et al. (2017).

The Galaxy And Mass Assembly survey (Driver et al., 2009, 2011; Liske et al., 2015) is a spectroscopic survey of $\sim$240 000 galaxies. vU18 used a subset of $\sim$180 000 galaxies fully overlapping with KiDS. vU18 selected two GAMA samples, a low redshift sample with $z_{\rm spec}<0.2$, and a high redshift sample with $0.2<z_{\rm spec}<0.5$. These two foreground samples are used to determine the galaxy-matter cross-correlation as well as the clustering power spectra; we adopt both from vU18.

4.3 Non-linear structure formation and baryonic feedback model

A proper inclusion of a series of astrophysical systematics and an appropriate modelling of the matter power spectrum on non-linear scales are important to derive accurate cosmological constraints. Here we will describe our choice for the non-linear power spectrum and in particular for the description of baryonic feedback, which modifies the matter power spectrum on small, non-linear scales (e.g. Semboloni et al., 2013).

The effect of non-linear structure formation and baryonic feedback can now be modelled in CLASS using a module called HMcode (Archidiacono, Brieden & Lesgourgues in prep.), based on the results of Mead et al. (2015), who account for baryonic effects (in particular, AGN feedback) by modifying parameters that describe the shape of dark matter haloes. This reflects, for example, what happens with AGN and supernova feedback, which make haloes less concentrated by pushing material out of them. In HMcode this is modelled by choosing a functional form for the mass-concentration relation of the type

[TABLE]

where $z_{\rm f}$ is the formation redshift of a halo, which depends on the halo mass. The free parameter in the fit, $c_{\mathrm{min}}$, modulates the amplitude of this mass-concentration relation. It also sets the amplitude of a ‘halo bloating’ parameter $\eta_{0}$, which changes the halo profile in a mass dependent way (see equation 26 of Mead et al., 2015). vU18 followed the recommendation of Mead et al. (2015) by fixing $\eta_{0}=1.03-0.11c_{\mathrm{min}}$; we follow the same choice. Setting $c_{\mathrm{min}}=3.13$ corresponds to a dark-matter-only model. Marginalizing over $c_{\mathrm{min}}$ allows us to account for the uncertainty on this parameter. This will be done in Section 5 every time we show contour plots for pairs of parameters (the ones not including $c_{\mathrm{min}}$), as these are marginalised over all other parameters except the two considered in each contour plot.

The problem of non-linear contributions in this modified gravity context is even more exacerbated than in $\Lambda$CDM, since clear prescriptions for modified gravity models on non-linear scales are highly model-dependent and currently not available at the required precision. We decided to produce our non-linear corrections using also in this modified gravity case HMcode, which takes into account baryonic effects. The main obstacle in this sense has been the absence of an implementation of HMcode in HiClass. To solve this problem, we developed our own version of HMcode for HiClass, modifying the one provided to us for Class. We remark here that the use of HMcode is not the most rigorous approach to non-linear corrections in modified gravity, because HMcode is a refined version of halofit (Smith et al., 2003), which has been developed for a $\Lambda$CDM model. However, while prescriptions for single modified gravity models exist (see e.g. Winther et al., 2015), there is no general one for the non-linear power spectrum for the whole Horndeski class of dark energy and modified gravity models. By choosing to deal with nonlinearities using HMcode we hope to give a flavour of the increased constraining power that comes from non-linear scales (as done e.g. in Alonso et al., 2017; Spurio Mancini et al., 2018; Reischke et al., 2019, for future Stage IV surveys), however we warn of the necessity of developing proper formalisms for the non-linear matter power spectrum in Horndeski gravity for the purpose of constraining these gravity models with future surveys.

Having said that, we notice that given the size of the error bars obtained in our analysis, and the fact that we already marginalise over the baryon feedback parameter, we are not sensitive to small modifications in the non-linear model. To substantiate this statement, in Appendix C we studied in detail the effect of non-linear corrections on our cosmological constraints, by introducing an additional phenomenological nuisance parameter that mimics the effect of different modified gravity models on the non-linear matter power spectrum. Our main finding is that our cosmological constraints are largely insensitive to this additional parameter, corroborating the strength of our analysis against different possible choices for the non-linear matter power spectrum. This will certainly need to be reassessed in future work, where the same analysis could be repeated with larger datasets, more sensitive to varying non-linear prescriptions for different modified gravity models.

4.4 Screening mechanisms

A screening mechanism is an important ingredient of modified gravity theories. Acting as a non-linear effect, it screens modifications of gravity at small scales to provide agreement with tests of General Relativity in short scales or high-density environments such as the Solar System (see e.g. Koyama, 2016). We model screening mechanisms in a phenomenological way, by applying a scale-dependent filter to the effective Newtonian coupling, the gravitational slip and the linear growth factor:

[TABLE]

This achieves the recovery of General Relativity predictions on scales $k>k_{s}$. A (conservative) value for $k_{s}$ is $0.1\,h\,\mathrm{Mpc}^{-1}$ (Alonso et al., 2017; Spurio Mancini et al., 2018; Reischke et al., 2019); we vary over this scale taking a uniform prior on the logarithm of $k_{s}$ over the interval $[-1,1]$, thus spanning three orders of magnitude for the screening scale.

4.5 Intrinsic alignment model

Intrinsic alignments affect both the cosmic shear power spectrum and the galaxy-matter power spectrum. For the cosmic shear power spectrum, there are two contributions, the intrinsic-intrinsic (II) and the shear-intrinsic (GI) terms, while the galaxy-matter power spectrum has a galaxy-intrinsic (gI) contribution (e.g. Joachimi & Bridle, 2010). We model the intrinsic alignment power spectrum following the non-linear modification of the linear alignment model (Catelan et al., 2001; Hirata & Seljak, 2004; Bridle & King, 2007; Hirata & Seljak, 2010; Blazek et al., 2017):

[TABLE]

with $P_{\delta}(k,z)$ the full non-linear matter power spectrum, $D_{+}(k,z)$ the linear growth factor, normalised to unity at $z=0$ and in general dependent on scale and redshift in modified gravity, $\rho_{\rm crit}$ the critical density, $C_{1}=5\times 10^{-14}h^{-2}M_{\odot}^{-1}$Mpc3 a normalization constant chosen such that $C_{1}\,\rho_{\rm crit}\approx 0.0134$ (Hirata & Seljak, 2004; Bridle & King, 2007; Joachimi et al., 2011), and $A_{\rm IA}$ the overall amplitude, which is a free parameter in our model.

As already mentioned in Secs.3.3 and 3.4, to model $P^{\rm gm}$ and $P^{\rm gg}$, we assume that the galaxy bias is constant and scale-independent, to be interpreted as an effective bias since we also consider non-linear scales. This effective bias is fitted separately for the low-redshift and high-redshift foreground sample. Motivation for not using a scale dependence for the bias comes from combined observations of galaxy-galaxy lensing and galaxy clustering, which found the scale dependence to be small (e.g. Hoekstra et al., 2002; Simon et al., 2007; Jullo et al., 2012; Cacciato et al., 2012). We refer to vU18 for details on how this approximation of scale-independent bias is benign in their analysis as well as in ours. For cosmic shear, in a modified gravity context the intrinsic alignments contributions are given by ($i\leftrightarrow j$ stands for the symmetric term obtained swapping indices $i$ and $j$):

[TABLE]

while for galaxy-galaxy lensing, the intrinsic alignments contribution is given by:

[TABLE]

We note that modifications to the linear growth factor $D_{+}$ and the modified gravity functions $\mu$ and $\eta$ are included for completeness, but their constraining power is limited due to the intrinsic alignment amplitude $A_{\rm IA}$ left free to vary.

4.6 Priors

The cosmological parameters in our parameter space are: the cold dark matter density $w_{\mathrm{cdm}}=\Omega_{\mathrm{cdm}}h^{2}$, the baryon density $w_{\mathrm{b}}=\Omega_{\mathrm{b}}h^{2}$, the amplitude of the primordial power spectrum $A_{s}$ (we vary over $\mathrm{ln}10^{10}\mathrm{A_{s}}$), the value $h$ of the Hubble parameter today divided by 100 km/s/Mpc and the exponent of the primordial power spectrum $n_{s}$. The parameters $\Omega_{\mathrm{m}},\sigma_{8}$ and $S_{8}=\sigma_{8}\sqrt{\Omega_{\mathrm{m}}/0.3}$, shown in the plots in the following Sections, are derived parameters in our analysis. When constraining Horndeski gravity, we vary additionally over either the proportionality coefficients $\hat{\alpha}_{B}$ and $\hat{\alpha}_{M}$, for $\alpha_{B}(\tau)=\hat{\alpha}_{B}\,\Omega_{\mathrm{DE}}(\tau)$ and $\alpha_{M}(\tau)=\hat{\alpha}_{M}\,\Omega_{\mathrm{DE}}(\tau)$, or $\Omega_{0}$, for the parameterization given by equation (5), where the only free function is $\Omega(\tau)=\Omega_{0}\,a(\tau)$.

The astrophysical nuisance parameters (discussed earlier in Secs. 4.3 and 4.5) considered in our parameter space are: the intrinsic alignment amplitude $A_{\rm IA}$, the HMcode parameter $c_{\rm min}$ and the effective bias values for the low-$z$ and high-$z$ foreground samples, denoted with $\mathrm{bias}_{z1}$ and $\mathrm{bias}_{z2}$, respectively. We also vary over a cut-off scale for the screening mechanism, as described in Section 4.4.

We adopt top-hat priors, with ranges specified in Table 1, on all cosmological and physical nuisance parameters. We fix $k_{\rm pivot}$, the pivot scale where the scalar spectrum has an amplitude of $A_{s}$, to 0.05 Mpc*-1*. While massive neutrinos may suppress structure formation in a way that is degenerate with modifications of gravity (Baldi et al., 2014; Spurio Mancini et al., 2018; Reischke et al., 2019; Peel et al., 2018), here we fix the sum of the neutrino mass to zero as done in Hildebrandt et al. (2017) and vU18, and adopt their same prior range. All other prior ranges are also set to reproduce the ones used in vU18.

For the Horndeski parameters, when we consider the proportionality to the dark energy density fraction as time dependence for the $\alpha(\tau)$ functions, we vary both $\hat{\alpha}_{B}$ and $\hat{\alpha}_{M}$ over the uniform range $[-2,+2]$. This is a symmetric range around the $\Lambda$CDM zero value, chosen such that on one edge, it reaches a point of theoretical singularity for the $\alpha(\tau)$ functions: $\hat{\alpha}_{B}=2$ (Noller & Nicola, 2018b).

When we consider the parameterization for the $\alpha(\tau)$ functions described by equation (5), the only free modified gravity parameter becomes $\Omega_{0}$, which parameterizes $\Omega(\tau)=\Omega_{0}\,a(\tau)$. For this parameter we consider the same prior range $\Omega_{0}\in[0,1]$ used in Planck Collaboration et al. (2016b).

5 Results

In this Section we present our cosmological constraints. All of them have been obtained using our new likelihood module for MontePython, which we tested thoroughly against the results obtained by vU18 in a $\Lambda$CDM context. A detailed comparison is described in Appendix D. Here we only comment on the excellent agreement in the parameter $S_{8}=\sigma_{8}\sqrt{\Omega_{m}/0.3}$, derived from $\Omega_{\mathrm{m}}$ and $\sigma_{8}$, i.e. the two parameters whose degenerate combination cosmic shear is most sensitive to (Hildebrandt et al., 2017): for example, the value for $S_{8}$ that we measure considering only cosmic shear is $0.760^{+0.039}_{-0.038}$, to be compared with $0.761^{+0.040}_{-0.038}$ obtained by vU18. We find excellent agreement with the results of vU18 on all parameters and considering all possible combinations of probes. The fact that our implementation is completely independent from the one in vU18 strengthens both analyses and allows for the use of our likelihood module to obtain constraints on modified gravity, after properly modifying the module so that the mathematical expressions implemented there for the power spectra of the three probes match those for modified gravity introduced in Section 3.

5.1 Modified gravity constraints

We investigate the constraining power of our likelihood on dark energy/modified gravity by studying two time parameterizations for the $\alpha(\tau)$ functions, namely the proportionality to the dark energy density fraction described by equation (4) and the planck_linear parameterization introduced in equation (5). Our numerical results for all probe combinations are always compatible with $\Lambda$CDM and are summarised for both parameterizations in Tables  2 and 3, respectively, showing mean values and 68 per cent credibility intervals.

We start with the proportionality to $\Omega_{\mathrm{DE}}(\tau)$, and set constraints on the $\hat{\alpha}$ coefficients. Specifically, we consider $\hat{\alpha}_{B}$ and $\hat{\alpha}_{M}$ on top of our usual cosmological parameters, while we fix $\alpha_{K}$ and $\alpha_{T}$ to zero as explained in Section 2. In Figs. 1, 2 and 3 we overplot the theoretical predictions for the cosmic shear, galaxy - galaxy lensing and galaxy clustering power spectra, respectively, obtained with the best fitting parameters from our combined analysis, in $\Lambda$CDM and in Horndeski gravity; we compare these predictions with the measured power spectra for each probe. We notice how the power spectra produced in Horndeski gravity are characterised by a suppression of the signal towards smaller angular scales, compared to the $\Lambda$CDM spectra. This suggests that our analysis was successful in modelling efficiently the bias factors; if our model had been incorrect (e.g. if the assumption of linear bias had turned out to be insufficient), we would have noticed a steep increase of the power spectra towards higher angular multipoles. As reported in Table 2, the level of agreement between theoretical predictions and observed data does not vary significantly going from $\Lambda$CDM to Horndeski gravity. While the $\chi^{2}$ is large in both cosmological scenarios given the relatively high number of degrees of freedom, indicating a bad fit, this is fully explained by a known under-prediction of shot noise contributions to the covariance at large scales (see Section 4.1 for details).

Fig. 4 shows 68 and 95 per cent marginalised contours on the cosmological parameters that our LSS analysis is most sensitive to, again for the parameterization of the $\alpha(\tau)$ functions that sets them proportional to $\Omega_{\mathrm{DE}}(\tau)$. We overplot results obtained considering either cosmic shear alone or all the three probes together. As we can see from the sharp cutoff in the $\hat{\alpha}_{B}$-$\hat{\alpha}_{M}$ plane, HiClass applies some stability checks that prevent the Markov Chain Monte Carlo (MCMC) from ending up in regions of parameter space where both of these Horndeski parameters are simultaneously negative (see also Fig.7). The investigation of these stability conditions is of the highest priority for future similar analyses (see also Denissenya & Linder, 2018; Noller & Nicola, 2018a). We notice that constraints with a multi-probe approach are tighter, as expected; in particular, the value of $S_{8}$ shrinks from $S_{8}=0.803_{-0.051}^{0.054}$ when measured with cosmic shear alone to $S_{8}=0.835_{-0.039}^{+0.041}$ considering the three probes combined. In both the single and joint analyses the contours for the Horndeski parameters $\hat{\alpha}_{B}$ and $\hat{\alpha}_{B}$ are consistent with $\Lambda$CDM values. Specifically, we find $\hat{\alpha}_{B}=0.05_{-0.36}^{+0.30}\,$ and $\,\hat{\alpha}_{M}=0.46_{-0.61}^{+0.25}$ when considering cosmic shear alone; we retrieve $\hat{\alpha}_{B}=0.20_{-0.33}^{+0.20}\,$ and $\,\hat{\alpha}_{M}=0.25_{-0.29}^{+0.19}$ in our joint analysis of the three probes. Using cosmic shear alone we find a 68$\%$ marginalised contour in the $\hat{\alpha}_{B}$-$\hat{\alpha}_{M}$ plane that covers 11 per cent of the parameter space that is allowed to be explored by the stability conditions applied by HiClass. The same contour obtained with the joint three-probe analysis shrinks down to 6 per cent. Similarly to vU18, we confirm that the increased constraining power with the multi-probe approach comes from the tighter constraints achievable on nuisance parameters. The inclusion of the cross-correlation with the other probes tightens the constraints in particular on the intrinsic alignment amplitude, mostly constrained by the inclusion of galaxy-galaxy lensing in the analysis, which shrinks the error on $A_{IA}$ from $0.44^{+1.16}_{-0.68}$ for cosmic shear alone, to $1.37^{+0.48}_{-0.44}$ with the multi-probe approach. This decrease in the error for $A_{\mathrm{IA}}$ is very similar to the one found by vU18 in their $\Lambda$CDM analysis; in both $\Lambda$CDM and Horndeski scenarios, the cross-correlation of the three probes shrinks the error on the intrinsic alignment amplitude to approximately half its cosmic shear-only value. At the same time, the best fit value for $A_{\mathrm{IA}}$ doubles in the multi-probe approach with respect to the cosmic shear analysis, making the measurement four times more significant. Our combined fit in Horndeski gravity constrains very well the bias parameters for the two foreground samples, similarly to what found in vU18; interestingly, however, our errors in Horndeski gravity are half of those reported by vU18 in $\Lambda$CDM. Constraints on $\sigma_{8}$ and $\Omega_{\mathrm{m}}$ become tighter in the combined analysis, and consequently the same happens for $S_{8}$.

The second parameterization for the $\alpha(\tau)$ functions that we consider is the one called planck_linear in HiClass, which we already introduced in Section 2. In this parameterization the four $\alpha$ functions can be specified through a single function of time $\Omega\left(\tau\right)$ affecting three of the $\alpha$ functions (see equation 5). In this case we do not observe tighter constraints on the modified gravity parameter $\Omega_{0}$ from the inclusion of the multi-probe cross-correlation, while we do find stronger constraints on $\Omega_{\mathrm{m}},\sigma_{8}$ and $S_{8}$, similarly to the case of Horndeski functions proportional to $\Omega_{\mathrm{DE}}(\tau)$. For both parameterizations considered for the $\alpha(\tau)$ functions we find that the screening scale $k_{s}$ is largely unconstrained.

5.2 Comparison and combination with Planck

We investigated the relevance of our analysis in the context of the tension between cosmological parameters estimated with weak gravitational lensing analyses and CMB measurements. Fiducial analyses of the CFHTLenS (Heymans et al., 2013; Joudaki et al., 2017) and KiDS (Hildebrandt et al., 2017; Köhlinger et al., 2017) collaborations have reported values for $\Omega_{\mathrm{m}}$ and $\sigma_{8}$ that are in mild tension with those measured by Planck Collaboration et al. (2016a). Here and throughout our analysis we will refer to the 2015 Planck results because in that case the likelihood modules, to be run with MontePython, are publicly available, thus allowing for re-running of the MCMC. We notice, however, that the latest Planck results (Planck Collaboration et al., 2018) do not differ significantly from the earlier release.

For comparison with Planck, we choose the combination of two likelihoods from the Planck 2015 release: the Plik_lite likelihood for the temperature-only power spectrum in the range $\ell=30-2508$, and the lowTEB likelihood for temperature and LFI polarization information in the range $\ell=2-29$. We choose this likelihood as it is the same used in vU18. For the high-$\ell$ temperature likelihood, we run the Plik lite likelihood, which is pre-marginalised over all nuisance parameters except for the Planck absolute calibration. Noller & Nicola (2018b) demonstrated that running the full high-$\ell$ temperature likelihood or its pre-marginalised version produces equivalent constraints on Horndeski parameters. However, the Plik_lite version produces a remarkable speed up in the analysis, as the number of nuisance parameters for the Planck part is reduced to only one. We have in common with the analysis of Noller & Nicola (2018b) a pure $\Lambda$CDM background which, as the authors of that paper suggest, is one of the reasons behind the excellent agreement between the Plik_lite and full Planck likelihood, since the Plik_lite likelihood is pre-marginalised assuming a $\Lambda$CDM cosmology. However, as also reported in Noller & Nicola (2018b), we remark that extensions of our analysis considering background cosmologies different from $\Lambda$CDM should arguably run the Planck likelihood in its full version.

We compare the results of our LSS analysis, which considers cross-correlations of our three probes, with Planck results. This comparison can be summarised as in Fig. 5, where the two subplots show results for different choices of cosmological model in this comparison. Inference results are presented for $\Omega_{\mathrm{m}}$ and $\sigma_{8}$, the two cosmological parameters that are best constrained by cosmic shear. In our analysis these two parameters are derived from the cold dark matter density $w_{\mathrm{cdm}}=\Omega_{\mathrm{cdm}}h^{2}$ and the rescaled amplitude of the primordial power spectrum $\ln(10^{10}A_{\rm s})$. In the left panel, we show the results of running our LSS probes in $\Lambda$CDM and overplotting the constraints in the $\Omega_{\mathrm{m}}-\sigma_{8}$ plane with those obtained by Planck Collaboration et al. (2016a) also in $\Lambda$CDM. This plot essentially reproduces Fig. 7 in vU18, but the LSS contours are obtained with our own, independently developed likelihood: we recognise the familiar mild ‘tension’ between LSS and CMB measurements. In the right panel of Fig. 5 we analyse the LSS probes in a Horndeski scenario, and also re-run the Planck likelihood in Horndeski gravity. The agreement between LSS and CMB contours is excellent. For comparison, in the same $\Omega_{\mathrm{m}}-\sigma_{8}$ plane we overplot the same $\Lambda$CDM contours for Planck presented in the left subplot, corresponding to the results published in Planck Collaboration et al. (2016a). We notice that there is already a mitigation of the tension between LSS and CMB results when one considers the LSS constraints in Horndeski gravity, and the CMB ones in $\Lambda$CDM. The alleviation of the tension, compared to the $\Lambda$CDM results, originates from the combined effect of the larger Horndeski parameter space, which widens the LSS contours in the $\Omega_{\mathrm{m}}-\sigma_{8}$ plane so that they better match with the $\Lambda$CDM contours from Planck, as well as a shift in the best fit value for the parameters.

To measure the agreement between the $S_{8}$ parameter retrieved by our Horndeski and $\Lambda$CDM gravity analyses, in comparison to the Planck values, we calculate the shift $\Delta S_{8}$ between this parameter measured by the CMB and LSS probes. Considering cosmic shear alone, we report a value $\Delta S_{8}=0.048_{-0.056}^{+0.059}$ in Horndeski, compared to $\Delta S_{8}=0.091_{-0.045}^{+0.046}$ in $\Lambda$CDM. In the joint three-probe analysis, we find $\Delta S_{8}=0.016_{-0.046}^{+0.048}$ in Horndeski gravity and $\Delta S_{8}=0.059_{-0.039}^{+0.040}$ in $\Lambda$CDM.

As a last step in the comparison of our LSS analysis with the Planck CMB results, we combine the LSS and CMB likelihoods. When doing this, we initially extended the Planck chains to the inclusion of CMB lensing; while this could in principle induce correlations between Planck and KiDS, we found that CMB lensing had no impact on the results presented; therefore, we decided to report results obtained without the addition of CMB lensing, as they do not differ significantly from the ones obtained with CMB lensing. Crucially, we remark here that we are allowed to combine the LSS and CMB likelihoods because, contrary to what happened when analysing LSS and CMB probes in $\Lambda$CDM (as in Hildebrandt et al. 2017 and vU18), our Horndeski analysis produces constraints on the cosmological parameters, for LSS and CMB, that are not in tension between them, as shown in Fig. 5.

The two panels of Fig. 6 show the increase in sensitivity of the combined LSS+CMB analysis to key cosmological parameters (left panel) including the Horndeski parameters $\hat{\alpha}_{B}$ and $\hat{\alpha}_{M}$, and astrophysical nuisance parameters (right panel) such as the bias coefficients, the intrinsic alignment amplitude and the HMcode parameter $c_{\mathrm{min}}$, resulting from running our LSS chains in combination with the Planck CMB ones. We find that, as expected, combining LSS and CMB datasets reduces uncertainties on all cosmological parameters. In particular, Planck measurements allow for better constraints on parameters that are weakly constrained by LSS probes, such as $n_{s}$ and $h$ (see Table 2). The Horndeski parameters $\hat{\alpha}_{B}$ and $\hat{\alpha}_{M}$ are also better constrained in the combined analysis, although we notice that most of the constraining power on these two parameters comes from the LSS probes. This can be seen better in Fig. 7, where we plot the region of parameter space in the $\hat{\alpha}_{B}-\hat{\alpha}_{M}$ plane that is allowed to be explored by the stability conditions applied by HiClass (which ensure that the choice of the Horndeski parameters corresponds to a theoretical model that is safe from ghost, gradient and tachyon instabilities; for details, see Zumalacárregui et al., 2017), as well as the marginalised 68 and 95 per cent contours that can be obtained on these parameters using either our joint LSS analysis, re-running the Planck likelihood in Horndeski gravity or running together the Planck and LSS likelihoods. The combination of LSS and CMB likelihoods increases the sensitivity to the Horndeski parameters, with the 68 per cent contour covering 5 per cent of the allowed parameter space in the $\hat{\alpha}_{B}$ - $\hat{\alpha}_{M}$ plane, although the contours obtained with the combined LSS+CMB analysis are comparable in size with the ones obtained from the LSS probes alone. We remark that the contours appear to be centered away from the $\Lambda$CDM point $\{0,0\}$ (while still being consistent with this value) due to the skewness in the marginal distributions of the parameters; this is due to stability conditions applied by HiClass, which artificially cut out a part of parameter space. A plot comparing the contours on all cosmological and nuisance parameters obtained with our LSS analysis and those produced running the LSS and Planck likelihood together can be found in Appendix B. The numerical values for the mean and 68 per cent credibility intervals for all cosmological and nuisance parameters in the combined LSS+CMB analysis are reported in Table 2, to be compared with the LSS-only and Planck-only results, also reported in Table 2. The exact number of degrees of freedom for the Planck likelihood is not readily available in the literature (see e.g. Planck Collaboration et al., 2016a) and its non-trivial calculation is beyond the scope of this work. We reproduced the value of $\chi^{2}$ obtained by the official Planck analysis in $\Lambda$CDM, which we found to be very similar to the one we obtained in Horndeski gravity. In the combined large-scale structure and CMB analysis we notice that the large-scale structure probes add a negligible contribution to the total number of degrees of freedom compared to Planck; since we find that the $\chi^{2}$ for the combined large-scale structure and CMB analysis increases by $\sim 1\%$ from the CMB-only value, we conclude that there is no indication that the goodness of fit degrades in the combined analysis with respect to considering Planck alone.

6 Discussion and conclusions

In this paper we analysed data from the KiDS and GAMA surveys, in a joint framework for tomographic cosmic shear, galaxy-galaxy lensing and angular clustering power spectra. Our numerical implementation reproduces the $\Lambda$CDM results of vU18 and extends their analysis to a Horndeski scenario: we set constraints on both standard cosmological parameters as well as functions that fully describe the evolution of linear perturbations in Horndeski gravity. We considered two different time parameterizations for these functions: one where they are proportional to the dark energy density fraction and another one where there is effectively only one free function which we set proportional to the scale factor; the latter represents a subset of models within the Horndeski class.

For all the combinations of probes considered, we found that our constraints are compatible with $\Lambda$CDM values: specifically, all Horndeski parameters are compatible with zero (their GR value) within their 68 per cent marginalised errors. We notice that while the inclusion of cross-correlations of the three probes tightens the constraints on both cosmological and nuisance parameters, as expected our constraints are rather large per se, for multiple reasons. Our parameter space is very large, because we do not fix any background values for the standard cosmological parameters; we vary over them as well. Additionally, while the KiDS data that we used represent state-of-the-art imaging data that are currently available from a cosmic shear survey, the survey volume considered is still relatively small; we do expect increased constraining power from larger data releases. We remark that, while the constraining power of the datasets considered in our analysis is not sufficient to produce very strong constraints on the values of the Horndeski $\alpha(\tau)$ functions (or more precisely, on the parameters that describe one of the possible time parameterizations for these functions), there are some important conclusions that we can draw from the constraints we found.

In particular, assuming the commonly studied relation of proportionality with $\Omega_{\mathrm{DE}}$ for the $\alpha(\tau)$ functions, we found that the values of both proportionality coefficients $\hat{\alpha}_{B}$ and $\hat{\alpha}_{M}$ are close to zero and smaller than unity. Following the line of reasoning presented in Bellini et al. (2016), if one assumes the time parameterization of proportionality to $\Omega_{\mathrm{DE}}(\tau)$, one would expect the coefficients $\hat{\alpha}$ to be $\mathcal{O}(1)$, because in this parameterization the $\alpha(\tau)$ functions are driven by the same functions of the scalar field as the energy density and its derivatives. The fact that we found values for the proportionality coefficients $\hat{\alpha}$ significantly smaller than unity can, in this sense, be regarded as a substantial observational constraint on dark energy/modified gravity derived from this analysis.

We compared and combined our LSS constraints on Horndeski gravity with the Planck analysis of the CMB. Running our LSS and the Planck CMB likelihood separately, we found that the enlarged Horndeski parameter space for the LSS probes helps reduce the tension in cosmological parameters with the CMB measurements obtained in $\Lambda$CDM; the agreement becomes excellent if we compare our LSS Horndeski results with the Planck chains re-run in Horndeski gravity. Since separately these LSS and Planck constraints in Horndeski gravity are not in tension, we are allowed to combine the LSS and CMB likelihoods. The constraints obtained by this combined analysis are tighter as expected, although we notice that most of the constraining power on the Horndeski parameters comes from the LSS likelihood; the Planck likelihood helps constrain much better other standard cosmological parameters, such as $n_{s}$ and $h$, that are weakly constrained by the LSS probes. When combining LSS and CMB likelihoods, we keep the assumption of flatness that we also imposed on our LSS-only analysis. While flatness is predicted by inflation, future analyses should allow for the possibility of non-vanishing curvature when combining LSS and Planck likelihoods.

For the non-linear matter power spectrum we followed the official KiDS prescription, using the Mead et al. (2015) correction implemented in HMcode, which in turn is based on halofit. A bespoke, accurate prescription for the non-linear matter power spectrum in the full model space covered by Horndeski gravity is not yet available. However, note that we are exploring constraints around the $\Lambda$CDM values, where halofit is fully calibrated, and that we apply the non-linear correction to the linear power spectrum rigorously calculated by HiClass in a modified gravity scenario. We argue that the physical mechanism implemented in halofit produces a reasonable degree of predictive power on the extended cosmological models considered, appropriate for the currently still fairly low constraining power of the data. This statement is corroborated by the analysis in Appendix C, which shows quantitatively the fairly low sensitivity of our data to the modelling uncertainty of the non-linear matter power spectrum in Horndeski gravity. Similarly, the screening mechanism has been implemented in a phenomenological way due to a lack of an evolved screening mechanism in the linear formalism of the $\alpha(\tau)$ functions. We remark that future work in these directions is needed (see e.g. Cataneo et al. 2019) to provide tighter constraints with either larger data releases of the current survey generation or future data from Stage IV surveys such as Euclid (Laureijs et al., 2011)141414https://www.euclid-ec.org/ and the Large Synoptic Survey Telescope (LSST Science Collaboration et al., 2009)151515https://www.lsst.org/.

One point to highlight is the difficult interpretation of parameter space as far as the Horndeski functions are concerned. Assuming proportionality to $\Omega_{\mathrm{DE}}$, in our analysis these parameters were let free to vary over negative values as well. However, we found a very distinctive cutoff in the $\hat{\alpha}_{\mathrm{B}}$-$\hat{\alpha}_{\mathrm{M}}$ plane, which prevents the chains to end up in regions of parameter space where both proportionality coefficients are significantly negative. This is due to some stability checks (in particular those concerning the positivity of the speed of sound) which reject points ending up in those regions of parameter space. Those regions are, however, interesting because they may refer to models of relevance, such as some $f(R)$ models for which the relation $\alpha_{B}=-\alpha_{M}$ holds. While we are aware that some of these stability conditions may in some cases be harmlessly bypassed, because they reject points in parameter space even when they produce a negative speed of sound only at very early times, the exploration of these stability conditions is an area of active research. The recent analysis of Denissenya & Linder (2018), for example (but see also a discussion on this in Kreisch & Komatsu, 2017), has investigated the relation between parameterization and stability in Horndeski gravity, showing that the relation is not trivial. In particular, they have shown how stability evolves with redshift, picking out different regions of parameter space that can have complex structure. These end up forming disconnected ’islands’ in parameter space, which are significantly sensitive to the time evolution assumed for the $\alpha(\tau)$ functions.

The very recent analysis of Noller & Nicola (2018a) suggests that radiative stability places strong constraints on the Horndeski parameter space; they used this condition to obtain tight constraints from the analysis of data from the Planck, SDSS/BOSS and 6dF surveys. Noller & Nicola (2018b) also presented cosmological constraints on Horndeski gravity using CMB, redshift space distortions, matter power spectrum and BAO measurements from the Planck, SDSS/BOSS and 6dF surveys. While our results are obtained considering different cosmological probes with respect to theirs, interestingly our constraints compare well with those presented in their analysis; this is true in particular in the $\hat{\alpha}_{B}-\hat{\alpha}_{M}$ plane, when considering the time dependence for the Horndeski parameters that makes them proportional to the dark energy density fraction. We remark that the investigation of the theoretical stability conditions for the $\alpha(\tau)$ functions seems a high priority issue to be addressed in future extensions of our analysis.

To conclude, we stress that, in view of future work, the importance of our KiDS+GAMA analysis lies in two main aspects. On the one hand, we have developed a first likelihood module for Horndeski gravity that can analyse all possible cross-combinations of cosmic shear, galaxy-galaxy lensing and galaxy clustering power spectra data, as can be produced by current and future surveys. On the other hand, our likelihood module can also be run in a standard $\Lambda$CDM scenario and in this regard our results are in excellent agreement with those from the fully independent analysis pipeline of vU18, representing a benchmark for our implementation as well as a further cross-check of their fiducial KiDS+GAMA analysis. In the future it will be interesting to carry out a similar analysis with larger data releases from the KiDS survey and/or with data from Stage IV surveys such as Euclid; we expect that the work presented in this paper and its publicly available implementation will represent a useful tool for these or similar analyses.

Acknowledgements

ASM acknowledges financial support from the graduate college Astrophysics of cosmological probes of gravity by Landesgraduiertenakademie Baden-Württemberg. ASM and BJ also acknowledge support by the UCL Cosmoparticle Initiative. We thank the German DFG Excellence Initiative for funding a mobility program, through the grant Gravity on the largest scales and the cosmic large-scale structure, within which this work originated. FK acknowledges support from the World Premier International Research Center Initiative (WPI), MEXT, Japan and from JSPS KAKENHI Grant Number JP17H06599. Most of the numerical results of this work have been derived with software run on the Hypatia cluster at University College London: we are grateful to Hiranya Peiris for making this possible and to Edd Edmondson for the technical support with it. We thank Carlos García García, Emilio Bellini and Miguel Zumalacárregui for their help with MontePython and Pedro Ferreira, Steven Gratton, Shahab Joudaki, Graça Rocha and Patricio Vielva for useful discussions.

Based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under programme IDs 177.A-3016, 177.A-3017 and 177.A-3018. We use cosmic shear measurements from the Kilo-Degree Survey (Kuijken et al., 2015; Hildebrandt et al., 2017; Fenech Conti et al., 2017), hereafter referred to as KiDS. The KiDS data are processed by THELI (Erben et al., 2013) and Astro-WISE (Begeman et al., 2013; de Jong et al., 2015). Shears are measured using lensfit (Miller et al., 2013), and photometric redshifts are obtained from PSF-matched photometry and calibrated using external overlapping spectroscopic surveys (Hildebrandt et al., 2017). GAMA is a joint European-Australasian project based around a spectroscopic campaign using the Anglo-Australian Telescope. The GAMA input catalogue is based on data taken from the Sloan Digital Sky Survey and the UKIRT Infrared Deep Sky Survey. Complementary imaging of the GAMA regions is being obtained by a number of independent survey programs including GALEX MIS, VST KiDS, VISTA VIKING, WISE, Herschel-ATLAS, GMRT, and ASKAP providing UV to radio coverage. GAMA is funded by the STFC (UK), the ARC (Australia), the AAO, and the participating institutions. The GAMA website is http://www.gama-survey.org/.

Appendix A Derivation of modified gravity power spectrum for cosmic shear

Here we present a detailed derivation of equation (8), i.e. the cosmic shear power spectrum for a generic modified gravity theory parameterised by the functions $\mu$ and $\eta$ introduced in Eqs. 6 and 7. Expressions for galaxy-galaxy lensing and galaxy clustering (corresponding to Eqs. 10 and 11) can be derived following a similar procedure.

Cosmic shear power spectra can be most elegantly derived from expressions for the lensing potential power spectra, since the shear is linearly defined in terms of the lensing potential, through its second derivatives. Moreover, to avoid complications arising from the spin-2 tensorial nature of the shear field, we can compute power spectra for the convergence, which is a scalar field similarly related in a linear way to the lensing potential. Convergence and shear power spectra are equal in the flat sky approximation and differ only as a result of a different multipole-dependent prefactor (Kilbinger et al., 2017), which for $\ell>45$ represents a difference of less than $0.1\%$ between the two power spectra; given our multipole range, we can safely ignore this difference.

As remarked in Section 3, a generic modified gravity scenario is characterised by Bardeen potentials $\Phi$ and $\Psi$ in equation (2), which can differ from each other. Thus, the lensing potential $\phi_{i}$ in tomographic bin $i$ can be defined as

[TABLE]

where $g_{i}$ is the geometric weight factor for bin $i$ introduced in equation (9), and the integration is carried out along the line of sight in Born approximation. The correlation of the lensing potential in tomographic bins $i$ and $j$ is thus

[TABLE]

where we made use of the extended Limber approximation, and $P_{\Phi+\Psi}$ is the power spectrum of the sum of the Bardeen potentials. Since the convergence $\kappa$ is related to the Laplacian of the lensing potential:

[TABLE]

rewriting equation (21) in terms of the convergence gives us

[TABLE]

where the $k^{-4}$ factor comes from equation (22) expressed in Fourier space. Note that this is a two dimensional Fourier vector on the flat sky. It can be augmented by a third dimension in the thin lens approximation (Bartelmann & Schneider, 2001). Using now the (modified) Poisson equation (6), to link the gravitational potential to the overdensity field, and recalling that the two Bardeen potentials $\Phi$ and $\Psi$ are related through $\eta$ (equation 7), we find that

[TABLE]

with $P_{\delta}$ the matter power spectrum. Thus,

[TABLE]

which leads to equation (8).

Appendix B Full posterior distribution

Fig. 8 shows the 68 and 95 per cent marginalised contours for all the cosmological and nuisance parameters listed in Table 1. We show the constraints that we obtain from our joint cosmic shear, galaxy-galaxy lensing and galaxy clustering analysis in Horndeski gravity with parameterization $\alpha(\tau)=\hat{\alpha}\Omega_{\mathrm{DE}}(\tau)$, as well as the contours produced from running the same LSS joint anlaysis together with Planck likelihoods, as discussed in Section 5.2. For parameters such as $n_{s}$ that are not well constrained by the LSS probes considered here, the combination of the LSS likelihood with the Planck CMB one essentially does not add much more information with respect to the one already contained in the CMB-only analysis. Conversely, for parameters such as the intrinsic alignment amplitude that are constrained only by LSS probes, the constraints do not benefit substantially from running the LSS and CMB likelihoods together. The biggest improvement that is determined by this joint LSS-CMB analysis is on parameters that are well constrained singularly by LSS and CMB datasets, such as $\Omega_{\mathrm{m}}$ and $\sigma_{8}$. Importantly for the analysis of this paper, this improvement also affects the Horndeski $\hat{\alpha}$ parameters.

Appendix C Effect of non-linearities on cosmological constraints

As remarked in the main text, currently there is no general prescription for the treatment of non-linear scales in cosmologies alternative to $\Lambda$CDM. Some of our data are, however, at scales affected by such non-linearities. In our analysis we decided to employ a non-linear prescription for the matter power spectrum which is the same as the one followed by the KiDS collaboration, i.e. based on the Mead et al. (2015) prescription implemented in HMcode, which is in turn a modification of halofit (Smith et al., 2003). We already warned the reader that this should be interpreted as a first order correction, meant to give an idea of the constraining power of the data, while being formally incorrect as it lacks of generality for the class of dark energy and modified gravity models considered in this analysis. To support the statement that the Mead et al. (2015) prescription is a sensitive correction to apply as a first order approach, we carried out a phenomenological study of the effect of modifications of the non-linear matter power spectrum. The idea is to parameterize our ignorance by introducing an additional nuisance parameter $\zeta_{NL}$ that quantifies the uncertainty on the non-linear matter power spectrum. This parameter only becomes ‘active’ when the Boltzmann code HiClass switches to HMcode to source the non-linear matter power spectrum; this happens at a certain redshift-dependent scale value $k_{\sigma}(z)$. When this happens, a phenomenological function of scale and redshift, $\beta(k,z)$, defined as

[TABLE]

modifies the non-linear matter power spectrum $P(k,z)$ sourced from HMcode into its ‘corrected’ version $P_{\rm{CORR}}(k,z)$

[TABLE]

to account for possible deviations from $\Lambda$CDM. We vary $\zeta_{NL}$ uniformly in the range $[-5,5]$, and marginalise over $\zeta_{NL}$ when presenting constraints on cosmological parameters. We also studied what happened when we varied the logarithm of $\zeta_{NL}$ uniformly in the range $[-5,0]$, and found similar results. Fig. 9 shows a comparison between the cosmological constraints obtained with our joint three-probe analysis, with and without the implementation of the additional parameter $\zeta_{NL}$. We notice that the constraints on all other parameters do not change significantly, suggesting that our data are only marginally affected by different non-linear prescriptions for the matter power spectrum. Hence, our choice of using HMcode appears sensible and major modifications will only be required with future data at even more non-linear scales, as coming from future releases of the KiDS survey or from Stage IV surveys such as Euclid and LSST.

Appendix D Comparison with vU18

In Figs. 10 and 11, we select two plots of comparison with the constraints obtained by vU18, which we overplot to ours. Specifically, in the two figures we show marginalised 68 and 95 per cent contours considering cosmic shear alone and the combination of all three probes, respectively. We find excellent agreement with the results of vU18 not only with these two choices of probes combination, but for all the possible combinations. Table 4 summarises this comparison: it shows the mean and 68 credibility intervals obtained from our analysis and those obtained from vU18. We notice in particular the excellent agreement in the parameter $S_{8}=\sigma_{8}\sqrt{\Omega_{\mathrm{m}}/0.3}$, derived from $\Omega_{\mathrm{m}}$ and $\sigma_{8}$, i.e. the two parameters whose degenerate combination cosmic shear is most sensitive to (Hildebrandt et al., 2017). Our results are obtained with a completely likelihood different implementation, which strengthens both the validity of the fiducial KiDS+GAMA analysis and the consistency of our likelihood.