© The Authors 2024

1 Introduction

We have entered a new era in gravitational physics in which it is now possible to test and exploit general relativity (GR) on a wide range of scales. The successes of Solar System constraints, and the precision measurements arising from observations of millisecond pulsars can now be combined with detections of gravitational waves and images of black hole shadows (Berti et al. 2015). To this battery of techniques must be added cosmological constraints using the large-scale structure of the Universe (Ferreira 2019).

There have been multiple attempts at constraining GR with cosmological surveys (see, for example Aghanim et al. 2020; Mueller et al. 2018; Joudaki et al. 2018; Lee et al. 2021; Raveri et al. 2021; Nguyen et al. 2023). With spectroscopic and imaging surveys of galaxies, combined with measurements of the cosmic microwave background (CMB), it has been possible to map out gravitational potentials over an appreciable part of the Universe as well as determine the growth rate of the structure and its morphology out to redshift z ~ 2. The resulting constraints have been an important first step in understanding gravity in an altogether untested regime, but they have been underwhelming. For example, the constraints on the simplest scalar-tensor modification to GR, Jordan-Brans-Dicke (JBD) gravity, are more than order of magnitude weaker (Sola et al. 2019, 2020; Joudaki et al. 2022; Ballardini et al. 2022) than those obtained from millisecond pulsar observations (Voisin et al. 2020), but on very different scales.

The relevance of cosmological measurements for gravitational physics is about to change with the upcoming generation of large-scale structure surveys (Ferreira 2019). By mapping out vast swathes of the Universe with exquisite precision, it is hoped that it will be possible to substantially tighten the constraints on gravitational physics on the largest observable scales. Of particular importance in this new vanguard is the Euclid mission. The Euclid satellite will undertake two key complementary surveys: a spectroscopic survey of galaxies and an imaging survey (targeting weak lensing; it can also be used to reconstruct galaxy clustering using photometric redshifts). The primary goal is to determine the nature of dark energy and it is ideally suited for cosmological constraints on gravity (Euclid Collaboration 2020, EC19 hereafter).

Given the potential of the Euclid mission it is imperative to assess its ability to constrain GR. One way of doing this is by assessing how well it will be able to constrain specific extensions of GR, in particular, modified gravity models (MG). MG models are particularly constrained by small-scale experiments such as those in the Solar System where fifth-force effects (Bertotti et al. 2003) and a violation of the equivalence principle (Williams et al. 2012; Touboul et al. 2017) have been thoroughly investigated. As a result, extensions of GR that preserve the equivalence principle in the Earth’s environment and exclude large fifth-force effects on Solar System scales must shield small-scales from potential deviations from GR on large cosmological distances. New physical effects would only reveal themselves on large scales where the Euclid mission would have the potential to unravel them. Screenings of GR extensions have been broadly classified into three families: the chameleon (Khoury & Weltman 2004), k-mouflage (Babichev et al. 2009) and Vainshtein (Vainshtein 1972) mechanisms (see, for a recent review, Brax et al. 2021). For chameleons, screening occurs in regions of space where Newton’s potential is large enough whilst for k-mouflage and Vainshtein, this takes place where the first or second spatial derivatives of Newton’s potential are also large enough. One particular example of an extensively studied theory with the chameleon screening property is f(R) gravity (Carroll et al. 2004; Hu & Sawicki 2007b). It has been shown that the Euclid mission, using the combination of spec-troscopic and photometric probes, will be able to distinguish this model from ΛCDM at more than 3σ confidence level, for realistic fiducial values of its free model parameter, f_R0 (Casas et al. 2023). Derivative screening mechanisms, that is, k-mouflage and Vainshtein, have not been investigated within the context of the Euclid mission, and this will be one of the outcomes of the present work.

In this paper, we forecast how well the surveys from the Euclid mission can be used to constrain a family of theories that modify the theory of GR, but retain one of its properties: a scale-independent linear growth rate. The three theories we consider are: JBD gravity (Brans & Dicke 1961), which is the simplest scalar-tensor theory and involves a non-minimal coupling between a scalar field and the metric; Dvali-Gabadadze-Porrati gravity (DGP, see Dvali et al. 2000), which is a braneworld model that introduces modifications on cosmological scales and screens with the Vainshtein mechanism; and a k-mouflage model (KM, see Babichev et al. 2009; Brax & Valageas 2014, 2016), that is a scalar-tensor theory with a non-canonical kinetic energy and the k-mouflage screening property. The JBD theory is not screened, and to be consistent with current astrophysical constraints on fifth forces, it requires that the Brans-Dicke (BD) parameter ω_BD > 4 × 10⁴ (Bertotti et al. 2003). This is larger than the range of parameters that can cosmologically be tested, as we show below. It implies that the JBD model investigated here must be taken as a template for large-scale deviations against which we compare the DGP and k-mouflage models. All these modifications of GR affect in one way or another the expansion rate and growth of structure and they are therefore prime candidates to be constrained by data from the Euclid mission.

We structure this paper as follows. In Sect. 2 we recapitulate the essential facts about linear cosmological perturbations in the context of extensions to GR, identifying the phenomenological time-dependent parameters, (µ, η, and Σ), that are fed into the linear evolution equations. We then describe the three candidate theories we explored by laying out their corresponding actions, background evolution, functional forms of {µ, η, and Σ}, and how all this is implemented numerically. In Sect. 3 we explain in detail how we calculated all aspects of the theoretical model predictions that go into the forecasting procedure. In Sect. 4 we describe the survey specifications and how they are integrated in the analysis method; in this paper, we use a Fisher forecasting approach. In Sect. 5 we present the results of our methods, and we conclude in Sect. 6.

2 Linearly scale-independent modified gravity

We followed the Bardeen formalism (Bardeen 1980; Ma & Bertschinger 1995) and defined the infinitesimal line element of the flat Friedmann-Lemaître-Robertson-Walker (FLRW) metric by $$\text{d}{s}^2=-(1+2\text{Ψ})c^2\text{d}{t}^2+a^2(t)(1-2\text{Φ}){\delta }_{ij}\text{d}{x}^i\text{d}{x}^j,$$(1)

where a(t) is the scale factor as a function of the cosmic time t, Ψ and Φ are the two scalar potentials and c is the speed of light. We work in Fourier space, where modes are functions of t and the comoving wave number kⁱ. At linear order, the matter energy-momentum tensor can be decomposed as $$T_0^0=-\overline{\rho }(1+\delta )c^2,T_i^0=\left(\overline{\rho }+\overline{p}/c^2\right)c{v}_i,T_j^i=(\overline{p}+\delta p){\delta }_j^i+{\text{Σ}}_j^i,$$(2)

where $\delta =\rho /\overline{\rho }-1$ is the energy density contrast with ρ being the matter density and $\overline{\rho }$ its background value, $p=\overline{p}+\delta p$ is the pressure with $\overline{p}$ the background value, v_i is the peculiar velocity and ${\Sigma }_j^i$ is the (traceless) anisotropic stress, ${\Sigma }_i^i=0$. In the following we work with the scalar component of the matter anisotropic stress and the comoving density perturbation, defined as $\left(\overline{\rho }{c}^2+\overline{p}\right)\left({\widehat{k}}^i\cdot {\widehat{k}}_j-1/3{\delta }_j^i\right)\sigma ={\Sigma }_j^i$ with ${\widehat{k}}_i=k_i/k$ and $\overline{\rho }\Delta \equiv \overline{\rho }\delta +3\left(aH/k/c^2\right)\left(\overline{\rho }+\overline{p}/c^2\right)v$, where v is the velocity potential defined through v_i = -ιkiv/k. Here H ≡ ȧ/a is the Hubble function and a dot stands for the derivative with respect to the coordinate time. In analogy to the ΛCDM model where curvature is found to be compatible with zero, we assumed a flat spatial geometry.

In models in which gravity is modified by a scalar field, the relations between the gravitational potentials and the matter perturbations are modified. These deviations from GR can be encoded into two functions, defined as follows: $$-k^2\text{Ψ}=\frac{4\pi G_{\text{N}}}{c^2}{a}^2\mu (a,k)\left[\overline{\rho }\text{Δ}+3\left(\overline{\rho }+\overline{p}/c^2\right)\sigma \right],$$(3) $$\begin{array}{c}-k^2(\text{Φ}+\text{Ψ})=\frac{8\pi G_{\text{N}}}{c^2}{a}^2\{\text{Σ}(a,k)\left[\overline{\rho }\text{Δ}+3\left(\overline{\rho }+\overline{p}/c^2\right)\sigma \right]\\ -\frac{3}{2}\mu (a,k)\left(\overline{\rho }+\overline{p}/c^2\right)\sigma \},\end{array}$$(4)

where G_N is Newton’s gravitational constant. Eqs. (3) and (4) can be obtained using the quasi-static approximation, which considers scales smaller than the horizon and the sound-horizon of the scalar field, where time derivatives become negligible with respect to spatial derivatives. In a given model, it allows us to determine the functional form of µ and Σ analytically (see, e.g. Silvestri et al. 2013; Bellini & Sawicki 2014; Gleyzes et al. 2015; Zucca et al. 2019; Pace et al. 2021). A further function that can be introduced is therefore the function that defines the ratio of the two potentials, η = Φ/Ψ. In the absence of anisotropic stress the three phenomenological functions are related by the following expression: $$\text{Σ}=\frac{1}{2}\mu (1+\eta ).$$(5)

The phenomenological functions µ, η and Σ are identically equal to 1 in the GR limit. In general, they are functions of time and scale. The models that we consider in this paper preserve the scale-independent (linear) growth pattern, that is, they have µ = µ(a), Σ = Σ(a) and η = η(a).

For a given theory, µ and Σ can be determined numerically after solving for the full dynamics of linear perturbations via an Einstein-Boltzmann solver. This can be achieved with hi_class (Zumalacárregui et al. 2017; Bellini et al. 2020) or EFTCAMB (Hu et al. 2014; Raveri et al. 2014), which implement the effective field theory formalism for dark energy into the standard CLASS (Blas et al. 2011; Lesgourgues 2011) and CAMB (Lewis et al. 2000) codes, respectively. These codes have been validated as part of an extended code comparison effort (Bellini et al. 2018). Alternatively, one can opt for the quasi-static (QS) limit, that is, scales sufficiently small to be well within the horizon and the sound-horizon of the scalar field, and derive the phenomenological functions analytically. In this case, one can use the MGCAMB patch to CAMB (Zhao et al. 2009; Hojjati et al. 2011; Zucca et al. 2019).

In the following, we introduce the three models under consideration, that is, JBD, DGP, and KM. We provide the background evolution equations and the expressions for µ, η, and Σ for each of these models.

2.1 Jordan-Brans-Dicke gravity

The JBD theory of gravity (Brans & Dicke 1961) is described by the following action, $$S_{\text{BD}}={{\displaystyle \int }}^{\text{}}{\text{d}}^{4}x\sqrt{-g}\left[\frac{c^4}{16\pi }\left(\varphi R-\frac{{\omega }_{\text{BD}}}{\varphi }{g}^{\mu \nu }{\partial }_{\mu }\varphi {\partial }_{\nu }\varphi -2\text{Λ}\right)+{ℒ}_{\text{m}}\right],$$(6)

where 𝑔_µν is the metric and R its associated Ricci scalar, ϕ is the JBD scalar field (which has the dimensions of the inverse of Newton’s constant), ω_BD is the dimensionless BD parameter, and 𝓛_m is the matter Lagrangian which is minimally coupled to the metric. We have included a cosmological constant, A. The model has only one free parameter, ω_BD; in the limit in which ω_BD → ∞, the scalar field is frozen and we recover Einstein gravity.

The JBD theory of gravity is remarkably simple in that it depends on so few parameters. However, cosmological constraints on JBD gravity can have wider implications when we take the view that it is the long-wavelength, low-energy limit of more general scalar-tensor theories (see Joudaki et al. 2022, for example). Furthermore, more general scalar theories may be endowed with gravitational screening, which appears on smaller scales, or alternatively, regions of high density, for example. As a consequence, local constraints will to some extent decouple from more global, or large scale, constraints. Thus, constraints on cosmological scales on JBD theory may cover a broad class of scalar-tensor theories, and furthermore, be independent from model-specific constraints on small-scales. Thus, while simple, JBD gravity is a powerful tool for constraining general classes of scalar-tensor theories using cosmological data.

