# Euclid preparation. Simulations and nonlinearities beyond ΛCDM. 2. Results from non-standard simulations

G Rácz, M.-A Breton, B Fiorini, A.M.C Le Brun, H.-A Winther, Z Sakr, L Pizzuti, A Ragagnin, T Gayoux, E Altamura, et al.

► To cite this version:

G Rácz, M.-A Breton, B Fiorini, A.M.C Le Brun, H.-A Winther, et al., Euclid preparation. Simulations and nonlinearities beyond ΛCDM. 2. Results from non-standard simulations. Astron.Astrophys., 2025, 695, pp.A232. 10.1051/0004-6361/202452185 . hal-04701125

HAL Id: hal-04701125

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

Submitted on 19 Jul 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Euclid preparation

LXIII. Simulations and non-linearities beyond Lambda cold dark matter.
2. Results from non-standard simulations

Euclid Collaboration: G. Rácz¹,²,*, M.-A. Breton³,⁴,⁵, B. Fiorini⁶,⁷, A. M. C. Le Brun⁵,⁸, H.-A. Winther⁹, Z. Sakr¹⁰,¹¹,¹², L. Pizzuti¹³, A. Ragagnin¹⁴,¹⁵,¹⁶,¹⁷, T. Gayoux⁵, E. Altamura¹⁸, E. Carella¹⁹,²⁰, K. Pardede²¹,²²,²³,¹⁵, G. Verza²⁴,²⁵, K. Koyama⁶, M. Baldi²⁶,¹⁴,²⁷, A. Pourtsidou²⁸,²⁹, F. Vernizzi³⁰, A. G. Adame³¹,³²,³³, J. Adamek³⁴, S. Avila³⁵, C. Carbone¹⁹, G. Despali¹⁶,¹⁴,²⁷, C. Giocoli¹⁴,³⁶, C. Hernández-Aguayo³⁷, F. Hassani⁹, M. Kunz³⁸, B. Li³⁹, Y. Rasera⁵,⁴⁰, G. Yepes³¹,³³, V. Gonzalez-Perez³¹, P.-S. Corasaniti⁵, J. García-Bellido³², N. Hamaus⁴¹,⁴², A. Kiessling¹, M. Marinucci⁴³,⁴⁴, C. Moretti²²,¹⁷,⁴⁵,¹⁵,⁴⁶, D. F. Mota⁹, L. Piga⁴⁷,²¹,¹⁹, A. Pisani⁴⁸,⁴⁹,²⁵,⁵⁰, I. Szapudi⁵¹, P. Tallada-Crespi⁵²,⁵³, N. Aghanim⁵⁴, S. Andreon⁵⁵, C. Baccigalupi¹⁵,⁴⁵,⁴⁶,²², S. Bardelli¹⁴, D. Bonino⁵⁶, E. Branchini⁵⁷,⁵⁸,⁵⁹, M. Brescia⁵⁹,⁶⁰,⁶¹, J. Brinchmann⁶²,⁶³, S. Camera⁶⁴,⁶⁵,⁵⁶, V. Capobianco⁵⁶, V. F. Cardone⁶⁶,⁶⁷, J. Carretero⁵²,⁵³, S. Casas⁶⁸, M. Castellano⁶⁶, G. Castignani¹⁴, S. Cavuoti⁶⁰,⁶¹, A. Cimatti⁶⁹, C. Colodro-Conde⁷⁰, G. Congedo²⁸, C. J. Conselice¹⁸, L. Conversi⁷¹,⁷², Y. Copin⁷³, F. Courbin⁷⁴, H. M. Courtois⁷⁵, A. Da Silva⁷⁶,⁷⁷, H. Degaudenzi⁷⁸, G. De Lucia⁴⁵, M. Douspis⁵⁴, F. Dubath⁷⁸, C. A. J. Duncan¹⁸, X. Dupac⁷², S. Dusini⁴⁴, A. Ealet⁷³, M. Farina⁷⁹, S. Farrens⁸⁰, S. Ferriol⁷³, P. Fosalba⁸¹,⁴, M. Frailis⁴⁵, E. Franceschi¹⁴, M. Fumana¹⁹, S. Galeotta⁴⁵, B. Gillis²⁸, P. Gómez-Alvarez⁸²,⁷², A. Grazian⁸³, F. Grupp⁸⁴,⁴¹, S. V. H. Haugan⁹, W. Holmes¹, F. Hormuth⁸⁵, A. Hornstrup⁸⁶,⁸⁷, S. Ilić⁸⁸,¹¹, K. Jahnke⁸⁹, M. Jhabvala⁹⁰, B. Joachimi⁹¹, E. Keihänen⁹², S. Kermiche⁴⁸, M. Kilbinger⁸⁰, T. Kitching⁹³, B. Kubik⁷³, H. Kurki-Suonio²,⁹⁴, P. B. Lilje⁹, V. Lindholm²,⁹⁴, I. Lloro⁹⁵, G. Mainetti⁹⁶, E. Maiorano¹⁴, O. Mansutti⁴⁵, O. Marggraf⁹⁷, K. Markovic¹, M. Martinelli⁶⁶,⁶⁷, N. Martinet⁹⁸, F. Marulli¹⁶,¹⁴,²⁷, R. Massey³⁹, E. Medinaceli¹⁴, S. Mei⁹⁹, Y. Mellier¹⁰⁰, M. Meneghetti¹⁴,²⁷, G. Meylan⁷⁴, M. Moresco¹⁶,¹⁴, L. Moscardini¹⁶,¹⁴,²⁷, S.-M. Niemi¹⁰¹, C. Padilla³⁵, S. Paltani⁷⁸, F. Pasian⁴⁵, K. Pedersen¹⁰², W. J. Percival¹⁰³,¹⁰⁴,¹⁰⁵, V. Pettorino¹⁰¹, S. Pires⁸⁰, G. Polenta¹⁰⁶, M. Poncet¹⁰⁷, L. A. Popa¹⁰⁸, F. Raison⁸⁴, R. Rebolo⁷⁰,¹⁰⁹, A. Renzi⁴³,⁴⁴, J. Rhodes¹, G. Riccio⁶⁰, E. Romelli⁴⁵, M. Roncarelli¹⁴, R. Saglia⁴¹,⁸⁴, J.-C. Salvignol¹⁰¹, A. G. Sánchez⁸⁴, D. Sapone¹¹⁰, B. Sartoris⁴¹,⁴⁵, M. Schirmer⁸⁹, T. Schrabback¹¹¹, A. Secroun⁴⁸, G. Seidel⁸⁹, S. Serrano⁸¹,¹¹²,³, C. Sirignano⁴³,⁴⁴, G. Sirri²⁷, L. Stanco⁴⁴, J. Steinwagner⁸⁴, A. N. Taylor²⁸, I. Tereno⁷⁶,¹¹³, R. Toledo-Moreo¹¹⁴, F. Torradeflot⁵³,⁵², I. Tutusaus¹¹, L. Valenziano¹⁴,¹¹⁵, T. Vassallo⁴¹,⁴⁵, G. Verdoes Kleijn¹¹⁶, Y. Wang¹¹⁷, J. Weller⁴¹,⁸⁴, E. Zucca¹⁴, A. Biviano⁴⁵,¹⁵, A. Boucaud⁹⁹, E. Bozzo⁷⁸, C. Burigana¹¹⁸,¹¹⁵, M. Calabrese¹¹⁹,¹⁹, D. Di Ferdinando²⁷, J. A. Escartin Vigo⁸⁴, G. Fabbian¹²⁰,¹²¹,²⁵, F. Finelli¹⁴,¹¹⁵, J. Gracia-Carpio⁸⁴, S. Matthew²⁸, N. Mauri⁶⁹,²⁷, A. Pezzotta⁸⁴, M. Pöntinen², C. Porciani⁹⁷, V. Scottez¹⁰⁰,¹²², M. Tenti²⁷, M. Viel¹⁵,⁴⁵,²²,⁴⁶,¹⁷, M. Wiesmann⁹, Y. Akrami³²,¹²³, V. Allevato⁶⁰, S. Anselmi⁴⁴,⁴³,⁵, M. Archidiacono²⁰,¹²⁴, F. Atrio-Barandela¹²⁵, A. Balaguera-Antolinez⁷⁰,¹⁰⁹, M. Ballardini¹²⁶,¹⁴,¹²⁷, D. Bertacca⁴³,⁸³,⁴⁴, L. Blot¹²⁸,⁵, S. Borgani¹²⁹,¹⁵,⁴⁵,⁴⁶, S. Bruton¹³⁰, R. Cabanac¹¹, A. Calabro⁶⁶, B. Camacho Quevedo⁸¹,³, A. Cappi¹⁴,¹³¹, F. Caro⁶⁶, C. S. Carvalho¹¹³, T. Castro⁴⁵,⁴⁶,¹⁵,¹⁷, K. C. Chambers⁵¹, S. Contarini⁸⁴, A. R. Cooray¹³², B. De Caro¹⁹, S. de la Torre⁹⁸, G. Desprez¹³³, A. Díaz-Sánchez¹³⁴, J. J. Diaz¹³⁵, S. Di Domizio⁵⁷,⁵⁸, H. Dole⁵⁴, S. Escoffier⁴⁸, A. G. Ferrari⁶⁹,²⁷, P. G. Ferreira¹³⁶, I. Ferrero⁹, A. Fontana⁶⁶, F. Fornari¹¹⁵, L. Gabarra¹³⁶, K. Ganga⁹⁹, T. Gasparetto⁴⁵, E. Gaztanaga³,⁸¹,⁶, F. Giacomini²⁷, F. Gianotti¹⁴, G. Gozalias¹³⁷,²⁹, C. M. Gutierrez¹³⁸, A. Hall²⁸, H. Hildebrandt¹³⁹, J. Hjorth¹⁴⁰, A. Jimenez Muñoz¹⁴¹, J. J. E. Kajava¹⁴²,¹⁴³, V. Kansal¹⁴⁴,¹⁴⁵, D. Karagiannis⁷,¹⁴⁶, C. C. Kirkpatrick⁹², F. Lacasa¹⁴⁷,⁵⁴, J. Le Graet⁴⁸, L. Legrand¹⁴⁸, J. Lesgourgues⁶⁸, T. I. Liaudat¹⁴⁹, A. Loureiro¹⁵⁰,¹⁵¹, J. Macias-Perez¹⁴¹, G. Maggio⁴⁵, M. Magliocchetti⁷⁹, F. Mannucci¹⁵²,

<footnote>* Corresponding author; gabor.racz@helsinki.fi</footnote>
R. Maoli¹⁵³,⁶⁶, C. J. A. P. Martins¹⁵⁴,⁶², L. Maurin⁵⁴, R. B. Metcalf¹⁶,¹⁴, M. Miluzio⁷²,¹⁵⁵,
P. Monaco¹²⁹,⁴⁵,⁴⁶,¹⁵, A. Montoro³,⁸¹, A. Mora¹⁵⁶, G. Morgante¹⁴, S. Nadathur⁶, L. Patrizii²⁷, V. Popa¹⁰⁸,
D. Potter³⁴, P. Reimberg¹⁰⁰, I. Risso¹⁵⁷, P.-F. Rocci⁵⁴, M. Sahlén¹⁵⁸, A. Schneider³⁴, M. Sereno¹⁴,²⁷,
A. Silvestri¹⁵⁹, A. Spurio Mancini¹⁶⁰,⁹³, J. Stadel³⁴, K. Tanidis¹³⁶, C. Tao⁴⁸, N. Tessore⁹¹, G. Testera⁵⁸,
R. Teyssier⁴⁹, S. Toft⁸⁷,¹⁶¹,¹⁶², S. Tosi⁵⁷,⁵⁸, A. Troja⁴³,⁴⁴, M. Tucci⁷⁸, C. Valieri²⁷, J. Valiviita²,⁹⁴,
D. Vergani¹⁴, P. Vielzeuf⁴⁸, N. A. Walton¹²⁰

(Affiliations can be found after the references)

Received 9 September 2024 / Accepted 11 November 2024

ABSTRACT

The *Euclid* mission will measure cosmological parameters with unprecedented precision. To distinguish between cosmological models, it is essential to generate realistic mock observables from cosmological simulations that were run in both the standard Λ-cold-dark-matter (ΛCDM) paradigm and in many non-standard models beyond ΛCDM. We present the scientific results from a suite of cosmological N-body simulations using non-standard models including dynamical dark energy, *k*-essence, interacting dark energy, modified gravity, massive neutrinos, and primordial non-Gaussianities. We investigate how these models affect the large-scale-structure formation and evolution in addition to providing synthetic observables that can be used to test and constrain these models with *Euclid* data. We developed a custom pipeline based on the *Rockstar* halo finder and the *nbodykit* large-scale structure toolkit to analyse the particle output of non-standard simulations and generate mock observables such as halo and void catalogues, mass density fields, and power spectra in a consistent way. We compare these observables with those from the standard ΛCDM model and quantify the deviations. We find that non-standard cosmological models can leave large imprints on the synthetic observables that we have generated. Our results demonstrate that non-standard cosmological N-body simulations provide valuable insights into the physics of dark energy and dark matter, which is essential to maximising the scientific return of *Euclid*.

**Key words.** methods: numerical – cosmology: theory – dark matter – dark energy – large-scale structure of Universe

# 1. Introduction

The concordance Λ-cold-dark-matter (ΛCDM) model is the simplest cosmological scenario that accounts for the cosmological observations thus far available. It is based on the assumption that in addition to baryonic matter and radiation, the Universe is filled with two invisible components: an exotic form of matter, dubbed dark energy and described by a cosmological constant (Λ) in Einstein's equations of general relativity, and a cold-dark-matter (CDM) component that is non-relativistic and only interacts through gravity. In this scenario, dark matter is primarily responsible for fostering the formation of the visible structures we observe today, while dark energy drives the accelerated expansion of the Universe at late times. This model has been remarkably successful in explaining a variety of cosmological observations, such as the Hubble diagram from luminosity distance measurements of type Ia supernovae (Riess et al. 1998; Perlmutter et al. 1999), temperature and polarisation anisotropy angular power spectra of the cosmic microwave background (CMB, de Bernardis et al. 2000; Spergel et al. 2003; Kovac et al. 2002; Planck Collaboration VI 2020), the galaxy power spectrum of the large-scale structure (LSS, Efstathiou et al. 2002; Colless et al. 2003; Tegmark et al. 2004, 2006), and the presence of baryonic acoustic oscillations (BAO) in the LSS (Eisenstein et al. 2005; Cole et al. 2005). Despite the great success of the ΛCDM model, the physical origin of dark energy and dark matter remains unknown. Unveiling the nature of these dark components is the primary motivation for many investigations in modern cosmology.

In the last decade, multiple tensions among different types of cosmological observations have emerged. As an example, while CMB measurements indicate a value of the Hubble constant of $H_0 = 67.7 \pm 0.4 \text{ km s}^{-1} \text{ Mpc}^{-1}$ (Planck Collaboration VI 2020), local measurements, often based on the observations of supernovae in nearby galaxies, suggest a higher value of $H_0 = 73.0 \pm 1.0 \text{ km s}^{-1} \text{ Mpc}^{-1}$ (Riess et al. 2022). This $5\sigma$ discrepancy

is called the Hubble tension. A similar tension has been identified in the $S_8 = \sigma_8 \sqrt{\Omega_m}/0.3$ parameter, which combines the amplitude of linear-matter density fluctuations on the $8 h^{-1} \text{ Mpc}$ scale, $\sigma_8$, and the cosmic matter density, $\Omega_m$. Measurements derived from the CMB (Planck Collaboration VI 2020) appear to yield a value of $S_8 = 2.9\sigma$ higher than that obtained from observations of the LSS (Joseph et al. 2023), such as measurements of the clustering of galaxies and weak gravitational lensing (Li et al. 2023; Abbott et al. 2022). Such tensions may result from systematic errors yet to be identified in the data. Alternatively, they may be a manifestation of the limits of the ΛCDM model, since modifications to the standard cosmological model can provide a solution to these tensions (see e.g. Martinelli & Tutusaus 2019; Di Valentino et al. 2021).

Ongoing and upcoming Stage-IV surveys, such as *Euclid* (Euclid Collaboration: Mellier et al. 2025; Laureijs et al. 2011; Euclid Collaboration: Scamarella et al. 2022), Dark Energy Spectroscopic Instrument (hereafter DESI, DESI Collaboration 2016), Vera C. Rubin Observatory Legacy Survey of Space and Time (LSST, Ivezić et al. 2019), Spectro-Photometer for the History of the Universe, Epoch of Reionization, and the Ices Explorer (SPHEREx, Doré et al. 2014), and the *Nancy Grace Roman* Space Telescope (Spergel et al. 2015), will collect unprecedented amounts of data on the LSS, which will enable detailed assessments of the Hubble and $S_8$ tensions in addition to shedding new light on the nature of the invisible components in the Universe.

