Euclid preparation. Simulations and nonlinearities beyond $Λ$CDM. 1. Numerical methods and validation
Euclid Collaboration J. Adamek B. Fiorini M. Baldi G. Brando M. -A. Breton F. Hassani K. Koyama A. M. C. Le Brun G. Rácz H. -A. Winther A. Casalino C. Hernández-Aguayo B. Li D. Potter E. Altamura C. Carbone C. Giocoli D. F. Mota A. Pourtsidou Z. Sakr F. Vernizzi A. Amara S. Andreon N. Auricchio C. Baccigalupi S. Bardelli P. Battaglia D. Bonino E. Branchini M. Brescia J. Brinchmann A. Caillat S. Camera V. Capobianco V. F. Cardone J. Carretero S. Casas F. J. Castander 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 X. Dupac S. Dusini M. Farina S. Farrens S. Ferriol P. Fosalba M. Frailis E. Franceschi M. Fumana S. Galeotta B. Gillis P. Gómez-Alvarez A. Grazian F. Grupp L. Guzzo S. V. H. Haugan W. Holmes F. Hormuth A. Hornstrup S. Ilić K. Jahnke M. Jhabvala B. Joachimi E. Keihänen S. Kermiche A. Kiessling M. Kilbinger B. Kubik M. Kümmel M. Kunz H. Kurki-Suonio S. Ligori P. B. Lilje V. Lindholm I. Lloro G. Mainetti 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 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 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í D. Tavagnacco I. Tereno R. Toledo-Moreo F. Torradeflot I. Tutusaus E. A. Valentijn L. Valenziano T. Vassallo G. Verdoes Kleijn A. Veropalumbo Y. Wang J. Weller G. Zamorani E. Zucca A. Biviano C. Burigana M. Calabrese D. Di Ferdinando J. A. Escartin Vigo G. Fabbian F. Finelli J. Gracia-Carpio S. Matthew N. Mauri A. Pezzotta M. Pöntinen V. Scottez M. Tenti M. Viel M. Wiesmann Y. Akrami V. Allevato S. Anselmi M. Archidiacono F. Atrio-Barandela A. Balaguera-Antolinez M. Ballardini A. Blanchard L. Blot H. Böhringer S. Borgani 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 A. G. Ferrari P. G. Ferreira I. Ferrero A. Finoguenov 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 S. Kruk J. Le Graet L. Legrand J. Lesgourgues T. I. Liaudat A. Loureiro G. Maggio M. Magliocchetti F. Mannucci R. Maoli C. J. A. P. Martins L. Maurin R. B. Metcalf M. Migliaccio M. Miluzio P. Monaco A. Montoro A. Mora C. Moretti G. Morgante S. Nadathur L. Patrizii V. Popa P. Reimberg I. Risso P. -F. Rocci M. Sahlén E. Sarpa A. Schneider M. Sereno A. Silvestri A. Spurio Mancini K. Tanidis C. Tao N. Tessore G. Testera R. Teyssier S. Toft S. Tosi A. Troja M. Tucci C. Valieri J. Valiviita D. Vergani G. Verza P. Vielzeuf N. A. Walton
 



Simulations and nonlinearities beyond ΛCDM. 1. Numerical methods and validation    
   Simulations and nonlinearities beyond ΛCDM – numerical methods and validation






										   

Department of Astrophysics, University of Zurich, Winterthurerstrasse 190, 8057 Zurich, Switzerland
Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth PO1 3FX, UK
Dipartimento di Fisica e Astronomia, Università di Bologna, Via Gobetti 93/2, 40129 Bologna, Italy
INAF-Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Piero Gobetti 93/3, 40129 Bologna, Italy
INFN-Sezione di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy
Max Planck Institute for Gravitational Physics (Albert Einstein Institute), Am Muhlenberg 1, D-14476 Potsdam-Golm, Germany
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
Laboratoire Univers et Théorie, Observatoire de Paris, Université PSL, Université Paris Cité, CNRS, 92190 Meudon, France
Institute of Theoretical Astrophysics, University of Oslo, P.O. Box 1029 Blindern, 0315 Oslo, Norway
Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA, 91109, USA
Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, 85748 Garching, Germany
Department of Physics, Institute for Computational Cosmology, Durham University, South Road, DH1 3LE, UK
Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, University of Manchester, Oxford Road, Manchester M13 9PL, UK
INAF-IASF Milano, Via Alfonso Corti 12, 20133 Milano, Italy
Istituto Nazionale di Fisica Nucleare, Sezione di Bologna, Via Irnerio 46, 40126 Bologna, Italy
Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK
Higgs Centre for Theoretical Physics, School of Physics and Astronomy, The University of Edinburgh, Edinburgh EH9 3FD, UK
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
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
INAF-Osservatorio Astrofisico di Torino, Via Osservatorio 20, 10025 Pino Torinese (TO), Italy
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
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 Astronomico di Roma, Via Frascati 33, 00078 Monteporzio Catone, Italy
INFN-Sezione di Roma, Piazzale Aldo Moro, 2 - c/o Dipartimento di Fisica, Edificio G. Marconi, 00185 Roma, 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
Institute for Theoretical Particle Physics and Cosmology (TTK), RWTH Aachen University, 52056 Aachen, Germany
Institut d'Estudis Espacials de Catalunya (IEEC),  Edifici RDIT, Campus UPC, 08860 Castelldefels, Barcelona, Spain
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
Departamento de Física, Faculdade de Ciências, Universidade de Lisboa, Edifício C8, Campo Grande, PT1749-016 Lisboa, Portugal
Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências, Universidade de Lisboa, Campo Grande, 1749-016 Lisboa, Portugal
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
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
Max Planck Institute for Extraterrestrial Physics, Giessenbachstr. 1, 85748 Garching, Germany
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
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
Universität Bonn, Argelander-Institut für Astronomie, Auf dem Hügel 71, 53121 Bonn, Germany
Dipartimento di Fisica e Astronomia "Augusto Righi" - Alma Mater Studiorum Università di Bologna, via Piero Gobetti 93/2, 40129 Bologna, Italy
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 d'Astrophysique de Paris, UMR 7095, CNRS, and Sorbonne Université, 98 bis 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
Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700 AV Groningen, The Netherlands
INFN-Bologna, Via Irnerio 46, 40126 Bologna, Italy
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
INFN-Sezione di Milano, Via Celoria 16, 20133 Milano, Italy
Departamento de Física Fundamental. Universidad de Salamanca. Plaza de la Merced s/n. 37008 Salamanca, Spain
Departamento de Astrofísica, Universidad de La Laguna, 38206, La Laguna, Tenerife, 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
Dipartimento di Fisica - Sezione di Astronomia, Università di Trieste, Via Tiepolo 11, 34131 Trieste, Italy
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
Ruhr University Bochum, Faculty of Physics and Astronomy, Astronomical Institute (AIRUB), German Centre for Cosmological Lensing (GCCL), 44780 Bochum, Germany
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
Centro de Astrofísica da Universidade do Porto, Rua das Estrelas, 4150-762 Porto, Portugal
Dipartimento di Fisica, Università di Roma Tor Vergata, Via della Ricerca Scientifica 1, Roma, Italy
INFN, Sezione di Roma 2, Via della Ricerca Scientifica 1, Roma, Italy
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
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
INAF-Osservatorio Astronomico di Brera, Via Brera 28, 20122 Milano, Italy, and INFN-Sezione di Genova, Via Dodecaneso 33, 16146, Genova, 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
Mullard Space Science Laboratory, University College London, Holmbury St Mary, Dorking, Surrey RH5 6NT, 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

   



 
  To constrain models beyond ΛCDM, the development of the analysis pipeline requires simulations that capture the nonlinear phenomenology of such models.
  We present an overview of numerical methods and N-body simulation codes developed to study the nonlinear regime of structure formation in alternative dark energy and modified gravity theories. We review a variety of numerical techniques and approximations employed in cosmological N-body simulations to model the complex phenomenology of scenarios beyond ΛCDM. This includes discussions on solving nonlinear field equations, accounting for fifth forces, and implementing screening mechanisms. Furthermore, we conduct a code comparison exercise to assess the reliability and convergence of different simulation codes across a range of models. Our analysis demonstrates a high degree of agreement among the outputs of different simulation codes, providing confidence in current numerical methods for modelling cosmic structure formation beyond ΛCDM. We highlight recent advances made in simulating the nonlinear scales of structure formation, which are essential for leveraging the full scientific potential of the forthcoming observational data from the mission.
   

   

   preparation
    Euclid Collaboration: J. Adamek0000-0002-0723-6740julian.adamek@uzh.ch<ref>