The modified Einstein field equations are (Clifton et al. 2012) $$\begin{array}{c}{G}_{\mu \nu }=\frac{8\pi }{c^4}\frac{1}{\varphi }{T}_{\mu \nu }+\frac{{\omega }_{\text{BD}}}{{\varphi }^2}\left({\nabla }_{\mu }\varphi {\nabla }_{\nu }\varphi -\frac{1}{2}{g}_{\mu \nu }{\nabla }_{\alpha }\varphi {\nabla }^{\alpha }\varphi \right)\\ +\frac{1}{\varphi }\left[{\nabla }_{\mu }{\nabla }_{\nu }\varphi -g_{\mu \nu }(\square \varphi +\text{Λ})\right].\end{array}$$(7)

Here, T_µν is the total matter stress-energy tensor, while □ denotes the d’Alembertian. On the background, these equations give $$3H^2=\frac{8\pi }{c^4}\frac{\overline{\rho }}{\varphi }-3H\frac{\dot{\varphi }}{\varphi }+\frac{{\omega }_{\text{BD}}}{2}\frac{{\dot{\varphi }}^2}{{\varphi }^2}+\frac{\text{Λ}}{\varphi },$$(8) $$2\dot{H}+3H^2=-\frac{8\pi }{c^4}\frac{\overline{p}}{\varphi }-\frac{{\omega }_{\text{BD}}}{2}\frac{{\dot{\varphi }}^2}{{\varphi }^2}-2H\frac{\dot{\varphi }}{\varphi }-\frac{\ddot{\varphi }}{\varphi }+\frac{\Lambda }{\varphi },$$(9)

and the scalar field equation of motion reads $$\square \varphi =\frac{8\pi }{c^4}\left(\frac{T}{3+2{\omega }_{\text{BD}}}\right)-\frac{4\Lambda }{3+2{\omega }_{\text{BD}}},$$(10)

where T ≡ 𝑔^µνT_µν is the trace of the stress-energy tensor on the background.

The phenomenological QS functions in this theory read $$\begin{array}{l}\text{Σ}=\frac{1}{G_{\text{N}}\varphi },\hfill \\ \mu =\frac{4+2{\omega }_{\text{BD}}}{3+2{\omega }_{\text{BD}}}\text{Σ},\hfill \\ \eta \equiv \frac{\Phi }{\Psi }=\frac{1+{\omega }_{\text{BD}}}{2+{\omega }_{\text{BD}}}.\hfill \end{array}$$(11)

Constraints on JBD gravity have been obtained by using a combination of different cosmological data sets and sampling over a different parameterisation of the BD parameter ω_BD. For example, Avilez & Skordis (2014) imposed a flat prior on − log₁₀ (ω_BD) to obtain a lower bound of ω_BD > 1900 at 95% CL with CMB information from Planck 2013 data. Ballardini et al. (2016) obtained log₁₀ (1 + 1/ω_BD) < 0.0030 at 95% CL combining Planck 2015 and BOSS DR10-11 data; this upper bound was subsequently updated in Ballardini et al. (2020) with a combination of Planck 2018 and Baryon Oscillation Spectroscopic Survey (BOSS) DR12 data to log₁₀ (1 + 1/ω_BD) < 0.0022 at 95% CL. Joudaki et al. (2022) used a combination of the Planck 2018 CMB data, the 3 × 2pt combination of the Kilo-Degree Survey (KiDS) and the 2 degree Field gravitational Lens survey (2dFLens) data, the Pantheon supernovae data and BOSS measurements of the BAO, to find the coupling constant ω_BD > {1540, 160, 160, 350} at 95% CL for the different choices of parametrization (or priors): {log₁₀ (1/ω_BD), log₁₀ (1 + 1/ω_BD), 1/ω_BD, 1/log 10 ω_BD}^. These constraints were obtained by fixing the value of the scalar field today to $$G_{\text{N}}\varphi (a=1)=\frac{2{\omega }_{\text{BD}}+4}{2{\omega }_{\text{BD}}+3},$$(12)

in order to guarantee that the effective gravitational constant at present corresponds to the one measured in a Cavendish-like experiment (Boisseau et al. 2000) (see also Avilez & Skordis (2014); Joudaki et al. (2022); Ballardini et al. (2022) for studies of the JBD model without imposing this condition and Ballardini et al. (2016); Ooba et al. (2017); Rossi et al. (2019); Braglia et al. (2020, 2021); Cheng et al. (2021) for a simple generalisation of these constraints extending the JBD action Eq. (6) with different potentials and couplings to the Ricci scalar).

The JBD model is implemented in CLASSig (Umiltà et al. 2015) and hi_class (Zumalacárregui et al. 2017; Bellini et al. 2020). For these codes, the agreement and validation of the background and linear perturbations was thoroughly studied in Bellini et al. (2018). In this paper, we use the results produced with hi_class.

2.2 Dvali-Gabadadze-Porrati braneworld gravity

The DGP model (Dvali et al. 2000) is a five-dimensional braneworld model defined by the action $$S=\frac{c^4}{16\pi G_5}\underset{ℳ}{{\displaystyle \int }}{\text{d}}^{5}x\sqrt{-\gamma }{R}_5+\underset{\partial ℳ}{{\displaystyle \int }}{\text{d}}^{4}x\sqrt{-g}\left[\frac{c^4}{16\pi G_{\text{N}}}R+{ℒ}_{\text{m}}\right],$$(13)

where γ is the five-dimensional metric and R₅ its Ricci curvature scalar. G₅ and G_N are the five- and four-dimensional Newton constants, respectively. The matter Lagrangian is denoted with 𝓛_m and is confined to a four-dimensional brane in a five-dimensional Minkowski spacetime. The induced gravity described by the usual four-dimensional Einstein-Hilbert action is responsible for the recovery of the four-dimensional gravity on the brane. The cross-over scale r_c = G₅/(2G_N) is the only parameter of the model and GR is recovered in the limit r_c → ∞.

The Friedmann equation on the brane is given by (Deffayet 2001) $$H^2=\pm c\frac{H}{r_{\text{c}}}+\frac{8\pi G_N}{3}\overline{\rho }.$$(14)

Two branches of solutions depend on the embedding of the brane: the self-accelerating branch (sDGP, Bowcock et al. 2000; Deffayet 2001) and the normal branch (nDGP, Bowcock et al. 2000; Deffayet 2001), corresponding to the + and – sign for the contribution from the five-dimensional gravity, respectively. The self-accelerating branch admits the late-time acceleration without dark energy but the solution is plagued by ghost instabilities (Luty et al. 2003; Gorbunov et al. 2006; Charmousis et al. 2006). We therefore focus on the normal branch. In order to separate the effect of MG on structure formation from that on the expansion history, it is often assumed that the background expansion is identical to that of ΛCDM. This can be achieved by introducing an additional dark energy contribution with an appropriate equation of state (Schmidt 2009; Bag et al. 2018). We adopted this approach.

The evolution of density and metric perturbations on the brane require solutions of the bulk metric equations. These bulk effects can be encapsulated in an effective 3 + 1 description that uses the combination of any two of the functions µ, η, and Σ. Using the results of Koyama & Maartens (2006); Hu & Sawicki (2007a); Song (2008); Cardoso et al. (2008); Lombriser et al. (2009); Seahra & Hu (2010), we have $$\eta =\frac{1+g}{1-g},$$(15)

where, using the quasi-static (QS) approximation, we have (Lombriser et al. 2013) $$g(a)=g_{\text{QS}}=-\frac{1}{3}{\left[1\mp \frac{2H{r}_c}{c}\left(1+\frac{\dot{H}}{3H^2}\right)\right]}^{-1},$$(16)

so that the effective modification introduced with η can be treated as scale-independent.

The effective 3 + 1 Poisson equation for the lensing potential in the QS approximation is $$-k^2(\text{Φ}+\text{Ψ})=\frac{8\pi G_{\text{N}}}{c^2}{a}^2\overline{\rho }\text{Δ},$$(17)

therefore (Lombriser et al. 2013) $$\text{Σ}=1,$$(18)

and hence $$\mu (a)=1+\frac{1}{3\beta },\beta (a)\equiv 1+\frac{H}{H_0}\frac{1}{\sqrt{{\text{Ω}}_{\text{rc}}}}\left(1+\frac{\dot{H}}{3H^2}\right).$$(19)

We chose to parametrise the modification to gravity by ${\mathrm{\Omega }}_{\text{rc}}\equiv c^2/\left(4r_{\text{c}}^2{H}_0^2\right)$.

Not many studies have constrained nDGP with an exact ΛCDM background. However, Raccanelli et al. (2013), using measurements of the zeroth- and second-order moments of the correlation function from SDSS DR7 data up to r_max = 120 Mpc h⁻¹, and marginalised bias, found an Ω_rc upper limit at 95% of ~ 40 (from r_c > 340 Mpc with fixed H₀ = 70 km s⁻¹ Mpc⁻¹). In addition, Barreira et al. (2016) used the clustering wedges statistic of the galaxy correlation function and the growth rate values estimated from more recent BOSS DR12 data to set [r_cH₀/c]⁻¹ < 0.97 at 95% C.L., corresponding to an upper bound of Ω_rc < 0.27.

There have also been recent attempts to forecast constraints on Ω_rc. Liu et al. (2021), using the galaxy cluster abundance from the Simons Observatory and galaxy correlation functions from a Dark Energy Spectroscopic Instrument (DESI)-like experiment and found δ(Ω_rc) ~ 0.038 around a fiducial value of 0.25. Cataneo et al. (2022) forecast for Euclid-like future constraints on a fiducial Ω_rc ~ 0.0625 a 1σ accuracy of 0.0125 from combining the 3D matter power spectrum and the probability distribution function of the smoothed three-dimensional matter density field probes. Bose et al. (2020) also forecast for a Large Synoptic Survey Telescope (LSST)-like survey a 1σ accuracy of 0.08 using cosmic shear alone on a fiducial Ω_rc ~ 0.

Constraints have also been inferred for an nDGP model with a cosmological constant rather than a constructed dark energy field, thus with an approximate ΛCDM background. For this model, Lombriser et al. (2009) and Xu (2014) inferred somewhat stronger constraints of Ω_rc < 0.020 (95% C.L.) and Ω_rc < 0.002 (68% C.L.) from the combination of CMB and large-scale structure data, where the constraints were mainly driven by the CMB data. Thus, high-precision CMB measurements such as by the Planck satellite constrain Ω_rc tightly. In this work, we focus on the model with an exact ΛCDM background.

nDGP has been implemented in MGCAMB and in QSA_CLASS (Pace et al. 2021). These codes solve for a ΛCDM-background evolution and Eqs. (18) and (19). The overall agreement in the linear matter power spectrum is never worse than 0.5%. In this paper, we use the results produced with MGCAMB.

2.3 k-mouflage gravity

KM theories are built complementing simple k-essence scenarios with a universal coupling of the scalar field φ to matter. They are defined by the action (Babichev et al. 2009; Brax & Valageas 2014, 2016) $$S={{\displaystyle \int }}^{\text{}}{\text{d}}^{4}x\sqrt{-\tilde{g}}\left[\frac{c^4}{16\pi {\tilde{G}}_{\text{N}}}\tilde{R}+c^2{ℳ}^{4}K(\tilde{\chi })\right]+S_{\text{m}}\left({\psi }_i,g_{\mu \nu }\right),$$(20)

where ${\tilde{G}}_{\text{N}}$ is the bare Newton constant, $\tilde{R}$ is the Ricci scalar in the Einstein frame, ℳ⁴ is the energy density scale of the scalar field, 𝑔_µν is the Jordan frame metric, ${\tilde{g}}_{\mu \nu }$ is the Einstein frame metric with $g_{\mu \nu }=A^2(\phi ){\tilde{g}}_{\mu \nu },S_{\text{m}}$ is the Lagrangian of the matter fields ${\psi }_m^{(i)},\tilde{\chi }$ is defined as $$\tilde{\chi }=-\frac{{\tilde{g}}^{\mu \nu }{\partial }_{\mu }\phi {\partial }_{\nu }\phi }{2{ℳ}^4},$$(21)

and ℳ⁴K is the non-standard kinetic term of the scalar field.

