11institutetext: 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 INFN-Sezione di Torino, Via P. Giuria 1, 10125 Torino, Italy
4 INAF-Osservatorio Astrofisico di Torino, Via Osservatorio 20, 10025 Pino
Torinese (TO), Italy
5 INAF-Osservatorio Astronomico di Roma, Via Frascati 33, 00078 Monteporzio
Catone, Italy
6 INFN-Sezione 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 Ernest-Ansermet, CH-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á 2-4,
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 INAF-Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via
Piero Gobetti 93/3, 40129 Bologna, Italy
15 INFN-Sezione 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é Paris-Saclay 91191
Gif-sur-Yvette Cedex, France
23 Dipartimento di Fisica e Astronomia ”G.Galilei”, Università di Padova, Via
Marzolo 8, 35131 Padova, Italy
24 INFN-Padova, Via Marzolo 8, 35131 Padova, Italy
25 INAF-Osservatorio 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 INFN-Bologna, 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, 1749-016 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, P.O. Box 1029
Blindern, 0315 Oslo, Norway
38 Université Paris-Saclay, 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 INFN-Sezione
di Genova, via Dodecaneso 33, 16146, Genova, Italy
41 INFN-Sezione di Roma Tre, Via della Vasca Navale 84, 00146, Roma, Italy
42 INAF-Osservatorio 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, PT4150-762 Porto, Portugal
44 INAF-IASF 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 INAF-Osservatorio 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 Michel-Ange, 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, PT1749-016 Lisboa, Portugal
57 Department of Astronomy, University of Geneva, ch. d’Ecogia 16, 1290
Versoix, Switzerland
58 Université Paris-Saclay, Université Paris Cité, CEA, CNRS, Astrophysique,
Instrumentation et Modélisation Paris-Saclay, 91191 Gif-sur-Yvette, France
59 INAF-Osservatorio 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äts-Sternwarte München, Fakultät für Physik, Ludwig-Maximilians-
Universitä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 INFN-Sezione di Milano, Via Celoria 16, 20133 Milano, Italy
65 INAF-Osservatorio 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 Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg,
Germany
70 Aix-Marseille 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 Argelander-Institut 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, DK-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 Aix-Marseille 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, 1349-018 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

#  Euclid: Constraining linearly scale-independent modifications of gravity
with the spectroscopic and photometric primary probes††thanks: This paper is
published on behalf of the Euclid Consortium.

N. Frusciante1<EMAIL_ADDRESS>F. Pace2,3,4 V.F. Cardone5,6 S.
Casas7 I. Tutusaus8,9,10,11 M. Ballardini12,13,14,15 E. Bellini 16,17,18,8 G.
Benevento19 B. Bose20,8,21 P. Valageas22 N. Bartolo23,24,25 P. Brax22,26 P. G.
Ferreira27 F. Finelli14,28 K. Koyama29 L. Legrand8 L. Lombriser8 D.
Paoletti14,15 M. Pietroni30,31 A. Rozas-Fernández32 Z. Sakr33,34,35 A.
Silvestri36 F. Vernizzi22 H.A. Winther37 N. Aghanim38 L. Amendola34 N.
Auricchio14 R. Azzollini39 M. Baldi13,14,15 D. Bonino4 E. Branchini40,41 M.
Brescia1,42 J. Brinchmann43 S. Camera2,3,4 V. Capobianco4 C. Carbone44 J.
Carretero45,46 M. Castellano5 S. Cavuoti42,47 A. Cimatti48,49 R.
Cledassou50,51 G. Congedo20 L. Conversi52,53 Y. Copin54 L. Corcione4 F.
Courbin55 M. Cropper39 A. Da Silva56,32 H. Degaudenzi57 J. Dinis32,56 F.
Dubath57 X. Dupac53 S. Dusini24 S. Farrens58 S. Ferriol54 P. Fosalba10,9 M.
Frailis59 E. Franceschi14 S. Galeotta59 B. Gillis20 C. Giocoli14,60 A.
Grazian25 F. Grupp61,62 L. Guzzo63,64,65 S.V.H. Haugan37 W. Holmes66 F.
Hormuth67 A. Hornstrup68 K. Jahnke69 S. Kermiche70 A. Kiessling66 M.
Kilbinger58 T. Kitching39 M. Kunz8 H. Kurki-Suonio71 S. Ligori4 P. B. Lilje37
I. Lloro72 E. Maiorano14 O. Mansutti59 O. Marggraf73 K. Markovic66 F.
Marulli74,14,15 R. Massey75 E. Medinaceli14 M. Meneghetti15,14 G. Meylan55 M.
Moresco74,14 L. Moscardini74,14,15 E. Munari59 S.M. Niemi76 J. Nightingale75
C. Padilla45 S. Paltani57 F. Pasian59 K. Pedersen77 W.J. Percival78,79,80 V.
Pettorino58 G. Polenta81 M. Poncet50 L. Popa82 F. Raison61 R. Rebolo83,84 A.
Renzi23,24 J. Rhodes66 G. Riccio42 E. Romelli59 R. Saglia61,62 D. Sapone85 B.
Sartoris62,59 A. Secroun70 G. Seidel69 C. Sirignano23,24 G. Sirri15 L.
Stanco24 C. Surace86 P. Tallada-Crespí87,46 A.N. Taylor20 I. Tereno56,88 R.
Toledo-Moreo89 F. Torradeflot87,46 E.A. Valentijn90 L. Valenziano14,28 T.
Vassallo59 G.A. Verdoes Kleijn90 Y. Wang91 A. Zacchei59 G. Zamorani14 J.
Zoubian70 V. Scottez92

###### Abstract

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 scale-independent growth of
perturbations while featuring different testable types of derivative screening
mechanisms at smaller nonlinear scales.

Aims. We consider three specific models, namely Jordan-Brans-Dicke (JBD), a
scalar-tensor theory with a flat potential, the normal branch of Dvali-
Gabadadze-Porrati (nDGP) gravity, a braneworld model in which our Universe is
a four-dimensional brane embedded in a five-dimensional Minkowski space-time,
and $k$-mouflage (KM) gravity, an extension of $k$-essence 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 extra parameters of the
models, respectively, $\omega_{\rm BD}$, $\Omega_{\rm rc}$ and
$\epsilon_{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 forecasts analysis employs the Fisher matrix method applied to
weak lensing (WL); photometric galaxy clustering (GCph); spectroscopic galaxy
clustering (GCsp) and the cross-correlation (XC) between GCph and WL. For the
Euclid survey specifications we define three scenarios, characterised by
different cuts in the maximum multipole and wavenumber, to assess the
constraining power of nonlinear scales. For each model we consider two
fiducial values for the corresponding model parameter.

Results. In an optimistic setting at 68.3% confidence interval, with Euclid
alone we find the following percentage relative errors: for
$\log_{10}{\omega_{\rm BD}}$, with a fiducial value of $\omega_{\rm BD}=800$,
35% using GCsp alone, 3.6% using GCph+WL+XC and 3.3% using GCph+WL+XC+GCsp;
for $\log_{10}{\Omega_{\rm rc}}$, with a fiducial value of $\Omega_{\rm
rc}=0.25$, we find respectively 90%, 20% and 17%; finally, for
$\epsilon_{2,0}=-0.04$ respectively 5%, 0.15% and 0.14%. From the relative
errors for fiducial values closer to their $\Lambda\mathrm{CDM}$ limits, we
find that most of the constraining power is lost. Our results highlight the
importance of the constraining power coming from nonlinear scales.

###### Key Words.:

Gravitational lensing: weak – large-scale structure of Universe – cosmological
parameters

##  1 Introduction

We have entered a new era in gravitational physics where 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 one must add cosmological constraints using the
large scale structure of the Universe (Ferreira, 2019).

There have been multiple attempts at constraining GR with cosmological surveys
(Aghanim et al., 2020; Mueller et al., 2018; Joudaki et al., 2018; Lee et al.,
2021; Raveri et al., 2021). Using 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 pin down the growth rate of structure and its
morphology out to redshift $z\sim 2$. The resulting constraints have been an
important first step in understanding gravity in an altogether untested regime
yet they have been underwhelming. For example, constraints on the simplest
scalar-tensor modification to GR – Jordan-Brans-Dicke (JBD) gravity – are more
than order of magnitude weaker (Joudaki et al., 2022; Ballardini et al., 2022)
than those obtained from millisecond pulsar observations (Voisin et al.,
2020), albeit on very different scales.

The relevance of cosmological measurements for gravitational physics is about
to change with the upcoming generation of large-scale structure surveys
(Ferreira, 2019). By mapping out vast swathes of the Universe with exquisite
precision, one hopes that it will be possible to tighten substantially
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 and which can
also be used to reconstruct galaxy clustering using photometric redshifts) –
with a primary goal of pinning down the nature of dark energy and ideally
suited for cosmological constraints on gravity (Euclid Collaboration:
Blanchard et al., 2020, EC19 hereafter).

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

In this paper we will forecast how well the surveys from the Euclid mission
can be used to constrain a family of theories that modify the theory of GR yet
retain one of its properties: a scale-independent linear growth rate. The
three theories we will consider are: JBD gravity (Brans & Dicke, 1961), the
simplest scalar-tensor theory which involves a non-minimal coupling between a
scalar field and the metric; Dvali-Gabadadze-Porrati gravity (DGP, see Dvali
et al., 2000), a braneworld model which introduces modifications on
cosmological scales and screens with the Vainshtein mechanism; a $k$-mouflage
model (KM, see Babichev et al., 2009; Brax & Valageas, 2014, 2016), i.e. a
scalar-tensor theory with a non-canonical kinetic energy and the $k$-mouflage
screening property. The JBD theory is not screened and, therefore, suffers
from large deviations from GR unless the Brans-Dicke (BD) parameter
$\omega_{\rm BD}>4\times 10^{4}$ (Bertotti et al., 2003). This is larger than
the range of parameters testable cosmologically as we shall see below implying
that the JBD model investigated here must be taken as a template for large
scale deviations against which we compare the DGP and $k$-mouflage models. All
these modifications of GR will, in one way or another, affect the expansion
rate and growth of structure and are, thus, prime candidates to be constrained
by data from the Euclid mission.

We structure this paper as follows. In Sect. 2 we recap the essential facts
about linear cosmological perturbations in the context of extensions to GR,
identifying the phenomenological time-dependent parameters, ($\mu$, $\eta$ and
$\Sigma$), that are fed into the linear evolution equations. We then delve
into the three candidate theories we are exploring, laying out their
corresponding actions, background evolution, functional forms of {$\mu$,
$\eta$ and $\Sigma$} and how all this is implemented numerically. In Sect. 3
we present a detailed explanation of how we calculate all aspects of the
theoretical model predictions that will go into the forecasting procedure. In
Sect. 4 we describe the survey specifications and how these are integrated in
the analysis method; in this paper we will be using a Fisher forecasting
approach. In Sect. 5 we present the results of our methods and in Sect. 6 we
conclude.

##  2 Linearly scale-independent modified gravity

We follow the Bardeen formalism (Bardeen, 1980; Ma & Bertschinger, 1995), and
define the infinitesimal line element of the flat Friedmann-Lemaître-
Robertson-Walker (FLRW) metric by

$\mathrm{d}s^{2}=-(1+2\Psi)\,c^{2}\,\mathrm{d}t^{2}+a^{2}(t)\,(1-2\Phi)\,\delta_{ij}\,\mathrm{d}x^{i}\,\mathrm{d}x^{j}\,,$
(1)

where $a(t)$ is the scale factor as a function of the cosmic time $t$, $\Psi$
and $\Phi$ are the two scalar potentials and $c$ the speed of light. We will
work in Fourier space, where modes are functions of $t$ and the comoving
wavenumber $k^{i}$. At linear order, the matter energy-momentum tensor can be
decomposed as follows,

$T^{0}_{0}=-\bar{\rho}(1+\delta)c^{2}\,\,,\,T^{0}_{i}=(\bar{\rho}+\bar{p}/c^{2})cv_{i}\,\,,\,T^{i}_{j}=(\bar{p}+\delta
p)\delta^{i}_{j}+\Sigma^{i}_{j}\,,$ (2)

where $\delta=\rho/\bar{\rho}-1$ is the energy density contrast with $\rho$
being the matter density and $\bar{\rho}$ its background value,
$p=\bar{p}+\delta p$ is the pressure with $\bar{p}$ the background value,
$v_{i}$ is the peculiar velocity and $\Sigma^{i}_{j}$ 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, respectively defined as
$(\bar{\rho}c^{2}+\bar{p})\left(\hat{k}^{i}\cdot\hat{k}_{j}-1/3\delta^{i}_{j}\right)\sigma=\Sigma^{i}_{j}$
with $\hat{k}_{i}=k_{i}/k$ and
$\bar{\rho}\Delta\equiv\bar{\rho}\delta+3(aH/k/c^{2})(\bar{\rho}+\bar{p}/c^{2})v$,
where $v$ is the velocity potential defined through $v_{i}=-\imath k_{i}v/k$.
Here $H\equiv\dot{a}/a$ is the Hubble function and a dot stands for the
derivative with respect to the cosmic time. In analogy to the
$\Lambda\mathrm{CDM}$ model where curvature is found to be compatible with
zero, in this work we will assume a flat spatial geometry.

In models where gravity is modified by the presence of 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:

$\displaystyle-k^{2}\Psi$ $\displaystyle=\frac{4\pi\,G_{\rm
N}}{c^{2}}\,a^{2}\mu(a,k)\left[\bar{\rho}\Delta+3\left(\bar{\rho}+\bar{p}/c^{2}\right)\sigma\right]\,,$
(3) $\displaystyle-k^{2}\left(\Phi+\Psi\right)$
$\displaystyle=\frac{8\pi\,G_{\rm
N}}{c^{2}}\,a^{2}\Big{\\{}\Sigma(a,k)\left[\bar{\rho}\Delta+3\left(\bar{\rho}+\bar{p}/c^{2}\right)\sigma\right]$
$\displaystyle-\frac{3}{2}\mu(a,k)\left(\bar{\rho}+\bar{p}/c^{2}\right)\sigma\Big{\\}}\,,$
(4)

where $G_{\rm N}$ is Newton’s gravitational constant and $c$ the speed of
light. As such, a further function that can be introduced is the one defining
the ratio of the two potentials, $\eta=\Phi/\Psi$. In the absence of
anisotropic stress the three phenomenological functions are related by the
following expression:

$\Sigma=\frac{1}{2}\mu(1+\eta)\,.$ (5)

The _phenomenological functions_ $\mu$, $\eta$ and $\Sigma$ are identically
equal to $1$ in the GR limit. In general they are functions of time and scale.
The models that we consider in this paper preserve the scale-independent
(linear) growth pattern, i.e. they have $\mu=\mu(a)$, $\Sigma=\Sigma(a)$ and
$\eta=\eta(a)$.

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