B. Fiorini0000-0002-0092-4321<ref>
M. Baldi0000-0003-4145-1943<ref>,<ref>,<ref>
G. Brando0000-0003-0805-1905<ref>
M.-A. Breton<ref>,<ref>,<ref>
F. Hassani0000-0003-2640-4460<ref>
K. Koyama0000-0001-6727-6915<ref>
A. M. C. Le Brun0000-0002-0936-4594<ref>
G. Rácz0000-0003-3906-5699<ref>
H.-A. Winther0000-0002-6325-2710<ref>
A. Casalino0000-0001-6709-5292<ref>
C. Hernández-Aguayo0000-0001-9921-8832<ref>
B. Li0000-0002-1098-9188<ref>
D. Potter0000-0002-0757-5195<ref>
E. Altamura0000-0001-6973-1897<ref>
C. Carbone0000-0003-0125-3563<ref>
C. Giocoli0000-0002-9590-7961<ref>,<ref>
D. F. Mota0000-0003-3141-142X<ref>
A. Pourtsidou0000-0001-9110-5550<ref>,<ref>
Z. Sakr0000-0002-4823-3757<ref>,<ref>,<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>
P. Battaglia0000-0002-7337-5909<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>
V. F. Cardone<ref>,<ref>
J. Carretero0000-0002-3130-0204<ref>,<ref>
S. Casas0000-0002-4751-5138<ref>
F. J. Castander0000-0001-7316-4573<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>
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>
M. Fumana0000-0001-6787-5950<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>
S. V. H. Haugan0000-0001-9648-7260<ref>
W. Holmes<ref>
F. Hormuth<ref>
A. Hornstrup0000-0002-3363-0936<ref>,<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. Kümmel0000-0003-2791-2117<ref>
M. Kunz0000-0002-3052-7394<ref>
H. Kurki-Suonio0000-0002-4618-3063<ref>,<ref>
S. Ligori0000-0003-4172-4606<ref>
P. B. Lilje0000-0003-4324-7794<ref>
V. Lindholm0000-0003-2317-5471<ref>,<ref>
I. Lloro<ref>
G. Mainetti0000-0003-2384-2377<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>
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>
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>
D. Tavagnacco0000-0001-7475-9894<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>
E. A. Valentijn<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>
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>
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>
V. Allevato0000-0001-7232-5152<ref>
S. Anselmi0000-0002-3579-9583<ref>,<ref>,<ref>
M. Archidiacono0000-0003-4952-9012<ref>,<ref>
F. Atrio-Barandela0000-0002-2130-2513<ref>
A. Balaguera-Antolinez0000-0001-5028-3035<ref>,<ref>
M. Ballardini0000-0003-4481-3559<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. Borgani0000-0001-6151-6439<ref>,<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>
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>
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>
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>
S. Kruk0000-0001-8010-8879<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>
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>
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. Migliaccio<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>
S. Nadathur0000-0001-9070-3102<ref>
L. Patrizii<ref>
V. Popa0000-0002-9118-8330<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>
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>
C. Valieri<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>
N. A. Walton0000-0003-3983-8778<ref>
    Received XXX; accepted ZZZ
=======================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================






§ INTRODUCTION

Significant progress in cosmological observations is expected
in the upcoming years, in particular from the survey <cit.>,  Vera Rubin Observatory’s Legacy Survey of Space and Time <cit.>, the Roman Space Telescope <cit.> and the Dark Energy Spectroscopic Instrument <cit.>. 
These surveys will offer precision observations to high redshifts, allowing
us to study the evolution of the Universe with unprecedented accuracy and potentially uncover the nature of dark matter and dark energy (DE). 
Gaining a deeper understanding of the nature of DE and addressing the long-standing question of whether the cosmological constant (Λ) is responsible for the late-time accelerated expansion of the Universe is indeed one of the primary goals of the Euclid survey <cit.>.

The space telescope was launched on July 1, 2023, and is going to observe billions of galaxies out to redshift z ≈ 2, covering more than a third of the sky in optical and near-infrared wavelengths. will deliver precise measurements of the shapes and redshifts of galaxies <cit.>, from which we will measure weak gravitational lensing <cit.> and galaxy clustering <cit.>. 
These primary probes can be used to rigorously investigate different cosmological scenarios, in particular those related to DE that go beyond the Λ-Cold-Dark-Matter (ΛCDM) concordance model.

Although the ΛCDM model is generally very successful in matching observations, the true identities of CDM and the cosmological constant Λ remain unknown. Additionally, some tensions have persisted in recent years, most notably the Hubble tension <cit.> where local measurements of the Hubble parameter today, H_0, appear to disagree with those inferred from high-redshift observations by around 5 σ. Further examples are the S_8 tension <cit.> and some anomalies found in measurements of the cosmic microwave background <cit.>.
The presence of these tensions
may hint at a breakdown of the ΛCDM model and further motivates the exploration of alternative scenarios.

Over the past few years, cosmologists have explored different possibilities to account for the late-time accelerating expansion of the Universe <cit.> either by introducing a new field, referred to as the DE field or by proposing a modified theory of gravity (MG). A wide range of MG or DE models is equivalent to adding a new light scalar degree of freedom to the theory of General Relativity (GR).

In these theories, the scalar degree of freedom exhibits
time evolution, sometimes accompanied by
spatial fluctuations within the cosmic horizon. Even in the absence of such fluctuations,
the background evolution may be different from ΛCDM,
leading to modifications in structure formation. 
Significant spatial
fluctuations in these models may arise due to various factors, including a low characteristic speed of sound in the theory 
<cit.>, or as a result of the non-minimal coupling of the scalar field to matter or gravity
<cit.>. 
MG and DE theories featuring a coupling of the scalar field to matter
can further affect perturbations at sub-horizon scales by mediating a fifth force. 
If the coupling is universal and includes baryons, a screening mechanism is essential to 
evade the precise constraints of local experiments <cit.>. Screening mechanisms are typically achieved through nonlinear phenomena in such theories.
If, on the other hand, the coupling to matter is non-universal and is confined entirely to the dark sector, local experiments have no constraining power, and cosmological observations provide the main constraints.

Given the diversity of possible DE or MG scenarios, a large information gain is expected from nonlinear scales in the cosmological large-scale structure. These scales must be studied using N-body simulations that capture the essential aspects of the DE or MG models under consideration. This usually means that at least one additional equation needs to be solved for the extra degree of freedom. In many cases, this leads to a difficult nonlinear problem that could require special techniques or approximations that need to be developed. This makes N-body simulations for models of DE or MG a challenging task.

