Issue 
A&A
Volume 652, August 2021



Article Number  A155  
Number of page(s)  25  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/202039123  
Published online  27 August 2021 
The mass function dependence on the dynamical state of dark matter haloes
^{1}
MaxPlanckInstitut für extraterrestrische Physik (MPE), Giessenbachstrasse 1, 85748 Garching bei München, Germany
email: rseppi@mpe.mpg.de
^{2}
Instituto de Astrofísica de Andalucía (CSIC), Glorieta de la Astronomía, 18080 Granada, Spain
^{3}
Astronomy Department, New Mexico State University, Las Cruces, NM, USA
^{4}
Department of Astronomy, University of Virginia, Charlottesville, VA, USA
Received:
7
August
2020
Accepted:
25
May
2021
Context. Galaxy clusters are luminous tracers of the most massive dark matter haloes in the Universe. To use them as a cosmological probe, a detailed description of the properties of dark matter haloes is required.
Aims. We characterize how the dynamical state of haloes impacts the dark matter halo mass function at the highmass end (i.e., for haloes hosting clusters of galaxies).
Methods. We used the dark matteronly MultiDark suite of simulations and the highmass objects M > 2.7 × 10^{13} M_{⊙} h^{−1} therein. We measured the mean relations of concentration, offset, and spin as a function of dark matter halo mass and redshift. We investigated the distributions around the mean relations. We measured the dark matter halo mass function as a function of offset, spin, and redshift. We formulated a generalized mass function framework that accounts for the dynamical state of the dark matter haloes.
Results. We confirm the recent discovery of the concentration upturn at high masses and provide a model that predicts the concentration for different values of mass and redshift with one single equation. We model the distributions around the mean values of concentration, offset, and spin with modified Schechter functions. We find that the concentration of lowmass haloes shows a faster redshift evolution compared to highmass haloes, especially in the highconcentration regime. We find that the offset parameter is systematically smaller at low redshift, in agreement with the relaxation of structures at recent times. The peak of its distribution shifts by a factor of ∼1.5 from z = 1.4 to z = 0. The individual models are combined into a comprehensive mass function model, which predicts the mass function as a function of spin and offset. Our model recovers the fiducial mass function with ∼3% accuracy at redshift 0 and accounts for redshift evolution up to z ∼ 1.5.
Results. This new approach accounts for the dynamical state of the halo when measuring the halo mass function. It offers a connection with dynamical selection effects in galaxy cluster observations. This is key toward precision cosmology using cluster counts as a probe.
Key words: cosmology: theory / galaxies: clusters: general / Xrays: galaxies: clusters / dark matter / cosmology: observations
© R. Seppi et al. 2021
Open 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.
Open Access funding provided by Max Planck Society.
1. Introduction
Galaxy clusters are the most massive virialized, gravitationally bound structures in the Universe. They grow hierarchically, starting from matter perturbations in the initial density field. This makes them good tracers of the underlying cosmic web and its densest regions. Clusters in cosmology are used to construct the halo mass function, which indicates the mass density of haloes in a specific volume, in a small mass interval between M and M + dM (Weinberg et al. 2013). An early theoretical description of the mass function was given by Press & Schechter (1974; PS) based on the assumption that Gaussian density perturbations overcoming a fixed density contrast collapse into haloes. This formally accounts for only half of the total halo mass in the Universe. An alternative approach, employing the excursion set theory, solved these shortcomings by considering the probability of crossing a given barrier with random walks (Bond et al. 1991). This provides a good prediction for highmass haloes, but it predicts too many lowmass objects. The introduction of ellipsoidal collapse corrected these differences between simulations and theory (Sheth & Tormen 1999, 2002).
The mass function has been extensively studied in more recent works in an attempt to find a universal model that is independent from cosmology (Jenkins et al. 2001; Tinker et al. 2008; Bhattacharya et al. 2011; Despali et al. 2016; Bocquet et al. 2016, 2020; Comparat et al. 2017). A robust way to build such a mass function model is using Nbody simulations (Kravtsov et al. 1997; Springel 2005), for example MultiDark (Prada et al. 2012; Klypin et al. 2016). The generalization of such models as a function of cosmological parameters is best handled by emulating the mass function based on large sets of simulations (e.g., McClintock et al. 2019; Nishimichi et al. 2019; Bocquet et al. 2020). It is important to precisely predict the halo mass function to fulfill the potential of current and future Xray, SZ, or optical cluster surveys such as eROSITA (Merloni et al. 2012; Predehl et al. 2021), Planck (Zubeldia & Challinor 2019), SPT3G (Benson et al. 2014), CMB S4 (Abazajian et al. 2019), SPIDERS/eBOSS (Dawson et al. 2016; Finoguenov et al. 2020), DESI (DESI Collaboration 2016), 4MOST (de Jong 2011; Finoguenov et al. 2019), Euclid (Laureijs et al. 2011), LSST (LSST Science Collaboration 2009), and WFIRST (Spergel et al. 2015). These future surveys will provide tighter constraints on massobservable scaling relations. This means that systematic uncertainties due to the accuracy of the halo mass function model and its evolution with redshift will contribute to the total error budget significantly more than in the previous cluster counts experiments (Salvati et al. 2020). Therefore, a detailed prediction of the theoretical dark matter halo statistic is required. In this work we calibrate a mass function model that includes the dynamical properties of dark matter haloes.
Relaxed dark matter haloes can be selected according to multiple diagnostics (Neto et al. 2007; Macciò et al. 2008; Prada et al. 2012; Klypin et al. 2016): (i) the virial parameter 2K/W−1, where K and W are respectively the kinetic and potential energy within the virial radius; (ii) the spin parameter (Peebles 1969), which traces the dynamical state of the halo; (iii) the fraction of substructures, which is higher in unrelaxed haloes; (iv) the offset parameter X_{off} = R_{peak} − R_{cm}/R_{vir} (Behroozi et al. 2013; Klypin et al. 2016), which is the difference between the position of the peak of the density profile and the center of mass, normalized by virial radius. If the halo is perfectly relaxed, the peak of the profile will correspond to the center of mass, and X_{off} will be small. On the other hand, higher X_{off} values will indicate an unrelaxed halo (e.g., merger, accretion). A combination of these quantities can be used.
To illustrate the respective contribution to the mass function of relaxed and disturbed haloes, we divided the sample of haloes in HugeMultiDark (hereafter HMD, see Sect. 3) at z = 0 using X_{off} and λ. We consider the offset parameter in physical scale X_{off, P} = R_{peak} − R_{cm} (i.e., not normalizing by the virial radius). We consider haloes with X_{off, P} < 100 kpc h^{−1} and λ < 0.007 as relaxed. In Fig. 1 we show the halo mass function of the complete halo population: the relaxed and disturbed halo population. The left panel shows the distribution of redshift zero haloes in the X_{off, P} and λ plane. The right panel shows the three multiplicity functions (see formalism in Sect. 2) sampled by all haloes (green), relaxed haloes (blue), and disturbed haloes (orange). The bottom panel shows their relative contribution as a function of halo mass. It is clear how at 10^{13.5} M_{⊙} h^{−1} the contribution from relaxed structures dominates by a factor of about 0.8 dex. The multiplicity functions of the two samples cross each other at 10^{14.5} M_{⊙} h^{−1}, then the unrelaxed structures take over at the highmass end. This approach offers a possible connection to selection effects in observations, such as the cool core bias in Xrays (Eckert et al. 2011; Käfer et al. 2019, 2020). It possibly offers a solution to mitigate biases in a cosmological interpretation of clusters abundance. This might improve cosmological constraints using Xray selected clusters (Ider Chitham et al. 2020).
Fig. 1. Contribution of relaxed and unrelaxed haloes to the total mass function. Left panel: distribution of redshift zero dark matter haloes in the X_{off, P} and λ plane. Cuts in X_{off, P} and λ are applied to divide relaxed (blue) and disturbed (orange) structures. The contours contain 1%, 5%, 10%, 30%, and 50% of the data. Right panel: Halo mass functions vs. mass (σ) at redshift 0, as defined in Sect. 2 and using Eqs. (15) and (16). This mass function is built with different subsets of haloes from HMD at z = 0 (see Sect. 3): the red line indicates the model from Comparat et al. (2017), the shaded areas represent the relaxed (blue), unrelaxed (orange), and the full sample (green) of haloes. The areas cover 1σ uncertainties. Lower panel: residuals fraction of each component compared to the red model in the upper panel. 
In this article, we investigate the variations of the dark matter halo mass function as a function of the dynamical state of the constituting haloes. To trace the dynamical state we use X_{off} and λ. This paper is structured as follows. We define the formalism in Sect. 2. We present the Nbody data used in Sect. 3. We present the average relations between concentration, X_{off}, spin, and mass, as well as their distribution around mean values, in Sect. 4. We define the generalized mass function framework in Sect. 5. We present the generalized model of the halo mass function as a function of X_{off} and λ in Sect. 6. We present the results of the fit and the bestfit parameters in Sect. 7. We summarize our findings in Sect. 8, and discuss their implications for cosmological studies.
2. Formalism and definitions
The growth of the density perturbations in the matter field is described by the evolution of the overdensity field and its variance as a function of scale. The variance of the smoothed density field is defined as
where k is the wavenumber, the Fourier transform of a tophat filtering function, and P(k, z) = D_{+}(z)^{2}P(k, 0) the linear matter power spectrum. Its redshift evolution is encoded in the growth factor D_{+}(z). The quantity
is the peak height, where δ_{cr} is the critical overdensity required for a structure to collapse in a dark matter halo (Press & Schechter 1974; Sheth & Tormen 2002). We report explicit values of mass, σ, and ν in Table 1 for z = 0 and z = 0.5.
Correspondence between mass, peak height, and variance of the linear density field at z = 0 and z = 0.5.
We write the mass function in its differential form as in Eq. (3) (see Allen et al. 2011, for a recent review),
where f(σ) is the multiplicity function. A comprehensive list of models of the multiplicity function is available in Table 1 of Murray et al. (2013).
In this paper as mass variable we use σ (Eq. (1)), peak height (Eq. (2)), and mass interchangeably.
3. Simulations
We describe the set of gravityonly simulations and the halo finding postprocess in this section.
3.1. MultiDark
We use the MultiDark simulations (Prada et al. 2012; Riebe et al. 2013; Klypin et al. 2016). They are computed in a flat ΛCDM Planck (Planck Collaboration XI 2014) cosmology (H_{0} = 67.77 km s^{−1} Mpc^{−1}, Ω_{m0} = 0.307115, Ω_{b0} = 0.048206, σ_{8} = 0.8228) with the GADGET2 code (Springel 2005). It is one of the largest sets of highresolution (∼4000^{3} particles) Nbody simulations.
We used three MultiDark simulations: HMD, BIGMD, MDPL2 (see details in Table 2). Alternative simulations that could be used for this project include MILLENNIUMXXL, DARKSKIES, Q CONTINUUM, V^{2}GC SIMULATION, described in Angulo et al. (2012), Skillman et al. (2014), Heitmann et al. (2015), and Ishiyama et al. (2015) respectively.
Nbody simulations used in this analysis.
The complete list of simulation outputs (snapshots) utilized are given in Table A.6, where the expansion parameter a and the corresponding redshifts are reported for each snapshot.
3.2. Halo finding
Finding haloes in dark matter simulations is not an easy task (see Knebe et al. 2013; Behroozi et al. 2015, for a review). In this study haloes are identified by the Robust Overdensity Calculation using KSpace Topologically Adaptive Refinement (ROCKSTAR) and CONSISTENTTREES algorithms (Behroozi et al. 2013). They are based on adaptive hierarchical refinement of friendsoffriends (FOF) groups. They work with six phasespace dimensions (halo positions and velocities), and one time dimension. This allows us to track relative motions and merging history between substructures in different snapshots. ROCKSTAR computes the halo mass of identified objects by removing the unbound particles inside the virial radius. Virial mass and virial radius are related by
where R_{vir} encompasses a mean halo density equal to the background matter density multiplied by Δ_{vir}, Ω_{M} is the matter density parameter, and ρ_{b} is the matter density of the Universe. The overdensity over the matter background Δ_{vir} is defined according to Bryan & Norman (1998) as
where x = Ω(z)−1, Ω_{M, 0} is the matter density parameter at the present day, and E(z) is the Hubble parameter as a function of redshift in units of H_{0}:
Given the reference cosmology adopted in this work, the virial overdensity is equal to 332.5 at z = 0 and asymptotically tends to 178 at high redshift.
The recovery of main properties such as position, mass, and circular velocity is consistent between different finders. However, derived properties such as spin show a ∼20% scatter (Knebe et al. 2013). This holds especially for lowmass haloes with fewer than 30–40 particles, where the identification of substructures is not straightforward (Knebe et al. 2011). The lowmass limits in Table 2 are set at more than ∼200 particles per halo in all boxes. This ensures accurate halo properties (Behroozi et al. 2013; Knebe et al. 2013).
In this work we use the virial overdensity (Eq. (5)). The virial mass function has been shown to be the one that comes closest to universality (Despali et al. 2016).
4. Concentration, offset, spin: Empirical relations with peak height and redshift
In the footsteps of Klypin et al. (2016) and RodriguezPuebla et al. (2016), we analyze the average relations linking concentration, λ, and X_{off} to the peak height, respectively in Sects. 4.1–4.3. We also analyze the distributions of these quantities around the mean relations. The mean relations are fitted by models that simultaneously account for the mass and redshift dependence of the relations (Eqs. (9), (11), and (13)). The probability density functions (PDFs) of concentration, spin, and X_{off} are fitted by modified Schechter models, respectively in Eqs. (10), (12), and (14). The PDFs at different redshifts are modeled independently.
4.1. Concentration–mass–redshift relation
We study the relation between concentration and mass. Numerical simulations (Navarro et al. 1996) show that to a good approximation (∼10−20%), the density of dark matter haloes is described by the profile in Eq. (6),
where R_{s} is the scale radius. The characteristic density of the halo ρ_{s} is equal to
where c is the concentration, Δ the overdensity, and ρ_{crit} is the critical density of the Universe. The concentration is a dimensionless quantity defined by
where Δ can refer to any threshold. In this work we consider the virial overdensity over the matter background Δ_{vir} (see Eq. (5)). Therefore, we study the virial concentration c_{vir} = R_{vir}/R_{s}. The concentration–mass relation is an important part of models describing galaxy clusters (e.g., reviews from Allen et al. 2011; Umetsu 2020) or gravitational lensing (e.g., reviews from Bartelmann 2010; Kilbinger 2015).
This relation has been extensively studied in simulations. The concentration anticorrelates with mass with negative redshift trend (e.g., Diemer & Joyce 2019; Ragagnin et al. 2019). Its detailed trend depends on the measurement method, in simulations (Meneghetti & Rasia 2013; Lang et al. 2015; PovedaRuiz et al. 2016) and in observations (Foëx et al. 2014; Phriksee et al. 2020; Du et al. 2015; Shan et al. 2017). This introduces possibles biases in the measure of concentration (Sereno et al. 2015; Cibirka et al. 2017; van Uitert et al. 2016). Diemer & Kravtsov (2015) found that the relation shows the smaller deviation from universality when adopting the definition c_{200c}. However, it is not completely universal, meaning that concentration is described not only by mass or ν, but also by assembly history. Leaving aside biases in definition and measurement, different models have been proposed to describe it: power laws (e.g., Duffy et al. 2008; Dutton & Macciò 2014); a combination of power laws to describe the highmass upturn (e.g., Klypin et al. 2016; Diemer & Joyce 2019); semianalytic models based on the PressSchechter theory (Correa et al. 2015). In this section we extend the models from Klypin et al. (2016). We adjust a global model that includes a redshift dependence for highmass haloes at relatively low redshift (z < 1.5), which is particularly interesting for haloes hosting galaxy clusters.
4.1.1. Model
We model the concentration c as a function of the rms of the overdensity field σ(M, z) (not as the function of halo mass M). Our parameterization of the concentration–σ relation is a generalization of the Klypin et al. (2016) relation and reads
The bestfit values are in Table A.1. We find bestfit values of a_{0} = 0.754091 ± 0.000004, b_{0} = 0.574413 ± 0.000002, in agreement with Klypin et al. (2016). This model is fitted using haloes with M > 10^{12.5} M_{⊙} h^{−1} (ν ∼ 0.95 at z = 0). The lowmass end is not sampled due to the particle mass resolution. We do not consider a highredshift regime where there is a statistical limitation of highmass haloes. We chose the binning limits following Klypin et al. (2016). The strength of our model resides in the ability to predict the concentration–σ relation for a variety of masses and redshifts with a single equation. We recover the hockey stick shape of the relation (see Fig. 2). This is consistent with Prada et al. (2012), Klypin et al. (2016), and Diemer & Joyce (2019) (the last also includes a cosmology dependence). This upturn feature is absent in Duffy et al. (2008) and Wang et al. (2020), due to the smaller volumes analyzed (400 Mpc h^{−1} and 500 Mpc h^{−1}, respectively). They are not large enough to obtain a significant number of highmass objects and probe the upturn.
Fig. 2. Concentration–σ relation (Eq. (9)). Circular dots, triangles, and squares represent HMD, BIGMD, MDPL2, respectively. They are colorcoded by redshift. Straight lines indicate our bestfit model, dotted lines show the model from Klypin et al. (2016), and the shaded blue area indicates the distribution from Wang et al. (2020) at z = 0. The upper xaxis converts peaks values into mass at z = 0. 
We describe the distribution of concentration around its mean value and we model its probability density function. We use snapshots from HMD, at four redshifts. Each snapshot is divided into six slices of mass. The PDFs obtained at each redshift are fitted simultaneously. It includes a σ dependence in our model (see Eq. (10)):
The bestfit parameters are included in Table A.1. Uncertainties are negligible for most parameters. The distribution of concentration around its mean value is well described by a modified Schechter function with massdependent terms (Eq. (10)). We use the bestfit model to predict the distribution of concentration in fixed mass slices as a function of redshift. The result is shown in Fig. 3. The full distribution of concentration (not sliced in different mass intervals) is well described by a modified Schechter as well.
Fig. 3. Probability density function of the concentration (Eq. (10)) at the different redshifts values given in each panel. Each set is divided into mass slices, and is colorcoded accordingly. The shaded areas represent the data with 1σ error, while the straight lines indicate the bestfit model. The blue points and line represent the total sample not sliced in mass. For clarity, each line and its fit is shifted by 0.1 dex on the yaxis. This means that the constant C_{0} assumes values of (+0.3, +0.2, +0.1, 0.0, −0.1, −0.2). The purple line is not shifted, and is thus the one with the proper normalization. 
4.1.2. Discussion
The average value of concentration increases at late times. It goes from 4 at z = 1.43 to 5.8 at z = 0 for haloes corresponding to peak height ν = 2. We confirm the recent discovery of the concentration upturn at high masses (Klypin et al. 2016; Diemer & Joyce 2019). Klypin et al. (2016) suggest that the highmass upturn is caused by the tendency of regions with smaller root mean square variance in the overdensity field to be more spherical. These regions are the ones that evolve into highmass dark matter haloes. Because of this aspect, gravitational accretion of matter toward the center is more efficient and results in higher concentration. In addition, very massive haloes correspond to the highest peaks of the density field. This means that the probability of finding them in a phase characterized by strong accretion and therefore higher concentration is higher, as shown by Ludlow et al. (2012).
We find that the concentration of haloes with different masses shows a different evolution with redshift (see Fig. 3). At z = 0 the concentration of highmass haloes gives a smaller contribution to the total PDF in the highconcentration tail. At higher redshift we see highconcentration haloes with both high and low mass. The opposite holds for lowmass haloes, which contribute more to the total PDF at higher z in the lowconcentration regime. At high redshift the distribution of concentration for different mass bins has the same shape. Conversely, at low redshift this statement does not hold. The lowmass distribution is much broader than the highmass one. There is also a general redshift trend: with decreasing redshift the number of highconcentration, lowmass haloes increases. More recently in time the distribution is flatter; in fact, the slope of the powerlaw is smaller, going from ∼4.5 at z = 1.43 to ∼1.4 in the present day. Moreover, the mass trend of the exponential decay at z = 0 is negative with sigma (e_{2} = −0.959), which translates into a faster decrease at high mass. This means that the distribution of concentration of lowmass haloes evolves more than that of highmass objects. For example, at concentration c = 8 the PDF of highmass haloes (1 × 10^{14} < M < 2 × 10^{14} M_{⊙} h^{−1}) changes by 0.32 dex between z = 1.43 and z = 0, while for lowmass objects (2 × 10^{13} < M < 3 × 10^{13} M_{⊙} h^{−1}) it varies by 0.49 dex. There is a net difference of 0.17 dex between the two, meaning that the high end of concentration evolves ∼1.5 times faster for lower mass haloes compared to highmass ones. This is in agreement with the fact that these haloes evolve in different environments. This causes a different redshift evolution of the distribution of concentration of haloes with different masses. Structures with the same mass but different formation histories present different properties. This aspect is related to the notion of assembly bias (Gao et al. 2005; Croton et al. 2007; Angulo et al. 2008). There is a correlation between halo concentration and the age of the universe when the progenitor reached a fixed fraction of its current mass (Zhao et al. 2009). Moreover, a faster mass accretion causes a slower increase of concentration. Very massive objects reside in highly populated environments, which slows down the evolution of concentration (Zhao et al. 2003). This can be used to predict concentration as a function of assembly history (Giocoli et al. 2012; Ludlow et al. 2013, 2014). The latter is also related to the evolution of cosmological parameters. This allows modeling of the concentration mass relation for cold and warm dark matter haloes (Ludlow et al. 2016). In general, highmass haloes cluster more than lowmass ones. Lowmass objects, some of which are found in isolated environments, will experience different histories. Isolated lowmass haloes will be highly concentrated, while the concentration of lowmass haloes in dense environments will stay low. Contreras et al. (2019) found that lowconcentration haloes cluster more than haloes with high concentration at high redshift. This is in agreement with our results.
4.2. Spin–mass–redshift relation
In this subsection we study the relation between the spin parameter and mass. An accurate theoretical description of the distribution of the halo spin is key to understanding the possible implications of the systemic rotation on cluster sample definitions. This rotation is induced by a combination of the initial spin of the halo, the infalling material, and the merging activity.
Measurements of cluster rotation were obtained at low redshift, using member galaxies to infer the rotation movement (Hwang & Lee 2007; Tovmassian 2015; Manolopoulou & Plionis 2017; Bilton et al. 2019). The spin in galaxy clusters has also been studied in Xrays (Bianconi et al. 2013; Eckert et al. 2019). Highresolution data are needed to explore this topic, in the optical (Song et al. 2018) and Xray (Hitomi Collaboration 2018) bands. Moreover, cluster cores can be analyzed by distortion of the 6.7 keV line (Sunyaev et al. 2003), fluctuations in the CMB due to rotation (Rephaeli 1995; Cooray & Chen 2002), and rotational kSZ effect (Baxter et al. 2019). Inferring motion with the SZ effect was also demonstrated for relaxed clusters with a significant spin on simulated (hydrodynamically) clusters (Baldi et al. 2018, 2019). Future SZ and Xrays surveys might enable such measurements on a large number of clusters. A statistical description of the halo population as a function of spin is thus of interest, and is developed in this section.
4.2.1. Model
We analyze the spin–mass–redshift relation and its PDF. The spin is defined as λ = JE^{1/2}/GM^{5/2} (Peebles 1969). We fit the mean relation with mass as a linear relation (Eq. (11)).
There is no noticeable dependence on redshift, so we do not consider a redshift evolution in this model. The bestfit parameters are reported in Table A.2. The relation is shown in Fig. 4. We find that the spin correlates weakly with halo mass, as Bett et al. (2007) and RodriguezPuebla et al. (2016) did. The PDF of the spin parameter is best fit by a modified Schechter law with no mass dependence:
Fig. 4. Spin–mass relation (λσ, Eq. (11)). The model is a linear relation with no redshift trend. Data points are colorcoded by redshift, while the different geometrical shapes refer to different simulations: squares for HMD, triangles for BigMD, and circles for MDPL. The straight line indicates the bestfit model, which considers all simulations and redshift at the same time. The bestfit parameters are given in Table A.2. The upper xaxis converts peaks values into mass at z = 0. 
Distributions for different redshift snapshots are shown in Fig. 5 (see Table A.2 for bestfit parameters). The result is consistent with previous findings, (e.g., RodriguezPuebla et al. 2016). The evolution with redshift of the PDF shows that with time haloes build up higher spins: the maximum of the PDF shifts to higher spin values when redshift decreases (Fig. 5). We also tried a lognormal distribution as a model, but it was not successful.
Fig. 5. PDF of λ (Eq. (12)). Individual points represent spin bins used to compute the distribution; straight lines refer to the bestfit modified Schechter model. They are colorcoded by redshift. We do not consider mass dependence in this relation. Each redshift slice is fitted independently. The bestfit parameters are given in Table A.2. 
4.2.2. Discussion
The small correlation with mass and the wellmodeled evolution with redshift makes it a rather simple dependence to account for in statistical studies of the halo population. From the perspective of the measurement of the halo mass function based on a cluster sample, marginalizing over the spin is possible, with limited complications. It also shows the spin cannot be considered a candidate for assembly bias.
4.3. Offset–mass–redshift relation
The offset parameter X_{off} traces the relaxation state of the halo (Thomas et al. 2001; Neto et al. 2007; Henson et al. 2017). Hollowood et al. (2019) analyzed the miscentering in SDSS galaxies followed up with Chandra Xray observations. Provided there is a link between X_{off} and the brightest cluster galaxy (BCG) to Xray displacement, an estimation of the bivariate mass and X_{off} function from observations of clusters is possible. To interpret it requires a detailed description of the link between X_{off} and mass, which is detailed in this section. We find that lowmass haloes have smaller offset than highmass ones at each redshift. Moreover, the offset parameter is reduced by a factor of ∼1.5 from z ∼ 1.5 to z = 0. This is in agreement with the fact that structures relax in time.
4.3.1. Model
We model the X_{off} − σ relation with a redshift dependent power law, see Eq. (13),
where E(z) is the dimensionless Hubble parameter. The bestfit parameters are given in Table A.3. We find a_{0} = −1.30418 ± 0.00001, b_{0} = 0.15084 ± 0.00001. We find a significant redshift evolution of the normalization and the slope of the relation. The data obtained from MultiDark simulations and the bestfit models are shown in Fig. 6. The bestfit parameters of the X_{off} − σ relation and the distribution of X_{off} are given in Table A.3. As for concentration, the strength of this equation relies upon its ability to predict an average X_{off} value given the mass and the redshift of a dark matter halo.
Fig. 6. Offset–mass relation (X_{off}–σ, Eq. (13)). Circular dots, triangles, and squares represent HMD, BIGMD, and MDPL2 respectively. They are colorcoded by redshift. Straight lines show the bestfit model. The bestfit parameters are given in Table A.3. The upper xaxis converts peaks values into mass at z = 0. 
We describe the distribution of X_{off} around its mean value with a modified Schechter function (Eq. (14)). The PDF of X_{off} does not show mass dependency. It is included in the normalization to the virial radius . Therefore, we do not consider any σ dependence:
We fit all haloes in each redshift snapshot together. Figure 7 shows distributions at different redshifts, fitted independently from one another. The parameters are given in Table A.3.
Fig. 7. Probability density function of X_{off}. Each panel shows the distribution at a specific redshift. Scatter points indicate the data, while straight lines represent the modified Schechter model. The samples are colorcoded by redshift. Each redshift slice is fitted independently by Eq. (14). The bestfit parameters are given in Table A.3. 
4.3.2. Discussion
We find a nonzero slope for the relation between X_{off} and mass. On average, highmass haloes have a larger offset parameter. This is described by the negative a_{0} parameter. This is in agreement with the hierarchical picture of structure formation. Highmass haloes formed recently and have not had time to dynamically relax. A further contribution is given by the environment surrounding these structures. Highmass haloes form in the knots of the largescale structure where more matter is available for inflows and mergers, making these objects more disturbed. In a picture where X_{off} possibly traces the cool core–noncool core classification of galaxy clusters in Xrays (Eckert et al. 2011), massive structures have a higher fraction of noncool cores. This influences the highmass tail of the mass function (see Fig. 1). Therefore, given an Xray flux limit, the mass function of a relaxed galaxy clusters sample will be complete to lower masses than an unrelaxed one. This effect also evolves with redshift. At high z, haloes show a larger offset. In this context it means that it is more difficult to detect structures at earlier times. More recently, structures have had time to relax and therefore show smaller values of the offset parameter. This is linked to the development of cool cores at low redshift (Ettori & Brighenti 2008). The PDF shows a powerlaw growth from low offset values (slope α = 3.71 at z = 0) and exponential cutoff at high X_{off}. It is described by a modified Schechter function (Eq. (14)). Its maximum shifts by a factor of ∼1.5, from X_{off} ∼ 0.09 to X_{off} ∼ 0.06, between redshift 1.4 and 0, confirming that haloes have more time to relax. The shape of the distribution does not show a significant redshift trend. The width of the probability density functions, measured at log_{10}P(X_{off}) = − 2.5, at z = 0 and z = 1.4 agree with 2.8% accuracy. At z = 0 these values span from X_{off} ∼ 0.01 to X_{off} ∼ 0.3. This precludes the offset from being a possible assembly bias candidate.
5. Generalized mass function
Generalizations of the mass function have been made in a number of directions: cosmological parameters, angular momentum, and friction (Achitouv & Corasaniti 2012; Achitouv et al. 2014; Del Popolo et al. 2017) via detailed modeling of the collapse barrier. Nevertheless, it is technically demanding to connect these parameters to observations. In this section we generalize the mass function formalism to include additional variables (X_{off},λ) in its formulation. These two quantities describe the properties of dark matter haloes and will hopefully become connectable to observational properties.
5.1. Definitions
We compute the mass function as follows:
Here ΔN_{M} is the number of haloes in each mass bin ΔlnM and V is the total volume of the simulation. It corresponds to Eq. (3). The M(σ) relation and its first derivative are computed with COLOSSUS (Diemer 2018). By convention, the relation between a certain mass and its corresponding scale is normalized with the matter density at z = 0. By combining the previous equations with Eq. (3), we estimate a multiplicity function:
5.2. Generalization
We include the offset parameter and spin by generalizing the approach described in the previous section. The number density of haloes as a function of mass and dynamical state is described by Eq. (17):
Here ΔN_{M, Xoff, λ} is the number of haloes in each mass, X_{off} and λ bin, V is the total volume of the simulated cube and s_{M}, s_{Xoff}, and s_{λ} are the natural (base 10) logarithm of mass (X_{off}, λ) binning. Equivalently to Eq. (16), we calculate
We consider a single integration of h(σ, X_{off}, λ), which results in the set of Eqs. (19). The notation g_{X} designates the marginalization of h(σ, X_{off}, λ) over the variable X:
Integrating again, we obtain
The functions f, g, h are thus linked by derivatives as
where X, Y, Z are permutations of the variables σ, X_{off}, λ. With this method we recover the multiplicity function f(σ), which in this notation is f_{Xoff, λ}(σ). This allows us to study the behavior of the dark matter halo mass function according to different variables, making sure that in the end our analysis provides an accurate multiplicity function.
5.3. Mass–offset–spin function h(σ, X_{off}, λ)
Here we present a model for the generalized mass function h(σ, X_{off}, λ) introduced in the previous subsection. The relaxation state of a dark matter halo is related to the values of X_{off} and λ parameters. We consider HMD, BIGMD, and MDPL2 to build a 3D histogram of halo counts in bins of σ, X_{off}, λ, according to Eqs. (17) and (18).
Mass functions are expressed as multiplicity functions f(σ). It allows the inclusion of part of the redshift evolution in σ. We focus on the highmass end of the mass function, using haloes with M > 2.7 × 10^{13} M_{⊙} h^{−1}.
We measure the halo number density in bins of log_{10}σ^{−1} instead of mass. We consider linear spaced bins from −0.09 to 0.6, with 0.01 width, corresponding to values close to 2.8 × 10^{13} and 1.3 × 10^{16} M_{⊙} h^{−1}. We consider 50 bins spanning logarithmically from 10^{−3.8} to 10^{−0.2} for X_{off} and 50 bins spanning logarithmically from 10^{−4.5} to 10^{−0.1} for λ. This is almost three orders of magnitude higher than the mass resolution in HMD (7.9 × 10^{10} M_{⊙} h^{−1}). Therefore, our results will not be impacted by the mass resolution of the simulations. The total sample consists of 8 051 654 haloes for HMD, 2 103 896 for BigMD and 142 527 for MDPL. First, we estimate directly f(σ) for each simulation, according to Eqs. (15) and (16). Then, we estimate h by computing a 3D histogram in the same mass bins and different X_{off} and λ bins.
Uncertainties on histogram values are computed considering a Poisson number count term and a cosmic variance term:
Here C is a term accounting for cosmic variance, which is set differently according to the type of simulation. Values for the cosmic variance are given in Table 3. These values are estimated by Comparat et al. (2017), using a jackknife method for variance at low masses. We fitted bins containing more than 50 haloes, which according to Eq. (22) means an uncertainty of around 15%; in this way our measurement is not dominated by Poisson uncertainty.
Cosmic variance in different MD simulations.
6. Model
We create a single model for the h(σ, X_{off}, λ) function. We describe in detail its features at z = 0 and its redshift evolution.
6.1. Redshift zero
We consider that both λ and X_{off} probability density functions are described by a modified Schechter function (Eqs. (12) and (14)), and combine these two functions with the multiplicity function along the mass axis. We obtain the following model,
with μ′ = 10^{1.83log10μ}, to disentangle the degeneracy between the two knees of the modified Schechter functions.
This model recalls the Bhattacharya et al. (2011) formulation along the mass axis and considers the combination of a power law and an exponential cutoff (i.e., a modified Schechter function) along the X_{off} and λ axes. The last exponential contains crossed terms between X_{off} and λ, which takes into account their correlation. Both X_{off} and λ modified Schechter functions do not have mass dependency (as suggested by Eqs. (12) and (14)). In the Bhattacharya et al. (2011) formulation, there is a double power law; here we consider a single σ power law. Additional σ dependencies are described by crossed massdependent terms in the exponential cutoff, which relates X_{off} and λ modified Schechter functions. The position of the knee of the X_{off} function (i.e., μ′) is the same in its two exponential cutoffs, but in the second one we introduce the scaling with mass through the parameter e. We correlate directly the knees of the modified Schechter functions for λ and X_{off} (μ, μ′ = 10^{1.83log10μ}), and the slopes of the power law and exponential cutoff of X_{off} (α, 0.05α), to write the model in the most compact way possible.
6.2. Evolution with redshift
Since structures accrete matter with time and grow, the mass function depends on redshift (Springel et al. 2005). Part of this redshift evolution is in the mass–σ relation through the matter power spectrum (Eq. (1)). However, this does not make the mass function completely universal at different times. Tinker et al. (2008) showed that a spherical overdensity mass function evolves up to 30% from z = 0 to z = 2.5. Despali et al. (2016) highlighted how only virial overdensity nears the universality for the mass function. Crocce et al. (2010) considered a FOF mass function and found a 10% evolution up to z = 2. In this work we use ROCKSTAR to identify haloes, its base is a FOF algorithm as well.
Departures from universality are partially explained by the cosmology dependence of the mass function on the power spectrum and growth rate (OndaroMallea et al. 2021). We find that the distribution functions of X_{off} and λ show a redshift trend as well (see Sect. 4, and Eqs. (12) and (14)). We provide a detailed description of the evolution of h(σ, X_{off}, λ) with redshift. Further investigation, using simulations in different cosmologies, is needed to assess a possible relation between parameters in Eq. (24) and cosmological parameters.
Given our model at z = 0, we use the latter as the benchmark model to fit the halo mass−X_{off} − λ function at a higher redshift. For this goal we concatenate again samples from HMD, BIGMD, and MDPL2 at redshifts 0.045, 0.117, 0.221, 0.425, 0.523, 0.702, 0.779, 1.032, and 1.425. We note that BIGMD is not tabulated at exactly the same redshift snapshots as the other two simulations. Nonetheless, we use snapshots that are as close as possible, resulting in a 1.3% difference for the worstcase scenario at z = 1.425. Further details about the snapshots are available in Appendix A. For all these snapshots we consider the same X_{off} and λ binning as we did for z = 0. However, we shift the σ binning slightly upward compared to the z = 0 case, which allows us to reach masses of 7 × 10^{12} M_{⊙} h^{−1} at z = 0.702 and 10^{12} M_{⊙} h^{−1} at z = 1.425.
We include a redshift evolution for all the parameters A, a, q, μ, α, β, γ, δ, and e. We note that we did not consider an evolving critical density contrast with redshift, fixing it at z = 0. So δ_{c} in Eq. (23) is fixed at the value of 1.68647. Considering its evolution, even if tiny, introduces the need for additional evolution of the parameters, as pointed out by Bhattacharya et al. (2011). We model the redshift evolution for these parameter using exponents k_{0}, k_{1}, k_{2}, k_{3}, k_{4}, k_{5}, k_{6}, k_{7}, and k_{8} as follows in Eq. (24):
7. Results
We present the result of the fits to the data (Sect. 3) and the parameters of the model (Sect. 6). We fit directly log_{10}h(σ, X_{off}, λ), which allows better modeling of the highmass end. We consider a Gaussian likelihood
where D is log_{10}h(σ, X_{off}, λ) computed from Eq. (17), M is the base 10 logarithm of the model (Eq. (23)), and E is the uncertainty of log_{10}h(σ, X_{off}, λ) (see log error in Eq. (22)).
7.1. Redshift zero
The bestfit parameters at z = 0 are obtained maximizing the likelihood in Eq. (25). We derive posterior probability distributions and the Bayesian evidence with the nested sampling Monte Carlo algorithm MLFriends (Buchner 2016, 2019), using the UltraNest^{1} software. The results are shown in Fig. A.1. We used flat priors. The description of the parameters, priors, and posteriors is summarized in Table A.4. We obtain log_{10}A = −22.004 ± 0.006, a = 0.885 ± 0.004, q = 2.284 ± 0.016, log_{10}μ= − 3.326 ± 0.001, α = 5.623 ± 0.002, β = −0.391 ± 0.001, γ = 3.024 ± 0.003, δ = 1.209 ± 0.001, e = −1.105 ± 0.005. The parameter a is in agreement with Bhattacharya et al. (2011). Our q parameter shows higher values, but this is expected because an additional mass trend is described by the knee of the exponential cutoff with mixed X_{off}, λ terms. The parameters α and γ describe the powerlaw increment from small X_{off} and λ values, respectively. The second is similar to the values computed for the spin distribution in Sect. 4, while the first is bigger than almost a factor of two. This is expected because an additional offset trend is described by the negative β parameter. Moreover, together with β, the parameter e accounts for the relation between offset and spin, including mass dependency as well. The negative e allows the shifting of the peak along the X_{off} axis to higher values with mass, according to the findings in Sect. 4.
Since this is a 3D model, we show in Fig. 8 all six combinations of h(σ, X_{off}, λ) integrated in 1D (Eq. (19)). To perform the integrals we exploit the Simpson’s rule method for numerical integration^{2}. We show each 2D distribution in five different slices of a single quantity. We recover the typical exponential cutoff along the mass axis, as well as the modified Schechter shapes for offset and spin. Our model describes the σ and λ evolution very well, g_{Xoff}(σ, λ) also agrees with the data in the tails of the distribution. We note the small deviations in the distribution of X_{off} when X_{off} values approach the spatial resolution limit. For lowredshift samples, each version of MultiDark has its own resolution limit: 25 kpc h^{−1} for HMD, 10 kpc h^{−1} for BIGMD, 5 kpc h^{−1} for MDPL2. The mass trend of the low X_{off} slice is slightly underpredicted by the model at low mass. The same holds for the spin trend; there is a 0.1 dex, 3.1σ tension between the peak values of model and data. This is expected within the resolution limit. An improvement toward efficient computation and future generation of Nbody simulations will be needed to probe the kiloparsec scales of Dark Matter haloes in simulation cubes of the gigaparsec scale. All the panels in Fig. 8 involving σ show different uncertainties between succeeding mass bins as a result of the concatenation of bins from different MultiDark versions. MDPL2 bins contain fewer haloes than BIGMD, which contain fewer haloes than HMD. This translates into smaller uncertainties for bins containing a higher number of haloes. Overall, the model provides an excellent representation of the data.
Fig. 8. Single integration of the 3D model. In each panel straight lines indicate the bestfit model, while shaded areas represent the data with 1σ uncertainties. Top left panel: g_{λ}(σ, X_{off}) as a function of X_{off} in different mass slices. Top right panel: g_{λ}(σ, X_{off}) as a function of σ in different X_{off} slices. Middle left panel: g_{Xoff}(σ, λ) as a function of λ in different mass slices. Middle right panel: g_{Xoff}(σ, λ) as a function of σ in different λ slices. Bottom left panel: g_{σ}(X_{off}, λ) as a function of λ in different X_{off} slices. Bottom right panel: g_{σ}(X_{off}, λ) as a function of X_{off} in different λ slices. The integrals are defined in Eq. (19). 
We further integrate the model in Eq. (23), obtaining the distributions of X_{off}f_{σ, λ}(X_{off}) and λf_{σ, Xoff}(λ) (Eq. (20)). Figure 9 shows the result. The distributions around the peaks are well described by the model. This is important because these functions are dominated by objects described by the peak of the PDF. The spin distribution is better behaved than the X_{off} in the tails. This is expected due to spatial resolution limits. Moreover, this is a further confirmation of the fact that X_{off}, λ scatter around their mean values with modified Schechter distributions (see Sect. 4).
Fig. 9. Comparison between data and model of f_{σ, λ}(X_{off}) and f_{σ, Xoff}(λ). Top panels: the straight red lines indicate the integral on the bestfit model, while the shaded blue areas represent the integral on the 3d h(σ, X_{off}, λ) data with 1σ uncertainties. Each bottom panel shows the residual trend with σ error; the straight black line represents the perfect match between data and model with null residual. Top left panel: f(X_{off}) as a function of X_{off}. Bottom left panel: residual between f_{σ, λ}(X_{off}) data and model in logarithmic scale. Top right panel: f(λ) as a function of λ. Bottom right panel: residual between f_{σ, Xoff}(λ) data, and model in logarithmic scale. 
We obtain the multiplicity function f_{Xoff, λ}(σ) marginalizing h(σ, X_{off}, λ) on X_{off}, λ (i.e., performing the double integral):
The result is shown in Fig. 10. In the top panel, we show a comparison between our model, the data obtained from simulations, and the Comparat et al. (2017) model, fitted on these same simulations at z = 0. The multiplicity functions computed directly on each simulation cube, without taking X_{off} and λ into account, are shown by three shaded regions in different colors: red for MDPL2, green for BIGMD, and orange for HMD. Bigger simulation boxes extend to higher mass values. The light blue shaded region represents the 2D integral computed on the concatenated sample of all three simulations. The solid blue line is the integral of our model and the dashed pink line is the Comparat et al. (2017) model. In the lower panel, we show the percentage difference between the multiplicity function f(σ) we recover and the Comparat et al. (2017) model, obtained on the same MultiDark simulations. This difference is always under 3.3% in the mass range of interest. It is also compatible with uncertainty on the data. Our model is able to recover the halo mass distribution in the simulations, with the advantage of taking into account parameters that describe the dynamical state as well. Once again, we note that our model is adapted to masses higher than 2.7 × 10^{13} M_{⊙} h^{−1} at z = 0 (10^{12} M_{⊙} h^{−1} at z = 1.4).
Fig. 10. Comparison of multiplicity functions. Top panel: the three shaded regions show the 1σ contours of f(σ) data directly computed on different simulations (orange for HMD, green for BIGMD, red for MDPL2). The light blue shaded region is the 1σ contour of the 2D integral computed on the concatenated sample containing all three simulations; the dashed pink line indicates the mass function from Comparat et al. (2017); the blue solid line is the f(σ) recovered by integrating our model along X_{off} and λ. Bottom panel: the blue thick line is the fractional difference between our f(σ) and that of Comparat et al. (2017). The light blue shaded area denotes the 1σ contours of the residual between the integrated data and our bestfit model; the black horizontal line indicates the perfect match with null residual. 
7.2. Evolution with redshift
To study the redshift evolution we start from the fiducial model at z = 0. We add the redshift dependence (Eq. (24)) to the bestfit parameters in Eq. (23) at z = 0. We concatenate samples for ten redshift values, as described in Sect. 6.2.
We obtain the values of the exponents in Eq. (24) fitting the z trend of each parameter for all the concatenated snapshots simultaneously. We obtain k_{0} = −0.0441 ± 0.0001, k_{1} = −0.161 ± 0.001, k_{2} = 0.041 ± 0.002, k_{3} = −0.1286 ± 0.0002, k_{4} = 0.1081 ± 0.0002, k_{5} = −0.311 ± 0.001, k_{6} = 0.0902 ± 0.0004, k_{7} = −0.0768 ± 0.0004, and k_{8} = 0.612 ± 0.002. The full result is shown by the triangular plot in Fig. A.2. Priors and posteriors for each parameter are given in Table A.5.
Redshift dependence is shown in Fig. 11. The shaded areas include the uncertainty on the bestfit parameter at z = 0 and on the z evolution, according to Eq. (26)
Fig. 11. Redshift evolution of the bestfit parameters of our model. Each panel shows a single parameter. The values at z = 0 are reported in Table A.4. The redshift evolution is described by Eq. (24); the slopes are given in Table A.5. 
where P is each parameter in Eq. (23), P_{0} is its value at z = 0, and k indicates each parameter describing the evolution in Eq. (24).
The parameters A, q, μ, α, γ, and β show an increasing redshift trend. On the other hand a, β, and e decrease with redshift. This means that with increasing redshift the modified Schechter functions need to increase more quickly and decrease more slowly. The knee describing the mass trend (parameter a) decreases with redshift, in agreement with Bhattacharya et al. (2011). They find no redshift dependence for the slope of mass trend (parameter q). This is not true for this work, where q increases with z. However, this is mitigated by the mass trend of the position of the knee in the crossed X_{off}, λ exponential cutoff, which decreases at early times. The fact that the position of the X_{off} knee (parameter μ) moves to higher values at high z confirms the results of Sect. 4, with the higher average value of X_{off} early in time (Figs. 6 and 7).
8. Summary and conclusions
In the context of the hierarchical model of structure formation, the evolution of the number density of galaxy clusters is a powerful cosmological probe. In order to achieve precision cosmology with the next generation of galaxy clusters samples (such as the eROSITA All Sky Survey; Merloni et al. 2012), precise modeling of the theoretical mass function is necessary. We calibrated a model that includes a dynamical description of dark matter haloes. Using the formalism described here and the MultiDark simulations, we quantified the impact on the mass function of the lack of unrelaxed structures (see Fig. 1). We explored relations between quantities that describe different aspects of dark matter haloes, including their dynamical state. We investigated the concentrationmass relation. We confirmed the recent discovery of the concentration upturn at high masses, in agreement with previous results from Prada et al. (2012) and Klypin et al. (2016), based on a similar set of MultiDark simulations. In addition, our model provides a prediction of concentration according to mass and redshift with one single equation (Eq. (9)). The probability density function of concentration is a modified Schechter law, with mass dependency (Eq. (10)). We find that the concentration of lowmass haloes has a faster redshift evolution than highmass objects, especially in the highconcentration regime. For concentration c = 8, the PDF for highmass haloes shifts by 0.32 dex from z = 1.43 to z = 0, while the lowmass value changes by 0.49 dex. We find the spin parameter λ to be modeled by a linear relation with mass and a probability density function well described by a modified Schechter function (Eq. (12)), in agreement with RodriguezPuebla et al. (2016). The offset parameter evolves with mass and redshift according to Eq. (13). The negative slope of the relation suggests that lowmass haloes are typically more relaxed compared to highmass objects. This is true at every redshift. The offset distribution around the mean value is well described by a modified Schechter function (Eq. (14)). The peak of the distribution shifts by a factor of 1.5 between z ∼ 1.4 and z = 0. This is in agreement with haloes relaxing with time and the recent formation of cool cores in galaxy clusters.
We define a general mass function framework, where dark matter haloes are not only described as a function of mass but also by the two additional variables X_{off}, λ. This approach considers the mass, offset parameter, and spin of each halo at the same time in a σ − X_{off} − λ function (Eq. (18)). We model it in Sect. 6 combining terms of a fiducial mass function (Bhattacharya et al. 2011) with modified Schechter functions for X_{off}, λ, as obtained in Sect. 4. This new approach accounts for the dynamical state of dark matter haloes directly in the context of the halo mass function, providing 2D and 1D distributions at the same time.
We describe the fitting procedure and results in Sect. 7. Our result at z = 0 recovers the Comparat et al. (2017) mass function, which is fitted on the same set of simulations, with 3.3% accuracy. This means that our model is able to account for the dynamical state of dark matter haloes simultaneously with mass and to describe the multiplicity function with great precision. In addition, our model includes the redshift evolution, according to Eq. (24). The result is shown in Fig. 11.
The link with observations will be explored further in future work.
Acknowledgments
We thank the anonymous referee for the constructive feedback. JC thanks Dominique Eckert for insightful discussions about this project. The CosmoSim database used in this paper is a service by the LeibnizInstitute for Astrophysics Potsdam (AIP). The MultiDark database was developed in cooperation with the Spanish MultiDark Consolider Project CSD200900064. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gausscentre.eu) and the Partnership for Advanced Supercomputing in Europe (PRACE, www.praceri.eu) for funding the MultiDark simulation project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (LRZ, www.lrz.de). We derive posterior probability distributions and the Bayesian evidence with the nested sampling Monte Carlo algorithm MLFriends (Buchner 2016, 2019) using the UltraNest software. To fit relations and distributions in Sect. 4, we use a χ^{2} minimization algorithm with CURVE_FIT python package (Virtanen et al. 2020) (https://scipy.org/). In order to perform cosmology calculation, we used the python toolkit COLOSSUS (Diemer 2018) (https://bdiemer.bitbucket.io/colossus/).
References
 Abazajian, K., Addison, G., Adshead, P., et al. 2019, ArXiv eprints [arXiv:1907.04473] [Google Scholar]
 Achitouv, I. E., & Corasaniti, P. S. 2012, JCAP, 2012, 002 [Google Scholar]
 Achitouv, I., Wagner, C., Weller, J., & Rasera, Y. 2014, JCAP, 2014, 077 [CrossRef] [Google Scholar]
 Allen, S. W., Evrard, A. E., & Mantz, A. B. 2011, ARA&A, 49, 409 [NASA ADS] [CrossRef] [Google Scholar]
 Angulo, R. E., Baugh, C. M., & Lacey, C. G. 2008, MNRAS, 387, 921 [Google Scholar]
 Angulo, R. E., Springel, V., White, S. D. M., et al. 2012, MNRAS, 426, 2046 [CrossRef] [Google Scholar]
 Baldi, A. S., De Petris, M., Sembolini, F., et al. 2018, MNRAS, 479, 4028 [Google Scholar]
 Baldi, A. S., De Petris, M., Sembolini, F., et al. 2019, J. Phys. Conf. Ser., 1226, 012003 [Google Scholar]
 Bartelmann, M. 2010, Class. Quant. Grav., 27, 233001 [Google Scholar]
 Baxter, E. J., Sherwin, B. D., & Raghunathan, S. 2019, JCAP, 2019, 001 [Google Scholar]
 Behroozi, P., Wechsler, R., & Wu, H.Y. 2013, ApJ, 762, 109 [Google Scholar]
 Behroozi, P., Knebe, A., Pearce, F. R., et al. 2015, MNRAS, 454, 3020 [Google Scholar]
 Benson, B. A., Ade, P. A. R., Ahmed, Z., et al. 2014, in SPT3G: A NextGeneration Cosmic Microwave Background Polarization Experiment on the South Pole Telescope, SPIE Conf. Ser., 9153, 91531P [Google Scholar]
 Bett, P., Eke, V., Frenk, C. S., et al. 2007, MNRAS, 376, 215 [NASA ADS] [CrossRef] [Google Scholar]
 Bhattacharya, S., Heitmann, K., White, M., et al. 2011, ApJ, 732, 122 [Google Scholar]
 Bianconi, M., Ettori, S., & Nipoti, C. 2013, MNRAS, 434, 1565 [Google Scholar]
 Bilton, L. E., Hunt, M., Pimbblet, K. A., & Roediger, E. 2019, MNRAS, 490, 5017 [Google Scholar]
 Bocquet, S., Saro, A., Dolag, K., & Mohr, J. J. 2016, MNRAS, 456, 2361 [Google Scholar]
 Bocquet, S., Heitmann, K., & Habib, S. 2020, ApJ, 901, 5 [Google Scholar]
 Bond, J. R., Cole, S., Efstathiou, G., & Kaiser, N. 1991, ApJ, 379, 440 [NASA ADS] [CrossRef] [Google Scholar]
 Bryan, G. L., & Norman, M. L. 1998, ApJ, 495, 80 [NASA ADS] [CrossRef] [Google Scholar]
 Buchner, J. 2016, Stat. Comput., 26, 383 [CrossRef] [Google Scholar]
 Buchner, J. 2019, PASP, 131, 108005 [Google Scholar]
 Cibirka, N., Cypriano, E. S., Brimioulle, F., et al. 2017, MNRAS, 468, 1092 [Google Scholar]
 Comparat, J., Prada, F., Yepes, G., & Klypin, A. 2017, MNRAS, 469, 4157 [Google Scholar]
 Contreras, S., Zehavi, I., Padilla, N., et al. 2019, MNRAS, 484, 1133 [Google Scholar]
 Cooray, A., & Chen, X. 2002, ApJ, 573, 43 [Google Scholar]
 Correa, C. A., Wyithe, J. S. B., Schaye, J., & Duffy, A. R. 2015, MNRAS, 452, 1217 [NASA ADS] [CrossRef] [Google Scholar]
 Crocce, M., Fosalba, P., Castand er, F. J., & Gaztañaga, E. 2010, MNRAS, 403, 1353 [Google Scholar]
 Croton, D. J., Gao, L., & White, S. D. M. 2007, MNRAS, 374, 1303 [Google Scholar]
 Dawson, K. S., Kneib, J.P., Percival, W. J., et al. 2016, AJ, 151, 44 [Google Scholar]
 de Jong, R. 2011, The Messenger, 145, 14 [NASA ADS] [Google Scholar]
 Del Popolo, A., Pace, F., & Le Delliou, M. 2017, JCAP, 3, 032 [Google Scholar]
 DESI Collaboration (Aghamousa, A., et al.) 2016, ArXiv eprints [arXiv:1611.00036] [Google Scholar]
 Despali, G., Giocoli, C., Angulo, R. E., et al. 2016, MNRAS, 456, 2486 [NASA ADS] [CrossRef] [Google Scholar]
 Diemer, B. 2018, ApJS, 239, 35 [Google Scholar]
 Diemer, B., & Joyce, M. 2019, ApJ, 871, 168 [Google Scholar]
 Diemer, B., & Kravtsov, A. V. 2015, ApJ, 799, 108 [NASA ADS] [CrossRef] [Google Scholar]
 Du, W., Fan, Z., Shan, H., et al. 2015, ApJ, 814, 120 [Google Scholar]
 Duffy, A. R., Schaye, J., Kay, S. T., & Dalla Vecchia, C. 2008, MNRAS, 390, L64 [NASA ADS] [CrossRef] [Google Scholar]
 Dutton, A. A., & Macciò, A. V. 2014, MNRAS, 441, 3359 [Google Scholar]
 Eckert, D., Molendi, S., & Paltani, S. 2011, A&A, 526, A79 [Google Scholar]
 Eckert, D., Ghirardini, V., Ettori, S., et al. 2019, A&A, 621, A40 [Google Scholar]
 Ettori, S., & Brighenti, F. 2008, MNRAS, 387, 631 [Google Scholar]
 Finoguenov, A., Merloni, A., Comparat, J., et al. 2019, The Messenger, 175, 39 [NASA ADS] [Google Scholar]
 Finoguenov, A., Rykoff, E., Clerc, N., et al. 2020, A&A, 638, A114 [Google Scholar]
 Foëx, G., Motta, V., Jullo, E., Limousin, M., & Verdugo, T. 2014, A&A, 572, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gao, L., Springel, V., & White, S. D. M. 2005, MNRAS, 363, L66 [Google Scholar]
 Giocoli, C., Tormen, G., & Sheth, R. K. 2012, MNRAS, 422, 185 [NASA ADS] [CrossRef] [Google Scholar]
 Heitmann, K., Frontiere, N., Sewell, C., et al. 2015, ApJS, 219, 34 [Google Scholar]
 Henson, M. A., Barnes, D. J., Kay, S. T., McCarthy, I. G., & Schaye, J. 2017, MNRAS, 465, 3361 [Google Scholar]
 Hitomi Collaboration (Aharonian, F., et al.) 2018, PASJ, 70, 9 [NASA ADS] [Google Scholar]
 Hollowood, D. L., Jeltema, T., Chen, X., et al. 2019, ApJS, 244, 22 [Google Scholar]
 Hwang, H. S., & Lee, M. G. 2007, ApJ, 662, 236 [Google Scholar]
 Ider Chitham, J., Comparat, J., Finoguenov, A., et al. 2020, MNRAS, 499, 4768 [Google Scholar]
 Ishiyama, T., Enoki, M., Kobayashi, M. A. R., et al. 2015, PASJ, 67, 61 [Google Scholar]
 Jenkins, A., Frenk, C. S., White, S. D. M., et al. 2001, MNRAS, 321, 372 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Käfer, F., Finoguenov, A., Eckert, D., et al. 2019, A&A, 628, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Käfer, F., Finoguenov, A., Eckert, D., et al. 2020, A&A, 634, A8 [CrossRef] [EDP Sciences] [Google Scholar]
 Kilbinger, M. 2015, Rep. Prog. Phys., 78, 086901 [Google Scholar]
 Klypin, A., Yepes, G., Gottlöber, S., Prada, F., & Heß, S. 2016, MNRAS, 457, 4340 [NASA ADS] [CrossRef] [Google Scholar]
 Knebe, A., Knollmann, S. R., Muldrew, S. I., et al. 2011, MNRAS, 415, 2293 [NASA ADS] [CrossRef] [Google Scholar]
 Knebe, A., Pearce, F. R., Lux, H., et al. 2013, MNRAS, 435, 1618 [Google Scholar]
 Kravtsov, A. V., Klypin, A. A., & Khokhlov, A. M. 1997, ApJS, 111, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Lang, M., HolleyBockelmann, K., & Sinha, M. 2015, ApJ, 811, 152 [Google Scholar]
 Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, ArXiv eprints [arXiv:1110.3193] [Google Scholar]
 LSST Science Collaboration (Abell, P. A., et al.) 2009, ArXiv eprints [arXiv:0912.0201] [Google Scholar]
 Ludlow, A. D., Navarro, J. F., Li, M., et al. 2012, MNRAS, 427, 1322 [Google Scholar]
 Ludlow, A. D., Navarro, J. F., BoylanKolchin, M., et al. 2013, MNRAS, 432, 1103 [NASA ADS] [CrossRef] [Google Scholar]
 Ludlow, A. D., Navarro, J. F., Angulo, R. E., et al. 2014, MNRAS, 441, 378 [NASA ADS] [CrossRef] [Google Scholar]
 Ludlow, A. D., Bose, S., Angulo, R. E., et al. 2016, MNRAS, 460, 1214 [Google Scholar]
 Macciò, A. V., Dutton, A. A., & van den Bosch, F. C. 2008, MNRAS, 391, 1940 [NASA ADS] [CrossRef] [Google Scholar]
 Manolopoulou, M., & Plionis, M. 2017, MNRAS, 465, 2616 [Google Scholar]
 McClintock, T., Rozo, E., Becker, M. R., et al. 2019, ApJ, 872, 53 [CrossRef] [Google Scholar]
 Meneghetti, M., & Rasia, E. 2013, ArXiv eprints [arXiv:1303.6158] [Google Scholar]
 Merloni, A., Predehl, P., Becker, W., et al. 2012, ArXiv eprints [arXiv:1209.3114] [Google Scholar]
 Murray, S. G., Power, C., & Robotham, A. S. G. 2013, Astron. Comput., 3, 23 [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563 [Google Scholar]
 Neto, A. F., Gao, L., Bett, P., et al. 2007, MNRAS, 381, 1450 [NASA ADS] [CrossRef] [Google Scholar]
 Nishimichi, T., Takada, M., Takahashi, R., et al. 2019, ApJ, 884, 29 [Google Scholar]
 OndaroMallea, L., Angulo, R. E., Zennaro, M., Contreras, S., & Aricò, G. 2021, ArXiv eprints [arXiv:2102.08958] [Google Scholar]
 Peebles, P. J. E. 1969, ApJ, 155, 393 [NASA ADS] [CrossRef] [Google Scholar]
 Phriksee, A., Jullo, E., Limousin, M., et al. 2020, MNRAS, 491, 1643 [Google Scholar]
 Planck Collaboration XI. 2014, A&A, 571, A11 [Google Scholar]
 PovedaRuiz, C. N., ForeroRomero, J. E., & MuñozCuartas, J. C. 2016, ApJ, 832, 169 [Google Scholar]
 Prada, F., Klypin, A. A., Cuesta, A. J., BetancortRijo, J. E., & Primack, J. 2012, MNRAS, 423, 3018 [Google Scholar]
 Predehl, P., Andritschke, R., Arefiev, V., et al. 2021, A&A, 647, A1 [EDP Sciences] [Google Scholar]
 Press, W. H., & Schechter, P. 1974, ApJ, 187, 425 [NASA ADS] [CrossRef] [Google Scholar]
 Ragagnin, A., Dolag, K., Moscardini, L., Biviano, A., & D’Onofrio, M. 2019, MNRAS, 486, 4001 [NASA ADS] [CrossRef] [Google Scholar]
 Rephaeli, Y. 1995, ARA&A, 33, 541 [Google Scholar]
 Riebe, K., Partl, A. M., Enke, H., et al. 2013, Astron. Nachr., 334, 691 [Google Scholar]
 RodriguezPuebla, A., Behroozi, P., Primack, J., et al. 2016, MNRAS, 462, 893 [Google Scholar]
 Salvati, L., Douspis, M., & Aghanim, N. 2020, A&A, 643, A20 [Google Scholar]
 Sereno, M., Giocoli, C., Ettori, S., & Moscardini, L. 2015, MNRAS, 449, 2024 [NASA ADS] [CrossRef] [Google Scholar]
 Shan, H., Kneib, J.P., Li, R., et al. 2017, ApJ, 840, 104 [Google Scholar]
 Sheth, R. K., & Tormen, G. 1999, MNRAS, 308, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Sheth, R. K., & Tormen, G. 2002, MNRAS, 329, 61 [Google Scholar]
 Skillman, S. W., Warren, M. S., Turk, M. J., et al. 2014, ArXiv eprints [arXiv:1407.2600] [Google Scholar]
 Song, H., Hwang, H. S., Park, C., Smith, R., & Einasto, M. 2018, ApJ, 869, 124 [Google Scholar]
 Spergel, D., Gehrels, N., Baltay, C., et al. 2015, ArXiv eprints [arXiv:1503.03757] [Google Scholar]
 Springel, V. 2005, MNRAS, 364, 1105 [Google Scholar]
 Springel, V., White, S. D. M., Jenkins, A., et al. 2005, Nature, 435, 629 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Sunyaev, R. A., Norman, M. L., & Bryan, G. L. 2003, Astron. Lett., 29, 783 [Google Scholar]
 Thomas, P. A., Muanwong, O., Pearce, F. R., et al. 2001, MNRAS, 324, 450 [Google Scholar]
 Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709 [NASA ADS] [CrossRef] [Google Scholar]
 Tovmassian, H. M. 2015, Astrophysics, 58, 328 [Google Scholar]
 Umetsu, K. 2020, A&ARv, 28, 1 [Google Scholar]
 van Uitert, E., Gilbank, D. G., Hoekstra, H., et al. 2016, A&A, 586, A43 [Google Scholar]
 Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
 Wang, J., Bose, S., Frenk, C. S., et al. 2020, Nature, 585, 39 [Google Scholar]
 Weinberg, D. H., Mortonson, M. J., Eisenstein, D. J., et al. 2013, Phys. Rep, 530, 87 [Google Scholar]
 Zhao, D. H., Mo, H. J., Jing, Y. P., & Börner, G. 2003, MNRAS, 339, 12 [Google Scholar]
 Zhao, D. H., Jing, Y. P., Mo, H. J., & Börner, G. 2009, ApJ, 707, 354 [Google Scholar]
 Zubeldia, Í., & Challinor, A. 2019, MNRAS, 489, 401 [Google Scholar]
Appendix A: Figures and tables
In this section we collect the figures and tables relative to this work. They describe the mean relations between concentration, offset parameter, spin, and mass, as well as the full probability density functions of these quantities. Moreover, we show additional plots describing the halo σ − X_{off} − λ function.
Bestfit parameters for concentration–σ relation and its PDF P(c).
Bestfit parameters for λ–σ relation and its PDF P(λ) at different redshifts.
Bestfit parameters for X_{off}σ relation and its PDF P(X_{off}).
Model parameters with priors and posterior constraints at redshift zero.
Model parameters with prior and posterior constraints for the redshift evolution of the halo σ − X_{off} − λ function.
Full list of snapshots used, available in HMD, BIGMD, MDPL2.
Fig. A.1. Marginalized posterior distributions of the bestfit parameters of the halo σ − X_{off} − λ function. The 0.68 and 0.95 confidence levels of the posteriors are shown as filled 2D contours. The 2.5th, 16th, 84th, and 97.5th percentiles of the onedimensional posterior distributions are indicated by the vertical lines on the diagonal plots. The model is given by Eq. (23). The parameters are also given in Table A.4. 
Fig. A.2. Marginalized posterior distributions of the bestfit parameters describing the redshift evolution of h(σ, X_{off}, λ). The 0.68 and 0.95 confidence levels of the posteriors are shown as filled 2D contours. The model is given by Eq. (24). The parameters are reported in Table A.5. 
Appendix B: Offset in physical units
Here we present the results of the same analysis elaborated in Sects. 4 and 7 to the offset parameter in physical units X_{off, P}, measured in kpc h^{−1}. This approach allows the comparison between the physics of dark matter simulations and observations.
B.1. Offset–mass–redshift relation
The relation between X_{off, P}, mass and redshift is modeled by
The distribution of X_{off, P} around its mean value is described by a modified Schechter function, but X_{off, P} is not normalized to the virial radius. Therefore, a mass dependence has to be included in the relation.
The bestfit parameters are given in Table B.1, and the results are shown in Figs. B.1 and B.2.
Fig. B.1. X_{off, P}–σ relation (Eq. (B.1)). Circular dots, triangles, and squares represent HMD, BIGMD, and MDPL2, respectively. They are colorcoded by redshift. Straight lines show the bestfit model. Parameters are given in Table B.1. 
Fig. B.2. Probability density function of X_{off, P} (Eq. (B.2)). Each panel shows the distribution at a specific redshift. Each set is divided in mass slices, identified by color. Scatter points indicate the data, while straight lines represent the modified Schechter model (14). For clarity, each line and its fit are shifted by 0.2 dex along both axes. This means that the coefficient C_{0} assumes values ( + 0.6, +0.4, +0.2, 0.0, −0.2, −0.4), while C_{1} is ( − 0.6, −0.4, −0.2, 0.0, +0.2, +0.4). The red line is not shifted, and thus is the one with the correct normalization. 
Bestfit parameters for X_{off, P}–σ relation and P(X_{off, P}).
B.2. Mass – X_{off, P} – λ function
In this section we collect figures and tables that describe h(σ, X_{off, P}, λ). The analysis is similar to h(σ, X_{off}, λ) explained in Sects. 5 and 6; the only difference is that the modified Schechter function that describes X_{off, P} needs a massdependent term. Therefore, we introduce an additional parameter and model the distribution according to
We recover the fiducial mass function at the ∼3.9% level.
Fig. B.3. Single integration of the 3D model. In each panel straight lines indicate the bestfit model, while shaded areas represent the data with 1σ uncertainties. Top left panel: g_{λ}(σ, X_{off, P}) as a function of X_{off, P} in different mass slices. Top right panel: g_{λ}(σ, X_{off, P}) as a function of σ in different X_{off, P} slices. Middle left panel: g_{Xoff, P}(σ, λ) as a function of λ in different mass slices. Middle right panel: g_{Xoff, P}(σ, λ) as a function of σ in different λ slices. Bottom left panel: g_{σ}(X_{off, P}, λ) as a function of λ in different X_{off, P} slices. Bottom right panel: g_{σ}(X_{off, P}, λ) as a function of X_{off, P} in different λ slices. 
Fig. B.4. Comparison of f(X_{off, P}) and f(λ) between data and model. Top panels: solid red lines indicate the integral on the bestfit model, while the shaded blue areas represent the integral on the 3D h(σ, X_{off, P}, λ) data with 1σ uncertainties. Each bottom panel shows the residual trend with σ error; the straight black line represents the perfect match between data and model, with null residual. Top left panel: f(X_{off, P}) as a function of X_{off, P}. Bottom left panel: residual between f(X_{off, P}) data and model in logarithmic scale. Top right panel: f(λ) as a function of λ. Bottom right panel: residual between f(λ) data and model in logarithmic scale. 
Fig. B.5. Comparison of multiplicty functions. Top panel: three shaded regions show the 1σ contours of f(σ) data directly computed on different simulations (orange for HMD, green for BIGMD, red for MDPL2); the light blue shaded region is the 1σ contour of the 2D integral computed on the concatenated sample containing all three simulations; the dashed pink line indicates the mass function from Comparat et al. (2017), while the blue solid line is the f(σ) we recover integrating our model along X_{off, P} and λ. Bottom panel: blue thick line is the fractional difference between our f(σ) and the Comparat et al. (2017) value. The light blue shaded area represent the 1σ contours of the residual between the integrated data and our bestfit model; the black horizontal line indicates the perfect match with null residual. 
Fig. B.6. Redshift evolution of the bestfit parameters for h(σ, X_{off, P}, λ). Each panel shows a single parameter. The values at z = 0 are reported in Table B.2. The slopes of the redshift trends are given in Table B.3. 
Fig. B.7. Marginalized posterior distributions of the h(σ, X_{off, P}, λ) bestfit parameters at redshift 0. The 0.68 and 0.95 confidence levels of the posteriors are shown as filled 2D contours. The 2.5th, 16th, 84th, and 97.5th percentiles of the onedimensional posterior distributions are indicated by the vertical lines on the diagonal plots. 
Fig. B.8. Marginalized posterior distributions of the bestfit parameters describing the redshift evolution of h(σ, X_{off, P}, λ). The 0.68 and 0.95 confidence levels of the posteriors are shown as filled 2D contours. The 2.5th, 16th, 84th, and 97.5th percentiles of the onedimensional posterior distributions are indicated by the vertical lines on the diagonal plots. 
h(σ, X_{off, P}, λ) model parameters with priors and posterior costraints.
Model parameters with prior and posterior constraints for the redshift evolution of h(σ, X_{off, P}, λ).
All Tables
Correspondence between mass, peak height, and variance of the linear density field at z = 0 and z = 0.5.
Model parameters with prior and posterior constraints for the redshift evolution of the halo σ − X_{off} − λ function.
Model parameters with prior and posterior constraints for the redshift evolution of h(σ, X_{off, P}, λ).
All Figures
Fig. 1. Contribution of relaxed and unrelaxed haloes to the total mass function. Left panel: distribution of redshift zero dark matter haloes in the X_{off, P} and λ plane. Cuts in X_{off, P} and λ are applied to divide relaxed (blue) and disturbed (orange) structures. The contours contain 1%, 5%, 10%, 30%, and 50% of the data. Right panel: Halo mass functions vs. mass (σ) at redshift 0, as defined in Sect. 2 and using Eqs. (15) and (16). This mass function is built with different subsets of haloes from HMD at z = 0 (see Sect. 3): the red line indicates the model from Comparat et al. (2017), the shaded areas represent the relaxed (blue), unrelaxed (orange), and the full sample (green) of haloes. The areas cover 1σ uncertainties. Lower panel: residuals fraction of each component compared to the red model in the upper panel. 

In the text 
Fig. 2. Concentration–σ relation (Eq. (9)). Circular dots, triangles, and squares represent HMD, BIGMD, MDPL2, respectively. They are colorcoded by redshift. Straight lines indicate our bestfit model, dotted lines show the model from Klypin et al. (2016), and the shaded blue area indicates the distribution from Wang et al. (2020) at z = 0. The upper xaxis converts peaks values into mass at z = 0. 

In the text 
Fig. 3. Probability density function of the concentration (Eq. (10)) at the different redshifts values given in each panel. Each set is divided into mass slices, and is colorcoded accordingly. The shaded areas represent the data with 1σ error, while the straight lines indicate the bestfit model. The blue points and line represent the total sample not sliced in mass. For clarity, each line and its fit is shifted by 0.1 dex on the yaxis. This means that the constant C_{0} assumes values of (+0.3, +0.2, +0.1, 0.0, −0.1, −0.2). The purple line is not shifted, and is thus the one with the proper normalization. 

In the text 
Fig. 4. Spin–mass relation (λσ, Eq. (11)). The model is a linear relation with no redshift trend. Data points are colorcoded by redshift, while the different geometrical shapes refer to different simulations: squares for HMD, triangles for BigMD, and circles for MDPL. The straight line indicates the bestfit model, which considers all simulations and redshift at the same time. The bestfit parameters are given in Table A.2. The upper xaxis converts peaks values into mass at z = 0. 

In the text 
Fig. 5. PDF of λ (Eq. (12)). Individual points represent spin bins used to compute the distribution; straight lines refer to the bestfit modified Schechter model. They are colorcoded by redshift. We do not consider mass dependence in this relation. Each redshift slice is fitted independently. The bestfit parameters are given in Table A.2. 

In the text 
Fig. 6. Offset–mass relation (X_{off}–σ, Eq. (13)). Circular dots, triangles, and squares represent HMD, BIGMD, and MDPL2 respectively. They are colorcoded by redshift. Straight lines show the bestfit model. The bestfit parameters are given in Table A.3. The upper xaxis converts peaks values into mass at z = 0. 

In the text 
Fig. 7. Probability density function of X_{off}. Each panel shows the distribution at a specific redshift. Scatter points indicate the data, while straight lines represent the modified Schechter model. The samples are colorcoded by redshift. Each redshift slice is fitted independently by Eq. (14). The bestfit parameters are given in Table A.3. 

In the text 
Fig. 8. Single integration of the 3D model. In each panel straight lines indicate the bestfit model, while shaded areas represent the data with 1σ uncertainties. Top left panel: g_{λ}(σ, X_{off}) as a function of X_{off} in different mass slices. Top right panel: g_{λ}(σ, X_{off}) as a function of σ in different X_{off} slices. Middle left panel: g_{Xoff}(σ, λ) as a function of λ in different mass slices. Middle right panel: g_{Xoff}(σ, λ) as a function of σ in different λ slices. Bottom left panel: g_{σ}(X_{off}, λ) as a function of λ in different X_{off} slices. Bottom right panel: g_{σ}(X_{off}, λ) as a function of X_{off} in different λ slices. The integrals are defined in Eq. (19). 

In the text 
Fig. 9. Comparison between data and model of f_{σ, λ}(X_{off}) and f_{σ, Xoff}(λ). Top panels: the straight red lines indicate the integral on the bestfit model, while the shaded blue areas represent the integral on the 3d h(σ, X_{off}, λ) data with 1σ uncertainties. Each bottom panel shows the residual trend with σ error; the straight black line represents the perfect match between data and model with null residual. Top left panel: f(X_{off}) as a function of X_{off}. Bottom left panel: residual between f_{σ, λ}(X_{off}) data and model in logarithmic scale. Top right panel: f(λ) as a function of λ. Bottom right panel: residual between f_{σ, Xoff}(λ) data, and model in logarithmic scale. 

In the text 
Fig. 10. Comparison of multiplicity functions. Top panel: the three shaded regions show the 1σ contours of f(σ) data directly computed on different simulations (orange for HMD, green for BIGMD, red for MDPL2). The light blue shaded region is the 1σ contour of the 2D integral computed on the concatenated sample containing all three simulations; the dashed pink line indicates the mass function from Comparat et al. (2017); the blue solid line is the f(σ) recovered by integrating our model along X_{off} and λ. Bottom panel: the blue thick line is the fractional difference between our f(σ) and that of Comparat et al. (2017). The light blue shaded area denotes the 1σ contours of the residual between the integrated data and our bestfit model; the black horizontal line indicates the perfect match with null residual. 

In the text 
Fig. 11. Redshift evolution of the bestfit parameters of our model. Each panel shows a single parameter. The values at z = 0 are reported in Table A.4. The redshift evolution is described by Eq. (24); the slopes are given in Table A.5. 

In the text 
Fig. A.1. Marginalized posterior distributions of the bestfit parameters of the halo σ − X_{off} − λ function. The 0.68 and 0.95 confidence levels of the posteriors are shown as filled 2D contours. The 2.5th, 16th, 84th, and 97.5th percentiles of the onedimensional posterior distributions are indicated by the vertical lines on the diagonal plots. The model is given by Eq. (23). The parameters are also given in Table A.4. 

In the text 
Fig. A.2. Marginalized posterior distributions of the bestfit parameters describing the redshift evolution of h(σ, X_{off}, λ). The 0.68 and 0.95 confidence levels of the posteriors are shown as filled 2D contours. The model is given by Eq. (24). The parameters are reported in Table A.5. 

In the text 
Fig. B.1. X_{off, P}–σ relation (Eq. (B.1)). Circular dots, triangles, and squares represent HMD, BIGMD, and MDPL2, respectively. They are colorcoded by redshift. Straight lines show the bestfit model. Parameters are given in Table B.1. 

In the text 
Fig. B.2. Probability density function of X_{off, P} (Eq. (B.2)). Each panel shows the distribution at a specific redshift. Each set is divided in mass slices, identified by color. Scatter points indicate the data, while straight lines represent the modified Schechter model (14). For clarity, each line and its fit are shifted by 0.2 dex along both axes. This means that the coefficient C_{0} assumes values ( + 0.6, +0.4, +0.2, 0.0, −0.2, −0.4), while C_{1} is ( − 0.6, −0.4, −0.2, 0.0, +0.2, +0.4). The red line is not shifted, and thus is the one with the correct normalization. 

In the text 
Fig. B.3. Single integration of the 3D model. In each panel straight lines indicate the bestfit model, while shaded areas represent the data with 1σ uncertainties. Top left panel: g_{λ}(σ, X_{off, P}) as a function of X_{off, P} in different mass slices. Top right panel: g_{λ}(σ, X_{off, P}) as a function of σ in different X_{off, P} slices. Middle left panel: g_{Xoff, P}(σ, λ) as a function of λ in different mass slices. Middle right panel: g_{Xoff, P}(σ, λ) as a function of σ in different λ slices. Bottom left panel: g_{σ}(X_{off, P}, λ) as a function of λ in different X_{off, P} slices. Bottom right panel: g_{σ}(X_{off, P}, λ) as a function of X_{off, P} in different λ slices. 

In the text 
Fig. B.4. Comparison of f(X_{off, P}) and f(λ) between data and model. Top panels: solid red lines indicate the integral on the bestfit model, while the shaded blue areas represent the integral on the 3D h(σ, X_{off, P}, λ) data with 1σ uncertainties. Each bottom panel shows the residual trend with σ error; the straight black line represents the perfect match between data and model, with null residual. Top left panel: f(X_{off, P}) as a function of X_{off, P}. Bottom left panel: residual between f(X_{off, P}) data and model in logarithmic scale. Top right panel: f(λ) as a function of λ. Bottom right panel: residual between f(λ) data and model in logarithmic scale. 

In the text 
Fig. B.5. Comparison of multiplicty functions. Top panel: three shaded regions show the 1σ contours of f(σ) data directly computed on different simulations (orange for HMD, green for BIGMD, red for MDPL2); the light blue shaded region is the 1σ contour of the 2D integral computed on the concatenated sample containing all three simulations; the dashed pink line indicates the mass function from Comparat et al. (2017), while the blue solid line is the f(σ) we recover integrating our model along X_{off, P} and λ. Bottom panel: blue thick line is the fractional difference between our f(σ) and the Comparat et al. (2017) value. The light blue shaded area represent the 1σ contours of the residual between the integrated data and our bestfit model; the black horizontal line indicates the perfect match with null residual. 

In the text 
Fig. B.6. Redshift evolution of the bestfit parameters for h(σ, X_{off, P}, λ). Each panel shows a single parameter. The values at z = 0 are reported in Table B.2. The slopes of the redshift trends are given in Table B.3. 

In the text 
Fig. B.7. Marginalized posterior distributions of the h(σ, X_{off, P}, λ) bestfit parameters at redshift 0. The 0.68 and 0.95 confidence levels of the posteriors are shown as filled 2D contours. The 2.5th, 16th, 84th, and 97.5th percentiles of the onedimensional posterior distributions are indicated by the vertical lines on the diagonal plots. 

In the text 
Fig. B.8. Marginalized posterior distributions of the bestfit parameters describing the redshift evolution of h(σ, X_{off, P}, λ). The 0.68 and 0.95 confidence levels of the posteriors are shown as filled 2D contours. The 2.5th, 16th, 84th, and 97.5th percentiles of the onedimensional posterior distributions are indicated by the vertical lines on the diagonal plots. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.