*Euclid* is a space mission led by the European Space Agency (ESA) with contributions from the National Aeronautics and Space Administration (NASA), aiming to study the nature and evolution of the dark Universe. The survey uses a 1.2-m-diameter telescope and two instruments, a visible-wavelength camera, and a near-infrared camera/spectrometer, to observe billions of galaxies over more than a third of the sky in optical and near-infrared wavelengths. *Euclid* measures the shapes (Euclid Collaboration: Bretonnière et al. 2022; 2023;
Euclid Collaboration: Merlin et al. 2023d) and redshifts (Euclid Collaboration: Desprez et al. 2020; Euclid Collaboration: Ilbert et al. 2021) of galaxies, in order to determine the weak gravitational lensing and clustering of galaxies, covering a period of cosmic history over which dark energy accelerated the expansion of the Universe. These measurements will provide detailed insights into the properties of dark energy, dark matter, and gravity by probing the expansion history of the Universe and the growth rate of structures over time (Martinelli et al. 2021; Nesseris et al. 2022; Euclid Collaboration: Castro et al. 2023). Euclid was launched on July 1 2023 and is designed to operate for six years. The survey will provide unprecedented constraints on cosmological parameters and tests of fundamental physics, as well as a rich catalogue of legacy data that can be used for a wide range of astrophysical research. The mission data will be publicly released within two years of acquisition. Euclid is one of the most ambitious and exciting space missions in the field of cosmology and will enable a thorough validation of a broad range of cosmological models.

Euclid observations will provide precise measurements of the clustering of matter over a wide range of scales, where effects due to the late-time non-linear gravitational collapse of matter need to be taken into account. A key tool in the preparation of the cosmological analyses and the interpretation of the Euclid data is the use of cosmological N-body simulations, which can follow the non-linear evolution of matter clustering. This is a numerical technique that calculates the evolution of the matter density field under the effect of gravity across cosmic time and predicts the LSS of the Universe for a given cosmological model (Press & Schechter 1974; Zeldovich 1978; Klypin & Shandarin 1983; Appel 1985; Potter et al. 2017; Angulo & Hahn 2022). In this method, the matter density field is sampled with discrete N-body particles, whose equations of motion are solved in the Newtonian limit in an expanding Friedmann–Lemaître–Robertson–Walker (FLRW) Universe. These simulations enable the study of the formation and growth of cosmic structures from linear to non-linear scales and predict the distribution of matter in galaxy clusters, filaments, and voids, for a range of cosmological models and parameters (Klypin et al. 2003; Dolag et al. 2004; Alimi et al. 2010; Li et al. 2012; Puchwein et al. 2013; Baldi & Simpson 2015), as well as initial conditions (Dalal et al. 2008). The cosmological models beyond the standard-ΛCDM paradigm are expected to have left imprints that should be detectable in the Euclid observables, such as the redshift-space power spectra of galaxies or the void-size functions.

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

1. Numerical methods and validation (Euclid Collaboration: Adamek et al. 2025).

2. Results from non-standard simulations (this work).

3. Constraints on $f(R)$ models from the photometric primary probes (Euclid Collaboration: Koyama et al. 2024).