In this paper, we first review the main features of the different classes of DE and MG models that have been proposed over the past years <cit.>. For each of them, we then discuss the
numerical methods implemented within a selection of existing N-body codes
(summarised in Table <ref>). 
Focusing on MG models with a universal coupling, we then compare the results of different N-body implementations for two well-studied theories, namely the Hu–Sawicki f(R) gravity <cit.> and the `normal branch' of the Dvali–Gabadadze–Porrati braneworld model <cit.>. 
We choose simulation parameters following the code comparison paper by <cit.> [W15 hereafter], allowing us to validate a number of new codes against existing results.

This article is part of a series that collectively explores simulations and nonlinearities beyond the ΛCDM model:

    
  * Numerical methods and validation (this work).
    
  * Results from non-standard simulations (Rácz et al. in prep.).
    
  * Cosmological constraints on non-standard cosmologies from simulated Euclid probes (D'Amico et al. in prep.).
    
  * Constraints on f(R) models from the photometric primary probes (Koyama et al. in prep.).

For further details, see our companion papers. The purpose of this first article in the series is to serve as a reference for models beyond ΛCDM and their existing implementations in various codes. This
paper is structured as follows. In Sect. <ref> we give a broad overview of different numerical approaches to treat the additional physics of models beyond ΛCDM. In Sect. <ref> we discuss a number of different codes that implement those approaches and carry out a validation exercise, comparing several recently developed codes with the existing state of the art. We conclude in Sect. <ref>. In an Appendix, we discuss some performance considerations.








§ METHODS 




 §.§ Non-standard background evolution


A wide range of models beyond the simplest cosmological constant scenario are based on an additional scalar degree of freedom – e.g. a classical scalar field ϕ – that evolves dynamically in the expanding Universe and whose background energy density ρ_ϕ provides the source for the observed DE abundance. To induce cosmic acceleration and to match existing constraints on the background expansion history, the equation-of-state parameter w of such an additional field must be sufficiently negative at recent epochs but is poorly constrained at earlier times, which allows for models where w also evolves dynamically as long as it converges to values close to w ≈ -1 in the late Universe. For these models, the DE component modifies the background expansion history of the Universe, which is encoded by the general expression of the Hubble function,

    H^2(a)/H^2_0 = m a^-3 + r a^-4 + k a^-2 + DEe^ -3∫ _1^a1+w(a')/a'da'  ,

where the equation-of-state parameter of DE can be obtained by solving the background field equations – including the evolution of the additional scalar degree of freedom ϕ – or can be parameterised. A common parameterisation suggested by Chevallier, Polarski & Linder <cit.> is based on the desired evolution of w at low redshifts,

    w(a) = w_0 + w_a (1 - a) .

Alternatively, one can set the desired relative abundance of DE at late (DE = 1-m) and early (Ω _ EDE) epochs as in the Early Dark Energy <cit.> parameterisation,

    w(a) = w_0/1 + bln (1/a) ,     
        b = 3w_0/ln1 - Ω _ EDE/Ω _ EDE + ln1 - m/m .

The modified expansion history expressed by Eq. (<ref>) will indirectly affect the evolution of matter density fluctuations and modify the formation process of collapsed structures by changing the Hubble friction term in the equation for linear matter perturbations, which in Newtonian gauge and in Fourier space for sub-horizon scales reads:

    δ̈_m + 2Hδ̇_m = 4π G ( mm + DEDE),

where m and DE are the density contrasts of matter and DE perturbations, respectively, G is  Newton's constant, and a dot represents a derivative with respect to cosmic time.

Besides the richer background dynamics that is endowed by an evolving field, whenever DE is promoted from a cosmological constant to a dynamical degree of freedom, the model also acquires an additional layer of complexity: the presence and evolution of DE fluctuations around the mean-field configuration. This corresponds to the situation where δ_DE in Eq. (<ref>) is non-negligible, whereas in ΛCDM it would vanish identically
at all scales. Like any other density perturbations, inhomogeneities in the DE would then contribute to the peculiar gravitational potential that governs
the evolution of matter perturbations and thus the
formation of cosmic structures as shown by Eq. (<ref>).

However, in many of the simplest scalar-field scenarios, such perturbations are negligible at sub-horizon scales because the speed of sound c_s of the scalar field is naturally close to the speed of light. Ignoring them for the purpose of numerical simulations, the only modification of N-body algorithms required to simulate these DE models is given by an appropriate calculation of the cosmic expansion rate. The most common approach amounts to tabulating the specific expansion rate of the universe for the model to be simulated according to Eq. (<ref>) and replacing the standard analytical calculation of the Hubble function within the N-body algorithm with an interpolated value from the tabulated solution that is provided to the code as an input.
This approach is implemented by most of the simulation codes employed within the Euclid Collaboration to perform cosmological simulations in homogeneous DE models beyond ΛCDM.



 §.§ Linearised DE perturbations


Although a wide range of DE models are characterised by negligible DE fluctuations as discussed above, some specific scenarios may not fulfil such a condition at all scales and/or at all times, either because they feature a lower value of the DE speed of sound, allowing DE perturbations to grow on scales above the associated Jeans length that then falls inside the cosmological horizon, or because additional interactions – besides gravity – can induce the growth of such perturbations. The former case corresponds to the class of clustering DE models, while the latter is known as coupled DE.



  §.§.§ Clustering DE
 
The clustering DE models are characterised by two time-dependent variables: the speed of sound c_s and the equation of state parameter w. In these theories, the DE component clusters on scales larger than the associated
sound-horizon/Jeans scale, λ_s = H/c_s and DE perturbations decay quickly below the sound-horizon scale. For a sufficiently small speed of sound, we may even
expect nonlinear DE structures to form. At a fundamental level, clustering DE models are analogous to the k-essence type of theories, so that the action reads <cit.>

    S = ∫d^4x √(-g)[ c^4 R/ 16π G  + P(X, ϕ)  + ℒ_m],

where P is a general function of the kinetic term X ≡ -1/2∇_μϕ∇^μϕ and the scalar field ϕ, and ℒ_m is the matter Lagrangian. For a given P(X, ϕ), the speed of sound and the equation of state are given by <cit.>

    w=P/P-2X P_,X ,      c_s^2 = P_,X/2XP_,XX+P_,X ,

where the subscript “,X” denotes the partial derivative with respect to X. We therefore need to specify the function P(X, ϕ) to derive the equations of motion for the k-essence scalar field. However, since there are
many possible choices, we can instead employ the effective field theory (EFT) approach to model the dynamics of k-essence DE. The EFT framework, although
not a fundamental theory, offers several advantages <cit.>, such as being a description of a wide range of theories within some scales. The EFT is a perturbative approach based on the assumption that the scalar field perturbations remain small over the scales of interest. It is worth noting that the regime of nonlinear matter clustering is accessible to the EFT framework as long as the scalar field perturbations remain small.
The k-essence theories or clustering DE models are implemented in several N-body and Einstein–Boltzmann codes. In  <cit.> and  <cit.>, these theories are implemented using the fluid picture.
In  <cit.>, the EFT equations are implemented and can be controlled using the EFT parameter α_K≡ 3(1+w)c_s^-2 within the code. On the other hand, in  <cit.>, which is an N-body code based on  <cit.>, nonlinear equations for clustering DE are implemented as an independent component, and the k-essence field for small c_s can form nonlinear structures. In some N-body codes, for example in  <cit.>, clustering DE is implemented through a linear solution from an Einstein–Boltzmann solver. This is a good assumption for large speeds of sound, but for small ones, this method does not allow for the response of DE to the nonlinear matter structures.



  §.§.§ Coupled quintessence


Moving to the case of coupled DE models, the interaction can be formulated at a fundamental level by introducing a direct coupling between the scalar field and the spatial curvature R in the so-called Jordan frame <cit.>, so that the action reads

    S = ∫d^4x √(-g)[c^4 f(ϕ, R)/16π G-1/2Z(ϕ)∇^μϕ∇_μϕ - V(ϕ) + ℒ_m] ,

where f(ϕ, R) is a function that couples the scalar field to the curvature, Z(ϕ) is a function that allows for non-standard kinetic terms, V(ϕ ) is the scalar field self-interaction potential, and the matter Lagrangian contains at least one cold species characterised by some rest mass m_0.

Alternatively, the interaction can be formulated by including source terms in the covariant conservation equations of the interacting species in the so-called Einstein frame,

    ∇^μ T^(c)_μν   =    -β_(c)(ϕ)/M_ Pl T^(c)∇_νϕ ,
    ∇^μ T^(b)_μν   =    -β_(b)(ϕ)/M_ PlT^(b)∇_νϕ , 
    ∇^μ T^(ϕ)_μν   =   1/M_ Pl[β_(c)(ϕ)T^(c) + β_(b)(ϕ)T^(b)]∇_νϕ ,

where T_μν^(Y) is the stress-energy tensor of a given species Y, T^(Y) is its trace, β_(Y)(ϕ) is the coupling function of species Y, the labels c, b, ϕ refer to the dark matter, baryon, and scalar field species, respectively, and M_ Pl≡ (ħ c)^1/2 (8π G)^-1/2 is the reduced Planck mass.

While in the former case the interaction will be universal (i.e. involving all matter species with the same strength), which goes under the name of Extended Quintessence,  the latter approach allows for non-universal couplings that may selectively involve individual species, for example by separately choosing the coupling functions for baryons and dark matter. 

In the case of a universal coupling (that is, if β_(b) = β_(c)), the two approaches can be related to one another through a Weyl transformation of the metric <cit.>, and are therefore equivalent.
On the other hand, the possibility to leave the baryonic component of the Universe only minimally coupled evades Solar System constraints <cit.> on the
deviations from standard gravity thereby avoiding the need for screening mechanisms. This is the case of Coupled Quintessence models <cit.>, where the direct coupling between the scalar field and massive (non-baryonic) particles can support stable perturbations of the DE field at sub-horizon scales <cit.>. In general, such perturbations may even become nonlinear in the presence of a sufficiently strong coupling <cit.>. Nonetheless, a large class of widely studied coupled DE models is known to feature scalar perturbations of the order of the standard Newtonian gravitational potential <cit.>, thereby remaining in the linear regime at all times and scales of cosmological interest. This allows us to linearise the corresponding field equations and derive modified equations of motions for massive particles, including the contribution of the additional force arising from the direct coupling with the scalar field <cit.>. 

In fact, a general feature of coupled DE models is the existence of a `fifth force' mediated by the
scalar field. The new force can be expressed as an additional acceleration
experienced by a massive coupled particle, which in comoving coordinates will be given by

    a⃗_Y,5th = -β_(Y)(ϕ )∇⃗δϕ ,

where Y identifies a coupled matter species, and δϕ is the scalar field fluctuation.
This extra acceleration term is added to the standard Newtonian acceleration acting on all massive particles,

    a⃗_N= -Hv⃗ -∇⃗Φ_N ,

