Euclid preparation. Simulations and nonlinearities beyond $Λ$CDM. 4. Constraints on $f(R)$ models from the photometric primary probes
Euclid Collaboration K. Koyama S. Pamuk S. Casas B. Bose P. Carrilho I. Sáez-Casares L. Atayde M. Cataneo B. Fiorini C. Giocoli A. M. C. Le Brun F. Pace A. Pourtsidou Y. Rasera Z. Sakr H. -A. Winther E. Altamura J. Adamek M. Baldi M. -A. Breton G. Rácz F. Vernizzi A. Amara S. Andreon N. Auricchio C. Baccigalupi S. Bardelli F. Bernardeau C. Bodendorf D. Bonino E. Branchini M. Brescia J. Brinchmann A. Caillat S. Camera V. Capobianco C. Carbone J. Carretero M. Castellano G. Castignani S. Cavuoti A. Cimatti C. Colodro-Conde G. Congedo C. J. Conselice L. Conversi Y. Copin F. Courbin H. M. Courtois A. Da Silva H. Degaudenzi G. De Lucia M. Douspis F. Dubath C. A. J. Duncan X. Dupac S. Dusini M. Farina S. Farrens S. Ferriol P. Fosalba M. Frailis E. Franceschi S. Galeotta B. Gillis P. Gómez-Alvarez A. Grazian F. Grupp L. Guzzo M. Hailey S. V. H. Haugan W. Holmes F. Hormuth A. Hornstrup P. Hudelot S. Ilić K. Jahnke M. Jhabvala B. Joachimi E. Keihänen S. Kermiche A. Kiessling M. Kilbinger B. Kubik M. Kunz H. Kurki-Suonio P. B. Lilje V. Lindholm I. Lloro G. Mainetti D. Maino E. Maiorano O. Mansutti O. Marggraf K. Markovic M. Martinelli N. Martinet F. Marulli R. Massey E. Medinaceli S. Mei M. Melchior Y. Mellier M. Meneghetti E. Merlin G. Meylan M. Moresco L. Moscardini E. Munari C. Neissner S. -M. Niemi C. Padilla S. Paltani F. Pasian K. Pedersen W. J. Percival V. Pettorino S. Pires G. Polenta M. Poncet L. A. Popa L. Pozzetti F. Raison A. Renzi J. Rhodes G. Riccio E. Romelli M. Roncarelli R. Saglia J. -C. Salvignol A. G. Sánchez D. Sapone B. Sartoris M. Schirmer T. Schrabback A. Secroun G. Seidel S. Serrano C. Sirignano G. Sirri L. Stanco J. Steinwagner P. Tallada-Crespí A. N. Taylor I. Tereno R. Toledo-Moreo F. Torradeflot I. Tutusaus L. Valenziano T. Vassallo G. Verdoes Kleijn A. Veropalumbo Y. Wang J. Weller G. Zamorani E. Zucca A. Biviano E. Bozzo C. Burigana M. Calabrese D. Di Ferdinando J. A. Escartin Vigo G. Fabbian R. Farinelli F. Finelli J. Gracia-Carpio S. Matthew N. Mauri A. Pezzotta M. Pöntinen V. Scottez M. Tenti M. Viel M. Wiesmann Y. Akrami S. Anselmi M. Archidiacono F. Atrio-Barandela M. Ballardini D. Bertacca A. Blanchard L. Blot H. Böhringer S. Bruton R. Cabanac A. Calabro B. Camacho Quevedo G. Cañas-Herrera A. Cappi F. Caro C. S. Carvalho T. Castro K. C. Chambers S. Contarini A. R. Cooray G. Desprez A. Díaz-Sánchez J. J. Diaz S. Di Domizio H. Dole S. Escoffier M. Ezziati A. G. Ferrari P. G. Ferreira I. Ferrero A. Finoguenov A. Fontana F. Fornari L. Gabarra K. Ganga J. García-Bellido T. Gasparetto V. Gautard E. Gaztanaga F. Giacomini F. Gianotti G. Gozaliasl C. M. Gutierrez A. Hall H. Hildebrandt J. Hjorth A. Jimenez Muñoz S. Joudaki J. J. E. Kajava V. Kansal D. Karagiannis C. C. Kirkpatrick J. Le Graet L. Legrand J. Lesgourgues T. I. Liaudat S. J. Liu A. Loureiro G. Maggio M. Magliocchetti F. Mannucci R. Maoli J. Martín-Fleitas C. J. A. P. Martins L. Maurin R. B. Metcalf M. Miluzio P. Monaco A. Montoro A. Mora C. Moretti G. Morgante C. Murray S. Nadathur Nicholas A. Walton L. Pagano L. Patrizii V. Popa D. Potter P. Reimberg I. Risso P. -F. Rocci M. Sahlén E. Sarpa A. Schneider M. Sereno A. Silvestri A. Spurio Mancini J. Stadel K. Tanidis C. Tao N. Tessore G. Testera R. Teyssier S. Toft S. Tosi A. Troja M. Tucci J. Valiviita D. Vergani G. Verza P. Vielzeuf