In these theories, the evolution of both the background and perturbations is affected by the universal coupling and by the scalar field dynamics. The degree of deviation from ΛCDM at the background level and in perturbation theory can be expressed in terms of two time-dependent functions, related to the coupling A and the kinetic function K, that is, $${ϵ}_2\equiv \frac{\text{d\hspace{0.17em}}\ln \overline{A}}{\text{d\hspace{0.17em}}\ln a},{ϵ}_1\equiv \frac{2}{{\overline{K}}^{\prime }}{\left[{ϵ}_2\frac{c}{\sqrt{8\pi {\tilde{G}}_{\text{N}}}}{\left(\frac{\text{d}\overline{\phi }}{\text{d\hspace{0.17em}}\ln a}\right)}^{-1}\right]}^2,$$(22)

where over-bars denote background quantities and a prime indicates derivatives with respect to $\overline{\tilde{\chi }}$. The KM Friedmann equation therefore reads $$H^2=\frac{8\pi {\tilde{G}}_{\text{N}}}{3}\frac{{\overline{A}}^2}{{\left(1-{ϵ}_2\right)}^2}\left[\overline{\rho }+\frac{{ℳ}^4}{{\overline{A}}^4}\left(2\overline{\tilde{\chi }}\frac{\text{d}\overline{K}}{\text{d}\overline{\tilde{\chi }}}-\overline{K}\right)\right].$$(23)

Considering linear scalar perturbations around an FLRW background, the phenomenological functions µ and Σ read (Ben-evento et al. 2019) $$\mu (a)=\left(1+{ϵ}_1\right){\overline{A}}^2,\Sigma (a)={\overline{A}}^2.$$(24)

We recover GR when ϵ₁ → 0, ϵ₂ → 0, and A and K are constants, with ℳ⁴K playing the role of the cosmological constant, see Eq. (20).

In addition to the six standard ΛCDM parameters, the KM model requires the specification of the kinetic function $K(\tilde{\chi })$ and of the coupling A(φ). Following Brax & Valageas (2016), this can also be expressed in terms of the time-dependent background values $\overline{K}(a)$ and $\overline{A}(a)$. Brax & Valageas (2016) proposed a simple parameterisation that satisfies self-consistency constraints. It involves five parameters, {ϵ_2,0, γ_A,m,α_U,γ_U}, where ϵ_2,0 is the value of ϵ₂ at redshift ɀ = 0 and must be negative to ensure that there are no ghosts in the theory. For the coupling function $\overline{A}(a)$, this reads $$\overline{A}(a)=1+{\alpha }_A-{\alpha }_A{\left[\frac{\left({\gamma }_A+1\right)a}{{\gamma }_A+a}\right]}^{{\nu }_A},$$(25)

with $${\nu }_A=\frac{3(m-1)}{2m-1},{\alpha }_A=-\frac{{ϵ}_{2,0}\left(1+{\gamma }_A\right)}{{\gamma }_A{\nu }_A}.$$(26)

Here, m is the exponent of the kinetic function for a large argument, $K\sim {\tilde{\chi }}^m$, which also specifies the high-redshift dependence of $\overline{A}(a)$, while γ_A sets the redshift at which $\overline{A}$ changes from the current unit value to the high-redshift value 1 + α_A, which is parameterised by ϵ_2,0. The background kinetic function $\overline{K}(a)$ is conveniently parameterised in terms of a function U(a), with $$\frac{\text{d}{\overline{K}}^{\prime }}{\text{d}\tilde{\chi }}=\frac{U(a)}{a^3\sqrt{\tilde{\chi }}},\sqrt{\tilde{\chi }}=-\frac{{\overline{\rho }}_0}{{ℳ}^4}\frac{{ϵ}_2{\overline{A}}^4}{2U\left(-3{ϵ}_2+\frac{\text{d\hspace{0.17em}}\ln U}{\text{d\hspace{0.17em}}\ln a}\right)},$$(27)

which we took to be of the form $$U(a)\propto \frac{a^2\ln \left({\gamma }_U+a\right)}{\left(\sqrt{a_{\text{eq}}}+\sqrt{a}\right)\ln \left({\gamma }_U+a\right)+{\alpha }_U{a}^2}.$$(28)

The various terms in this expression allowed us to follow the radiation, matter, and dark-energy eras. The parameters α_U and γ_U, of order unity, set the shape of the transition to the dark-energy era. The main parameter is ϵ_2,0, which measures the amplitude of the deviation from GR and ΛCDM at ɀ = 0. Other parameters, of order unity, mostly describe the shapes of the transitions between different cosmological regimes.

The KM model has been tested against a set of complementary cosmological data sets, including CMB temperature, polarisation and lensing, type Ia supernovae, baryon acoustic oscillations (BAO) and local measurements of H₀ (Benevento et al. 2019). This gives the 95% C.L. bounds –0.04 ≤ ϵ_2,0 ≤ 0, while the other parameters are unconstrained (as long as they do not become very large). These results are consistent with earlier more qualitative CMB studies (Barreira et al. 2015). The X-ray cluster multiplicity function could provide bounds of the same order (Brax et al. 2015), but more detailed analyses are needed to derive robust constraints.

The fact that ϵ_2,0 is the main parameter constrained by the data can be understood from Eqs. (22) and (24) and the property that at leading order over ϵ₂ we have ϵ₁ ≃ –ϵ₂. This shows that the running of Newton’s constant and the impact on the linear metric and matter density fluctuations are set by ϵ₂. Because most of the deviations from ΛCDM occur at low redshift, this is mostly determined by the value ϵ_2,0 at ɀ = 0. Therefore, we here focused on the dependence on ϵ_2,0, which was left free to vary, while we fixed the other parameters to γ_A = 0.2, m = 3, α_U = 0.2, and γ_U = 1, as some representative values for the other parameters. Thus, this gives one additional parameter, ϵ_2,0, in addition to the six standard ΛCDM parameters. The ΛCDM model is recovered when ϵ_2,0 = 0, independently of the value of the other parameters; for ϵ_2,0 → 0, we have: $$\overline{A}(a)\to 1,{ϵ}_2(a)\to 0,{ϵ}_1(a)\to 0,\forall a,$$(29)

and the kinetic function in Eq. (20) reduces to a cosmological constant.

To produce linear predictions for this theory we used its implementation in the EFTCAMB code that is described in Benevento et al. (2019).

3 Theoretical predictions for Euclid observables

Both the forecasting method and the tools used in this paper are the same as those in EC19 except for the changes needed to account for the use of MG rather than GR in the predictions of Euclid observables. The motivations and the relevant changes in the recipe have been described in a previous paper on forecasts for f (R) theories. We therefore refer to Casas et al. (2023), while here we only recall the main steps.

3.1 Photometric survey

The observables we considered for the Euclid photometric survey are the angular power spectra $C_{ij}^{XY}(\ell )$ between the probe X in the i - th redshift bin, and the probe Y in the j - th bin, where X refers to weak lensing (WL), photometric galaxy clustering (GC_ph), or their cross correlation, XC. These are still calculated as in EC19 relying on the Limber approximation and setting to unity the ℓ - dependent factor in the flat-sky limit. The spectra are then given by $$C_{ij}^{XY}(\ell )=c{{{\displaystyle \int }}^{\text{}}}_{{ɀ}_{\min }}^{{ɀ}_{\max }}\text{d}ɀ\frac{W_i^X(ɀ)W_j^Y(ɀ)}{H(ɀ)r^2(ɀ)}{P}_{\delta \delta }\left(k_{\ell },ɀ\right),$$(30)

with k_ℓ = (ℓ + 1/2)/r(z), r(z) the comoving distance to redshift z = 1/a – 1, and P_δδ (k_{ℓ, z}) the non-linear power spectrum of matter density fluctuations, δ, at wave number k_ℓ and redshift ɀ. We set (ɀ_min, ɀ_max) = (0.001, 4), which spans the full range within which the source redshift distributions n_i(ɀ) are non-vanishing. The GC_ph and WL window functions read (Spurio Mancini et al. 2019) $$W_i^{\text{G}}(k,ɀ)=b_i(k,ɀ)\frac{n_i(ɀ)}{{\overline{n}}_i}\frac{H(ɀ)}{c},$$(31) $$\begin{array}{l}{W}_i^{\text{L}}(k,ɀ)=\frac{3}{2}{\text{Ω}}_{\text{m},0}\frac{H_0^2}{c^2}(1+ɀ)r(ɀ)\text{Σ}(ɀ){{{\displaystyle \int }}^{\text{}}}_{ɀ}^{ɀ\max }\text{d}{ɀ}^{\prime }\frac{n_i\left({ɀ}^{\prime }\right)}{{\overline{n}}_i}\frac{r\left({ɀ}^{\prime }-ɀ\right)}{r\left({ɀ}^{\prime }\right)}\\ +W_i^{\text{IA}}(k,ɀ),\end{array}$$(32)

where $n_i(ɀ)/{\overline{n}}_i$ and b_i(k, ɀ) are the normalised galaxy distribution and the galaxy bias in the i-th redshift bin, respectively, and $W_i^{\text{IA}}(k,ɀ)$ encodes the contribution of intrinsic alignments (IA) to the WL power spectrum. The function Σ(z) in the WL window function explicitly accounts for the changes in the lensing potential due to the particular MG theory of interest. Its explicit form for the cases under consideration can be found in Section 2. The impact on the background quantities H(ɀ), r(ɀ) and the matter power spectrum P_δδ(k, ɀ), in contrast, are already taken into account by the dedicated Boltzmann solver so that the GC_ph window function remains unchanged.

The IA contribution was computed following the eNLA model adopted in EC19 so that the corresponding window function is $$W_i^{\text{IA}}(k,ɀ)=-\frac{{\mathcal{A}}_{\text{IA}}{C}_{\text{IA}}{\text{Ω}}_{\text{m},0}{ℱ}_{\text{IA}}(ɀ)}{\delta (k,ɀ)/\delta (k,ɀ=0)}\frac{n_i(ɀ)}{{\overline{n}}_i(ɀ)}\frac{H(ɀ)}{c},$$(33)

where $${ℱ}_{\text{IA}}(ɀ)={(1+ɀ)}^{{\eta }_{\text{IA}}}{\left[\frac{\langle L\rangle (ɀ)}{L_{\star }(ɀ)}\right]}^{{\beta }_{\text{IA}}},$$(34)

with 〈L〉(z) and L_*(ɀ) redshift-dependent mean and the characteristic luminosity of source galaxies as computed from the luminosity function, 𝒜_IA, β_IA and η_IA are the nuisance parameters of the model, and 𝒞_IA is a constant accounting for dimensional units. This model is the same as was used in EC19 since IA takes place on astrophysical scales that are unaffected by modifications to gravity. However, MG has an impact on the growth factor, introducing a possible scale dependence. This is explicitly taken into account in Eq. (33) through the matter perturbation δ(k, ɀ), which is considered to be scale dependent in this case. This allows us to also consider the scale dependence introduced by massive neutrinos, which was assumed to be negligible in EC19. We nevertheless stress that for the models we considered, the scale dependence is quite small so that the IA is essentially the same as in the GR case.

3.2 Spectroscopic survey

We now discuss the modelling of the power spectrum to analyse the data from the Euclid spectroscopic survey.

For the models considered in this paper, the Compton wavelength of the scalar field is much larger than the scales probed by Euclid. Moreover, it is assumed that the speed of propagation of the scalar fluctuations is of the order of the speed of light, so that the sound horizon is of the order of the Hubble scale. Under these assumptions, we can apply the quasi-static approximation and relate the scalar field perturbation to the gravitational potential. Since in all these models the weak equivalence principle holds, the modelling of the bias as an expansion in the derivatives of the gravitational potential remains unchanged with respect to the ΛCDM one. To be consistent with the official forecast, we used the same modelling for galaxy clustering as in EC19.

The observed galaxy power spectrum is given by $$\begin{array}{cc}{P}_{\text{obs}}\left(k,{\mu }_{\theta };ɀ\right)=\hfill & \frac{1}{q_{\perp }^2(ɀ)q_{\parallel }(ɀ)}\left\{\frac{{\left[b{\sigma }_8(ɀ)+f{\sigma }_8(ɀ){\mu }_{\theta }^2\right]}^2}{1+{\left[f(ɀ)k{\mu }_{\theta }{\sigma }_{\text{p}}(ɀ)\right]}^2}\right\}\hfill \\ & \times \frac{P_{\text{dw}}\left(k,{\mu }_{\theta };ɀ\right)}{{\sigma }_8^2(ɀ)}{F}_{ɀ}\left(k,{\mu }_{\theta };ɀ\right)+P_{\text{s}}(ɀ),\hfill \end{array}$$(35)

