11institutetext: 1 Institute for Theoretical Particle Physics and Cosmology
(TTK), RWTH Aachen University, 52056 Aachen, Germany
2 INAF-Osservatorio Astronomico di Roma, Via Frascati 33, 00078 Monteporzio
Catone, Italy
3 INFN-Sezione di Roma, Piazzale Aldo Moro, 2 - c/o Dipartimento di Fisica,
Edificio G. Marconi, 00185 Roma, Italy
4 Departamento de Física, FCFM, Universidad de Chile, Blanco Encalada 2008,
Santiago, Chile
5 Department of Physics ”E. Pancini”, University Federico II, Via Cinthia 6,
80126, Napoli, Italy
6 Dipartimento di Fisica, Universitá degli Studi di Torino, Via P. Giuria 1,
10125 Torino, Italy
7 INFN-Sezione di Torino, Via P. Giuria 1, 10125 Torino, Italy
8 INAF-Osservatorio Astrofisico di Torino, Via Osservatorio 20, 10025 Pino
Torinese (TO), Italy
9 Dipartimento di Fisica, Universitá degli studi di Genova, and INFN-Sezione
di Genova, via Dodecaneso 33, 16146, Genova, Italy
10 INAF-Osservatorio Astronomico di Trieste, Via G. B. Tiepolo 11, 34143
Trieste, Italy
11 SISSA, International School for Advanced Studies, Via Bonomea 265, 34136
Trieste TS, Italy
12 IFPU, Institute for Fundamental Physics of the Universe, via Beirut 2,
34151 Trieste, Italy
13 Dipartimento di Fisica ”Aldo Pontremoli”, Universitá degli Studi di Milano,
Via Celoria 16, 20133 Milano, Italy
14 INFN-Sezione di Milano, Via Celoria 16, 20133 Milano, Italy
15 Institute of Cosmology and Gravitation, University of Portsmouth,
Portsmouth PO1 3FX, UK
16 Institut de Recherche en Astrophysique et Planétologie (IRAP), Université
de Toulouse, CNRS, UPS, CNES, 14 Av. Edouard Belin, 31400 Toulouse, France
17 Université de Genève, Département de Physique Théorique and Centre for
Astroparticle Physics, 24 quai Ernest-Ansermet, CH-1211 Genève 4, Switzerland
18 Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans,
s/n, 08193 Barcelona, Spain
19 Institut d’Estudis Espacials de Catalunya (IEEC), Carrer Gran Capitá 2-4,
08034 Barcelona, Spain
20 Université Paris-Saclay, Université Paris Cité, CEA, CNRS, Astrophysique,
Instrumentation et Modélisation Paris-Saclay, 91191 Gif-sur-Yvette, France
21 Institut für Theoretische Physik, University of Heidelberg, Philosophenweg
16, 69120 Heidelberg, Germany
22 Université St Joseph; Faculty of Sciences, Beirut, Lebanon
23 Institute Lorentz, Leiden University, PO Box 9506, Leiden 2300 RA, The
Netherlands
24 Dipartimento di Scienze Matematiche, Fisiche e Informatiche, Università di
Parma, Viale delle Scienze 7/A 43124 Parma, Italy
25 INFN Gruppo Collegato di Parma, Viale delle Scienze 7/A 43124 Parma, Italy
26 Institut de Physique Théorique, CEA, CNRS, Université Paris-Saclay 91191
Gif-sur-Yvette Cedex, France
27 Mullard Space Science Laboratory, University College London, Holmbury St
Mary, Dorking, Surrey RH5 6NT, UK
28 Institute for Astronomy, University of Edinburgh, Royal Observatory,
Blackford Hill, Edinburgh EH9 3HJ, UK
29 Higgs Centre for Theoretical Physics, School of Physics and Astronomy, The
University of Edinburgh, Edinburgh EH9 3FD, UK
30 Université Paris-Saclay, CNRS, Institut d’astrophysique spatiale, 91405,
Orsay, France
31 INAF-IASF Milano, Via Alfonso Corti 12, 20133 Milano, Italy
32 Instituto de Física Teórica UAM-CSIC, Campus de Cantoblanco, 28049 Madrid,
Spain
33 ESAC/ESA, Camino Bajo del Castillo, s/n., Urb. Villafranca del Castillo,
28692 Villanueva de la Cañada, Madrid, Spain
34 INAF-Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via
Piero Gobetti 93/3, 40129 Bologna, Italy
35 Dipartimento di Fisica e Astronomia, Universitá di Bologna, Via Gobetti
93/2, 40129 Bologna, Italy
36 INFN-Sezione di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy
37 Max Planck Institute for Extraterrestrial Physics, Giessenbachstr. 1, 85748
Garching, Germany
38 Dipartimento di Fisica, Universitá di Genova, Via Dodecaneso 33, 16146,
Genova, Italy
39 INFN-Sezione di Genova, Via Dodecaneso 33, 16146, Genova, Italy
40 INAF-Osservatorio Astronomico di Capodimonte, Via Moiariello 16, 80131
Napoli, Italy
41 Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP,
Rua das Estrelas, PT4150-762 Porto, Portugal
42 Institut de Física d’Altes Energies (IFAE), The Barcelona Institute of
Science and Technology, Campus UAB, 08193 Bellaterra (Barcelona), Spain
43 Port d’Informació Científica, Campus UAB, C. Albareda s/n, 08193 Bellaterra
(Barcelona), Spain
44 INFN section of Naples, Via Cinthia 6, 80126, Napoli, Italy
45 Dipartimento di Fisica e Astronomia ”Augusto Righi” - Alma Mater Studiorum
Universitá di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy
46 Centre National d’Etudes Spatiales – Centre spatial de Toulouse, 18 avenue
Edouard Belin, 31401 Toulouse Cedex 9, France
47 Institut national de physique nucléaire et de physique des particules, 3
rue Michel-Ange, 75794 Paris Cédex 16, France
48 Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy,
University of Manchester, Oxford Road, Manchester M13 9PL, UK
49 European Space Agency/ESRIN, Largo Galileo Galilei 1, 00044 Frascati, Roma,
Italy
50 University of Lyon, Univ Claude Bernard Lyon 1, CNRS/IN2P3, IP2I Lyon, UMR
5822, 69622 Villeurbanne, France
51 Institute of Physics, Laboratory of Astrophysics, Ecole Polytechnique
Fédérale de Lausanne (EPFL), Observatoire de Sauverny, 1290 Versoix,
Switzerland
52 UCB Lyon 1, CNRS/IN2P3, IUF, IP2I Lyon, 4 rue Enrico Fermi, 69622
Villeurbanne, France
53 Departamento de Física, Faculdade de Ciências, Universidade de Lisboa,
Edifício C8, Campo Grande, PT1749-016 Lisboa, Portugal
54 Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências,
Universidade de Lisboa, Campo Grande, 1749-016 Lisboa, Portugal
55 Department of Astronomy, University of Geneva, ch. d’Ecogia 16, 1290
Versoix, Switzerland
56 Department of Physics, Oxford University, Keble Road, Oxford OX1 3RH, UK
57 INFN-Padova, Via Marzolo 8, 35131 Padova, Italy
58 Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191, Gif-
sur-Yvette, France
59 Aix-Marseille Université, CNRS/IN2P3, CPPM, Marseille, France
60 INAF-Osservatorio Astronomico di Padova, Via dell’Osservatorio 5, 35122
Padova, Italy
61 University Observatory, Faculty of Physics, Ludwig-Maximilians-Universität,
Scheinerstr. 1, 81679 Munich, Germany
62 INAF-Osservatorio Astronomico di Brera, Via Brera 28, 20122 Milano, Italy
63 Institute of Theoretical Astrophysics, University of Oslo, P.O. Box 1029
Blindern, 0315 Oslo, Norway
64 von Hoerner & Sulger GmbH, SchloßPlatz 8, 68723 Schwetzingen, Germany
65 Technical University of Denmark, Elektrovej 327, 2800 Kgs. Lyngby, Denmark
66 Cosmic Dawn Center (DAWN), Denmark
67 Institut d’Astrophysique de Paris, UMR 7095, CNRS, and Sorbonne Université,
98 bis boulevard Arago, 75014 Paris, France
68 Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg,
Germany
69 Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak
Grove Drive, Pasadena, CA, 91109, USA
70 Department of Physics, P.O. Box 64, 00014 University of Helsinki, Finland
71 Helsinki Institute of Physics, Gustaf Hällströmin katu 2, University of
Helsinki, 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 CEA Saclay, DFR/IRFU, Service d’Astrophysique, Bat. 709, 91191 Gif-sur-
Yvette, France
77 European Space Agency/ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The
Netherlands
78 Department of Physics and Astronomy, University of Aarhus, Ny Munkegade
120, DK-8000 Aarhus C, Denmark
79 Centre for Astrophysics, University of Waterloo, Waterloo, Ontario N2L 3G1,
Canada
80 Department of Physics and Astronomy, University of Waterloo, Waterloo,
Ontario N2L 3G1, Canada
81 Perimeter Institute for Theoretical Physics, Waterloo, Ontario N2L 2Y5,
Canada
82 Space Science Data Center, Italian Space Agency, via del Politecnico snc,
00133 Roma, Italy
83 Institute of Space Science, Str. Atomistilor, nr. 409 Măgurele, Ilfov,
077125, Romania
84 Dipartimento di Fisica e Astronomia ”G. Galilei”, Universitá di Padova, Via
Marzolo 8, 35131 Padova, Italy
85 Universitäts-Sternwarte München, Fakultät für Physik, Ludwig-Maximilians-
Universität München, Scheinerstrasse 1, 81679 München, Germany
86 Institut d’Astrophysique de Paris, 98bis Boulevard Arago, 75014, Paris,
France
87 Junia, EPA department, 41 Bd Vauban, 59800 Lille, France
88 Satlantis, University Science Park, Sede Bld 48940, Leioa-Bilbao, Spain
89 AIM, CEA, CNRS, Université Paris-Saclay, Université de Paris, 91191 Gif-
sur-Yvette, France
90 Aix-Marseille Université, CNRS, CNES, LAM, Marseille, France
91 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas
(CIEMAT), Avenida Complutense 40, 28040 Madrid, Spain
92 Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências,
Universidade de Lisboa, Tapada da Ajuda, 1349-018 Lisboa, Portugal
93 Universidad Politécnica de Cartagena, Departamento de Electrónica y
Tecnología de Computadoras, Plaza del Hospital 1, 30202 Cartagena, Spain
94 Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700
AV Groningen, The Netherlands
95 INFN-Bologna, Via Irnerio 46, 40126 Bologna, Italy
96 Infrared Processing and Analysis Center, California Institute of
Technology, Pasadena, CA 91125, USA

#  Euclid: Constraints on f(R) cosmologies from the spectroscopic and
photometric primary probes††thanks: This paper is published on behalf of the
Euclid Consortium.

S. Casas1 0000-0002-4751-5138<EMAIL_ADDRESS>V. F. Cardone2,3 D.
Sapone4 0000-0001-7089-4503 N. Frusciante5 F. Pace6,7,8 0000-0001-8039-0480 G.
Parimbelli9,10,11,12 0000-0002-2539-2472 M. Archidiacono13,14 K. Koyama15 I.
Tutusaus16,17,18,19 0000-0002-3199-0399 S. Camera6,7,8 0000-0003-3399-3574 M.
Martinelli2,3 0000-0002-6943-7732 V. Pettorino20 Z. Sakr21,16,22
0000-0002-4823-3757 L. Lombriser17 A. Silvestri23 0000-0001-6904-5061 M.
Pietroni24,25 0000-0001-5480-5996 F. Vernizzi26 M. Kunz17 0000-0002-3052-7394
T. Kitching27 0000-0002-4061-4598 A. Pourtsidou28,29 0000-0001-9110-5550 F.
Lacasa17,30 0000-0002-7268-3440 C. Carbone31 0000-0003-0125-3563 J. Garcia-
Bellido32 0000-0002-9370-8360 N. Aghanim30 B. Altieri33 0000-0003-3936-0284 A.
Amara15 N. Auricchio34 0000-0003-4444-8651 M. Baldi35,34,36
0000-0003-4145-1943 C. Bodendorf37 E. Branchini38,39 0000-0002-0808-6908 M.
Brescia5,40 0000-0001-9506-5680 J. Brinchmann41 0000-0003-4359-8797 V.
Capobianco8 0000-0002-3309-7692 J. Carretero42,43 0000-0002-3130-0204 M.
Castellano2 0000-0001-9875-8263 S. Cavuoti40,44 0000-0002-3787-4196 A.
Cimatti45 R. Cledassou46,47 0000-0002-8313-2230 G. Congedo28
0000-0003-2508-0046 C. J. Conselice48 L. Conversi49,33 0000-0002-6710-8476 Y.
Copin50 0000-0002-5317-7518 L. Corcione8 0000-0002-6497-5881 F. Courbin51
0000-0003-0758-6510 H. M. Courtois52 0000-0003-0509-1776 A. Da Silva53,54
0000-0002-6385-1609 H. Degaudenzi55 0000-0002-5887-6799 F. Dubath55
0000-0002-6533-2810 C. A. J. Duncan56,48 X. Dupac33 S. Dusini57
0000-0002-1128-0664 S. Farrens58 0000-0002-9594-9387 S. Ferriol50 P.
Fosalba18,19 0000-0002-1510-5214 M. Frailis10 0000-0002-7400-2135 E.
Franceschi34 0000-0002-0585-6591 M. Fumana31 0000-0001-6787-5950 S. Galeotta10
0000-0002-3748-5115 B. Garilli31 0000-0001-7455-8750 W. Gillard59
0000-0003-4744-9748 B. Gillis28 0000-0002-4478-1270 C. Giocoli34,36 A.
Grazian60 0000-0002-5688-0663 F. Grupp37,61 L. Guzzo13,62,14
0000-0001-8264-5192 S. V. H. Haugan63 0000-0001-9648-7260 F. Hormuth64 A.
Hornstrup65,66 0000-0002-3363-0936 P. Hudelot67 K. Jahnke68
0000-0003-3804-2137 S. Kermiche59 0000-0002-0302-5735 A. Kiessling69
0000-0002-2590-1273 M. Kilbinger20 0000-0001-9513-7138 H. Kurki-Suonio70,71
0000-0002-4618-3063 S. Ligori8 0000-0003-4172-4606 P. B. Lilje63
0000-0003-4324-7794 I. Lloro72 E. Maiorano34 0000-0003-2593-4355 O. Mansutti10
0000-0001-5758-4658 O. Marggraf73 0000-0001-7242-3852 F. Marulli74,34,36
0000-0002-8850-0303 R. Massey75 0000-0002-6085-3780 E. Medinaceli34
0000-0002-4040-7783 Y. Mellier67,76 M. Meneghetti34,36 0000-0003-1225-7084 E.
Merlin2 0000-0001-6870-8900 G. Meylan51 M. Moresco74,34 0000-0002-7616-7136 L.
Moscardini74,34,36 0000-0002-3473-6716 E. Munari10 0000-0002-1751-5946 S.-M.
Niemi77 C. Padilla42 0000-0001-7951-0166 S. Paltani55 F. Pasian10 K.
Pedersen78 W. J. Percival79,80,81 S. Pires58 0000-0002-0249-2104 G. Polenta82
0000-0003-4067-9196 M. Poncet46 L. A. Popa83 F. Raison37 A. Renzi84,57
0000-0001-9856-1970 J. Rhodes69 G. Riccio40 E. Romelli10 0000-0003-3069-9222
M. Roncarelli34 0000-0001-9587-7822 E. Rossetti35 R. Saglia85,37
0000-0003-0378-7032 B. Sartoris85,10 V. Scottez86,87 A. Secroun59
0000-0003-0505-3710 G. Seidel68 0000-0003-2907-353X S. Serrano19,18,88
0000-0002-0211-2861 C. Sirignano84,57 0000-0002-0995-7146 G. Sirri36
0000-0003-2626-2853 L. Stanco57 0000-0002-9706-5104 J.-L. Starck89
0000-0003-2177-7794 C. Surace90 P. Tallada-Crespí91,43 0000-0002-1336-8328 A.
N. Taylor28 I. Tereno53,92 R. Toledo-Moreo93 0000-0002-2997-4859 F.
Torradeflot43,91 0000-0003-1160-1517 E. A. Valentijn94 L. Valenziano34,95
0000-0002-1170-0104 T. Vassallo85,10 0000-0001-6512-6358 Y. Wang96
0000-0002-4749-2984 J. Weller85,37 0000-0002-8282-2010 J. Zoubian59