Simulations and nonlinearities beyond ΛCDM. 4. Constraints on f(R) models from the photometric primary probes    



 



Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth PO1 3FX, UK
	Institute for Theoretical Particle Physics and Cosmology (TTK), RWTH Aachen University, 52056 Aachen, Germany
	Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK
	Laboratoire Univers et Théorie, Observatoire de Paris, Université PSL, Université Paris Cité, CNRS, 92190 Meudon, France
	Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências, Universidade de Lisboa, Campo Grande, 1749-016 Lisboa, Portugal
	Departamento de Física, Faculdade de Ciências, Universidade de Lisboa, Edifício C8, Campo Grande, PT1749-016 Lisboa, Portugal
	Ruhr University Bochum, Faculty of Physics and Astronomy, Astronomical Institute (AIRUB), German Centre for Cosmological Lensing (GCCL), 44780 Bochum, Germany
	Universität Bonn, Argelander-Institut für Astronomie, Auf dem Hügel 71, 53121 Bonn, Germany
	INAF-Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Piero Gobetti 93/3, 40129 Bologna, Italy
	Istituto Nazionale di Fisica Nucleare, Sezione di Bologna, Via Irnerio 46, 40126 Bologna, Italy
	Dipartimento di Fisica, Università degli Studi di Torino, Via P. Giuria 1, 10125 Torino, Italy
	INFN-Sezione di Torino, Via P. Giuria 1, 10125 Torino, Italy
	INAF-Osservatorio Astrofisico di Torino, Via Osservatorio 20, 10025 Pino Torinese (TO), Italy
	Higgs Centre for Theoretical Physics, School of Physics and Astronomy, The University of Edinburgh, Edinburgh EH9 3FD, UK
	Institut universitaire de France (IUF), 1 rue Descartes, 75231 PARIS CEDEX 05, France
	Institut für Theoretische Physik, University of Heidelberg, Philosophenweg 16, 69120 Heidelberg, Germany
	Institut de Recherche en Astrophysique et Planétologie (IRAP), Université de Toulouse, CNRS, UPS, CNES, 14 Av. Edouard Belin, 31400 Toulouse, France
	Université St Joseph; Faculty of Sciences, Beirut, Lebanon
	Institute of Theoretical Astrophysics, University of Oslo, P.O. Box 1029 Blindern, 0315 Oslo, Norway
	Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, University of Manchester, Oxford Road, Manchester M13 9PL, UK
	Department of Astrophysics, University of Zurich, Winterthurerstrasse 190, 8057 Zurich, Switzerland
	Dipartimento di Fisica e Astronomia, Università di Bologna, Via Gobetti 93/2, 40129 Bologna, Italy
	INFN-Sezione di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy
	Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s/n, 08193 Barcelona, Spain
	Institut de Ciencies de l'Espai (IEEC-CSIC), Campus UAB, Carrer de Can Magrans, s/n Cerdanyola del Vallés, 08193 Barcelona, Spain
	Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA, 91109, USA
	Institut de Physique Théorique, CEA, CNRS, Université Paris-Saclay 91191 Gif-sur-Yvette Cedex, France
	School of Mathematics and Physics, University of Surrey, Guildford, Surrey, GU2 7XH, UK
	INAF-Osservatorio Astronomico di Brera, Via Brera 28, 20122 Milano, Italy
	IFPU, Institute for Fundamental Physics of the Universe, via Beirut 2, 34151 Trieste, Italy
	INAF-Osservatorio Astronomico di Trieste, Via G. B. Tiepolo 11, 34143 Trieste, Italy
	INFN, Sezione di Trieste, Via Valerio 2, 34127 Trieste TS, Italy
	SISSA, International School for Advanced Studies, Via Bonomea 265, 34136 Trieste TS, Italy
	Institut d'Astrophysique de Paris, UMR 7095, CNRS, and Sorbonne Université, 98 bis boulevard Arago, 75014 Paris, France
	Max Planck Institute for Extraterrestrial Physics, Giessenbachstr. 1, 85748 Garching, Germany
	Dipartimento di Fisica, Università di Genova, Via Dodecaneso 33, 16146, Genova, Italy
	INFN-Sezione di Genova, Via Dodecaneso 33, 16146, Genova, Italy
	Department of Physics "E. Pancini", University Federico II, Via Cinthia 6, 80126, Napoli, Italy
	INAF-Osservatorio Astronomico di Capodimonte, Via Moiariello 16, 80131 Napoli, Italy
	INFN section of Naples, Via Cinthia 6, 80126, Napoli, Italy
	Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, PT4150-762 Porto, Portugal
	Faculdade de Ciências da Universidade do Porto, Rua do Campo de Alegre, 4150-007 Porto, Portugal
	Aix-Marseille Université, CNRS, CNES, LAM, Marseille, France
	INAF-IASF Milano, Via Alfonso Corti 12, 20133 Milano, Italy
	Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Avenida Complutense 40, 28040 Madrid, Spain
	Port d'Informació Científica, Campus UAB, C. Albareda s/n, 08193 Bellaterra (Barcelona), Spain
	INAF-Osservatorio Astronomico di Roma, Via Frascati 33, 00078 Monteporzio Catone, Italy
	Dipartimento di Fisica e Astronomia "Augusto Righi" - Alma Mater Studiorum Università di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy
	Instituto de Astrofísica de Canarias, Calle Vía Láctea s/n, 38204, San Cristóbal de La Laguna, Tenerife, Spain
	European Space Agency/ESRIN, Largo Galileo Galilei 1, 00044 Frascati, Roma, Italy
	ESAC/ESA, Camino Bajo del Castillo, s/n., Urb. Villafranca del Castillo, 28692 Villanueva de la Cañada, Madrid, Spain
	Université Claude Bernard Lyon 1, CNRS/IN2P3, IP2I Lyon, UMR 5822, Villeurbanne, F-69100, France
	Institute of Physics, Laboratory of Astrophysics, Ecole Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, 1290 Versoix, Switzerland
	Institut de Ciències del Cosmos (ICCUB), Universitat de Barcelona (IEEC-UB), Martí i Franquès 1, 08028 Barcelona, Spain
	Institució Catalana de Recerca i Estudis Avançats (ICREA), Passeig de Lluís Companys 23, 08010 Barcelona, Spain
	UCB Lyon 1, CNRS/IN2P3, IUF, IP2I Lyon, 4 rue Enrico Fermi, 69622 Villeurbanne, France
	Department of Astronomy, University of Geneva, ch. d'Ecogia 16, 1290 Versoix, Switzerland
	Université Paris-Saclay, CNRS, Institut d'astrophysique spatiale, 91405, Orsay, France
	INFN-Padova, Via Marzolo 8, 35131 Padova, Italy
	INAF-Istituto di Astrofisica e Planetologia Spaziali, via del Fosso del Cavaliere, 100, 00100 Roma, Italy
	Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191, Gif-sur-Yvette, France
	Institut d'Estudis Espacials de Catalunya (IEEC),  Edifici RDIT, Campus UPC, 08860 Castelldefels, Barcelona, Spain
	FRACTAL S.L.N.E., calle Tulipán 2, Portal 13 1A, 28231, Las Rozas de Madrid, Spain
	INAF-Osservatorio Astronomico di Padova, Via dell'Osservatorio 5, 35122 Padova, Italy
	Universitäts-Sternwarte München, Fakultät für Physik, Ludwig-Maximilians-Universität München, Scheinerstrasse 1, 81679 München, Germany
	Dipartimento di Fisica "Aldo Pontremoli", Università degli Studi di Milano, Via Celoria 16, 20133 Milano, Italy
	Mullard Space Science Laboratory, University College London, Holmbury St Mary, Dorking, Surrey RH5 6NT, UK
	Felix Hormuth Engineering, Goethestr. 17, 69181 Leimen, Germany
	Technical University of Denmark, Elektrovej 327, 2800 Kgs. Lyngby, Denmark
	Cosmic Dawn Center (DAWN), Denmark
	Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France
	Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany
	NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA
	Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK
	Department of Physics and Helsinki Institute of Physics, Gustaf Hällströmin katu 2, 00014 University of Helsinki, Finland
	Aix-Marseille Université, CNRS/IN2P3, CPPM, Marseille, France
	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
	Department of Physics, P.O. Box 64, 00014 University of Helsinki, Finland
	Helsinki Institute of Physics, Gustaf Hällströmin katu 2, University of Helsinki, Helsinki, Finland
	NOVA optical infrared instrumentation group at ASTRON, Oude Hoogeveensedijk 4, 7991PD, Dwingeloo, The Netherlands
	Centre de Calcul de l'IN2P3/CNRS, 21 avenue Pierre de Coubertin 69627 Villeurbanne Cedex, France
	INFN-Sezione di Milano, Via Celoria 16, 20133 Milano, Italy
	INFN-Sezione di Roma, Piazzale Aldo Moro, 2 - c/o Dipartimento di Fisica, Edificio G. Marconi, 00185 Roma, Italy
	Dipartimento di Fisica e Astronomia "Augusto Righi" - Alma Mater Studiorum Università di Bologna, via Piero Gobetti 93/2, 40129 Bologna, Italy
	Department of Physics, Institute for Computational Cosmology, Durham University, South Road, DH1 3LE, UK
	Université Paris Cité, CNRS, Astroparticule et Cosmologie, 75013 Paris, France
	University of Applied Sciences and Arts of Northwestern Switzerland, School of Engineering, 5210 Windisch, Switzerland
	Institut d'Astrophysique de Paris, 98bis Boulevard Arago, 75014, Paris, France
	Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, 08193 Bellaterra (Barcelona), Spain
	European Space Agency/ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands
	DARK, Niels Bohr Institute, University of Copenhagen, Jagtvej 155, 2200 Copenhagen, Denmark
	Waterloo Centre for Astrophysics, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada
	Department of Physics and Astronomy, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada
	Perimeter Institute for Theoretical Physics, Waterloo, Ontario N2L 2Y5, Canada
	Space Science Data Center, Italian Space Agency, via del Politecnico snc, 00133 Roma, Italy
	Centre National d'Etudes Spatiales – Centre spatial de Toulouse, 18 avenue Edouard Belin, 31401 Toulouse Cedex 9, France
	Institute of Space Science, Str. Atomistilor, nr. 409 Măgurele, Ilfov, 077125, Romania
	Dipartimento di Fisica e Astronomia "G. Galilei", Università di Padova, Via Marzolo 8, 35131 Padova, Italy
	Departamento de Física, FCFM, Universidad de Chile, Blanco Encalada 2008, Santiago, Chile
	Universität Innsbruck, Institut für Astro- und Teilchenphysik, Technikerstr. 25/8, 6020 Innsbruck, Austria
	Satlantis, University Science Park, Sede Bld 48940, Leioa-Bilbao, Spain
	Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências, Universidade de Lisboa, Tapada da Ajuda, 1349-018 Lisboa, Portugal
	Universidad Politécnica de Cartagena, Departamento de Electrónica y Tecnología de Computadoras,  Plaza del Hospital 1, 30202 Cartagena, Spain
	INFN-Bologna, Via Irnerio 46, 40126 Bologna, Italy
	Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700 AV Groningen, The Netherlands
	Dipartimento di Fisica, Università degli studi di Genova, and INFN-Sezione di Genova, via Dodecaneso 33, 16146, Genova, Italy
	Infrared Processing and Analysis Center, California Institute of Technology, Pasadena, CA 91125, USA
	INAF, Istituto di Radioastronomia, Via Piero Gobetti 101, 40129 Bologna, Italy
	Astronomical Observatory of the Autonomous Region of the Aosta Valley (OAVdA), Loc. Lignan 39, I-11020, Nus (Aosta Valley), Italy
	Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
	School of Physics and Astronomy, Cardiff University, The Parade, Cardiff, CF24 3AA, UK
	Junia, EPA department, 41 Bd Vauban, 59800 Lille, France
	ICSC - Centro Nazionale di Ricerca in High Performance Computing, Big Data e Quantum Computing, Via Magnanelli 2, Bologna, Italy
	Instituto de Física Teórica UAM-CSIC, Campus de Cantoblanco, 28049 Madrid, Spain
	CERCA/ISO, Department of Physics, Case Western Reserve University, 10900 Euclid Avenue, Cleveland, OH 44106, USA
	Departamento de Física Fundamental. Universidad de Salamanca. Plaza de la Merced s/n. 37008 Salamanca, Spain
	Dipartimento di Fisica e Scienze della Terra, Università degli Studi di Ferrara, Via Giuseppe Saragat 1, 44122 Ferrara, Italy
	Istituto Nazionale di Fisica Nucleare, Sezione di Ferrara, Via Giuseppe Saragat 1, 44122 Ferrara, Italy
	Center for Data-Driven Discovery, Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan
	Ludwig-Maximilians-University, Schellingstrasse 4, 80799 Munich, Germany
	Max-Planck-Institut für Physik, Boltzmannstr. 8, 85748 Garching, Germany
	Minnesota Institute for Astrophysics, University of Minnesota, 116 Church St SE, Minneapolis, MN 55455, USA
	Institute Lorentz, Leiden University, Niels Bohrweg 2, 2333 CA Leiden, The Netherlands
	Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Bd de l'Observatoire, CS 34229, 06304 Nice cedex 4, France
	Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu, HI 96822, USA
	Department of Physics & Astronomy, University of California Irvine, Irvine CA 92697, USA
	Department of Astronomy & Physics and Institute for Computational Astrophysics, Saint Mary's University, 923 Robie Street, Halifax, Nova Scotia, B3H 3C3, Canada
	Departamento Física Aplicada, Universidad Politécnica de Cartagena, Campus Muralla del Mar, 30202 Cartagena, Murcia, Spain
	Instituto de Astrofísica de Canarias (IAC); Departamento de Astrofísica, Universidad de La Laguna (ULL), 38200, La Laguna, Tenerife, Spain
	Department of Physics, Oxford University, Keble Road, Oxford OX1 3RH, UK
	CEA Saclay, DFR/IRFU, Service d'Astrophysique, Bat. 709, 91191 Gif-sur-Yvette, France
	Department of Computer Science, Aalto University, PO Box 15400, Espoo, FI-00 076, Finland
	Instituto de Astrofísica de Canarias, c/ Via Lactea s/n, La Laguna E-38200, Spain. Departamento de Astrofísica de la Universidad de La Laguna, Avda. Francisco Sanchez, La Laguna, E-38200, Spain
	Univ. Grenoble Alpes, CNRS, Grenoble INP, LPSC-IN2P3, 53, Avenue des Martyrs, 38000, Grenoble, France
	Department of Physics and Astronomy, Vesilinnantie 5, 20014 University of Turku, Finland
	Serco for European Space Agency (ESA), Camino bajo del Castillo, s/n, Urbanizacion Villafranca del Castillo, Villanueva de la Cañada, 28692 Madrid, Spain
	ARC Centre of Excellence for Dark Matter Particle Physics, Melbourne, Australia
	Centre for Astrophysics & Supercomputing, Swinburne University of Technology,  Hawthorn, Victoria 3122, Australia
	School of Physics and Astronomy, Queen Mary University of London, Mile End Road, London E1 4NS, UK
	Department of Physics and Astronomy, University of the Western Cape, Bellville, Cape Town, 7535, South Africa
	ICTP South American Institute for Fundamental Research, Instituto de Física Teórica, Universidade Estadual Paulista, São Paulo, Brazil
	IRFU, CEA, Université Paris-Saclay 91191 Gif-sur-Yvette Cedex, France
	Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, Stockholm, SE-106 91, Sweden
	Astrophysics Group, Blackett Laboratory, Imperial College London, London SW7 2AZ, UK
	INAF-Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125, Firenze, Italy
	Dipartimento di Fisica, Sapienza Università di Roma, Piazzale Aldo Moro 2, 00185 Roma, Italy
	Aurora Technology for European Space Agency (ESA), Camino bajo del Castillo, s/n, Urbanizacion Villafranca del Castillo, Villanueva de la Cañada, 28692 Madrid, Spain
	Centro de Astrofísica da Universidade do Porto, Rua das Estrelas, 4150-762 Porto, Portugal
	HE Space for European Space Agency (ESA), Camino bajo del Castillo, s/n, Urbanizacion Villafranca del Castillo, Villanueva de la Cañada, 28692 Madrid, Spain
	Dipartimento di Fisica - Sezione di Astronomia, Università di Trieste, Via Tiepolo 11, 34131 Trieste, Italy
	Theoretical astrophysics, Department of Physics and Astronomy, Uppsala University, Box 515, 751 20 Uppsala, Sweden
	Department of Physics, Royal Holloway, University of London, TW20 0EX, UK
	Department of Astrophysical Sciences, Peyton Hall, Princeton University, Princeton, NJ 08544, USA
	Cosmic Dawn Center (DAWN)
	Niels Bohr Institute, University of Copenhagen, Jagtvej 128, 2200 Copenhagen, Denmark
	Center for Cosmology and Particle Physics, Department of Physics, New York University, New York, NY 10003, USA
	Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, 10010, New York, NY, USA    




Euclid Collaboration: K. Koyama et al.

Constraints on f(R) models from the photometric primary probes





  
  
We study the constraint on f(R) gravity that can be obtained by photometric primary probes of the mission. Our focus is the dependence of the constraint on the theoretical modelling of the nonlinear matter power spectrum. In the Hu–Sawicki f(R) gravity model, we consider four different predictions for the ratio between the power spectrum in f(R) and that in : a fitting formula, the halo model reaction approach,  and two emulators based on dark matter only N-body simulations,  and . These predictions are added to the  implementation to predict the angular power spectra for weak lensing (WL), photometric galaxy clustering and their cross-correlation. By running Markov Chain Monte Carlo, we compare constraints on parameters and investigate the bias of the recovered f(R) parameter if the data are created by a different model. For the pessimistic setting of WL, one dimensional bias for the f(R) parameter, , is found to be 0.5 σ when  is used to create the synthetic data with =-5.301 and fitted by . The impact of baryonic physics on WL is studied by using a baryonification emulator . For the optimistic setting, the f(R) parameter and two main baryon parameters are well constrained despite the degeneracies among these parameters. However, the difference in the nonlinear dark matter prediction can be compensated by the adjustment of baryon parameters, and the one-dimensional marginalised constraint on  is biased. This bias can be avoided in the pessimistic setting at the expense of weaker constraints. For the pessimistic setting, using the  synthetic data for WL, we obtain the prior-independent upper limit of < -5.6. Finally, we implement a method to include theoretical errors to avoid the bias due to inaccuracies in the nonlinear matter power spectrum prediction. 



