Open Access
Issue
A&A
Volume 691, November 2024
Article Number A62
Number of page(s) 20
Section Cosmology (including clusters of galaxies)
DOI https://doi.org/10.1051/0004-6361/202451230
Published online 30 October 2024

© The Authors 2024

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 structure-formation process in the Universe is hierarchical, with smaller structures collapsing and merging to form larger ones. Galaxy clusters, the most massive virialized objects in the Universe, lie at the apex of this hierarchy. They serve as valuable cosmological probes, offering insights into the growth of density perturbations and the geometry of the Universe (see, for instance, Allen et al. 2011; Kravtsov & Borgani 2012, for reviews).

The cosmological exploitation of cluster surveys is primarily based on number count analysis. This method involves comparing the observed number of clusters in a survey, as a function of redshift and of a given observable quantity, to the theoretical prediction of the halo mass function (HMF) within a cosmological model, thus enabling the derivation of constraints on cosmological parameters. Numerous studies have been conducted in this area (e.g., Borgani et al. 2001; Holder et al. 2001; Rozo et al. 2010; Hasselfield et al. 2013; Planck Collaboration XX 2014; Bocquet et al. 2015; Mantz et al. 2015; Planck Collaboration XXIV 2016; Bocquet et al. 2019; Abbott et al. 2020; Costanzi et al. 2021; Lesci et al. 2022a).

Complementing cluster count analysis is cluster clustering statistics, which examines the spatial distribution of clusters in the Universe (Mana et al. 2013; Castro et al. 2016; Baxter et al. 2016; To et al. 2021; Lesci et al. 2022b; Sunayama et al. 2023; Romanello et al. 2024; Fumagalli et al. 2024). The halo bias (HB) is a fundamental concept in this analysis since it reflects the ratio between the number overdensity of the cluster and of the matter distribution. This relationship is expected to bring cosmological information through the mass- and redshift-dependence of the HB and to be linear on large scales (i.e., ≳ 100 Mpc) to guarantee the scale-independence of the HB.

The Euclid mission (Laureijs et al. 2011; Euclid Collaboration 2022, 2024c) is projected to provide significant advancements to cluster cosmology. Sartoris et al. (2016) have forecasted that the combined cluster count and clustering analysis by the Euclid mission will provide constraints on the amplitude of the matter power spectrum and the mass density parameters independent and competitive with other cosmological probes, underlining the potential of galaxy clusters as cosmological probes for ongoing and future missions.

At the heart of cluster cosmology are the theoretical models for the HMF and HB. Simplified models based on linear perturbation theory and spherical collapse have provided invaluable insights into the potential of cluster counts and clustering as cosmological probes (see, e.g., Press & Schechter 1974; Bond et al. 1991). However, given the complexity and strongly non-linear nature of cluster formation dynamics, a refinement of these models to the precision level required by available and forthcoming surveys has to rely on cosmological simulations as the primary method to capture such complexity. Several studies have been dedicated to calibrating semi-analytical models for the HMF and HB, aiming to align these models’ predictions with the results from extensive sets of simulations (see, for instance, Sheth & Tormen 1999; Sheth et al. 2001; Jenkins et al. 2001; Warren et al. 2006; Tinker et al. 2008, 2010; Bhattacharya et al. 2011; Watson et al. 2013; Despali et al. 2016; Comparat et al. 2017; Euclid Collaboration 2023). These simulations not only accurately describe the gravitational interactions that predominantly drive structure formation, they but also attempt to account for the effects of baryonic matter.

The influence of baryons, albeit a minor component in the Universe’s overall composition, plays a significant role in the formation of structures, particularly in the context of these simulations (Cui et al. 2014; Velliscig et al. 2014; Bocquet et al. 2016; Castro et al. 2020). Given the sensitivity of baryon evolution to the inclusion and modeling of astrophysical processes occurring at scales much smaller than those resolved in simulations, the modeling of baryonic feedback in hydrodynamical simulations remains a subject of active debate. At the scale of galaxy clusters, for instance, baryonic feedback is known to reorganize the mass density profile of halos without disrupting the structures, thereby altering the mass enclosed within a given radius compared to predictions from collisionless N-body simulations. Owing to the substantially greater computational demands of hydrodynamical simulations, it has become standard practice to derive the HMF from gravitational N -body simulations, with subsequent post-processing to account for baryonic effects (see, e.g. Schneider & Teyssier 2015; Aricò et al. 2021). In this paper, we concentrate on the initial step of calibrating the HB using collisionless simulations. This approach is intended to be a foundational phase, with the baryonic effects being integrated later, and we employed a methodology akin to that used for the HMF (see, Castro et al. 2020; Euclid Collaboration 2024a). This strategy underscores our commitment to systematically exploring the cosmological parameter space, acknowledging the importance of baryonic effects while methodically building toward their inclusion in our analysis.

Systematic errors in the calibration of the HMF and HB can significantly impact the final cosmological constraints. Studies such as those by Salvati et al. (2020), Artis et al. (2021), and Euclid Collaboration (2023) have highlighted how inaccuracies in theoretical models can propagate biases into cosmological parameter inferences. In response to these challenges, Euclid Collaboration (2023) presented a new, rigorously studied calibration of the HMF based on a suitably designed set of N -body simulations, offering the required accuracy to analyze Euclid cluster count data.

Semi-analytical modeling typically starts with a simplified physical model, such as the peak-background split (PBS; Mo & White 1996), which is then extended and refined by adding more degrees of freedom. These additional degrees of freedom are subsequently fitted to simulations. Conceptually, PBS links the HB to the HMF by decomposing the density field into high- and low-frequency modes. The high-frequency modes that cross the collapsing barrier describe the collapse of structures. In contrast, the low-frequency modes modulate the density field fluctuations, thereby enhancing the number of peaks that cross the collapse threshold, therefore linking the clustering of collapsed objects with the local density field. Despite its qualitative consistency with simulations, the quantitative precision of PBS must be enhanced, especially in the context of the Euclid mission’s requirements.

In this paper, we address the challenge of enhancing the accuracy of HB predictions to the level required to fully exploit the cosmological potential of the two-point clustering statistics from the Euclid photometric cluster survey. Our approach involves calibrating a semi-analytical model to quantify the discrepancies between PBS predictions and simulation results. This calibration aims to refine HB predictions, improving the reliability of cosmological parameter estimation derived from cluster counts and clustering.

This paper is organized as follows. We revisit the theoretical aspects used in this paper in Sect. 2. In Sect. 3, we describe the methodology used in our analysis. We present the HB model and its calibration in Sect. 4 along with an assessment of our model’s impact in a forecast Euclid cluster cosmology analysis. Final remarks are made in Sect. 5. The implementation of our model is publicly available at https://github.com/TiagoBsCastro/CCToolkit and presented in Sect. 5.

2 Theory

2.1 The halo mass function

The differential HMF is given by dn dM dM=ρmMvf(v)d ln v,${{{\rm{d}}n} \over {{\rm{d}}M}}{\rm{d}}M = {{{\rho _{\rm{m}}}} \over M}vf(v){\rm{d}}\ln v,$(1)

where n is the comoving number density of halos with mass in the range [M, M + dM], v is the peak height, and the function vf (v) is known as the multiplicity function. The term pm is the comoving cosmic mean matter density, ρm=3H02Ωm,08πG,${\rho _{\rm{m}}} = {{3H_0^2{\Omega _{{\rm{m}},0}}} \over {8\pi G}},$(2)

where H0 and Ωm,0 are the current value of the Hubble parameter and the matter density parameter, and G is the gravitational constant. The peak height is defined as v = δc/σ(M, z), where δc is the critical density for spherical collapse (Peebles 2020) and σ2(M, z) is the filtered mass variance at redshift z so that it measures how rare a halo is. The mass variance is expressed in terms of the linear matter power spectrum Pm(k, z) as σ2(M,ɀ)=12π20dkk2Pm(k,ɀ)W2(kR),${\sigma ^2}(M,) = {1 \over {2{\pi ^2}}}\mathop \smallint \limits_0^\infty {\rm{d}}k\,{k^2}{P_{\rm{m}}}(k,){W^2}(kR),$(3)