4. Cosmological constraints on non-standard cosmologies from simulated Euclid probes (D'Amico et al. in prep.).

(For further details, see our companion papers.) In this work, we consistently analyse large numbers of N-body simulations over a wide range of non-standard cosmological scenarios, to generate catalogues of synthetic observables for Euclid. This analysis is achieved using a pipeline that was specifically written for that task. We calculate reconstructed density fields, halo and void catalogues, halo mass functions, dark matter, and halo power spectra in real and redshift space, as well as halo bias functions. The paper is organised as follows: in Sect. 2, we introduce the analysed non-standard models; then, in Sect. 3, we present an overview of the analysed cosmological N-body sim-

ulations. In Sect. 4, we describe the analysis pipeline and the calculated quantities. We demonstrate the imprints of the non-standard models in the computed observables in Sect. 5 and finally, we summarise our results in Sect. 6.

## 2. Cosmological models beyond the standard ΛCDM paradigm

To address the tensions and anomalies in the ΛCDM model, various non-standard cosmological models have been proposed that extend or modify the standard model in different ways. Some examples of non-standard cosmological models are dark energy models, such as quintessence and phantom energy, modified-gravity theories, such as $f(R)$ gravity, and massive-neutrino models, such as sterile neutrinos and self-interacting neutrinos. These models introduce new degrees of freedom or new mechanisms that can affect the dynamics and observables of the Universe at different scales and epochs. In this section, we discuss the main features, motivations, and challenges of these non-standard cosmological models.

### 2.1. Dark energy models

#### 2.1.1. wCDM

A simple generalisation of the cosmological constant assumes that dark energy is a fluid with a constant equation-of-state, $w \equiv p_{de}/(\rho_{de} c^2)$, where $p_{de}$ and $\rho_{de}$ are, respectively, the pressure and density of the fluid, and $c$ is the speed of light. To trigger an accelerated phase of cosmic expansion, the dark energy equation-of-state parameter must be $w < -1/3$. The ΛCDM model corresponds to the $w = -1$ specific case, while $w < -1$ corresponds to so-called phantom dark energy models (Caldwell et al. 2003), though such values may also result from an unaccounted interaction between dark energy and dark matter (Das et al. 2006).

#### 2.1.2. Dynamical dark energy

The dark energy equation-of-state could be a function of redshift. Chevallier, Polarski (Chevallier & Polarski 2001) and Linde (Linde 2003) proposed a simple parameterisation of

$$w_{de}(z) = w_0 + w_a \frac{z}{1+z} = w_0 + w_a (1-a), \quad (1)$$

where the $w_0$ parameter represents the value of the equation-of-state at the present time, and $w_a$ defines the rate of change with redshift. This model is also called the CPL parametrisation of dark energy, after the initials of the authors who proposed it.

This dark energy parametrisation is a fitting function of a general $w_{de}(z)$ around $z = 0$, assuming that $w(z)$ is smooth and slowly changing with the scale factor. As a consequence, this model can closely follow the expansion history of a wide range of other models with $w_{de}(z)$ at late times. Despite its simple form, it shows a wide range of interesting properties (Linde 2008; Linden & Virey 2008). The cosmological constant corresponds to $w_0 = -1$ and $w_a = 0$ in the CPL parametrisation.

#### 2.1.3. K-essence

The k-essence model is characterised by an action for the scalar field of the following form:

$$S = \int d^4 x \sqrt{-g} p(\phi, X), \quad (2)$$
where $X = (1/2)g^{\mu\nu}\nabla_{\mu}\phi\nabla_{\nu}\phi$. The energy density of the scalar field is given by

$$u_{\phi} = 2X \frac{dp}{dX} - p, \qquad (3)$$

and the pressure is $p_{\phi} = p(\phi, X)$. This pressure gives an effective fluid equation-of-state parameter as

$$w_{\phi} = \frac{p_{\phi}}{u_{\phi}} = -\frac{p}{p - 2Xp_{,X}}, \qquad (4)$$

where the subscript ,X indicates a derivative with respect to X, and a dimensionless speed-of-sound parameter for the k-essence fluctuations as

$$c_s^2 = \frac{p_{,X}}{p_{,X} + 2Xp_{,XX}}. \qquad (5)$$

The k-essence field satisfies the continuity equation

$$\dot{u}_{\phi} = -3H(u_{\phi} + p_{\phi}), \qquad (6)$$

which results in the scalar equation of motion

$$G^{\mu\nu}\nabla_{\mu}\nabla_{\nu}\phi + 2X \frac{\partial^2 p}{\partial X \partial \phi} - \frac{\partial p}{\partial \phi} = 0, \qquad (7)$$

where

$$G^{\mu\nu} = \frac{\partial p}{\partial X} g^{\mu\nu} + \frac{\partial^2 p}{\partial X^2} \nabla^{\mu} \phi \nabla^{\nu} \phi. \qquad (8)$$

The k-essence was first proposed by Armendariz-Picon et al. (2000, 2001), who showed that there exist tracking attractor solutions to the equation of motion during the radiation and matter-dominated eras of the Universe, and that with a suitably chosen $p$, the scalar can have an appropriate equation of state that allows it to act as dark energy for the background accelerated expansion. In addition, whenever the kinetic terms for the scalar field are not linear in $X$, the speed of sound of the fluctuations differs from unity, allowing for the clustering of the dark energy field at sub-horizon scales, which should be modelled at the perturbations level.

### 2.1.4. Interacting dark energy

In the interacting dark energy (IDE) models (Amendola 2000; Farrar & Peebles 2004; Baldi et al. 2010), dark energy and cold dark matter are allowed to interact through an exchange of energy-momentum in order to keep the total stress-energy tensor, $T_{\mu\nu}$, conserved:

$$\nabla_{\mu} T_{\nu}^{(c)\mu} = C_{\nu}(\phi) = -\nabla_{\mu} T_{\nu}^{(\phi)\mu}, \qquad (9)$$

where $C_{\nu}(\phi)$ is a conformal coupling function expressed in the form:

$$C_{\nu}(\phi) = \kappa\beta(\phi)u_{c}\nabla_{\nu}\phi, \qquad (10)$$

where $\kappa \equiv \frac{8\pi G_N}{c^2}$, $G_N$ is Newton's gravitational constant, $u_c$ is the cold dark matter energy density in the IDE model<fnref n="1" />, and $\beta(\phi)$

is a coupling function. The dark energy scalar field, $\phi$, has an intrinsic energy density and pressure given by

$$u_{\phi} = \frac{1}{2} g^{\mu\nu} \partial_{\mu} \phi \partial_{\nu} \phi + V(\phi), \qquad (11)$$

$$p_{\phi} = \frac{1}{2} g^{\mu\nu} \partial_{\mu} \phi \partial_{\nu} \phi - V(\phi), \qquad (12)$$

where $V(\phi)$ is a self-interaction potential. The conservation equations then translate into the following set of background-dynamic equations under the assumption of a constant coupling function $\beta(\phi) = \beta$:

$$\ddot{\phi} + 3H\dot{\phi} + \frac{dV}{d\phi} = \kappa\beta u_c, \qquad (13)$$

$$\dot{u}_c + 3Hu_c = -\kappa\beta u_c. \qquad (14)$$

In the standard approach, a theoretically motivated analytical form for the self-interaction potential function, $V(\phi)$, is chosen. However, the simulations that are considered in the present work implement the alternative approach proposed by Barros (2019) that consists of imposing a standard $\Lambda$CDM background expansion history by setting

$$H^2 = H_{\Lambda\text{CDM}}^2, \qquad (15)$$

where $H_{\Lambda\text{CDM}}$ is the standard Hubble function defined by

$$H_{\Lambda\text{CDM}}^2 = \frac{8\pi G_N}{3}(\rho_r + \rho_b + \rho_{\text{CDM}} + \rho_{\Lambda}), \qquad (16)$$

where $\rho_r$, $\rho_b$, $\rho_{\text{CDM}}$, and $\rho_{\Lambda}$ are the mass densities of the radiation, baryon, CDM, and $\Lambda$ components of the background $\Lambda$CDM model. This will determine an effective potential, $V(\phi)$, according to the resulting evolution of the scalar field, $\phi$. Taking the time derivative of Eq. (15) and using the continuity Eqs. (13 and 14), one gets the scalar-field energy density and pressure as

$$u_{\phi} = \rho_{\text{CDM}} c^2 + \rho_{\Lambda} c^2 - u_c, \qquad (17)$$

$$p_{\phi} = p_{\Lambda} = -u_{\Lambda}, \qquad (18)$$

which can be combined with Eqs. (11 and 12) to obtain the dynamics of the scalar field:

$$\dot{\phi}^2 = \rho_{\text{CDM}} c^2 - u_{\phi}. \qquad (19)$$

The scalar-field potential, $V(\phi)$, can then be reconstructed using Eqs. (17 and 18) as:

$$V(\phi) = \frac{1}{2}\dot{\phi}^2 + \rho_{\Lambda} c^2, \qquad (20)$$

and taking the time derivative of Eq. (19), one can derive the scalar-field equation of motion

$$2\ddot{\phi} + 3H\dot{\phi} - \kappa\beta u_c = 0, \qquad (21)$$

which can be numerically solved for the dynamical evolution of the system. With this choice, the $\beta$ coupling remains the only free parameter of this model. Observational constraints on the model were computed in Barros et al. (2023), which found that the model can alleviate the $\sigma_8$ tension, but that CMB prefers the $\Lambda$CDM limit. In particular, they find that the CMB constrains $|\beta| \lesssim 0.02$, RSD constraints $|\beta| \lesssim 0.10$, while weak lensing data from the Kilo-Degree Survey actually prefers a non-zero value $|\beta| \sim 0.1$.

<footnote n="1">¹ Note that we choose $u_c$ here for our IDE model to better distinguish it from $\rho_{\text{CDM}}$ used right afterwards to describe the $\Lambda$CDM background evolution.</footnote>
## 2.2. Modified gravity models

### 2.2.1. nDGP gravity

The Dvali–Gabadadze–Porrati (DGP) model (Dvali et al. 2000) assumes that our Universe is described by a five-dimensional bulk, while the visible matter component is confined to the four-dimensional brane described by the Minkowski metric, $\gamma$. This model’s action is

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

where $G_5$ and $G_N$ are the five- and four-dimensional Newton’s constants, respectively, and $\mathcal{L}_m$ is the matter Lagrangian. At small scales, four-dimensional gravity is recovered due to an intrinsic Einstein-Hilbert term sourced by brane curvature causing a gravitational force that scales as $r^{-2}$, while, at large scales, the gravity behaves as a five-dimensional force. The transition between the five-dimensional modifications and the four-dimensional gravity is given by the cross-over scale, $r_c = G_5/(2G_N)$, from which we construct the dimensionless parameter $\Omega_{rc} \equiv c^2/(4r_c^2H_0^2)$. The modified Friedmann equation on the brane (Deffayet 2001) becomes

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

The model we investigate in this paper is the normal branch with the – sign (Bowcock et al. 2000) characterised by a $\Lambda$CDM background achieved by introducing an additional dark energy contribution with an appropriate equation-of-state (Schmidt 2009),

$$\rho_{de}(a) = \rho_{cr,0} (\Omega_{\Lambda} + 2 \sqrt{\Omega_{rc}} \sqrt{\Omega_{\Lambda} + \Omega_{m} a^{-3}}), \quad (24)$$

where $\rho_{cr,0}$ is the critical density. The observational constraints on the model require the cross-over scale, $r_c$, to be larger than the size of the horizon $H_0^{-1}$ today. For example, Solar System constraints require $r_cH_0 \gtrsim 1.6$ (Battat et al. 2008), and galaxy clustering in the BOSS survey constraints $r_cH_0 \gtrsim 4.5$ (Piga et al. 2023).

### 2.2.2. $f(R)$ gravity

The $f(R)$ theory of gravity (Buchdahl 1970) is characterised by the following action:

$$S = \frac{c^4}{16\pi G_N} \int d^4 x \sqrt{-g} [R + f(R)], \quad (25)$$

where $g_{\mu\nu}$ is the metric tensor and $f(R)$ is a functional form of the Ricci scalar, $R$. Here, we consider the Hu–Sawicki model (Hu & Sawicki 2007) with $n = 1$, where in the limit of $f_R = df/dR \ll 1$ we have

$$f(R) = -6\Omega_{\Lambda} \frac{H_0^2}{c^2} + |f_{R0}| \frac{\bar{R}_0^2}{R}, \quad (26)$$

where $f_{R0}$ is the free parameter of the model, $\bar{R}_0$ is the Ricci scalar evaluated at background at present time, $H_0$ is the Hubble constant, and $\Omega_{\Lambda}$ is the energy-density parameter of the cosmological constant. $|f_{R0}|$ characterises the magnitude of the deviation from $\Lambda$CDM, with smaller values corresponding to weaker departures from general relativity until we recover $\Lambda$CDM in

the limit of $f_{R0} \to 0$, but for the small $|f_{R0}|$ values still allowed by observations, the background expansion history approximates that of $\Lambda$CDM and

$$\bar{R}_0 = 3\Omega_m \frac{H_0^2}{c^2} \left( 1 + 4 \frac{\Omega_{\Lambda}}{\Omega_m} \right), \quad (27)$$

with the matter energy density parameter $\Omega_m = 1 - \Omega_{\Lambda}$. However, though the background expansion could mimic that of a cosmological-constant model, it still differs at the level of cosmological perturbations where the growth of structure is driven by a modification of gravity following the above adopted model of $f(R)$.

The observational constraints on the model parameter $|f_{R0}|$ vary from $|f_{R0}| \lesssim 10^{-6}$ in the Solar System, $|f_{R0}| \lesssim 10^{-8}$ from galaxy scales (Burrage et al. 2024) to $|f_{R0}| \lesssim 10^{-6}-10^{-4}$ from various cosmological probes (see, e.g., Fig. 28 in Koyama 2016, for a summary). The parameter values of the simulations presented in this paper are similar to the current cosmological constraints.

## 2.3. Massive relativistic neutrinos

Neutrinos are mainly characterised by two properties: their mass, $M_\nu$, and the number of neutrino species, $N_{\text{eff}}$. More in general, $N_{\text{eff}}$ parametrises the contribution of relativistic species to the background density of radiation, $\rho_r$, as

$$\rho_r = \left[ 1 + \frac{7}{8} \left( \frac{4}{11} \right)^{4/3} N_{\text{eff}} \right] \rho_{\gamma}, \quad (28)$$

where $\rho_\gamma$ is the photon background density. In the standard model, $N_{\text{eff}}$ is expected to be ~3.045 (Cielo et al. 2023) for three families of active neutrinos that thermalised in the early Universe and decoupled well before electron-positron annihilation. The calculation of $N_{\text{eff}}$ involves the complete treatment of neutrino decoupling, which incorporates non-instantaneous decoupling. A deviation from the fiducial value serves to account for the presence of non-standard neutrino features, or additional relativistic relics contributing to the energy budget (Mangano et al. 2002). Here we focus on standard neutrino families only.

In addition, oscillation experiments (Maltoni et al. 2004; Kajita 2016) showed that at least two neutrinos are massive by measuring two squared-mass differences. It can be shown that the minimum value of the neutrino mass sum is either 0.06 eV in the normal or 0.10 eV in the inverted hierarchy. This value can be well constrained through cosmological observations since neutrinos are known to impact the expansion history and suppress the clustering of cold dark matter, which can be observed in the large-scale distribution of galaxies (Sakr 2022). Neutrinos with mass $\lesssim 0.6$ eV become non-relativistic after the epoch of recombination probed by the CMB, and this mechanism allows massive neutrinos to alter the matter-radiation equality for a fixed $\Omega_m h^2$ (Lesgourgues & Pastor 2006). Massive neutrinos act as non-relativistic particles on scales $k > k_{\text{nr}} = 0.018(m_\nu/1\text{eV})^{1/2}\Omega_m^{1/2} h^{-1} \text{Mpc}$, where $k_{\text{nr}}$ is the wavenumber corresponding to the Hubble horizon size at the epoch $z_{\text{nr}}$ when the given neutrino species becomes non-relativistic following $1 + z_{\text{nr}} \approx 1900 (m_\nu / 1\text{eV})$, $\Omega_m$ is the matter density parameter, and $h = H_0 / 100 \text{ km s}^{-1} \text{Mpc}^{-1}$. The large velocity dispersion of non-relativistic neutrinos suppresses the formation of neutrino perturbations in a way that depends on $m_\nu$ and redshift $z$, leaving an
imprint on the matter power spectrum at scales $k > k_{fs}(z)$, with

$$k_{\text{fs}} = \frac{0.82H(z)}{H_0(1+z)^2} \left(\frac{m_{\nu}}{\text{leV}}\right) h \, \text{Mpc}^{-1}, \qquad (29)$$

where neutrinos cannot cluster and do not contribute to the gravitational potential wells produced by cold dark matter and baryons (Takada et al. 2006; Lesgourgues & Pastor 2006). This modifies the shape of the matter power spectrum and the correlation function on these scales.

## 2.4. Primordial non-Gaussianities

The simplest inflation models predict that primordial curvature perturbations follow a distribution that is close to Gaussian (Maldacena 2003; Creminelli & Zaldarriaga 2004). However, there are many alternative inflation models that predict certain amounts of primordial non-Gaussianity (PNG). One of the simplest cases is that of the so-called local primordial non-Gaussianities (Salopek & Bond 1990; Komatsu & Spergel 2001). For this case, the primordial potential $\phi$ is given by

$$\phi(\mathbf{x}) = \phi_G(\mathbf{x}) + f_{\text{NL}}^{\text{local}}(\phi^2(\mathbf{x}) - \langle\phi^2(\mathbf{x})\rangle), \qquad (30)$$

where $\phi_G(\mathbf{x})$ is the Gaussian potential, while $\phi$ is the non-Gaussian potential. $f_{\text{NL}}^{\text{local}}$ measures the level of deviations from Gaussianity.

The perturbations in the primordial potential produce perturbations in the density field and they are related through Poisson's equation. Therefore, in Fourier space, the density field is given by

$$\delta(k, z) = \alpha(k, z) \phi(k, z), \qquad (31)$$

where

$$\alpha(k, z) = \frac{2D(z)}{3\Omega_{\text{m}}} \frac{c^2}{H_0^2} \frac{g(0)}{g(z_{\text{rad}})} k^2 T(k), \qquad (32)$$

$T(k)$ is the transfer function normalised at $T(k \to 0) = 1$, and $D(z)$ is the growth factor normalised at $D(z = 0) = 1$. The factor $g(0)/g(z_{\text{rad}})$, where $g(z) = (1+z)D(z)$, takes into account the difference between our normalisation of $D(z)$ and the early-time normalisation where $D(z) \propto 1/(1+z)$ during matter-domination. This factor is $g(z_{\text{rad}})/g(0) \sim 1.3$, with a small dependency on the cosmology.

This type of non-Gaussianity characteristically affects the clustering of biased tracers, inducing a scale-dependent bias (Dalal et al. 2008; Slosar et al. 2008; Matarrese & Verde 2008). To linear order, the power spectrum of galaxies can be given as

$$P_{\text{t,t}}(k, z) = \left[ b_1 + \frac{b_\phi f_{\text{NL}}^{\text{local}}}{\alpha(k, z)} \right]^2 P_{\text{m,m}}(k, z), \qquad (33)$$

where $P_{\text{t,t}}(k, z)$ is the power spectrum of the tracer, $P_{\text{m,m}}(k, z)$ is the power spectrum of the matter, $b_1$ is the linear bias, and $b_\phi$ is the response of the tracer to the presence of the local-PNG. Now, $P_{\text{t,t}}(k, z)$ has a dependency with $k$ which scales as $k^{-2}$ at leading order due to the $\alpha(k, z)$ term. The $b_\phi$ is usually parametrised as

$$b_{\phi} = 2\delta_{c}(b_{1} - p). \qquad (34)$$

Although it is possible to make a theoretical prediction for $p$ (by assuming a universal mass function, $p = 1$, Dalal et al. 2008), several studies using numerical simulations have shown that the prediction may be different depending on the type of galaxy or tracer under consideration (Slosar et al. 2008; Desjacques et al. 2009; Hamaus et al. 2011; Biagetti et al. 2017; Barreira et al. 2020; Gutiérrez Adame et al. 2024).

## 3. Simulations

This section summarises the simulations used for this project and gives a very brief description of each setup. The analysed simulations followed the evolution of the matter field with discrete N-body method in the models described in Sect. 2. Baryonic and hydrodynamical effects are neglected in this paper. For a comprehensive description of each of the simulation suites, we refer the reader to the main references given in Table 1 along with the volumes, resolutions, initial redshifts, and the used order of the Lagrangian perturbation theory (LPT) during the initial-condition generation.

### 3.1. The COMPLEMENTARY simulations

The COMPLEMENTARY simulation series is a set of 4 cosmological N-body simulations in wCDM and ΛCDM cosmologies. This suite used the complementary-simulation method (Rácz et al. 2023), which is a novel technique in which cosmological N-body simulations are run in phase-shifted matching pairs. One simulation starts from a regular random Gaussian initial condition, while the second simulation has modified initial amplitudes of the Fourier modes to ensure that the average power spectrum of the pair is equal to the cosmic mean power spectrum from linear theory at the initial time. The average statistical properties of a pair of such simulations have greatly suppressed variance. In this paper, we have analysed two complementary pairs using ΛCDM and wCDM cosmologies. The ΛCDM simulation pair used the best-fit Planck2018 (Planck Collaboration VI 2020) cosmological parameters: $\Omega_{\text{m}} = 1 - \Omega_{\Lambda} = 0.3111$, $\Omega_{\text{b}} = 0.04897$, $H_0 = 67.66 \, \text{km s}^{-1} \, \text{Mpc}^{-1}$, $n_s = 0.9665$, and $\sigma_8 = 0.8102$. The wCDM pair had the following parameters: $w_0 = -1.04$, $\Omega_{\text{m}} = 0.3096$, $\Omega_{\text{b}} = 0.04899$, $\Omega_{\text{de}} = 0.6904$, $H_0 = 67.66 \, \text{km s}^{-1} \, \text{Mpc}^{-1}$, $n_s = 0.9331$, and $\sigma_8 = 0.8438$. The cosmological simulations of this series were run using the cosmological N-body code GIZMO (Hopkins 2015). All simulations in the series contained $2160^3$ dark matter particles in a $(1.5 \, h^{-1} \, \text{Gpc})^3$ volume, with $\varepsilon = 13.8 \, h^{-1} \, \text{kpc}$ softening length. The initial conditions (ICs) were generated by a modified version of the N-GenIC code (Springel 2015) by using the Zel'dovich approximation and initial linear power spectra from the Boltzmann code CAMB (Lewis & Challinor 2011). The simulations started from redshift $z_{\text{init}} = 127$, with a total of 48 output times. In this project, 31 particle snapshots were analysed in the $0.5 \le z \le 2.0$ redshift range for each simulation.

### 3.2. The DEMNUNI simulation suite

The ‘Dark Energy and Massive Neutrino Universe’ (DEMNUNI) simulations (Carbone et al. 2016; Parimbelli et al. 2022) have been produced with the aim of investigating the LSS in the presence of massive neutrinos and dynamical dark energy, and they were conceived for the non-linear analysis and modelling of different probes, including dark matter, halo, and galaxy clustering (see Castorina et al. 2015; Zennaro et al. 2018; Parimbelli et al. 2022; Gouyou Beauchamps et al. 2024), weak lensing, CMB lensing, Sunyaev-Zeldovich, and integrated Sachs-Wolfe (ISW) effects (Roncarelli et al. 2015; Carbone et al. 2016), cosmic void statistics (Kreisch et al. 2019), and cross-correlations among these probes (Cuozzo et al. 2024). The DEMNUNI simulations were run using the tree particle mesh-smoothed particle hydrodynamics (TreePM-SPH) code p-GADGET3 (Springel 2005), specifically modified as in Viel et al. (2010) to account for the presence of massive neutrinos. This modified version of
<table><caption>**Table 1.** Overview of the simulation suites analysed for this project.</caption><thead><tr><th>Name/Reference</th><th>Code</th><th>N<sub>sim</sub></th><th>L [h<sup>-1</sup> Gpc]</th><th>N<sub>DM</sub></th><th>m<sub>DM</sub> [h<sup>-1</sup> M<sub>☉</sub>]</th><th>z<sub>init</sub></th><th>LPT order</th><th>Model</th></tr></thead><tbody><tr><td rowspan="2">COMPLEMENTARY 1</td><td rowspan="2">GIZMO</td><td>2</td><td>1.5</td><td>2160<sup>3</sup></td><td>2.89 × 10<sup>10</sup></td><td rowspan="2">127</td><td rowspan="2">1LPT</td><td>ΛCDM</td></tr><tr><td>2</td><td>1.5</td><td>2160<sup>3</sup></td><td>2.88 × 10<sup>10</sup></td><td>wCDM</td></tr><tr><td rowspan="3">DEMNUNI</td><td>p-GADGET3</td><td>100</td><td>1.0</td><td>1024<sup>3</sup></td><td>8 × 10<sup>10</sup></td><td rowspan="3">99</td><td rowspan="3">1LPT</td><td>ΛCDM+m<sub>ν</sub></td></tr><tr><td>p-GADGET3</td><td>15</td><td>2.0</td><td>2048<sup>3</sup></td><td>8 × 10<sup>10</sup></td><td>CPL + m<sub>ν</sub></td></tr><tr><td>p-GADGET3</td><td>2</td><td>0.5</td><td>2048<sup>3</sup></td><td>1.25 × 10<sup>9</sup></td><td>ΛCDM+m<sub>ν</sub></td></tr><tr><td rowspan="2">RAYGAL 3</td><td rowspan="2">RAMSES</td><td>1</td><td>2.625</td><td>4096<sup>3</sup></td><td>1.88 × 10<sup>10</sup></td><td rowspan="2">49</td><td rowspan="2">2LPT</td><td>ΛCDM</td></tr><tr><td>1</td><td>2.625</td><td>4096<sup>3</sup></td><td>2.0 × 10<sup>10</sup></td><td>wCDM</td></tr><tr><td rowspan="2">ELEPHANT 4</td><td rowspan="2">ECOSMOG</td><td>11</td><td>1.024</td><td>1024<sup>3</sup></td><td>8.85 × 10<sup>10</sup></td><td rowspan="2">49</td><td rowspan="2">2LPT</td><td>ΛCDM</td></tr><tr><td>11</td><td>1.024</td><td>1024<sup>3</sup></td><td>8.85 × 10<sup>10</sup></td><td>nDGP</td></tr><tr><td rowspan="2">COLA HIRES 5</td><td rowspan="2">MG-COLA</td><td>5</td><td>1.024</td><td>2048<sup>3</sup></td><td>1.07 × 10<sup>10</sup></td><td rowspan="2">127</td><td rowspan="2">2LPT</td><td>nDGP</td></tr><tr><td>2</td><td>1.024</td><td>2048<sup>3</sup></td><td>1.07 × 10<sup>10</sup></td><td>ΛCDM</td></tr><tr><td rowspan="2">DUSTGRAIN 6</td><td rowspan="2">MG-GADGET</td><td>3</td><td>2.0</td><td>2048<sup>3</sup></td><td>8.27 × 10<sup>10</sup></td><td rowspan="2">99</td><td rowspan="2">1LPT</td><td>nDGP</td></tr><tr><td>11</td><td>0.75</td><td>768<sup>3</sup></td><td>8.1 × 10<sup>10</sup></td><td>f(R) + m<sub>ν</sub></td></tr><tr><td>CIDER 7</td><td>c-GADGET</td><td>4</td><td>1.0</td><td>1024<sup>3</sup></td><td>8.1 × 10<sup>10</sup></td><td>99</td><td>2LPT</td><td>cDE</td></tr><tr><td>DAKAR (1&2) 8</td><td>c-GADGET</td><td>5</td><td>1.0</td><td>1024<sup>3</sup></td><td>8.1 × 10<sup>10</sup></td><td>99</td><td>1LPT</td><td>DS</td></tr><tr><td rowspan="2">CLUSTERING DE 9, 10</td><td rowspan="2">k-evolution</td><td>1</td><td>2.0</td><td>2400<sup>3</sup></td><td>4.2 × 10<sup>10</sup></td><td rowspan="2">100</td><td rowspan="2">1LPT</td><td>w<sub>0</sub>c<sub>s</sub><sup>2</sup>CDM</td></tr><tr><td>1</td><td>2.0</td><td>1200<sup>3</sup></td><td>3.3 × 10<sup>11</sup></td><td>w<sub>0</sub>c<sub>s</sub><sup>2</sup>CDM</td></tr><tr><td rowspan="2">FORGE 11</td><td rowspan="2">MG-Arepo</td><td>100</td><td>0.5</td><td>1024<sup>3</sup></td><td>~10<sup>10</sup></td><td rowspan="2">127</td><td rowspan="2">2LPT</td><td>ΛCDM</td></tr><tr><td>98</td><td>0.5</td><td>1024<sup>3</sup></td><td>~10<sup>10</sup></td><td>f(R)</td></tr><tr><td>BRIDGE 12</td><td>MG-Arepo</td><td>98</td><td>0.5</td><td>1024<sup>3</sup></td><td>~10<sup>10</sup></td><td>127</td><td>2LPT</td><td>nDGP</td></tr><tr><td>PNG-UNIT 13</td><td>GADGET-2</td><td>1</td><td>1.0</td><td>4096<sup>3</sup></td><td>1.2 × 10<sup>9</sup></td><td>99</td><td>2LPT</td><td>PNG</td></tr></tbody></table>

**Notes.** All simulations were dark matter only, except the CIDER and DAKAR2 suites, which used a two-component collisionless approximation to follow the baryonic and dark matter components separately. Hydrodynamic simulations were not analysed in this project. (1) Rácz et al. (2023); (2) Carbone et al. (2016); (3) Rasera et al. (2022); (4) Cautun et al. (2018); (5) Fiorini et al. (2023); (6) Giocoli et al. (2018); (7) Baldi (2023); (8) Baldi & Simpson (2017); (9) Hassani et al. (2019); (10) Hassani et al. (2020); (11) Arnold et al. (2021); (12) Harnois-Déraps et al. (2023); (13) Gutiérrez Adame et al. (2024).

p-GADGET3 follows the evolution of CDM and neutrino particles, treating them as two distinct collisionless components. The reference cosmological parameters were chosen to be close to the baseline Planck 2013 cosmology (Planck Collaboration XVI 2014): $\Omega_b = 0.05$, $\Omega_m = 0.32$, $H_0 = 67.0 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $n_s = 0.96$, and $A_s = 2.127 \times 10^{-9}$. Given these values, the reference (i.e. the massless neutrino case) CDM-particle mass resolution is $m_{\text{CDM}}^p = 8.27 \times 10^{10} h^{-1} M_{\odot}$, which is decreased according to the mass of neutrino particles, in order to keep the same $\Omega_m$ among all the DEMNUNI simulations. In fact, massive neutrinos are assumed to come as a particle component in a three-mass-degenerate scenario, therefore, to keep $\Omega_m$ fixed, an increase in the massive neutrino density fraction yields a decrease in the CDM density fraction. The DEMNUNI simulations balance mass resolution and volume to include perturbations at both large and small scales. The simulations are characterised by a softening length of $\epsilon = 20 h^{-1} \text{ kpc}$, a comoving volume of $8 h^{-3} \text{ Gpc}^3$ filled with 2048$^3$ dark matter particles and, when present, 2048$^3$ neutrino particles. The simulations are initialised at $z_{\text{init}} = 99$ with Zeldovich initial conditions. The initial power spectrum is rescaled to the initial redshift via the rescaling method developed in Zennaro et al. (2017). Initial conditions are then generated with a modified version of the N-GenIC software, assuming Rayleigh random amplitudes and uniform random phases.

### 3.3. The RAYGAL simulations

The RAYGAL simulations (Breton et al. 2019; Rasera et al. 2022) are a set of two dark-matter only simulations in wCDM and ΛCDM cosmologies. The simulations were performed with the adaptive-mesh refinement (AMR) N-body code RAMSES (Teyssier 2002; Guillet & Teyssier 2011). These simulations have a box size of 2625 h$^{-1}$ Mpc for 4096$^3$ particles, which results in a smoothing scale of 5 h$^{-1}$ kpc at the maximum refinement level. Both simulations share the parameters $H_0 = 72.0 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $n_s = 0.963$, $\Omega_b = 0.04356$ and $\Omega_r = 8.076 \times 10^{-5}$. The flat ΛCDM simulation has a WMAP7 cosmology (Komatsu et al. 2011): $\Omega_m = 0.25733$, and $\sigma_8 = 0.80101$, while the flat wCDM simulation is consistent at the 1σ-level with a WMAP7 cosmology with $\Omega_m = 0.27508$, $\sigma_8 = 0.85205$, and $w = -1.2$. In both cases, Gaussian initial conditions are generated using a modified version of the code MPGRAFIC (Prunet et al. 2008) with the displacement field computed using second-order Lagrangian perturbation theory (2LPT) to minimise the effect of transients (Crocce et al. 2006). The initial redshift has been set to $z_{\text{init}} \sim 46$ such as to ensure that the maximum displacement is of the order of one coarse cell. Such a late start guarantees smaller discreteness errors (see Michaux et al. 2021, for more details). For the present work, we focus on the snapshots at $z = 0, 1$, and 2.
### 3.4. The ELEPHANT simulation suite

The Extended LEnsing Physics using ANalYTic ray Tracing (ELEPHANT) cosmological simulation suite was run using the ECOSMOG simulation code (Li et al. 2012, 2013b; Barreira et al. 2015; Bose et al. 2017), which is based on the dark matter and hydrodynamic AMR simulation code RAMSES and includes various types of modified gravity models (e.g. Li et al. 2012; Brax et al. 2012, 2013; Li et al. 2013b,a; Becker et al. 2020). It is particularly designed to solve for a non-linear scalar field using AMR. New simulations were run for the purpose of testing the effective field theory of large-scale structure (EFTofLSS) pipeline for spectroscopic galaxy clustering (Cautun et al. 2018; Fiorini et al. 2021; Casas et al. 2023; Euclid Collaboration: Koyama et al. 2024). For this purpose, 11 simulations were carried out using the Euclid reference cosmology without massive neutrinos for ΛCDM and the nDGP model (Table 2 of Euclid Collaboration: Knabenhans et al. 2021). The cosmological parameters of the ΛCDM simulations are: $\Omega_m = 0.319$, $\Omega_b = 0.049$, $\Omega_\Lambda = 0.681$, $H_0 = 67.0 \text{ km s}^{-1}\text{Mpc}^{-1}$, $A_s = 2.1 \times 10^{-9}$, and $n_s = 0.96$. The nDGP simulations used the same parameters as the ΛCDM simulations with the cross-over scale $r_c = 1.2 c/H_0$. All of the simulations in this simulation suite had a box size of $1024 h^{-1} \text{ Mpc}$ and $1024^3$ particles. The initial conditions were generated at $z_{\text{init}} = 49$ with 2LPT using the FML code<fnref n="2" /> with fixed initial amplitudes. The phases of 10 realisations were extracted with different random seeds, while one realisation shares the same random seed as one of the other simulations, but with opposite phases to have a single paired-and-fixed simulation pair with suppressed cosmic variance (Angulo & Pontzen 2016). Output redshifts were selected from the Euclid Collaboration forecast paper for galaxy clustering (Euclid Collaboration: Blanchard et al. 2020, $z = 1.0, 1.2, 1.4$ and 1.65).

### 3.5. The COLA HIRES simulations

This simulation series contains overall seven simulations in ΛCDM and nDGP cosmologies that were run with MG-COLA, a modified gravity extension of the COmoving Lagrangian Acceleration (COLA) algorithm as implemented in the FML code. The COLA method uses a combination of analytic 2LPT displacement and particle mesh (PM) simulations to perform fast approximate simulations (Tassev et al. 2013). These techniques are extended to modified-gravity models using approximate screening methods to preserve the speed advantage of COLA simulations (Winther et al. 2017). The downside of PM simulations is that the internal structure of dark matter haloes is not well resolved due to limited resolution. This has an important implication for dark matter halo statistics. To mitigate this problem, the COLA simulations were run with an increased mass resolution (Fiorini et al. 2023). All simulations in this suite have a box size of $1024 h^{-1} \text{ Mpc}$, with $1024^3$ particles. The base cosmological parameters of the simulations are the Planck 2015 parameters (Planck Collaboration XIII 2016): $\Omega_m = 1 - \Omega_\Lambda = 0.3089$, $\Omega_b = 0.0486$, $H_0 = 67.74 \text{ km s}^{-1}\text{Mpc}^{-1}$, $n_s = 0.9667$, and $\sigma_8 = 0.8159$. This simulation series focusses on nDGP gravity and tested four cases: $r_c = \{0.5, 1, 2, 5\} c/H_0$. The series contains paired-and-fixed simulations (Angulo & Pontzen 2016) to suppress cosmic variance in ΛCDM and in the nDGP model for $r_c = 1 c/H_0$, while for the others they were only run for a single fixed amplitude realisation. The initial conditions were generated at $z_{\text{init}} = 127$ using 2LPT. Full particle snapshots were stored at four redshift values, $z = 1.0, 1.2, 1.4$, and 1.65, motivated by the expected Hα-emitters redshifts in the Euclid spectroscopic survey (Euclid Collaboration: Blanchard et al. 2020).

### 3.6. The DUSTGRAIN and DUSTGRAIN-PF simulations

The DUSTGRAIN (Dark Universe Simulations to Test GRAvity In the presence of Neutrinos) project is an initiative aimed at investigating the degeneracy between $f(R)$ gravity and massive neutrinos at the level of non-linear cosmological observables, which was first pointed out in Baldi et al. (2014). More specifically, the project includes two suites of cosmological dark-matter-only simulations named the DUSTGRAIN-pathfinder (DUSTGRAIN-PF, Giocoli et al. 2018) and the DUSTGRAIN-fullscale simulations that have been run by joining the MG-GADGET (Puchwein et al. 2013) solver for $f(R)$ gravity and the massive neutrinos implementation (Viel et al. 2010) available within the p-GADGET3 code. The former has been described and validated in Winther et al. (2015) and Euclid Collaboration: Adamek et al. (2025), while the latter has been compared with other methods in Adamek et al. (2023).

The DUSTGRAIN-PF simulations have been developed to sample the joint ($f_{R0}, m_\gamma$) parameter space to identify the most degenerate combinations of parameters with respect to some basic LSS statistics. These include the non-linear matter power spectrum, the halo mass function, weak-lensing-convergence power spectrum, various higher-order statistics, cosmic voids, velocity fields (see Peel et al. 2018, 2019; Merten et al. 2019; Contarini et al. 2021; García-Farieta et al. 2019; Hagstotz et al. 2019a,b; Boyle et al. 2021). This series includes in total 13 simulations in $f(R) + m_\gamma$ cosmology, plus an additional suite of 12 standard ΛCDM simulations for varying one single standard cosmological parameter at a time that have been specifically run for the Higher-Order Weak Lensing Statistics (HOWLS) project (Euclid Collaboration: Ajani et al. 2023). These simulations have a box size of $750 h^{-1} \text{ Mpc}$ per side, used a softening length of $\epsilon = 20 h^{-1} \text{ kpc}$, and include $(2\times)768^3$ particles (for the CDM and neutrinos components). The cosmological parameters (for the reference ΛCDM cosmology with massless neutrinos) have been set to $\Omega_m = 1 - \Omega_\Lambda = 0.31345$, $\sigma_8 = 0.842$, $H_0 = 67.31 \text{ km s}^{-1}\text{Mpc}^{-1}$, $n_s = 0.9658$, and the total matter density has been kept constant when varying the neutrino mass. Full snapshots have been stored at 34 output times between $z = 99$ (corresponding to the starting redshift of the simulation) and $z = 0$.

The DUSTGRAIN-fullscale simulations include only three runs (a reference ΛCDM cosmology and two $f(R)$ gravity models with $f_{R0} = -10^{-5}$ and different values of the total neutrino mass, namely $m_\gamma = \{0.1, 0.16\} \text{ eV}$) simulated in a $8 h^{-3} \text{ Gpc}^3$ volume containing $(2\times)2048^3$ particles. In order to allow for a direct comparison with the DEMNUNI simulations described above, and to produce an extension to the latter for $f(R)$ gravity with massive neutrino cosmologies, the DUSTGRAIN-fullscale simulations share the same initial conditions with DEMNUNI for each of the values of the neutrino mass. Therefore, the two sets of simulations have the same statistical realisations of the Universe and identical cosmological parameters. Full snapshots have been stored for 73 output times between $z = 99$ (i.e. the initial conditions) and $z = 0$.

### 3.7. The CIDER simulations

The Constrained Interacting Dark Energy scenario (or CIDER, Barros 2019) is a particular type of coupled Quintessence models

<footnote n="2">² https://github.com/HAWinther/FML</footnote>
characterised by a background cosmic expansion that is fixed by construction to be identical to a standard ΛCDM cosmology. As is discussed in Sect. 2.1.4, this implies refraining from choosing a priori any specific functional form for the scalar self-interaction potential and letting the dynamic evolution of the field sample the potential shape required to match the imposed expansion history. The main feature of the CIDER models is that they show a suppressed growth of structures compared to a standard ΛCDM model with the same expansion history, thereby possibly easing the $\sigma_8$ tension without further exacerbating the tension on $H_0$. For these reasons, the model has received some attention even though – at least in its original form – it may already be quite tightly constrained by CMB observations (Barros et al. 2023). The CIDER simulations have been run with the c-GADGET code (Baldi et al. 2010, see also Euclid Collaboration: Adamek et al. 2025) that implements all the relevant features of interacting dark energy models, and includes three values of the coupling $\beta = 0.03, 0.05, 0.08$ besides a reference ΛCDM cosmology corresponding to the case $\beta = 0$. All simulations clearly share the same expansion history, consistent with the following cosmological parameters: $\Omega_m = 1 - \Omega_\Lambda = 0.311$, $\Omega_b = 0.049$, $H_0 = 67.7 \text{ km s}^{-1}\text{Mpc}^{-1}$, $n_s = 0.9665$, $A_s = 1.992 \times 10^{-9}$, corresponding to a value of $\sigma_8 = 0.788$ at $z = 0$ in the reference ΛCDM model. The simulations follow the evolution of $2 \times 10^{24}^3$ particles for the (coupled) dark matter and (uncoupled) baryon components in a cosmological volume of $1 h^{-3} \text{Gpc}^3$ with a softening length of $\epsilon = 25 h^{-1} \text{ kpc}$. The baryonic species are treated as a separate family of collisionless particles, in other words, neither hydrodynamic forces nor radiative processes are considered in the simulations, and its inclusion is required in order to consistently represent the effects of the non-universal coupling characterising these models. Therefore, baryonic particles will interact with other massive particles according to standard Newtonian forces, while the interaction between pairs of CDM particles will be governed by an effective gravitational constant $G_{\text{eff}} = G_N [1 + (4/3)\beta^2]$ (see e.g. Amendola 2004; Baldi et al. 2010). Full snapshots have been stored for 25 output times between $z = 99$ and $z = 0$.

### 3.8. The DAKAR and DAKAR2 simulations

The dark scattering (DS) scenario (Simpson 2010) is another particular class of coupled Quintessence models where a non-universal interaction between dark matter particles and a classical scalar field playing the role of dark energy is characterised by a pure momentum exchange between the two species, with no transfer of rest-frame energy (see e.g. Pourtsidou et al. 2013; Skordis et al. 2015). In this respect, this interaction resembles a process of elastic scattering of massive particles (i.e. the dark matter) moving in a homogeneous fluid with an equation-of-state parameter, $w$ (i.e. the dark energy field), which can be simulated by introducing a velocity-dependent force acting on dark matter particles that will depend on the evolution of the dark energy equation-of-state parameter, $w$, and on the cross-section, $\sigma$, characterising the interaction strength (Baldi & Simpson 2015).

The DAKAR (Baldi & Simpson 2017) and DAKAR2 simulations have been run with the c-GADGET code and cover various combinations of the shape of $w(z)$, including the CPL parametrisation given by Eq. (1) and hyperbolic tangent shapes, and of the cross-section, $\sigma$, giving rise to a diverse phenomenology at both linear and non-linear scales. In particular, DS models have been shown to suppress the linear growth of perturbations for equation-of-state parameters $w > -1$ (Pourtsidou & Tram 2016;

Bose et al. 2018; Carrilho et al. 2022), thereby possibly addressing the $\sigma_8$ tension, but such suppression is typically paired with a substantial enhancement of structure growth at deeply non-linear scales.

The DAKAR simulations are subject to the approximation of considering the entirety of matter in the Universe is in the form of dark matter, thereby slightly overestimating the effect of the interaction as well as not capturing the segregation effects between dark matter and baryons due to the non-universality of the coupling. These have been run for a cosmology with $\Omega_m = 1 - \Omega_\Lambda = 0.308$, $H_0 = 67.8 \text{ km s}^{-1}\text{Mpc}^{-1}$, $n_s = 0.966$, $A_s = 2.215 \times 10^{-9}$, in a simulation box with a volume of $1 h^{-3} \text{Gpc}^3$ filled with $10^{24}^3$ dark matter particles and using a softening length of $\epsilon = 12 h^{-1} \text{ kpc}$.

The DAKAR2 simulations, instead, share the same cosmology and the same statistical realisation as the CIDER simulations described above (i.e. the two sets of simulations share exactly the same reference ΛCDM run) and include collisionless baryons as a separate family of uncoupled particles, thereby consistently capturing the non-universality of the DS interaction. As for the CIDER simulations, a collection of 25 full snapshots for redshifts between $z = 99$ and $z = 0$ has been stored.

### 3.9. The CLUSTERING DE simulations

The CLUSTERING DARK ENERGY simulations are run using the k-evolution code, a relativistic N-body code (Hassani et al. 2019, 2020) based on gevolution-1.2 (Adamek et al. 2016). In k-evolution, the field equations for k-essence type theories (Eq. 2) are solved using the effective field theory (EFT) framework. We have two free parameters in the EFT framework of these theories: the equation-of-state parameter, $w(\tau)$, appearing at the background level and kineticity, $\alpha_K(\tau)$, at the perturbation level. In the fluid picture of these theories, the relevant parameters are the speed of sound, $c_s(\tau)$, and the equation-of-state parameter, $w(\tau)$, which in general are time-dependent. The term ‘clustering dark energy’ refers to the fact that these theories include a sound-horizon scale, beyond which scalar-field perturbations can grow. In the analysed simulations, constant $w_0$ and $c_s^2$ are used, with cosmological parameters based on the Euclid reference cosmology (Euclid Collaboration: Knabenhans et al. 2021). The suite contains one ΛCDM simulation and four clustering dark energy simulations: $(w_0, c_s^2) = (-0.9, 1 c^2)$, $(-0.9, 10^{-4} c^2)$, $(-0.9, 10^{-7} c^2)$, and $(-0.8, 10^{-7} c^2)$. In these simulations, the box size was set to $2 h^{-1} \text{Gpc}$ with $N = 1200^3$ particles. Moreover, two sets of simulations with different resolutions were considered to study the convergence of the results. In this high-resolution simulation set, the box size was set to $2 h^{-1} \text{Gpc}$ with $N = N_{\text{grid}} = 2400^3$. In this series, the particle snapshots were saved in GADGET-2 format at five different redshifts $z \in \{2, 1.5, 1, 0.5, 0\}$.

### 3.10. The FORGE and BRIDGE simulation suites

The FORGE simulation suite (Arnold et al. 2021) is a set of 198 dark matter only simulations for $f(R)$ gravity and ΛCDM run with the Arepo cosmological simulation code (Springel 2010; Weinberger et al. 2020) using its MG module (Arnold et al. 2019). The simulations explore the cosmological and $f(R)$ parameter space spanned by $\Omega_m (\Omega_\Lambda = 1 - \Omega_m)$, $h$, $\sigma_8$, and $f_{R0}$ through 50 combinations (nodes) of these parameters sampled in a Latin-hypercube. All other cosmological parameters are fixed to a Planck cosmology ($n_s = 0.9652$, $\Omega_b = 0.049199$, Planck Collaboration VI 2020). For each node, FORGE consists
of a pair of large box simulations with 512³ particles in a 1.5 h⁻¹ Gpc side-length box and a pair of high-resolution runs with 1024³ particles in a 500 h⁻¹ Mpc box. For each pair, the initial conditions were chosen such that the large-scale variance in the 3D matter power spectrum approximately cancels when averaged over the two simulations (see Arnold et al. 2021; Ruan et al. 2024; Harnois-Déraps et al. 2023, for further details and some applications of these simulations). All simulations in this suite started from zₐₗₐₗ = 127, with initial conditions generated using the 2LPTic (Crocce et al. 2006) code.

The BRIDGE simulation suite (Harnois-Déraps et al. 2023) uses the same base setup and ICs as FORGE, but for the nDGP gravity model implemented into the Arepo code (Hernández-Aguayo et al. 2021). Accordingly, instead of varying the $f_{R0}$ parameter, the nDGP parameter, $H_0 r_c$ is varied to explore the cosmological parameter space, with all the other parameters identical to those in FORGE for corresponding nodes.

The fact that the FORGE and BRIDGE simulations have the same cosmological parameters, node by node, allows a third suite of simulations to be done as a control set to quantify the effect of modified gravity and how it correlates to the effect of varying cosmological parameters. This additional suite of simulations, FORGE-ΛCDM, uses the same setup but runs for the ΛCDM counterparts of the corresponding $f(R)$ and nDGP models.

As the FORGE, FORGE-ΛCDM, and BRIDGE simulations have different cosmologies, the mass resolution differs amongst them. In Table 1, we thus quoted an order of magnitude for the mass resolution.

### 3.11. The PNG-UNIT simulation

The PNG-UNIT (Gutiérrez Adame et al. 2024) is a twin of one of the existing UNITsims (Universe N-body simulations for the Investigation of Theoretical models, Chuang et al. 2019), but with local primordial non-Gaussianities given by $f_{NL} = 100$. The simulation assumes the following ΛCDM parameters: $\Omega_m = 1 - \Omega_\Lambda = 0.3089$, $H_0 = 67.74 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $n_s = 0.9667$, $\sigma_8 = 0.8147$. It consists of $N = 4096^3$ particles in $L = 1 h^{-1} \text{ Gpc}$ evolved with L-GADGET, which is a version of GADGET-2 optimised for massive parallelisation, using a tree-PM algorithm with a softening length of $\epsilon = 6 h^{-1} \text{ kpc}$. The initial conditions are run with the 2LPT implementation in the FastPM (Feng et al. 2016) code at $z = 99$. Both the UNIT and PNG-UNIT are run with fixed initial conditions (Angulo & Pontzen 2016), which set the amplitude of the ICs to their expected value. Whereas there are four UNIT simulations in two sets of pairs (within each pair, each simulation has the inverted phases one with respect to another, following Angulo & Pontzen 2016), we only have one simulation for the PNG-UNIT. The PNG-UNITsim is run with the phases of the ICs matched to one of the UNITsims, which is labelled in the databases as ‘Ampl1’. The usage of fixed ICs with local PNG was validated in Avila & Adame (2023), where it was also shown how to increment the precision of the statistics measured from matched simulations. Overall, 129 snapshots were stored during the simulation, and 32 in the $0.5 < z < 2.0$ range.

## 4. Analysis

We have developed a cosmological analysis pipeline to generate mock observables from non-standard cosmological simulations

in a consistent and rapid way. The pipeline is a SLURM<fnref n="3" /> script that runs in parallel on multiple nodes on the machines where the simulations are stored. The pipeline consists of several modules that can be activated or deactivated independently. The modules are controlled by a configuration file that specifies the input and output parameters, as well as the options for each module. The input of the pipeline is the particle snapshots of the non-standard cosmological simulations. The supported input formats are GADGET binary and Hierarchical Data Format version 5 (HDF5, Springel 2005), Arepo HDF5 (Weinberger et al. 2020), RAMSES HDF5 format from RAYGAL (Roy et al. 2014; Rasera et al. 2022), and GIZMO HDF5 (Hopkins 2015). The main steps of the analysis are summarised in Fig. 1. In this section, we describe the quantities generated by this pipeline.

### 4.1. Dark matter density field

The pipeline uses nbodykit (Hand et al. 2018) to read and analyse the dark matter particle data of the input-simulation snapshots. This Python package is an open-source, massively parallel toolkit that provides a set of LSS algorithms useful in the analysis of cosmological data sets from N-body simulations and observational surveys. During the dark matter density-field analysis, nbodykit generates a reconstructed density field from the input-particle distribution with the triangular-shaped-cloud (TSC) density-assignment function. We chose to use a

$$N_{\text{grid}} = 2^{\lfloor \log_2(\sqrt[3]{N_{\text{part}}}) \rfloor - 1} \quad (35)$$

linear grid size for every analysed snapshot, where $N_{\text{part}}$ is the number of stored particles in the snapshot. With this choice, there will always be at least eight particles on average in each cubic density cell. The reconstructed density fields were saved in bigfile format (Feng et al. 2017) for future analysis.

#### 4.1.1. Real-space power spectrum

The real-space matter power spectrum is defined via

$$\langle \tilde{\delta}(\mathbf{k}) \tilde{\delta}^*(\mathbf{k}') \rangle = (2\pi)^3 P(\mathbf{k}) \delta_D^3(\mathbf{k} - \mathbf{k}') \quad (36)$$

where $\tilde{\delta}(\mathbf{k})$ is the Fourier-transform of the matter overdensity field

$$\delta(\mathbf{r}) = \frac{\rho(\mathbf{r})}{\bar{\rho}} - 1 \quad (37)$$

and $\mathbf{k}$ is the wavevector. We estimated the power spectrum using nbodykit. The density field was created by binning the particles into a grid using a TSC-density-assignment function, with the linear-grid size defined in Eq. (35). The density field was Fourier-transformed and the power spectrum was computed by binning $|\tilde{\delta}(\mathbf{k})|^2$, deconvolving the window function and subtracting shot noise. We also used the interlacing technique to reduce aliasing (Sefusatti et al. 2016). The bin size of the power spectrum was set to

$$\Delta k = k_{\text{f}} = \frac{2\pi}{L_{\text{box}}} \quad (38)$$

where $k_{\text{f}}$ is the fundamental wavenumber, and $L_{\text{box}}$ is the linear size of the simulation. The pipeline saves the power spectrum of every calculated bin below the

$$k_{\text{Ny}} = \frac{\pi N_{\text{grid}}}{L_{\text{box}}} \quad (39)$$

<footnote n="3">³ https://slurm.schedmd.com/</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Flowchart summarising the main steps of the analysis pipeline. In the first step, the pipeline executes the Rockstar halo finder to generate the halo catalogues and additional BGC2 particle data containing all relevant information for the halo profile calculation. Then, the 2D and 3D halo profiles, the halo mass function, and void catalogues are calculated by the corresponding modules. The real- and redshift-space halo power spectra and the triangular shaped cloud (TSC) reconstructed halo density fields on a regular cubic grid are computed for a predefined mass bin. After this, the dark matter power spectrum calculator module computes the real- and redshift-space dark matter power spectra with the reconstructed TSC dark matter density field. By using the real-space dark matter and halo power spectra, the linear-halo-bias estimator module calculates the halo-bias table using only the linear scales. Finally, by using this halo-bias table, the linear matter power spectrum, and the cosmological parameters, the halo redshift-space power spectrum Gaussian covariances are computed. This process is repeated for all selected particle snapshots.</figcaption></figure>

Nyquist wavenumber with the number of modes into a simple ASCII format file.

### 4.1.2. Redshift-space power spectrum

The real-space matter power spectrum is not directly measurable in galaxy surveys because we cannot probe the real-space positions of galaxies. What we can directly measure is the redshift-space power spectrum $P^s(k, \mu)$ where $\mu = \hat{n}_{\text{LOS}} \cdot \hat{k}$ and $\hat{n}_{\text{LOS}}$ is a unit vector in the line-of-sight (LOS) direction. This can be expanded in multipoles, $P^s(k, \mu) = \sum_{l=0}^{\infty} P_l^s(k) \mathcal{L}_l(\mu)$, where $\mathcal{L}_l(\mu)$ are the Legendre polynomials. The multipoles are then computed from the redshift-space power spectrum as

$$P_l^s(k) = \frac{2l + 1}{2} \int_{-1}^{1} P^s(k, \mu) \mathcal{L}_l(\mu) d\mu. \quad (40)$$

We computed the redshift-space power spectrum in 25 $\mu$ bins and the redshift-space multipoles (the monopole $P_0$, the quadrupole $P_2$, and the hexadecapole $P_4$) using nbodykit from the input dark matter density field. For this, we used the distant-observer approximation

$$s_i = r_i + \left( \hat{n}_{\text{LOS}} \cdot \frac{v_i}{aH} \right) \cdot \hat{n}_{\text{LOS}}, \quad (41)$$

to add the redshift-space distortions using the three coordinate axes as the LOS directions (observables were computed as the mean over these three individual axes). Here, $r_i$ and $v_i$ are the real-space particle coordinates and peculiar velocities inside the periodic simulation box, $s_i$ is the corresponding redshift-space position we computed, $a = 1/(z + 1)$ is the scale factor, and $H$ is the Hubble parameter at the redshift of the snapshot. We deconvolved the window function for the density assignment, applying interlacing, and the resulting density power spectrum
was finally shot-noise subtracted. The saved wavenumber bins are the same as in Sect. 4.1.1.

### 4.1.3. Linear dark matter power spectrum

During the analysis of scale-independent models, multiple modules use the linear dark matter power spectrum as an input. To make the analysis more transparent, we also generated and saved the linear power spectrum for all analysed redshifts. For this, the pipeline needs the $P_{\text{lin}}(k, z_{\text{start}})$ input linear power spectrum of the simulation that was used during the initial condition generation. This can be defined at any $z_{\text{start}}$ redshift. Then, this linear power spectrum is renormalised with the cosmological parameter $\sigma_8$ at $z=0$. The normalised $P_{\text{lin}}(k, z=0)$ power spectrum is rescaled to $z_{\text{snap}}$ redshifts of all analysed snapshots as

$$P_{\text{lin}}(k, z_{\text{snap}}) = \left[ \frac{D(z_{\text{snap}})}{D(z=0)} \right]^2 P_{\text{lin}}(k, z=0), \quad (42)$$

where $D(z)$ is the linear growth function. The pipeline uses this back-scaling since the linear growth in these Newtonian simulations follows this scale-independent evolution. For ΛCDM reference simulations, we used the

$$D(a) = \frac{5\Omega_m H_0^2}{2} H(a) \int_0^a \frac{da'}{\dot{a}'^3} \quad (43)$$

linear growth to scale the linear spectrum (Peebles 1993). This growth function only describes the linear growth in the ΛCDM framework. In the case of wCDM and CPL models, we solve the

$$G'' + \left[ \frac{7}{2} + \frac{3}{2} \frac{w(a)}{1+X(a)} \right] \frac{G'}{a} + \frac{3}{2} \frac{1-w(a)}{1+X(a)} \frac{G}{a^2} = 0 \quad (44)$$

ordinary differential equation (Linder & Jenkins 2003) with the COLOSSUS python package (Diemer 2018), where $G(a) = D(a)/a$ and

$$X(a) = \frac{\Omega_m}{1 - \Omega_m} e^{-\frac{1}{2} \int_a^1 da' \frac{w(a')}{a'}}, \quad (45)$$

For every other model, we use tabulated linear growth functions. The linear power spectra are calculated in the same wavenumber bins as the non-linear real-space matter power spectra.

### 4.2. Halo catalogues

Rockstar (Behroozi et al. 2013) is a friends-of-friends (FoF) halo-finder algorithm that uses information from the full 6D phase space (positions and the velocities) of the particles. The code initially creates FoF groups in real space, with a large linking length ($b \approx 0.28$). It then does a new FoF search using the phase-space metric

$$d = \sqrt{\frac{|\mathbf{x}_1 - \mathbf{x}_2|^2}{\sigma_x^2} + \frac{|\mathbf{v}_1 - \mathbf{v}_2|^2}{\sigma_v^2}}, \quad (46)$$

where $\sigma_x$ and $\sigma_v$ are the particle-position and velocity dispersions for the given FoF group. Finally, it links particles into subgroups and this is done iteratively on each subgroup creating a hierarchical set of structures. By default, the algorithm calculates halo and subhalo masses using dark matter particles from the spherical regions around the friends-of-friends group with

<figure><img src="image_3.png" /><figcaption>Fig. 2. Calculated halo mass function of the ELEPHANT simulation suite. The vertical and horizontal dash-dotted lines indicate the mass relative to haloes with 50 particles and the number density relative to a 1% shot-noise error, respectively. The shaded region highlights the mass bin used to calculate the halo power spectra shown in Fig. 3.</figcaption></figure>

gravitationally unbound particles removed. The halo masses calculated this way are called bound-only (BO) masses. If the unbound particles are not removed during the mass calculation, the calculated masses are strict spherical-overdensity (SO) masses.

We made a custom version of the publicly available Rockstar code (Behroozi et al. 2012) to analyse our non-standard simulations. We added new input formats for simulations, options to read tabulated expansion histories, and already internally computed quantities in the outputs such as halo minor axis vectors and radii at different mass definitions. None of these modifications impact the halo-finding algorithm. In our pipeline, we use the following mass-definitions: $M_{200c}$ (SO & BO), $M_{500c}$ (SO), $M_{1000c}$ (SO), $M_{2500c}$ (SO), $M_{200b}$ (SO). The $M_{\text{vir}}$ masses are not calculated by the pipeline, since this mass definition is dependent on the cosmological parameters and on the laws of gravity. Many non-standard cosmological models are changing the dynamics of the dark matter component, and this choice simplifies the future expansion of the database without the need of implementing new cosmologies in Rockstar. After the catalogue (in ASCII format) is produced, we run a post-processing script to find parent haloes for subhaloes and store the information as an additional index column. Extra information is saved in the header such as scale factor, box length, and particle mass. Additional particle data for each halo are also saved by Rockstar in a custom BGC2 binary data format. During the execution of our pipeline, these BGC2 files are temporarily stored to provide additional input for other analysis modules.

#### 4.2.1. Halo mass function and power spectra

By default, we compute the halo mass function (HMF) in the range $11 < \log_{10}[M/(h^{-1} M_{\odot})] < 14$ with 32 logarithmic bins using the main mass definition (200c) and excluding
<figure><img src="image_4.png" /><figcaption>Fig. 3. Calculated matter and halo power spectra of the ELEPHANT simulation suite in the mass bin $10^{12.7}h^{-1} M_{\odot} < M_{\text{halo}} < 10^{13.2}h^{-1} M_{\odot}$. The solid lines represent the reference ΛCDM simulations, while the dashed lines the results of the nDGP simulations. Top left: Real-space power spectra for dark matter. The dots above the solid lines highlight the locations where the power spectrum is estimated. Top right: Real-space power spectra for haloes. Bottom left: Monopole of the halo power spectrum in redshift space. The shaded region shows the calculated variance of the monopole of the ΛCDM halo power spectrum at redshift $z = 1$. Bottom right: Quadrupole of the halo power spectrum in redshift space. The shaded region represents the calculated variance of the quadrupole of the ΛCDM halo power spectrum at redshift $z = 1$.</figcaption></figure>

substructures. The pipeline allows the user to use different mass definitions (see Sect. 4.2), include substructures, and vary the HMF range and binning.

In practice, many simulations produce very large ASCII files that are not practical to read using standard libraries such as NumPy<fnref n="4" /> or Pandas<fnref n="5" />. To speed up the analysis pipeline, we therefore use Polars<fnref n="6" />, a fast multi-threaded dataframe library.

The halo real-space and redshift-space power spectra are computed using the same tools as in Sects. 4.1.1 and 4.1.2. The user can specify the halo mass range, SO or BO for the main mass definition, and whether or not to include substructures.

### 4.2.2. Halo bias

For each catalogue, we infer the linear halo bias with the estimator

$$b = \left\langle \sqrt{\frac{P_h(k)}{P_m(k)}} \right\rangle_{k<k_{\text{max}}}, \qquad (47)$$

<footnote n="4">⁴ https://numpy.org/</footnote>

<footnote n="5">⁵ https://pandas.pydata.org/</footnote>

<footnote n="6">⁶ https://pola.rs/</footnote>
with $P_m(k)$ and $P_h(k)$ the matter and halo real-space power spectra, respectively estimated in Sects. 4.1.1 and 4.2.1. This estimator calculates the bias by taking the square root of the ratio of the halo power spectrum $P_h(k)$ to the matter power spectrum $P_m(k)$, and then averaging this ratio over all $k$ bins where $k < k_{max}$ with uniform weighting. We only use this computed quantity to be as model-independent as possible and to remove cosmic variance (since matter and halo both share the same sample and cosmic variance). We compute the mean power spectra ratio up to a conservative value of $k_{max} = 0.1 h \text{ Mpc}^{-1}$ to mitigate the effects of non-linear clustering. This method works reliably for scale-independent bias with sub-percent accuracy, but cannot be used for models that have scale-dependent bias at $k < k_{max}$ wavenumber.

### 4.2.3. Redshift-space Gaussian covariance

To produce Gaussian covariances of the power-spectrum multipoles in redshift space, we need the linear bias and power-spectrum multipoles including the shot-noise contributions as inputs. The former is estimated numerically in Sect. 4.2.2, while the latter can be internally computed from an input linear power spectrum; in this case, the covariance is estimated as in Taruya et al. (2010), or with the EFT model using the COMET emulator (Eggemeier et al. 2023) and the covariance formulae from Grieb et al. (2016).

When analysing snapshots, it may be interesting to compute the power-spectrum multipoles averaged over the three box directions to significantly suppress variance. However, this procedure also has to be carefully accounted for in the covariance (Smith et al. 2021) since the LOS-averaged covariance is not equal to the single-LOS covariance divided by three (as one might naively expect). Thanks to the LOS-averaged covariance implemented in COMET, it can also be part of the outputs of our pipeline.

### 4.2.4. 2D and 3D halo profiles

The generation of binned three-dimensional and two-dimensional projected profiles was performed using a custom analysis module that reads the halo catalogues as well as the BGC2 particle data, and stores the resulting profiles of each halo in a separated HDF5 file.

All the profiles are obtained considering 50 log-spaced bins in a fixed radial range [0.001, 5], in units of $r_{500c}$. This analysis module provides cumulative mass, density, number density, and cumulative number density profiles, as well as velocity and velocity dispersion profiles for the cartesian- and spherical-coordinates components of the velocities. The 2D profiles correspond to projecting the LOS along each of the cartesian coordinates in a cylinder of length $5r_{500c}$. In an upcoming version of the pipeline, the profiles will also be available for the projections along the axes of the inertia ellipsoid *a*, *b*, and *c*.

### 4.3. Cosmic voids

The Void Identification and Examination toolkit, VIDE (Sutter et al. 2015), is a parameter-free topological void finder, conceived for galaxy-redshift surveys and N-body simulations. The VIDE pipeline is an open-source Python/C++ code, based on the ZOBOV (Neyrinck 2008) software and can be launched on any tracer distribution. The algorithm follows the following main steps: i) estimation of the density field of a tracer distribution using the Voronoi tessellation (Schaap & van de Weygaert

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Calculated halo mass function from the FORGE simulation suite. The vertical and horizontal dash-dotted lines indicate the mass relative to halos with 50 particles and the number density relative to a 1% shot-noise error, respectively. The shaded region highlights the mass bin used to calculate the halo power spectra shown in Fig. 5.</figcaption></figure>

2000); ii) detection of all the relative minima; iii) merging of nearby Voronoi cells into zones via the watershed transform (Platen et al. 2007), cells correspond to local catchment ‘basins’, which are identified as voids. VIDE can also merge adjacent voids to construct a nested hierarchy of voids if a merging threshold is provided. In this case, when two adjacent voids have at least one Voronoi cell on the ridge separating them lower than the threshold, they are merged into a parent void. In this work, in order to leave the algorithm parameter-free, and for consistency with other Euclid void analyses (Hamaus et al. 2022; Contarini et al. 2022), we do not explore this possibility.