where v⃗ is the peculiar particle velocity in comoving coordinates, and Φ_N is the peculiar Newtonian potential obeying the standard Poisson equation 

    ∇ ^2Φ_N = 4π G a^2 ∑_Yρ_Yδ_Y ,

where the sum runs over all clustering species in the Universe.

Therefore, solving for the dynamical evolution of massive coupled particles requires solving for the scalar field perturbation δϕ entering in Eq. (<ref>), which in the most general case follows a 
nonlinear elliptic equation,

    ∇ ^2δϕ = F(δϕ ) + ∑_Y 8π Ga^2 β_(Y)(ϕ)δ_Y ,

with F a function of the scalar field fluctuation δϕ, and where the sum runs over all the coupled matter species with their respective couplings β _(Y)(ϕ ). 

For the particular case of a coupled DE model with a non-universal interaction <cit.> involving only dark matter and leaving baryons uncoupled (i.e. β_(b) = 0) the function F(δϕ ) in Eq. (<ref>) is negligible compared to the term associated with matter density perturbations <cit.> and can be safely discarded. As a result, the scalar-field equation reduces to

    ∇ ^2δϕ≈ 8π Ga^2 β_(c)(ϕ)δ_c = 2β_(c)(ϕ)∇ ^2Φ_c ,

where Φ_c is the Newtonian potential generated by the distribution of the coupled dark matter particles, that is

    ∇ ^2Φ_c = 4π G a^2  ρ_cδ_c .

Therefore, the solution for the scalar field perturbations will be directly proportional to the potential Φ_c according to the relation

    δϕ≈ 2β_(c)Φ_c .

The acceleration equation for a coupled particle can then be rewritten as 

    a⃗_c = -Hv⃗_c - ∇⃗Φ_b - ∇⃗(1 + 2β_(c)^2(ϕ))Φ_c ,

assuming here for simplicity that other clustering species (such as massive neutrinos) give a negligible contribution to the total Newtonian potential such that Φ_N = Φ_c + Φ_b.
This modified acceleration equation introduces a further modification to be implemented in N-body simulation codes for Coupled Quintessence cosmologies besides the specific expansion history of each particular model. This often requires substantial modifications in the gravity solvers of conventional N-body codes, as the algorithms need to evolve coupled and uncoupled massive particles (typically dark matter and baryons, respectively) with different equations and should therefore treat these components separately. Even under the approximation of a purely collisionless treatment (i.e., ignoring the hydrodynamical and astrophysical processes that affect standard baryonic matter leading to the formation of stars and galaxies) that is often employed for large-volume simulations targeted at galaxy surveys such as Euclid, both coupled and uncoupled matter species must be included in the simulation to provide a consistent representation of the dynamics at all scales: as baryons and dark matter evolve differently, assuming that all matter is dark would lead to an overestimation of the effects of the coupling and, thus, biased results. This approach is implemented in the  code <cit.> which is employed for
Coupled DE simulations performed within the Euclid Collaboration.

Distinguishing between coupled and uncoupled particle types in simulations of Coupled Quintessence is also crucial for proper treatment of two other effects that characterise these cosmological models beyond the fifth force described by Eq. (<ref>).
The first is the mass variation of coupled particles due to the exchange of rest-frame energy with the DE scalar field, which arises as a direct consequence of the modified continuity equations (<ref>) and of the assumption of particle number conservation. More specifically, the mass of coupled particles evolves as a result of the evolution of the background scalar field according to

    m_Y(a) =  m_0exp(- ∫_ϕ _0^ϕ (a)β_(Y) (ϕ )dϕ/M_ Pl) .

Such a mass variation, which involves only particle species with a non-vanishing coupling to the scalar field, must be taken into account in N-body algorithms by changing the mass of individual simulation particles at every time step. This is normally done by tabulating the mass as a function of scale factor a by numerically integrating Eq. (<ref>) along with the background dynamics of the scalar field ϕ, and interpolating from that table as the simulation progresses.

The second effect is an additional force (on top of the fifth force) acting on coupled particles as a consequence of momentum conservation due to the particles'  mass variation described by Eq. (<ref>), which takes the form of a velocity-dependent extra acceleration behaving either as a friction or as a drag, depending on the relative signs of the coupling function β (ϕ ) and of the scalar field velocity ϕ̇ <cit.>,

    a⃗_Y,v = -β _(Y)(ϕ)/M_ Plϕ̇v⃗ .

Such a velocity-dependent acceleration is responsible for a very rich phenomenology characterising Coupled Quintessence models, especially on highly nonlinear scales <cit.>, and must be included in N-body simulations as well for a fully consistent treatment of these scenarios. This is done by adding the extra acceleration described in Eq. (<ref>) to the total acceleration (i.e. Newtonian plus fifth force) of all coupled particles in each time step,

    a⃗_Y = a⃗_N + a⃗_Y,5th + a⃗_Y,v .

The relevant quantities β (ϕ) and ϕ̇ can again be interpolated from a table obtained by integrating the background dynamics of the system. This is the approach implemented in the  code that has been used to run Coupled Quintessence simulations within the Euclid Collaboration.



  §.§.§ Momentum exchange and dark scattering


A further example of interacting DE cosmologies characterised by scalar-field perturbations that always remain linear is given by models of pure momentum exchange <cit.> between the DE field and massive particles like dark matter or baryons. A limiting case is given by the Dark Scattering scenario <cit.> where the momentum transfer between the two components is modelled as the elastic scattering of massive particles moving through a homogeneous DE fluid with equation of state w. This results in an extra force acting on the moving massive particles which is proportional to their comoving velocity, similar to the velocity-dependent force described by Eq. (<ref>) for Coupled Quintessence models. However, the origin of this force is completely different in this case, as it does not originate from the mass variation of particles but rather from the momentum transfer with the DE field. As a result, the scattering acceleration can be expressed as

    a⃗_s = -(1+w)3H^2Ω_DE/8π Gξv⃗ ,

where the parameter ξ is defined as

    ξ≡σ/m ,

with σ denoting the scattering cross section and m the typical mass of the scattering particle species.

This type of interaction can be implemented in N-body algorithms <cit.> in a very similar way as the velocity-dependent acceleration in Coupled Quintessence scenarios, as the factors entering Eq. (<ref>) are all either constants or background quantities that can be interpolated at every timestep from tabulated data. This is the approach implemented in the  code that has been used to run the  and  simulations <cit.>.

Although Dark Scattering represents a limiting case of the more general class of pure momentum-exchange models between matter and DE <cit.>, for which further modifications to the standard particle dynamics are expected besides the drag force of Eq. (<ref>), recent works <cit.> have shown that such additional modifications are generally subleading with respect to the drag force so that their effect on structure formation can be neglected. This ensures that the current implementation of Dark Scattering within the simulations used in the Euclid Collaboration can be considered representative of the general class of momentum-exchange cosmologies.





 §.§ Nonlinear scalar field perturbations



In models where a scalar field couples to matter universally, or at least to baryons in a relevant way, some mechanism to suppress the coupling is required to satisfy the stringent local tests of gravity. This is commonly referred to as `screening'. Screening mechanisms are achieved by nonlinearity in the scalar-field equation coupled to matter. 
The equation determining the evolution of the scalar field is typically a wave equation of the form 

    □ϕ = S(ϕ, ∇_μϕ, ∇_μ∇_νϕ, m) .

Here, □≡∇_μ∇^μ represents the d'Alembertian operator and S is a nonlinear function that depends on the matter density, the scalar field, and its derivatives. 
Various methods have been developed to solve this nonlinear scalar field equation in N-body simulations, where the nonlinear density m is modelled by collisionless particles .

Several approximations are often used to solve these nonlinear equations. The most common one is the quasi-static approximation. The scalar field can be split into a background part, ϕ̅, and a perturbation, δϕ, as
ϕ = ϕ̅ + δϕ.
The quasi-static approximation amounts to ignoring the time dependence of the scalar field perturbation, i.e. assuming ϕ̇≃ϕ̇̅̇.
The partial differential equation (PDE) of the field perturbation, which in its original form may have been of the hyperbolic or parabolic type, is therefore cast into an elliptic form so that the scalar field solution at any given time depends solely on the matter configuration at that time. This is a good approximation whenever the speed of sound of the scalar field is small
<cit.>, which is the case for the MG models considered here.
Non-quasistatic cosmological simulations have been conducted for several MG models using different techniques, such as the explicit leap-frog method and the implicit Newton–Gauss–Seidel method <cit.>.

The scalar-field solution is required for the computation of the total gravitational potential Φ that acts on the matter particles,

    ∇^2 Φ = 4 π G a^2 (δm + δeff (ϕ) ) ,

where the effective density depends on the scalar field. There are two common ways of solving for this total gravitational force. The first option is to solve first for ϕ and then use this solution to compute the source term in Eq. (<ref>) and solve for Φ using a standard Poisson solver to get the total force ∇Φ. The other option is to 
apply the total force ∇Φ_N + ∇ϕ to the particles.

Under the quasi-static approximation, the scalar-field equation assumes the same form as the usual Poisson equation, 

    ∇^2 ϕ = S(ϕ, ∇_i ϕ, ∇_i ∇_j ϕ, δm) .

The main difference is that the scalar-field equation is generally nonlinear. This nonlinear behaviour implies that conventional techniques, such as using Fourier analysis, cannot be used to solve the equation. Numerous approaches have been developed to address this challenge, and we refer the reader to  for details. For computational methods that aim to accurately solve a nonlinear equation on refined grids, the approach typically involves discretising the equation in a suitable way and employing an iterative algorithm, such as the Newton–Raphson method, to successively refine solutions based on an initial guess. To speed up convergence, many of these methods incorporate so-called `multigrid' acceleration techniques which we quickly review here.



  §.§.§ Nonlinear multigrid algorithm



