© The Authors 2025

1. Introduction

Structure formation in the Universe proceeds hierarchically. Small-scale perturbations collapse and merge to form progressively larger structures over cosmic time. Galaxy clusters, being the most massive virialized objects, occupy a prominent position in this hierarchy and serve as powerful probes of cosmology (see reviews by Allen et al. 2011; Kravtsov & Borgani 2012).

The halo mass function (HMF), which provides the comoving number density of matter halos as a function of mass and redshift, is a fundamental theoretical tool for interpreting observations of galaxy clusters. Due to the strong nonlinearity behind galaxy cluster formation and evolution, analytical modeling of the HMF is limited to a qualitative description that quantitatively does not fulfill the requirements needed for precision cosmology (e.g., Press & Schechter 1974; Bond et al. 1991). To circumvent this, the community started developing physically motivated semi-analytical models (e.g., Sheth & Tormen 1999; Sheth et al. 2001; Tinker et al. 2008; Despali et al. 2016; Ondaro-Mallea et al. 2021; Euclid Collaboration: Castro et al. 2023) that, when calibrated on simulations, provide a much more accurate prediction of the halo abundances.

The most comprehensive simulations to calibrate HMF models are cosmological hydrodynamic simulations, as they simulate gravity and baryonic effects associated with galaxy formation. Baryons, although subdominant to dark matter, are known to impact the cluster abundance in a significant way (e.g., Cui et al. 2014; Velliscig et al. 2014; Bocquet et al. 2016; Castro et al. 2018, 2020; Schaye et al. 2023; Euclid Collaboration: Castro et al. 2024). However, their cost limits the number of simulations that can be used for a robust exploratory model calibration. Therefore, the community took the strategy of calibrating the HMF model on gravity-only simulations and adding the baryonic impact as post-processing. We followed this strategy too and focus on gravity-only simulations in the rest of the paper.

Recent and forthcoming observational campaigns, such as those conducted by the Vera C. Rubin Observatory’s Legacy Survey of Space and Time (LSST; LSST Science Collaboration 2009)¹, the third generation of the South Pole Telescope (SPT-3G; Benson 2014)², the extended ROentgen Survey with an Imaging Telescope Array (eROSITA; Predehl et al. 2021)³, the Square Kilometre Array (SKA; Maartens et al. 2015)⁴, the Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration 2016)⁵, the Nancy Grace Roman Space Telescope (Spergel et al. 2015)⁶, and Euclid (Euclid Collaboration: Mellier et al. 2025)⁷, will provide measurements with unprecedented statistical precision. To fully exploit these data, theoretical models must keep pace and achieve sub-percent accuracy so that systematic uncertainties remain subdominant.

In the standard cosmological model, the Universe’s accelerated expansion is driven by a cosmological constant, Λ, which dominates the current energy density composition of the Universe, followed by cold dark matter (CDM) and baryons. Altogether, these components correspond to more than 99 percent of the current Universe’s energetic composition, with massive neutrinos and relativistic species making up the remaining 1 percent. Euclid Collaboration: Castro et al. (2023) present a calibration of the HMF that reproduces with sub-percent accuracy the outcome of simulations assuming a standard ΛCDM cosmological model. A straightforward extension of the cosmological constant (see, e.g., Peebles & Ratra 2003; Copeland et al. 2006, for reviews) is to consider a smooth dark energy (DE) component with a variable equation of state (EOS) and a constant sound speed, c_s = 1. Within this phenomenological approach, DE remains effectively homogeneous on sub-horizon scales, and its perturbations do not significantly affect structure formation. Castro et al. (2025, hereafter Paper I) extended the HMF model of Euclid Collaboration: Castro et al. (2023) to dynamical DE models described by the Chevallier–Polarski–Linder (CPL) parametrization (Chevallier & Polarski 2001; Linder 2003).

Clustering dark energy (CDE) becomes significant on sub-horizon scales when the DE sound speed (c_s) drops below unity, permitting its perturbations to grow rather than remain suppressed. This regime was first quantified by Weller & Lewis (2003) and Bean & Dore (2004), who analyzed its imprints on the cosmic microwave background and the matter power spectrum. Explicit field-theory realizations – such as k-essence and tachyon models – were formulated by Garriga & Mukhanov (1999), Armendariz-Picon et al. (2000), and Sen (2002), who demonstrated that c_s can deviate substantially from unity. To assess nonlinear effects, Mota & van de Bruck (2004) applied a spherical collapse model to inhomogeneous DE, finding consequential shifts in the collapse threshold. These analytic approaches were refined by Abramo et al. (2009) and Chang et al. (2018). More recently, cosmological simulation codes, such as Dakin et al. (2019a), Hassani et al. (2019), and Blot et al. (2023), were developed. In Dakin et al. (2019a), a linear-realization ansatz was adopted, capturing the leading impact of DE perturbations on dark matter dynamics while omitting the backreaction of nonlinear matter clustering on the DE field. In Hassani et al. (2019), the authors employed an effective field theory (EFT; see Gubitosi et al. 2013) framework to evolve DE perturbations and dark matter in a coupled nonlinear scheme. Lastly, in Blot et al. (2023), DE was modeled as an effective fluid, and the associated continuity and Euler equations were solved using a hydrodynamic scheme.

The nonlinear modeling of CDE remains an open research frontier (Creminelli et al. 2009; Frusciante & Papadomanolakis 2017; Cusin et al. 2018a,b). For instance, Hassani et al. (2022, 2023) and Eckmann et al. (2023) report potential instabilities in the EFT description. Due to the unresolved issues with fully nonlinear schemes, we adopted the linear realization approach proposed by Dakin et al. (2019a) hereafter and refer to Batista (2021) for a deeper discussion on CDE models and their effects on structure formation.

Batista & Marra (2017) show that a null sound speed introduces an additional degree of freedom that alters the growth of density perturbations. In scenarios where c_s < 1, DE perturbations are less suppressed on small scales, modifying the gravitational instability that drives the collapse of dark matter halos. This change affects key quantities in the spherical collapse model – such as the critical density threshold (e.g., Press & Schechter 1974; Bond et al. 1991) and the virial overdensity (e.g., Sheth & Tormen 1999; Sheth et al. 2001) – and ultimately alters the shape of the HMF. Therefore, an accurate calibration of the HMF in CDE cosmologies must explicitly account for the impact of a variable sound speed on the nonlinear evolution of density perturbations.

The impending arrival of high-precision data underscores the necessity for improved theoretical modeling. As noted by Salvati et al. (2020) and Artis et al. (2021), ongoing and forthcoming surveys of galaxy clusters will probe the large-scale structure of the Universe with an accuracy that demands theoretical predictions be controlled at the sub-percent level. Extending the HMF’s calibration and incorporating the effects of DE clustering with a variable sound speed is imperative to meet these requirements.