VIDE provides some fundamental properties of voids. The void size is measured by the effective radius, defined as the radius of a sphere with the same volume as the void, $R_{\text{eff}} = [(3/4\pi) \sum_i V_i]^{1/3}$, where $V_i$ is the volume of the $i^{\text{th}}$ Voronoi cell belonging to the void. The void centre is defined as the volume-weighted barycentre, $X_v = \sum_i x_i V_i / V_{\text{tot}}$, with $V_{\text{tot}} = \sum_i V_i$. We note that this corresponds to the geometric centre of the void. In addition, VIDE also provides the position of the tracer sitting in the lowest-density Voronoi cell, that is, the minimum. The void’s depth is estimated via the central density, defined as the mean density in a sphere centred in the barycenter, $X_v$, with radius $R_{\text{eff}}/4$. VIDE also computes void shapes via the inertia tensor as well as the corresponding eigenvalues and eigenvectors. The ellipticity is then computed as $\epsilon = 1 - (J_1/J_3)^{1/4}$, where $J_1$ and $J_3$ are the smallest and largest eigenvalues.

We detect voids in the distribution of Rockstar haloes. After the void catalogues are produced, we post-process them to measure the void-size function, which is the number density of voids as a function of their size, $R_{\text{eff}}$. The void-size function is a sensitive probe for cosmology, strongly complementary to the galaxy 2pt-statistics (Pisani et al. 2015; Massara et al. 2015; Kreisch et al. 2019; Verza et al. 2019; Contarini et al. 2021, 2022, 2023; Verza et al. 2022, 2023, 2024). Additionally, albeit not computed for this paper, the void catalogues allow one
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Calculated matter and halo power spectra from the FORGE simulation suite in the mass bin $10^{12.7}h^{-1}M_{\odot} < M_{\text{halo}} < 10^{13.2}h^{-1}M_{\odot}$. The solid lines represent the reference ΛCDM simulations, while the dashed lines are the results of the nDGP simulations. **Top left:** Real-space power spectra for dark matter. The dots above the solid lines highlight the locations where the power spectrum is estimated. **Top right:** Real-space power spectra for haloes. **Bottom left:** Monopole of the halo power spectrum in redshift space. The shaded region shows the calculated variance of the monopole of the ΛCDM halo power spectrum at redshift $z=0$. **Bottom right:** Quadrupole of the halo power spectrum in redshift space. The shaded region represents the calculated variance of the quadrupole of the ΛCDM halo power spectrum at redshift $z=0$.</figcaption></figure>