where P_dw(k,µ_θ; z), the de-wiggled power spectrum, includes the correction that accounts for the smearing of the BAO features, $$P_{\text{dw}}\left(k,{\mu }_{\theta };ɀ\right)=P_{\delta \delta }^{\mathrm{lin}}(k;ɀ){\text{e}}^{-g_{\mu }{k}^2}+P_{\text{nw}}(k;ɀ)\left(1-{\text{e}}^{-g_{\mu }{k}^2}\right),$$(36)

and $P_{\delta \delta }^{\text{lin}}(k;ɀ)$ stands for the linear matter power spectrum. P_nw(k; ɀ) is a no-wiggle power spectrum with the same broadband shape as $P_{\delta \delta }^{\text{lin}}(k;ɀ)$ but without BAO features (see the discussion below). The function $$g_{\mu }\left(k,{\mu }_{\theta },ɀ\right)={\sigma }_{\text{v}}^2(ɀ)\left\{1-{\mu }_{\theta }^2+{\mu }_{\theta }^2{[1+f(ɀ)]}^2\right\},$$(37)

is the non-linear damping factor of the BAO signal derived by Eisenstein et al. (2007), with $${\sigma }_{\text{v}}^2(ɀ)=\frac{1}{6{\pi }^2}{\displaystyle \int \text{d}}k{P}_{\delta \delta }^{\mathrm{lin}}(k,ɀ).$$(38)

The curly bracket in Eq. (35) is the redshift-space-distortion (RSD) contribution correcting for the non-linear finger-of-God (FoG) effect, where we defined bσ₈(ɀ) as the product of the effective scale-independent bias of galaxy samples and the r.m.s. matter density fluctuation σ₈(ɀ) (we marginalized over bσ₈(ɀ)), while µ_θ is the cosine of the angle θ between the wave vector k and the line-of-sight direction $\widehat{r}$ and ${\sigma }_{\text{p}}^2(ɀ)={\sigma }_{\text{v}}^2(ɀ)$. Although these parameters were assumed to be the same, they come from two different physical effects, namely large-scale bulk flow for the former and virial motion for the latter.

The factor F_z accounts for the smearing of the galaxy density field along the line of sight due to redshift uncertainties. It is given by $$F_{ɀ}\left(k,{\mu }_{\theta };ɀ\right)={\text{e}}^{-k^2{\mu }_{\theta }^2{\sigma }_r^2(ɀ)},$$(39)

where σ_r(ɀ) = c (1 + ɀ)σ_0,z/H(ɀ) and σ_0,z is the error on the measured redshifts.

The factor in front of the curly bracket in Eq. (35) describes the Alcock-Paczynski effect, which is parametrised in terms of the angular diameter distance D_A(ɀ) and the Hubble parameter H(ɀ) as $$q_{\perp }(ɀ)=\frac{D_{\text{A}}(ɀ)}{D_{\text{A},\mathrm{ref}}(ɀ)},$$(40) $$q_{\Vert }(ɀ)=\frac{H_{\text{ref}}(ɀ)}{H(ɀ)}.$$(41)

Due to the Alcock-Paczynski effect, µ_θ and k were rescaled in a cosmology-dependent way as a function of the projection along and perpendicular to the line of sight. In the previous Eq. (35), the arguments µ and k are themselves functions of the true µ_θ,ref and k_ref at the reference cosmology. This relation, which for each argument, is a function of q|| and q_⊥, can be found in Sect. 3.2.1 of EC19. Finally, the P_s(ɀ) is a scale-independent shot-noise term that enters as a nuisance parameter (see EC19).

Due to the scale independence of σ_v(= σ_p), we evaluated it at each redshift bin but we kept it fixed in the Fisher matrix analysis. This method corresponds to the optimistic settings in EC19. We would like to highlight that we directly took the derivatives of the observed galaxy power spectrum with respect to the final parameters, in contrast to EC19 where first a Fisher matrix analysis was performed for the redshift-dependent parameters H(ɀ), D_A(ɀ), and fσ₈(ɀ) and then projected to the final cosmological parameters of interest. However, we verified that both approaches lead to consistent results when considering the ΛCDM and w₀w_aCDM models. The other term appearing in Eq. (36) is the non-wiggle matter power spectrum, and this was obtained applying a Savitzky-Golay filter to the matter power spectrum $P_{\delta \delta }^{\text{lin}}(k;ɀ)$. For more details on this implementation, we refer to EC19.

3.3 non-linear modelling for weak lensing

The galaxy power spectrum on mildly non-linear scales can be modelled by a modified version of the implementation in EC19. This is no longer the case when we move to the deeply non-linear regime in which no an analytical description of the matter power spectrum is available for the models of our interest.

The wide window functions (in particular, the WL function) entering the prediction of the photometric observables $C_{ij}^{XY}(\ell )$ require an accurate description of non-linear scales. In particular, the impact of baryons needs to be explicitely accounted for and nuisance parameters must be included, which control the baryonic feedback in order to avoid biasing the parameter estimation (see e.g. Schneider et al. 2020b,a). However, we currently do not have accurate Euclid-like simulations including baryonic effects, especially for MG cosmologies. Therefore, we ignored these effects in our analysis, and left their inclusion for a future work. We note that our scale cuts, in particular for the pessimistic scenario (ℓ < 1500), partially mitigate the fact of neglecting baryonic effects and are in agreement with recent studies on the impact of baryonic effects on constraints for modified gravity (see, for instance Spurio Mancini & Bose 2023). In the following subsections, we comment on the individual prescriptions for each MG model, wthat is used to compute the non-linear matter power spectrum. While having a common prescription at non-linear scales for all models would be desirable, given the different behaviours of each model at non-linear scales, this remains an ambitious goal. Currently, we have at our disposal codes implementing the specific features of each model that have already been tested and used in the literature. This reinforces the validity of our procedure.

3.3.1 Jordan-Brans-Dicke gravity

For the non-linear prescription, we used a modified version of HMCODE (Mead et al. 2015, 2016). HMCODE is an augmented halo model that can be used to accurately predict the non-linear matter power spectrum over a wide range of cosmologies. A brief summary of how this works is given below.

In the halo model (Cooray & Sheth 2002), the non-linear matter power spectrum can be written as a sum P_NL(k, ɀ) = P_1H(k, ɀ) + P_2H(k, ɀ) where $$P_{2\text{H}}(k,ɀ)=P_{\text{lin\hspace{0.17em}}}(k,ɀ){\left[{\displaystyle {\int }_0^{\infty }b}(M,ɀ)W(M,k,ɀ)n(M,ɀ)\text{d}M\right]}^2,$$(42)

is the so-called two-halo term (correlation between different haloes), and $$P_{1\text{H}}(k,ɀ)={\displaystyle {\int }_0^{\infty }{W}^2}(M,k,ɀ)n(M,ɀ)\text{d}M,$$(43)

is the so-called one-halo term (correlations between mass-elements within each halo). Above M is the halo mass, P_lin(k, ɀ) is the linear matter power spectrum, n(M, ɀ) is the halo mass function and b(M, ɀ) is the linear halo bias.

The window function W is the Fourier transform of the halo matter density profiles: $$W(M,k,ɀ)=\frac{1}{\overline{\rho }}{\displaystyle {\int }_0^{r_{\text{v}}}4}\pi r^2\frac{\sin (kr)}{kr}\rho (M,r,ɀ)\text{d}r,$$(44)

where ρ(M, r, ɀ) is the radial matter density profile in a host halo with a mass M, and $\overline{\rho }$ is the mean matter density. The halo mass is related to the virial radius, r_v, via $M=4\pi r_{\text{v}}^3{\Delta }_v(ɀ)\overline{\rho }/3$, where Δ_v(ɀ) is the virial halo overdensity. The halo profiles, the halo definition, and the halo mass function can either be computed from excursion set models, extracted from numerical simulations and/or parametrised as functions with free parameters that are then fit to data. HMCODE combines this, and the resulting fitting function has, for ΛCDM, been shown to be accurate to 2.5% for scales k < 10 h Mpc^–1 and redshifts z < 2 (Mead et al. 2021).

Joudaki et al. (2022) modifed the HMC0DE, which is able to include the effects of JBD. This was done using a suite of N-body simulations obtained with modified versions of both the COmoving Lagrangian Acceleration (COLA, Tassev et al. 2013; Winther et al. 2017) and RAMSES (Teyssier 2002) codes. COLA solves for the perturbations around paths predicted from second-order Lagrangian perturbation theory (2LPT), and it has proven to be fast and accurate on large scales. This was used for scales k < 0.5 h Mpc^–1 to generate a large enough ensemble of simulations and to reduce sample variance. On very small-scales, up to k < 10 h Mpc^–1, the RAMSES grid-based hydrodynamical solver with adaptive mesh refinement was used.

The spectra generated by these simulations were then used to calibrate HMC0DE. While we did not consider the effects of baryons here, the advantage of HMC0DE, over Halofit is that the former is able to capture baryonic feedback at an accuracy level of ~5% up to k ≃ 10 h Mpc^–1 (Mead et al. 2021). This will become an important consideration (Schneider et al. 2020b,a) in future more detailed analyses. To take JBD into account in addition to modifying the expansion history and growth of perturbations, the virialised halo overdensity was modified as $$\begin{array}{l}{\Delta }_{\text{v}}={\mathrm{\Omega }}_{\text{m}}{(ɀ)}^{-0.352}\\ \times \left\{d_0+\left(418.0-d_0\right)\arctan \left[{\left(0.001\left|{\omega }_{\text{BD}}-50.0\right|\right)}^{0.2}\right]\frac{2}{\pi }\right\},\end{array}$$(45)

where d₀ = 320.0 + 40.0z^0.26. This modification has the feature that it reduces to the original ΛCDM HMC0DE as ω_BD → ∞ We stress that this modification is not a physical claim about JBD causing this particular change in the virialised halo overdensity, but rather that this change in the virialised halo overdensity within the HMC0DE machinery is able to accurately reproduce the JBD power spectrum to the given accuracy. The resulting power spectrum with this modification was found to be accurate to 10% for the fitted range 10⁴ ≳ ω_BD ≥ 50 and on scales up to k ≃ 10 h Mpc^–1. This non-linear prescription was made with past weak-lensing surveys in mind and fitted to a range of ω_BD and scales that we can realistically constrain at the present time. For a parameter inference with actual Euclid data the accuracy here might not be good enough.

We show in Fig. 1 the matter power spectrum and the lensing angular power spectrum of JBD with ω_BD = 800, which we refer to as JBD1, and their comparison to ΛCDM. In the top left panel we plot the linear (dashed) and non-linear (solid) matter power spectrum as a function of scale for ɀ = 0. In the top right panel we show the ratio of these power spectra with respect to their corresponding ΛCDM cases, using the same σ₈ normalisation today.

As was shown in detail in Joudaki et al. (2022), the linear growth rate will undergo a scale-independent enhancement (without altering the shape of the linear power spectrum) while the non-linear growth will be mildly suppressed on smaller scales due to the presence of the scalar field. This is clear for JBD1, where the linear power spectrum is almost identical to its standard model counterpart, while the non-linear power spectrum shows a small suppression on small-scales that sets in at about k = 0.3 h Mpc⁻¹. It is important to note that, unlike nDGP, the JBD theory is not endowed with gravitational screening, so that any modification on small-scales in the non-linear power spectrum is inherited from the change in the primordial amplitude, A_s, in the linear power spectrum. In order for σ₈ to be the same in ΛCDM as in JBD we need a higher primordial amplitude As . A higher A_s in ΛCDM means that non-linear structures form faster and become more massive, which in the matter power spectrum then translates into an even higher amplitude on nonlinear scales (i.e. it increases more than linear theory predicts) which gives rise to the small-scale suppression seen in Fig. 1. This effect can also be found in pure ΛCDM when we consider the ratio of the matter power spectrum for two different values of A_s . The result is then constant on large scales, but shows an enhancement or suppression on small-scales. When we instead show results with the same primordial amplitude in JBD as in ΛCDM then the JBD power spectrum would show a certain enhancement to ΛCDM on linear scales and an ever stronger enhancement on non-linear scales, which is more in line with naive expectations of a larger gravitational constant leading to an enhancement in the matter power spectrum. The signature in the matter power spectrum naturally impacts the lensing spectrum, which shows a similar suppression at small scales because for a given multipole and redshift bin, it is just proportional to the matter power spectrum as defined in Eq. (30).