In the following, we will introduce the three models under consideration,
i.e., JBD, DGP and KM. We will provide the background evolution equations and
the expressions for $\mu$, $\eta$ and $\Sigma$ for each of these models.

###  2.1 Jordan-Brans-Dicke gravity

The Jordan-Brans-Dicke theory of gravity (Brans & Dicke, 1961) is described by
the following action,

$S_{\rm
BD}=\int\mathrm{d}^{4}x\sqrt{-g}\left[\frac{c^{4}}{16\pi}\left({\phi}R-\frac{\omega_{\rm
BD}}{\phi}g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi-2\Lambda\right)+\mathcal{L}_{\rm
m}\right]\,,$ (6)

where $g_{\mu\nu}$ is the metric and $R$ its associated Ricci scalar, $\phi$
is the JBD scalar field (which has the dimensions of the inverse of Newton’s
constant), $\omega_{\rm BD}$ is the dimensionless BD parameter, and
$\mathcal{L}_{\rm m}$ is the matter Lagrangian which is minimally coupled to
the metric. We have included a cosmological constant, $\Lambda$. The model has
only one free parameter, $\omega_{\rm BD}$; in the limit in which $\omega_{\rm
BD}\rightarrow\infty$, 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. Yet, cosmological constraints on JBD gravity can have wider
implications if one takes the view that it is the long-wavelength, low-energy
limit of more general scalar-tensor theories (see Joudaki et al., 2022, for
example). Furthermore, more general scalar theories may be endowed with
gravitational screening which will kick in on smaller scales or,
alternatively, regions of high density, for example. As a consequence local
constraints will, to some extent, decouple from more global, or large scale,
constraints. Thus, constraints on cosmological scales on JBD theory may cover
a broad class of scalar-tensor theories and, furthermore, be independent from
model specific constraints on small scales. Thus, while simple, JBD gravity is
a powerful tool for constraining general classes of scalar-tensor theories
using cosmological data.

The modified Einstein field equations are

$\displaystyle G_{\mu\nu}={\frac{8\pi}{c^{4}}}\frac{1}{\phi}T_{\mu\nu}+$
$\displaystyle\,\frac{\omega_{\rm{BD}}}{\phi^{2}}\left(\nabla_{\mu}\phi\nabla_{\nu}\phi-\frac{1}{2}g_{\mu\nu}\nabla_{\alpha}\phi\nabla^{\alpha}\phi\right)$
$\displaystyle+$
$\displaystyle\,\frac{1}{\phi}\left[\nabla_{\mu}\nabla_{\nu}\phi-
g_{\mu\nu}\left(\Box\phi+\Lambda\right)\right]\,.$ (7)

Here $T_{\mu\nu}$ is the total matter stress-energy tensor while $\Box$
denotes the d’Alembertian. On the background, these equations give

$\displaystyle 3H^{2}=$
$\displaystyle\,{\frac{8\pi}{c^{4}}}\frac{\bar{\rho}}{\phi}-3H\frac{\dot{\phi}}{\phi}+\frac{\omega_{\rm{BD}}}{2}\frac{\dot{\phi}^{2}}{\phi^{2}}+\frac{\Lambda}{\phi}\,,$
(8) $\displaystyle 2\dot{H}+3H^{2}=$
$\displaystyle\,-{\frac{8\pi}{c^{4}}}\frac{\bar{p}}{\phi}-\frac{\omega_{\rm{BD}}}{2}\frac{\dot{\phi}^{2}}{\phi^{2}}-2H\frac{\dot{\phi}}{\phi}-\frac{\ddot{\phi}}{\phi}+\frac{\Lambda}{\phi}\,,$
(9)

and the scalar field equation of motion reads

$\Box\phi={\frac{8\pi}{c^{4}}}\left({\frac{T}{3+2\omega_{\rm{BD}}}}\right)-\frac{4\Lambda}{3+2\omega_{\rm{BD}}}\,,$
(10)

where $T\equiv g^{\mu\nu}T_{\mu\nu}$ is the trace of the stress-energy tensor
on the background.

The phenomenological QS functions in this theory read

$\displaystyle\Sigma=$ $\displaystyle\,\frac{1}{G_{\rm N}\phi}\,,$
$\displaystyle\mu=$ $\displaystyle\,\frac{4+2\omega_{\rm BD}}{3+2\omega_{\rm
BD}}\Sigma\,,$ (11) $\displaystyle\eta\equiv$
$\displaystyle\,\frac{\Phi}{\Psi}=\,\frac{1+\omega_{\rm BD}}{2+\omega_{\rm
BD}}\,.$

Constraints on JBD gravity have been obtained by using a combination of
different cosmological data sets and sampling over different parameterisation
of the BD parameter $\omega_{\rm BD}$. For example, Avilez & Skordis (2014)
imposed a flat prior on $-\log_{10}\left(\omega_{\rm BD}\right)$ to obtain a
lower bound of $\omega_{\rm BD}>1900$ at 95% CL with CMB information from
Planck 2013 data. Ballardini et al. (2016) obtained
$\log_{10}\left(1+1/\omega_{\rm BD}\right)<0.0030$ at 95% CL combining Planck
2015 and BOSS DR10-11 data; this upper bound was subsequently updated in
Ballardini et al. (2020) with a combination of Planck 2018 and BOSS DR12 data
to $\log_{10}\left(1+1/\omega_{\rm BD}\right)<0.0022$ at 95% CL. Joudaki et
al. (2022) used a combination of the Planck 2018 CMB data, the $3\times 2$pt
combination of the KiDS and 2dFLens data, the Pantheon supernovae data, BOSS
measurements of the BAO, to find the coupling constant $\omega_{\rm
BD}>\\{1540,\,160,\,160,\,350\\}$ at 95% CL for the different choices of
parametrization (or priors): $\left\\{\log_{10}\left(1/\omega_{\rm
BD}\right),\,\log_{10}\left(1+1/\omega_{\rm BD}\right),\,1/\omega_{\rm
BD},\,1/\log_{10}\omega_{\rm BD}\right\\}$. These constraints were obtained
fixing the value of the scalar field today to

$G_{\rm N}\phi\left(a=1\right)\,=\,\frac{2\omega_{\rm BD}+4}{2\omega_{\rm
BD}+3}\,,$ (12)

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

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

###  2.2 Dvali-Gabadadze-Porrati braneworld gravity

The Dvali-Gabadadze-Porrati model (Dvali et al., 2000) is a five-dimensional
braneworld model defined by the action

$S=\frac{c^{4}}{16\pi G_{5}}\int_{\cal
M}\mathrm{d}^{5}x\sqrt{-\gamma}R_{5}+\int_{\partial{\cal
M}}\mathrm{d}^{4}x\sqrt{-g}\left[\frac{c^{4}}{16\pi G_{\rm N}}R+{\cal L}_{\rm
m}\right]\,,$ (13)

where $\gamma$ is the five-dimensional metric and $R_{5}$ its Ricci curvature
scalar. $G_{5}$ and $G_{\rm N}$ are the five- and four-dimensional Newton’s
constants, respectively. The matter Lagrangian is denoted with ${\cal L}_{\rm
m}$ and is confined to a four-dimensional brane in a five-dimensional
Minkowski spacetime. The induced gravity described by the usual four-
dimensional Einstein-Hilbert action is responsible for the recovery of the
four-dimensional gravity on the brane. The cross-over scale $r_{\rm
c}=G_{5}/(2G_{\rm N})$ is the only parameter of the model and GR is recovered
in the limit $r_{\rm c}\to\infty$.

The Friedmann equation on the brane is given by (Deffayet, 2001)

$H^{2}=\pm\,c\frac{H}{r_{\rm c}}+\frac{8\pi G}{3}\bar{\rho}\,.$ (14)

There are two branches of solutions depending on the embedding of the brane:
the self-accelerating branch (sDGP, Bowcock et al., 2000; Deffayet, 2001) and
the normal branch (nDGP, Bowcock et al., 2000; Deffayet, 2001), corresponding
to the $+$ and $-$ sign for the contribution from the five-dimensional
gravity, respectively. The self-accelerating branch admits the late-time
acceleration without dark energy but the solution is plagued by ghost
instabilities (Luty et al., 2003; Gorbunov et al., 2006; Charmousis et al.,
2006). Thus we 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
$\Lambda\mathrm{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 will adopt 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 $\mu$, $\eta$, $\Sigma$. Using the results of Koyama & Maartens
(2006); Hu & Sawicki (2007a); Song (2008); Cardoso et al. (2008); Lombriser et
al. (2009); Seahra & Hu (2010), we have

$\eta=\frac{1+g}{1-g}\,,$ (15)

where, using the quasi-static (QS) approximation, we have (Lombriser et al.,
2013)

$g(a)=g_{\rm
QS}=-\frac{1}{3}\left[1\mp\frac{2Hr_{c}}{c}\left(1+\frac{\dot{H}}{3H^{2}}\right)\right]^{-1}\,,$
(16)

such that the effective modification introduced with $\eta$ can be treated as
scale-independent.

The effective $3+1$ Poisson equation for the lensing potential in the QS
approximation is

$-k^{2}(\Phi+\Psi)=\frac{8\pi G_{\rm N}}{c^{2}}\,a^{2}\,\bar{\rho}\,\Delta\,,$
(17)

therefore (Lombriser et al., 2013)

$\Sigma=1\,,$ (18)

and hence

$\mu(a)=1+\frac{1}{3\beta}\,,\qquad\beta(a)\equiv
1+\frac{H}{H_{0}}\frac{1}{\sqrt{\Omega_{\rm
rc}}}\left(1+\frac{\dot{H}}{3H^{2}}\right)\,.$ (19)

We have chosen to parameterise the modification to gravity by $\Omega_{\rm
rc}\equiv c^{2}/(4r_{\rm c}^{2}H_{0}^{2})$.

There are not many studies on constraining nDGP with an exact
$\Lambda\mathrm{CDM}$ background. However, Raccanelli et al. (2013), using
measurements of the zeroth and second order moments of the correlation
function from SDSS DR7 data up to $r_{\rm max}=120\,{\rm Mpc}\,h^{-1}$, and
marginalized bias, found an $\Omega_{\rm rc}$ upper limit at 95% of $\sim 40$
(from $r_{\rm c}>340$ Mpc with fixed $H_{0}=70\,{\rm km}\,\rm{s}^{-1}\,{\rm
Mpc}^{-1}$). In addition, Barreira et al. (2016) used clustering wedges
statistic of the galaxy correlation function and the growth rate values
estimated from more recent BOSS DR12 data to set $[r_{\rm
c}H_{0}/c]^{-1}<0.97$ at 95% C.L., corresponding to an upper bound of
$\Omega_{\rm rc}<0.27$.

There have also been recent attempts to forecast $\Omega_{\rm rc}$. Liu et al.
(2021), using galaxy cluster abundance from the Simons Observatory and galaxy
correlation functions from a DESI-like experiment, found $\delta(\Omega_{\rm
rc})\sim$ 0.038 around a fiducial value of 0.25. Cataneo et al. (2022)
forecast for Euclid-like future constraints on a fiducial $\Omega_{\rm rc}\sim
0.0625$ a $1\sigma$ accuracy of 0.0125 from combining the 3D matter power
spectrum and the probability distribution function of the smoothed three-
dimensional matter density field probes. Also Bose et al. (2020) forecast for
a LSST-like survey, using cosmic shear alone, on a fiducial $\Omega_{\rm
rc}\sim 0$ a $1\sigma$ accuracy of 0.08.

Note that constraints have also been inferred for a nDGP model with
cosmological constant rather than a constructed dark energy field, thus with
an approximate $\Lambda\mathrm{CDM}$ background. For this model, Lombriser et
al. (2009) and Xu (2014) inferred somewhat stronger constraints of
$\Omega_{\rm rc}<0.020$ (95% C.L.) and $\Omega_{\rm rc}<0.002$ (68% C.L.) from
the combination of CMB and large-scale structure data, where the constraints
are mainly driven by the CMB data. Thus, high-precision CMB measurements such
as by the Planck satellite tightly constrain $\Omega_{\rm rc}$. In this work
we focus on the model with the exact $\Lambda\mathrm{CDM}$ background as our
interest is to see how Euclid can constrain the modification in structure
formation that is not coming from a different background.

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

###  2.3 $k$-mouflage gravity

$k$-mouflage theories are built complementing simple $k$-essence scenarios
with a universal coupling of the scalar field $\varphi$ to matter. They are
defined by the action (Babichev et al., 2009; Brax & Valageas, 2014, 2016)

$\displaystyle S$
$\displaystyle=\int\mathrm{d}^{4}x\sqrt{-\tilde{g}}\left[\frac{c^{4}}{16\pi\tilde{G}_{\rm
N}}\tilde{R}+c^{2}{\cal M}^{4}K(\tilde{\chi})\right]+S_{\rm
m}(\psi_{i},g_{\mu\nu})\,,$ (20)

where $\tilde{G}_{\rm N}$ is the bare Newton constant, ${\cal M}^{4}$ is the
energy density scale of the scalar field, $g_{\mu\nu}$ is the Jordan frame
metric, $\tilde{g}_{\mu\nu}$ is the Einstein frame metric with
$g_{\mu\nu}=A^{2}(\varphi)\tilde{g}_{\mu\nu}$, $S_{\rm 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}\varphi\partial_{\nu}\varphi}{2{\cal
M}^{4}}\,,$ (21)

and ${\cal M}^{4}K$ is the non-standard kinetic term of the scalar field.

In these theories the evolution of both the background and perturbations is
affected by the universal coupling and by the scalar field dynamics. The
amount of deviation from $\Lambda\mathrm{CDM}$ at the background level and in
perturbation theory can be expressed in terms of two time-dependent functions,
related to the coupling $A$ and the kinetic function $K$, i.e.,

$\epsilon_{2}\equiv\frac{\mathrm{d}\ln\bar{A}}{\mathrm{d}\ln
a}\,,\quad\epsilon_{1}\equiv\frac{2}{\bar{K}^{\prime}}\left[\epsilon_{2}\frac{c}{\sqrt{8\pi\tilde{G}_{\rm
N}}}\left(\frac{\mathrm{d}\bar{\varphi}}{\mathrm{d}\ln
a}\right)^{-1}\right]^{2}\,,$ (22)

where over-bars denote background quantities and a prime indicates derivatives
w.r.t. $\bar{\tilde{\chi}}$. The KM Friedmann equation therefore reads

$H^{2}=\frac{8\pi\tilde{G}_{\rm
N}}{3}\frac{\bar{A}^{2}}{(1-\epsilon_{2})^{2}}\;\left[\bar{\rho}+\frac{{\cal
M}^{4}}{\bar{A}^{4}}\left(2\bar{\tilde{\chi}}\frac{\mathrm{d}\bar{K}}{\mathrm{d}\bar{\tilde{\chi}}}-\bar{K}\right)\right]\,.$
(23)