where R(M) = [3 M / (4 πρm)1/3 is the Lagrangian radius of a sphere containing the mass M, and W(kR) is the Fourier transform of a top-hat filter of radius R.

The multiplicity function is considered universal if its cosmological dependence is solely through the peak height. However, numerous studies based on N-body simulations have challenged this assumption. These analyses reveal that, while the initial approximation of HMF universality is generally valid, systematic deviations from this universality become evident at late times in the universe’s evolution. This deviation has been demonstrated in various independent investigations, each indicating a nuanced understanding of the HMF’s behavior (e.g., Crocce et al. 2010; Courtin et al. 2011; Watson et al. 2013; Diemer 2020; Ondaro-Mallea et al. 2021; Euclid Collaboration 2023).

The non-universality of the HMF is affected by both the halo definition and the residual dependence of Ac on cosmology. Various studies have shown this dependence, including Watson et al. (2013), Despali et al. (2016), Diemer (2020), and Ondaro- Mallea et al. (2021) for the dependence on the halo definition, and Courtin et al. (2011) for the cosmology dependence of δc. In our study, we define halos as spherical overdensities (SO) with an average enclosed mean density equal to Δvir(z) times the background density, where Δvir(z) is the non-linear density contrast of virialized structures as predicted by spherical collapse (Eke et al. 1996; Bryan & Norman 1998). The multiplicity function for halo masses computed at the virial radius has been shown to preserve universality better than other commonly assumed definition of halo radii (Despali et al. 2016; Diemer 2020; Ondaro-Mallea et al. 2021). As for δc, we use the fitting formula introduced by Kitayama & Suto (1996) that ignores the effect of massive neutrinos; however, for the adopted values for total neutrino masses in this work, the fitting formula is still percent level accurate (LoVerde 2014).

In this paper, we use the HMF presented in Euclid Collaboration (2023): vf(v)=A(p,q)2av2πeav2/2(1+1(av2))(va)q1,$v\,f(v) = A(p,q)\sqrt {{{2a{v^2}} \over \pi }} {{\rm{e}}^{ - a{v^2}/2}}\left( {1 + {1 \over {\left( {a{v^2}} \right)}}} \right){(v\sqrt a )^{q - 1}},$(4)

where the parameters {a, p, q} depend on background evolution and power spectrum shape as a=aRΩmaɀ(ɀ),$a = {a_R}{\rm{\Omega }}_{\rm{m}}^{{a_}}(),$(5) p=p1+p2(d ln σd ln R+0.5),$p = {p_1} + {p_2}\left( {{{{\rm{d}}\ln \sigma } \over {{\rm{d}}\ln R}} + 0.5} \right),$(6) q=qRΩmqɀ(ɀ),$q = {q_R}{\rm{\Omega }}_{\rm{m}}^{{q_}}(),$(7)

and where Ωm is the fractional density of matter in the Universe as a function of redshift, encompassing both baryonic and dark matter contributions aR=a1+a2(d ln σd ln R+0.6125)2,${a_R} = {a_1} + {a_2}{\left( {{{{\rm{d}}\ln \sigma } \over {{\rm{d}}\ln R}} + 0.6125} \right)^2},$(8) qR=q1+q2(d ln σd ln R+0.5).${q_R} = {q_1} + {q_2}\left( {{{{\rm{d}}\ln \sigma } \over {{\rm{d}}\ln R}} + 0.5} \right)$(9)

Lastly, the normalization parameter A is not a free parameter but a function of the other parameters: A(p,q)={ 21/2p+q/2π[ 2pΓ(q2)+Γ(p+q2) ] }1,$A(p,q) = {\left\{ {{{{2^{ - 1/2 - p + q/2}}} \over {\sqrt \pi }}\left[ {{2^p}\Gamma \left( {{q \over 2}} \right) + \Gamma \left( { - p + {q \over 2}} \right)} \right]} \right\}^{ - 1}},$(10)

where Γ denotes the Gamma function. The adopted values for the HMF parameters are presented in Table 4 of Euclid Collaboration (2023) and depend on the halo-finder used. In this work, we mostly use the ROCKSTAR calibration. The SUBFIND calibration is also used in Sect. 4.1.2 to assess the impact of the halo-finder in our model.

2.2 The linear halo bias

The overdensity of halos of mass M at the position r at redshift z, δh(r,M,ɀ)=n(r,M,ɀ)/n¯(M,ɀ)1,${\delta _{\rm{h}}}({\bf{r}},M,) = n({\bf{r}},M,)/\bar n(M,) - 1,$(11)

is expressed in terms of the corresponding local number halo number density, n(r, M, z), and of the cosmic mean number density of such halos, n¯(M,z)$\bar n(M,z)$. In linear theory, it is related to the matter density contrast δm(r, z) as δh(r,M,ɀ)=b(M,ɀ)δm(r,ɀ)+ϵ(r,M,ɀ),${\delta _{\rm{h}}}({\bf{r}},M,) = b(M,){\delta _{\rm{m}}}({\bf{r}},) + ({\bf{r}},M,),$(12)

where b(M, z) is the linear halo bias and ϵ is a stochastic term that in the following we assume to be associated with shot-noise.

It follows from Eq. (12) that the halo-halo, Ph , and halomatter power spectrum, Phm , are written as a function of the linear matter power spectrum, Pm , for sufficiently large scales as Ph(k,M,ɀ)=b2(M,ɀ)Pm(k,ɀ)+PSN,${P_{\rm{h}}}(k,M,) = {b^2}(M,){P_{\rm{m}}}(k,) + {P_{{\rm{SN}}}},$(13) Phm(k,M,ɀ)=b(M,ɀ)Pm(k,ɀ),${P_{{\rm{hm}}}}(k,M,) = b(M,){P_{\rm{m}}}(k,),$(14)

where PSN represents the shot-noise component. Under the assumption that halos offer a discrete Poisson sample of the underlying continuous matter density field, PSN denotes a shotnoise component commonly assumed to be equivalent to the Poisson term, PSN = 1/N, where N is the mean number density of tracers. On the other hand, halos are known to correspond with high-density peaks of the underlying matter distribution. Therefore, they are expected not to provide a purely Poissonian sampling of this continuous density field. In fact, Casas-Miranda et al. (2002) and Hamaus et al. (2010) showed that positive and negative corrections to Poisson shot-noise are expected for low- and high-mass halos. In this paper, we parameterize PSN as PSN=1αn¯,${P_{{\rm{SN}}}} = {{1 - \alpha } \over {\bar n}},$(15)

where α, a fitting parameter that we calibrate through simulations, controls the deviation from the assumption of Poisson noise.

Assuming the universality of the HMF, Mo & White (1996) derived the HB b(M, z) directly from the HMF by following the PBS framework. The PBS prediction for the HB as a function of the peak height ν reads bPBS(v)=11δc d ln vf(v)d ln v.${b_{{\rm{PBS}}}}(v) = 1 - {1 \over {{\delta _{\rm{c}}}}}{{{\rm{d}}\ln v\,f(v)} \over {{\rm{d}}\ln v}}.$(16)

Although the PBS provides an estimate of the bias that presents the correct qualitative behavior, Tinker et al. (2010) claims a relatively poor performance of the PBS in reproducing results from N -body simulations, with an accuracy of about 10–20%.

Given the correct qualitative behavior of the PBS prescription, in this paper, we aim to improve the prediction of the bias by calibrating a model for the bias as a function of bPBS – that is, we assume Eq. (16) to be valid, with the HMF (4) – and model its difference with respect to the simulations.

Table 1

Cosmological parameters of the PICCOLO set of simulations.

3 Methodology

3.1 Simulations

3.1.1 N-body simulations

Table 1 presents the adopted values for the matter density parameter and baryonic density parameter at redshift 0 (Ωm,0 and Ωb,0), the dimensionless Hubble parameter h, the spectral index of the primordial power spectrum ns, and the amplitude of matter density fluctuations on scales of 8h1 Mpc σ8 for the N-body simulations used in this work. We extended the set of PICCOLO simulations introduced and used by Euclid Collaboration (2023) to calibrate the HMF model. We maintain the same technical configurations as the original PICCOLO simulations and refer to the above-mentioned HMF paper for further details while summarizing the main aspects. The set comprises 69 cosmological boxes, each with a comoving size of 2000 h1 Mpc, and 4 × 12803 dark matter particles. The simulations were conducted using OpenGADGET3, with initial conditions generated by monofonIC(Michaux et al. 2021), based on third-order Lagrangian Perturbation Theory (3LPT) at a redshift of z = 24. The adopted gravitational softening is equivalent to one-fortieth of the mean inter-particle distance.

The original PICCOLO set of simulations included nine different choices for cosmological parameters, which were randomly chosen from the 95% confidence level hyper-volume of the joint SPT and DES cluster abundance constraints (Costanzi et al. 2021). Those represent the cosmologies C0 to C8 in Table 1. To further stress our modeling and guarantee its robustness, we also add the cosmologies C9 and C10, which sample the (Ωm,0, σ8) plane in the direction orthogonal to the degeneracy direction of the constraints from Costanzi et al. (2021), and in significant tension with such constraints. For each cosmology, two white-noise realizations were created to generate initial conditions. For each noise realization, a pair was generated by fixing the amplitudes of the Fourier modes of the density fluctuation field and pairing the phases (Angulo & Pontzen 2016), except for the reference C0 model, which had 20 realizations.

We further added three simulations with Einstein–de Sitter cosmology (EdS; i.e., Ωm = 1 and ΩΛ = 0) with powerlaw initial matter power spectrum Pm(k)kns${P_{\rm{m}}}(k) \propto {k^{{n_{\rm{s}}}}}$ with ns ∈ {–1.5, –2.0, –2.5}. Those simulations, which have the same box size and the same number of particles as the PICCOLO set, are only instrumental for the modeling and are not used for the calibration, as they are far from the regime used to calibrate the model of Euclid Collaboration (2023).

Lastly, we carried out three pairs of simulations with massive neutrinos, again using the same box size as the PICCOLO set. Each pair has a total neutrino mass Σmν ∈ {0.15, 0.30, 0.60} eV. The simulation set-up for the neutrino simulations is the same used for the OpenGADGET3 simulations extensively validated in Adamek et al. (2023). The baseline cosmological parameters is the C0 set, with Ων,0 subtracted from Ωm,0. The simulations have the same primordial amplitude As, resulting in a lower σ8 for increasing neutrino mass (see Table 1).

The initial conditions (ICs) for the neutrino simulations were generated using the FastDF (Elbers 2022) implementation to monofonIC (Michaux et al. 2021)1. The forked repository with the FastDF integration can be found here2. We employed the same number of neutrino particles as the number of grid resolution elements used for the cold dark matter particles. The total neutrino mass specified is attributed to a single massive neutrino species.

3.1.2 Approximate methods: PINOCCHIO

In this paper, we also analyzed 200 (100 pairs with each pair having fixed amplitudes and paired phases) halo catalogs simulated with the approximate LPT-based method implemented in the PINOCCHIO code (Monaco et al. 2002, 2013; Munari et al. 2017). All these simulations have been carried out under the assumption of the C0 cosmological parameters. The rationale for this extra set of simulated catalogs is to model the impact of fixing and pairing the Fourier mode amplitudes in the ICs on the cluster clustering.

3.2 Halo finders

Euclid Collaboration (2023) showed that the halo-finder adopted for the analysis of the N -body simulations can significantly alter the HMF. To understand if the halo-finder also impacts the HB, we selected two algorithms to extract halo catalogs: ROCKSTAR (Behroozi et al. 2013a)3 and SUBFIND (Springel et al. 2001; Dolag et al. 2009; Springel et al. 2021). Although all these algorithms rely on the SO method to define halo boundaries, they differ in the method used to identify the center from which the spheres are grown and the criteria to classify between structures and sub-structures.

ROCKSTAR divides the simulation volume into 3D friend-of- friends (FOF; see, for instance, Davis et al. 1985) groups and runs a recursive 6D FOF algorithm on each group to create a hierarchy of FOF subgroups. Halo centers are determined by averaging the positions of particles in the innermost subgroup. To improve consistency, we apply the CONSISTENT algorithm, which dynamically tracks halo progenitors, to the extracted ROCKSTAR catalogs as demonstrated in Behroozi et al. (2013b). SUBFIND also determines halo centers using a parallel implementation of the 3D FOF algorithm but directly assigns it to the particle with the lowest gravitational potential.

Among the halo-finders studied by Euclid Collaboration (2023), ROCKSTAR and SUBFIND are good representatives of the heterogeneity of possible results as they are close to the extremes, with SUBFIND suppressing the abundance of objects more massive than 1013 M h1 by roughly 10%. See Knebe et al. (2011) for a more detailed comparison between the halo-finding algorithms.

3.3 Measuring the halo bias

To measure the HB, we bin the halo distribution in log10(Mvir/M h1) with equispaced width intervals of 0.1 dex at each redshift. If the number of halos inside a bin was less than 10 000, we merged it with its neighbor to avoid having bins where the power-spectrum measurements were primarily dominated by shot noise. We measured the cross-spectrum Phm between the halos in each mass bin and the matter distribution traced by particles. The bias is then obtained as the ratio between this cross-spectrum and the matter power spectrum Pm(k). The matter density field was computed at the initial conditions and rescaled according to the linear growth factor for the simulations without massive neutrinos. For the simulations with massive neutrinos, the matter density field was built from the particle data from the same snapshot from which the halo catalog was extracted to consider the scale-dependent growth factor. We used the PYLIANS4 Python libraries to construct the density field and compute power spectra on a 10243 piecewise cubic spline mesh grid. PYLIANS averages the power spectra in k-space shells with the width given by the fundamental mode of the box, kf ≡ 2π/L. Following Castro et al. (2020), we only considered modes with k values smaller than 0.05 h Mpc–1 to measure the HB, to ensure the validity of the linear approximation. The maximum k used for the measurements corresponds to the 16th harmonic of the box and is much smaller than the Nyquist frequency of the grid used to compute the power spectrum.

To calculate the bias for each mass bin, we used the ratio of the halo-matter cross-spectrum and the matter power spectrum bi,jsim=Phm(kj)i/Pm(kj),$b_{i,j}^{{\rm{sim}}} = {P_{{\rm{hm}}}}{\left( {{k_j}} \right)_i}/{P_{\rm{m}}}\left( {{k_j}} \right),$(17)

where i and j are the mass bin and the k-shell indexes, respectively.

3.4 Halo bias calibration

We used a Bayesian approach with uninformative uniform priors on all parameters to fit our model for the linear HB parameters to simulation results. The best fits were obtained using the Dual Annealing method to find the posterior maximum as implemented by Virtanen et al. (2020), and covariance between parameters was estimated using PyMC (Salvatier et al. 2016)5. The No-U-Turn sampler (NUTS) (Hoffman et al. 2014) was automatically assigned internally by PyMC to sample the likelihood.

We assumed a Gaussian likelihood for the bias, since the power spectrum estimation for a Gaussian field realization follows a χ2-distribution when averaged over a shell. Since the number of modes Nk inside the shell increases rapidly with k, this distribution approaches a Gaussian distribution by the central limit theorem. However, the number of modes is small for the first bins, and deviation from the Gaussian approximation is evident. To avoid this issue, we re-binned the measurement of the first 3 k-bins by merging them, ensuring that the bin with the fewest modes still contains 117 modes to recover the Gaussian approximation’s validity.

We note that differently than Castro et al. (2020), we used ICs with fixed amplitudes. Therefore, the distribution of the simulated bias is not approximately a ratio of two Gaussian distributions but approximately Gaussian itself since the denominator of Eq. (17) is not a random variable. The variance of the shell-average halo-matter cross-spectrum on simulations with random Gaussian initial conditions is given by (σPhmPm)2=1Nk(b2+1αn¯Pm),${\left( {{{{\sigma _{{P_{{\rm{hm}}}}}}} \over {{P_{\rm{m}}}}}} \right)^2} = {1 \over {{N_k}}}\left( {{b^2} + {{1 - \alpha } \over {\bar n\,{P_{\rm{m}}}}}} \right),$(18)

where we assumed that the shot-noise contribution to the halo power spectrum follows Eq. (15) and a linear HB b. However, Zhang et al. (2022) showed that the predictions for random Gaussian initial conditions overestimate the variance observed for biased tracers on simulations with fixed amplitudes. Therefore, we modified Eq. (18) as follows: (σPhmPm)2σb2=1Nk(βb2+1αn¯Pm)+b2σsys2.${\left( {{{{\sigma _{{P_{{\rm{hm}}}}}}} \over {{P_{\rm{m}}}}}} \right)^2} \equiv \sigma _b^2 = {1 \over {{N_k}}}\left( {\beta \,{b^2} + {{1 - \alpha } \over {\bar n\,{P_{\rm{m}}}}}} \right) + {b^2}\sigma _{{\rm{sys}}}^2.$(19)

Here, β and σsys are parameters we marginalize over, which control the variance suppression and the relative error due to the limited accuracy of the HMF used as the backbone for the PBS prescription. We note that, based on the halo model (see, for instance Cooray & Sheth 2002), one should expect a value for β close to zero as in the limit where all halos are considered, the shot-noise term on the right-hand side of Eq. (19) tends to zero so that one should recover the matter power spectrum that by construction has zero variance. Furthermore, the HMF presented in Euclid Collaboration (2023) was shown to have percent-level accuracy; thus, σsys is expected to assume similar values during the calibration, presuming the PBS framework is valid. However, it is crucial to note that should σsys significantly deviate from zero, such an occurrence could indicate a potential violation or limitation within the PBS framework, underscoring the necessity for careful interpretation of these parameters.

We obtained the total log-likelihood by summing up all mass bins, modes, redshifts, and simulations. Following Euclid Collaboration (2023), we use the redshifts z ∊{2.00, 1.25, 0.90, 0.52, 0.29, 0.14, 0.0}, translating to a timespacing of about 1.7 Gyr. This spacing is larger than the characteristic dynamical time of galaxy clusters and effectively suppresses the correlation between the results of different snapshots.

Similarly, we assume that the correlations between different mass bins and modes are negligible. This assumption is justified by linear theory, which posits that different modes evolve independently in the linear regime, thus minimizing their mutual influence. In our analysis, we have three fitting parameters: α, β, and σsys , in addition to the parameters of the bias model that are subject to calibration (to be discussed in Sect. 4.1.3). This approach allows for a comprehensive calibration of the bias model, taking into account the shot-noise correction α, the suppression of variance β, and the systematic uncertainties σsys inherent to our method.

3.5 Forecasting Euclid’s cluster counts and cluster clustering observations

To understand the impact of the HB calibration on cosmological constraints, it is important to realistically forecast the cosmological information to be extracted from the Euclid photometric cluster survey. For this purpose, we first quantify the impact of the HB on cluster counts and cluster clustering analyses. More precisely, the HB enters the modeling of the cluster counts covariance, which we compute analytically following the model of Hu & Kravtsov (2003), as validated in Fumagalli et al. (2021). Regarding the cluster clustering, the HB enters both in the computation of the mean value (power spectrum or two-point correlation function), and in the associated covariance matrix; in this work, we test the effect on the real-space two-point correlation function and its covariance, following the model presented and validated in Euclid Collaboration (2024b).

After assessing the impact on the two statistics, we forecast how the accuracy of the HB calibration propagates on the cosmological constraints obtained by cluster counts and cluster clustering experiments. We generate synthetic cluster abundance data as described in Section 2.5 of Euclid Collaboration (2023), assuming the HMF calibrated and the HB of this work as benchmarks. Through a likelihood analysis, we constrain the cosmological parameters Ωm,0 and σ8, and the mass-observable relation (MOR) parameters Aλ, Bλ, Cλ, Dλ (see Section 2.5 of Euclid Collaboration 2024a), assuming flat priors for all the parameters. The MOR parameters describe the optical richness λ distribution as a function of the halo mass according to lnλMvir,ɀ =lnAλ+Bλln(Mvir3×1014h1M)                                   +Cλln(H(ɀ)H(ɀ=0.6)),$\matrix{ {\langle \ln \lambda \mid {M_{{\rm{vir}}}},\rangle = \ln {A_\lambda } + {B_\lambda }\ln \left( {{{{M_{{\rm{vir}}}}} \over {3 \times {{10}^{14}}{h^{ - 1}}{M_ \odot }}}} \right)} \hfill \cr {\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, + {C_\lambda }\ln \left( {{{H()} \over {H( = 0.6)}}} \right),} \hfill \cr } $(20)

where H(z) denotes the Hubble parameter at redshift z. The range for richness λ is considered to be between 20 and 2000, with the variation in logarithmic richness for a given virial mass and redshift expressed as a log-normal scatter: σlnλMvir,ɀ2=Dλ2.$\sigma _{\ln \lambda \mid {M_{{\rm{vir}},}}}^2 = D_\lambda ^2.$(21)

The reference values for the parameters are Aλ = 37.8, Bλ = 1.16, Cλ = 0.91, and Dλ = 0.15. These values have been obtained refitting the parameters presented in Saro et al. (2015) to the virial mass definition, with the assumption that halo profiles follow an NFW profile with a mass-concentration relationship as outlined in Diemer & Joyce (2019). The parameter values adopted in this work are consistent with the model presented by Castignani & Benoist (2016) to assign cluster membership probabilities to galaxies from photometric surveys.

We perform the analysis on the synthetic catalogs, comparing them with the predictions made by using our bias and the one of Tinker et al. (2010), and compare the resulting posteriors with two estimators: we quantify the broadening/tightening of the posterior’s amplitude by computing the difference of the figure of merit (∆FOM; see Huterer & Turner 2001; Albrecht et al. 2006) in the Ωm,0σ8 plane, and the shift of the posterior’s position by computing the posterior agreement (Bocquet et al. 2019), which determines whether the difference between two posterior distributions is consistent with zero.

thumbnail Fig. 1

Constraints on the parameters in Eq. (19) fitted to the unbiased standard deviation of the 200 PINOCCHIO mocks with C0 cosmology. We fitted α and β considering k ≤ 0.05 h Mpc−1, assuming a Gaussian likelihood with error bar estimated from the measurements for 0.05 ≤ k/h Mpc−1 ≤ 0.2 and assuming it to be constant in k and equivalent to the unbiased standard deviation.

4 Results

4.1 Calibration of the halo bias

4.1.1 Biased tracers statistics on fixed and paired simulations

Before modeling and calibrating the HB to the simulations, we investigate the impact of the variance suppression technique on the halo matter cross-spectrum.

In Fig.1, we present the constraints on the parameters in Eq. (19) fitted to the unbiased standard deviation of the 200 PINOCCHIO mocks with C0 cosmology. We fixed σsys to zero, as this exercise does not involve modeling errors on the bias. The parameters α and β were fitted using only modes with k ≤ 0.05 h Mpc−1, under the assumption of a Gaussian likelihood. We estimated the error bars from the measurements within the range 0.05 ≤ k/h Mpc−1 ≤ 0.2, treating it as constant and equivalent to the unbiased standard deviation. This approach prevents the estimation of the expectation of the mean and the error from using the same data. Figure 1 reveals a positive correlation between the parameters, with both assuming positive values. As expected, the posterior for the β parameter peaks close to zero, validating the effectiveness of the fixed amplitudes technique in reducing the variance of biased tracer statistics. Notably, the small value of α found in our analysis suggests that the shot noise is well modeled as Poissonian to within 1–3%.

In Fig. 2, we present the relative difference between Eq. (19) computed with the best-fit values of α and β, and the unbiased standard deviation of the measurements, for different mass bins and redshift values. We present the results for three values z ∊ 0.0, 1.0, 2.0 spanning the redshift range of interest, while for the masses, we present the first, the intermediate, and the last occupied bin for each redshift. We observe that the residuals of the fit always oscillate around zero, with no statistically significant mass, redshift, or k dependence.

In Fig. 3, we present the Pearson correlation coefficient ρ between the measurements of the halo-matter cross-power spectrum on a given simulation and its paired realization. The correlation coefficient between two random variables X and Y is defined as ρ(X,Y)=(XX)(YY)σXσY,$\rho (X,Y) = {{\langle (X - \langle X\rangle )(Y - \langle Y\rangle )\rangle } \over {{\sigma _X}{\sigma _Y}}},$(22)

where σX and σY are the standard deviation of the random variables X and Y. The cross-power spectrum was measured for k ≤ 0.05 h Mpc−1, for different ranges of halo masses (as reported in each panel) and redshift values (different columns). We also present the correlation coefficient between simulations that assumes uncorrelated white-noise realizations for comparison. We note that paired simulations do not show a statistically significant difference in their correlation with respect to simulations that assume uncorrelated noise realizations. The same conclusion is obtained by running a p-value test on all mass and redshift bins. This result justifies the assumption that different simulations are, in fact, independent, even if they have the same amplitudes and paired phases. This conclusion aligns with the claims of Villaescusa-Navarro et al. (2018) that variance suppression techniques are unlikely to affect the halo abundance distribution. As the shot-noise term in Eq. (19) dominates over the other terms for our sample selection, one could anticipate the independence of the bias results from fixed and paired simulations due to the independence of abundance fluctuations and modes paring.

Lastly, we use the PINOCCHIO catalogs to assess the impact of neglecting the correlation between different mass bins and Fourier modes on the calibration likelihood presented in Sect. 3.4. We measured the bias on the PINOCCHIO catalogs by applying the same mass and modes binning we used in our calibration and measured the correlation matrix between different simulations. In line with the results presented in Euclid Collaboration (2024b) for the two-point correlation function, we explicitly verified in our analysis that the off-diagonal terms are sub-dominant and of the order of 10%, validating our calibration likelihood.

4.1.2 Impact of the halo finder on the peak-background split performance

In Fig. 4, the impact of the halo finder on the PBS is examined through the bias ratio measured in halo catalogs generated from the same simulation using either ROCKSTAR or SUBFIND, alongside the corresponding PBS prescription. For the ROCKSTAR catalogs, the standard error of the mean is further illustrated using an additional 19 realizations, with the assumed cosmology being C0. The analysis spans different redshifts within the z ∊ [0, 2] range. It is observed that the PBS tends to underestimate the simulation-derived bias by approximately 10% at higher redshifts, though it shows improved accuracy at z = 0. Given the minimal impact of the choice of halo finder on PBS’s overall accuracy, subsequent results focus exclusively on the ROCKSTAR halo catalogs. Despite PBS’s tendency to underpredict the HB compared to simulations, it is noteworthy that the deviation remains consistently within 5–15% across all explored values of v and redshifts. Our future efforts will aim to refine the PBS model by addressing these discrepancies, with the objective of achieving a simulation-calibrated HB model that is precise to within a few percent.

thumbnail Fig. 2

Relative difference between Eq. (19) best fit and the unbiased standard deviation of the PINOCCHIO measurements. Different columns are for different redshifts, and the corresponding mass bins are shown in each panel. The vertical dotted line demarcates two distinct sets of measurements: Those to the left of the line were utilized as data points in the parameter fitting process, while the scatter of the points to the right was analyzed to estimate the variance. The gray regions highlight areas within a 5% deviation from the expected values.

4.1.3 Modeling the halo bias

In Fig. 5, we present the mean of the ratio of the measured bias with respect to the PBS prescription for different simulations. In the left panel, we use all the CO runs and show the ratio of the bias as a function of v/(1 + z) for redshifts 0, 0.29, and 1.25. The factor (1 + z) is only used to scale the results from different red-shifts of the C0 model to the same range. In the middle panel, we show the mean ratio as a function of v for the three EdS cosmologies. Lastly, in the right panel, we present the mean ratio as a function of the background evolution Ωm(z) for the C0, C9, and C10 cosmologies.

The left panel of Fig. 5 illustrates that the performance of PBS is influenced by the cosmological background evolution, encapsulated by Ωm(z), yet appears unaffected by variations in v. This contrasts with the observations in the central panel, where PBS performance varies with both v and changes in the power spectrum’s shape. This sensitivity to v is attributed to the limitations of the HMF calibration by Euclid Collaboration (2023) when applied to an EdS cosmology that is far from its calibration regime, introducing artificial dependency.

However, it is important to note that, although these extrapolations to EdS scenarios are beyond the initial calibration range of the HMF model, the model’s accuracy is not disproportionately affected across different values of ns. Indeed, in Euclid Collaboration (2023), it has been demonstrated that the HMF model retains a consistent level of precision across various EdS cosmologies characterized by scale-free linear power spectra. Therefore, the dependence on the shape of the power spectrum is more likely related to the varying degree of accuracy of the PBS bias model as the shape of the power spectrum changes.

We interpret the dependence of the PBS performance as a function of the shape of the power spectrum as follows. The extrapolation of the results on the central panel of Fig. 5 indicates that the PBS performance on EdS cosmologies with a steeper power spectrum (ns < −2.5) degrades with decreasing ns. Within the PBS framework, the mass variance σ(RL) smoothed on the scale of the Lagrangian patch RL is assumed to be dominated by the contribution of scales RLSSRL For a power-law power spectrum, it is d ln σ d ln R=(n+3)2.${{{\rm{d}}\ln \sigma } \over {{\rm{d}}\ln R}} = - {{(n + 3)} \over 2}.$(23)

Therefore, the ratio σ (RL) / σ(RLSS) tends to unity as ns tends to −3. On the one hand, this explains why ns = −2.5 presents better performance than the other cases as one of the PBS assumptions is better suited. On the other hand, for ns ≡ −3, perturbations on all scales reach the collapse at the same time, and it is no longer possible to distinguish between peaks and a large-scale modulation of a background perturbation, thus breaking PBS’s fundamental assumption.

The right panel of Fig. 5 shows that the residuals with respect to the PBS prediction increase linearly with the value of the density parameter Ωm(z). While the slope of this linear dependence is similar for the three cosmologies, the normalization is a decreasing function of the clustering amplitude S8. In fact, C9 is the simulation with the lowest S8=(Ωm,0/0.3  σ8)=0.438${S_8} = \left( {\sqrt {{{\rm{\Omega }}_{{\rm{m}},0}}/0.3} \,\,{\sigma _8}} \right) = 0.438$ while C10 has the highest clustering amplitude S8 = 1.07. The C1 to C8 simulations are not shown in this panel for better readability, but they cluster around C0 as they have similar S8 values.

The better performance of the PBS in cosmologies with more clustering suggests that the difference between this model and the simulation results is related to the connection between Lagrangian patches in the initial density field and the collapsed structures identified by the halo-finder. Collapsed structures stand out more clearly from the non-linear density field in more evolved and clustered cosmologies. For less clustered models, large halos are still forming and overlapping due to ongoing mergers. This makes it more challenging to identify and link them to their corresponding Lagrangian patches clearly. Not surprisingly, in the EdS cosmologies, the PBS best performance is for ns = −2.5, where the evolution of the power spectrum is the steepest.

The above line of reasoning suggests then that an SO algorithm could not be accurate in providing a one-to-one mapping between Lagrangian patches, destined to form virialized halos according to spherical collapse, and for which the PBS method predicts the bias, and halos identified in the non-linearly evolved density field. In this vein, since both ROCKSTAR and SUBFIND are based on an SO algorithm, it is not surprising that they predict similar deviations from PBS (see Fig. 4).

On the other hand, we expect that collapsed structures have had more time to relax in cosmological models characterized by a higher value of S8. As a consequence, they are more likely to be spherical. Again, this is in line with the better performance of the PBS on evolved cosmologies, as shown in the right panel of Fig. 5. Although suggestive, this interpretation of the deviations of PBS predictions would require a dedicated analysis to track their origin in detail, which goes beyond the scope of the analysis presented here.

From the results shown in Fig. 5, it emerges that deviations from the PBS should depend on cosmic evolution, parameterized by Ωm(z), on the slope of the linear power spectrum, and on the clustering amplitude S86. To capture such dependencies, we adopted the following description of the correction to the PBS prediction for the linear HB: bbPBS:=f(Ωm(ɀ),d ln σd ln R,S8)           =A0f0(Ωm(ɀ))f1(d ln σd ln R)f2(S8).$\matrix{ {{b \over {{b_{{\rm{PBS}}}}}}: = f\left( {{{\rm{\Omega }}_{\rm{m}}}(),{{{\rm{d}}\ln \sigma } \over {{\rm{d}}\ln R}},{S_8}} \right)} \hfill \cr {\,\,\,\,\,\,\,\,\,\,\, = {A_0}{f_0}\left( {{{\rm{\Omega }}_{\rm{m}}}()} \right){f_1}\left( {{{{\rm{d}}\ln \sigma } \over {{\rm{d}}\ln R}}} \right){f_2}\left( {{S_8}} \right).} \hfill \cr } $(24)

In the above expression, we assumed the following functional forms for the three above dependencies: f0(x)=1+a1x,${f_0}(x) = 1 + {a_1}x,$(25) f1(x)=1+b1x+b2x2,${f_1}(x) = 1 + {b_1}x + {b_2}{x^2},$(26) f2(x)=1+c1x.${f_2}(x) = 1 + {c_1}x.$(27)

The parameters A0, a1, b1, b2, and c1 are calibrated in the next section through a detailed comparison with simulations. A balance between simplicity and empirical accuracy drives the parameterization chosen for these contributions. Specifically, we opted for a linear relationship for Ωm(z) and S8 while modeling the shape of the power spectrum using a quadratic function. In order to assess the potential parameter redundancy in using an extra parameter for the shape of the power spectrum, we performed Watanabe-Akaike Information Criterion (WAIC) (Watanabe 2010) and Pareto Smoothed Importance Sampling Leave-One-Out Cross-Validation (PSIS-LOO) (Vehtari et al. 2017) analyses comparing the model with b2 free to vary with respect to b2 fixed to 0. Both analyses confirmed that the fewer degrees of freedom adopted by the surrogate model do not compensate for the decrease in the model’s predictability.

Lastly, we report that we do not observe any significant correlation between the model prediction residual and other cosmological parameters assumed in the simulations. Thus, we conclude that the three fi components in Eq. (24) used in our analysis are sufficient to achieve our goal.

thumbnail Fig. 3

Correlation coefficient ρ between the measurements of the halo-matter power spectrum, Phm(k), in different simulations for k ≤ 0.05 h Mpc−1 as a function of halo mass bin and redshift. The red histograms show the distribution of the correlation coefficient between a simulation and its paired realization, while blue histograms are for the correlation coefficient between simulations with uncorrelated white-noise realizations.

thumbnail Fig. 4

Relative difference between the bias measured in halo catalogs and the bias predicted by the PBS model of Eq. (16) at different red-shifts in the range z ∊ [0, 2], Results refer to simulations carried out for the C0 cosmology. In each panel, blue and red lines refer to the results obtained for the halo catalogs based on the application of ROCKSTAR and SUBFIND, respectively. For the ROCKSTAR catalogs, we also show the standard error of the mean using other 19 realizations of the same cosmology.

4.2 Calibration of the halo bias correction

In Fig. 6, we present the marginalized two-dimensional and unidimensional constraints on the model parameters presented in Eqs. (24) to (27). The best-fit and 95% limits are reported in Table 2. We calibrate our model using the subset of 60 PICCOLO simulations covering the cosmological parameters from C0 to C10.

In Fig. 7, we present the ratio between the mean of the observations on the 20 C0 simulations with respect to our model predictions. Different rows correspond to different redshifts, while each panel corresponds to a different mass bin. The presented mass bins were selected as before to span from the least massive to the most massive occupied bin in that redshift. The shaded region in red corresponds to the error on the mean, assuming that each measurement follows Eq. (19). The shaded regions in gray correspond to 2% and 4% regions. As can be seen, our model’s prediction presents a performance below 2% for different masses and redshift regimes when not primarily dominated by the sample variance, as it happens, for instance, for k ≲ 10−2 h Mpc−1. This accuracy holds over the k range up to which the onset of non-linearity occurs and the approximation of scale-independent bias breaks down, that is, k ≳ 0.05 h Mpc−1 (marked by a vertical line).

We note that at large k values, non-linearity effects cause the bias measured from the simulations to take a scale dependence in exceeding the model prediction. As expected, this effect is smaller at higher redshift, consistent with the effect of non-linearity shifting to larger k values.

Similarly to Fig. 7, in Fig. 8, we present the ratio between the mean of the observations on the C9 and C10 simulations with respect to our model predictions. C9 and C10 are the simulations with lowest and highest S8, respectively. Even for such extreme scenarios, our model performs well, thus confirming that our linear bias model, with the previously described calibration, can reproduce results from simulations with an accuracy of a few percent for ΛCDM cosmologies.

thumbnail Fig. 5

Mean of the ratio of the bias measured in simulations with respect to the PBS prescription. Left: ratio of the bias as a function of v/(1 + z) for different z, labeled by Ωm(z), for all the CO runs. Center, mean ratio as a function of v for the three EdS cosmologies of pure power-law shapes of the linear power spectrum at z = 0. We report the value of the three spectral indexes in the inset. Right: mean ratio as a function of the background evolution Ωm(z) for the C0, C9, and C10 cosmologies with varying S8.

thumbnail Fig. 6

Marginalized 68% and 95% confidence level contours on the model parameters presented in Eqs. (24) to (27). We calibrated our model using the subset of PICCOLO simulations C0–C10. (See Table 2 for the best fit and confidence levels.)

Table 2

Best fit and 95% limits on the model parameters presented in Eqs. (24)(27).

thumbnail Fig. 7

Ratio between the mean of the observations on the 20 C0 simulations with respect to our model predictions. Different rows correspond to different redshifts, while each panel corresponds to a different mass bin. The shaded region in red corresponds to the error on the mean, assuming that each measurement follows Eq. (19). The shaded regions in gray correspond to 2% and 4% regions.

4.3 Cosmologies with massive neutrinos

We present our model’s performance when considering simulations with massive neutrinos in Fig. 9. This allowed us to assess the performance of our HB calibration for this minimal extension of ΛCDM. In this case, the simulation’s bias has been computed by comparing it to the linear power spectrum of the corresponding cosmological model, which includes only the contributions from cold dark matter and baryons. For consistency, the same choice of considering only CDM and baryon contribution is also made for the computation of the HMF entering in our model for the HB (see, for instance, Castorina et al. 2014; Costanzi et al. 2013). From Fig. 9, it is evident that our model also precisely describes the bias in Λ(v)CDM models, despite the fact that such models were not used during the HB calibration.

We note that, unlike for the pure ACDM models, in this case, the measured bias underpredicts the model bias at large k. We also note that the dependence of this effect on redshift, if any, goes in the direction of being larger at higher z. Also, there is some hint for it to be slightly smaller for smaller values of mv and therefore of Ωv. These effects align with the expectation that such deviations are not dominated by non-linear evolution but rather by the effect of neutrino-free streaming (Castorina et al. 2014).

thumbnail Fig. 8

Similar to Fig. 7 but for the C9 and C10 cosmological parameters. Among the PICCOLO set, C9 and C10 correspond to the cosmologies with the lowest and the highest S8, respectively.

thumbnail Fig. 9

Similar to Fig. 7 but for simulations with massive neutrinos. For better plot readability, we only show the uncertainties (red shaded regions) for the simulation with total neutrino masses equal to 0.15 eV.

thumbnail Fig. 10

Comparison between the HB predicted by our model with predictions from other models presented in the literature: Cole & Kaiser (1989). Sheth et al. (2001), Tinker et al. (2010), and Comparat et al. (2017). We present both our benchmark model as well as the PBS predictions based on the HMF model of Euclid Collaboration (2023) used as a baseline of our model. Different columns correspond to different redshifts. The relative difference with respect to our benchmark model is presented in the panels in the second row. We adopted a composite scale for the residual plot to show the dynamic range of differences between the models: The scale is linear for values between [−10, 10] % and symmetric log outside. For reference, we show the zero line in black. The predictions of the models from the literature have been computed using the COLOSSUS toolkit (Diemer 2018).

4.4 Comparison with previous models

In Fig. 10, we compare our model prediction with other models in the literature: Cole & Kaiser (1989), Sheth et al. (2001), Tinker et al. (2010), and Comparat et al. (2017). We present both our benchmark model as well as the PBS predictions based on the HMF model of Euclid Collaboration (2023), which we use as a baseline of our model. Different columns correspond to different redshifts. The relative difference to our benchmark model is presented in the panels in the second row. The predictions of the external models have been computed using the COLOSSUS toolkit (Diemer 2018)7.

To ensure a fair comparison, we adopted the Planck-like CO cosmology, where the compared models have their peak performance among the PICCOLO cosmologies. All compared models have degraded performance as we move from this benchmark cosmology while we have shown the robustness of our model in Figs. 7, 8, and 9. That is due to the models assuming either the universality of the bias relation or a redshift-only dependence, while our method explicitly models the cosmology dependence of this relation. As such, comparisons between these models should be interpreted cautiously, as the underlying cosmology influences the exact figures.

As for the comparison with the PBS prediction based on the HMF calibration by Euclid Collaboration (2023), the results shown here confirm those shown in Fig. 5: the PBS-based predictions underestimate our simulations-based calibration by about 5–10%, almost independently of v. The models of Cole & Kaiser (1989) and Sheth et al. (2001) over- and under-estimate the bias by ~10%, respectively. Their relatively poor performance is not surprising. The prediction by Cole & Kaiser (1989) corresponds to the PBS prediction when using the HMF from Press & Schechter (1974). As the Press & Schechter (1974) HMF only qualitatively explains the abundance of halos, it is expected that the bias from the PBS will not perform much better. On the other hand, the Sheth et al. (2001) model was calibrated on simulations. However, such simulations had a resolution that allowed those authors to cover a dynamic range significantly narrower than that accessible to our simulations. The prediction by Cole & Kaiser (1989) differs from ours by an amount almost independent of y and redshift. On the other hand, the HB by Sheth et al. (2001) differs from ours in a v- and z-dependent way.

Notably, the model of Tinker et al. (2010) is only superior to the abovementioned models for low peak-height. The differences with respect to our model grow with redshift and peak-height. This could be due to the heterogeneity of the simulations used to calibrate the model of Tinker et al. (2010) and the possible limitation of the model itself to capture the cosmological dependence of the HB. In this paper, we calibrate to a set of simulations that have been run with the same code and setup. On the other hand, Tinker et al. (2010) based their analysis on a collection of simulations carried out with different codes and configurations in terms of resolution and box sizes. Also, their model assumes a redshift dependency for the evolution, while from Fig. 6 we note that a parametrization of this evolution through Ωm(z) is more universal.

Lastly, the model by Comparat et al. (2017) shows a good agreement with our model. The most significant differences are at low-redshift where the model of Comparat et al. (2017) predicts a bias that is smaller than ours by about 5%. This difference reduces to a sub-percent at high redshift. The agreement is unsurprising as the model of Comparat et al. (2017) was also calibrated on ROCKSTAR catalogs.

thumbnail Fig. 11

Percentage residuals of cluster counts (left panel) and cluster clustering covariance matrices (central and right panels), computed with the bias from Tinker et al. (2010) in comparison to the one calculated using the bias calibrated in this study (Eq. (24)). We show the full covariance matrix for number counts in mass and redshift bins. For the two-point correlation function of galaxy clusters, we show two blocks of the full covariance (low and high redshift bins) as a function of the radial separation.

4.5 Impact on cluster cosmology analysis

In this section, we forecast the impact of the HB calibration on cosmological analysis of cluster counts and cluster clustering from Euclid. We present the results for the bias model of Tinker et al. (2010) and the model calibrated in this study. The rationale for assuming Tinker et al. (2010) is to use a widely used model in cluster cosmology representative of the difference in the bias models presented in Fig. 10. Nonetheless, we do not expect the results to change significantly if we assumed the model of Comparat et al. (2017) that presents a better concordance with our model at high redshift but worse at low redshift. Therefore, the overall impact on cosmological constraints will compensate partially as the clustering cosmological signal for Euclid peaks at lower redshifts.

As described in Sect. 3.5, we assess the effect of the HB calibration in a more realistic scenario, performing a likelihood analysis of both the individual analysis of number counts and cluster clustering statistics and the combination of these probes. In all scenarios, the observable-mass relation (Eq. (20)) is calibrated by combining the probes with weak lensing (WL) mass estimates, assuming a constant error of 1%. The mass calibration is the primary source of systematic uncertainty in cluster cosmology studies, and a 1% calibration represents the goal for stage IV surveys. Therefore, the chosen setting offers a forecast of the utmost cosmological bias resulting from inaccurate modeling of the HB. Lastly, we assume three independent Gaussian likelihoods (Fumagalli et al. 2024) for number counts, clustering, and WL masses.

In the left panel of Fig. 11, we start by presenting the percentage residuals of the number counts covariance. We show the full covariance matrix for the number counts analysis, with the mass dependence within each redshift bin. Notably, the impact of the HB model is minimal at low redshift but becomes significant, reaching up to 20%, at higher redshifts, especially in the high-mass bins. However, the impact of a different bias calibration is mitigated by the shot-noise contribution when the latter becomes dominant along the diagonal, as the HB only plays a role in the computation of sample variance. To quantify the impact of such a discrepancy on parameter posteriors, we perform the likelihood analysis for a number counts experiment, as described in Sect. 3.5. From the comparison of the two posteriors, we obtain ΔFOM = −0.67 and a perfect agreement between the positions of the two contours, meaning that the impact of the HB calibration is below other systematics.

As for the analysis of cluster counts, in the central and right panels of Fig. 11, we present the percent residuals for the clustering covariance as a function of the radial separation in both a low-redshift (central panel) and high redshift (right panel) bin. Similar to our findings for the number counts, the most significant impact is observed on the off-diagonal elements, particularly at high redshift. However, in contrast to number counts covariance, the inclusion of shot-noise in cluster clustering contributes to the off-diagonal elements, helping mitigate the effect of different HB calibrations across all scales. This results in a difference that never exceeds 10%. Importantly, in the case of cluster clustering, the HB also affects the expected signal – either the two-point correlation function or power spectrum – leading to a difference independent of the radial separation but increasing with redshift, reaching a 10% level in the high-redshift interval.

The cosmological forecasts from the clustering experiment show that the minimal variation in covariance terms produces a negligible difference in the posterior amplitude, equal to ΔFOM = 1.03. However, comparing the two posteriors reveals an agreement at only 0.68 σ. This implies that the differences in the two-point correlation function translate into a sizeable shift in the cosmological constraints. Notably, the difference in the posteriors induced by the different calibrations of the HB alone surpasses the 0.25 σ threshold commonly employed in other studies (Deshpande et al. 2024) to flag systematic errors that, if exceeded, could accumulate and lead to a collectively significant difference.

The combined analysis (cluster counts + cluster clustering) posteriors are shown in Fig. 12. We notice that assuming the HB calibration by Tinker et al. (2010) still causes a shift in the posteriors with respect to the HB calibration presented in this paper. This aligns with the forecast results for the cluster clustering analysis presented above. Although the difference is reduced to 0.39 σ, the combination with number counts and weak lensing masses cannot compensate for the impact of the HB calibration that affects mostly the cluster clustering.

thumbnail Fig. 12

Parameter posteriors at 68% and 95% confidence levels obtained by analyzing number counts, weak lensing masses, and cluster clustering computed with the HB calibrated in this work (cyan contours) and the bias from Tinker et al. (2010) (orange contours). The error associated with weak lensing mass is set at 1% of the mass.

5 Conclusions

This paper presents a calibrated semi-analytical model for the HB in view of the joint cosmological exploitation of number counts and clustering of galaxy clusters from the Euclid survey. Our approach began with the PBS model, based on the HMF of Euclid Collaboration (2023), and we extended it by introducing a novel parametric correction. This correction was designed to align the PBS prediction with the results from an extended and homogenous set of N-body simulations that we carried out for vanilla ACDM models, varying cosmological parameters, and Λ(v)CDM models by varying the sum of the neutrino masses. The simulations employed fixed and paired initial conditions (see, Angulo & Pontzen 2016), providing a robust, reduced variance framework for our analysis. We measured the HB by examining the ratio of the matter-halo cross-spectrum. Additionally, we modeled the covariance of these measurements using 200 mock catalogs of the Euclid cluster survey, based on the approximate LPT-based PINOCCHIO code. This ensured a thorough understanding of the uncertainties involved in our calibration of the HB. The key findings and implications of our study are summarized in the following paragraphs.

The use of fixed and paired initial conditions for the simulations analyzed in our study proved highly advantageous for estimating the bias of tracers. By parametrizing the covariance of the bias measurements with two parameters – one controlling the shot-noise contribution and the other for suppression due to fixing, respectively α and β in Eq. (19) – we observed significant effectiveness in the variance suppression term. This was demonstrated in the constraints on the terms describing the variance in the halo-matter power spectrum, Phm(k), shown in Fig. 1. Furthermore, our analysis of the measurements of Phm(k) between paired simulations, as illustrated in Fig. 3, revealed no significant correlation between them. This finding underscores the efficacy of the fixed and paired simulation approach in providing reliable estimates of the bias factor characterizing the distribution of tracers (i.e., halos), which is free from the influences of inherent correlations that could affect the results.

The impact of the choice of the halo finder used in the analysis of the N-body simulations on the performance of the PBS is shown in Fig. 4. While comparing the ROCKSTAR and SUBFIND halo finders, we observed that the PBS generally underestimates the measured bias from simulations, particularly at higher red-shifts. However, the impact of the halo finder choice on the PBS prescription’s performance is almost negligible, thus reinforcing the robustness of our approach in assessing the PBS performance across different redshift ranges.

Our modeling of the HB, as illustrated in Fig. 5, reveals significant insights into the cosmological dependency of the performance of the PBS. The background cosmological evolution influences the PBS performance more than the peak height parameter v. This was particularly evident in different cosmologies, where PBS’s effectiveness varied with the shape of the power spectrum and the degree of clustering evolution, as described by the S8 parameter. Notably, in more clustered cosmologies, PBS improved its performance. This suggests a possible link between the ease of identifying collapsed structures in cosmologies with more evolved clustering and their corresponding Lagrangian patches. This result led us to develop a refined model for the PBS correction, expressed in Eq. (24), which incorporates terms depending on Ωm(z), the local slope of the power spectrum, and S 8.

The calibration of our model parameters, with the best fit presented in Fig. 6 and Table 2, demonstrated its robust performance across a range of cosmological conditions. Figures 7 and 8 illustrate our model prediction performance on the reference C0 simulations and on the C9 and C10 simulations. The accuracy of our model is particularly noteworthy, as it always remains below a 2% deviation for different masses and redshift regimes, with the possible exception of unrealistic cases largely influenced by sample variance. Quite remarkably, this level of precision is maintained even in extreme scenarios represented by the C9 and C10 simulations, which have the lowest and highest S8 values, respectively.

The robustness of our HB calibration is further demonstrated in scenarios involving massive neutrinos, as showcased in Fig. 9. Despite not incorporating massive neutrino simulations during the calibration phase, our model accurately predicts the HB in these cosmologies. Neutrinos are treated according to the model presented by Castorina et al. (2014, see also Costanzi et al. 2013) and the measurements of the bias with respect to the matter power spectrum of cold dark matter and baryons, as was done for the HMF in simulations with massive neutrinos in Euclid Collaboration (2023). The ability of our model to adapt and perform reliably in such scenarios without the need for recalibration highlights its robustness and versatility.

As for the comparison with HB models already introduced in the literature (see Figure 10), the models by Cole & Kaiser (1989) and Sheth et al. (2001) show significant deviations from our results, likely due to their calibration in simulations covering narrower dynamic ranges and a variety of cosmological models. The HB model by Tinker et al. (2010) shows increasing discrepancies with our model at higher redshifts and peak heights. Such differences could be attributed to its calibration on a heterogeneous set of simulations and inadequately accounting for the cosmological dependence of the HB. In contrast, the model by Comparat et al. (2017) aligns more closely with our findings, particularly at higher redshifts. This agreement is expected, as their model was also calibrated using ROCKSTAR catalogs.

As for the impact of changing the calibration of the HB on the derived cosmological posteriors, we showed in Fig. 11 the differences of the covariance matrices for a Euclid cluster count and cluster clustering analysis using both our calibration and the one provided by Tinker et al. (2010). While the impact on number counts covariance is minimal at low redshifts, it becomes substantial, up to 20%, at higher redshifts. However, the presence of shot-noise in the analysis helps mitigate this effect. In cluster clustering, we observed that the HB calibration could lead to differences in the two-point correlation function, particularly at high redshifts. This difference can potentially bias cosmological constraints beyond the 0.25 σ threshold commonly used to flag significant systematic errors (Adamek et al. 2023). Moreover, the combined analysis of number counts, cluster clustering, and weak lensing masses demonstrates that even with these additional data, the calibration of HB cannot be entirely compensated for. This highlights the importance of precise HB calibration in cluster cosmology, especially for a survey, such as the one being provided by Euclid, which will reach an unprecedented sensitivity and level of statistics.

In summary, the analysis presented in this paper has systematically calibrated and tested the HB for a range of cosmological scenarios, demonstrating its critical impact on the precision of cosmological analysis based on galaxy clusters for the Euclid mission. The resilience of our HB model against variations of cosmological models, including the presence of massive neutrinos and different degrees of clustering amplitude, highlights its robustness and adaptability. Importantly, our model is robust against the halo finder definition, inheriting its dependence through the HMF only. This is a remarkable feature, as the correspondence between halos in N-body simulations and real clusters in surveys remains a complex issue, with uncertainties in halo identification and characterization potentially influencing the extraction of cosmological parameters. Future research should focus on understanding and quantifying these uncertainties, especially concerning observational challenges, such as projection effects and the mass-observable relation. As we move forward, extending this precision to departures from the standard Λ(v)CDM framework will be crucial to fully harnessing the capabilities of next-generation cosmological surveys.

Data availability

In Castro & Fumagalli (2024), we implement the model presented in this paper, together with the models for the HMF presented in Euclid Collaboration (2023) and for the impact of baryonic feedback on cluster masses presented in Euclid Collaboration (2024a). The source code can be accessed in https://github.com/TiagoBsCastro/CCToolkit

Acknowledgements

It is a pleasure to thank Valerio Marra for constructive comments during the production of this work, Fabio Pitari and Caterina Caravita for support with the CINECA environment, Peter Berhoozi for the support with ROCKSTAR, Oliver Hahn for the support with monofonIC, and Luca Di Mascolo for the support with PyMC. TC is supported by the Agenzia Spaziale Italiana (ASI) under - Euclid-FASE D Attivita’ scientifica per la missione – Accordo attuativo ASI-INAF no. 2018-23-HH.0. SB, TC, PM, and AS are supported by the PRIN 2022 PNRR project “Space-based cosmology with Euclid: the role of High-Performance Computing” (code no. P202259YAF), by the Italian Research Center on High-Performance Computing Big Data and Quantum Computing (ICSC), project funded by European Union – NextGenerationEU – and National Recovery and Resilience Plan (NRRP) – Mission 4 Component 2, within the activities of Spoke 3, Astrophysics and Cosmos Observations, and by the INFN INDARK PD51 grant. TC and AS are also supported by the FARE MIUR grant ‘ClustersXEuclid’ R165SBKTMA. AS is also supported by the ERC ‘ClustersXCosmo’ grant agreement 716762. MC and TC are supported by the PRIN 2022 project EMC2 – Euclid Mission Cluster Cosmology: unlock the full cosmological utility of the Euclid photometric cluster catalog (code no. J53D23001620006). KD acknowledges support by the DFG (EXC-2094 – 390783311) as well as support through the COMPLEX project from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program grant agreement ERC-2019-AdG 882679. We acknowledge the computing centre of CINECA and INAF, under the coordination of the “Accordo Quadro (MoU) per lo svolgimento di attività congiunta di ricerca Nuove frontiere in Astrofisica: HPC e Data Exploration di nuova generazione”, for the availability of computing resources and support. We acknowledge the use of the HOTCAT computing infrastructure of the Astronomical Observatory of Trieste – National Institute for Astrophysics (INAF, Italy) (see Bertocco et al. 2020; Taffoni et al. 2020). The Euclid Consortium acknowledges the European Space Agency and a number of agencies and institutes that have supported the development of Euclid, in particular the Agenzia Spaziale Italiana, the Austrian Forschungsförderungsgesellschaft funded through BMK, the Belgian Science Policy, the Canadian Euclid Consortium, the Deutsches Zentrum für Luft- und Raumfahrt, the DTU Space and the Niels Bohr Institute in Denmark, the French Centre National d’Etudes Spatiales, the Fundação para a Ciência e a Tecnologia, the Hungarian Academy of Sciences, the Ministerio de Ciencia, Innovación y Universidades, the National Aeronautics and Space Administration, the National Astronomical Observatory of Japan, the Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the Research Council of Finland, the Romanian Space Agency, the State Secretariat for Education, Research, and Innovation (SERI) at the Swiss Space Office (SSO), and the United Kingdom Space Agency. A complete and detailed list is available on the Euclid web site (www.euclid-ec.org).

References

  1. Abbott, T. M. C., Aguena, M., Alarcon, A., et al. 2020, PRD, 102, 023509 [NASA ADS] [CrossRef] [Google Scholar]
  2. Adamek, J., Angulo, R. E., Arnold, C., et al. 2023, JCAP, 06, 035 [NASA ADS] [Google Scholar]
  3. Albrecht, A., Bernstein, G., Cahn, R., et al. 2006, arXiv e-prints [arXiv:astro-ph/0609591] [Google Scholar]
  4. Allen, S. W., Evrard, A. E., & Mantz, A. B. 2011, ARA&A, 49, 409 [Google Scholar]
  5. Angulo, R. E., & Pontzen, A. 2016, MNRAS, 462, L1 [NASA ADS] [CrossRef] [Google Scholar]
  6. Aricò, G., Angulo, R. E., Contreras, S., et al. 2021, MNRAS, 506, 4070 [CrossRef] [Google Scholar]
  7. Artis, E., Melin, J.-B., Bartlett, J. G., & Murray, C. 2021, A&A, 649, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Baxter, E. J., Rozo, E., Jain, B., Rykoff, E., & Wechsler, R. H. 2016, MNRAS, 463, 205 [NASA ADS] [CrossRef] [Google Scholar]
  9. Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013a, ApJ, 762, 109 [NASA ADS] [CrossRef] [Google Scholar]
  10. Behroozi, P. S., Wechsler, R. H., Wu, H.-Y., et al. 2013b, ApJ, 763, 18 [NASA ADS] [CrossRef] [Google Scholar]
  11. Bertocco, S., Goz, D., Tornatore, L., et al. 2020, in Astronomical Society of the Pacific Conference Series, 527, 303 [NASA ADS] [Google Scholar]
  12. Bhattacharya, S., Heitmann, K., White, M., et al. 2011, ApJ, 732, 122 [NASA ADS] [CrossRef] [Google Scholar]
  13. Bocquet, S., Saro, A., Mohr, J. J., et al. 2015, ApJ, 799, 214 [NASA ADS] [CrossRef] [Google Scholar]
  14. Bocquet, S., Saro, A., Dolag, K., & Mohr, J. J. 2016, MNRAS, 456, 2361 [Google Scholar]
  15. Bocquet, S., Dietrich, J. P., Schrabback, T., et al. 2019, ApJ, 878, 55 [Google Scholar]
  16. Bond, J. R., Cole, S., Efstathiou, G., & Kaiser, N. 1991, ApJ, 379, 440 [NASA ADS] [CrossRef] [Google Scholar]
  17. Borgani, S., Rosati, P., Tozzi, P., et al. 2001, ApJ, 561, 13 [NASA ADS] [CrossRef] [Google Scholar]
  18. Bryan, G. L., & Norman, M. L. 1998, ApJ, 495, 80 [NASA ADS] [CrossRef] [Google Scholar]
  19. Casas-Miranda, R., Mo, H. J., Sheth, R. K., & Boerner, G. 2002, MNRAS, 333, 730 [NASA ADS] [CrossRef] [Google Scholar]
  20. Castignani, G., & Benoist, C. 2016, A&A, 595, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Castorina, E., Sefusatti, E., Sheth, R. K., Villaescusa-Navarro, F., & Viel, M. 2014, JCAP, 02, 049 [CrossRef] [Google Scholar]
  22. Castro, T., & Fumagalli, A. 2024, https://github.com/TiagoBsCastro/CCToolkit [Google Scholar]
  23. Castro, T., Marra, V., & Quartin, M. 2016, MNRAS, 463, 1666 [NASA ADS] [CrossRef] [Google Scholar]
  24. Castro, T., Borgani, S., Dolag, K., et al. 2020, MNRAS, 500, 2316 [Google Scholar]
  25. Cole, S., & Kaiser, N. 1989, MNRAS, 237, 1127 [NASA ADS] [Google Scholar]
  26. Comparat, J., Prada, F., Yepes, G., & Klypin, A. 2017, MNRAS, 469, 4157 [NASA ADS] [CrossRef] [Google Scholar]
  27. Cooray, A., & Sheth, R. K. 2002, Phys. Rept., 372, 1 [Google Scholar]
  28. Costanzi, M., Villaescusa-Navarro, F., Viel, M., et al. 2013, JCAP, 12, 012 [Google Scholar]
  29. Costanzi, M., Saro, A., Bocquet, S., et al. 2021, PRD, 103, 043522 [NASA ADS] [CrossRef] [Google Scholar]
  30. Courtin, J., Rasera, Y., Alimi, J. M., et al. 2011, MNRAS, 410, 1911 [NASA ADS] [Google Scholar]
  31. Crocce, M., Fosalba, P., Castander, F. J., & Gaztanaga, E. 2010, MNRAS, 403, 1353 [NASA ADS] [CrossRef] [Google Scholar]
  32. Cui, W., Borgani, S., & Murante, G. 2014, MNRAS, 441, 1769 [NASA ADS] [CrossRef] [Google Scholar]
  33. Davis, M., Efstathiou, G., Frenk, C. S., & White, S. D. M. 1985, ApJ, 292, 371 [Google Scholar]
  34. Deshpande, A. C., et al. 2024, A&A, 684, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Despali, G., Giocoli, C., Angulo, R. E., et al. 2016, MNRAS, 456, 2486 [NASA ADS] [CrossRef] [Google Scholar]
  36. Diemer, B. 2018, ApJS, 239, 35 [NASA ADS] [CrossRef] [Google Scholar]
  37. Diemer, B. 2020, ApJ, 903, 87 [NASA ADS] [CrossRef] [Google Scholar]
  38. Diemer, B., & Joyce, M. 2019, ApJ, 871, 168 [NASA ADS] [CrossRef] [Google Scholar]
  39. Dolag, K., Borgani, S., Murante, G., & Springel, V. 2009, MNRAS, 399, 497 [Google Scholar]
  40. Eke, V. R., Cole, S., & Frenk, C. S. 1996, MNRAS, 282, 263 [CrossRef] [Google Scholar]
  41. Elbers, W. 2022, JCAP, 11, 058 [CrossRef] [Google Scholar]
  42. Euclid Collaboration (Scaramella, R., et al.) 2022, A&A, 662, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Euclid Collaboration (Castro, T., et al.) 2023, A&A, 671, A100 [CrossRef] [EDP Sciences] [Google Scholar]
  44. Euclid Collaboration (Castro, T., et al.) 2024a, A&A, 685, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Euclid Collaboration (Fumagalli, A., et al.) 2024b, A&A, 683, A253 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Euclid Collaboration (Mellier, Y., et al.) 2024c, A&A, in press https://doi.org/10.1051/0004-6361/202450810 [Google Scholar]
  47. Fumagalli, A., Saro, A., Borgani, S., et al. 2021, A&A, 652, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. Fumagalli, A., Costanzi, M., Saro, A., Castro, T., & Borgani, S. 2024, A&A, 682, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Hamaus, N., Seljak, U., Desjacques, V., Smith, R. E., & Baldauf, T. 2010, PRD, D82, 043515 [Google Scholar]
  50. Hasselfield, M., Hilton, M., Marriage, T. A., et al. 2013, JCAP, 07, 008 [CrossRef] [Google Scholar]
  51. Hoffman, M. D., Gelman, A., et al. 2014, J. Mach. Learn. Res., 15, 1593 [Google Scholar]
  52. Holder, G., Haiman, Z., & Mohr, J. 2001, ApJ, 560, L111 [NASA ADS] [CrossRef] [Google Scholar]
  53. Hu, W., & Kravtsov, A. V. 2003, ApJ, 584, 702 [Google Scholar]
  54. Huterer, D., & Turner, M. S. 2001, PRD, 64, 123527 [NASA ADS] [CrossRef] [Google Scholar]
  55. Jenkins, A., Frenk, C. S., White, S. D. M., et al. 2001, MNRAS, 321, 372 [Google Scholar]
  56. Kitayama, T., & Suto, Y. 1996, ApJ, 469, 480 [CrossRef] [Google Scholar]
  57. Knebe, A., Knollmann, S. R., Muldrew, S. I., et al. 2011, MNRAS, 415, 2293 [Google Scholar]
  58. Kravtsov, A. V., & Borgani, S. 2012, ARA&A, 50, 353 [Google Scholar]
  59. Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv e-prints [arXiv:1110.3193] [Google Scholar]
  60. Lesci, G. F., Marulli, F., Moscardini, L., et al. 2022a, A&A, 659, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Lesci, G. F., Nanni, L., Marulli, F., et al. 2022b, A&A, 665, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. LoVerde, M. 2014, Phys. Rev. D, 90, 083518 [NASA ADS] [CrossRef] [Google Scholar]
  63. Mana, A., Giannantonio, T., Weller, J., et al. 2013, MNRAS, 434, 684 [NASA ADS] [CrossRef] [Google Scholar]
  64. Mantz, A. B., von der Linden, A., Allen, S. W., et al. 2015, MNRAS, 446, 2205 [Google Scholar]
  65. Michaux, M., Hahn, O., Rampf, C., & Angulo, R. E. 2021, MNRAS, 500, 663 [Google Scholar]
  66. Mo, H. J., & White, S. D. M. 1996, MNRAS, 282, 347 [Google Scholar]
  67. Monaco, P., Theuns, T., & Taffoni, G. 2002, MNRAS, 331, 587 [Google Scholar]
  68. Monaco, P., Sefusatti, E., Borgani, S., et al. 2013, MNRAS, 433, 2389 [NASA ADS] [CrossRef] [Google Scholar]
  69. Munari, E., Monaco, P., Sefusatti, E., et al. 2017, MNRAS, 465, 4658 [Google Scholar]
  70. Ondaro-Mallea, L., Angulo, R. E., Zennaro, M., Contreras, S., & Aricò, G. 2021, MNRAS, 509, 6077 [NASA ADS] [CrossRef] [Google Scholar]
  71. Peebles, P. J. E. 2020, The Large-scale Structure of the Universe, 98 (Princeton University Press) [Google Scholar]
  72. Planck Collaboration XX. 2014, A&A, 571, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  73. Planck Collaboration XXIV. 2016, A&A, 594, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  74. Press, W. H., & Schechter, P. 1974, ApJ, 187, 425 [Google Scholar]
  75. Romanello, M., Marulli, F., Moscardini, L., et al. 2024, A&A, 682, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Rozo, E., Wechsler, R. H., Rykoff, E. S., et al. 2010, ApJ, 708, 645 [Google Scholar]
  77. Salvati, L., Douspis, M., & Aghanim, N. 2020, A&A, 643, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  78. Salvatier, J., Wiecki, T. V., & Fonnesbeck, C. 2016, PeerJ Comput. Sci., 2, e55 [Google Scholar]
  79. Saro, A., Bocquet, S., Rozo, E., et al. 2015, MNRAS, 454, 2305 [Google Scholar]
  80. Sartoris, B., Biviano, A., Fedeli, C., et al. 2016, MNRAS, 459, 1764 [Google Scholar]
  81. Schneider, A., & Teyssier, R. 2015, JCAP, 12, 049 [Google Scholar]
  82. Sheth, R. K., & Tormen, G. 1999, MNRAS, 308, 119 [Google Scholar]
  83. Sheth, R. K., Mo, H. J., & Tormen, G. 2001, MNRAS, 323, 1 [NASA ADS] [CrossRef] [Google Scholar]
  84. Springel, V., White, S. D. M., Tormen, G., & Kauffmann, G. 2001, MNRAS, 328, 726 [Google Scholar]
  85. Springel, V., Pakmor, R., Zier, O., & Reinecke, M. 2021, MNRAS, 506, 2871 [NASA ADS] [CrossRef] [Google Scholar]
  86. Sunayama, T., Miyatake, H., Sugiyama, S., More, S., et al. 2023, arXiv e-prints [arXiv:2309.13025] [Google Scholar]
  87. Taffoni, G., Becciani, U., Garilli, B., et al. 2020, in Astronomical Society of the Pacific Conference Series, 527, 307 [NASA ADS] [Google Scholar]
  88. Tinker, J. L., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709 [NASA ADS] [CrossRef] [Google Scholar]
  89. Tinker, J. L., Robertson, B. E., Kravtsov, A. V., et al. 2010, ApJ, 724, 878 [NASA ADS] [CrossRef] [Google Scholar]
  90. To, C., Krause, E., Rozo, E., et al. 2021, PRL, 126, 141301 [NASA ADS] [CrossRef] [Google Scholar]
  91. Vehtari, A., Gelman, A., & Gabry, J. 2017, Statist. Comput., 27, 1413 [CrossRef] [Google Scholar]
  92. Velliscig, M., van Daalen, M. P., Schaye, J., et al. 2014, MNRAS, 442, 2641 [Google Scholar]
  93. Villaescusa-Navarro, F., Naess, S., Genel, S., et al. 2018, ApJ, 867, 137 [NASA ADS] [CrossRef] [Google Scholar]
  94. Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261 [CrossRef] [Google Scholar]
  95. Warren, M. S., Abazajian, K., Holz, D. E., & Teodoro, L. 2006, ApJ, 646, 881 [NASA ADS] [CrossRef] [Google Scholar]
  96. Watanabe, S. 2010, J. Mach. Learn. Res., 11, 3571 [Google Scholar]
  97. Watson, W. A., Iliev, I. T., D’Aloisio, A., et al. 2013, MNRAS, 433, 1230 [Google Scholar]
  98. Zhang, T., Chuang, C.-H., Wechsler, R. H., et al. 2022, MNRAS, 518, 3737 [NASA ADS] [CrossRef] [Google Scholar]

6

There is no fundamental reason for choosing S8 over other param֊ eterizations, such as σ8m,0/0.3)α with α free to vary. However, we found empirically that S8 works well within our analysis, and thus, we opted to adopt this widely used variable in the literature.