Fig. 1 Large-scale structure observables for the different MG models evaluated at the fiducial values of the model parameters used for the Fisher matrix analysis. Top left: linear (dashed) and non-linear (solid lines) matter power spectrum Pδδ(z, k) entering Eqs. (36) and (30), respectively, evaluated at redshift ζ = 0, for KM1 (ϵ2,0 = –0.04, blue), JBD1 (ωBD = 800, red) and nDGPl (Ωrc = 0.25, green). Top right: ratio of the matter power spectra for the above mentioned models with respect to their ΛCDM counterpart with the same value of σ8 today, for the linear (dashed) and non-linear (solid) cases. Bottom left: cosmic shear (WL) angular power spectra for the auto-correlation of the first photometric bin $C_{1,1}^{\text{WL}}(\ell )$ (solid line) and the cross-correlation of the first and last bin $C_{1,10}^{\text{WL}}(\ell )$ (dotted line), defined in Eq. (30) for the three models considered above, with the same colour labelling. Bottom right: ratio of the cosmic shear Ci,j(ℓ) to their ΛCDM counterpart, for the bin combination i = j = 1 and i = 1; j = 10.

3.3.2 Dvali-Gabadadze-Porrati gravity

We modelled the non-linear power spectrum using the halo model reaction (Cataneo et al. 2019) which has been shown to agree with N-body simulations at the 2% level down to k = 3h Mpc⁻¹ with small variation depending on redshift, degree of modification to GR and mass of neutrinos (Bose et al. 2021). The approach attempts to model non-linear corrections to the power spectrum from modified gravity through the so-called reaction 𝓡(k, ɀ), which employs both one-loop perturbation theory and the halo model. The non-linear power spectrum in nDGP is given by the product $$P_{\text{NL}}(k,ɀ)=ℛ(k,ɀ)P_{\text{NL}}^{\text{pseudo}}(k,ɀ),$$(46)

where the pseudo-power spectrum is a spectrum in which all non-linear physics are modelled using GR but the initial conditions are tuned in such a way as to replicate the modified linear clustering at the target redshift. We used the halofit formula of Takahashi et al. (2012) to model $P_{\text{NL}}^{\text{pseudo\hspace{0.17em}}}$(k, ɀ) by providing the halofit formula a linear matter power spectrum modelled within nDGP gravity as input. This ensures that we have the correct modified linear clustering at ɀ while keeping a non-linear clustering as described by GR, in line with our definition of the pseudo-power spectrum.

The halo model reaction, 𝓡(k, ɀ), is given by a corrected ratio of target to pseudo-halo model spectra, $$ℛ(k,ɀ)=\frac{\left\{[1-ℰ(ɀ)]e^{-k/k_{\star }(ɀ)}+ℰ(ɀ)\right\}{P}_{2\text{H}}(k,ɀ)+P_{1\text{H}}(k,ɀ)}{P_{\text{hm}}^{\text{pesudo\hspace{0.17em}}}(k,ɀ)}.$$(47)

The components are given explicitly as $$P_{\text{hm}}^{\text{pseudo\hspace{0.17em}}}(k,ɀ)=P_{2\text{H}}(k,ɀ)+P_{1\text{H}}^{\text{pseudo\hspace{0.17em}}}(k,ɀ),$$(48) $$ℰ(ɀ)=\underset{k\to 0}{\lim }\frac{P_{1\text{H}}(k,ɀ)}{P_{1\text{H}}^{\text{pseudo\hspace{0.17em}}}(k,ɀ)},$$(49) $$k_{\star }(ɀ)=-\overline{k}{\left\{\ln \left[\frac{A(\overline{k},ɀ)}{P_{2\text{H}}(\overline{k},ɀ)}-ℰ(ɀ)\right]-\ln [1-ℰ(ɀ)]\right\}}^{-1},$$(50)

where $$A(k,ɀ)=\frac{P_{1\text{−loop\hspace{0.17em}}}(k,ɀ)+P_{1\text{H}}(k,ɀ)}{P_{1\text{−loop\hspace{0.17em}}}^{\text{pseudo\hspace{0.17em}}}(k,ɀ)+P_{1\text{H}}^{\text{pseudo\hspace{0.17em}}}(k,ɀ)}{P}_{\text{hm}}^{\text{pseudo\hspace{0.17em}}}(k,ɀ)-P_{1\text{H}}(k,ɀ)\text{.\hspace{0.17em}}$$(51)

P_շH(k, ɀ) is the two-halo contribution, which we can approximate with the (nDGP) linear power spectrum, P_L (k, ɀ) (see for a review on the halo model Cooray & Sheth 2002). P_1H(k, ɀ) and $P_{1\text{H}}^{\text{pseudo\hspace{0.17em}}}$(k, ɀ) are the one-halo contributions to the power spectrum as predicted by the halo model, with and without modifications to the standard ΛCDM spherical collapse equations, respectively. We recall that by definition, the pseudo-cosmology has no non-linear beyond-ΛCDM modifications. Similarly, P _1–loop(k, ɀ) and $P_{1-\text{\hspace{0.17em}loop\hspace{0.17em}}}^{\text{pseudo\hspace{0.17em}}}$ (k, ɀ) are the one-loop predictions with and without non-linear modifications to ΛCDM, respectively. As in previous works, the limit in Eq. (49) was taken to be at k = 0.01 h Mpc⁻¹, and we computed k_⋆ using $\overline{k}=0.06h{\text{Mpc}}^{-1}$. This scale was chosen such that the one-loop predictions are sufficiently accurate at all the redshifts considered.

We note that the correction to the halo-model ratio in Eq. (47) has been shown to improve this ratio when there are modifications to gravity that invoke some sort of screening mechanism (Cataneo et al. 2019).

As in previous works, in our halo-model calculations, we used a Sheth-Tormen mass function (Sheth & Tormen 1999, 2002), a power-law concentration-mass relation (see for example Bullock et al. 2001) and an NFW halo-density profile (Navarro et al. 1997). These calculations were performed using the publicly available ReACT code (Bose et al. 2020). We refer to this reference for further computational and theoretical details.

Finally, as previously discussed in the JBD case, we show in Fig. 1 the comparison of the matter power spectrum and the lensing angular power spectrum of nDGP versus ΛCDM as a function of scale for ɀ = 0. We chose for Ω_rc the value Ω_rc = 0.25, which is the largest nDGP modification to ΛCDM used in this work. We refer to this choice as the nDGP1 model. The ratio of these power spectra with respect to their corresponding ΛCDM cases was determined using the same σ₈ normalisation today. For nDGP1, the linear matter power spectrum is identical to that of the ΛCDM. However, the non-linear matter power spectrum shows a suppression starting from scales k = 0.1 h Mpc⁻¹ , which reaches about 10% at scales of about k = 1.0 h Mpc⁻¹. This impacts the lensing spectrum, which shows a similar suppression at small angular scales because for a given multipole and redshift bin, it is just proportional to the matter power spectrum, as defined in Eq. (30).

As the σ₈ normalisation is the same for the ΛCDM reference spectrum, it is essentially the $P_{\text{NL}}^{\text{pseudo\hspace{0.17em}}}$ by definition. This means the quantity associated with the solid green curve in the top right panel of Fig. 1 is approximately 𝓡(k, ɀ = 0). This quantity is largely governed by the ratio of the nDGP to pseudo one-halo terms at k > 0.1 h Mpc⁻¹ (Cataneo et al. 2019). In a nDGP cosmology, there will be more high-mass halos when compared to its GR counterpart due to the supplemental fifth force sourced by the additional degree of freedom. High-mass contributions to the one-halo term are relatively more suppressed by the NFW density profile (see for example Fig. 9 of Cooray & Sheth 2002, and Eq. (52)). This will cause the 1-halo GR power spectrum to be larger than the DGP one-halo spectrum if the linear clustering in the two cosmologies is the same at the target redshift.

We note that when we instead fix the same primordial amplitude of perturbations, we obtain the reverse: more halos will have formed by the target redshift in the nDGP cosmology, which means an enhancement of the one-halo term over GR. We refer to Figs. 6 and 8 in Cataneo et al. (2019), which highlight these two scenarios, GR and nDGP, with the same late-time amplitude of linear perturbations and one with the same early-time amplitude, respectively.

3.3.3 k-mouflage gravity

For the KM non-linear matter power spectrum, we used a similar analytical approach that combined one-loop perturbation theory with a halo-model. This method was introduced in Valageas et al. (2013) for the standard ΛCDM cosmology. As in usual halo models, it splits the matter power spectrum over two- halo and one-halo contributions, as in Eq. (48), but it uses a Lagrangian framework to evaluate these two contributions in terms of the pair-separation probability distribution. This provides a Lagrangian-space regularisation of perturbation theory, which matches standard perturbation theory up to one-loop order and includes a partial resummation of higher-order terms, such that the pair-separation probability distribution is positive and normalised to unity at all scales. Thus, the two-halo term goes beyond the Zeldovich approximation by including a non-zero skewness, which enables the consistency with standard perturbation theory up to one-loop order, as well as a simple Ansatz for higher-order cumulants inspired by the adhesion model. The one-halo term includes a counter-term that ensures its falloff at low k, in agreement with the conservation of matter and momentum, $$P_{1\text{H}}(k)={\displaystyle {\int }_0^{\infty }\frac{\text{d}\nu }{\nu }}f(\nu )\frac{M}{\overline{\rho }{(2\pi )}^3}{\left[{\tilde{u}}_M(k)-\tilde{W}\left(k{q}_M\right)\right]}^2,$$(52)

where ${\tilde{u}}_M$ (k) is the normalised Fourier transform of the halo radial profile, $\tilde{W}\left(k{q}_M\right)$ is the normalised Fourier transform of the top-hat of Lagrangian radius q_M, and f (ν) is the normalised halo mass function. This counter-term is not introduced by hand, as it arises directly from the splitting of the matter power spectrum over two-halo and one-halo contributions within a Lagrangian framework, which by construction satisfies the conservation of matter. Additionally, ν is defined as δ_c/σ(M), where δ_c is the linear density contrast associated with a non-linear density contrast of 200, and σ(M) is the root mean square of the linear density contrast at mass scale M. The threshold δ_c is sensitive to the modification of gravity. In principle, it also depends on M, as the shells in the spherical collapse are coupled and evolve differently depending on their masses due to non-linear screening. However, for the models that we considered here, screening is negligible beyond 1 h⁻¹ Mpc (clusters are not screened), so that δ_c is independent of M.

This non-linear modelling was extended in Brax & Valageas (2013) to several modified-gravity scenarios [ f (R) theories, dilaton, and symmetron models], and in Brax & Valageas (2014); Brax et al. (2015) to KM models. This involved the computation of the impact of modified gravity on both the linear and one-loop contributions to the matter power spectrum, which enter the two- halo term, and on the spherical collapse dynamics, which enter the one-halo term through the halo mass function. Therefore, this non-linear modelling exactly captures the modification of gravity both up to one-loop order and at the level of the fully non-linear spherical collapse. The Ansatz for the partial resummation of higher-order terms is the same as for the standard ΛCDM cosmology (this implies that their values are also modified as they depend on the lower orders to ensure the positivity and the normalisation of the underlying pair distribution). A comparison with numerical simulations (Brax & Valageas 2013) for f(R) theories (with $\left|f_{R_0}\right|={10}^{-4}$, 10⁻⁵ and 10⁻⁶), shows that this approach captures the relative deviation from the ΛCDM power spectrum at ɀ = 0 up to k = 3 h Mpc⁻¹. We have not yet compared this recipe with numerical simulations of the KM theories, which have only recently been performed (Hernández-Aguayo et al. 2022). However, we expect at least the same level of agreement because the KM models are simpler and closer to the ΛCDM cosmology than the f (R) models (as in ΛCDM, the linear growth rate is scale independent).