###### Abstract

Context. Euclid will provide a powerful compilation of data including
spectroscopic redshifts, the angular clustering of galaxies, weak lensing
cosmic shear, and the cross-correlation of these last two photometric
observables. This will lead to very stringent constraints on the
$\Lambda\mathrm{CDM}$ concordance cosmological model and models beyond it, in
which for instance standard gravity is modified.

Aims. In this study we extend recently presented Euclid forecasts into the Hu-
Sawicki $f(R)$ cosmological model, a popular extension of the Hilbert-Einstein
action that introduces an universal modified gravity force in a scale-
dependent way. This scale-dependent modification requires a generalisation of
our previous recipes for spectroscopic and photometric galaxy clustering, and
for weak lensing, both in the linear and in the non-linear regimes. Our aim is
to estimate how well future Euclid data will be able to constrain the extra
parameter of the theory, $f_{R0}$, for the range in which this parameter is
still allowed by current observations.

Methods. For the spectroscopic probe, we use a phenomenological approach to
account for the scale dependence of the growth of perturbations in the terms
related to baryon acoustic oscillations and redshift-space distortions. For
the photometric observables, which probe deeper into the non-linear regime, we
use a fitting formula developed in the literature that captures the
modifications in the non-linear matter power spectrum caused by the $f(R)$
model.

Results. We show that, in an optimistic setting, and for a fiducial value of
$|f_{R0}|=5\times 10^{-6}$, Euclid alone will be able to constrain the
additional parameter $\operatorname{log_{10}}|f_{R0}|$ at the $3\%$ level,
using spectroscopic galaxy clustering alone; at the $1.4\%$ level, using the
combination of photometric probes on their own; and at the $1\%$ level, using
the combination of spectroscopic and photometric observations. This last
constraint corresponds to an error of the order of $6\times 10^{-7}$ at the
$1\sigma$ level on the model parameter $|f_{R0}|=5\times 10^{-6}$. We report
also forecasted constraints on a model with a large value of $|f_{R0}|$,
namely $|f_{R0}|=5\times 10^{-5}$ and a model closer to standard gravity with
$|f_{R0}|=5\times 10^{-7}$ and show that in the optimistic scenario, Euclid
will be able to distinguish these models from $\Lambda\mathrm{CDM}$ at more
than 3$\sigma$.

Conclusions. We find that with a good control of systematic effects and
modelling of the matter power spectrum in the mildly and deeply non-linear
regime, Euclid will be a powerful probe for $f(R)$ models. It will constrain
the scale-dependence of the perturbations as well as a substantial part of the
non-linear regime of structure formation, discerning these models from
$\Lambda\mathrm{CDM}$ .

###### Key Words.:

Cosmology: theory; large-scale structure of Universe; cosmological parameters;
dark energy. Gravitational lensing: weak

##  1 Introduction

The origin of the accelerated expansion of the Universe is still challenging
our understanding of late-time cosmology. A cosmological constant, $\Lambda$,
remains in agreement with current data but its value, when considered as
vacuum energy, does not correspond to theoretical predictions and is rather
considered as a phenomenological parameter that fits the data. An appealing
proposal for an alternative modelling is that of modifying gravitational
interactions felt by particles, either in a universal (same interaction for
all particles) or non-universal way (acting differently on different
particles). In this paper, we investigate one popular scenario which belongs
to the first class, in which the theory of general relativity is modified by
extending the Ricci scalar $R$ in the Hilbert-Einstein action with a general
function of it, $R\rightarrow R+f(R)$. We forecast the ability that the
forthcoming Euclid satellite will have to constrain this scenario, via galaxy
clustering (GC, photometric GCph and spectroscopic GCsp), weak lensing (WL),
or their combination, either of the photometric probes alone (XCph) or all
together. In particular, we produce for the first time validated forecasts on
the Hu-Sawicki $f(R)$ model (Hu & Sawicki, 2007), whose background expansion
mimics that of a cosmological constant model, while differing at the level of
cosmological perturbations: the growth of structure is driven here by a
modification of gravity (MG).

Euclid 111http://www.euclid-ec.org/ is a European Space Agency medium-class
space mission due for launch in 2023. It will carry on-board a near-infrared
spectrophotometric instrument (Costille et al., 2018) and a visible imager
(Cropper et al., 2018) that will allow it to perform both a spectroscopic and
a photometric survey over $15\,000\,\deg^{2}$ of extra-Galactic sky (Laureijs
et al., 2011). The main aim of the mission is to measure the geometry of the
Universe and the growth of structures up to redshift $z\sim 2$ and beyond.

Euclid will include a photometric survey, measuring positions and shapes of
over a billion galaxies, enabling the analysis of WL and GCph. Given the
relatively large redshift uncertainties that we expect from photometric
measurements (compared to spectroscopic observations), these analyses will be
performed via a tomographic approach, in which galaxies are binned into
redshift slices that are considered as two-dimensional (projected) data sets.
On the other hand, the spectroscopic survey will provide very precise radial
measurements of the position of galaxies. Even if the number density will be
lower – compared to the photometric survey – it will allow us to perform a
galaxy clustering analysis in three dimensions, GCsp. The combination of
photometric and spectroscopic surveys will enable a powerful test of the two
independent gravitational potentials that are predicted to be different within
$f(R)$ cosmologies.

We want to quantify the effect of combining the complementary information
obtained from the two probes. It has been shown in Euclid Collaboration:
Blanchard et al. (2020, EC19 hereafter) that combining GCph and its cross-
correlation (XCph) with WL is able to improve the figure of merit by a factor
of three, for dynamical dark energy models. Our goal here is to explore the
impact of cross-correlation on the additional parameter
$f_{R0}\equiv\mathrm{d}f/\mathrm{d}R\,(z=0)$ describing the standard model
extensions within the $f(R)$ model, as will be defined in the next section.

Several studies have tried to constrain the Hu-Sawicki model. Among them, Hu
et al. (2016) used Planck15 cosmic microwave background (CMB) data, baryon
acoustic oscillations (BAO), Supernovae Ia from JLA, WiggleZ and CFHTLenS data
sets and obtained an upper bound of $|f_{R0}|<6.3\times 10^{-4}$ at the 95%
confidence level (C.L.). Similar constraints from cosmological data were
obtained for instance in Nunes et al. (2017), Okada et al. (2013) and Pérez-
Romero & Nesseris (2018). A review of local and astrophysical constraints on
$f_{R0}$ is given by Lombriser (2014). Notably, a bound of $|f_{R0}|<10^{-6}$
is obtained assuming that the Milky Way can be treated as an isolated system
in the cosmological background with no environmental screening. Under similar
assumptions, astrophysical tests are capable of constraining $|f_{R0}|$ at the
level of $10^{-7}$ (Baker et al., 2019). Furthermore Desmond & Ferreira (2020)
constrained $|f_{R0}|<1.4\times 10^{-8}$ using galaxy morphology.

After reviewing the $f(R)$ formalism in Sect. 2, we present in Sect. 3 the
Euclid primary probes, WL, $\text{GC}_{\mathrm{ph}}$ and XCph, for the
photometric part and $\text{GC}_{\mathrm{sp}}$ for the spectroscopic part. We
then present the survey specifications and analysis scheme in Sect. 4.
Finally, we present our results for the considered fiducial models in Sect. 5
and conclude in Sect. 6.

##  2 Hu-Sawicki $f(R)$ gravity

A modification of Einstein’s theory of general relativity (GR) can be obtained
by promoting the linear dependence of the Hilbert-Einstein action, $S$, on the
Ricci scalar $R$ to a non-linear function $R+f(R)$ (Buchdahl, 1970),

$\displaystyle S=\frac{c^{4}}{16\pi G_{\rm
N}}\int{\mathrm{d}^{4}x\sqrt{-g}\left[R+f(R)\right]}\,,$ (1)

where $g_{\mu\nu}$ is the metric tensor, $G_{\rm N}$ is Newton’s gravitational
constant and we have expressed explicitly the speed of light $c$, to allow for
consistency with the choice of units in the observable quantities below.

The $f(R)$ family of cosmologies implies a universal coupling with all matter
species, inducing an additional ‘fifth force’. Therefore, an important
attribute that a viable late-time $f(R)$ modification must possess to leave a
detectable signature in the cosmic structure formation while complying with
stringent constraints on gravity in the Solar System, is that the functional
form $f(R)$ gives rise to a screening mechanism, the so-called ‘chameleon
mechanism’ (Khoury & Weltman, 2004).

The fifth force has a range determined by the Compton wavelength $\lambda_{C}$
which has a very direct relation with the parameter $f_{R0}$. For cosmological
densities one has $\lambda_{C}=32\sqrt{|f_{R0}|/10^{-4}}$ Mpc (Hu & Sawicki
(2007) and Cabre et al. (2012)). This relation is important because it is the
screening scale that prevents us to have $f_{R0}=0$ as the fiducial model and
thus being able to forecast what is the minimum $f_{R0}$ which is detectable.

Figure 1: Relative difference between EFTCAMB and MGCAMB for the linear matter
power spectrum ($\Delta
P/P\equiv(P_{\texttt{EFTCAMB}}-P_{\texttt{MGCAMB}})/P_{\texttt{MGCAMB}}$) at
three different redshifts, $z=0$ (solid orange line), $z=1$ (dashed blue line)
and $z=3$ (dot-dashed purple line), for the three fiducial models with
$|f_{R0}|=5\times 10^{-5}$ (HS5), $|f_{R0}|=5\times 10^{-6}$ (HS6) and
$|f_{R0}|=5\times 10^{-7}$ (HS7).

As a specific example of this class of theories, we consider here the model
proposed by Hu & Sawicki (2007). The functional form of $f(R)$, adopting $n=1$
for simplicity and the limit $|f_{R}|\ll 1$, is given by

$f(R)=-6\Omega_{\mathrm{DE},0}\frac{H_{0}^{2}}{c^{2}}+|f_{R0}|\frac{\bar{R}_{0}^{2}}{R}\,,$
(2)

where $f_{R0}<0$, $\bar{R}_{0}$ denoting the Ricci scalar in the cosmological
background today, $H_{0}$ the Hubble constant, and $\Omega_{\mathrm{DE},0}$ is
the current fractional energy density attributed to a cosmological constant.
The only additional free parameter of the model over $\Lambda\mathrm{CDM}$ is
therefore $f_{R0}$.

For the $|f_{R0}|\ll 1$ values of interest here, the background expansion
history approximates that of $\Lambda\mathrm{CDM}$ and

$\bar{R}_{0}=3\Omega_{\mathrm{m},0}\frac{H_{0}^{2}}{c^{2}}\left(1+4\frac{\Omega_{\mathrm{DE},0}}{\Omega_{\mathrm{m},0}}\right)\,,$
(3)

with matter energy density parameter
$\Omega_{\mathrm{m},0}=1-\Omega_{\mathrm{DE},0}$. It characterises the
magnitude of the deviation from $\Lambda\mathrm{CDM}$ , with smaller
$|f_{R0}|$ values corresponding to weaker departures from GR.
$\Lambda\mathrm{CDM}$ is recovered in the limit of $f_{R0}\rightarrow 0$.

At the perturbation level, deviations from GR can be encoded in
phenomenological functions of the metric (Zhang et al., 2007; Amendola et al.,
2008; Planck Collaboration, 2016). Using the Bardeen formalism (Ma &
Bertschinger, 1995), we can define the conformal metric of the infinitesimal
line element, $\mathrm{d}s$, in an expanding Universe as

$\mathrm{d}s^{2}=a^{2}(\tau)\left[-(1+2\Psi)\,c^{2}\,\mathrm{d}\tau^{2}+(1-2\Phi)\,\mathrm{d}x^{i}\,\mathrm{d}x_{i}\right]\,,$
(4)

where $a(\tau)$ is the scale factor in conformal time, $\tau$,
$\mathrm{d}x^{i}$ is the three dimensional infinitesimal spatial element, and
$\Psi$ and $\Phi$ are the two scalar potentials. Then the phenomenological
functions describe the modifications to the Poisson equations, namely

$\displaystyle-k^{2}\Psi$ $\displaystyle=\frac{4\pi\,G_{\rm
N}}{c^{2}}\,a^{2}\mu\left[\bar{\rho}\Delta+3\left(\bar{\rho}+\frac{\bar{p}}{c^{2}}\right)\sigma\right]\,,$
(5) $\displaystyle k^{2}\left(\Phi-\eta\Psi\right)$
$\displaystyle=\frac{12\pi\,G_{\rm
N}}{c^{2}}\,a^{2}\mu\left(\bar{\rho}+\frac{\bar{p}}{c^{2}}\right)\sigma\,,$
(6) $\displaystyle-k^{2}\left(\Phi+\Psi\right)$
$\displaystyle=\frac{8\pi\,G_{\rm
N}}{c^{2}}\,a^{2}\left\\{\Sigma\left[\bar{\rho}\Delta+3\left(\bar{\rho}+\frac{\bar{p}}{c^{2}}\right)\sigma\right]-\right.$
$\displaystyle\left.\qquad\qquad\quad\frac{3}{2}\mu\left(\bar{\rho}+\frac{\bar{p}}{c^{2}}\right)\sigma\right\\}\,,$
(7)

where the background quantities $\bar{\rho}$ and $\bar{p}$ are respectively
the density and pressure of the matter species and are only a function of
time, whereas perturbations are functions of time and scale, $\sigma$ is the
matter anisotropic stress,
$\bar{\rho}\Delta=\bar{\rho}\delta+3(aH/k)(\bar{\rho}+\bar{p}/c^{2})v$ is the
comoving density perturbation, with $\delta=\rho/\bar{\rho}-1$ the density
contrast, and $v$ is the velocity potential. The phenomenological functions
$\mu(a,k),\,\eta(a,k),\,\Sigma(a,k)$ are identically equal to 1 in the GR
limit. Notice that only two of them are independent from each other, the third
one being a combination of the other two; in the limit of negligible
anisotropic stress from matter, the relation reduces to

$\Sigma(a,k)=\frac{\mu(a,k)}{2}\left[1+\eta(a,k)\right]\,.$ (8)

These phenomenological functions can be determined analytically considering
the quasi-static limit (i.e. scales sufficiently small to be well within the
horizon and the sound-horizon of the scalar field). In the case of $f(R)$
gravity, the expressions reflect the presence of an additional fifth force
with a characteristic mass scale

$m_{f_{R}}^{2}\sim\frac{1+f_{R}}{3f_{RR}}\sim\frac{1}{3f_{RR}}\,.$ (9)

For negligible matter anisotropic stress, one finds (Pogosian & Silvestri,
2008)

$\displaystyle\mu(a,k)$
$\displaystyle=\frac{1}{1+f_{R}(a)}\frac{1+4k^{2}a^{-2}m_{f_{R}}^{-2}(a)}{1+3k^{2}a^{-2}m_{f_{R}}^{-2}(a)}\,,$
(10)

and for the Hu-Sawicki model under consideration, $m_{f_{R}}$ is given by
(Brax & Valageas, 2013)

