# Euclid: Constraining linearly scale-independent modifications of gravity with the spectroscopic and photometric primary probes

N Frusciante, F Pace, V.F Cardone, S Casas, I Tutusaus, M Ballardini, E Bellini, G Benevento, B Bose, P Valageas, et al.

► To cite this version:

N Frusciante, F Pace, V.F Cardone, S Casas, I Tutusaus, et al.. Euclid: Constraining linearly scale-independent modifications of gravity with the spectroscopic and photometric primary probes. Astron.Astrophys., 2024, 690, pp.A133. 10.1051/0004-6361/202347526. hal-04152083

HAL Id: hal-04152083

https://hal.science/hal-04152083v1

Submitted on 26 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Euclid: Constraining linearly scale-independent modifications of gravity with the spectroscopic and photometric primary probes*

N. Frusciante¹,**, D. F. Pace²,³,⁴, V. F. Cardone⁵,⁶, S. Casas⁷, I. Tutusaus⁸,⁹,¹⁰,¹¹, M. Ballardini¹²,¹³,¹⁴,¹⁵,
E. Bellini¹⁶,¹⁷,¹⁸,⁸, G. Benevento¹⁹, B. Bose²⁰,⁸,²¹, D. P. Valageas²², N. Bartolo²³,²⁴,²⁵, P. Brax²²,²⁶, P. G. Ferreira²⁷,
F. Finelli¹⁴,²⁸, K. Koyama²⁹, L. Legrand⁸, L. Lombriser⁸, D. Paoletti¹⁴,¹⁵, M. Pietroni³⁰,³¹, A.
Rozas-Fernández³², Z. Sakr³³,³⁴,³⁵, A. Silvestri³⁶, F. Vernizzi²², H. A. Winther³⁷, N. Aghanim³⁸,
L. Amendola³⁴, N. Auricchio¹⁴, R. Azzollini³⁹, M. Baldi¹³,¹⁴,¹⁵, D. Bonino⁴, E. Branchini⁴⁰,⁴¹, M. Brescia¹,⁴²,
J. Brinchmann⁴³, S. Camera²,³,⁴, V. Capobianco⁴, C. Carbone⁴⁴, J. Carretero⁴⁵,⁴⁶, M. Castellano⁵,
S. Cavuoti⁴²,⁴⁷, A. Cimatti⁴⁸,⁴⁹, R. Cledassou⁵⁰,⁵¹, G. Congedo²⁰, L. Conversi⁵²,⁵³, Y. Copin⁵⁴, L. Corcione⁴,
F. Courbin⁵⁵, M. Cropper³⁹, A. Da Silva⁵⁶,³², H. Degaudenzi⁵⁷, J. Dinis³²,⁵⁶, F. Dubath⁵⁷, X. Dupac⁵³,
S. Dusini²⁴, S. Farrens⁵⁸, S. Ferriol⁵⁴, P. Fosalba¹⁰,⁹, M. Frailis⁵⁹, E. Franceschi¹⁴, S. Galeotta⁵⁹,
B. Gillis²⁰, C. Giocoli¹⁴,⁶⁰, A. Grazian²⁵, F. Grupp⁶¹,⁶², L. Guzzo⁶³,⁶⁴,⁶⁵, S. V. H. Haugan³⁷, W. Holmes⁶⁶,
F. Hormuth⁶⁷, A. Hornstrup⁶⁸, K. Jahnke⁶⁹, S. Kermiche⁷⁰, A. Kiessling⁶⁶, M. Kilbinger⁵⁸, T. Kitching³⁹,
M. Kunz⁸, H. Kurki-Suonio⁷¹, S. Ligori⁴, P. B. Lilje³⁷, I. Lloro⁷², E. Maiorano¹⁴, O. Mansutti⁵⁹,
O. Marggraf⁷³, K. Markovic⁶⁶, F. Marulli⁷⁴,¹⁴,¹⁵, R. Massey⁷⁵, E. Medinaceli¹⁴, M. Meneghetti¹⁵,¹⁴,
G. Meylan⁵⁵, M. Moresco⁷⁴,¹⁴, L. Moscardini⁷⁴,¹⁴,¹⁵, E. Munari⁵⁹, S. M. Niemi⁷⁶, J. Nightingale⁷⁵,
C. Padilla⁴⁵, S. Paltani⁵⁷, F. Pasian⁵⁹, K. Pedersen⁷⁷, W. J. Percival⁷⁸,⁷⁹,⁸⁰, V. Pettorino⁵⁸, G. Polenta⁸¹,
M. Poncet⁵⁰, L. Popa⁸², F. Raison⁶¹, R. Rebolo⁸³,⁸⁴, A. Renzi²³,²⁴, J. Rhodes⁶⁶, G. Riccio⁴², E. Romelli⁵⁹,
R. Saglia⁶¹,⁶², D. Sapone⁸⁵, B. Sartoris⁶², A. Secroun⁷⁰, G. Seidel⁶⁹, C. Sirignano²³,²⁴, G. Sirri¹⁵,
L. Stanco²⁴, C. Surace⁸⁶, P. Tallada-Crespi⁸⁷,⁴⁶, A. N. Taylor²⁰, I. Tereno⁵⁶,⁸⁸, R. Toledo-Moreo⁸⁹,
F. Torradeflot⁸⁷,⁴⁶, E. A. Valentijn⁹⁰, L. Valenziano¹⁴,²⁸, T. Vassallo⁵⁹, G. A. Verdoes Kleijn⁹⁰, Y. Wang⁹¹,
A. Zacchei⁵⁹, G. Zamorani¹⁴, J. Zoubian⁷⁰, and V. Scottez⁹²

*(Affiliations can be found after the references)*

Received 21 July 2023 / Accepted 29 January 2024

## 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 non-linear scales.

