Cluster Cosmology Redux: A Compact Representation for the Halo Mass Function

Massive halos hosting groups and clusters of galaxies imprint coherent, arcminute-scale features across the spectrophotometric sky, especially optical-IR clusters of galaxies, distortions in the sub-mm CMB, and extended sources of X-ray emission. Statistical modeling of such features often rely upon the evolving space-time density of dark matter halos – the halo mass function (HMF) – as a common theoretical ground for cosmological, astrophysical and fundamental physics studies. We propose a compact (eight parameter) representation of the HMF with readily interpretable parameters that stem from polynomial expansions, first in terms of log-mass, then expanding those coefficients similarly in redshift. We demonstrate good ( ∼ 5%) agreement of this form, referred to as the dual-quadratic (DQ-HMF), with Mira-Titan N-body emulator estimates for halo masses above 10 13 . 7 ℎ − 1 M ⊙ over the redshift range 0 . 1 < 𝑧 < 1 . 5, present best-fit parameters for a Planck 2018 cosmology, and present parameter variation in the 𝜎 8 − Ω m plane. Convolving with a minimal mass–observable relation (MOR) yields closed-form expressions for counts, mean mass, and mass variance of cluster samples characterized by some observable property. Performing information-matrix forecasts of potential parameter constraints from existing and future surveys under different levels of systematic uncertainties, we demonstrate the potential for percent-level constraints on model parameters by an LSST-like optical cluster survey of 300,000 clusters and a richness–mass variance of 0 . 3 2 . Even better constraints could potentially be achieved by a survey with one-tenth the sample size but with a reduced selection property variance of 0 . 1 2 . Potential benefits and extensions to the basic MOR parameterization are discussed.