$m_{f_{R}}(a)=\frac{H_{0}}{c\sqrt{2|f_{R0}|}}\frac{\left(4\Omega_{\mathrm{DE},0}+\Omega_{\mathrm{m},0}a^{-3}\right)^{3/2}}{4\Omega_{\mathrm{DE},0}+\Omega_{\mathrm{m},0}}\,.$
(11)

Since $f(R)$ models have a conformal coupling, light deflection is weakly
affected as follows

$\Sigma(a)=\frac{1}{1+f_{R}(a)}\,,$ (12)

and weak lensing is affected in the same way as matter growth, but with a
different weight in time and scale.

This approach is at the basis of the Einstein-Boltzmann solver MGCAMB (Zhao et
al., 2009; Hojjati et al., 2011; Zucca et al., 2019), or MGCLASS Baker & Bull
(2015); Sakr & Martinelli (2022), each a modification of the standard
Einstein-Boltzmann solver CAMB or CLASS, respectively. Note that here we
assume the pre-factor $1/(1+f_{R})$ in Eq. 10 to be unity (Hojjati et al.,
2016). This approximation is valid for viable values of $f_{R0}$. Given our
choice of the fiducial values, $|f_{R0}|\ll 1$, the deviation of $\Sigma$ from
unity is also negligible.

Alternatively, the phenomenological functions, $\mu$, $\eta$, and $\Sigma$,
can be determined numerically, after solving for the full dynamics of linear
perturbations via EFTCAMB (Hu et al., 2014; Raveri et al., 2014), which
implements the effective field theory formalism for dark energy into the
standard CAMB code (Lewis et al., 2000); see Hu et al. (2016) for an
application to Hu-Sawicki $f(R)$ gravity. This code has been validated as part
of an extended code comparison effort (Bellini et al., 2018).

For the model under consideration, we have compared predictions of the angular
power spectrum up to $\ell=5000$ of the CMB and the matter power spectrum up
to $k=10\,h\,\mathrm{Mpc}^{-1}$ both from MGCAMB (quasi-static) and EFTCAMB
(full evolution). Both codes lead to a sub-percent agreement, well within the
desired level of accuracy. For the range of values of $|f_{R0}|$ considered in
this work the agreement of the angular power spectra is never worse than
$0.25\%$ for the temperature-temperature power spectrum (for $\ell<10^{3}$ it
is below 0.1%) and 0.1% for the lensing power spectrum, for the matter power
spectrum the two codes agree extremely well up to
$k=0.02\,h\,\mathrm{Mpc}^{-1}$ ($<0.1\%$) and for larger $k$ the relative
difference is always below $0.25\%$. We show in Fig. 1 for the matter power
spectrum the relative difference between EFTCAMB and MGCAMB at three different
redshifts. The agreement of the two codes has been tested against the choice
of GR transition time, i.e. the time at which a MG model starts to deviate
significantly from its GR limit. We find that the level of agreement is not
affected by this parameter (once this is the same in both codes). For the
present analysis we set the GR transition time at $a=10^{-3}$. This choice is
justified by considering that in the Hu-Sawicki model the growth function is
scale dependent and large-$k$ modes show a significant deviation from GR
already at high redshift. Moreover in this case an early GR transition time
guarantees a smooth transition between the MG and GR regimes. We have also
verified that this is the case using MGCLASS. Given the agreement of the
codes, we conclude that the quasi-static approximation for the Hu-Sawicki
$f(R)$ model is a valid assumption, and we proceed with the forecasts using
the inputs from MGCAMB.

Figure 2: Ratio of the scale-dependent matter growth rate $f(k,z)$ in $f(R)$
gravity (for the fiducial value $|f_{R0}|=5\times 10^{-6}$) with respect to
$\Lambda\mathrm{CDM}$, for three different wavenumbers $5\times 10^{-3}$ (blue
solid line), $5\times 10^{-2}$ (dashed purple line) and $5\times
10^{-1}\,h\,\mathrm{Mpc}^{-1}$ (dot-dashed orange line), as a function of
redshift from $z=0$ to $z=2.5$. The smaller the spatial scales and the lower
the redshifts the larger is the enhancement of the growth rate compared to
$\Lambda\mathrm{CDM}$.

Finally, the scale-dependent $\mu$ function in Eq. (10) introduces a scale-
dependent growth of structures as shown by Zhang (2006) and Song et al.
(2007). In Fig. 2, we plot the growth rate of perturbations,
$f(k,a)\equiv\mathrm{d}\ln\delta/\mathrm{d}\ln a$, in $f(R)$ Hu-Sawicki
gravity with respect to the one in $\Lambda\mathrm{CDM}$, for three different
scales $k$, namely $5\times 10^{-3}$ (blue solid line), $5\times 10^{-2}$
(dashed purple line) and $5\times 10^{-1}$ $h\,\mathrm{Mpc}^{-1}$ (dot-dashed
orange line), as a function of redshift from $z=0$ to $2.5$. This shows that
the growth of perturbations at large scales is very similar to the one in
standard GR, at the redshifts of interests for large-scale structure
formation. However, at smaller scales (larger $k$) the growth of perturbations
is enhanced in $f(R)$ at low redshifts (see Bueno Belloso et al., 2011, for a
parametrization of the growth rate in generic $f(R)$ models, in terms of the
growth index $\gamma$). These very features also complicate the observational
modelling, which moreover serves as a stress test of the forecasting pipeline
for constraints on theories beyond $\Lambda\mathrm{CDM}$.

##  3 Theoretical predictions for Euclid observables

As it will be described in the next section, the forecasting methods and tools
used in this paper are the same of EC19. However, we must notice here that the
change in the theory of gravity introduced through the Hu-Sawicki model
implies significant modifications of the recipes used to compute theoretical
predictions for Euclid observables. We discuss in this Section how moving away
from the standard GR assumption impacts the predictions for the angular power
spectra $C(\ell)$ that will be compared with the photometric survey data, and
on the power spectra $P_{\rm obs}$ compared with data of the spectroscopic
survey.

###  3.1 Photometric survey

For the Euclid photometric survey, the observables that need to be computed
and compared with the data are the angular power spectra for WL,
$\text{GC}_{\mathrm{ph}}$ and their cross correlation, XCph.

In EC19 these were calculated using the Limber approximation plus the flat-sky
approximation with pre-factor set to unity in a flat $\Lambda\mathrm{CDM}$
Universe, as

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

with $k_{\ell}=(\ell+1/2)/r(z)$, $r(z)$ the comoving distance to redshift
$z=1/a-1$, $P_{\delta\delta}(k_{\ell},z)$ the non-linear power spectrum of
matter density fluctuations, $\delta$, at wave number $k_{\ell}$ and redshift
$z$, in the redshift range of the integral from $z_{\rm min}=0.001$ to $z_{\rm
max}=4$. The dimensionless Hubble function is defined as $E(z)=H(z)/H_{0}$ and
in all subsequent equations $H_{0}$ is expressed in units of ${\rm km}\,{\rm
s}^{-1}\,{\rm Mpc}^{-1}$.

For each tomographic redshift bin $i$, the window functions $W^{X}_{i}(z)$
with $X=\\{{\rm L,G}\\}$ (corresponding to WL and GCph, respectively) need to
be computed differently with respect to what was done in EC19 as, when
abandoning the assumption of a GR gravity theory, one has to account for
changes in the evolution of both the homogeneous background and of
cosmological perturbations. However, in the case of the $f(R)$ model
considered in this work, the background is $\Lambda\mathrm{CDM}$ up to a high
precision. In general, one also has to account for both the modified evolution
of the Bardeen potentials, $\Phi$ and $\Psi$, and for the fact that in MG the
GR relation $\Phi=\Psi$ is not necessarily satisfied. Using the modified
Poisson equation for $\Phi+\Psi$ of Eq. 7, this combination can be related to
$P_{\delta\delta}$ as

$P_{\Phi+\Psi}(k,z)=\left[-3\Omega_{\mathrm{m},0}\left(\frac{H_{0}}{c}\right)^{2}(1+z)\Sigma(k,z)\right]^{2}P_{\delta\delta}(k,z)\,,$
(14)

where we assume a standard background evolution of the matter component, i.e.
$\rho_{\rm m}(z)=\rho_{{\rm m},0}(1+z)^{3}$, and $P_{\delta\delta}$ is
computed accounting for the MG effects introduced through Eq. 5.

We can, therefore, use the recipe of Eq. 13 accounting for the effects of
these modifications of gravity, not general modifications, e.g. if matter
coupling shifts $\Omega_{\rm m}(a)$. We can calculate $H$, $r$ and
$P_{\delta\delta}$, provided by dedicated Boltzmann solvers, but with the new
window functions (Spurio Mancini et al., 2019)

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

where $n_{i}(z)/\bar{n}_{i}$ and $b_{i}(k,z)$ are, respectively, the
normalised galaxy distribution and the galaxy bias in the $i$-th redshift bin,
and $W^{\rm IA}_{i}(k,z)$ encodes the contribution of intrinsic alignments
(IA) to the WL power spectrum. We follow EC19 in assuming an effective scale-
independent galaxy bias. The main reason for this choice is to be able to
compare to the standard analysis with the modified gravity model as the only
variable. Accounting for a scale-dependent galaxy bias would introduce further
degrees of freedom that could confuse the comparison between the different
cosmological models. A detailed analysis of both the concordance model and
modified gravity theories accounting for scale-dependent galaxy bias is beyond
the scope of this work (see e.g. Tutusaus et al., 2020, for an analysis on the
concordance model with a local, non-linear galaxy bias model).

The IA contribution is computed following the eNLA model from EC19, in which

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

where

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

with $\langle L\rangle(z)$ and $L_{\star}(z)$ redshift-dependent mean and the
characteristic luminosity of source galaxies as computed from the luminosity
function, $\mathcal{A}_{\rm IA}$, $\beta_{\rm IA}$ and $\eta_{\rm IA}$ are the
nuisance parameters of the model, and $\mathcal{C}_{\rm IA}$ is a constant
accounting for dimensional units.

Changes in the theory of gravity impact the IA contribution introducing a
scale dependence through the modified perturbations growth. This is explicitly
taken into account in Eq. 17 through the matter perturbation $\delta(k,z)$,
which is considered to be scale dependent in this case. This allows us to
consider also the scale dependence introduced by massive neutrinos, which was
assumed to be negligible in EC19.

Notice that while the window function for WL includes the MG function $\Sigma$
in order to properly account for the modifications to $\Phi+\Psi$, the
$\text{GC}_{\mathrm{ph}}$ one does not have any explicit MG contribution, as
the modifications on the clustering of matter are accounted for in the new
$P_{\delta\delta}(k_{\ell},z)$.

We can, therefore, apply this recipe to the Hu-Sawicki $f(R)$ model. Given our
choice of the fiducial $f_{R0}$, discussed in Sect. 4, the background
modifications with respect to $\Lambda\mathrm{CDM}$ are negligible, while this
model affects the evolution of perturbations, and therefore
$P_{\delta\delta}$, through Eq. 5, with the $\mu$ function given by Eq. 10.

The function $\Sigma$ needs to relate the $\Phi+\Psi$ and matter power
spectra, as in Eq. 14, and its is given by Eq. 12. As previously discussed for
our fiducial choice, $|f_{R0}|\ll 1$, the deviations of $\Sigma$ from unity
are negligible and therefore the geometrical part of lensing kernel entering
Eq. 16 reduces to the standard one.

###  3.2 Spectroscopic survey

In order to exploit data from the Euclid spectroscopic survey, we need to
compute the theoretical prediction for the observed galaxy power spectrum in
the extended model considered here.

The full non-linear model for the observed galaxy power spectrum is given by

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

where the $P_{\rm dw}(k,\mu;\,z)$ is the de-wiggled power spectrum which
models the smearing of the BAO features due to the displacement field of
wavelengths smaller than the BAO scale,

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

where the $P_{\rm nw}(k;z)$ is a ‘no-wiggle’ power spectrum with the same
broad band shape as $P_{\delta\delta}(k;z)$ but without BAO features (see
below for details on how we compute it).

In Eq. 19, $k$ is the modulus of the wave vector $\bm{k}$ and $\mu_{\theta}$
is the cosine of the angle $\theta$ between this vector and the line-of-sight
direction $\hat{\bm{r}}$. These quantities on the right hand side are
functions of their counterparts at a reference cosmology, i.e. $k\equiv
k(k_{\rm ref})$, $\mu_{\theta}\equiv\mu_{\theta,\text{ref}}$, which are
transformed due to the Alcock-Paczynski effect, see EC19 and Casas et al.
(2023) for the explicit formula. This transform, which also scales the overall
$P_{\text{obs}}$ is parameterised in terms of the angular diameter distance
$D_{\rm A}(z)$ and the Hubble parameter $H(z)$ as

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

The term in the curly brackets in Eq. 19 is the contribution of redshift space
distortions (RSD) corrected for the non-linear Finger-of-God (FoG) effect,
where we defined $b\sigma_{8}(k,z)$ as the product of the effective scale-
dependent bias of galaxy samples and the r.m.s. matter density fluctuation
$\sigma_{8}(z)$; similarly, $f\sigma_{8}(k,z)$ is the product of the scale-
dependent growth rate and $\sigma_{8}(z)$. As in the photometric survey, we
follow EC19 in using an effective scale-independent galaxy bias. An analysis
of modified gravity with scale-dependent galaxy bias models is left for future
work.

The observed galaxy power spectrum is modulated by the redshift uncertainties
which is manifested as a smearing of the galaxy density field along the line-
of-sight, hence the factor $F_{z}$ in Eq. 19 reads

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

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

Finally, the $P_{\rm s}(z)$ is a scale-independent shot noise term, which
enters as a nuisance parameter (see EC19, ).

The change in the gravity model affects the way to compute the theoretical
predictions, as these expressions need to account for the possibility for the
growth rate $f(z)$ to also depend on the wave number $k$. In general, this is
the case for any MG model, and also when perturbation in the dark energy
sector are considered. For the model considered in this paper, the only terms
affected are those directly related to the growth rate, in brief: $f(k,z)$
itself, and the two phenomenological parameters related to the velocity
dispersion, $\sigma_{\rm v}$, and the pairwise velocity dispersion,
$\sigma_{\rm p}$. These are

$\displaystyle\sigma^{2}_{\rm v}(z,\mu_{\theta})$
$\displaystyle=\frac{1}{6\pi^{2}}\int\mathrm{d}k\,P_{\delta\delta}(k,z)\left\\{1-\mu_{\theta}^{2}+\mu_{\theta}^{2}\left[1+f(k,z)\right]^{2}\right\\},$
(24) $\displaystyle\sigma_{\rm p}^{2}(z)$
$\displaystyle=\frac{1}{6\pi^{2}}\int\mathrm{d}k\,P_{\delta\delta}(k,z)f^{2}(k,z)\,.$
(25)

The phenomenological parameters $\sigma_{\rm v}(z,\mu_{\theta})$ and
$\sigma_{\rm p}(z)$ account for the damping of the BAO features and the FoG
effect, respectively. The smearing of the BAO peak is due to the bulk motion
of scales smaller than the BAO scale. For the power spectrum, this can be
modeled in the Zeldovich approximation by a multiplicative damping term of the
form $\exp\left[-k^{i}k^{j}\langle d^{i}(z)d^{j}(z)\rangle\right]$, where
$\langle d^{i}(z)d^{j}(z)\rangle$ is the correlation function of the
displacement field $d^{i}$ evaluated at zero distance (see for instance the
Appendix C of Peloso et al. 2015 for more details). Finally, we would like to
clarify that in Eq. 20 we used the function $g_{\mu}$ to express the damping
of the BAO features in the matter power spectrum to keep the recipe closer to
the EC19; in this work, $g_{\mu}=\sigma_{\rm v}(z,\,\mu_{\theta})$.

