EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 695 (March 2025) A&A, 695 (2025) A232 Full HTML
Open Access
Issue
A&A
Volume 695, March 2025
Article Number A232
Number of page(s) 22
Section Cosmology (including clusters of galaxies)
DOI https://doi.org/10.1051/0004-6361/202452185
Published online 27 March 2025

© The Authors 2025

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

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 H0 = 67.7  ±  0.4 km s−1 Mpc−1 (Planck Collaboration VI 2020), local measurements, often based on the observations of supernovae in nearby galaxies, suggest a higher value of H0 = 73.0  ±  1.0 km s−1Mpc−1 (Riess et al. 2022). This 5σ discrepancy is called the Hubble tension. A similar tension has been identified in the S 8 = σ 8 Ω m / 0.3 $ S_8=\sigma_8\,\sqrt{\Omega_{\mathrm{m}}/0.3} $ parameter, which combines the amplitude of linear-matter density fluctuations on the 8 h−1 Mpc scale, σ8, and the cosmic matter density, Ωm. Measurements derived from the CMB (Planck Collaboration VI 2020) appear to yield a value of S8 2.9σ 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 S8 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 simulations. 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 ≡ pde/(ρde c2), where pde and ρ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 Linder (Linder 2003) proposed a simple parameterisation of

w de ( z ) = w 0 + w a z 1 + z = w 0 + w a ( 1 − a ) , $$ \begin{aligned} w_{\mathrm{de}}(z) = w_0 + w_a \frac{z}{1+z} = w_0 + w_a \left(1-a\right) , \end{aligned} $$(1)