Considering linear scalar perturbations around an FLRW background, the
phenomenological functions $\mu$ and $\Sigma$ read (Benevento et al., 2019):

$\mu(a)=(1+\epsilon_{1})\bar{A}^{2}\,,\quad\Sigma(a)=\bar{A}^{2}\,.$ (24)

We recover GR when $\epsilon_{1}\to 0$, $\epsilon_{2}\to 0$, and $A$ and $K$
are constants, with ${\cal M}^{4}K$ playing the role of the cosmological
constant, see Eq. (20).

In addition to the six standard $\Lambda\mathrm{CDM}$ parameters, the KM model
requires the specification of the kinetic function $K(\tilde{\chi})$ and of
the coupling $A(\varphi)$. Following Brax & Valageas (2016), this can also be
expressed in terms of the time-dependent background values $\bar{K}(a)$ and
$\bar{A}(a)$. Brax & Valageas (2016) proposed a simple parameterization that
satisfies self-consistency constraints. It involves five parameters,
$\\{\epsilon_{2,0},\gamma_{A},m,\alpha_{U},\gamma_{U}\\}$, where
$\epsilon_{2,0}$ is the value of $\epsilon_{2}$ at redshift $z=0$ and must be
negative to ensure there are no ghosts in the theory. For the coupling
function $\bar{A}(a)$, this reads

$\bar{A}(a)=1+\alpha_{A}-\alpha_{A}\left[\frac{(\gamma_{A}+1)\,a}{\gamma_{A}+a}\right]^{\nu_{A}}\,,$
(25)

with

$\nu_{A}=\frac{3(m-1)}{2m-1}\,,\quad\alpha_{A}=-\frac{\epsilon_{2,0}\,(1+\gamma_{A})}{\gamma_{A}\nu_{A}}\,.$
(26)

Here, $m$ is the exponent of the kinetic function for large argument,
$K\sim\tilde{\chi}^{m}$, which also specifies the high-redshift dependence of
$\bar{A}(a)$, while $\gamma_{A}$ sets the redshift where $\bar{A}$ goes from
the current unit value to the high-redshift value $1+\alpha_{A}$, which is
parameterized by $\epsilon_{2,0}$. The background kinetic function
$\bar{K}(a)$ is conveniently parameterized in terms of a function $U(a)$, with

$\frac{\mathrm{d}\bar{K}^{\prime}}{\mathrm{d}\tilde{\chi}}=\frac{U(a)}{a^{3}\sqrt{\tilde{\chi}}}\,,\quad\sqrt{\tilde{\chi}}=-\frac{\bar{\rho}_{0}}{{\cal
M}^{4}}\frac{\epsilon_{2}\bar{A}^{4}}{2U\left(-3\epsilon_{2}+\frac{\mathrm{d}\ln
U}{\mathrm{d}\ln a}\right)}\,,$ (27)

which we take to be of the form

$U(a)\propto\frac{a^{2}\ln(\gamma_{U}+a)}{(\sqrt{a_{\rm
eq}}+\sqrt{a})\ln(\gamma_{U}+a)+\alpha_{U}a^{2}}\,.$ (28)

The various terms in this expression allow us to follow the radiation, matter
and dark-energy eras. The parameters $\alpha_{U}$ and $\gamma_{U}$, of order
unity, set the shape of the transition to the dark-energy era. The main
parameter is $\epsilon_{2,0}$, which measures the amplitude of the deviation
from GR and $\Lambda\mathrm{CDM}$ at $z=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, polarization 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\leq\epsilon_{2,0}\leq 0$, while the other parameters are unconstrained
(as long as they do not become very large). These results are consistent with
earlier more qualitative CMB studies (Barreira et al., 2015). The X-ray
cluster multiplicity function could provide bounds of the same order (Brax et
al., 2015), but more detailed analysis are needed to derive robust
constraints.

The fact that $\epsilon_{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 $\epsilon_{2}$, we have $\epsilon_{1}\simeq-\epsilon_{2}$. Indeed,
this shows that the running of Newton’s constant and the impact on the linear
metric and matter density fluctuations are set by $\epsilon_{2}$. Because most
of the deviations from $\Lambda\mathrm{CDM}$ occur at low redshift, this is
mostly determined by the value $\epsilon_{2,0}$ at $z=0$. Therefore, in this
work we focus on the dependence on $\epsilon_{2,0}$, which is left free to
vary, while we fix the other parameters to $\gamma_{A}=0.2$, $m=3$,
$\alpha_{U}=0.2$, and $\gamma_{U}=1$, as some representative values for the
other parameters. Thus, this gives one additional parameter, $\epsilon_{2,0}$,
in addition to the six standard $\Lambda\mathrm{CDM}$ parameters. The
$\Lambda\mathrm{CDM}$ model is recovered when $\epsilon_{2,0}=0$,
independently of the value of the other parameters; indeed for
$\epsilon_{2,0}\rightarrow 0$ we have:

$\bar{A}(a)\rightarrow 1\,,\quad\epsilon_{2}(a)\rightarrow
0\,,\quad\epsilon_{1}(a)\rightarrow 0\,,\quad\forall\ a\,,$ (29)

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

To produce linear predictions for this theory we use its implementation in the
EFTCAMB code, 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
the ones in EC19 apart from 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 the reader to Casas et al.
(2023), while here we only remind the reader of the main steps.

###  3.1 Photometric survey

The observables we consider 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 (GCph), or their cross
correlation, XC. These are still calculated as in EC19 relying on the Limber
approximation and setting to unity the $\ell$ \- dependent factor in the flat
sky limit. The spectra are then given by

$C^{XY}_{ij}(\ell)=c\int_{z_{\rm min}}^{z_{\rm
max}}\mathrm{d}z\,{\frac{W_{i}^{X}(z)W_{j}^{Y}(z)}{H(z)r^{2}(z)}P_{\delta\delta}(k_{\ell},z)}\,,$
(30)

with $k_{\ell}=(\ell+1/2)/r(z)$, $r(z)$ the comoving distance to redshift
$z=1/a-1$, $P_{\delta\delta}(k_{\ell},z)$ the nonlinear power spectrum of
matter density fluctuations, $\delta$, at wave number $k_{\ell}$ and redshift
$z$. We set $(z_{\rm min},z_{\rm max})=(0.001,4)$ which span the full range
where the source redshift distributions $n_{i}(z)$ are non vanishing. The GCph
and WL window functions read (Spurio Mancini et al., 2019)

$\displaystyle W_{i}^{\rm G}(k,z)=$
$\displaystyle\;b_{i}(k,z)\,\frac{n_{i}(z)}{\bar{n}_{i}}\frac{H(z)}{c}\,,$
(31) $\displaystyle W_{i}^{\rm L}(k,z)=$
$\displaystyle\;\frac{3}{2}\Omega_{\mathrm{m},0}\frac{H_{0}^{2}}{c^{2}}(1+z)\,r(z)\,\Sigma(z)\int_{z}^{z_{\rm
max}}{\mathrm{d}z^{\prime}\frac{n_{i}(z^{\prime})}{\bar{n}_{i}}\frac{r(z^{\prime}-z)}{r(z^{\prime})}}$
$\displaystyle+W^{\rm IA}_{i}(k,z)\,,$ (32)

where $n_{i}(z)/\bar{n}_{i}$ and $b_{i}(k,z)$ are, respectively, the
normalised galaxy distribution and the galaxy bias in the $i$-th redshift bin,
and $W^{\rm IA}_{i}(k,z)$ encodes the contribution of intrinsic alignments
(IA) to the WL power spectrum. Notice the presence of the function $\Sigma(z)$
in the WL window function which 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(z)$, $r(z)$ and the matter power spectrum
$P_{\delta\delta}(k,z)$ are, on the contrary, already taken into account by
the dedicated Boltzmann solver so that the GCph window function remains
unchanged.

The IA contribution is computed following the eNLA model adopted in EC19 so
that the corresponding window function is

$W^{\rm IA}_{i}(k,z)=-\frac{\mathcal{A}_{\rm IA}\,\mathcal{C}_{\rm
IA}\,\Omega_{\rm m,0}\,\mathcal{F}_{\rm
IA}(z)}{\delta(k,z)/\delta(k,z=0)}\frac{n_{i}(z)}{\bar{n}_{i}(z)}\frac{H(z)}{c}\,,$
(33)

where

$\mathcal{F}_{\rm IA}(z)=(1+z)^{\eta_{\rm IA}}\left[\frac{\langle
L\rangle(z)}{L_{\star}(z)}\right]^{\beta_{\rm IA}}\,,$ (34)

with $\langle L\rangle(z)$ and $L_{\star}(z)$ redshift-dependent mean and the
characteristic luminosity of source galaxies as computed from the luminosity
function, $\mathcal{A}_{\rm IA}$, $\beta_{\rm IA}$ and $\eta_{\rm IA}$ are the
nuisance parameters of the model, and $\mathcal{C}_{\rm IA}$ is a constant
accounting for dimensional units. This model is the same as the one used in
EC19 since IA takes place on astrophysical scales which 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 $\delta(k,z)$, which is considered
to be scale dependent in this case. This allows us to consider also the scale
dependence introduced by massive neutrinos, which was assumed to be negligible
in EC19. We nevertheless stress that, for the models we are considering in
this paper, 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 scalar fluctuations is of the order of the
speed of light, so that the sound horizon is of the order of the Hubble scale.
Under these assumptions, we can apply the quasi-static approximation and
relate the scalar field perturbation to the gravitational potential. Since in
all these models the weak equivalence principle holds, the modelling of the
bias as an expansion in the derivatives of the gravitational potential remains
unchanged with respect to the $\Lambda\mathrm{CDM}$ one. To be consistent with
the official forecast, we will use the same modelling for galaxy clustering as
in EC19.

The observed galaxy power spectrum is given by

$P_{\text{obs}}(k,\mu_{\theta};z)=\frac{1}{q_{\perp}^{2}(z)\,q_{\parallel}(z)}\left\\{\frac{\left[b\sigma_{8}(z)+f\sigma_{8}(z)\mu_{\theta}^{2}\right]^{2}}{1+\Big{[}f(z)k\mu_{\theta}\sigma_{\rm
p}(z)\Big{]}^{2}}\right\\}\\\
\times\frac{P_{\text{dw}}(k,\mu_{\theta};z)}{\sigma_{8}^{2}(z)}F_{z}(k,\mu_{\theta};z)+P_{\text{s}}(z)\,,$
(35)

where $P_{\rm dw}(k,\mu_{\theta};\,z)$, the de-wiggled power spectrum,
includes the correction that accounts for the smearing of the BAO features,

$P_{\text{dw}}(k,\mu_{\theta};z)=P_{\delta\delta}^{\rm
lin}(k;z)\,\text{e}^{-g_{\mu}k^{2}}+P_{\text{nw}}(k;z)\left(1-\text{e}^{-g_{\mu}k^{2}}\right)\,,$
(36)

and $P_{\delta\delta}^{\rm lin}(k;z)$ stands for the linear matter power
spectrum. $P_{\rm nw}(k;z)$ is a ‘no-wiggle’ power spectrum with the same
broad-band shape as $P_{\delta\delta}^{\rm lin}(k;z)$ but without BAO features
(see discussion below). The function

$g_{\mu}(k,\mu_{\theta},z)=\sigma_{\rm
v}^{2}(z)\left\\{1-\mu_{\theta}^{2}+\mu_{\theta}^{2}[1+f(z)]^{2}\right\\}\,,$
(37)

is the nonlinear damping factor of the BAO signal derived by Eisenstein et al.
(2007), with

$\sigma_{\rm v}^{2}(z)=\frac{1}{6\pi^{2}}\int\mathrm{d}kP_{\delta\delta}^{\rm
lin}(k,z)\,.$ (38)

The curly bracket in Eq. (35) is the redshift-space-distortion (RSD)
contribution correcting for the nonlinear finger-of-God (FoG) effect, where we
defined $b\sigma_{8}(z)$ as the product of the effective scale-independent
bias of galaxy samples and the r.m.s. matter density fluctuation
$\sigma_{8}(z)$ (let us note that we marginalize over $b\sigma_{8}(z)$), while
$\mu_{\theta}$ is the cosine of the angle $\theta$ between the wave vector
$\bm{k}$ and the line-of-sight direction $\hat{\bm{r}}$ and $\sigma_{\rm
p}^{2}(z)=\sigma_{\rm v}^{2}(z)$. Although these parameters are assumed to be
the same, they come from two different physical effects, namely large-scale
bulk flow for the former and virial motion for the latter.

The factor $F_{z}$ accounts for the smearing of the galaxy density field along
the line-of-sight due to redshift uncertainties. It is given by

$F_{z}(k,\mu_{\theta};z)=\text{e}^{-k^{2}\mu_{\theta}^{2}\sigma_{r}^{2}(z)}\,,$
(39)

where $\sigma_{r}(z)=c\,(1+z)\sigma_{0,z}/H(z)$ and $\sigma_{0,z}$ is the
error on the measured redshifts.

The factor in front of the curly bracket in Eq. (35) describes the Alcock-
Paczynski effect, which is parameterised in terms of the angular diameter
distance $D_{\rm A}(z)$ and the Hubble parameter $H(z)$ as

$\displaystyle q_{\perp}(z)$ $\displaystyle=\frac{D_{\rm A}(z)}{D_{\rm
A,\,ref}(z)},$ (40) $\displaystyle q_{\parallel}(z)$
$\displaystyle=\frac{H_{\text{ref}}(z)}{H(z)}\,.$ (41)

Due to the Alcock-Paczynski effect, $\mu_{\theta}$ and $k$ are re-scaled in a
cosmology-dependent way as a function of the projection along and
perpendicular to the line-of-sight. In the previous Eq. (35), the arguments
$\mu$ and $k$ are themselves functions of the “true” $\mu_{\theta,\rm ref}$
and $k_{\rm ref}$ at the reference cosmology. This relation, which for each
argument, is a function of $q_{\parallel}$ and $q_{\perp}$, can be found in
Sect. 3.2.1 of EC19. Finally, the $P_{\rm s}(z)$ is a scale-independent shot-
noise term, which enters as a nuisance parameter (see EC19, ).

Due to the scale independence of $\sigma_{v}(=\sigma_{p})$, we evaluate 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 in this work we take directly the derivatives of the observed
galaxy power spectrum with respect to the final parameters, contrary to EC19
where first a Fisher matrix analysis is performed for the redshift dependent
parameters $H(z)$, $D_{\rm A}(z)$ and $f\sigma_{8}(z)$ and then projected to
the final cosmological parameters of interest. However, we checked that both
approaches lead to consistent results when considering the $\Lambda$CDM and
$w_{0}w_{a}$CDM models. The other term appearing in Eq. (36) is the non-wiggle
matter power spectrum and this has been obtained applying a Savitzky-Golay
filter to the matter power spectrum $P_{\delta\delta}^{\rm lin}(k;z)$. For
more details on this implementation, see EC19.

###  3.3 Nonlinear modelling for WL

The galaxy power spectrum on mildly nonlinear scales can be modelled by a
modified version of the implementation in EC19. This is no more the case when
we move to the deeply nonlinear regime where an analytical description of the
matter power spectrum is unavailable for the models of our interest.

The wide window functions (in particular, the WL one) entering the prediction
of the photometric observables $C_{ij}^{XY}(\ell)$ require an accurate
description of nonlinear scales. In particular, one should explicitly account
for the impact of baryons, and include nuisance parameters which control the
baryonic feedback in order to avoid biasing the parameter estimation (see e.g.
Schneider et al., 2020b, a). However, at the moment we do not have accurate
Euclid-like simulations including baryonic effects, especially in the case of
MG cosmologies. Therefore, we ignore these effects in our analysis, leaving
their inclusion for a future work. In the following subsections, we comment on
the individual prescriptions for each MG model, which will be 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
behaviors 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 which have already been tested and used in literature thus
reinforcing the validity of our procedure.

####  3.3.1 Jordan-Brans-Dicke gravity

For the nonlinear prescription, we use a modified version of HMCODE (Mead et
al., 2015, 2016). HMCODE is an augmented halo model that can be used to make
accurate predictions of 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_{\rm
NL}(k,z)=P_{1\mathrm{H}}(k,z)+P_{2\mathrm{H}}(k,z)$ where

$P_{2\mathrm{H}}(k,z)=P_{\mathrm{lin}}(k,z)\left[\int_{0}^{\infty}b(M,z)W(M,k,z)\,n(M,z)\;\mathrm{d}M\right]^{2}\,,$
(42)

is the so-called 2-halo term (correlation between different halos) and

$P_{1\mathrm{H}}(k,z)=\int_{0}^{\infty}W^{2}(M,k,z)\,n(M,z)\;\mathrm{d}M\,,$
(43)

is the so-called 1-halo term (correlations between mass-elements within each
halo). Above $M$ is the halo mass, $P_{\mathrm{lin}}(k,z)$ is the linear
matter power spectrum, $n(M,z)$ is the halo mass function and $b(M,z)$ is the
linear halo bias.

The window function $W$ is the Fourier transform of the halo matter density
profiles:

$W(M,k,z)=\frac{1}{\bar{\rho}}\int_{0}^{r_{\mathrm{v}}}4\pi
r^{2}\frac{\sin(kr)}{kr}\,\rho(M,r,z)\;\mathrm{d}r\,,$ (44)

where $\rho(M,r,z)$ is the radial matter density profile in a host halo of
mass $M$ and $\bar{\rho}$ is the mean matter density. The halo mass is related
to the virial radius, $r_{\mathrm{v}}$, via $M=4\pi
r_{\mathrm{v}}^{3}\Delta_{\mathrm{v}}(z)\bar{\rho}/3$, where
$\Delta_{\mathrm{v}}(z)$ 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 which are then fit to data. HMCODE does a
combination of this and the resulting fitting function has, for
$\Lambda\mathrm{CDM}$, been shown to be accurate to $2.5\%$ for scales
$k<10\,h\,{\rm Mpc}^{-1}$ and redshifts $z<2$ (Mead et al., 2021).

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

The spectra generated by these simulations were then used to calibrate HMCODE.
While we do not consider the effects of baryons here, the advantage of HMCODE,
rather than Halofit is that the former is able to capture baryonic feedback at
$\sim 5\%$ accuracy level up to $k\simeq 10\,h\,{\rm Mpc}^{-1}$ (Mead et al.,
2021). This will become an important consideration (Schneider et al., 2020b,
a) in future, more detailed analyses. To take into account JBD, on top of
modifying the expansion history and growth of perturbations, the virialised
halo overdensity was modified as

$\displaystyle\Delta_{\rm v}=$ $\displaystyle\Omega_{\rm m}(z)^{-0.352}$ (45)
$\displaystyle\times\left\\{d_{0}+(418.0-d_{0})\arctan\left[(0.001|\omega_{\rm
BD}-50.0|)^{0.2}\right]\frac{2}{\pi}\right\\}\,,$

where $d_{0}=320.0+40.0z^{0.26}$. This modification has the feature that it
reduces to the original $\Lambda\mathrm{CDM}$ HMCODE as $\omega_{\rm
BD}\to\infty$. 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 HMCODE
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}\gtrsim\omega_{\rm BD}\geq 50$
and on scales up to $k\simeq 10\,h\,{\rm Mpc}^{-1}$. This nonlinear
prescription was made with past weak-lensing surveys in mind and fitted to a
range of $\omega_{\rm BD}$ and scales that we can realistically constrain at
the present time. Let us note that if doing parameter inference with actual
Euclid data the accuracy here might not be good enough.