**Aims.** We considered three specific models, namely Jordan-Brans-Dicke, 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 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 additional parameters of the models, respectively, $\omega_{BD}$, $\Omega_{\text{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 forecast analysis employs the Fisher matrix method applied to weak lensing (WL); photometric galaxy clustering ($GC_{\text{ph}}$), spectroscopic galaxy clustering ($GC_{\text{sp}}$) and the cross-correlation (XC) between $GC_{\text{ph}}$ and WL. For the Euclid survey specifications, we define three scenarios that are characterised by different cuts in the maximum multipole and wave number, to assess the constraining power of non-linear scales. For each model we considered two fiducial values for the corresponding model parameter.

**Results.** In an optimistic setting at 68.3% confidence interval, we find the following percentage relative errors with Euclid alone: for $\log_{10} \omega_{BD}$, with a fiducial value of $\omega_{BD} = 800$, 27.1% using $GC_{\text{sp}}$ alone, 3.6% using $GC_{\text{ph}}+WL+XC$ and 3.2% using $GC_{\text{ph}}+WL+XC+GC_{\text{sp}}$; for $\log_{10} \Omega_{\text{rc}}$, with a fiducial value of $\Omega_{\text{rc}} = 0.25$, we find 93.4, 20 and 15% respectively; and finally, for $\epsilon_{2,0} = -0.04$, we find 3.4%, 0.15%, and 0.14%. From the relative errors for fiducial values closer to their $\Lambda$CDM limits, we find that most of the constraining power is lost. Our results highlight the importance of the constraining power from non-linear scales.

**Key words.** cosmology: theory – large-scale structure of Universe

<footnote>* This paper is published on behalf of the Euclid Consortium.</footnote>

<footnote>** Corresponding author; noemi.frusciante@unina.it</footnote>
# 1. Introduction

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

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

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

Given the potential of the *Euclid* mission it is imperative to assess its ability to constrain GR. One way of doing this is by assessing how well it will be able to constrain specific extensions of GR, in particular, modified gravity models (MG). MG models are particularly constrained by small-scale experiments such as those in the Solar System where fifth-force effects (Bertotti et al. 2003) and a violation of the equivalence principle (Williams et al. 2012; Touboul et al. 2017) have been thoroughly investigated. As a result, extensions of GR that preserve the equivalence principle in the Earth's environment and exclude large fifth-force effects on Solar System scales must shield small-scales from potential deviations from GR on large cosmological distances. New physical effects would only reveal themselves on large scales where the *Euclid* mission would have the potential to unravel them. Screenings of GR extensions have been broadly classified into three families: the chameleon (Khoury & Weltman 2004), k-mouflage (Babichev et al. 2009) and Vainshtein (Vainshtein 1972) mechanisms (see, for a recent review, Brax et al. 2021). For chameleons, screening occurs in regions of space where Newton's potential is large enough whilst for k-mouflage and Vainshtein, this takes place where the first or second spatial derivatives of Newton's potential are also large enough. One particular example of an extensively

studied theory with the chameleon screening property is $f(R)$ gravity (Carroll et al. 2004; Hu & Sawicki 2007b). It has been shown that the *Euclid* mission, using the combination of spectroscopic and photometric probes, will be able to distinguish this model from $\Lambda$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, that is, k-mouflage and Vainshtein, have not been investigated within the context of the *Euclid* mission, and this will be one of the outcomes of the present work.

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

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

# 2. Linearly scale-independent modified gravity

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

$$ds^2 = -(1 + 2\Psi)c^2 dt^2 + a^2(t)(1 - 2\Phi)\delta_{ij}dx^i dx^j, \quad (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$ is the speed of light. We work in Fourier space, where modes are functions of $t$ and the comoving wave number $k^i$. At linear order, the matter energy-momentum tensor can be decomposed as

$$T_0^0 = -\bar{\rho}(1 + \delta)c^2, \quad T_i^0 = (\bar{\rho} + \bar{\rho}/c^2)cv_i, \quad T_j^i = (\bar{\rho} + \delta\rho)\delta_{ij} + \Sigma_j^i, \quad (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{\rho} + \delta p$ is the pressure with $\bar{\rho}$ the background value, $v_i$ is the peculiar velocity and $\Sigma_i^i$ is the (traceless) anisotropic stress, $\Sigma_i^i = 0$. In the following we work with the scalar component of the matter anisotropic stress and the comoving density perturbation, defined as $(\bar{\rho}c^2 + \bar{\rho})(\hat{k}^i \cdot \hat{k}_j - 1/3\delta_i^j)\sigma = \Sigma_i^i$ with $\hat{k}_i = k_i/k$ and $\bar{\rho}\Delta \equiv \bar{\rho}\delta + 3(aH/k/c^2)(\bar{\rho} + \bar{\rho}/c^2)v$, where $v$ is the velocity potential defined through $v_i = -ik_i v/k$. Here $H \equiv \dot{a}/a$ is the Hubble function and a dot stands for the derivative with respect to the coordinate time. In analogy to the $\Lambda$CDM model where curvature is found to be compatible with zero, we assumed a flat spatial geometry.

In models in which gravity is modified by a scalar field, the relations between the gravitational potentials and the matter perturbations are modified. These deviations from GR can be encoded into two functions, defined as follows:

$$-k^2 \Psi = \frac{4\pi G_N}{c^2} a^2 \mu(a, k) [\bar{\rho}\Delta + 3(\bar{\rho} + \bar{\rho}/c^2)\sigma], \quad (3)$$

$$-k^2 (\Phi + \Psi) = \frac{8\pi G_N}{c^2} a^2 \{\Sigma(a, k) [\bar{\rho}\Delta + 3(\bar{\rho} + \bar{\rho}/c^2)\sigma] \\ - \frac{3}{2}\mu(a, k) (\bar{\rho} + \bar{\rho}/c^2)\sigma\}, \quad (4)$$

where $G_N$ is Newton's gravitational constant. Eqs. (3) and (4) can be obtained using the quasi-static approximation, which considers scales smaller than the horizon and the sound-horizon of the scalar field, where time derivatives become negligible with respect to spatial derivatives. In a given model, it allows us to determine the functional form of $\mu$ and $\Sigma$ analytically (see, e.g. Silvestri et al. 2013; Bellini & Sawicki 2014; Gleyzes et al. 2015; Zucca et al. 2019; Pace et al. 2021). A further function that can be introduced is therefore the function that defines the ratio of the two potentials, $\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). \quad (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, that is, 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 implement the effective field theory formalism for dark energy into the standard CLASS (Blas et al. 2011; Lesgourgues 2011) and CAMB (Lewis et al. 2000) codes, respectively. These codes have been validated as part of an extended code comparison effort (Bellini et al. 2018). Alternatively, one can opt for the quasi-static (QS) limit, that is, scales sufficiently small to be well within the horizon and the sound-horizon of the scalar field, and derive the phenomenological functions analytically. In this case, one can use the MGCAMB patch to CAMB (Zhao et al. 2009; Hojjati et al. 2011; Zucca et al. 2019).

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

## 2.1. Jordan-Brans-Dicke gravity

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

$$S_{\text{BD}} = \int d^4 x \sqrt{-g} \left[ \frac{c^4}{16\pi} \left( \phi R - \frac{\omega_{\text{BD}}}{\phi} g^{\mu\nu} \partial_{\mu} \phi \partial_{\nu} \phi - 2\Lambda \right) + \mathcal{L}_{\text{m}} \right], \quad (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_{\text{BD}}$ is the dimensionless BD parameter, and $\mathcal{L}_{\text{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_{\text{BD}}$; in the limit in which $\omega_{\text{BD}} \to \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. However, cosmological constraints on JBD gravity can have wider implications when we take the view that it is the long-wavelength, low-energy limit of more general scalar-tensor theories (see Joudaki et al. 2022, for example). Furthermore, more general scalar theories may be endowed with gravitational screening, which appears on smaller scales, or alternatively, regions of high density, for example. As a consequence, local constraints will to some extent decouple from more global, or large scale, constraints. Thus, constraints on cosmological scales on JBD theory may cover a broad class of scalar-tensor theories, and furthermore, be independent from model-specific constraints on small-scales. Thus, while simple, JBD gravity is a powerful tool for constraining general classes of scalar-tensor theories using cosmological data.

The modified Einstein field equations are (Clifton et al. 2012)

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

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

$$3H^2 = \frac{8\pi \bar{\rho}}{c^4 \phi} - 3H \frac{\dot{\phi}}{\phi} + \frac{\omega_{\text{BD}}}{2} \frac{\dot{\phi}^2}{\phi^2} + \frac{\Lambda}{\phi}, \quad (8)$$

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

and the scalar field equation of motion reads

$$\Box\phi = \frac{8\pi}{c^4} \left( \frac{T}{3 + 2\omega_{\text{BD}}} \right) - \frac{4\Lambda}{3 + 2\omega_{\text{BD}}}, \quad (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

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

$$G_N \phi (a = 1) = \frac{2\omega_{BD} + 4}{2\omega_{BD} + 3}, \qquad (12)$$

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

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

## 2.2. Dvali-Gabadadze-Porrati braneworld gravity

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

$$S = \frac{c^4}{16\pi G_5} \int_M d^5 x \sqrt{-\gamma} R_5 + \int_{\partial M} d^4 x \sqrt{-g} \left[ \frac{c^4}{16\pi G_N} R + \mathcal{L}_m \right], \qquad (13)$$

where $\gamma$ is the five-dimensional metric and $R_5$ its Ricci curvature scalar. $G_5$ and $G_N$ are the five- and four-dimensional Newton constants, respectively. The matter Lagrangian is denoted with $\mathcal{L}_m$ and is confined to a four-dimensional brane in a five-dimensional Minkowski spacetime. The induced gravity described by the usual four-dimensional Einstein-Hilbert action is responsible for the recovery of the four-dimensional gravity on the brane. The cross-over scale $r_c = G_5/(2G_N)$ is the only parameter of the model and GR is recovered in the limit $r_c \to \infty$.

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

$$H^2 = \pm c \frac{H}{r_c} + \frac{8\pi G_N}{3} \bar{\rho}. \qquad (14)$$

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

The evolution of density and metric perturbations on the brane require solutions of the bulk metric equations. These bulk effects can be encapsulated in an effective 3 + 1 description that uses the combination of any two of the functions $\mu$, $\eta$, and $\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}, \qquad (15)$$

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

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

so 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_N}{c^2} a^2 \bar{\rho} \Delta, \qquad (17)$$

therefore (Lombriser et al. 2013)

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

and hence

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

We chose to parametrise the modification to gravity by $\Omega_{rc} \equiv c^2/(4r_c^2H_0^2)$.

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

There have also been recent attempts to forecast constraints on $\Omega_{rc}$. Liu et al. (2021), using the galaxy cluster abundance from the Simons Observatory and galaxy correlation functions from
a Dark Energy Spectroscopic Instrument (DESI)-like experiment and found $\delta(\Omega_{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_{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. Bose et al. (2020) also forecast for a Large Synoptic Survey Telescope (LSST)-like survey a $1\sigma$ accuracy of 0.08 using cosmic shear alone on a fiducial $\Omega_{rc} \sim 0$.

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

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

## 2.3. k-mouflage gravity

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

$$S = \int d^4 x \sqrt{-\tilde{g}} \left[ \frac{c^4}{16\pi\tilde{G}_N} \tilde{R} + c^2 M^4 K(\tilde{\chi}) \right] + S_m(\psi_i, g_{\mu\nu}), \quad (20)$$

where $\tilde{G}_N$ is the bare Newton constant, $\tilde{R}$ is the Ricci scalar in the Einstein frame, $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_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 M^4}, \quad (21)$$

and $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 degree of deviation from $\Lambda$CDM at the background level and in perturbation theory can be expressed in terms of two time-dependent functions, related to the coupling $A$ and the kinetic function $K$, that is,

$$\epsilon_2 \equiv \frac{d \ln \tilde{A}}{d \ln a}, \quad \epsilon_1 \equiv \frac{2}{\tilde{K}'} \left[ \epsilon_2 \frac{c}{\sqrt{8\pi\tilde{G}_N}} \left( \frac{d\tilde{\varphi}}{d\ln a} \right)^{-1} \right]^2, \quad (22)$$

where over-bars denote background quantities and a prime indicates derivatives with respect to $\tilde{\chi}$. The KM Friedmann equation therefore reads

$$H^2 = \frac{8\pi\tilde{G}_N}{3} \frac{\tilde{A}^2}{(1 - \epsilon_2)^2} \left[ \bar{\rho} + \frac{M^4}{\tilde{A}^4} \left( 2\bar{\chi} \frac{d\tilde{K}}{d\tilde{\chi}} - \tilde{K} \right) \right]. \quad (23)$$

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

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

We recover GR when $\epsilon_1 \to 0$, $\epsilon_2 \to 0$, and $A$ and $K$ are constants, with $M^4 K$ playing the role of the cosmological constant, see Eq. (20).

In addition to the six standard $\Lambda$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 $\tilde{K}(a)$ and $\tilde{A}(a)$. Brax & Valageas (2016) proposed a simple parameterisation 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 that there are no ghosts in the theory. For the coupling function $\tilde{A}(a)$, this reads

$$\tilde{A}(a) = 1 + \alpha_A - \alpha_A \left[ \frac{(\gamma_A + 1) a}{\gamma_A + a} \right]^{\nu_A}, \quad (25)$$

with

$$\nu_A = \frac{3(m-1)}{2m-1}, \quad \alpha_A = -\frac{\epsilon_{2,0}(1+\gamma_A)}{\gamma_A \nu_A}. \quad (26)$$

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

$$\frac{d\tilde{K}'}{d\tilde{\chi}} = \frac{U(a)}{a^3 \sqrt{\tilde{\chi}}}, \quad \sqrt{\tilde{\chi}} = -\frac{\bar{\rho}_0}{M^4} \frac{\epsilon_2 \tilde{A}^4}{2U(-3\epsilon_2 + \frac{d\ln U}{d\ln a})}, \quad (27)$$

which we took to be of the form

$$U(a) \propto \frac{a^2 \ln(\gamma_U + a)}{(\sqrt{a}_{eq} + \sqrt{a}) \ln(\gamma_U + a) + \alpha_U a^2}. \quad (28)$$

The various terms in this expression allowed 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$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, polarisation and lensing, type Ia supernovae, baryon acoustic oscillations (BAO) and local measurements of $H_0$ (Benevento et al. 2019). This gives the 95% C.L. bounds $-0.04 \le \epsilon_{2,0} \le 0$, while the other parameters are unconstrained (as long as they do not become very large). These results are consistent with earlier more qualitative CMB studies (Barreira et al. 2015). The X-ray cluster multiplicity function could provide bounds of the same order (Brax et al. 2015), but more detailed analyses are needed to derive robust constraints.

The fact that $\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 \approx -\epsilon_2$. 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 ΛCDM occur at low redshift, this is mostly determined by the value $\epsilon_{2,0}$ at $z=0$. Therefore, we here focused on the dependence on $\epsilon_{2,0}$, which was left free to vary, while we fixed 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 ΛCDM parameters. The ΛCDM model is recovered when $\epsilon_{2,0} = 0$, independently of the value of the other parameters; for $\epsilon_{2,0} \to 0$, we have:

$$ \bar{A}(a) \to 1, \quad \epsilon_2(a) \to 0, \quad \epsilon_1(a) \to 0, \quad \forall a, \qquad (29) $$

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

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

### 3. Theoretical predictions for Euclid observables

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

#### 3.1. Photometric survey

The observables we considered for the Euclid photometric survey are the angular power spectra $C_{ij}^{XY}(\ell)$ between the probe X in the i-th redshift bin, and the probe Y in the j-th bin, where X refers to weak lensing (WL), photometric galaxy clustering (GC$_{ph}$), or their cross correlation, XC. These are still calculated as in EC19 relying on the Limber approximation and setting to unity the $\ell$-dependent factor in the flat-sky limit. The spectra are then given by

$$ C_{ij}^{XY}(\ell) = c \int_{z_{\min}}^{z_{\max}} dz \frac{W_i^X(z) W_j^Y(z)}{H(z) r^2(z)} P_{\delta\delta}(k_\ell, z), \qquad (30) $$

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

$$ W_i^G(k, z) = b_i(k, z) \frac{n_i(z)}{\bar{n}_i} \frac{H(z)}{c}, \qquad (31) $$

$$ W_i^L(k, z) = \frac{3}{2} \Omega_{m,0} \frac{H_0^2}{c^2} (1+z) r(z) \Sigma(z) \int_z^{z_{\max}} dz' \frac{n_i(z')}{\bar{n}_i} \frac{r(z' - z)}{r(z')} + W_i^{IA}(k, z), \qquad (32) $$

where $n_i(z)/\bar{n}_i$ and $b_i(k, z)$ are the normalised galaxy distribution and the galaxy bias in the i-th redshift bin, respectively, and $W_i^{IA}(k, z)$ encodes the contribution of intrinsic alignments (IA) to the WL power spectrum. The function $\Sigma(z)$ in the WL window function explicitly accounts for the changes in the lensing

potential due to the particular MG theory of interest. Its explicit form for the cases under consideration can be found in Section 2. The impact on the background quantities $H(z)$, $r(z)$ and the matter power spectrum $P_{\delta\delta}(k, z)$, in contrast, are already taken into account by the dedicated Boltzmann solver so that the GC$_{ph}$ window function remains unchanged.

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

$$ W_i^{IA}(k, z) = - \frac{\mathcal{A}_{\text{IA}} C_{\text{IA}} \Omega_{m,0} \mathcal{F}_{\text{IA}}(z) n_i(z) H(z)}{\delta(k, z)/\delta(k, z=0) \bar{n}_i(z) c}, \qquad (33) $$

where

$$ \mathcal{F}_{\text{IA}}(z) = (1+z)^{\eta_{\text{IA}}} \left[ \frac{\langle L \rangle(z)}{L_{\star}(z)} \right]^{\beta_{\text{IA}}}, \qquad (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}_{\text{IA}}$, $\beta_{\text{IA}}$ and $\eta_{\text{IA}}$ are the nuisance parameters of the model, and $C_{\text{IA}}$ is a constant accounting for dimensional units. This model is the same as was used in EC19 since IA takes place on astrophysical scales that are unaffected by modifications to gravity. However, MG has an impact on the growth factor, introducing a possible scale dependence. This is explicitly taken into account in Eq. (33) through the matter perturbation $\delta(k, z)$, which is considered to be scale dependent in this case. This allows us to also consider the scale dependence introduced by massive neutrinos, which was assumed to be negligible in EC19. We nevertheless stress that for the models we considered, the scale dependence is quite small so that the IA is essentially the same as in the GR case.

#### 3.2. Spectroscopic survey

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

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

The observed galaxy power spectrum is given by

$$ P_{\text{obs}}(k, \mu_{\theta}; z) = \frac{1}{q_{\perp}^2(z) q_{\parallel}(z)} \left\{ \frac{[b\sigma_8(z) + f\sigma_8(z)\mu_{\theta}^2]^2}{1 + [f(z)k\mu_{\theta}\sigma_p(z)]^2} \right\} \times \frac{P_{\text{dw}}(k, \mu_{\theta}; z)}{\sigma_8^2(z)} F_z(k, \mu_{\theta}; z) + P_s(z), \qquad (35) $$

where $P_{\text{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}^{\text{lin}}(k; z) e^{-g_{\mu}k^2} + P_{\text{nw}}(k; z) (1 - e^{-g_{\mu}k^2}), \qquad (36) $$
and $P_{\delta\delta}^{\text{lin}}(k; z)$ stands for the linear matter power spectrum. $P_{\text{nw}}(k; z)$ is a no-wiggle power spectrum with the same broad-band shape as $P_{\delta\delta}^{\text{lin}}(k; z)$ but without BAO features (see the discussion below). The function

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

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

$$\sigma_{v}^{2}(z) = \frac{1}{6\pi^{2}} \int dk P_{\delta\delta}^{\text{lin}}(k, z). \quad (38)$$

The curly bracket in Eq. (35) is the redshift-space-distortion (RSD) contribution correcting for the non-linear finger-of-God (FoG) effect, where we defined $b\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)$ (we marginalized over $b\sigma_8(z)$), while $\mu_\theta$ is the cosine of the angle $\theta$ between the wave vector $k$ and the line-of-sight direction $\hat{r}$ and $\sigma_p^2(z) = \sigma_v^2(z)$. Although these parameters were assumed to be the same, they come from two different physical effects, namely large-scale bulk flow for the former and virial motion for the latter.

The factor $F_z$ accounts for the smearing of the galaxy density field along the line of sight due to redshift uncertainties. It is given by

$$F_z(k, \mu_{\theta}; z) = e^{-k^2 \mu_{\theta}^2 \sigma_r^2(z)}, \quad (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 parametrised in terms of the angular diameter distance $D_A(z)$ and the Hubble parameter $H(z)$ as

$$q_{\perp}(z) = \frac{D_A(z)}{D_{A, \text{ref}}(z)}, \quad (40)$$

$$q_{\parallel}(z) = \frac{H_{\text{ref}}(z)}{H(z)}. \quad (41)$$

Due to the Alcock-Paczynski effect, $\mu_\theta$ and $k$ were rescaled in a cosmology-dependent way as a function of the projection along and perpendicular to the line of sight. In the previous Eq. (35), the arguments $\mu$ and $k$ are themselves functions of the true $\mu_{\theta, \text{ref}}$ and $k_{\text{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_s(z)$ is a scale-independent shot-noise term that enters as a nuisance parameter (see EC19).

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

### 3.3. non-linear modelling for weak lensing

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

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

#### 3.3.1. Jordan-Brans-Dicke gravity

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

In the halo model (Cooray & Sheth 2002), the non-linear matter power spectrum can be written as a sum $P_{\text{NL}}(k, z) = P_{\text{1H}}(k, z) + P_{\text{2H}}(k, z)$ where

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

is the so-called two-halo term (correlation between different haloes), and

$$P_{\text{1H}}(k, z) = \int_{0}^{\infty} W^{2}(M, k, z) n(M, z) \, dM, \quad (43)$$

is the so-called one-halo term (correlations between mass-elements within each halo). Above $M$ is the halo mass, $P_{\text{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_{\text{v}}} 4\pi r^2 \frac{\sin(kr)}{kr} \rho(M, r, z) \, dr, \quad (44)$$

where $\rho(M, r, z)$ is the radial matter density profile in a host halo with a mass $M$, and $\bar{\rho}$ is the mean matter density. The halo mass
<figure><img src="image_2.png" /><figcaption>Fig. 1. Large-scale structure observables for the different MG models evaluated at the fiducial values of the model parameters used for the Fisher matrix analysis. Top left: linear (dashed) and non-linear (solid lines) matter power spectrum $P_{\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_{BD} = 800$, red) and nDGP1 ($\Omega_{rc} = 0.25$, green). Top right: ratio of the matter power spectra for the above mentioned models with respect to their ΛCDM counterpart with the same value of $\sigma_8$ today, for the linear (dashed) and non-linear (solid) cases. Bottom left: cosmic shear (WL) angular power spectra for the auto-correlation of the first photometric bin $C_{1,1}^{WL}(\ell)$ (solid line) and the cross-correlation of the first and last bin $C_{1,10}^{WL}(\ell)$ (dotted line), defined in Eq. (30) for the three models considered above, with the same colour labelling. Bottom right: ratio of the cosmic shear $C_{i,j}(\ell)$ to their ΛCDM counterpart, for the bin combination $i = j = 1$ and $i = 1; j = 10$.</figcaption></figure>

is related to the virial radius, $r_v$, via $M = 4\pi r_v^3 \Delta_v(z) \bar{\rho}/3$, where $\Delta_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 that are then fit to data. HM CODE combines this, and the resulting fitting function has, for ΛCDM, been shown to be accurate to 2.5% for scales $k < 10 h \text{ Mpc}^{-1}$ and redshifts $z < 2$ (Mead et al. 2021).

Joudaki et al. (2022) modified the HM CODE, which is able to include the effects of JBD. This was done using a suite of N-body simulations obtained with modified versions of both the CO moving 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 \text{ Mpc}^{-1}$ to generate a large enough ensemble of simulations and to reduce sample variance. On very small-scales, up to $k < 10 h \text{ Mpc}^{-1}$, the RAMSES grid-based hydrodynamical solver with adaptive mesh refinement was used.

The spectra generated by these simulations were then used to calibrate HM CODE. While we did not consider the effects of baryons here, the advantage of HM CODE, over Halofit is that the former is able to capture baryonic feedback at an accuracy level of ~5% up to $k \approx 10 h \text{ Mpc}^{-1}$ (Mead et al. 2021). This will

become an important consideration (Schneider et al. 2020b,a) in future more detailed analyses. To take JBD into account in addition to modifying the expansion history and growth of perturbations, the virialised halo overdensity was modified as

$$ \Delta_v = \Omega_m(z)^{-0.352} \times \left\{ d_0 + (418.0 - d_0) \arctan \left[ (0.001|\omega_{\text{BD}} - 50.0|)^{0.2} \right] \frac{2}{\pi} \right\}, \quad (45) $$

where $d_0 = 320.0 + 40.0z^{0.26}$. This modification has the feature that it reduces to the original ΛCDM HM CODE as $\omega_{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 HM CODE 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_{BD} \gtrsim 50$ and on scales up to $k \approx 10 h \text{ Mpc}^{-1}$. This non-linear prescription was made with past weak-lensing surveys in mind and fitted to a range of $\omega_{BD}$ and scales that we can realistically constrain at the present time. For a parameter inference with actual Euclid data the accuracy here might not be good enough.

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

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

### 3.3.2. Dvali-Gabadadze-Porrati gravity

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

$$P_{\text{NL}}(k, z) = \mathcal{R}(k, z) P_{\text{NL}}^{\text{pseudo}}(k, z), \qquad (46)$$

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

modified linear clustering at $z$ while keeping a non-linear clustering as described by GR, in line with our definition of the pseudo-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_{2\text{H}}(k, z) + P_{1\text{H}}(k, z)}{P_{\text{hm}}^{\text{pseudo}}(k, z)}. \qquad (47)$$

The components are given explicitly as

$$P_{\text{hm}}^{\text{pseudo}}(k, z) = P_{2\text{H}}(k, z) + P_{1\text{H}}^{\text{pseudo}}(k, z), \qquad (48)$$

$$\mathcal{E}(z) = \lim_{k \to 0} \frac{P_{1\text{H}}(k, z)}{P_{1\text{H}}^{\text{pseudo}}(k, z)}, \qquad (49)$$

$$k_{\star}(z) = -\bar{k} \left\{ \ln \left[ \frac{A(\bar{k}, z)}{P_{2\text{H}}(\bar{k}, z)} - \mathcal{E}(z) \right] - \ln[1 - \mathcal{E}(z)] \right\}^{-1}, \qquad (50)$$

where

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

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

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

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

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

As the $\sigma_8$ normalisation is the same for the $\Lambda$CDM reference spectrum, it is essentially the $P_{\text{NL}}^{\text{pseudo}}$ by definition. This means the quantity associated with the solid green 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 one-halo terms at $k > 0.1 h \text{ 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 one-halo term are relatively more suppressed by the NFW density profile (see for example Fig. 9 of Cooray & Sheth 2002, and Eq. (52)). This will cause the 1-halo GR power spectrum to be larger than the DGP one-halo spectrum if the linear clustering in the two cosmologies is the same at the target redshift.

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

### 3.3.3. k-mouflage gravity

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

$$P_{\text{1H}}(k) = \int_{0}^{\infty} \frac{dv}{v} f(v) \frac{M}{\bar{\rho}(2\pi)^3} \left[ \tilde{u}_{M}(k) - \tilde{W}(kq_{M}) \right]^2, \quad (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(v)$ is the normalised halo mass function. This counter-term is not introduced by hand, as it arises directly from the splitting of the matter power spectrum over two-halo and one-halo contributions within a Lagrangian framework, which by construction satisfies the conservation of matter. Additionally, $v$ is defined as $\delta_c/\sigma(M)$, where $\delta_c$ is the linear density contrast associated with a non-linear density contrast of 200, and $\sigma(M)$ is the root mean square of the linear density contrast at mass scale $M$. The threshold $\delta_c$ is sensitive to the

modification of gravity. In principle, it also depends on $M$, as the shells in the spherical collapse are coupled and evolve differently depending on their masses due to non-linear screening. However, for the models that we considered here, screening is negligible beyond $1 h^{-1}$ Mpc (clusters are not screened), so that $\delta_c$ is independent of $M$.

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

Fig. 1 shows the comparison between the KM1 model defined with $\epsilon_{2,0} = -0.04$ and its $\Lambda$CDM counterpart, normalised to the same $\sigma_8$ and cosmological parameters at $z = 0$, for the matter and cosmic shear power spectra. The ratio of the linear matter power spectra increases on intermediate scales ($10^{-2} < k < 1 h \text{ Mpc}^{-1}$): This is due to the change of the background evolution, which slightly affects the scales at matter-radiation equality. $H(z)$ is indeed greater for KM1 at high redshift (Benevento et al. 2019), which implies a smaller Hubble radius at matter-radiation equality. This in turns means that the change in the slope of the linear power spectrum, from $P_L(k) \sim k^{n_s}$ to $P_L(k) \propto k^{n_s-4}$, is shifted to smaller scales, and hence to higher $k$. This implies that the ratio $P_{\text{KM1}}^{\text{mm}}/P_{\Lambda\text{CDM}}^{\text{mm}}$ grows on these intermediate scales, where the linear power spectrum changes slope, because of this delay to higher $k$ of the transition. Because the two spectra have the same $\sigma_8$ and the ratio $P_{\text{KM1}}^{\text{mm}}/P_{\Lambda\text{CDM}}^{\text{mm}}$ grows with $k$, the ratio is below unity at low $k$ and above unity at high $k$. Moreover the linear growth rates are scale independent for both KM and $\Lambda$CDM, and therefore the ratio $P_{\text{KM1}}^{\text{mm}}/P_{\Lambda\text{CDM}}^{\text{mm}}$ reaches a constant at very low $k$ and at very high $k$, where the slopes of both linear power spectra are equal to $n_s$ or $n_s - 4$ (i.e. outside of the intermediate scales where the slope of the linear power spectrum slowly runs with $k$). This flat ratio is clearly shown in Fig. 1 at low $k$. The amplitude of $P_{\text{KM1}}^{\text{mm}}/P_{\Lambda\text{CDM}}^{\text{mm}}$ is set by the parameter $\epsilon_{2,0}$, and it deviates substantially from 1, by about 30%, for a value of $\epsilon_{2,0} = -0.04$. The increase in $P_{\text{KM1}}^{\text{mm}}/P_{\Lambda\text{CDM}}^{\text{mm}}$ with $k$ is partly due to non-linear effects in the scalar field sector, both in the background and in the linear growth rate, which are not completely irrelevant. At $z = 0$, even though $\epsilon_2 = -0.04$ we have $\epsilon_1 \approx 0.12$, whereas at linear order in the scalar field dynamics we would have expected that $\epsilon_1 = -\epsilon_2$. As the growth of structures is 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 behaviour of the ratio $P_{\text{KM1}}^{\text{mm}}/P_{\text{\Lambda CDM}}^{\text{mm}}$, another impact of the slight modification of the background is that the baryon acoustic oscillations are shifted. This leads to the oscillations seen at BAO scales, $k \sim 0.01 h \text{ Mpc}^{-1}$. Non-linear effects lead to a decrease in the ratio $P_{\text{KM1}}^{\text{mm}}/P_{\text{\Lambda CDM}}^{\text{mm}}$ as compared with the linear prediction of mildly non-linear scales. This is due to the loss of information in the power spectrum in the non-linear regime, which can be related to the universal NFW profile of virialised haloes within the halo model. Then, the shape of the non-linear power spectrum at those scales is mostly set by the scale where the non-linear transition takes place, and this occurs at somewhat larger scales for the \Lambda CDM counterpart (i.e. the higher linear power at higher $k$ in the KM1 scenario is mostly erased by the non-linear dynamics).

As usual, the integration along the line of sight leads to a much smoother curve for the cosmic shear angular power spectrum. Even though the 3D non-linear matter power spectrum is reduced at $z=0$ for KM1 as compared with \Lambda CDM, as shown in the upper right panel, the lensing power is increased, as shown in the lower right panel. This is due to the coupling function $A(\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.

## 4. Survey specifications and analysis method

We varied as base cosmological parameters:

$$ \mathbf{\Theta} = \{\Omega_{m,0}, \Omega_{b,0}, h, n_s, \sigma_8\}, \qquad (53) $$

for which we chose the following fiducial values:

$$ \begin{aligned} \Omega_{m,0} &= 0.315, & \Omega_{b,0} &= 0.049, \\ h &= 0.6737, & n_s &= 0.966, \end{aligned} \qquad (54) $$

where $\Omega_{m,0} = \Omega_{c,0} + \Omega_{b,0} + \Omega_{\gamma,0}$. To fix the fiducial value for $\sigma_8$, we used the same initial amplitude of primordial perturbations for all models, namely $A_s = 2.09681 \times 10^{-9}$. The fiducial cosmology includes massive neutrinos with a total mass of $\Sigma m_\nu = 0.06$ eV, but we kept $\Sigma m_\nu$ fixed in the following Fisher matrix analysis. In the following, we specify the $\sigma_8$ value for each model along with the fiducial values of the model parameters,

### 1. JBD

$$ \mathbf{\Theta}_{\text{fid},1} = \{\sigma_8 = 0.816, \log_{10} \omega_{\text{BD}} = 2.90309\} \quad \text{(JBD1);} $$

$$ \mathbf{\Theta}_{\text{fid},2} = \{\sigma_8 = 0.812, \log_{10} \omega_{\text{BD}} = 3.39794\} \quad \text{(JBD2).} \quad (55) $$

These values correspond to $\omega_{\text{BD}} = 800$ and $\omega_{\text{BD}} = 2500$ respectively. In these cases the initial value of the scalar field was fixed to $\phi_{\text{ini}} = 1/G_N$.

### 2. nDGP

$$ \mathbf{\Theta}_{\text{fid},1} = \{\sigma_8 = 0.8690, \log_{10} \Omega_{\text{rc}} = -0.60206\} \quad \text{(nDGP1);} $$

$$ \mathbf{\Theta}_{\text{fid},2} = \{\sigma_8 = 0.8105, \log_{10} \Omega_{\text{rc}} = -6\} \quad \text{(nDGP2).} \quad (56) $$

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

### 3. KM

$$ \begin{aligned} \mathbf{\Theta}_{\text{fid},1} &= \{\sigma_8 = 0.994, \epsilon_2 = -0.04\} \quad \text{(KM1);} \\ \mathbf{\Theta}_{\text{fid},2} &= \{\sigma_8 = 0.813, \epsilon_2 = -0.0001\} \quad \text{(KM2).} \quad (57) \end{aligned} $$

We followed EC19 to set up the specifics of the photometric probes as we briefly summarise here. The sources were split into ten equi-populated redshift bins whose limits were obtained from the redshift distribution

$$ n(z) \propto \left(\frac{z}{z_0}\right)^2 \exp\left[-\left(\frac{z}{z_0}\right)^{3/2}\right], \qquad (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 \text{ arcmin}^{-2}$. This was then convolved with the sum of two Gaussians to account for the effect of photometric redshift. Galaxy bias was assumed to be constant within each redshfit bin, with fiducial values $b_i = \sqrt{1 + \bar{z}_i}$, and $\bar{z}_i$ the bin centre. Any possible scale dependence of the bias introduced by GR was taken to be negligible.

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

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

$$ \qquad + [C_{il}^{\text{AD}}(\ell) + N_{il}^{\text{AD}}(\ell)] [C_{jk}^{\text{BC}}(\ell') + N_{jk}^{\text{BC}}(\ell')] \qquad (61) $$

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

$$ N_{ij}^{\text{LL}}(\ell) = \frac{\delta_{ij}^{\text{K}}}{\bar{n}_i} \sigma_{\epsilon}^2, \qquad (62) $$

$$ N_{ij}^{\text{GG}}(\ell) = \frac{\delta_{ij}^{\text{K}}}{\bar{n}_i}, \qquad (63) $$

$$ N_{ij}^{\text{GL}}(\ell) = 0, \qquad (64) $$

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

We continued to follow EC19 to evaluate 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 $\mathbf{\Theta}$, and the index $i$ labels the redshift bin, each centred in $z_i = \{1.0, 1.2, 1.4, 1.65\}$, whose widths were $\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 adopted the direct derivative approach and directly varied the observed galaxy power spectrum with respect to the cosmological parameters. We included two additional redshift-dependent parameters $\ln b \sigma_8(z_i)$ and $P_s(z_i)$ over which we marginalised. The spectroscopic galaxy bias, $b(z)$, and the expected number density of H\alpha emitters, $n(z)$, are as reported in Table 3 of EC19.

For all probes, we considered an optimistic and a pessimistic scenario defined according to the specifications for WL, GCph and GCsp in Table 1. For GCsp, we added a third scenario, referred to as quasi-linear, for which we fixed the maximum wavenumber to $k_{\text{max}} = 0.15 h \text{ Mpc}^{-1}$. We explored this more
<table><caption>**Table 1.** *Euclid* survey specifications for WL, GC<sub>ph</sub> and GC<sub>sp</sub>.</caption><thead><tr><th>Survey area</th><th>A<sub>survey</sub></th><th>15 000 deg<sup>2</sup></th></tr></thead><tbody><tr><td colspan="3"><strong>WL</strong></td></tr><tr><td>Number of photo-z bins</td><td>N<sub>z</sub></td><td>10</td></tr><tr><td>Galaxy number density</td><td>n&#773;<sub>gal</sub></td><td>30 arcmin<sup>-2</sup></td></tr><tr><td>Intrinsic ellipticity dispersion</td><td>&sigma;<sub>&epsilon;</sub></td><td>0.30</td></tr><tr><td>Minimum multipole</td><td>&ell;<sub>min</sub></td><td>10</td></tr><tr><td>Maximum multipole</td><td>&ell;<sub>max</sub></td><td></td></tr><tr><td>- Pessimistic</td><td></td><td>1500</td></tr><tr><td>- Optimistic</td><td></td><td>5000</td></tr><tr><td colspan="3"><strong>GC<sub>ph</sub></strong></td></tr><tr><td>Number of photo-z bins</td><td>N<sub>z</sub></td><td>10</td></tr><tr><td>Galaxy number density</td><td>n&#773;<sub>gal</sub></td><td>30 arcmin<sup>-2</sup></td></tr><tr><td>Minimum multipole</td><td>&ell;<sub>min</sub></td><td>10</td></tr><tr><td>Maximum multipole</td><td>&ell;<sub>max</sub></td><td></td></tr><tr><td>- Pessimistic</td><td></td><td>750</td></tr><tr><td>- Optimistic</td><td></td><td>3000</td></tr><tr><td colspan="3"><strong>GC<sub>sp</sub></strong></td></tr><tr><td>Number of spectro-z bins</td><td>n<sub>z</sub></td><td>4</td></tr><tr><td>Centres of the bins</td><td>z<sub>i</sub></td><td>{1.0, 1.2, 1.4, 1.65}</td></tr><tr><td>Error on redshift</td><td>&sigma;<sub>0,z</sub></td><td>0.001</td></tr><tr><td>Minimum scale</td><td>k<sub>min</sub></td><td>0.001 h Mpc<sup>-1</sup></td></tr><tr><td>Maximum scale</td><td>k<sub>max</sub></td><td></td></tr><tr><td>- Quasi-linear</td><td></td><td>0.15 h Mpc<sup>-1</sup></td></tr><tr><td>- Pessimistic</td><td></td><td>0.25 h Mpc<sup>-1</sup></td></tr><tr><td>- Optimistic</td><td></td><td>0.30 h Mpc<sup>-1</sup></td></tr></tbody></table>

conservative case since the underlying matter power spectrum of Eq. (36) that we used in our observed galaxy power spectrum recipe is a linear one. Non-linear corrections begin to play a role for wave numbers larger than $k = 0.1 h \text{ Mpc}^{-1}$ for the redshifts we considered (see Taruya et al. 2010). Cutting at a lower $k_{\text{max}}$ should avoid any bias induced by neglecting non-linear corrections below this scale. However, this severe cut also removes part of the information encoded in GC<sub>sp</sub> so that we wished to quantify how this affects the constraints. In all scenarios, we fixed the $\sigma_v = \sigma_p$ nuisance parameter for GC<sub>sp</sub> to the values directly calculated from Eq. (38) for the fiducial cosmological parameters.

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

# 5. Results

In this section, we discuss the results of the Fisher matrix analysis for the three models under investigation, JBD, nDGP, and

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

$$X^{(\pm)} = X_{\text{fid}} \times 10^{\pm \sigma_{\log_{10} X}}. \qquad (65)$$

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

In order to have a global view on the results for the three MG models analysed, we summarise in Fig. 2 the 68.3% marginal errors on the model parameters, $\log_{10} \omega_{\text{BD}}$ for JBD, $\log_{10} \Omega_{\text{rc}}$ for nDGP, and $\epsilon_{2,0}$ for KM in all the scenarios and combinations of probes we used. We note that we did not consider WL probe alone because we verified that WL provides constraints of about the same order as the GC<sub>sp</sub> ones. When we combined this with GC<sub>ph</sub> and XC, the constraints were significantly boosted. Additionally, we note that the constraining power of the spectroscopic sample is very weak because we did not go deep enough ($k \gtrsim 1 h \text{ Mpc}^{-1}$) in the non-linear regime and also because when we used GC probes, we introduced degeneracies among the MG parameters, amplitudes, and bias parameters. These features characterise all the models.

## 5.1. Jordan-Brans-Dicke gravity

There are remarkably strong non-cosmological constraints on JBD gravity, with current bounds of about $\omega_{\text{BD}} > 10^5$. As pointed out above, however, it is useful to consider JBD gravity as a long-wavelength limit of more general scalar-tensor theories with a wide range of small-scale limits that may be endowed with gravitational screening. It then makes sense to focus on what we already know from cosmology, where the constraints are about $\omega_{\text{BD}} > 10^3$, independent of what we know from non-cosmological systems. Two further practical reasons led us to choose $\omega_{\text{BD}} \sim 10^3$ as a fiducial value. First of all, this value is at the limits of the range for which we calibrated the non-linear prescription we used. Second, the likelihood is very flat
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Marginalised $1\sigma$ errors on the model parameters relative to their corresponding fiducial values, for JBD (case 1 in blue and case 2 in cyan), nDGP (case 1 in green and case 2 in light green), and KM (case 1 in red and case 2 in light red). We show the marginalised $1\sigma$ errors for $GC_{sp}$ in the optimistic (dotted area), pessimistic (dashed area) and quasi-linear (filled area) scenarios (top panel), for WL+XC+GC$_{ph}$ in the optimistic and pessimistic cases (central panel) and for the $GC_{sp}+WL+XC+GC_{ph}$ in the optimistic and pessimistic cases.</figcaption></figure>

for high values of $\omega_{BD}$ which can lead to numerical errors in the finite differences of the Fisher matrix evaluation for some of the observables. We thus considered two cases for JBD. The first fiducial case with $\omega_{BD} = 800$ (JBD1), and the second case with $\omega_{BD} = 2500$ (JBD2). Whereas the first case is nearly compatible with the most recent constraints from publicly available

CMB and LSS data, the second case is well within the current bounds (Ballardini et al. 2020; Joudaki et al. 2022; Ballardini et al. 2022). Table 2 summarises the relative (with respect to the fiducial values) uncertainties for the cosmological and model parameters in the quasi-linear, pessimistic and optimistic *Euclid* cases. Fig. 3 shows the 68.3% and 95.4% joint marginal error contours on the cosmological model parameters for the optimistic settings for both JBD1 (left panel) and JBD2 (right panel). As in Casas et al. (2023) for $f(R)$, we performed the Fisher matrix analysis on the parameter $\log_{10} \omega_{BD}$, and not directly on $\omega_{BD}$, because for large numbers and differences in the order of magnitude, the Fisher matrix derivatives might become unstable.

In the optimistic setting, we find that *Euclid* will be able to constrain $\log_{10} \omega_{BD}$ as follows:

- JBD1:

• at 27% for the $GC_{sp}$ alone;

• at 3.6% when considering the combination WL+XC+GC$_{ph}$;

• at 3.2% for $GC_{sp}+WL+XC+GC_{ph}$.

- JBD2:

• at 72.2% for the $GC_{sp}$ alone;

• at 5.1% when considering the combination WL+XC+GC$_{ph}$;

• at 4.6% for $GC_{sp}+WL+XC+GC_{ph}$.

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

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

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

## 5.2. Dvali-Gabadadze-Porrati gravity

The two cases considered for the nDGP model are $\log_{10} \Omega_{rc} = -0.60206$ (nDGP1) and $\log_{10} \Omega_{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. A.2 we show the same but for the pessimistic setting.

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

- nDGP1:

• at 93.4% for $GC_{sp}$ alone;

• at 20% when considering WL+XC+GC$_{ph}$;
<table><caption>**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 GC<sub>sp</sub>, WL+XC+GC<sub>ph</sub> and GC<sub>sp</sub>+WL+XC+GC<sub>ph</sub>.</caption><thead><tr><th></th><th>Ω<sub>m,0</sub></th><th>Ω<sub>b,0</sub></th><th>log<sub>10</sub> ω<sub>BD</sub></th><th>h</th><th>n<sub>s</sub></th><th>σ<sub>8</sub></th></tr></thead><tbody><tr><td colspan="7"><strong>JBD1 ω<sub>BD</sub> = 800</strong></td></tr><tr><td><strong>Pessimistic setting</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>GC<sub>sp</sub> (k<sub>max</sub> = 0.15 h Mpc<sup>-1</sup>) (quasi-linear)</td><td>2.13%</td><td>4.49%</td><td>43.04%</td><td>1.72%</td><td>2.62%</td><td>1.11%</td></tr><tr><td>GC<sub>sp</sub> (k<sub>max</sub> = 0.25 h Mpc<sup>-1</sup>)</td><td>1.55%</td><td>2.48%</td><td>29.62%</td><td>1.72%</td><td>1.53%</td><td>0.83%</td></tr><tr><td>WL+XC+GC<sub>ph</sub></td><td>0.37%</td><td>5.10%</td><td>10.89%</td><td>3.20%</td><td>1.25%</td><td>0.34%</td></tr><tr><td>GC<sub>sp</sub>+WL+XC+GC<sub>ph</sub></td><td>0.36%</td><td>1.72%</td><td>9.21%</td><td>0.91%</td><td>0.66%</td><td>0.28%</td></tr><tr><td><strong>Optimistic setting</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>GC<sub>sp</sub> (k<sub>max</sub> = 0.3 h Mpc<sup>-1</sup>)</td><td>1.49%</td><td>1.21%</td><td>27.13%</td><td>1.41%</td><td>1.33%</td><td>0.77%</td></tr><tr><td>WL+XC+GC<sub>ph</sub></td><td>0.27%</td><td>3.43%</td><td>3.55%</td><td>1.69%</td><td>0.52%</td><td>0.11%</td></tr><tr><td>GC<sub>sp</sub>+WL+XC+GC<sub>ph</sub></td><td>0.24%</td><td>1.43%</td><td>3.21%</td><td>0.51%</td><td>0.28%</td><td>0.10%</td></tr><tr><td colspan="7"><strong>JBD2 ω<sub>BD</sub> = 2500</strong></td></tr><tr><td><strong>Pessimistic setting</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>GC<sub>sp</sub> (k<sub>max</sub> = 0.15 h Mpc<sup>-1</sup>) (quasi-linear)</td><td>2.14%</td><td>4.52%</td><td>114.54%</td><td>3.80%</td><td>2.64%</td><td>1.11%</td></tr><tr><td>GC<sub>sp</sub> (k<sub>max</sub> = 0.25 h Mpc<sup>-1</sup>)</td><td>1.55%</td><td>2.49%</td><td>78.82%</td><td>1.74%</td><td>1.53%</td><td>0.83%</td></tr><tr><td>WL+XC+GC<sub>ph</sub></td><td>0.34%</td><td>5.09%</td><td>17.32%</td><td>3.21%</td><td>1.24%</td><td>0.27%</td></tr><tr><td>GC<sub>sp</sub>+WL+XC+GC<sub>ph</sub></td><td>0.33%</td><td>1.68%</td><td>15.51%</td><td>0.78%</td><td>0.54%</td><td>0.24%</td></tr><tr><td><strong>Optimistic setting</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>GC<sub>sp</sub> (k<sub>max</sub> = 0.3 h Mpc<sup>-1</sup>)</td><td>1.50%</td><td>2.24%</td><td>72.20%</td><td>1.42%</td><td>1.34%</td><td>0.77%</td></tr><tr><td>WL+XC+GC<sub>ph</sub></td><td>0.25%</td><td>3.40%</td><td>5.12%</td><td>1.70%</td><td>0.56%</td><td>0.10%</td></tr><tr><td>GC<sub>sp</sub>+WL+XC+GC<sub>ph</sub></td><td>0.22%</td><td>1.42%</td><td>4.56%</td><td>0.50%</td><td>0.28%</td><td>0.09%</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** 68.3 and 95.4% joint marginal error contours on the cosmological parameters for JBD1 (left panel) and JBD2 (right panel) in the optimistic case. In red, we plot GC<sub>sp</sub>, in blue, we plot WL+XC+GC<sub>ph</sub> and in green we plot GC<sub>sp</sub>+WL+XC+GC<sub>ph</sub>.</figcaption></figure>

* at 15% for GC<sub>sp</sub>+WL+XC+GC<sub>ph</sub>.

– nDGP2:

* GC<sub>sp</sub> has no power in constraining the parameter.

* at 81% for both WL+XC+GC<sub>ph</sub> and GC<sub>sp</sub>+WL+XC+GC<sub>ph</sub>.

Regardless of the specific case considered (nDGP1 or nDGP2), the power in constraining Ω<sub>m,0</sub> and Ω<sub>b,0</sub> is the same, while n<sub>s</sub> and σ<sub>8</sub> are better constrained in the nDGP2 case and h in the nDGP1 case. For example, in the optimistic setting, n<sub>s</sub> for

WL+XC+GC<sub>ph</sub> is constrained at 0.52% for nDGP1 and 0.37% for nDGP2 and they further improve when considering the full combination of probes, when they are 0.27 and 0.19%, respectively. In Fig. 4 we note in the log<sub>10</sub> Ω<sub>rc</sub>-n<sub>s</sub> panel that the two parameters are anti-correlated in the nDGP1 case. This makes the forecast errors on n<sub>s</sub> larger than the nDGP2 case where the anti-correlation disappears. A similar discussion can be made for the σ<sub>8</sub> parameter. In the optimistic setting, for nDGP1, we have 0.37% for WL+XC+GC<sub>ph</sub> and 0.21% for the full combination of
<table><caption>**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 GC<sub>sp</sub>, WL+XC+GC<sub>ph</sub> and GC<sub>sp</sub>+WL+XC+GC<sub>ph</sub>.</caption><thead><tr><th></th><th>Ω<sub>m,0</sub></th><th>Ω<sub>b,0</sub></th><th>log<sub>10</sub> Ω<sub>rc</sub></th><th>h</th><th>n<sub>s</sub></th><th>σ<sub>8</sub></th></tr></thead><tbody><tr><td colspan="7"><strong>nDGP1 Ω<sub>rc</sub> = 0.25</strong></td></tr><tr><td><strong>Pessimistic setting</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>GC<sub>sp</sub> (k<sub>max</sub> = 0.15 h Mpc<sup>-1</sup>) (quasi-linear)</td><td>2.14%</td><td>4.59%</td><td>145.65%</td><td>3.23%</td><td>1.82%</td><td>1.21%</td></tr><tr><td>GC<sub>sp</sub> (k<sub>max</sub> = 0.25 h Mpc<sup>-1</sup>)</td><td>1.54%</td><td>2.51%</td><td>101.27%</td><td>1.51%</td><td>0.99%</td><td>0.89%</td></tr><tr><td>WL+XC+GC<sub>ph</sub></td><td>1.13%</td><td>5.64%</td><td>49.59%</td><td>2.65%</td><td>0.92%</td><td>1.08%</td></tr><tr><td>GC<sub>sp</sub>+WL+XC+GC<sub>ph</sub></td><td>0.66%</td><td>2.18%</td><td>33.98%</td><td>1.26%</td><td>0.49%</td><td>0.42%</td></tr><tr><td><strong>Optimistic setting</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>GC<sub>sp</sub> (k<sub>max</sub> = 0.3 h Mpc<sup>-1</sup>)</td><td>1.46%</td><td>2.24%</td><td>93.39%</td><td>1.25%</td><td>0.87%</td><td>0.83%</td></tr><tr><td>WL+XC+GC<sub>ph</sub></td><td>0.30%</td><td>5.08%</td><td>19.59%</td><td>2.33%</td><td>0.52%</td><td>0.37%</td></tr><tr><td>GC<sub>sp</sub>+WL+XC+GC<sub>ph</sub></td><td>0.26%</td><td>1.62%</td><td>14.78%</td><td>0.59%</td><td>0.27%</td><td>0.21%</td></tr><tr><td colspan="7"><strong>nDGP2 Ω<sub>rc</sub> = 10<sup>-6</sup></strong></td></tr><tr><td><strong>Pessimistic setting</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>GC<sub>sp</sub> (k<sub>max</sub> = 0.15 h Mpc<sup>-1</sup>) (quasi-linear)</td><td>2.19%</td><td>5.01%</td><td>3053.42%</td><td>3.53%</td><td>1.94%</td><td>2.53%</td></tr><tr><td>GC<sub>sp</sub> (k<sub>max</sub> = 0.25 h Mpc<sup>-1</sup>)</td><td>1.56%</td><td>2.65%</td><td>1969.23%</td><td>1.65%</td><td>1.01%</td><td>1.68%</td></tr><tr><td>WL+XC+GC<sub>ph</sub></td><td>0.86%</td><td>5.61%</td><td>398.64%</td><td>2.65%</td><td>0.81%</td><td>0.58%</td></tr><tr><td>GC<sub>sp</sub>+WL+XC+GC<sub>ph</sub></td><td>0.77%</td><td>1.89%</td><td>379.52%</td><td>1.13%</td><td>0.77%</td><td>0.39%</td></tr><tr><td><strong>Optimistic setting</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>GC<sub>sp</sub> (k<sub>max</sub> = 0.3 h Mpc<sup>-1</sup>)</td><td>1.49%</td><td>2.24%</td><td>1842.43%</td><td>1.35%</td><td>0.89%</td><td>1.54%</td></tr><tr><td>WL+XC+GC<sub>ph</sub></td><td>0.29%</td><td>5.08%</td><td>81.02%</td><td>2.32%</td><td>0.37%</td><td>0.25%</td></tr><tr><td>GC<sub>sp</sub>+WL+XC+GC<sub>ph</sub></td><td>0.26%</td><td>1.39%</td><td>80.77%</td><td>0.56%</td><td>0.19%</td><td>0.14%</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** 68.3 and 95.4% joint marginal error contours on the cosmological parameters for nDGP1 (left panel) and nDGP2 (right panel) in the optimistic case. In red, we plot GC<sub>sp</sub>, in blue, we plot WL+XC+GC<sub>ph</sub> and in green we plot GC<sub>sp</sub>+WL+XC+GC<sub>ph</sub>.</figcaption></figure>

probes, while for the nDGP2 model we obtain 0.25 and 0.14%, respectively. Instead *h* is better constrained in nDGP1 being 1.25 and 1.35% for nDGP2 using GC<sub>sp</sub>.

As previously done for the JBD model, we translated the forecast errors on the parameter log<sub>10</sub> Ω<sub>rc</sub> into the parameter Ω<sub>rc</sub> [see Eq. (65)]. For the optimistic scenario, we obtain for nDGP1 Ω<sub>rc</sub> = 0.25<sup>+0.66</sup><sub>-0.18</sub> for GC<sub>sp</sub>, Ω<sub>rc</sub> = 0.25<sup>+0.08</sup><sub>-0.06</sub> for WL+XC+GC<sub>ph</sub> and Ω<sub>rc</sub> = 0.25<sup>+0.06</sup><sub>-0.05</sub> for GC<sub>sp</sub>+WL+XC+GC<sub>ph</sub>. These correspond to relative errors that can reach 264, 32, and 24%, respectively.

For nDGP2, we obtain an upper bound that is Ω<sub>rc</sub> < 0.07 in the optimistic setting, while for the pessimistic case, Ω<sub>rc</sub> is unconstrained regardless of the probe we consider because the fiducial parameter is very close to ΛCDM and we do not have enough non-linear information. To examine this, we considered the linear power spectrum: We find that the derivative with respect to log<sub>10</sub> Ω<sub>rc</sub> approaches 0 in the ΛCDM limit, and P'<sub>L</sub>|<sub>10<sup>-6</sup></sub> ≈ 0.0002 compared to P'<sub>L</sub>|<sub>0.25</sub> ≈ 0.08, where a prime here denotes a derivative with respect to log<sub>10</sub> Ω<sub>rc</sub>. This leads us to conclude that
<figure><img src="image_6.png" /><figcaption>Fig. 5. 68.3 and 95.4% joint marginal error contours on the cosmological parameters for KM1 (left panel) and KM2 (right panel) in the optimistic case. In red, we plot GC<sub>sp</sub>, in blue, we plot WL+XC+GC<sub>ph</sub> and in green, we plot GC<sub>sp</sub>+WL+XC+GC<sub>ph</sub>.</figcaption></figure>

the pessimistic case is still too dominated by linear structure
information to constrain Ω<sub>rc</sub> when we are close to the ΛCDM
limit. It should also be noted that Bose et al. (2020) reported
meaningful constraints using cosmic shear with the same scale
cuts we employed in our pessimistic case, but they used a full
Markov chain Monte Carlo (MCMC) approach to a ΛCDM fiduc-
ial data vector, sampling in Ω<sub>rc</sub>, as well as employing a different
covariance and binning scheme. Additionally, it is known that
the results obtained with a Fisher approach can be different from
those computed with the MCMC method even though the sam-
pled parameter is the same (Perotto et al. 2006; Wolz et al. 2012).
Therefore a direct comparison of the results is not straightfor-
ward. Lastly, our loss in sensitivity in the ΛCDM limit justifies
the tight upper limit in nDGP1 when compared to nDGP2.

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

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

5.3. k-mouflage

We present our results for the KM model with the two fiducial choices of the ε₂₀ parameter (KM1 ε₂₀ = -0.04 and KM2 ε₂₀ = -0.0001). The KM1 model represents an extreme case that is allowed within the 95.4% confidence interval by current cosmological data, while KM2 is practically indistinguishable

from ΛCDM. We show in Fig. 5 the forecast 68.3 and 95.4%
joint marginal error contours on cosmological and model param-
eters for the optimistic scenario. Fig. A.3 shows forecast results
for the same models in the pessimistic scenario. The 68.3%
errors for both the optimistic and pessimistic configurations are
reported in Table 4.

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

- KM1:

* at 3.4% for the GC<sub>sp</sub> alone;

* at 0.15% when considering the combination WL+XC+GC<sub>ph</sub>;

* at 0.14% for GC<sub>sp</sub>+WL+XC+GC<sub>ph</sub>.

- KM2:

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

* upper bound at 10⁻³ for the GC<sub>sp</sub> alone;

* upper bound at 7 × 10⁻⁴ when considering the combina-
tion WL+XC+GC<sub>ph</sub>;

* upper bound at 4 × 10⁻⁴ for GC<sub>sp</sub>+WL+XC+GC<sub>ph</sub>.

For the KM1 fiducial model, the optimistic configuration will
allow us to constrain ε₂₀ to ~ 3.4% accuracy using GC<sub>sp</sub> alone.
In the worst-case scenario of cutting the observed galaxy power
spectrum at k<sub>max</sub> = 0.15 h Mpc<sup>-1</sup>, the error on ε₂₀ increases to
~ 5.5%. The constraining power greatly improves for the combi-
nation WL+XC+GC<sub>ph</sub>, and the percentage error on ε₂₀ decreases
to 0.15% in the optimistic case, without a substantial worsening
in the pessimistic scenario. Of the cosmological parameters we
varied for the KM1 model, ε₂₀ is most tightly constrained by
the combination WL+XC+GC<sub>ph</sub>. This shows that this combina-
tion of data is well suited to capturing the modification induced
by k-mouflage to the Poisson and lensing equations, probing μ
and Σ independently. Moreover, the *Euclid* survey will probe the
redshift range z ≥ 1, where the largest effect determined by the
running of Newton's constant is expected to manifest (see, e.g.,
<table><caption>**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 GC<sub>sp</sub>, WL+XC+GC<sub>ph</sub> and GC<sub>sp</sub>+WL+XC+GC<sub>ph</sub>.</caption><thead><tr><th></th><th>Ω<sub>m,0</sub></th><th>Ω<sub>b,0</sub></th><th>ε<sub>2,0</sub></th><th>h</th><th>n<sub>s</sub></th><th>σ<sub>8</sub></th></tr></thead><tbody><tr><td colspan="7"><strong>KM1 ε<sub>2,0</sub> = -0.04</strong></td></tr><tr><td><strong>Pessimistic setting</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>GC<sub>sp</sub> (k<sub>max</sub> = 0.15 h Mpc<sup>-1</sup>) (quasi-linear)</td><td>6.53%</td><td>10.43%</td><td>5.48%</td><td>3.88%</td><td>1.95%</td><td>2.39%</td></tr><tr><td>GC<sub>sp</sub> (k<sub>max</sub> = 0.25 h Mpc<sup>-1</sup>)</td><td>2.23%</td><td>4.02%</td><td>3.75%</td><td>1.52%</td><td>1.10%</td><td>1.31%</td></tr><tr><td>WL+XC+GC<sub>ph</sub></td><td>0.36%</td><td>7.36%</td><td>0.22%</td><td>1.40%</td><td>1.93%</td><td>1.10%</td></tr><tr><td>GC<sub>sp</sub>+WL+XC+GC<sub>ph</sub></td><td>0.31%</td><td>1.77%</td><td>0.22%</td><td>0.44%</td><td>0.63%</td><td>0.46%</td></tr><tr><td><strong>Optimistic setting</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>GC<sub>sp</sub> (k<sub>max</sub> = 0.3 h Mpc<sup>-1</sup>)</td><td>1.82%</td><td>3.50%</td><td>3.40%</td><td>1.23%</td><td>0.89%</td><td>1.19%</td></tr><tr><td>WL+XC+GC<sub>ph</sub></td><td>0.20%</td><td>4.30%</td><td>0.15%</td><td>0.59%</td><td>0.98%</td><td>0.84%</td></tr><tr><td>GC<sub>sp</sub>+WL+XC+GC<sub>ph</sub></td><td>0.17%</td><td>1.57%</td><td>0.14%</td><td>0.32%</td><td>0.43%</td><td>0.37%</td></tr><tr><td colspan="7"><strong>KM2 ε<sub>2,0</sub> = -0.0001</strong></td></tr><tr><td><strong>Pessimistic setting</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>GC<sub>sp</sub> (k<sub>max</sub> = 0.15 h Mpc<sup>-1</sup>) (quasi-linear)</td><td>7.85%</td><td>12.28%</td><td>1720.00%</td><td>4.48%</td><td>2.24%</td><td>2.84%</td></tr><tr><td>GC<sub>sp</sub> (k<sub>max</sub> = 0.25 h Mpc<sup>-1</sup>)</td><td>3.19%</td><td>5.41%</td><td>1190.00%</td><td>1.99%</td><td>1.13%</td><td>1.47%</td></tr><tr><td>WL+XC+GC<sub>ph</sub></td><td>1.18%</td><td>5.22%</td><td>939.23%</td><td>1.03%</td><td>1.51%</td><td>1.17%</td></tr><tr><td>GC<sub>sp</sub>+WL+XC+GC<sub>ph</sub></td><td>0.73%</td><td>2.23%</td><td>642.65%</td><td>0.79%</td><td>0.82%</td><td>0.80%</td></tr><tr><td><strong>Optimistic setting</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>GC<sub>sp</sub> (k<sub>max</sub> = 0.3 h Mpc<sup>-1</sup>)</td><td>2.57%</td><td>4.41%</td><td>1120.00%</td><td>0.62%</td><td>1.00%</td><td>1.34%</td></tr><tr><td>WL+XC+GC<sub>ph</sub></td><td>0.40%</td><td>3.43%</td><td>658.67%</td><td>0.66%</td><td>0.87%</td><td>0.75%</td></tr><tr><td>GC<sub>sp</sub>+WL+XC+GC<sub>ph</sub></td><td>0.33%</td><td>1.73%</td><td>396.60%</td><td>0.43%</td><td>0.55%</td><td>0.42%</td></tr></tbody></table>

Benevento et al. 2019). The fact that with the photometric probes we obtain the tightest relative constraints for the KM1 model, compared to the other two models studied in this work, is consistent with the observation of the lower panels of Figure 1, which show that the lensing angular power spectrum is more distant to its ΛCDM counterpart than in the other two models. We note that the addition of GC<sub>sp</sub> to WL+XC+GC<sub>ph</sub> does not improve the constraints on ε<sub>2,0</sub> in either the pessimistic and optimistic scenarios.

The percentage error on ε<sub>2,0</sub> in the ΛCDM-proximate KM2 fiducial model decreases from 1120% using GC<sub>sp</sub> alone to a value of 397% with the full combination of GC<sub>sp</sub>+WL+XC+GC<sub>ph</sub> in the best-case scenario. In the optimistic case, the *Euclid* survey will therefore be able to detect a deviation from ΛCDM only when |ε<sub>2,0</sub>| ≥ 4 × 10<sup>-4</sup>, which improves over the constraining power of present CMB and LSS data by somewhat more than one order of magnitude.

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

## 6. Conclusions

We have explored the constraining power of the future *Euclid* mission for linearly scale-independent extensions of the concordance cosmological model, that is, models that induce modifications to the linear growth of perturbations that are solely time-dependent while featuring different testable types of screening mechanisms at smaller non-linear scales. We considered three

specific models, namely JBD, a scalar-tensor theory with a standard kinetic term and a flat potential (Sect. 2.1); the nDGP gravity, a braneworld model in which our Universe is a four-dimensional brane embedded in a five-dimensional Minkowski space-time (Sect. 2.2); and KM gravity, an extension of k-essence scenarios with a universal coupling of the scalar field to matter (Sect. 2.3).

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

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

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

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

For the KM scenario, we obtain constraints for the full combination of probes in the fiducial KM1 case of $\epsilon_{2,0} = -0.04 \pm 5.6 \times 10^{-5}$. Out of the full set of cosmological plus model parameters, $\epsilon_{2,0}$ is the most tightly constrained parameter for KM1 by the full combination because of the independent constraints on clustering and lensing, which probe the running of Newton's constant in the range $z \gtrsim 1$ where the strongest 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 ΛCDM limit, the nDGP2 and KM2 models remain unconstrained in the pessimistic setting, it is remarkable that *Euclid* alone can significantly reduce the allowed space for these models. For the full probe combination (GC<sub>sp</sub>+WL+XC+GC<sub>ph</sub>) the nDGP2 upper bound at $\Omega_{rc} < 0.07$ for the optimistic scenario significantly improves the current constraints derived with the same approach and similar type of data. This is consistent with the results found in Bose et al. (2020).

The forecasts for KM2 give $|\epsilon_{2,0}| < 4 \times 10^{-4}$ in the optimistic case and $|\epsilon_{2,0}| < 2 \times 10^{-3}$ in the pessimistic case. This improves the current constraints by more than one order of magnitude, even in the pessimistic case, by combining different cosmological datasets (Benevento et al. 2019). In addition, *Euclid* alone will provide constraints on $\omega_{BD}$ that will be tighter than those obtained by a combination of current CMB and low-redshift data Joudaki et al. (2022).

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

(Euclid Collaboration 2022), the improvements in the constraints when *Euclid* is cross-correlated with CMB data can be of the order of two to three and in some cases even larger. A large effort is currently ongoing to extend this analysis to other extended models, including those in this work.

*Acknowledgements.* The Euclid Consortium acknowledges the European Space Agency and a number of agencies and institutes that have supported the development of *Euclid*, in particular the Academy of Finland, the Agenzia Spaziale Italiana, the Belgian Science Policy, the Canadian Euclid Consortium, the French Centre National d'Etudes Spatiales, the Deutsches Zentrum für Luft- und Raumfahrt, the Danish Space Research Institute, the Fundação para a Ciência e a Tecnologia, the Ministerio de Ciencia e Innovación, the National Aeronautics and Space Administration, the National Astronomical Observatory of Japan, the Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the Romanian Space Agency, the State Secretariat for Education, Research and Innovation (SERI) at the Swiss Space Office (SSO), and the United Kingdom Space Agency. A complete and detailed list is available on the *Euclid* web site (http://www.euclid-ec.org). B.B. and L.L. were supported by a Swiss National Science Foundation Professorship grant (Nos. 170547 & 202671). E.B. has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie 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. G.B. acknowledges support from ASI/INFN grant n. 202one-43-HH.0. M.P. acknowledges support by the MIUR 'Progetti di Ricerca di Rilevante Interesse Nazionale' (PRIN) Bando 2022 - grant 20228RMX4A.

## References

Aghanim, N., Akrami, Y., Ashdow, Y., et al. 2020, A&A, 641, A6; Erratum: 2021, A&A, 652, C4

Avilez, A., & Skordis, C. 2014, Phys. Rev. Lett., 113, 011101

Babichev, E., Deffayet, C., & Ziour, R. 2009, Int. J. Mod. Phys. D, 18, 2147

Bag, S., Mishra, S. S., & Sahni, V. 2018, Phys. Rev. D, 97, 123537

Ballardini, M., Finelli, F., Umiltà, C., & Paoletti, D. 2016, JCAP, 05, 067

Ballardini, M., Sapone, D., Umiltà, C., et al. 2019, JCAP, 05, 049

Ballardini, M., Braglia, M., Finelli, F., et al. 2020, JCAP, 10, 044

Ballardini, M., Finelli, F., & Sapone, D. 2022, JCAP, 06, 004

Bardeen, J. M. 1980, Phys. Rev. D, 22, 1882

Barreira, A., Brax, P., Clesse, S., et al. 2015, Phys. Rev. D, 91, 063528

Barreira, A., Sánchez, A. G., & Schmidt, F. 2016, Phys. Rev. D, 94, 084022

Bellini, E., & Sawicki, I. 2014, JCAP, 07, 050

Bellini, E., Barreira, A., Frusciante, N., et al. 2018, Phys. Rev. D, 97, 023520

Bellini, E., Sawicki, I., & Zumalacárregui, M. 2020, JCAP, 02, 008

Benevento, G., Raveri, M., Lazanu, A., et al. 2019, JCAP, 05, 027

Berti, E., Barausse, E., Cardoso, V., et al. 2015, Class. Quant. Grav., 32, 243001

Bertotti, B., Iess, L., & Tortora, P. 2003, Nature, 425, 374

Blas, D., Lesgourgues, J., & Tram, T. 2011, JCAP, 07, 034

Boisseau, B., Esposito-Farese, G., Polarski, D., & Starobinsky, A. A. 2000, Phys. Rev. Lett., 85, 2236

Bose, B., Cataneo, M., Tröster, T., et al. 2020, MNRAS, 498, 4650

Bose, B., Wright, B. S., Cataneo, M., et al. 2021, MNRAS, 508, 2479

Bowcock, P., Charmousis, C., & Gregory, R. 2000, Class. Quant. Grav., 17, 4745

Braglia, M., Ballardini, M., Emond, W. T., et al. 2020, Phys. Rev. D, 102, 023529

Braglia, M., Ballardini, M., Finelli, F., & Koyama, K. 2021, Phys. Rev. D, 103, 043528
Brans, C., & Dicke, R. H. 1961, Phys. Rev., 124, 925

Brax, P., & Valageas, P. 2013, Phys. Rev. D, 88, 023527

Brax, P., & Valageas, P. 2014, Phys. Rev. D, 90, 023508

Brax, P., & Valageas, P. 2016, JCAP, 01, 020

Brax, P., Rizzo, L. A., & Valageas, P. 2015, Phys. Rev. D, 92, 043519

Brax, P., Casas, S., Desmond, H., & Elder, B. 2021, Universe, 8, 11

Bullock, J. S., Kolatt, T. S., Sigad, Y., et al. 2001, MNRAS, 321, 559

Cardoso, A., Koyama, K., Seahra, S. S., & Silva, F. P. 2008, Phys. Rev. D, 77, 083512

Carroll, S. M., Duvvuri, V., Trodden, M., & Turner, M. S. 2004, Phys. Rev. D, 70, 043528

Casas, S. J., Cardone, V. F., Sapone, D., et al. 2023, A&A, submitted [arXiv:2306.11053]

Cataneo, M., Lombriser, L., Heymans, C., et al. 2019, MNRAS, 488, 2121

Cataneo, M., Uhlemann, C., Arnold, C., et al. 2022, MNRAS, 513, 1623

Charmousis, C., Gregory, R., Kaloper, N., & Padilla, A. 2006, J. High Energy Phys., 10, 066

Chen, S.-F., Vlah, Z., & White, M. 2022, JCAP, 02, 008

Cheng, G., Wu, F., & Chen, X. 2021, Phys. Rev. D, 103, 103527

Clifton, T., Ferreira, P. G., Padilla, A., & Skordis, C. 2012, Phys. Rept., 513, 1

Cooray, A., & Sheth, R. K. 2002, Phys. Rev., 372, 1

D'Amico, G., Gleyzes, J., Kokron, N., et al. 2020, JCAP, 05, 005

Deffayet, C. 2001, Phys. Lett. B, 502, 199

Dvali, G. R., Gabadadze, G., & Porrati, M. 2000, Phys. Lett. B, 485, 208

Eisenstein, D. J., Seo, H.-j., & White, M. J. 2007, ApJ, 664, 660

Euclid Collaboration (Blanchard, A., et al.) 2020, A&A, 642, A191

Euclid Collaboration (Ilić, S., et al.) 2022, A&A, 657, A91

Ferreira, P. G. 2019, ARA&A, 57, 335

Gleyzes, J., Langlois, D., & Vernizzi, F. 2015, Int. J. Mod. Phys. D, 23, 1443010

Gorbunov, D., Koyama, K., & Sibiryakov, S. 2006, Phys. Rev. D, 73, 044016

Hernández-Aguayo, C., Ruan, C.-Z., Li, B., et al. 2022, JCAP, 01, 048

Hojjati, A., Pogosian, L., & Zhao, G.-B. 2011, JCAP, 1108, 005

Hu, B., Raveri, M., Frusciante, N., & Silvestri, A. 2014, Phys. Rev. D, 89, 103530

Hu, W., & Sawicki, I. 2007a, Phys. Rev. D, 76, 104043

Hu, W., & Sawicki, I. 2007b, Phys. Rev. D, 76, 064004

Ivanov, M. M., Simonović, M., & Zaldarriaga, M. 2020, JCAP, 05, 042

Joudaki, S., Blake, C., Johnson, A., et al. 2018, MNRAS, 474, 4894

Joudaki, S., Ferreira, P. G., Lima, N. A., & Winther, H. A. 2022, Phys. Rev. D, 105, 043522

Khoury, J., & Weltman, A. 2004, Phys. Rev. Lett, 93, 171104

Koyama, K., & Maartens, R. 2006, JCAP, 01, 016

Lee, S., Huff, E. M., Choi, A., et al. 2021, MNRAS, 509, 4982

Lesgourgues, J. 2011, arXiv e-prints [arXiv:1104.2932]

Lewis, A., Challinor, A., & Lasenby, A. 2000, ApJ, 538, 473

Liu, R., Valogiannis, G., Battaglia, N., & Bean, R. 2021, Phys. Rev. D, 104, 103519

Lombriser, L., Hu, W., Fang, W., & Seljak, U. 2009, Phys. Rev. D, 80, 063536

Lombriser, L., Yoo, J., & Koyama, K. 2013, Phys. Rev. D, 87, 104019

Luty, M. A., Porrati, M., & Rattazzi, R. 2003, J. High Energy Phys., 09, 029

Ma, C.-P., & Bertschinger, E. 1995, ApJ, 455, 7

Mead, A., Peacock, J., Heymans, C., Joudaki, S., & Heavens, A. 2015, MNRAS, 454, 1958

Mead, A., Heymans, C., Lombriser, L., et al. 2016, MNRAS, 459, 1468

Mead, A., Brieden, S., Tröster, T., & Heymans, C. 2021, MNRAS, 502, 1401

Müller, E.-M., Percival, W., Linder, E., et al. 2018, MNRAS, 475, 2122

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493

Nguyen, N.-M., Huterer, D., & Wen, Y. 2023, Phys. Rev. Lett., 131, 111001

Ooba, J., Ichiki, K., Chiba, T., & Sugiyama, N. 2017, Prog. Theor. Exp. Phys., 2017, 043E03

Pace, F., Battye, R., Bellini, E., et al. 2021, JCAP, 06, 017

Perotto, L., Lesgourgues, J., Hannestad, S., Tu, H., & Wong, Y. Y. Y. 2006, JCAP, 10, 013

Piga, L., Marinucci, M., D'Amico, G., et al. 2023, JCAP, 04, 038

Raccanelli, A., Bertacca, D., Pietrobon, D., et al. 2013, MNRAS, 436, 89

Raveri, M., Hu, B., Frusciante, N., & Silvestri, A. 2014, Phys. Rev. D, 90, 043513

Raveri, M., Pogosian, L., Martinelli, M., et al. 2021, arXiv e-prints [arXiv:2107.12990]

Rossi, M., Ballardini, M., Braglia, M., et al. 2019, Phys. Rev. D, 100, 103524

Schmidt, F. 2009, Phys. Rev. D, 80, 123003

Schneider, A., Refregier, A., Grandis, S., et al. 2020a, JCAP, 04, 020

Schneider, A., Stoira, N., Refregier, A., et al. 2020b, JCAP, 04, 019

Seahra, S. S., & Hu, W. 2010, Phys. Rev. D, 82, 124015

Sheth, R. K., & Tormen, G. 1999, MNRAS, 308, 119

Sheth, R. K., & Tormen, G. 2002, MNRAS, 329, 61

Silvestri, A., Pogosian, L., & Buniy, R. V. 2013, Phys. Rev. D, 87, 104015

Sola, J., Gomez-Valent, A., de Cruz Perez, J., et al. 2019, Astrophys. J. Lett. 1, L6

Sola, J., Gomez-Valent, A., de Cruz Perez, J., et al. 2020, Class. Quant. Grav. 24, 245003

Song, Y.-S. 2008, Phys. Rev. D, 77, 124031

Spurio Mancini, A., & Bose, B. 2023, arXiv e-prints [arXiv:2305.06350]

Spurio Mancini, A., Köhlinger, F., Joachimi, B., et al. 2019, MNRAS, 490, 2155

Takahashi, R., Sato, M., Nishimichi, T., Taruya, A., & Oguri, M. 2012, ApJ, 761, 152

Taruya, A., Nishimichi, T., & Saito, S. 2010, Phys. Rev. D, 82, 063522

Tassev, S., Zaldarriaga, M., & Eisenstein, D. 2013, JCAP, 06, 036

Teysier, R. 2002, A&A, 385, 337

Touboul, P., Métris, G., Rodrigues, M., et al. 2017, Phys. Rev. Lett., 119, 231101

Umiltà, C., Ballardini, M., Finelli, F., & Paoletti, D. 2015, JCAP, 08, 017

Vainshtein, A. I. 1972, Phys. Lett. B, 39, 393

Valageas, P., Nishimichi, T., & Taruya, A. 2013, Phys. Rev. D, 87, 083522

Voisin, G., Cognard, I., Freire, P. C. C., et al. 2020, A&A, 638, A24

Williams, J. G., Turyshev, S. G., & Boggs, D. 2012, Class. Quant. Grav., 29, 184004

Winther, H. A., Koyama, K., Manera, M., Wright, B. S., & Zhao, G.-B. 2017, JCAP, 08, 006

Wolz, L., Kilbinger, M., Weller, J., & Giannantonio, T. 2012, JCAP, 09, 009

Xu, L. 2014, JCAP, 02, 048

Zhao, G.-B., Pogosian, L., Silvestri, A., & Zylberberg, J. 2009, Phys. Rev. D, 79, 083513

Zucca, A., Pogosian, L., Silvestri, A., & Zhao, G.-B. 2019, JCAP, 05, 001

Zumalacárregui, M., Bellini, E., Sawicki, I., et al. 2017, JCAP, 08, 019

<footnote n="1">¹ Department of Physics "E. Pancini", University Federico II, Via Cinthia 6, 80126 Napoli, Italy</footnote>

<footnote n="2">² Dipartimento di Fisica, Università degli Studi di Torino, Via P. Giuria 1, 10125 Torino, Italy</footnote>

<footnote n="3">³ INFN-Sezione di Torino, Via P. Giuria 1, 10125 Torino, Italy</footnote>

<footnote n="4">⁴ INAF-Osservatorio Astrofisico di Torino, Via Osservatorio 20, 10025 Pino Torinese (TO), Italy</footnote>

<footnote n="5">⁵ INAF-Osservatorio Astronomico di Roma, Via Frascati 33, 00078 Monteporzio Catone, Italy</footnote>

<footnote n="6">⁶ INFN-Sezione di Roma, Piazzale Aldo Moro, 2 - c/o Dipartimento di Fisica, Edificio G. Marconi, 00185 Roma, Italy</footnote>

<footnote n="7">⁷ Institute for Theoretical Particle Physics and Cosmology (TTK), RWTH Aachen University, 52056 Aachen, Germany</footnote>

<footnote n="8">⁸ Université de Genève, Département de Physique Théorique and Centre for Astroparticle Physics, 24 quai Ernest-Ansermet, 1211 Genève 4, Switzerland</footnote>

<footnote n="9">⁹ Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s/n, 08193 Barcelona, Spain</footnote>

<footnote n="10">¹⁰ Institut d'Estudis Espacials de Catalunya (IEEC), Carrer Gran Capità 2-4, 08034 Barcelona, Spain</footnote>

<footnote n="11">¹¹ Institut de Recherche en Astrophysique et Planétologie (IRAP), Université de Toulouse, CNRS, UPS, CNES, 14 Av. Edouard Belin, 31400 Toulouse, France</footnote>

<footnote n="12">¹² Dipartimento di Fisica e Scienze della Terra, Università degli Studi di Ferrara, Via Giuseppe Saragat 1, 44122 Ferrara, Italy</footnote>

<footnote n="13">¹³ Dipartimento di Fisica e Astronomia, Università di Bologna, Via Gobetti 93/2, 40129 Bologna, Italy</footnote>

<footnote n="14">¹⁴ INAF-Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Piero Gobetti 93/3, 40129 Bologna, Italy</footnote>

<footnote n="15">¹⁵ INFN-Sezione di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy</footnote>

<footnote n="16">¹⁶ INFN, Sezione di Trieste, Via Valerio 2, 34127 Trieste TS, Italy</footnote>

<footnote n="17">¹⁷ IFPU, Institute for Fundamental Physics of the Universe, via Beirut 2, 34151 Trieste, Italy</footnote>

<footnote n="18">¹⁸ SISSA, International School for Advanced Studies, Via Bonomea 265, 34136 Trieste TS, Italy</footnote>

<footnote n="19">¹⁹ Johns Hopkins University 3400 North Charles Street Baltimore, MD 21218, USA</footnote>

<footnote n="20">²⁰ Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="21">²¹ Institute for Computational Science, University of Zurich, Winterthurerstrasse 190, 8057 Zurich, Switzerland</footnote>

<footnote n="22">²² Institut de Physique Théorique, CEA, CNRS, Université Paris-Saclay 91191 Gif-sur-Yvette Cedex, France</footnote>

<footnote n="23">²³ Dipartimento di Fisica e Astronomia “G.Galilei”, Università di Padova, Via Marzolo 8, 35131 Padova, Italy</footnote>
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, PO 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 Instituto 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 Astrofísico 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, 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
## Appendix A: Further results

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

<figure><img src="image_7.png" /><figcaption>Fig. A.1. 68.3% and 95.4% joint marginal error contours on the cosmological parameters for JBD1 (left panel) and JBD2 (right panel) in the pessimistic case. In red, we plot GC<sub>sp</sub>, in blue, we plot WL+XC+GC<sub>ph</sub> and in green, we plot GC<sub>sp</sub>+WL+XC+GC<sub>ph</sub>.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. A.2. 68.3% and 95.4% joint marginal error contours on the cosmological parameters for nDGP1 (left panel) and nDGP2 (right panel) in the pessimistic case. In red, we plot GC<sub>sp</sub>, in blue, we plot WL+XC+GC<sub>ph</sub> and in green, we plot GC<sub>sp</sub>+WL+XC+GC<sub>ph</sub>.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. A.3.** 68.3% and 95.4% joint marginal error contours on the cosmological parameters for KM1 (left panel) and KM2 (right panel) in the pessimistic case. In red, we plot GC<sub>sp</sub>, in blue, we plot WL+XC+GC<sub>ph</sub> and in green, we plot GC<sub>sp</sub>+WL+XC+GC<sub>ph</sub>.</figcaption></figure>