where the w0 parameter represents the value of the equation-of-state at the present time, and wa 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 wde(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 wde(z) at late times. Despite its simple form, it shows a wide range of interesting properties (Linder 2008; Linden & Virey 2008). The cosmological constant corresponds to w0 = −1 and wa = 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 = ∫ d 4 x − g p ( ϕ , X ) , $$ \begin{aligned} S=\int \mathrm{d}^4 x\;\sqrt{-g}\;p(\phi ,X), \end{aligned} $$(2)

where X = (1/2)gμν∇μϕ∇νϕ. The energy density of the scalar field is given by

u ϕ = 2 X d p d X − p , $$ \begin{aligned} {u}_{\phi }=2X\;\frac{\mathrm{d}p}{\mathrm{d} X}-p, \end{aligned} $$(3)

and the pressure is pϕ = p(ϕ, X). This pressure gives an effective fluid equation-of-state parameter as

w ϕ = p ϕ u ϕ = − p p − 2 X p , X , $$ \begin{aligned} w_{\phi }=\frac{p_{\phi }}{{u}_{\phi }}=-\frac{p}{p-2X p,_{X}}\;, \end{aligned} $$(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 = p , X p , X + 2 X p , XX . $$ \begin{aligned} c_{\rm s}^2 =\frac{p,_{X}}{p,_X+2Xp,_{XX}} \;. \end{aligned} $$(5)

The k-essence field satisfies the continuity equation

u ˙ ϕ = − 3 H ( u ϕ + p ϕ ) , $$ \begin{aligned} \dot{{u}}_{\phi }=-3H\,({u}_{\phi }+p_{\phi }), \end{aligned} $$(6)

which results in the scalar equation of motion

G μ ν ∇ μ ∇ ν ϕ + 2 X ∂ 2 p ∂ X ∂ ϕ − ∂ p ∂ ϕ = 0 , $$ \begin{aligned} G^{\mu \nu }\nabla _{\mu }\nabla _{\nu }\phi +2X\frac{\partial ^2p}{\partial X \partial \phi }-\frac{\partial p}{\partial \phi }=0, \end{aligned} $$(7)

where

G μ ν = ∂ p ∂ X g μ ν + ∂ 2 p ∂ X 2 ∇ μ ϕ ∇ ν ϕ . $$ \begin{aligned} G^{\mu \nu }=\frac{\partial p}{\partial X}g^{\mu \nu }+\frac{\partial ^2 p}{\partial X^2}\nabla ^{\mu }\phi \nabla ^{\nu }\phi \;. \end{aligned} $$(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μν, conserved:

∇ μ T ν ( c ) μ = C ν ( ϕ ) = − ∇ μ T ν ( ϕ ) μ , $$ \begin{aligned} \nabla _{\mu }T^{(c)\mu }_{\nu } = C_{\nu }(\phi ) = - \nabla _{\mu }T^{(\phi )\mu }_{\nu } , \end{aligned} $$(9)

where Cν(ϕ) is a conformal coupling function expressed in the form:

C ν ( ϕ ) = κ β ( ϕ ) u c ∇ ν ϕ , $$ \begin{aligned} C_{\nu }(\phi ) = \kappa \, \beta (\phi )\,{u} _{\rm c}\nabla _{\nu }\phi , \end{aligned} $$(10)

where κ ≡ 8 π G N c 2 $ \kappa \equiv \frac{8\pi G_{\mathrm{N}}}{c^2} $, GN is Newton’s gravitational constant, uc is the cold dark matter energy density in the IDE model1, and β(ϕ) is a coupling function. The dark energy scalar field, ϕ, has an intrinsic energy density and pressure given by

u ϕ = 1 2 g μ ν ∂ μ ϕ ∂ ν ϕ + V ( ϕ ) , $$ \begin{aligned} {u} _{\phi }&= \frac{1}{2}\,g^{\mu \nu }\,\partial _{\mu }\phi \,\partial _{\nu }\phi + V(\phi ),\end{aligned} $$(11)

p ϕ = 1 2 g μ ν ∂ μ ϕ ∂ ν ϕ − V ( ϕ ) , $$ \begin{aligned} p_{\phi }&= \frac{1}{2}\,g^{\mu \nu }\,\partial _{\mu }\phi \,\partial _{\nu }\phi - V(\phi ), \end{aligned} $$(12)

where V(ϕ) 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 β(ϕ) = β:

ϕ ¨ + 3 H ϕ ˙ + d V d ϕ = κ β u c , $$ \begin{aligned} \ddot{\phi } + 3H\dot{\phi } + \frac{\mathrm{d} V}{\mathrm{d}\phi }&=\kappa \, \beta \, {u} _{\rm c} ,\end{aligned} $$(13)

u ˙ c + 3 H u c = − κ β u c . $$ \begin{aligned} \dot{{u} }_{\rm c} + 3H{u} _{\rm c}&=-\kappa \,\beta \, {u} _{\rm c}. \end{aligned} $$(14)

In the standard approach, a theoretically motivated analytical form for the self-interaction potential function, V(ϕ), 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 ΛCDM background expansion history by setting

H 2 = H Λ CDM 2 , $$ \begin{aligned} H^{2}=H^{2}_{\Lambda \mathrm{CDM}} , \end{aligned} $$(15)

where HΛCDM is the standard Hubble function defined by

H Λ CDM 2 = 8 π G N 3 ( ρ r + ρ b + ρ CDM + ρ Λ ) , $$ \begin{aligned} H^{2}_{\Lambda \mathrm{CDM}} = \frac{8 \pi G_{\rm N}}{3}(\rho _{\rm r} + \rho _{\rm b} + \rho _{\rm CDM} + \rho _{\Lambda }), \end{aligned} $$(16)

where ρr, ρb, ρCDM, and ρΛ are the mass densities of the radiation, baryon, CDM, and Λ components of the background ΛCDM model. This will determine an effective potential, V(ϕ), according to the resulting evolution of the scalar field, ϕ. 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 ϕ = ρ CDM c 2 + ρ Λ c 2 − u c , $$ \begin{aligned} {u} _{\phi }&= \rho _{\rm CDM}\,c^2 + \rho _{\Lambda }\,c^2 - {u} _{\rm c}, \end{aligned} $$(17)

p ϕ = p Λ = − u Λ , $$ \begin{aligned} p_{\phi }&=p_{\Lambda }=-{u} _{\Lambda } , \end{aligned} $$(18)

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

ϕ ˙ 2 = ρ CDM c 2 − u ϕ . $$ \begin{aligned} \dot{\phi }^{2} = \rho _{\rm CDM}\,c^2 - {u} _{\phi }. \end{aligned} $$(19)

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

V ( ϕ ) = 1 2 ϕ ˙ 2 + ρ Λ c 2 , $$ \begin{aligned} V(\phi ) = \frac{1}{2}\dot{\phi }^{2} + \rho _{\Lambda }\,c^2, \end{aligned} $$(20)

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

2 ϕ ¨ + 3 H ϕ ˙ − κ β u c = 0 , $$ \begin{aligned} 2\ddot{\phi } + 3H\dot{\phi } -\kappa \beta {u} _{c} = 0 , \end{aligned} $$(21)

which can be numerically solved for the dynamical evolution of the system. With this choice, the β 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 σ8 tension, but that CMB prefers the ΛCDM limit. In particular, they find that the CMB constrains |β|≲0.02, RSD constraints |β|≲0.10, while weak lensing data from the Kilo-Degree Survey actually prefers a non-zero value |β|∼0.1.

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, γ. This model’s action is

S = c 4 16 π G 5 ∫ M d 5 x − γ R 5 + ∫ ∂ M d 4 x − g ( c 4 16 π G N R + L m ) , $$ \begin{aligned} S = \frac{c^4}{16\pi G_5}\int _\mathcal{M} \mathrm{d}^5 x \sqrt{-\gamma }\,R_5 + \int _{\partial \mathcal{M}} \mathrm{d}^4 x \,\sqrt{-g}\, \left(\frac{c^4}{16\pi G_{\rm N}}\,R + \mathcal{L}_{\rm m} \right), \end{aligned} $$(22)

where G5 and GN are the five- and four-dimensional Newton’s constants, respectively, and ℒ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, rc = G5/(2GN), from which we construct the dimensionless parameter Ωrc ≡ c2/(4rc2H02). The modified Friedmann equation on the brane (Deffayet 2001) becomes

H 2 = ± c H r c + 8 π G N 3 ρ ¯ . $$ \begin{aligned} H^2 = \pm \, c \frac{H}{r_{\rm c}} + \frac{8 \pi G_{\rm N}}{3} \bar{\rho }. \end{aligned} $$(23)

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

ρ de ( a ) = ρ cr , 0 ( Ω Λ + 2 Ω rc Ω Λ + Ω m a − 3 ) , $$ \begin{aligned} \rho _{\mathrm{de}}(a) = \rho _{\rm cr,0} \, \left( {\Omega _\Lambda } + 2 \sqrt{\Omega _{\rm rc}} \sqrt{{\Omega _\Lambda } + \Omega _\mathrm{m} \, a^{-3}} \, \right), \end{aligned} $$(24)

where ρcr, 0 is the critical density. The observational constraints on the model require the cross-over scale, rc, to be larger than the size of the horizon H0−1 today. For example, Solar System constraints require rcH0 ≳ 1.6 (Battat et al. 2008), and galaxy clustering in the BOSS survey constraints rcH0 ≳ 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 = c 4 16 π G N ∫ d 4 x − g [ R + f ( R ) ] , $$ \begin{aligned} S = \frac{c^4}{16\pi G_{\rm N}} \int {\mathrm{d}^4 x \,\sqrt{-g}\, \left[\,R+f(R)\,\right]} , \end{aligned} $$(25)

where gμν 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 fR = df/dR ≪ 1 we have

f ( R ) = − 6 Ω Λ H 0 2 c 2 + | f R 0 | R ¯ 0 2 R , $$ \begin{aligned} f(R) = - 6 {\Omega _\Lambda }\,\frac{H_0^2}{c^2} + |f_{R0}|\,\frac{\bar{R}_0^2}{R}, \end{aligned} $$(26)

where fR0 is the free parameter of the model, R ¯ 0 $ \bar R_0 $ is the Ricci scalar evaluated at background at present time, H0 is the Hubble constant, and ΩΛ is the energy-density parameter of the cosmological constant. |fR0| characterises the magnitude of the deviation from ΛCDM, with smaller values corresponding to weaker departures from general relativity until we recover ΛCDM in the limit of fR0 → 0, but for the small |fR0| values still allowed by observations, the background expansion history approximates that of ΛCDM and

R ¯ 0 = 3 Ω m H 0 2 c 2 ( 1 + 4 Ω Λ Ω m ) , $$ \begin{aligned} \bar{R}_0 = 3 \Omega _\mathrm{m} \frac{H_0^2}{c^2} \left(1+ 4 \frac{{\Omega _\Lambda }}{\Omega _\mathrm{m} } \right), \end{aligned} $$(27)

with the matter energy density parameter Ωm = 1 − ΩΛ. 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 |fR0| vary from |fR0|≲10−6 in the Solar System, |fR0|≲10−8 from galaxy scales (Burrage et al. 2024) to |fR0|≲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ν, and the number of neutrino species, Neff. More in general, Neff parametrises the contribution of relativistic species to the background density of radiation, ρr, as

ρ r = [ 1 + 7 8 ( 4 11 ) 4 / 3 N eff ] ρ γ , $$ \begin{aligned} \rho _{\rm r} = \left[1 + \frac{7}{8}\left(\frac{4}{11}\right)^{4/3}{N_{\rm eff}}\right]\rho _{\gamma }, \end{aligned} $$(28)

where ργ is the photon background density. In the standard model, Neff 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 Neff 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 ≲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 Ωmh2 (Lesgourgues & Pastor 2006). Massive neutrinos act as non-relativistic particles on scales k > knr = 0.018(mν/1eV)1/2Ωm1/2  h−1 Mpc, where knr is the wavenumber corresponding to the Hubble horizon size at the epoch znr when the given neutrino species becomes non-relativistic following 1 + z nr ≃ 1900 ( m ν 1 eV ) $ 1+z_{\mathrm{nr}} \simeq 1900 \left( \frac{m_\nu}{\mathrm{{1 eV}}} \right) $, Ωm is the matter density parameter, and h = H0/100 km s−1 Mpc−1. The large velocity dispersion of non-relativistic neutrinos suppresses the formation of neutrino perturbations in a way that depends on mν and redshift z, leaving an imprint on the matter power spectrum at scales k > kfs(z), with

k fs = 0.82 H ( z ) H 0 ( 1 + z ) 2 ( m ν 1 eV ) h Mpc − 1 , $$ \begin{aligned} k_{\rm fs}=\frac{0.82 H(z)}{H_0(1+z)^2} \left( \frac{m_\nu }{\mathrm{{1 eV}}} \right)\, h\,\mathrm{Mpc} ^{-1}, \end{aligned} $$(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 ϕ is given by

ϕ ( x ) = ϕ G ( x ) + f NL local ( ϕ 2 ( x ) − ⟨ ϕ 2 ( x ) ⟩ ) , $$ \begin{aligned} \phi (\mathbf x ) = \phi _G(\mathbf x ) + f_{\rm NL}^\mathrm{local} (\phi ^2(\mathbf x ) - \langle \phi ^2(\mathbf x ) \rangle ), \end{aligned} $$(30)

where ϕG(x) is the Gaussian potential, while ϕ is the non-Gaussian potential. f NL local $ f_{\mathrm{NL}}^{\mathrm{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

δ ( k , z ) = α ( k , z ) ϕ ( k , z ) , $$ \begin{aligned} \delta (k,z)\; = \;\alpha (k,z)\;\phi (k,z), \end{aligned} $$(31)

where

α ( k , z ) = 2 D ( z ) 3 Ω m c 2 H 0 2 g ( 0 ) g ( z rad ) k 2 T ( k ) , $$ \begin{aligned} \alpha (k,z) = \frac{2 D(z)}{3 \Omega _\mathrm{m} } \frac{c^2}{H_0^2} \frac{g(0)}{g(z_{\rm rad})} k^2 T(k), \end{aligned} $$(32)

T(k) is the transfer function normalised at T(k → 0) = 1, and D(z) is the growth factor normalised at D(z = 0) = 1. The factor g(0)/g(zrad), 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)∝1/(1 + z) during matter-domination. This factor is g ( z rad ) g ( 0 ) ∼ 1.3 $ \frac{g(z_{\mathrm{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 t , t ( k , z ) = [ b 1 + b ϕ f NL local α ( k , z ) ] 2 P m , m ( k , z ) , $$ \begin{aligned} P_{\rm t,t}(k,z) = \left[b_1 + \frac{b_\phi f_{\rm NL}^\mathrm{local}}{\alpha (k,z)}\right]^2 P_{\rm m,m}(k,z), \end{aligned} $$(33)

where Pt, t(k, z) is the power spectrum of the tracer, Pm, m(k, z) is the power spectrum of the matter, b1 is the linear bias, and bϕ is the response of the tracer to the presence of the local-PNG. Now, Pt, t(k, z) has a dependency with k which scales as k−2 at leading order due to the α(k, z) term. The bϕ is usually parametrised as

b ϕ = 2 δ c ( b 1 − p ) . $$ \begin{aligned} b_\phi = 2 \delta _c (b_1 -p). \end{aligned} $$(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.

Table 1.

Overview of the simulation suites analysed for this project.

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: Ωm = 1 − ΩΛ = 0.3111, Ωb = 0.04897, H0 = 67.66 km s−1 Mpc−1, ns = 0.9665, and σ8 = 0.8102. The wCDM pair had the following parameters: w0 = −1.04, Ωm = 0.3096, Ωb = 0.04899, Ωde = 0.6904, H0 = 67.66 km s−1 Mpc−1, ns = 0.9331, and σ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 21603 dark matter particles in a (1.5 h−1 Gpc)3 volume, with ε = 13.8 h−1 kpc softening length. The initial conditions (ICs) were generated by a modified version of the N-GenIC code (Springel 2015) by using the Zeldovich approximation and initial linear power spectra from the Boltzmann code CAMB (Lewis & Challinor 2011). The simulations started from redshift zinit = 127, with a total of 48 output times. In this project, 31 particle snapshots were analysed in the 0.5 ≤ z ≤ 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 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): Ωb = 0.05, Ωm = 0.32, H0 = 67.0 km s−1 Mpc−1, ns = 0.96, and As = 2.127 × 10−9. Given these values, the reference (i.e. the massless neutrino case) CDM-particle mass resolution is m CDM p = 8.27 × 10 10 h − 1 M ⊙ $ m^{\mathrm{p}}_{\mathrm{CDM}} = 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 Ω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 Ω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 ε = 20 h−1 kpc, a comoving volume of 8 h−3 Gpc3 filled with 20483 dark matter particles and, when present, 20483 neutrino particles. The simulations are initialised at zinit = 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 40963 particles, which results in a smoothing scale of 5 h−1 kpc at the maximum refinement level. Both simulations share the parameters H0 = 72.0 km s−1Mpc−1, ns = 0.963, Ωb = 0.04356 and Ωr = 8.076 × 10−5. The flat ΛCDM simulation has a WMAP7 cosmology (Komatsu et al. 2011): Ωm = 0.25733, and σ8 = 0.80101, while the flat wCDM simulation is consistent at the 1σ-level with a WMAP7 cosmology with Ωm = 0.27508, σ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 zinit ∼ 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 ANalaytic 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: Ωm = 0.319, Ωb = 0.049, ΩΛ = 0.681, H0 = 67.0 km s−1Mpc−1, As = 2.1 × 10−9, and ns = 0.96. The nDGP simulations used the same parameters as the ΛCDM simulations with the cross-over scale rc = 1.2 c/H0. All of the simulations in this simulation suite had a box size of 1024 h−1 Mpc and 10243 particles. The initial conditions were generated at zinit = 49 with 2LPT using the FML code2 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 H IRES 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 Mpc, with 20483 particles. The base cosmological parameters of the simulations are the Planck 2015 parameters (Planck Collaboration XIII 2016): Ωm = 1 − ΩΛ = 0.3089, Ωb = 0.0486, H0 = 67.74 km s−1 Mpc−1, ns = 0.9667, and σ8 = 0.8159. This simulation series focusses on nDGP gravity and tested four cases: rc = {0.5,  1,  2,  5} c/H0. The series contains paired-and-fixed simulations (Angulo & Pontzen 2016) to suppress cosmic variance in ΛCDM and in the nDGP model for rc = 1 c/H0, while for the others they were only run for a single fixed amplitude realisation. The initial conditions were generated at zinit = 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 (fR0, mν) 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ν 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 Mpc per side, used a softening length of ε = 20 h−1 kpc, and include (2 × )7683 particles (for the CDM and neutrinos components). The cosmological parameters (for the reference ΛCDM cosmology with massless neutrinos) have been set to Ωm = 1 − ΩΛ = 0.31345, σ8 = 0.842, H0 = 67.31 km s−1 Mpc−1, ns = 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 fR0 = −10−5 and different values of the total neutrino mass, namely mν = {0.1,  0.16} eV) simulated in a 8 h−3 Gpc3 volume containing (2 × )20483 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 C IDER simulations

The Constrained Interacting Dark EneRgy scenario (or C IDER, Barros 2019) is a particular type of coupled Quintessence models 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 C IDER 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 σ8 tension without further exacerbating the tension on H0. 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 C IDER 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 β = 0.03, 0.05, 0.08 besides a reference ΛCDM cosmology corresponding to the case β = 0. All simulations clearly share the same expansion history, consistent with the following cosmological parameters: Ωm = 1 − ΩΛ = 0.311, Ωb = 0.049, H0 = 67.7 km s−1Mpc−1, ns = 0.9665, As = 1.992 × 10−9, corresponding to a value of σ8 = 0.788 at z = 0 in the reference ΛCDM model. The simulations follow the evolution of 2 × 10243 particles for the (coupled) dark matter and (uncoupled) baryon components in a cosmological volume of 1 h−3 Gpc3 with a softening length of ε = 25 h−1 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 Geff = GN[1+(4/3)β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, σ, 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, σ, 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 σ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 Ωm = 1 − ΩΛ = 0.308, H0 = 67.8 km s−1 Mpc−1, ns = 0.966, As = 2.215 × 10−9, in a simulation box with a volume of 1 h−3 Gpc3 filled with 10243 dark matter particles and using a softening length of ε = 12 h−1 kpc.

The DAKAR2 simulations, instead, share the same cosmology and the same statistical realisation as the C IDER 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 C IDER 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(τ), appearing at the background level and kineticity, αK(τ), at the perturbation level. In the fluid picture of these theories, the relevant parameters are the speed of sound, cs(τ), and the equation-of-state parameter, w(τ), 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 w0 and cs2 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: (w0,  cs2) = (−0.9, 1 c2), ( − 0.9, 10−4 c2), ( − 0.9, 10−7 c2), and ( − 0.8, 10−7 c2). In these simulations, the box size was set to 2 h−1 Gpc with N = 12003 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 Gpc with N = Ngrid = 24003. In this series, the particle snapshots were saved in GADGET-2 format at five different redshifts z ∈ {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 Ωm (ΩΛ = 1 − Ωm), h, σ8, and f R 0 ¯ $ {\overline{f_{\mathrm{R0}}}} $ through 50 combinations (nodes) of these parameters sampled in a Latin-hypercube. All other cosmological parameters are fixed to a Planck cosmology (ns = 0.9652, Ωb = 0.049199, Planck Collaboration VI 2020). For each node, FORGE consists of a pair of large box simulations with 5123 particles in a 1.5 h−1 Gpc side-length box and a pair of high-resolution runs with 10243 particles in a 500 h−1 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 zinit = 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 R 0 ¯ $ {\overline{f_{\mathrm{R0}}}} $ parameter, the nDGP parameter, H0rc 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 fNL = 100. The simulation assumes the following ΛCDM parameters: Ωm = 1 − ΩΛ = 0.3089,H0 = 67.74 km s−1 Mpc−1, ns = 0.9667,σ8 = 0.8147. It consists of N = 40963 particles in L = 1 h−1 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 ε = 6 h−1 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 SLURM3 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.

thumbnail 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.

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 grid = 2 { floor [ log 2 ( N part 3 ) ] − 1 } $$ \begin{aligned} N_{\rm grid} = 2^{\left\{ \mathrm{floor} \left[\log _2\left(\root 3 \of {N_{\rm part}}\right)\right]-1\right\} } \end{aligned} $$(35)

linear grid size for every analysed snapshot, where Npart 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

⟨ δ ∼ ( k ) δ ∼ ∗ ( k ′ ) ⟩ = ( 2 π ) 3 P ( k ) δ D ( 3 ) ( k − k ′ ) , $$ \begin{aligned} \langle \tilde{\delta }(\boldsymbol{k})\tilde{\delta }^*(\boldsymbol{k}^{\prime }) \rangle = (2\pi )^3\,P(k)\,\delta _{\rm D}^{(3)}(\boldsymbol{k} - \boldsymbol{k}^{\prime }), \end{aligned} $$(36)

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

δ ( r ) = ρ ( r ) ρ ¯ − 1 , $$ \begin{aligned} \delta (\boldsymbol{r})\,=\,\frac{\rho (\boldsymbol{r})}{\overline{\rho }} - 1 , \end{aligned} $$(37)

and 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 | δ ∼ ( k ) | 2 $ |\tilde{\delta}({\boldsymbol 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

Δ k = k f = 2 π L box , $$ \begin{aligned} \Delta k = k_{\rm f} = \frac{2\pi }{L_{\rm box}}, \end{aligned} $$(38)

where kf is the fundamental wavenumber, and Lbox is the linear size of the simulation. The pipeline saves the power spectrum of every calculated bin below the

k Ny = π N grid L box $$ \begin{aligned} k_{\rm Ny} = \frac{\pi N_{\rm grid}}{L_{\rm box}} \end{aligned} $$(39)

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 Ps(k, μ) where μ = n ̂ LOS · k ̂ $ \mu = \hat{\boldsymbol{n}}_{\mathrm{LOS}} \cdot \hat{\boldsymbol{k}} $ and n ̂ LOS $ \hat{\boldsymbol{n}}_{\mathrm{LOS}} $ is a unit vector in the line-of-sight (LOS) direction. This can be expanded in multipoles, P s ( k , μ ) = ∑ ℓ = 0 ∞ P ℓ s ( k ) L ℓ ( μ ) $ P^s(k,\mu) = \sum\nolimits_{\ell=0}^\infty P_\ell^s(k)\mathcal{L}_\ell(\mu) $, where ℒℓ(μ) are the Legendre polynomials. The multipoles are then computed from the redshift-space power spectrum as

P ℓ s ( k ) = 2 ℓ + 1 2 ∫ − 1 1 P s ( k , μ ) L ℓ ( μ ) d μ . $$ \begin{aligned} P_\ell ^s(k) = \frac{2\ell +1}{2}\int _{-1}^1 P^s(k,\mu ) \mathcal{L} _\ell (\mu )\,\mathrm{d}\mu . \end{aligned} $$(40)

We computed the redshift-space power spectrum in 25 μ bins and the redshift-space multipoles (the monopole P0, the quadrupole P2, and the hexadecapole P4) using nbodykit from the input dark matter density field. For this, we used the distant-observer approximation

s i = r i + ( n ̂ LOS · v i aH ) · n ̂ LOS , $$ \begin{aligned} \boldsymbol{s}_{i} = \boldsymbol{r}_{i} + \left( \hat{\boldsymbol{n}}_{\rm LOS} \cdot \frac{\boldsymbol{v}_i}{a H} \right) \cdot \hat{\boldsymbol{n}}_{\rm LOS}, \end{aligned} $$(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, ri and vi are the real-space particle coordinates and peculiar velocities inside the periodic simulation box, si 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 Plin(k, zstart) input linear power spectrum of the simulation that was used during the initial condition generation. This can be defined at any zstart redshift. Then, this linear power spectrum is renormalised with the cosmological parameter σ8 at z = 0. The normalised Plin(k, z = 0) power spectrum is rescaled to zsnap redshifts of all analysed snapshots as

P lin ( k , z snap ) = [ D ( z snap ) D ( z = 0 ) ] 2 P lin ( k , z = 0 ) , $$ \begin{aligned} P_{\rm lin}(k,z_{\rm snap}) = \left[\frac{D(z_{\rm snap})}{D(z=0)}\right]^2 P_{\rm lin}(k,z=0) \;, \end{aligned} $$(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 ) = 5 Ω m H 0 2 2 H ( a ) ∫ 0 a d a ′ a ˙ ′ 3 $$ \begin{aligned} D(a) = \frac{5\Omega _\mathrm{m} H_0^2}{2}H(a)\int \limits _{0}^{a}\frac{\mathrm{d} a^{\prime }}{\dot{a}^{\prime 3}} \end{aligned} $$(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 ″ + [ 7 2 + 3 2 w ( a ) 1 + X ( a ) ] G ′ a + 3 2 1 − w ( a ) 1 + X ( a ) G a 2 = 0 $$ \begin{aligned} G^{\prime \prime } + \left[\frac{7}{2} + \frac{3}{2}\frac{w(a)}{1+X(a)} \right]\frac{G^{\prime }}{a} + \frac{3}{2} \frac{1-w(a)}{1+X(a)}\frac{G}{a^2} = 0 \end{aligned} $$(44)

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

X ( a ) = Ω m 1 − Ω m e − 3 ∫ a 1 d a ′ w ( a ′ ) a ′ . $$ \begin{aligned} X(a) = \frac{\Omega _\mathrm{m} }{1-\Omega _\mathrm{m} } \mathrm{e}^{-3\int \limits _{a}^{1}\mathrm{d} a^{\prime }\frac{w(a^{\prime })}{a^{\prime }}}. \end{aligned} $$(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 ≃ 0.28). It then does a new FoF search using the phase-space metric

d = | x 1 − x 2 | 2 σ x 2 + | v 1 − v 2 | 2 σ v 2 , $$ \begin{aligned} d = \sqrt{\frac{|\mathbf{x}_1 - \mathbf{x}_2|^2}{\sigma _x^2} + \frac{|\mathbf{v}_1 - \mathbf{v}_2|^2}{\sigma _v^2}} , \end{aligned} $$(46)

where σx and σ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 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: M200c (SO & BO), M500c (SO), M1000c (SO), M2500c (SO), M200b (SO). The Mvir 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 < log10[M/(h−1 M⊙)] < 14 with 32 logarithmic bins using the main mass definition (200c) and excluding 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 NumPy4 or Pandas5. To speed up the analysis pipeline, we therefore use Polars6, 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 = 〈 P h ( k ) P m ( k ) 〉 k < k max , $$ \begin{aligned} b = \left\langle \sqrt{\frac{P_{\rm h}(k)}{P_{\rm m}(k)}}\,\right\rangle _{k < k_{\rm max}}, \end{aligned} $$(47)

with Pm(k) and Ph(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 Ph(k) to the matter power spectrum Pm(k), and then averaging this ratio over all k bins where k < kmax 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 kmax = 0.1 h 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 < kmax 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 r500c. 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 5r500c. 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 http://skysrv.pha.jhu.edu/ neyrinck/voboz/ (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 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, Reff = [(3/4π)∑iVi]1/3, where Vi is the volume of the ith Voronoi cell belonging to the void. The void centre is defined as the volume-weighted barycentre, Xv = ∑ixiVi/Vtot, with Vtot = ∑iVi. 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, Xv, with radius Reff/4. VIDE also computes void shapes via the inertia tensor as well as the corresponding eigenvalues and eigenvectors. The ellipticity is then computed as ϵ = 1 − (J1/J3)1/4, where J1 and J3 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, Reff. 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 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 C IDER 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 cosmology7 corresponding to | f ¯ R 0 | = 10 − 5.34219 $ \left|\bar{f}_{\mathrm{R}0}\right| = 10^{-5.34219} $. In the case of the C IDER simulations, we focus on the IDE model with β = 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.

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 β = 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.

thumbnail 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.

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 ≃3%, and this difference increases to ≃5% for z = 0.55 compared to the ΛCDM model.

thumbnail Fig. 3.

Calculated matter and halo power spectra of the ELEPHANT simulation suite in the mass bin 1012.7h−1 M⊙ < Mhalo < 1013.2h−1 M⊙. 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.

thumbnail 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.

thumbnail Fig. 5.

Calculated matter and halo power spectra from the FORGE simulation suite in the mass bin 1012.7h−1 M⊙ < Mhalo < 1013.2h−1 M⊙. 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.

thumbnail Fig. 6.

Calculated halo mass function from the C IDER 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.

thumbnail Fig. 7.

Calculated matter and halo power spectra from the C IDER simulation suite in the mass bin 1012.7h−1 M⊙ < Mhalo < 1013.2h−1 M⊙. 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.

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 1012.7h−1 M⊙ < Mhalo < 1013.2h−1 M⊙ 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 C IDER 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 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 CosmoHub8 platform (Tallada et al. 2020; Carretero et al. 2017) designed for interactive exploration and distribution of massive cosmological datasets.

1

Note that we choose uc here for our IDE model to better distinguish it from ρCDM used right afterwards to describe the ΛCDM background evolution.

7

In particular we use node 10 of the Latin-Hypercube sampling described in Table 1 of Arnold et al. (2021).

Acknowledgments

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"orderungsgesellschaft 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 (http://www.euclid-ec.org). GR’s research was supported 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 partition. 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 collaboration with ICE-CSIC. CosmoHub received funding from the Spanish government (MCIN/AEI/10.13039/501100011033), the EU NextGeneration/PRTR (PRTRC17. I1), and the Generalitat de Catalunya. ZS acknowledges funding from DFG project 456622116 and support from the IRAP and IN2P3 Lyon computing centers. During part of this work, AMCLB was supported by a fellowship of PSL University-Paris Observatory. CG thanks the support from INAF theory Grant 2022: Illuminating Dark Matter using Weak Lensing by Cluster Satellites, PI: Carlo Giocoli. VGP is supported by the Atracción de Talento 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 Supercomputació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 significantly influenced the development of this project.

References

  1. Abbott, T. M. C., Aguena, M., Alarcon, A., et al. 2022, Phys. Rev. D, 105, 023520 [CrossRef] [Google Scholar]
  2. Adamek, J., Daverio, D., Durrer, R., & Kunz, M. 2016, JCAP, 7, 053 [CrossRef] [Google Scholar]
  3. Adamek, J., Angulo, R. E., Arnold, C., et al. 2023, JCAP, 06, 035 [NASA ADS] [Google Scholar]
  4. Alimi, J. M., Füzfa, A., Boucher, V., et al. 2010, MNRAS, 401, 775 [CrossRef] [Google Scholar]
  5. Amendola, L. 2000, Phys. Rev. D, 62, 043511 [NASA ADS] [CrossRef] [Google Scholar]
  6. Amendola, L. 2004, Phys. Rev. D, 69, 103524 [Google Scholar]
  7. Angulo, R. E., & Hahn, O. 2022, Liv. Rev. Comput. Astrophys., 8, 1 [CrossRef] [Google Scholar]
  8. Angulo, R. E., & Pontzen, A. 2016, MNRAS, 462, L1 [NASA ADS] [CrossRef] [Google Scholar]
  9. Appel, A. W. 1985, SIAM J. Sci. Stat. Comput., 6, 85 [NASA ADS] [CrossRef] [Google Scholar]
  10. Armendariz-Picon, C., Mukhanov, V., & Steinhardt, P. J. 2000, Phys. Rev. Lett., 85, 4438 [Google Scholar]
  11. Armendariz-Picon, C., Mukhanov, V., & Steinhardt, P. J. 2001, Phys. Rev. D, 63, 103510 [NASA ADS] [CrossRef] [Google Scholar]
  12. Arnold, C., Leo, M., & Li, B. 2019, Nat. Astron., 3, 945 [Google Scholar]
  13. Arnold, C., Li, B., Giblin, B., Harnois-Déraps, J., & Cai, Y. C. 2022, MNRAS, 515, 4161 [CrossRef] [Google Scholar]
  14. Avila, S., & Adame, A. G. 2023, MNRAS, 519, 3706 [NASA ADS] [CrossRef] [Google Scholar]
  15. Baldi, M. 2023, MNRAS, 521, 613 [NASA ADS] [Google Scholar]
  16. Baldi, M., & Simpson, F. 2015, MNRAS, 449, 2239 [NASA ADS] [CrossRef] [Google Scholar]
  17. Baldi, M., & Simpson, F. 2017, MNRAS, 465, 653 [CrossRef] [Google Scholar]
  18. Baldi, M., Pettorino, V., Robbers, G., & Springel, V. 2010, MNRAS, 403, 1684 [CrossRef] [Google Scholar]
  19. Baldi, M., Villaescusa-Navarro, F., Viel, M., et al. 2014, MNRAS, 440, 75 [NASA ADS] [CrossRef] [Google Scholar]
  20. Barreira, A., Li, B., Hellwing, W. A., et al. 2014, JCAP, 4, 029 [CrossRef] [Google Scholar]
  21. Barreira, A., Bose, S., & Li, B. 2015, JCAP, 12, 059 [CrossRef] [Google Scholar]
  22. Barreira, A., Cabass, G., Schmidt, F., Pillepich, A., & Nelson, D. 2020, JCAP, 12, 013 [NASA ADS] [CrossRef] [Google Scholar]
  23. Barros, B. J. 2019, Phys. Rev. D, 99, 064051 [NASA ADS] [Google Scholar]
  24. Barros, B. J., Castelão, D., da Fonseca, V., et al. 2023, JCAP, 1, 013 [NASA ADS] [Google Scholar]
  25. Bartelmann, M., & Schneider, P. 2001, Phys. Rep., 340, 291 [Google Scholar]
  26. Battat, J. B. R., Stubbs, C. W., & Chandler, J. F. 2008, Phys. Rev. D, 78, 022003 [NASA ADS] [Google Scholar]
  27. Becker, C., Arnold, C., Li, B., & Heisenberg, L. 2020, JCAP, 10, 055 [CrossRef] [Google Scholar]
  28. Behroozi, P., Wechsler, R., & Wu, H.-Y. 2012, Rockstar: Phase-space halo finder, Astrophysics Source Code Library [record ascl:1210.008] [Google Scholar]
  29. Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013, ApJ, 762, 109 [NASA ADS] [CrossRef] [Google Scholar]
  30. Berlind, A. A., Weinberg, D. H., Benson, A. J., et al. 2003, ApJ, 593, 1 [NASA ADS] [CrossRef] [Google Scholar]
  31. Biagetti, M., Lazeyras, T., Baldauf, T., Desjacques, V., & Schmidt, F. 2017, MNRAS, 468, 3277 [NASA ADS] [CrossRef] [Google Scholar]
  32. Bose, S., Li, B., Barreira, A., et al. 2017, JCAP, 2, 050 [CrossRef] [Google Scholar]
  33. Bose, B., Baldi, M., & Pourtsidou, A. 2018, JCAP, 4, 032 [CrossRef] [Google Scholar]
  34. Bowcock, P., Charmousis, C., & Gregory, R. 2000, CQG, 17, 4745 [NASA ADS] [Google Scholar]
  35. Boyle, A., Uhlemann, C., Friedrich, O., et al. 2021, MNRAS, 505, 2886 [NASA ADS] [CrossRef] [Google Scholar]
  36. Brax, P., Davis, A.-C., Li, B., Winther, H. A., & Zhao, G.-B. 2012, JCAP, 10, 002 [CrossRef] [Google Scholar]
  37. Brax, P., Davis, A.-C., Li, B., Winther, H. A., & Zhao, G.-B. 2013, JCAP, 4, 029 [CrossRef] [Google Scholar]
  38. Breton, M.-A., Rasera, Y., Taruya, A., Lacombe, O., & Saga, S. 2019, MNRAS, 483, 2671 [NASA ADS] [CrossRef] [Google Scholar]
  39. Buchdahl, H. A. 1970, MNRAS, 150, 1 [NASA ADS] [Google Scholar]
  40. Burrage, C., March, B., & Naik, A. P. 2024, JCAP, 04, 004 [NASA ADS] [Google Scholar]
  41. Caldwell, R. R., Kamionkowski, M., & Weinberg, N. N. 2003, Phys. Rev. Lett., 91, 071301 [NASA ADS] [CrossRef] [Google Scholar]
  42. Carbone, C., Petkova, M., & Dolag, K. 2016, JCAP, 7, 034 [CrossRef] [Google Scholar]
  43. Carretero, J., Tallada, P., Casals, J., et al. 2017, Proceedings of the European Physical Society Conference on High Energy Physics. 5-12 July, 488 [Google Scholar]
  44. Carrilho, P., Carrion, K., Bose, B., et al. 2022, MNRAS, 512, 3691 [Google Scholar]
  45. Casas, S., Cardone, V. F., Sapone, D., et al. 2023, A&A, submitted [arXiv:2306.11053] [Google Scholar]
  46. Castorina, E., Carbone, C., Bel, J., Sefusatti, E., & Dolag, K. 2015, JCAP, 7, 043 [CrossRef] [Google Scholar]
  47. Cautun, M., Paillas, E., Cai, Y.-C., et al. 2018, MNRAS, 476, 3195 [NASA ADS] [CrossRef] [Google Scholar]
  48. Chevallier, M., & Polarski, D. 2001, Int. J. Mod. Phys. D, 10, 213 [Google Scholar]
  49. Chuang, C.-H., Yepes, G., Kitaura, F.-S., et al. 2019, MNRAS, 487, 48 [NASA ADS] [CrossRef] [Google Scholar]
  50. Cielo, M., Escudero, M., Mangano, G., & Pisanti, O. 2023, Phys. Rev. D, 108, L121301 [CrossRef] [Google Scholar]
  51. Cole, S., Percival, W. J., Peacock, J. A., et al. 2005, MNRAS, 362, 505 [Google Scholar]
  52. Colless, M., Peterson, B. A., Jackson, C., et al. 2003, arXiv e-prints [arXiv:astro-ph/0306581] [Google Scholar]
  53. Contarini, S., Marulli, F., Moscardini, L., et al. 2021, MNRAS, 504, 5021 [NASA ADS] [CrossRef] [Google Scholar]
  54. Contarini, S., Verza, G., Pisani, A., et al. 2022, A&A, 667, A162 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Contarini, S., Pisani, A., Hamaus, N., et al. 2023, ApJ, 953, 46 [NASA ADS] [CrossRef] [Google Scholar]
  56. Creminelli, P., & Zaldarriaga, M. 2004, JCAP, 12, 006 [CrossRef] [Google Scholar]
  57. Crocce, M., Pueblas, S., & Scoccimarro, R. 2006, MNRAS, 373, 369 [NASA ADS] [CrossRef] [Google Scholar]
  58. Cuozzo, V., Carbone, C., Calabrese, M., Carella, E., & Migliaccio, M. 2024, JCAP, 4, 073 [NASA ADS] [Google Scholar]
  59. Dalal, N., Doré, O., Huterer, D., & Shirokov, A. 2008, Phys. Rev. D, 77, 123514 [CrossRef] [Google Scholar]
  60. Das, S., Corasaniti, P. S., & Khoury, J. 2006, Phys. Rev. D, 73, 083509 [NASA ADS] [CrossRef] [Google Scholar]
  61. de Bernardis, P., Ade, P. A. R., Bock, J. J., et al. 2000, Nature, 404, 955 [Google Scholar]
  62. Deffayet, C. 2001, Phys. Lett. B, 502, 199 [NASA ADS] [CrossRef] [Google Scholar]
  63. DESI Collaboration (Aghamousa, A., et al.) 2016, arXiv e-prints [arXiv:1611.00036] [Google Scholar]
  64. Desjacques, V., Seljak, U., & Iliev, I. T. 2009, MNRAS, 396, 85 [NASA ADS] [CrossRef] [Google Scholar]
  65. Di Valentino, E., Mena, O., Pan, S., et al. 2021, CQG, 38, 153001 [NASA ADS] [CrossRef] [Google Scholar]
  66. Diemer, B. 2018, ApJS, 239, 35 [NASA ADS] [CrossRef] [Google Scholar]
  67. Dolag, K., Bartelmann, M., Perrotta, F., et al. 2004, A&A, 416, 853 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  68. Doré, O., Bock, J., Ashby, M., et al. 2014, arXiv e-prints [arXiv:1412.4872] [Google Scholar]
  69. Dvali, G., Gabadadze, G., & Porrati, M. 2000, Phys. Lett. B, 485, 208 [Google Scholar]
  70. Efstathiou, G., Moody, S., Peacock, J. A., et al. 2002, MNRAS, 330, L29 [NASA ADS] [CrossRef] [Google Scholar]
  71. Eggemeier, A., Camacho-Quevedo, B., Pezzotta, A., et al. 2023, MNRAS, 519, 2962 [Google Scholar]
  72. Eisenstein, D. J., Zehavi, I., Hogg, D. W., et al. 2005, ApJ, 633, 560 [Google Scholar]
  73. Euclid Collaboration (Blanchard, A., et al.) 2020, A&A, 642, A191 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  74. Euclid Collaboration (Desprez, G., et al.) 2020, A&A, 644, A31 [EDP Sciences] [Google Scholar]
  75. Euclid Collaboration (Ilbert, O., et al.) 2021, A&A, 647, A117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Euclid Collaboration (Knabenhans, M., et al.) 2021, MNRAS, 505, 2840 [NASA ADS] [CrossRef] [Google Scholar]
  77. Euclid Collaboration (Bretonnière, H., et al.) 2022, A&A, 657, A90 [CrossRef] [EDP Sciences] [Google Scholar]
  78. Euclid Collaboration (Scaramella, R., et al.) 2022, A&A, 662, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  79. Euclid Collaboration (Ajani, V., et al.) 2023, A&A, 675, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  80. Euclid Collaboration (Bretonnière, H., et al.) 2023, A&A, 671, A102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  81. Euclid Collaboration (Castro, T., et al.) 2023, A&A, 671, A100 [CrossRef] [EDP Sciences] [Google Scholar]
  82. Euclid Collaboration (Merlin, E., et al.) 2023d, A&A, 671, A101 [Google Scholar]
  83. Euclid Collaboration (Koyama, K., et al.) 2024, A&A, submitted [arXiv:2409.03524] [Google Scholar]
  84. Euclid Collaboration (Adamek, J., et al.) 2025, A&A, 695, A230 [Google Scholar]
  85. Euclid Collaboration (Mellier, Y., et al.) 2025, A&A, in press, https://doi.org/10.1051/0004-6361/202450810 [Google Scholar]
  86. Farrar, G. R., & Peebles, P. J. E. 2004, ApJ, 604, 1 [Google Scholar]
  87. Feng, Y., Chu, M.-Y., Seljak, U., & McDonald, P. 2016, MNRAS, 463, 2273 [NASA ADS] [CrossRef] [Google Scholar]
  88. Feng, Y., Bird, S., & Lanusse, F. 2017, https://doi.org/10.5281/zenodo.1051252 [Google Scholar]
  89. Fiorini, B., Koyama, K., Izard, A., et al. 2021, JCAP, 9, 021 [NASA ADS] [Google Scholar]
  90. Fiorini, B., Koyama, K., & Baker, T. 2023, JCAP, 12, 045 [NASA ADS] [Google Scholar]
  91. García-Farieta, J. E., Marulli, F., Veropalumbo, A., et al. 2019, MNRAS, 488, 1987 [Google Scholar]
  92. Giannantonio, T., Scranton, R., Crittenden, R. G., et al. 2008, Phys. Rev. D, 77, 123520 [Google Scholar]
  93. Giocoli, C., Baldi, M., & Moscardini, L. 2018, MNRAS, 481, 2813 [Google Scholar]
  94. Gouyou Beauchamps, S., Baratta, P., Escoffier, S., et al. 2024, A&A, 693, A226 [Google Scholar]
  95. Grieb, J. N., Sánchez, A. G., Salazar-Albornoz, S., & Dalla Vecchia, C. 2016, MNRAS, 457, 1577 [NASA ADS] [CrossRef] [Google Scholar]
  96. Guillet, T., & Teyssier, R. 2011, J. Comput. Phys., 230, 4756 [Google Scholar]
  97. Gutiérrez Adame, A., Avila, S., Gonzalez-Perez, V., et al. 2024, A&A, 689, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  98. Hagstotz, S., Costanzi, M., Baldi, M., & Weller, J. 2019a, MNRAS, 486, 3927 [Google Scholar]
  99. Hagstotz, S., Gronke, M., Mota, D. F., & Baldi, M. 2019b, A&A, 629, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  100. Hamaus, N., Seljak, U., & Desjacques, V. 2011, Phys. Rev. D, 84, 083509 [NASA ADS] [CrossRef] [Google Scholar]
  101. Hamaus, N., Aubert, M., Pisani, A., et al. 2022, A&A, 658, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  102. Hand, N., Feng, Y., Beutler, F., et al. 2018, AJ, 156, 160 [Google Scholar]
  103. Harnois-Déraps, J., Hernandez-Aguayo, C., Cuesta-Lazaro, C., et al. 2023, MNRAS, 525, 6336 [Google Scholar]
  104. Hassani, F., Adamek, J., Kunz, M., & Vernizzi, F. 2019, JCAP, 12, 011 [CrossRef] [Google Scholar]
  105. Hassani, F., L’Huillier, B., Shafieloo, A., Kunz, M., & Adamek, J. 2020, JCAP, 4, 039 [CrossRef] [Google Scholar]
  106. Hernández-Aguayo, C., Arnold, C., Li, B., & Baugh, C. M. 2021, MNRAS, 503, 3867 [Google Scholar]
  107. Hopkins, P. F. 2015, MNRAS, 450, 53 [Google Scholar]
  108. Hu, W., & Sawicki, I. 2007, Phys. Rev. D, 76, 064004 [CrossRef] [Google Scholar]
  109. Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111 [Google Scholar]
  110. Jaroszynski, M., Park, C., Paczynski, B., Gott, J., & Richard, I. 1990, ApJ, 365, 22 [Google Scholar]
  111. Joseph, M., Aloni, D., Schmaltz, M., Sivarajan, E. N., & Weiner, N. 2023, Phys. Rev. D, 108, 023520 [NASA ADS] [CrossRef] [Google Scholar]
  112. Kajita, T. 2016, Rev. Mod. Phys., 88, 030501 [Google Scholar]
  113. Klypin, A. A., & Shandarin, S. F. 1983, MNRAS, 204, 891 [NASA ADS] [Google Scholar]
  114. Klypin, A., Macciò, A. V., Mainini, R., & Bonometto, S. A. 2003, ApJ, 599, 31 [Google Scholar]
  115. Komatsu, E., & Spergel, D. N. 2001, Phys. Rev. D, 63, 063002 [NASA ADS] [CrossRef] [Google Scholar]
  116. Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18 [Google Scholar]
  117. Kovac, J. M., Leitch, E. M., Pryke, C., et al. 2002, Nature, 420, 772 [NASA ADS] [CrossRef] [Google Scholar]
  118. Koyama, K. 2016, Rept. Prog. Phys., 79, 046902 [NASA ADS] [CrossRef] [Google Scholar]
  119. Kreisch, C. D., Pisani, A., Carbone, C., et al. 2019, MNRAS, 488, 4413 [NASA ADS] [CrossRef] [Google Scholar]
  120. Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv e-prints [arXiv:1110.3193] [Google Scholar]
  121. Le Brun, A. M. C., Arnaud, M., Pratt, G. W., & Teyssier, R. 2018, MNRAS, 473, L69 [Google Scholar]
  122. Lesgourgues, J., & Pastor, S. 2006, Phys. Rep., 429, 307 [Google Scholar]
  123. Lewis, A., & Challinor, A. 2011, CAMB: Code for Anisotropies in the Microwave Background, Astrophysics Source Code Library [record ascl:1102.026] [Google Scholar]
  124. Li, B., Zhao, G.-B., Teyssier, R., & Koyama, K. 2012, JCAP, 1, 051 [Google Scholar]
  125. Li, B., Barreira, A., Baugh, C. M., et al. 2013a, JCAP, 11, 012 [CrossRef] [Google Scholar]
  126. Li, B., Zhao, G.-B., & Koyama, K. 2013b, JCAP, 5, 023 [CrossRef] [Google Scholar]
  127. Li, S.-S., Hoekstra, H., Kuijken, K., et al. 2023, A&A, 679, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  128. Linden, S., & Virey, J.-M. 2008, Phys. Rev. D, 78, 023526 [NASA ADS] [Google Scholar]
  129. Linder, E. V. 2003, Phys. Rev. Lett., 90, 091301 [Google Scholar]
  130. Linder, E. V. 2008, Gen. Rel. Grav., 40, 329 [NASA ADS] [CrossRef] [Google Scholar]
  131. Linder, E. V., & Jenkins, A. 2003, MNRAS, 346, 573 [Google Scholar]
  132. Maldacena, J. 2003, J. High Energy Phys., 2003, 013 [CrossRef] [Google Scholar]
  133. Maltoni, M., Schwetz, T., Tórtola, M., & Valle, J. W. F. 2004, New J. Phys., 6, 122 [Google Scholar]
  134. Mangano, G., Miele, G., Pastor, S., & Peloso, M. 2002, Phys. Lett. B, 534, 8 [CrossRef] [Google Scholar]
  135. Martinelli, M., & Tutusaus, I. 2019, Symmetry, 11, 986 [Google Scholar]
  136. Martinelli, M., Martins, C. J. A. P., Nesseris, S., et al. 2021, A&A, 654, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  137. Massara, E., Villaescusa-Navarro, F., Viel, M., & Sutter, P. M. 2015, JCAP, 11, 018 [Google Scholar]
  138. Matarrese, S., & Verde, L. 2008, ApJ, 677, L77 [NASA ADS] [CrossRef] [Google Scholar]
  139. Merson, A. I., Baugh, C. M., Helly, J. C., et al. 2013, MNRAS, 429, 556 [NASA ADS] [CrossRef] [Google Scholar]
  140. Merten, J., Giocoli, C., Baldi, M., et al. 2019, MNRAS, 487, 104 [Google Scholar]
  141. Michaux, M., Hahn, O., Rampf, C., & Angulo, R. E. 2021, MNRAS, 500, 663 [Google Scholar]
  142. Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563 [Google Scholar]
  143. Nesseris, S., Sapone, D., Martinelli, M., et al. 2022, A&A, 660, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  144. Neyrinck, M. C. 2008, MNRAS, 386, 2101 [CrossRef] [Google Scholar]
  145. Parimbelli, G., Carbone, C., Bel, J., et al. 2022, JCAP, 11, 041 [Google Scholar]
  146. Peebles, P. J. E. 1993, Principles of Physical Cosmology (Princeton University Press) [Google Scholar]
  147. Peel, A., Pettorino, V., Giocoli, C., Starck, J.-L., & Baldi, M. 2018, A&A, 619, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  148. Peel, A., Lalande, F., Starck, J.-L., et al. 2019, Phys. Rev. D, 100, 023508 [Google Scholar]
  149. Perlmutter, S., Aldering, G., Goldhaber, G., et al. 1999, ApJ, 517, 565 [Google Scholar]
  150. Piga, L., Marinucci, M., D’Amico, G., et al. 2023, JCAP, 04, 038 [Google Scholar]
  151. Pisani, A., Sutter, P. M., Hamaus, N., et al. 2015, Phys. Rev. D, 92, 083531 [NASA ADS] [CrossRef] [Google Scholar]
  152. Planck Collaboration XVI. 2014, A&A, 571, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  153. Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  154. Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  155. Platen, E., van de Weygaert, R., & Jones, B. J. T. 2007, MNRAS, 380, 551 [NASA ADS] [CrossRef] [Google Scholar]
  156. Potter, D., Stadel, J., & Teyssier, R. 2017, Comput. Astrophys. Cosmol., 4, 2 [NASA ADS] [CrossRef] [Google Scholar]
  157. Pourtsidou, A., & Tram, T. 2016, Phys. Rev. D, 94, 043518 [CrossRef] [Google Scholar]
  158. Pourtsidou, A., Skordis, C., & Copeland, E. J. 2013, arXiv e-prints [arXiv:1307.0458] [Google Scholar]
  159. Press, W. H., & Schechter, P. 1974, ApJ, 187, 425 [Google Scholar]
  160. Prunet, S., Pichon, C., Aubert, D., et al. 2008, ApJS, 178, 179 [NASA ADS] [CrossRef] [Google Scholar]
  161. Puchwein, E., Baldi, M., & Springel, V. 2013, MNRAS, 436, 348 [NASA ADS] [CrossRef] [Google Scholar]
  162. Rácz, G., Kiessling, A., Csabai, I., & Szapudi, I. 2023, A&A, 672, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  163. Rasera, Y., Breton, M. A., Corasaniti, P. S., et al. 2022, A&A, 661, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  164. Riess, A. G., Filippenko, A. V., Challis, P., et al. 1998, AJ, 116, 1009 [Google Scholar]
  165. Riess, A. G., Yuan, W., Macri, L. M., et al. 2022, ApJ, 934, L7 [NASA ADS] [CrossRef] [Google Scholar]
  166. Roncarelli, M., Carbone, C., & Moscardini, L. 2015, MNRAS, 447, 1761 [NASA ADS] [CrossRef] [Google Scholar]
  167. Roy, F., Bouillot, V. R., & Rasera, Y. 2014, A&A, 564, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  168. Ruan, C.-Z., Cuesta-Lazaro, C., Eggemeier, A., et al. 2024, MNRAS, 527, 2490 [Google Scholar]
  169. Sakr, Z. 2022, Universe, 8, 284 [Google Scholar]
  170. Salopek, D. S., & Bond, J. R. 1990, Phys. Rev. D, 42, 3936 [NASA ADS] [CrossRef] [Google Scholar]
  171. Schaap, W. E., & van de Weygaert, R. 2000, A&A, 363, L29 [Google Scholar]
  172. Schmidt, F. 2009, Phys. Rev. D, 80, 123003 [NASA ADS] [CrossRef] [Google Scholar]
  173. Sefusatti, E., Crocce, M., Scoccimarro, R., & Couchman, H. M. P. 2016, MNRAS, 460, 3624 [NASA ADS] [CrossRef] [Google Scholar]
  174. Simpson, F. 2010, Phys. Rev. D, 82, 083505 [NASA ADS] [CrossRef] [Google Scholar]
  175. Skordis, C., Pourtsidou, A., & Copeland, E. J. 2015, Phys. Rev. D, 91, 083537 [NASA ADS] [Google Scholar]
  176. Slosar, A., Hirata, C., Seljak, U., Ho, S., & Padmanabhan, N. 2008, JCAP, 8, 031 [CrossRef] [Google Scholar]
  177. Smith, A., de Mattia, A., Burtin, E., Chuang, C.-H., & Zhao, C. 2021, MNRAS, 500, 259 [Google Scholar]
  178. Spergel, D. N., Verde, L., Peiris, H. V., et al. 2003, ApJS, 148, 175 [Google Scholar]
  179. Spergel, D., Gehrels, N., Baltay, C., et al. 2015, arXiv e-prints [arXiv:1503.03757] [Google Scholar]
  180. Springel, V. 2005, MNRAS, 364, 1105 [Google Scholar]
  181. Springel, V. 2010, MNRAS, 401, 791 [Google Scholar]
  182. Springel, V. 2015, N-GenIC: Cosmological structure initial conditions, Astrophysics Source Code Library [record ascl:1502.003] [Google Scholar]
  183. Sutter, P. M., Lavaux, G., Hamaus, N., et al. 2015, Astron. Comput., 9, 1 [NASA ADS] [CrossRef] [Google Scholar]
  184. Takada, M., Komatsu, E., & Futamase, T. 2006, Phys. Rev. D, 73, 083520 [NASA ADS] [Google Scholar]
  185. Tallada, P., Carretero, J., Casals, J., et al. 2020, Astron. Comput., 32, 100391 [Google Scholar]
  186. Taruya, A., Nishimichi, T., & Saito, S. 2010, Phys. Rev. D, 82, 063522 [NASA ADS] [CrossRef] [Google Scholar]
  187. Tassev, S., Zaldarriaga, M., & Eisenstein, D. J. 2013, JCAP, 6, 036 [CrossRef] [Google Scholar]
  188. Tegmark, M., Blanton, M. R., Strauss, M. A., et al. 2004, ApJ, 606, 702 [Google Scholar]
  189. Tegmark, M., Eisenstein, D. J., Strauss, M. A., et al. 2006, Phys. Rev. D, 74, 123507 [Google Scholar]
  190. Teyssier, R. 2002, A&A, 385, 337 [CrossRef] [EDP Sciences] [Google Scholar]
  191. Verza, G., Pisani, A., Carbone, C., Hamaus, N., & Guzzo, L. 2019, JCAP, 12, 040 [CrossRef] [Google Scholar]
  192. Verza, G., Carbone, C., & Renzi, A. 2022, ApJ, 940, L16 [NASA ADS] [CrossRef] [Google Scholar]
  193. Verza, G., Carbone, C., Pisani, A., & Renzi, A. 2023, JCAP, 12, 044 [Google Scholar]
  194. Verza, G., Carbone, C., Pisani, A., Porciani, C., & Matarrese, S. 2024, JCAP, 10, 079 [NASA ADS] [Google Scholar]
  195. Viel, M., Haehnelt, M. G., & Springel, V. 2010, JCAP, 6, 015 [CrossRef] [Google Scholar]
  196. Weinberger, R., Springel, V., & Pakmor, R. 2020, ApJS, 248, 32 [Google Scholar]
  197. Winther, H. A., Schmidt, F., Barreira, A., et al. 2015, MNRAS, 454, 4208 [NASA ADS] [CrossRef] [Google Scholar]
  198. Winther, H. A., Koyama, K., Manera, M., Wright, B. S., & Zhao, G.-B. 2017, JCAP, 8, 006 [CrossRef] [Google Scholar]
  199. Zeldovich, I. B. 1978, in Large Scale Structures in the Universe, eds. M. S. Longair, & J. Einasto, 79, 409 [Google Scholar]
  200. Zennaro, M., Bel, J., Villaescusa-Navarro, F., et al. 2017, MNRAS, 466, 3244 [NASA ADS] [CrossRef] [Google Scholar]
  201. Zennaro, M., Bel, J., Dossett, J., Carbone, C., & Guzzo, L. 2018, MNRAS, 477, 491 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1.

Overview of the simulation suites analysed for this project.

In the text

All Figures

thumbnail 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.
In the text
thumbnail 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.
In the text
thumbnail Fig. 3.

Calculated matter and halo power spectra of the ELEPHANT simulation suite in the mass bin 1012.7h−1 M⊙ < Mhalo < 1013.2h−1 M⊙. 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.
In the text
thumbnail 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.
In the text
thumbnail Fig. 5.

Calculated matter and halo power spectra from the FORGE simulation suite in the mass bin 1012.7h−1 M⊙ < Mhalo < 1013.2h−1 M⊙. 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.
In the text
thumbnail Fig. 6.

Calculated halo mass function from the C IDER 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.
In the text
thumbnail Fig. 7.

Calculated matter and halo power spectra from the C IDER simulation suite in the mass bin 1012.7h−1 M⊙ < Mhalo < 1013.2h−1 M⊙. 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.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.