Fig. 1 shows the comparison between the KM1 model defined with є_2,0 = −0.04 and its ΛCDM counterpart, normalised to the same σ₈ and cosmological parameters at ɀ = 0, for the matter and cosmic shear power spectra. The ratio of the linear matter power spectra increases on intermediate scales (10⁻² < k < 1 h Mpc⁻¹): This is due to the change of the background evolution, which slightly affects the scales at matterradiation equality. H(ɀ) is indeed greater for KM1 at high redshift (Benevento et al. 2019), which implies a smaller Hubble radius at matter-radiation equality. This in turns means that the change in the slope of the linear power spectrum, from $P_{\text{L}}(k)~k^{n_{\text{s}}}$ to $P_{\text{L}}(k)\propto k^{n_{\text{s}}-4}$ , is shifted to smaller scales, and hence to higher k. This implies that the ratio $P_{\text{KM}1}^{\text{mm}}/P_{\Lambda \text{CDM}}^{\text{mm}}$grows on these intermediate scales, where the linear power spectrum changes slope, because of this delay to higher k of the transition. Because the two spectra have the same σ₈ and the ratio $P_{\text{KM}1}^{\text{mm}}/P_{\Lambda \text{CDM}}^{\text{mm}}$ grows with k, the ratio is below unity at low k and above unity at high k. Moreover the linear growth rates are scale independent for both KM and ΛCDM, and therefore the ratio $P_{\text{KM}1}^{\text{mm}}/P_{\Lambda \text{CDM}}^{\text{mm}}$ reaches a constant at very low k and at very high k, where the slopes of both linear power spectra are equal to n_s or n_s − 4 (i.e. outside of the intermediate scales where the slope of the linear power spectrum slowly runs with k). This flat ratio is clearly shown in Fig. 1 at low k. The amplitude of $P_{\text{KM}1}^{\text{mm}}/P_{\Lambda \text{CDM}}^{\text{mm}}$ is set by the parameter є_2,0, and it deviates substantially from 1, by about 30%, for a value of є_2,0 = −0.04. The increase in $P_{\text{KM}1}^{\text{mm}}/P_{\Lambda \text{CDM}}^{\text{mm}}$ with k is partly due to non-linear effects in the scalar field sector, both in the background and in the linear growth rate, which are not completely irrelevant. At ɀ = 0, even though є₂ = −0.04 we have є₁ ≃ 0.12, whereas at linear order in the scalar field dynamics we would have expected that є₁ = −є₂. As the growth of structures is sensitive to є₁ , this leads to a stronger effect than could be expected from the value of є₂ itself. Then, the power spectrum being a quadratic functional of the density field gives an additional amplification factor of 2 for relative perturbations.

Beyond this global behaviour of the ratio $P_{\text{KM}1}^{\text{mm}}/P_{\Lambda \text{CDM}}^{\text{mm}}$, another impact of the slight modification of the background is that the baryon acoustic oscillations are shifted. This leads to the oscillations seen at BAO scales, k ~ 0.01 h Mpc⁻¹. Non-linear effects lead to a decrease in the ratio $P_{\text{KM}1}^{\text{mm}}/P_{\Lambda \text{CDM}}^{\text{mm}}$ as compared with the linear prediction of mildly non-linear scales. This is due to the loss of information in the power spectrum in the non-linear regime, which can be related to the universal NFW profile of virialised haloes within the halo model. Then, the shape of the non-linear power spectrum at those scales is mostly set by the scale where the non-linear transition takes place, and this occurs at somewhat larger scales for the ΛCDM counterpart (i.e. the higher linear power at higher k in the KM1 scenario is mostly erased by the non-linear dynamics).

As usual, the integration along the line of sight leads to a much smoother curve for the cosmic shear angular power spectrum. Even though the 3D non-linear matter power spectrum is reduced at ɀ = 0 for KM1 as compared with ΛCDM, as shown in the upper right panel, the lensing power is increased, as shown in the lower right panel. This is due to the coupling function A(φ), which is greater than unity at ɀ > 0 and leads to a higher value of the phenomenological function Σ, see Eq. (24), or in other words, to a greater Newton constant at high redshift.

4 Survey specifications and analysis method

We varied as base cosmological parameters: $$\Theta =\left\{{\mathrm{\Omega }}_{\text{m},0},{\mathrm{\Omega }}_{\text{b},0},h,n_{\text{s}},{\sigma }_8\right\},$$(53)

for which we chose the following fiducial values: $$\begin{array}{l}{\text{Ω}}_{\text{m},0}=0.315,{\text{Ω}}_{\text{b},0}=0.049,\hfill \\ h=0.6737,n_{\text{s}}=0.966,\hfill \end{array}$$(54)

where Ω_m,0 = Ω_c,0 + Ω_b,0 + Ω_ν,0. To fix the fiducial value for σ₈, we used the same initial amplitude of primordial perturbations for all models, namely A_s = 2.09681 × 10⁻⁹. The fiducial cosmology includes massive neutrinos with a total mass of Σm_ν = 0.06 eV, but we kept Σm_ν fixed in the following Fisher matrix analysis. In the following, we specify the σ₈ value for each model along with the fiducial values of the model parameters,

1. JBD $$\begin{array}{ll}{\Theta }_{\text{fid},1}=\left\{{\sigma }_8=0.816,{\log }_{10}{\omega }_{\text{BD}}=2.90309\right\}\hfill & \text{\hspace{0.17em}(JBD1);\hspace{0.17em}}\hfill \\ {\Theta }_{\text{fid},2}=\left\{{\sigma }_8=0.812,{\log }_{10}{\omega }_{\text{BD}}=3.39794\right\}\hfill & \text{\hspace{0.17em}(JBD2)}.\hfill \end{array}$$(55)

   These values correspond to ω_BD = 800 and ω_BD = 2500 respectively. In these cases the initial value of the scalar field was fixed to ϕ_ini = 1 / G_N.

2. nDGP $$\begin{array}{l}{\Theta }_{\text{fid},1}=\left\{{\sigma }_8=0.8690,{\log }_{10}{\mathrm{\Omega }}_{\text{rc}}=-0.60206\right\}(\text{nDGP}1);\\ {\Theta }_{\text{fid},2}=\left\{{\sigma }_8=0.8105,{\log }_{10}{\mathrm{\Omega }}_{\text{rc}}=-6\right\}\text{\hspace{0.17em}\hspace{0.17em}(nDGP2)}.\end{array}$$(56)

   These values correspond to Ω_rc = 0.25 and Ω_rc = 10⁻⁶, for nDGP1 and nDGP2, respectively.

3. KM $$\begin{array}{l}{\Theta }_{\text{fid},1}=\left\{{\sigma }_8=0.994,{ϵ}_2=-0.04\right\}(\text{KM}1);\hfill \\ {\Theta }_{\text{fid},2}=\left\{{\sigma }_8=0.813,{ϵ}_2=-0.0001\right\}(\text{KM}2).\hfill \end{array}$$(57)

We followed EC19 to set up the specifics of the photometric probes as we briefly summarise here. The sources were split into ten equi-populated redshift bins whose limits were obtained from the redshift distribution $$n(ɀ)\propto {\left(\frac{ɀ}{{ɀ}_0}\right)}^2\exp \left[-{\left(\frac{ɀ}{{ɀ}_0}\right)}^{3/2}\right],$$(58)

with ${ɀ}_0=0.9/\sqrt{2}$ and the normalisation set by the requirement that the surface density of galaxies is ${\overline{n}}_g=30{\text{arcmin}}^{-2}$. This was then convolved with the sum of two Gaussians to account for the effect of photometric redshift. Galaxy bias was assumed to be constant within each redshfit bin, with fiducial values $b_i=\sqrt{1+{\overline{ɀ}}_i}$, and ${\overline{ɀ}}_i$ the bin centre. Any possible scale dependence of the bias introduced by GR was taken to be negligible.

As in EC19, we considered a Gaussian-only covariance, whose elements are given by $$\mathrm{cov}\left[C_{ij}^{\text{AB}}(\ell ),C_{kl}^{\text{CD}}\left({\ell }^{\prime }\right)\right]=\frac{{\delta }_{\ell {\ell }^{\prime }}^{\text{K}}}{(2\ell +1)f_{\text{sky}}\text{Δ}\ell }$$(59) $$\times \{\left[C_{ik}^{\text{AC}}(\ell )+N_{ik}^{\text{AC}}(\ell )\right]\left[C_{jl}^{\text{BD}}\left({\ell }^{\prime }\right)+N_{jl}^{\text{BD}}\left({\ell }^{\prime }\right)\right]$$(60) $$+\left[C_{il}^{\text{AD}}(\ell )+N_{il}^{\text{AD}}(\ell )\right]\left[C_{jk}^{\text{BC}}\left({\ell }^{\prime }\right)+N_{jk}^{\text{BC}}\left({\ell }^{\prime }\right)\right]\},$$(61)

where the upper-case (lower-case) Latin indexes run over WL, GC_ph (all tomographic bins), ${\delta }_{\ell {\ell }^{\prime }}^{\text{K}}$ is the Kronecker delta symbol coming from the lack of correlation between different multipoles (ℓ,ℓ′), f_sky ≃ 0.36 is the survey sky fraction, and Δℓ denotes the width of the 100 logarithmic equi-spaced multipole bins. For the observables of interest here, we assumed a white noise, so that it is $$N_{ij}^{\text{LL}}(\ell )=\frac{{\delta }_{ij}^{\text{K}}}{{\overline{n}}_i}{\sigma }_{ϵ}^2,$$(62) $$N_{ij}^{\text{GG}}(\ell )=\frac{{\delta }_{ij}^{\text{K}}}{{\overline{n}}_i},$$(63) $$N_{ij}^{\text{GL}}(\ell )=0,$$(64)

where ${\sigma }_{ϵ}^2={0.3}^2$ is the variance of the observed ellipticities.

We continued to follow EC19 to evaluate the Fisher matrix F_αβ (ɀ_i) for the observed galaxy power spectrum. Here, α and β run over the cosmological parameters of the set Θ, and the index i labels the redshift bin, each centred in ɀ_i = {1.0, 1.2, 1.4, 1.65}, whose widths were Δɀ = 0.2 for the first three bins and Δɀ = 0.3 for the last bin. As a difference with respect to EC19, we adopted the direct derivative approach and directly varied the observed galaxy power spectrum with respect to the cosmological parameters. We included two additional redshift-dependent parameters ln bσ₈(ɀ_i) and P_s (ɀ_i) over which we marginalised. The spectroscopic galaxy bias, b(ɀ), and the expected number density of Hα emitters, n(ɀ), are as reported in Table 3 of EC19.

For all probes, we considered an optimistic and a pessimistic scenario defined according to the specifications for WL, GC_ph and GC_sp in Table 1. For GC_sp , we added a third scenario, referred to as quasi-linear, for which we fixed the maximum wavenumber to k_max = 0.15 h Mpc⁻¹. We explored this more conservative case since the underlying matter power spectrum of Eq. (36) that we used in our observed galaxy power spectrum recipe is a linear one. Non-linear corrections begin to play a role for wave numbers larger than k = 0.1 h Mpc⁻¹ for the redshifts we considered (see Taruya et al. 2010). Cutting at a lower k_max should avoid any bias induced by neglecting non-linear corrections below this scale. However, this severe cut also removes part of the information encoded in GC_sp so that we wished to quantify how this affects the constraints. In all scenarios, we fixed the σ_v = σ_p nuisance parameter for GC_sp to the values directly calculated from Eq. (38) for the fiducial cosmological parameters.

We considered both GC_sp alone and the combination of all photometric probes: GC_ph, WL, and their XC. Moreover, we also considered the full combination of all Euclid main probes: GC_sp , GC_ph , WL, and the XC between GC_ph and WL. It is important to mention that we accounted for all the mixed terms in the covariance matrix for the photometric probes such as, for instance, $\mathrm{Cov}\left[C_{ij}^{\text{LL}}(\ell ),C_{kl}^{\text{GL}}(\ell )\right]$. However, we followed EC19 in neglecting any correlation between GC_sp and the photometric probes. In the pessimistic scenario, we further imposed a cut ɀ < 0.9 to GC_ph and XC to remove any overlap between these probes and GC_sp . This cut was applied only when we added GC_sp to the photometric probes, while was not applied when we used the photometric probes alone in the pessimistic case. No cut was ever applied in the optimistic case.

5 Results

In this section, we discuss the results of the Fisher matrix analysis for the three models under investigation, JBD, nDGP, and KM, and for each of them we considered the two fiducial cases presented in Sect. 4. Additionally, for each case we considered a quasi-linear, pessimistic, and optimistic scenario as discussed in Sect. 4. In the following we present the 68.3 and 95.4% joint marginal error contours on the cosmological and model parameters for the optimistic setting, and for completeness, we include the contour plots for the pessimistic scenario in Appendix A. We stress that we are interested in detecting a non-zero value for the model parameters, and therefore we present the relative errors on the corresponding fiducial values. In the specific cases of JBD and nDGP, we performed the Fisher analysis on the logarithm of the additional parameter, namely log₁₀ X with X = ω_BD and X = Ω_rc for JBD and nDGP, respectively. This choice was dictated by the fact that we cannot perform a Fisher analysis on a parameter that varies by some orders of magnitude. Additionally, we need to use as parameter a quantity that is of order one in the analysis in order to avoid large differences between the highest and lowest eigenvalues of the Fisher matrix. For the specific case of nDGP, the use of the logarithm of the model parameter allowed us to consider only its positive viable range, which would not be possible if we performed a Fisher analysis directly on Ω_rc. We note indeed that a prior such as Ω_rc > 0 cannot be imposed a priori in the Fisher analysis. In addition to discussing the results on log₁₀ X we also report on the uncertainties on the parameter itself, X. We note that we cannot use a Jacobian transformation to convert between the Fisher matrices because the transformation between the two parametrisations is non-linear and the assumption of Gaussianity is only valid for the logarithmic parametrisation. In order to obtain the constraints on the parameter X, we used $$X^{(\pm )}=X_{\text{fid}}\times {10}^{\pm {\sigma }_{{\log }_{10}X}}.$$(65)