INTRODUCTION
The evolving population of galaxy clusters on the sky is a cosmological diagnostic whose value has been recognized since the era when 4-m class telescopes opened the study of clusters at redshifts above 0.5 (Gunn et al. 1986;Peebles et al. 1989;Evrard 1989).The massive halos that host groups and clusters of galaxies represent a rare event tail of hierarchical structure formation that is sensitive to both the growth rate of linear structure (White et al. 1993) and the nature of the initial fluctuation power spectrum (Dalal et al. 2008).
Cosmological constraints from the aforementioned studies are sometimes inconsistent.Dark Energy Survey Year One (DES-Y1) analysis, based on counts and mean lensing masses in four richness and three redshift bins with a total sample size of 6500 clusters, find  8 = 0.65 ± 0.04, significantly (4) below the 0.830 ± 0.013 value from Planck 2018 CMB analysis (Planck Collaboration et al. 2020).In contrast KIDS-DR3 cluster population analysis (Lesci et al. 2022), based on a data vector similar to that of DES-Y1 derived from an optical sample of nearly 3700 clusters, yields  8 = 0.78 ± 0.04, 1 consistent with the Planck CMB value.
Within a given cosmology, formulating expectations for cluster counts and aggregate lensing masses of samples selected on some observable property is challenged by several sources of systematic uncertainty.The physical extent of massive halos and their preference to form in large-scale overdense regions of the cosmic web creates source confusion; the virial regions of  > 3 × 10 13 M ⊙ halos hosting groups and clusters of galaxies cover one-third of the ΛCDM sky within  < 1.5 (Voit et al. 2001).Projection tends to boost intrinsic properties (e.g., White et al. 2002;Cohn et al. 2007;Costanzi et al. 2019), but the fact that the effects of projected structure on optical, X-ray, and SZ measurements will generally differ reinforces the value of multi-wavelength cluster sample analysis.
The statistical relationship between the bulk observable properties of a halo, such as its X-ray temperature, gas or stellar mass, or galaxy richness 1 , and its true total mass is another source of uncertainty (e.g., Salvati et al. 2020;Wu et al. 2021).This relationship connects the sky+redshift-space abundance of a property-selected cluster population to the space-time density of massive halos.The differential form of the latter point density, known as halo mass function (HMF), is now well characterized in the space of standard ΛCDM cosmological parameters by large N-body simulation campaigns (e.g., Bocquet et al. 2020, and references therein).
A convolution of the HMF with the mass-observable relation (MOR) is the basis of survey statistical expectations, and a powerlaw mean with log-normal variance is a canonical MOR form motivated by cosmological hydrodynamics simulations (e.g., Bryan & Norman 1998;Angulo et al. 2012;Farahi et al. 2018;Anbajagane et al. 2020, and references therein).Over a wide dynamic range in mass, a single power law mean may be insufficient, especially for hot gas properties (Farahi et al. 2018;Pop et al. 2022), and the variance may also be mass-dependent (Anbajagane et al. 2020).Extensions to accommodate such behavior are discussed in §5.5.
The convolution naturally joins cluster astrophysics, encapsulated by the MOR, to the (primarily) cosmology-driven HMF, and the parameter couplings of these spaces have been explored previously (Evrard et al. 2014, hereafter, E14).The model we present here extends previous work by letting the HMF shape parameters vary continuously with redshift.Essentially, E14 introduced approximate HMF forms at a fixed epoch in order to develop expressions for conditional statistics of samples selected by an observable property.This present work develops a continuous space-time representation of the differential space density at high halo masses with the goal of constructing a compact, interpretable form of the HMF.
In this paper, we first show that a simple, eight parameter representation captures the near-field ( < 1.5) group/cluster HMF derived by the Mira-Titan universe ensemble (Bocquet et al. 2020).Coupled with a log-normally distributed MOR, we derive closedform expressions for both the evolving space density and the logmean selected mass of the group/cluster population as a function of the selection property and redshift.
With those ingredients, we then perform an information matrix (IM) 2 analysis to explore the ability of current and future cluster surveys -specifically, those based on counts and mean gravitational lensing mass as a function of a chosen selection property and redshift -to constrain the parameters of this HMF form.The forecasts require information on the expected uncertainty in lensing mass measurements as well as the uncertainty in the MOR variance, and we consider both current estimates and future advances in our projections.
Another reason is that HMF-centric analyses can exploit increasingly tight constraints on the differential comoving volume element, /, from baryon acoustic oscillations (BAO, Alam et al. 2021;Abbott et al. 2022) and Type Ia supernovae (SN, Guy et al. 2010;Abbott et al. 2019;Brout et al. 2022;Mitra et al. 2023).In current methods of cluster survey analysis, / is left free to vary in the space of ΛCDM parameters.
From a practical perspective, a compact representation for the HMF at galaxy cluster scales can serve as a consistency check among cluster samples selected at different wavelengths.As the common ground that underlies all cluster samples, the inferred HMF needs to be consistent across surveys and independent of the chosen sample selection property.An important feature of our model is that it naturally incorporates multiple, intrinsically correlated physical properties (Mulroy et al. 2019;Farahi et al. 2019a).
We are far from the first to emphasize the HMF shape.The original study of Bahcall & Cen (1993) used counts of nearby clusters selected by optical and X-ray properties to directly estimate the HMF of the low-redshift universe.That work benefited from the insensitivity of nearby volume to cosmological mean density parameters.Subsequent studies derived HMF estimates from X-ray samples (Reiprich & Böhringer 2002;Böhringer et al. 2017) or optical cluster samples using galaxy richness (Bahcall et al. 2003) or velocity dispersion (Pisani et al. 2003;Rines et al. 2007Rines et al. , 2008) ) as a proxy for mass.The statistical power of these samples was limited by their moderate sample sizes, typically several hundred systems.
Compact HMF representations already exist, but historically they have been expressed in terms of a similarity variable,  2 (), the rms amplitude of linear perturbations smoothed on a Lagrangian scale  ∝  1/3 (Press & Schechter 1974).An assumption about cosmology is required to convert these forms to a function of mass.The Sheth-Tormen (ST) form (Sheth et al. 2001) is a popular example, and constraints on the parameters of this model have been published from analysis of magnified images of sub-mm galaxies (Cueli et al. 2022) and from counts of GAMA groups and clusters (Driver et al. 2022).
Because the ST model represents a non-linear function of the similarity variable rather than mass directly, its parameters are difficult to interpret.By expressing the HMF directly in terms of halo mass and redshift, the eight free parameters of our model (see Table 1 below) have clear interpretations as coefficients of polynomial expansions.
The aims of this paper are twofold.We first demonstrate the model's ability to reproduce LCDM sky counts from the Mira-Titan emulator 3 in the space of fluctuation amplitude,  8 , and matter density parameter, Ω m (Bocquet et al. 2020), and provide Planck 2018 model parameters.We then apply an information matrix approach to estimate potential constraints on DQ-HMF parameters from idealized cluster surveys patterned after existing (DES-Y1, Abbott et al. 2020) and future (LSST, Chisari et al. 2019) galaxy cluster surveys.The parameter forecasts employ cluster counts and mean weak lensing masses, each derived within finite ranges of selection property and redshift, along with an additional input on the degree of scatter in log-mass at fixed value of the selection property.
In §2 we detail the model's structure, demonstrate its utility at capturing emulator predictions in the space of { 8 , Ω m }.Expressions for counts and mean mass as a function of an observable property are presented in §3, and the IM elements we employ for survey analysis are also defined there.Forecasts of parameter constraints from existing and planned cluster surveys are presented in §4.In §5 we discuss ideas for implementing the model and review how massive halos tie to many non-Gaussian LSS signatures.Benefits of selecting a sub-sample with reduced property variance are made explicit in §5.3.An appendix offers a three-parameter toy model that helps illustrate the key role of MOR variance.
We employ a cosmology with matter density Ω  = 0.311, baryon density Ω  = 0.0489, Hubble constant  0 = 67.7 km s −1 Mpc −1 , primordial spectral index   = 0.967, and power spectrum normalization  8 = 0.810, values derived by Planck2018 CMB+BAO analysis.Our measure of halo mass is  200c , the mass defined by a mean interior spherical overdensity of 200 times the critical density,   (), and we express this mass in units of ℎ −1 M ⊙ , where ℎ =  0 /100 km s −1 Mpc −1 .Our spatial density unit for the HMF is ℎ 3 Mpc −3 .The IM analysis is patterned after optical survey samples but is generalizable to samples selected by other properties.Relative to X-ray and SZ selection, optical samples have the benefit of distance estimation from spectroscopic or photometric redshifts.We ignore distance uncertainties, as the redshift bins we employ are much wider than typical uncertainties (Rykoff et al. 2014(Rykoff et al. , 2016;;Maturi et al. 2023).

METHODS
A key component of ΛCDM structure formation is an initially Gaussian random density field whose amplitude grows due to gravity.A spherical collapse model (Gunn & Gott 1972) argues for a linearly evolved perturbation amplitude threshold at which halos form.Combining these elements, the HMF was originally derived by Press & Schechter (1974) as a derivative with respect to scale of the fraction of mass in the universe that satisfies the collapse condition.At the highest masses, for which only extreme peaks in the density field can have collapsed, this fraction is an error function with large argument, and its derivative leads to a steeply falling HMF with mass.At lower masses, where more modest-sized perturbations can collapse, the HMF transitions to nearly a power-law form.
The model presented in §2.1 represents the high mass portion by the tail of a Gaussian in log-mass, meaning the log of the HMF scales as a negative quadratic with mass.These three coefficients are themselves expanded as polynomials with redshift.While E14 included a cubic log-HMF representation, we defer that approach to future work as the quadratic form captures much of the information available in cluster counts, as we show in § 2.2 below.HMF parameter values in the space of { 8 , Ω m } are presented in §2.3.

A compact form for the cluster-scale HMF
The HMF describes the comoving spatial number density of halos as a function of mass and redshift.Considering a small volume, , at some redshift, , the probability that the center of a halo of mass, , lies within that volume defines the differential HMF  ≡ (, ) ln ln  .
(1) The convention of number density per logarithmic unit of mass used above implies that the HMF has dimension of inverse volume per logarithmic unit of mass.We express the HMF amplitude in units of ℎ 3 Mpc −3 .We introduce an eight parameter model that employs low-order polynomial forms in log-mass and redshift.Letting  ≡ ln(/  ), where   is a characteristic (pivot) mass scale, we begin with the E14 quadratic form for the log of the HMF, ln (, ) The characteristic mass is essentially the pivot scale of a quadratic expansion of the log HMF, with  0 () is the normalization,  1 () the local slope, and  2 () the curvature of ln[/].We choose a pivot mass of 10 14.3 ℎ −1 M ⊙ and apply this form for  ≥ 10 13.7 ℎ −1 M ⊙ .Below this mass scale the HMF transitions to a pure power-law form (e.g., Sheth et al. 2001, and references therein).
The explicit negative sign on the RHS of equation ( 2) is used so that the   () parameters take on positive values.We choose this form over an explicit Gaussian representation because the latter would imply a global representation over a very wide mass range.Instead, we are operating on a relatively narrow mass range, roughy 1.5 decades wide, out on the Gaussian's tail, so a description using canonical location and width of a normal distribution is not as useful or meaningful.
The E14 analysis used only   values defined at a few specific redshifts.We extend that work by allowing the first pair of coefficients to run as quadratic functions of (1 + ), where   is a pivot redshift which we take to be 0.5.The mass curvature evolves linearly with redshift, Hereafter, we refer to equations (2), ( 3) and (4) as the dual-quadratic (DQ-HMF) model.Table 2 summarizes the DQ-HMF parameters for the default ΛCDM Planck2018 cosmology.The following section describes how we obtained these values using the Mira-Titan emulator.

Fitting to Mira-Titan expectations
We evaluate the model using ΛCDM expectations based on the Mira-Titan emulator (Bocquet et al. 2020) using a process guided by expectations for the LSST survey sky area of 18, 000 deg 2 (Ivezić et al. 2019).We use fourteen redshift bins ranging from  min = 0.1 to  max = 1.5, each of width Δ = 0.1.At the central redshift,   , of each bin, we evaluate the Mira-Titan HMF and convert it to a differential number density function for an LSST-like sky area, where Δ  is the volume of an 18,000 sq degree survey between redshifts   −0.05 and   +0.05.We then integrate this form to obtain counts per bin in twenty -bins between masses of 10 13.7 ℎ −1 M ⊙ and 10 15 ℎ −1 M ⊙ .We assign a Poisson uncertainty to each bin, and obtain best-fit parameters by minimizing  2 across the combined set of 280 sampled count values.This approach emphasizes fitting at lower masses, the range that provides the majority of information in cosmological surveys. 4igure 1 compares the eight-parameter DQ-HMF differential model counts to the Mira-Titan emulator values.For clarity, we show four redshifts selected from the full range used in the fit; other redshifts behave similarly.The differences between the DQ form and the emulator expectations are below 5% at masses < 5×10 14 ℎ −1 M ⊙ , increasing to tens of percent at the highest masses.
To contextualize the differences in the lower panel we first note that the fits are best at the lowest masses, where the information content is highest.In addition, the finite volumes of the Mira-Titan Nbody ensemble yield an HMF uncertainty of ∼ 10% at 10 15 ℎ −1 M ⊙ (Bocquet et al. 2020).Finally, the emulator is based on universes in which the clustered matter is purely collisionless matter (dark matter with an optional minority neutrino component).The gravitational back-reaction effects of baryons cycling through the process of compact object formation can drive HMF deviations larger than 5% over the mass and redshift range shown (Stanek et al. 2009;Cui et al. 2012Cui et al. , 2014;;Martizzi et al. 2014;Cusworth et al. 2014;Castro et al. 2021;Schaye et al. 2023).Given these uncertainties, the DQ-HMF model can be considered a sufficient representation of the cluster population in the late universe.
The model parameters resulting from the Mira-Titan emulator fits for a Planck 2018 cosmology are listed in Table 2. Points in Figure 2 show   () values determined by fitting the HMF at each sampled redshift, while lines show the redshift-continuous DQ fit with quadratic behavior of the HMF normalization and slope and linear behavior of the HMF curvature.Halos at the pivot mass scale become increasingly rare with increasing redshift -the normalization varies by roughly a factor of 100 over the redshift range shown -and the HMF shape at the pivot mass becomes both steeper and more strongly curved at earlier times.
In mapping to observable properties, the value of the slope,  1 (), is particularly important as it controls the amplitude of a convolution-induced bias (often referred to as Eddington bias) discussed below.The local slope at the pivot mass of 10 14.3 ℎ −1 M ⊙ steepens from −2 at  = 0.2 to −4 at  = 1.5, implying that the magnitude of this bias will grow by a factor of two over this redshift range.
In the IM analysis below we explore two idealized cases, one patterned after the existing DES-Y1 cluster sample, which covers roughly 5000 deg 2 over redshifts, 0.2 <  < 0.65 and another patterned after the wider, 18000 deg 2 , and deeper LSST survey.Figure 3 shows the DQ model expectations for LSST survey counts of massive halos with  200c ≥ 10 13.7 ℎ −1 M ⊙ in 0.1-wide redshift bins.A quarter million such halos should lie in the range 0.1 <  < 1.5, with the population strongly peaked near our chosen pivot redshift of 0.5.

Relating HMF shape to cosmological parameters
The DQ-HMF shape is generic in ΛCDM cosmologies.Here we use the same Mira-Titan sky expectation fitting process to map how DQ-HMF parameters vary in the canonical cluster cosmology plane of {Ω m ,  8 }.All other cosmological parameters are held constant in this exercise.Figure 4 shows the resultant behavior of the HMF shape parameters, with the top row showing normalizations,  , , at the pivot redshift, the middle the gradients with redshift,  , , and the bottom row the redshift curvature values for the HMF normalization and mass gradient,  0,2 and  1,2 , respectively.(Recall that the HMF curvature evolves only linearly with redshift, meaning  2,2 = 0.) Since the amplitude at the pivot redshift,  0, , is the primary controller of counts, it is not surprising that its contours tend to follow loci of  8 √ Ω m ≃ const. in the top-left panel.The negative of the HMF log-mass slope at   (top middle panel) is sensitive only to  8 , reducing from 2.8 to 1.9 as  8 increases from 0.7 to 0.9.The negative of the mass curvature at the pivot redshift (upper right) behaves somewhat orthogonal to  0, , with values ranging from 0.6 to 0.9 over the range shown.
The rate at which the HMF shape shifts over time depends is also dependent on cosmology.The middle row of Figure 4 shows redshift gradients of the mass expansion terms.The gradient of the normalization,  0, , is highly sensitive to  8 , scaling inversely from a low of 1.8 at  8 = 0.9 to a high of 3.5 at  8 = 0.7.The redshift evolution of the local HMF slope,  1, , exhibits sensitivity similar to that the curvature normalization,  2, , with an amplitude variation of nearly a factor of two.
The three terms of the highest order,  0,2 ,  1,2 and  2, , display mildly non-linear behaviors in the space of  8 and Ω m .Both the redshift gradient of the HMF curvature,  2, , and the second redshift derivative of the HMF slope,  1,2 , primarily depend on Ω m , but the latter shifts behavior at high  8 .Not surprisingly, the highest-order parameters are anticipated to be the least well constrained in our IM analysis below.
Note that the features displayed in Figure 4 emerge from a model employing a fixed pivot mass and redshift.Alternative choices, such as scaling the pivot mass with Ω m , would lead to slightly different outcomes.Given that Planck CMB+BAO analysis limits the matter density to within a few percent (Ω m = 0.3111 ± 0.0065) (Planck Collaboration et al. 2020), the shifts in the practical regions of these panels would be quite modest in size.

Massive neutrinos
The Mira-Titan emulator allows for a non-zero neutrino mass.Using a total neutrino mass   = 0.5 eV while keeping all other parameters constant, we find reductions in the HMF curvature and in several redshift gradients terms of order −0.05.The largest shift of −0.1 occurs in the rate of change of the HMF slope,  1, .In terms of the effect on the linear growth rate of perturbations, this is roughly equivalent to reducing Ω m by 0.02.
Recent CMB lensing analysis by the ACT collaboration (Madhavacheril et al. 2023) limits the neutrino family total mass to 0.12 eV at 95% confidence, equivalent to Ω  ℎ 2 = 0.001.For this smaller neutrino mass, DQ-HMF parameters shift at the level of 0.01 or smaller.This level of error is potentially achievable in future surveys, albeit under optimistic conditions that would take many years to develop, as discussed in §5.

OBSERVABLE FEATURES OF CLUSTER SAMPLES
Because the true 3D mass measure of the theoretical HMF is not directly observable, proxies that correlate with that mass measure are required.We use a minimal MOR based on a power-law relation with log-normal scatter, a common assumption of many survey analysis models (Rozo et al. 2010;Sehgal et al. 2011;de Haan et al. 2016;Planck Collaboration et al. 2016;Bocquet et al. 2019;Abdullah et al. 2020;Chiu et al. 2022;Lesci et al. 2022).Generically motivated by central limit theorem arguments (see, e.g., Adams & Fatuzzo 1996, for an application to star formation), this form is also measured in total gas and stellar mass statistics of halos realized by cosmological hydrodynamics simulations (Farahi et al. 2018;Truong et al. 2018;Anbajagane et al. 2020).X-ray scaling relations (Pratt et al. 2009) and lensing analysis of CAMIRA clusters (Chiu et al. 2020) support this form empirically.Generalizations of this approach are discussed in §5.
The MOR model is described in §3.1, followed by expressions for counts, mean mass and mass variance for samples selected by an observed property ( §3.2).The ingredients of our IM analysis are presented in §3.3.The first three rows of   1), and all use the same pivot mass and redshift given in Table 1.The HMF normalization,  0, (top left panel), follows the familiar negative slope traditionally derived from cluster counts.The mass gradient,  1, (top middle), is mainly sensitive to  8 while the curvature,  2, (top right), adds information somewhat orthogonal to that of  0, .The redshift gradient and curvature terms in the middle and lower rows display a range of behaviors, with the normalization terms sensitive only to  8 and the curvature's redshift derivative is sensitive primarily to Ω m .The highest-order parameter,  1,2 , shows non-monotonic behavior within a narrow range of values.A 20 × 20 sampling grid was used in the domain shown.fractional error in mass variance see Table 4 Table 3. Non-HMF parameters.The first three rows describe the massobservable relation (MOR), equations ( 6) and ( 7).The variance,  2 , is in the observed property variance at fixed mass; its inverse, the mass variance at fixed observed property,  2  , is given by equation ( 11).The two bottom rows are assumed fractional errors in mean mass and mass variance used in the IM analysis.