to compute the void-galaxy cross-correlation function, another powerful statistic to constrain cosmology (see e.g. Hamaus et al. 2022).

**5. Interpretation**

The mock observables that we compute contain several interesting signatures of non-standard models. Due to a large number of analysed models, in this paper, we focus on showing results from the nDGP, $f(R)$, and interacting-dark-energy models, which we obtained thanks to the analysis of the ELEPHANT, FORGE, and CIDER simulations suites, respectively.

To mitigate the noise due to sample variance in the figures below, we average the signals over the available realisations of the ELEPHANT simulations. In the case of the FORGE suite, we focus on a single cosmology<fnref n="7" /> corresponding to $|\tilde{f}_{R0}| = 10^{-5.34219}$. In the case of the CIDER simulations, we focus on the IDE model with $\beta = 0.03$ coupling. The nDGP, $f(R)$, and IDE simulations were run with the same initial conditions as their ΛCDM counterpart so the effects of non-standard cosmologies can be studied directly by comparing the generated observables.

<footnote n="7">⁷ In particular we use node 10 of the Latin-Hypercube sampling described in Table 1 of Arnold et al. (2021).</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Calculated halo mass function from the CIDER simulation suite. The vertical and horizontal dash-dotted lines indicate the mass relative to halos with 50 particles and the number density relative to a 1% shot-noise error, respectively. The shaded region highlights the mass bin used to calculate the halo power spectra shown in Fig. 7.</figcaption></figure>