A generic way to solve nonlinear elliptic PDEs is to couple the multigrid algorithm to the Newton–Raphson method,

    u^ new = u^ old - ℒ(u^ old)/∂ℒ/∂ u^ old ,

where u is the discretised field, ℒ is the differential operator (which is a Laplacian for Newtonian gravity) and the superscripts refer to the new or old estimate of the solution in one Newton–Raphson iteration. The Newton–Raphson method produces linear equations for the correction terms, which are solved by the Full-Approximation-Storage Multigrid algorithm <cit.>. For a review of these methods applied to MG simulations, see e.g. <cit.>, <cit.> and .

A simple sketch of the algorithm goes as follows. One starts with a guess for the solution on a grid; this could be anything from a constant value across the grid to using the solution from the previous timestep in the simulation. One then loops a few times over all cells in the grid, updating the solution using Eq. (<ref>). This solution is then restricted to a grid with half the resolution, the solution is updated again, and this process is repeated recursively up to the coarsest grid (one with only 2^3 cells). The solution is then interpolated to the finer grid, updated once more, and this is done recursively until one reaches the finest grid we started with. One such cycle is called a V-cycle, and one repeats such V-cycles until convergence is achieved. The advantage of having this stack of coarser grids is that it helps to accelerate the convergence of the largest modes in the solution.



  §.§.§ Screening with nonlinearity in potentials
 
In models where screening is achieved by nonlinearity in a potential or coupling function, the equation for the scalar field becomes 

    ∇^2 ϕ = 4 π G a^2 β(ϕ) δm + V(ϕ) ,

where β(ϕ) and V(ϕ) are nonlinear functions of ϕ. A typical example is f(R) gravity. In this class of models, the value of the scalar field changes by orders of magnitude. To enhance numerical stability, a common technique involves redefining the scalar field in terms of a new variable. The redefinition to choose depends on the specific model under consideration. It is typically chosen to prevent the occurrence of unphysical values of the scalar field during Newton–Raphson iterations. For example, for f(R) models the scalar field ϕ = f_,R will be driven towards zero in high-density regions, but at the same time f_,R cannot cross zero, as the potential becomes singular in this scenario. To avoid this issue, a commonly used field redefinition is u ≡ln[f_,R/f̅_,R(a)] <cit.>. However, this transformation introduces additional nonlinearity and in some models, such as Hu–Sawicki f(R) models, this transformation is not necessary and might even lead to considerable performance losses in a simulation.

In some cases, this can be avoided. For example, <cit.> noticed that for Hu–Sawicki f(R) gravity with n = 1 <cit.>, when making the change of variable u = √(-f_,R), the field equation could be recast as a depressed cubic equation,

    u^3 + pu + q = 0 ,

which possesses analytical solutions <cit.>. Although the Gauss-Seidel smoothing procedure is still needed (because p depends on the values of the field u in neighbouring cells), this removes the Newton–Raphson part and expensive exponential/logarithmic operations from the method of <cit.>, therefore leading to significant performance gains. <cit.> also generalised this improved relaxation approach to the cases of n=0 (which strictly speaking is not a variation of the Hu–Sawicki model) and n=2.

For other models like the symmetron, which has a Higgs-like potential, the scalar field is free to cross zero, and no field redefinition is needed (apart from a simple rescaling).



  §.§.§ Screening with nonlinearity in kinetic terms

In another class of models, nonlinearity emerges within the kinetic term.
For example, in models with the Vainshtein screening mechanism <cit.>, the equation exhibits nonlinearity in the second derivatives of the scalar field, 

    ∇^2 ϕ = 4 π G a^2β(a) δm + g(∇_i ϕ, ∇_i ∇_j ϕ) ,

where β(a) is a time-dependent coupling function. The simplest example here is the DGP model which was first simulated by <cit.>. In such cases, the operator-splitting trick <cit.> can be employed.  This approach can simplify the equations, avoiding potential issues associated with imaginary square roots, and improving code performance. This trick is particularly useful for the DGP braneworld models and other Vainshtein screening models, such as Cubic and Quartic Galileons .



  §.§.§ Approximate treatments of screening


Some models allow linearisation of the nonlinear equation using some approximation.
One approach ( ;  ; see also the appendix of  ) is to introduce the screening factor for the matter density perturbation 

    ∇^2 ϕ = c^2m^2a^2/ħ^2ϕ  + 4π Ga^2 δmϵ_ screen(Φ_N, |∇Φ_N|, ∇^2 Φ_N) ,

where the screening function depends on the Newtonian potential Φ_N. This type of parameterised modified gravity is referred to as `type 1' in Table <ref>.
One specific method, developed in <cit.>, starts from linearising Klein–Gordon's equation. In this formalism, one solves the Poisson equation in Fourier space,

    - k^2 Φ= 4 π G_ eff(a,k) a^2 δm,,

where the function G_ eff(a,k) approximates the effective Newton's constant introduced by the screening effect of the scalar field on small scales. This function is given by 

    G_eff(a,k) = G + Δ G_eff(a,k) = G + (G_eff^lin(a) - G)M(a,k) ,

with G_ eff^ lin(a) being the asymptotic linear effective Newton's constant that depends only on time, and M(a,k) is a function that approximately captures the nonlinear corrections introduced by the scalar field on small scales. This function allows Eq. (<ref>) to transition from G_ eff(a,k)→ G_ eff^ lin(a) on large scales to G_ eff(a,k)→ G on small scales. This type of parameterised modified gravity is referred to as `type 2' in Table <ref>. A procedure to fix M(a,k) is described in <cit.>, which has the advantage of avoiding additional parameters to tune the screening efficiency.

One can also choose to parameterise the nonlinear contribution using an effective Newton's constant at both small and large scales. If the modifications of gravity are encoded in a scale-dependent function Δ G_eff(a,k) as in Eq. (<ref>),
then we can propose a similar equation in real space,

    G̃_eff(a, r) = G +ΔG̃_eff(a, r) ,

where ΔG̃_eff(a,r) is the Fourier transform of Δ G_eff(a,k). In practice, an additional approximation is made, namely ΔG̃_eff(a,r) ≈Δ G_eff(a,k → 1/r).








This approach allows the encoding of nonlinear contributions over the whole range of scales modelled by N-body algorithms through real-space equations, for instance, the Tree Particle-Mesh (TreePM) method implemented in codes like .
Provided the parameterisation is effective, this is expected to increase the accuracy of the estimation of the nonlinear effects.

Several parameterisations have been proposed for this kind of approach, either with additional tuning parameters, such as in <cit.>, or based on local small-scale environmental properties to avoid the need for any extra parameters, such as in <cit.>.



§ ADDITIONAL CODE VALIDATION


A code comparison for simulations that solve nonlinear scalar field perturbations is presented in . Since then, various new codes have been developed.  In this section, we show a comparison of the predictions for the power spectrum and the halo mass function using the simulations of  as a reference and starting from the same initial conditions. These were generated using second-order Lagrangian perturbation theory in a ΛCDM cosmology with m = 0.269, Ω_Λ = 0.731, 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. As in , we compare simulations for f(R) and nDGP models. In these models <cit.>, the background expansion history is closely approximated by that of ΛCDM.
Furthermore, the effect of modified gravity can be ignored at z = 49, thus it is justified to use the initial conditions of ΛCDM. The measurements of the power spectrum and mass function are performed by the pipeline developed in the second article of this series, Rácz et al. (in prep.). Based on two models only, our comparison does not encompass the full diversity of numerical methods discussed in the previous section. In many cases some validation of the various implementations can be found in the corresponding references.



 §.§ Summary of codes used in the validation