In this work we present a new HMF model that explicitly includes the impact of a variable sound speed of the DE fluid. We achieve this by extending the Dynamic Universe Cosmological Analysis (DUCA) N-body suite of simulations presented in Paper I to include models with c_s < 1. Remarkably, our new HMF model reproduces the results from the new CDE simulations without reducing the accuracy achieved before.

This paper is organized as follows: Section 2 describes the theoretical framework for CDE with a variable sound speed. Section 3 details the simulation setup employed to model the extended HMF model. Section 4 presents the main results, including the modeling, comparisons with previous HMF calibrations, and tests of model robustness. We discuss the implications of our results in Sect. 5. Finally, in Sect. 6 we draw our conclusions.

2. Theoretical framework

In this section we present a short overview of the main concepts of the CDE model adopted. We present the background evolution in Sect. 2.1, the perturbations in Sect. 2.2, and the HMF formalism in Sect. 2.3.

2.1. Background evolution

We adopted the CPL phenomenological framework to characterize the evolution of the DE EOS, w_DE, which links the DE pressure (p_DE) and energy density (ρ_DE):

$$\begin{array}{c}\hfill w_{\mathrm{DE}}(a)=\frac{p_{\mathrm{DE}}}{{\rho }_{\mathrm{DE}}}.\end{array}$$(1)

In this approach, w_DE is modeled as a linear function of the scale factor (a):

$$\begin{array}{c}\hfill w_{\mathrm{DE}}(a)=w_0+w_a(1-a),\end{array}$$(2)

where w₀ represents the present-day value of the EOS, and w_a quantifies its rate of change with the scale factor.

In a spatially flat universe that contains matter, radiation, massive neutrinos, and DE, the Hubble parameter H(z) is given by

$$\begin{array}{cc}\hfill \frac{H^2(z)}{H_0^2}& ={\mathrm{\Omega }}_{\mathrm{m},0}{(1+z)}^3+{\mathrm{\Omega }}_{\mathrm{r},0}{(1+z)}^4+{\mathrm{\Omega }}_{\nu }(z)\frac{{\rho }_{\mathrm{c}}(z)}{{\rho }_{\mathrm{c},0}}\hfill \\ \hfill & +{\mathrm{\Omega }}_{\mathrm{DE},0}{(1+z)}^{3(1+w_0+w_a)}exp[-\frac{3w_{a}z}{1+z}],\hfill \end{array}$$(3)

where H₀ denotes the present-day Hubble constant; Ω_m,0, Ω_r,0, and Ω_DE,0 are the current density parameters for matter (both baryonic and CDM), radiation, and DE, respectively. Here, Ω_ν(z) is the neutrino density parameter at redshift z, while ρ_c(z) and ρ_c,0 are the critical densities at redshift z and at present, respectively. While massive neutrinos also significantly impact cosmic evolution by influencing both the expansion history and the growth of structure (Lesgourgues & Pastor 2006), we do not consider them in the rest of the paper as their impact on the HMF has already been studied (Castorina et al. 2014; Costanzi et al. 2013; Vagnozzi et al. 2018) and validated (Euclid Collaboration: Castro et al. 2023, Paper I).

A key aspect of DE models is their ability to account for the observed accelerated expansion of the Universe. The acceleration is quantified by the deceleration parameter, q_dec, defined as (Weinberg 2008; Peebles 2020)

$$\begin{array}{c}\hfill q_{\mathrm{dec}}=-\frac{\ddot{a}}{a{H}^2},\end{array}$$(4)

with the overdots indicating derivatives with respect to cosmic time. By invoking the Friedmann equations, q_dec can be expressed in terms of the density parameters and their associated EOSs. At the current epoch (z = 0), the deceleration parameter reduces approximately to

$$\begin{array}{c}\hfill q_{\mathrm{dec},0}\approx \frac{1}{2}[2{\mathrm{\Omega }}_{\mathrm{m},0}+{\mathrm{\Omega }}_{\mathrm{DE},0}(1+3w_0)].\end{array}$$(5)

A negative q_dec signifies that the Universe is accelerating, which requires that w₀ < −1/3.

2.2. Dark energy perturbations

As the nature of the DE component that drives the accelerated expansion remains unknown, a compelling description of DE is achieved by parameterizing its perturbations. Two widely used approaches that have been implemented in both CLASS (Blas et al. 2011) and CAMB (Lewis et al. 2000) are the fluid formalism and the parametrized post-Friedmann (PPF) formalism (Hu 2008; Fang et al. 2008; Baker et al. 2013). The fluid formalism evolves DE density and velocity perturbations directly from an assumed EOS and sound speed, offering transparent physical interpretation. The PPF framework, by contrast, enforces a stable interpolation across regimes where the DE EOS crosses the phantom divide (w = −1), thus avoiding the divergences that plague simple fluid treatments in that case. In our simulations, we adopted one of these two prescriptions, both of which are available within the CONncept code and encompass models such as quintessence-like scenarios. In what follows, we briefly review both formalisms.

2.2.1. Dark energy as a fluid

In the fluid formalism, DE is described by an EOS, w(a), and a constant rest-frame sound speed, c_s. The EOS adopted in this paper is described in Sect. 2.1; however, the discussion below is independent of the particular form of w(a).

In the synchronous gauge, the continuity and Euler equations for the DE fluid are written in Fourier space as (see, e.g., Ballesteros & Lesgourgues 2010)

$$\begin{array}{cc}\hfill {\delta }_{\mathrm{DE}}^{\prime \mathrm{s}}& =-(1+w)({\theta }_{\mathrm{DE}}^{\mathrm{s}}+\frac{{\mathfrak{H}}^{\prime }}{2})-3(c_{\mathrm{s}}^2-w)\mathcal{H}{\delta }_{\mathrm{DE}}^{\mathrm{s}}\hfill \\ \hfill & -9(1+w)(c_{\mathrm{s}}^2-c_{\mathrm{a}}^2){\mathcal{H}}^2\frac{{\theta }_{\mathrm{DE}}^{\mathrm{s}}}{k^2},\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill {\theta }_{\mathrm{DE}}^{\prime \mathrm{s}}& =-(1-3c_{\mathrm{s}}^2)\mathcal{H}{\theta }_{\mathrm{DE}}^{\mathrm{s}}+\frac{c_{\mathrm{s}}^2{k}^2}{1+w}{\delta }_{\mathrm{DE}}^{\mathrm{s}}-k^2{\sigma }_{\mathrm{DE}},\hfill \end{array}$$(7)

where a prime superscript denotes a derivative with respect to conformal time τ, ${\delta }_{\mathrm{DE}}^{\mathrm{s}}$ and ${\theta }_{\mathrm{DE}}^{\mathrm{s}}$ are the DE density contrast and velocity divergence, respectively, and ℋ ≡ a′/a, ℌ is one of the synchronous gauge metric perturbations⁸. The standard assumption is vanishing anisotropic stress, i.e., σ_DE = 0. Lastly, the adiabatic sound speed is defined as