Mass-conditioned property likelihood
Let  be the observable property used for sample selection, made dimensionless by the choice of a convenient reference unit, and let  ≡ ln().The MOR kernel is assumed to be Normal, where  2 is the variance in  at fixed halo mass.
The mean selection property scales as a power-law in mass, meaning linearly in log-space, () =  +  . (7) While carrying value as a mass proxy,  is not a perfect indicator of mass.We consider only cases where  ≠ 0 and  2 > 0. At fixed variance, steeper proxies are better at selecting mass (see equation ( 11) below).In most practical cases of bulk observable properties, such as galaxy count or velocity dispersion or X-ray gas mass, the simple maxim that "bigger is bigger" holds, and so we generally expect that  > 0.5 For survey forecasting purposes, we assign values to the MOR parameters given in the right column of Table 3.The normalization,  = 3.1, is equivalent to an optical richness,  = 22, at the pivot mass scale of 10 14.3 ℎ −1 M ⊙ , and we assume a slope of unity.Both values are consistent with the mass-richness relation of HSC clusters (Murata et al. 2019).Other studies have found somewhat different values (see Abdullah et al. 2022, and references therein) but we do not seek to resolve those differences here.The variance of 0.3 2 is consistent with estimates derived from X-ray observations of DES-Y1 clusters (Farahi et al. 2019b).
Although, in general, , , and  2 could be functions of redshift, and the latter two also functions of mass, we consider them to be constants for the purpose of this work.In the analytic expressions below, one may simply replace these constants with appropriate functions.As this would introduce more degrees of freedom, and more sources of degeneracy, into the model, we defer such extensions to future work.Our focus here is to establish a baseline model for cluster sample statistics derived with the minimum of astrophysical complications.