In $\Lambda\mathrm{CDM}$ model the growth rate is scale-independent and both
Eqs. 25 and 24 are the same (see EC19, ). Notice that even if these parameters
are assumed to be the same, they come from two different physical effects,
namely large-scale bulk flow for the former and virial motion for the latter.
Finally, due to the scale dependence of $\sigma_{\rm p}$ and $\sigma_{\rm v}$,
we evaluated both parameters at each redshift bin but we kept them fixed in
the Fisher matrix analysis. This method corresponds to the optimistic settings
in EC19. We would like to highlight that in this work we take directly the
derivatives of the observed galaxy power spectrum with respect to the final
parameters, contrary to EC19 where first we performed the Fisher matrix
analysis for the redshift dependent parameters $H(z)$, $D_{\rm A}(z)$ and
$f\sigma_{8}(z)$ and then projected to the final cosmological parameters of
interest.

The no-wiggle matter power spectrum $P_{\rm nw}(k;z)$ entering Eq. 20 has been
obtained using a Savitzky-Golay filter to the matter power spectrum
$P_{\delta\delta}(k;z)$. The Savitzky-Golay filter is usually applied to noisy
data in order to smooth their behavior. This convolution method consists in
fitting successive sub-sets of adjacent data points with a low-order
polynomial. If the data are equally spaced (as it is in our case, equally
spaced in $\log_{10}k$), then an analytic solution to the least-squares can be
found as series of coefficients that can be applied to all the sub-sets. In
practice, using the Savitzky-Golay filter we recover exactly the same shape
and amplitude of the matter power spectrum without the BAO wiggles. While in
EC19 we used the Eisenstein-Hu fitting formula (Eisenstein & Hu, 1998) for the
no-wiggle power spectrum, this is a fitting formula that only applies
approximately to $\Lambda\textrm{CDM}$ models and therefore cannot be
straightforwardly applied in our case. We find that the Savitzky-Golay (SG)
method is more accurate for models where the growth of matter density field
depends on the scale $k$. The aforementioned smoothing filter has also been
used in previous works (Boyle & Komatsu, 2018) to reconstruct the neutrinos
masses from galaxy redshift survey. In Fig. 8 we plot our reconstruction of
the wiggles computed with the Eisenstein-Hu formula compared to the SG method
used in this work, where we can see that it performs very well also in the
case of $f(R)$.

###  3.3 Non-linear modelling

While for the galaxy power spectrum on mildly non-linear scales we use a
modified version of the model in EC19, we do not have, in general, an
analytical solution for the deeply non-linear power spectrum in a $f(R)$
cosmology. In this work, we will, therefore, use a fitting formula designed in
Winther et al. (2019) that captures the enhancement in the power spectrum
compared to a $\Lambda\mathrm{CDM}$ non-linear power spectrum, as a function
of the parameter $f_{R0}$. This fitting function has been calibrated using the
DUSTGRAIN (Giocoli et al., 2018) and the ELEPHANT (Cautun et al., 2018)
$N$-body simulations (see Winther et al., 2015, for a comparison of different
$N$-body codes for $f(R)$ cosmologies).

The fitting function we use is given by

$\displaystyle\frac{P^{\rm fit}_{f(R)}(k,z)}{P_{\Lambda\mathrm{CDM}}(k,z)}$
$\displaystyle=\Xi(k,z)$ $\displaystyle\equiv
1+X_{1}\frac{1+X_{2}\,k}{1+X_{3}\,k}\arctan(X_{4}\,k)^{X_{5}+X_{6}\,k}\,,$
(26)

where the $\\{X_{i}\\}$ are themselves functions of $f_{R0}$ and redshift,
$X\equiv X(z;f_{R0})$. The redshift dependence is a polynomial relation to the
scale factor $a=1/(1+z)$ given by

$\displaystyle X_{i}(z;f_{R0})=X_{i0}(y)+X_{i1}(y)(a-1)+X_{i2}(y)(a-1)^{2}\,,$
(27)

with each $X_{ij}$ coefficient in itself defined as a polynomial in
$y\equiv\log(f_{R0}/f_{R0}^{\rm fid})$, given by

$\displaystyle X_{ij}(y)=X_{ij0}+X_{ij1}y+X_{ij2}y^{2}\,.$ (28)

With three indices for $X_{ijk}$ this gives in total $6\times 3\times 3=54$
free parameters for the full scale, redshift and $f_{R0}$ dependence. The
response function $\Xi(k,z)$ is the ratio of the non-linear matter power
spectrum in $f(R)$ theory to the non-linear spectrum calculated within the
$\Lambda\mathrm{CDM}$ model. The fitting formula Sect. 3.3 found in Winther et
al. (2019) is a direct fit to this ratio, $\Xi(k,z)$. This fitting formula is
not defined outside the range $10^{-4}<|f_{R0}|<10^{-7}$, therefore this will
limit our smallest $|f_{R0}|$ fiducial value across probes to be $5\times
10^{-7}$, since we need to be far from the lower limit to be able to compute
the numerical derivatives accurately. In Fig. 3 we plot the function
$\Xi(k,z)$ for each of the fiducial $|f_{R0}|$ values chosen in this model,
namely $5\times 10^{-5}$ (blue line), $5\times 10^{-6}$ (purple line) and
$5\times 10^{-7}$ (orange line), as a function of scale $k$. To compare the
enhancement with respect to $\Lambda\mathrm{CDM}$ at our redshifts of
interest, we evaluate it at $z=0.25$ (solid lines) and $z=1.75$ (dashed
lines). These redshifts correspond approximately to the means of the first and
last lensing tomographic bins of our survey, respectively (see Sect. 4). As
can already be seen from this figure, the enhancement with respect to
$\Lambda\mathrm{CDM}$ decreases rapidly with a smaller $|f_{R0}|$ value,
therefore we expect to have worst constraints for smaller values of $|f_{R0}|$
when testing this model against probes that are sensitive to the deeply non-
linear power spectrum. Notice, however, that the Fisher matrix sensitivity is
dominated by the response of this function to small changes in $|f_{R0}|$
(i.e. its first derivatives), which can become large at small scales and
therefore, we can get relatively tight constraints, even if the absolute
enhancement over $\Lambda\mathrm{CDM}$ is of just a few percent evaluated at
the fiducial. In Sect. 4 we will specify our choice of fiducial parameters for
each model.

Figure 3: Ratio of the non-linear power spectrum in $f(R)$ gravity to
$\Lambda\mathrm{CDM}$ from the fitting function in Sect. 3.3 for three
different values of the $|f_{R0}|$ parameter, namely $5\times 10^{-5}$ (blue
line), $5\times 10^{-6}$ (purple line) and $5\times 10^{-7}$ (orange line), as
a function of scale $k$, evaluated at two different redshifts $z=0.25$ (solid
lines) and $z=1.75$ (dashed lines). These redshifts correspond approximately
to the means of the first and last lensing tomographic bins of our survey,
respectively (see Sect. 4). The fitting formula designed in Winther et al.
(2019) is not defined outside $10^{-4}<|f_{R0}|<10^{-7}$.

We implement $\Xi(k,z)$ into the Boltzmann codes MGCAMB and EFTCAMB and in
order to obtain the non-linear $f(R)$ matter power spectrum, we then multiply
$\Xi(k,z)$ by a $\Lambda\mathrm{CDM}$ non-linear spectrum, $P^{\rm
fit}_{f(R)}(k,z)=\Xi(k,z)\,P_{\Lambda\mathrm{CDM}}(k,z)$. For the
$\Lambda\mathrm{CDM}$ power spectrum $P_{\Lambda\mathrm{CDM}}(k,z)$, we use
the Halofit ‘Takahashi’ prescription (Takahashi et al., 2012), since this is
the prescription most readily available in MGCAMB, and also the prescription
used in Winther et al. (2019) to test the fitting formula against $N$-body
simulations. During the preparation of this work, an emulator for the deeply
non-linear matter power spectrum has been developed by the Euclid
collaboration, the EuclidEmulator (see Euclid Collaboration: Knabenhans et
al., 2019, for details on its implementation), calibrated on the Euclid
Flagship simulation (Potter et al., 2017). Whilst it might be interesting to
use it in the future, the current available version, EuclidEmulator2 (see
Euclid Collaboration: Knabenhans et al., 2021), offers a Python wrapper to the
Boltzmann code CLASS (Lesgourgues, 2011; Blas et al., 2011) rather than to
MGCAMB and EFTCAMB we use for this paper; we are therefore not using the
EuclidEmulator in the current analysis (see Sect. 2). Also, during the
preparation of this work, an emulator for the non-linear matter power spectrum
in $|f_{R0}|$ was developed by Arnold et al. (2021) and goes under the acronym
FORGE. It has been checked independently by the authors of this paper that
this emulator agrees well with the fitting formula by Winther et al. (2019)
around the fiducial values of interest. However, this FORGE emulator has been
calibrated with massless neutrino simulations and allows only variations of
the $\Omega_{\mathrm{m},0}$, $\sigma_{8}$ and $h$ cosmological parameters.
Since in this work we want to have the flexibility to vary all other
cosmological parameters and also the ability to connect to accurate modified
gravity Boltzmann codes, we leave the application of this emulator for future
work. The fitting formula in Winther et al. (2019) is also checked to be valid
in the presence of non-zero neutrino masses, using data from Baldi et al.
(2014). With the value $\sum m_{\nu}=0.06\,\mathrm{eV}$ chosen in this work,
this fitting formula is accurate enough across our ranges of scales and
redshifts of interest.

Note that in the present analysis, for the spectroscopic probe, we adopt the
non-linear modelling described above in Eq. 19 as it represents the minimal
modification to what used in the $\Lambda\mathrm{CDM}$ forecasts of EC19. We
choose to follow this more simplistic route, which allows a clean comparison
with EC19 with fixed theoretical systematics. However, it is worth noting that
the state of the art, based e.g. on the Taruya-Nishimichi-Saito model (2010)
or on perturbation theory prescriptions with bias expansion, has been
formulated and recently implemented in the analysis of real data from the BOSS
survey (Song et al., 2015; Colas et al., 2020). Furthermore, eBOSS analyses
(Beutler et al., 2014, 2017; de Mattia et al., 2021), as well as forecasts for
unbiased parameter estimation for Stage IV cosmological surveys, show that the
choice of $k_{\rm max}$ is not universal (see e.g. Markovic et al., 2019), and
different $k_{\rm max}$ for the monopole, quadrupole, and hexadecapole are
required (the $k_{\rm max}$ for the latter being considerably smaller).
Therefore, we study several $k_{\rm max}$ choices, for the full shape of the
power spectrum of spectroscopically-observed galaxies.

For the photometric observables, on the other hand, we are probing up to
smaller scales. For this reason, including nuisance parameters about baryonic
feedback on the matter power spectrum would be necessary for unbiased
parameter estimation, and it would also possibly entail a degradation of the
constraints from WL (see e.g. Schneider et al., 2020, 2020). However, at the
moment we do not have accurate Euclid-like simulations including baryonic
effects, especially in the case of MG cosmologies. Therefore, we ignore these
effects in our analysis, leaving their inclusion for a future work.

We are aware that there are important degeneracies among the effect of $f(R)$
and neutrino masses, especially at non-linear scales (see e.g. Hu et al.,
2015; Baldi et al., 2014). While on the one hand massive neutrinos suppress
the power spectrum at small scales, $f(R)$ will increase clustering at a
similar range of scales, therefore partially cancelling the former effect
(Baldi et al., 2014). Ignoring these degeneracies might indeed artificially
tighten our constraints on $|f_{R0}|$ since we would be considering a much
higher signal than the one actually present under a large sum of neutrino
masses (He, 2013; Motohashi et al., 2013; Harnois-Déraps et al., 2015).
Breaking these degeneracies is possible, either by using statistics on the
cosmic web (Shim et al., 2014), higher than second order statistics in weak
lensing (Peel et al., 2018; Giocoli et al., 2018) or machine learning (Peel et
al., 2019; Merten et al., 2019). These techniques are however beyond the scope
of this work. Galaxy clusters and voids also offer a possibility to
distinguish between these two possible scenarios, as investigated in Hagstotz
et al. (2019a, b); Ryu et al. (2020); Contarini et al. (2021). For our
spectroscopic observable, $\text{GC}_{\mathrm{sp}}$, there is also the
possibility of breaking this degeneracy, by using information on redshift
space distortions. This involves extending the perturbation theory modelling
(Wright et al., 2019) or performing a detailed comparison with galaxy mocks
from simulations (García-Farieta et al., 2019). An exploration of these
extensions on our non-linear modelling of RSD will be left for a future
investigation.

##  4 Survey specifications and analysis method

In order to forecast constraints on this specific model, we follow the same
approach of EC19. We adopt the same Fisher matrix formalism, as well as the
codes validated therein. Given that the theoretical model considered is
crucially different, we update the forecast recipe and the corresponding codes
as described in the previous section. The cosmological parameters here
considered, and their fiducial values, for which we again follow EC19, read

$\displaystyle\bm{\Theta}\phantom{{}_{\rm fid,1}}=\\{\Omega_{\rm
m,0},\,\Omega_{\rm b,0},\,h,\,n_{\rm
s},\,\sigma_{8},\,\operatorname{log_{10}}|f_{R0}|\\}\,,$
$\displaystyle\mathrm{HS5}:\bm{\Theta}_{\rm
fid,HS5}=\\{0.32,\,0.05,\,0.67,\,0.96,\,0.911,\,-4.301\\}\,,\,$
$\displaystyle\mathrm{HS6}:\bm{\Theta}_{\rm
fid,HS6}=\\{0.32,\,0.05,\,0.67,\,0.96,\,0.853,\,-5.301\\}\,,\,$
$\displaystyle\mathrm{HS7}:\bm{\Theta}_{\rm
fid,HS7}=\\{0.32,\,0.05,\,0.67,\,0.96,\,0.823,\,-6.301\\}\,.$ (29)

The fiducial values of $\sigma_{8}$ shown above are obtained keeping the
scalar amplitude of the primordial power spectrum $A_{\rm s}$ fixed for all
three cosmologies. Given the impact of the $f(R)$ model under examination,
this leads to $\sigma_{8}$ values that appear in tension with currently
available results. However, it is important to stress that the comparison
should be done with modified gravity analyses of present data, rather than
with results obtained assuming $\Lambda$CDM. Indeed, the fiducial values of
$\sigma_{8}$ we quote are compatible with such analyses, as degeneracies
between $\sigma_{8}$ and modified gravity parameters increase the mean value
of the former while also broadening the constraints (see e.g. Abbott et al.,
2023).

In the following text and in the figures we will refer to the model with the
corresponding $|f_{R0}|=5\times 10^{-5}$ fiducial as HS5, the model with
$|f_{R0}|=5\times 10^{-6}$ as HS6 and the model with a fiducial value of
$|f_{R0}|=5\times 10^{-7}$ as HS7. The baseline fiducial used for the rest of
this work will be the HS6 model, since it corresponds to a value of $|f_{R0}|$
still allowed by observations and with enough distinctive signatures with
respect to $\Lambda\mathrm{CDM}$ as to be detected by future observations.
Note that our fiducial cosmology includes massive neutrinos with total mass of
$\sum m_{\nu}=0.06\,\mathrm{eV}$, but we keep $\sum m_{\nu}$ fixed in the
following Fisher matrix analysis. We also use the same initial amplitude of
primordial perturbations for both fiducial models, namely $A_{\rm
s}=2.12605\times 10^{-9}$. As discussed previously in Sect. 3.3, fixing
neutrino masses ignores the degeneracies between their effect on the power
spectrum and the increase of clustering at small scales, therefore our
constraints might be tighter than in a scenario in which also $\sum m_{\nu}$
is varied.