$$\begin{array}{c}\hfill c_{\mathrm{a}}^2\equiv \frac{p_{\mathrm{DE}}^{\prime }}{{\rho }_{\mathrm{DE}}^{\prime }}.\end{array}$$(8)

This fluid parameterization is physically intuitive and adequately describes a range of DE models; however, in scenarios where 1 + w crosses zero (dubbed phantom crossing), Eq. (7) becomes ill-behaved. In such cases, an alternative prescription is required.

2.2.2. Parametrized post-Friedmann dark energy

The PPF formalism was developed to address the shortcomings of the fluid approach near the phantom crossing. In the PPF approach, DE perturbations are recast in terms of a new variable, Γ, defined in the DE rest frame by

$$\begin{array}{c}\hfill k^2\mathrm{\Gamma }\equiv -4\pi G{a}^2\delta {\rho }_{\mathrm{DE}}^{\mathrm{rest}}.\end{array}$$(9)

In the conformal Newtonian gauge, the Poisson equation is written as

$$\begin{array}{c}\hfill k^2\varphi =-4\pi G{a}^2[\delta {\rho }_{\mathrm{total}}^{\mathrm{N}}-3\mathcal{H}({\rho }_{\mathrm{total}}+p_{\mathrm{total}})\frac{{\theta }_{\mathrm{total}}^{\mathrm{N}}}{k^2}],\end{array}$$(10)

where the subscript “total” tautologically indicates the contributions from all species, not only from DE, and the superscript N denotes the Newtonian gauge. One can split the contributions of DE from the other species using that

$$\begin{array}{c}\hfill (1+w_{\mathrm{total}}){\theta }_{\mathrm{total}}={\displaystyle \sum _i}(1+w_i)\frac{{\rho }_i}{{\rho }_{\mathrm{total}}}{\theta }_i,\end{array}$$(11)

where

$$\begin{array}{c}\hfill w_{\text{total}}={\displaystyle \frac{{\sum }_i{w}_i{\rho }_i}{{\sum }_i{\rho }_i}}\text{and}{\rho }_{\text{total}}={\displaystyle \sum _i}{\rho }_i,\end{array}$$(12)

and the summation is made over DE and the other components. Therefore, Eq. (10) can be rewritten as

$$\begin{array}{c}\hfill k^2\varphi =-4\pi G{a}^2[\delta {\rho }_{\mathrm{t}}^{\mathrm{N}}-3\mathcal{H}({\rho }_{\mathrm{t}}+p_{\mathrm{t}})\frac{{\theta }_{\mathrm{t}}^{\mathrm{N}}}{k^2}]+k^2\mathrm{\Gamma },\end{array}$$(13)

where we use the same notation as Dakin et al. (2019a) and the subscript “t” represents the contributions of all components but DE. Imposing the physical requirements of energy conservation on super-horizon scales and recovering the standard Newtonian Poisson equation on small scales implies that Γ satisfies

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}^{\prime }=\mathcal{H}(S-\mathrm{\Gamma }),\end{array}$$(14)

with

$$\begin{array}{c}\hfill S\equiv \frac{4\pi G{a}^2}{\mathcal{H}}({\rho }_{\mathrm{DE}}+p_{\mathrm{DE}})\frac{{\theta }_{\mathrm{t}}^{\mathrm{N}}}{k^2}.\end{array}$$(15)

To interpolate between the super-horizon regime and the small-scale limit, an effective parameter c_Γ is introduced so that the evolution of Γ is governed by

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}^{\prime }=\mathcal{H}[\frac{S}{1+c_{\mathrm{\Gamma }}^2{k}^2/{\mathcal{H}}^2}-\mathrm{\Gamma }(1+\frac{c_{\mathrm{\Gamma }}^2{k}^2}{{\mathcal{H}}^2})].\end{array}$$(16)

A value of c_Γ ∼ 0.4 c_s reproduces the behavior of quintessence models and is widely adopted. This PPF approach provides a robust description of DE perturbations that remain well behaved even when the fluid description fails.

Dakin et al. (2019a) present a careful study of DE perturbations in N-body simulations by comparing the fluid and the PPF implementations within the CONncept code (Dakin et al. 2019b,a, 2022)⁹ They show that the two methods yield similar results on scales k ≳ 10⁻³h/Mpc. In contrast, the methods deviate significantly on larger scales. In addition, for a sound speed equal to unity, both parameterizations agree with Newtonian predictions at k ≳ 10⁻²h/Mpc, while on larger scales, general relativistic corrections due to DE perturbations cause deviations of tens of percent. Furthermore, percent level differences in the matter power spectrum at k ∼ 10⁻²h/Mpc are only achieved for c_s² ≲ 10⁻².

2.3. The halo mass function

The HMF describes the comoving number density of matter halos as a function of mass and redshift. It is conventionally written in differential form as

$$\begin{array}{c}\hfill \frac{\mathrm{d}n}{\mathrm{d}M}\mathrm{d}M=\frac{{\rho }_{\mathrm{m}}}{M}\nu f(\nu )\mathrm{d}ln\nu ,\end{array}$$(17)

where ρ_m represents the mean comoving matter density of the Universe, and ν is the peak height parameter. This parameter is defined as

$$\begin{array}{c}\hfill \nu =\frac{{\delta }_{\mathrm{c}}}{\sigma (M,z)},\end{array}$$(18)

with δ_c being the critical linear overdensity for spherical collapse extrapolated to redshift z (Kitayama & Suto 1996) and σ(M,z) the mass variance. The mass variance is computed from the linear matter power spectrum P_m(k,z) through

$$\begin{array}{c}\hfill {\sigma }^2(M,z)=\frac{1}{2{\pi }^2}{\displaystyle {\int }_0^{\infty }}dk{k}^2{P}_{\mathrm{m}}(k,z)W^2(kR),\end{array}$$(19)

where R is the Lagrangian radius corresponding to mass M (i.e., $R={[3M/(4\pi {\rho }_{\mathrm{m}})]}^{1/3}$) and W(kR) is the Fourier transform of a top-hat window function.

In Euclid Collaboration: Castro et al. (2023), we demonstrated that the multiplicity function νf(ν) given by

$$\begin{array}{c}\hfill \nu f(\nu )=A(p,q)\sqrt{\frac{2a{\nu }^2}{\pi }}exp(-\frac{a{\nu }^2}{2})[1+\frac{1}{{(a{\nu }^2)}^p}]{(\sqrt{a}\nu )}^{q-1},\end{array}$$(20)

provides flexibility to reproduce the outcomes of several scale-free simulations accurately. In that work, the challenge of capturing the non-universality of the HMF was reduced mainly to describing the cosmological evolution of the free parameters that appear in Eq. (20). It was shown that both the background evolution and the power spectrum shape impact the HMF universality. Their cosmological dependence was modeled as