We can see in Fig. 1 the comparison of the matter power spectrum and the
lensing angular power spectrum of JBD with $\omega_{\rm BD}=800$, which we
will refer to as JBD1, and their comparison to $\Lambda\mathrm{CDM}$. In the
top left panel we plot the linear (dashed) and nonlinear (solid) matter power
spectrum as a function of scale for $z=0$. In the top right panel we show the
ratio of these power spectra with respect to their corresponding
$\Lambda\mathrm{CDM}$ cases, using the same $\sigma_{8}$ normalization today.

We can see that for JBD1, the linear power spectrum is almost identical to its
standard model counterpart, however the nonlinear power spectrum shows a small
suppression on small scales, kicking off at about $k=0.3\,h\,{\rm Mpc}^{-1}$.
It is important to note that, unlike nDGP, the JBD theory is not endowed with
gravitational screening so any modification on small scales in the nonlinear
power spectrum is inherited from the change in the primordial amplitude,
$A_{\rm s}$, in the linear power spectrum. In order for $\sigma_{8}$ to be the
same in $\Lambda\mathrm{CDM}$ as in JBD we need a higher primordial amplitude
$A_{\rm s}$. A higher $A_{\rm s}$ in $\Lambda\mathrm{CDM}$ means nonlinear
structures form faster and get more massive, which in the matter power-
spectrum then translates into an even higher amplitude on nonlinear scales
(i.e. it increases more than linear theory predicts) which gives rise to the
small-scale suppression seen in Fig. 1. This effect can also be found in pure
$\Lambda\mathrm{CDM}$ if one considers the ratio of the matter power spectrum
for two different values of $A_{\rm s}$. The result will then be constant on
large scales, but shows an enhancement or suppression on small scales. If we
instead showed results with the same primordial amplitude in JBD as in
$\Lambda\mathrm{CDM}$ then the JBD power-spectrum would show a certain
enhancement to $\Lambda\mathrm{CDM}$ on linear scales and an ever bigger
enhancement on nonlinear scales, which is more in line with naive expectations
of having a larger gravitational constant leading to an enhancement in the
matter power-spectrum. The signature we see in the matter power-spectrum
naturally impacts the lensing spectrum, which shows a similar suppression at
small scales, given that for a given multipole and redshift bin, it is just
proportional to the matter power spectrum as defined in Eq. (30).

####  3.3.2 Dvali-Gabadadze-Porrati gravity