All Tables

Table 1

Cosmological parameters of the PICCOLO set of simulations.

Table 2

Best fit and 95% limits on the model parameters presented in Eqs. (24)(27).

All Figures

thumbnail Fig. 1

Constraints on the parameters in Eq. (19) fitted to the unbiased standard deviation of the 200 PINOCCHIO mocks with C0 cosmology. We fitted α and β considering k ≤ 0.05 h Mpc−1, assuming a Gaussian likelihood with error bar estimated from the measurements for 0.05 ≤ k/h Mpc−1 ≤ 0.2 and assuming it to be constant in k and equivalent to the unbiased standard deviation.

In the text
thumbnail Fig. 2

Relative difference between Eq. (19) best fit and the unbiased standard deviation of the PINOCCHIO measurements. Different columns are for different redshifts, and the corresponding mass bins are shown in each panel. The vertical dotted line demarcates two distinct sets of measurements: Those to the left of the line were utilized as data points in the parameter fitting process, while the scatter of the points to the right was analyzed to estimate the variance. The gray regions highlight areas within a 5% deviation from the expected values.

In the text
thumbnail Fig. 3

Correlation coefficient ρ between the measurements of the halo-matter power spectrum, Phm(k), in different simulations for k ≤ 0.05 h Mpc−1 as a function of halo mass bin and redshift. The red histograms show the distribution of the correlation coefficient between a simulation and its paired realization, while blue histograms are for the correlation coefficient between simulations with uncorrelated white-noise realizations.