$$\begin{array}{cc}\hfill a& =a_R{\mathrm{\Omega }}_{\mathrm{m}}^{a_z}(z),\hfill \end{array}$$(21)

$$\begin{array}{cc}\hfill a_R& =a_1+a_2{(\frac{\text{d}ln\sigma }{\text{d}lnR}+0.6125)}^2,\hfill \end{array}$$(22)

$$\begin{array}{cc}\hfill p& =p_1+p_2(\frac{\text{d}ln\sigma }{\text{d}lnR}+0.5),\hfill \end{array}$$(23)

$$\begin{array}{cc}\hfill q& =q_R{\mathrm{\Omega }}_{\mathrm{m}}^{q_z}(z),\hfill \end{array}$$(24)

$$\begin{array}{cc}\hfill q_R& =q_1+q_2(\frac{\text{d}ln\sigma }{\text{d}lnR}+0.5),\hfill \end{array}$$(25)

and we refer the reader to Euclid Collaboration: Castro et al. (2023) for the values of the free parameters.

For the DUCA simulation suite, which spans a considerably broader dynamic range, we found that the previous baseline function lacked the flexibility to capture the simulation data fully and proposed in Paper I the following parametrization for the HMF on a scenario with dynamic DE:

$$\begin{array}{c}\hfill \nu f(\nu )=A(p,q,r)\sqrt{\frac{2{\nu }_{\ast }^2}{\pi }}exp(-\frac{{\nu }_{\ast }^2}{2})\{{\nu }_{\ast }^r+{\left(\frac{{\nu }_p}{{\nu }_{\ast }}\right)}^{2p}\}{\nu }_{\ast }^{q-1},\end{array}$$(26)

where we define ${\nu }_{\ast }=\sqrt{a}\nu$. In this formulation, the normalization constant A is determined by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\left(\frac{A}{\mathrm{\Xi }}\right)}^{-1}=& -2^{p+\frac{r}{2}}q\mathrm{\Gamma }(\frac{q}{2}+\frac{r}{2}+1)+2^{p+\frac{r}{2}+1}p\mathrm{\Gamma }(\frac{q}{2}+\frac{r}{2}+1)\hfill \\ \hfill [1ex]& -q{\nu }_p^{2p}\mathrm{\Gamma }(-p+\frac{q}{2}+1)-r{\nu }_p^{2p}\mathrm{\Gamma }(-p+\frac{q}{2}+1),\hfill \end{array}\end{array}$$(27)

with

$$\begin{array}{c}\hfill \mathrm{\Xi }=\frac{\sqrt{\pi }(2p-q)(q+r)}{2^{-p+\frac{q}{2}+\frac{1}{2}}}.\end{array}$$(28)

In the limits r → 0 and ν_p → 1, Eq. (26) naturally reverts to the multiplicity function given in Eq. (20). The dependence on DE is introduced by allowing the parameters of the multiplicity function, collectively denoted as x_i ∈ {a, p, q, r}, to vary with cosmology while ν_p is cosmology independent. The cosmological dependence of the parameters (x_i) was modeled according to

$$\begin{array}{c}\hfill x_i=x_{\mathrm{\Lambda },i}\{1+\frac{3{\alpha }_i}{2}[w_{\mathrm{DE}}(z_{\mathrm{ta}})+1]{\mathrm{\Omega }}_{\mathrm{DE}}(z_{\mathrm{ta}})\},\end{array}$$(29)

where x_Λ, i is given by Eqs. (21)–(25) for x_i ∈ {a, p, q}, while

$$\begin{array}{c}\hfill r=r_1+r_2(\frac{\text{d}ln\sigma }{\text{d}lnR}+0.5).\end{array}$$(30)

In Eq. (29), z_ta is the turnaround redshift at which a spherical top-hat perturbation, collapsing at redshift z, reaches its maximum expansion. The free parameters α_i modulate the impact of the evolving DE component on the HMF. This formulation links the modifications in the multiplicity function to the epoch when the effects of DE and gravitational collapse are in balance. We adopted the HMF presented in Paper I as our baseline model and, for conciseness, we refer to Paper I for further details on the implementation and values of the free parameters in Eq. (26).

3. Methodology

In this section we present our methodology. We introduce the simulations produced in Sect. 3.1, the halo catalogs in Sect. 3.2, and present the model calibration setup in Sect. 3.3.

3.1. Simulations

We performed 16 cosmological simulations using the CONncept code, employing the same integration accuracy settings as in the DUCA simulation set presented in Paper I. All runs were carried out in a cubic volume with a comoving side length of 2 h⁻¹ Gpc and evolved with 2048³ particles, ensuring a particle mass of roughly 8 × 10¹⁰ h⁻¹ M_⊙, which provides the resolution needed for our analyses in the regime of objects more massive than 2 × 10¹³ h⁻¹ M_⊙. The initial conditions were generated on the fly using third-order Lagrangian perturbation theory (3LPT) at a starting redshift of z = 24.

The background cosmology, with the exception of the DE parameters, is kept fixed in all simulations with the following parameters: a Hubble constant of H₀ = 67 km s⁻¹ Mpc⁻¹, a baryon density parameter Ω_b = 0.049, and a CDM density parameter Ω_cdm = 0.27. The primordial power spectrum is specified by an amplitude A_s = 2.1 × 10⁻⁹, a spectral index n_s = 0.96, zero running (α_s = 0), and a pivot scale of k_pivot = 0.05 Mpc⁻¹.

We varied the DE parameters across the simulation suite to investigate the impact of DE dynamics. In particular, the DE sound speed squared, c_s², is set to one of four values (10⁻⁷, 10⁻⁵, 10⁻³, or 1), thereby spanning a wide range of clustering behaviors. For each choice of c_s², we explored three combinations of the DE EOS parameters by varying w₀ and w_a. Specifically, w₀ takes the values −1.0, −0.9, and −0.8, while the corresponding w_a values are 0.3, 0.2, and 0.1, respectively. This systematic variation results in 12 simulations using the PPF prescription, which enables a comprehensive exploration of the phenomenology associated with dynamical DE.

Table 1 summarizes our simulation suite’s DE parameter combinations, including four additional simulations: simulations 3 and 11 have also been run using the fluid description and its phantom twin – a simulation with the DE EOS reflected with respect to −1. We selected the sound speed values in our simulations to examine various CDE levels while working within the constraints of the linear-species implementation in CONncept. These constraints impose a lower bound on c_s as the realization method cannot accurately simulate values that are too low. In Hassani et al. (2019), the linear-realization scheme was compared with an alternative EFT approach for DE perturbations. The resulting difference in the nonlinear matter power spectrum is below 1 percent for c_s² = 10⁻⁷ at z = 0 and decreases further at higher redshift. While this finding motivates our choice of c_s, it is essential to remember that our implementation remains at the level of linear DE perturbations, capturing the impact of these perturbations on the nonlinear collapse of dark matter but, by construction, it excludes any backreaction of nonlinear matter clustering on the DE field. Furthermore, the DE EOS can also influence the CDE nonlinear behavior (Batista et al. 2023; Nouri-Zonoz et al. 2025). Therefore, our lowest-c_s results should be regarded as a baseline, pending more advanced treatments of nonlinear CDE.