Counts, mean masses, and mass variance
Motivated by DES-Y1 (Abbott et al. 2020) and similar analysis, we now consider two key observable quantities: i) the counts and; ii) mean masses of galaxy clusters disaggregated into bins of redshift and .
Convolving the HMF, equation ( 2), with the MOR kernel, equation ( 6), results in an analytic form for the space density of clusters as a function of the selection property.While originally derived in E14, that paper did not write the form explicitly in terms of  but rather implicitly in terms of the mean selected mass (see equations ( 5), ( 10) and ( 11) of that work).The explicit expression is ln where The logarithmic space density is quadratic in the observable, as expected.The last term in the second row of the expression reflects the so-called Eddington bias in mean selected mass, a topic to which we now turn.Following E14, Bayes' theorem implies that the mass distribution of clusters selected at fixed observable property, , is lognormally distributed with mean The first term in the numerator is simply the inverse of the mean MOR scaling.This value is lowered by the second term, which is approximately the HMF slope,  1 (), times the mass variance.This is the same mathematics as Eddington bias, but the source of variance differs.Eddington's scatter arose from flux measurement errors, which can be reduced by better observations.The scatter we are dealing with here is intrinsic to the population, driven by stochastic processes within a coeval population of equal-mass halos, and so cannot be reduced by improved measurement.We suggest convolution bias is a more appropriate label for this effect.Note that the log-mean mass is actually linear in the logobservable, , rather than quadratic.This behavior arises from an exact cancellation in the quadratic terms of the Bayes' theorem derivation.The result has an important implication ; additional information must be added to our IM analysis in order to invert the information matrix.We take this extra constraint to be the mass variance conditioned on the observable, , which is related to the MOR variance by Values of  2 () are of order unity, so when the MOR scatter is small then the mass scatter can be approximated by the simpler expectation,   = /||.In our IM analysis below, we impose a fractional uncertainty,  Var , on empirical estimates of this mass variance.
Rather than log-mean mass, what is directly measured via cumulative, or stacked, analysis of a cluster ensemble is a mean mass, with lensing mass derived from stacking weak lensing galaxy shear patterns (McClintock et al. 2019) or virial mass derived from an ensemble velocity likelihood (Farahi et al. 2016) being two viable methods.The log of this mean mass is shifted6 high by  2  /2 from equation ( 10), leading to the result A systematic error floor on this quantity,  ⟨  ⟩ , is employed in the following IM analysis.

Information Matrix Analysis
We use an information matrix approach to forecast DQ-HMF model parameter uncertainties anticipated from current and future cluster surveys.Such analysis, while necessarily idealized, is helpful in guiding intuition and exposing parameter degeneracies.The observable measures we consider are traditional elements (Payerne et al. 2023) of counts and mean system masses in the observable property (e.g.,richness) and redshift bins as well as estimates of the mass variance at fixed property within each redshift bin.The last two rows of Table 3 list control parameters that characterize sample data quality for mean mass and mass variance measurements.

Counts and mean masses
To derive counts,  , , within richness and redshift bins, we integrate the differential form, equation ( 8), in our chosen cosmology Here the {, } subscript denotes bins defined by the chosen limits of integration.Values for these limits are given in the relevant survey application sections of §4.An exact form for the expected mean mass in each bin requires a volume-weighted integral of the exponential of equation ( 12).Motivated by the mean value theorem, we take a simpler approach by evaluating equation ( 12) at the median property value and midpoint redshift of each bin, where the median value of  is determined by integrating the counts in each bin.

Degrees of freedom
We note that surveys with limited dynamic range in either redshift or selection property will be incapable of returning significant constraints on higher order terms of the DQ form.The cases examined in § 4 are progressively more ambitious in terms of sample size and data quality.
A survey limited to a single redshift shell centered on the pivot redshift, for example, returns no redshift gradient information, so the HMF parameters  , and  ,2 are irrelevant, leaving only the three normalizations,  , .These parameters, joined with the three MOR parameters, make a total of six.
At noted above, the forms of the observable counts, equation (13), and lensing mass measurements, equation ( 12), would return only five independent quantities: three from the quadratic counts and two from the linear log-mean mass.A data vector consisting only of counts and mean lensing mass is thus insufficient to uniquely constrain all six model parameters.To produce an soluble system, we add an empirical constraint on  2  , equation ( 11).For optically-selected DES-Y1 clusters this quantity has been derived by Farahi et al. (2019b) using X-ray temperatures of roughly 200 systems, finding   = 0.30 ± 0.04 (stat) ± 0.09 (sys).This estimate motivates our use of 0.3 2 for the default MOR variance.As improved mass estimates from lensing and dynamics become available for larger numbers of clusters, the uncertainty on this constraint is bound to improve.

Information matrix
We assume Poisson uncertainties in binned counts, a fractional error,  ⟨  ⟩ , in each mean mass measurement and a fractional error,  Var in the mass variance measurement.Using p to represent the set of model parameters, the information matrix for a single redshift bin takes the form The first term assumes Poisson variance in the count within each observable property bin and the second assumes a constant fractional uncertainty of the mean mass measured in each bin.The final term accounts for uncertainty in the mass variance, which is taken to be property-independent but depends on redshift through the  2 () term in equation ( 11).We evaluate this term at the midpoint of the redshift bin under consideration.
For the survey-specific expectations, the full information matrix is determined by a sum over all redshift bins For the DES-Y1 case, the three redshift bins used in the Abbott et al. (2020) analysis,  ∈ [0.2, 0.35), [0.35, 0.5) and [0.5, 0.65) are employed.For LSST, we use seven equally spaced redshift bins covering the interval 0.1 <  < 1.5.Appendix A provides explicit analysis of a reduced toy model based on a single redshift and property bin and having only three free parameters.This example helps illustrate the coupling of HMF and MOR parameters but the simplified scenario (two of the three MOR parameters are known) limits generalization of the results to the more complex survay applications.

