Issue 
A&A
Volume 690, October 2024



Article Number  A133  
Number of page(s)  22  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/202347526  
Published online  04 October 2024 
Euclid: Constraining linearly scaleindependent modifications of gravity with the spectroscopic and photometric primary probes^{★}
^{1}
Department of Physics “E. Pancini”, University Federico II,
Via Cinthia 6,
80126
Napoli,
Italy
^{2}
Dipartimento di Fisica, Università degli Studi di Torino,
Via P. Giuria 1,
10125
Torino,
Italy
^{3}
INFNSezione di Torino,
Via P. Giuria 1,
10125
Torino,
Italy
^{4}
INAFOsservatorio Astrofisico di Torino,
Via Osservatorio 20,
10025
Pino Torinese (TO),
Italy
^{5}
INAFOsservatorio Astronomico di Roma,
Via Frascati 33,
00078
Monteporzio Catone,
Italy
^{6}
INFNSezione di Roma, Piazzale Aldo Moro,
2  c/o Dipartimento di Fisica, Edificio G. Marconi,
00185
Roma,
Italy
^{7}
Institute for Theoretical Particle Physics and Cosmology (TTK), RWTH Aachen University,
52056
Aachen,
Germany
^{8}
Université de Genève, Département de Physique Théorique and Centre for Astroparticle Physics,
24 quai ErnestAnsermet,
1211
Genève 4,
Switzerland
^{9}
Institute of Space Sciences (ICE, CSIC), Campus UAB,
Carrer de Can Magrans, s/n,
08193
Barcelona,
Spain
^{10}
Institut d’Estudis Espacials de Catalunya (IEEC),
Carrer Gran Capitá 24,
08034
Barcelona,
Spain
^{11}
Institut de Recherche en Astrophysique et Planétologie (IRAP), Université de Toulouse, CNRS, UPS, CNES,
14 Av. Edouard Belin,
31400
Toulouse,
France
^{12}
Dipartimento di Fisica e Scienze della Terra, Universitá degli Studi di Ferrara,
Via Giuseppe Saragat 1,
44122
Ferrara,
Italy
^{13}
Dipartimento di Fisica e Astronomia, Universitá di Bologna,
Via Gobetti 93/2,
40129
Bologna,
Italy
^{14}
INAFOsservatorio di Astrofisica e Scienza dello Spazio di Bologna,
Via Piero Gobetti 93/3,
40129
Bologna,
Italy
^{15}
INFNSezione di Bologna,
Viale Berti Pichat 6/2,
40127
Bologna,
Italy
^{16}
INFN, Sezione di Trieste,
Via Valerio 2,
34127
Trieste TS,
Italy
^{17}
IFPU, Institute for Fundamental Physics of the Universe,
via Beirut 2,
34151
Trieste,
Italy
^{18}
SISSA, International School for Advanced Studies,
Via Bonomea 265,
34136
Trieste TS,
Italy
^{19}
Johns Hopkins University
3400 North Charles Street
Baltimore,
MD
21218,
USA
^{20}
Institute for Astronomy, University of Edinburgh,
Royal Observatory, Blackford Hill,
Edinburgh
EH9 3HJ,
UK
^{21}
Institute for Computational Science, University of Zurich,
Winterthurerstrasse 190,
8057
Zurich,
Switzerland
^{22}
Institut de Physique Théorique, CEA, CNRS, Université ParisSaclay
91191
GifsurYvette Cedex,
France
^{23}
Dipartimento di Fisica e Astronomia “G.Galilei”, Università di Padova,
Via Marzolo 8,
35131
Padova,
Italy
^{24}
INFNPadova,
Via Marzolo 8,
35131
Padova,
Italy
^{25}
INAFOsservatorio Astronomico di Padova,
Via dell’Osservatorio 5,
35122
Padova,
Italy
^{26}
CERN, Theoretical Physics Department,
Geneva,
Switzerland
^{27}
Department of Physics, Oxford University,
Keble Road,
Oxford
OX1 3RH,
UK
^{28}
INFNBologna,
Via Irnerio 46,
40126
Bologna,
Italy
^{29}
Institute of Cosmology and Gravitation, University of Portsmouth,
Portsmouth
PO1 3FX,
UK
^{30}
Dipartimento di Scienze Matematiche, Fisiche e Informatiche, Università di Parma,
Viale delle Scienze 7/A
43124
Parma,
Italy
^{31}
INFN Gruppo Collegato di Parma,
Viale delle Scienze 7/A
43124
Parma,
Italy
^{32}
Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências, Universidade de Lisboa,
Campo Grande,
1749016
Lisboa,
Portugal
^{33}
Institut d’Astrophysique de Paris, UMR 7095, CNRS, and Sorbonne Université,
98 bis boulevard Arago,
75014
Paris,
France
^{34}
Institut für Theoretische Physik, University of Heidelberg,
Philosophenweg 16,
69120
Heidelberg,
Germany
^{35}
Université St Joseph; Faculty of Sciences,
Beirut,
Lebanon
^{36}
Institute Lorentz, Leiden University,
PO Box 9506,
Leiden
2300 RA,
The Netherlands
^{37}
Institute of Theoretical Astrophysics, University of Oslo,
PO Box 1029
Blindern,
0315
Oslo,
Norway
^{38}
Université ParisSaclay, CNRS, Institut d’astrophysique spatiale,
91405,
Orsay,
France
^{39}
Mullard Space Science Laboratory, University College London,
Holmbury St Mary, Dorking,
Surrey
RH5 6NT,
UK
^{40}
Dipartimento di Fisica, Università degli studi di Genova, and INFNSezione di Genova,
via Dodecaneso 33,
16146
Genova,
Italy
^{41}
INFNSezione di Roma Tre,
Via della Vasca Navale 84,
00146
Roma,
Italy
^{42}
INAFOsservatorio Astronomico di Capodimonte,
Via Moiariello 16,
80131
Napoli,
Italy
^{43}
Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto,
CAUP, Rua das Estrelas,
PT4150762
Porto,
Portugal
^{44}
INAFIASF Milano,
Via Alfonso Corti 12,
20133
Milano,
Italy
^{45}
Institut de Física d’Altes Energies (IFAE), The Barcelona Institute of Science and Technology,
Campus UAB,
08193
Bellaterra (Barcelona),
Spain
^{46}
Port d’Informació Científica,
Campus UAB, C. Albareda s/n,
08193
Bellaterra (Barcelona),
Spain
^{47}
INFN section of Naples,
Via Cinthia 6,
80126,
Napoli,
Italy
^{48}
Dipartimento di Fisica e Astronomia “Augusto Righi”  Alma Mater Studiorum Università di Bologna,
Viale Berti Pichat 6/2,
40127
Bologna,
Italy
^{49}
INAFOsservatorio Astrofisico di Arcetri,
Largo E. Fermi 5,
50125
Firenze,
Italy
^{50}
Centre National d’Etudes Spatiales,
Toulouse,
France
^{51}
Institut national de physique nucléaire et de physique des particules,
3 rue MichelAnge,
75794
Paris Cédex 16,
France
^{52}
European Space Agency/ESRIN,
Largo Galileo Galilei 1,
00044
Frascati, Roma,
Italy
^{53}
ESAC/ESA, Camino Bajo del Castillo,
s/n., Urb. Villafranca del Castillo,
28692
Villanueva de la Cañada, Madrid,
Spain
^{54}
Univ Lyon, Univ Claude Bernard Lyon 1, CNRS/IN2P3,
IP2I Lyon, UMR 5822,
69622
Villeurbanne,
France
^{55}
Institute of Physics, Laboratory of Astrophysics, Ecole Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny,
1290
Versoix,
Switzerland
^{56}
Departamento de Física, Faculdade de Ciências, Universidade de Lisboa,
Edifício C8, Campo Grande,
PT1749016
Lisboa,
Portugal
^{57}
Department of Astronomy, University of Geneva,
ch. d’Ecogia 16,
1290
Versoix,
Switzerland
^{58}
Université ParisSaclay, Université Paris Cité, CEA, CNRS, Astrophysique, Instrumentation et Modélisation ParisSaclay,
91191
GifsurYvette,
France
^{59}
INAFOsservatorio Astronomico di Trieste,
Via G. B. Tiepolo 11,
34143
Trieste,
Italy
^{60}
Istituto Nazionale di Fisica Nucleare, Sezione di Bologna,
Via Irnerio 46,
40126
Bologna,
Italy
^{61}
Max Planck Institute for Extraterrestrial Physics,
Giessenbachstr. 1,
85748
Garching,
Germany
^{62}
UniversitätsSternwarte München, Fakultät für Physik, LudwigMaximiliansUniversität München,
Scheinerstrasse 1,
81679
München,
Germany
^{63}
Dipartimento di Fisica “Aldo Pontremoli”, Università degli Studi di Milano,
Via Celoria 16,
20133
Milano,
Italy
^{64}
INFNSezione di Milano,
Via Celoria 16,
20133
Milano,
Italy
^{65}
INAFOsservatorio Astronomico di Brera,
Via Brera 28,
20122
Milano,
Italy
^{66}
Jet Propulsion Laboratory, California Institute of Technology,
4800 Oak Grove Drive,
Pasadena,
CA,
91109,
USA
^{67}
von Hoerner & Sulger GmbH,
SchloßPlatz 8,
68723
Schwetzingen,
Germany
^{68}
Technical University of Denmark,
Elektrovej 327,
2800
Kgs. Lyngby,
Denmark
^{69}
MaxPlanckInstitut für Astronomie,
Königstuhl 17,
69117
Heidelberg,
Germany
^{70}
AixMarseille Université, CNRS/IN2P3, CPPM,
Marseille,
France
^{71}
Department of Physics and Helsinki Institute of Physics,
Gustaf Hällströmin katu 2, 00014 University of Helsinki,
Finland
^{72}
NOVA optical infrared instrumentation group at ASTRON,
Oude Hoogeveensedijk 4,
7991PD
Dwingeloo,
The Netherlands
^{73}
ArgelanderInstitut für Astronomie, Universität Bonn,
Auf dem Hügel 71,
53121
Bonn,
Germany
^{74}
Dipartimento di Fisica e Astronomia “Augusto Righi”  Alma Mater Studiorum Università di Bologna,
via Piero Gobetti 93/2,
40129
Bologna,
Italy
^{75}
Department of Physics, Institute for Computational Cosmology, Durham University,
South Road,
DH1 3LE,
UK
^{76}
European Space Agency/ESTEC,
Keplerlaan 1,
2201 AZ
Noordwijk,
The Netherlands
^{77}
Department of Physics and Astronomy, University of Aarhus,
Ny Munkegade 120,
8000
Aarhus C,
Denmark
^{78}
Centre for Astrophysics, University of Waterloo,
Waterloo,
Ontario
N2L 3G1,
Canada
^{79}
Department of Physics and Astronomy, University of Waterloo,
Waterloo,
Ontario
N2L 3G1,
Canada
^{80}
Perimeter Institute for Theoretical Physics,
Waterloo,
Ontario
N2L 2Y5,
Canada
^{81}
Space Science Data Center, Italian Space Agency,
via del Politecnico snc,
00133
Roma,
Italy
^{82}
Institute of Space Science,
Bucharest,
077125,
Romania
^{83}
Instituto de Astrofísica de Canarias,
Calle Vía Láctea s/n,
38204,
San Cristóbal de La Laguna,
Tenerife,
Spain
^{84}
Departamento de Astrofísica, Universidad de La Laguna,
38206,
La Laguna,
Tenerife,
Spain
^{85}
Departamento de Física, FCFM, Universidad de Chile,
Blanco Encalada 2008,
Santiago,
Chile
^{86}
AixMarseille Université, CNRS, CNES, LAM,
Marseille,
France
^{87}
Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT),
Avenida Complutense 40,
28040
Madrid,
Spain
^{88}
Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências, Universidade de Lisboa,
Tapada da Ajuda,
1349018
Lisboa,
Portugal
^{89}
Universidad Politécnica de Cartagena, Departamento de Electrónica y Tecnología de Computadoras,
30202
Cartagena,
Spain
^{90}
Kapteyn Astronomical Institute, University of Groningen,
PO Box 800,
9700 AV
Groningen,
The Netherlands
^{91}
Infrared Processing and Analysis Center, California Institute of Technology,
Pasadena,
CA
91125,
USA
^{92}
Institut d’Astrophysique de Paris,
98bis Boulevard Arago,
75014
Paris,
France
^{★★} Corresponding author; noemi.frusciante@unina.it
Received:
21
July
2023
Accepted:
29
January
2024
Context. The future Euclid space satellite mission will offer an invaluable opportunity to constrain modifications to Einstein’s general relativity at cosmic scales. In this paper, we focus on modified gravity models characterised, at linear scales, by a scaleindependent growth of perturbations while featuring different testable types of derivative screening mechanisms at smaller nonlinear scales.
Aims. We considered three specific models, namely JordanBransDicke, a scalartensor theory with a flat potential, the normal branch of DvaliGabadadzePorrati (nDGP) gravity, a braneworld model in which our Universe is a fourdimensional brane embedded in a fivedimensional Minkowski spacetime, and kmouflage gravity, an extension of kessence scenarios with a universal coupling of the scalar field to matter. In preparation for real data, we provide forecasts from spectroscopic and photometric primary probes by Euclid on the cosmological parameters and the additional parameters of the models, respectively, ω_{BD}, Ω_{гc} and ϵ_{2,0}, which quantify the deviations from general relativity. This analysis will improve our knowledge of the cosmology of these modified gravity models.
Methods. The forecast analysis employs the Fisher matrix method applied to weak lensing (WL); photometric galaxy clustering (GC_{ph}), spectroscopic galaxy clustering (GC_{sp}) and the crosscorrelation (XC) between GC_{ph} and WL. For the Euclid survey specifications, we define three scenarios that are characterised by different cuts in the maximum multipole and wave number, to assess the constraining power of nonlinear scales. For each model we considered two fiducial values for the corresponding model parameter.
Results. In an optimistic setting at 68.3% confidence interval, we find the following percentage relative errors with Euclid alone: for log_{10} ω_{BD}, with a fiducial value of ω_{BD} = 800, 27.1% using GC_{sp} alone, 3.6% using GC_{ph}+WL+XC and 3.2% using GC_{ph}+WL+XC+GC_{sp}; for log_{10} Ω_{гc}, with a fiducial value of Ω_{гc} = 0.25, we find 93.4, 20 and 15% respectively; and finally, for ϵ_{2,0} = −0.04, we find 3.4%, 0.15%, and 0.14%. From the relative errors for fiducial values closer to their ΛCDM limits, we find that most of the constraining power is lost. Our results highlight the importance of the constraining power from nonlinear scales.
Key words: cosmology: theory / largescale structure of Universe
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
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 largescale 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 scalartensor modification to GR, JordanBransDicke (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 largescale 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 smallscale experiments such as those in the Solar System where fifthforce 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 fifthforce effects on Solar System scales must shield smallscales 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), kmouflage (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 kmouflage 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 spectroscopic 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, kmouflage 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 scaleindependent linear growth rate. The three theories we consider are: JBD gravity (Brans & Dicke 1961), which is the simplest scalartensor theory and involves a nonminimal coupling between a scalar field and the metric; DvaliGabadadzePorrati 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 kmouflage model (KM, see Babichev et al. 2009; Brax & Valageas 2014, 2016), that is a scalartensor theory with a noncanonical kinetic energy and the kmouflage screening property. The JBD theory is not screened, and to be consistent with current astrophysical constraints on fifth forces, it requires that the BransDicke (BD) parameter ω_{BD} > 4 × 10^{4} (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 largescale deviations against which we compare the DGP and kmouflage 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 timedependent 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 scaleindependent modified gravity
We followed the Bardeen formalism (Bardeen 1980; Ma & Bertschinger 1995) and defined the infinitesimal line element of the flat FriedmannLemaîtreRobertsonWalker (FLRW) metric by $$\text{d}{s}^{2}=(1+2\text{\Psi}){c}^{2}\text{d}{t}^{2}+{a}^{2}(t)(12\text{\Phi}){\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^{i}. At linear order, the matter energymomentum 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{\Sigma}}_{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{\Psi}=\frac{4\pi {G}_{\text{N}}}{{c}^{2}}{a}^{2}\mu (a,k)\left[\overline{\rho}\text{\Delta}+3\left(\overline{\rho}+\overline{p}/{c}^{2}\right)\sigma \right],$$(3) $$\begin{array}{c}{k}^{2}(\text{\Phi}+\text{\Psi})=\frac{8\pi {G}_{\text{N}}}{{c}^{2}}{a}^{2}\{\text{\Sigma}(a,k)\left[\overline{\rho}\text{\Delta}+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 quasistatic approximation, which considers scales smaller than the horizon and the soundhorizon 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{\Sigma}=\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 scaleindependent (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 EinsteinBoltzmann 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 quasistatic (QS) limit, that is, scales sufficiently small to be well within the horizon and the soundhorizon 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 JordanBransDicke 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{\Lambda}\right)+{\mathcal{L}}_{\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 longwavelength, lowenergy limit of more general scalartensor 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 scalartensor theories, and furthermore, be independent from modelspecific constraints on smallscales. Thus, while simple, JBD gravity is a powerful tool for constraining general classes of scalartensor 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{\Lambda})\right].\end{array}$$(7)
Here, T_{µν} is the total matter stressenergy tensor, while □ denotes the d’Alembertian. On the background, these equations give $$3{H}^{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{\Lambda}}{\varphi},$$(8) $$2\dot{H}+3{H}^{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 stressenergy tensor on the background.
The phenomenological QS functions in this theory read $$\begin{array}{l}\text{\Sigma}=\frac{1}{{G}_{\text{N}}\varphi},\hfill \\ \mu =\frac{4+2{\omega}_{\text{BD}}}{3+2{\omega}_{\text{BD}}}\text{\Sigma},\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_{10} (ω_{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_{10} (1 + 1/ω_{BD}) < 0.0030 at 95% CL combining Planck 2015 and BOSS DR1011 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_{10} (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 KiloDegree 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_{10} (1/ω_{BD}), log_{10} (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 Cavendishlike 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 DvaliGabadadzePorrati braneworld gravity
The DGP model (Dvali et al. 2000) is a fivedimensional braneworld model defined by the action $$S=\frac{{c}^{4}}{16\pi {G}_{5}}\underset{\mathcal{M}}{{\displaystyle \int}}{\text{d}}^{5}x\sqrt{\gamma}{R}_{5}+\underset{\partial \mathcal{M}}{{\displaystyle \int}}{\text{d}}^{4}x\sqrt{g}\left[\frac{{c}^{4}}{16\pi {G}_{\text{N}}}R+{\mathcal{L}}_{\text{m}}\right],$$(13)
where γ is the fivedimensional metric and R_{5} its Ricci curvature scalar. G_{5} and G_{N} are the five and fourdimensional Newton constants, respectively. The matter Lagrangian is denoted with 𝓛_{m} and is confined to a fourdimensional brane in a fivedimensional Minkowski spacetime. The induced gravity described by the usual fourdimensional EinsteinHilbert action is responsible for the recovery of the fourdimensional gravity on the brane. The crossover scale r_{c} = G_{5}/(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 selfaccelerating 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 fivedimensional gravity, respectively. The selfaccelerating branch admits the latetime 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}{1g},$$(15)
where, using the quasistatic (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}}{3{H}^{2}}\right)\right]}^{1},$$(16)
so that the effective modification introduced with η can be treated as scaleindependent.
The effective 3 + 1 Poisson equation for the lensing potential in the QS approximation is $${k}^{2}(\text{\Phi}+\text{\Psi})=\frac{8\pi {G}_{\text{N}}}{{c}^{2}}{a}^{2}\overline{\rho}\text{\Delta},$$(17)
therefore (Lombriser et al. 2013) $$\text{\Sigma}=1,$$(18)
and hence $$\mu (a)=1+\frac{1}{3\beta},\text{\hspace{0.17em}}\beta (a)\equiv 1+\frac{H}{{H}_{0}}\frac{1}{\sqrt{{\text{\Omega}}_{\text{rc}}}}\left(1+\frac{\dot{H}}{3{H}^{2}}\right).$$(19)
We chose to parametrise the modification to gravity by ${\mathrm{\Omega}}_{\text{rc}}\equiv {c}^{2}/\left(4{r}_{\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 secondorder moments of the correlation function from SDSS DR7 data up to r_{max} = 120 Mpc h^{−1}, and marginalised bias, found an Ω_{rc} upper limit at 95% of ~ 40 (from r_{c} > 340 Mpc with fixed H_{0} = 70 km s^{−1} Mpc^{−1}). 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_{c}H_{0}/c]^{−1} < 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 Euclidlike 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 threedimensional 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 largescale structure data, where the constraints were mainly driven by the CMB data. Thus, highprecision 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 ΛCDMbackground 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 kmouflage gravity
KM theories are built complementing simple kessence 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}{\mathcal{M}}^{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, ℳ^{4} 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{\mathcal{M}}^{4}},$$(21)
and ℳ^{4}K is the nonstandard 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 timedependent functions, related to the coupling A and the kinetic function K, that is, $${\u03f5}_{2}\equiv \frac{\text{d\hspace{0.17em}}\mathrm{ln}\overline{A}}{\text{d\hspace{0.17em}}\mathrm{ln}a},\text{\hspace{0.17em}}{\u03f5}_{1}\equiv \frac{2}{{\overline{K}}^{\prime}}{\left[{\u03f5}_{2}\frac{c}{\sqrt{8\pi {\tilde{G}}_{\text{N}}}}{\left(\frac{\text{d}\overline{\phi}}{\text{d\hspace{0.17em}}\mathrm{ln}a}\right)}^{1}\right]}^{2},$$(22)
where overbars 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{\u03f5}_{2}\right)}^{2}}\left[\overline{\rho}+\frac{{\mathcal{M}}^{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 (Benevento et al. 2019) $$\mu (a)=\left(1+{\u03f5}_{1}\right){\overline{A}}^{2},\text{\hspace{0.17em}}\Sigma (a)={\overline{A}}^{2}.$$(24)
We recover GR when ϵ_{1} → 0, ϵ_{2} → 0, and A and K are constants, with ℳ^{4}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 timedependent background values $\overline{K}(a)$ and $\overline{A}(a)$. Brax & Valageas (2016) proposed a simple parameterisation that satisfies selfconsistency constraints. It involves five parameters, {ϵ_{2,0}, γ_{A},m,α_{U},γ_{U}}, where ϵ_{2,0} is the value of ϵ_{2} 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(m1)}{2m1},\text{\hspace{0.17em}}{\alpha}_{A}=\frac{{\u03f5}_{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 highredshift dependence of $\overline{A}(a)$, while γ_{A} sets the redshift at which $\overline{A}$ changes from the current unit value to the highredshift 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}}},\text{\hspace{0.17em}}\sqrt{\tilde{\chi}}=\frac{{\overline{\rho}}_{0}}{{\mathcal{M}}^{4}}\frac{{\u03f5}_{2}{\overline{A}}^{4}}{2U\left(3{\u03f5}_{2}+\frac{\text{d\hspace{0.17em}}\mathrm{ln}U}{\text{d\hspace{0.17em}}\mathrm{ln}a}\right)},$$(27)
which we took to be of the form $$U(a)\propto \frac{{a}^{2}\mathrm{ln}\left({\gamma}_{U}+a\right)}{\left(\sqrt{{a}_{\text{eq}}}+\sqrt{a}\right)\mathrm{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 darkenergy eras. The parameters α_{U} and γ_{U}, of order unity, set the shape of the transition to the darkenergy 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_{0} (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 Xray 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 ϵ_{2} we have ϵ_{1} ≃ –ϵ_{2}. This shows that the running of Newton’s constant and the impact on the linear metric and matter density fluctuations are set by ϵ_{2}. 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,\text{\hspace{0.17em}}{\u03f5}_{2}(a)\to 0,\text{\hspace{0.17em}}{\u03f5}_{1}(a)\to 0,\text{\hspace{0.17em}}\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 flatsky limit. The spectra are then given by $${C}_{ij}^{XY}(\ell )=c{{{\displaystyle \int}}^{\text{}}}_{{\u0240}_{\mathrm{min}}}^{{\u0240}_{\mathrm{max}}}\text{d}\u0240\frac{{W}_{i}^{X}(\u0240){W}_{j}^{Y}(\u0240)}{H(\u0240){r}^{2}(\u0240)}{P}_{\delta \delta}\left({k}_{\ell},\u0240\right),$$(30)
with k_{ℓ} = (ℓ + 1/2)/r(z), r(z) the comoving distance to redshift z = 1/a – 1, and P_{δδ} (k_{ℓ, z}) the nonlinear 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 nonvanishing. The GC_{ph} and WL window functions read (Spurio Mancini et al. 2019) $${W}_{i}^{\text{G}}(k,\u0240)={b}_{i}(k,\u0240)\frac{{n}_{i}(\u0240)}{{\overline{n}}_{i}}\frac{H(\u0240)}{c},$$(31) $$\begin{array}{l}{W}_{i}^{\text{L}}(k,\u0240)=\frac{3}{2}{\text{\Omega}}_{\text{m},0}\frac{{H}_{0}^{2}}{{c}^{2}}(1+\u0240)r(\u0240)\text{\Sigma}(\u0240){{{\displaystyle \int}}^{\text{}}}_{\u0240}^{\u0240\mathrm{max}}\text{d}{\u0240}^{\prime}\frac{{n}_{i}\left({\u0240}^{\prime}\right)}{{\overline{n}}_{i}}\frac{r\left({\u0240}^{\prime}\u0240\right)}{r\left({\u0240}^{\prime}\right)}\\ \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}+{W}_{i}^{\text{IA}}(k,\u0240),\end{array}$$(32)
where ${n}_{i}(\u0240)/{\overline{n}}_{i}$ and b_{i}(k, ɀ) are the normalised galaxy distribution and the galaxy bias in the ith redshift bin, respectively, and ${W}_{i}^{\text{IA}}(k,\u0240)$ 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,\u0240)=\frac{{\mathcal{A}}_{\text{IA}}{C}_{\text{IA}}{\text{\Omega}}_{\text{m},0}{\mathcal{F}}_{\text{IA}}(\u0240)}{\delta (k,\u0240)/\delta (k,\u0240=0)}\frac{{n}_{i}(\u0240)}{{\overline{n}}_{i}(\u0240)}\frac{H(\u0240)}{c},$$(33)
where $${\mathcal{F}}_{\text{IA}}(\u0240)={(1+\u0240)}^{{\eta}_{\text{IA}}}{\left[\frac{\langle L\rangle (\u0240)}{{L}_{\star}(\u0240)}\right]}^{{\beta}_{\text{IA}}},$$(34)
with 〈L〉(z) and L_{*}(ɀ) redshiftdependent 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 quasistatic 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};\u0240\right)=\hfill & \frac{1}{{q}_{\perp}^{2}(\u0240){q}_{\parallel}(\u0240)}\left\{\frac{{\left[b{\sigma}_{8}(\u0240)+f{\sigma}_{8}(\u0240){\mu}_{\theta}^{2}\right]}^{2}}{1+{\left[f(\u0240)k{\mu}_{\theta}{\sigma}_{\text{p}}(\u0240)\right]}^{2}}\right\}\hfill \\ & \times \frac{{P}_{\text{dw}}\left(k,{\mu}_{\theta};\u0240\right)}{{\sigma}_{8}^{2}(\u0240)}{F}_{\u0240}\left(k,{\mu}_{\theta};\u0240\right)+{P}_{\text{s}}(\u0240),\hfill \end{array}$$(35)
where P_{dw}(k,µ_{θ}; z), the dewiggled power spectrum, includes the correction that accounts for the smearing of the BAO features, $${P}_{\text{dw}}\left(k,{\mu}_{\theta};\u0240\right)={P}_{\delta \delta}^{\mathrm{lin}}(k;\u0240){\text{e}}^{{g}_{\mu}{k}^{2}}+{P}_{\text{nw}}(k;\u0240)\left(1{\text{e}}^{{g}_{\mu}{k}^{2}}\right),$$(36)
and ${P}_{\delta \delta}^{\text{lin}}(k;\u0240)$ stands for the linear matter power spectrum. P_{nw}(k; ɀ) is a nowiggle power spectrum with the same broadband shape as ${P}_{\delta \delta}^{\text{lin}}(k;\u0240)$ but without BAO features (see the discussion below). The function $${g}_{\mu}\left(k,{\mu}_{\theta},\u0240\right)={\sigma}_{\text{v}}^{2}(\u0240)\left\{1{\mu}_{\theta}^{2}+{\mu}_{\theta}^{2}{[1+f(\u0240)]}^{2}\right\},$$(37)
is the nonlinear damping factor of the BAO signal derived by Eisenstein et al. (2007), with $${\sigma}_{\text{v}}^{2}(\u0240)=\frac{1}{6{\pi}^{2}}{\displaystyle \int \text{d}}k{P}_{\delta \delta}^{\mathrm{lin}}(k,\u0240).$$(38)
The curly bracket in Eq. (35) is the redshiftspacedistortion (RSD) contribution correcting for the nonlinear fingerofGod (FoG) effect, where we defined bσ_{8}(ɀ) as the product of the effective scaleindependent bias of galaxy samples and the r.m.s. matter density fluctuation σ_{8}(ɀ) (we marginalized over bσ_{8}(ɀ)), while µ_{θ} is the cosine of the angle θ between the wave vector k and the lineofsight direction $\widehat{r}$ and ${\sigma}_{\text{p}}^{2}(\u0240)={\sigma}_{\text{v}}^{2}(\u0240)$. Although these parameters were assumed to be the same, they come from two different physical effects, namely largescale 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}_{\u0240}\left(k,{\mu}_{\theta};\u0240\right)={\text{e}}^{{k}^{2}{\mu}_{\theta}^{2}{\sigma}_{r}^{2}(\u0240)},$$(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 AlcockPaczynski effect, which is parametrised in terms of the angular diameter distance D_{A}(ɀ) and the Hubble parameter H(ɀ) as $${q}_{\perp}(\u0240)=\frac{{D}_{\text{A}}(\u0240)}{{D}_{\text{A},\mathrm{ref}}(\u0240)},$$(40) $${q}_{\Vert}(\u0240)=\frac{{H}_{\text{ref}}(\u0240)}{H(\u0240)}.$$(41)
Due to the AlcockPaczynski effect, µ_{θ} and k were rescaled in a cosmologydependent 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 scaleindependent shotnoise 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 redshiftdependent parameters H(ɀ), D_{A}(ɀ), and fσ_{8}(ɀ) 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_{0}w_{a}CDM models. The other term appearing in Eq. (36) is the nonwiggle matter power spectrum, and this was obtained applying a SavitzkyGolay filter to the matter power spectrum ${P}_{\delta \delta}^{\text{lin}}(k;\u0240)$. For more details on this implementation, we refer to EC19.
3.3 nonlinear modelling for weak lensing
The galaxy power spectrum on mildly nonlinear 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 nonlinear 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 nonlinear 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 Euclidlike 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 nonlinear matter power spectrum. While having a common prescription at nonlinear scales for all models would be desirable, given the different behaviours of each model at nonlinear 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 JordanBransDicke gravity
For the nonlinear 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 nonlinear 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 nonlinear matter power spectrum can be written as a sum P_{NL}(k, ɀ) = P_{1H}(k, ɀ) + P_{2H}(k, ɀ) where $${P}_{2\text{H}}(k,\u0240)={P}_{\text{lin\hspace{0.17em}}}(k,\u0240){\left[{\displaystyle {\int}_{0}^{\infty}b}(M,\u0240)W(M,k,\u0240)n(M,\u0240)\text{d}M\right]}^{2},$$(42)
is the socalled twohalo term (correlation between different haloes), and $${P}_{1\text{H}}(k,\u0240)={\displaystyle {\int}_{0}^{\infty}{W}^{2}}(M,k,\u0240)n(M,\u0240)\text{d}M,$$(43)
is the socalled onehalo term (correlations between masselements 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,\u0240)=\frac{1}{\overline{\rho}}{\displaystyle {\int}_{0}^{{r}_{\text{v}}}4}\pi {r}^{2}\frac{\mathrm{sin}(kr)}{kr}\rho (M,r,\u0240)\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}(\u0240)\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 Nbody 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 secondorder 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 smallscales, up to k < 10 h Mpc^{–1}, the RAMSES gridbased 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}}{(\u0240)}^{0.352}\\ \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}\times \left\{{d}_{0}+\left(418.0{d}_{0}\right)\mathrm{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_{0} = 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^{4} ≳ ω_{BD} ≥ 50 and on scales up to k ≃ 10 h Mpc^{–1}. This nonlinear prescription was made with past weaklensing 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 nonlinear (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 σ_{8} normalisation today.
As was shown in detail in Joudaki et al. (2022), the linear growth rate will undergo a scaleindependent enhancement (without altering the shape of the linear power spectrum) while the nonlinear 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 nonlinear power spectrum shows a small suppression on smallscales that sets in at about k = 0.3 h Mpc^{−1}. It is important to note that, unlike nDGP, the JBD theory is not endowed with gravitational screening, so that any modification on smallscales in the nonlinear power spectrum is inherited from the change in the primordial amplitude, A_{s}, in the linear power spectrum. In order for σ_{8} to be the same in ΛCDM as in JBD we need a higher primordial amplitude As . A higher A_{s} in ΛCDM means that nonlinear 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 smallscale 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 smallscales. 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 nonlinear 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 Largescale 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 nonlinear (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 nonlinear (solid) cases. Bottom left: cosmic shear (WL) angular power spectra for the autocorrelation of the first photometric bin ${C}_{1,1}^{\text{WL}}(\ell )$ (solid line) and the crosscorrelation 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 C_{i,j}(ℓ) to their ΛCDM counterpart, for the bin combination i = j = 1 and i = 1; j = 10. 
3.3.2 DvaliGabadadzePorrati gravity
We modelled the nonlinear power spectrum using the halo model reaction (Cataneo et al. 2019) which has been shown to agree with Nbody simulations at the 2% level down to k = 3h Mpc^{−1} with small variation depending on redshift, degree of modification to GR and mass of neutrinos (Bose et al. 2021). The approach attempts to model nonlinear corrections to the power spectrum from modified gravity through the socalled reaction 𝓡(k, ɀ), which employs both oneloop perturbation theory and the halo model. The nonlinear power spectrum in nDGP is given by the product $${P}_{\text{NL}}(k,\u0240)=\mathcal{R}(k,\u0240){P}_{\text{NL}}^{\text{pseudo}}(k,\u0240),$$(46)
where the pseudopower spectrum is a spectrum in which all nonlinear 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 nonlinear clustering as described by GR, in line with our definition of the pseudopower spectrum.
The halo model reaction, 𝓡(k, ɀ), is given by a corrected ratio of target to pseudohalo model spectra, $$\mathcal{R}(k,\u0240)=\frac{\left\{[1\mathcal{E}(\u0240)]{e}^{k/{k}_{\star}(\u0240)}+\mathcal{E}(\u0240)\right\}{P}_{2\text{H}}(k,\u0240)+{P}_{1\text{H}}(k,\u0240)}{{P}_{\text{hm}}^{\text{pesudo\hspace{0.17em}}}(k,\u0240)}.$$(47)
The components are given explicitly as $${P}_{\text{hm}}^{\text{pseudo\hspace{0.17em}}}(k,\u0240)={P}_{2\text{H}}(k,\u0240)+{P}_{1\text{H}}^{\text{pseudo\hspace{0.17em}}}(k,\u0240),$$(48) $$\mathcal{E}(\u0240)=\underset{k\to 0}{\mathrm{lim}}\frac{{P}_{1\text{H}}(k,\u0240)}{{P}_{1\text{H}}^{\text{pseudo\hspace{0.17em}}}(k,\u0240)},$$(49) $${k}_{\star}(\u0240)=\overline{k}{\left\{\mathrm{ln}\left[\frac{A(\overline{k},\u0240)}{{P}_{2\text{H}}(\overline{k},\u0240)}\mathcal{E}(\u0240)\right]\mathrm{ln}[1\mathcal{E}(\u0240)]\right\}}^{1},$$(50)
where $$A(k,\u0240)=\frac{{P}_{1\text{\u2212loop\hspace{0.17em}}}(k,\u0240)+{P}_{1\text{H}}(k,\u0240)}{{P}_{1\text{\u2212loop\hspace{0.17em}}}^{\text{pseudo\hspace{0.17em}}}(k,\u0240)+{P}_{1\text{H}}^{\text{pseudo\hspace{0.17em}}}(k,\u0240)}{P}_{\text{hm}}^{\text{pseudo\hspace{0.17em}}}(k,\u0240){P}_{1\text{H}}(k,\u0240)\text{.\hspace{0.17em}}$$(51)
P_{շH}(k, ɀ) is the twohalo 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 onehalo 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 pseudocosmology has no nonlinear 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 oneloop predictions with and without nonlinear modifications to ΛCDM, respectively. As in previous works, the limit in Eq. (49) was taken to be at k = 0.01 h Mpc^{−1}, and we computed k_{⋆} using $\overline{k}=0.06h\text{\hspace{0.17em}}{\text{Mpc}}^{1}$. This scale was chosen such that the oneloop predictions are sufficiently accurate at all the redshifts considered.
We note that the correction to the halomodel 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 halomodel calculations, we used a ShethTormen mass function (Sheth & Tormen 1999, 2002), a powerlaw concentrationmass relation (see for example Bullock et al. 2001) and an NFW halodensity 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 σ_{8} normalisation today. For nDGP1, the linear matter power spectrum is identical to that of the ΛCDM. However, the nonlinear matter power spectrum shows a suppression starting from scales k = 0.1 h Mpc^{−1} , which reaches about 10% at scales of about k = 1.0 h Mpc^{−1}. 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 σ_{8} 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 onehalo terms at k > 0.1 h Mpc^{−1} (Cataneo et al. 2019). In a nDGP cosmology, there will be more highmass halos when compared to its GR counterpart due to the supplemental fifth force sourced by the additional degree of freedom. Highmass contributions to the onehalo 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 1halo GR power spectrum to be larger than the DGP onehalo 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 onehalo 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 latetime amplitude of linear perturbations and one with the same earlytime amplitude, respectively.
3.3.3 kmouflage gravity
For the KM nonlinear matter power spectrum, we used a similar analytical approach that combined oneloop perturbation theory with a halomodel. 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 onehalo contributions, as in Eq. (48), but it uses a Lagrangian framework to evaluate these two contributions in terms of the pairseparation probability distribution. This provides a Lagrangianspace regularisation of perturbation theory, which matches standard perturbation theory up to oneloop order and includes a partial resummation of higherorder terms, such that the pairseparation probability distribution is positive and normalised to unity at all scales. Thus, the twohalo term goes beyond the Zeldovich approximation by including a nonzero skewness, which enables the consistency with standard perturbation theory up to oneloop order, as well as a simple Ansatz for higherorder cumulants inspired by the adhesion model. The onehalo term includes a counterterm 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 tophat of Lagrangian radius q_{M}, and f (ν) is the normalised halo mass function. This counterterm is not introduced by hand, as it arises directly from the splitting of the matter power spectrum over twohalo and onehalo 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 nonlinear 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 nonlinear screening. However, for the models that we considered here, screening is negligible beyond 1 h^{−1} Mpc (clusters are not screened), so that δ_{c} is independent of M.
This nonlinear modelling was extended in Brax & Valageas (2013) to several modifiedgravity 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 oneloop contributions to the matter power spectrum, which enter the two halo term, and on the spherical collapse dynamics, which enter the onehalo term through the halo mass function. Therefore, this nonlinear modelling exactly captures the modification of gravity both up to oneloop order and at the level of the fully nonlinear spherical collapse. The Ansatz for the partial resummation of higherorder 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^{−5} and 10^{−6}), shows that this approach captures the relative deviation from the ΛCDM power spectrum at ɀ = 0 up to k = 3 h Mpc^{−1}. We have not yet compared this recipe with numerical simulations of the KM theories, which have only recently been performed (HernándezAguayo 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 σ_{8} 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^{−2} < k < 1 h Mpc^{−1}): 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 matterradiation 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 σ_{8} 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 nonlinear 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 є_{2} = −0.04 we have є_{1} ≃ 0.12, whereas at linear order in the scalar field dynamics we would have expected that є_{1} = −є_{2}. As the growth of structures is sensitive to є_{1} , this leads to a stronger effect than could be expected from the value of є_{2} 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^{−1}. Nonlinear 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 nonlinear scales. This is due to the loss of information in the power spectrum in the nonlinear regime, which can be related to the universal NFW profile of virialised haloes within the halo model. Then, the shape of the nonlinear power spectrum at those scales is mostly set by the scale where the nonlinear 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 nonlinear 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 nonlinear 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{\Omega}}_{\text{m},0}=0.315,\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{\text{\Omega}}_{\text{b},0}=0.049,\hfill \\ \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}h=0.6737,\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{n}_{\text{s}}=0.966,\hfill \end{array}$$(54)
where Ω_{m,0} = Ω_{c,0} + Ω_{b,0} + Ω_{ν,0}. To fix the fiducial value for σ_{8}, we used the same initial amplitude of primordial perturbations for all models, namely A_{s} = 2.09681 × 10^{−9}. 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 σ_{8} value for each model along with the fiducial values of the model parameters,

JBD $$\begin{array}{ll}{\Theta}_{\text{fid},1}=\left\{{\sigma}_{8}=0.816,{\mathrm{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,{\mathrm{log}}_{10}{\omega}_{\text{BD}}=3.39794\right\}\hfill & \text{\hspace{0.17em}(JBD2)}.\text{\hspace{0.17em}}\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}.

nDGP $$\begin{array}{l}{\Theta}_{\text{fid},1}=\left\{{\sigma}_{8}=0.8690,{\mathrm{log}}_{10}{\mathrm{\Omega}}_{\text{rc}}=0.60206\right\}\text{\hspace{0.17em}}(\text{nDGP}1);\\ {\Theta}_{\text{fid},2}=\left\{{\sigma}_{8}=0.8105,{\mathrm{log}}_{10}{\mathrm{\Omega}}_{\text{rc}}=6\right\}\text{\hspace{0.17em}\hspace{0.17em}(nDGP2)}.\text{\hspace{0.17em}}\end{array}$$(56)
These values correspond to Ω_{rc} = 0.25 and Ω_{rc} = 10^{−6}, for nDGP1 and nDGP2, respectively.

KM $$\begin{array}{l}{\Theta}_{\text{fid},1}=\left\{{\sigma}_{8}=0.994,{\u03f5}_{2}=0.04\right\}\text{\hspace{0.17em}}(\text{KM}1);\hfill \\ {\Theta}_{\text{fid},2}=\left\{{\sigma}_{8}=0.813,{\u03f5}_{2}=0.0001\right\}\text{\hspace{0.17em}}(\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 equipopulated redshift bins whose limits were obtained from the redshift distribution $$n(\u0240)\propto {\left(\frac{\u0240}{{\u0240}_{0}}\right)}^{2}\mathrm{exp}\left[{\left(\frac{\u0240}{{\u0240}_{0}}\right)}^{3/2}\right],$$(58)
with ${\u0240}_{0}=0.9/\sqrt{2}$ and the normalisation set by the requirement that the surface density of galaxies is ${\overline{n}}_{g}=30\text{\hspace{0.17em}}{\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{\u0240}}_{i}}$, and ${\overline{\u0240}}_{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 Gaussianonly 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{\Delta}\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 uppercase (lowercase) 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 equispaced 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}_{\u03f5}^{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}_{\u03f5}^{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 redshiftdependent parameters ln bσ_{8}(ɀ_{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 quasilinear, for which we fixed the maximum wavenumber to k_{max} = 0.15 h Mpc^{−1}. 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. Nonlinear corrections begin to play a role for wave numbers larger than k = 0.1 h Mpc^{−1} for the redshifts we considered (see Taruya et al. 2010). Cutting at a lower k_{max} should avoid any bias induced by neglecting nonlinear 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.
Euclid survey specifications for WL, GC_{ph} and GC_{sp} .
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 quasilinear, 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 nonzero 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_{10} 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_{10} 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 nonlinear 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}_{{\mathrm{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_{10} ω_{BD} for JBD, log_{10} Ω_{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^{−1}) in the nonlinear 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 JordanBransDicke gravity
There are remarkably strong noncosmological constraints on JBD gravity, with current bounds of about ω_{BD} > 10^{5}. As pointed out above, however, it is useful to consider JBD gravity as a longwavelength limit of more general scalartensor theories with a wide range of smallscale 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^{3}, independent of what we know from non cosmological systems. Two further practical reasons led us to choose ω_{BD} ~ 10^{3} 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 quasilinear, 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_{10} ω_{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_{10} ω_{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_{10} ω_{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 GC_{sp} in the optimistic (dotted area), pessimistic (dashed area) and quasilinear (filled area) scenarios (top panel), for WL+XC+GC_{ph} in the optimistic and pessimistic cases (central panel) and for the GC_{Sp}+WL+XC+GC_{ph} in the optimistic and pessimistic cases. 
5.2 DvaliGabadadzePorrati gravity
The two cases considered for the nDGP model are log_{10} Ω_{rc} = −0.60206 (nDGPl) and log_{10} Ω_{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 quasilinear 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_{10} Ω_{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 σ_{8} 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_{10} Ω_{rc}−n_{s} panel that the two parameters are anticorrelated in the nDGP1 case. This makes the forecast errors on n_{s} larger than the nDGP2 case where the anticorrelation disappears. A similar discussion can be made for the σ_{8} 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_{10} Ω_{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 nonlinear information. To examine this, we considered the linear power spectrum: We find that the derivative with respect to log_{10} Ω_{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_{10} Ω_{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 quasilinear) scenario, leads to larger uncertainties on the model parameter (log_{10} Ω_{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 nonlinear 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.
Forecast 1σ marginal errors on the cosmological and model parameters relative to their corresponding fiducial value for JBD1 and JBD2 in the pessimistic, quasilinear and optimistic cases, using Euclid observations GC_{sp}, WL+XC+GC_{ph} and GC_{sp}+WL+XC+GC_{ph}.
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 GC_{sp}, in blue, we plot WL+XC+GC_{ph} and in green we plot GC_{sp}+WL+XC+GC_{ph}. 
Forecast 1σ marginal errors on the cosmological and model parameters relative to their corresponding fiducial value for nDGPl and nDGP2 in the pessimistic, quasilinear and optimistic cases, using Euclid observations of GC_{sp}, WL+XC+GC_{ph} and GC_{sp}+WL+XC+GC_{ph}.
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 GC_{sp}, in blue, we plot WL+XC+GC_{ph} and in green we plot GC_{sp}+WL+XC+GC_{ph}. 
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 GC_{sp}, in blue, we plot WL+XC+GC_{ph} and in green, we plot GC_{sp}+WL+XC+GC_{ph}. 
5.3 kmouflage
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 modifiedgravity 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^{−3} for the GC_{sp} alone;
upper bound at 7 × 10^{−4} when considering the combination WL+XC+GC_{ph};
upper bound at 4 × 10^{−4} 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 worstcase scenario of cutting the observed galaxy power spectrum at k_{max} = 0.15 h Mpc^{−1}, 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 kmouflage 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 ϵ_{20} in either the pessimistic and optimistic scenarios.
The percentage error on ϵ_{20} in the ΛCDMproximate 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 bestcase scenario. In the optimistic case, the Euclid survey will therefore be able to detect a deviation from λCDM only when ϵ_{20} ≳ 4 × 10^{−4}, 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 anticorrelated 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}.
Forecast 1σ marginal errors on the cosmological and model parameters, relative to their corresponding fiducial value, for KM1 and KM2 in the pessimistic, quasilinear and optimistic cases, using Euclid observations of GC_{sp}, WL+XC+GC_{ph} and GC_{sp}+WL+XC+GC_{ph}.
6 Conclusions
We have explored the constraining power of the future Euclid mission for linearly scaleindependent extensions of the concordance cosmological model, that is, models that induce modifications to the linear growth of perturbations that are solely timedependent while featuring different testable types of screening mechanisms at smaller nonlinear scales. We considered three specific models, namely JBD, a scalartensor 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 fourdimensional brane embedded in a fivedimensional Minkowski spacetime (Sect. 2.2); and KM gravity, an extension of kessence 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_{10} ω_{BD} for JBD, log_{10}Ω_{rc} for nDGP, and ϵ_{2} 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 crosscorrelation (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 nonlinear matter power spectrum using different prescriptions for each MG model: The HMCODE (Mead et al. 2015, 2016) calibrated on a suite of Nbody 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 oneloop 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 nonlinear scales: the quasilinear scenario with k_{max} = 0.15 h Mpc^{−1} for GC_{sp}; the pessimistic scenario with k_{max} = 0.25 h Mpc^{−1} for GC_{sp}, ℓ_{max} = 1500 for WL, and ℓ_{max} = 750 for GC_{ph}; and the optimistic scenario with k_{max} = 0.3 h Mpc^{−1} 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 nonzero modifiedgravity 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^{−6}) 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^{−5}. 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 ϵ_{20} = −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^{−4} in the optimistic case and ϵ_{20} < 2 × 10^{−3} 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 lowredshift data Joudaki et al. (2022).
Finally, our analysis clearly showed that most of the constraining power of Euclid comes from the nonlinear 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 nonlinearities 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 scaleindependent 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 crosscorrelation 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 crosscorrelated 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.euclidec.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 SklodowskaCurie 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/FISAST/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 largescale structure observables” and acknowledges support from the Spanish Ministry of Science, Innovation and Universities through grant ESP201789838, 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. 202one43HH.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 GC_{sp}, in blue, we plot WL+XC+GC_{ph} and in green, we plot GC_{sp}+WL+XC+GC_{ph}. 
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 GC_{sp}, in blue, we plot WL+XC+GC_{ph} and in green, we plot GC_{sp}+WL+XC+GC_{ph}. 
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 GC_{sp}, in blue, we plot WL+XC+GC_{ph} and in green, we plot GC_{sp}+WL+XC+GC_{ph}. 
References
 Aghanim, N., Akrami, Y., Ashdow, Y., et al. 2020, A&A, 641, A6; Erratum: 2021, A&A, 652, C4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Avilez, A., & Skordis, C. 2014, Phys. Rev. Lett., 113, 011101 [CrossRef] [Google Scholar]
 Babichev, E., Deffayet, C., & Ziour, R. 2009, Int. J. Mod. Phys. D, 18, 2147 [NASA ADS] [CrossRef] [Google Scholar]
 Bag, S., Mishra, S. S., & Sahni, V. 2018, Phys. Rev. D, 97, 123537 [NASA ADS] [CrossRef] [Google Scholar]
 Ballardini, M., Finelli, F., Umiltà, C., & Paoletti, D. 2016, JCAP, 05, 067 [CrossRef] [Google Scholar]
 Ballardini, M., Sapone, D., Umiltà, C., et al. 2019, JCAP, 05, 049 [CrossRef] [Google Scholar]
 Ballardini, M., Braglia, M., Finelli, F., et al. 2020, JCAP, 10, 044 [CrossRef] [Google Scholar]
 Ballardini, M., Finelli, F., & Sapone, D. 2022, JCAP, 06, 004 [CrossRef] [Google Scholar]
 Bardeen, J. M. 1980, Phys. Rev. D, 22, 1882 [NASA ADS] [CrossRef] [Google Scholar]
 Barreira, A., Brax, P., Clesse, S., et al. 2015, Phys. Rev. D, 91, 063528 [NASA ADS] [CrossRef] [Google Scholar]
 Barreira, A., Sánchez, A. G., & Schmidt, F. 2016, Phys. Rev. D, 94, 084022 [NASA ADS] [CrossRef] [Google Scholar]
 Bellini, E., & Sawicki, I. 2014, JCAP, 07, 050 [CrossRef] [Google Scholar]
 Bellini, E., Barreira, A., Frusciante, N., et al. 2018, Phys. Rev. D, 97, 023520 [NASA ADS] [CrossRef] [Google Scholar]
 Bellini, E., Sawicki, I., & Zumalacárregui, M. 2020, JCAP, 02, 008 [CrossRef] [Google Scholar]
 Benevento, G., Raveri, M., Lazanu, A., et al. 2019, JCAP, 05, 027 [CrossRef] [Google Scholar]
 Berti, E., Barausse, E., Cardoso, V., et al. 2015, Class. Quant. Grav., 32, 243001 [NASA ADS] [CrossRef] [Google Scholar]
 Bertotti, B., Iess, L., & Tortora, P. 2003, Nature, 425, 374 [Google Scholar]
 Blas, D., Lesgourgues, J., & Tram, T. 2011, JCAP, 07, 034 [CrossRef] [Google Scholar]
 Boisseau, B., EspositoFarese, G., Polarski, D., & Starobinsky, A. A. 2000, Phys. Rev. Lett., 85, 2236 [NASA ADS] [CrossRef] [Google Scholar]
 Bose, B., Cataneo, M., Tröster, T., et al. 2020, MNRAS, 498, 4650 [Google Scholar]
 Bose, B., Wright, B. S., Cataneo, M., et al. 2021, MNRAS, 508, 2479 [NASA ADS] [CrossRef] [Google Scholar]
 Bowcock, P., Charmousis, C., & Gregory, R. 2000, Class. Quant. Grav., 17, 4745 [NASA ADS] [CrossRef] [Google Scholar]
 Braglia, M., Ballardini, M., Emond, W. T., et al. 2020, Phys. Rev. D, 102, 023529 [NASA ADS] [CrossRef] [Google Scholar]
 Braglia, M., Ballardini, M., Finelli, F., & Koyama, K. 2021, Phys. Rev. D, 103, 043528 [NASA ADS] [CrossRef] [Google Scholar]
 Brans, C., & Dicke, R. H. 1961, Phys. Rev., 124, 925 [NASA ADS] [CrossRef] [Google Scholar]
 Brax, P., & Valageas, P. 2013, Phys. Rev. D, 88, 023527 [NASA ADS] [CrossRef] [Google Scholar]
 Brax, P., & Valageas, P. 2014, Phys. Rev. D, 90, 023508 [NASA ADS] [CrossRef] [Google Scholar]
 Brax, P., & Valageas, P. 2016, JCAP, 01, 020 [CrossRef] [Google Scholar]
 Brax, P., Rizzo, L. A., & Valageas, P. 2015, Phys. Rev. D, 92, 043519 [NASA ADS] [CrossRef] [Google Scholar]
 Brax, P., Casas, S., Desmond, H., & Elder, B. 2021, Universe, 8, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Bullock, J. S., Kolatt, T. S., Sigad, Y., et al. 2001, MNRAS, 321, 559 [Google Scholar]
 Cardoso, A., Koyama, K., Seahra, S. S., & Silva, F. P. 2008, Phys. Rev. D, 77, 083512 [NASA ADS] [CrossRef] [Google Scholar]
 Carroll, S. M., Duvvuri, V., Trodden, M., & Turner, M. S. 2004, Phys. Rev. D, 70, 043528 [NASA ADS] [CrossRef] [Google Scholar]
 Casas, S., Cardone, V. F., Sapone, D., et al. 2023, A&A, submitted [arXiv:2306.11053] [Google Scholar]
 Cataneo, M., Lombriser, L., Heymans, C., et al. 2019, MNRAS, 488, 2121 [Google Scholar]
 Cataneo, M., Uhlemann, C., Arnold, C., et al. 2022, MNRAS, 513, 1623 [CrossRef] [Google Scholar]
 Charmousis, C., Gregory, R., Kaloper, N., & Padilla, A. 2006, J. High Energy Phys., 10, 066 [CrossRef] [Google Scholar]
 Chen, S.F., Vlah, Z., & White, M. 2022, JCAP, 02, 008 [NASA ADS] [CrossRef] [Google Scholar]
 Cheng, G., Wu, F., & Chen, X. 2021, Phys. Rev. D, 103, 103527 [NASA ADS] [CrossRef] [Google Scholar]
 Clifton, T., Ferreira, P. G., Padilla, A., & Skordis, C. 2012, Phys. Rept., 513, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Cooray, A., & Sheth, R. K. 2002, Phys. Rev., 372, 1 [NASA ADS] [Google Scholar]
 D’Amico, G., Gleyzes, J., Kokron, N., et al. 2020, JCAP, 05, 005 [Google Scholar]
 Deffayet, C. 2001, Phys. Lett. B, 502, 199 [NASA ADS] [CrossRef] [Google Scholar]
 Dvali, G. R., Gabadadze, G., & Porrati, M. 2000, Phys. Lett. B, 485, 208 [NASA ADS] [CrossRef] [Google Scholar]
 Eisenstein, D. J., Seo, H.j., & White, M. J. 2007, ApJ, 664, 660 [NASA ADS] [CrossRef] [Google Scholar]
 Euclid Collaboration (Blanchard, A., et al.) 2020, A&A, 642, A191 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Euclid Collaboration (Ilic, S., et al.) 2022, A&A, 657, A91 [CrossRef] [EDP Sciences] [Google Scholar]
 Ferreira, P. G. 2019, ARA&A, 57, 335 [Google Scholar]
 Gleyzes, J., Langlois, D., & Vernizzi, F. 2015, Int. J. Mod. Phys. D, 23, 1443010 [Google Scholar]
 Gorbunov, D., Koyama, K., & Sibiryakov, S. 2006, Phys. Rev. D, 73, 044016 [NASA ADS] [CrossRef] [Google Scholar]
 HernándezAguayo, C., Ruan, C.Z., Li, B., et al. 2022, JCAP, 01, 048 [CrossRef] [Google Scholar]
 Hojjati, A., Pogosian, L., & Zhao, G.B. 2011, JCAP, 1108, 005 [Google Scholar]
 Hu, B., Raveri, M., Frusciante, N., & Silvestri, A. 2014, Phys. Rev. D, 89, 103530 [Google Scholar]
 Hu, W., & Sawicki, I. 2007a, Phys. Rev. D, 76, 104043 [Google Scholar]
 Hu, W., & Sawicki, I. 2007b, Phys. Rev. D, 76, 064004 [CrossRef] [Google Scholar]
 Ivanov, M. M., Simonovic, M., & Zaldarriaga, M. 2020, JCAP, 05, 042 [CrossRef] [Google Scholar]
 Joudaki, S., Blake, C., Johnson, A., et al. 2018, MNRAS, 474, 4894 [NASA ADS] [CrossRef] [Google Scholar]
 Joudaki, S., Ferreira, P. G., Lima, N. A., & Winther, H. A. 2022, Phys. Rev. D, 105, 043522 [NASA ADS] [CrossRef] [Google Scholar]
 Khoury, J., & Weltman, A. 2004, Phys. Rev. Lett, 93, 171104 [NASA ADS] [CrossRef] [Google Scholar]
 Koyama, K., & Maartens, R. 2006, JCAP, 01, 016 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, S., Huff, E. M., Choi, A., et al. 2021, MNRAS, 509, 4982 [CrossRef] [Google Scholar]
 Lesgourgues, J. 2011, arXiv eprints [arXiv:1104.2932] [Google Scholar]
 Lewis, A., Challinor, A., & Lasenby, A. 2000, ApJ, 538, 473 [Google Scholar]
 Liu, R., Valogiannis, G., Battaglia, N., & Bean, R. 2021, Phys. Rev. D, 104, 103519 [NASA ADS] [CrossRef] [Google Scholar]
 Lombriser, L., Hu, W., Fang, W., & Seljak, U. 2009, Phys. Rev. D, 80, 063536 [CrossRef] [Google Scholar]
 Lombriser, L., Yoo, J., & Koyama, K. 2013, Phys. Rev. D, 87, 104019 [NASA ADS] [CrossRef] [Google Scholar]
 Luty, M. A., Porrati, M., & Rattazzi, R. 2003, J. High Energy Phys., 09, 029 [CrossRef] [Google Scholar]
 Ma, C.P., & Bertschinger, E. 1995, ApJ, 455, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Mead, A., Peacock, J., Heymans, C., Joudaki, S., & Heavens, A. 2015, MNRAS, 454, 1958 [Google Scholar]
 Mead, A., Heymans, C., Lombriser, L., et al. 2016, MNRAS, 459, 1468 [NASA ADS] [CrossRef] [Google Scholar]
 Mead, A., Brieden, S., Tröster, T., & Heymans, C. 2021, MNRAS, 502, 1401 [Google Scholar]
 Mueller, E.M., Percival, W., Linder, E., et al. 2018, MNRAS, 475, 2122 [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [Google Scholar]
 Nguyen, N.M., Huterer, D., & Wen, Y. 2023, Phys. Rev. Lett., 131, 111001 [CrossRef] [Google Scholar]
 Ooba, J., Ichiki, K., Chiba, T., & Sugiyama, N. 2017, Prog. Theor. Exp. Phys., 2017, 043E03 [CrossRef] [Google Scholar]
 Pace, F., Battye, R., Bellini, E., et al. 2021, JCAP, 06, 017 [CrossRef] [Google Scholar]
 Perotto, L., Lesgourgues, J., Hannestad, S., Tu, H., & Wong, Y. Y. Y. 2006, JCAP, 10, 013 [CrossRef] [Google Scholar]
 Piga, L., Marinucci, M., D’Amico, G., et al. 2023, JCAP, 04, 038 [Google Scholar]
 Raccanelli, A., Bertacca, D., Pietrobon, D., et al. 2013, MNRAS, 436, 89 [NASA ADS] [CrossRef] [Google Scholar]
 Raveri, M., Hu, B., Frusciante, N., & Silvestri, A. 2014, Phys. Rev. D, 90, 043513 [Google Scholar]
 Raveri, M., Pogosian, L., Martinelli, M., et al. 2021, arXiv eprints [arXiv:2107.12990] [Google Scholar]
 Rossi, M., Ballardini, M., Braglia, M., et al. 2019, Phys. Rev. D, 100, 103524 [CrossRef] [Google Scholar]
 Schmidt, F. 2009, Phys. Rev. D, 80, 123003 [NASA ADS] [CrossRef] [Google Scholar]
 Schneider, A., Refregier, A., Grandis, S., et al. 2020a, JCAP, 04, 020 [Google Scholar]
 Schneider, A., Stoira, N., Refregier, A., et al. 2020b, JCAP, 04, 019 [CrossRef] [Google Scholar]
 Seahra, S. S., & Hu, W. 2010, Phys. Rev. D, 82, 124015 [NASA ADS] [CrossRef] [Google Scholar]
 Sheth, R. K., & Tormen, G. 1999, MNRAS, 308, 119 [Google Scholar]
 Sheth, R. K., & Tormen, G. 2002, MNRAS, 329, 61 [Google Scholar]
 Silvestri, A., Pogosian, L., & Buniy, R. V. 2013, Phys. Rev. D, 87, 104015 [NASA ADS] [CrossRef] [Google Scholar]
 Sola, J., GomezValent, A., de Cruz Perez, J., et al. 2019, Astrophys. J. Lett. 1, L6 [Google Scholar]
 Sola, J., GomezValent, A., de Cruz Perez, J., et al. 2020, Class. Quant. Grav. 24, 245003 [Google Scholar]
 Song, Y.S. 2008, Phys. Rev. D, 77, 124031 [NASA ADS] [CrossRef] [Google Scholar]
 Spurio Mancini, A., & Bose, B. 2023, arXiv eprints [arXiv:2305.06350] [Google Scholar]
 Spurio Mancini, A., Köhlinger, F., Joachimi, B., et al. 2019, MNRAS, 490, 2155 [Google Scholar]
 Takahashi, R., Sato, M., Nishimichi, T., Taruya, A., & Oguri, M. 2012, ApJ, 761, 152 [Google Scholar]
 Taruya, A., Nishimichi, T., & Saito, S. 2010, Phys. Rev. D, 82, 063522 [NASA ADS] [CrossRef] [Google Scholar]
 Tassev, S., Zaldarriaga, M., & Eisenstein, D. 2013, JCAP, 06, 036 [CrossRef] [Google Scholar]
 Teyssier, R. 2002, A&A, 385, 337 [CrossRef] [EDP Sciences] [Google Scholar]
 Touboul, P., Métris, G., Rodrigues, M., et al. 2017, Phys. Rev. Lett., 119, 231101 [NASA ADS] [CrossRef] [Google Scholar]
 Umiltà, C., Ballardini, M., Finelli, F., & Paoletti, D. 2015, JCAP, 08, 017 [CrossRef] [Google Scholar]
 Vainshtein, A. I. 1972, Phys. Lett. B, 39, 393 [NASA ADS] [CrossRef] [Google Scholar]
 Valageas, P., Nishimichi, T., & Taruya, A. 2013, Phys. Rev. D, 87, 083522 [NASA ADS] [CrossRef] [Google Scholar]
 Voisin, G., Cognard, I., Freire, P. C. C., et al. 2020, A&A, 638, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Williams, J. G., Turyshev, S. G., & Boggs, D. 2012, Class. Quant. Grav., 29, 184004 [NASA ADS] [CrossRef] [Google Scholar]
 Winther, H. A., Koyama, K., Manera, M., Wright, B. S., & Zhao, G.B. 2017, JCAP, 08, 006 [CrossRef] [Google Scholar]
 Wolz, L., Kilbinger, M., Weller, J., & Giannantonio, T. 2012, JCAP, 09, 009 [CrossRef] [Google Scholar]
 Xu, L. 2014, JCAP, 02, 048 [CrossRef] [Google Scholar]
 Zhao, G.B., Pogosian, L., Silvestri, A., & Zylberberg, J. 2009, Phys. Rev. D, 79, 083513 [CrossRef] [Google Scholar]
 Zucca, A., Pogosian, L., Silvestri, A., & Zhao, G.B. 2019, JCAP, 05, 001 [Google Scholar]
 Zumalacárregui, M., Bellini, E., Sawicki, I., et al. 2017, JCAP, 08, 019 [CrossRef] [Google Scholar]
All Tables
Forecast 1σ marginal errors on the cosmological and model parameters relative to their corresponding fiducial value for JBD1 and JBD2 in the pessimistic, quasilinear and optimistic cases, using Euclid observations GC_{sp}, WL+XC+GC_{ph} and GC_{sp}+WL+XC+GC_{ph}.
Forecast 1σ marginal errors on the cosmological and model parameters relative to their corresponding fiducial value for nDGPl and nDGP2 in the pessimistic, quasilinear and optimistic cases, using Euclid observations of GC_{sp}, WL+XC+GC_{ph} and GC_{sp}+WL+XC+GC_{ph}.
Forecast 1σ marginal errors on the cosmological and model parameters, relative to their corresponding fiducial value, for KM1 and KM2 in the pessimistic, quasilinear and optimistic cases, using Euclid observations of GC_{sp}, WL+XC+GC_{ph} and GC_{sp}+WL+XC+GC_{ph}.
All Figures
Fig. 1 Largescale 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 nonlinear (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 nonlinear (solid) cases. Bottom left: cosmic shear (WL) angular power spectra for the autocorrelation of the first photometric bin ${C}_{1,1}^{\text{WL}}(\ell )$ (solid line) and the crosscorrelation 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 C_{i,j}(ℓ) to their ΛCDM counterpart, for the bin combination i = j = 1 and i = 1; j = 10. 

In the text 
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 GC_{sp} in the optimistic (dotted area), pessimistic (dashed area) and quasilinear (filled area) scenarios (top panel), for WL+XC+GC_{ph} in the optimistic and pessimistic cases (central panel) and for the GC_{Sp}+WL+XC+GC_{ph} in the optimistic and pessimistic cases. 

In the text 
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 GC_{sp}, in blue, we plot WL+XC+GC_{ph} and in green we plot GC_{sp}+WL+XC+GC_{ph}. 

In the text 
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 GC_{sp}, in blue, we plot WL+XC+GC_{ph} and in green we plot GC_{sp}+WL+XC+GC_{ph}. 

In the text 
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 GC_{sp}, in blue, we plot WL+XC+GC_{ph} and in green, we plot GC_{sp}+WL+XC+GC_{ph}. 

In the text 
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 GC_{sp}, in blue, we plot WL+XC+GC_{ph} and in green, we plot GC_{sp}+WL+XC+GC_{ph}. 

In the text 
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 GC_{sp}, in blue, we plot WL+XC+GC_{ph} and in green, we plot GC_{sp}+WL+XC+GC_{ph}. 

In the text 
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 GC_{sp}, in blue, we plot WL+XC+GC_{ph} and in green, we plot GC_{sp}+WL+XC+GC_{ph}. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.