3.2. Halo identification

In our study, we made use of the ROCKSTAR¹⁰ halo finder in combination with the CONSISTENT Trees algorithm¹¹. ROCKSTAR is an advanced phase-space halo finder that employs both spatial and velocity information to detect dark matter halos (Behroozi et al. 2013a).

The procedure commences by dividing the simulation volume into friends-of-friends (FoF) groups using a relatively large linking length of 0.28 (expressed in units of the mean inter-particle separation), which exceeds the conventional value of approximately 0.20 typically used in other FoF algorithms. This initial grouping facilitates efficient parallel processing. Within each FoF group, ROCKSTAR performs an adaptive hierarchical refinement in six-dimensional phase space (three dimensions for position and three for velocity), resulting in a nested hierarchy of subgroups. This approach permits the precise identification of halos and their substructures, even in regions with high particle density.

To improve both the temporal consistency and the accuracy of halo properties, we processed the ROCKSTAR outputs with the CONSISTENT trees algorithm (Behroozi et al. 2013b). This algorithm tracks halos across successive simulation snapshots, ensuring that properties such as halo masses and positions evolve smoothly over time. Such dynamic tracking is essential for reducing scatter and enhancing the reliability of the halo catalogs in our analysis.

We adopted the spherical overdensity criterion for defining halos, whereby the mass of a halo is determined by the mean density within a sphere that is equal to the virial overdensity Δ_vir(z). ROCKSTAR uses the fitting formula of Bryan & Norman (1998) for Δ_vir and considers only the contribution of matter. In the case of CDE, including the DE contribution within the overdensity radius could, in principle, improve the universality of the HMF, as suitable choices of the collapse threshold are known to preserve universality in homogeneous DE scenarios (Diemer 2020). However, a self-consistent CDE spherical collapse calculation would introduce significant computational overhead and complicate the HMF evaluation, and we therefore retained the standard matter-only Bryan & Norman (1998) fit as implemented in ROCKSTAR.

ROCKSTAR defines the spherical overdensity center as the average positions of the particles belonging to the innermost subgroup in the hierarchy. All particles within the virial radius are included in the halo mass, irrespective of their gravitational binding.

Finally, the halo catalogs generated by the combined use of ROCKSTAR and CONSISTENT Trees are post-processed, counting the number of objects as a function of halo mass in logarithmically spaced bins based on the number of particles per halo. This binning strategy minimizes the effects of mass discretization and ensures a statistically robust sampling across the halo mass range. Only halos containing more than 300 particles are considered. This selection cut in mass is sufficient to ensure sub-percent agreement with higher-resolution simulations (see, for instance, Euclid Collaboration: Castro et al. 2023).

3.3. Model calibration

We adopted a Bayesian approach to calibrate the HMF model. The calibration procedure follows that of Euclid Collaboration: Castro et al. (2023) and in Paper I; the reader is referred to those works for a complete discussion. Here, we summarize the main aspects.

The calibration fits the theoretical prediction as a function of the HMF parameters θ

$$\begin{array}{c}\hfill N_i^{\mathrm{th}}=N_i(\theta ,z),\end{array}$$(31)

obtained by integrating the HMF model over the mass bin and multiplying by the simulation volume, to the halo counts $N_i^{\mathrm{sim}}$ extracted from our simulations (see Sect. 3.2). The number of halos in each mass bin is assumed to follow a Poisson distribution as motivated by the Press–Schechter formalism (Press & Schechter 1974). The Poisson log-likelihood is given by

$$\begin{array}{c}\hfill ln\mathcal{L}(N_i^{\mathrm{sim}}\mid \theta ,z)=N_i^{\mathrm{sim}}ln{N}_i^{\mathrm{th}}-N_i^{\mathrm{th}}-ln(N_i^{\mathrm{sim}}!).\end{array}$$(32)

To account for systematic effects and numerical uncertainties, a composite likelihood is adopted. In bins with a large number of halos (i.e., $N_i^{\mathrm{sim}}>25$), the Poisson likelihood is approximated by a Gaussian distribution (Euclid Collaboration: Castro et al. 2023). The composite log-likelihood is expressed as