Table <ref> shows an overview of the simulation codes considered in this section and provides a quick reference of their capabilities and limitations. For each of them, a short summary is presented here. In the Appendix we comment on the trade-off between accuracy and computational cost of the implementations.



  §.§.§ 


First presented in <cit.>, this code is based on the moving-mesh N-body and hydrodynamical simulation code  <cit.>, which uses a TreePM algorithm to calculate gravitational forces. The additional modified gravity force (fifth force) is calculated with a relaxation solver <cit.> that is accelerated by the multigrid method and uses adaptive mesh refinement (AMR). It currently also supports simulations for the nDGP model <cit.>, as well as massive neutrinos implemented using the δ f method <cit.>.

To solve the modified gravity equations, the density field is projected onto the AMR grid, constructed in such a way that each cell on the highest refinement level contains at most one particle (except if a pre-set maximum refinement level is reached; the cell size at this level is of the order of the smoothing length of the standard gravity solver). Once the field equation is solved to obtain the scalar field configuration, the modified gravity force can be computed from its gradient using finite differencing. 

Since  computes the forces on all particles simultaneously, and the modified gravity field equations are generally highly nonlinear (with a poor convergence rate of the relaxation algorithm), this is computationally expensive compared to 's Newtonian gravity solver. However, the maximum acceleration of the modified gravity force is smaller than that of Newtonian gravity, mainly because the latter occurs in regions with high density where screening occurs. This allows the modified gravity solver to run using larger time steps <cit.>, resulting in significantly reduced computational cost. Together with 's efficient MPI parallelisation and lean memory footprint, these have made it possible to run the large number of f(R) simulations used in various recent works, such as the - <cit.> simulation suite of 200 f(R) and nDGP models. This has allowed accurate emulators of various physical quantities or observables to be constructed.

Another highlight of  is its capability to run realistic galaxy formation simulations in a cosmological box <cit.>, thanks to the use of the Illustris-TNG subgrid physics model <cit.>. More recently, it has been used for larger-box hydrodynamical simulations with a realistic recalibrated Illustris-TNG model <cit.>, enabling the study of galaxy clusters in modified gravity.

The  simulations used in this paper were run using a residual criterion of ϵ = 10^-2 and a maximum refinement level () of 10 for the nDGP simulations and 18 for the f(R) models with a gravitational softening of 0.01 h^-1 Mpc.



  §.§.§ 


 <cit.> is a modified version of the TreePM N-body code  <cit.> implementing an AMR solver for the scalar degree of freedom f_,R characterising the widely-studied Hu–Sawicki f(R) gravity model.
In , the same tree structure that is employed to solve for standard Newtonian gravity is also used as an adaptive grid to solve for the scalar field configuration through an iterative Newton–Gauss–Seidel (NGS) relaxation scheme (see Sect. <ref>) with the Full-Approximation-Storage Multigrid method (see Sect. <ref>) and with the field redefinition u ≡ln[f_,R/f̅_,R(a)] (see Sect. <ref>).  also allows MG simulations to be run with massive neutrinos <cit.> using the neutrino particle method <cit.>.

For the simulations presented here, the relative tree opening criterion was used with an acceleration relative error threshold of 0.0025, and a uniform grid with 512^3 cells was employed to compute long-range Newtonian forces. Concerning the MG field solver, a residual tolerance of ϵ = 10^-2 was set for the V-cycle iteration and a maximum refinement level of 18 was used for the AMR grid, corresponding to a spatial resolution of 1 h^-1 kpc at the finest grid level, compared to a gravitational softening of 18 h^-1 kpc, following the setup adopted in . 



  §.§.§ 


 <cit.> is a generic modified gravity simulation code based on the publicly-available N-body and hydrodynamical simulation code  <cit.>. Originally developed for f(R) gravity, this code takes advantage of the adaptive mesh refinement of  to achieve the high resolution needed to solve the scalar field and hence the fifth force in high-density regions. The nonlinear f(R) field equation is solved with the standard Gauss-Seidel approach as first applied by <cit.>, but it was later replaced by the more efficient algorithm of <cit.>. The code has since been extended for simulations for the generalised chameleon <cit.>, symmetron and dilaton <cit.>, nDGP <cit.>, cubic Galileon <cit.>, quartic Galileon <cit.>, vector Galileon <cit.> and nonlocal gravity <cit.>.

For the  simulations used in this paper, we have used a domain grid (the uniform mesh that covers the whole simulation domain) with 2^9=512 cells per dimension, and the cells are hierarchically refined if they contain 8 or more effective[Since the code uses a cloud-in-cell mass assignment scheme, it is tricky to count particles in cells because a particle can contribute to the densities of 8 nearby cells. Here “effective” is used to mean that the density in a given cell, multiplied by the volume of the cell, is equivalent to the total mass of a specified number of particles.] particles. The highest refined levels have effectively 2^16 cells, leading to a force resolution of about 0.0075 h^-1 Mpc. 



  §.§.§ ISIS


The  code <cit.>, like the  code above, is based on  <cit.>. It contains a scalar field solver that can be used to simulate generic MG models with nonlinear equation of motion and has been used to simulate models such as f(R) gravity, the symmetron model, nDGP and disformal coupled models <cit.>. It also allows for hydrodynamical simulations that have been used to study the interplay between baryonic physics and modified gravity <cit.>.  
Furthermore, the code has the capability to go beyond the quasistatic limit and study the full time dependence of the scalar field <cit.>.
The scalar field solver used in the code is a Gauss–Seidel relaxation method with multigrid acceleration, very similar to the one in  described above.

For the  simulations presented in this paper, we have used a domain grid (the uniform mesh that covers the whole simulation domain) with 2^9 = 512 cells per dimension, and the cells are hierarchically refined if they contain 8 or more effective particles.



  §.§.§ 