HMF PARAMETER FORECASTS
We now explore potential parameter constraints from current and future optical cluster surveys under two conditions for the quality of derived mean mass per bin and mass variance.We refer to these conditions as Weak and Strong, with the latter having a factor two or better levels of uncertainty compared to the former.The choices of fractional errors in mean mass per bin,  ⟨  ⟩ , and mass variance,  Var , are summarized in Table 4.The Weak choices for DES-Y1 of (0.1, 0.6) are based on current systematic uncertainties estimates (McClintock et al. 2019;Farahi et al. 2019b), while the Strong assumption improves each by a factor of two.The LSST Weak quality are slightly improved over DES-Y1 Strong, and the LSST Strong values represent a further improvement of a factor four, to  ⟨  ⟩ = 0.01 and  Var = 0.05.
The latter constraints are certainly aspirational and will require substantial effort to achieve.For example, the fractional error in mass scatter, Δ  /  =  Var /2, is only 2.5% in the LSST Strong case, an uncertainty of only 0.0075 on a central value of 0.300.For the LSST Weak case, the uncertainty in mass scatter would be a less stringent value of 0.015.
In common with much IM analysis, the spirit of this work is to reveal best-case DQ-HMF+MOR parameter constraints from analysis of counts, mean mass and mass variance assuming no prior knowledge.Our example applications are tuned to optical cluster surveys with galaxy richness, , as the observable selection property but the model can be generalized to searches at other wavelengths.For example, on the mass scales investigated here, hydrodynamic simulations suggest that hot gas mass has a smaller intrinsic variance, ∼ 0.1 2 , at the pivot mass scale for redshifts,  < 1 (Farahi et al. 2018).Benefits of a sharper proxy are discussed in § 5.3.

Current survey application: DES-Y1
The DES-Y1 survey identified 6500 optical clusters with  ≥ 20 lying at redshifts 0.2 <  < 0.65 within roughly 5000 square degrees of the southern sky using the redMaPPer algorithm (Rykoff et al. 2016) equations ( 3) and ( 4).The model thus has eight degrees of freedom consisting of five HMF and three MOR parameters.
Applying the IM analysis using the counts and mean masses in these twelve bins, along with the uncertainty in mass variance within each redshift bin, yields the parameter constraints listed in Table 5 for the Weak or Strong quality assumptions.Figure 5 plots these parameter uncertainties, with the orange line offering a reference value of 0.1 for future reference to LSST sample expectations.
Under the Weak quality case, the normalization at the pivot mass and redshift,  0, , is forecast to have an uncertainty of 0.23, implying a fractional uncertainty of 26% in the number density.The redshift gradient of the normalization,  0, , is slightly better constrained, at 0.17, which represents a 7% fractional uncertainty on a central value of 2.38 in ΛCDM.This result is helped by the fact that the MOR is independent of redshift; the change in counts across redshift in each richness bin feeds information primarily to  0, .
The mass slope of the HMF at the pivot redshift,  1, , is forecast to have an uncertainty of 0.29, which is 14% of its ΛCDM central value of 2.38.The highest order terms, those describing the redshift gradient of the slope with mass,  1, , and the mass curvature at   ,  2, , are forecast to be weakly constrained with errors > 0.5 on central values of 1.33 and 0.75, respectively.
Under the Weak quality case, forecast errors on MOR normalization and slope are ∼ 0.1.The intrinsic variance of the observable conditioned on mass,  2 is anticipated to be returned within an error of 0.035 on a central value of 0.09.Note that these constraints come entirely from the sample itself.In practice, one might imagine external priors on these parameters being imposed in informative ways.
The couplings between parameters for the Weak quality case are displayed in Figure 6.From the analytic expressions for the space density, equation ( 8), and log-mean mass, equation ( 10), we can anticipate significant degeneracies among MOR and HMF parameters.Unsurprisingly, the two normalization parameters,  0, and , are strongly coupled, as are the two slope measures,  1, and .The MOR intrinsic variance,  2 , couples strongly to all of these parameters.
The fact that the DES-Y1 richness threshold of 20 lies close to the MOR normalization,   = 22, means that the limiting mass scale lies close to the pivot mass   .Because counts and mean masses in higher richness bins provide leverage to only one side of   , there is non-zero covariance between the HMF normalization, slope and curvature at the pivot redshift.These correlations are somewhat weaker than those associated with the MOR.Because the MOR is assumed to be a pure power-law, with zero curvature, there is very weak coupling between MOR parameters and the HMF curvature,  2, .
Forecasts for the Strong quality case are shown as filled circles in Figure 5.The MOR sector receives the primary benefit of improved quality in mean mass and mass variance, with improvements close to a factor of two.In the HMF sector, the pivot normalization and mass slope,  0, and  1, , see significant improvement while the remaining terms improve only modestly.Because the counts in each bin are the same for the two cases, the redshift gradient parameters,  0, and  1, , are little improved.
A picture that emerges is that improved measurements of mean mass and mass variance tighten the MOR sector, and these improvements filter primarily into the HMF pivot normalization and slope and secondarily to higher-order HMF parameters.This behavior is repeated in the LSST analysis below.
Note that the forecast uncertainty in HMF pivot normalization does not include any contribution from errors in distance measurements.At the forecast level of 0.12 for the Strong quality case, current volume uncertainties are sub-dominant, though still contribute at the level of ∼ 0.08 (Alam et al. 2017;Abbott et al. 2022).

Future survey application: LSST
The Rubin Observatory Legacy Survey of Space and Time (LSST) will be both wider and deeper than DES (Ivezić et al. 2019;Chisari et al. 2019).The increase in depth will yield improved measurements of galaxy shapes and colors, and this improvement should translate to more precise estimates of weak lensing mass.For the quality of mean mass estimates, we employ systematic error levels of 0.04 (Weak) and 0.01 (Strong).For the quality level of mass variance we assume values of 0.20 and 0.05, respectively.The LSST Weak values represent modest improvements over the DES Strong case.
To push this idealized case further, we anticipate that improvements in optical cluster finding will allow for a factor of two reduction in the sample richness limit, to a value of 10.The overall number of clusters expected using this observable threshold is 380,000, and the redshift distribution of the counts is similar to that shown for the mass-limited case of Figure 3. Note, however, that in the IM analysis we use seven redshift bins, each of width 0.2, between 0.1 <  ≤ 1.5, as well as five richness bins comprised of the four DES-Y1 bins joined with  ∈ [10, 20).The total number of terms  4). in the IM is 77, consisting of 35 counts, 35 mean masses, and seven mass variance measures.
The larger counts and better quality assumptions increase the volume of the IM matrix determinant relative to the DES-Y1 case (see the toy model in Appendix A).We thus employ the full set of eight DQ-HMF parameters.
Figure 7 and listed in Table 6 show that, despite the increase in model dimension, the larger information content improves the constraints on all parameters relative to DES-Y1.In the left panel of Figure 7, the horizontal dashed line reproduces the 0.1 amplitude in Figure 5.While none of the forecast uncertainties fall below this value for DES-Y1, in the LSST-Weak case all but two high-order parameters,  1,2 and  2, , lie below it.
In the LSST-Weak case, all three HMF shape parameters at the pivot redshift,  0, ,  1, and  2, are forecast to have uncertainties of 4 percent or better.For the Strong quality case, the pivot normalization and mass slope have forecast errors of one percent.Tight constraints on the three MOR parameters emergy, shown in the right panel of Figure 7.For the Weak quality case, the MOR normalization and slope are forecast to have uncertainties of 0.015 and 0.012, respectively.In the Strong case, the forecast errors below 0.004 for both.Such sub-percent errors reflect the powerful potential of LSST-era samples, but achieving such tight constraints will be difficult in practice, as discussed in §5 below.
The property variance uncertainties translate to errors on the richness scatter (square root of variance), of ±0.01 and ±0.003, in the Weak and Strong cases, respectively.These values, roughly three and one percent fractional errors on the central value  = 0.3, will again be quite challenging to achieve in practice.
The full IM covariance structure under Strong data quality assumptions for LSST is shown in Figure 8.As in the DES-Y1 example, the MOR variance and normalization parameters ( 2 and ) remain strongly coupled, as are the MOR and HMF normalizations ( and  0, ) and slopes ( and  1, ).
For the LSST case, our choice of pivot redshift,   = 0.5, lies below both the median sample redshift and the redshift midpoint of 0.8.As a result, redshift evolution parameters of the DQ-HMF are coupled.For example, the redshift gradient of the HMF normalization,  0, , is mildly correlated with its redshift curvature,  0,2 , and the redshift curvature of the local slope,  1,2 .While a more optimal choice of pivot redshift would reduce these correlations, we maintain a common pivot redshift for both samples analyzed here in order to provide a fair comparison of potential gains.
Relative to the DES-Y1 analysis, the lower richness threshold assumed for the case of LSST offers leverage below the pivot mass scale of 10 14.3 ℎ −1 M ⊙ .Unlike the DES-Y1 case, the HMF shape parameters at the pivot redshift ( 0, ,  1, ,  2, ) are largely uncorrelated for the case of LSST.
Increasing the richness limit to 20, the overall sample size drops to 60,000, a factor of roughly 2 2.5 lower than the richness 10 counts.Parameter constraints are degraded accordingly, with forecast errors in the Strong case of 0.04 in  0, and 0.08 in  1, and  1, .