$$\begin{array}{c}\hfill ln\mathcal{L}(N_i^{\mathrm{sim}}\mid \theta ,z)=\{\begin{array}{cc}\begin{array}{cc}& N_i^{\mathrm{sim}}ln{N}_i^{\mathrm{th}}-N_i^{\mathrm{th}}\hfill \\ \hfill & -ln(N_i^{\mathrm{sim}}!)\hfill \end{array}\hfill & \text{if}{N}_i^{\mathrm{sim}}\le 25,\hfill \\ [5ex]\begin{array}{cc}& -\frac{1}{2}{\left(\frac{N_i^{\mathrm{sim}}-N_i^{\mathrm{th}}}{{\sigma }_i}\right)}^2\hfill \\ \hfill & -\frac{1}{2}ln(2\pi {\sigma }_i^2)\hfill \end{array}\hfill & \text{if}{N}_i^{\mathrm{sim}}>25,\hfill \end{array}\end{array}$$(33)

with the standard deviation defined as

$$\begin{array}{c}\hfill {\sigma }_i^2=N_i^{\mathrm{th}}(1+N_i^{\mathrm{th}}{\sigma }_{\mathrm{sys}}^2).\end{array}$$(34)

In this formulation, σ_sys represents the additional variance component to account for systematic uncertainties and is fixed to the best fit presented in Paper I, namely 8.82 × 10⁻². This correction is not included in the Poissonian part of the likelihood.

The overall log-likelihood was computed by summing over all mass bins, redshifts, and simulation outputs:

$$\begin{array}{c}\hfill ln{\mathcal{L}}_{\mathrm{total}}(\theta )={\displaystyle \sum _s\sum _z\sum _i}ln\mathcal{L}(N_{i,z,s}^{\mathrm{sim}}\mid \theta ,z).\end{array}$$(35)

In the above equation, the index s labels the simulation outputs. Although outputs from the same simulation at different redshifts are not strictly independent, this issue is mitigated by selecting snapshots separated by intervals exceeding the typical dynamical time of galaxy clusters (approximately 1.7 Gyr) in accordance with Bocquet et al. (2016). The best-fit value for θ is determined by the global maximum identified with the DUAL ANNEALING algorithm (Tsallis 1988; Tsallis & Stariolo 1996; Xiang et al. 1997; Xiang & Gong 2000) as implemented in SCIPY (Virtanen et al. 2020).

4. Results

4.1. Impact of a variable sound speed

Figure 1 illustrates how modifications in the DE sound speed affect the HMF at z = 0. Not surprisingly, the impact of CDE on the HMF is more substantial the lower the sound speed, as, in this case, its clustering is higher. While in the case for c_s² = 10⁻³ we only observe an increasing scatter of the number counts at higher masses, for c_s² = 10⁻⁵ and c_s² = 10⁻⁷ we observe a slight boost in the number of objects on the order of a few percent for the most massive objects.

Fig. 1. Relative differences in HMF at redshift zero for varying DE sound speeds (cs2), normalized to the reference case cs2 = 1. Different colors correspond to different DE EOS parameters, while different panels correspond to different sound speeds. From left to right, cs2 is 10−7, 10−5, and 10−3. Gray bands mark ±1% (dark) and ±2% (light) differences.

From Fig. 1, we also observe that the larger the deviation from w = −1, the higher the impact; again, as in this condition, DE perturbations are more important. To verify that, we present in Fig. 2 the DE power spectrum computed using CONncept. Therefore, the higher the clustering in DE, the higher the impact on the HMF.

Fig. 2. DE power spectrum at redshift zero and sound speeds cs2 = 10−7 and cs2 = 10−7 for the different DE EOS parameters.

4.2. Extended model

In Fig. 3 we compare the relative differences in the HMF at z = 0 for c_s² = 10⁻⁷ with respect to the reference case c_s² = 1, alongside the prediction of our baseline model from Paper I. The baseline model accounts for the impact of CDE on clustering via Eq. (19), where P_m is computed directly from simulations with c_s² = 10⁻⁷. Although it qualitatively reproduces the simulation results, it overestimates the effect by a factor of a few.

Fig. 3. Relative differences in the HMF at z = 0 for the case where cs2 = 10−7, compared to the reference case cs2 = 1, for different DE EOS parameters (filled lines in different colors). The figure also shows the prediction of this quantity from our baseline model, presented in Paper I, as indicated by the dotted lines.

To understand the reason for the mismatch between the baseline model and the simulations, we have to keep in mind the key role played by the mapping between linear theory and nonlinear clustering within the HMF formalism. Relativistic species, which couple to matter only through gravity, cluster significantly less than matter. As a result, linear theory overestimates the actual impact that CDE has on nonlinear matter clustering. Consequently, using the linear matter power spectrum for c_s² < 1 to compute the mass variance in Eq. (19) and to predict the number of collapsed structures leads to an overestimate of the effect of CDE. A similar conclusion was shown by Nouri-Zonoz et al. (2025) studying the linear and nonlinear coupling between DE and matter perturbations.

To reconcile this issue with the sensitivity of the HMF to CDE clustering, we propose the following empirical modification to the computation of the peak height:

$$\begin{array}{c}\hfill \nu (M,z|c{}_{\mathrm{s}}^2<1)\to \frac{\nu (M,z|c{}_{\mathrm{s}}^2=1)}{1+\lambda \mathrm{sgn}(1+w(z_{\mathrm{ta}}))\frac{{\mathrm{\Omega }}_{\mathrm{DE}}(z_{\mathrm{ta}}){\sigma }_{\mathrm{DE}}(R_{\mathrm{ta}},z_{\mathrm{ta}}|c{}_{\mathrm{s}}^2<1)}{{\mathrm{\Omega }}_{\mathrm{m}}(z_{\mathrm{ta}}){\sigma }_{\mathrm{m}}(R,z_{\mathrm{ta}}|c{}_{\mathrm{s}}^2<1)}\frac{{\delta }_{\mathrm{ta}}^{\mathrm{EdS}}}{{\mathrm{\Delta }}_{\mathrm{ta}}^{\mathrm{EdS}}}}.\end{array}$$(36)

Here σ_DE and σ_m are the variance of DE and nonrelativistic matter perturbations, respectively, sgn is the sign function, and the subscript “ta” indicates that the quantities are computed at the turnaround. As in Paper I, the Einstein–de Sitter (EdS) values are used for the turnaround redshift as well as for the linear and nonlinear density contrasts, ${\delta }_{\mathrm{ta}}^{\mathrm{EdS}}$ and ${\mathrm{\Delta }}_{\mathrm{ta}}^{\mathrm{EdS}}$, respectively. Numerically, they are

$$\begin{array}{cc}\hfill z_{\mathrm{ta}}& =\frac{1+z}{{(1/2)}^{2/3}}-1,\hfill \end{array}$$(37)

$$\begin{array}{cc}\hfill {\delta }_{\mathrm{ta}}^{\mathrm{EdS}}& =\frac{9{\pi }^2}{16},\hfill \end{array}$$(38)

$$\begin{array}{cc}\hfill {\mathrm{\Delta }}_{\mathrm{ta}}^{\mathrm{EdS}}& =\frac{3}{5}{\left(\frac{3\pi }{4}\right)}^{2/3}.\hfill \end{array}$$(39)

The following considerations motivate the terms in Eq. (36). The sign function, $\mathrm{sgn}(1+w(z_{\mathrm{ta}}))$, is introduced because when the DE EOS crosses −1, the pressure contribution overtakes the density contribution in affecting the gravitational potential; this leads to a de-boost in matter clustering. This is illustrated in Fig. 4, where we note that a region of overdensity of matter corresponds to an overdensity of DE for simulation 3, but to an underdensity of DE on its phantom-twin counterpart, in agreement with the results of Batista et al. (2023). In addition, the ratio of the density parameters, Ω_DE(z_ta)/Ω_m(z_ta), is included to weight the DE perturbations relative to those of matter, as the gravitational impact is governed by the absolute density perturbations–as reflected in the Euler equation in the N-body gauge–rather than by the density contrast alone. The ratio of variances, σ_DE(R_ta, z_ta|c_s² < 1)/σ_m(R, z_ta|c_s² < 1), quantifies the relative amplitude of DE fluctuations compared to matter fluctuations. In this context, evaluating the variance for DE fluctuations at the turnaround scale is crucial because the late onset of DE clustering implies that, when DE becomes dynamically important, the collapsing matter is already confined to a smaller volume than the original Lagrangian patch R. As CDE perturbations are well described by linear theory virtually at all scales and epochs, the difference between Eulerian and Lagrangian coordinates of a CDE fluid element is always small. Therefore, to assess the impact that CDE has on the dark matter collapse, one has to compare the DE fluctuations at the scale when DE is dynamically and kinematically relevant for the collapse – in our model, the turnaround radius – to the fluctuations of matter perturbations linearly extrapolated to the epoch of turnaround. Note also, from Eq. (6), that DE perturbations go to zero when the EOS parameter approaches −1. Finally, the factor ${\delta }_{\mathrm{ta}}^{\mathrm{EdS}}/{\mathrm{\Delta }}_{\mathrm{ta}}^{\mathrm{EdS}}$ is introduced to rescale the matter variance computed from the linear power spectrum, correcting for the mismatch that arises because DE perturbations remain in the linear regime at the turnaround. In contrast, matter perturbations have already entered the nonlinear regime. Finally, the empirical constant λ calibrates the overall strength of the DE clustering effect. Together, these modifications yield a peak height that more accurately reflects the interplay between DE clustering and matter perturbations in structure formation.

Fig. 4. Density fields in a comoving volume of size 1 Gpc/h at z = 0. Left: Matter distribution for simulation 3, i.e., w0 = −0.8 and wa = 0.1. Middle and right: DE distributions for the same simulation and its phantom twin (w0 = −1.2,  wa = −0.1), respectively. A zoomed-in view of the region inside the red rectangles is displayed in an inset at the lower right.

4.3. Model calibration and validation

In Fig. 5 we present the model calibrated on simulations 1–12 using the PPF description presented in Table 1. The best fit for λ obtained by maximizing the likelihood presented in Sect. 3.3 is λ = 0.34. We have not observed a statistically significant evolution of λ with redshift or c_s. We compared the model prediction against the simulations for c_s² ∈ {10⁻⁷, 10⁻⁵} as larger sound speeds resulted in an effect that is too small, as we show in Fig. 1. Different rows present the comparison at various redshifts.

Fig. 5. Comparison between the model calibrated on simulations 1–12 using the PPF description presented in Table 1 (dotted) and the simulations for cs2 ∈ {10−7, 10−5} (solid lines). Different DE EOSs are presented in different colors, and the rows present the comparison and respective relative residuals at the selected redshifts, z ∈ {0.0, 0.3, 0.9, 2.0}.

From Fig. 5, we observe that the relative impact of CDE is below one percent for all scenarios considered for redshift z ≳ 1. From the plots of the residuals, we note that the calibrated model reproduces the simulation with sub-percent accuracy, except in the high mass end, where the count errors become evident. Since the simulations share the same seed, the ratio concerning the c_s² = 1 case suppresses the Poisson noise, making the small impact of the CDE assessable. Still, when the number of halos inside a mass bin becomes small, the ratio of the simulations presented in Fig. 5 starts fluctuating more significantly around the model.

In Fig. 6 we present our model prediction for simulation 3 and its phantom twin. While the former has been used to calibrate our model, the latter is a real prediction of our model. We observe that our model’s performance is consistent for simulation 3 and its twin. We also observe that the amplitude of the overall impact is not simply reflected about zero, with the phantom impact being smaller by a factor of two. This can be understood as the contribution of DE perturbations to the Poisson equation is driven by

Fig. 6. Comparison between the best-fit model prediction for the cosmological parameters from simulation 3 (w0 = −0.8, wa = 0.1) and its phantom twin (w0 = −1.2, wa = −0.1) for cs2 = 10−7. Different columns present the comparison at the selected redshifts, z ∈ {0.0, 0.3, 0.9}.

$$\begin{array}{c}\hfill \delta {\rho }_{\mathrm{DE}}^{\mathrm{Nb}}=\delta {\rho }_{\mathrm{DE}}^{\mathrm{S}/\mathrm{N}}+3\mathcal{H}(1+w_{\mathrm{DE}})\frac{{\theta }_{\mathrm{total}}^{\mathrm{S}/\mathrm{N}}}{k^2}{\overline{\rho }}_{\mathrm{DE}},\end{array}$$(40)

where the super-scripts refer to N-body, synchronous, and Newtonian gauges, respectively. Thus, the contribution of the velocity divergence field ${\theta }_{\mathrm{total}}^{S/N}$ flips its sign when phantom crossing. Starting from adiabatic initial conditions, where density fluctuations of all species couple to the same curvature fluctuations, the phantom density perturbations δρ_DE and divergence of the velocity field initially have different signs, causing a slower evolution with respect to the non-phantom twin. Moreover, in phantom models, the background energy density decays faster with redshift than in non-phantom models, which further limits the impact of phantom DE perturbations.

Lastly, in Fig. 7 we compare the differences in the HMF between the PPF and fluid descriptions. We present the predictions of the best-fit model for c_s² = 10⁻⁷ at redshift z = 0 where the differences are the largest. Still, the predicted relative difference is the order of 10⁻⁴ below both the accuracy and precision our simulations can probe and the overall effect of CDE.

Fig. 7. Relative difference in the HMF between the fluid and PPF approaches plotted as a function of the halo mass. Results are shown at redshift z = 0 for the cosmologies corresponding to simulations 1–3. The solid lines denote the original cosmology, while dotted lines indicate the corresponding phantom twin.

5. Discussion

An interesting result of our analysis is that previous models used to assess the impact of CDE on the HMF tended to predict a larger effect than what is observed in our simulations. For instance, Batista & Marra (2017) predicted that CDE could affect the abundance of clusters by up to 30 percent. In contrast, our simulations reveal a maximum relative difference of only a few percent at the mass scale of 10¹⁵ h⁻¹ M_⊙.

To reconcile these divergent results, it is important to note that Batista & Marra (2017) employed the spherical collapse model framework to describe the nonlinear evolution of DE, assuming a null sound speed. This approach presumes that the nonlinearities between DE and matter remain correlated throughout their coevolution. However, as demonstrated in Nouri-Zonoz et al. (2025), the correlation between the two fields decays rapidly once DE becomes nonlinear. Moreover, as shown in Batista et al. (2023), a value of c_s² ≤ 10⁻⁷ can be considered effectively as a null sound speed at cluster scales, leading to significant DE contrasts near halo virialization overdensities. Therefore, given that our current analysis assumes linear DE perturbations, our results should be regarded as a lower limit on the impact of DE fluctuations. Likely, the accurate effect for c_s² ∼ 10⁻⁷ lies between our method and the predictions based on the spherical collapse approach.

Another interesting effect is that CDE is enhanced in non-phantom DE scenarios. In our analysis, models with an EOS that remains above the phantom divide exhibit a more pronounced impact on the HMF. This behavior is expected because DE perturbations contribute more effectively to the gravitational potential when w_DE > −1.

It is also noteworthy that CDE could potentially induce an effect analogous to the tension observed in σ₈ between cosmic microwave background and large-scale structure surveys (Battye et al. 2015; Douspis et al. 2018). In fact, DE perturbations can either boost or de-boost the matter fluctuations, making the collapse of structures easier or harder. The signal of the impact, as depicted in Eq. (36) and Fig. 4, depends on the DE EOS.

However, the overall impact remains very modest for the scenarios simulated with our current numerical setup. For instance, in our simulation 3, which has the largest CDE impact, at redshift z = 0, the bias between the fiducial value of A_s and the one that reproduces the number of halos more massive than 10¹⁴ h⁻¹ M_⊙ ignoring CDE is below one percent. The bias reaches roughly one percent if we extrapolate our model to c_s² = 0, but the behavior of our model beyond its validity is not easily assessable and, therefore, not encouraged. A different approach would be necessary to reliably explore the extreme case of c_s² = 0.

Finally, the variations stemming from assuming the fluid or PPF description for DE are much smaller than the already limited impact of CDE. This finding suggests that uncertainties in DE microphysics are unlikely to introduce significant biases in the predictions for halo abundances when CDE’s dominant contributions are adequately considered.

6. Conclusions

We have presented an improved HMF model that explicitly incorporates the impact of CDE with a variable sound speed. Our analysis is based on an extended suite of N-body simulations that cover a wide range of DE parameters and sound speeds. The main conclusions of our study are as follows:

• The effect of CDE on the HMF depends sensitively on the DE EOS. In models with w_DE > −1, DE perturbations boost gravitational collapse, whereas in phantom scenarios (w_DE < −1) the opposite effect is observed.

• The impact of CDE is more pronounced in non-phantom DE scenarios, where DE perturbations contribute more effectively to the gravitational potential.

• When the matter power spectrum computed for c_s² < 1 is applied directly to the baseline model of Paper I, the qualitative behavior agrees with the simulation outcomes; however, the impact is over-predicted by roughly a factor of two.

• We argue that this overprediction arises because DE clusters much less efficiently than matter. Therefore, extrapolating the impact of CDE from the linear power spectrum exaggerates its effect on the nonlinear collapse of structures.

• To include the effect of CDE on halo collapse, we introduced an effective peak height [Eq. (36)] that is defined in terms of the matter power spectrum for the c_s² = 1 case and modulated by the relative contribution of CDE perturbations at turnaround. This formulation more accurately captures the interplay between DE and matter during collapse.

• Although CDE could, in principle, generate an effect analogous to the σ₈ tension observed between cosmic microwave background and large-scale structure measurements, the overall impact remains modest for the range of sound speeds explored in our simulations. Extrapolation to c_s² = 0 suggests a larger bias; however, this extreme regime would require an alternative numerical treatment as the linear realization methodology will not be valid if CDE starts to deviate strongly from linearity.

• Finally, the difference between the PPF and fluid approaches for DE is much smaller than the limited impact of CDE itself. Hence, uncertainties in DE microphysics are unlikely to introduce significant biases in halo abundance predictions when the dominant effects of CDE are properly modeled.

Our extended HMF model maintains sub-percent accuracy in reproducing simulation outcomes and provides a robust framework for interpreting forthcoming galaxy cluster surveys. Future work should focus on alternative simulation approaches to reliably probe the extreme c_s² = 0 regime and further refine the model in light of next-generation observational data.

Data availability

The baseline model used in this paper is presented in Castro & Fumagalli (2024) (https://github.com/TiagoBsCastro/CCToolkit). The code to compute the DE power spectrum is presented in Castro & Dakin (2025) (https://github.com/TiagoBsCastro/ConceptSpectra).

Acknowledgments

It is a pleasure to thank Duca dos Anjos. We warmly thank Francesco Pace for sharing solutions to the spherical collapse with us. We thank Isabella Baccarelli, Fabio Pitari, and Caterina Caravita for their support with the CINECA environment. We acknowledge the CINECA award for the availability of high-performance computing resources and support as part of the Leonardo Early Access Program (LEAP) and under the ISCRA initiative. RCB, TC, VM also thank Conselho Nacional de Desenvolvimento Científico e Tecnológico – CNPq (process 444368/2024-8) for partially supporting this work. TC and SB are supported by the Agenzia Spaziale Italiana (ASI) under – Euclid-FASE D Attivita’ scientifica per la missione – Accordo attuativo ASI-INAF n. 2018-23-HH.0, by the National Recovery and Resilience Plan (NRRP), Mission 4, Component 2, Investment 1.1, Call for tender No. 1409 published on 14.9.2022 by the Italian Ministry of University and Research (MUR), funded by the European Union – NextGenerationEU– Project Title “Space-based cosmology with Euclid: the role of High-Performance Computing” – CUP J53D23019100001 – Grant Assignment Decree No. 962 adopted on 30/06/2023 by the Italian Ministry of University and Research (MUR); by the Italian Research Center on High-Performance Computing Big Data and Quantum Computing (ICSC), a project funded by European Union – NextGenerationEU – and National Recovery and Resilience Plan (NRRP) – Mission 4 Component 2, by the INFN INDARK PD51 grant, and by the PRIN 2022 project EMC2 – Euclid Mission Cluster Cosmology: unlock the full cosmological utility of the Euclid photometric cluster catalog (code no. J53D23001620006). VM acknowledges partial financial support FAPES (Brazil).

References

All Tables

All Figures

	Fig. 1. Relative differences in HMF at redshift zero for varying DE sound speeds (cs2), normalized to the reference case cs2 = 1. Different colors correspond to different DE EOS parameters, while different panels correspond to different sound speeds. From left to right, cs2 is 10−7, 10−5, and 10−3. Gray bands mark ±1% (dark) and ±2% (light) differences.
In the text

	Fig. 2. DE power spectrum at redshift zero and sound speeds cs2 = 10−7 and cs2 = 10−7 for the different DE EOS parameters.
In the text

	Fig. 3. Relative differences in the HMF at z = 0 for the case where cs2 = 10−7, compared to the reference case cs2 = 1, for different DE EOS parameters (filled lines in different colors). The figure also shows the prediction of this quantity from our baseline model, presented in Paper I, as indicated by the dotted lines.
In the text

	Fig. 4. Density fields in a comoving volume of size 1 Gpc/h at z = 0. Left: Matter distribution for simulation 3, i.e., w0 = −0.8 and wa = 0.1. Middle and right: DE distributions for the same simulation and its phantom twin (w0 = −1.2,  wa = −0.1), respectively. A zoomed-in view of the region inside the red rectangles is displayed in an inset at the lower right.
In the text

	Fig. 5. Comparison between the model calibrated on simulations 1–12 using the PPF description presented in Table 1 (dotted) and the simulations for cs2 ∈ {10−7, 10−5} (solid lines). Different DE EOSs are presented in different colors, and the rows present the comparison and respective relative residuals at the selected redshifts, z ∈ {0.0, 0.3, 0.9, 2.0}.
In the text

	Fig. 6. Comparison between the best-fit model prediction for the cosmological parameters from simulation 3 (w0 = −0.8, wa = 0.1) and its phantom twin (w0 = −1.2, wa = −0.1) for cs2 = 10−7. Different columns present the comparison at the selected redshifts, z ∈ {0.0, 0.3, 0.9}.
In the text

	Fig. 7. Relative difference in the HMF between the fluid and PPF approaches plotted as a function of the halo mass. Results are shown at redshift z = 0 for the cosmologies corresponding to simulations 1–3. The solid lines denote the original cosmology, while dotted lines indicate the corresponding phantom twin.
In the text