Concerning the photometric probes, the galaxy distribution is binned into $10$
equi-populated redshift bins with an overall distribution following

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

with $z_{0}=0.9/\sqrt{2}$ and the normalisation set by the requirement that
the surface density of galaxies is $\bar{n}_{g}=30\,\mathrm{arcmin}^{-2}$. The
redshift distribution is then convolved with a sum of two Gaussian
distributions to account for the photometric redshift uncertainties (see EC19,
for details). The galaxy bias is assumed to be constant within each redshift
bin, and its values $b_{i}$ are introduced as nuisance parameters in our
analysis, with their fiducial values determined by
$b_{i}=\sqrt{1+\bar{z}_{i}}$, where $\bar{z}_{i}$ is the mean redshift of each
redshift bin. Even though deviations from GR introduce in principle a scale
dependence also in the galaxy bias, we assume that this is negligible in our
case.

Moreover, we follow EC19 in accounting for a Gaussian covariance between the
different photometric probes:

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

where upper-case Latin indexes $A,\ldots=\\{{\rm
WL},\,\text{GC}_{\mathrm{ph}}\\}$, lower-case Latin indexes $i,\ldots$ run
over all tomographic bins, $\delta_{\ell\ell^{\prime}}^{\rm K}$ is the
Kronecker delta symbol, $f_{\rm sky}\simeq 0.36$ represents the fraction of
the sky observed by Euclid, and $\Delta\ell$ denotes the width of the
multipole bins, where we use $100$ equi-spaced bins in log-space. The noise
terms, which for the observables considered here are in fact white noise,
namely $N_{ij}^{AB}(\ell)\equiv N_{ij}^{AB}$, read

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

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

For the spectroscopic probe, we evaluate the Fisher matrix
$F_{\alpha\beta}(z_{i})$ for the observed galaxy power spectrum according to
the recipe outlined in EC19 (see Section 3.2). Here, $\alpha$ and $\beta$ run
over the cosmological parameters of the set $\bm{\Theta}$, the index $i$
labels the redshift bin, each respectively centred in
$z_{i}=\\{1.0,\,1.2,\,1.4,\,1.65\\}$, whose widths are $\Delta z=0.2$ for the
first three bins and $\Delta z=0.3$ for the last bin. In this paper, in
comparison to EC19, we adopt the direct derivative approach, i.e. we vary the
observed galaxy power spectrum with respect to the cosmological parameters of
$\bm{\Theta}$ directly, plus two additional redshift-dependent parameters $\ln
b\sigma_{8}(z_{i})$ and $P_{\rm s}(z_{i})$ that we marginalise over. We
consider the numerical values for the galaxy bias, $b(z)$, and the expected
number density of the observed H$\alpha$ emitters, $n(z)$, reported in table 3
of EC19.

For both probes, we consider two different scenarios: an optimistic and a
pessimistic case. In the optimistic case, we consider $k_{\rm
max}=0.30\,h\,\mathrm{Mpc}^{-1}$ for GCsp, $\ell_{\rm max}=5000$ for WL, and
$\ell_{\rm max}=3000$ for $\text{GC}_{\mathrm{ph}}$ and XCph. Instead, in the
pessimistic scenario, we consider $k_{\rm max}=0.25\,h\,\mathrm{Mpc}^{-1}$ for
GCsp, $\ell_{\rm max}=1500$ for WL, and $\ell_{\rm max}=750$ for
$\text{GC}_{\mathrm{ph}}$ and XCph. At the smallest photometric redshift bin,
the galaxy number density distribution $n(z)$ peaks at around $z=0.25$, which
means that under the Limber approximation and for our fiducial cosmology, the
corresponding maximum values of $k$ evaluated in the power spectrum
corresponding to the pessimistic and optimistic scenarios for GCph are $k_{\rm
max}=[0.7,2.9]\,h\,\mathrm{Mpc}^{-1}$, respectively and for WL the maximum
wavenumbers probed are $k_{\rm max}=[1.4,4.8]\,h\,\mathrm{Mpc}^{-1}$ for
pessimistic and optimistic, respectively. For smaller values of $z$, the
values of $k$ at a given $\ell$ increase monotonically, but there the window
functions in Eq. 15 and Eq. 16 suppress the power spectrum and we set it to
zero after a fixed $k_{\rm max}=30\,h\,\mathrm{Mpc}^{-1}$. In both scenarios,
we fix the $\sigma_{\rm p}$ and $\sigma_{\rm v}$ nuisance parameters for
$\text{GC}_{\mathrm{sp}}$, which we calculate directly from Eqs. 25 and 24 for
the fiducial value of the cosmological parameters. We also perform a second
pessimistic forecast for GCsp only, where we set our maximum wave number at
$k_{\rm max}=0.15\,h\,{\rm Mpc}^{-1}$ in order to have a more conservative
estimate of the constraining power of the GCsp probe. The reason for this is
that the underlying matter power spectrum of Eq. 19 that we are using in our
observed galaxy power spectrum recipe is a linear one, as we detailed in Sect.
3.2. It is known that non-linear corrections start playing a role above scales
of around $k=0.1\,h\,{\rm Mpc}^{-1}$ for the redshifts under consideration
(see Taruya et al., 2010) and, therefore, the use of a linear power spectrum
beyond these scales can bias our constraints. Hence, we aim to estimate what
would happen if we used just quasi-linear scales in the analysis. As shown in
Sect. 5, the effect of these two different scale cuts on the constraining
power on the $\operatorname{log_{10}}|f_{R0}|$ parameter is minimal for our
GCsp recipe. As a reference for the reader, we list the specific choices of
scales and settings used for each observable in Table 1.

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

In this work, as in EC19, we show the results for most of these single probes,
but also for their combinations. It is important to mention that when we
consider the combination of $\text{GC}_{\mathrm{ph}}$ with WL, we neglect any
cross-correlation. However, when we add their cross-correlation XCph, we
include it both in the data vector and in the covariance, i.e. we perform a
full analysis taking into account the cross-covariances between
$\text{GC}_{\mathrm{ph}}$, WL, and their cross-correlation. This combination
of three distinct two-point correlation functions is also known in the
literature as 3x2pt and we will use this terminology interchangeably in this
work. Moreover, again following EC19, we do not present the values for
$\text{GC}_{\mathrm{ph}}$ alone. The main reason for this choice is that we
consider both $\sigma_{8}$ and the galaxy biases as parameters, which are
degenerate in the linear regime. Even if this degeneracy might be partially
broken when adding non-linear information, the Fisher formalism can still
manifest numerical instabilities for this single probe alone. Therefore, we
always show the constraints from GCph in combination with other probes. In the
optimistic scenario, we assume GCsp to be uncorrelated to photometric probes.
In the pessimistic setting, we neglect any correlation between GCsp and WL,
and also apply a redshift cut at $z<0.9$ for GCph and XCph, in order to
minimise the overlap between the different galaxy clustering probes. Note
however, that this redshift cut is only applied when combining spectroscopic
and photometric data. Even in the pessimistic case, we do not apply any
redshift cut when considering photometric data alone.

Finally, in this work we have re-validated the codes used in EC19 to account
for the modified recipe outlined in the previous sections. We have compared,
for the constraint on each parameter of $\bm{\Theta}$ and the nuisance
parameters of each probe, the performance of each independent code to the
median of the constraints obtained by the available codes. We have verified
that the relative difference in percentage between a given code and such a
median is always below 10 %, both for the marginalized and unmarginalized
parameters, which was the threshold for code consistency adopted in EC19 and
more recently in Casas et al. (2023).

Table 2: Forecast $1\sigma$ marginal relative errors on the cosmological
parameters for a flat $f(R)$ model with $|f_{R0}|=5\times 10^{-6}$
($\operatorname{log_{10}}|f_{R0}|=-5.301$) in the pessimistic and optimistic
cases, using Euclid observations of spectroscopic Galaxy Clustering
($\text{GC}_{\mathrm{sp}}$), Weak Lensing (WL), photometric Galaxy Clustering
($\text{GC}_{\mathrm{ph}}$) and the cross-correlation among the photometric
probes XCph. $|f_{R0}|=5\times 10^{-6}$
---
| $\Omega_{\rm m,0}$ | $\Omega_{\rm b,0}$ | $\operatorname{log_{10}}|f_{R0}|$ | $h$ | $n_{\rm s}$ | $\sigma_{8}$
Fiducial values | $0.32$ | $0.05$ | $-5.301$ | $0.67$ | $0.96$ | $0.853$
Pessimistic setting
GCsp$(k_{\rm max}=0.15\,h\,{\rm Mpc}^{-1})$ | $2.3\%$ | $4.6\%$ | $6.3\%$ | $0.6\%$ | $1.7\%$ | $1.9\%$
GCsp$(k_{\rm max}=0.25\,h\,{\rm Mpc}^{-1})$ | $1.4\%$ | $2.6\%$ | $3.6\%$ | $0.3\%$ | $1.1\%$ | $1.1\%$
WL | $2.3\%$ | $47\,\,\%$ | $8.3\%$ | $21\,\,\%$ | $4.8\%$ | $1.5\%$
GCph | $2.4\%$ | $5.8\%$ | $9.4\%$ | $4.1\%$ | $3.4\%$ | $1.9\%$
GCsp+WL | $0.7\%$ | $1.7\%$ | $2.3\%$ | $0.2\%$ | $0.8\%$ | $0.3\%$
GCph+WL | $1.0\%$ | $5.1\%$ | $3.8\%$ | $3.4\%$ | $1.9\%$ | $0.5\%$
GCsp+WL+GCph | $0.6\%$ | $1.6\%$ | $2.2\%$ | $0.2\%$ | $0.7\%$ | $0.3\%$
WL+GCph+XCph | $0.8\%$ | $5.0\%$ | $2.7\%$ | $3.3\%$ | $1.7\%$ | $0.4\%$
GCsp+WL+GCph+XCph | $0.6\%$ | $1.6\%$ | $1.7\%$ | $0.2\%$ | $0.7\%$ | $0.3\%$
Optimistic setting
GCsp$(k_{\rm max}=0.30\,h\,{\rm Mpc}^{-1})$ | $1.3\%$ | $2.3\%$ | $3.0\%$ | $0.3\%$ | $1.0\%$ | $1.0\%$
WL | $1.5\%$ | $45\,\,\%$ | $4.7\%$ | $19\,\,\%$ | $3.2\%$ | $0.9\%$
GCph | $1.5\%$ | $4.4\%$ | $1.7\%$ | $2.9\%$ | $0.7\%$ | $0.2\%$
GCsp+WL | $0.5\%$ | $1.5\%$ | $2.0\%$ | $0.2\%$ | $0.6\%$ | $0.3\%$
GCph+WL | $0.4\%$ | $4.3\%$ | $1.6\%$ | $2.1\%$ | $0.7\%$ | $0.2\%$
GCsp+WL+GCph | $0.3\%$ | $1.2\%$ | $1.2\%$ | $0.1\%$ | $0.3\%$ | $0.2\%$
WL+GCph+XCph | $0.3\%$ | $4.3\%$ | $1.4\%$ | $2.1\%$ | $0.7\%$ | $0.2\%$
GCsp+WL+GCph+XCph | $0.3\%$ | $1.2\%$ | $1.0\%$ | $0.1\%$ | $0.3\%$ | $0.2\%$

##  5 Results

Figure 4:  $1$ and $2\sigma$ joint marginal error contours on the cosmological
parameters for a flat $f(R)$ model with $|f_{R0}|=5\times 10^{-6}$ in the
optimistic scenario. Purple is for GCsp, blue for WL, orange for the
combination GCsp+WL, and yellow for all the photometric probes including their
cross-correlation XCph, combined with GCsp, namely GCsp+WL+GCph+XCph. While
the WL probe is unable to properly constrain the Hubble parameter $h$ and the
primordial slope $n_{\rm s}$ on its own, the orthogonality of the contours for
GCsp and WL in the subspaces involving $h$ and $n_{\rm s}$, helps to lift
degeneracies and further improves the fully-marginalized constraints on
$\operatorname{log_{10}}|f_{R0}|$, when probe combinations are used.

### Optimistic scenario

As shown in Eq. 29, we have chosen three different fiducial values of the Hu-
Sawicki $f(R)$ model parameter, namely $|f_{R0}|=5\times 10^{-7}$ (HS7),
$|f_{R0}|=5\times 10^{-6}$ (HS6) and $|f_{R0}|=5\times 10^{-5}$ (HS5), based
on our discussion of the current observational constraints in Sect. 1. We
employ HS7, which contains a very small value of $|f_{R0}|$ as our GR-limit
test, since we cannot correctly perform forecasts at lower values of
$|f_{R0}|$ due to the limitations we have in the non-linear modelling with the
Winther fitting formula mentioned in Sect. 3.3. We remind the reader of the
two considered scenarios, pessimistic and optimistic, as explained in Sect. 4
plus the ‘quasi-linear’ scenario for GCsp, defined by $k_{\rm
max}=0.15\,h\,{\rm Mpc}^{-1}$.

For our baseline fiducial (HS6) and in an optimistic setting, Euclid alone
will be able to constrain the additional parameter
$\operatorname{log_{10}}|f_{R0}|$, which has a value of
$\operatorname{log_{10}}|f_{R0}|=-5.301$ at the $1\sigma$ level with an
absolute error of

  * •

$\sigma_{\operatorname{log_{10}}|f_{R0}|}=0.16$ with spectroscopic GCsp alone
(corresponding to a relative $3.0\%$ error);

  * •

$\sigma_{\operatorname{log_{10}}|f_{R0}|}=0.25$ with WL alone
(corresponding to a relative $4.7\%$ error);

  * •

$\sigma_{\operatorname{log_{10}}|f_{R0}|}=0.07$
combining WL, GCph, and XCph
(corresponding to a relative $1.4\%$ error);

  * •

$\sigma_{\operatorname{log_{10}}|f_{R0}|}=0.05$
using the full combination GCsp+WL+GCph+XCph
(corresponding to a relative $1.0\%$ error).

Figure 5: Marginalised $1\sigma$ errors on cosmological parameters, relative
to their corresponding fiducial value for the optimistic case. Corresponding
to the value of $|f_{R0}|=5\times 10^{-6}$. We show results for GCsp (purple),
WL (blue), GCsp+WL (orange), the combination of all photometric probes,
including cross correlations, WL+GCph+XCph (red) and the combination of all
spectroscopic and photometric probes GCsp+WL+GCph+XCph (yellow). We can see
that the constraints on $\operatorname{log_{10}}|f_{R0}|$ and $\sigma_{8}$ are
very similar among GCsp and WL, but for other parameters like the Hubble
parameter $h$ or the fraction of baryons $\Omega_{\rm b,0}$ the constraints
coming from GCsp alone are more stringent than the constraints of all
photometric probes combined on their own.

In Table 2 and in Fig. 5, we list the forecasted $1\sigma$ fully-marginalized
errors (relative errors to its fiducial) on all the cosmological parameters
considered for our model with $|f_{R0}|=5\times 10^{-6}$ (HS6), for the
individual probes and their combinations in the optimistic scenario. The
probes shown are $\text{GC}_{\mathrm{sp}}$ (in purple), WL (in blue),
$\text{GC}_{\mathrm{sp}}$+WL (orange), the 3x2pt combination of all
photometric probes, including cross correlations, WL+GCph+XCph (red) and the
combination of all spectroscopic and photometric probes GCsp+WL+GCph+XCph
(yellow). We keep the same color convention in all figures of the paper, when
showing constraints from different probes. Table 2 contains this information
also for the pessimistic survey settings.