DISCUSSION
In § 2.2 we showed that a compact form sufficiently captures the near-field, space-time density of high mass halos derived from an N-body emulator, where sufficiency here is in relation to current systematic uncertainties associated with the effects of galaxy formation feedback.The eight DQ-HMF parameters have straightforward interpretations as polynomial coefficients in log-mass and redshift.
A key benefit of the model is that convolution with a log-normal MOR produces closed-form expressions for observable features of group and cluster samples: counts, mean mass, and mass variance as a function of an observable property and redshift.Information matrix analysis designed around existing and planned optical cluster surveys indicate that potential constraints from an LSST-scale survey could be percent-level on many DQ-HMF and MOR parameters.
Achieving such precise constraints will be challenging.We begin by discussing the role of volume uncertainties, projection effects and more complex MOR forms.We then briefly touch on sample selection, focusing on the potential benefits of joining multiwavelength samples.In the case of LSST, we show that constraints on all model parameters could be improved using a one-tenth subsample of clusters having a tighter mass proxy, with intrinsic property variance of 0.1 2 rather than 0.3 2 .Machine learning techniques that employ all available measurements could provide a pathway to classifying such a sample, particularly if tuned accurately by synthetic data from cosmological hydrodynamics simulations.

Comoving volume uncertainties
In the IM forecasts above, we have ignored uncertainties in the comoving comoving cosmic volume.Uncertainties in cosmic volume will introduce additional error to the normalization terms,  0, ∈ { 0, ,  0, ,  0,2 }.At the chosen pivot redshift of 0.5, SDSS-III measurements of baryon acoustic oscillation (BAO) and galaxy clustering (Alam et al. 2017) constrain the local cosmic volume to within 4%, and the 2.7% distance measurement to  = 0.835 from DES BAO analysis (Abbott et al. 2022) implies a roughly 8% error in volume.These uncertainties are subdominant to the DES-Y1 IM errors in  0,  (Table 5) but achieving future LSST constraints (Table 6) will require more precise distance measurements.
Increased precision will almost certainly come.For example, current forecasts for distances derived from Type Ia SN in LSST suggest comoving distance errors of roughly 0.25% in the redshift range 0.5 <  < 1.2 (Mitra et al. 2023), meaning volume errors below one percent.This is similar to the 1% LSST-Strong constraint on  0, from the IM analysis.Alternatively, the HMF could simply be redefined in terms of the directly observable volume element, in units of number per square degree per unit redshift rather than cubic megaparsecs.The above IM forecasts apply directly to this alternative HMF framing.

Projection and more complex MOR forms
A galaxy cluster sample defines a discrete population in which each member is minimally defined by a location on the sky and (ideally) redshift, an angular size, and one or more aggregate observable properties, preferably measured within that aperture.The inverse mapping of a given cluster sample onto the underlying space-time population of massive halos is complicated by several effects arising from projection and other factors such as halo orientation.In addition, the minimal MOR used above may require extensions for practical application to specific cluster surveys.
These non-trivial issues pose a challenge to precise modeling of the cluster-halo connection, especially for optically-selected samples (e.g., Wetzell et al. 2022;Zhou et al. 2023;Varga et al. 2022;Giles et al. 2022;Upsdell et al. 2023;Zhang et al. 2023).We sketch here some ideas for how to incorporate them into DQ-HMF-focused analysis.

Projection
Typically, a single massive halo subtends a few arcminutes of sky and, due to its origin as a peak in an initially Gaussian noise field (Kaiser 1984), tend to be more strongly clustered than the general dark matter distribution.The intrinsic properties of a given halo are thus superposed with projected contributions from other halos along the same line-of-sight.A general way to accommodate this effect on measured properties is by adding another statistical factor, ( obs |, , ), that accounts for projection-induced distortions (e.g., Mulroy et al. 2019).This function will introduce additional parameters, prior values of which can be estimated by survey-specific simulations (Costanzi et al. 2019;Chiu et al. 2020; LSST Dark Energy Science Collaboration (LSST DESC) et al. 2021).
Projection will generally boost aperture-based signals (White et al. 2002;Cohn et al. 2007;Costanzi et al. 2019), driving positive skewness into the observed property kernel, ( obs |, ).The inverse kernel, (| obs , ), will lean toward lower halo masses, and this implies a similar lean in potential well depth measures such as X-ray temperature (Ge et al. 2019).
In terms of the model, kernel skew can be accommodated multiple ways, including by a Gaussian mixture where the first term represents a majority fraction,  , of clear sightlines with mean  obs () and variance  2 , and the second term represents a highly-projected subset boosted in the mean by Δ  with variance  2  .This form, which is supported by red sequence cluster finding using Millennium Simulation galaxies (Cohn et al. 2007), brings the benefit of retaining the analytical forms in §3.2 which would be fast to compute in survey analysis.A downside is the introduction of three additional parameters, but these dimensions could be coupled and reduced to a simple skewness measure implemented by Markov Chain Monte Carlo chains.

Intrinsic MOR Complexity
Intrinsic property statistics are sensitive to both cosmology (through environmentally-sensitive formation histories) and astrophysics related to galaxy formation and plasma evolution.The minimal MOR form used above, with three parameters, is likely to require some extensions for precise survey likelihood application.Based largely on the behavior of halos in cosmological hydrodynamics simulations, we briefly outline modifications that may apply to different observable properties.
MOR shapes from cosmological hydrodynamics simulations.Large samples of high-mass halos from cosmological hydrodynamics simulations provide the means to test the MOR kernel for multiple observable properties.In BAHAMAS+MACSIS simulations, the hot gas mass and the total stellar mass within  200c follow log-normal kernel shapes (Farahi et al. 2018, hereafter F18).
The existence of a log-normal PDF for the total stellar mass of halos was confirmed using three independent cosmological hydrodynamics simulations by Anbajagane et al. (2020).That work also finds slight skewness in halo mass-conditioned statistics for the total number of satellite galaxies,  sat , and the BCG stellar mass,  ★,BCG .A common Gaussian mixture fit is derived for the normalized  sat kernel, with 79 ± 1 percent of halos in a dominant component with mean, 0.28 ± 0.01, and scatter, 0.68 ± 0.01, and the remaining 21% component having mean −1.04 ± 0.05 and scatter 1.13 ± 0.02.More work is needed to understand intrinsic MOR shapes for other observable properties, such as X-ray luminosity and temperature or thermal SZ decrement amplitude, and efforts to verify statistic forms from different cosmological hydrodynamics methods are also warranted.
Running of MOR parameters with redshift and/or halo mass.The property normalization, , is likely to evolve with redshift.A self-similarity assumption (Kaiser 1986) that ties physical properties to the evolving critical density is often used to express, (), in terms of powers of  () ≡  ()/ 0 .7Under strict self-similarity, the total stellar or gas mass fractions are independent of redshift.In the BAHAMAS+MACSIS simulations, F18 find modest (several percent) redshift dependence in both measures, with the gas mass fraction declining, and stellar mass fraction increasing, slightly from  = 1 to  = 0.These shifts are mildly mass-dependent, being larger at lower halo masses that are more strongly influenced by galaxy evolution.Free parameters introduced to capture deviations in normalization from self-similarity would couple most strongly to the DQ-HMF normalization parameters,  0,  .The intrinsic property variance,  2 , of hot gas and stellar mass was also found to run weakly with mass and redshift by F18.
The constancy of the MOR slope, , is also a simplification that may require modification for some properties.For example, F18 find that the slopes of hot gas mass and stellar mass vary modestly with both halo mass scale and, for the former, redshift.At lower halo masses, the hot gas mass slope steepens to values above unity, and the stellar mass scaling becomes shallower than unity.A parameter introduced to describe an MOR slope gradient, /, would couple most strongly to the MOR curvature,  2, .Extending further to allow for this parameter to run linearly with (1 + ) would then couple to  2, .