We will model the nonlinear power spectrum using the halo model reaction
(Cataneo et al., 2019) which has been shown to agree with N-body simulations
at the $2\%$ level down to $k=3h\,{\rm 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 coming from modified gravity through the so-called reaction
$\mathcal{R}(k,z)$, which employs both 1-loop perturbation theory and the halo
model. The nonlinear power spectrum in nDGP is given by the product

$P_{\rm NL}(k,z)=\mathcal{R}(k,z)P^{\rm pseudo}_{\rm NL}(k,z)\,,$ (46)

where the _pseudo_ power spectrum is a spectrum where 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. In this work
we use the halofit formula of Takahashi et al. (2012) to model $P^{\rm
pseudo}_{\rm NL}(k,z)$ by supplying the formula with a nDGP input linear
matter power spectrum.

The halo model reaction, $\mathcal{R}(k,z)$, is given by a corrected ratio of
target-to-pseudo halo model spectra

$\mathcal{R}(k,z)=\frac{\\{[1-\mathcal{E}(z)]e^{-k/k_{\star}(z)}+\mathcal{E}(z)\\}P_{\rm
2H}(k,z)+P_{\rm 1H}(k,z)}{P_{\rm hm}^{\rm pseudo}(k,z)}\,.$ (47)

The components are given explicitly as

$\displaystyle P_{\rm hm}^{\rm pseudo}(k,z)=$ $\displaystyle P_{\rm
2H}(k,z)+P_{\rm 1H}^{\rm pseudo}(k,z),$ (48) $\displaystyle\mathcal{E}(z)=$
$\displaystyle\lim_{k\rightarrow 0}\frac{P_{\rm 1H}(k,z)}{P_{\rm 1H}^{\rm
pseudo}(k,z)},$ (49) $\displaystyle k_{\rm\star}(z)=$
$\displaystyle-\bar{k}\left\\{\ln\left[\frac{A(\bar{k},z)}{P_{\rm
2H}(\bar{k},z)}-\mathcal{E}(z)\right]-\ln\left[1-\mathcal{E}(z)\right]\right\\}^{-1}\,,$
(50)

where

$A(k,z)=\frac{P_{\rm 1-loop}(k,z)+P_{\rm 1H}(k,z)}{P^{\rm pseudo}_{\rm
1-loop}(k,z)+P_{\rm 1H}^{\rm pseudo}(k,z)}P_{\rm hm}^{\rm pseudo}(k,z)-P_{\rm
1H}(k,z)\,.$ (51)

$P_{\rm 2H}(k,z)$ is the 2-halo contribution which we can approximate with the
(nDGP) linear power spectrum, $P_{\rm L}(k,z)$ (see for a review on the halo
model Cooray & Sheth, 2002). $P_{\rm 1H}(k,z)$ and $P_{\rm 1H}^{\rm
pseudo}(k,z)$ are the 1-halo contributions to the power spectrum as predicted
by the halo model, with and without modifications to the standard
$\Lambda\mathrm{CDM}$ spherical collapse equations, respectively. We remind
the reader that by definition, the pseudo cosmology has no nonlinear
beyond-$\Lambda\mathrm{CDM}$ modifications. Similarly, $P_{\rm 1-loop}(k,z)$
and $P_{\rm 1-loop}^{\rm pseudo}(k,z)$ are the 1-loop predictions with and
without nonlinear modifications to $\Lambda\mathrm{CDM}$, respectively. As in
previous works, the limit in Eq. (49) is taken to be at $k=0.01\,h\,{\rm
Mpc}^{-1}$ and we compute $k_{\star}$ using $\bar{k}=0.06\,h\,{\rm Mpc}^{-1}$.
This scale is chosen such that the 1-loop predictions are sufficiently
accurate at all the redshifts considered.

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

As in previous works, in our halo model calculations we use a Sheth-Tormen
mass function (Sheth & Tormen, 1999, 2002), a power-law concentration-mass
relation (see for example Bullock et al., 2001) and an NFW halo density
profile (Navarro et al., 1997). These calculations are performed using the
publicly available ReACT code (Bose et al., 2020). We refer the interested
reader 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 vs. $\Lambda\mathrm{CDM}$ as a function of scale for $z=0$. Here we
have chosen for $\Omega_{\rm rc}$ the value $\Omega_{\rm rc}=0.25$ that is
farther away from $\Lambda\mathrm{CDM}$. We will refer to this choice as nDGP1
model. The ratio of these power spectra with respect to their corresponding
$\Lambda\mathrm{CDM}$ cases has been done using the same $\sigma_{8}$
normalization today. We can see that for nDGP1, the linear matter power
spectrum is identical to the $\Lambda\mathrm{CDM}$ one. However, the nonlinear
matter power spectrum shows a suppression starting from scales $k=0.1\,h\,{\rm
Mpc}^{-1}$, which reaches about $10\%$ at scales of about $k=1.0\,h\,{\rm
Mpc}^{-1}$. This impacts of course the lensing spectrum, which shows a similar
suppression at small angular scales, given that for a given multipole and
redshift bin, it is just proportional to the matter power spectrum as defined
in Eq. (30).

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

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

####  3.3.3 $k$-mouflage gravity

For the KM nonlinear matter power spectrum, we use a similar analytical
approach that combines 1-loop perturbation theory with a halo model. This
method was introduced in Valageas et al. (2013) for the standard
$\Lambda\mathrm{CDM}$ cosmology. As in usual halo models, it splits the matter
power spectrum over 2-halo and 1-halo contributions, as in Eq. (48), but it
uses a Lagrangian framework to evaluate these two contributions in terms of
the pair-separation probability distribution. This provides a Lagrangian-space
regularisation of perturbation theory, which matches standard perturbation
theory up to 1-loop order and includes a partial resummation of higher-order
terms, such that the pair-separation probability distribution is positive and
normalised to unity at all scales. Thus, the two-halo term goes beyond the
Zeldovich approximation by including a non-zero skewness, which enables the
consistency with standard perturbation theory up to 1-loop order, as well as a
simple Ansatz for higher-order cumulants inspired by the adhesion model. The
1-halo term includes a counter-term that ensures its falloff at low $k$, in
agreement with the conservation of matter and momentum,

$P_{\rm
1H}(k)=\int_{0}^{\infty}\frac{\mathrm{d}\nu}{\nu}f(\nu)\frac{M}{\bar{\rho}(2\pi)^{3}}\left[\tilde{u}_{M}(k)-\tilde{W}(kq_{M})\right]^{2}\,,$
(52)

where $\tilde{u}_{M}(k)$ is the normalised Fourier transform of the halo
radial profile, $\tilde{W}(kq_{M})$ is the normalised Fourier transform of the
top-hat of Lagrangian radius $q_{M}$, and $f(\nu)$ 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 2-halo and 1-halo
contributions within a Lagrangian framework, which by construction satisfies
the conservation of matter. Additionally, $\nu$ is defined as
$\delta_{c}/\sigma(M)$, where $\delta_{\rm c}$ is the linear density contrast
associated to a nonlinear density contrast of 200, and $\sigma(M)$ is the root
mean square of the linear density contrast at mass scale $M$. The threshold
$\delta_{\rm c}$ is sensitive to the modification of gravity. In principles 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 consider here, screening is negligible beyond
$1\,h^{-1}$ Mpc (clusters are not screened), so that $\delta_{\rm c}$ is
independent of $M$.

This nonlinear modelling was extended in Brax & Valageas (2013) to several
modified-gravity scenarios ($f(R)$ theories, dilaton and symmetron models),
and in Brax & Valageas (2014); Brax et al. (2015) to KM models. This involves
the computation of the impact of modified gravity on both the linear and
1-loop contributions to the matter power spectrum, which enter the 2-halo
term, and on the spherical collapse dynamics, which enter the 1-halo term
through the halo mass function. Therefore, this nonlinear modelling exactly
captures the modification of gravity both up to 1-loop order and at the level
of the fully nonlinear spherical collapse. The Ansatz for the partial
resummation of higher-order terms is the same as for the standard
$\Lambda\mathrm{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
$|f_{R_{0}}|=10^{-4}$, $10^{-5}$ and $10^{-6}$), shows that this approach
captures the relative deviation from the $\Lambda\mathrm{CDM}$ power spectrum
at $z=0$ up to $k=3\,h\,{\rm Mpc}^{-1}$. We have not compared this recipe with
numerical simulations of the KM theories yet, which have only recently been
performed (Hernández-Aguayo et al., 2022). However, we expect at least the
same level of agreement, as the KM models are simpler and closer to the
$\Lambda\mathrm{CDM}$ cosmology than the $f(R)$ models (as in
$\Lambda\mathrm{CDM}$, the linear growth rate is scale independent).

We can see in Fig. 1 the comparison between the KM1 model defined with
$\epsilon_{2,0}=-0.04$ and its $\Lambda\mathrm{CDM}$ counterpart, normalised
to the same $\sigma_{8}$ and cosmological parameters at $z=0$, for the matter
and cosmic shear power spectra. We can notice that the ratio of the linear
matter power spectra increases on intermediate scales ($10^{-2}<k<1\,h\,{\rm
Mpc}^{-1}$): this is due to the change of the background evolution, which
slightly affects the scales at matter-radiation equality. Indeed, $H(z)$ is
greater for KM1 at high redshift (Benevento et al., 2019), which implies a
smaller Hubble radius at matter-radiation equality. This in turns means that
the change of the slope of the linear power spectrum, from $P_{\rm L}(k)\sim
k^{n_{\rm s}}$ to $P_{\rm L}(k)\propto k^{n_{\rm s}-4}$, is shifted to smaller
scales, and whence to higher $k$. This implies that the ratio $P^{\rm mm}_{\rm
KM1}/P^{\rm mm}_{\Lambda\mathrm{CDM}}$ 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 $\sigma_{8}$ and
the ratio $P^{\rm mm}_{\rm KM1}/P^{\rm mm}_{\Lambda\mathrm{CDM}}$ 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
$\Lambda\mathrm{CDM}$, therefore the ratio $P^{\rm mm}_{\rm KM1}/P^{\rm
mm}_{\Lambda\mathrm{CDM}}$ goes to a constant at very low $k$ and at very high
$k$, where the slopes of both linear power spectra are equal to $n_{\rm s}$ or
$n_{\rm s}-4$ (i.e., outside of the intermediate scales where the slope of the
linear power spectrum slowly runs with $k$). This flat ratio can be clearly
seen in Fig. 1 at low $k$. The amplitude of the deviation is set by the
parameter $\epsilon_{2,0}$. There is a significant multiplicative factor, of
about $7$, as it reaches $30\%$ for $\epsilon_{2,0}=-0.04$. This 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 $z=0$, even
though $\epsilon_{2}=-0.04$ we have $\epsilon_{1}\simeq 0.12$, whereas at
linear order in the scalar field dynamics we should have expected that
$\epsilon_{1}=-\epsilon_{2}$. As the growth of structures is actually
sensitive to $\epsilon_{1}$, this leads to a stronger effect than could be
expected from the value of $\epsilon_{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 behavior of the ratio $P^{\rm mm}_{\rm KM1}/P^{\rm
mm}_{\Lambda\mathrm{CDM}}$, another impact of the slight modification of the
background is to shift the baryon acoustic oscillations. This leads to the
oscillations seen at BAO scales, $k\sim 0.01\,h\,{\rm Mpc}^{-1}$. Nonlinear
effects lead to a decrease of the ratio $P^{\rm mm}_{\rm KM1}/P^{\rm
mm}_{\Lambda\mathrm{CDM}}$ 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
virialized halos 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 happens to occur at somewhat larger scales for
the $\Lambda\mathrm{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 $z=0$ for KM1 as compared with
$\Lambda\mathrm{CDM}$, as seen in the upper right panel, the lensing power is
increased, as seen in the lower right panel. This is due to the coupling
function $A(\varphi)$, which is greater than unity at $z>0$ and leads to a
higher value of the phenomenological function $\Sigma$, see Eq. (24), or in
other words to a greater Newton constant at high redshift.

Figure 1: Large-scale structure observables for the different MG models
considered in this work evaluated at the fiducial values of the model
parameters used for the Fisher matrix analysis. Top Left: The linear (dashed)
and nonlinear (solid lines) matter power spectrum $P_{\delta\delta}(z,k)$
entering Eqs. (36) and (30), respectively, evaluated at redshift $z=0$, for
KM1 ($\epsilon_{2,0}=-0.04$, blue), JBD1 ($\omega_{\rm BD}=800$, red) and
nDGP1 ($\Omega_{\rm rc}=0.25$, green). Top Right: Ratio of the matter power
spectra for the above mentioned models with respect to their
$\Lambda\mathrm{CDM}$ counterpart with the same value of $\sigma_{8}$ today,
both for the linear (dashed) and nonlinear (solid) cases. While for nDGP1 and
JBD1 the linear ratio remains constant, in the case of KM1 the ratio shows
large-scale damping and oscillatory features that appear at scales around
$k\approx 0.1\,h\,\mathrm{Mpc}^{-1}$. Bottom Left: Cosmic shear (WL) angular
power spectra for the first bin $i=j=1$, defined in Eq. (30) for the three
models considered above, with the same legend labels. Bottom Right: Ratio of
the cosmic shear $C_{i,j}(\ell)$ to their $\Lambda\mathrm{CDM}$ counterpart,
for the bin combination $i=j=1$. For JBD1 and nDGP1 the suppression in the
nonlinear power spectrum can also be seen as a suppression of the lensing
angular spectra. For KM1, however, although the nonlinear power spectrum is
also decreased at small scales the lensing enhancement function $\Sigma$
defined in Eq. (24) is larger than 1 at $z>0$, giving an overall increase in
the lensing $C_{i,j}(\ell)$ when compared to $\Lambda\mathrm{CDM}$.

##  4 Survey specifications and analysis method

We will vary as base cosmological parameters:

$\bm{\Theta}=\\{\Omega_{\rm m,0},\,\Omega_{\rm b,0},\,h,\,n_{\rm
s},\,\sigma_{8}\\}\,,$ (53)

for which we choose the following fiducial values:

$\displaystyle\Omega_{\rm m,0}$ $\displaystyle=0.315\,,$
$\displaystyle\Omega_{\rm b,0}$ $\displaystyle=0.049\,,$ $\displaystyle h$
$\displaystyle=0.6737\,,$ $\displaystyle n_{\rm s}$ $\displaystyle=0.966\,,$
(54)

where $\Omega_{\rm m,0}=\Omega_{\rm c,0}+\Omega_{\rm b,0}+\Omega_{\rm\nu,0}$.
To fix the fiducial value for $\sigma_{8}$ we will use the same initial
amplitude of primordial perturbations for all models, namely $A_{\rm
s}=2.09681\times 10^{-9}$. Note that the fiducial cosmology includes massive
neutrinos with total mass of $\Sigma m_{\nu}=0.06\,\mathrm{eV}$, but we keep
$\Sigma m_{\nu}$ fixed in the following Fisher matrix analysis. In the
following we will specify the $\sigma_{8}$ value for each model along with the
fiducial values of the model parameters:

  1. 1.

JBD

$\displaystyle\bm{\Theta}_{\rm fid,1}$
$\displaystyle=\\{\sigma_{8}=0.816,\,\log_{10}\omega_{\rm
BD}=2.90309\\}\,\quad({\rm JBD1});\,$ $\displaystyle\bm{\Theta}_{\rm fid,2}$
$\displaystyle=\\{\sigma_{8}=0.812,\,\log_{10}\omega_{\rm
BD}=3.39794\\}\,\quad({\rm JBD2}).$ (55)

These values correspond to $\omega_{\rm BD}=800$ and $\omega_{\rm BD}=2500$
respectively. In these cases the initial value of the scalar field has been
fixed to $\phi_{\rm ini}=1/G_{\rm N}$.

  2. 2.

nDGP

$\displaystyle\bm{\Theta}_{\rm fid,1}$
$\displaystyle=\\{\sigma_{8}=0.8690,\,\log_{10}\Omega_{\rm
rc}=-0.60206\\}\,\quad({\rm nDGP1});\,$ $\displaystyle\bm{\Theta}_{\rm fid,2}$
$\displaystyle=\\{\sigma_{8}=0.8105,\,\log_{10}\Omega_{\rm
rc}=-6\\}\,\quad({\rm nDGP2}).$ (56)

These values correspond to $\Omega_{\rm rc}=0.25$ and $\Omega_{\rm
rc}=10^{-6}$, respectively for nDGP1 and nDGP2.

  3. 3.

KM

$\displaystyle\bm{\Theta}_{\rm fid,1}$
$\displaystyle=\\{\sigma_{8}=0.994,\,\epsilon_{2}=-0.04\\}\,\quad({\rm
KM1});\,$ $\displaystyle\bm{\Theta}_{\rm fid,2}$
$\displaystyle=\\{\sigma_{8}=0.813,\,\epsilon_{2}=-0.0001\\}\,\quad({\rm
KM2}).$ (57)

We follow EC19 for setting up the specifics of the photometric probes as we
briefly summarise here. The sources are split in 10 equi-populated redshift
bins whose limits are obtained from the redshift distribution

$n(z)\propto\left(\frac{z}{z_{0}}\right)^{2}\,\text{exp}\left[-\left(\frac{z}{z_{0}}\right)^{3/2}\right]\,,$
(58)

with $z_{0}=0.9/\sqrt{2}$ and the normalisation set by the requirement that
the surface density of galaxies is $\bar{n}_{g}=30\,\mathrm{arcmin}^{-2}$.
This is then convolved with the sum of two Gaussians to account for the effect
of photometric redshift. Galaxy bias is assumed constant within each redshfit
bin, with fiducial values $b_{i}=\sqrt{1+\bar{z}_{i}}$, and $\bar{z}_{i}$ the
bin centre. Any possible scale dependence of the bias introduced by GR is
taken to be negligible.

As in EC19, we consider a Gaussian only covariance whose elements are given by

$\text{Cov}\left[C_{ij}^{\rm AB}(\ell),C_{kl}^{\rm
CD}(\ell^{\prime})\right]=\frac{\delta_{\ell\ell^{\prime}}^{\rm
K}}{(2\ell+1)f_{\rm sky}\Delta\ell}\\\ \times\left\\{\left[C_{ik}^{\rm
AC}(\ell)+N_{ik}^{\rm AC}(\ell)\right]\left[C_{jl}^{\rm
BD}(\ell^{\prime})+N_{jl}^{\rm BD}(\ell^{\prime})\right]\right.\\\
+\left.\left[C_{il}^{\rm AD}(\ell)+N_{il}^{\rm
AD}(\ell)\right]\left[C_{jk}^{\rm BC}(\ell^{\prime})+N_{jk}^{\rm
BC}(\ell^{\prime})\right]\right\\}\,,$ (59)

where upper-case (lower-case) Latin indexes run over WL, GCph (all tomographic
bins), $\delta_{\ell\ell^{\prime}}^{\rm K}$ is the Kronecker delta symbol
coming from the lack of correlation between different multipoles
$(\ell,\ell^{\prime})$, $f_{\rm sky}\simeq 0.36$ is the survey sky fraction,
and $\Delta\ell$ denotes the width of the $100$ logarithmic equi-spaced
multipole bins. Note that, for the observables of interest here, we can take a
white noise so that it is

$\displaystyle N_{ij}^{\rm LL}(\ell)$ $\displaystyle=\frac{\delta_{ij}^{\rm
K}}{\bar{n}_{i}}\sigma_{\epsilon}^{2}\,,$ (60) $\displaystyle N_{ij}^{\rm
GG}(\ell)$ $\displaystyle=\frac{\delta_{ij}^{\rm K}}{\bar{n}_{i}}\,,$ (61)
$\displaystyle N_{ij}^{\rm GL}(\ell)$ $\displaystyle=0\,,$ (62)

where $\sigma_{\epsilon}^{2}=0.3^{2}$ is the variance of observed
ellipticities.

We keep following EC19 for evaluating the Fisher matrix
$F_{\alpha\beta}(z_{i})$ for the observed galaxy power spectrum. Here,
$\alpha$ and $\beta$ run over the cosmological parameters of the set
$\bm{\Theta}$, the index $i$ labels the redshift bin, each respectively
centred in $z_{i}=\\{1.0,\,1.2,\,1.4,\,1.65\\}$, whose widths are $\Delta
z=0.2$ for the first three bins and $\Delta z=0.3$ for the last bin. As a
difference with respect to EC19, we adopt here the direct derivative approach
varying the observed galaxy power spectrum with respect to the cosmological
parameters directly. We include two additional redshift-dependent parameters
$\ln b\sigma_{8}(z_{i})$ and $P_{\rm s}(z_{i})$ that we marginalise over.
Spectroscopic galaxy bias, $b(z)$, and expected number density of H$\alpha$
emitters, $n(z)$, are as reported in Table 3 of EC19.

For all probes, we consider an optimistic and a pessimistic scenario defined
according to the specifications for WL, GCph and GCsp in Table 1. Note that,
for GCsp, we have added a third scenario, referred to as quasi-linear, fixing
the maximum wave number to $k_{\rm max}=0.15\,h\,{\rm Mpc}^{-1}$. We explore
this more conservative case since the underlying matter power spectrum of Eq.
(36) that we are using in our observed galaxy power spectrum recipe is a
linear one. Nonlinear corrections actually start playing a role for
wavenumbers larger than $k=0.1\,h\,{\rm Mpc}^{-1}$ for the redshifts in
consideration (see Taruya et al., 2010). Cutting at a lower $k_{\rm max}$
should avoid any bias induced by neglecting nonlinear corrections below this
scale. However, such a severe cut also removes part of the information encoded
in GCsp so that we want to quantify how this affects the constraints. In all
scenarios, we fix the $\sigma_{\rm v}$ = $\sigma_{\rm p}$ nuisance parameter
for GCsp to the values directly calculated from Eq. (38) for the fiducial
cosmological parameters.

Table 1: Euclid survey specifications for WL, GCph and GCsp. Survey area | $A_{\rm survey}$ | $15\,000\,\deg^{2}$
---|---|---
WL
Number of photo-$z$ bins | $N_{z}$ | 10
Galaxy number density | $\bar{n}_{\rm gal}$ | $30\,\mathrm{arcmin}^{-2}$
Intrinsic ellipticity dispersion | $\sigma_{\epsilon}$ | 0.30
Minimum multipole | $\ell_{\rm min}$ | 10
Maximum multipole | $\ell_{\rm max}$ |
Maximum– Pessimistic |  | $1500$
Maximum– Optimistic |  | $5000$
GCph
Number of photo-$z$ bins | $N_{z}$ | 10
Galaxy number density | $\bar{n}_{\rm gal}$ | $30\,\mathrm{arcmin}^{-2}$
Minimum multipole | $\ell_{\rm min}$ | 10
Maximum multipole | $\ell_{\rm max}$ |
Maximum– Pessimistic |  | $750$
Maximum– Optimistic |  | $3000$
GCsp
Number of spectro-$z$ bins | $n_{z}$ | 4
Centres of the bins | $z_{i}$ | $\\{1.0,\,1.2,\,1.4,\,1.65\\}$
Error on redshift | $\sigma_{0,z}$ | 0.001
Minimum scale | $k_{\rm min}$ | $0.001\,h\,{\rm Mpc}^{-1}$
Maximum scale | $k_{\rm max}$ |
Maximum– Quasi-linear |  | $0.15\,h\,{\rm Mpc}^{-1}$
Maximum– Pessimistic |  | $0.25\,h\,{\rm Mpc}^{-1}$
Maximum– Optimistic |  | $0.30\,h\,{\rm Mpc}^{-1}$

We consider both GCsp alone and the combination of all photometric probes:
GCph, WL, and their XC. Moreover, we also consider the full combination of all
Euclid main probes: GCsp, GCph, WL, and the XC between GCph and WL. It is
important to mention that we account for all the mixed terms in the covariance
matrix for the photometric probes such as, for instance,
$\text{Cov}\left[C_{ij}^{\rm LL}(\ell),C_{kl}^{\rm GL}(\ell)\right]$. However,
we follow EC19 in neglecting any correlation between GCsp and the photometric
probes. In the pessimistic scenario, we further impose a cut $z<0.9$ to GCph
and XC to remove any overlap between these probes and GCsp. Note that such a
cut is applied only when we add GCsp to the photometric probes, while it is
not when using photometric probes alone in the pessimistic case. No cut is
ever applied in the optimistic case.

##  5 Results

In this section we discuss the results of the Fisher matrix analysis for the
three models under investigation, JBD, nDGP and KM, and for each of them we
consider two fiducial cases presented in Sect. 4. Additionally, for each case
we consider a quasi-linear, pessimistic and optimistic scenario as discussed
in Sect. 4. In the following we present the 68.3% and 95.4% joint marginal
error contours on the cosmological and model parameters for the optimistic
setting and for completeness we include the contour plots for the pessimistic
scenario in the Appendix A. We stress that we are interested in detecting a
non zero value for the model’s parameters, 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 extra parameter,
namely $\log_{10}X$ with $X=\omega_{\rm BD}$ and $X=\Omega_{\rm rc}$ for JBD
and nDGP respectively. This choice is dictated by the fact that we cannot
perform a Fisher analysis on a parameter varying by order of magnitudes.
Additionally, we need to use as parameter a quantity that is of order one in
the analysis in order to avoid having large differences between the largest
and lowest eigenvalues of the Fisher matrix. For the specific case of nDGP,
the use of the logarithm of the model’s parameter allows us to consider only
its positive viable range, which would not be possible if we performed a
Fisher analysis directly on $\Omega_{\rm rc}$. We note indeed that a prior
like $\Omega_{\rm rc}>0$ cannot be imposed a priori in the Fisher analysis.
Beside discussing the results on $\log_{10}X$ we will 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 parameterisations is nonlinear and the
assumption of Gaussianity is valid for the logarithmic parameterisation only.
In order to obtain the constraints on the parameter $X$ we use:

$X^{(\pm)}=X_{\rm fid}\times 10^{\pm\sigma_{\log_{10}X}}\,.$ (63)

As a consequence the uncertainties on $X$ will be asymmetrical. We now proceed
to discuss in detail each model.

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 models
parameters, $\log_{10}\omega_{\rm BD}$ for JBD, $\log_{10}\Omega_{\rm rc}$ for
nDGP and $\epsilon_{2,0}$ for KM in all the scenarios and combinations of
probes we used. Let us note that we did not consider WL probe alone because we
verified that WL provides constraints roughly of the same order as the GCsp
ones. When we combine with GCph and XC the constraints are significantly
boosted. Additionally, we note that the constraining power of the
spectroscopic sample is very weak and this is due to the fact that we do not
go deep enough in the nonlinear regime and also because when using GC probes
we introduce degeneracies among the MG parameters, amplitudes and bias
parameters. These features characterise all the models.

Figure 2: Marginalised $1\sigma$ errors on the models’ 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\sigma$ errors for GCsp
in the optimistic (dotted area), pessimistic (dashed area) and quasi-linear
(filled area) scenarios (top panel), for WL+XC+GCph in the optimistic and
pessimistic cases (central panel) and for the GCsp+WL+XC+GCph in the the
optimistic and pessimistic cases.

###  5.1 Jordan-Brans-Dicke gravity

There are remarkably strong non-cosmological constraints on JBD gravity, with
current bounds of order $\omega_{\rm BD}>10^{5}$. But, as pointed out above,
it is useful to look at JBD gravity as a long wavelength limit of more general
scalar-tensor theories with a wide range of small scale limits which may be
endowed with gravitational screening. It then makes sense to focus on what we
already know from cosmology, where constraints are of order $\omega_{\rm
BD}>10^{3}$, independently from what we know from non-cosmological systems.
Two further, practical reasons lead us to choose $\omega_{\rm BD}\sim 10^{3}$.
First of all this value is at limits of the range for which we have calibrated
the nonlinear prescription we are using. Second, the likelihood is very flat
for high values of $\omega_{\rm BD}$ which can lead to numerical errors in the
finite differences of the Fisher matrix evaluation for some of the
observables. We thus consider two cases for JBD. The first fiducial with
$\omega_{\rm BD}=800$ (JBD1) and the second one with $\omega_{\rm BD}=2500$
(JBD2). Whereas the first one is nearly compatible with the most recent
constraints from publicly available CMB and LSS data, the second one is well
within the current bounds (Ballardini et al., 2020; Joudaki et al., 2022;
Ballardini et al., 2022). Table 2 summarises the relative (with respect to the
fiducial values) uncertainties for the cosmological and model parameters in
the quasi-linear, pessimistic and optimistic Euclid cases. Fig. 3 shows the
68.3% and 95.4% joint marginal error contours on the cosmological model
parameters for the optimistic settings for both JBD1 (left panel) and JBD2
(right panel). Note that as in Casas et al. (2023) for $f(R)$, we have
performed the Fisher matrix analysis on the parameter $\log_{10}\omega_{\rm
BD}$, instead of directly on $\omega_{\rm BD}$, since 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}\omega_{\rm BD}$ as follows:

  * •

JBD1:

    * –

at 35% for the GCsp alone;

    * –

at 3.6% when considering the combination WL+XC+GCph;

    * –

at 3.3% for GCsp+WL+XC+GCph.

  * •

JBD2:

    * –

at 92% for the GCsp alone;

    * –

at 5.1% when considering the combination WL+XC+GCph;

    * –

at 4.6% for GCsp+WL+XC+GCph.

For JBD1 the parameter $\log_{10}\omega_{\rm BD}$ can be measured at high
statistical significance even in the pessimistic case, whereas this does not
hold for JBD2.

We now discuss the forecast errors on the parameter $\omega_{\rm BD}$
propagating the uncertainties according to Eq. (63). For JBD1 we obtain
$\omega_{\rm BD}=800^{+8000}_{-730}$ for GCsp, $\omega_{\rm
BD}=800^{+210}_{-170}$ for WL+XC+GCph and $\omega_{\rm BD}=800^{+200}_{-160}$
for GCsp+WL+XC+GCph; these results correspond to relative errors of $1000\%$,
$27\%$ and $25\%$, respectively. For JBD2, always in the case of optimistic
settings, we obtain $\omega_{\rm BD}=2500^{+1200}_{-820}$ for WL+XC+GCph and
$\omega_{\rm BD}=2500^{+1070}_{-750}$ for GCsp+WL+XC+GCph; these correspond to
relative errors of $48\%$ and $40\%$, respectively (JBD2 is unconstrained by
GC alone).

These results for the optimistic settings demonstrate how Euclid alone could
potentially detect at statistically significant level values of $\omega_{\rm
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 which are compatible with Solar System constraints. Although Euclid
specifications are not exactly the same, these results we find here are
consistent with previous ones in Ballardini et al. (2019), once rewritten in
the coupling to the curvature $\gamma=(4\omega_{\rm BD})^{-1}$.

Finally, let us note that in the pessimistic setting we find that, for some
parameters, the combination of the photometric probes has more power in
constraining than the full combination with spectroscopic probes, see Table 2
and Fig. 6. This is because for the pessimistic setting we adopt a cut in
redshift as described in Sect. 4 to avoid possible overlaps in $z$. Therefore,
we lose the power from GCph at high redshift. This contingency is present only
for the JBD case and not for the other models we study in this work. We
identify the physical mechanism behind it in the fact that JBD does not have a
proper screening mechanism, while, on the contrary, nDGP and KM include it
such that at nonlinear scales they recover $\Lambda\mathrm{CDM}$. This reduces
the importance of the constraining power coming from nonlinearities (which is
captured by photometric probes) for the latter models, being instead important
for the former. We show further results in Appendix A where we present the
results without this $z$-cut.

Figure 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 GCsp, in blue WL+XC+GCph and in green GCsp+WL+XC+GCph.  Table 2:
Forecast $1\sigma$ marginal errors on the cosmological and model parameters,
relative to their corresponding fiducial value, for JBD1 and JBD2 in the
pessimistic, quasi-linear and optimistic cases, using Euclid observations
GCsp, WL+XC+GCph and GCsp+WL+XC+GCph.  JBD1 $\omega_{\rm BD}=800$
---
| $\Omega_{\rm m,0}$ | $\Omega_{\rm b,0}$ | $\log_{10}\omega_{\rm BD}$ | $h$ | $n_{\rm s}$ | $\sigma_{8}$
Pessimistic setting
GCsp $(k_{\rm max}=0.15\,h\,{\rm Mpc}^{-1})$ (quasi-linear) | 3.87% | 6.68% | 87.72% | 1.52% | 2.46% | 1.95%
GC${}_{\rm sp}\,(k_{\rm max}=0.25\,h\,{\rm Mpc}^{-1})$ | 1.59% | 2.35% | 39.49% | 0.53% | 1.24% | 0.98%
WL+XC+GCph | 0.37% | 5.10% | 10.89% | 3.20% | 1.25% | 0.34%
GCsp+WL+XC+GCph | 0.55% | 1.44% | 15.01% | 0.27% | 0.52% | 0.43%
Optimistic setting
GCsp $(k_{\rm max}=0.3\,h\,{\rm Mpc}^{-1})$ | 1.39% | 1.91% | 34.60% | 0.42% | 1.09% | 0.85%
WL+XC+GCph | 0.27% | 3.43% | 3.55% | 1.69% | 0.52% | 0.11%
GCsp+WL+XC+GCph | 0.23% | 1.08% | 3.29% | 0.11% | 0.21% | 0.10%
JBD2 $\omega_{\rm BD}=2500$
Pessimistic setting
GCsp $(k_{\rm max}=0.15\,h\,{\rm Mpc}^{-1})$ (quasi-linear) | 3.88% | 6.74% | 233.99% | 1.53% | 2.45% | 2.01%
GCsp $(k_{\rm max}=0.25\,h\,{\rm Mpc}^{-1})$ | 1.59% | 2.36% | 104.99% | 0.53% | 1.24% | 0.99%
WL+XC+GCph | 0.34% | 5.09% | 17.32% | 3.21% | 1.24% | 0.27%
GCsp+WL+XC+GCph | 1.42% | 1.96% | 74.17% | 0.47% | 1.01% | 0.97%
Optimistic setting
GCsp $(k_{\rm max}=0.3\,h\,{\rm Mpc}^{-1})$ | 1.39% | 1.92% | 92.30% | 0.43% | 1.09% | 0.87%
WL+XC+GCph | 0.25% | 3.40% | 5.12% | 1.70% | 0.56% | 0.10%
GCsp+WL+XC+GCph | 0.20% | 1.09% | 4.55% | 0.11% | 0.21% | 0.09%

###  5.2 Dvali-Gabadadze-Porrati gravity

The two cases considered for the nDGP model are $\log_{10}{\Omega_{\rm
rc}}=-0.60206$ (nDGP1) and $\log_{10}{\Omega_{\rm rc}=-6}$ (nDGP2). We
summarise the 68.3% marginalised errors (relative to the fiducial values) in
Table 3 on the cosmological and model parameters in the optimistic,
pessimistic and quasi-linear case. We show in Fig. 4 the 68.3% and 95.4% joint
marginal error contours on cosmological and model parameters for the
optimistic scenario for both nDGP1 (left panel) and nDGP2 (right panel). In
Fig. 7 we show the same but for the pessimistic setting.

In the optimistic setting, we find that Euclid will be able to constrain
$\log_{10}{\Omega_{\rm rc}}$ as follows:

  * •

nDGP1:

    * –

at 90% for GCsp alone;

    * –

at 20% when considering WL+XC+GCph;

    * –

at 17% for GCsp+WL+XC+GCph.

  * •

nDGP2:

    * –

GCsp has no power in constraining the parameter.

    * –

at 81% for both WL+XC+GCph and GCsp+WL+XC+GCph.

Regardless of the specific case considered (nDGP1 or nDGP2), the power in
constraining $\Omega_{\rm m,0}$, $\Omega_{\rm b,0}$ and $h$ is the same, while
$n_{\rm s}$ and $\sigma_{8}$ are better constrained in the nDGP2 case. For
example, in the optimistic setting, $n_{\rm s}$ for WL+XC+GCph is constrained
at $0.52\%$ for nDGP1 and $0.37\%$ for nDGP2 and they further improve when
considering the full combination of probes, being, respectively, $0.37\%$ and
$0.17\%$. In Fig. 4, in the $\log_{10}{\Omega_{\rm rc}}$–$n_{\rm s}$ panel, we
notice that the two parameters are anti-correlated in the nDGP1 case. This
makes the forecast errors on $n_{\rm s}$ larger than the nDGP2 case where the
anti-correlation disappears. A similar discussion can be done for the
$\sigma_{8}$ parameter. Indeed, in the optimistic setting, for nDGP1 we have
$0.37\%$ for WL+XC+GCph and $0.26\%$ for the full combination of probes, while
for the nDGP2 model we obtain $0.25\%$ and $0.14\%$, respectively.

As previously done for the JBD model, we translate the forecast errors on the
parameter $\log_{10}{\Omega_{\rm rc}}$ to the parameter $\Omega_{\rm rc}$, see
Eq. (63). Quantitatively for the optimistic scenario, for nDGP1 we get
$\Omega_{\rm rc}=0.25^{+0.62}_{-0.18}$ for GCsp, $\Omega_{\rm
rc}=0.25^{+0.08}_{-0.06}$ for WL+XC+GCph and $\Omega_{\rm
rc}=0.25^{+0.07}_{-0.05}$ for GCsp+WL+XC+GCph. These correspond to relative
errors that can reach $247\%$, $32\%$ and $26\%$ respectively. For nDGP2 we
obtain an upper bound that is $\Omega_{\rm rc}<0.072$ in the optimistic
setting. While for the pessimistic case $\Omega_{\rm rc}$ is unconstrained
regardless of the probe we consider because the fiducial parameter is very
close to $\Lambda\mathrm{CDM}$ and we do not have enough nonlinear
information. To examine this, we have considered the linear power spectrum: we
find that the derivative with respect to $\log_{10}\Omega_{\rm rc}$ approaches
0 in the $\Lambda\mathrm{CDM}$ limit, and $P_{\rm
L}^{\prime}|_{10^{-6}}\approx 0.0002$ compared to $P_{\rm
L}^{\prime}|_{0.25}\approx 0.08$, where a prime here denotes a derivative with
respect to $\log_{10}\Omega_{\rm rc}$. This leads us to conclude that the
pessimistic case is still too dominated by linear structure information to
constrain $\Omega_{\rm rc}$ when we are close to the $\Lambda\mathrm{CDM}$
limit. It should also be noted that in Bose et al. (2020) the authors do find
meaningful constraints using cosmic shear with the same scale cuts we employ
in our pessimistic case, but they use a full Markov chain Monte Carlo (MCMC)
approach to a $\Lambda\mathrm{CDM}$ fiducial data vector, sampling in
$\Omega_{\rm rc}$, as well as employing a different covariance and binning
scheme. Additionally it is known that the results obtained by using a Fisher
approach can be different from those computed with MCMC method though the
sampled parameter is the same (Perotto et al., 2006; Wolz et al., 2012).
Therefore a direct comparison of results is not straightforward. Lastly, our
loss in sensitivity in the $\Lambda\mathrm{CDM}$ limit justifies the tight
upper limit in nDGP1 when compared to the nDGP2 one.

Moreover a cut in the maximum multipole and scale at lower values, as it is
for the pessimistic (or quasi-linear) scenario, leads to larger uncertainties
on the model parameter ($\log_{10}\Omega_{\rm rc}$) compared to the optimistic
setting as shown in Table 3 and in Figure 2, where the impact of the different
scales adopted is highlighted. On the contrary, the relative errors on the
cosmological parameters only slightly deteriorate. This means that most of the
power of Euclid in constraining this model comes from the nonlinear scales.

Finally, in this analysis we have performed the forecasts on the logarithm of
$\Omega_{\rm rc}$ following the motivation discussed in Sect. 5, however we
tested that a similar constraining power can be obtained by sampling
$\Omega_{\rm rc}$ directly. For the $3x2$pt in the optimistic setting we
obtain for the nDGP1: a 24% constraint with respect to the 32% obtained on
$\Omega_{\rm rc}$ indirectly from the log sampling.

Figure 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 GCsp, in blue WL+XC+GCph and in green GCsp+WL+XC+GCph. Table 3:
Forecast $1\sigma$ marginal errors on the cosmological and model parameters,
relative to their corresponding fiducial value, for nDGP1 and nDGP2 in the
pessimistic, quasi-linear and optimistic cases, using Euclid observations of
GCsp, WL+XC+GCph and GCsp+WL+XC+GCph.  nDGP1 $\Omega_{\rm rc}=0.25$
---
| $\Omega_{\rm m,0}$ | $\Omega_{\rm b,0}$ | $\log_{10}\Omega_{\rm rc}$ | $h$ | $n_{\rm s}$ | $\sigma_{8}$
Pessimistic setting
GCsp $(k_{\rm max}=0.15\,h\,{\rm Mpc}^{-1})$ (quasi-linear) | 2.23% | 4.30% | 142.20% | 0.57% | 1.69% | 1.55%
GCsp $(k_{\rm max}=0.25\,h\,{\rm Mpc}^{-1})$ | 1.41% | 2.43% | 97.63% | 0.30% | 1.12% | 1.04%
WL+XC+GCph | 1.13% | 5.64% | 49.59% | 2.65% | 0.92% | 1.08%
GCsp+WL+XC+GCph | 0.55% | 1.95% | 40.27% | 0.26% | 0.68% | 0.46%
Optimistic setting
GCsp $(k_{\rm max}=0.3\,h\,{\rm Mpc}^{-1})$ | 1.33% | 2.15% | 90.06% | 0.25% | 0.97% | 0.95%
WL+XC+GCph | 0.30% | 5.08% | 19.59% | 2.33% | 0.52% | 0.37%
GCsp+WL+XC+GCph | 0.25% | 1.28% | 16.99% | 0.15% | 0.33% | 0.26%
nDGP2 $\Omega_{\rm rc}=10^{-6}$
Pessimistic setting
GCsp $(k_{\rm max}=0.15\,h\,{\rm Mpc}^{-1})$ (quasi-linear) | 2.36% | 4.51% | 3020.58% | 0.59% | 1.68% | 2.52%
GCsp $(k_{\rm max}=0.25\,h\,{\rm Mpc}^{-1})$ | 1.46% | 2.51% | 1923.03% | 0.31% | 1.10% | 1.68%
WL+XC+GCph | 0.86% | 5.61% | 398.64% | 2.65% | 0.81% | 0.58%
GCsp+WL+XC+GCph | 0.89% | 1.84% | 378.82% | 0.22% | 1.01% | 0.39%
Optimistic setting
GCsp $(k_{\rm max}=0.3\,h\,{\rm Mpc}^{-1})$ | 1.38% | 2.21% | 1801.03% | 0.26% | 0.97% | 1.54%
WL+XC+GCph | 0.29% | 5.08% | 81.02% | 2.32% | 0.37% | 0.25%
GCsp+WL+XC+GCph | 0.25% | 1.12% | 80.83% | 0.12% | 0.17% | 0.14%

###  5.3 $k$-mouflage

Here we presents our results for the KM model with the two fiducial choices of
the $\epsilon_{2,0}$ parameter (KM1 $\epsilon_{2,0}=-0.04$ and KM2
$\epsilon_{2,0}=-0.0001$). The KM1 model represents an extreme case allowed
within 95.4% confidence interval by current cosmological data, while KM2 is
practically indistinguishable from $\Lambda\mathrm{CDM}$. Shown in Fig. 5 are
the forecast 68.3% and 95.4% joint marginal error contours on cosmological and
model parameters for the optimistic scenario. Fig. 8 shows forecast results
for the same models in the pessimistic scenario. The 68.3% errors for both the
optimistic and pessimistic configurations are reported in Table 4.

In the optimistic setting, we find that Euclid will be able to constrain the
modified-gravity parameter $\epsilon_{2,0}$ as follows:

  * •

KM1 :

    * –

at 5% for the GCsp alone;

    * –

at 0.15% when considering the combination WL+XC+GCph;

    * –

at 0.14% for GCsp+WL+XC+GCph.

  * •

KM2:

    * –

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

    * –

upper bound at $8\times 10^{-4}$ for the GCsp alone;

    * –

upper bound at $7\times 10^{-4}$ when considering the combination WL+XC+GCph;

    * –

upper bound at $6\times 10^{-4}$ for GCsp+WL+XC+GCph.

Considering the KM1 fiducial, the optimistic configuration will allow us to
constrain $\epsilon_{2,0}$ to $\sim 5\%$ accuracy, using GCsp alone. In the
worst case scenario of cutting the observed galaxy power spectrum at $k_{\rm
max}=0.15\,h\,{\rm Mpc}^{-1}$, the error on $\epsilon_{2,0}$ increases to
$\sim 11\%$. The constraining power greatly improves when considering the
combination WL+XC+GCph, with percentage error on $\epsilon_{2,0}$ lowering to
$0.15\%$ in the optimistic case, with no substantial worsening in the
pessimistic scenario. Among the cosmological parameters varied for the KM1
model, $\epsilon_{2,0}$ is the most tightly constrained by the combination
WL+XC+GCph. This shows that this combination of data is well suited to capture
the modification induced by $k$-mouflage to the Poisson and lensing equations,
probing $\mu$ and $\Sigma$ independently. Moreover the Euclid survey will
probe the redshift range $z\gtrsim 1$ where the largest effect determined by
the running of Newton’s constant is expected to manifest (see, e.g., Benevento
et al., 2019). We note that the addition of GCsp to WL+XC+GCph does not
improve constraints on $\epsilon_{2,0}$ in both the pessimistic and optimistic
scenarios.

The percentage error on $\epsilon_{2,0}$ in the
$\Lambda\mathrm{CDM}$-proximate KM2 fiducial model goes from $760\%$ using
GCsp alone to a value of $464\%$ with the full combination of GCsp+WL+XC+GCph,
in the best case scenario. This tells us that, in the optimistic case, the
Euclid survey will be able to detect a deviation from $\Lambda\mathrm{CDM}$
only if $|\epsilon_{2,0}|\gtrsim 6\times 10^{-4}$, improving over the
constraining power of present CMB and LSS data by somewhat more than one order
of magnitude.

From Fig. 5 we notice that $\epsilon_{2,0}$ is anti-correlated with the $h$
parameter. This correlation can be exploited to reduce the Hubble tension, as
noticed in Benevento et al. (2019). This effect determines a lower forecast
error for the $h$ parameter in the KM1 fiducial model, for which
$\epsilon_{2,0}$ is tightly constrained by Euclid probes, especially when the
full data combination is considered. A similar argument applies for the
$\Omega_{\rm m,0}$ parameter which is also better constrained in the KM1 case.
The uncertainties for the other parameters do not vary sensibly with the
fiducial choice of $\epsilon_{2,0}$.

Figure 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 GCsp, in blue WL+XC+GCph and in green GCsp+WL+XC+GCph. Table 4:
Forecast $1\sigma$ marginal errors on the cosmological and model parameters,
relative to their corresponding fiducial value, for KM1 and KM2 in the
pessimistic, quasi-linear and optimistic cases, using Euclid observations of
GCsp, WL+XC+GCph and GCsp+WL+XC+GCph.  KM1 $\epsilon_{2,0}=-0.04$
---
| $\Omega_{\rm m,0}$ | $\Omega_{\rm b,0}$ | $\epsilon_{2,0}$ | $h$ | $n_{\rm s}$ | $\sigma_{8}$
Pessimistic setting
GCsp $(k_{\rm max}=0.15\,h\,{\rm Mpc}^{-1})$ (quasi-linear) | 2.68% | 5.34% | 11.14% | 1.42% | 1.48% | 1.48%
GCsp $(k_{\rm max}=0.25\,h\,{\rm Mpc}^{-1})$ | 1.57% | 3.22% | 5.77% | 0.52% | 1.02% | 1.08%
WL+XC+GCph | 0.36% | 7.36% | 0.22% | 1.40% | 1.93% | 1.10%
GCsp+WL+XC+GCph | 0.34% | 1.73% | 0.22% | 0.14% | 0.60% | 0.52%
Optimistic setting
GCsp $(k_{\rm max}=0.3\,h\,{\rm Mpc}^{-1})$ | 1.48% | 3.08% | 4.99% | 0.41% | 0.89% | 1.04%
WL+XC+GCph | 0.20% | 4.30% | 0.15% | 0.59% | 0.98% | 0.84%
GCsp+WL+XC+GCph | 0.17% | 1.47% | 0.14% | 0.10% | 0.34% | 0.36%
KM2 $\epsilon_{2,0}=-0.0001$
Pessimistic setting
GCsp $(k_{\rm max}=0.15\,h\,{\rm Mpc}^{-1})$ (quasi-linear) | 2.99% | 5.61% | 1264.58% | 0.77% | 1.63% | 1.15%
GCsp $(k_{\rm max}=0.25\,h\,{\rm Mpc}^{-1})$ | 1.60% | 3.08% | 824.38% | 0.38% | 1.07% | 0.98%
WL+XC+GCph | 1.18% | 5.22% | 939.23% | 1.03% | 1.51% | 1.17%
GCsp+WL+XC+GCph | 0.68% | 1.90% | 645.16% | 0.28% | 0.79% | 0.67%
Optimistic setting
GCsp $(k_{\rm max}=0.3\,h\,{\rm Mpc}^{-1})$ | 1.45% | 2.69% | 759.45% | 0.31% | 0.94% | 0.93%
WL+XC+GCph | 0.40% | 3.43% | 658.67% | 0.66% | 0.87% | 0.75%
GCsp+WL+XC+GCph | 0.32% | 1.52% | 464.38% | 0.18% | 0.34% | 0.41%

##  6 Conclusions

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

We derived forecasts from Euclid spectroscopic and photometric primary probes
on the cosmological parameters and the extra parameters of the models,
respectively, $\log_{10}\omega_{\rm BD}$ for JBD, $\log_{10}\Omega_{\rm rc}$
for nDGP and $\epsilon_{2}$ for KM. In order to do so, we applied the Fisher
matrix method to weak lensing (WL); photometric galaxy clustering (GCph);
spectroscopic galaxy clustering (GCsp) and the cross-correlation (XC) between
GCph and WL. For each MG model we considered two fiducial values for the
corresponding model parameter, following the rationale of having a case which
is representative of the $\Lambda\mathrm{CDM}$ limit and another which 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 N-body simulations from modified versions of COLA (Tassev et al., 2013;
Winther et al., 2017) and RAMSES (Teyssier, 2002) for JBD; the halo model
reaction (Cataneo et al., 2019) for nDGP; and an analytical approach that
combined 1-loop perturbation theory with a halo model for KM following Brax &
Valageas (2014); Brax et al. (2015).

When setting the Euclid survey specifications, we defined three scenarios,
characterized by different cuts in the maximum multipole and wavenumber, to
assess the constraining power of nonlinear scales: the quasi-linear with
$k_{\rm max}=0.155\,h\,{\rm Mpc}^{-1}$ for GCsp; the pessimistic with $k_{\rm
max}=0.25\,h\,{\rm Mpc}^{-1}$ for GCsp, $\ell_{\rm max}=1500$ for WL and
$\ell_{\rm max}=750$ for GCph; the optimistic with $k_{\rm max}=0.3\,h\,{\rm
Mpc}^{-1}$ for GCsp, $\ell_{\rm max}=5000$ for WL and $\ell_{\rm max}=3000$
for GCph. In Sect. 5 we have reported in great detail, and discussed, the
results corresponding to these different scenarios; for each case, we
considered three different combinations of data: GCsp only, WL+XC+GCph and the
full GCsp+WL+XC+GCph. For each of such combinations we presented the relative
errors on the fiducial values of the three models as they can give us a
glimpse into possible future detections of non-zero modified gravity
parameters.

With the full set of probes, we found that, in the optimistic scenario, Euclid
alone will be able to constrain the JBD parameter $\omega_{\rm
BD}=800^{+200}_{-160}$ in the fiducial case JBD1 for the full combination of
probes, and in the JBD2 case $\omega_{\rm BD}=2500^{+1070}_{-750}$. These
indicate that Euclid alone will be capable of detecting at statistically
significant level only values of $\omega_{\rm BD}$ which are compatible with
current bounds from CMB and LSS. While we found that in the pessimistic case
the combination of all probes can be less constraining than that of the
photometric ones; this is to be attributed to the cut $z<0.9$ that we apply to
GCph and XC in the full data combination for the pessimistic scenario and
which removes the constraining power of the photometric part at lower
redshift. This cut impacts the JBD model, given that the latter does not have
a screening mechanism and nonlinear scales (which are captured by GCph)
provide therefore quite stringent constraints.

For the nDGP model, we get for the optimistic scenario that $\Omega_{\rm
rc}=0.25^{+0.07}_{-0.05}$ with the full combination in the nDGP1 fiducial
case, while for the nDGP2 fiducial ($\Omega_{\rm rc}=10^{-6}$) we find that
$\Omega_{\rm rc}<0.072$. Furthermore, nDGP2 case is closer to the
$\Lambda\mathrm{CDM}$ model than the KM2 one, which is why nDGP2 has worse
constraints compared to KM2.

For the KM scenario, we get constraints for the full combination of probes in
the fiducial KM1 case of $\epsilon_{2,0}=-0.04\pm 0.0056$. Out of the full set
of cosmological plus model parameters, $\epsilon_{2,0}$ is the most tightly
constrained one for KM1 by the full combination, thanks to the independent
constraints on clustering and lensing, probing the running of Newton’s
constant in the range $z\gtrsim 1$ where the largest effect is expected. The
KM2 fiducial model (with $\epsilon_{2,0}=-0.0001$) is unconstrained for all
the different combinations of probes.

Although it has been shown that for values close to the $\Lambda$CDM limit,
the nDGP2 and KM2 models remain unconstrained in the pessimistic setting, it
is remarkable how Euclid alone can significantly reduce the allowed space for
these models. For the full probe combination (GCsp+WL+XC+GCph) nDGP2 upper
bound at $\Omega_{\rm rc}<0.072$ 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 $|\epsilon_{2,0}|<5.6\times 10^{-4}$ in the
optimistic case and $|\epsilon_{2,0}|<7.5\times 10^{-4}$ in the pessimistic,
improving by more than one order of magnitude, even in the pessimistic case,
the current constraints from a combination of different cosmological datasets
(Benevento et al., 2019). In addition, Euclid alone will provide constraints
on $\omega_{\rm BD}$ which will be tighter than those obtained by a
combination of current CMB and low-redshift data Joudaki et al. (2022).

Finally, our analysis clearly showed that most of the constraining power of
Euclid comes from the 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, like
the one used in this analysis, seen from the substantial differences in
constraining power between the optimistic and pessimistic case. Additionally
we would like to stress that the combination and cross-correlation of future
Euclid data with CMB measurements will be of crucial importance for
cosmological constraints of extended models. As already shown within the
Euclid Collaboration (Euclid Collaboration: Ilić et al., 2022) the
improvements in the constraints when cross-correlation of Euclid with CMB data
is considered can be of the order of two to three and in some cases even
larger. A large effort is currently ongoing to extend this analysis to other
extended models including those in this work.

###### Acknowledgements.

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

##  Appendix A Further results

In this Appendix we show additional results to complement what has been
discussed in Sect. 5. In Figs. 6, 7 and 8 we provide the 68.3% and 95.4% joint
marginal error contours in the pessimistic case on the cosmological parameters
for the JBD, nDGP and KM respectively. We also include Fig. 9 for the JBD case
in the pessimistic scenario without the cut at $z<0.9$. As discussed in Sect.
5.1 the combination of all probes in the pessimistic scenario might result in
worse constraints for some parameters with respect to the combination
WL+XC+GCph because in the former we perform a cut in redshift to avoid overlap
between XC and GCph with GCsp. Comparing Figs. 6 and 5.1 we conclude that the
z-cut removes the constraining power of the photometric part. In detail,
without z-cut we have the following forecast 68.3% relative error for
GCsp+WL+XC+GCph:

  * •

JBD1: $\\{\Omega_{{\rm m,0}},\Omega_{{\rm b,0}},\log_{10}\omega_{{\rm
BD}},h,n_{\rm
s},\sigma_{8}\\}=\\{0.36\%,1.28\%,9.95\%,0.20\%,0.44\%,0.29\%\\}$ ,

  * •

JBD2: $\\{\Omega_{{\rm m,0}},\Omega_{{\rm b,0}},\log_{10}\omega_{{\rm
BD}},h,n_{\rm
s},\sigma_{8}\\}=\\{0.33\%,1.26\%,16.55\%,0.17\%,0.39\%,0.24\%\\}$.

Figure 6: 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 GCsp, in blue WL+XC+GCph and in green GCsp+WL+XC+GCph.

Figure 7: 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 GCsp, in blue WL+XC+GCph and in green GCsp+WL+XC+GCph.

Figure 8: 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 GCsp, in blue WL+XC+GCph and in green GCsp+WL+XC+GCph.

Figure 9: 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 without the redshift cut at $z<0.9$ (see the text for details). In red
GCsp, in blue WL+XC+GCph and in green GCsp+WL+XC+GCph.

## References

  * Aghanim et al. (2020) Aghanim, N., Akrami, Y., Ashdow, Y., et al. 2020, A&A, 641, A6, [Erratum: A&A, 652, C4 (2021)]
  * Avilez & Skordis (2014) Avilez, A. & Skordis, C. 2014, PhRvL, 113, 011101
  * Babichev et al. (2009) Babichev, E., Deffayet, C., & Ziour, R. 2009, IJMPD, 18, 2147
  * Bag et al. (2018) Bag, S., Mishra, S. S., & Sahni, V. 2018, PhRvD, 97, 123537
  * Ballardini et al. (2020) Ballardini, M., Braglia, M., Finelli, F., et al. 2020, JCAP, 10, 044
  * Ballardini et al. (2022) Ballardini, M., Finelli, F., & Sapone, D. 2022, JCAP, 06, 004
  * Ballardini et al. (2016) Ballardini, M., Finelli, F., Umiltà, C., & Paoletti, D. 2016, JCAP, 05, 067
  * Ballardini et al. (2019) Ballardini, M., Sapone, D., Umiltà, C., et al. 2019, JCAP, 05, 049
  * Bardeen (1980) Bardeen, J. M. 1980, PhRvD, 22, 1882
  * Barreira et al. (2015) Barreira, A., Brax, P., Clesse, S., et al. 2015, PhRvD, 91, 063528
  * Barreira et al. (2016) Barreira, A., Sánchez, A. G., & Schmidt, F. 2016, PhRvD, 94, 084022
  * Bellini et al. (2018) Bellini, E., Barreira, A., Frusciante, N., et al. 2018, PhRvD, 97, 023520
  * Bellini et al. (2020) Bellini, E., Sawicki, I., & Zumalacárregui, M. 2020, JCAP, 02, 008
  * Benevento et al. (2019) Benevento, G., Raveri, M., Lazanu, A., et al. 2019, JCAP, 05, 027
  * Berti et al. (2015) Berti, E., Barausse, E., Cardoso, V., et al. 2015, CQGra, 32, 243001
  * Bertotti et al. (2003) Bertotti, B., Iess, L., & Tortora, P. 2003, Natur, 425, 374
  * Blas et al. (2011) Blas, D., Lesgourgues, J., & Tram, T. 2011, JCAP, 07, 034
  * Boisseau et al. (2000) Boisseau, B., Esposito-Farese, G., Polarski, D., & Starobinsky, A. A. 2000, PhRvL, 85, 2236
  * Bose et al. (2020) Bose, B., Cataneo, M., Tröster, T., et al. 2020, MNRAS, 498, 4650
  * Bose et al. (2021) Bose, B., Wright, B. S., Cataneo, M., et al. 2021, MNRAS, 508, 2479
  * Bowcock et al. (2000) Bowcock, P., Charmousis, C., & Gregory, R. 2000, CQGra, 17, 4745
  * Braglia et al. (2020) Braglia, M., Ballardini, M., Emond, W. T., et al. 2020, PhRvD, 102, 023529
  * Braglia et al. (2021) Braglia, M., Ballardini, M., Finelli, F., & Koyama, K. 2021, PhRvD, 103, 043528
  * Brans & Dicke (1961) Brans, C. & Dicke, R. H. 1961, PhRv, 124, 925
  * Brax et al. (2021) Brax, P., Casas, S., Desmond, H., & Elder, B. 2021, Univ, 8, 11
  * Brax et al. (2015) Brax, P., Rizzo, L. A., & Valageas, P. 2015, PhRvD, 92, 043519
  * Brax & Valageas (2013) Brax, P. & Valageas, P. 2013, PhRvD, 88, 023527
  * Brax & Valageas (2014) Brax, P. & Valageas, P. 2014, PhRvD, 90, 023508
  * Brax & Valageas (2016) Brax, P. & Valageas, P. 2016, JCAP, 01, 020
  * Bullock et al. (2001) Bullock, J. S., Kolatt, T. S., Sigad, Y., et al. 2001, MNRAS, 321, 559
  * Cardoso et al. (2008) Cardoso, A., Koyama, K., Seahra, S. S., & Silva, F. P. 2008, PhRvD, 77, 083512
  * Carroll et al. (2004) Carroll, S. M., Duvvuri, V., Trodden, M., & Turner, M. S. 2004, PhRvD, 70, 043528
  * Casas et al. (2023) Casas, S. et al. 2023 [arXiv:2306.11053]
  * Cataneo et al. (2019) Cataneo, M., Lombriser, L., Heymans, C., et al. 2019, MNRAS, 488, 2121
  * Cataneo et al. (2022) Cataneo, M., Uhlemann, C., Arnold, C., et al. 2022, MNRAS, 513, 1623
  * Charmousis et al. (2006) Charmousis, C., Gregory, R., Kaloper, N., & Padilla, A. 2006, JHEP, 10, 066
  * Cheng et al. (2021) Cheng, G., Wu, F., & Chen, X. 2021, PhRvD, 103, 103527
  * Cooray & Sheth (2002) Cooray, A. & Sheth, R. K. 2002, PhR, 372, 1
  * Deffayet (2001) Deffayet, C. 2001, PhLB, 502, 199
  * Dvali et al. (2000) Dvali, G. R., Gabadadze, G., & Porrati, M. 2000, PhLB, 485, 208
  * Eisenstein et al. (2007) Eisenstein, D. J., Seo, H.-j., & White, M. J. 2007, ApJ, 664, 660
  * Euclid Collaboration: Blanchard et al. (2020) Euclid Collaboration: Blanchard, A., Camera, S., Carbone, C., et al. 2020, A&A, 642, A191
  * Euclid Collaboration: Ilić et al. (2022) Euclid Collaboration: Ilić, S., Aghanim, N., Baccigalupi, C., et al. 2022, A&A, 657, A91
  * Ferreira (2019) Ferreira, P. G. 2019, ARA&A, 57, 335
  * Gorbunov et al. (2006) Gorbunov, D., Koyama, K., & Sibiryakov, S. 2006, PhRvD, 73, 044016
  * Hernández-Aguayo et al. (2022) Hernández-Aguayo, C., Ruan, C.-Z., Li, B., et al. 2022, JCAP, 01, 048
  * Hojjati et al. (2011) Hojjati, A., Pogosian, L., & Zhao, G.-B. 2011, JCAP, 1108, 005
  * Hu et al. (2014) Hu, B., Raveri, M., Frusciante, N., & Silvestri, A. 2014, PhRvD, 89, 103530
  * Hu & Sawicki (2007a) Hu, W. & Sawicki, I. 2007a, PhRvD, 76, 104043
  * Hu & Sawicki (2007b) Hu, W. & Sawicki, I. 2007b, PhRvD, 76, 064004
  * Joudaki et al. (2018) Joudaki, S., Blake, C., Johnson, A., et al. 2018, MNRAS, 474, 4894
  * Joudaki et al. (2022) Joudaki, S., Ferreira, P. G., Lima, N. A., & Winther, H. A. 2022, PhRvD, 105, 043522
  * Khoury & Weltman (2004) Khoury, J. & Weltman, A. 2004, PhRvL, 93, 171104
  * Koyama & Maartens (2006) Koyama, K. & Maartens, R. 2006, JCAP, 01, 016
  * Lee et al. (2021) Lee, S. et al. 2021, MNRAS, 509, 4982
  * Lesgourgues (2011) Lesgourgues, J. 2011 [arXiv:1104.2932]
  * Lewis et al. (2000) Lewis, A., Challinor, A., & Lasenby, A. 2000, ApJ, 538, 473
  * Liu et al. (2021) Liu, R., Valogiannis, G., Battaglia, N., & Bean, R. 2021, PhRvD, 104, 103519
  * Lombriser et al. (2009) Lombriser, L., Hu, W., Fang, W., & Seljak, U. 2009, PhRvD, 80, 063536
  * Lombriser et al. (2013) Lombriser, L., Yoo, J., & Koyama, K. 2013, PhRvD, 87, 104019
  * Luty et al. (2003) Luty, M. A., Porrati, M., & Rattazzi, R. 2003, JHEP, 09, 029
  * Ma & Bertschinger (1995) Ma, C.-P. & Bertschinger, E. 1995, ApJ, 455, 7
  * Mead et al. (2021) Mead, A., Brieden, S., Tröster, T., & Heymans, C. 2021, MNRAS, 502, 1401
  * Mead et al. (2016) Mead, A., Heymans, C., Lombriser, L., et al. 2016, MNRAS, 459, 1468
  * Mead et al. (2015) Mead, A., Peacock, J., Heymans, C., Joudaki, S., & Heavens, A. 2015, MNRAS, 454, 1958
  * Mueller et al. (2018) Mueller, E.-M., Percival, W., Linder, E., et al. 2018, MNRAS, 475, 2122
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493
  * Ooba et al. (2017) Ooba, J., Ichiki, K., Chiba, T., & Sugiyama, N. 2017, PTEP, 2017, 043E03
  * Pace et al. (2021) Pace, F., Battye, R., Bellini, E., et al. 2021, JCAP, 06, 017
  * Perotto et al. (2006) Perotto, L., Lesgourgues, J., Hannestad, S., Tu, H., & Wong, Y. Y. Y. 2006, JCAP, 10, 013
  * Raccanelli et al. (2013) Raccanelli, A., Bertacca, D., Pietrobon, D., et al. 2013, MNRAS, 436, 89
  * Raveri et al. (2014) Raveri, M., Hu, B., Frusciante, N., & Silvestri, A. 2014, PhRvD, 90, 043513
  * Raveri et al. (2021) Raveri, M., Pogosian, L., Martinelli, M., et al. 2021 [arXiv:2107.12990]
  * Rossi et al. (2019) Rossi, M., Ballardini, M., Braglia, M., et al. 2019, PhRvD, 100, 103524
  * Schmidt (2009) Schmidt, F. 2009, PhRvD, 80, 123003