The main difference between the modified-gravity simulation and ΛCDM is the inclusion of a fifth force. For nDGP, this fifth force acts on all scales and increases in strength toward redshift zero (and goes to zero as we go to higher and higher redshifts). This is also true for $f(R)$, with the exception that the fifth force has only a finite range, so it does not affect the clustering on the largest scales.

In the IDE simulations, the dark energy component interacts with the dark matter, resulting in a transfer of energy between the two. This interaction affects the growth of cosmic structures by modifying the gravitational potential. This cosmological scenario is expected to suppress structure formation at late time compared to the standard ΛCDM model.

These differences lead to a number of different observable signatures, a few of which we highlight below.

**Abundance of dark matter haloes** – In Fig. 2, we compare the cumulative halo mass function of the nDGP and ΛCDM models. The inclusion of the fifth force means structures will form more rapidly than in ΛCDM and this is indeed what we see. This is most pronounced at the high-mass end where the abundance is up to 50% larger. In Fig. 4, we compare the cumulative halo mass function of the $f(R)$ and ΛCDM model. We see roughly the same qualitative features as for nDGP in that the halo abundance generally increases with halo mass and with redshift compared to the reference ΛCDM results. However, as opposed to nDGP, we see an over-abundance of ‘small’ haloes at earlier times in the $f(R)$ simulations. This is a consequence of the fact that the fifth force only acts on ‘small’ scales and the fact that the screening mechanism is more effective at suppressing the fifth force in and around the most massive haloes. The comparison of the halo mass function between the standard ΛCDM and the $\beta = 0.03$ IDE model can be seen in Fig. 6. In the IDE model, the interaction between the dark energy and dark matter caused a substantial reduction

in the HMF. This is a straightforward consequence of the suppressed growth rate of the matter fluctuations.

**Clustering of dark matter** – In Figs. 3, 5, and 7, we show the calculated real- and redshift-space power-spectrum multipoles for the haloes and dark matter for nDGP, $f(R)$, and IDE, respectively. For the modified gravity models, the effect of the fifth force is again clearly in the dark matter power spectrum and shows two different effects: for nDGP, we have a scale-independent growth rate causing the power spectrum to be boosted on all scales displayed, while for $f(R)$ we have a scale-dependent growth-rate where $f(R)$ agrees with ΛCDM on the largest scales, but is boosted below a critical scale that is related to the range of the fifth force. For both models, the difference with respect to ΛCDM increases in strength as we get closer to the present time. In the case of the interacting-dark-energy model, the energy transfer between the dark energy and dark matter caused a scale-independent suppression in the dark matter power spectrum. At redshift 2, this is ≈5%, and this difference increases to ≈5% for $z = 0.55$ compared to the ΛCDM model.

**Halo bias** – When it comes to halo clustering in real space, we see the opposite effect as for the dark matter power spectrum in Figs. 3 and 5, with nDGP and $f(R)$ being less clustered than ΛCDM. This comes from a smaller halo bias in these modified-gravity models (see e.g. Barreira et al. 2014 for a theoretical explanation for nDGP). In the interacting-dark-energy scenario, the halo bias in real space is 6% higher than for ΛCDM. As a consequence, the real-space clustering of the dark matter haloes is more prominent in the IDE simulation.

**Redshift space distortions** – For the redshift-space halo power spectra, the boost in the ratio with respect to ΛCDM is seen to be larger than in real space, which comes from the larger velocities in the modified-gravity simulations, leading to enhanced redshift-space distortions. The monopole redshift-space power spectra of the haloes in the mass bin $10^{12.7}h^{-1}M_{\odot} < M_{\text{halo}} < 10^{13.2}h^{-1}M_{\odot}$ in the IDE simulations are showing a 5–10% excess power compared to the ΛCDM counterpart, similarly to the real-space clustering. On the other hand, the quadrupole only shows a power increase at smaller, non-linear scales.

## 6. Summary

In this paper, we describe a new pipeline based on the Rockstar halo finder and the nbodykit LSS toolkit to post-process cosmological simulations with modified gravity, non-standard expansion history, modified dark matter or dark energy components, or altered initial conditions. We used this pipeline to analyse 474 cosmological N-body simulations in various ΛCDM and non-standard cosmological scenarios in a consistent way. With this pipeline, we generated halo catalogues, halo mass functions, reconstructed density fields, real- and redshift-space power spectra, Gaussian covariances, halo biases, and void catalogues. This generated data will serve as a theoretical prediction and reference for Euclid as well as other Stage-IV cosmology projects. Using the calculated quantities, we identified distinctive signatures of non-standard behaviour in nDGP and $f(R)$ modified-gravity models, and in the CIDER interacting-dark-energy scenario.

The synthetic halo catalogues are crucial in the production of additional observables, which can be used for a direct comparison with cosmological observations of Euclid. In the near future, we shall extend the generated database with halo density profiles (Navarro et al. 1996; Le Brun et al. 2018), synthetic galaxy catalogues (Berlind et al. 2003), weak
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Calculated matter and halo power spectra from the CIDER simulation suite in the mass bin $10^{12.7}h^{-1}M_{\odot} < M_{\text{halo}} < 10^{13.2}h^{-1}M_{\odot}$. The solid lines represent the reference ΛCDM simulations, while the dashed lines are the results of the IDE simulations. **Top left:** Real-space power spectra for dark matter. The dots above the solid lines highlight the locations where the power spectrum is estimated. **Top right:** Real-space power spectra for haloes. **Bottom left:** Monopole of the halo power spectrum in redshift space. The shaded region shows the calculated variance of the monopole of the redshift space ΛCDM halo power spectrum at redshift $z = 0.55$. **Bottom right:** Quadrupole of the halo power spectrum in redshift space. The shaded region represents the calculated variance of the quadrupole of the redshift space ΛCDM halo power spectrum at redshift $z = 0.55$.</figcaption></figure>

lensing (Jaroszynski et al. 1990; Bartelmann & Schneider 2001) and ISW (Giannantonio et al. 2008) maps, and lightcones (Merson et al. 2013).

We have generated overall more than 100 TB of post-processed data from the available non-standard simulations. During the analysis, the pipeline used 66 CPU hours and 60GB of memory per billion particles per snapshot on average. The data are available on request on the CosmoHub<fnref n="8" /> platform (Tallada et al. 2020; Carretero et al. 2017) designed for interactive exploration and distribution of massive cosmological datasets.

**Acknowledgements.** The Euclid Consortium acknowledges the European Space Agency and a number of agencies and institutes that have supported the development of Euclid, in particular the Agenzia Spaziale Italiana, the Austrian Forschungsförderungsgesellschaft funded through BMK, the Belgian Science Policy, the Canadian Euclid Consortium, the Deutsches Zentrum für Luft- und Raumfahrt, the DTU Space and the Niels Bohr Institute in Denmark, the French Centre National d'Etudes Spatiales, the Fundação para a Ciência e a Tecnologia, the Hungarian Academy of Sciences, the Ministerio de Ciencia, Innovación y Universidades, the National Aeronautics and Space Administration, the National Astronomical Observatory of Japan, the Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the Research Council of Finland, 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 (www.euclid-ec.org). GR's research was supported

<footnote n="8">⁸ https://cosmohub.pic.es/home</footnote>
by an appointment to the NASA Postdoctoral Program administered by Oak
Ridge Associated Universities under contract with NASA. GR and AK were
supported by JPL, which is run under contract by the California Institute of
Technology for NASA (80NM0018D0004). GR acknowledges the support of
the Research Council of Finland grant 354905. The authors acknowledge the
Texas Advanced Computing Center (TACC) at The University of Texas at Austin
for providing HPC and visualization resources that have contributed to the
research results reported within this paper. URL: http://www.tacc.utexas.
edu. This project was provided with computer and storage resources by GENCI
at TGCC thanks to the grant 2023-A0150402287 on Joliot Curie's SKL parti-
tion. DFM thanks the Research Council of Norway for their support and the
resources provided by UNINETT Sigma2 – the National Infrastructure for High-
Performance Computing and Data Storage in Norway. This work has made use
of Cosmo- Hub, developed by PIC (maintained by IFAE and CIEMAT) in col-
laboration with ICE-CSIC. CosmoHub received funding from the Spanish gov-
ernment (MCIN/AEI/10.13039/501100011033), the EU NextGeneration/PRTR
(PRTRC17. 11), and the Generalitat de Catalunya. ZS acknowledges funding
from DFG project 456622116 and support from the IRAP and IN2P3 Lyon com-
puting centers. During part of this work, AMCLB was supported by a fellow-
ship of PSL University-Paris Observatory. CG thanks the support from INAF
theory Grant 2022: Illuminating Dark Matter using Weak Lensing by Clus-
ter Satellites, PI: Carlo Giocoli. VGP is supported by the Attraction de Tal-
ento Contract no. 2019-T1/TIC-12702 granted by the Comunidad de Madrid in
Spain. VGP, and by the Ministerio de Ciencia e Innovación (MICINN) under
research grant PID2021-122603NB-C21. PNG-UNITsim was run thanks to
the MareNostrum supercomputer in Spain and the Red Española de Supercom-
putación through grants: RES-AECT-2021-3-0004, RES-AECT-1-0007 & RES-
AECT-2022-3-0030. We extend our sincere gratitude to Christian Arnold and
Claudio Llinares for their valuable contributions to this research. Their work sig-
nificantly influenced the development of this project.

References

Abbott, T. M. C., Aguena, M., Alarcon, A., et al. 2022, Phys. Rev. D, 105, 023520

Adamek, J., Daverio, D., Durrer, R., & Kunz, M. 2016, JCAP, 7, 053

Adamek, J., Angulo, R. E., Arnold, C., et al. 2023, JCAP, 06, 035

Alimi, J. M., Füzfa, A., Boucher, V., et al. 2010, MNRAS, 401, 775

Amendola, L. 2000, Phys. Rev. D, 62, 043511

Amendola, L. 2004, Phys. Rev. D, 69, 103524

Angulo, R. E., & Hahn, O. 2022, Liv. Rev. Comput. Astrophys., 8, 1

Angulo, R. E., & Pontzen, A. 2016, MNRAS, 462, L1

Appel, A. W. 1985, SIAM J. Sci. Stat. Comput., 6, 85

Armendariz-Picon, C., Mukhanov, V., & Steinhardt, P. J. 2000, Phys. Rev. Lett., 85, 4438

Armendariz-Picon, C., Mukhanov, V., & Steinhardt, P. J. 2001, Phys. Rev. D, 63, 103510

Arnold, C., Leo, M., & Li, B. 2019, Nat. Astron., 3, 945

Arnold, C., Li, B., Giblin, B., Harnois-Déraps, J., & Cai, Y. C. 2022, MNRAS, 515, 4161

Avila, S., & Adame, A. G. 2023, MNRAS, 519, 3706

Baldi, M. 2023, MNRAS, 521, 613

Baldi, M., & Simpson, F. 2015, MNRAS, 449, 2239

Baldi, M., & Simpson, F. 2017, MNRAS, 465, 653

Baldi, M., Pettorino, V., Robbers, G., & Springel, V. 2010, MNRAS, 403, 1684

Baldi, M., Villaescusa-Navarro, F., Viel, M., et al. 2014, MNRAS, 440, 75

Barreira, A., Li, B., Hellwing, W. A., et al. 2014, JCAP, 4, 029

Barreira, A., Bose, S., & Li, B. 2015, JCAP, 12, 059

Barreira, A., Cabass, G., Schmidt, F., Pillepich, A., & Nelson, D. 2020, JCAP, 12, 013

Barros, B. J. 2019, Phys. Rev. D, 99, 064051

Barros, B. J., Castelão, D., da Fonseca, V., et al. 2023, JCAP, 1, 013