"Gold Sample" Selection with Machine Learning using Multiple Properties
Cluster samples are generally defined by a threshold in a single observed selection property.The DES-Y1 sample, for example, is limited by red galaxy richness,  ≥ 20.The mapping between a set of observed clusters and their underlying host halos is assumed to be bĳective; a chosen halo maps uniquely to a single cluster, and viceversa.This is not always the case8 , and multi-wavelength studies are critical to understanding how frequently this assumption is violated.
In a recent joint study of cluster samples identified independently by X-ray and optical observations in roughly 60 deg 2 of sky, Upsdell et al. (2023) find that only one of 178 X-ray sources has two optical clusters identified along the same line of sight.Such effects, as well as more prosaic issues such as survey masking (e.g., Rykoff et al. 2016), will affect cluster selection and require calibration by multi-wavelength observations and simulations.Joint property analysis of large cluster samples can improve cosmological parameter constraints (Cunha 2009) because combining multiple observable properties can substantially reduce mass variance relative to single-property characterization Ho et al. (2023).The anti-correlation of hot gas and stellar mass contents observed in the LoCuSS sample (Farahi et al. 2019a) is an important feature; selecting on just these two intrinsic properties in the Magneticum simulation yields a variance in halo mass of 0.05 2 (Ho et al. 2023).