preparation
    Euclid Collaboration: K. Koyama0000-0001-6727-6915kazuya.koyama@port.ac.uk<ref>
	S. Pamuk0009-0004-0852-8624<ref>
	S. Casas0000-0002-4751-5138<ref>
	B. Bose0000-0003-1965-8614<ref>
	P. Carrilho0000-0003-1339-0194<ref>
	I. Sáez-Casares0000-0003-0013-5266<ref>
	L. Atayde0000-0001-6373-9193<ref>,<ref>
	M. Cataneo0000-0002-7992-0656<ref>,<ref>
	B. Fiorini0000-0002-0092-4321<ref>
	C. Giocoli0000-0002-9590-7961<ref>,<ref>
	A. M. C. Le Brun0000-0002-0936-4594<ref>
	F. Pace0000-0001-8039-0480<ref>,<ref>,<ref>
	A. Pourtsidou0000-0001-9110-5550<ref>,<ref>
	Y. Rasera0000-0003-3424-6941<ref>,<ref>
	Z. Sakr0000-0002-4823-3757<ref>,<ref>,<ref>
	H.-A. Winther0000-0002-6325-2710<ref>
	E. Altamura0000-0001-6973-1897<ref>
	J. Adamek0000-0002-0723-6740<ref>
	M. Baldi0000-0003-4145-1943<ref>,<ref>,<ref>
	M.-A. Breton<ref>,<ref>,<ref>
	G. Rácz0000-0003-3906-5699<ref>
	F. Vernizzi0000-0003-3426-2802<ref>
	A. Amara<ref>
	S. Andreon0000-0002-2041-8784<ref>
	N. Auricchio0000-0003-4444-8651<ref>
	C. Baccigalupi0000-0002-8211-1630<ref>,<ref>,<ref>,<ref>
	S. Bardelli0000-0002-8900-0298<ref>
	F. Bernardeau<ref>,<ref>
	C. Bodendorf<ref>
	D. Bonino0000-0002-3336-9977<ref>
	E. Branchini0000-0002-0808-6908<ref>,<ref>,<ref>
	M. Brescia0000-0001-9506-5680<ref>,<ref>,<ref>
	J. Brinchmann0000-0003-4359-8797<ref>,<ref>
	A. Caillat<ref>
	S. Camera0000-0003-3399-3574<ref>,<ref>,<ref>
	V. Capobianco0000-0002-3309-7692<ref>
	C. Carbone0000-0003-0125-3563<ref>
	J. Carretero0000-0002-3130-0204<ref>,<ref>
	M. Castellano0000-0001-9875-8263<ref>
	G. Castignani0000-0001-6831-0687<ref>
	S. Cavuoti0000-0002-3787-4196<ref>,<ref>
	A. Cimatti<ref>
	C. Colodro-Conde<ref>
	G. Congedo0000-0003-2508-0046<ref>
	C. J. Conselice0000-0003-1949-7638<ref>
	L. Conversi0000-0002-6710-8476<ref>,<ref>
	Y. Copin0000-0002-5317-7518<ref>
	F. Courbin0000-0003-0758-6510<ref>,<ref>,<ref>
	H. M. Courtois0000-0003-0509-1776<ref>
	A. Da Silva0000-0002-6385-1609<ref>,<ref>
	H. Degaudenzi0000-0002-5887-6799<ref>
	G. De Lucia0000-0002-6220-9104<ref>
	M. Douspis0000-0003-4203-3954<ref>
	F. Dubath0000-0002-6533-2810<ref>
	C. A. J. Duncan<ref>
	X. Dupac<ref>
	S. Dusini0000-0002-1128-0664<ref>
	M. Farina0000-0002-3089-7846<ref>
	S. Farrens0000-0002-9594-9387<ref>
	S. Ferriol<ref>
	P. Fosalba0000-0002-1510-5214<ref>,<ref>
	M. Frailis0000-0002-7400-2135<ref>
	E. Franceschi0000-0002-0585-6591<ref>
	S. Galeotta0000-0002-3748-5115<ref>
	B. Gillis0000-0002-4478-1270<ref>
	P. Gómez-Alvarez0000-0002-8594-5358<ref>,<ref>
	A. Grazian0000-0002-5688-0663<ref>
	F. Grupp<ref>,<ref>
	L. Guzzo0000-0001-8264-5192<ref>,<ref>
	M. Hailey<ref>
	S. V. H. Haugan0000-0001-9648-7260<ref>
	W. Holmes<ref>
	F. Hormuth<ref>
	A. Hornstrup0000-0002-3363-0936<ref>,<ref>
	P. Hudelot<ref>
	S. Ilić0000-0003-4285-9086<ref>,<ref>
	K. Jahnke0000-0003-3804-2137<ref>
	M. Jhabvala<ref>
	B. Joachimi0000-0001-7494-1303<ref>
	E. Keihänen0000-0003-1804-7715<ref>
	S. Kermiche0000-0002-0302-5735<ref>
	A. Kiessling0000-0002-2590-1273<ref>
	M. Kilbinger0000-0001-9513-7138<ref>
	B. Kubik0009-0006-5823-4880<ref>
	M. Kunz0000-0002-3052-7394<ref>
	H. Kurki-Suonio0000-0002-4618-3063<ref>,<ref>
	P. B. Lilje0000-0003-4324-7794<ref>
	V. Lindholm0000-0003-2317-5471<ref>,<ref>
	I. Lloro<ref>
	G. Mainetti0000-0003-2384-2377<ref>
	D. Maino<ref>,<ref>,<ref>
	E. Maiorano0000-0003-2593-4355<ref>
	O. Mansutti0000-0001-5758-4658<ref>
	O. Marggraf0000-0001-7242-3852<ref>
	K. Markovic0000-0001-6764-073X<ref>
	M. Martinelli0000-0002-6943-7732<ref>,<ref>
	N. Martinet0000-0003-2786-7790<ref>
	F. Marulli0000-0002-8850-0303<ref>,<ref>,<ref>
	R. Massey0000-0002-6085-3780<ref>
	E. Medinaceli0000-0002-4040-7783<ref>
	S. Mei0000-0002-2849-559X<ref>
	M. Melchior<ref>
	Y. Mellier<ref>,<ref>
	M. Meneghetti0000-0003-1225-7084<ref>,<ref>
	E. Merlin0000-0001-6870-8900<ref>
	G. Meylan<ref>
	M. Moresco0000-0002-7616-7136<ref>,<ref>
	L. Moscardini0000-0002-3473-6716<ref>,<ref>,<ref>
	E. Munari0000-0002-1751-5946<ref>,<ref>
	C. Neissner0000-0001-8524-4968<ref>,<ref>
	S.-M. Niemi<ref>
	C. Padilla0000-0001-7951-0166<ref>
	S. Paltani0000-0002-8108-9179<ref>
	F. Pasian0000-0002-4869-3227<ref>
	K. Pedersen<ref>
	W. J. Percival0000-0002-0644-5727<ref>,<ref>,<ref>
	V. Pettorino<ref>
	S. Pires0000-0002-0249-2104<ref>
	G. Polenta0000-0003-4067-9196<ref>
	M. Poncet<ref>
	L. A. Popa<ref>
	L. Pozzetti0000-0001-7085-0412<ref>
	F. Raison0000-0002-7819-6918<ref>
	A. Renzi0000-0001-9856-1970<ref>,<ref>
	J. Rhodes0000-0002-4485-8549<ref>
	G. Riccio<ref>
	E. Romelli0000-0003-3069-9222<ref>
	M. Roncarelli0000-0001-9587-7822<ref>
	R. Saglia0000-0003-0378-7032<ref>,<ref>
	J.-C. Salvignol<ref>
	A. G. Sánchez0000-0003-1198-831X<ref>
	D. Sapone0000-0001-7089-4503<ref>
	B. Sartoris0000-0003-1337-5269<ref>,<ref>
	M. Schirmer0000-0003-2568-9994<ref>
	T. Schrabback0000-0002-6987-7834<ref>
	A. Secroun0000-0003-0505-3710<ref>
	G. Seidel0000-0003-2907-353X<ref>
	S. Serrano0000-0002-0211-2861<ref>,<ref>,<ref>
	C. Sirignano0000-0002-0995-7146<ref>,<ref>
	G. Sirri0000-0003-2626-2853<ref>
	L. Stanco0000-0002-9706-5104<ref>
	J. Steinwagner0000-0001-7443-1047<ref>
	P. Tallada-Crespí0000-0002-1336-8328<ref>,<ref>
	A. N. Taylor<ref>
	I. Tereno<ref>,<ref>
	R. Toledo-Moreo0000-0002-2997-4859<ref>
	F. Torradeflot0000-0003-1160-1517<ref>,<ref>
	I. Tutusaus0000-0002-3199-0399<ref>
	L. Valenziano0000-0002-1170-0104<ref>,<ref>
	T. Vassallo0000-0001-6512-6358<ref>,<ref>
	G. Verdoes Kleijn0000-0001-5803-2580<ref>
	A. Veropalumbo0000-0003-2387-1194<ref>,<ref>,<ref>
	Y. Wang0000-0002-4749-2984<ref>
	J. Weller0000-0002-8282-2010<ref>,<ref>
	G. Zamorani0000-0002-2318-301X<ref>
	E. Zucca0000-0002-5845-8132<ref>
	A. Biviano0000-0002-0857-0732<ref>,<ref>
	E. Bozzo0000-0002-8201-1525<ref>
	C. Burigana0000-0002-3005-5796<ref>,<ref>
	M. Calabrese0000-0002-2637-2422<ref>,<ref>
	D. Di Ferdinando<ref>
	J. A. Escartin Vigo<ref>
	G. Fabbian0000-0002-3255-4695<ref>,<ref>
	R. Farinelli<ref>
	F. Finelli0000-0002-6694-3269<ref>,<ref>
	J. Gracia-Carpio<ref>
	S. Matthew0000-0001-8448-1697<ref>
	N. Mauri0000-0001-8196-1548<ref>,<ref>
	A. Pezzotta0000-0003-0726-2268<ref>
	M. Pöntinen0000-0001-5442-2530<ref>
	V. Scottez<ref>,<ref>
	M. Tenti0000-0002-4254-5901<ref>
	M. Viel0000-0002-2642-5707<ref>,<ref>,<ref>,<ref>,<ref>
	M. Wiesmann0009-0000-8199-5860<ref>
	Y. Akrami0000-0002-2407-7956<ref>,<ref>
	S. Anselmi0000-0002-3579-9583<ref>,<ref>,<ref>
	M. Archidiacono0000-0003-4952-9012<ref>,<ref>
	F. Atrio-Barandela0000-0002-2130-2513<ref>
	M. Ballardini0000-0003-4481-3559<ref>,<ref>,<ref>
	D. Bertacca0000-0002-2490-7139<ref>,<ref>,<ref>
	A. Blanchard0000-0001-8555-9003<ref>
	L. Blot0000-0002-9622-7167<ref>,<ref>
	H. Böhringer0000-0001-8241-4204<ref>,<ref>,<ref>
	S. Bruton0000-0002-6503-5218<ref>
	R. Cabanac0000-0001-6679-2600<ref>
	A. Calabro0000-0003-2536-1614<ref>
	B. Camacho Quevedo0000-0002-8789-4232<ref>,<ref>
	G. Cañas-Herrera0000-0003-2796-2149<ref>,<ref>
	A. Cappi<ref>,<ref>
	F. Caro<ref>
	C. S. Carvalho<ref>
	T. Castro0000-0002-6292-3228<ref>,<ref>,<ref>,<ref>
	K. C. Chambers0000-0001-6965-7789<ref>
	S. Contarini0000-0002-9843-723X<ref>
	A. R. Cooray0000-0002-3892-0190<ref>
	G. Desprez0000-0001-8325-1742<ref>
	A. Díaz-Sánchez0000-0003-0748-4768<ref>
	J. J. Diaz<ref>
	S. Di Domizio0000-0003-2863-5895<ref>,<ref>
	H. Dole0000-0002-9767-3839<ref>
	S. Escoffier0000-0002-2847-7498<ref>
	M. Ezziati0009-0003-6065-1585<ref>
	A. G. Ferrari0009-0005-5266-4110<ref>,<ref>
	P. G. Ferreira0000-0002-3021-2851<ref>
	I. Ferrero0000-0002-1295-1132<ref>
	A. Finoguenov0000-0002-4606-5403<ref>
	A. Fontana0000-0003-3820-2823<ref>
	F. Fornari0000-0003-2979-6738<ref>
	L. Gabarra0000-0002-8486-8856<ref>
	K. Ganga0000-0001-8159-8208<ref>
	J. García-Bellido0000-0002-9370-8360<ref>
	T. Gasparetto0000-0002-7913-4866<ref>
	V. Gautard<ref>
	E. Gaztanaga0000-0001-9632-0815<ref>,<ref>,<ref>
	F. Giacomini0000-0002-3129-2814<ref>
	F. Gianotti0000-0003-4666-119X<ref>
	G. Gozaliasl0000-0002-0236-919X<ref>,<ref>
	C. M. Gutierrez0000-0001-7854-783X<ref>
	A. Hall0000-0002-3139-8651<ref>
	H. Hildebrandt0000-0002-9814-3338<ref>
	J. Hjorth0000-0002-4571-2306<ref>
	A. Jimenez Muñoz0009-0004-5252-185X<ref>
	S. Joudaki0000-0001-8820-673X<ref>
	J. J. E. Kajava0000-0002-3010-8333<ref>,<ref>
	V. Kansal0000-0002-4008-6078<ref>,<ref>
	D. Karagiannis0000-0002-4927-0816<ref>,<ref>
	C. C. Kirkpatrick<ref>
	J. Le Graet0000-0001-6523-7971<ref>
	L. Legrand0000-0003-0610-5252<ref>
	J. Lesgourgues0000-0001-7627-353X<ref>
	T. I. Liaudat0000-0002-9104-314X<ref>
	S. J. Liu0000-0001-7680-2139<ref>
	A. Loureiro0000-0002-4371-0876<ref>,<ref>
	G. Maggio0000-0003-4020-4836<ref>
	M. Magliocchetti0000-0001-9158-4838<ref>
	F. Mannucci0000-0002-4803-2381<ref>
	R. Maoli0000-0002-6065-3025<ref>,<ref>
	J. Martín-Fleitas0000-0002-8594-569X<ref>
	C. J. A. P. Martins0000-0002-4886-9261<ref>,<ref>
	L. Maurin0000-0002-8406-0857<ref>
	R. B. Metcalf0000-0003-3167-2574<ref>,<ref>
	M. Miluzio<ref>,<ref>
	P. Monaco0000-0003-2083-7564<ref>,<ref>,<ref>,<ref>
	A. Montoro0000-0003-4730-8590<ref>,<ref>
	A. Mora0000-0002-1922-8529<ref>
	C. Moretti0000-0003-3314-8936<ref>,<ref>,<ref>,<ref>,<ref>
	G. Morgante<ref>
	C. Murray<ref>
	S. Nadathur0000-0001-9070-3102<ref>
	Nicholas A. Walton0000-0003-3983-8778<ref>
	L. Pagano0000-0003-1820-5998<ref>,<ref>
	L. Patrizii<ref>
	V. Popa0000-0002-9118-8330<ref>
	D. Potter0000-0002-0757-5195<ref>
	P. Reimberg0000-0003-3410-0280<ref>
	I. Risso0000-0003-2525-7761<ref>
	P.-F. Rocci<ref>
	M. Sahlén0000-0003-0973-4804<ref>
	E. Sarpa0000-0002-1256-655X<ref>,<ref>,<ref>
	A. Schneider0000-0001-7055-8104<ref>
	M. Sereno0000-0003-0302-0325<ref>,<ref>
	A. Silvestri0000-0001-6904-5061<ref>
	A. Spurio Mancini0000-0001-5698-0990<ref>,<ref>
	J. Stadel0000-0001-7565-8622<ref>
	K. Tanidis<ref>
	C. Tao0000-0001-7961-8177<ref>
	N. Tessore0000-0002-9696-7931<ref>
	G. Testera<ref>
	R. Teyssier0000-0001-7689-0933<ref>
	S. Toft0000-0003-3631-7176<ref>,<ref>
	S. Tosi0000-0002-7275-9193<ref>,<ref>
	A. Troja0000-0003-0239-4595<ref>,<ref>
	M. Tucci<ref>
	J. Valiviita0000-0001-6225-3693<ref>,<ref>
	D. Vergani0000-0003-0898-2216<ref>
	G. Verza0000-0002-1886-8348<ref>,<ref>
	P. Vielzeuf0000-0003-2035-9339<ref>
    September 9, 2024