As a consequence, the uncertainties on X are asymmetrical. We now proceed to discuss each model in detail.

In order to have a global view on the results for the three MG models analysed, we summarise in Fig. 2 the 68.3% marginal errors on the model parameters, log₁₀ ω_BD for JBD, log₁₀ Ω_rc for nDGP, and є_2,0 for KM in all the scenarios and combinations of probes we used. We note that we did not consider WL probe alone because we verified that WL provides constraints of about the same order as the GC_sp ones. When we combined this with GC_ph and XC, the constraints were significantly boosted. Additionally, we note that the constraining power of the spectroscopic sample is very weak because we did not go deep enough (k ≳ 1 h Mpc⁻¹) in the non-linear regime and also because when we used GC probes, we introduced degeneracies among the MG parameters, amplitudes, and bias parameters. These features characterise all the models.

5.1 Jordan-Brans-Dicke gravity

There are remarkably strong non-cosmological constraints on JBD gravity, with current bounds of about ω_BD > 10⁵. As pointed out above, however, it is useful to consider JBD gravity as a long-wavelength limit of more general scalar-tensor theories with a wide range of small-scale limits that may be endowed with gravitational screening. It then makes sense to focus on what we already know from cosmology, where the constraints are about ω_BD > 10³, independent of what we know from non- cosmological systems. Two further practical reasons led us to choose ω_BD ~ 10³ as a fiducial value. First of all, this value is at the limits of the range for which we calibrated the nonlinear prescription we used. Second, the likelihood is very flat for high values of ω_BD which can lead to numerical errors in the finite differences of the Fisher matrix evaluation for some of the observables. We thus considered two cases for JBD. The first fiducial case with ω_BD = 800 (JBD1), and the second case with ω_BD = 2500 (JBD2). Whereas the first case is nearly compatible with the most recent constraints from publicly available CMB and LSS data, the second case is well within the current bounds (Ballardini et al. 2020; Joudaki et al. 2022; Ballardini et al. 2022). Table 2 summarises the relative (with respect to the fiducial values) uncertainties for the cosmological and model parameters in the quasi-linear, pessimistic and optimistic Euclid cases. Fig. 3 shows the 68.3% and 95.4% joint marginal error contours on the cosmological model parameters for the optimistic settings for both JBD1 (left panel) and JBD2 (right panel). As in Casas et al. (2023) for f(R), we performed the Fisher matrix analysis on the parameter log₁₀ ω_BD, and not directly on ω_BD, because for large numbers and differences in the order of magnitude, the Fisher matrix derivatives might become unstable.

In the optimistic setting, we find that Euclid will be able to constrain log₁₀ ω_BD as follows:

• JBD1:

  • at 27% for the GC_sp alone;

  • at 3.6% when considering the combination WL+XC+GC_ph;

  • at 3.2% for GC_sp+WL+XC+GC_ph.

• JBD2:

  • at 72.2% for the GC_sp alone;

  • at 5.1% when considering the combination WL+XC+GC_ph;

  • at 4.6% for GC_Sp + WL+XC+GC_ph.

For JBD1 the parameter log₁₀ ω_BD can be measured at high statistical significance even in the pessimistic case, but this does not hold for JBD2.

We now discuss the forecast errors on the parameter ω_BD propagating the uncertainties according to Eq. (65). For JBD1 we obtain ${\omega }_{\text{BD}}={800}_{-670}^{+4100}$ for GC_sp, ${\omega }_{\text{BD}}={800}_{-170}^{+210}$ for WL+XC+GC_ph and ${\omega }_{\text{BD}}={800}_{-160}^{+200}$ for GC_sp+WL+XC+GC_ph; these results correspond to relative errors of 513, 27 and 24%, respectively. For JBD2, always in the case of optimistic settings, we obtain ${\omega }_{\text{BD}}={2500}_{-820}^{+1200}$ for WL+XC+GC_ph and ${\omega }_{\text{BD}}={2500}_{-750}^{+1070}$ for GC_Sp+WL+XC+GC_ph. This corresponds to relative errors of 49 and 43%, respectively (JBD2 is unconstrained by GC alone).

These results for the optimistic settings demonstrate that Euclid alone could potentially detect at statistically significant level values of ω_BD which are compatible with current publicly available CMB and LSS data (Ballardini et al. 2020; Joudaki et al. 2022; Ballardini et al. 2022), but not those that are compatible with Solar System constraints. Although the Euclid specifications we used here are different from those used in Ballardini et al. (2019), the results we find are consistent with theirs.

Fig. 2 Marginalised 1σ errors on the model parameters relative to their corresponding fiducial values, for JBD (case 1 in blue and case 2 in cyan), nDGP (case 1 in green and case 2 in light green), and KM (case 1 in red and case 2 in light red). We show the marginalised 1σ errors for GCsp in the optimistic (dotted area), pessimistic (dashed area) and quasi-linear (filled area) scenarios (top panel), for WL+XC+GCph in the optimistic and pessimistic cases (central panel) and for the GCSp+WL+XC+GCph in the optimistic and pessimistic cases.

5.2 Dvali-Gabadadze-Porrati gravity

The two cases considered for the nDGP model are log₁₀ Ω_rc = −0.60206 (nDGPl) and log₁₀ Ω_rc = −6 (nDGP2). We summarise the 68.3% marginalised errors (relative to the fiducial values) in Table 3 on the cosmological and model parameters in the optimistic, pessimistic and quasi-linear case. We show in Fig. 4 the 68.3 and 95.4% joint marginal error contours on cosmological and model parameters for the optimistic scenario for both nDGPl (left panel) and nDGP2 (right panel). In Fig. A.2 we show the same but for the pessimistic setting.

In the optimistic setting, we find that Euclid will be able to constrain log₁₀ Ω_rc as follows:

• nDGPl:

  • at 93.4% for GC_sp alone;

  • at 20% when considering WL+XC+GC_ph;

  • at 15% for GC_sp+WL+XC+GC_ph.

• nDGP2:

  • GC_sp has no power in constraining the parameter.

  • at 81% for both WL+XC+GC_ph and GC_sp+WL+XC+GC_ph.

Regardless of the specific case considered (nDGP1 or nDGP2), the power in constraining Ω_m,0 and Ω_b,0 is the same, while n_s and σ₈ are better constrained in the nDGP2 case and h in the nDGP1 case. For example, in the optimistic setting, n_s for WL+XC+GC_ph is constrained at 0.52% for nDGP1 and 0.37% for nDGP2 and they further improve when considering the full combination of probes, when they are 0.27 and 0.19%, respectively. In Fig. 4 we note in the log₁₀ Ω_rc−n_s panel that the two parameters are anti-correlated in the nDGP1 case. This makes the forecast errors on n_s larger than the nDGP2 case where the anti-correlation disappears. A similar discussion can be made for the σ₈ parameter. In the optimistic setting, for nDGP1, we have 0.37% for WL+XC+GC_ph and 0.21% for the full combination of probes, while for the nDGP2 model we obtain 0.25 and 0.14%, respectively. Instead h is better constrained in nDGP1 being 1.25 and 1.35% for nDGP2 using GC_sp.

As previously done for the JBD model, we translated the forecast errors on the parameter log₁₀ Ω_rc into the parameter Ω_rc [see Eq. (65)]. For the optimistic scenario, we obtain for nDGP1 ${\mathrm{\Omega }}_{\text{rc}}={0.25}_{-0.18}^{+0.66}$ for GC_sp, ${\mathrm{\Omega }}_{\text{rc}}={0.25}_{-0.06}^{+0.08}$ for WL+XC+GC_ph and ${\mathrm{\Omega }}_{\text{rc}}={0.25}_{-0.05}^{+0.06}$ for GC_sp+WL+XC+GC_ph. These correspond to relative errors that can reach 264, 32, and 24%, respectively.

For nDGP2, we obtain an upper bound that is Ω_rc < 0.07 in the optimistic setting, while for the pessimistic case, Ω_rc is unconstrained regardless of the probe we consider because the fiducial parameter is very close to ACDM and we do not have enough non-linear information. To examine this, we considered the linear power spectrum: We find that the derivative with respect to log₁₀ Ω_rc approaches 0 in the ΛCDM limit, and ${{P^{\prime }}_{\text{L}}|}_{{10}^{-6}}\approx 0.0002$ compared to ${{P^{\prime }}_{\text{L}}|}_{0.25}\approx 0.08$, where a prime here denotes a derivative with respect to log₁₀ Ω_rc. This leads us to conclude that the pessimistic case is still too dominated by linear structure information to constrain Ω_rc when we are close to the ACDM limit. It should also be noted that Bose et al. (2020) reported meaningful constraints using cosmic shear with the same scale cuts we employed in our pessimistic case, but they used a full Markov chain Monte Carlo (MCMC) approach to a ΛCDM fiducial data vector, sampling in Ω_rc, as well as employing a different covariance and binning scheme. Additionally, it is known that the results obtained with a Fisher approach can be different from those computed with the MCMC method even though the sampled parameter is the same (Perotto et al. 2006; Wolz et al. 2012). Therefore a direct comparison of the results is not straightforward. Lastly, our loss in sensitivity in the ΛCDM limit justifies the tight upper limit in nDGP1 when compared to nDGP2.

Moreover, a cut in the maximum multipole and scale at lower values, as for the pessimistic (or quasi-linear) scenario, leads to larger uncertainties on the model parameter (log₁₀ Ω_rc) compared to the optimistic setting, as shown in Table 3 and in Fig. 2, where the impact of the different scales is highlighted. Conversely, the relative errors on the cosmological parameters deteriorate only slightly. This means that most of the power of Euclid in constraining this model comes from the non-linear scales.

Finally, in this analysis we forecast the logarithm of Ω_rc following the motivation discussed in Sect. 5, but we tested that a similar constraining power can be directly obtained by sampling Ω_rc. For the 3 × 2pt in the optimistic setting we obtain for nDGP1 a 24% constraint with respect to the 32% obtained on Ω_rc indirectly from the log sampling.

Fig. 3 68.3 and 95.4% joint marginal error contours on the cosmological parameters for JBD1 (left panel) and JBD2 (right panel) in the optimistic case. In red, we plot GCsp, in blue, we plot WL+XC+GCph and in green we plot GCsp+WL+XC+GCph.

Fig. 4 68.3 and 95.4% joint marginal error contours on the cosmological parameters for nDGP1 (left panel) and nDGP2 (right panel) in the optimistic case. In red, we plot GCsp, in blue, we plot WL+XC+GCph and in green we plot GCsp+WL+XC+GCph.

Fig. 5 68.3 and 95.4% joint marginal error contours on the cosmological parameters for KM1 (left panel) and KM2 (right panel) in the optimistic case. In red, we plot GCsp, in blue, we plot WL+XC+GCph and in green, we plot GCsp+WL+XC+GCph.

5.3 k-mouflage

We present our results for the KM model with the two fiducial choices of the ϵ_2,0 parameter (KM1 ϵ_2,0 = −0.04 and KM2 ϵ_2,0 = −0.0001). The KM1 model represents an extreme case that is allowed within the 95.4% confidence interval by current cosmological data, while KM2 is practically indistinguishable from ΛCDM. We show in Fig. 5 the forecast 68.3 and 95.4% joint marginal error contours on cosmological and model parameters for the optimistic scenario. Fig. A.3 shows forecast results for the same models in the pessimistic scenario. The 68.3% errors for both the optimistic and pessimistic configurations are reported in Table 4.

In the optimistic setting, we find that Euclid will be able to constrain the modified-gravity parameter ϵ_2,0 as follows:

• KM1:

  • at 3.4% for the GC_sp alone;

  • at 0.15% when considering the combination WL+XC+GC_ph;

  • at 0.14% forGC_sp+WL+XC+GC_ph.

• KM2:

  • none of the combinations of probes we consider are able to constrain KM2;

  • upper bound at 10⁻³ for the GC_sp alone;

  • upper bound at 7 × 10⁻⁴ when considering the combination WL+XC+GC_ph;

  • upper bound at 4 × 10⁻⁴ for GC_sp +WL+XC+GC_ph.

For the KM1 fiducial model, the optimistic configuration will allow us to constrain ϵ_2,0 to ~ 3.4% accuracy using GC_sp alone. In the worst-case scenario of cutting the observed galaxy power spectrum at k_max = 0.15 h Mpc⁻¹, the error on ϵ_2,0 increases to ~ 5.5%. The constraining power greatly improves for the combination WL+XC+GC_ph, and the percentage error on ϵ_2,0 decreases to 0.15% in the optimistic case, without a substantial worsening in the pessimistic scenario. Of the cosmological parameters we varied for the KM1 model, ϵ_2,0 is most tightly constrained by the combination WL+XC+GC_ph. This shows that this combination of data is well suited to capturing the modification induced by k-mouflage to the Poisson and lensing equations, probing µ and ∑ independently. Moreover, the Euclid survey will probe the redshift range ɀ ≳ 1, where the largest effect determined by the running of Newton’s constant is expected to manifest (see, e.g., Benevento et al. 2019). The fact that with the photometric probes we obtain the tightest relative constraints for the KM1 model, compared to the other two models studied in this work, is consistent with the observation of the lower panels of Figure 1, which show that the lensing angular power spectrum is more distant to its λCDM counterpart than in the other two models. We note that the addition of GC_sp to WL+XC+GC_ph does not improve the constraints on ϵ₂₀ in either the pessimistic and optimistic scenarios.

The percentage error on ϵ₂₀ in the ΛCDM-proximate KM2 fiducial model decreases from 1120% using GC_sp alone to a value of 397% with the full combination of GC_sp+WL+XC+GC_ph in the best-case scenario. In the optimistic case, the Euclid survey will therefore be able to detect a deviation from λCDM only when |ϵ₂₀| ≳ 4 × 10⁻⁴, which improves over the constraining power of present CMB and LSS data by somewhat more than one order of magnitude.

Fig. 5 shows that ϵ_2,0 is anti-correlated with the h parameter. This correlation can be exploited to reduce the Hubble tension, as noted in Benevento et al. (2019). This effect determines a lower forecast error for the h parameter in the KM1 fiducial model, for which ϵ_2,0 is tightly constrained by Euclid probes, especially when the full data combination is considered. A similar argument applies to the Ω_m,0 parameter, which is also better constrained for KM1. The uncertainties for the other parameters do not vary much with the fiducial choice of ϵ_2,0.

6 Conclusions

We have explored the constraining power of the future Euclid mission for linearly scale-independent extensions of the concordance cosmological model, that is, models that induce modifications to the linear growth of perturbations that are solely time-dependent while featuring different testable types of screening mechanisms at smaller non-linear scales. We considered three specific models, namely JBD, a scalar-tensor theory with a standard kinetic term and a flat potential (Sect. 2.1); the nDGP gravity, a braneworld model in which our Universe is a four-dimensional brane embedded in a five-dimensional Minkowski space-time (Sect. 2.2); and KM gravity, an extension of k-essence scenarios with a universal coupling of the scalar field to matter (Sect. 2.3).

We derived forecasts from Euclid spectroscopic and photometric primary probes on the cosmological parameters and the additional parameters of the models, log₁₀ ω_BD for JBD, log₁₀Ω_rc for nDGP, and ϵ₂ for KM. In order to do this, we applied the Fisher matrix method to weak lensing (WL), photometric galaxy clustering (GC_ph), spectroscopic galaxy clustering (GC_sp), and the cross-correlation (XC) between GC_ph and WL. For each MG model, we considered two fiducial values for the corresponding model parameter, following the rationale of a case that is representative of the λCDM limit and another that differs more significantly while still being (nearly) compatible with current bounds. We modelled the non-linear matter power spectrum using different prescriptions for each MG model: The HMCODE (Mead et al. 2015, 2016) calibrated on a suite of N-body simulations from modified versions of COLA (Tassev et al. 2013; Winther et al. 2017) and RAMSES (Teyssier 2002) for JBD; the halo model reaction (Cataneo et al. 2019) for nDGP; and an analytical approach that combined one-loop perturbation theory with a halo model for KM following (Brax & Valageas 2014; Brax et al. 2015).

When setting the Euclid survey specifications, we defined three scenarios that were characterized by different cuts in the maximum multipole and wavenumber to assess the constraining power of non-linear scales: the quasi-linear scenario with k_max = 0.15 h Mpc⁻¹ for GC_sp; the pessimistic scenario with k_max = 0.25 h Mpc⁻¹ for GC_sp, ℓ_max = 1500 for WL, and ℓ_max = 750 for GC_ph; and the optimistic scenario with k_max = 0.3 h Mpc⁻¹ for GC_sp, ℓ_max = 5000 for WL, and ℓ_max = 3000 for GC_ph. In Sect. 5 we discussed and reported in great detail the results corresponding to these different scenarios. For each case, we considered three different combinations of data: GC_sp alone, WL+XC+GC_ph and the full GC_sp +WL+XC+GC_ph. For each of these combinations we presented the relative errors on the fiducial values of the three models as they can provide indications about possible future detections of non-zero modified-gravity parameters.

With the full set of probes, we found that in the optimistic scenario, Euclid alone will be able to constrain the JBD parameter ${\omega }_{\text{BD}}={800}_{-160}^{+200}$ in the fiducial case JBD1 for the full combination of probes, and in the JBD2 case ${\omega }_{\text{BD}}={2500}_{-750}^{+1070}$. This indicates that Euclid alone will be capable of detecting at a statistically significant level only values of ω_BD that are compatible with the current bounds from CMB and LSS.

For the nDGP model, we obtain for the optimistic scenario that ${\mathrm{\Omega }}_{\text{rc}}={0.25}_{-0.05}^{+0.06}$ of with the full combination in the nDGP1 fiducial case, while for the nDGP2 fiducial (ω_rc = 10⁻⁶) we find that ω_rc < 0.07. Furthermore, the nDGP2 case is closer to the λCDM model than the KM2 model. The constraints of nDGP2 are therefore worse than those of KM2.

For the KM scenario, we obtain constraints for the full combination of probes in the fiducial KM1 case of ϵ_2,0 = −0.04 ± 5.6 × 10⁻⁵. Out of the full set of cosmological plus model parameters, ϵ_2,0 is the most tightly constrained parameter for KM1 by the full combination because of the independent constraints on clustering and lensing, which probe the running of Newton’s constant in the range ɀ ≳ 1 where the strongest effect is expected. The KM2 fiducial model (with ϵ₂₀ = −0.0001) is unconstrained for all the different combinations of probes.

Although it has been shown that for values close to the λCDM limit, the nDGP2 and KM2 models remain unconstrained in the pessimistic setting, it is remarkable that Euclid alone can significantly reduce the allowed space for these models. For the full probe combination (GC_sp+WL+XC+GC_ph) the nDGP2 upper bound at ω_rc < 0.07 for the optimistic scenario significantly improves the current constraints derived with the same approach and similar type of data. This is consistent with the results found in Bose et al. (2020).

The forecasts for KM2 give |ϵ_2,0| < 4 × 10⁻⁴ in the optimistic case and |ϵ₂₀| < 2 × 10⁻³ in the pessimistic case. This improves the current constraints by more than one order of magnitude, even in the pessimistic case, by combining different cosmological datasets (Benevento et al. 2019). In addition, Euclid alone will provide constraints on ω_BD that will be tighter than those obtained by a combination of current CMB and low-redshift data Joudaki et al. (2022).

Finally, our analysis clearly showed that most of the constraining power of Euclid comes from the non-linear scales. We conclude that Euclid will be able to provide outstanding constraints on extensions beyond the concordance model given a good modelling of our theoretical observables at these scales, such as we used in this analysis, seen from the substantial differences in constraining power between the optimistic and pessimistic case. The current dedicated modelling of non-linearities for spectroscopic galaxy clustering, developed and tested for λCDM and applied to BOSS data by D’Amico et al. (2020); Ivanov et al. (2020); Chen et al. (2022), for instance, can be straightforwardly extended to the class of scale-independent models considered in this paper, as was already shown, in the nDGP model in Piga et al. (2023), for example. Additionally we would like to stress that the combination and cross-correlation of future Euclid data with CMB measurements will be of crucial importance for the cosmological constraints of extended models. As already shown within the Euclid Collaboration (Euclid Collaboration 2022), the improvements in the constraints when Euclid is cross-correlated with CMB data can be of the order of two to three and in some cases even larger. A large effort is currently ongoing to extend this analysis to other extended models, including those in this work.

Acknowledgements

The Euclid Consortium acknowledges the European Space Agency and a number of agencies and institutes that have supported the development of Euclid, in particular the Academy of Finland, the Agenzia Spaziale Italiana, the Belgian Science Policy, the Canadian Euclid Consortium, the French Centre National d’Etudes Spatiales, the Deutsches Zentrum für Luft- und Raumfahrt, the Danish Space Research Institute, the Fundação para a Ciência e a Tecnologia, the Ministerio de Ciencia e Innovación, the National Aeronautics and Space Administration, the National Astronomical Observatory of Japan, the Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the Romanian Space Agency, the State Secretariat for Education, Research and Innovation (SERI) at the Swiss Space Office (SSO), and the United Kingdom Space Agency. A complete and detailed list is available on the Euclid web site (http://www.euclid-ec.org). B.B. and L.L. were supported by a Swiss National Science Foundation Professorship grant (Nos. 170547 & 202671). E.B. has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 754496. B.B was supported by a UK Research and Innovation Stephen Hawking Fellowship (EP/W005654/2). P.G.F. acknowledges funding from STFC, the Beecroft Trust and the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme (grant agreement No 693024). N.F. is supported by the Italian Ministry of University and Research (MUR) through Rita Levi Montalcini project “Tests of gravity at cosmological scales” with reference PGR19ILFGP. N.F., F.P. and A.R.F. also acknowledge the FCT project with ref. number PTDC/FIS-AST/0054/2021. K.K. is supported by the UK STFC grant ST/S000550/1 and ST/W001225/1. For the purpose of open access, the author(s) has applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising. F.P. acknowledges partial support from the INFN grant InDark and the Departments of Excellence grant L.232/2016 of the Italian Ministry of University and Research (MUR). A.R.F. acknowledges support from DL 57/2016 from the ‘Departamento de Física, Faculdade de Ciências, Universidade de Lisboa’. I.T. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 863929; project title “Testing the law of gravity with novel large-scale structure observables” and acknowledges support from the Spanish Ministry of Science, Innovation and Universities through grant ESP2017-89838, and the H2020 programme of the European Commission through grant 776247. Z.S. acknowledges funding from DFG project 456622116. G.B. acknowledges support from ASI/INFN grant n. 202one-43-HH.0. M.P. acknowledges support by the MIUR ‘Progetti di Ricerca di Rilevante Interesse Nazionale’ (PRIN) Bando 2022 - grant 20228RMX4A.

Appendix A Further results

In this appendix we show additional results to complement what we discussed in Sect. 5. In Figs. A.1, A.2, and A.3 we provide the 68.3% and 95.4% joint marginal error contours in the pessimistic case on the cosmological parameters for JBD, nDGP, and KM respectively.

Fig. A.1 68.3% and 95.4% joint marginal error contours on the cosmological parameters for JBD1 (left panel) and JBD2 (right panel) in the pessimistic case. In red, we plot GCsp, in blue, we plot WL+XC+GCph and in green, we plot GCsp+WL+XC+GCph.

Fig. A.2 68.3% and 95.4% joint marginal error contours on the cosmological parameters for nDGP1 (left panel) and nDGP2 (right panel) in the pessimistic case. In red, we plot GCsp, in blue, we plot WL+XC+GCph and in green, we plot GCsp+WL+XC+GCph.

Fig. A.3 68.3% and 95.4% joint marginal error contours on the cosmological parameters for KM1 (left panel) and KM2 (right panel) in the pessimistic case. In red, we plot GCsp, in blue, we plot WL+XC+GCph and in green, we plot GCsp+WL+XC+GCph.

References