Bartelmann, M., & Schneider, P. 2001, Phys. Rep., 340, 291

Battat, J. B. R., Stubbs, C. W., & Chandler, J. F. 2008, Phys. Rev. D, 78, 022003

Becker, C., Arnold, C., Li, B., & Heisenberg, L. 2020, JCAP, 10, 055

Behroozi, P., Wechsler, R., & Wu, H.-Y. 2012, Rockstar: Phase-space halo finder, Astrophysics Source Code Library [record ascl:1210.0008]

Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013, ApJ, 762, 109

Berlind, A. A., Weinberg, D. H., Benson, A. J., et al. 2003, ApJ, 593, 1

Biagetti, M., Lazeyras, T., Baldauf, T., Desjacques, V., & Schmidt, F. 2017, MNRAS, 468, 3277

Bose, S., Li, B., Barreira, A., et al. 2017, JCAP, 2, 050

Bose, S., Baldi, M., & Pourtsidou, A. 2018, JCAP, 4, 032

Bowcock, P., Charmousis, C., & Gregory, R. 2000, CQG, 17, 4745

Boyle, A., Uhlemann, C., Friedrich, O., et al. 2021, MNRAS, 505, 2886

Brax, P., Davis, A.-C., Li, B., Winther, H. A., & Zhao, G.-B. 2012, JCAP, 10, 002

Brax, P., Davis, A.-C., Li, B., Winther, H. A., & Zhao, G.-B. 2013, JCAP, 4, 029
Breton, M.-A., Rasera, Y., Taruya, A., Lacombe, O., & Saga, S. 2019, MNRAS,
483, 2671

Buchdahl, H. A. 1970, MNRAS, 150, 1

Burrage, C., March, B., & Naik, A. P. 2024, JCAP, 04, 004

Caldwell, R. R., Kamionkowski, M., & Weinberg, N. N. 2003, Phys. Rev. Lett., 91, 071301

Carbone, C., Petkova, M., & Dolag, K. 2016, JCAP, 7, 034

Carretero, J., Tallada, P., Casals, J., et al. 2017, Proceedings of the European Physical Society Conference on High Energy Physics. 5-12 July, 488

Carrilho, P., Carrion, K., Bose, B., et al. 2022, MNRAS, 512, 3691

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

Castorina, E., Carbone, C., Bel, J., Sefusatti, E., & Dolag, K. 2015, JCAP, 7, 043

Cautun, M., Paillas, E., Cai, Y.-C., et al. 2018, MNRAS, 476, 3195

Chevallier, M., & Polarski, D. 2001, Int. J. Mod. Phys. D, 10, 213

Chuang, C.-H., Yepes, G., Kitaura, F.-S., et al. 2019, MNRAS, 487, 48

Cielo, M., Escudero, M., Mangano, G., & Pisanti, O. 2023, Phys. Rev. D, 108, L121301

Cole, S., Percival, W. J., Peacock, J. A., et al. 2005, MNRAS, 362, 505

Colless, M., Peterson, B. A., Jackson, C., et al. 2003, arXiv e-prints [arXiv:astro-ph/0306581]

Contarini, S., Marulli, F., Moscardini, L., et al. 2021, MNRAS, 504, 5021

Contarini, S., Verza, G., Pisani, A., et al. 2022, A&A, 667, A162

Contarini, S., Pisani, A., Hamaus, N., et al. 2023, ApJ, 953, 46

Creminelli, P., & Zaldarriaga, M. 2004, JCAP, 12, 006

Crocce, M., Pueblas, S., & Scoccimarro, R. 2006, MNRAS, 373, 369

Cuozzo, V., Carbone, C., Calabrese, M., Carella, E., & Migliaccio, M. 2024, JCAP, 4, 073

Dalal, N., Doré, O., Huterer, D., & Shirokov, A. 2008, Phys. Rev. D, 77, 123514

Das, S., Corasaniti, P. S., & Khoury, J. 2006, Phys. Rev. D, 73, 083509

de Bernardis, P., Ade, P. A. R., Bock, J. J., et al. 2000, Nature, 404, 955

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

DESI Collaboration (Aghamousa, A., et al.) 2016, arXiv e-prints [arXiv:1611.00036]

Desjacques, V., Seljak, U., & Iliev, I. T. 2009, MNRAS, 396, 85

Di Valentino, E., Mena, O., Pan, S., et al. 2021, CQG, 38, 153001

Diemer, B. 2018, ApJS, 239, 35

Dolag, K., Bartelmann, M., Perrotta, F., et al. 2004, A&A, 416, 853

Doré, O., Bock, J., Ashby, M., et al. 2014, arXiv e-prints [arXiv:1412.4872]

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

Efstathiou, G., Moody, S., Peacock, J. A., et al. 2002, MNRAS, 330, L29

Eggemeier, A., Camacho-Quevedo, B., Pezzotta, A., et al. 2023, MNRAS, 519, 2962

Eisenstein, D. J., Zehavi, I., Hogg, D. W., et al. 2005, ApJ, 633, 560

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

Euclid Collaboration (Desprez, G., et al.) 2020, A&A, 644, A31

Euclid Collaboration (Ilbert, O., et al.) 2021, A&A, 647, A117

Euclid Collaboration (Knabenhans, M., et al.) 2021, MNRAS, 505, 2840

Euclid Collaboration (Bretonnière, H., et al.) 2022, A&A, 657, A90

Euclid Collaboration (Scaramella, R., et al.) 2022, A&A, 662, A112

Euclid Collaboration (Ajani, V., et al.) 2023, A&A, 675, A120

Euclid Collaboration (Bretonnière, H., et al.) 2023, A&A, 671, A102

Euclid Collaboration (Castro, T., et al.) 2023, A&A, 671, A100

Euclid Collaboration (Merlin, E., et al.) 2023d, A&A, 671, A101

Euclid Collaboration (Koyama, K., et al.) 2024, A&A, submitted [arXiv:2409.03524]

Euclid Collaboration (Adamek, J., et al.) 2025, A&A, 695, A230

Euclid Collaboration (Mellier, Y., et al.) 2025, A&A, in press, https://doi.org/10.1051/0004-6361/202450810

Farrar, G. R., & Peebles, P. J. E. 2004, ApJ, 604, 1

Feng, Y., Chu, M.-Y., Seljak, U., & McDonald, P. 2016, MNRAS, 463, 2273

Feng, Y., Bird, S., & Lanusse, F. 2017, https://doi.org/10.5281/zenodo.1051252

Fiorini, B., Koyama, K., Izard, A., et al. 2021, JCAP, 9, 021

Fiorini, B., Koyama, K., & Baker, T. 2023, JCAP, 12, 045

García-Farieta, J. E., Marulli, F., Veropalumbo, A., et al. 2019, MNRAS, 488, 1987

Giannantonio, T., Scranton, R., Crittenden, R. G., et al. 2008, Phys. Rev. D, 77, 123520

Giocoli, C., Baldi, M., & Moscardini, L. 2018, MNRAS, 481, 2813

Gouyou Beauchamps, S., Baratta, P., Escoffier, S., et al. 2024, A&A, 693, A226

Grieb, J. N., Sánchez, A. G., Salazar-Albornoz, S., & Dalla Vecchia, C. 2016, MNRAS, 457, 1577

Guillet, T., & Teysnier, R. 2011, J. Comput. Phys., 230, 4756

Gutiérrez Adame, A., Avila, S., Gonzalez-Perez, V., et al. 2024, A&A, 689, A69

Hagstotz, S., Costanzi, M., Baldi, M., & Weller, J. 2019a, MNRAS, 486, 3927
Hagstotz, S., Gronke, M., Mota, D. F., & Baldi, M. 2019b, A&A, 629, A46

Hamaus, N., Seljak, U., & Desjacques, V. 2011, Phys. Rev. D, 84, 083509

Hamaus, N., Aubert, M., Pisani, A., et al. 2022, A&A, 658, A20

Hand, N., Feng, Y., Beutler, F., et al. 2018, AJ, 156, 160

Harnois-Dérap, J., Hernandez-Aguayo, C., Cuesta-Lazaro, C., et al. 2023, MNRAS, 525, 6336

Hassani, F., Adamek, J., Kunz, M., & Vernizzi, F. 2019, JCAP, 12, 011

Hassani, F., L'Huillier, B., Shafieloo, A., Kunz, M., & Adamek, J. 2020, JCAP, 4, 039

Hernández-Aguayo, C., Arnold, C., Li, B., & Baugh, C. M. 2021, MNRAS, 503, 3867

Hopkins, P. F. 2015, MNRAS, 450, 53

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

Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111

Jaroszynski, M., Park, C., Paczynski, B., Gott, J., & Richard, I. 1990, ApJ, 365, 22

Joseph, M., Aloni, D., Schmaltz, M., Sivarajan, E. N., & Weiner, N. 2023, Phys. Rev. D, 108, 023520

Kajita, T. 2016, Rev. Mod. Phys., 88, 030501

Klypin, A. A., & Shandarin, S. F. 1983, MNRAS, 204, 891

Klypin, A., Macciò, A. V., Mainini, R., & Bonometto, S. A. 2003, ApJ, 599, 31

Komatsu, E., & Spergel, D. N. 2001, Phys. Rev. D, 63, 063002

Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18

Kovac, J. M., Leitch, E. M., Pryke, C., et al. 2002, Nature, 420, 772

Koyama, K. 2016, Rept. Prog. Phys., 79, 046902

Kreisch, C. D., Pisani, A., Carbone, C., et al. 2019, MNRAS, 488, 4413

Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv e-prints [arXiv:1110.3193]

Le Brun, A. M. C., Arnaud, M., Pratt, G. W., & Teyssier, R. 2018, MNRAS, 473, L69

Lesgourgues, J., & Pastor, S. 2006, Phys. Rep., 429, 307

Lewis, A., & Challinor, A. 2011, CAMB: Code for Anisotropies in the Microwave Background, Astrophysics Source Code Library [record ascl:1102.026]

Li, B., Zhao, G.-B., Teyssier, R., & Koyama, K. 2012, JCAP, 1, 051

Li, B., Barreira, A., Baugh, C. M., et al. 2013a, JCAP, 11, 012

Li, B., Zhao, G.-B., & Koyama, K. 2013b, JCAP, 5, 023

Li, S.-S., Hoekstra, H., Kuijken, K., et al. 2023, A&A, 679, A133

Linden, S., & Virey, J.-M. 2008, Phys. Rev. D, 78, 023526

Linder, E. V. 2003, Phys. Rev. Lett., 90, 091301

Linder, E. V. 2008, Gen. Rel. Grav., 40, 329

Linder, E. V., & Jenkins, A. 2003, MNRAS, 346, 573

Maldacena, J. 2003, J. High Energy Phys., 2003, 013

Maltoni, M., Schwetz, T., Tórtola, M., & Valle, J. W. F. 2004, New J. Phys., 6, 122

Mangano, G., Miele, G., Pastor, S., & Peloso, M. 2002, Phys. Lett. B, 534, 8

Martinelli, M., & Tutusaus, I. 2019, Symmetry, 11, 986

Martinelli, M., Martins, C. J. A. P., Nesseris, S., et al. 2021, A&A, 654, A148

Massara, E., Villaescusa-Navarro, F., Viel, M., & Sutter, P. M. 2015, JCAP, 11, 018

Matarrese, S., & Verde, L. 2008, ApJ, 677, L77

Merson, A. I., Baugh, C. M., Helly, J. C., et al. 2013, MNRAS, 429, 556

Merten, J., Giocoli, C., Baldi, M., et al. 2019, MNRAS, 487, 104

Michaux, M., Hahn, O., Rampf, C., & Angulo, R. E. 2021, MNRAS, 500, 663

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563

Nesseris, S., Sapone, D., Martinelli, M., et al. 2022, A&A, 660, A67

Neyrinck, M. C. 2008, MNRAS, 386, 2101

Parimbelli, G., Carbone, C., Bel, J., et al. 2022, JCAP, 11, 041

Peebles, P. J. E. 1993, Principles of Physical Cosmology (Princeton University Press)

Peel, A., Pettorino, V., Giocoli, C., Starck, J.-L., & Baldi, M. 2018, A&A, 619, A38

Peel, A., Lalande, F., Starck, J.-L., et al. 2019, Phys. Rev. D, 100, 023508

Perlmutter, S., Aldering, G., Goldhaber, G., et al. 1999, ApJ, 517, 565

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

Pisani, A., Sutter, P. M., Hamaus, N., et al. 2015, Phys. Rev. D, 92, 083531

Planck Collaboration XVI. 2014, A&A, 571, A16

Planck Collaboration XIII. 2016, A&A, 594, A13

Planck Collaboration VI. 2020, A&A, 641, A6

Platen, E., van de Weygaert, R., & Jones, B. J. T. 2007, MNRAS, 380, 551

Potter, D., Stadel, J., & Teyssier, R. 2017, Comput. Astrophys. Cosmol., 4, 2

Pourtsidou, A., & Tram, T. 2016, Phys. Rev. D, 94, 043518

Pourtsidou, A., Skordis, C., & Copeland, E. J. 2013, arXiv e-prints [arXiv:1307.0458]

Press, W. H., & Schechter, P. 1974, ApJ, 187, 425

Prunet, S., Pichon, C., Aubert, D., et al. 2008, ApJS, 178, 179

Puchwein, E., Baldi, M., & Springel, V. 2013, MNRAS, 436, 348

Rácz, G., Kiessling, A., Csabai, I., & Szapudi, I. 2023, A&A, 672, A59

Rasera, Y., Breton, M. A., Corasaniti, P. S., et al. 2022, A&A, 661, A90

Riess, A. G., Filippenko, A. V., Challis, P., et al. 1998, AJ, 116, 1009

Riess, A. G., Yuan, W., Macri, L. M., et al. 2022, ApJ, 934, L7

Roncarelli, M., Carbone, C., & Moscardini, L. 2015, MNRAS, 447, 1761

Roy, F., Bouillot, V. R., & Rasera, Y. 2014, A&A, 564, A13

Ruan, C.-Z., Cuesta-Lazaro, C., Eggemeier, A., et al. 2024, MNRAS, 527, 2490

Sakr, Z. 2022, Universe, 8, 284

Salopek, D. S., & Bond, J. R. 1990, Phys. Rev. D, 42, 3936

Schaap, W. E., & van de Weygaert, R. 2000, A&A, 363, L29

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

Sefusatti, E., Crocce, M., Scoccimarro, R., & Couchman, H. M. P. 2016, MNRAS, 460, 3624

Simpson, F. 2010, Phys. Rev. D, 82, 083505

Skordis, C., Pourtsidou, A., & Copeland, E. J. 2015, Phys. Rev. D, 91, 083537

Slosar, A., Hirata, C., Seljak, U., Ho, S., & Padmanabhan, N. 2008, JCAP, 8, 031

Smith, A., de Mattia, A., Burtin, E., Chuang, C.-H., & Zhao, C. 2021, MNRAS, 500, 259

Spergel, D. N., Verde, L., Peiris, H. V., et al. 2003, ApJS, 148, 175

Spergel, D., Gehrels, N., Baltay, C., et al. 2015, arXiv e-prints [arXiv:1503.03757]

Springel, V. 2005, MNRAS, 364, 1105

Springel, V. 2010, MNRAS, 401, 791

Springel, V. 2015, N-GenIC: Cosmological structure initial conditions, Astrophysics Source Code Library [record ascl:1502.0003]

Sutter, P. M., Lavaux, G., Hamaus, N., et al. 2015, Astron. Comput., 9, 1

Takada, M., Komatsu, E., & Futamase, T. 2006, Phys. Rev. D, 73, 083520

Tallada, P., Carretero, J., Casals, J., et al. 2020, Astron. Comput., 32, 100391

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

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

Tegmark, M., Blanton, M. R., Strauss, M. A., et al. 2004, ApJ, 606, 702

Tegmark, M., Eisenstein, D. J., Strauss, M. A., et al. 2006, Phys. Rev. D, 74, 123507

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

Verza, G., Pisani, A., Carbone, C., Hamaus, N., & Guzzo, L. 2019, JCAP, 12, 040

Verza, G., Carbone, C., & Renzi, A. 2022, ApJ, 940, L16

Verza, G., Carbone, C., Pisani, A., & Renzi, A. 2023, JCAP, 12, 044

Verza, G., Carbone, C., Pisani, A., Porciani, C., & Matarrese, S. 2024, JCAP, 10, 079

Viel, M., Haehnelt, M. G., & Springel, V. 2010, JCAP, 6, 015

Weinberger, R., Springel, V., & Pakmor, R. 2020, ApJS, 248, 32

Winther, H. A., Schmidt, F., Barreira, A., et al. 2015, MNRAS, 454, 4208

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

Zeldovich, I. B. 1978, in Large Scale Structures in the Universe, eds. M. S. Longair, & J. Einasto, 79, 409

Zennaro, M., Bel, J., Villaescusa-Navarro, F., et al. 2017, MNRAS, 466, 3244

Zennaro, M., Bel, J., Dossett, J., Carbone, C., & Guzzo, L. 2018, MNRAS, 477, 491

<footnote n="1">¹ Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA</footnote>

<footnote n="2">² Department of Physics, P.O. Box 64, 00014 University of Helsinki, Finland</footnote>

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

<footnote n="4">⁴ Institut de Ciencies de l'Espai (IEEC-CSIC), Campus UAB, Carrer de Can Magrans, s/n Cerdanyola del Vallés, 08193 Barcelona, Spain</footnote>