[https://github.com/mianbreton/pyscogithub https://github.com/mianbreton/pysco] is a particle-mesh (PM) code written in  and accelerated with  which currently supports Newtonian and f(R) gravity <cit.>. While multiple flavours of solvers based on Fast Fourier Transforms (FFT) are available, for the present paper, we use a multigrid solver to propose something different from other codes in this comparison project (other codes that also use multigrid are AMR-based). uses a triangular-shaped cloud mass assignment scheme and solves the linear Poisson equation using multigrid V-cycles with a tolerance threshold of the residual of 10^-3, and two F-cycles <cit.> to solve the additional field in f(R) gravity with the nonlinear multigrid method described in Sect. <ref> and Eq. (<ref>). Our convergence threshold is very conservative since we do not intend to conduct a convergence study in this paper (a less conservative threshold could still give reasonable results at much lower computational cost). Furthermore, to resolve the small scales we use a coarse grid with 2048^3 cells, resulting in roughly 500 time steps to complete the simulations.



  §.§.§ 


 simulations were performed using the Fourier-Multigrid Library .[https://github.com/HAWinther/FMLgithub https://github.com/HAWinther/FML] These simulations are based on the COmoving Lagrangian Acceleration (COLA) method <cit.>, which combines Lagrangian perturbation theory with the PM method to reduce the number of time steps that are required to recover clustering on large scales. The  N-body solver in the  library contains implementations of various DE and MG models like the DGP model, the symmetron model, f(R) gravity and the Jordan–Brans–Dicke model <cit.>.
For the PM part, we used N_mesh^1/3 = 5 N_p^1/3, i.e. a mesh discretisation five times smaller than the mean particle separation, and the total of 150 time-steps linearly spaced along the scale factor to achieve a good agreement of the mass function with AMR simulations. We also ran low-resolution simulations with N_mesh^1/3 = 3 N_p^1/3 with 100 time steps, and 
checked that the nonlinear enhancement of the power spectrum from these low-resolution simulations agrees well with the one from the high-resolution ones. Screening is included using approximate treatments described in Sect. <ref>. 
For f(R) gravity models, the  code has one parameter to tune the strength of chameleon screening called . For the model with f̅_,R = 10^-6, we used the default value  while we used  for f̅_,R = 10^-5. For nDGP, we used Gaussian smoothing with a smoothing radius of 1 h^-1 Mpc to compute the density field for screening and did not use an option to enforce the linear solution at small wavenumbers k. We also ran simulations based on the screening approximation using G_eff(a, k) given by Eq. (<ref>) for nDGP, and in these simulations, we have used the same  settings as in the other screening approximation implementation. The biggest advantage of this screening implementation is that it does not require any additional tuning parameter related to screening, that is, the function G_eff(a,k) is completely defined by the theoretical model one wants to simulate.



  §.§.§ 


 is an extension of the TreePM code . It introduces modifications at large scales in the PM part with an effective Newton's constant according to Eq. (<ref>), while the force at small scales is modified in the tree part with Eq. (<ref>). 
In this respect,  implements an approximate solver for the extra force induced by different possible MG theories. On the other hand, differently from other approximate methods, the dynamics of matter particles under the effect of the (dominant) standard gravity force is treated with the full TreePM solver of  of which it retains the accuracy in modeling the nonlinear density field.
For nDGP, the functional form of G̃_eff(a,r) is described in <cit.>, while for f(R) it is introduced as in <cit.>. In the latter case, no additional parameters are needed to describe the screening. For the first nDGP case instead, in addition to the theoretically defined parameters N_0=1, B = 1/[3β(a)], b=2 and a=3, we consider the screening scale k_th = 0.4 h Mpc^-1 defined by the parameterisation

    r_V≅2/3(3 Ω_m H^2 r_c^2/β^2)^1/3 r_th .

where r_ V is the Vainshtein radius of the model for a spherically symmetric density perturbation <cit.>.





 §.§ Comparison of the power spectrum




Before discussing the effect of extra degrees of freedom, we compare the different codes within a ΛCDM cosmology. For reference and only for the case of ΛCDM, we include results from the code  that was also used for the Flagship simulations in <cit.>. Figure <ref> shows the matter power spectra at two different redshifts, z=1 (left panels) and z=0.667 (right panels). The lower panels show the relative difference to the simulation carried out with  that we use as a reference throughout. Since  and  have both been developed from the  code, the agreement between these two codes is better than 1%. These AMR codes tend to underestimate the power on small scales when compared to tree-based codes, mainly due to the mesh refinement criterion used in our simulations. Fine-tuning some precision parameters, we expect that a better agreement can be achieved, cf. <cit.>. As we can see, the codes that perform closely to our benchmark at small scales are the ones that also share the same tree-PM gravity solver,  and . , which uses a tree structure and the fast multipole method to compute gravitational forces, is in excellent agreement with  up to k ∼ 1 h Mpc^-1, yielding slightly more power at higher k.

Also shown in the figure are the results from , a pure PM code that uses a fixed grid to solve the Poisson equation. The fixed resolution of the PM grid explains why  consistently suffers from low force resolution at k ≳ 1 h Mpc^-1. 
All our simulations use the exact same initial particle data, such that our comparisons should not be contaminated by cosmic variance. However, for  simulations, we also reconstructed the initial density field from the initial particle distribution to estimate the displacement fields required for the 2LPT calculations in those simulations.
While incorporates a full N-body solver, its small-scale accuracy is constrained by resolution limitations that arise from the absence of AMR. Moreover, within the scope of this code comparison project, employs a multigrid algorithm, resulting in reduced small-scale clustering compared to FFT (on a regular grid), thereby explaining the observed deficiency in power at wavenumbers k ≳ 1 h Mpc^-1.

For modified gravity models, we compare the ratio between the power spectrum computed for that model and the one found in the ΛCDM reference cosmology. Taking such ratios largely cancels out the differences seen in ΛCDM between the codes, making the differences due to the treatments of the extra scalar field more apparent. Following , we consider two MG models that exhibit two different types of growth dependence and screening mechanisms. The first one is the so-called Hu–Sawicki f(R) which exhibits a scale-dependent growth factor at linear order, and a screening mechanism realised by nonlinearities in the potential, see Sect. <ref>. The strength of the modification of gravity is characterised by the background value of the scalar degree of freedom, and we study the cases f̅_,R = 10^-5 and f̅_,R = 10^-6, labelled `F5' and `F6', respectively, the latter being closer to GR. The second one is nDGP, a braneworld theory defined in a five-dimensional spacetime. The growth factor in this theory is scale independent, and screening is realised through nonlinearities in the kinetic terms, as discussed in Sect. <ref>. Here, the strength of the modification of gravity is characterised by the value of the cross-over scale, and we study the cases r_c = 1.2 H_0^-1 and r_c = 5.6 H_0^-1, labelled `N1.2' and `N5.6', respectively, the latter being closer to GR.





Figure <ref> shows the relative change (with respect to ΛCDM) of the matter power spectrum due to modifications of gravity, sometimes called `boost', for the two f(R) scenarios at redshifts z=1 (left panels) and z=0.667 (right panels). Since f(R) already exhibits a scale-dependent growth factor at linear order, we can see that on scales of k≈ 0.1 h Mpc^-1 this effect is already present, and gets enhanced at large wavenumbers due to nonlinearities of the density field. As we can see the codes , , , , and all roughly agree to better than one percent down to scales of k≈ 10 h Mpc^-1, as expected. These minor discrepancies can be caused by the refinement criterion used in the latter code or by slight variations in the redshift of the particle snapshot output. The present results are largely in agreement with a similar analysis performed in .

In the same plots, we can also see how approximation schemes to introduce screening from nonlinearities in the potential perform in contrast to the exact solutions. The results that use these methods are showcased by the examples of  and , where each uses a different scheme to approximate the effects of the dynamics of the scalar field in dense environments. As expected, the two codes do not exhibit the same level of agreement down to scales of k≈ 10 h Mpc^-1 as their counterparts using full solvers. However, in both cases, we can see that the deviations from the  reference results are limited to 2% even in the most extreme regime of departure from GR. Through closer inspection, we can see small differences in the agreement between the approximate schemes and full solvers at different values of k. These deviations are caused by different implementations of the approximate MG solvers in the two codes. In fact, while  is an approximate method that uses a PM algorithm to solve the Poisson equation on a fixed grid with the use of the screening approximation to linearise the Poisson equation,  is a new implementation that exploits the TreePM structure of the baseline  code to solve for the small-scale particle dynamics by incorporating the MG effects (including the screening) through a scale-dependent Newton's constant in both real and Fourier space.



Figure <ref> shows results for the matter power spectra in the two nDGP scenarios we consider, characterised by the two values of the cross-over scale, r_c = 1.2 H_0^-1 (N1.2) and r_c = 5.6 H_0^-1 (N5.6). As before, we compare simulation results at two different redshifts, z=1 (left panels) and z=0.667 (right panels). Since the linear growth function has a scale-independent enhancement compared to ΛCDM, we see a constant amplification of power
at small values of k. The Vainshtein mechanism suppresses the deviation from ΛCDM on small scales, leading to a diminishing boost at large k. The agreement between different codes at k < 1 h Mpc^-1 is better than 1% for all codes, including the approximate simulations with  and . At k > 1 h Mpc^-1, the AMR-based codes,  and , show larger deviations at the level of 2%, where the suppression of the deviation from ΛCDM is slightly underestimated compared with . This could be caused by the difference in the ΛCDM power spectrum shown in Fig. <ref>, even though this difference in the baseline code is largely cancelled out in the boost. On the other hand, despite their approximations in the MG solvers,  and  agree with  at the level of 1% even on these scales.



 §.§ Comparison of the halo mass function


To gain further insight into the nonlinear dynamics of the different models and their respective implementations we also compare the cumulative halo mass functions measured in our simulations. For this purpose, halo catalogues are obtained with the  halo finder by running a pipeline described in the second paper of this series (Rácz et al. in prep.). To establish the baseline for the comparison, in the spirit of the previous section, Fig. <ref> shows a comparison of the cumulative halo mass function for the different codes in a ΛCDM cosmology. Since  and  are based on the same TreePM gravity solver, their agreement is excellent. They also agree very well with results from  which, as we like to remind the reader, are only available for the case of ΛCDM and are shown for reference here. On the other hand,  and  are based on the AMR method. They agree with each other, but these simulations underestimate the abundance of low-mass halos below 10^13 h^-1 M_⊙.  uses a fixed-grid PM method, thus it also underestimates the abundance of low-mass halos. With relatively high resolution in these simulations (N_mesh^1/3 = 5 N_p^1/3),  agrees well with  and . We note that exhibits a similar behaviour to ,  and , albeit with a lower amplitude. This discrepancy can be attributed to the fact that has the least small-scale clustering (see Fig. <ref>), resulting in a smaller number of halos within the simulation.



The chameleon screening in f(R) depends on the halo mass such that high-mass halos are typically screened. The critical mass above which screening is effective is determined by the parameter f̅_R. 
As we can see from Fig. <ref>, the ratio of the halo mass function between f(R) and ΛCDM is enhanced for low-mass halos but approaches unity at the high-mass end where all halos are effectively screened. The agreement between full simulations ( and , , ) is around 4% for all halo masses.  uses an approximation for screening based on the thin-shell condition for a spherically symmetric object. Although this captures an overall effect of screening, it fails for low-mass halos in F6 and high-mass halos in F5, leading to larger deviations.



In the case of the Vainshtein mechanism, there is no halo-mass dependence on the screening. Despite screening being effective inside dark matter halos, these halos still feel enhanced gravitational attraction. This increases the merger rate and ultimately leads to larger enhancements of the halo mass function for halos of larger masses. In nDGP, shown in Fig. <ref>, all codes agree within 4%, with a sub-per cent agreement seen in the intermediate mass range 
10^12 h^-1 M_⊙ < M < 10^13 h^-1 M_⊙.  shows a large deviation for M ≲ 5 × 10^11 h^-1 M_⊙. For these masses, there are less than 50 dark matter particles assigned to each halo, which indicates that the results could be affected by the refinement criteria of the simulations. The deviations are larger for the most massive halos, but the number of these halos is low and the halo mass function therefore becomes very noisy in this regime.




§ CONCLUSIONS




In this work, we have presented a comprehensive review of numerical methods for cosmological N-body simulations in scenarios extending beyond the standard ΛCDM model. Our exploration spanned a variety of alternative DE and MG theories, highlighting the critical role of N-body simulations in connecting theoretical models with observational data. Through the detailed examination of numerical solvers and approximations tailored to these extended theories, we have showcased the state of the art of modelling the nonlinear scales of cosmic structure formation under a wide range of cosmological scenarios. Our code comparison exercise, based on the simulations from  and extended by incorporating new codes and approximation techniques, has demonstrated a fair consensus among different numerical implementations. This validation is particularly important for the mission, as the forthcoming observational data will require precise nonlinear modelling to constrain cosmological parameters effectively.

This article is part of a series that explores simulations and nonlinearities in models beyond ΛCDM. The simulation codes that we have considered in this article are used to generate simulation products in the companion paper by Rácz et al. (in preparation) and are crucial for generating the simulations needed to create the nonlinear modelling in the companion paper by Koyama et al. (in preparation), which forecasts constraints for f(R) gravity from photometric primary probes of the mission. The validation checks performed in this paper are therefore critical for being able to trust these results.

The main outcomes of this work can be summarised as follows:

    
  ⋆ N-body simulation codes have been developed for a wide range of extended cosmological scenarios, ranging from simple DE models to more complex interacting scalar field models and MG theories, to non-standard dark matter and initial conditions; the availability of such codes will be a crucial asset for the future developments of large galaxy surveys such as ; we have provided a concise yet comprehensive overview of several such codes, their main features, implementation methods, assumptions and approximations.
    
  ⋆ As a matter of fact, among these simulation codes the ones involving algorithms for the solution of nonlinear differential equations of some additional degrees of freedom, as e.g. for the case of MG theories, are the most challenging in terms of implementation and numerical convergence; we have therefore performed a thorough validation effort of these methods through a code comparison study, extending the approach adopted in  to more recent and diverse algorithms.
    
  ⋆ As a result of our validation effort, we found  agreement in the power spectrum boost at ≲ 1% up to k ≲ 1 h Mpc^-1 and at ≲ 3% up to k ≲ 10 h Mpc^-1 among all the codes implementing full field solvers (), while approximate methods () display slightly larger deviations not exceeding 3% up to k ≲ 7 h Mpc^-1.
    
  ⋆ The halo mass function shows larger deviations among the codes, also due to larger differences in the outcomes for the baseline ΛCDM simulations; nonetheless, all codes agree within less than 5% on the relative change of the halo mass function except for the very largest and the smallest mass ranges where poor statistics and insufficient resolution, respectively, may impact the results.
    
  ⋆ We also compared the computational requirements of the different codes by measuring the CPU time needed to complete a reference MG simulation starting from identical initial conditions; we found that while full field solvers generally imply a substantial increase – up to a factor of ten – of the total CPU time relative to a ΛCDM simulation, approximate solvers are not significantly more demanding for MG simulations compared to standard ΛCDM runs.



Looking forward, the continued evolution of simulation techniques will be paramount in leveraging the full potential of upcoming large-scale structure surveys such as Euclid. N-body simulations therefore continue to set a solid foundation for the next generation of cosmological inquiries. By persistently pushing the boundaries of computational astrophysics, we are poised to uncover the underlying physics driving the accelerated expansion of the Universe, thereby opening new windows onto the fundamental nature of DE, dark matter, and gravity itself.


The lead authors thank A. Schneider and R.E. Smith for their diligent work as internal referees and P. Schneider for carefully proofreading the manuscript.
      The work of JA is supported by the Swiss National Science Foundation.
      During part of this work, AMCLB was supported by a fellowship of PSL University hosted by the Paris Observatory.
      This project was provided with computer and storage resources by GENCI at CINES thanks to the grant 2023-A0150402287 on the supercomputer Adastra's GENOA partition.
      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).
      

aa




§ COMPUTATIONAL COST




One of the obstacles in producing accurate predictions in the nonlinear regime of structure formation is presented by the computational cost of running large N-body simulations. This issue becomes even more pronounced in the case of MG 
simulations, due to the additional computational cost associated with the solution of the Klein–Gordon equation for the scalar field describing the additional degree of freedom in these theories.
To provide some insights into the trade-offs between accuracy and time-to-solution, 
we attempt here a comparative analysis of the computational cost of the simulations run for this paper. Given that these simulations were run on different machines and with different parallelisation settings, we cannot conduct a precise assessment of their computational cost. Instead, we limit this analysis to an order-of-magnitude comparison of the simulations run with the various codes and models discussed in this paper.

For this exercise, we use the information on the wall-clock time, T_ real, and the number of cores, N_ cores, as recorded in the log files of the simulation runs. This information was available for all simulations except for the ones run with . We estimate the computational cost of the simulations, C, as the product of the wall-clock time and the number of cores:

    C ≡ N_ cores  T_ real  .

The estimates of the computational cost for the simulations are compared in Fig. <ref>.
We can see that the cost of ΛCDM simulations of Tree-PM and AMR codes is C ∼ 10^3   CPUh, while for  and  the cost is about one order of magnitude lower at C ∼ 10^2   CPUh. For f(R) gravity models instead, the computational cost increases significantly (approximately by a factor of ten) for the codes that solve the full Klein–Gordon equation of the scalar field, namely , , , and , while the overhead is smaller for the codes that adopt screening approximations, namely  and . Finally, in nDGP gravity, only  has a significant overhead compared to ΛCDM, while  and the approximate codes,  and , have just a small overhead.

The performance of multigrid codes depends on the convergence criterion chosen. For , we chose an extremely conservative approach (see Sect. <ref>) with a very low tolerance threshold, resulting in more V-cycles and almost double the CPU time needed to solve the linear Poisson equation compared to a more standard setup. Regarding the nonlinear solver, we use two F-cycles instead of a single one (which should in principle be enough, but it is not the goal of the present paper to provide a convergence study), therefore roughly doubling the CPU time needed for the f(R) gravity models.

We stress that a thorough assessment of the efficiency of the codes is beyond the scope of this paper and would have required a much more methodical effort including (but not limited to)

    
  * running the simulations in a controlled environment,
    
  * conducting convergence tests for the various hyper-parameters,
    
  * testing the scaling performance of each code.


In fact, when focusing only on predictions of the amplification factors, it is possible to achieve a similar level of accuracy with lower force, mass or time resolution, since resolution effects mostly cancel out when taking ratios of quantities affected by the same inaccuracies <cit.>.
This has been shown to be the case for  simulations in <cit.>, where the use of a lower resolution allowed accurate predictions of power spectrum boosts in nDGP gravity with a theoretical gain of ∼ 300 with respect to the computational cost of the COLA simulations presented here. 

Such large speed-ups have paved the way for creating emulators for the nonlinear amplification of the power spectrum in models beyond ΛCDM in a cost-effective way, i.e. without the need for supercomputers. This has already been done for some of the models we consider in this paper <cit.>. For instance, <cit.> found that an emulator for the nDGP model can be constructed with as little as a few thousand CPUh worth of computational time. Likewise, <cit.> who presented a generic pipeline for using COLA to create such emulators, used f(R) gravity as an example and found similar numbers for the required computational time. <cit.> described a simulation setup that can also be used for emulating the full power spectrum (up to a reasonable high wavenumber k∼ 1 h Mpc^-1), requiring around ∼ 100 CPUh per simulation on a modern CPU. Emulators have also been constructed for beyond-ΛCDM models using high-resolution direct simulations in the same way as has been done for ΛCDM. This approach generally gives more accurate emulators than those created with approximate methods, but this comes at a much higher cost. For example, both <cit.> and <cit.> have each presented a high-fidelity emulator for the f(R) model considered in this paper, but at a higher cost of about 3.5-4 million CPUh.

