The mass function dependence on the dynamical state of dark matter haloes

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. We characterize how the dynamical state of haloes impacts the halo mass function at the high-mass end. We used the dark matter-only MultiDark suite of simulations and the high-mass objects M>2.7e13 M/h therein. We measured mean relations of concentration, offset, and spin as a function of halo mass and redshift. We investigated the distributions around the mean relations. We measured the 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. We confirm the 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 concentration, offset, and spin with modified Schechter functions. The concentration of low-mass haloes shows a faster redshift evolution compared to high-mass haloes, especially in the high-concentration regime. The offset parameter is 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, 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. This 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.


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 (PS 1974) 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 E-mail: rseppi@mpe.mpg.de (Bond et al. 1991). This provides a good prediction for highmass haloes, but it predicts too many low-mass objects. The introduction of ellipsoidal collapse corrected these differences between simulations and theory (Sheth & Tormen 1999. 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;Comparat et al. 2017;Bocquet et al. 2020). A robust way to build such a mass function model is using N-body simulations (Kravtsov et al. 1997;Springel 2005), for example Mul-tiDark (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 X-ray, SZ, or optical cluster surveys such as eROSITA (Mer-Article number, page 1 of 26 arXiv:2008.03179v2 [astro-ph.CO] 10 Jun 2021 A&A proofs: manuscript no. 39123corr_rseppi loni et al. 2012;Predehl et al. 2020), Planck (Zubeldia & Challinor 2019), SPT-3G (Benson et al. 2014), CMB S4 (Abazajian et al. 2019), SPIDERS/eBOSS (Dawson et al. 2016;Finoguenov et al. 2019b), DESI (DESI Collaboration et al. 2016), 4MOST (de Jong 2011;Finoguenov et al. 2019a), Euclid (Laureijs et al. 2011), LSST (LSST Science Collaboration et al. 2009), and WFIRST (Spergel et al. 2015). These future surveys will provide tighter constraints on mass-observable 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 λ = J √ E/GM 5/2 (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 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 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, 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 X-rays (Eckert et al. 2011;Käfer et al. 2019Käfer et al. , 2020. It possibly offers a solution to mitigate biases in a cosmological interpretation of clusters abundance. This might improve cosmological constraints using X-ray selected clusters (Ider Chitham et al. 2020).
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 N-body 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 func- Table 1: Correspondence between mass, peak height, and variance of the linear density field at z=0 and z=0.5.

Mass
z=0 z=0.5 log 10 M /h ν = δ c /σ log 10 (1/σ) ν = δ c /σ log 10 (1/σ) 10 0. tion as a function of X off and λ in Sect. 6. We present the results of the fit and the best-fit parameters in Sect. 7. We summarize our findings in Sect. 8, and discuss their implications for cosmological studies.

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 top-hat 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. 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 σ (Equation 1), peak height (Equation 2), and mass interchangeably.
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.

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 K-Space Topologically Adaptive Refinement (rockstar) and consistentTrees algorithms (Behroozi et al. 2013). They are based on adaptive hierarchical refinement of friends-of-friends (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 low-mass haloes with fewer than 30-40 particles, where the identification of substructures is not straightforward (Knebe et al. 2011). The low-mass 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 (Equation 5). The virial mass function has been shown to be the one that comes closest to universality (Despali et al. 2016).

Concentration, offset, spin: Empirical relations with peak height and redshift
In the footsteps of Klypin et al. (2016) and Rodriguez-Puebla et al. (2016), we analyze the average relations linking concentration, λ, and X off to the peak height, respectively in sections 4.1, 4.2, and 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 equations 10, 12, and 14. The PDFs at different redshifts are modeled independently.

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 Equation  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 anti-correlates 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;Poveda-Ruiz 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 high-mass upturn (e.g., Klypin et al. 2016;Diemer & Joyce 2019); semi-analytic models based on the Press-Schechter 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 high-mass haloes at relatively low redshift (z < 1.5), which is particularly interesting for haloes hosting galaxy clusters.

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 best-fit values are in Table A.1. We find best-fit 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 (ν ∼ 0.95 at z=0). The low-mass end is not sampled due to the particle mass resolution. We do not consider a high-redshift regime where there is a statistical limitation of high-mass 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. (2019), due to the smaller volumes analyzed (400 Mpc/h and 500 Mpc/h, respectively). They are not large enough to obtain a significant number of high-mass objects and probe the upturn.
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 Equation 10): The best-fit 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 mass-dependent terms (Equation 10). We use the best-fit 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.

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 high-mass 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 high-mass 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 high-mass haloes gives a smaller con-tribution to the total PDF in the high-concentration tail. At higher redshift we see high-concentration haloes with both high and low mass. The opposite holds for low-mass haloes, which contribute more to the total PDF at higher z in the low-concentration 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 low-mass distribution is much broader than the high-mass one. There is also a general redshift trend: with decreasing redshift the number of highconcentration, low-mass haloes increases. More recently in time the distribution is flatter; in fact, the slope of the power-law 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 low-mass haloes evolves more than that of highmass objects. For example, at concentration c = 8 the PDF of high-mass haloes (1 × 10 14 < M < 2 × 10 14 M /h) changes by 0.32 dex between z = 1.43 and z = 0, while for low-mass objects (2 × 10 13 < M < 3 × 10 13 M /h) 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 high-mass 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. 2013Ludlow et al. , 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, high-mass haloes cluster more than low-mass ones. Low-mass objects, some of which are found in isolated environments, will experience different histories. Isolated low-mass haloes will be highly concentrated, while the concentration of low-mass haloes in dense environments will stay low. Contreras et al. (2019) found that low-concentration haloes cluster more than haloes with high concentration at high redshift. This is in agreement with our results.

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 X-rays (Bianconi et al. 2013;Eckert et al. 2019). Highresolution data are needed to explore this topic, in the optical (Song et al. 2018) and X-ray (Hitomi Collaboration et al. 2018) bands. Moreover, cluster cores can be analyzed by distortion of Each set is divided into mass slices, and is color-coded accordingly. The shaded areas represent the data with 1σ error, while the straight lines indicate the best-fit 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 y-axis. 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. 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(Baldi et al. , 2019. Future SZ and X-rays 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.

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 best-fit 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 Rodriguez-Puebla et al. (2016) did.The PDF of The best-fit parameters are given in Table A.2. The upper x-axis converts peaks values into mass at z=0. the spin parameter is best fit by a modified Schechter law with no mass dependence: Distributions for different redshift snapshots are shown in Fig. 5 (see Table A.2 for best-fit parameters). The result is consistent with previous findings, (e.g., Rodriguez-Puebla 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.

Discussion
The small correlation with mass and the well-modeled 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.  12). Individual points represent spin bins used to compute the distribution; straight lines refer to the best-fit modified Schechter model. They are color-coded by redshift. We do not consider mass dependence in this relation. Each redshift slice is fitted independently. The best-fit parameters are given in Table A.2.

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 X-ray observations. Provided there is a link between X off and the brightest cluster galaxy (BCG) to Xray displacement, an estimation of the bi-variate 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 low-mass haloes have smaller offset than high-mass 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.

Model
We model the X off − σ relation with a redshift dependent power law, see Equation 13, where E(z) is the dimensionless Hubble parameter. The best-fit 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 best-fit models are shown in Fig. 6. The best-fit 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. ). Circular dots, triangles, and squares represent HMD, BigMD, and MDPL2 respectively. They are color-coded by redshift. Straight lines show the best-fit model. The best-fit parameters are given in Table A.3.
The upper x-axis converts peaks values into mass at z=0.
We describe the distribution of X off around its mean value with a modified Schechter function (Equation 14). The PDF of X off does not show mass dependency. It is included in the normalization to the virial radius R vir αM 1/3 vir . 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.

Discussion
We find a non-zero slope for the relation between X off and mass. On average, high-mass 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. High-mass haloes form in the knots of the large-scale 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-non-cool core classification of galaxy clusters in X-rays (Eckert et al. 2011), massive structures have a higher fraction of non-cool cores. This influences the high-mass tail of the mass function (see Fig. 1). . The best-fit parameters are given in Table A.3. Therefore, given an X-ray 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 power-law 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 (Equation 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.

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.

Definitions
We compute the mass function as follows: Here ∆N M is the number of haloes in each mass bin ∆ ln M 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: (16)

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 Equation 17: Here ∆N M,X off ,λ 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 X off , and s λ are the natural (base 10) logarithm of mass (X off , λ) binning. Equivalently to equation 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 X off ,λ (σ). 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.

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 equations 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 high-mass end of the mass function, using haloes with M > 2.7 × 10 13 M /h.
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. 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). 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 Equations 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 equation 22 means an uncertainty of around 15%; in this way our measurement is not dominated by Poisson uncertainty.

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.

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, h(σ, X off , λ, z, A, a, q, µ, α, β, γ, δ, e) = ...
with µ = 10 1.83 log 10 µ , 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 Equations 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 mass-dependent 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.83 log 10 µ ), 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.

Evolution with redshift
Since structures accrete matter with time and grow, the mass function depends on redshift . Part of this redshift evolution is in the mass-σ relation through the matter power spectrum (Equation 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 (Ondaro-Mallea et al. 2021). We find that the distribution functions of X off and λ show a redshift trend as well (see Sect. 4,and Equations 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 Equation 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 worst-case 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 at z=0.702 and 10 12 M /h 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 Equation 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 Equation 24:

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 high-mass end. We consider a Gaussian likelihood where D is log 10 h(σ, X off , λ) computed from Equation 17, M is the base 10 logarithm of the model (Equation 23), and E is the uncertainty of log 10 h(σ, X off , λ) (see log error in Equation 22).

Redshift zero
The best-fit 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 2014(Buchner , 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 power-law 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 (Equation 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 X off (σ, λ) 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 low-redshift samples, each version of MultiDark has its own resolution limit: 25 kpc/h for HMD, 10 kpc/h for BigMD, 5 kpc/h 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 N-body 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.
We further integrate the model in Equation 23, obtaining the distributions of X off f σ,λ (X off ) and λ f σ,X off (λ) (Equation 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).
We obtain the multiplicity function f X off ,λ (σ) 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.
2 https://docs.scipy.org/doc/ (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 Multi-Dark 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 at z=0 (10 12 M /h at z=1.4).

Evolution with redshift
To study the redshift evolution we start from the fiducial model at z=0. We add the redshift dependence (Equation 24) to the bestfit parameters in Equation 23 at z=0. We concatenate samples for ten redshift values, as described in 6.2.
Redshift dependence is shown in Fig. 11. The shaded areas include the uncertainty on the best-fit parameter at z=0 and on the z evolution, according to equation 26 where P is each parameter in equation 23, P 0 is its value at z=0, and k indicates each parameter describing the evolution in equation 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 (Figures 6 and 7).

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 Fig. 8: Single integration of the 3D model. In each panel straight lines indicate the best-fit model, while shaded areas represent the data with 1σ uncertainties. Top left: g λ (σ, X off ) as a function of X off in different mass slices. Top right: g λ (σ, X off ) as a function of σ in different X off slices. Middle left: g X off (σ, λ) as a function of λ in different mass slices. Middle right: g X off (σ, λ) as a function of σ in different λ slices. Bottom left: g σ (X off , λ) as a function of λ in different X off slices. Bottom right: g σ (X off , λ) as a function of X off in different λ slices. The integrals are defined in Equations 19. Fig. 9: Comparison between data and model of f σ,λ (X off ) and f σ,X off (λ). In the top panels the straight red lines indicate the integral on the best-fit 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: f (X off ) as a function of X off . Bottom left: Residual between f σ,λ (X off ) data and model in logarithmic scale. Top right: f (λ) as a function of λ. Bottom right: Residual between f σ,X off (λ) data, and model in logarithmic scale.
concentration-mass 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 (Equation 9). The probability density function of concentration is a modified Schechter law, with mass dependency (Equation 10). We find that the concentration of low-mass haloes has a faster redshift evolution than high-mass objects, especially in the high-concentration regime. For concentration c = 8, the PDF for high-mass haloes shifts by 0.32 dex from z = 1.43 to z = 0, while the low-mass 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 (Equation 12), in agreement with Rodriguez-Puebla et al. (2016). The offset parameter evolves with mass and redshift according to Equation 13. The negative slope of the relation suggests that low-mass haloes are typically more relaxed compared to high-mass objects. This is true at every redshift. The offset distribution around the mean value is well described by a modified Schechter function (Equation 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 (Equation 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 Equation 24. The result is shown in Fig. 11. The link with observations will be explored further in future work.  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 best-fit model; the black horizontal line indicates the perfect match with null residual.
In order to perform cosmology calculation, we used the python toolkit Colossus (Diemer 2018) (https://bdiemer. bitbucket.io/colossus/).  Table A.4. The redshift evolution is described by Equation 24; the slopes are given in Table A.5.    The full distribution of the posteriors is shown in Fig. A.1. In the posteriors, when the 1σ uncertainty on the one-dimensional distribution is smaller than 3 order of magnitudes with respect to the parameter, we write null error. The redshift evolution of each parameter is shown in Fig. 11 Table A.5.   best-fit model, while shaded areas represent the data with 1σ uncertainties. Top left: g λ (σ, X off,P ) as a function of X off,P in different mass slices. Top right: g λ (σ, X off,P ) as a function of σ in different X off,P slices. Middle left: g X off,P (σ, λ) as a function of λ in different mass slices. Middle right: g X off,P (σ, λ) as a function of σ in different λ slices. Bottom left: g σ (X off,P , λ) as a function of λ in different X off,P slices. Bottom right: g σ (X off,P , λ) as a function of X off,P in different λ slices.    Table B.2. The slopes of the redshift trends are given in Table B.3.