=========================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================





§ INTRODUCTION

In 1998, astronomers made a surprising discovery that the expansion of the Universe is accelerating, not slowing down <cit.>.
This late-time acceleration of the Universe has become the most challenging problem in theoretical physics. The main aim of ongoing and future cosmological surveys is to address the key questions about the origin of the late-time acceleration. Is the acceleration driven by a cosmological constant or dark energy that evolves with the expansion of the Universe? Alternatively, there could be no dark energy if General Relativity (GR) itself is in error on cosmological scales. There has been significant progress in developing modified theories of gravity and these have been developed into tests of GR itself via cosmological observations  <cit.>. Testing gravity is one of the main objectives of stage-IV dark energy surveys <cit.>.  

Of particular importance in these surveys is the mission <cit.>. The satellite undertakes a spectroscopic survey of galaxies and an imaging survey (targeting weak lensing which can also be used to reconstruct galaxy clustering using photometric redshifts). The combination of these two probes is essential for cosmological tests of gravity <cit.>.

Modified gravity models typically include an additional scalar degree of freedom that gives rise to a fifth-force. To satisfy the stringent constraints on deviations from GR in the solar system, many modified gravity models utilise screening mechanisms to hide modifications of gravity on small scales <cit.>. This is accomplished by nonlinearity in the equation that governs the dynamics of the scalar degree of freedom. This significantly complicates the nonlinear modelling of matter clustering in these models as the nonlinear equation for the scalar mode coupled to nonlinear density fields needs to be solved. A systematic comparison of N-body simulations in f(R) gravity and normal-branch Dvali–Gabadadze–Porratti (nDGP) models was done in <cit.>. Since then, new simulations have been developed, including those using approximate methods to accelerate simulations <cit.>. A fitting formula <cit.> and emulators for the nonlinear matter power spectrum have been developed <cit.>. At the same time, a semi-analytic method based on the halo model to predict the nonlinear matter power spectrum for general dark energy and modified gravity models has been developed <cit.>. These nonlinear predictions were used to study their modifications of the Weak Lensing (WL) observables <cit.>, cross-correlation of galaxies and cosmic microwave background <cit.> for example in f(R) gravity. 

In Euclid:2023tqw, Fisher Matrix forecasts were performed to predict 's ability to constrain f(R) gravity models. In the Hu–Sawicki f(R) model, it was shown that in the optimistic setting defined in , and for a fiducial value of = 5 × 10^-6, alone will be able to constrain the additional parameter  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. The forecast for photometric probes used the fitting formula for the nonlinear matter power spectrum obtained in Winther:2019mus. Further Fisher Matrix forecasts have been done for other modified gravity models with linearly scale-independent growth <cit.>. 