In the text
thumbnail Fig. 4

Relative difference between the bias measured in halo catalogs and the bias predicted by the PBS model of Eq. (16) at different red-shifts in the range z ∊ [0, 2], Results refer to simulations carried out for the C0 cosmology. In each panel, blue and red lines refer to the results obtained for the halo catalogs based on the application of ROCKSTAR and SUBFIND, respectively. For the ROCKSTAR catalogs, we also show the standard error of the mean using other 19 realizations of the same cosmology.

In the text
thumbnail Fig. 5

Mean of the ratio of the bias measured in simulations with respect to the PBS prescription. Left: ratio of the bias as a function of v/(1 + z) for different z, labeled by Ωm(z), for all the CO runs. Center, mean ratio as a function of v for the three EdS cosmologies of pure power-law shapes of the linear power spectrum at z = 0. We report the value of the three spectral indexes in the inset. Right: mean ratio as a function of the background evolution Ωm(z) for the C0, C9, and C10 cosmologies with varying S8.

In the text
thumbnail Fig. 6

Marginalized 68% and 95% confidence level contours on the model parameters presented in Eqs. (24) to (27). We calibrated our model using the subset of PICCOLO simulations C0–C10. (See Table 2 for the best fit and confidence levels.)

In the text
thumbnail Fig. 7