<footnote n="5">⁵ Laboratoire Univers et Théorie, Observatoire de Paris, Université PSL, Université Paris Cité, CNRS, 92190 Meudon, France</footnote>

<footnote n="6">⁶ Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth PO1 3FX, UK</footnote>

<footnote n="7">⁷ School of Physics and Astronomy, Queen Mary University of London, Mile End Road, London E1 4NS, UK</footnote>

<footnote n="8">⁸ Institut d'Astrophysique de Paris, UMR 7095, CNRS, and Sorbonne Université, 98 bis boulevard Arago, 75014 Paris, France</footnote>

<footnote n="9">⁹ Institute of Theoretical Astrophysics, University of Oslo, P.O. Box 1029, Blindern 0315, Oslo, Norway</footnote>

<footnote n="10">¹⁰ Institut für Theoretische Physik, University of Heidelberg, Philosophenweg 16, 69120 Heidelberg, Germany</footnote>

<footnote n="11">¹¹ Institut de Recherche en Astrophysique et Planétologie (IRAP), Université de Toulouse, CNRS, UPS, CNES, 14 Av. Edouard Belin, 31400 Toulouse, France</footnote>
12 Université St Joseph; Faculty of Sciences, Beirut, Lebanon

13 Dipartimento di Fisica “G. Occhiali”, Università degli Studi di Milano Bicocca, Piazza della Scienza 3, 20126 Milano, Italy

14 INAF-Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Piero Gobetti 93/3, 40129 Bologna, Italy

15 IFPU, Institute for Fundamental Physics of the Universe, via Beirut 2, 34151 Trieste, Italy

16 Dipartimento di Fisica e Astronomia “Augusto Righi” – Alma Mater Studiorum Università di Bologna, via Piero Gobetti 93/2, 40129 Bologna, Italy

17 ICSC – Centro Nazionale di Ricerca in High Performance Computing, Big Data and Quantum Computing, Via Magnanelli 2, Bologna, Italy

18 Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, University of Manchester, Oxford Road, Manchester M13 9PL, UK

19 INAF-IASF Milano, Via Alfonso Corti 12, 20133 Milano, Italy

20 Dipartimento di Fisica “Aldo Pontremoli”, Università degli Studi di Milano, Via Celoria 16, 20133 Milano, Italy

21 INFN Gruppo Collegato di Parma, Viale delle Scienze 7/A, 43124 Parma, Italy

22 SISSA, International School for Advanced Studies, Via Bonomea 265, 34136 Trieste, TS, Italy

23 International Centre for Theoretical Physics (ICTP), Strada Costiera 11, 34151 Trieste, Italy

24 Center for Cosmology and Particle Physics, Department of Physics, New York University, New York, NY 10003, USA

25 Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, 10010 New York, NY, USA

26 Dipartimento di Fisica e Astronomia, Università di Bologna, Via Gobetti 93/2, 40129 Bologna, Italy

27 INFN-Sezione di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy

28 Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK

29 Higgs Centre for Theoretical Physics, School of Physics and Astronomy, The University of Edinburgh, Edinburgh EH9 3FD, UK

30 Institut de Physique Théorique, CEA, CNRS, Université Paris-Saclay, 91191 Gif-sur-Yvette Cedex, France

31 Departamento de Física Teórica, Facultad de Ciencias, Universidad Autónoma de Madrid, 28049 Cantoblanco, Madrid, Spain

32 Instituto de Física Teórica UAM-CSIC, Campus de Cantoblanco, 28049 Madrid, Spain

33 Centro de Investigación Avanzada en Física Fundamental (CIAFF), Facultad de Ciencias, Universidad Autónoma de Madrid, 28049 Madrid, Spain

34 Department of Astrophysics, University of Zurich, Winterthur-strasse 190, 8057 Zurich, Switzerland

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

36 Istituto Nazionale di Fisica Nucleare, Sezione di Bologna, Via Irnerio 46, 40126 Bologna, Italy

37 Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, 85748 Garching, Germany

38 Université de Genève, Département de Physique Théorique and Centre for Astroparticle Physics, 24 quai Ernest-Ansermet, CH-1211 Genève 4, Switzerland

39 Department of Physics, Institute for Computational Cosmology, Durham University, South Road DH1 3LE, UK

40 Institut universitaire de France (IUF), 1 rue Descartes, 75231 Paris Cedex 05, France

41 Universitäts-Sternwarte München, Fakultät für Physik, Ludwig-Maximilians-Universität München, Scheinerstrasse 1, 81679 München, Germany

42 Excellence Cluster ORIGINS, Boltzmannstrasse 2, 85748 Garching, Germany

43 Dipartimento di Fisica e Astronomia “G. Galilei”, Università di Padova, Via Marzolo 8, 35131 Padova, Italy

44 INFN-Padova, Via Marzolo 8, 35131 Padova, Italy

45 INAF-Osservatorio Astronomico di Trieste, Via G. B. Tiepolo 11, 34143 Trieste, Italy

46 INFN, Sezione di Trieste, Via Valerio 2, 34127 Trieste, TS, Italy

47 Dipartimento di Scienze Matematiche, Fisiche e Informatiche, Università di Parma, Viale delle Scienze 7/A, 43124 Parma, Italy

48 Aix-Marseille Université, CNRS/IN2P3, CPPM, Marseille, France

49 Department of Astrophysical Sciences, Peyton Hall, Princeton University, Princeton, NJ 08544, USA

50 The Cooper Union for the Advancement of Science and Art, 41 Cooper Square, New York, NY 10003, USA

51 Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu, HI 96822, USA

52 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Avenida Complutense 40, 28040 Madrid, Spain

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

54 Université Paris-Saclay, CNRS, Institut d'astrophysique spatiale, 91405 Orsay, France

55 INAF-Osservatorio Astronomico di Brera, Via Brera 28, 20122 Milano, Italy

56 INAF-Osservatorio Astrofisico di Torino, Via Osservatorio 20, 10025 Pino Torinese, (TO), Italy

57 Dipartimento di Fisica, Università di Genova, Via Dodecaneso 33, 16146 Genova, Italy

58 INFN-Sezione di Genova, Via Dodecaneso 33, 16146 Genova, Italy

59 Department of Physics “E. Pancini”, University Federico II, Via Cinthia 6, 80126 Napoli, Italy

60 INAF-Osservatorio Astronomico di Capodimonte, Via Moiariello 16, 80131 Napoli, Italy

61 INFN section of Naples, Via Cinthia 6, 80126 Napoli, Italy

62 Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, PT4150-762 Porto, Portugal

63 Faculdade de Ciências da Universidade do Porto, Rua do Campo de Alegre, 4150-007 Porto, Portugal

64 Dipartimento di Fisica, Università degli Studi di Torino, Via P. Giuria 1, 10125 Torino, Italy

65 INFN-Sezione di Torino, Via P. Giuria 1, 10125 Torino, Italy

66 INAF-Osservatorio Astronomico di Roma, Via Frascati 33, 00078 Monteporzio Catone, Italy

67 INFN-Sezione di Roma, Piazzale Aldo Moro, 2 – c/o Dipartimento di Fisica, Edificio G. Marconi, 00185 Roma, Italy

68 Institute for Theoretical Particle Physics and Cosmology (TTK), RWTH Aachen University, 52056 Aachen, Germany

69 Dipartimento di Fisica e Astronomia “Augusto Righi” – Alma Mater Studiorum Università di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy

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

71 European Space Agency/ESRIN, Largo Galileo Galilei 1, 00044 Frascati, Roma, Italy

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

73 Université Claude Bernard Lyon 1, CNRS/IN2P3, IP2I Lyon, UMR 5822, Villeurbanne F-69100, France

74 Institute of Physics, Laboratory of Astrophysics, Ecole Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, 1290 Versoix, Switzerland

75 UCB Lyon 1, CNRS/IN2P3, IUF, IP2I Lyon, 4 rue Enrico Fermi, 69622 Villeurbanne, France

76 Departamento de Física, Faculdade de Ciências, Universidade de Lisboa, Edifício C8, Campo Grande PT1749-016 Lisboa, Portugal

77 Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências, Universidade de Lisboa, Campo Grande 1749-016, Lisboa, Portugal

78 Department of Astronomy, University of Geneva, ch. d'Ecogia 16, 1290 Versoix, Switzerland

79 INAF-Istituto di Astrofisica e Planetologia Spaziali, via del Fosso del Cavaliere, 100, 00100 Roma, Italy
80 Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191 Gif-sur-Yvette, France

81 Institut d'Estudis Espacials de Catalunya (IEEC), Edifici RDIT, Campus UPC, 08860 Castelldefels, Barcelona, Spain

82 FRACTAL S.L.N.E., calle Tulipán 2, Portal 13 1A, 28231 Las Rozas de Madrid, Spain

83 INAF-Osservatorio Astronomico di Padova, Via dell'Osservatorio 5, 35122 Padova, Italy

84 Max Planck Institute for Extraterrestrial Physics, Giessenbachstr. 1, 85748 Garching, Germany

85 Felix Hormuth Engineering, Goethestr. 17, 69181 Leimen, Germany

86 Technical University of Denmark, Elektrovej 327, 2800 Kgs. Lyn-gby, Denmark

87 Cosmic Dawn Center (DAWN), Denmark

88 Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France

89 Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany

90 NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

91 Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK

92 Department of Physics and Helsinki Institute of Physics, Gustaf Hällströmin katu 2, 00014 University of Helsinki, Finland

93 Mullard Space Science Laboratory, University College London, Holmbury St Mary, Dorking, Surrey RH5 6NT, UK

94 Helsinki Institute of Physics, Gustaf Hällströmin katu 2, University of Helsinki, Helsinki, Finland

95 NOVA optical infrared instrumentation group at ASTRON, Oude Hoogeveensedijk 4, 7991PD Dwingeloo, The Netherlands

96 Centre de Calcul de l'IN2P3/CNRS, 21 avenue Pierre de Coubertin, 69627 Villeurbanne Cedex, France

97 Universität Bonn, Argelander-Institut für Astronomie, Auf dem Hügel 71, 53121 Bonn, Germany

98 Aix-Marseille Université, CNRS, CNES, LAM, Marseille, France

99 Université Paris Cité, CNRS, Astroparticule et Cosmologie, 75013 Paris, France

100 Institut d'Astrophysique de Paris, 98bis Boulevard Arago, 75014 Paris, France

101 European Space Agency/ESTEC, Keplerlaan 1, 2201, AZ Noordwijk, The Netherlands

102 Department of Physics and Astronomy, University of Aarhus, Ny Munkegade 120, DK-8000 Aarhus C, Denmark

103 Waterloo Centre for Astrophysics, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada

104 Department of Physics and Astronomy, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada

105 Perimeter Institute for Theoretical Physics, Waterloo, Ontario N2L 2Y5, Canada

106 Space Science Data Center, Italian Space Agency, via del Politecnico snc, 00133 Roma, Italy

107 Centre National d'Etudes Spatiales – Centre spatial de Toulouse, 18 avenue Edouard Belin, 31401 Toulouse Cedex 9, France

108 Institute of Space Science, Str. Atomistilor, nr. 409 Măgurele, Ilfov 077125, Romania

109 Departamento de Astrofísica, Universidad de La Laguna, 38206, La Laguna, Tenerife, Spain

110 Departamento de Física, FCFM, Universidad de Chile, Blanco Encalada 2008, Santiago, Chile

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

112 Satlantis, University Science Park, Sede Bld 48940, Leioa-Bilbao, Spain

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

114 Universidad Politécnica de Cartagena, Departamento de Eléctronica y Tecnología de Computadoras, Plaza del Hospital 1, 30202 Cartagena, Spain

115 INFN-Bologna, Via Irnerio 46, 40126 Bologna, Italy

116 Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700 AV Groningen, The Netherlands

117 Infrared Processing and Analysis Center, California Institute of Technology, Pasadena, CA 91125, USA

118 INAF, Istituto di Radioastronomia, Via Piero Gobetti 101, 40129 Bologna, Italy

119 Astronomical Observatory of the Autonomous Region of the Aosta Valley (OAVdA), Loc. Lignan 39, I-11020 Nus (Aosta Valley), Italy

120 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK

121 School of Physics and Astronomy, Cardiff University, The Parade, Cardiff CF24 3AA, UK

122 Junia, EPA department, 41 Bd Vauban, 59800 Lille, France

123 CERCA/ISO, Department of Physics, Case Western Reserve University, 10900 Euclid Avenue, Cleveland, OH 44106, USA

124 INFN-Sezione di Milano, Via Celoria 16, 20133 Milano, Italy

125 Departamento de Física Fundamental. Universidad de Salamanca, Plaza de la Merced s/n, 37008 Salamanca, Spain

126 Dipartimento di Fisica e Scienze della Terra, Università degli Studi di Ferrara, Via Giuseppe Saragat 1, 44122 Ferrara, Italy

127 Istituto Nazionale di Fisica Nucleare, Sezione di Ferrara, Via Giuseppe Saragat 1, 44122 Ferrara, Italy

128 Kavli Institute for the Physics and Mathematics of the Universe (WPI), University of Tokyo, Kashiwa, Chiba 277-8583, Japan

129 Dipartimento di Fisica – Sezione di Astronomia, Università di Trieste, Via Tiepolo 11, 34131 Trieste, Italy

130 Minnesota Institute for Astrophysics, University of Minnesota, 116 Church St SE, Minneapolis, MN 55455, USA

131 Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Bd de l'Observatoire, CS 34229, 06304 Nice Cedex 4, France

132 Department of Physics & Astronomy, University of California Irvine, Irvine, CA 92697, USA

133 Department of Astronomy & Physics and Institute for Computational Astrophysics, Saint Mary's University, 923 Robie Street, Halifax, Nova Scotia B3H 3C3, Canada

134 Departamento Física Aplicada, Universidad Politécnica de Cartagena, Campus Muralla del Mar, 30202 Cartagena, Murcia, Spain

135 Instituto de Astrofísica de Canarias (IAC); Departamento de Astrofísica, Universidad de La Laguna (ULL), 38200 La Laguna, Tenerife, Spain

136 Department of Physics, Oxford University, Keble Road, Oxford OX1 3RH, UK

137 Department of Computer Science, Aalto University, PO Box 15400, Espoo FI-00 076, Finland

138 Instituto de Astrofísica de Canarias, c/ Via Lactea s/n, La Laguna E-38200, Spain. Departamento de Astrofísica de la Universidad de La Laguna, Avda. Francisco Sanchez, La Laguna E-38200, Spain

139 Ruhr University Bochum, Faculty of Physics and Astronomy, Astronomical Institute (AIRUB), German Centre for Cosmological Lensing (GCCL), 44780 Bochum, Germany

140 DARK, Niels Bohr Institute, University of Copenhagen, Jagtvej 155, 2200 Copenhagen, Denmark

141 Univ. Grenoble Alpes, CNRS, Grenoble INP, LPSC-IN2P3, 53, Avenue des Martyrs, 38000 Grenoble, France

142 Department of Physics and Astronomy, Vesilinnantie 5, 20014 University of Turku, Finland

143 Serco for European Space Agency (ESA), Camino bajo del Castillo, s/n, Urbanizacion Villafranca del Castillo, Villanueva de la Cañada 28692, Madrid, Spain

144 ARC Centre of Excellence for Dark Matter Particle Physics, Melbourne, Australia

145 Centre for Astrophysics & Supercomputing, Swinburne University of Technology, Hawthorn, Victoria 3122, Australia

146 Department of Physics and Astronomy, University of the Western Cape, Bellville, Cape Town 7535, South Africa
147 Université Libre de Bruxelles (ULB), Service de Physique Théorique CP225, Boulevard du Triophe, 1050 Bruxelles, Belgium

148 ICTP South American Institute for Fundamental Research, Instituto de Física Teórica, Universidade Estadual Paulista, São Paulo, Brazil

149 IRFU, CEA, Université Paris-Saclay, 91191 Gif-sur-Yvette Cedex, France

150 Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, Stockholm SE-106 91, Sweden

151 Astrophysics Group, Blackett Laboratory, Imperial College London, London SW7 2AZ, UK

152 INAF-Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy

153 Dipartimento di Fisica, Sapienza Università di Roma, Piazzale Aldo Moro 2, 00185 Roma, Italy

154 Centro de Astrofísica da Universidade do Porto, Rua das Estrelas, 4150-762 Porto, Portugal

155 HE Space for European Space Agency (ESA), Camino bajo del Castillo, s/n, Urbanizacion Villafranca del Castillo, Villanueva de la Cañada 28692, Madrid, Spain

156 Aurora Technology for European Space Agency (ESA), Camino bajo del Castillo, s/n, Urbanizacion Villafranca del Castillo, Villanueva de la Cañada 28692, Madrid, Spain

157 Dipartimento di Fisica, Università degli studi di Genova, and INFN-Sezione di Genova, via Dodecaneso 33, 16146 Genova, Italy

158 Theoretical astrophysics, Department of Physics and Astronomy, Uppsala University, Box 515, 751 20 Uppsala, Sweden

159 Institute Lorentz, Leiden University, Niels Bohrweg 2, 2333 CA Leiden, The Netherlands

160 Department of Physics, Royal Holloway, University of London, London TW20 0EX, UK

161 Cosmic Dawn Center (DAWN)

162 Niels Bohr Institute, University of Copenhagen, Jagtvej 128, 2200 Copenhagen, Denmark