# Euclid: Modelling massive neutrinos in cosmology – a code comparison

J. Adamek, R.E. Angulo, C. Arnold, M. Baldi, M. Biagetti, B. Bose, C.

Carbone, T. Castro, J. Dakin, K. Dolag, et al.

► To cite this version:

J. Adamek, R.E. Angulo, C. Arnold, M. Baldi, M. Biagetti, et al. Euclid: Modelling massive neutrinos in cosmology – a code comparison. JCAP, 2023, 06, pp.035. 10.1088/1475-7516/2023/06/035. hal-03892198

HAL Id: hal-03892198

https://hal.science/hal-03892198

Submitted on 9 Nov 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Euclid: modelling massive neutrinos in cosmology — a code comparison

To cite this article: On behalf of Euclid consortium JCAP06(2023)035

View the [article online](#) for updates and enhancements.

You may also like

- [Cosmological Constant and Stability of the Cosmological Models](#)
  Wang Yong-Jiu and Tang Zhi-Ming

- [The GIGANTES Data Set: Precision Cosmology from Voids in the Machine-Learning Era](#)
  Christina D. Kreisch, Alice Pisani, Francisco Villaescusa-Navarro et al.

- [THE ATACAMA COSMOLOGY TELESCOPE: COSMOLOGICAL PARAMETERS FROM THE 2008 POWER SPECTRUM](#)
  J. Dunkley, R. Hlozek, J. Sievers et al.
# Euclid: modelling massive neutrinos in cosmology — a code comparison

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

## On behalf of Euclid consortium

J. Adamek,¹ R.E. Angulo,²,³ C. Arnold,⁴ M. Baldi,⁵,⁶,⁷

M. Biagetti,⁸,⁹,¹⁰,¹¹ B. Bose,¹²,¹³,¹ C. Carbone,¹⁴

T. Castro,⁹,¹⁰,⁸ J. Dakin,¹,¹⁵ K. Dolag,¹⁶,¹⁷ W. Elbers,⁴

C. Fidler,¹⁸ C. Giocoli,⁶,¹⁹ S. Hannestad,¹⁵ F. Hassani,²⁰

C. Hernández-Aguayo,¹⁷ K. Koyama,²¹ B. Li,⁴

R. Mauland,²⁰ P. Monaco,²²,⁹,¹⁰,⁸ C. Moretti,¹²,⁹

D.F. Mota,²⁰ C. Partmann,¹⁷ G. Parimbelli,²³,⁹,¹¹,⁸ D. Potter,¹

A. Schneider,¹ S. Schulz,¹ R.E. Smith,²⁴ V. Springel,¹⁷

J. Stadel,¹ T. Tram,¹⁵ M. Viel,¹¹,⁸,⁹,¹⁰

F. Villaescusa-Navarro,²⁵,²⁶ H.A. Winther,²⁰ B.S. Wright,²⁷

M. Zennaro,² N. Aghanim,²⁸ L. Amendola,²⁹ N. Auricchio,⁶

D. Bonino,³⁰ E. Branchini,³¹,³² M. Brescia,³³,³⁴

S. Camera,³⁵,³⁶,³⁰ V. Capobianco,³⁰ V.F. Cardone,³⁷,³⁸

J. Carretero,³⁹,⁴⁰ F.J. Castander,⁴¹,⁴² M. Castellano,³⁷

S. Cavuoti,³⁴,⁴³ A. Cimatti,⁴⁴,⁴⁵ R. Cledassou,⁴⁶,⁴⁷

G. Congedo,¹² L. Conversi,⁴⁸,⁴⁹ Y. Copin,⁵⁰

A. Da Silva,⁵¹,⁵² H. Degaudenzi,⁵³ M. Douspis,²⁸

F. Dubath,⁵³ C.A.J. Duncan,⁵⁴,⁵⁵ X. Dupac,⁴⁸ S. Dusini,⁵⁶

S. Farrens,⁵⁷ S. Ferriol,⁵⁰ P. Fosalba,⁴²,⁴¹ M. Frailis,⁹

E. Franceschi,⁶ S. Galeotta,⁹ B. Garilli,¹⁴ W. Gillard,⁵⁸

B. Gillis,¹² A. Grazian,⁵⁹ S.V. Haugan,²⁰ W. Holmes,⁶⁰

A. Hornstrup,⁶¹,⁶² K. Jahnke,⁶³ S. Kermiche,⁵⁸

A. Kiessling,⁶⁰ M. Kilbinger,⁵⁷ T. Kitching,⁶⁴ M. Kunz,¹³

H. Kurki-Suonio,⁶⁵,⁶⁶ P.B. Lilje,²⁰ I. Lloro,⁶⁷ O. Mansutti,⁹

O. Marggraf,⁶⁸ F. Marulli,⁵,⁶,⁷ R. Massey,⁴ E. Medinaceli,⁶

M. Meneghetti,⁶,⁷ G. Meylan,⁶⁹ M. Moresco,⁵,⁶

L. Moscardini,⁵,⁶,⁷ E. Munari,⁹ S.-M. Niemi,⁷⁰ C. Padilla,³⁹

<figure><img src="image_3.png" /></figure>
S. Paltani,⁵³ F. Pasian,<fnref n="9" /> K. Pedersen,<fnref n="15" /> W.J. Percival ID,⁷¹,⁷²,⁷³

V. Pettorino,⁵⁷ G. Polenta ID,⁷⁴ M. Poncet,⁴⁶ L.A. Popa,⁷⁵

F. Raison ID,⁷⁶ R. Rebolo,⁷⁷,⁷⁸ A. Renzi ID,⁷⁹,⁵⁶ J. Rhodes,⁶⁰

G. Riccio,³⁴ E. Romelli ID,<fnref n="9" /> M. Roncarelli,<fnref n="6" /> R. Saglia ID,<fnref n="16" />,⁷⁶

D. Sapone ID,⁸⁰ B. Sartoris,<fnref n="16" />,<fnref n="9" /> P. Schneider ID,⁶⁸

T. Schrabback ID,⁶⁸,⁸¹ A. Secroun ID,⁵⁸ G. Seidel ID,⁶³

C. Sirignano ID,⁷⁹,⁵⁶ G. Sirri ID,<fnref n="7" /> L. Stanco ID,⁵⁶ J.-L. Starck,⁸²

P. Tallada-Crespí,⁸³,⁴⁰ A.N. Taylor,<fnref n="12" /> I. Tereno,⁵¹,⁸⁴

R. Toledo-Moreo ID,⁸⁵ F. Torradeflot ID,⁸³,⁴⁰ I. Tutusaus ID,⁸⁶,<fnref n="13" />

L. Valenziano ID,<fnref n="6" />,<fnref n="7" /> T. Vassallo ID,<fnref n="9" /> Y. Wang ID,⁸⁷ J. Weller ID,<fnref n="16" />,⁷⁶

A. Zacchei ID,<fnref n="9" />,<fnref n="8" /> G. Zamorani ID,<fnref n="6" /> J. Zoubian,⁵⁸ G. Fabbian ID,²⁵,⁸⁸
and V. Scottez⁸⁹,⁹⁰

<footnote n="1">¹Institute for Computational Science, Universität Zürich,
Winterthurerstrasse 190, 8057 Zürich, Switzerland</footnote>

<footnote n="2">²Donostia International Physics Center (DIPC),
Paseo Manuel de Lardizabal 4, 20018 Donostia-San Sebastián, Guipuzkoa, Spain</footnote>

<footnote n="3">³IKERBASQUE, Basque Foundation for Science, 48013 Bilbao, Spain</footnote>

<footnote n="4">⁴Department of Physics, Institute for Computational Cosmology, Durham University,
South Road, Durham DH1 3LE, U.K.</footnote>

<footnote n="5">⁵Dipartimento di Fisica e Astronomia “Augusto Righi” — Alma Mater Studiorum Università
di Bologna, via Piero Gobetti 93/2, 40129 Bologna, Italy</footnote>

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

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

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

<footnote n="9">⁹INAF-Osservatorio Astronomico di Trieste, Via G. B. Tiepolo 11, 34143 Trieste, Italy</footnote>

<footnote n="10">¹⁰INFN, Sezione di Trieste, Via Valerio 2, 34127 Trieste, Italy</footnote>

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

<footnote n="12">¹²Institute for Astronomy, University of Edinburgh, Royal Observatory,
Blackford Hill, Edinburgh EH9 3HJ, U.K.</footnote>

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

<footnote n="14">¹⁴INAF-IASF Milano, Via Alfonso Corti 12, 20133 Milano, Italy</footnote>

<footnote n="15">¹⁵Department of Physics and Astronomy, University of Aarhus,
Ny Munkegade 120, 8000 Aarhus C, Denmark</footnote>

<footnote n="16">¹⁶Universitäts-Sternwarte München, Fakultät für Physik,
Ludwig-Maximilians-Universität München, Scheinerstrasse 1, 81679 München, Germany</footnote>

<footnote n="17">¹⁷Max-Planck-Institut für Astrophysik, Karl-Schwarzschild Str. 1, 85741 Garching, Germany</footnote>

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

<footnote n="19">¹⁹Istituto Nazionale di Fisica Nucleare, Sezione di Bologna,
Via Irnerio 46, 40126 Bologna, Italy</footnote>
<sup>20</sup>Institute of Theoretical Astrophysics, University of Oslo,
P.O. Box 1029 Blindern, 0315 Oslo, Norway

<sup>21</sup>Institute of Cosmology and Gravitation, University of Portsmouth,
Portsmouth PO1 3FX, U.K.

<sup>22</sup>Dipartimento di Fisica — Sezione di Astronomia, Università di Trieste,
Via Tiepolo 11, 34131 Trieste, Italy

<sup>23</sup>Dipartimento di Fisica, Università degli studi di Genova, and INFN-Sezione di Genova,
via Dodecaneso 33, 16146 Genova, Italy

<sup>24</sup>Department of Physics & Astronomy, University of Sussex, Brighton BN1 9QH, U.K.

<sup>25</sup>Center for Computational Astrophysics, Flatiron Institute,
162 5th Avenue, New York, NY 10010, U.S.A.

<sup>26</sup>Department of Astrophysical Sciences, Peyton Hall, Princeton University,
Princeton, NJ 08544, U.S.A.

<sup>27</sup>School of Physics and Astronomy, Queen Mary University of London,
Mile End Road, London E1 4NS, U.K.

<sup>28</sup>Université Paris-Saclay, CNRS, Institut d'astrophysique spatiale, 91405 Orsay, France

<sup>29</sup>Institut für Theoretische Physik, University of Heidelberg,
Philosophenweg 16, 69120 Heidelberg, Germany

<sup>30</sup>INAF-Osservatorio Astrofisico di Torino, Via Osservatorio 20, 10025 Pino Torinese, Italy

<sup>31</sup>Dipartimento di Fisica, Università di Genova, Via Dodecaneso 33, 16146 Genova, Italy

<sup>32</sup>INFN-Sezione di Roma Tre, Via della Vasca Navale 84, 00146 Roma, Italy

<sup>33</sup>Department of Physics “E. Pancini”, University Federico II,
Via Cinthia 6, 80126 Napoli, Italy

<sup>34</sup>INAF-Osservatorio Astronomico di Capodimonte, Via Moiariello 16, 80131 Napoli, Italy

<sup>35</sup>Dipartimento di Fisica, Università degli Studi di Torino,
Via P. Giuria 1, 10125 Torino, Italy

<sup>36</sup>INFN-Sezione di Torino, Via P. Giuria 1, 10125 Torino, Italy

<sup>37</sup>INAF-Osservatorio Astronomico di Roma,
Via Frascati 33, 00078 Monteporzio Catone, Italy

<sup>38</sup>INFN-Sezione di Roma, Piazzale Aldo Moro, 2 — c/o Dipartimento di Fisica,
Edificio G. Marconi, 00185 Roma, Italy

<sup>39</sup>Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science
and Technology, Campus UAB, 08193 Bellaterra (Barcelona), Spain

<sup>40</sup>Port d'Informació Científica, Campus UAB,
C. Albareda s/n, 08193 Bellaterra (Barcelona), Spain

<sup>41</sup>Institut d'Estudis Espacials de Catalunya (IEEC),
Carrer Gran Capitá 2-4, 08034 Barcelona, Spain

<sup>42</sup>Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s/n, 08193
Barcelona, Spain

<sup>43</sup>INFN section of Naples, Via Cinthia 6, 80126 Napoli, Italy

<sup>44</sup>Dipartimento di Fisica e Astronomia “Augusto Righi” — Alma Mater Studiorum Università
di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy

<sup>45</sup>INAF-Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy

<sup>46</sup>Centre National d'Etudes Spatiales, Toulouse, France
<sup>47</sup>Institut national de physique nucléaire et de physique des particules,
3 rue Michel-Ange, 75794 Paris CEDEX 16, France

<sup>48</sup>ESAC/ESA, Camino Bajo del Castillo, s/n., Urb. Villafranca del Castillo,
28692 Villanueva de la Cañada, Madrid, Spain

<sup>49</sup>European Space Agency/ESRIN, Largo Galileo Galilei 1, 00044 Frascati, Roma, Italy

<sup>50</sup>Université Lyon, Université Claude Bernard Lyon 1, CNRS/IN2P3, IP2I Lyon,
UMR 5822, 69622 Villeurbanne, France

<sup>51</sup>Departamento de Física, Faculdade de Ciências, Universidade de Lisboa,
Edifício C8, Campo Grande, 1749-016 Lisboa, Portugal

<sup>52</sup>Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências,
Universidade de Lisboa, Campo Grande, 1749-016 Lisboa, Portugal

<sup>53</sup>Department of Astronomy, University of Geneva, ch. d'Ecogia 16, 1290 Versoix, Switzerland

<sup>54</sup>Department of Physics, Oxford University, Keble Road, Oxford OX1 3RH, U.K.

<sup>55</sup>Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy,
University of Manchester, Oxford Road, Manchester M13 9PL, U.K.

<sup>56</sup>INFN-Padova, Via Marzolo 8, 35131 Padova, Italy

<sup>57</sup>Université Paris-Saclay, Université Paris Cité, CEA, CNRS, Astrophysique,
Instrumentation et Modélisation Paris-Saclay, 91191 Gif-sur-Yvette, France

<sup>58</sup>Aix-Marseille Université, CNRS/IN2P3, CPPM, Marseille, France

<sup>59</sup>INAF-Osservatorio Astronomico di Padova, Via dell'Osservatorio 5, 35122 Padova, Italy

<sup>60</sup>Jet Propulsion Laboratory, California Institute of Technology,
4800 Oak Grove Drive, Pasadena, CA 91109, U.S.A.

<sup>61</sup>Technical University of Denmark, Elektrovej 327, 2800 Kgs. Lyngby, Denmark

<sup>62</sup>Cosmic Dawn Center (DAWN), Denmark

<sup>63</sup>Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany

<sup>64</sup>Mullard Space Science Laboratory, University College London,
Holmbury St Mary, Dorking, Surrey RH5 6NT, U.K.

<sup>65</sup>Department of Physics, P.O. Box 64, 00014 University of Helsinki, Finland

<sup>66</sup>Helsinki Institute of Physics, Gustaf Hällströmin katu 2, University of Helsinki,
Helsinki, Finland

<sup>67</sup>NOVA optical infrared instrumentation group at ASTRON,
Oude Hoogeveensedijk 4, 7991PD, Dwingeloo, The Netherlands

<sup>68</sup>Argelander-Institut für Astronomie, Universität Bonn,
Auf dem Hügel 71, 53121 Bonn, Germany

<sup>69</sup>Institute of Physics, Laboratory of Astrophysics, Ecole Polytechnique Fédérale de Lausanne
(EPFL), Observatoire de Sauverny, 1290 Versoix, Switzerland

<sup>70</sup>European Space Agency/ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands

<sup>71</sup>Centre for Astrophysics, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada

<sup>72</sup>Department of Physics and Astronomy, University of Waterloo,
Waterloo, Ontario N2L 3G1, Canada

<sup>73</sup>Perimeter Institute for Theoretical Physics, Waterloo, Ontario N2L 2Y5, Canada

<sup>74</sup>Space Science Data Center, Italian Space Agency,
via del Politecnico snc, 00133 Roma, Italy

<sup>75</sup>Institute of Space Science, Bucharest 077125, Romania
<sup>76</sup>Max Planck Institute for Extraterrestrial Physics,
Giessenbachstr. 1, 85748 Garching, Germany

<sup>77</sup>Instituto de Astrofísica de Canarias,
Calle Vía Láctea s/n, 38204 San Cristóbal de La Laguna, Tenerife, Spain

<sup>78</sup>Departamento de Astrofísica, Universidad de La Laguna, 38206 La Laguna, Tenerife, Spain

<sup>79</sup>Dipartimento di Fisica e Astronomia “G. Galilei”, Università di Padova,
Via Marzolo 8, 35131 Padova, Italy

<sup>80</sup>Departamento de Física, FCFM, Universidad de Chile,
Blanco Encalada 2008, Santiago, Chile

<sup>81</sup>Institut für Astro- und Teilchenphysik, Universität Innsbruck,
Technikerstr. 25/8, 6020 Innsbruck, Austria

<sup>82</sup>AIM, CEA, CNRS, Université Paris-Saclay, Université de Paris,
91191 Gif-sur-Yvette, France

<sup>83</sup>Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT),
Avenida Complutense 40, 28040 Madrid, Spain

<sup>84</sup>Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências,
Universidade de Lisboa, Tapada da Ajuda, 1349-018 Lisboa, Portugal

<sup>85</sup>Universidad Politécnica de Cartagena, Departamento de Electrónica y Tecnología
de Computadoras, 30202 Cartagena, Spain

<sup>86</sup>Institut de Recherche en Astrophysique et Planétologie (IRAP), Université de Toulouse,
CNRS, UPS, CNES, 14 Av. Edouard Belin, 31400 Toulouse, France

<sup>87</sup>Infrared Processing and Analysis Center, California Institute of Technology,
Pasadena, CA 91125, U.S.A.

<sup>88</sup>School of Physics and Astronomy, Cardiff University, The Parade, Cardiff CF24 3AA, U.K.

<sup>89</sup>Institut d'Astrophysique de Paris, 98bis Boulevard Arago, 75014 Paris, France

<sup>90</sup>Junia, EPA department, 59000 Lille, France

E-mail: julian.adamek@uzh.ch

Received November 23, 2022

Accepted April 27, 2023

Published June 16, 2023
**Abstract.** The measurement of the absolute neutrino mass scale from cosmological large-scale clustering data is one of the key science goals of the *Euclid* mission. Such a measurement relies on precise modelling of the impact of neutrinos on structure formation, which can be studied with *N*-body simulations. Here we present the results from a major code comparison effort to establish the maturity and reliability of numerical methods for treating massive neutrinos. The comparison includes eleven full *N*-body implementations (not all of them independent), two *N*-body schemes with approximate time integration, and four additional codes that directly predict or emulate the matter power spectrum. Using a common set of initial data we quantify the relative agreement on the nonlinear power spectrum of cold dark matter and baryons and, for the *N*-body codes, also the relative agreement on the bispectrum, halo mass function, and halo bias. We find that the different numerical implementations produce fully consistent results. We can therefore be confident that we can model the impact of massive neutrinos at the sub-percent level in the most common summary statistics. We also provide a code validation pipeline for future reference.

**Keywords:** cosmological neutrinos, cosmological simulations, neutrino masses from cosmology, power spectrum

ArXiv ePrint: 2211.12457
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>Neutrino physics</td><td>3</td></tr><tr><td>3</td><td>Numerical methods</td><td>4</td></tr><tr><td>3.1</td><td>Particle-based methods</td><td>4</td></tr><tr><td>3.2</td><td>Mesh-based methods</td><td>6</td></tr><tr><td>3.3</td><td>Approximations and other methods</td><td>7</td></tr><tr><td>3.3.1</td><td>COLA</td><td>8</td></tr><tr><td>3.3.2</td><td>PINOCCHIO</td><td>8</td></tr><tr><td>3.3.3</td><td>Newtonian motion gauge</td><td>10</td></tr><tr><td>3.4</td><td>Halo-model reaction</td><td>10</td></tr><tr><td>3.5</td><td>Power-spectrum emulation</td><td>12</td></tr><tr><td>4</td><td>Simulations</td><td>13</td></tr><tr><td>4.1</td><td>Initial conditions</td><td>14</td></tr><tr><td>4.2</td><td>Post-processing pipeline</td><td>15</td></tr><tr><td>4.2.1</td><td>Power spectra</td><td>15</td></tr><tr><td>4.2.2</td><td>Bispectra</td><td>16</td></tr><tr><td>4.2.3</td><td>Halo catalogues</td><td>16</td></tr><tr><td>5</td><td>Results</td><td>17</td></tr><tr><td>5.1</td><td>Power spectra</td><td>17</td></tr><tr><td>5.1.1</td><td>CDM and baryons</td><td>17</td></tr><tr><td>5.1.2</td><td>Convergence tests</td><td>21</td></tr><tr><td>5.1.3</td><td>Contributions from neutrinos</td><td>22</td></tr><tr><td>5.1.4</td><td>Total matter</td><td>23</td></tr><tr><td>5.2</td><td>Bispectra</td><td>23</td></tr><tr><td>5.3</td><td>Halo mass function</td><td>29</td></tr><tr><td>5.4</td><td>Halo bias</td><td>31</td></tr><tr><td>6</td><td>Discussion</td><td>34</td></tr></table>

## 1 Introduction

The upcoming *Euclid* mission [1] will provide very detailed observations of the large-scale structure of our Universe, making it possible to probe physics related to dark energy and neutrinos at an unprecedented level of precision. The analysis and interpretation of these data require a very accurate modelling of the process of structure formation. This is of particular relevance since a precise modelling of the mass-dependent effect neutrinos have on various summary statistics will allow a cosmological measurement of the absolute neutrino mass scale, which is one of the key science objectives of *Euclid*.

Here we focus on the treatment of massive neutrinos in cosmological N-body simulations and investigate the convergence of a number of different codes over a variety of different scales
and redshifts. Oscillation experiments have established a firm lower bound on the sum of neutrino masses of around 0.06 eV. Using the well-known relation between neutrino mass and cosmological energy density, $\sum m_\nu \approx \Omega_\nu h^2 \times 94$ eV, this lower bound on the sum of neutrino masses corresponds to a lower bound of $\Omega_\nu h^2 > 6 \times 10^{-4}$, or approximately 0.5% of the total matter density. As usual, the cosmological density of any component X can be given in terms of $\Omega_X$, which is its present-day energy density in units of the critical density, and physical density parameters are then denoted as $\omega_X = \Omega_X h^2$, where h is the reduced Hubble parameter.

The sum of the neutrino masses is already constrained using information from the cosmic microwave background (CMB) combined with observations of large-scale structure like baryon acoustic oscillations [2], redshift-space distortions [3], and the Lyman-α forest [4]. While currently providing only upper bounds, these constraints are expected to improve significantly with upcoming surveys like *Euclid* which will be able to measure the neutrino mass fraction even if it is close to the lower bound. The matter power spectrum is still affected at the level of 4% in this scenario, which is well within the sensitivity of the *Euclid* main probes.

For instance, the weak-lensing signal probed by *Euclid* is sensitive to the matter power spectrum up to $k \approx 7 h \text{ Mpc}^{-1}$ [5]. A linear model would be completely inadequate at such short scales and we therefore need robust nonlinear models. The forecast for galaxy clustering in *Euclid* typically assumes that the matter power spectrum and the galaxy bias is well understood at least up to $k \approx 0.3 h \text{ Mpc}^{-1}$ which also requires some nonlinear prescription [6]. Many codes considered here are used in *Euclid* preparation papers and reference simulations. For instance, *PKDGRAV3* ran *Euclid*'s flagship simulations and the simulations used by Knabenhans et al. [7, 8] to train the *EuclidEmulator2*; *openGADGET3* was used to calibrate the halo mass function for the *Euclid* cluster abundance analysis by Castro et al. [9], and *PINOCCHIO* was used to create the synthetic catalogues for the validation of the covariance matrix of cluster abundance and the clustering of clusters by Fumagalli et al. [10, 11]. Our objective is therefore to establish a reliable calibration baseline for the measurement of the neutrino mass scale within the cosmological analysis of *Euclid* data. We expect that our results are also relevant in the context of other so-called “stage IV surveys” like the *Vera C. Rubin Observatory* or the *Nancy Grace Roman Space Telescope*.

Over the past decade, a variety of different methods have been employed to incorporate massive neutrinos in N-body simulations [12]. As discussed in more detail later, they broadly fall into two categories which we may refer to as “particle based” and “mesh based,” respectively. They follow different philosophies of keeping track of the evolving neutrino phase-space distribution function. Particle-based methods sample the six-dimensional phase-space directly, see e.g. refs. [13–22]. Mesh-based methods, on the other hand, work under the approximation that neutrino perturbations remain small and can be treated with perturbation theory [23–26]. In the simplest case that is sufficient for low neutrino masses one works with a linear realisation of the neutrino density field on a grid [23]. The linear theory for neutrinos may also be solved using the full nonlinear gravitational potential calculated in the simulation [24–26]. A different approach is to treat the massive neutrinos as a fluid and then solve the corresponding fluid equations, employing some approximation scheme to close the set of equations [27]. One can also attempt to integrate the Vlasov-Poisson equations on a six-dimensional phase-space grid [28].

Generally speaking, the mesh-based schemes work best for relatively small neutrino masses where neutrino perturbations remain linear or quasi-linear at all times. There are also hybrid schemes that use elements from both approaches. For instance, one may use a linear mesh-based representation at early times which is then converted to a particle representation
at late time [29, 30]. The so-called “$\delta f$ method” introduced in Elbers et al. [31] is another hybrid approach that uses a particle ensemble to estimate perturbations $\delta f$ to the smooth background phase-space distribution function $f$. Finally, a coordinate (gauge) transformation can be used to include linear neutrino perturbations without modifications to the N-body simulation code [32]. Cosmology-rescaling algorithms that are applied in post-processing have been shown to provide accurate results as well [33]. In this work, we aim to compare these different numerical approaches by employing them to run cosmological N-body simulations starting from the same initial conditions, and comparing the properties of the resulting matter and halo distributions using a controlled post-processing pipeline.

This paper is structured as follows. We begin with a brief review of neutrino physics and its impact on cosmology in section 2. In section 3, we describe the numerical methods that can be used to account for the cosmological effects of neutrinos. Our simulations are described in section 4 and in section 5 we present our numerical results. We conclude with a discussion in section 6.

## 2 Neutrino physics

From oscillation experiments it is firmly established that at least two of the standard-model neutrino mass states have non-zero mass, but the absolute mass scale is unknown and two mass orderings remain possible: normal and inverted. The current best-fit values for the mass-square differences measured in oscillation experiments are given by [34]

$$ \Delta m_{21}^2 = 7.42_{-0.20}^{+0.21} \times 10^{-5} \text{ eV}^2, \qquad (2.1) $$

$$ \Delta m_{31}^2 = 2.517_{-0.028}^{+0.026} \times 10^{-3} \text{ eV}^2 \quad (\text{NO}), \qquad (2.2) $$

$$ \Delta m_{32}^2 = -2.498_{-0.028}^{+0.028} \times 10^{-3} \text{ eV}^2 \quad (\text{IO}), \qquad (2.3) $$

where “NO” denotes the normal mass ordering and “IO” the inverted one. This leads to lower bounds on the sum of neutrino masses of $\sum m_\nu \gtrsim 0.06 \text{ eV}$ (NO) and $\sum m_\nu \gtrsim 0.1 \text{ eV}$ (IO), respectively. The best current experimental upper bound on the neutrino masses comes from the KATRIN experiment, which measures an incoherent sum of mass states using beta decay of tritium [35, 36]. This bound approximately translates to $\sum m_\nu \lesssim 2.4 \text{ eV}$.

However, cosmology already provides much more stringent bounds, typically around $\sum m_\nu \lesssim 0.1 - 0.2 \text{ eV}$. Assuming a minimal $\Lambda$CDM cosmology with massive neutrinos, a joint analysis of cosmological probes currently obtains $\sum m_\nu < 0.09 \text{ eV}$ at 95% confidence, already putting pressure on the inverted mass ordering scenario [37]. These constraints are based on several physical effects affecting the CMB and large-scale structure in different ways, see Archidiacono et al. [38] for a review. Some of the effects are simply related to the change in the expansion history, others to the explicit coupling of neutrinos to cosmological perturbations. More specifically, massive neutrinos modify the shape of the matter power spectrum both in the linear and the nonlinear regimes. First, as neutrinos behave like radiation in the early Universe, they move the radiation-matter equality to slightly later times, therefore shifting the peak of the power spectrum towards smaller wavenumbers. Second, after the non-relativistic transition, they slow down the linear growth of perturbations at scales smaller than the free-streaming length, leading to a scale-dependent growth rate. The small-scale suppression in the linear power spectra of cold dark matter (CDM) and baryons, $P_{\text{cb}}$, or total matter (which includes massive neutrinos), $P_{\text{m}}$, with respect to a model with massless
neutrinos, can be quantified as [18, 39, 40]

$$ \frac{\Delta P_{cb}}{P_{cb}} \approx 6 f_{\nu}, \qquad \frac{\Delta P_{m}}{P_{m}} \approx 8 f_{\nu}, \qquad (2.4) $$

respectively, where $f_{\nu} = \Omega_{\nu}/\Omega_{m}$ is the neutrino mass fraction. In the nonlinear regime, this suppression is even more prominent and exhibits a dip at $k \approx 1 h \text{ Mpc}^{-1}$ for low redshift, giving rise to the well-known “spoon-like” feature [41]. In the context of the halo model, this feature appears at the transition region of the two-halo and the one-halo terms. In particular, the dip is caused by the small-scale suppression of two-halo clustering that is induced by free-streaming, while the subsequent rise reflects the fact that the number of the most massive halos is rather independent of the neutrino masses. All these effects can be accurately predicted by modelling the neutrino component in cosmological N-body simulations. Assessing the relative accuracy and convergence of such modelling over a range of different numerical methods and simulation codes is the main goal of this paper.

## 3 Numerical methods

In this section, we give an overview of the various methods that have been developed for the treatment of massive neutrinos in N-body simulations and other numerical models. The most accurate results are expected when the local density of neutrinos in configuration space is accounted for within a simulation itself. This is technically challenging because of the large phase-space volume that is populated by neutrino particles. The methods to deal with this broadly fall into two categories that shall be discussed in turn: particle-based and mesh-based. Hybrid methods that use concepts from both categories have also been developed. Apart from full N-body simulations that may try to incorporate (as much as possible) the neutrino physics, there also exist approximate methods to generate realisations of large-scale structure. These can be augmented with recipes to account for the effect of massive neutrinos. Finally, if one is only interested in summary statistics like the power spectrum, emulators are a powerful tool that can be calibrated to include the sum of the neutrino masses as a free parameter. An overview of the various numerical codes employed in this work is given in table 1.

### 3.1 Particle-based methods

Conceptually, the most straightforward method of accounting for cosmic neutrinos in a simulation is to represent them by a separate N-body ensemble. However, this method faces several challenges that need to be addressed carefully.

The first challenge is posed by the aforementioned phase-space volume that needs to be sampled. The phase-space distribution of neutrino particles typically has a very large velocity dispersion that is orders of magnitude larger than the bulk velocity. Therefore, representing the neutrinos with a small number of N-body particles that simply track the bulk velocities, which is essentially the method of choice for CDM or baryons, would completely miss the fact that most neutrinos are unbound and easily free stream out of gravity wells. Thus, the common practice is to sample the N-body particles from the true phase-space distribution, effectively performing a Monte-Carlo integration of the evolution equations. The main drawback of this method is that a poor sampling usually introduces significant shot noise, while high sampling rates quickly become very expensive as both the memory requirement and computations become completely dominated by the neutrino particle load. This is undesirable since neutrinos are just a tiny fraction of the matter after all.
<table><caption>**Table 1.** Overview of the numerical codes employed in this code comparison. N-body codes typically use a particle-mesh (PM) method coupled to some scheme to increase the force resolution in high-density regions. Methods featured here include Tree-PM, adaptive mesh refinement (AMR), fast multipole method (FMM), particle-particle/particle-mesh (P<sup>3</sup>M), and moving mesh. No adaptive force computation is used in *gevolution*, *COLA*, and *PINOCCHIO*, which all work with a uniform mesh.</caption><thead><tr><th>Code</th><th>type</th><th>neutrino method(s)</th><th>reference(s)</th></tr></thead><tbody><tr><td>GADGET-3</td><td>N-body (Tree-PM)</td><td>particle</td><td>[42, 43]</td></tr><tr><td>L-GADGET3</td><td>N-body (Tree-PM)</td><td>mesh</td><td>[43–45]</td></tr><tr><td>openGADGET3</td><td>N-body (Tree-PM)</td><td>particle</td><td>[46, 47]</td></tr><tr><td>GADGET-4</td><td>N-body (Tree-PM)</td><td>particle</td><td>[48]</td></tr><tr><td>NM-GADGET4</td><td>N-body (Tree-PM)</td><td>Newtonian motion gauge</td><td>[32, 49]</td></tr><tr><td>AREPO</td><td>N-body (Tree-PM)</td><td>particle</td><td>[50, 51]</td></tr><tr><td>CONCEPT</td><td>N-body (P<sup>3</sup>M)</td><td>mesh</td><td>[27, 52]</td></tr><tr><td>PKDGRAV3</td><td>N-body (Tree + FMM)</td><td>mesh</td><td>[53]</td></tr><tr><td>SWIFT</td><td>N-body (PM + FMM)</td><td>particle / δf</td><td>[31, 54, 55]</td></tr><tr><td>ANUBIS</td><td>N-body (PM + AMR)</td><td>particle</td><td>[56, 57]</td></tr><tr><td>gevolution</td><td>N-body (uniform PM)</td><td>particle / mesh</td><td>[20, 58, 59]</td></tr><tr><td>COLA</td><td>N-body surrogate</td><td>mesh</td><td>[60, 61]</td></tr><tr><td>PINOCCHIO</td><td>N-body surrogate</td><td>linear growth factor</td><td>[62, 63]</td></tr><tr><td>ReACT</td><td>P(k) prediction</td><td>halo-model reaction</td><td>[64, 65]</td></tr><tr><td>BACCOemulator</td><td>P(k) prediction</td><td>emulation</td><td>[45]</td></tr><tr><td>EuclidEmulator2</td><td>P(k) prediction</td><td>emulation</td><td>[8]</td></tr><tr><td>CosmicEmu</td><td>P(k) prediction</td><td>emulation</td><td>[66, 67]</td></tr></tbody></table>

Shot noise affects all moments of the distribution function, and in particular the density. This means that shot noise will also propagate into the gravitational field computed from the neutrino perturbations. While this could in principle severely degrade the accuracy of the gravitational evolution as a whole, the impact is actually mitigated by the fact that neutrinos only account for a very small fraction of the total matter, and the gravitational fields are therefore dominated by cold species. Nevertheless, some shot noise does propagate into the other matter species, particularly on large scales where the contribution of neutrino perturbations is largest. Various strategies have been developed to reduce the impact of shot noise, e.g. by filtering small-scale fluctuations. One effective strategy is to implement a statistical weighting of the neutrino particles, as is done in the *δf* method [31]. This method works by decomposing the distribution function *f* into an analytical background component *f̃* and a perturbation *δf* computed from the particle ensemble. The weights are given at each time step by requiring phase-space density conservation. These weights are typically negligible, except for particles that are significantly perturbed, such as those captured by halos. Shot noise is thereby minimised as particles only contribute to the gravitational potential when needed. The *δf* method has been implemented in a number of codes, but is only used by *SWIFT* in this comparison (see table 1). Other strategies aimed at reducing shot noise while still using particles include alternate sampling of neutrino momenta [21] and various hybrid methods [29, 30, 68].
The second challenge pertains to the kinematics of neutrino N-body particles. If one were to apply a similar time-stepping criterion as for cold matter species, the high velocities would typically result in extremely small integration time steps, making the simulations considerably more expensive. This is often solved by relaxing the time-stepping criterion and allowing the neutrino particles to travel a larger distance in each integration step than what would be considered “safe” for cold species. High-velocity neutrino particles then respond poorly to small-scale fluctuations in the gravitational forces. Given that the neutrino distribution at small scales is plagued by the shot-noise problem anyway, this additional problem is often considered to be of little concern.

Still regarding kinematics, further issues arise due to the relativistic nature of neutrinos. In the weak-field limit, the propagation of a collisionless massive particle is governed by the Hamiltonian equations of motion [69]

$$v = \frac{pc}{\sqrt{p^2 + m^2c^2a^2}}, \qquad (3.1)$$

$$p' = -\frac{2p^2 + m^2c^2a^2}{c\sqrt{p^2 + m^2c^2a^2}} \nabla\Psi, \qquad (3.2)$$

where $v$ is the peculiar velocity, $p$ is the canonical momentum, $m$ is the particle's rest-mass, $a$ is the scale factor, and $\Psi$ is the gravitational potential, assuming as usual that gravitational slip can be neglected, i.e. that non-relativistic and ultra-relativistic particles essentially see the same potential. Here, a prime denotes the derivative with respect to conformal time and $c$ denotes the speed of light. The canonical momentum is conserved in the absence of a gravitational force. For non-relativistic particles one usually considers the limit $p^2 \ll m^2c^2a^2$ in which the equations simplify to

$$v = \frac{p}{ma}, \qquad (3.3)$$

$$p' = -ma\nabla\Psi. \qquad (3.4)$$

This simpler set of equations has the advantage that it is easy to find integration methods that are symplectic, i.e. that preserve the phase-space volume exactly as demanded by Hamiltonian time evolution. Note also that, in the absence of a gravitational force, the peculiar velocity scales exactly as $\propto a^{-1}$ in this case.

Even though the evolution of high-momentum particles suffers from severe errors, including the breakdown of causality for $p^2 > m^2c^2a^2$, some implementations might still use the simplified equations. The propagation of these errors into the clustering amplitude of matter is limited by the fact that high-momentum particles barely contribute to clustering in the first place. However, using the more accurate eqs. (3.1) and (3.2) is a more common choice. The issue of symplectic time integration in this case is discussed e.g. in appendix A of Adamek et al. [20] and appendix D of Elbers et al. [31].

## 3.2 Mesh-based methods

The main alternative to particle-based methods is to represent the distribution function on a spatial mesh. Yoshikawa et al. [28] discretise the distribution function on a six-dimensional mesh in phase-space, however, a brute-force approach like this is rather expensive and cannot easily be applied to very large simulations where memory requirements are a particular concern. A possible way out is to take moments of the distribution function (density, bulk velocity, and so on) where the momentum coordinates are integrated out so that a discretisation in the three spatial dimensions is sufficient.
Experience from solving the hierarchy of moments (the so-called Boltzmann hierarchy) in linear perturbation theory shows that a considerable number of moments must be taken into account in order to reach good accuracy for the evolution of the lowest moments. This concerns in particular the density that also affects the clustering of other matter components through gravitational coupling. On the other hand, the numerical solutions are readily available in the linear regime where they can be expressed in terms of linear transfer functions. Throughout this paper, we follow the convention from standard linear cosmological perturbation theory where the transfer function $T_X$ of any perturbation variable $X$ in Fourier space is defined through the relation

$$X(\mathbf{k}, z) = T_X(k, z)\zeta(\mathbf{k}), \qquad (3.5)$$

where $\zeta(\mathbf{k})$ is the comoving curvature perturbation of the mode $\mathbf{k}$ before it enters the horizon. Those transfer functions provide a deterministic factor by which any given initial random perturbation mode needs to be multiplied in order to obtain, for instance, the density perturbation at any given time. Using these transfer functions that can be calculated at the outset, a simulation code can therefore construct the linear density field of neutrinos at any point in time for any given realisation of the random initial conditions. This is precisely what basic mesh-based methods do: they use the density field of neutrinos extrapolated from linear perturbation theory, which is often a reasonable approximation because neutrinos do not cluster strongly. The method is free of shot noise and is by construction exact in the limit of linear perturbations. However, it obviously lacks any response to nonlinear gravitational potentials that develop in a simulation.

While the linear method produces results that are sufficiently accurate for many purposes, some more advanced approaches have been developed in attempts to address the shortcomings. Ali-Haïmoud & Bird [24] solve for the transfer function of the neutrino perturbations using the nonlinear matter power spectrum of the simulation to construct an effective source term, assuming that the phase correlation between neutrino and matter perturbations remain largely intact even at nonlinear scales. Dakin et al. [27] employ the coupled evolution equations for the lowest moments in their nonlinear form. Then, to avoid having to calculate a large Boltzmann hierarchy for every wavevector represented in the simulation, the hierarchy is truncated by assuming that a “scaling” holds approximately for ratios of higher moments, where the scaling coefficients are taken from linear theory. This approximation is then used to close the system of equations using only a small number of nonlinear moments. By construction this method agrees with the simpler method in the limit of linear perturbations. While the resulting nonlinear neutrino density is somewhat more realistic, the distribution function still has some residual errors that cannot easily be reduced without including further moments in the nonlinear computation.

## 3.3 Approximations and other methods

For some purposes, such as the computation of covariance matrices for different cosmological probes, it is useful to have methods for making cosmological predictions that are faster, although less accurate, than the traditional N-body methods discussed so far. Here we present two such methods that can be used as surrogates for N-body simulations: the COmoving Lagrangian Acceleration (COLA) approach and the PINpointing Orbit-Crossing Collapsed Hierarchical Objects (PINOCCHIO) approach. In both cases, a speed-up is achieved by drastically simplifying the time integration in the particle evolution. Finally, we also present a method
that avoids the need to include any neutrino physics in the actual *N*-body simulation altogether, apart from in the background solution. This method employs the so-called *Newtonian motion gauge* and can be used with virtually any numerical scheme that solves the Newtonian gravity problem.

### 3.3.1 COLA

The **COLA** approach by Tassev et al. [70] produces fast, approximate simulations of cosmological structure formation. Essentially, instead of solving for a full particle trajectory $\boldsymbol{x}(t)$, in this approach we solve for the deviations of the full trajectory about the trajectory predicted by second-order Lagrangian perturbation theory (2LPT) $\delta\boldsymbol{x}(t) = \boldsymbol{x}(t) - \boldsymbol{x}_{2\text{LPT}}(t)$. Since the evolution of the particles on large scales will be very close to that predicted by 2LPT, we can decrease the number of time steps of the simulation to trade accuracy at small scales for overall simulation speed while maintaining good accuracy at large scales. For a large number of time steps, the method effectively converges to a standard PM *N*-body method.

Adding massive neutrinos to the **COLA** method was described by Wright et al. [61], which also included an implementation in the **MG-PICOLA** simulation code<fnref n="1" /> by Winther et al. [60]. This implementation was carried over to the **COLA** solver within **FML**<fnref n="2" /> which succeeded **MG-PICOLA**. It is this implementation of the **COLA** solver within **FML** that we use in this paper. These implementations rely on the linear mesh-based method described above, i.e. we use the density field of neutrinos extrapolated from linear perturbation theory on a mesh for the PM part. For the 2LPT part of the **COLA** code, we make a further approximation to the 2LPT equation and use the $\Lambda$CDM kernel to speed up the computation.

To demonstrate the key advantage of **COLA** over traditional *N*-body codes, we use only 50 time steps linearly distributed in scale factor for the **COLA** simulations in this paper. However, the **COLA** method does not work well for simulations starting from high initial redshifts when using a relatively small number of time steps; for a discussion on how to optimise initial redshift and number of time steps in **COLA** simulations see sections 4.1 and 4.3 of Izard et al. [71]. Therefore, we use a slightly modified<fnref n="3" /> version of FML's built-in generator of initial conditions to generate initial particle data at $z = 19$ instead of $z = 127$ as is described in section 4.1 and used for the other methods in this paper. In addition, we use the CAMB Boltzmann solver by Lewis et al. [72, 73] to generate the density transfer function for massive neutrinos. Finally, we note that for all **COLA** simulations in this paper we use a force grid that is a factor of three finer than the mean inter-particle distance; for a thorough investigation of the impact of varying this factor in **COLA** simulations see section 4.4 of Izard et al. [71].

### 3.3.2 PINOCCHIO

The **PINOCCHIO** code<fnref n="4" /> [62, 74, 75] is an approximate method to generate halo catalogues in a very small fraction (of the order of 1/1000) of the time taken by an equivalent *N*-body simulation. Starting from a linear density field generated in Lagrangian space over a regular grid, its main goal is to construct halo catalogues by predicting which particles will end up in dark matter halos. To achieve this goal the algorithm first smoothes the linear density on a grid of

<footnote n="1">¹https://github.com/HAWinther/MG-PICOLA-PUBLIC.</footnote>

<footnote n="2">²https://github.com/HAWinther/FML/tree/master/FML/COLASolver.</footnote>

<footnote n="3">³The modifications are to the order in which (pseudo-)random numbers are drawn for the phases and amplitudes, such that we now do the same as the version of N-GenIC used for the other methods in this paper, see section 4.1.</footnote>

<footnote n="4">⁴https://github.com/pigimonaco/Pinocchio.</footnote>
smoothing radii, then uses ellipsoidal collapse to compute the collapse time of each particle. In the second step, it proceeds to group the collapsed particles into halos, using an algorithm that mimics their hierarchical clustering and distinguishes between halos and filaments.

As opposed to the other N-body methods employed in this work, **PINOCCHIO** does not integrate particle orbits but places halos at their final position using a single 2LPT or 3LPT displacement. Indeed, once the displacement fields are averaged over the multi-stream region that corresponds to a dark matter halo, LPT is very effective in predicting halo positions [76]. Another difference is that **PINOCCHIO** does not start from a set of displaced particles but generates the linear density field internally. Having implemented the same sequence for populating modes in k-space as the **N-GenIC** code that is used for the purpose of generating initial data in this work (see section 4.1), it can reproduce the same large-scale structure if the same seed for random numbers is provided.

The extension of **PINOCCHIO** to massive neutrinos is presented by Rizzo et al. [63] and is based on the result of Castorina et al. [18, 77] that the nonlinear clustering of massive neutrinos is negligible. We use **CAMB** to compute linear power spectra in massive neutrino cosmologies, and compute the scale-dependent growth rate of matter by taking ratios of power spectra of CDM and baryons (i.e. without neutrinos) at different times. We also adapt the code to incorporate a scale-dependent growth rate. With respect to the original implementation of Rizzo et al. [63], which was limited to 2LPT, we extend here the computation to third order: as shown by Munari et al. [75] this results in a significant improvement at mildly nonlinear scales.

Although the code has been conceived to predict the properties of dark matter halos, it can produce a full nonlinear density field as follows: particles that do not belong to halos are moved to their final position using 3LPT, halo particles are distributed around their halo center of mass following a Navarro-Frenk-White (NFW) profile [78] with Maxwellian velocity distributions. This allows us to construct snapshots like an N-body simulation, representing density fields that are far more accurate than a straight LPT implementation. Because we have only one type of particle, to compute the power spectrum of CDM and baryons needed below (section 5) we subtract the linear neutrino contribution from the total matter power spectrum obtained from the snapshot. To this end we also assume that the neutrino-matter cross-power spectrum, $P_{\nu,m}(k)$, can be approximated by $P_{\nu,m}(k) = \sqrt{P_{\nu}^{\text{L}}(k)P_{m}(k)}$, where the superscript “L” denotes a power spectrum from linear theory. This approximation is strictly only true in the linear regime but we apply it at all scales.

We do not expect **PINOCCHIO** to be competitive with N-body codes in predicting the matter power spectrum: taking it as a sophisticated implementation of 3LPT, we expect it to lose power on scales smaller than $k = 0.3 h \text{ Mpc}^{-1}$ for the halo power spectrum and $k = 0.2 h \text{ Mpc}^{-1}$ for the matter power spectrum. It will not be competitive with **COLA** as well which, being a PM code, can converge to the solution (on scales larger than the mesh) if a sufficient number of time steps is used. This better accuracy comes at a higher cost in computing time, by approximately a factor of eight in the configuration used in this paper (see also Blot et al. [79] for a similar benchmark), as well as in memory since our **COLA** runs use a grid three times finer than the mean particle separation. The **PINOCCHIO** code is widely used especially to characterise the covariance of galaxy clustering measurements, thanks to its low computational cost and its ability to generate halo catalogues.
### 3.3.3 Newtonian motion gauge

The Newtonian motion gauge approach for massive neutrinos was developed by Partmann et al. [32] and Heuschling et al. [49]. It allows for a simulation of nonlinear CDM in an ordinary Newtonian N-body simulation while accounting for the impact of linear neutrinos via a modification of the dark matter initial conditions and by employing a dynamically evolving coordinate system. The method is agnostic towards the implementation of the N-body simulation and for this paper we choose to employ GADGET-4 [48]. However, any method solving Newtonian nonlinear gravity is compatible, even methods other than N-body simulations. We would like to stress that our method is exact in the weak-field limit of general relativity (see Fidler et al. [80]) and therefore captures the full effect of linear neutrino perturbations on the nonlinear matter clustering.

The Newtonian motion approach allows for a very simple inclusion of massive neutrinos, requiring only three additional steps applied to a simulation without any neutrinos. First, we start from a set of “back-scaled” initial conditions based on the present-day power spectrum of CDM and baryons in the Newtonian motion coordinates, excluding neutrino perturbations. In contrast to other neutrino methods presented in this paper, these initial conditions do not assume a scale-dependent growth, i.e. the rescaling of the power spectrum is done using the scale-independent growth factor $D_+$. The residual effect of decaying modes due to neutrinos is included in the construction of the present-day matter power spectrum in the Newtonian motion coordinate system. This also has the added benefit that ordinary generators of initial conditions can be used without modifications, provided that the correctly back-scaled Newtonian motion gauge power spectrum is used. We then evolve the initial data with the Newtonian solver, taking into account the impact of the massive neutrinos on the background evolution via the Hubble rate. Finally, we obtain the output in the Newtonian motion gauge. To make it comparable to the output of other methods, we need to transform the result to the gauge employed therein (usually the “N-boisson” gauge [80]). This step is realised by a displacement field acting on the particle positions that is implemented in a similar way to how the initial conditions are set. The transformation accounts for the residual impact of neutrinos and other relativistic effects on the evolution of the CDM and baryon particles. However, by construction, the transformation vanishes exactly at $z = 0$ (or another chosen target redshift) while it is in general very small at late times, for small neutrino masses and on small scales. Therefore, it can often be neglected as it will only lead to small corrections in which case the output of the simulation can be used as-is.

For this work, we include only the first two steps, while omitting the final particle displacement. This leads to a small mismatch in the results shown for $z = 1$ at large scales for the cases of the highest neutrino masses. By leaving this correction out, we demonstrate that for most neutrino masses, box sizes, and redshifts the method is already sufficiently accurate in its simplest form. For more details on the transformation we refer the reader to the original work by Partmann et al. [32].

## 3.4 Halo-model reaction

The halo-model reaction approach provides the nonlinear corrections caused by massive neutrinos to a ΛCDM power spectrum through a ratio of halo-model predictions. Following Cataneo et al. [81], the nonlinear power spectrum is then given by

$$P^{\text{NL}}(k, z) = \mathcal{R}(k, z) P_{\text{pseudo}}^{\text{NL}}(k, z), \qquad (3.6)$$
where $\mathcal{R}(k, z)$ is the halo-model reaction and $P_{\text{pseudo}}^{\text{NL}}(k, z)$ is the nonlinear pseudo power spectrum.

The pseudo spectrum is a nonlinear $\Lambda$CDM power spectrum but with the initial conditions tuned such that its linear clustering exactly matches the linear clustering in the non-standard cosmology at the target redshift. For example, if the non-standard physics introduces a simple rescaling of the linear clustering amplitude, one could just rescale the amplitude of any $\Lambda$CDM power spectrum to produce the pseudo spectrum. In the case of scale-dependent modifications, this becomes a bit trickier in practice. We approximate this quantity as in previous works by Cataneo et al. [64, 81] and pass the modified linear spectrum as produced by CAMB to HMcode developed by Mead et al. [82], with $\Lambda$CDM presets, i.e. no baryonic feedback nor massive neutrinos. The benefit of using the pseudo rather than $\Lambda$CDM cosmology is that it guarantees the mass functions in target and pseudo cosmologies are similar as they have the same linear clustering. This produces a smoother transition between the two-halo and one-halo terms. This transition was one of the previous issues in calculating this nonlinear response using the halo model [83–85].

Following Cataneo et al. [64], the halo-model reaction for massive neutrinos is given by

$$ \mathcal{R}(k) = \frac{(1 - f_{\nu})^2 P_{\text{cb}}^{\text{HM}}(k) + 2f_{\nu}(1 - f_{\nu}) P_{\nu, \text{cb}}^{\text{HM}}(k) + f_{\nu}^2 P_{\nu}^{\text{L}}(k)}{P_{\text{m}}^{\text{L}}(k) + P_{\text{pseudo}}^{\text{1h}}(k)}, \quad (3.7) $$

with “cb” denoting the CDM and baryon component and “$\nu$” denoting massive neutrinos. We include the effects of massive neutrinos at the linear level in the numerator via the weighted sum of the nonlinear halo-model (cb) spectrum and the massive neutrino linear spectrum [15]. The components of the reaction are

$$ P_{\nu, \text{cb}}^{\text{HM}}(k) \approx \sqrt{P_{\text{cb}}^{\text{HM}}(k) P_{\nu}^{\text{L}}(k)}, \quad (3.8) $$

$$ P_{\text{cb}}^{\text{HM}}(k) = P_{\text{cb}}^{\text{L}}(k) + P_{\text{cb}}^{\text{1h}}(k). \quad (3.9) $$

Explicitly, the one-halo terms are given as integrals over the Fourier space halo density profile $u(k, M)$ and the halo mass function $n(M)$,

$$ P_{\text{cb}}^{\text{1h}}(k) = \int d\ln M n_{\text{cb}}(M) \left(\frac{M}{\bar{\rho}_{\text{cb}}}\right)^2 |u_{\text{cb}}(k, M)|^2, \quad (3.10) $$

$$ P_{\text{pseudo}}^{\text{1h}}(k) = \int d\ln M n_{\text{pseudo}}(M) \left(\frac{M}{\bar{\rho}_{\text{m}}}\right)^2 |u_{\text{pseudo}}(k, M)|^2, \quad (3.11) $$

where $\bar{\rho}$ is the background density for the specified matter species. The halo mass functions are given as

$$ n_{\text{cb}}(M) = \frac{\bar{\rho}_{\text{cb}}}{M} [\nu' f(\nu')] \frac{d \ln \nu'}{d \ln M}, \quad (3.12) $$

$$ n_{\text{pseudo}}(M) = \frac{\bar{\rho}_{\text{m}}}{M} [\nu'' f(\nu'')] \frac{d \ln \nu''}{d \ln M}. \quad (3.13) $$

The peak heights are defined as $\nu' = \delta_{\text{sc,cb}}(M)/\sigma_{\text{cb}}[R_{\text{cb}}(M)]$ and $\nu'' = \delta_{\text{sc,m}}/\sigma_{\text{m}}[R_{\text{m}}(M)]$, where the subscript “sc” indicates this quantity is calculated by solving the standard $\Lambda$CDM spherical-collapse equations but with the indicated matter density. The mass fluctuation
variances are given by

$$[\sigma_{\text{cb}}(R)]^2 = \int \frac{d^3 k}{(2\pi)^3} |\tilde{W}(kR)|^2 P_{\text{cb}}^{\text{L}}(k), \qquad (3.14)$$

$$[\sigma_{\text{m}}(R)]^2 = \int \frac{d^3 k}{(2\pi)^3} |\tilde{W}(kR)|^2 P_{\text{m}}^{\text{L}}(k). \qquad (3.15)$$

In all predictions from halo-model reaction, we employ a Sheth-Tormen halo mass function [86, 87], a power-law concentration-mass relation (see for example the work by Bullock et al. [88]), and an NFW halo density profile [78]. The predictions are computed numerically using the public code ReACT<fnref n="5" /> by Bose et al. [65, 89].

## 3.5 Power-spectrum emulation

Fast predictors of the matter power spectrum are an essential ingredient for many inference pipelines in cosmology. Since numerical simulations are too costly to be directly applied in this context, different approaches based on approximate methods or elaborate fitting techniques have been used in the past. Apart from the halo-model reaction discussed in the previous section, well-known examples are the *halofit* predictor developed by Smith et al. [90] and later improved by Takahashi et al. [91], and the *HMcode* predictor by Mead et al. [82], the latter being based on the halo model. In terms of the power-suppression signal of neutrinos, fitting routines by Bird et al. [16] have been used in the past.

More recently, the emulation technique has become a popular alternative to obtain fast predictions of the matter power spectrum within the cosmological parameter space. Broadly speaking, emulators are interpolation routines based on a suite of numerical simulations that sample the cosmological parameter space and act as a training set. There are different surrogate techniques currently used for cosmological emulators, such as Gaussian process regression [66, 92–94], polynomial chaos expansion [7, 8], or neural network approaches [95].

In addition to the predictions from *halofit* and *HMcode* mentioned above, we focus in this paper on the *CosmicEmu* [66], the *EuclidEmulator2* [8], and the *BACCOemulator* [45]. These emulators provide predictions of the matter power spectrum and include a free parameter for the sum of the neutrino masses. We will now summarise the particularities of these three emulators, specifically focusing on the neutrino implementation.

* The *CosmicEmu-2022* is built upon the *Mira-Titan* simulations [67, 96], a suite of 111 simulations run with the *HACC* code [97]. The simulations are distributed over an eight-dimensional cosmological parameter space comprising ($\omega_m, \omega_b, \omega_\nu, \sigma_8, h, n_s, w_0, w_a$), where $\sigma_8$ is the present-day amplitude of linear matter density fluctuations at the scale of $8h^{-1}$ Mpc, $n_s$ is the scalar spectral index, and $w_0, w_a$ parameterise the effective equation of state of dark energy in terms of the first two coefficients of the Taylor series expansion around $a = 1$. The emulator achieves an absolute precision of about four percent for modes of $k < 5h^{-1}$ Mpc$^{-1}$ within the redshift range $z \in [0, 2]$. Neutrinos are not incorporated in the simulations and are effectively treated as a smooth background component. The power spectra are then corrected on large scales for enhanced growth beyond the neutrino free-streaming scale using the scale-dependent linear growth factor.

* The *EuclidEmulator2* is trained on 200 paired and fixed simulations that were run with *PKDGRAV3*. It emulates the nonlinear boost factor that is then multiplied by the results

<footnote n="5">⁵https://github.com/nebblu/ReACT.</footnote>
of a linear Boltzmann solver. The emulator covers eight cosmological parameters ($\Omega_m$, $\Omega_b$, $\sum m_\nu$, $A_s$, $h$, $n_s$, $w_0$, $w_a$), where $A_s$ is the amplitude of primordial perturbations at the scale $k_p = 0.05$ Mpc$^{-1}$, and includes redshifts of $z \in [0, 3]$ and modes up to $k = 10$ h Mpc$^{-1}$. It claims an error of below one percent which is better than the other emulators discussed here. Note, however, that the EuclidEmulator2 covers a somewhat narrower parameter space motivated by the results of the Planck mission [98]. Within the training set the neutrinos are modelled using the mesh-based method implemented in PKDGRAV3.

* The BACCOemulator is trained on a very large suite of simulations based on the cosmology-rescaling technique [99]. More specifically, four high-resolution simulations with judiciously chosen cosmologies are rescaled to more than 800 cosmologies at different redshifts. Whenever the target cosmologies included massive neutrinos, their effect is added following the extension of the cosmology-rescaling technique presented by Zennaro et al. [33]. This emulator varies eight cosmological parameters ($\Omega_{cb}$, $\Omega_b$, $\sum m_\nu$, $\sigma_8$, $h$, $n_s$, $w_0$, $w_a$) and covers a redshift range of $z \in [0, 1.5]$ for modes up to $k = 5$ h Mpc$^{-1}$. The claimed precision is better than three percent.

We refer to the original references for more information about the emulators.

# 4 Simulations

To compare different numerical methods, we carry out a large suite of N-body simulations where we employ different codes to run the same set of ten simulations summarised in table 2. These ten simulations cover different choices of total neutrino mass $\sum m_\nu$ (including the massless case), different box sizes $L_{\text{box}}$, and different mass resolutions to check for numerical convergence with respect to finite-volume and discretisation effects. $N_{\text{part}}$ denotes the number of particles used for CDM and baryons, as well as the number of particles for neutrinos if a particle-based method is employed. For simplicity, we assume degenerate neutrino mass eigenstates because cosmology is mainly sensitive to the total neutrino mass scale [100]. We keep the total matter density (at redshift zero) fixed at $\Omega_m = 0.319$ by adjusting the CDM density parameter together with the neutrino mass. The baryon density is fixed at $\Omega_b = 0.049$, and the remaining cosmological parameters are $A_s = 2.215 \times 10^{-9}$ at the pivot scale $k_p = 0.05$ Mpc$^{-1}$, $n_s = 0.9619$, and $h = 0.67$, which are based on the Euclid Flagship 2 simulation. Dark energy is modelled as a cosmological constant that provides a spatially flat background solution, and the CMB temperature is set to 2.7255 K and the effect of radiation is taken into account in the simulations at the linear level, either by using carefully tailored initial conditions as detailed below or by including the radiation component on the mesh if a mesh-based method is used.

Using identical initial data (see section 4.1 below for details) in each case, the ten simulations are run with each of the thirteen N-body methods listed in table 1 to produce particle snapshots at redshifts $z = 1$ and $z = 0$. For AREPO, due to resource constraints, only the four simulations with $N_{\text{part}} = 512^3$ are run, precluding the possibility of conducting numerical convergence tests in this case. Therefore, a total of 248 individual snapshots are analysed in this code comparison. In addition, nonlinear power spectra are predicted for each distinct choice of neutrino masses using the remaining methods in table 1, as well as using the HMcode and halofit fitting methods.
<table><caption>**Table 2.** Overview of the basic parameters used in our simulation suite. The cases with $\sum m_\nu = 0$ eV and $\sum m_\nu = 0.15$ eV are the two main baselines for our comparison, but we include some cases with larger masses, up to $\sum m_\nu = 0.6$ eV, to probe more “extreme” regions of parameter space.</caption><thead><tr><th>Simulation</th><th>L<sub>box</sub></th><th>N<sub>part</sub></th><th>mass resolution</th><th>Σ m<sub>ν</sub></th></tr></thead><tbody><tr><td>0.0eV</td><td>512 h<sup>-1</sup> Mpc</td><td>512<sup>3</sup></td><td>8.85×10<sup>10</sup> h<sup>-1</sup> M<sub>☉</sub></td><td>0.0 eV</td></tr><tr><td>0.15eV</td><td>512 h<sup>-1</sup> Mpc</td><td>512<sup>3</sup></td><td>8.75×10<sup>10</sup> h<sup>-1</sup> M<sub>☉</sub></td><td>0.15 eV</td></tr><tr><td>0.3eV</td><td>512 h<sup>-1</sup> Mpc</td><td>512<sup>3</sup></td><td>8.65×10<sup>10</sup> h<sup>-1</sup> M<sub>☉</sub></td><td>0.3 eV</td></tr><tr><td>0.6eV</td><td>512 h<sup>-1</sup> Mpc</td><td>512<sup>3</sup></td><td>8.45×10<sup>10</sup> h<sup>-1</sup> M<sub>☉</sub></td><td>0.6 eV</td></tr><tr><td>0.0eV_HR</td><td>512 h<sup>-1</sup> Mpc</td><td>1024<sup>3</sup></td><td>1.11×10<sup>10</sup> h<sup>-1</sup> M<sub>☉</sub></td><td>0.0 eV</td></tr><tr><td>0.15eV_HR</td><td>512 h<sup>-1</sup> Mpc</td><td>1024<sup>3</sup></td><td>1.09×10<sup>10</sup> h<sup>-1</sup> M<sub>☉</sub></td><td>0.15 eV</td></tr><tr><td>0.0eV_1024Mpc</td><td>1024 h<sup>-1</sup> Mpc</td><td>1024<sup>3</sup></td><td>8.85×10<sup>10</sup> h<sup>-1</sup> M<sub>☉</sub></td><td>0.0 eV</td></tr><tr><td>0.15eV_1024Mpc</td><td>1024 h<sup>-1</sup> Mpc</td><td>1024<sup>3</sup></td><td>8.75×10<sup>10</sup> h<sup>-1</sup> M<sub>☉</sub></td><td>0.15 eV</td></tr><tr><td>0.3eV_1024Mpc</td><td>1024 h<sup>-1</sup> Mpc</td><td>1024<sup>3</sup></td><td>8.65×10<sup>10</sup> h<sup>-1</sup> M<sub>☉</sub></td><td>0.3 eV</td></tr><tr><td>0.6eV_1024Mpc</td><td>1024 h<sup>-1</sup> Mpc</td><td>1024<sup>3</sup></td><td>8.45×10<sup>10</sup> h<sup>-1</sup> M<sub>☉</sub></td><td>0.6 eV</td></tr></tbody></table>

## 4.1 Initial conditions

The initial conditions of all simulations are generated at redshift $z = 127$.<fnref n="6" /> The linear matter power spectra and transfer functions are obtained by running either CAMB or the CLASS Boltzmann code by Blas et al. [101]. These files are then used by the REPS<fnref n="7" /> code to compute the rescaled power spectra and transfer functions at $z = 127$ by solving the multi-fluid linear equations as outlined by Zennaro et al. [102]. This procedure, known as rescaling, guarantees that the power spectrum of the output of the simulation on linear scales at low redshift will match the correct linear power spectra. A realisation of initial data is then generated by drawing random phases for all perturbation modes and fixing their amplitudes according to the initial transfer functions. This approach of “fixing” the amplitudes effectively removes cosmic variance at linear scales and has been shown to generally produce less noisy summary statistics [103]. It introduces a specific type of non-Gaussianity that is not expected to affect any of our results. Given the density field of CDM and baryons, the initial positions and velocities of the N-body particles are computed from the Zeldovich approximation [104]. We employ a modified version of the N-GenIC code<fnref n="8" /> that accounts for the scale-dependence present in both the growth rate and growth factor in cosmologies with massive neutrinos.

For neutrinos, the different implementations make use of distinct methods. For particle-based implementations, the positions and velocities of the massive neutrino particles are generated in a similar fashion as CDM. This means we effectively make use of the Zeldovich approximation to set the first two moments of the phase-space distribution function which correspond to the density perturbation and bulk velocity, respectively. The initial velocities of neutrino particles are then offset by a random thermal component drawn from their Fermi-Dirac distribution [105, 106]. We assume the standard Big Bang scenario where this distribution is set in equilibrium before the weak interaction freezes out — when the Universe was about 1 second old — and after freeze-out simply redshifts as the Universe expands. Note that this means that the typical thermal velocities are much larger than those of a dis-

<footnote n="6">⁶COLA exceptionally uses $z = 19$ as initial redshift (see section 3.3.1).</footnote>

<footnote n="7">⁷https://github.com/matteozennaro/reps.</footnote>

<footnote n="8">⁸https://github.com/franciscovillaescusa/N-GenIC_growth.</footnote>
tribution that is in thermal equilibrium at low redshift. For SWIFT, neutrino particles are instead set up with the FastDF code,<fnref n="9" /> using geodesic integration from high redshift [20, 69], to reproduce the full distribution function and to prevent the initial perturbations from being erased by thermal motions [55]. For mesh-based implementations, on the other hand, the density field of each neutrino species is directly computed using the phases from the random field realisation of the linear initial conditions.

## 4.2 Post-processing pipeline

In order to quantify differences in our numerical schemes as precisely as possible, we analyse the snapshots of all our N-body simulations in a common pipeline. We compute the power spectra and bispectra of the CDM and baryon component, and produce halo catalogues from which we measure the halo mass functions and halo bias. In cases where the simulations provide a neutrino distribution, we also compute the cross-power spectra of neutrinos with the CDM and baryon component, as well as the neutrino auto-power spectra.

### 4.2.1 Power spectra

The power spectra of the different snapshots have been estimated using Pylians3.<fnref n="10" /> The routine first deposits particle masses into a regular 3D grid with $N^3$ voxels using the cloud-in-cell mass-assignment scheme. In this work, we always use a mesh with $N^3 = N_{\text{part}}$ such that the Nyquist scales match between particles and mesh. Although using larger grids may improve measurements on smaller scales, we recommend caution due to potential systematics and advise against relying on results near or beyond the Nyquist scale set by the mean particle separation. The constructed field is then Fourier transformed and the effects of the mass-assignment scheme are corrected. Next, for each mode the square of its amplitude is computed, $|\delta(\mathbf{k})|^2$. The modes are then binned in intervals of width equal to the fundamental frequency $k_{\text{F}} = 2\pi L_{\text{box}}^{-1}$ and the power spectrum is finally estimated as

$$P(k_i) = \frac{1}{N_i} \sum_{\mathbf{k} \in B_i} |\delta(\mathbf{k})|^2, \qquad (4.1)$$

where $N_i$ is the number of independent modes in the considered bin $B_i = \{\mathbf{k} | i k_{\text{F}} \le |\mathbf{k}| < (i+1)k_{\text{F}}\}$ and $k_i$ is computed as

$$k_i = \frac{1}{N_i} \sum_{\mathbf{k} \in B_i} |\mathbf{k}|. \qquad (4.2)$$

To compute the cross-power spectrum of two fields instead, the estimator is generalised in the most straightforward way,

$$P_{X,Y}(k_i) = \frac{1}{2N_i} \sum_{\mathbf{k} \in B_i} \left[ \delta_X(\mathbf{k}) \delta_Y^*(\mathbf{k}) + \delta_X^*(\mathbf{k}) \delta_Y(\mathbf{k}) \right], \qquad (4.3)$$

where the subscripts “X” and “Y” denote the two different fields. When presenting the results, we combine the measurements into larger bins logarithmically spaced in $k$. This reduces the noise at large $k$ and makes our plots more readable.

Some codes do not produce snapshots at exactly the desired redshift, but at redshifts that deviate by less than ± 0.01 from the target redshift. These differences can be visible when

<footnote n="9">⁹https://github.com/wullm/fastdf.</footnote>

<footnote n="10">¹⁰https://pylians3.readthedocs.io.</footnote>
comparing power spectra on a sub-percent level. For those cases, we rescale the power spectra by the square of the ratio of the linear growth factors at the respective redshift values. Such a rescaling is applied to ANUBIS for all snapshots and to L-GADGET3 and PKDGRAV3 at $z = 1$.

### 4.2.2 Bispectra

We measure the bispectrum of CDM and baryons (ccc) using the estimator

$$ \hat{B}_{\text{ccc}}(k_l, k_m, k_n) = \frac{k_F^3}{N_{\text{tr}}(k_l, k_m, k_n)} \sum_{q_1 \in B_l} \sum_{q_2 \in B_m} \sum_{q_3 \in B_n} \delta_K(\mathbf{q}_1 + \mathbf{q}_2 + \mathbf{q}_3) \delta(\mathbf{q}_1) \delta(\mathbf{q}_2) \delta(\mathbf{q}_3), \quad (4.4) $$

where $N_{\text{tr}}$ is the number of “fundamental triangles”,

$$ N_{\text{tr}}(k_l, k_m, k_n) = \sum_{\mathbf{q}_1 \in B_l} \sum_{\mathbf{q}_2 \in B_m} \sum_{\mathbf{q}_3 \in B_n} \delta_K(\mathbf{q}_1 + \mathbf{q}_2 + \mathbf{q}_3), \quad (4.5) $$

formed by the vectors $\mathbf{q}_i$ satisfying the triangle condition $\mathbf{q}_1 + \mathbf{q}_2 + \mathbf{q}_3 = 0$ that are included within the “triangle bin” defined by the triplet of centers $(k_l, k_m, k_n)$ and corresponding bins $B_l, B_m, B_n$.

We use a Python code implementing the fourth-order density interpolation and the interlacing scheme described by Sefusatti et al. [107]. In order to compare the large-volume simulations ($L_{\text{box}} = 1024 h^{-1} \text{Mpc}$) more easily with the small-volume ones, we use the same $k$-space binning in both cases, fixing the bin width to $k_F$ of the small box. Just like for the power spectra, to account for inaccuracies in the redshift of some snapshots, we rescale some of the resulting bispectra by the cube of the ratio of the linear growth factors at the respective redshift values. Such a rescaling is applied to ANUBIS for all snapshots and to L-GADGET3 and PKDGRAV3 at $z = 1$.

### 4.2.3 Halo catalogues

For the considered snapshots of the various simulations, we identify halos with the code Denhf [108–111] which uses a “spherical overdensity” criterion. The algorithm does not rely on any pre-identification method. Only CDM and baryon particles are considered in the characterisation of halos; neutrino particles (if present at all in the simulation) are considered as a background component [17, 77].

Denhf estimates the local density at the position of each $N$-body particle by calculating the distance to its $10^{\text{th}}$-nearest neighbour $d_{10}$, and assigning to each particle a density that is proportional to $d_{10}^{-3}$. Centered on the particle with the highest density value, the algorithm grows a sphere and stops when the mean density within the sphere falls below a desired overdensity threshold, set to 200 times the background density of CDM and baryons for the purpose of this work. All particles assigned to this spherical overdensity halo are then removed from the global list of particles, and the algorithm proceeds recursively until none of the remaining particles has a local density large enough to be the center of a 10-particle halo. Particles not assigned to halos will be part of the field.

Only in the case of PINOCCHIO we use the halo catalogues as produced by the code itself instead of Denhf. Because PINOCCHIO is calibrated on the friends-of-friends halo mass function, we translate its masses to spherical overdensity ones by applying the rescaling of halo masses that translates the halo mass function of Watson et al. [112], that has been used to calibrate the code, to the one of Tinker et al. [113]. Such a rescaling has been used, e.g., by Fumagalli et al. [10] to force the halo mass function averaged over 1000 realisations to
<figure><img src="image_4.png" /><figcaption>Figure 1. Power spectrum of CDM and baryons as measured from different codes, relative to GADGET-3 at $z=0$ and $z=1$ for a neutrino mass of $\sum m_\nu = 0.15$eV in the simulation with $L_{\text{box}} = 512 h^{-1} \text{ Mpc}$ and $N_{\text{part}} = 512^3$ CDM and baryon particles. The corresponding particle Nyquist wavenumber is indicated by the grey dash-dotted line. The grey bands highlight the interval of $\pm 0.01$.</figcaption></figure>

follow a target one. We compute the rescaling only once, in the case of massless neutrinos, and use it for all neutrino masses.

For the codes that do not produce snapshots at the exact values of the desired redshifts, we apply no further corrections here. The error in the redshift is less than $\pm 0.01$ while our halo properties typically display disagreements larger than 1% between different codes. We therefore assume that the error due to mismatching redshift values is subdominant.

# 5 Results

## 5.1 Power spectra

A key prediction from neutrino simulations is a suppression of the matter power spectrum that exceeds the maximum linear theory prediction of $\Delta P_m/P_m \approx 8f_\nu$. The matter power spectrum can be decomposed as follows

$$P_m(k) = (1 - f_\nu)^2 P_{cb}(k) + 2f_\nu(1 - f_\nu)P_{\nu,cb}(k) + f_\nu^2 P_\nu(k), \quad (5.1)$$

where $P_{cb}$ is the power spectrum of CDM and baryons, $P_{\nu,cb}$ the cross-power spectrum of neutrinos with CDM and baryons, and $P_\nu$ the neutrino power spectrum. Various methods treat these components differently or make predictions for only some of them, so we discuss each component in turn. Finally, we will also compare the total $P_m(k)$ with various power spectrum emulators, including EuclidEmulator2 which only predicts this quantity.

### 5.1.1 CDM and baryons

The leading contribution to $P_m(k)$ is $P_{cb}(k)$, which is suppressed in massive neutrino models. Figure 1 shows the ratio of $P_{cb}(k)$ for models with a neutrino mass of $\sum m_\nu = 0.15$ eV relative
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Suppression of the power spectrum of CDM and baryons at $z=0$ and $z=1$ for three different neutrino masses, $\sum m_\nu \in \{0.15, 0.3, 0.6\}$ eV, when compared to the massless case. Results are from the simulations with $L_{\text{box}} = 512 \, h^{-1} \text{Mpc}$ and $N_{\text{part}} = 512^3$. The corresponding particle Nyquist wavenumber is indicated by the grey dash-dotted line.</figcaption></figure>

to the GADGET-3 simulation which we arbitrarily pick as the reference. In all our figures, results from codes where massive neutrinos are represented through an N-body ensemble are plotted using solid lines, while other N-body methods, including surrogates, use dashed lines. Any additional predictions use dotted lines. In figure 1, the linear prediction, computed with CLASS, is shown by the rose dotted line that drops off sharply at $k \approx 0.1 \, h \, \text{Mpc}^{-1}$ beyond which the error quickly exceeds 10%. At $z=0$, on the largest scales, all codes deviate less than 1% — the fluctuations seen in the dotted lines are largely due to the lack of cosmic variance in the codes that predict $P_{\text{cb}}(k)$ directly. On smaller scales, some of the codes start to deviate from the GADGET-3 reference. PINOCCHIO is in agreement within 1% up to $k \approx 0.1 \, h \, \text{Mpc}^{-1}$, ANUBIS up to $k \approx 0.3 \, h \, \text{Mpc}^{-1}$, gevolution up to $k \approx 0.5 \, h \, \text{Mpc}^{-1}$, and COLA up to $k \approx 0.7 \, h \, \text{Mpc}^{-1}$. The other codes stay within a 1% deviation from GADGET-3 for all scales down to the particle Nyquist scale of $k_{\text{Nyq}} = \pi N_{\text{part}}^{1/3} L_{\text{box}}^{-1} \approx 3 \, h \, \text{Mpc}^{-1}$. This scale is indicated by a vertical dash-dotted line. The ReACT and BACCOemulator codes stay accurate within 1% to 5% on all scales down to the Nyquist scale. At $z=1$ the qualitative behaviour is similar, but most notably the codes disagree more on large scales while still staying within a 1% agreement. Here, PINOCCHIO is accurate up to $k \approx 0.2 \, h \, \text{Mpc}^{-1}$.

Figure 2 shows the ratio of $P_{\text{cb}}(k)$ for models with massive neutrinos relative to the massless case for $\sum m_\nu \in \{0.15, 0.3, 0.6\}$ eV. For later convenience, we denote this suppression ratio by $S_{\text{cb}}(k)$, or in general

$$S_X(k) = \frac{P_X^{(\text{massive})}(k)}{P_X^{(\text{massless})}(k)}, \qquad (5.2)$$

where the subscript “X” denotes any component in question, here $X \to \text{cb}$. Linear calculations (taken from CLASS) predict that, on the largest scales, the growth of structure is mostly
<figure><img src="image_6.png" /><figcaption>**Figure 3.** The suppression of the power spectrum of CDM and baryons relative to the one measured in the GADGET-3 reference runs at $z = 0$ and $z = 1$ for $\sum m_\nu = 0.15$ eV for the simulations with $L_{\text{box}} = 512 \, h^{-1} \text{Mpc}$ and $N_{\text{part}} = 512^3$. The corresponding particle Nyquist wavenumber is indicated by the grey dash-dotted line. The grey bands highlight the interval of $\pm 0.01$.</figcaption></figure>

unaffected so that $S_{\text{cb}}(k)$ approaches unity, while on small scales $S_{\text{cb}}(k)$ reaches a plateau. Compared to this linear expectation, all codes in the comparison predict slightly less suppression around $k = 0.1 \, h \, \text{Mpc}^{-1}$ and a much greater suppression for $k > 0.3 \, h \, \text{Mpc}^{-1}$, followed by an upturn on nonlinear scales. This upturn has been repeatedly demonstrated and results from the reduced sensitivity of the one-halo contribution [41, 114]. At $z = 0$, we obtain excellent agreement between all simulations and most approximate methods up to the scale of maximum suppression at $k_{\text{max}} \approx 1 \, h \, \text{Mpc}^{-1}$, where the suppression is 20% greater than the linear prediction. At $z = 1$, the scale of maximum suppression shifts to $k_{\text{max}} \approx 2 \, h \, \text{Mpc}^{-1}$ and the differences are greater, both with the linear prediction and between the codes, with the exception of PINOCCHIO$^{11}$ which fares significantly better compared to $z = 0$.

To study these relative differences in greater detail, we show $S_{\text{cb}}(k)$ for the smallest neutrino mass, $\sum m_\nu = 0.15$ eV, relative to the GADGET-3 prediction for this quantity in figure 3. With the exception of PINOCCHIO and CLASS, all codes agree to better than 1% at $z = 0$ up to the particle Nyquist scale. Near $k_{\text{Nyq}}$, the approximate COLA method and the BACCOemulator differ by more than 1% from the bulk of the simulations, while gevolution differs by slightly less than 1%. Beyond this scale, the predictions diverge and should be compared to higher-resolution runs since our estimator of the power spectrum is computed on a mesh with a matching Nyquist scale. The measured power spectra therefore cannot be used beyond $k_{\text{Nyq}}$ where they become strongly biased. This particularly affects the comparison between simulations — where the power spectrum estimator is employed — and other methods to predict $P_{\text{cb}}(k)$. At $z = 1$, nonlinearities are smaller and the agreement between the simulations is better. Here the snapshot produced by PINOCCHIO achieves percent accuracy to $k \approx 0.3 \, h \, \text{Mpc}^{-1}$. Some of the approximate methods fare slightly worse at this earlier time

<footnote>$^{11}$The reason for this disagreement is due to the fact that the accuracy of LPT-based PINOCCHIO depends on the level of nonlinearity that varies with neutrino mass.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** The suppression of the power spectrum of CDM and baryons relative to the one measured in the GADGET-3 reference runs at $z = 0$ for $\sum m_\nu \in \{0.3, 0.6\}$ eV for the simulations with $L_{\text{box}} = 512 \, h^{-1} \text{ Mpc}$ and $N_{\text{part}} = 512^3$. The corresponding particle Nyquist wavenumber is indicated by the grey dash-dotted line. The grey bands highlight the interval of $\pm 0.01$.</figcaption></figure>

compared to $z = 0$, with the difference between CLASS and the simulations increasing by 50% and COLA diverging beyond $k = 2 \, h \, \text{Mpc}^{-1}$. The NM-GADGET4 method requires an additional post-processing coordinate transformation at any redshift except $z = 0$. Because this additional step is omitted in this work for simplicity, a small error at low wavenumbers remains in the $z = 1$ data. This explains why the error is larger at that redshift than at the final time.

ANUBIS notably drops off for $k > 1 \, h \, \text{Mpc}^{-1}$ at $z = 1$ and also shows a small excess on linear scales both at $z = 0$ and $z = 1$. This excess originates from the massless case and is a result of ANUBIS being run with a coarser base grid than the other codes ($512^3$ for all simulations) due to limited resources. A finer base grid requires more memory but the lack of it can be somewhat compensated by using a smaller time step. For the ANUBIS massive neutrino runs, this is done automatically but for the massless case the time step has been set to half of that originally calculated by the code. Tests indicate that further reducing the time step or ideally using a finer base grid should lessen the excess at large scales, but finding the optimal choice of code settings is not the aim of this work. The drop-off observed for ANUBIS at $z = 1$ for $k > 1 \, h \, \text{Mpc}^{-1}$ is due to differences in resolution between the various simulations. As ANUBIS is an AMR-code, a modified version of the RAMSES code originally written by Teyssier [56], the inclusion of massive neutrinos, which suppresses clustering on small scales, also reduces the number of refinements reached in the simulation compared to the massless case. This effect increases with the neutrino mass and is more noticeable for higher redshifts where there is also less refinement due to less clustering. This can be solved by a higher particle density which automatically leads to more refinement. Generally this is also necessary to find a better agreement between ANUBIS and GADGET-3 for smaller scales, as can be observed from the high-resolution runs and also in a previous comparison conducted by Schneider et al. [115] between RAMSES, GADGET-3, and PKDGRAV3.

The observed differences are greater when the neutrino mass is increased, especially for the approximate methods. Figure 4 shows $S_{\text{cb}}(k)$ relative to the GADGET-3 prediction for
models with $\sum m_\nu = 0.3$ eV (left panel) and $\sum m_\nu = 0.6$ eV (right panel) at $z=0$. While the **BACCO emulator** agrees to better than 1% with **GADGET-3** up to $k_{Nyq}$ for 0.3 eV, it makes no prediction for 0.6 eV because that value lies far outside the range covered by the training set of the emulator. We deliberately include a case of such a large neutrino mass to exacerbate the differences between the various numerical implementations. The reaction method differs by slightly more than 1% and 2% at $k = 0.7$ h Mpc$^{-1}$ for total neutrino masses of 0.3 eV and 0.6 eV, respectively. The differences with the linear prediction (**CLASS**) and with **PINOCCHIO** are large, as noted above for 0.15 eV.

The agreement is excellent for the other codes, but some subtle differences can be discerned. On large scales, we observe that **SWIFT**, **CONCEPT**, **COLA**, **gevolution**, **PKDGRAV3**, **NM-GADGET4**, and **L-GADGET3** show the same coherent scatter relative to **GADGET-3**, especially for a neutrino mass of 0.6 eV. This is due to the contamination from shot noise in the neutrino particle implementation used by the **GADGET-3** run. The other mentioned codes have implementations that do not suffer from shot noise or take measures to limit the contamination. For instance, in **gevolution** the neutrino N-body ensemble is evolved throughout the simulation, but it is only used as source of gravitational fields from redshift $z=7$ and below. At higher redshifts, the code uses the linear grid-based density instead. The reasoning behind this strategy is that shot noise is constant over time and hence more problematic at high redshift where cosmological perturbations are smaller in comparison. On the other hand, the linear prediction for neutrinos is expected to be very accurate at high redshift. One can therefore reduce the total error by judiciously choosing the time at which the code switches from linear to fully nonlinear neutrino treatment. It is nonetheless reassuring that even without mitigating against shot noise the scatter remains far below 1%.

On small scales the differences are also larger. The lines for **COLA** and **gevolution** track each other closely, but differ by more than 1% from the other codes for $k > 1$ h Mpc$^{-1}$. **SWIFT**, **PKDGRAV3**, and **CONCEPT** are low compared to **GADGET-3** for $k > 1$ h Mpc$^{-1}$, unlike what was seen in figure 3. **ANUBIS** diverges from **GADGET-3** by more than 1% for $k > 1$ h Mpc$^{-1}$ for the 0.6 eV neutrino mass case. As mentioned earlier, this is due to the fact that the AMR scheme has a lower effective resolution as the neutrino mass increases, simply because refinement is triggered by clustering. Finally, **NM-GADGET4** is slightly higher than **GADGET-3**. However, these differences remain below 1% well beyond $k_{Nyq}$.

### 5.1.2 Convergence tests

To study the numerical convergence of our results, we consider the effects of finite box size and resolution. Figure 5 shows the relative suppression for the runs with larger volume (left panel) and higher resolution (right panel). When $L_{box}$ is doubled at fixed resolution, the agreement remains excellent on linear scales and is sometimes even slightly better around $k_{Nyq}$, providing an important consistency check for most codes. Increasing instead the mass resolution by doubling $k_{Nyq}$, the agreement between the simulations improves significantly on nonlinear scales. Including more scales in either direction, most codes remain within 1% of the reference runs done with **GADGET-3**. The excess on large scales for the case of **ANUBIS** persists as a result of the coarse base grid. For the runs with $N_{part} = 1024^3$, this base grid is even less suited and the time-steppings for the massless cases are set to 0.1 and 0.15 times the original time step for the larger-box and high-resolution runs, respectively.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** The suppression of the power spectrum of CDM and baryons relative to the one measured in the GADGET-3 reference runs at $z=0$ for $\sum m_\nu = 0.15$ eV for the simulations with a larger volume, $L_{\text{box}} = 1024 \, h^{-1} \text{Mpc}$ and $N_{\text{part}} = 1024^3$ (left panel), and at a higher resolution in the small volume, $L_{\text{box}} = 512 \, h^{-1} \text{Mpc}$ and $N_{\text{part}} = 1024^3$ (right panel). The respective particle Nyquist wavenumbers are indicated by the grey dash-dotted lines. The grey bands highlight the interval of $\pm 0.01$.</figcaption></figure>

### 5.1.3 Contributions from neutrinos

The subdominant contributions to $P_m(k)$ are the cross-power spectrum $P_{\nu,cb}(k)$ between neutrinos and the CDM and baryon component, and the auto-power spectrum $P_\nu(k)$ of neutrinos. As can be seen from eq. (5.1), these are suppressed by the small factors $f_\nu$ and $f_\nu^2$, respectively, and they are themselves additionally strongly suppressed with respect to $P_{cb}(k)$ on scales smaller than the neutrino free-streaming scale. While both of these contributions will be exceedingly hard to constrain individually from observations, it is nonetheless interesting to study them in the context of our code comparison in order to highlight some more subtle differences in the numerical schemes. Figure 6 (left panel) shows $P_{\nu,cb}(k)$ for various codes relative to the result from the SWIFT code, for the smallest neutrino mass, $\sum m_\nu = 0.15$ eV, computed from the high-resolution simulations. We use SWIFT as the reference here because it has a very low level of shot noise in the neutrino component, yet is able to track the nonlinear evolution of neutrinos very accurately. As is the case with CDM and baryons, linear theory cannot describe neutrino clustering on nonlinear scales and therefore CLASS significantly underestimates the cross-power spectrum for $k > 0.1 \, h \, \text{Mpc}^{-1}$. Most other codes use a particle implementation of neutrinos and scatter about the SWIFT prediction partially due to shot noise. The results from CONCEPT have no shot noise, but depart from the other codes for $k > 0.2 \, h \, \text{Mpc}^{-1}$. The relative difference to SWIFT is 8% at $k = 1 \, h \, \text{Mpc}^{-1}$. Much the same applies to the neutrino auto-power spectrum, $P_\nu(k)$, shown relative to the neutrino auto-power spectrum of SWIFT in the right panel of figure 6. Here, for the codes with a neutrino particle ensemble, the dominant contribution to the shot noise is removed. Except for SWIFT this is done by subtracting the inverse of the average neutrino particle density, $\bar{n}^{-1}$, from the measured power spectrum, where $\bar{n}^{-1} = 0.125 \, h^{-3} \, \text{Mpc}^3$ in our high-resolution simulations. For SWIFT the subtracted values are derived from the $\delta f$ method. The difference
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Relative cross-power spectrum of neutrinos with CDM and baryons (left) and neutrino auto-power spectrum (right) with respect to SWIFT at $z=0$ for $\sum m_\nu = 0.15$ eV for the higher-resolution simulation with $L_{\text{box}} = 512 \, h^{-1} \text{Mpc}$ and $N_{\text{part}} = 1024^3$. The Nyquist wavenumber is indicated by the grey dash-dotted line. The grey bands highlight the interval of $\pm 0.01$.</figcaption></figure>

between SWIFT and CONCEPT is 19% (off the chart) at $k = 1 \, h \, \text{Mpc}^{-1}$ and the effects of shot noise are even more evident in the other codes. The agreement between gevolution and SWIFT is quite remarkable, and on the largest scales these results are also more consistent with CONCEPT than with the other N-body codes.

### 5.1.4 Total matter

Some of the codes, in particular the EuclidEmulator2, only provide predictions for the power spectrum of total matter, $P_m(k)$. In figure 7, we show the relative agreement of different emulators and other codes predicting $P_m(k)$ to our high-resolution reference run with GADGET-3. The left panel shows the result for the matter power spectrum itself, $P_m(k)$, at $\sum m_\nu = 0.15$ eV. The correlated fluctuations on large scales are due to sample variance which is only present in the reference simulation and not in the predicted spectra. We note that emulators, HMcode, and the halo-model reaction method perform slightly better than the fitting recipe of *halofit*. The results of the EuclidEmulator2 are marginally consistent with the claimed accuracy of 1%, but the neutrino mass lies close to the boundary of parameter space the emulator was trained for. The right panel of figure 7 shows corresponding results for the power suppression factor, $S_m(k)$, with respect to the massless scenario. CosmicEmu shows a disagreement larger than 1% around $k \approx 1 \, h \, \text{Mpc}^{-1}$ where the power suppression is largest. Also *halofit* performs poorly, worse even than linear theory (CLASS). The other codes remain within 1% of the GADGET-3 result up to the particle Nyquist scale. Overall our results are in fair agreement with the detailed comparison carried out by Parimbelli et al. [116].

## 5.2 Bispectra

The nonlinear evolution of matter fluctuations generates a non-vanishing bispectrum, the three-point correlation function of matter in Fourier space, even if non-Gaussianity is negligible.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Total matter power spectrum $P_m(k)$ at $z=0$ for $\sum m_\nu = 0.15$ eV from emulators and fitting methods (left panel), and the respective suppression $S_m(k)$ with respect to the massless case (right panel), compared to the higher-resolution reference GADGET-3 simulation with $L_{\text{box}} = 512 h^{-1} \text{ Mpc}$ and $N_{\text{part}} = 1024^3$. The corresponding particle Nyquist wavenumber of the GADGET-3 run is indicated by the grey dash-dotted line and marks the limit beyond which the estimator of the power spectrum becomes unreliable. Disagreements near and beyond this line are therefore not indicative of errors in the emulators. The grey bands highlight the interval of $\pm 0.01$.</figcaption></figure>

gibly small in the initial conditions. This represents an opportunity for the measurement of neutrino masses, as we expect that the suppression predicted by linear theory is enhanced at the nonlinear level [117–120]. The total matter bispectrum in the presence of massive neutrinos can be schematically defined as

$$B_{mmm}(k_1, k_2, k_3) = (1 - f_\nu)^3 B_{ccc}(k_1, k_2, k_3) + f_\nu (1 - f_\nu)^2 B_{cc\nu}^{(\text{sym})}(k_1, k_2, k_3) + f_\nu^2 (1 - f_\nu) B_{c\nu\nu}^{(\text{sym})}(k_1, k_2, k_3) + f_\nu^3 B_{\nu\nu\nu}(k_1, k_2, k_3), \quad (5.3)$$

where “ccc” denotes the CDM and baryon bispectrum (we do not write “cbcbcb” to avoid clutter) and we note the presence of cross cold-cold-neutrino “cc\nu” and cold-neutrino-neutrino “c\nu\nu” terms that are symmetrised, as indicated by the superscript “(sym)”. In this work, we focus on the bispectrum of CDM and baryons only, but cross terms have also been investigated in the literature, e.g. by Ruggeri et al. [121]. As for the power spectrum case, the leading term is the one of CDM and baryons.

Figure 8 shows a comparison of the bispectrum $B_{ccc}$ as measured at redshift $z=1$ by all simulation codes<fnref n="12" /> for all triangle configurations and all neutrino masses considered, for the set of runs with $N_{\text{part}} = 512^3$. Triangles are plotted as a function of $k_{\text{max}} = \max(k_1, k_2, k_3)$. The comparison shows the maximum percentage difference with respect to GADGET-3 of all triangle configurations at each value of $k_{\text{max}}$. For L-GADGET3, openGADGET3, GADGET-4, PKDGRAV3, CONCEPT, and NM-GADGET4, discrepancies are mostly within 5% for $\sum m_\nu = 0.0$ eV, while steadily growing up to around 10% for massive neutrino cosmologies. Other simulation

<footnote n="12">¹²PINOCCHIO is not used in this test that is mostly focused on nonlinear scales beyond its range of validity.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** The coloured bands show the scatter of $B_{\text{ccc}}$ for all triangle configurations with respect to the measurement from GADGET-3 at redshift $z = 1$. The horizontal axis indicates the maximum wavenumber in each triangle configuration, $k_{\text{max}} = \max(k_1, k_2, k_3)$, in order to present the results in a simple plot. The four panels show different neutrino masses. Each plot uses three subpanels to make the results of the many different codes more visible.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Squeezed bispectrum of CDM and baryons as measured from different codes, relative to GADGET-3 at $z=0$ and $z=1$ for a neutrino mass of $\sum m_\nu = 0.15$eV in the simulation with $L_{\text{box}} = 512 \, h^{-1} \text{Mpc}$ and $N_{\text{part}} = 512^3$ CDM and baryon particles. Results shown are for the squeezed configuration where $k_1 = k_2 \equiv k$, $k_3 = 0.012 \, h \, \text{Mpc}^{-1}$. The grey bands highlight the interval of $\pm 0.01$.</figcaption></figure>

codes are within 10% already at $\sum m_\nu = 0.0$ eV. ANUBIS and gevolution show some of the strongest deviations at large $k_{\text{max}}$, but this is mainly a result of finite resolution as we find much better agreement in the higher-resolution runs. At low $k_{\text{max}}$, strong fluctuations can be observed where measurements can cross zero because of sampling variance, which in turn leads to numerical issues when taking ratios.

In figure 9, we consider specifically a squeezed configuration for which $k_1 = k_2$ and $k_3 = 0.012 \, h \, \text{Mpc}^{-1}$ and show the agreement between different codes for a total neutrino mass of $\sum m_\nu = 0.15$ eV at redshift $z=0$ (left panel) and $z=1$ (right panel). We use measurements from our simulations with $L_{\text{box}} = 512 \, h^{-1} \text{Mpc}$ and $N_{\text{part}} = 512^3$. The relative agreement is better at low redshift, partially due to the fact that the signal amplitude is larger there.

Figure 10 is a comparison of simulation codes for different triangle configurations, for all choices of neutrino masses at redshift $z=1$. We consider four triangle configurations: squeezed, equilateral and two different scalene configurations. Squeezed configurations refer to triangles where one of the side is much shorter than the other two (in Fourier space) which corresponds to looking at the correlation of a distant point with two points close to each other. Equilateral configurations, instead, refer to the correlation of three points at equal distance. Scalene triangles do not have any specific symmetry. In detail, we consider

* squeezed configurations, for which $k_1 = k_2$ and $k_3 = 0.012 \, h \, \text{Mpc}^{-1}$, plotted as a function of $k_1$ as in figure 9;

* equilateral configurations, for which $k_1 = k_2 = k_3$, plotted as a function of $k_1$;

* scalene configurations A, for which $k_1 = 0.79 \, h \, \text{Mpc}^{-1}$, $k_2 = 0.56 \, h \, \text{Mpc}^{-1}$, plotted as a function of $k_3$;
<figure><img src="image_13.png" /><figcaption>**Figure 10.** The four panels show bispectrum measurements at redshift $z = 1$ in the simulations with $L_{\text{box}} = 512 \, h^{-1} \text{ Mpc}$ and $N_{\text{part}} = 512^3$ for different triangle configurations: squeezed (top left panel), equilateral (top right panel), and scalene configurations A and B (bottom panels). In each panel, the top subpanel shows the suppression ratio of the bispectrum of CDM and baryons for three different neutrino masses $\sum m_{\nu} \in \{0.15, 0.3, 0.6\} \text{ eV}$ with respect to the massless case, and the three bottom subpanels show the respective relative differences of the various codes when compared to GADGET-3.</figcaption></figure>
* scalene configurations **B**, for which $k_1 = 0.39 h \, \text{Mpc}^{-1}$, $k_3 = 0.2 h \, \text{Mpc}^{-1}$, plotted as a function of $k_2$.

In each of the four panels of figure 10, the top subpanel shows the ratio between the CDM and baryon bispectrum for massive neutrino cosmologies over massless ones, which we define in analogy to the case of the power spectrum as

$$T_{\text{ccc}}(k_1, k_2, k_3) = \frac{B_{\text{ccc}}^{(\text{massive})}(k_1, k_2, k_3)}{B_{\text{ccc}}^{(\text{massless})}(k_1, k_2, k_3)} \quad (5.4)$$

The three bottom subpanels show the relative differences of the measurements of the suppression ratio in the various codes with respect to GADGET-3 at each of the three neutrino masses, i.e. $\sum m_\nu \in \{0.15, 0.3, 0.6\} \, \text{eV}$ (from top to bottom). For all configurations considered, discrepancies fall broadly within the 5% range. As expected, massive neutrinos suppress the bispectrum of CDM and baryons at all scales, with a stronger effect at smaller scales. For comparison, we also show the tree-level prediction from perturbation theory, using

$$B_{\text{ccc}}^{(\text{tree-level})}(k_1, k_2, k_3) = 2F_2(\mathbf{k}_1, \mathbf{k}_2) P_{\text{cb}}^{\text{L}}(k_1) P_{\text{cb}}^{\text{L}}(k_2) + 2 \text{ permutations}, \quad (5.5)$$

where

$$F_2(\mathbf{k}_1, \mathbf{k}_2) = \frac{5}{7} + \frac{1}{2} \frac{\mathbf{k}_1 \cdot \mathbf{k}_2}{k_1 k_2} \left( \frac{k_1}{k_2} + \frac{k_2}{k_1} \right) + \frac{2}{7} \frac{(\mathbf{k}_1 \cdot \mathbf{k}_2)^2}{k_1^2 k_2^2}, \quad (5.6)$$

and $P_{\text{cb}}^{\text{L}}$ is the linear power spectrum of CDM and baryons generated by CLASS. We can see in figure 10 that the suppression ratio is in good agreement with this prediction when all the three scales $k_1, k_2, k_3$ have a small wavenumber, and that the measured suppression is generally stronger if some of the wavenumbers are large. This is in line with the results seen in the power spectrum, where strong nonlinearities lead to additional suppression.

At the largest scales, i.e. when $k_1, k_2, k_3 \lesssim 0.1 h \, \text{Mpc}^{-1}$, a distinct feature can be discerned which is particularly prominent in the equilateral configuration and for larger neutrino masses: the measurements from the various codes separate into two groups, GADGET-3, openGADGET3, GADGET-4, AREPO, and ANUBIS on the one side, and L-GADGET3, NM-GADGET4, CONCEPT, PKDGRAV3, SWIFT, gevolution, and COLA on the other side. The latter group includes all the codes that employ a mesh-based method, and all of them have means to mitigate against shot noise. We therefore suspect that this dichotomy originates from shot noise in the particle method. This could be tested, e.g. by increasing the number of neutrino particles until convergence is achieved.

As with the case of the power spectrum, we conducted various checks concerning numerical convergence with respect to finite-volume and resolution effects. These show a consistent picture that is in line with what we discussed in section 5.1.2. As an example, figure 11 presents the results for the case of the squeezed configuration in simulations with $N_{\text{part}} = 1024^3$, with a larger volume (left panel) or a higher resolution (right panel) than our simulations with $N_{\text{part}} = 512^3$. In both cases, the agreement on smaller scales is improved: for the larger volume this happens because more independent triangles contribute to each measurement, while for the higher resolution this is due to the better numerical convergence of the density field on small scales.

Overall we may conclude that the N-body methods which produce highly consistent two-point statistics also tend to agree very well on the three-point statistics presented here. NM-GADGET4 appears to be an outlier, showing considerable deviations for squeezed configurations when the sum of the neutrino masses is large. This is however not unexpected since the
<figure><img src="image_14.png" /><figcaption>**Figure 11.** The suppression of the squeezed bispectrum of CDM and baryons relative to the one measured in the GADGET-3 reference runs at $z=1$ for $\sum m_\nu = 0.15$ eV for the simulations with a larger volume, $L_{\text{box}} = 1024\,h^{-1}\,\text{Mpc}$ and $N_{\text{part}} = 1024^3$ (left panel), and at a higher resolution in the small volume, $L_{\text{box}} = 512\,h^{-1}\,\text{Mpc}$ and $N_{\text{part}} = 1024^3$ (right panel). Results shown are for the squeezed configuration where $k_1 = k_2 \equiv k$, $k_3 = 0.012\,h\,\text{Mpc}^{-1}$. The grey bands highlight the interval of $\pm 0.01$.</figcaption></figure>

final gauge transformation from Newtonian motion gauge has been neglected here. This transformation mainly acts at large scales and would therefore affect the squeezed configurations. At low neutrino mass, where the method works best, this effect is almost negligible though. It also becomes minimal at redshift $z=0$ which was set as the target redshift for this method.

## 5.3 Halo mass function

From the halo catalogues produced by **Denhf**, we estimate the halo mass functions (considering only the contribution from CDM and baryons) and compare them to the predictions by Tinker et al. [113], hereafter Tinker10, as well as Despali et al. [111], hereafter Despali16. For these predictions, we use the linear power spectra of CDM and baryons calculated by **CLASS** for the respective neutrino cosmologies in the modelling of the theoretical halo mass functions. It has been shown by Costanzi et al. [122] that this approach reproduces the halo mass function well for neutrino cosmologies. Figure 12 shows the ratio of the halo mass functions relative to the halo mass function of GADGET-3 at a neutrino mass of $\sum m_\nu = 0.15$ eV for the runs with $N_{\text{part}} = 1024^3$, in the large volume where $L_{\text{box}} = 1024\,h^{-1}\,\text{Mpc}$ (left panel), as well as for the higher-resolution setup with $L_{\text{box}} = 512\,h^{-1}\,\text{Mpc}$ (right panel). At the high-mass end, the agreement between different codes is generally very good, and fluctuations are smaller in the larger volume due to better statistics. At low masses $M_{200b} < 10^{13}\,h^{-1}\,\text{M}_\odot$, the number density of halos is underestimated by **COLA**, **gevolution**, and **ANUBIS** by up to 50% for the lower resolution. At higher resolution the agreement improves. The relatively poor performance of **gevolution** in predicting the halo mass function can be understood from the fact that the code uses a uniform mesh. This leads to a smoothing of small-scale structures and generally to a mass estimate of halos that is poorly converged at the low-mass end. **COLA** suffers from the same limitation, but the simulations used a mesh with significantly higher
<figure><img src="image_15.png" /><figcaption>Figure 12. Halo mass functions relative to the one of GADGET-3, for $N_{\text{part}} = 1024^3$ at $z = 0$ and neutrino mass $\sum m_{\nu} = 0.15$ eV. The result for the larger-volume runs with $L_{\text{box}} = 1024 h^{-1} \text{ Mpc}$ is shown in the left panel while the right panel shows the result for the higher-resolution runs with $L_{\text{box}} = 512 h^{-1} \text{ Mpc}$.</figcaption></figure>

resolution in this case. Specifically, in the runs with $N_{\text{part}} = 1024^3$, COLA used a mesh of 3072³ grid points, while gevolution used a mesh of 2048³ grid points.<fnref n="13" />

Figure 13 shows the suppression of the halo mass function due to neutrinos with masses $\sum m_{\nu} \in \{0.15, 0.3, 0.6\}$ eV at redshift $z = 0$ (left panel) and $z = 1$ (right panel). At low halo masses, there is little suppression, while going to higher masses the number density of halos is more and more suppressed. The higher the neutrino masses, and the higher the redshift, the stronger the suppression: at $z = 0$ and $\sum m_{\nu} = 0.15$ eV the suppression goes down to a factor of 0.9 at halo masses of $10^{14} h^{-1} M_{\odot}$, while at $z = 1$ and $\sum m_{\nu} = 0.6$ eV the suppression goes down to 0.4 at the same halo mass. In analogy to the case of the power spectrum, we define the suppression ratio with respect to the massless case as

$$R = \frac{\frac{dn^{\text{(massive)}}}{d \ln M}}{\frac{dn^{\text{(massless)}}}{d \ln M}}. \qquad (5.7)$$

Figure 14 shows this suppression ratio relative to the one measured from GADGET-3 for a sum of neutrino masses of $\sum m_{\nu} = 0.15$ eV. The different codes generally agree within 3% on the suppression ratio, even in cases where the halo mass function was poorly converged in figure 12. Our data show a trend that the suppression is about 1% stronger than predicted by the models of Tinker10 and Despali16. The COLA results are a slight outlier, agreeing more with these models than with the other simulations.

<footnote n="13">¹³Running gevolution on such a fine mesh, i.e. with less than one particle per cell on average, was not supported in the public release of the code. In such a situation, the evolution would become unstable due to the low order of the finite-difference gradients used in the particle update. A second-order gradient computation was therefore implemented for this work, a feature that is made available in a recent patch of the code.</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure 13.** The halo mass function for different neutrino masses, relative to the case with massless neutrinos for each code, respectively. We show the results for the simulations with $L_{\text{box}} = 512 \ h^{-1} \text{ Mpc}$ and $N_{\text{part}} = 512^3$ at $z = 0$ (left panel) and at $z = 1$ (right panel).</figcaption></figure>

## 5.4 Halo bias

Finally, we study the halo bias for a fixed selection of halos defined by a mass threshold of $M_{200b} > 10^{13} h^{-1} M_{\odot}$. However, since the halo mass function shows considerable differences between the different simulations, sometimes due to the fact that the mass estimate is not well converged at the low-mass end (certainly for gevolution, COLA, and ANUBIS), we apply the halo selection as follows. First, we select the halos above the mass threshold in the reference runs done with GADGET-3. We may call the size of the selected population $N_h$. Then, for each other code, we generate the sample by selecting the $N_h$ most massive halos. The reasoning for this approach is that, while the estimated masses of individual halos may differ significantly between different codes, we still expect there to be a tight correlation that largely preserves the mass ordering. Another way to think about this is to consider a simple abundance matching of $N_h$ sources, assigned to the centers of the most massive halos. We compare the bias measurements to the prediction by Tinker et al. [113] (Tinker10). Here the large-scale bias is estimated from the mass-dependent peak height of halos in the linear density field. Given the GADGET-3 halo masses, we model the peak heights using the linear power spectrum of CDM and baryons calculated by CLASS for the respective neutrino cosmology. We then take the average of all biases obtained for each halo mass to get the final prediction. Note that the prediction of Tinker10 is only expected to work in the linear regime.

Figure 15 shows the bias measurements from the simulations with larger volume, $L_{\text{box}} = 1024 h^{-1} \text{ Mpc}$. We define the scale-dependent halo bias $b(k)$ as the ratio of the cross-power spectrum of halos with CDM and baryons and the auto-power spectrum of CDM and baryons, i.e.

$$b(k) = \frac{P_{h,cb}(k)}{P_{cb}(k)}. \qquad (5.8)$$
<figure><img src="image_17.png" /><figcaption>Figure 14. Suppression of the halo mass functions relative to the one measured in GADGET-3 at redshift $z = 0$ and for a total neutrino mass of $\sum m_\nu = 0.15$ eV. We show the results for the simulations with $L_{\text{box}} = 1024 h^{-1} \text{ Mpc}$ and $N_{\text{part}} = 1024^3$ in the left panel, and $L_{\text{box}} = 512 h^{-1} \text{ Mpc}$ and $N_{\text{part}} = 1024^3$ in the right panel. The grey bands highlight the interval of $\pm 0.01$.</figcaption></figure>

It has been shown by Castorina et al. [77] that defining the bias factor with respect to cold species gives closer-to-universal and less scale-dependent results than using the same definition with respect to total matter. As can be seen in the left panel of figure 15, the bias measurements agree reasonably well on large scales except for gevolution where the bias is measured to be about 4% larger, and for PINOCCHIO where the bias is measured to be about 10% smaller. In analogy to the case of the power spectrum, we again define a bias ratio with respect to the massless case, which in this situation will quantify the increase (rather than suppression) of the bias in the presence of massive neutrinos,

$$Q(k) = \frac{b^{\text{(massive)}}(k)}{b^{\text{(massless)}}(k)}. \qquad (5.9)$$

Results for this bias ratio are shown in the right panel of figure 15. Here the agreement between different codes is excellent, well within 1% over almost the entire range of scales probed. For PINOCCHIO the bias ratio is about 1% accurate up to $k \simeq 0.3 h \text{ Mpc}^{-1}$.

To study the robustness of our results with respect to the mass resolution of the simulations we repeat the bias measurements in the runs with higher resolution, i.e. with $L_{\text{box}} = 512 h^{-1} \text{ Mpc}$ and $N_{\text{part}} = 1024^3$. The smaller simulation volume leads to a higher level of shot noise in the halo counts, which incurs somewhat larger fluctuations when compared to the larger volume. As can be seen in figure 16, the agreement between the different codes is improved significantly, in particular for codes that have difficulties in predicting the halo mass function accurately (gevolution, COLA, and ANUBIS). The results for PINOCCHIO do not improve significantly, as the discrepancy is mainly due to the approximate nature of the method rather than lack of resolution.

It is worth pointing out that the bias ratio $Q$ has a shape similar to the inverted power spectrum ratio $S_{\text{cb}}$. This is actually expected, and was recently discussed by Hassani et
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Halo bias with respect to CDM and baryons at $z=0$ and neutrino mass $\sum m_\nu = 0.15$ eV for the simulations with $L_{\text{box}} = 1024 \, h^{-1} \text{Mpc}$ and $N_{\text{part}} = 1024^3$. In the **GADGET-3** reference simulations, all halos with $M_{200b} > 10^{13} \, h^{-1} M_{\odot}$ are selected, providing a sample of size $N_{\text{h}}$. For the other simulations, we then select the most massive $N_{\text{h}}$ halos. The grey bands in the lower panels highlight the interval of $\pm 0.01$.</figcaption></figure>

al. [123]. It means that the power spectrum of halos, when selected at fixed mass threshold, is much less sensitive to the neutrino mass than the power spectrum of CDM and baryons. On the other hand, synthetic catalogues are often created in such a way that the observed abundance of a certain type of object is reproduced (abundance matching). In such a situation, it may be more appropriate to study the dependence of the bias on the neutrino mass at fixed number count. We therefore repeat our measurements, but keeping the size of all samples fixed at $N_h^{\text{(massless)}}$, which is the number of halos with $M_{200b} > 10^{13} \, h^{-1} M_{\odot}$ in the **GADGET-3** reference simulation at zero neutrino mass. In other words, when selecting the halo sample for non-zero neutrino mass, we still select the $N_h^{\text{(massless)}}$ most massive halos in all simulations. This effectively reduces the mass threshold of the selection for the massive neutrino case when compared to the previous procedure.

Figure 17 shows the results of the bias measurement obtained through this procedure, using the higher-resolution simulations. We observe that the bias still increases with neutrino mass, but not quite as much as in figure 16 where a fixed halo mass threshold was used. The large-scale bias indicated by the dotted line (Tinker10) is about 0.01, or one percentage point, lower so that the corresponding bias ratio $Q$ drops by 0.7%. The largest effect in $Q$ appears around $k \approx 1 \, h \, \text{Mpc}^{-1}$ were the change can be as much as 2% due to the different halo selection. The interplay between bias enhancement and the suppression of the matter power
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Halo bias with respect to CDM and baryons at $z=0$ and neutrino mass $\sum m_\nu = 0.15$ eV for the higher-resolution simulations with $L_{\text{box}} = 512 \, h^{-1} \text{Mpc}$ and $N_{\text{part}} = 1024^3$. In the GADGET-3 reference simulations, all halos with $M_{200b} > 10^{13} \, h^{-1} \text{M}_\odot$ are selected, providing a sample of size $N_h$. For the other simulations, we then select the most massive $N_h$ halos. The grey bands in the lower panels highlight the interval of $\pm 0.01$.</figcaption></figure>

spectrum has another layer of complexity due to the way in which the sample is selected. This needs to be studied carefully in the context of the specific numerical recipes that are employed in the production of synthetic catalogues.

# 6 Discussion

Accurate and reliable modelling of the signatures that the neutrino mass imprints on the observables used to test the cosmological model is an essential ingredient for the data analysis of all upcoming galaxy surveys, and in particular for *Euclid*. Such modelling necessarily requires a self-consistent description of the linear and possibly nonlinear clustering of neutrinos along with the nonlinear evolution of dark matter and baryonic structures. By comparing results across different implementations, including eleven full N-body implementations, two N-body schemes with fast time integration based on Lagrangian perturbation theory, and a further four codes that predict the nonlinear matter power spectra directly, we establish that current numerical techniques are in sub-percent agreement with regards to modelling the impact of massive neutrinos on the most common summary statistics of cosmological large-scale structure. We identify several specific situations where larger modelling errors can occur, but such shortfalls are generally well understood in terms of approximations or other compromises that were made in these situations. Our results can therefore be used as
<figure><img src="image_20.png" /><figcaption>**Figure 17.** Halo bias with respect to CDM and baryons at $z=0$ and neutrino mass $\sum m_\nu = 0.15$ eV for the higher-resolution simulations with $L_{\text{box}} = 512 \, h^{-1} \text{Mpc}$ and $N_{\text{part}} = 1024^3$. In the GADGET-3 reference simulation for $\sum m_\nu = 0$ eV, all halos with $M_{200b} > 10^{13} h^{-1} M_{\odot}$ are selected, providing a sample of size $N_h^{(\text{massless})}$. We then select the most massive $N_h^{(\text{massless})}$ halos in all the other simulations, even for those with non-zero neutrino mass. The grey bands in the lower panels highlight the interval of $\pm 0.01$.</figcaption></figure>

a detailed guide for choosing the preferred modelling techniques for any application given its requirements in terms of resources, accuracy, and quantities that need to be modelled. The validation presented here is the crucial first step for building up confidence in the numerical tools employed in the data analysis pipeline of *Euclid*. It is particularly vital when considering the actual measurement of the neutrino mass scale from the data, which is one of the key science goals of the mission.

The fastest method of predicting simple summary statistics like the nonlinear matter power spectrum are emulators, and they will therefore play a crucial role in the cosmological likelihood analysis of *Euclid* and other large-scale structure surveys. They are of course many orders of magnitude faster than simulations but tend to outperform even semi-analytic models which often have some bottlenecks in their numerical evaluation. However, emulators can only be as accurate as the simulations they are trained on, and it is therefore important to understand the modelling errors of simulations too. We find that the current state of the art for emulators yield an absolute precision on the power spectrum of total matter better than 2% and can predict the relative change due to the neutrino mass to better than 1% on all scales considered in this work. Interestingly, the best semi-analytic fitting methods available, in particular *ReACT* and *HMcode*, can achieve similar performance.
Overall our results demonstrate that we are in a fairly comfortable position, with several independent numerical techniques at our disposal that produce consistent results at the sub-percent level if we employ them diligently. The community has also implemented such techniques in a large number of different N-body codes, such that there is no shortage in choice of which code one wants to use. Moreover, our detailed comparison of particle-based and mesh-based techniques shows that the assumption of linear neutrinos is clearly sufficient to reach percent accuracy, even up to scales of the order of $k \approx 7 h \text{ Mpc}^{-1}$ relevant for predicting the weak-lensing signal in *Euclid*. We note that some codes have inherent difficulties reaching such levels of absolute accuracy due to effects of finite resolution. This is obvious in the cases where a uniform mesh is employed in the computation of gravitational interactions (*gevolution* and *COLA*), but AMR does not solve the issue entirely as the example of *ANUBIS* illustrates. The relative impact of massive neutrinos can nonetheless be predicted very accurately with those codes. Here we do of course not attempt to address the additional challenge of modelling baryonic effects, i.e. astrophysical processes, down to such scales as this can be treated separately, see e.g. Martinelli et al. [124] for a discussion. On mildly nonlinear scales and in particular at redshifts $z \gtrsim 1$ relevant for *Euclid*, a “sophisticated 3LPT” realisation like the one produced with *PINOCCHIO*, based on a scale-dependent linear growth rate computed from *CAMB* and propagated to second- and third-order LPT with standard techniques, can be useful to produce a large number of halo catalogues in a limited amount of computing time, see e.g. Fumagalli et al. [10].

The summary statistics considered in our analysis include auto- and cross-power spectra of the CDM and baryon component and neutrinos, bispectra of the CDM and baryon component, halo mass functions, and halo bias. We present results for redshifts $z=0$ and $z=1$, relevant for galaxy surveys like *Euclid*. We do not consider redshift-space distortions or other effects that occur due to taking observations on our past light cone, and leave a detailed investigation of these effects to future work. However, we expect that no big surprises would appear given our level of confidence in the modelling of the summary statistics presented here.

In order to aid future code development, we make our reference simulations and analysis pipelines available via a public repository (see data availability statement below). This provides a reliable baseline against which further numerical methods can be validated, and it showcases the current state-of-the-art in modelling massive neutrinos in cosmology.

## Acknowledgments

JA acknowledges financial support from the Swiss National Science Foundation. BB was supported by a U.K. Research and Innovation Stephen Hawking Fellowship (EP/W005654/1). MBa acknowledges support by the project “Combining Cosmic Microwave Background and Large Scale Structure data: an Integrated Approach for Addressing Fundamental Questions in Cosmology”, funded by the MIUR Progetti di Ricerca di Rilevante Interesse Nazionale (PRIN) Bando 2017 — grant 2017YJYZAH. CH-A acknowledges support from the Excellence Cluster ORIGINS which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy — EXC-2094 — 390783311. BSW is supported by a Royal Society Enhancement Award (grant no. RGF\EA\181023). RM, DFM, HAW, FH thank the Research Council of Norway for their support and our computations were performed on resources provided by UNINETT Sigma2 — the National Infrastructure for High Performance Computing and Data Storage in Norway. CM acknowledges support from a U.K. Research and Innovation Future Leaders Fellowship [grant
MR/S016066/1]. The PINOCCHIO simulations were performed on Queen Mary's Apocrita HPC facility, supported by QMUL Research-IT. CG acknowledges the support from the grant PRIN-MIUR 2017 WSCC32 ZOOMING and from the Italian National Institute of Astrophysics under the grant "Bando PRIN 2019," PI: Viola Allevato. TC is supported by the INFN INDARK PD51 grant and the FARE MIUR grant 'ClustersXEuclid' R165SBKTMA. KD acknowledges support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy — EXC-2094 — 390783311 and the COMPLEX project from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program grant agreement ERC-2019-AdG 882679. The openGADGET3 simulations were carried out at the Leibniz Supercomputer Center (LRZ) under the project pr86re. CA and BL are supported by an European Research Council Starting Grant (ERC-StG-716532), and BL also acknowledges support by the U.K. Science and Technology Facilities Council Consolidated Grants No. ST/I00162X/1 and ST/P000541/1. KK is supported by the U.K. Science and Technology Facilities Council (grant numbers ST/S000550/1 and ST/W001225/1). MV and GP are supported by the INFN INDARK PD 51 grant and by the ASI-INAF n. 2017-14-H.0 agreement. REA and MZ acknowledge the support of the ERC-STG number 716151 (BACCO). This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (https://dirac.ac.uk/). The equipment was funded by BEIS via STFC capital grants ST/K00042X/1, ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and STFC operation grant ST/R000832/1. DiRAC is part of the U.K. National e-Infrastructure.

The Euclid Consortium acknowledges the European Space Agency and a number of agencies and institutes that have supported the development of Euclid, in particular the Academy of Finland, the Agenzia Spaziale Italiana, the Belgian Science Policy, the Canadian Euclid Consortium, the French Centre National d'Etudes Spatiales, the Deutsches Zentrum für Luft- und Raumfahrt, the Danish Space Research Institute, the Fundação para a Ciência e a Tecnologia, the Ministerio de Ciencia e Innovación, the National Aeronautics and Space Administration, the National Astronomical Observatory of Japan, the Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the Romanian Space Agency, the State Secretariat for Education, Research and Innovation (SERI) at the Swiss Space Office (SSO), and the United Kingdom Space Agency. A complete and detailed list is available on the Euclid web site (http://www.euclid-ec.org).

For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this work and upload the accepted version to arXiv with CC BY licence.

**Data availability.** Raw data from our reference runs (GADGET-3), and all reduced data (summary statistics) presented in this work are available at https://doi.org/10.5281/zenodo.7868793. We also provide initial data, parameter specifications for CLASS, and a documented analysis pipeline to facilitate further validations against our results.

## References

[1] R. Laureijs et al., Euclid definition study report, arXiv:1110.3193 [INSPIRE].

[2] E. Di Valentino, A. Melchiorri and J. Silk, Cosmological constraints in extended parameter space from the Planck 2018 legacy release, JCAP 01 (2020) 013 [arXiv:1908.01391] [INSPIRE].
[3] EBOSS collaboration, *Completed SDSS-IV extended Baryon Oscillation Spectroscopic Survey: cosmological implications from two decades of spectroscopic surveys at the Apache Point Observatory*, Phys. Rev. D 103 (2021) 083533 [arXiv:2007.08991] [INSPIRE].

[4] N. Palanque-Delabrouille et al., *Hints, neutrino bounds and WDM constraints from SDSS DR14 Lyman-α and Planck full-survey data*, JCAP 04 (2020) 038 [arXiv:1911.09073] [INSPIRE].

[5] P.L. Taylor, T.D. Kitching and J.D. McEwen, *Preparing for the cosmic shear data flood: optimal data extraction and simulation requirements for stage IV dark energy experiments*, Phys. Rev. D 98 (2018) 043532 [arXiv:1804.03667] [INSPIRE].

[6] EUCLID collaboration: A. Blanchard et al., *Euclid preparation: VII. Forecast validation for Euclid cosmological probes*, Astron. Astrophys. 642 (2020) A191 [arXiv:1910.09273] [INSPIRE].

[7] EUCLID collaboration: M. Knabenhans et al., *Euclid preparation: II. The EuclidEmulator — a tool to compute the cosmology dependence of the nonlinear matter power spectrum*, Mon. Not. Roy. Astron. Soc. 484 (2019) 5509 [arXiv:1809.04695] [INSPIRE].

[8] EUCLID collaboration: M. Knabenhans et al., *Euclid preparation: IX. EuclidEmulator2 — power spectrum emulation with massive neutrinos and self-consistent dark energy perturbations*, Mon. Not. Roy. Astron. Soc. 505 (2021) 2840 [arXiv:2010.11288] [INSPIRE].

[9] EUCLID collaboration: T. Castro et al., *Euclid preparation: XXIV. Calibration of the halo mass function in $\Lambda(\nu)$ CDM cosmologies*, Astron. Astrophys. 671 (2023) A100 [arXiv:2208.02174] [INSPIRE].

[10] A. Fumagalli et al., *Euclid: effects of sample covariance on the number counts of galaxy clusters*, Astron. Astrophys. 652 (2021) A21 [arXiv:2102.08914] [INSPIRE].

[11] EUCLID collaboration: A. Fumagalli et al., *Euclid preparation. XXVII. Covariance model validation for the 2-point correlation function of galaxy clusters*, arXiv:2211.12965 [INSPIRE].

[12] R.E. Angulo and O. Hahn, *Large-scale dark matter simulations*, arXiv:2112.05165 [DOI:10.1007/s41115-021-00013-z] [INSPIRE].

[13] J. Brandbyge, S. Hannestad, T. Haugbølle and B. Thomsen, *The effect of thermal neutrino motion on the non-linear cosmological matter power spectrum*, JCAP 08 (2008) 020 [arXiv:0802.3700] [INSPIRE].

[14] M. Viel, M.G. Haehnelt and V. Springel, *The effect of neutrinos on the matter distribution as probed by the intergalactic medium*, JCAP 06 (2010) 015 [arXiv:1003.2422] [INSPIRE].

[15] S. Agarwal and H.A. Feldman, *The effect of massive neutrinos on the matter power spectrum*, Mon. Not. Roy. Astron. Soc. 410 (2011) 1647 [arXiv:1006.0689] [INSPIRE].

[16] S. Bird, M. Viel and M.G. Haehnelt, *Massive neutrinos and the non-linear matter power spectrum*, Mon. Not. Roy. Astron. Soc. 420 (2012) 2551 [arXiv:1109.4416] [INSPIRE].

[17] F. Villaescusa-Navarro et al., *Cosmology with massive neutrinos I: towards a realistic modeling of the relation between matter, haloes and galaxies*, JCAP 03 (2014) 011 [arXiv:1311.0866] [INSPIRE].

[18] E. Castorina et al., *DEMNUni: the clustering of large-scale structures in the presence of massive neutrinos*, JCAP 07 (2015) 043 [arXiv:1505.07148] [INSPIRE].

[19] J.D. Emberson et al., *Cosmological neutrino simulations at extreme scale*, Res. Astron. Astrophys. 17 (2017) 085 [arXiv:1611.01545] [INSPIRE].

[20] J. Adamek, R. Durrer and M. Kunz, *Relativistic N-body simulations with massive neutrinos*, JCAP 11 (2017) 004 [arXiv:1707.06938] [INSPIRE].

[21] A. Banerjee, D. Powell, T. Abel and F. Villaescusa-Navarro, *Reducing noise in cosmological N-body simulations with neutrinos*, JCAP 09 (2018) 028 [arXiv:1801.03906] [INSPIRE].

[22] J. Brandbyge, S. Hannestad and T. Tram, *Momentum space sampling of neutrinos in N-body simulations*, JCAP 03 (2019) 047 [arXiv:1806.05874] [INSPIRE].
[23] J. Brandbyge and S. Hannestad, Grid based linear neutrino perturbations in cosmological N-body simulations, JCAP 05 (2009) 002 [arXiv:0812.3149] [INSPIRE].

[24] Y. Ali-Haimoud and S. Bird, An efficient implementation of massive neutrinos in non-linear structure formation simulations, Mon. Not. Roy. Astron. Soc. 428 (2012) 3375 [arXiv:1209.0461] [INSPIRE].

[25] J. Liu et al., MassiveNuS: cosmological massive neutrino simulations, JCAP 03 (2018) 049 [arXiv:1711.10524] [INSPIRE].

[26] J.Z. Chen, A. Upadhye and Y.Y.Y. Wong, One line to run them all: SuperEasy massive neutrino linear response in N-body simulations, JCAP 04 (2021) 078 [arXiv:2011.12504] [INSPIRE].

[27] J. Dakin et al., νCONCEPT: cosmological neutrino simulations from the non-linear Boltzmann hierarchy, JCAP 02 (2019) 052 [arXiv:1712.03944] [INSPIRE].

[28] K. Yoshikawa, S. Tanaka, N. Yoshida and S. Saito, Cosmological Vlasov-Poisson simulations of structure formation with relic neutrinos: nonlinear clustering and the neutrino mass, Astrophys. J. 904 (2020) 159 [arXiv:2010.00248] [INSPIRE].

[29] J. Brandbyge and S. Hannestad, Resolving cosmic neutrino structure: a hybrid neutrino N-body scheme, JCAP 01 (2010) 021 [arXiv:0908.1969] [INSPIRE].

[30] S. Bird, Y. Ali-Haïmoud, Y. Feng and J. Liu, An efficient and accurate hybrid method for simulating non-linear neutrino structure, Mon. Not. Roy. Astron. Soc. 481 (2018) 1486 [arXiv:1803.09854] [INSPIRE].

[31] W. Elbers et al., An optimal non-linear method for simulating relic neutrinos, Mon. Not. Roy. Astron. Soc. 507 (2021) 2614 [arXiv:2010.07321] [INSPIRE].

[32] C. Partmann, C. Fidler, C. Rampf and O. Hahn, Fast simulations of cosmic large-scale structure with massive neutrinos, JCAP 09 (2020) 018 [arXiv:2003.07387] [INSPIRE].

[33] M. Zennaro et al., How to add massive neutrinos to your ΛCDM simulation — extending cosmology rescaling algorithms, Mon. Not. Roy. Astron. Soc. 489 (2019) 5938 [arXiv:1905.08696] [INSPIRE].

[34] I. Esteban et al., The fate of hints: updated global analysis of three-flavor neutrino oscillations, JHEP 09 (2020) 178 [arXiv:2007.14792] [INSPIRE].

[35] KATRIN collaboration, Improved upper limit on the neutrino mass from a direct kinematic method by KATRIN, Phys. Rev. Lett. 123 (2019) 221802 [arXiv:1909.06048] [INSPIRE].

[36] KATRIN collaboration, Direct neutrino-mass measurement with sub-electronvolt sensitivity, Nature Phys. 18 (2022) 160 [arXiv:2105.08533] [INSPIRE].

[37] E. Di Valentino, S. Gariazzo and O. Mena, Most constraining cosmological neutrino mass bounds, Phys. Rev. D 104 (2021) 083504 [arXiv:2106.15267] [INSPIRE].

[38] M. Archidiacono, T. Brinckmann, J. Lesgourgues and V. Poulin, Physical effects involved in the measurements of neutrino masses with future cosmological data, JCAP 02 (2017) 052 [arXiv:1610.09852] [INSPIRE].

[39] J.R. Bond, G. Efstathiou and J. Silk, Massive neutrinos and the large scale structure of the universe, Phys. Rev. Lett. 45 (1980) 1980 [INSPIRE].

[40] W. Hu, D.J. Eisenstein and M. Tegmark, Weighing neutrinos with galaxy surveys, Phys. Rev. Lett. 80 (1998) 5255 [astro-ph/9712057] [INSPIRE].

[41] S. Hannestad, A. Upadhye and Y.Y.Y. Wong, Spoon or slide? The non-linear matter power spectrum in the presence of massive neutrinos, JCAP 11 (2020) 062 [arXiv:2006.04995] [INSPIRE].

[42] V. Springel, The cosmological simulation code GADGET-2, Mon. Not. Roy. Astron. Soc. 364 (2005) 1105 [astro-ph/0505010] [INSPIRE].
[43] V. Springel et al., *The Aquarius project: the subhalos of galactic halos*, Mon. Not. Roy. Astron. Soc. 391 (2008) 1685 [arXiv:0809.0898] [INSPIRE].

[44] R.E. Angulo et al., *Scaling relations for galaxy clusters in the Millennium-XXL simulation*, Mon. Not. Roy. Astron. Soc. 426 (2012) 2046 [arXiv:1203.3216] [INSPIRE].

[45] R.E. Angulo et al., *The BACCO simulation project: exploiting the full power of large-scale structure for cosmology*, Mon. Not. Roy. Astron. Soc. 507 (2021) 5869 [arXiv:2004.06245] [INSPIRE].

[46] A.M. Beck et al., *An improved SPH scheme for cosmological simulations*, Mon. Not. Roy. Astron. Soc. 455 (2016) 2110 [arXiv:1502.07358] [INSPIRE].

[47] T. Marin-Gilabert, M. Valentini, U.P. Steinwandel and K. Dolag, *The role of physical and numerical viscosity in hydrodynamical instabilities*, Mon. Not. Roy. Astron. Soc. 517 (2022) 5971 [arXiv:2205.09135] [INSPIRE].

[48] V. Springel, R. Pakmor, O. Zier and M. Reinecke, *Simulating cosmic structure formation with the gadget-4 code*, Mon. Not. Roy. Astron. Soc. 506 (2021) 2871 [arXiv:2010.03567] [INSPIRE].

[49] P. Heuschling, C. Partmann and C. Fidler, *A minimal model for massive neutrinos in Newtonian N-body simulations*, JCAP 09 (2022) 068 [arXiv:2201.13186] [INSPIRE].

[50] V. Springel, *E pur si muove: Galilean-invariant cosmological hydrodynamical simulations on a moving mesh*, Mon. Not. Roy. Astron. Soc. 401 (2010) 791 [arXiv:0901.4107] [INSPIRE].

[51] R. Weinberger, V. Springel and R. Pakmor, *The Arepo public code release*, Astrophys. J. Suppl. 248 (2020) 32 [arXiv:1909.04667] [INSPIRE].

[52] J. Dakin, S. Hannestad and T. Tram, *The cosmological simulation code CONCEPT 1.0*, Mon. Not. Roy. Astron. Soc. 513 (2022) 991 [arXiv:2112.01508] [INSPIRE].

[53] D. Potter, J. Stadel and R. Teyssier, *PKDGRAV3: beyond trillion particle cosmological simulations for the next era of galaxy surveys*, arXiv:1609.08621 [INSPIRE].

[54] M. Schaller et al., *SWIFT: SPH with inter-dependent fine-grained tasking*, Astrophysics source code library https://ascl.net/1805.020, May 2018.

[55] W. Elbers, *Geodesic motion and phase-space evolution of massive neutrinos*, JCAP 11 (2022) 058 [arXiv:2207.14256] [INSPIRE].

[56] R. Teyssier, *Cosmological hydrodynamics with adaptive mesh refinement: a new high resolution code called Ramses*, Astron. Astrophys. 385 (2002) 337 [astro-ph/0111367] [INSPIRE].

[57] R. Mauland, O. Elgarøy, D.F. Mota and H.A. Winther, *The void-galaxy cross-correlation function with massive neutrinos and modified gravity*, arXiv:2303.05820 [INSPIRE].

[58] J. Adamek, D. Daverio, R. Durrer and M. Kunz, *General relativity and cosmic structure formation*, Nature Phys. 12 (2016) 346 [arXiv:1509.01699] [INSPIRE].

[59] J. Adamek, D. Daverio, R. Durrer and M. Kunz, *Gevolution: a cosmological N-body code based on general relativity*, JCAP 07 (2016) 053 [arXiv:1604.06065] [INSPIRE].

[60] H.A. Winther et al., *COLA with scale-dependent growth: applications to screened modified gravity models*, JCAP 08 (2017) 006 [arXiv:1703.00879] [INSPIRE].

[61] B.S. Wright, H.A. Winther and K. Koyama, *COLA with massive neutrinos*, JCAP 10 (2017) 054 [arXiv:1705.08165] [INSPIRE].

[62] P. Monaco et al., *An accurate tool for the fast generation of dark matter halo catalogs*, Mon. Not. Roy. Astron. Soc. 433 (2013) 2389 [arXiv:1305.1505] [INSPIRE].

[63] L.A. Rizzo et al., *Simulating cosmologies beyond ΛCDM with PINOCCHIO*, JCAP 01 (2017) 008 [arXiv:1610.07624] [INSPIRE].

[64] M. Cataneo et al., *On the road to per cent accuracy — III. Non-linear reaction of the matter power spectrum to massive neutrinos*, Mon. Not. Roy. Astron. Soc. 491 (2020) 3101 [arXiv:1909.02561] [INSPIRE].
[65] B. Bose et al., On the road to per cent accuracy — V. The non-linear power spectrum beyond $\Lambda$CDM with massive neutrinos and baryonic feedback, Mon. Not. Roy. Astron. Soc. 508 (2021) 2479 [arXiv:2105.12114] [INSPIRE].

[66] E. Lawrence et al., The Mira-Titan universe II: matter power spectrum emulation, Astrophys. J. 847 (2017) 50 [arXiv:1705.03388] [INSPIRE].

[67] K.R. Moran et al., The Mira-Titan universe — IV. High-precision power spectrum emulation, Mon. Not. Roy. Astron. Soc. 520 (2023) 3443 [arXiv:2207.12345] [INSPIRE].

[68] A. Banerjee and N. Dalal, Simulating nonlinear cosmological structure formation with massive neutrinos, JCAP 11 (2016) 015 [arXiv:1606.06167] [INSPIRE].

[69] C.-P. Ma and E. Bertschinger, A calculation of the full neutrino phase space in cold + hot dark matter models, Astrophys. J. 429 (1994) 22 [astro-ph/9308006] [INSPIRE].

[70] S. Tassev, M. Zaldarriaga and D. Eisenstein, Solving large scale structure in ten easy steps with COLA, JCAP 06 (2013) 036 [arXiv:1301.0322] [INSPIRE].

[71] A. Izard, M. Crocce and P. Fosalba, ICE-COLA: towards fast and accurate synthetic galaxy catalogues optimizing a quasi N-body method, Mon. Not. Roy. Astron. Soc. 459 (2016) 2327 [arXiv:1509.04685] [INSPIRE].

[72] A. Lewis, A. Challinor and A. Lasenby, Efficient computation of CMB anisotropies in closed FRW models, Astrophys. J. 538 (2000) 473 [astro-ph/9911177] [INSPIRE].

[73] A. Lewis and S. Bridle, Cosmological parameters from CMB and other data: a Monte Carlo approach, Phys. Rev. D 66 (2002) 103511 [astro-ph/0205436] [INSPIRE].

[74] P. Monaco, T. Theuns and G. Taffoni, Pinocchio: pinpointing orbit-crossing collapsed hierarchical objects in a linear density field, Mon. Not. Roy. Astron. Soc. 331 (2002) 587 [astro-ph/0109323] [INSPIRE].

[75] E. Munari et al., Improving fast generation of halo catalogues with higher order Lagrangian perturbation theory, Mon. Not. Roy. Astron. Soc. 465 (2017) 4658 [arXiv:1605.04788] [INSPIRE].

[76] E. Munari et al., Testing approximate predictions of displacements of cosmological dark matter halos, JCAP 07 (2017) 050 [arXiv:1704.00920] [INSPIRE].

[77] E. Castorina et al., Cosmology with massive neutrinos II: on the universality of the halo mass function and bias, JCAP 02 (2014) 049 [arXiv:1311.1212] [INSPIRE].

[78] J.F. Navarro, C.S. Frenk and S.D.M. White, A universal density profile from hierarchical clustering, Astrophys. J. 490 (1997) 493 [astro-ph/9611107] [INSPIRE].

[79] L. Blot et al., Comparing approximate methods for mock catalogues and covariance matrices II: power spectrum multipoles, Mon. Not. Roy. Astron. Soc. 485 (2019) 2806 [arXiv:1806.09497] [INSPIRE].

[80] C. Fidler et al., A new approach to cosmological structure formation with massive neutrinos, JCAP 01 (2019) 025 [arXiv:1807.03701] [INSPIRE].

[81] M. Cataneo et al., On the road to percent accuracy: non-linear reaction of the matter power spectrum to dark energy and modified gravity, Mon. Not. Roy. Astron. Soc. 488 (2019) 2121 [arXiv:1812.05594] [INSPIRE].

[82] A. Mead, S. Brieden, T. Tröster and C. Heymans, HMcode-2020: improved modelling of non-linear cosmological power spectra with baryonic feedback, arXiv:2009.01858 [DOI:10.1093/mnras/stab082] [INSPIRE].

[83] A. Cooray and R.K. Sheth, Halo models of large scale structure, Phys. Rept. 372 (2002) 1 [astro-ph/0206508] [INSPIRE].

[84] M. Cacciato et al., Galaxy clustering & galaxy-galaxy lensing: a promising union to constrain cosmological parameters, Mon. Not. Roy. Astron. Soc. 394 (2009) 929 [arXiv:0807.4932] [INSPIRE].
[85] C. Giocoli, M. Bartelmann, R.K. Sheth and M. Cacciato, Halo model description of the non-linear dark matter power spectrum at $k \gg 1 \text{ Mpc}^{-1}$, Mon. Not. Roy. Astron. Soc. 408 (2010) 300 [arXiv:1003.4740] [INSPIRE].

[86] R.K. Sheth and G. Tormen, Large scale bias and the peak background split, Mon. Not. Roy. Astron. Soc. 308 (1999) 119 [astro-ph/9901122] [INSPIRE].

[87] R.K. Sheth and G. Tormen, An excursion set model of hierarchical clustering: ellipsoidal collapse and the moving barrier, Mon. Not. Roy. Astron. Soc. 329 (2002) 61 [astro-ph/0105113] [INSPIRE].

[88] J.S. Bullock et al., Profiles of dark haloes. Evolution, scatter, and environment, Mon. Not. Roy. Astron. Soc. 321 (2001) 559 [astro-ph/9908159] [INSPIRE].

[89] B. Bose et al., On the road to per cent accuracy IV: ReACT — computing the non-linear power spectrum beyond $\Lambda$CDM, Mon. Not. Roy. Astron. Soc. 498 (2020) 4650 [arXiv:2005.12184] [INSPIRE].

[90] VIRGO CONSORTIUM collaboration, Stable clustering, the halo model and nonlinear cosmological power spectra, Mon. Not. Roy. Astron. Soc. 341 (2003) 1311 [astro-ph/0207664] [INSPIRE].

[91] R. Takahashi et al., Revising the halofit model for the nonlinear matter power spectrum, Astrophys. J. 761 (2012) 152 [arXiv:1208.2701] [INSPIRE].

[92] K. Heitmann et al., The Coyote universe II: cosmological models and precision emulation of the nonlinear matter power spectrum, Astrophys. J. 705 (2009) 156 [arXiv:0902.0429] [INSPIRE].

[93] T. Nishimichi et al., Dark quest. I. Fast and accurate emulation of halo clustering statistics and its application to galaxy clustering, Astrophys. J. 884 (2019) 29 [arXiv:1811.09504] [INSPIRE].

[94] S.K. Giri and A. Schneider, Emulation of baryonic effects on the matter power spectrum and constraints from galaxy cluster data, JCAP 12 (2021) 046 [arXiv:2108.08863] [INSPIRE].

[95] G. Aricò et al., The BACCO simulation project: a baryonification emulator with neural networks, Mon. Not. Roy. Astron. Soc. 506 (2021) 4070 [arXiv:2011.15018] [INSPIRE].

[96] K. Heitmann et al., The Mira-Titan universe: precision predictions for dark energy surveys, Astrophys. J. 820 (2016) 108 [arXiv:1508.02654] [INSPIRE].

[97] S. Habib et al., HACC: simulating sky surveys on state-of-the-art supercomputing architectures, New Astron. 42 (2016) 49 [arXiv:1410.2805] [INSPIRE].

[98] PLANCK collaboration, Planck 2015 results. XIII. Cosmological parameters, Astron. Astrophys. 594 (2016) A13 [arXiv:1502.01589] [INSPIRE].

[99] R.E. Angulo and S.D.M. White, One simulation to fit them all — changing the background parameters of a cosmological $N$-body simulation, Mon. Not. Roy. Astron. Soc. 405 (2010) 143 [arXiv:0912.4277] [INSPIRE].

[100] M. Archidiacono, S. Hannestad and J. Lesgourgues, What will it take to measure individual neutrino mass states using cosmology?, JCAP 09 (2020) 021 [arXiv:2003.03354] [INSPIRE].

[101] D. Blas, J. Lesgourgues and T. Tram, The Cosmic Linear Anisotropy Solving System (CLASS) II: approximation schemes, JCAP 07 (2011) 034 [arXiv:1104.2933] [INSPIRE].

[102] M. Zennaro et al., Initial conditions for accurate $N$-body simulations of massive neutrino cosmologies, Mon. Not. Roy. Astron. Soc. 466 (2017) 3244 [arXiv:1605.05283] [INSPIRE].

[103] R.E. Angulo and A. Pontzen, Cosmological $N$-body simulations with suppressed variance, Mon. Not. Roy. Astron. Soc. 462 (2016) L1 [arXiv:1603.05253] [INSPIRE].

[104] Y.B. Zeldovich, Gravitational instability: an approximate theory for large density perturbations, Astron. Astrophys. 5 (1970) 84 [INSPIRE].
[105] M. Davis, F.J. Summers and D. Schlegel, Large scale structure in a universe with mixed hot and cold dark matter, *Nature* **359** (1992) 393 [INSPIRE].

[106] A. Klypin, J. Holtzman, J. Primack and E. Regos, *Structure formation with cold plus hot dark matter*, *Astrophys. J.* **416** (1993) 1 [astro-ph/9305011] [INSPIRE].

[107] E. Sefusatti, M. Crocce, R. Scoccimarro and H. Couchman, *Accurate estimators of correlation functions in Fourier space*, *Mon. Not. Roy. Astron. Soc.* **460** (2016) 3624 [arXiv:1512.07295] [INSPIRE].

[108] G. Tormen, *The assembly of matter in galaxy clusters*, *Mon. Not. Roy. Astron. Soc.* **297** (1998) 648 [astro-ph/9802290] [INSPIRE].

[109] G. Tormen and L. Moscardini, *Properties of cluster satellites in hydrodynamical simulations*, *Mon. Not. Roy. Astron. Soc.* **350** (2004) 1397 [astro-ph/0304375] [INSPIRE].

[110] C. Giocoli, G. Tormen and F.C. Bosch, *The population of dark matter subhaloes: mass functions and average mass loss rates*, *Mon. Not. Roy. Astron. Soc.* **386** (2008) 2135 [arXiv:0712.1563] [INSPIRE].

[111] G. Despali et al., *The universality of the virial halo mass function and models for non-universality of other halo definitions*, *Mon. Not. Roy. Astron. Soc.* **456** (2016) 2486 [arXiv:1507.05627] [INSPIRE].

[112] W.A. Watson et al., *The halo mass function through the cosmic ages*, *Mon. Not. Roy. Astron. Soc.* **433** (2013) 1230 [arXiv:1212.0095] [INSPIRE].

[113] J.L. Tinker et al., *The large scale bias of dark matter halos: numerical calibration and model tests*, *Astrophys. J.* **724** (2010) 878 [arXiv:1001.3162] [INSPIRE].

[114] E. Massara, F. Villaescusa-Navarro and M. Viel, *The halo model in a massive neutrino cosmology*, *JCAP* **12** (2014) 053 [arXiv:1410.6813] [INSPIRE].

[115] A. Schneider et al., *Matter power spectrum and the challenge of percent accuracy*, *JCAP* **04** (2016) 047 [arXiv:1503.05920] [INSPIRE].

[116] G. Parimbelli et al., *DEMNUni: comparing nonlinear power spectra prescriptions in the presence of massive neutrinos and dynamical dark energy*, *JCAP* **11** (2022) 041 [arXiv:2207.13677] [INSPIRE].

[117] S. Saito, M. Takada and A. Taruya, *Impact of massive neutrinos on nonlinear matter power spectrum*, *Phys. Rev. Lett.* **100** (2008) 191301 [arXiv:0801.0607] [INSPIRE].

[118] Y.Y.Y. Wong, *Higher order corrections to the large scale matter power spectrum in the presence of massive neutrinos*, *JCAP* **10** (2008) 035 [arXiv:0809.0693] [INSPIRE].

[119] B. Audren et al., *Neutrino masses and cosmological parameters from a Euclid-like survey: Markov chain Monte Carlo forecasts including theoretical errors*, *JCAP* **01** (2013) 026 [arXiv:1210.2194] [INSPIRE].

[120] R. de Belsunce and L. Senatore, *Tree-level bispectrum in the effective field theory of large-scale structure extended to massive neutrinos*, *JCAP* **02** (2019) 038 [arXiv:1804.06849] [INSPIRE].

[121] R. Ruggeri, E. Castorina, C. Carbone and E. Sefusatti, *DEMNUni: massive neutrinos and the bispectrum of large scale structures*, *JCAP* **03** (2018) 003 [arXiv:1712.02334] [INSPIRE].

[122] M. Costanzi et al., *Cosmology with massive neutrinos III: the halo mass function and an application to galaxy clusters*, *JCAP* **12** (2013) 012 [arXiv:1311.1514] [INSPIRE].

[123] F. Hassani, J. Adamek, R. Durrer and M. Kunz, *Biased tracers as a probe of beyond-ΛCDM cosmologies*, *Astron. Astrophys.* **668** (2022) A56 [arXiv:2206.14179] [INSPIRE].

[124] M. Martinelli et al., *Euclid: impact of non-linear and baryonic feedback prescriptions on cosmological parameter estimation from weak lensing cosmic shear*, *Astron. Astrophys.* **649** (2021) A100 [arXiv:2010.12382] [INSPIRE].