In Fig. 4, we plot the elliptical $1\sigma$ and $2\sigma$ contours for the
probes $\text{GC}_{\mathrm{sp}}$, WL, the combination
$\text{GC}_{\mathrm{sp}}$+WL, and all the Euclid probes combined
GCsp+WL+GCph+XCph, using the same color convention mentioned before. From the
parameters used in the Fisher matrix, we leave out of this plot
$\Omega_{\mathrm{b},0}$ since it does not provide any extra insight. As it can
be seen from this figure, GCsp is always better at constraining $h$ and
$n_{\rm s}$, compared to the cosmic shear probe (WL) on its own. However, due
to the orthogonality of the contours, especially in the subspaces combining
$\sigma_{8}$ with $\Omega_{\mathrm{b},0}$, $n_{\rm s}$ and $h$, there is an
important lifting of degeneracies, which makes the combination of GCsp and WL
(shown in orange) much more constraining. The relative constraint on
$\operatorname{log_{10}}|f_{R0}|$ coming from WL alone (blue) is of the order
of $4.7\%$, this reduces by more than a factor of two when combining with the
spectroscopic probe (orange) and another factor two when adding all Euclid
probes together (yellow), yielding in total a relative constraint on
$\operatorname{log_{10}}|f_{R0}|$ of the order of $1.0\%$ in this optimistic
case. Throughout this work, we indicate the $1\sigma$ constraints rounded to
the nearest significant digit, since our Fisher matrix method has been
validated at the 10% level on the discrepancy between $1\sigma$ marginalized
and unmarginalized errors on the cosmological and nuisance parameters, as
mentioned in Sect. 4.

Contour ellipses comparing the constraints from the combination of
spectroscopic GC and WL to the 3x2pt photometric combination
(WL+$\text{GC}_{\mathrm{ph}}$+$\text{XC}_{\mathrm{ph}}$) can be found in Fig.
10 in the Appendix (using the same global color convention), where it can be
seen that it is the spectroscopic probe that helps to break degeneracies in
the $h$ and $\Omega_{\mathrm{b},0}$ planes, mainly due to the sensitivity of
the BAO wiggles on these two parameters. The cosmic shear probe (WL) is
relatively insensitive to $\Omega_{\mathrm{b},0}$ and $h$ and the full
combination of photometric probes is also not good at constraining $h$. It is
the breaking of degeneracies when combining GCsp and WL probes, that improves
considerably the constraints on all parameters, showcasing the power of Euclid
to measure parameters in and beyond the standard model of cosmology.

### Constraints on the fundamental model parameter $|f_{R0}|$

Note that we perform the Fisher matrix analysis on the parameter
$\operatorname{log_{10}}|f_{R0}|$, instead of directly on $|f_{R0}|$, since
for very small numbers and for large order of magnitude differences, the
Fisher matrix derivatives might become unstable (see, e.g., Camera et al.,
2018, Appendix A1). Therefore, it is recommended to have all the involved
parameters in the Fisher matrix being of the same order of magnitude. Since
the transformation between $\operatorname{log_{10}}|f_{R0}|$ and $|f_{R0}|$ is
non-linear and the parameter constraints are not small in some cases, we
cannot simply use a Jacobian transformation to convert between the Fisher
matrices in this case. Our assumption of Gaussianity is only true for the
logarithmic parametrization $\operatorname{log_{10}}|f_{R0}|$, therefore the
posterior contours for $|f_{R0}|$ will be non-Gaussian.

We can, nevertheless, obtain the fully marginalized constraints on $|f_{R0}|$
by transforming the $\operatorname{log_{10}}|f_{R0}|$ symmetric bounds

$\operatorname{log_{10}}|f_{R0}|^{(\pm)}=\operatorname{log_{10}}|f_{R0}|_{\rm
fid}\pm\sigma_{\operatorname{log_{10}}|f_{R0}|},$ (35)

into the upper and lower bounds for the linearly parameterized parameter
$|f_{R0}|$

${|f_{R0}|}^{(\pm)}=10^{\left(\operatorname{log_{10}}|f_{R0}|_{\rm
fid}\pm\sigma_{\operatorname{log_{10}}|f_{R0}|}\right)}=|f_{R0}|_{\rm
fid}\times 10^{\left(\pm\sigma_{\operatorname{log_{10}}|f_{R0}|}\right)}\,.$
(36)

This will yield asymmetrical errors in $|f_{R0}|$ since the upper and lower
bounds will be given by the exponentiation of the symmetrical $1\sigma$
bounds.

Using these formulas we can obtain the upper and lower $1\sigma$-bounds for
our fiducial parameter $|f_{R0}|$ for the different cases

  * •

$|f_{R0}|=(5.0^{+2.2}_{-1.5}\times 10^{-6})$ with spectroscopic GCsp alone;

  * •

$|f_{R0}|=(5.0^{+3.9}_{-2.2}\times 10^{-6})$ with WL alone;

  * •

$|f_{R0}|=(5.0^{+0.91}_{-0.77}\times 10^{-6})$ combining WL, GCph, and XCph

  * •

$|f_{R0}|=(5.0^{+0.62}_{-0.55}\times 10^{-6})$
with the combination GCsp+WL+GCph+XCph.

As one can clearly see from these numbers, the stronger the constraint on the
$\operatorname{log_{10}}|f_{R0}|$ parameter, the more symmetric the upper and
lower bounds on $|f_{R0}|$ become, simply due to the central limit theorem and
the fact that for a very peaked likelihood, a Gaussian approximation is always
possible around the maximum of the posterior distribution.

In order to visualize this, we take our Fisher matrices computed for each of
the cases and probes and assign them to a multivariate Gaussian distribution,
by using the Fisher matrix as the inverse covariance. We then sample from this
distribution and transform the samples of $\operatorname{log_{10}}|f_{R0}|$
into samples on $|f_{R0}|$ using the inverse of the logarithmic
transformation, see Eq. 36. Using this technique, we plot on the left side of
Fig. 7 the marginalized posterior 1-d and 2-d distributions for the original
model parameter $|f_{R0}|$ in the HS6 case, for the optimistic scenario, using
the same color convention as before. Compared to Fig. 4 we show now, for
illustration purposes, the 3x2pt combination of photometric probes
(WL+GCph+XCph), instead of the $\text{GC}_{\mathrm{sp}}$+WL combination, and
we focus only on the parameters $|f_{R0}|$, $\Omega_{\mathrm{m},0}$ and
$\sigma_{8}$. For the subspaces involving the $|f_{R0}|$ parameter, we see now
non-symmetrical posterior distributions, which are clearly non-Gaussian, and
noticing the marginalized 1-d posterior distributions in the top of each sub-
panel, we see that they have a longer tail on the upper side of the
constraints. This asymmetry is also due to the fact that using the logarithmic
parameterisation automatically prevents us from going into negative regions of
the parameter space for $|f_{R0}|$, therefore the probability distributions
are skewed towards the positive side. However, the stronger the constraints on
this parameter (for example in the case of the full combination of the Euclid
probes, in yellow), the better the Gaussian approximation.

Figure 6:  1 and 2$\sigma$ joint marginal error contours on the cosmological
parameters for a flat $f(R)$ model with $|f_{R0}|=5\times 10^{-6}$ (HS6) vs.
$|f_{R0}|=5\times 10^{-7}$ (HS7) in the optimistic scenario. Lighter colors
and solid contours correspond to HS6, while darker colors and dashed empty
contours correspond to HS7. In purple the spectroscopic GCsp, in blue the WL
probe, in yellow the photometric and spectroscopic probes combined together
(GCsp+WL+GCph+XCph). One can see that for certain parameter combinations like
$\Omega_{\mathrm{m},0}$ vs $h$ the contours are very similar, in size and
orientation, since the fiducial values for the standard cosmological
parameters in these models only differ by their value of $\sigma_{8}$.
However, for $\operatorname{log_{10}}|f_{R0}|$ in HS7 the ellipses become much
more elongated, due to this model being very close to $\Lambda\mathrm{CDM}$
and therefore, being less constrained in its extra model parameter compared to
HS6.

### Comparison of the HS5, HS6 and HS7 Hu-Sawicki models

Figure 7:  Left: Same as Fig. 4 (HS6) for a subset of the cosmological
parameters when transforming the contours back into the fundamental parameter
$|f_{R0}|$. While in the logarithmic case the marginalized Fisher contours are
simply Gaussians by definition, this is not the case for the contours on the
linearly parameterized case $|f_{R0}|$. For the probes
$\text{GC}_{\mathrm{sp}}$ (purple) and WL (blue) on their own which are less
constraining, the errors on $|f_{R0}|$ are very unsymmetric with heavy tails
on the right. However, the smaller the $1\sigma$ marginalized error on
$\operatorname{log_{10}}|f_{R0}|$ is, the better can the marginalized
posterior on $|f_{R0}|$ be approximated by a Gaussian. This can be well seen
for the full combination of Euclid probes (yellow). Right:  1-dimensional
fully-marginalized posteriors on the $|f_{R0}|$ parameter for
$|f_{R0}|=5\times 10^{-7}$ (HS7, dotted lines), $|f_{R0}|=5\times 10^{-6}$
(HS6, solid lines) and $|f_{R0}|=5\times 10^{-5}$ (HS5, dashed lines), after
applying the transformation in Eq. 36. Here we show the spectroscopic probe
$\text{GC}_{\mathrm{sp}}$ on its own in purple, and the 3x2pt combination of
photometric probes (WL+$\text{GC}_{\mathrm{ph}}$+$\text{XC}_{\mathrm{ph}}$) in
red, using the same color convention throughout the whole paper. The shaded
areas correspond to the $1\sigma$ confidence intervals, while the lines extend
up to $3\sigma$. We can see that all three models can be differentiated at
$1\sigma$ from each other, even using the spectroscopic probe on its own.
However, the larger tail of HS7 on the right side, makes it impossible to
discriminate this model from HS6 at $2\sigma$ and from HS5 at $3\sigma$. Using
the full combination of photometric probes, all models can be distinguished at
more than $3\sigma$ from each other.

Besides our baseline fiducial for $f(R)$ gravity, $|f_{R0}|=5\times 10^{-6}$,
we also perform forecasts on a second and third fiducial, namely
$|f_{R0}|=5\times 10^{-5}$ (HS5) and $|f_{R0}|=5\times 10^{-7}$ (HS7), the
latter is closer to $\Lambda\mathrm{CDM}$ in terms of perturbations and
background evolution. A closer limit to $\Lambda\mathrm{CDM}$ is currently not
possible using non-linear predictions, since our non-linear fitting formula
discussed in Sect. 3.3 cannot be extrapolated for values of the $f(R)$ model
parameter smaller than $|f_{R0}|=5\times 10^{-7}$. We also perform forecasts
on a flat $\Lambda\mathrm{CDM}$ model in order to have a comparison scenario
within the assumptions of our recipe for the large scale structure probes.

In Table 3 we report the $1\sigma$ fully marginalized forecasted errors for
these models in the optimistic and pessimistic scenarios for the most
important Euclid probe-combinations discussed above. For the HS5 model we find
in general very similar constraints to HS6, but it is important to remark that
due to the increased signal at non-linear scales and also the larger
$k$-dependence of the growth, due to a higher value of $|f_{R0}|$, we obtain
better constraints on $\operatorname{log_{10}}|f_{R0}|$ for the full
combination of probes, both in the pessimistic and optimistic scenarios, with
relative errors of $1.4\%$ and $0.6\%$, respectively. In the case of the HS7
model, the fiducial value of $|f_{R0}|$ is yet another order of magnitude
closer to $\Lambda\mathrm{CDM}$ as compared to HS6 and therefore, we expect
the extra signal on the power spectrum coming from modified gravity to be much
smaller. We indeed find that this is the case, and compared to HS6, the
relative percentage constraints on $\operatorname{log_{10}}|f_{R0}|$ increase
by a factor 4.8 in GCsp, a factor 3.5 in WL and a factor 2.2 when considering
the full Euclid spectroscopic and photometric combination. The final relative
constraint on $\operatorname{log_{10}}|f_{R0}|$ for this latter case is
$1.8\%$, which corresponds to a determination of
$|f_{R0}|=(5.0^{+1.5}_{-1.1}\times 10^{-7})$. In Fig. 6 we plot the confidence
contours for this model HS7 in the optimistic scenario compared to the HS6
baseline model. In dark blue/blue we plot the constraints for GCsp for
HS7/HS6, in dark purple/violet the constraints on WL for HS7/HS6 and in
yellow/orange the constraints on the full spectroscopic+photometric Euclid
combination. For the standard cosmological parameters, such as in the subspace
of $\Omega_{\mathrm{m},0}$ vs. $h$ or $\sigma_{8}$, we can see that we recover
very similar contours in size and in orientation. On the other hand, one can
see in the panel of $\operatorname{log_{10}}|f_{R0}|$ vs. $\sigma_{8}$ that
the GCsp and WL probes alone are not sufficient to distinguish these two
models from each other, but the full spectroscopic+photometric combination is
able to discriminate between them at several $\sigma$. As discussed
previously, this is under the assumption that both models have the same
primordial amplitude and, therefore, due to a larger growth in the case of
HS6, they end up with different values of $\sigma_{8}$ today. To be certain
that the choice of $\sigma_{8}$ values does not affect our constraints, we
performed a Fisher matrix forecast of the HS6 model, changing its primordial
amplitude $A_{\rm s}$, so that it has the same amplitude of fluctuations
$\sigma_{8}=0.816$ today as the baseline $\Lambda\mathrm{CDM}$ model. We found
that none of the parameter constraints changed considerably with respect to
the baseline HS6 case.

We have also tested that when fixing the extra model parameter
$\operatorname{log_{10}}|f_{R0}|$ in the HS7 model, which is the closest we
have to $\Lambda\mathrm{CDM}$ and just deviates in the non-linear power
spectrum and the growth rate function by a few percent (see Figs. 2 and 3), we
then recover the same constraints on the five remaining $\Lambda\mathrm{CDM}$
parameters. This comparison is shown in the right panel of Fig. 10, where in
dark yellow we show the full combination of Euclid probes for the pessimistic
HS7 case and in cyan the corresponding 1-d and 2-d contours for the
$\Lambda\mathrm{CDM}$ scenario as described in EC19.

On the right panel of Fig. 7, we compare now the constraints across the three
models studied in this work, HS5 (in dashed lines), HS6 (in solid lines) and
HS7 (in dotted lines) for the fundamental parameter $|f_{R0}|$. We plot the
1-d fully-marginalized posterior on $|f_{R0}|$ for the spectroscopic probe
$\text{GC}_{\mathrm{sp}}$ on its own and the 3x2pt photometric combination of
probes, using again the same color codes as before. The shaded areas under the
curve represent the $1\sigma$ confindence intervals and the extent of the
lines corresponds to a $3\sigma$ deviation from the mean. Due to the
logarithmic scale of the $x$-axis, the distributions look symmetric again,
nevertheless we can see that in absolute $|f_{R0}|$ numbers, the constraints
are much larger on the right side of the constraints. In the HS7 model for the
optimistic scenario, the forecasted upper bound on $|f_{R0}|$ would be
$2.9\times 10^{-6}$ at $1\sigma$, $1.7\times 10^{-5}$ at $2\sigma$ and
$9.6\times 10^{-5}$ at $3\sigma$. For instance, assuming data would follow the
HS7 model, using only the spectroscopic probe, it would not be possible to
distinguish it from HS6 at the $2\sigma$ confidence level. Discriminating HS7
from HS5 would be impossible at a $3\sigma$ confidence level. Logically, this
situation degrades even more in the pessimistic scenario. However, when using
the full combination of photometric probes (red lines and shaded regions), all
three models would be distinguishable from each other at more than $3\sigma$.
The 1-, 2- and $3\sigma$ level upper bounds for HS7 are in this combination of
probes of the order of $8.4\times 10^{-7}$, $1.4\times 10^{-6}$ and $2.4\times
10^{-6}$, respectively. Using Eq. 36, the rest of the non-symmetrical bounds
on $|f_{R0}|$ can be recovered for all the numbers listed in Tables 2 and 3.
Since we performed the Fisher matrix analysis on the logarithm of the original
model parameter, we opt to report here only the forecasted Fisher matrix
constraints which are Gaussian and therefore, symmetric.