For real data analysis, it is imperative to check the effect of the accuracy of the theoretical modelling. This article is part of a series that collectively explores simulations and nonlinearities beyond the Λ Cold Dark Matter () model:

  * Numerical methods and validation <cit.>
    
  * Results from non-standard simulations <cit.>
    
  * Cosmological constraints on non-standard cosmologies from simulated probes (D'Amico et al. in prep.)
    
  * Constraints on f(R) models from the photometric primary probes (this work)
In , the comparison of N-body simulations performed in Winther:2015wla in the Hu–Sawicki f(R) and nDGP models was extended to add more simulations. The comparison was done for the matter power spectrum and the halo mass function. The measurements of these quantities were done using the dedicated pipeline developed by . In , additional simulations have been analysed in addition to those used in . In this paper, we utilise these developments and compare several predictions for the nonlinear matter power spectrum in f(R) gravity. We will perform Markov Chain Monte Carlo (MCMC) simulations using synthetic data for photometric probes, and assess the impact of using different nonlinear models for the matter power spectrum at the level of parameter constraints. In addition, we will add baryonic effects using a baryonification method, which was not included in the Fisher Matrix forecast. Although this paper focuses on f(R) gravity, for which multiple public codes are available to predict the nonlinear matter power spectrum, the methodology and code developed in this paper are readily applicable to other modified gravity models. The validation of the nonlinear models for spectroscopic probes has been done in Euclid:2023bgs and D'Amico et al. (paper 3) will perform a similar analysis to this paper's for the spectroscopic probes. 

This paper is organised as follows. In <ref>, we summarise theoretical predictions for observables based on  and  Euclid:2023tqw. In <ref>, we introduce the Hu–Sawicki f(R) gravity model and summarise the four nonlinear models for the matter power spectrum. We then discuss the implementation of these models in the  code introduced in Euclid:2023pxu. In <ref>, we compare predictions for the nonlinear matter power spectrum with N-body simulations and study their impact on the angular power spectra for photometric probes. Forecasts for errors from the combination of photometric probes considered in  will be presented based on the synthetic data created by four different nonlinear models. We also compare the result with the Fisher Matrix forecast. We then study the bias in the recovered f(R) gravity parameter when the synthetic data are created by a different nonlinear model. <Ref> is devoted to the study of baryonic effects using the  baryonification method. We show how the bias is affected by baryonic effects and obtain the upper bound on |f_R0| using the  synthetic data. In <ref>, we implement theoretical errors to take into account the uncertainties of theoretical predictions. We conclude in <ref>. 



§ THEORETICAL PREDICTIONS FOR EUCLID OBSERVABLES


In this section, we discuss how moving away from the standard GR assumption impacts the predictions for the angular power spectra C(ℓ) that will be compared with the photometric data survey. The observables that need to be computed and compared with the data are the angular power spectra for weak lensing (WL), photometric galaxy clustering () and their cross-correlation ().
In , these were calculated using the Limber and flat-sky approximations in a flat  Universe, as 

    C^XY_ij(ℓ) = c ∫_z_ min^z_ max z W_i^X(z) W_j^Y(z)/H(z)  r^2(z)P_δδ(k_ℓ,z) ,

where W_i^X(z) is the window function in each tomographic redshift bin i with X={ L,G} corresponding to WL and , respectively, k_ℓ=(ℓ+1/2)/r(z), r(z) is the comoving distance to redshift z,  P_δδ(k_ℓ,z) is the nonlinear power spectrum of matter density fluctuations, δ, at wave number k_ℓ and redshift z, in the redshift range of the integral from z_ min=0.001 to z_ max=4 and H(z) is  the Hubble function. 

However, when abandoning the assumption of GR, one has to account for changes in the evolution of both the homogeneous background and cosmological perturbations. For WL, it is the Weyl potential ϕ_W that determines the angular power spectrum. The power spectrum of the Weyl potential is related to P_δδ as <cit.> 

    P_ϕ_W(k,z) = [-3  (H_0/c)^2 (1+z) Σ(k,z)]^2P_δδ(k,z) ,

where Σ(k,z) is a phenomenological parameterisation to account for deviations from the standard  lensing prediction.  is the density parameter of matter and H_0 is the Hubble parameter – both at the present time. Here we assumed a standard background evolution of the matter component, i.e. ρ_ m(z)=ρ_ m,0(1+z)^3. 

We can, therefore, use the recipe of <ref> with H, r and P_δδ provided by a Boltzmann solver, but with the new window functions <cit.>

    W_i^ G(k,z) =     1/c b_i(k,z) n_i(z) H(z) , 
      
     W_i^ L(k,z) =     3/2  (H_0/c)^2(1+z) r(z) Σ(k,z)  
       ×∫_z^z_ max z'  n_i(z) r(z')-r(z)/r(z') + W^ IA_i(k,z)  ,

where n_i(z) is the normalised galaxy number-density distribution in a tomographic redshift bin i such that ∫_z_ min^z_ max n_i(z) dz = 1, and W^ IA_i(k,z) encodes the contribution of intrinsic alignments (IA) to the WL power spectrum. We follow  in assuming an effective scale-independent galaxy bias, 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 = √(1+z̅_i), where z̅_i is the mean redshift of each redshift bin.

The IA contribution is computed following the nonlinear alignment model with a redshift dependent amplitude ,
in which
    W^ IA_i(k,z)=-𝒜_ IA 𝒞_ IA  ℱ_ IA(z)/δ(k,z) / δ(k,0) n_i(z) H(z)/c ,

where 

    ℱ_ IA(z)=(1+z)^η_ IA.

The parameters 𝒜_ IA and η_ IA are the nuisance parameters of the model, and 𝒞_ IA is a constant accounting for dimensional units. The galaxy distribution is binned into 10 equipopulated redshift bins with an overall distribution following

    n(z)∝(z/z_0)^2 exp[-(z/z_0)^3/2] ,

with z_0=0.9/√(2) and the normalisation set by the requirement that the surface density of galaxies is n̅_ g=30 arcmin^-2 .

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 <ref> through the matter perturbation δ(k,z), while the modifications on the clustering of matter in the  case are accounted for in the new P_δδ(k_ℓ,z).

Finally, we present the likelihood that we will use. We follow the approach presented in Euclid:2023pxu. We first construct an (N^G+N^L)× (N^G+N^L) angular power spectrum matrix for each multipole, where the different N correspond to the number of redshift bins for each probe (WL and ):


    C_ℓ = [ [ C^LL_ij(ℓ) C^GL_ij(ℓ); C^LG_ij(ℓ) C^GG_ij(ℓ) ]],

where lower-case Latin indexes i,… run over all tomographic bins. Similarly, the noise contributions can also be condensed into one noise matrix N:

    N_ℓ = [ [ N^LL_ij(ℓ) N^GL_ij(ℓ); N^LG_ij(ℓ) N^GG_ij(ℓ) ]],

where the noise terms N_ij^AB(ℓ) are given by

    N_ij^ LL(ℓ)    = δ_ij^ K/n̅_iσ_ϵ^2 ,
    
        N_ij^ GG(ℓ)    = δ_ij^ K/n̅_i ,
    
        N_ij^ GL(ℓ)    = 0 ,

where σ_ϵ^2=0.3^2 is the variance of observed ellipticities. We can further define Ĉ_ℓC_ℓ + N_ℓ. This is the covariance of the spherical multipole moments a_lm.

In the Gaussian approximation, it can be shown that the description in , using the covariance of the angular power spectra, is equivalent to the description found in Euclid:2023pxu, using the covariance of the a_ℓ m. The likelihood then can be expressed in terms of the observed covariance 

    Ĉ^obs_ij(ℓ) = 1/2 ℓ+1 ∑_m=-ℓ^ℓ(a_ℓ m)_i (a_ℓ m)^*_j ,
 and the theoretically predicted one Ĉ^th_ij(ℓ) as 

    χ^2 = f^sky ∑_ℓ = ℓ_min^ℓ_max (2ℓ+1)  [d^mix/d^th + ln(d^th/d^obs) - N] ,

where the determinants d are defined as 

    d^th(ℓ)    = [Ĉ^th_ij(ℓ)],
        
    
        d^obs(ℓ)    = [Ĉ^obs_ij(ℓ)], 
         
    
        d^mix(ℓ)    = ∑_k^N[ {[ Ĉ^obs_ij(ℓ)  for  k = j;  Ĉ^th_ij(ℓ)  for  k ≠ j ]. ]  .

Here N is the number of bins and thus either (N^G+N^L) for multipoles where we include the cross correlation, or the respective N for multipoles where we treat the two probes separately. The additional factor f^ sky comes from an approximation to account for having less available independent ℓ modes due to partial sky coverage.
In this paper, we set the observed covariance Ĉ^ obs to the theoretically predicted one computed at the fiducial cosmology.
Following  and Euclid:2023pxu, we do not include super-sample covariance (SSC) in this work. The SSC impact was shown to be non-negligible and will need to be included in future analyses as shown in Euclid:2023ove. 

We consider two different scenarios: an optimistic and a pessimistic case. In the optimistic case, we consider ℓ_ max=5000 for WL, and ℓ_ max=3000 for  and . Instead, in the pessimistic scenario, we consider ℓ_ max=1500 for WL, and ℓ_ max=750 for  and . 

In the smallest photometric redshift bin, the galaxy number density distribution n(z) peaks at around z=0.25. Under the Limber approximation for our fiducial cosmology, the corresponding maximum values of k evaluated in the power spectrum corresponding to the pessimistic and optimistic scenario for  are k_ max=[0.7, 2.9] h Mpc^-1, respectively, while for the WL, maximum wavenumbers probed are k_ max=[1.4, 4.8] h Mpc^-1 at the peak redshift z=0.25 <cit.>. Here h denotes the dimensionless Hubble parameter h  H_0 / (100 ). For smaller values of z, the values of k at a given ℓ increase, but the window functions in <ref> suppress the power spectrum and we set it to zero after a fixed k_ max=30 h Mpc^-1.
We list the specific choices of scales and settings used for each observable in <Ref>.
Although the currently foreseen specification of the Euclid Wide Survey differs from
that assumed in , for example in terms of the survey area, we use their results to allow for
comparison with earlier forecasts.
 
 



§  IMPLEMENTATION


In this section, we describe the implementation of the likelihood discussed in <ref> in the  code developed in Euclid:2023pxu, in the Hu–Sawicki f(R) gravity model.



 §.§ Hu–Sawicki fR gravity


Modified gravity f(R) models <cit.> are models constructed by promoting the Ricci scalar R in the Einstein–Hilbert action to a generic function of R→ R + f(R), i.e.

    S = c^4/16π G∫^4 x √(-g) [R+f(R)] + S_ m[g_μν, Ψ_ m]  ,

where g_μν is the metric tensor, g its determinant and S_ m represents the matter sector with its matter fields Ψ_ m. For further discussions we refer to Sotiriou_2010,Clifton_2012 and <cit.>.

How this modification changes gravity is more easily seen by formulating the theory in the so-called Einstein frame by performing a conformal transformation g_μν = g̃_μνA^2(ϕ) where A(ϕ) = √(1 + f_R) = e^ϕ/√(6) to obtain

    S =     c^4/16π G∫^4 x √(-g̃)[R̃ + 1/2(∂ϕ)^2 - V(ϕ)]
         + S_ m[A^2(ϕ) g̃_μν, Ψ_ m] ,

with the potential V = [f_R R - f(R)]/2(1+f_R)^2 and f_R f(R)/ R. This demonstrates that the theory reduces to standard GR together with an extra scalar field that is coupled to the matter sector giving rise to a fifth-force. This fifth-force has a finite range λ_C and within this range, it mediates a force that has a strength that, in the linear regime, corresponds to 1/3 of the usual gravitational force:

    F_ fifth = 1/3 G Mm/r^2  (1 + r/λ_C) e^-r/λ_C .

This effectively changes G→43G on small scales (r≪λ_C) while keeping the usual G on large scales. Such a large modification would be ruled out by observations if it was not for the fact that the theory possesses a screening mechanism <cit.> which hides the modifications in high-density regions. This screening mechanism effectively suppresses the fifth-force by a factor ∝ f_R / Φ_ N where Φ_ N is the Newtonian gravitational potential. 

Not all f(R) models one can write down possess such a screening mechanism as this places some restrictions on their functional form <cit.>. One concrete example of a model that has all the right ingredients is the model proposed by Hu:2007nk which in the large-curvature limit is defined by

    f(R) = - 6 H_0^2/c^2 + |f_R0| R̅_0^2/R ,

where 

    R̅_0 = 3 H_0^2/c^2(1+ 4 /) ,

is the Ricci scalar in the cosmological background and  is the density parameter of dark energy at the present time. The first term in <ref> corresponds to a cosmological constant and the only free parameter is |f_R0|. In the limit |f_R0| → 0, we recover GR and the  model. This parameter controls the range of the fifth-force and, in the cosmological background, we have λ_C≃ 32√(|f_R0|/10^-4)  Mpc at the present time. Solar System constraints require |f_R0| ≲ 10^-6, cosmological constraints currently lie around 10^-6–10^-4 depending on the probe in question <cit.> while astrophysical constraints at the galaxy scale can be as tight as ≲ 10^-8 <cit.>, but see <cit.> for a recent note of caution on such galactic-scale constraints.

The energy density of the scalar field contributes in general to the expansion of the Universe, however for viable models, like the one we consider here, this contribution is tiny (of the order |f_R0| Ω_ DE, 0) apart from the constant part of the potential which is indistinguishable from a cosmological constant. The background evolution of such models is therefore very close to . Since f(R) models have a conformal coupling, light deflection is weakly affected as follows

    Σ(z)=1/1+f_R(z) .

Since the maximum value of |f_R(z)| is given by |f_R0|, for the values of |f_R0| we consider in this paper, we can ignore this effect. Thus gravitational lensing is also not modified in the sense that the lensing potential is sourced by matter in the same way as in GR (though the underlying density field will of course be different). The main cosmological signatures of the model therefore come from having a fifth-force, acting only on small scales r ≲λ_C, in the process of structure formation. The main effect of the screening mechanism is that the prediction for the amount of clustering will generally be much smaller than what naive linear perturbation theory predicts.



 §.§ Nonlinear modelling

We implement Ξ(k,z) defined as 

    Ξ(k,z) P_f(R)(k,z)/P_(k,z) ,

to obtain the nonlinear f(R) matter power spectrum. For the  power spectrum P_(k,z), we use the  `Takahashi' prescription <cit.> following Euclid:2023tqw.  
It includes the minimum mass for massive neutrinos in P_(k,z), but ignores the effect of massive neutrinos on Ξ(k,z). This approximation was shown to be well justified for the minimum mass of neutrinos in <cit.> using data from <cit.>. We describe below four models for  Ξ(k,z) used in this paper. We note that we can use any  power spectrum prediction in our approach such as  <cit.> and  <cit.>. 
The exercise we perform in this paper is a comparison of the different prescriptions for Ξ(k, z). Given this,  the  nonlinear spectrum prescription does not matter and it is common to all nonlinear models.



  §.§.§ Fitting formula

A fitting formula for Ξ(k,z) was developed in <cit.> and describes the enhancement in the power spectrum compared to a  nonlinear power spectrum as a function of the parameter |f_R0|. This fitting function has been calibrated using the  <cit.> N-body simulations run by  <cit.> and the  N-body simulations <cit.> run by  <cit.>. The main approximation is that the cosmological parameter dependence of Ξ(k,z) is ignored. In <cit.>, this assumption was checked using simulations with different σ_8, , as well as the mass of massive neutrinos. <cit.> also corrected the fitting formula to account for additional dependence on these parameters. In this paper, we do not include these corrections as the previous forecast paper <cit.> used the fitting formula without these corrections. 

The fitting function has in total 54 parameters for the full scale, redshift and |f_R0| dependence. 
This fitting formula is not defined outside the range 10^-7 <  < 10^-4. The code is publicly available (https://github.com/HAWinther/FofrFittingFunctiongithub). 





  §.§.§ Halo model reaction

We further consider the halo model reaction approach of Cataneo:2018cic which combines the halo model and perturbation theory frameworks to model corrections coming from non-standard physics. The nonlinear power spectrum is given by


    P_ NL(k,z) = ℛ(k,z)   P^ pseudo_ NL(k,z) ,


where P^ pseudo_ NL(k,z) is called the pseudo-power spectrum and is defined as a nonlinear  spectrum with initial conditions tuned such that the linear clustering at the target redshift z matches the beyond- case. This choice ensures the halo mass functions in the beyond- and pseudo-universes are similar, giving a smoother transition of the power spectrum over inter- and intra-halo scales. We model the pseudo-cosmology nonlinear power spectrum using  <cit.> by supplying the code with the linear f(R) power spectrum. 

The halo model reaction, ℛ(k,z), models all the corrections to the pseudo spectrum coming from nonlinear beyond- physics. We refer the reader to Cataneo:2018cic,Bose:2021mkz and <cit.> for the exact expressions for this term. The halo model reaction can be computed efficiently using the publicly available  <cit.> code.

Despite  being highly efficient, having been used in previous real cosmic shear analyses <cit.>, it is still too computationally expensive for the number of tests we wish to perform. To accelerate our inference pipeline, we create a neural network emulator using the  package <cit.> for the halo model reaction-based boost


    Ξ(k,z) = ℛ(k,z)  P_ HMCode2020^ pseudo(k,z)/P_(k,z) ,


where P_(k,z) and P_ HMCode2020^ pseudo(k,z) are calculated using  <cit.> while ℛ(k,z) is calculated using .
We chose the  to model the pseudo-power spectrum as it has been shown to have improved accuracy and does not show suppression of power with respect to  (Ξ <1) at high redshifts, which is not expected. 



To train the emulator, we follow the procedure of SpurioMancini:2023mpt but widen the parameter priors. We produce ∼10^5 boost predictions in the range  k ∈ [0.01, 3] h  Mpc^-1 and z ∈ [0,2.5]. We take cosmologies from the Latin hypercube given by the ranges in <ref> and <ref>, with the massive neutrino  density parameter today's range being ∈ [0.0, 0.00317].  Emulation of the boost speeds up the computation by 4 orders of magnitude and we find similar accuracy of our emulator as found in SpurioMancini:2023mpt. The emulator is publicly available (https://github.com/nebblu/ReACT-emusgithub). Finally, we note another small difference between our emulator and that of SpurioMancini:2023mpt. In this work, we assume P_ in <ref> with  equal to the total of the true cosmology, whereas in SpurioMancini:2023mpt they assume  =  +,  and  being the baryon and cold dark matter density parameters today respectively. The emulators give the same output for = 0, which is what we assume in this work for Ξ. 





  §.§.§ 

The  emulator was introduced in Arnold:2021xtm. It is based on simulations for 50 combinations of |f_R0|, , σ_8^ and h with all other parameters fixed. We note that σ_8^ is the σ_8 we would obtain in a  model with the same cosmological parameters and initial amplitude A_ s and not σ_8 in an f(R) gravity model. The emulator accuracy is better than 2.5 % around the centre of the explored parameter space, up to scales of k = 10  h  Mpc^-1. f(R) simulations are performed by a modified version of the  code,  <cit.> that solves the nonlinear scalar field equation using a relaxation method. See  for more details. 

The emulation was made for the ratio between the power spectrum in f(R) and the  prediction in . We note that this is different from Ξ(k,z) as the power spectrum in  models in these simulations can have slight deviations from the  prediction. To account for this effect, the authors provided the ratio of the power spectrum in a reference  model to the  prediction. This  simulation uses =0.31315, h=0.6737, σ_8^=0.82172.
This can be used to obtain Ξ(k,z). However, this correction is provided only in this  model. Thus, the assumption here is that this correction is independent of cosmological parameters. The latest version of the  simulations has  counterparts to the f(R) simulations using the same seed. These simulations were analysed in . Hence, it is in principle possible to emulate directly Ξ(k,z) from these simulations. However, as this is not publicly available, we opted for using the original  emulator (https://bitbucket.org/arnoldcn/forge_emulator/src/master/bitbucket) for this paper. We will use one of these simulations for validation in <ref>. 





  §.§.§ 


We also consider the predictions given by the  emulator presented in <cit.>, which can predict the ratio Ξ(k,z) between the nonlinear matter power spectrum in f(R) gravity and in . 
The predictions are calibrated from a large suite of N-body simulations run with the code  <cit.>, a modified gravity version of the Adaptive Mesh Refinement (AMR) N-body code  <cit.>.
The  simulation suite covers the {|f_R_0|, , σ_8^} parameter space with 110 cosmological models sampled from a Latin hypercube <cit.>. 
The remaining  parameters, h, n_ s and  remain fixed, which means that their impact on Ξ(k,z) is ignored.
In saez-casares_2023, it was shown that the error made by this approximation is at the sub-percent level.
The quantity Ξ(k,z) is measured from pairs of f(R) and  simulations, both run with the same initial conditions, which leads to a large cancellation of cosmic variance and numerical resolution errors.
The emulation is done through a Gaussian Process regression <cit.>.
The emulator can give predictions for redshifts z ∈[0, 2], wavenumbers k ∈[0.03, 10] h Mpc^-1 and cosmological parameters in the following range: |f_R_0| ∈[10^-7, 10^-4], ∈[0.2365, 0.3941] and σ_8^∈[0.6083, 1.0140].
The estimated accuracy of , including emulation errors and systematic errors in the training data, is better than 3 % for scales k≲ 7 h Mpc^-1 and across the whole parameter space although in most cases, the accuracy is of order 1 %. The code is publicly available (https://gitlab.obspm.fr/e-mantis/e-mantisgitlab).






 §.§  implementation

Our implementation of the likelihood is based on the  implementation described in Euclid:2023pxu. For a more detailed explanation of the numerical implementation, we refer to that work. For our purposes, we have modified this code to include our prediction for Ξ(k,z) to obtain the nonlinear power spectrum in the f(R) gravity model. As the different implementations have different ranges in cosmological parameters, wavenumbers, and redshifts, we have chosen an extrapolation scheme to unify the ranges. We note that the exact implementation does not affect the final results strongly. We checked that the contributions from high redshifts and high wavenumbers which used extrapolations were sub-dominant.

For the modified gravity parameter |f_R0|, we use flat priors in terms of  to stay within the validity range of each model. For the largest scales beyond the range of the emulators, we set Ξ = 1. This is because the effect of the fifth-force has to vanish on these scales. For small scales, we do a power-law extrapolation. To obtain the spectral index for the extrapolation, we proceed as follows: we calculate Ξ on a grid close to the edges of the emulators; we then fit, for fixed z, a linear function in the log-log space onto this grid. This is done to average out any numerical noise at the edge. We do the same for high redshifts by fitting the power law for fixed wavenumbers. For regions of both high k and z, we do a constant extrapolation of the spectral index. To keep the extrapolated function Ξ(k, z) from going to non-physical values, we set a hard lower bound of Ξ = 1 and an upper bound of Ξ = 2. The ranges of the different emulators are found in <Ref>. We do a constant extrapolation in the  parameters. This can be done because during a typical MCMC the majority of the suggested points are within the ranges of our emulators. 
The validity range of  cosmological parameters is summarised in <ref>.







We also add the effect of baryonic feedback in the form of an additional correction Ξ^BFM to our power spectrum prediction. The physics of the baryonic feedback effects is discussed in <ref>. We refer the reader to Schneider:2019xpf and Mead:2020vgs for further details. We obtain the correction from  by Giri_2021. This emulator is only trained for redshifts below z=2 and wavenumbers k<12.5 h Mpc^-1. In this case, we use constant extrapolation for both k and z. We checked that this had little effect on our results. This is because the main probes are mostly sensitive to redshifts around z∼ 1. For these redshifts, the extrapolation only affects very high multipoles ℓ≳ 4500. Thus, the extrapolation has negligible effect on the angular power spectrum. We estimate the overall effect this has on our result to be at most at the percent level. In the absence of modified gravity, to obtain the nonlinear baryonic feedback power spectrum, we multiply the  nonlinear power spectrum by  Ξ^BFM. When adding the effect of modified gravity, we combine both boosts as

    P_f(R)^BFM(k,z) = Ξ^BFM(k, z) Ξ(k,z) P_(k,z) .

We can make this approximation if both effects are independent. This was shown to be the case for small deviations from  in f(R) gravity models considered in this paper using hydrodynamical simulations by Arnold:2019vpg and <cit.>. 

Our final addition to the code is the inclusion of theoretical errors. For this, we have adjusted the prescription of Audren_2013. The theory and numerical implementation is further discussed in <ref>.



§ COMPARISON OF DIFFERENT NONLINEAR MODELS

In this section, we compare four different predictions for the nonlinear dark matter power spectrum in terms of Ξ(k, z) defined in <ref> introduced in <ref>. We start from a comparison of the matter power spectrum with N-body simulations. We then compare the angular power spectra in the  reference cosmology. We perform MCMC simulations to compare errors and investigate bias due to the difference in the prediction of the nonlinear matter power spectrum using the settings defined in .



 §.§ Comparison of predictions



  §.§.§ Comparison with N-body simulations



<Ref> shows a comparison of the ratio of the power spectra between f(R) gravity and  measured from N-body simulations with the theoretical predictions for Ξ. 
These simulations use the same initial conditions and the ratio removes the cosmic variance and the effect of mass resolution. The left-hand side plot shows a comparison using the measurements from . This is based on the comparison project in Winther:2015wla. These simulations were run in a  cosmology with = 0.269, h = 0.704, n_s = 0.966 and σ_8 = 0.801. The simulations have N_p = 512^3 particles of mass M_ p≃ 8.756 × 10^9 h^-1 M_⊙ in a box of size B = 250 h^-1 Mpc and start at redshift z = 49. We picked a model called F5 with |f_R0| =  10^-5 and showed the result at z=0.667 that is presented in . In the comparison, we included the measurements from  and , as  is based on , while  is based on . The prediction of  agrees with  very well. On the other hand, the  prediction deviates from  as well as . We note that the  prediction is corrected using the  simulation (Node 0) in the  simulation suite with =0.31315, h=0.6737, σ_8^=0.82172 to obtain Ξ. We find that if we use =0.31315 and h=0.6737 in the  prediction, the agreement with  is much better. 

To further investigate this issue with , we use the measurement of the power spectrum in one of the nodes in the  simulation suite (Node 13) run by  that is closest to the  reference cosmology that we will use in this paper with non-zero |f_R0|. This simulation has the following cosmological parameters: = 0.34671, h = 0.70056 and σ_8 = 0.78191 and we show a comparison at z=0.652. Both  and f(R) simulations with |f_R0| = 10^-4.90056 are available so that we can measure Ξ directly. We note that the pipeline developed in  measures the power spectrum only up k =3 h  Mpc^-1. In this case, both   and  agree with  within 1 % up to k =3 h  Mpc^-1. The fitting formula and  agree with  within 1 % up to k =1 h  Mpc^-1. Given this result, the large discrepancy between  and  is likely to be caused by emulation errors as well as calibrations using the   power spectra predictions.  



  §.§.§ Comparison in the  reference cosmology

In this paper, we consider the model called HS6 in Euclid:2023tqw, which has the following parameters 
    Θ_ fid   ={, ,  h,  n_ s, σ_8, } ,
    HS6:Θ_ fid   ={ 0.32,  0.05,  0.67,  0.96,  0.853,  -5.301} .

The cosmological parameters are the same ones adopted in . 
Euclid:2023tqw showed that this value of |f_R0| = 5 × 10^-6 can be well constrained by the photometric probes. Also the range of |f_R0| covered by the four models is wider than the errors predicted in Euclid:2023tqw. Our fiducial cosmology includes massive neutrinos with a total mass of ∑ m_ν=0.06 eV, but we keep ∑ m_ν fixed in the following analysis. 

<Ref> shows a comparison of the power spectrum and the angular power spectrum for WL,  and their cross-correlation . In these plots, we show the ratio to the  prediction and error bars from the diagonal part of the covariance matrix. For the power spectrum comparison at z=0.5, we see that  and  agree best. This is not surprising as these two predictions are based on N-body simulations ran by the same code . On the other hand,  overestimates Ξ at k= 0.1 h  Mpc^-1. This is similar to the deviation we find in the comparison with the N-body simulation from  although the deviation is smaller at the 2 % level. 
This is likely due to the fact that  in HS6 is closer to  in the fiducial   simulation (=0.31315). On the other hand,  underestimates Ξ at k < 0.1 h  Mpc^-1 and even predicts Ξ <1. As we discussed in the previous section, we enforce Ξ≥  1.  





 §.§ Forecasting errors and biases
 


  §.§.§ Forecasting errors for WL

We first compare errors obtained by running MCMC using the synthetic data created by one of the four nonlinear models and fitting it by the same model. In this case, by definition, we recover the input parameters that were used to create the synthetic data. Our interests are constraints on the |f_R0| parameter and cosmological parameters. We first consider the WL-only case. In this case, we impose a tight Gaussian prior on the spectral index, n_ s, taken from the Planck results <cit.> and on the baryon density parameter, using Big Bang nucleosynthesis constraints <cit.>: 

    n_ s = 0.96 ± 0.004,    h^2 = 0.022445 ± 0.00036 ,

as we do not expect to obtain strong constraints on these parameters from WL alone. The parameters that are used in the MCMC runs are summarised in <Ref>, including their fiducial values and prior ranges. As a convergence criterion, we use a Gelman–Rubin <cit.> value of R-1 < 0.01 for each individual sampling parameter using . For post-processing chains, we use 
 
<cit.>. 



<Ref> shows the 2D contours of the constraints on parameters,  
and <Ref>
and <Ref> 
summarise constraints on  for the optimistic and pessimistic settings. The constraints on cosmological parameters are consistent among the four different nonlinear models, while we see some notable differences in the constraints on |f_R0|. In the case of the optimistic setting,  gives the tightest constraints, which is also closer to Gaussian. The fitting formula agrees with  for small  but has a weaker constraint for larger . Constraints from  agree with the fitting formula for large , but give a weaker constraint for small . The degeneracy between , h^2 and ln(10^10 A_ s) also presents some notable differences. The degeneracy for larger  is different for the fitting formula when compared with  and . This could be attributed to the fact that the fitting formula does not include any cosmological parameter dependence in the prediction for Ξ(k,z). We observe similar agreements and disagreements for the pessimistic setting, but the agreement among ,  and the fitting formula is better, particularly for small .  gives a weaker constraint on , but the constraints on cosmological parameters are consistent among the four different models. This is due to the weak cosmology dependence of Ξ(k,z), as discussed in <cit.>, and the fact that the constraint on cosmological parameters is coming from the  power spectrum, which is common in all these four models. 

We compare the constraints from the MCMC analysis with the Fisher Matrix forecast and show this in <ref>. The  pipeline can be used as a Fisher Matrix forecast tool <cit.>,  which was shown to agree very well with previous Fisher Matrix forecasts given by . We note that in Euclid:2023tqw, no prior was imposed on n_ s and h^2, thus no direct comparison is possible. Instead, we use  as a Fisher Matrix forecast tool and compare the result with the MCMC analysis to validate the Fisher Matrix forecast. To be consistent with Euclid:2023tqw, we use the fitting formula as the nonlinear model. The 1 σ error is very consistent: the Fisher Matrix forecast gives 0.111 while the MCMC analysis gives 0.116. The constraint from MCMC is non-Gaussian and the 1D posterior has a slightly longer tail for large . The constraints on cosmological parameters agree very well between the Fisher Matrix forecast and the MCMC result. 



  §.§.§ Assessing biases for WL

Next, we create the synthetic data using  and we fit it by different nonlinear models to assess the biases in the recovered parameters due to the difference in the nonlinear modelling. We selected  to create the data because it has a relatively narrow prior range for |f_R0| and we encountered a problem with that range when using  as the model to fit when including baryonic effects. We should note that 
has a larger discrepancy from other nonlinear models. As we discussed before, this could be attributed to emulation errors and calibration with . In particular, this choice is disadvantageous to  as the difference of the prediction for Ξ from  is the largest. 
Thus, the estimation of bias in recovered parameters presented in this section is conservative, particularly for  

<Ref> shows the 2D contours of the constraints on parameters,  
and <Ref>
and <Ref> 
summarise constraints on  for the optimistic and pessimistic settings.

We first start from the optimistic setting. Since  is valid only up to 
k = 3 h  Mpc^-1, we will not include  in this case. 
The fitting formula and  recover the input parameters within 1σ. In the case of , cosmological parameters are well recovered, but there is a slight bias in the recovered . On the other hand, a slight bias appears in h in the case of the fitting formula. This may be attributed to the fact that the fitting formula does not include cosmological parameter dependence in the prediction of Ξ. To quantify the bias, we define a 1D bias as 

    B_ 1D= μ - μ_/σ_ ,

where μ and σ are the mean and 1 σ error computed from the 1D marginalised posteriors. If μ > μ_ (μ < μ_), we use the upper (lower) 68.3% confidence interval to obtain σ_. 
The 1D bias for  is B_ 1D=0.273 and 0.602 for the fitting formula and , respectively. 

We see a similar result in the pessimistic setting for  and the fitting formula. In this case, the input  is well within 1σ although h is again slightly biased. The 1D bias for  is B_ 1D=0.441 and 0.518 for the fitting formula and , respectively while 
the 1D bias for h is 0.988 and 0.812 for the fitting formula and , respectively. 
On the other hand, we find a bias in the recovered  when  is used. As mentioned above, we should bear in mind that the choice of  as a fiducial is disadvantageous for . In the case of  , the difference from  leads to biases also in the cosmological parameters leading to larger h^2 and smaller 
ln(10^10 A_ s) and h. This is due to the different k dependence of Ξ between  and  and this is compensated by adjusting cosmological parameters as well as  leading to a stronger bias. The 1D bias for  reaches B_ 1D=3.11. We will discuss how to mitigate this bias by correcting the prediction for the fiducial cosmology and including theoretical errors in <ref>. 

Since the  prediction deviates from the other three models, we also performed the same analysis using  as fiducial data for the pessimistic setting. We show the results in Appendix A. Qualitatively, we obtain consistent results. As we can see from <ref>, the 2D contours overlap in the same way as the case where FORGE is used as fiducial, while the mean of  obtained by  is closer to the input value. However, we still find that the 1D bias is at 3 sigma level due to the smaller errors of  compared with . 













  §.§.§ 3x2pt  analysis

We consider the constraints from  statistics by adding  and its cross-correlation  to WL. In this paper, we use a simple model of the scale-independent linear bias as in . This assumption needs to be reexamined in f(R) gravity as the scale-dependent growth will lead to a scale-dependent bias. However, the linear bias assumption also needs to be relaxed even in  and it is beyond the scope of this paper to implement a more complete bias description. For this reason, we consider the pessimistic case only and focus our attention on a comparison with the Fisher Matrix forecast and check if the bias when  is used to create the data becomes worse with the increasing statistical power. We also do not include  as we already observe a significant bias with WL only. 

We have 10 scale-invariant bias parameters, one for each redshift bin. For the  analysis, we vary n_ s, but we still impose a tight Gaussian prior on h^2
as it is unlikely to get a tighter constraint than the Big Bang nucleosynthesis constraints.

The parameters used in addition to the WL analysis are listed in <Ref>.
We did not change the fiducial bias parameters from   as there is no prediction of the linear bias in  f(R) gravity. The fiducial bias adopted in  will need to be improved even in . Our prime focus is the study of the effect of nonlinear models and we vary the linear bias in the MCMC analysis resulting in different constraints depending on the nonlinear model used. The observed covariance is built from synthetic data. Thus the effect of f(R) gravity is included in the covariance. We independently checked that the constraints on parameters do not change by using the covariance. Thus we expect that the effect of f(R) on galaxy bias has also negligible effects on the covariance. 

The left-hand side of <ref> shows the 2D contours of the constraints on parameters, and <Ref>
summarises constraints on  where the same nonlinear model is used for the data and the model. 
The right-hand side of <ref> shows the 2D contours of the constraints on parameters, and <Ref>
summarises constraints on  as well as 1D bias for  when  is used to create the data. 

In the case where the same nonlinear model is used for the data and the model, errors are consistent between the fitting formula and , although we still observe the same longer tail for larger  for the fitting formula as we observe for WL.  gives slightly tighter constraints on . Constraints on cosmological parameters are very consistent among the three different nonlinear models. The Fisher Matrix forecast 
using the fitting formula is also very consistent with the MCMC results. We note that errors from  analysis in the pessimistic setting are comparable to or better than WL alone in the optimistic setting. This demonstrates the strength of the  analysis although we should bear in mind the limitation of the bias model used in this analysis. 

We find that the increased statistical power does not degrade the agreement between  and . 
When  is used to create the data, the input parameters are well recovered by , although h is slightly biased as in the WL-only case. The 1D bias for  is given by 0.150. In the case of the fitting formula, the bias in cosmological parameters becomes worse compared with WL. This could be attributed to the fact that the fitting formula does not take into account the cosmological parameter dependence of Ξ as mentioned before. The 1D bias for  also becomes slightly larger, with a value of B_ 1D =0.667. 











§ BARYONIC EFFECTS
 





 §.§ Adding baryonic effects using 

In this section, we study the impact of baryonic effects on the constraints on the f(R) parameter. We use the seven-parameter emulator of baryonic effects called  <cit.> and we assume that the baryonic effects and the modified gravity effects can be treated independently as discussed in <ref> <cit.>. The  parameters govern the gas profiles and stellar abundances in haloes. It is not our intention to study in detail the baryonic effects. We use the default values for the baryonic parameters provided with  and use the full prior range. We choose not to include any redshift dependence. This is because we need to impose a tight prior range for these parameters at z=0 to keep them within the prior range and we may miss important degeneracies between baryonic parameters and . 
We found that two baryonic parameters, M_ c^ Giri et al. and θ_ ej, have the strongest effects on the matter power spectrum and they are well constrained in the presence of the f(R) parameter. Thus, we only vary these two parameters.
In , the gas profile is modelled as a cored double power law. M_ c^ Giri et.al. controls the dark matter halo mass dependence of the logarithmic slope of the first-cored power law. It allows the profile to become less steep than the Navarro_96 one for M<M_ c^ Giri et.al.. The parameter θ_ ej determines the scale radius (with respect to the virial radius) of the second-cored power law. <Ref> summarises the fiducial values and priors for baryon parameters. We will use the dimensionless parameter M_ c that is related to the one defined in <cit.> via M_ c M_ c^ Giri et.al. / M_⊙ where M_⊙ is the solar mass. 


















 §.§ Effects of adding baryons

<Ref> shows the effect of changing , h and two baryon parameters on the ratio of the WL angular power spectrum to the  one. Baryonic effects introduce scale-dependent modifications to the angular power spectrum that are similar to the effect of f(R) at ℓ > 100. On the other hand, h changes the overall amplitude of the angular power spectrum as it changes . We will see that the interplay between these parameters leads to interesting degeneracies. 

<Ref> shows the 2D contours of the WL constraints on parameters for the optimistic settings with baryons and <Ref> summarises constraints on . 

In the case of the optimistic setting, when  is used as the data and the model, the two baryonic parameters degrade the constraints on : the 1σ error becomes 0.281 while it was 0.088 without baryons. However, it is important to stress that  and the two baryonic parameters are still constrained well within the prior range of these parameters. It means that the impact of modified gravity on the high-k tail of the matter power spectrum is not washed out by the baryonic feedback, and we can still distinguish between the effect of f(R) and baryons. 

Nevertheless, if  is used as the data, the difference between  and  in the matter power spectrum can be absorbed by the shift in M_ c and h, and the constraint on  is shifted along the degeneracy direction between  and M_ c as well as h. This can be understood from <ref>. The difference of the scale dependence between  and  at ℓ >100 can be adjusted by decreasing  and decreasing M_ c. This leads to a lower amplitude that can be adjusted by decreasing h. Due to the combination of these effects, the best-fit  becomes smaller. We note that 
 does not include the h dependence in the prediction for Ξ(k ,z) explicitly, but it depends on h implicitly through σ_8^. However, this dependence of h on Ξ(k,z) is very weak. To compute the power spectrum in f(R) we use 
P_f(R)
= Ξ(k,z) P_(k,z). 
The  emulator provides us with the first factor, Ξ(k,z), and a  emulator provides the last factor. The h dependence comes from the  power spectrum not from Ξ(k,z).

As a result, we obtain a 95.5% upper limit of  as < -5.477, which is incompatible with the input value of = -5.301. This is partly due to prior volume effects shifting the contour to lower values of , caused by the strong degeneracy between  and M_ c as well as h.  

To confirm this, we obtain a profile likelihood for  by fixing  and finding a minimum χ^2 by varying other parameters. The profile likelihood was obtained by  
<cit.>. This is shown in <ref>.  
The global best-fit value for  is obtained as =-5.765, which is  larger than the mean of the 1D marginalised posterior =-6.087. <Ref> also shows the prediction of   with the best-fit values, which agrees well with . The Δχ^2 for the input value of = -5.301 is found as Δχ^2 = 2.657, thus it is still within 2 σ. We also observe that the Δχ^2 curve becomes flat for smaller values of . 

We also tested the fitting formula as a model to fit to the data generated by . In this case, the posterior is highly non-Gaussian and the chains did not converge well. 

This implies that we can break the degeneracy between f(R) gravity and baryonic effects, but we need to have an accurate nonlinear model for the f(R) gravity model. Otherwise, we could obtain a significantly biased result in terms of the 1D marginalised constraint. We note that if we have a physical prior on M_ c from baryonic physics, we could break the degeneracy between  and M_ c. For example, in the case where  was used as the data and they were modelled by , the mean of M_ c is significantly shifted to a smaller value as we can see in <ref>. If we had a prior on M_ c to prevent this, this bias could be avoided. 
To test this idea, we also ran an analysis imposing a prior on M_ c and showed the result in <ref>. We used a Gaussian prior with the width of 0.2, which was estimated from the weak lensing informed gas and stellar mass fraction measurements of massive haloes by <cit.>. This prior is consistent with the error on M_ c that we obtained by using e-Mantis both for data and model. We observe that the bias is relaxed slightly, however, the degeneracy between M_ c and  is quite strong and this prior is not enough to alleviate the bias in the marginalised constraints. An improved prior from external data is needed to fully break this degeneracy.

We assumed that the effect of f(R) gravity and baryonic effect can be treated independently. 
Recently, the coupling between baryonic feedback and cosmology has been studied showing that the combined effect of baryonic and non-baryonic suppression mechanisms is greater than the sum of its parts for decaying dark matter <cit.>. The effect of this coupling on the degeneracy needs to be studied in the future. 


Next, we consider the pessimistic setup. <Ref> shows the 2D contours of the WL constraints on parameters for the optimistic settings with baryons and <Ref> summarises constraints on . In this case, due to the weaker constraining power, the 99.7% confidence level upper bound on  reaches the prior boundary of =-7. Also, the 95.5% confidence level upper bound of θ_ ej is bounded by the prior. If the data are generated by , strong degeneracies appear among  and two baryonic parameters. The input  is consistent within 1σ for the fitting formula and  although the prior bound  =-7 is reached at 95.5% confidence level lower bound. On the other hand, for , the input value is outside the 95.5% confidence level upper bound. As we already observed without baryons, h is biased to a lower value. 

We also checked the case where fiducial data in generated by . As we can see in <ref>, unlike the case without baryons, the agreements between ,  and  are improved significantly compared with the case where  is used for fiducial data. Particularly, for , the mean of  is now consistent with the input value within 1σ. 



 §.§ Constraints on fR0 with the  data

Based on these results, we consider the pessimistic setting and use  to obtain an upper limit of  in the presence of baryonic effects to be conservative. 

<Ref> shows the 2D contours of the WL constraints on parameters for the pessimistic settings where  is used to create the data. 
We find that the recovered cosmological parameters are consistent with the input values. The mean of 1D marginalised constraints of baryonic parameters are slightly biased due to the strong degeneracy between them, but it is still consistent within 1 σ. We obtain the 95.5% confidence level upper bound on  as 

    < - 5.21 .

We note that this bound depends on the prior =-7. To obtain the prior-independent bound, we follow the approach presented in Gordon:2007xm, Piga:2022mge and <cit.>. 
The ratio of the marginalised posterior and prior is given by 

    b(x | d, p) = 𝒫(x | d, p)/p(x) ,

where x is the parameter we are interested in (i.e. ), d is the data and p is the prior, and 𝒫 is the Posterior. The Bayes factor B(x_1, x_2), which quantifies the support of the models with x=x_1 over the models with x=x_2 is given by 

    B(x_1, x_2) =  b(x_1 | d, p) /b(x_2 | d, p) 
    = ℒ(d | x_1)/ℒ(d| x_2),

where ℒ(d| x) is the marginalised likelihood of the data for parameter x. 
For our purpose, we choose x = and fix x_1 to be the upper bound of the prior, x_1 = -7. 
Following Gordon:2007xm, we then use B(x_1, x_2) =2.5  to find x_2 so that the model with x=x_1 is favoured compared to the model with x=x_2 at 95.5% confidence level. This gives the 95.5% confidence interval of  that does not depend on the prior. We note that this method applies only to a ℒ(d| x) that is a monotonic function of x. Employing this technique, we obtain

    < - 5.58 .


Finally, in <ref>, we present the constraints on  and   so that we can compare these with those in the literature <cit.>. It is not possible to make a direct comparison due to various differences in the setting. Nonetheless, our constraint is comparable to the one presented in Schneider:2019xpf [< -5.3] where a similar analysis was done by combining  and the fitting formula. We cannot also make a comparison with Harnois-Deraps:2022bie as baryonic effects are not included in their analysis, but again the constraint is comparable [< -5.24]. 
On the other hand, our constraint is much weaker than the one found in SpurioMancini:2023mpt that used  with . We refer the readers to SpurioMancini:2023mpt for possible explanations. In our analysis,  tended to prefer a lower value of  when  is used as data and this might be one of the reasons for this difference. Finally, a recent study by Tsedrik:2024cdi used similar settings to forecast cosmic shear constraints on model-independent parametrisations of modified gravity theories with scale-independent linear growth. They found that a much better understanding of baryonic feedback is needed in order to detect a screening transition.








§ THEORETICAL ERROR


Any theoretical prediction will in general come with an associated error and should be included in the likelihood if it is not subdominant to the statistical error. In this section, we discuss the implementation of theoretical errors in our pipeline. To be conservative, we will estimate the theoretical errors using  and  as these two give the most discrepant results. We will apply this method to 
 for the pessimistic case of WL to check if we can remove the parameter biases that we observe for the data created by  at the expense of enlarging error bars.



 §.§ Adding uncorrelated theoretical errors to the likelihood

The implementation presented in this work is based on the work of Audren_2013. We have adjusted the recipe to the full  likelihood, although we only show the result for the WL analysis in this paper. Here we will discuss the idea behind the formulation. We calculate the angular power spectrum from the power spectrum using <ref>. If there is some uncertainty in the modelling of the nonlinear power spectrum, this propagates to the angular power spectrum. As <ref> is a linear functional of the nonlinear matter power spectrum, ℱ[P_δδ], we can propagate the error on the power spectrum, Δ P_δδ, and find 
    ℱ[P_δδ+Δ P_δδ] = ℱ[P_δδ] + ℱ[Δ P_δδ]  C_ij^XY(ℓ) + E_ij^XY(ℓ) ,

where we have defined the angular power spectrum error E_ij^XY(ℓ). Following Audren_2013, we define the relative error on the power spectrum

    Δ P_δδ(k,z) α(k,z)   P_δδ(k,z) .

We take the most conservative approach and assume that the error is uncorrelated between different values of ℓ.
To account for this uncorrelated theoretical error, we add a new nuisance parameter ε_ℓ for each multipole and define the shifted covariance as 
    C^XY_ij(ℓ) Ĉ^XY_ij(ℓ) + ε_ℓ  L^1/2  E^XY_ij(ℓ) .

The free parameters ε_ℓ have the following meaning: they quantify how many standard deviations the shifted covariance matrix is from the theoretically predicted one. We treat them as a random Gaussian variable with zero mean and a standard deviation of one. The normalisation 
factor L^1/2 where L=ℓ_max-ℓ_min+1 will become clear soon. The log-likelihood function becomes 
    χ̃^2 (ε_ℓ) =  ∑_ℓ =ℓ_min^ℓ_max[ (2ℓ+1) f_ sky( d̃^ mix_ℓ(ε_ℓ)/d̃^ th_ℓ(ε_ℓ) + lnd̃^ th_ℓ(ε_ℓ)/d^ obs_ℓ-N) + ε_ℓ^2 ].

The quantities d̃^ th and d̃^ mix are constructed in the same way as d^ th and d^ mix 
[see eq:ds1eq:ds3] using the shifted covariance matrix C^XY_ij(ℓ). To include the theoretical error, we vary ε_ℓ and marginalise over them. To a very good approximation, this is equivalent to minimising the χ̃^2 with respect to ε_ℓ at the level of the likelihood.  Thus, we define our new χ^2 as the minimum of the χ̃^2 with respect to ε_ℓ:

    χ^2 min_ε_ℓ∈ℝ^Lχ̃^2.

The normalisation factor L^1/2 can be explained as follows. If we were to measure Ĉ^obs(ℓ) = Ĉ^th(ℓ) + E(ℓ), the minimisation would find ε_ℓ = L^-1/2. The resulting χ^2 = ∑_ℓε_ℓ^2 = ∑_ℓ L^-1 =1 would match our expectation that a one-sigma theoretical error for each ℓ results in an increase of the χ^2 by one. 

The main ingredient of this formulation is the relative error function α(k,z) defined in <ref>. 
The construction of α(k,z) using  and  will be discussed in <Ref>. 



 §.§ Numerical implementation

The numerical computation of the theory error covariance E^XY_ij(ℓ) follows the prescription presented by Euclid:2023pxu. For the minimisation, we can use the fact that all the free ε_ℓ are independent of each other, and we can do the minimisation for each multipole separately. We use Newton's method to find the minimum. For this, we have to compute the first and second derivatives of the likelihood with respect to ε_ℓ. For any single multipole, we find 
    dχ̃/dε_ℓ = (2 ℓ+1) f_ sky ((d̃^mix_ℓ)^' + (d̃^th_ℓ)^'/d̃^th_ℓ - d̃^mix_ℓ (d̃^th_ℓ)^'/(d̃^th_ℓ)^2) + 2 ε_ℓ .

The derivatives of the determinants are computed using Jacobi's formula. This gives, for example

    (d̃^th_ℓ)^' = ( Ĉ_ℓ^th) L^1/2 [ (Ĉ_ℓ^th)^-1  E_ℓ] ,

and a similar expression for (d̃^mix_ℓ)^'. We compute the inverse of the covariance numerically and obtain the second derivatives numerically from the first derivative by doing a double-sided three-point stencil. The minimisation would then need to be done for each multipole. For the pessimistic settings, this would correspond to a minimisation in a 1500-dimensional parameter space. To save time, we only do the minimisation on a logarithmically-spaced grid with 100 discrete values. The other values are obtained from an interpolating function. We checked that the effect of the interpolation does not change the results by more than 1%. The obtained ε_ℓ tend to vary continuously with ℓ, as they try to mimic the effects of changing other theory parameters to the observed power spectrum.








 §.§ Application to 

We first apply this method to the case where the data are created by  and the parameter fitting is done by . Due to the significant bias in the recovered parameters for both |f_R0| and cosmological parameters, we find that the inclusion of theoretical errors is not enough to mitigate the bias. Therefore, we additionally correct the prediction of  by  for the fiducial cosmology as 

    Ξ_
        = Ξ_(
    Ξ_/Ξ_)_ fiducial.

The constraints on cosmological parameters are shown in the left panel of <ref>. 
In this case, as expected, we recover the input parameters in an unbiased way when the data are created by . We see that theoretical errors affect mainly  and h.  In order to check a non-trivial case, next we create the data with . The result is shown in the right panel of <ref>. 
Without theoretical errors, the means of 1D marginalised constraints for , h^2 and ln(10^10 A_ s) are slightly biased compared with the input values. The inclusion of the theoretical errors largely resolves these biases not only by enlarging the error bars but also by making the means closer to the input values. We note that the inclusion of theoretical errors in this case is important to justify the rescaling of Ξ in <ref>. We also note that the theoretical error included here overestimates the errors significantly now that we corrected the prediction of  by  for the fiducial cosmology. 

Finally, we study the impact of adding baryonic effects in the presence of theoretical errors. We again use  as the data. The result is shown in <ref>. Due to the enlarged errors, constraints on parameters are not affected significantly by baryonic effects. Also, the means of the 1D marginalised constraint remain consistent with the input parameters. However, the inclusion of the theoretical errors changes the degeneracies between cosmological parameters, baryon parameters and . This leads to a tighter lower bound on  compared with the case where the same data are fitted by  itself. This reinforces our conclusion that the 1D marginalised constraint on  is sensitive to degeneracies among parameters, and the difference in the theoretical predictions strongly affects the constraint. 

The result shown here is just an illustration of the inclusion of theoretical errors and their effects on the parameter constraints. Our implementation is flexible and it can be applied to any theoretical error described by the relative error function α(k,z) defined in <ref>. 



§ CONCLUSION


In this paper, we studied the effect of using different nonlinear predictions for the dark matter power spectrum on the parameter constraints in the Hu–Sawikci f(R) gravity model obtained from primary photometric probes. We implemented four different models in the  pipeline to predict angular power spectra for weak lensing (WL), photometric galaxy clustering () and their cross-correlation (). Comparing with the N-body simulation data obtained in , we found that  agreed very well with  that was used to run simulations to construct the emulator, while  had larger errors compared with the  simulation. The agreement is better for one of the N-body simulations used to construct  obtained in  run by . This indicates that the difference between  and  is larger than the one in the baseline N-body simulations (i.e.  and ) mainly due to the way  was constructed. In the  reference cosmology,  gives a larger Ξ, the ratio between the power spectrum in f(R) and in , on all scales compared with  and the fitting formula at the 2 % level.  underestimates Ξ compared with  more than  and the fitting formula at intermediate k. 

We used the fiducial value of |f_R0| = 5 × 10^-6 (= - 5.301) and ran MCMC in the  fiducial cosmology defined in . For the fitting formula, the Fisher Matrix forecast and MCMC results generally agree well although the MCMC result is non-Gaussian. This is partly caused by a lack of cosmological parameter dependence of Ξ, which affects the degeneracy between  and cosmological parameters for large  .  gives more Gaussian constraints with smaller errors. When  is used to create the data, the 1D mean of  is not strongly biased in the case of  and the fitting formula and the 1D bias is at most 0.6 σ. Even for the  analysis including all the probes and their cross-correlations, the 1D bias is 0.15 σ for  in the case of . 

The impact of baryonic physics on WL was studied by using a baryonification emulator . For the optimistic setting, the f(R) parameter and two main baryon parameters are well constrained despite the degeneracies among these parameters. However, the difference in the nonlinear dark matter prediction can be compensated by the adjustment of baryon parameters as well as cosmological parameters, and the 1D marginalised constraint on  is biased. This bias can be avoided in the pessimistic setting at the expense of weaker constraints. For the pessimistic setting, using the  synthetic data for WL, we obtained the prior-independent bound of < -5.6 using . 

 shows a large bias in  as well as cosmological parameters when  was used to create the data. 
This is because the prediction of  is furthest away from . We implemented a method to include uncorrelated theoretical errors proposed in Audren_2013 to address this issue. The method is based on the relative error function for the nonlinear dark matter power spectrum. We estimated this using the difference between  and . We found that the inclusion of theoretical errors alone was not enough to mitigate the bias. We then corrected the prediction of  with  for the fiducial model. We applied this model to the data created by . We found that theoretical errors, in this case, helped reduce the bias not only by enlarging errors but also by making the means of 1D marginalised constraint closer to the input values. When we added baryonic effects with , errors were not significantly affected and the input values were still recovered. However, the lower bound on  is tighter than the one obtained by applying  itself. This reinforces our conclusion that the 1D marginalised constraint on  is sensitive to the degeneracies among parameters. 

Based on the result of this paper, we draw the following conclusions:



  * It is important to check the agreement of different N-body codes that are used to create theoretical predictions. This is not only the code itself, but also accuracy settings. In , we found that the accuracy setting such as the refinement criteria in the adaptive mesh refinement method has a large effect on the power spectrum. With the controlled accuracy setting, it is possible to realise 1 % agreements between different N-body codes in terms of Ξ(k, z) in the Hu–Sawicki f(R) gravity model. 


  * We then need to check the accuracy of emulators. We found that  suffers from larger emulation errors and this leads to a larger difference between  and  compared with the difference in their baseline N-body code  and  in some cosmologies. Improving the emulation technique will be able to make the agreement better. 


  * Including baryonic effects can worsen the bias in the 1D marginalised constraint on . This is because the difference in the nonlinear dark matter power spectrum prediction can be compensated by the adjustment of baryon parameters, and the best-fit values are biased. In addition, the degeneracy between baryon parameters, cosmological parameters and  leads to a stronger volume effect. This bias can be avoided if there is a prior on baryon parameters from external data sets for example. 


  * To account for the uncertainty in the theoretical prediction for the nonlinear power spectrum, it is safer to include theoretical errors. In this paper, we used a conservative error estimation using  and . An improvement in emulators will make the theoretical error smaller leading to better constraints. However, we still need to check whether or not the inclusion of baryons will change this conclusion. 



The pipeline developed in this paper can be used to test the readiness of the nonlinear power spectrum prediction for the application to real data from in other extended cosmologies. For the Data Release 1 of and future data releases, we plan to improve the 
 emulator and make the pipeline ready for the data analysis. We also plan to extend the analysis to the models studied in Euclid:2023rjj. 

Finally, we note that our forecasts have ignored observational systematics such as shear and redshift measurement biases. For Stage-IV surveys like , it is known that these need to be very well characterised in order to get accurate cosmological parameter estimates. This is a challenge that needs to be addressed for both ΛCDM and exotic cosmologies. A detailed study of this issue is beyond the scope of our paper, but we refer to EuclidSkyOverview for details. In summary, the mean galaxy redshifts within the bins need to be known with an accuracy better than ∼ 0.002(1+z) if errors in cosmological parameters are not to be degraded. On the other hand, <cit.> incldued a parameter that shifts the mean of the redshift distribution in each redshift bin and found that these parameters are not strongly degenerate with the f(R) parameter.


K. K. is supported by STFC grant ST/W001225/1.

B. B. is supported by a UK Research and Innovation Stephen Hawking Fellowship (EP/W005654/2). 
A. P. is a UKRI Future Leaders Fellow [grant MR/X005399/1]. P. C. is supported by grant RF/ERE/221061.
M. C. acknowledges the financial support provided by the Alexander von Humboldt Foundation through the Humboldt Research Fellowship program, as well as support from the Max Planck Society and the Alexander von Humboldt Foundation in the framework of the Max Planck-Humboldt Research Award endowed by the Federal Ministry of Education and Research. 
F. P. acknowledges partial support from the INFN grant InDark and the Departments of Excellence grant L.232/2016 of the Italian Ministry of University and Research (MUR). FP also acknowledges the FCT project with ref. number PTDC/FIS-AST/0054/2021.

GR’s research was supported by an appointment to the NASA Postdoctoral Program administered by Oak Ridge Associated Universities under contract with NASA. GR was supported by JPL, which is run under contract by the California Institute of Technology for NASA (80NM0018D0004).

Numerical computations were done on the Sciama High Performance Compute (HPC) cluster which is supported by the ICG, SEPNet and the University of Portsmouth.
This work has made use of the Infinity Cluster hosted by Institut d'Astrophysique de Paris.







We acknowledge open libraries support  <cit.>,  <cit.>,  <cit.>, and  <cit.>.

For the purpose of open access, we have applied a Creative Commons Attribution (CC BY) licence to any
Author Accepted Manuscript version arising.

Supporting research data are available on reasonable request from the corresponding author.











§  FIDUCIAL DATA


In this Appendix, we show the result of the cases where we use  as fiducial data for the pessimistic setting of WL with and without baryons. <Ref> and <ref> show the results without baryons while the <Ref> and <ref> show the results with baryons. We observe that the agreements of the three models are generally better as expected although without baryons, 1D bias is still at the 3σ level for  due to smaller errors of . On the other hand, in the case with baryons, 1D bias is reduced to less than 1σ. 













§ CONSTRUCTION OF THE RELATIVE ERROR
 
To construct the relative error function, α(k,z), we can use the fact that Ξ(k, z) is constructed as the ratio to the same  nonlinear power spectrum calculated with the Halofit `Takahashi' prescription in all models. This means that any deviations of the  power spectrum from N-body simulations is also absorbed in Ξ. We can thus write the relative error as

    α(k,z) = Δ P_δδ(k,z)/ P_δδ(k,z) = ΔΞ(k,z)/Ξ(k,z) .

To construct Δ P_δδ we follow the basic idea that all predictions for Ξ are equally accurate in calculating the modified-gravity power spectrum. The true power spectrum is thus only known up to the spread of the predictions. To be conservative with our forecast, we choose to calculate the error from the difference between the  predictions and the  predictions, which is the largest among these predictions. This is shown in <ref>. The  prescription underpredicts the power spectrum at intermediate scales and overpredicts slightly at smaller scales. This is most likely due to  handling the one- and two-halo power spectra separately.  predictions for the pseudo-cosmology non-linear power spectrum also contribute to these inaccuracies as we can see by comparing the left and right panels of Figure 4 in Cataneo:2018cic. 
To simplify the exact difference between these two codes, we argue that the theoretical error should plateau at the maximum deviation. This is because at nonlinear scales, the modes of different wave numbers are no longer independent of each other. By plateauing, we essentially say that we can no longer be more precise in the computation of the power spectrum on smaller scales. 

Finally, we construct α using only the difference of the predictions at a fiducial of |f_R0| = 5 × 10^-6 in order not to bias our results. Our final fit has the form

    α(k,z)    = A(z) x^2+x/x^2+x+1 ,
    
        A(z)    = A_1/exp(z-A_2/A_3)+1+A_4 , 
    
        k_ p(z)     = B_1  exp[tanh(z-B_2/B_3)] ,
 
where we use x  k / k_ p. We separately fit the amplitude function A(z) to the maximum deviation and the plateau wavenumber k_ p(z) to the wavenumber of maximum deviation. The best-fit values can be found in <Ref> and are shown as the shaded region in <ref>. We see that the zero line is within the 68.3% confidence bounds most of the time. This is in accordance with our construction that the two emulators should differ by the theoretical error.