Ratio between the mean of the observations on the 20 C0 simulations with respect to our model predictions. Different rows correspond to different redshifts, while each panel corresponds to a different mass bin. The shaded region in red corresponds to the error on the mean, assuming that each measurement follows Eq. (19). The shaded regions in gray correspond to 2% and 4% regions.

In the text
thumbnail Fig. 8

Similar to Fig. 7 but for the C9 and C10 cosmological parameters. Among the PICCOLO set, C9 and C10 correspond to the cosmologies with the lowest and the highest S8, respectively.

In the text
thumbnail Fig. 9

Similar to Fig. 7 but for simulations with massive neutrinos. For better plot readability, we only show the uncertainties (red shaded regions) for the simulation with total neutrino masses equal to 0.15 eV.

In the text
thumbnail Fig. 10

Comparison between the HB predicted by our model with predictions from other models presented in the literature: Cole & Kaiser (1989). Sheth et al. (2001), Tinker et al. (2010), and Comparat et al. (2017). We present both our benchmark model as well as the PBS predictions based on the HMF model of Euclid Collaboration (2023) used as a baseline of our model. Different columns correspond to different redshifts. The relative difference with respect to our benchmark model is presented in the panels in the second row. We adopted a composite scale for the residual plot to show the dynamic range of differences between the models: The scale is linear for values between [−10, 10] % and symmetric log outside. For reference, we show the zero line in black. The predictions of the models from the literature have been computed using the COLOSSUS toolkit (Diemer 2018).

In the text
thumbnail Fig. 11

Percentage residuals of cluster counts (left panel) and cluster clustering covariance matrices (central and right panels), computed with the bias from Tinker et al. (2010) in comparison to the one calculated using the bias calibrated in this study (Eq. (24)). We show the full covariance matrix for number counts in mass and redshift bins. For the two-point correlation function of galaxy clusters, we show two blocks of the full covariance (low and high redshift bins) as a function of the radial separation.

In the text
thumbnail Fig. 12

Parameter posteriors at 68% and 95% confidence levels obtained by analyzing number counts, weak lensing masses, and cluster clustering computed with the HB calibrated in this work (cyan contours) and the bias from Tinker et al. (2010) (orange contours). The error associated with weak lensing mass is set at 1% of the mass.

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.