Under our given assumptions (and the fact that we have fixed for both models
the sum of neutrino masses and the same initial amplitude of perturbations
$A_{\rm s}$), we can say that the full combination of Euclid probes would be
capable of distinguishing between these three $f(R)$ models with a high degree
of certainty. Since most of the constraints in the photometric probe come from
small scales, this highlights again the importance of using the deeply non-
linear regime of structure formation to test this $f(R)$ theory accurately and
in an ideal scenario, discern it from the $\Lambda\mathrm{CDM}$ model,
provided the data indicates a strong preference for the standard cosmological
model with a large statistical significance.

### Pessimistic survey scenario

For the baseline HS6 model and in the pessimistic scenario, the errors on
$\operatorname{log_{10}}|f_{R0}|$, using GCsp, do not degrade that much,
increasing to $3.6\%$ when scales of $k_{\rm max}=0.25\,h\,{\rm Mpc}^{-1}$ are
considered, while increasing to $6.3\%$ when the cut is performed at a smaller
wavenumber $k_{\rm max}=0.15\,h\,{\rm Mpc}^{-1}$, which we label as the quasi-
linear scenario in this work. On the other hand, for photometric probes the
pessimistic scenario has a much larger impact on the errors, because scales
are decreased from $\ell_{\rm max}=5000$ to $1500$ in the case of WL and from
$\ell_{\rm max}=3000$ to $750$ in the case of GCph (see Table 1). For WL
alone, the errors degrade by a factor 1.7, for the combination of all
photometric probes, the errors become 2 times larger, and since these probes
dominate the total constraining power on this parameter, for the combination
of spectroscopic and photometric probes, the errors on
$\operatorname{log_{10}}|f_{R0}|$ are a factor of 1.7 larger than compared to
the optimistic scenario. For HS5 the degradation between the optimistic and
pessimistic scenarios are of the same order, of about a factor 2, while for
HS7 the degradation is stronger and it is of a factor 3, when photometric
probes are involved.

Our tests demonstrate that the non-linear scales have a much larger impact on
the photometric probes, since in order to probe multipoles of about
$\ell=5000$, one needs to have an accurate determination of the power spectrum
in the deeply non-linear regime, at least to scales of about $k=5\,h\,{\rm
Mpc}^{-1}$ at $z=0$, before noise and other systematic effects, like baryonic
contributions start dominating (for an analysis of the effect of $k$-cuts on
the information content in cosmic shear, see Taylor et al., 2018). In the
right panel of Fig. 9, we compare the Fisher matrix contours for the
optimistic scenario (yellow) versus the pessimistic scenario (blue), within
our baseline model and using all Euclid primary probes combined.

Table 3: Forecast $1\sigma$ marginal relative errors on the cosmological
parameters for a flat $f(R)$ model with $|f_{R0}|=5\times 10^{-5}$ ,
$|f_{R0}|=5\times 10^{-7}$ and the flat $\Lambda\mathrm{CDM}$ model in the
pessimistic and optimistic scenarios. $|f_{R0}|=5\times 10^{-5}$
---
| $\Omega_{\rm m,0}$ | $\Omega_{\rm b,0}$ | $\operatorname{log_{10}}|f_{R0}|$ | $h$ | $n_{\rm s}$ | $\sigma_{8}$
Fiducial values | $0.32$ | $0.05$ | $-4.301$ | $0.67$ | $0.96$ | $0.911$
Pessimistic setting
GCsp | 1.5% | 2.9% | 5.2% | 0.4% | 1.3% | 1.0%
WL+GCph+XCph | 0.8% | 4.9% | 1.9% | 2.8% | 0.9% | 0.5%
GCsp+WL+GCph+XCph | 0.6% | 1.6% | 1.4% | 0.2% | 0.5% | 0.3%
Optimistic setting
GCsp | 1.5% | 2.7% | 5.0% | 0.3% | 1.2% | 0.9%
WL+GCph+XCph | 0.3% | 4.3% | 0.9% | 2.1% | 0.5% | 0.1%
GCsp+WL+GCph+XCph | 0.2% | 1.2% | 0.6% | 0.1% | 0.2% | 0.1%
$|f_{R0}|=5\times 10^{-7}$
| $\Omega_{\rm m,0}$ | $\Omega_{\rm b,0}$ | $\operatorname{log_{10}}|f_{R0}|$ | $h$ | $n_{\rm s}$ | $\sigma_{8}$
Fiducial values | $0.32$ | $0.05$ | $-6.301$ | $0.67$ | $0.96$ | $0.823$
Pessimistic setting
GCsp | 1.4% | 2.8% | 14% | 0.4% | 1.0% | 1.8%
WL+GCph+XCph | 0.8% | 5.4% | 8.9% | 4.0% | 2.0% | 0.7%
GCsp+WL+GCph+XCph | 0.6% | 1.6% | 5.4% | 0.2% | 0.6% | 0.4%
Optimistic setting
GCsp | 1.3% | 2.4% | 12% | 0.3% | 0.9% | 1.5%
WL+GCph+XCph | 0.6% | 4.3% | 3.6% | 2.3% | 0.7% | 0.5%
GCsp+WL+GCph+XCph | 0.3% | 1.3% | 1.8% | 0.1% | 0.2% | 0.2%
flat $\Lambda\mathrm{CDM}$
| $\Omega_{\rm m,0}$ | $\Omega_{\rm b,0}$ | - | $h$ | $n_{\rm s}$ | $\sigma_{8}$
Fiducial values | $0.32$ | $0.05$ | - | $0.67$ | $0.96$ | $0.816$
Pessimistic setting
GCsp | 1.3% | 2.1% | - | 0.3% | 0.9% | 0.9%
WL+GCph+XCph | 0.8% | 5.2% | - | 2.7% | 0.9% | 0.4%
GCsp+WL+GCph+XCph | 0.5% | 1.3% | - | 0.2% | 0.5% | 0.2%
Optimistic setting
GCsp | 1.2% | 1.9% | - | 0.2% | 0.8% | 0.8%
WL+GCph+XCph | 0.3% | 4.6% | - | 2.0% | 0.4% | 0.1%
GCsp+WL+GCph+XCph | 0.2% | 0.9% | - | 0.1% | 0.1% | 0.1%

##  6 Conclusions

In this work we have focused on the ability of the future Euclid mission to
constrain extensions of the concordance cosmological model. The large
literature on $f(R)$ cosmologies, as well as their properties (screening
mechanism, scale-dependent growth of structures, possibility to rewrite it as
a scalar-tensor theory embedded in the Horndeski action), make $f(R)$ an ideal
test case of modified gravity models. In this work we have considered Hu &
Sawicki (2007) $f(R)$ cosmologies and predicted forecasts on their peculiar
parameter, $f_{R0}$, and other cosmological parameters.

As done in EC19, forecasts are validated and computed for spectroscopic galaxy
clustering (GCsp), photometric galaxy clustering (GCph), weak lensing cosmic
shear (WL), and the cross-correlation between the last two probes (XCph). The
extension presented in this paper over the forecasting pipeline developed in
EC19 is threefold: first, for the class of cosmologies investigated here,
which goes beyond $\Lambda\mathrm{CDM}$; secondly, in the encoded equations
for all probes, as they now include a growth factor which is scale dependent;
and finally, in the software development of Fisher matrix codes, which are now
validated also for this extended recipe. Validation was pursued as follows:

  1. 1.

We have first compared the input quantities obtained from different Boltzmann
solvers and compared the quasi-static limit with the full evolution, finding
equivalent results for the redshifts and scales of interest—in particular,
predictions on the matter power spectrum up to $k=10\,h\,\mathrm{Mpc}^{-1}$
lead to sub-percent agreement between quasi-static and full evolution.

  2. 2.

We have then used the forecasting pipeline presented in EC19 and modified it
to account for the scale-dependent quantities and the other MG modifications
arising in $f(R)$.

  3. 3.

As done in EC19, we have implemented independent forecast codes and considered
their output validated when the discrepancy on each parameter uncertainty,
compared to the median, was less than $10\%$, both before and after
marginalisation.

Hence, we have computed forecasts in an optimistic and a pessimistic scenario,
depending on the range in scales considered for each probe and on the level of
systematics included. This choice follows the one in EC19, for easier
comparison, and it is summarised in Table 1. In the specific case of GCsp, we
have also shown the impact of a quasi-linear, more conservative choice, due to
the additional uncertainty in the non-linear predictions for a theoretical
model beyond $\Lambda\mathrm{CDM}$.

In an optimistic scenario, combining all Euclid primary probes, and
considering our baseline fiducial $|f_{R0}|=5\times 10^{-6}$, we have obtained
that Euclid alone will be able to constrain $\operatorname{log_{10}}|f_{R0}|$
at the $1.0$% level using the full combination GCsp+WL+GCph+XCph. This
uncertainty increases to $1.4$%, when considering only the photometric probes,
which highlights the importance of combining spectroscopic and photometric
observables. The constraint considering the spectroscopic probe alone
increases to $3\%$. We have discussed in the text how these constraints on the
logarithmic parameter transform into asymmetrical constraints of the original
model parameter $|f_{R0}|$. The $\operatorname{log_{10}}|f_{R0}|$
parametrization has been used several times in the literature, not only for
theoretical considerations but also in Markov-Chain-Monte-Carlo explorations
of cosmological likelihoods, such as in Dossett et al. (2014); Schneider et
al. (2020) and Tröster et al. (2021). However, for future analysis with real
data, one should reconsider if other compact re-parametrizations of
$|f_{R0}|$, for example the one used in Terukina et al. (2014) would be
physically more meaningful.

The reported bounds in Table 2 of the optimistic HS6 scenario, correspond to
determinations of $|f_{R0}|$ that are $|f_{R0}|=(5.0^{+2.2}_{-1.5}\times
10^{-6})$ with spectroscopic GCsp alone, $|f_{R0}|=(5.0^{+0.91}_{-0.77}\times
10^{-6})$ combining WL, GCph, and XCph and
$|f_{R0}|=(5.0^{+0.62}_{-0.55}\times 10^{-6})$ with the full Euclid
combination GCsp+WL+GCph+XCph. For the fiducial model HS5, $|f_{R0}|=5\times
10^{-5}$, which is already at the margins of the possible observational
constraints, we have obtained constraints of the same order. This hints at the
fact that Euclid should be able to distinguish these two models from each
other once data becomes available. For the model HS7 we find that the bounds
get significantly less constraining, since it is a model very close to
$\Lambda\mathrm{CDM}$ with a similar growth of perturbations as compared to
the standard cosmological model. Nevertheless we also find that under the
optimistic scenario and using the full combination of Euclid photometric and
spectroscopic probes, this model can as well be distinguished from HS6 and
from a $\Lambda\mathrm{CDM}$ cosmology.

This highlights the significant amount of information that can be extracted
from cross-correlating photometric probes, and by including a larger range in
non-linear scales. In conclusion, Euclid will be able to provide outstanding
constraints on extensions beyond the concordance model thanks to its ability
to probe the non-linear scales, where most of the information is contained.
However, a good knowledge of the modelling of our theoretical observables at
these scales, like the one used in this analysis, and a good control of the
systematic uncertainties is required in order to reach reliable final results
when data will be available.

###### Acknowledgements.

We thank Bill Wright for useful comments on the degeneracy between $f(R)$
models and neutrino masses. S. Cas. thanks Sabarish V. M. and Sefa Pamuk for
improvements in the development of the code CosmicFish used in this work. S.
Cas. acknowledges support from CNRS and ‘Centre National d’Études Spatiales’
(CNES) grants during the early stages of this work. DS acknowledges financial
support from the Fondecyt Regular project number 1200171. IT 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. S. Cam. acknowledges support from the
‘Departments of Excellence 2018-2022’ Grant (L. 232/2016) awarded by the
Italian Ministry of Education, University and Research (miur), as well as
support by miur Rita Levi Montalcini project ‘prometheus – Probing and
Relating Observables with Multi-wavelength Experiments To Help Enlightening
the Universe’s Structure’ for the early stages of this project. P. N. is
financially supported by ‘Centre National d’Études Spatiales’ (CNES). The
research of N. F. is supported by the Italian Ministry of University and
Research (miur) through the Rita Levi Montalcini project “Tests of gravity on
cosmic scales” with reference PGR19ILFGP. N. F. and F. P. acknowledges the FCT
project with ref. number PTDC/FIS-AST/0054/2021 F. P. acknowledges support
from the INFN grant InDark and the Departments of Excellence grant L.232/2016
of the Italian Ministry of Education, University and Research (miur). K. K. is
supported by the UK STFC grant ST/S000550/1, and the European Research Council
under the European Union’s Horizon 2020 programme (grant agreement No. 646702
“CosTesGrav”). V. Y. acknowledges funding from the European Research Council
(ERC) under the European Union’s Horizon 2020 research and innovation
programme (grant agreement No. 769130). L. L. was supported by a Swiss
National Science Foundation (SNSF) Professorship grant (No. 170547). A. P. is
a UK Research and Innovation Future Leaders Fellow [grant MR/S016066/2]. We
acknowledge open libraries support IPython (Perez & Granger, 2007), Matplotlib
(Hunter, 2007), Numpy (Walt et al., 2011), and SciPy (Virtanen et al., 2020).
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).