Potential Gains of a Gold Sample
There is potential to improve DQ-HMF parameter constraints using a selection approach that identifies a Gold Sample of clusters with reduced intrinsic MOR variance.For this example, we imagine a classifier returning 10% of the overall population with intrinsic MOR variance, 0.1 2 .While a significant improvement over the 0.3 2 value used in our default analysis, we note that, for high halos masses, the hot gas mass is seen to have such a small variance (Truong et al. 2018;Farahi et al. 2018;Pop et al. 2022;Farahi et al. 2022;Pellissier et al. 2023).
Using the reduced, three-parameter model of Appendix A as a guide, the information volume scaling of  −4 for low-scatter proxies (other parameters held fixed), equation (A15), would imply that the improvement in MOR variance wins over the decrease in sample size.Figure 9 confirms this to be the case.The filtered cluster subsample with 10 percent of the counts but 0.1 2 variance yields improvements in all HMF parameters, with the biggest gains occurring for the highest order quantities,  0,2 ,  1,2 and  2, .As discussed in §2.3.1, the shifts in such higher-order terms caused by massive neutrinos are of the order 0.01, potentially within reach of Gold Sample analysis.
Machine learning (ML) techniques have been demonstrated to yield improved estimates of galaxy cluster masses from noisy observations derived from simulations of massive halos (Ntampaka et al. 2016(Ntampaka et al. , 2019;;Cohn & Battaglia 2020;Krippendorf et al. 2023;Ho et al. 2023), and sample selection in the low signal-to-noise regime has been explored by Kosiba et al. (2020).Symbolic regression has been used to identify property combinations that minimize mass variance (Wadekar et al. 2023) and random forest techniques have been used to classify galaxies into orbit classes using projected phase space information (Aung et al. 2023;Farid et al. 2022).
We encourage other researchers to explore whether ML methods can be trained to identify a Gold Sample with characteristics similar to that assumed above.Synthetic sky maps and catalogs  (Han et al. 2021).As multiple synthetic skies that jointly meet the requirements of surveys in optical/IR, sub-millimeter and X-ray become available, methods for sample selection can be cross-verified, trained on one simulation methodology and tested on another.

Lensing and Correlated LSS Measures
Massive halos impose peaks in weak lensing maps on arcminute scales, and tangential shear analysis has long been a staple method of estimating the underlying true halo masses of galaxy clusters (Tyson et al. 1990;Miralda-Escude 1991;Kaiser & Squires 1993;Luppino & Kaiser 1997), see the review of Hoekstra et al. (2013).
Weak lensing peaks contain information on cosmological parameters including neutrino mass (Ajani et al. 2020;Zürcher et al. 2022;Liu et al. 2023).In addition, the spatial auto-and cross-correlations of galaxies, gravitational lensing and both thermal and kinetic SZ maps contain some degree of information about massive halos, and higher-order statistical signatures at non-linear scales are even more strongly connected. 9The spatial clustering of the cluster population itself is a signal that improves cosmological inference (Majumdar & Mohr 2004;Euclid Collaboration et al. 2022), and the power spectrum and bispectrum of massive halos contains potentially powerful information on primordial non-Gaussianities (Coulton et al. 2023).
Cluster counts offer complementary information to other cosmological probes, especially as the population is sensitive to both 9 For example, this Snowmass2021 Letter of Interest.cosmic geometry and the gravitational growth of structure (Frieman et al. 2008;Cunha et al. 2009).A recent study that combines DES redMaPPeR cluster counts with spatial correlations of galaxy and lensing demonstrates the value of this approach (To et al. 2021).Clusters could be used to independently assess a recent CMB+LSS finding of a 4.2 larger than ΛCDM growth factor index (Nguyen et al. 2023).
These types of studies could potentially benefit from a compact mass function form, as DQ-HMF parameters could be used either as informative priors or as part of the focus of posterior likelihood evaluation.

Other Caveats and Extensions
We mention here a few additional caveats and potential extensions.
Alternative Mass Conventions.In N-body simulations, the mass of a halo is typically defined by percolation or spherical overdensity approaches (White 2001).spherical (see, e.g., Diemer 2020, and references therein).For the spherical overdensity approach, common choices for the interior mean density threshold and/or the reference density (critical or mean mass are typical choices) induce scale-dependent shifts in mass.The resultant HMF forms are follow similar forms, however, remain similar and can be converted using mean mass density profile shapes (see Appendix B of Evrard et al. 2002).We suspect, but do not attempt to prove here, that a compact representation would be valid for most, if not all, existing conventions for true halo mass.
Alternative Formulations for Extended Dynamic Range.Our model aims at near-field studies of groups and clusters.To extend to model to lower-mass halos, one could include a transition mass scale below which the HMF would become pure power-law.We note that the pure power law form at low masses ignores effects of baryon feedback during galaxy formation.A recent internal structure study of halos across nearly six orders of magnitude in mass in the IllustrisTNG simulations (Anbajagane et al. 2022) finds wiggles in dark matter halo scaling relations near the Milky Way mass of 10 12 ℎ −1 M ⊙ , where star formation efficiency in the late universe peaks (Behroozi et al. 2013).This finding suggests that the HMF may also have a localized deviation from a pure power-law form at that scale.
The near-field halos above our chosen limiting mass of 10 13.7 ℎ −1 M ⊙ comprise several percent of the overall matter density at  < 1.5, but this fraction becomes negligible at much higher redshifts.The mass scale associated with the most extreme few percent of the halo population declines with redshift, reaching Milky Way-scale halos that host bright galaxies at  > 8, as seen in JWST observations (Boylan-Kolchin 2023).
To span a wider range in redshift, one could redesign the model by reframing the normalization.Instead of the number density at fixed mass,  0 (), one could employ a mass scale at fixed number density parameter, for example,  −6 () to represent the mass scale at which the comoving space density is 10 −6 ℎ 3 Mpc −3 .To avoid cosmic volume uncertainties, the space density itself could be reframed in observable terms, in units of number per square degree per unit redshift.
Beyond binning.The IM forecasts employ binned values for key sample characteristics of counts and mean mass.As multiple observable properties become available for larger population ensembles, a likelihood analysis that considers each system's true mass as additional model parameters (Mulroy et al. 2019) could prove powerful.
Multi-property statistics.The expressions derived in E14 for selection property-conditioned statistics still apply.We emphasize above only the mean mass and mass variance conditioned on the selection property,   , but expressions for one or more additional properties,   (see equations ( 12) through ( 14) of E14) remain applicable, except now the HMF mass-shape parameters are explicitly redshift dependent,   →   ().

SUMMARY
We introduce a compact representation for the differential space density of high mass halos that host groups and clusters and demonstrate its utility to match well the output of the Mira-Titan emulator of purely collisionless universes for masses > 10 13.7 ℎ −1 M ⊙ in the near cosmic field of redshifts  < 1.5.Convolving with a minimal MOR yields analytic forms for the space density and propertyselected statistics that explicitly expose parameter degeneracies and that are fast to compute.Such a compact representation offers a common ground for cluster sample analysis independent of selection method.With roughly one million halos above 10 14 M ⊙ available on the full sky (Allen et al. 2011), and studies of protoclusters at moderate redshifts in ascendancy (Alberts & Noble 2022), there is abundant information available from galaxy cluster surveys.Unlocking that information will require careful modeling of sample selection, an endeavor that will be aided by sophisticated sky maps (e.g., Schaye et al. 2023).Near-term, more efforts to empirically study the MOR using high quality multi-wavelength data are needed.As the sample size of clusters with multiple well-measured properties grows from tens (e.g., Mulroy et al. 2019) to hundreds (e.g., Giles et al. 2022;Upsdell et al. 2023) to thousands, the detailed form of the multi-property MOR will come into focus, which can unlock more precise estimates of the underlying true mass of each system and, via collective study, the HMF and its behavior over cosmic time.

Figure 1 .
Figure1.The upper panel shows DQ-HMF fits (solid) to the Mira-Titan emulator expectations (dashed) for counts of halos above 10 13.7 ℎ −1 M ⊙ centered at redshifts shown in the legend over the Rubin-LSST area of 18, 000 deg 2 .A total of 280 sampled counts -20 mass bins in each of fourteen redshift shells covering the interval 0.1 <  < 1.5 -are used to fit the eight parameters of the model (see Table2); we show only a subset for clarity.Poisson uncertainties applied to each binned count yield a model that best fits lower masses.The lower panel displays the fractional deviation of the fits,  DQ / MiraTitan − 1, with the grey band highlighting 5% agreement.Baryon effects associated with galaxy formation drive deviations at this level or larger (e.g.,Castro et al. 2021), and the emulator itself is uncertain at the 10%-level at 10 15 ℎ −1 M ⊙(Bocquet et al. 2020).

Figure 2 .Figure 3 .
Figure2.DQ-HMF model parameters,   (), derived from fitting Mira-Titan sky count expectations in 0.1-wide redshift shells centered at the redshifts given by the points.A Planck 2018 ΛCDM cosmology is assumed.Lines show the fits to the redshift-continuous forms, equations (3) and (4), with parameter values given in Table2.

Figure 4 .
Figure 4. Contours showing DQ-HMF parameter values in the  8 and Ω m plane.Values are made positive by definition, equation (1), and all use the same pivot mass and redshift given in Table1.The HMF normalization,  0, (top left panel), follows the familiar negative slope traditionally derived from cluster counts.The mass gradient,  1, (top middle), is mainly sensitive to  8 while the curvature,  2, (top right), adds information somewhat orthogonal to that of  0, .The redshift gradient and curvature terms in the middle and lower rows display a range of behaviors, with the normalization terms sensitive only to  8 and the curvature's redshift derivative is sensitive primarily to Ω m .The highest-order parameter,  1,2 , shows non-monotonic behavior within a narrow range of values.A 20 × 20 sampling grid was used in the domain shown.

Figure 6 .
Figure 6.Parameter covariance forecast for a survey patterned after DES-Y1 under Weak quality assumptions (see Table4).

Figure 7 .
Figure 7. IM-forecasted uncertainties on DQ-HMF (left) and MOR (right) parameters from an LSST-like optical survey of  > 10 clusters under the quality assumptions listed in Table 4. Five richness bins and seven redshift bins across 0.1 <  < 1.5 are employed.The reference value of 0.1 is repeated from Figure 5. Parameter correlations are displayed in Figure 8.

Figure 8 .
Figure 8. Parameter covariance for an LSST-like survey under Strong quality constraints for weak lensing mass and mass variance measurements.

Table 1 .
Summary of DQ-HMF model parameters.The   () terms represent negatives of the normalization, slope and curvature of the log-space HMF at redshift  using  = ln( /  ), equation (2), and units of ℎ 3 Mpc −3 .Rows two through four list the eight core HMF parameters; elements of the polynomial redshift expansions of the   () terms around   , equations (3) and (4).Values for these parameters are determined by fitting to ΛCDM Mira-Titan emulator predictions listed in Table2.The final three rows list our choices of pivot locations in mass and redshift as well as the minimum mass scale of the Mira-Titan fits.

Table 2 .
DQ-HMF parameters of the Mira-Titan ΛCDM model.The  0 normalization at   = 0.5 is equivalent to a space density of 10 −5.85 Mpc −3 for Hubble parameter ℎ = 0.677.

Table 2
Table 3 list the parameters needed to describe the MOR: a slope, normalization and variance.The last three rows introduce data quality measures used in the IM analysis.

Table 4 .
Assumed fractional errors mean mass and mass variance for the IM analysis, equation 15.For each sample, two levels of quality, Weak and Strong, are used, with the latter improving over the former by a factor of two (DES-Y1) or four (LSST).

Table 5 .
DQ-HMF and MOR parameter constraints anticipated from a sample patterned after  > 20 DES-Y1 clusters, shown in Figure5.
Our reference model, which uses a different cosmology and MOR, yields a similar total count but with slightly different counts per redshift bin(1498, 2286, and 2752), a level of agreement acceptable for the purpose of illustration.We employ the same four richness bins as DES-Y1 in the IM analysis.Due to the limited redshift range in the DES-Y1 sample, we ignore the highest-order terms from the   redshift expansions,

Table 6 .
DQ-HMF and MOR parameter constraints anticipated from a sample patterned after  > 10 LSST clusters, shown in Figure7.
HMF parameter constraints for the LSST-Strong case with MOR variance,  2 = 0.3 2 (filled circles, same as Figure7), are compared to those from a cleaner subset ("Gold Sample") consisting of 10% of the former sample with a reduced MOR variance of 0.1 2 (open circles).The clean subset yields improvements, particularly in the higher-order quantities such as  0,2 ,  1,2 , and  2, .Note the logarithmic scale on the constraint amplitude.