## References

  * Abbott et al. (2023) Abbott, T. M. C. et al. 2023, Phys. Rev. D, 107, 083504
  * Amendola et al. (2008) Amendola, L., Kunz, M., & Sapone, D. 2008, JCAP, 04, 013
  * Arnold et al. (2021) Arnold, C., Li, B., Giblin, B., Harnois-Déraps, J., & Cai, Y.-C. 2021, arXiv preprint arXiv:2109.04984
  * Baker et al. (2019) Baker, T., Barreira, A., Desmond, H., et al. 2019 [arXiv:1908.03430]
  * Baker & Bull (2015) Baker, T. & Bull, P. 2015, ApJ, 811, 116
  * Baldi et al. (2014) Baldi, M., Villaescusa-Navarro, F., Viel, M., et al. 2014, MNRAS, 440, 75
  * Bellini et al. (2018) Bellini, E. et al. 2018, Phys. Rev. D, 97, 023520
  * Beutler et al. (2014) Beutler, F., Saito, S., Seo, H.-J., et al. 2014, MNRAS, 443, 1065
  * Beutler et al. (2017) Beutler, F., Seo, H.-J., Saito, S., et al. 2017, MNRAS, 466, 2242
  * Blas et al. (2011) Blas, D., Lesgourgues, J., & Tram, T. 2011, JCAP, 07, 034
  * Boyle & Komatsu (2018) Boyle, A. & Komatsu, E. 2018, JCAP, 03, 035
  * Brax & Valageas (2013) Brax, P. & Valageas, P. 2013, Phys. Rev. D, 88, 023527
  * Buchdahl (1970) Buchdahl, H. A. 1970, MNRAS, 150, 1
  * Bueno Belloso et al. (2011) Bueno Belloso, A., Garcia-Bellido, J., & Sapone, D. 2011, JCAP, 10, 010
  * Cabre et al. (2012) Cabre, A., Vikram, V., Zhao, G.-B., Jain, B., & Koyama, K. 2012, JCAP, 07, 034
  * Camera et al. (2018) Camera, S., Fonseca, J., Maartens, R., & Santos, M. G. 2018, MNRAS, 481, 1251
  * Casas et al. (2023) Casas, S. et al. 2023 [arXiv:2303.09451]
  * Cautun et al. (2018) Cautun, M., Paillas, E., Cai, Y.-C., et al. 2018, MNRAS, 476, 3195
  * Colas et al. (2020) Colas, T., D’amico, G., Senatore, L., Zhang, P., & Beutler, F. 2020, JCAP, 06, 001
  * Contarini et al. (2021) Contarini, S., Marulli, F., Moscardini, L., et al. 2021, MNRAS, 504, 5021
  * Costille et al. (2018) Costille, A., Caillat, A., Rossin, C., et al. 2018, in Space Telescopes and Instrumentation 2018: Optical, Infrared, and Millimeter Wave, ed. M. Lystrup, H. A. MacEwen, G. G. Fazio, N. Batalha, N. Siegler, & E. C. Tong, Vol. 10698, International Society for Optics and Photonics (SPIE), 730
  * Cropper et al. (2018) Cropper, M., Pottinger, S., Azzollini, R., et al. 2018, in Space Telescopes and Instrumentation 2018: Optical, Infrared, and Millimeter Wave, ed. M. Lystrup, H. A. MacEwen, G. G. Fazio, N. Batalha, N. Siegler, & E. C. Tong, Vol. 10698, International Society for Optics and Photonics (SPIE), 709
  * de Mattia et al. (2021) de Mattia, A., Ruhlmann-Kleider, V., Raichoor, A., et al. 2021, MNRAS, 501, 5616
  * Desmond & Ferreira (2020) Desmond, H. & Ferreira, P. G. 2020, Phys. Rev. D, 102, 104060
  * Dossett et al. (2014) Dossett, J., Hu, B., & Parkinson, D. 2014, JCAP, 03, 046
  * Eisenstein & Hu (1998) Eisenstein, D. J. & Hu, W. 1998, ApJ, 496, 605
  * Euclid Collaboration: Blanchard et al. (2020) Euclid Collaboration: Blanchard, A., Camera, S., Carbone, C., et al. 2020, A&A, 642, A191
  * Euclid Collaboration: Knabenhans et al. (2021) Euclid Collaboration: Knabenhans, Knabenhans, M., Stadel, J., et al. 2021, MNRAS, 505, 2840
  * Euclid Collaboration: Knabenhans et al. (2019) Euclid Collaboration: Knabenhans, M. et al. 2019, MNRAS, 484, 5509
  * García-Farieta et al. (2019) García-Farieta, J. E., Marulli, F., Veropalumbo, A., et al. 2019, MNRAS, 488, 1987
  * Giocoli et al. (2018) Giocoli, C., Baldi, M., & Moscardini, L. 2018, MNRAS, 481, 2813
  * Hagstotz et al. (2019a) Hagstotz, S., Costanzi, M., Baldi, M., & Weller, J. 2019a, MNRAS, 486, 3927
  * Hagstotz et al. (2019b) Hagstotz, S., Gronke, M., Mota, D. F., & Baldi, M. 2019b, A&A, 629, A46
  * Harnois-Déraps et al. (2015) Harnois-Déraps, J., Munshi, D., Valageas, P., et al. 2015, MNRAS, 454, 2722
  * He (2013) He, J.-h. 2013, Phys. Rev. D, 88, 103523
  * Hojjati et al. (2016) Hojjati, A., Plahn, A., Zucca, A., et al. 2016, Phys. Rev. D, 93, 043531
  * Hojjati et al. (2011) Hojjati, A., Pogosian, L., & Zhao, G.-B. 2011, JCAP, 1108, 005
  * Hu et al. (2014) Hu, B., Raveri, M., Frusciante, N., & Silvestri, A. 2014, Phys. Rev. D, 89, 103530
  * Hu et al. (2016) Hu, B., Raveri, M., Rizzato, M., & Silvestri, A. 2016, MNRAS, 459, 3880
  * Hu et al. (2015) Hu, B., Raveri, M., Silvestri, A., & Frusciante, N. 2015, Phys. Rev. D, 91, 063524
  * Hu & Sawicki (2007) Hu, W. & Sawicki, I. 2007, Phys. Rev., D76, 064004
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * Khoury & Weltman (2004) Khoury, J. & Weltman, A. 2004, Phys. Rev. Lett., 93, 171104
  * Laureijs et al. (2011) Laureijs, R., Amiaux, J., Arduini, S., et al. 2011 [arXiv:1110.3193]
  * Lesgourgues (2011) Lesgourgues, J. 2011 [arXiv:1104.2932]
  * Lewis et al. (2000) Lewis, A., Challinor, A., & Lasenby, A. 2000, ApJ, 538, 473
  * Lombriser (2014) Lombriser, L. 2014, Annalen Phys., 526, 259
  * Ma & Bertschinger (1995) Ma, C.-P. & Bertschinger, E. 1995, ApJ, 455, 7
  * Markovic et al. (2019) Markovic, K., Bose, B., & Pourtsidou, A. 2019, The Open Journal of Astrophysics, 2, 13
  * Merten et al. (2019) Merten, J., Giocoli, C., Baldi, M., et al. 2019, MNRAS, 487, 104
  * Motohashi et al. (2013) Motohashi, H., Starobinsky, A. A., & Yokoyama, J. 2013, Phys. Rev. Lett., 110, 121302
  * Nunes et al. (2017) Nunes, R. C., Pan, S., Saridakis, E. N., & Abreu, E. M. C. 2017, JCAP, 01, 005
  * Okada et al. (2013) Okada, H., Totani, T., & Tsujikawa, S. 2013, Phys. Rev. D, 87, 103002
  * Peel et al. (2019) Peel, A., Lalande, F., Starck, J.-L., et al. 2019, Phys. Rev. D, 100
  * Peel et al. (2018) Peel, A., Pettorino, V., Giocoli, C., Starck, J.-L., & Baldi, M. 2018, A&A, 619, A38
  * Peloso et al. (2015) Peloso, M., Pietroni, M., Viel, M., & Villaescusa-Navarro, F. 2015, JCAP, 1507, 001
  * Perez & Granger (2007) Perez, F. & Granger, B. E. 2007, Computing in Science Engineering, 9, 21
  * Pérez-Romero & Nesseris (2018) Pérez-Romero, J. & Nesseris, S. 2018, Phys. Rev. D, 97, 023525
  * Planck Collaboration (2016) Planck Collaboration. 2016, A&A, 594, A14
  * Pogosian & Silvestri (2008) Pogosian, L. & Silvestri, A. 2008, Phys. Rev. D, 77, 023503, [Erratum: Phys.Rev.D 81, 049901 (2010)]
  * Potter et al. (2017) Potter, D., Stadel, J., & Teyssier, R. 2017, Comput. Astrophys. Cosmol., 4, 2
  * Raveri et al. (2014) Raveri, M., Hu, B., Frusciante, N., & Silvestri, A. 2014, Phys. Rev. D, 90, 043513
  * Ryu et al. (2020) Ryu, S., Lee, J., & Baldi, M. 2020, ApJ, 904, 93
  * Sakr & Martinelli (2022) Sakr, Z. & Martinelli, M. 2022, JCAP, 05, 030
  * Schneider et al. (2020) Schneider, A., Refregier, A., Grandis, S., et al. 2020, JCAP, 04, 020
  * Schneider et al. (2020) Schneider, A., Stoira, N., Refregier, A., et al. 2020, JCAP, 2020, 019
  * Shim et al. (2014) Shim, J., Lee, J., & Baldi, M. 2014 [arXiv:1404.3639]
  * Song et al. (2007) Song, Y.-S., Hu, W., & Sawicki, I. 2007, Phys. Rev. D, 75, 044004
  * Song et al. (2015) Song, Y.-S., Taruya, A., Linder, E., et al. 2015, Phys. Rev. D, 92, 043522
  * Spurio Mancini et al. (2019) Spurio Mancini, A., Köhlinger, F., Joachimi, B., et al. 2019, MNRAS, 490, 2155
  * Takahashi et al. (2012) Takahashi, R., Sato, M., Nishimichi, T., Taruya, A., & Oguri, M. 2012, ApJ, 761, 152
  * Taruya et al. (2010) Taruya, A., Nishimichi, T., & Saito, S. 2010, Phys. Rev. D, 82, 063522
  * Taylor et al. (2018) Taylor, P. L., Bernardeau, F., & Kitching, T. D. 2018, Phys. Rev. D, 98, 083514
  * Taylor & Markovič (2022) Taylor, P. L. & Markovič, K. 2022, Phys. Rev. D, 106, 063536
  * Terukina et al. (2014) Terukina, A., Lombriser, L., Yamamoto, K., et al. 2014, JCAP, 04, 013
  * Tröster et al. (2021) Tröster, T. et al. 2021, Astron. Astrophys., 649, A88
  * Tutusaus et al. (2020) Tutusaus, I., Martinelli, M., Cardone, V. F., et al. 2020, A&A, 643, A70
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Walt et al. (2011) Walt, S. v. d., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engg., 13, 22
  * Winther et al. (2019) Winther, H., Casas, S., Baldi, M., et al. 2019, Phys. Rev. D, 100, 123540
  * Winther et al. (2015) Winther, H. A. et al. 2015, MNRAS, 454, 4208
  * Wright et al. (2019) Wright, B. S., Koyama, K., Winther, H. A., & Zhao, G.-B. 2019, JCAP, 2019, 040–040
  * Zhang (2006) Zhang, P. 2006, Phys. Rev., D73, 123504
  * Zhang et al. (2007) Zhang, P., Liguori, M., Bean, R., & Dodelson, S. 2007, Phys. Rev. Lett., 99, 141302
  * Zhao et al. (2009) Zhao, G.-B., Pogosian, L., Silvestri, A., & Zylberberg, J. 2009, Phys. Rev., D79, 083513
  * Zucca et al. (2019) Zucca, A., Pogosian, L., Silvestri, A., & Zhao, G. B. 2019, J. Cosmology Astropart. Phys., 2019, 001

##  Appendix A Further results

Figure 8: Comparison between the Eisenstein-Hu fitting method (in purple) and
the Savitzky-Golay filter (in orange) for calculating the non-wiggle matter
power spectrum. The curves show the ratio of the total linear matter power
spectrum from CAMB divided by the corresponding smooth non-wiggle power
spectrum.

In this Appendix we show further results to highlight different aspects of our
analysis. In Fig. 8 we show the difference between our determination of the
non-wiggle matter power spectrum using the Savitzky-Golay (SG) filter (in
blue) as compared to the standard Eisenstein-Hu method (in red) used
previously in EC19. As discussed in Sect. 3.2 the SG method performs well and
it is easily generalizable to modified gravity theories and models in which
there is a scale-dependent growth of perturbations.

In Fig. 9 we show the comparison between the optimistic and pessimistic
scenarios of the Euclid probes as described in Table 1, for the HS6 model. In
the left panel we show for the spectroscopic GC probe the pessimistic quasi-
linear (dotted-dashed lines), the pessimistic (green) and the optimistic
(purple) cases. One can see that, as expected, the difference between the
pessimistic scale-cut at $0.25\,h\,{\rm Mpc}^{-1}$ and the optimistic one at
$0.3\,h\,{\rm Mpc}^{-1}$, is not that visible in the final contours. However,
the quasi-linear case ($k_{\rm max}=0.15\,h\,{\rm Mpc}^{-1}$), does discard
important information in the quasi-linear regime and this has a major impact
on the contours, degrading them by a factor 2 or more. In the right panel of
Fig. 9 we show the pessimistic (in cyan, empty contours) and optimistic (in
orange, filled contours) scenarios for the full combination of photometric
probes. As discussed in Sect. 5 and shown in Table 2 the pessimistic errors on
$\operatorname{log_{10}}|f_{R0}|$ are about a factor 2 to 3 larger than in the
optimistic scenario, mainly due to the fact that non-linear scales contain
crucial information on the scale-dependent growth of perturbations and the
additional clustering caused by the fifth force present in $f(R)$. This
highlights the importance in the future for real-data analysis to model these
scales correctly and accurately to avoid introducing systematic errors in the
analysis that could severely bias the parameter determination.

To study the effect of the probe combinations in more detail, we show in the
left panel of Fig. 10, for the $|f_{R0}|=5\times 10^{-6}$ model (HS6) and in
the optimistic case, the main combinations possible with Euclid primary
probes. Using the same color convention of the paper, we show in orange the
elliptical 1$\sigma$ error contours for GCsp+WL only, and in red the effect of
the full combination of photometric probes (including its cross-correlation,
$\text{XC}_{\mathrm{ph}}$). This is under the assumption that we can neglect
the cross-correlation between the spectroscopic and photometric probes from
the Euclid observations, which has been shown to be a good assumption for
Stage-IV surveys (see Taylor & Markovič (2022) and Paganin 2023, in prep.).
For the $h$ and the $\Omega_{\mathrm{b},0}$ parameter, as discussed in Sect.
5, $\text{GC}_{\mathrm{sp}}$ is the probe that gives the tighter constraints
and the breaking of degeneracies in these parameter combinations, and it is
what yields the improved constraints on $\operatorname{log_{10}}|f_{R0}|$ when
combining GCsp with the photometric probes.

Finally, in the right panel of Fig. 10 we show in dark yellow the fully
marginalized contours for the HS7 model and in dashed cyan the contours for a
$\Lambda\mathrm{CDM}$ equivalent, after fixing (maximizing) the extra
parameter $\operatorname{log_{10}}|f_{R0}|$. We do this to prove that the HS7
model is close enough to $\Lambda\mathrm{CDM}$ and that we recover very
similar constraints when comparing the same number of free parameters, even if
the HS7 model has a larger $\sigma_{8}$. In Table 3 the full list of $1\sigma$
marginalized errors can be found for HS5, HS7 and the $\Lambda\mathrm{CDM}$
counterpart.

Figure 9: Left panel: 1 and 2$\sigma$ joint marginal error contours on the
cosmological parameters for a flat $f(R)$ model with $|f_{R0}|=5\times
10^{-6}$ (HS6) in the optimistic case for the GC spectroscopic probe. In black
dot-dashed lines the quasi-linear case ($k_{\rm max}=0.15$), in solid green
the pessimistic case ($k_{\rm max}=0.25$) and in purple the optimistic case
($k_{\rm max}=0.3$). Right panel: The full combination GCsp+WL+GCph+XCphfor
the HS6 case, comparing the optimistic (yellow, solid) and pessimistic (cyan,
empty) specifications. See Table 1 for details on these scenarios.

Figure 10: Left panel: 1 and 2$\sigma$ joint marginal error contours on the
cosmological parameters for a flat $f(R)$ model with $|f_{R0}|=5\times
10^{-6}$ (HS6) in the optimistic case. In orange the combination of the
spectroscopic probe GCsp and cosmic shear WL (same color code as in Fig. 5),
in red the photometric probes with their cross-correlation (WL+GCph+XCph). As
explained in the main text, it is only the spectroscopic probe that is capable
to break degeneracies with the Hubble parameter $h$ and
$\Omega_{\mathrm{b},0}$, such that the combination of spectroscopic and
photometric yields indeed much better constraints on the model parameter
$\operatorname{log_{10}}|f_{R0}|$. The full combination of the spectroscopic
and all the photometric probes including its cross-correlation is shown in
Fig. 4. Right panel: Comparison of the posterior contours when fixing (i.e.
maximizing) the $\operatorname{log_{10}}|f_{R0}|$ parameter from the Fisher
matrices for the full combination GCsp+WL+GCph+XCph for the HS7 case. In dark
yellow the HS7 case and in cyan dashed lines the GR baseline case. The contour
shapes match very well in all parameter subspaces, considering that in
$\sigma_{8}$ their fiducial values are different. This is expected from the
fact that the HS7 model is very close to $\Lambda\mathrm{CDM}$.