Halo Counts-in-cells for Cosmological Models with Different Dark Energy

The nature of dark energy may not be as simple as a vacuum energy represented by constant $\Lambda$ and $w$. Scalar-field models introduce quintessence (see Wetterich, 1988; Ratra & Peebles, 1988; Frieman et al., 1995; Copeland et al., 2006, for reviews) and allow a time-varying equation of state parametrized by a Taylor expansion $w(a)=w_{0}+(1-a)w_{a}$, where $a$ is the scale factor (Planck Collaboration et al., 2018). A cosmological constant $\Lambda$ corresponds to $w_{0}=-1$ and $w_{a}=0$. A constant $\Lambda$ and $w=-1$ is consistent with CMB+BAO+SNe+$H_{0}$ measurements within approximately 10% uncertainty with 95% confidence (Planck Collaboration et al., 2018). However, $w_{0}$ and $w_{a}$ for a dynamical dark energy equation of state are less constrained by combinations of different probes (Fig. 4 in Planck Collaboration et al., 2016). The case when $w<-1$ is ”phantom” dark energy, in which the sum of the pressure and energy density is negative, causing a ”big rip” universe (Caldwell et al., 2003). For a weakly coupled scalar field, $w$ can be parametrised by the slope parameter, $\epsilon_{s}$, whose negative values correspond to phantom models. The value of $\epsilon_{s}$ is not well constrained even combining multiple probes.

Statistics of large scale structure in the CMB, in galaxy surveys and large cosmological simulations can test theoretical models of cosmic evolution. The description of the CMB largely relies on the measurement of angular power spectra (e.g. Planck Collaboration et al., 2018). Power spectra are also widely used to validate results of dark-matter-only cosmological N-body simulations (e.g. Springel et al., 2005; Alimi et al., 2010; Klypin et al., 2011; Skillman et al., 2014). To compare simulations with galaxy surveys, halo finders (Davis et al., 1985; Behroozi et al., 2013; Knebe et al., 2013) generate halo catalogs and the halos are then populated with galaxies using theoretical models, such as halo abundance matching (Kravtsov et al., 2004; Tasitsiomi et al., 2004; Vale & Ostriker, 2004; Conroy et al., 2006). Another method to quantify large scale structure is the two-point correlation function; this is the Fourier transform of the power spectrum (Peebles, 1980; Saslaw, 2000). The power spectrum is often more convenient for continuous matter density, and the two-point correlation functions for discrete distributions of objects. Higher order spectra or correlation functions are usually more difficult to estimate accurately and are computationally challenging (Smith, & Zaldarriaga, 2011; Munshi et al., 2011).

Nonlinear clustering produces higher order correlations (e.g. Peebles, 1980; Bernardeau et al., 2002) for galaxies, halos and the underlying density field. The hierarchy of correlation functions, in principle, provides complete information about the statistical distribution of structures, but it is impractical to measure accurate correlation functions of high orders. Moreover, two-point correlation functions and power spectra also do not contain any phase information. Two spatial density distributions can have the same power spectrum, but have very different clustering morphology caused by phase differences (Chiang, 2001). Primordial non-Gaussianity in the initial density field from inflation can leave signatures of spatial phases. Nonlinear evolution can produce non-random phases (Coles, 2003). Useful statistical measures of the spatial phase distribution are yet to be discovered.

The number counts of discrete objects in an ensemble of cells of a fixed size and shape gives a CiC probability distribution function (PDF) (Saslaw, 2000). We compare the CiC PDFs in N-body dark matter cosmological simulations for different dark energy models. This method is an extension of galaxy cluster abundance, which can be considered the first moment, i.e. mean, of the CiC distribution with a high mass cut at $\sim 10^{14}M_{\odot}$. If halos greater than galaxy masses are included, the sample is much larger and nonlinear clustering is also probed. CiC distributions contain much more information from higher moments and higher order correlations (Saslaw, 2000). Galaxy-mass halo distributions can be related to the galaxy distributions in galaxy surveys through a mass-luminosity relation. Halos and galaxies are both biased tracers of the underlying dark matter density, but simulations have not found any qualitative difference between the mean bias relations for galaxies and halos (Uhlemann et al., 2018). Measuring the count-in-cells distributions of halos is a crucial first step to measure the differences of different dark energy models, before similar measurements are conducted with observed galaxy survey samples with greater uncertainties.

It is possible that the halo and galaxy distributions are very similar if each halo is usually occupied by one galaxy. Even so, tests over a wide range of scales and redshifts are necessary to confirm the validity of the analytical PDF for structure formation and clustering. If these distribution functions are roughly similar at low redshift, then high redshift samples could provide tighter constraints on their functional form. The effect of dark energy on the distribution functions of halo and galaxy clustering has yet to be explored. The analytical models of CiC statistics considered in this paper (Section 2.4) have not explicitly included evolving dark energy. To distinguish the difference between galaxy and halo distribution, the low mass halos should be included, although they may not host any galaxies. However, this may be difficult because the resolution of simulations across a wide range of halo masses is limited. The scale-dependence of halo bias (Dalal et al., 2008) and galaxy bias (Weinberg et al., 2004) indicate that the distribution functions of halos and galaxies due to clustering are scale-dependent. While the scale-dependence of the moments of the CiC distribution has been measured in surveys and simulations, the analytical form of the CiC PDF as a function of scale and redshift has not been extensively studied. The halo number density and halo CiC PDF can be influenced differently by different dark energy models, but difficult to measure due to limited simulation volume (e.g. Shi & Baugh, 2016). The CiC distributions for different equations of states with dark energy might show measurable differences for halos at certain mass ranges, length scales and redshifts.

The main goals of this paper are:

a) to demonstrate the utility of CiC for distinguishing different dark energy models;

b) to find the best-fit analytical model of the halo CiC PDF; and,

c) to make predictions of the scales and redshifts of CiC statistics that show the largest deviations from $\Lambda$CDM.

In Section 2, we describe our CiC algorithm and the analytical models used for model fitting. In Section 3, we give results for the residuals of CiC PDFs, for the moments of the CiC PDFs, and for the best-fit analytical models as functions of cell size and redshift for various dark energy models. In Section 4, we discuss how our method describes the differences in dark energy models and the sources of uncertainties. In Section 5, we summarize our findings.

CiC statistics contain more information than the two-point correlation function. The CiC PDF shows the ensemble average of clustering properties and combines information from a range of length scales up to the cell size (Wall & Jenkins, 2012). The PDF moments can be readily computed (in both two and three dimensions), and be related to the correlation functions (Peebles, 1980). It has been shown that the CiC statistics do not depend on cell shapes strongly (Saslaw & Crane, 1991; Szapudi & Colombi, 1996; Szapudi, 1998).

The mean of the CiC at a given redshift (e.g. galaxy cluster abundance) constrains models of structure formation and clustering (Mana et al., 2013). It depends on the evolution of the halo mass function (Tinker et al., 2008). For a Gaussian-distribution density field or Poisson-distributed halos or galaxies, i.e. without clustering, the mean contains complete information about the statistical properties (Wall & Jenkins, 2012). A clustered distribution has an increased variance and a non-zero skewness compared to a Poisson distribution. These indicate nonlinear gravitational clustering (Wall & Jenkins, 2012). The second order moment $\sigma_{8}$ is the matter density fluctuation at a scale of 8 Mpc/h. It has been readily measured for the CMB (Planck Collaboration et al., 2018) and used as a cosmological parameter in simulations (e.g. Klypin et al., 2011). To constrain cosmological models, higher order moments of the CiC as a function of redshift and scale are necessary.

The moments of the galaxy CiC PDF have been used to quantify galaxy clustering in the CFHTLS-Wide survey (Wolk et al., 2013) and the spectroscopic Sloan Digital Sky Survey (SDSS) (Bel & Marinoni, 2014). CiC has also been used to measure galaxy bias for submillimeter galaxies observed with the Atacama Large Millimeter Array (ALMA) (Ono et al., 2014), galaxy samples from the ALHAMBRA survey (López-Sanjuan et al., 2015), the VIPERS survey (Di Porto et al., 2016) and the Dark Energy Survey (DES) (Salvador et al., 2019). The third and fourth order moments of CiC have been measured for simulated galaxies in smoothed particle hydrodynamic simulations (Weinberg et al., 2004) and dark matter halos in different simulated cosmologies as a function of scale and redshift (Casas-Miranda et al., 2003; Angulo et al., 2008). All these CiC measurements only assumed the $\Lambda$CDM model, except for Bel & Marinoni (2014) who also included a varied $w$.

There are advantages in fitting an analytic form to the measured discrete CiC PDF rather than computing the first few moments numerically. An analytical form inherently contains moments of all orders and sheds light on the stochastic processes underlying the clustering statistics and it provides hints about the physical processes behind structure formation and clustering. The CiC PDF from galaxy surveys and N-body simulations have been compared with theoretical models of structure formation and clustering, but there is no clear consensus on the detailed functional form of the CiC PDF of galaxies, halos or dark matter densities, other than clear evidence for non-Poisson distributions. The analytical forms of CiC PDF models are described in Section 2.4. Some examples in the literature of fitting an analytic CiC PDF to survey data or simulations include studies using the 2dF Galaxy Redshift Survey (2dFGRS) and dark matter in the Hubble Volume simulation (Croton et al., 2004), the DEEP2 Galaxy Redshift Survey and mock galaxies from simulations (Conroy et al., 2005), SDSS galaxies (Yang & Saslaw, 2011; Hurtado-Gil et al., 2017), VIPERS galaxies (Bel et al., 2016), DES galaxies and MICE mock galaxies (Clerkin et al., 2017) and dark matter halos in other simulations (Neyrinck et al., 2014; Ahn et al., 2015).

In this study, we measure the CiC statistics of simulated dark matter halos in three cosmologies and find the best-fit parameters of model probability distribution functions in a range of length scales and redshifts. We make use of the publicly available Dark Energy Universe Simulations halo position catalogs (DEUS; Alimi et al., 2012). A set of three simulations in comoving space with box-lengths of 648 $h^{-1}$Mpc and $2048^{3}$ dark matter particles (Alimi et al., 2012) are selected for comparison between three dark energy models, namely, the cosmological constant ($\Lambda$CDM), Ratra-Peebles potential for quintessence scalar field (RPCDM) and phantom dark energy model ($w$CDM) with constant equation of state. Table 1 lists the cosmological parameters used for these simulations (see Bouillot et al., 2015, for details on the choice of cosmologies in DEUS). In summary, the cosmological parameters for RPCDM and $w$CDM were chosen by DEUS to be statistically consistent with the best-fit flat $\Lambda$CDM model to the WMAP-7yr data while being indistinguishable from $\Lambda$CDM at the significance level of the WMAP-7yr data. In our work we use the simulations as provided by DEUS. It is worth noting that the values of $\Omega_{m}$, $w_{0}$ and $w_{1}$ chosen at the time of DEUS are not fully consistent with recent constraints combining multiple observational probes (Abbott et al., 2019a, b). New simulations adopting more tightly constrained cosmologies would produce smaller differences in CiC between different dark energy models, but should still be distinguishable with a subset of similar CiC measurements. We plan to use a wider range of simulations that systematically vary the cosmological parameters in future work to test our proposed framework for distinguishing dark energy models. Here we use the halo catalogs from DEUS with smaller box-lengths than the DEUS-FUR, but with the same cosmological parameters. The smaller box-length and fewer number of particles allow a higher mass resolution so that galaxy-sized halos are included in the halo catalogs. Each halo contains at least 100 particles, found using the Friends-of-Friends algorithm (Roy et al., 2014), which gives a mass cut at $M_{halo}>2.4\times 10^{11}M_{\odot}$ for all three cosmologies. The halos are from 9 snapshots at redshifts spaced between $0\leqslant z\leqslant 4$.

The CiC distribution sampled here is the PDF of a discrete random variable, $N$, the number count of dark matter halos with $M_{halo}>2.4\times 10^{11}M_{\odot}$ inside spherical cells of a given radius. We make 216 CiC measurements using the 3 halo catalogs for 8 different cell sizes of radii $R=$ 2, 4, 6, 8, 10, 15, 20, 25 $h^{-1}$Mpc, at 9 different redshifts $z$ = 0, 0.1, 0.25, 0.4, 0.65, 1, 1.5, 2.3 and 4. Each side of the cubic simulation box is divided into 512 sections and the centers of the $512^{3}$ cubes are the centers of the $512^{3}$ spherical cells used for the CiC measurements (see Appendix A). Since the cells can be overlapping, our CiC is oversampling the simulation box. The oversampling method is inspired by the infinite sampling CiC method (Szapudi, 1998), which transforms counting discrete objects in sampling cells into calculating the ratios of overlapping area/volume for cells around discrete objects. The infinite sampling CiC method is argued to be free of measurement errors. The densely populated cell centers in our method can be viewed as sampling the overlapping regions in the infinite sampling CiC method. When the number of cells tends to infinity, the ratios of overlapping volumes are then exactly recovered. Our choice of the cell number is based on a resolution study (see Appendix A). The number of halos in each cell is found by comparing the cell radius and the distance between the cell center and halo center. The occurrence of each count $N$ is then counted for cells that are completely within the simulation box, excluding incomplete cells on the edges of the simulation box. Lastly, the histogram is normalized to produce a PDF, called the CiC distribution. The CiC algorithm is parallelized with Message Passing Interface (MPI; Gropp et al., 2014) in C and is run on the Blue Waters supercomputer at the National Center for Supercomputing Applications (NCSA).

Four models are fitted to the resulting CiC distributions. The first is the gravitational quasi-equilibrium distribution (GQED) given by

where $\overline{N}$ is the average number of halos in a given cell volume and $b$ is the ensemble average ratio of the gravitational correlation energy to twice the kinetic energy of peculiar velocities (Saslaw & Hamilton, 1984).

The second is the negative binomial distribution (NBD) derived by Elizalde & Gaztanaga (1992). We use a formulation with the Gamma function (Yang & Saslaw, 2011),

where $n$ is the average number of halos in a given cell volume and $g$ is a clustering parameter. The third model is a log-normal distribution, commonly used to describe the matter density distribution. The Poisson sampled log-normal distribution includes the shot noise as halos are discretely sampled from a log-normally distributed continuous matter density field (Clerkin et al., 2017),

where

$x=\mathrm{ln}(1+\delta_{g})+w^{2}/2$ and $w^{2}$ is the variance of the corresponding normal distribution $f[\mathrm{ln}(1+\delta_{g})]$. The fourth model is the log-normal distribution plus a halo bias parameter $b$ and a matter density variance parameter $C_{b}$ (Hurtado-Gil et al., 2017),

where

Appendix A shows two examples of the best-fit models for different cell radii and cosmologies. Of these four models, the GQED model has a physical basis in the thermodynamics (Saslaw, 2000, and earlier references) and statistical mechanics (Saslaw & Ahmad, 2010) of cosmological gravitational many-body systems. The physics behind the other three models seems more obscure. All four models are Lagrangian distributions. The NBD is essentially a Poisson sampled gamma distribution and the GQED is a Poisson sampled truncated Borel distribution (also called a convolved or compound distribution in some literature) (Saslaw, 2000). One interpretation for the Poisson component of these distribution functions is that the formation of galaxies or dark matter halos around galaxies out of the underlying density field is a Poisson process. The question then is what statistical distribution the underlying dark matter density field follows. Another interpretation is a halo model, where the distribution of dark matter halos that host galaxy clusters follows a Poisson distribution, and the distribution of galaxies in galaxy clusters can be modeled by a halo occupation distribution, essentially governed by the halo mass and mass profile as well as biasing effects (Sheth & Saslaw, 1994; Fry et al., 2011). With suitable values of the bias factors, this full halo model shows good agreement with the first few moments of the dark matter halos and subhalos CiC from an adaptive mesh refinement cosmological simulation of dark matter particles (Fry et al., 2011). However, the physics behind this phenomenological halo model is not clear. The assumption that the validity of the emergence of cluster-sized halos is a Poisson process requires verification.

We find the best-fit model parameter with the least squares weighted uniformly for all counts $N$

where $E=0.0005$ and $f_{CiC}(N)$ is the measured CiC PDF. The uniform weighting factor $E$ allows the cells with underdensity and overdensity to have the same weighting as the peak of the CiC distribution function. The least-squares is scaled larger with a small $E$ for easier comparison and smaller round-off errors during fitting. We use MPFIT¹¹1https://pages.physics.wisc.edu/$\sim$craigm/idl/cmpfit.html (Markwardt, 2009), a least squares fitting library in C using the Levenberg-Marquardt technique (Press et al., 2002), to iteratively search for the best-fit parameters and the least squares for all four models. Both PLN and PLNB are integral functions. The CQUAD doubly-adaptive integration in the GNU Scientific Library²²2https://www.gnu.org/software/gsl/doc/html/integration.html (Galassi et al., 2009) is used to evaluate the integral functions. Due to the presence of factorials of large integers as well as very small or large exponential functions in all four models, double precision floating-point numbers in C/C++ programming language are insufficient for evaluating our model over the desired range of parameters. We make use of the GNU MPFR library³³3https://www.mpfr.org/ (Fousse et al., 2007), a C library for multiple-precision floating-point computations with correct rounding, for function evaluations. At least 80-bits significand must be kept to ensure the smoothness of model functions, the success of integration and the convergence of the least squares fitting.

The jackknife method (Shao & Tu, 1995) is used to estimate the uncertainties of the CiC distributions for two test cases, cell radii $R=2h^{-1}$ Mpc and $R=10h^{-1}$ Mpc at $z=0$ in $\Lambda$CDM. The sample of $512^{3}$ cells are divided into $n$ equal volume subsamples by coordinates, then one subsample is deleted at a time to obtain a jackknife CiC subsample. For a given $n$, there are therefore $n$ different jackknife CiC subsamples, and every cell is deleted once. The jackknife error is given by

where

and $\delta_{i}$ denotes the rms error when all the cells except the $i$-th are used. The sum of the squared jackknife errors for a CiC distribution may vary with delete fraction. For the two test cases, we choose $n$ between 2 and 2048.

For two test cases, cells of radii $R=2h^{-1}$ Mpc and $R=10h^{-1}$ Mpc at $z=0$ in $\Lambda$CDM, the sums of jackknife errors for a given delete fraction are similar (Figure 1). Except when half or 1/4 of the cells are deleted, the sum of jackknife errors decreases as the delete fraction decreases. The jackknife errors are slightly larger than the uncertainties due to uncertain counts near the boundaries of cells in the resolution study, but are much smaller than the least squares of best-fit analytical models. The jackknife errors for the CiC distributions of other cell sizes, redshifts and cosmologies are expected to be similar and well below the deviations between the analytical models and CiC distributions. The origin of the uncertainties in our CiC measurement is discussed further in Section 4.

The measured CiC distributions $f(N)$ shown in Figure 2 and 3 are smooth unimodal probability distribution functions. In Figure 3, at a given redshift the peak of the distribution moves to larger halo counts as the cell radius increases, because larger cells contain more halos. On a log scale, $f(N)$ is less smooth in larger cells at higher redshifts. For a fixed cell radius (Figure 2), the peak of $f(N)$ shifts towards larger $N$ at higher redshifts, because the same physical volume encloses a larger comoving volume. The peak also decreases at higher redshifts, as the probability spreads into a wider range of $N$. The differences between distributions in different dark energy models in smaller cells and lower redshifts are smaller, but grow very significantly as a function of both cell radius and redshift.

In Figure 4 we examine the lowest-order statistical moments of the residuals between quintessence and $\Lambda$CDM and phantom dark energy and $\Lambda$CDM, denoted by $\triangle f(N)_{RP-\Lambda}=f(N)_{RPCDM}-f(N)_{\Lambda CDM}$ and $\triangle f(N)_{w-\Lambda}=f(N)_{wCDM}-f(N)_{\Lambda CDM}$ respectively. The mean, variance, skewness and kurtosis of the residuals of $f(N)$ between both alternative dark energy models and $\Lambda$CDM are larger in smaller cells. The mean of the residuals $\bar{\triangle}f(N)_{RP-\Lambda}$ and $\bar{\triangle}f(N)_{w-\Lambda}$ decrease at higher redshifts. The residual means are computed over all $N$ (from 0 to the largest N with non-zero count), so their values are vanishingly small and do not distinguish well between dark energy models due to this averaging. The variance of of the residuals $\triangle f(N)_{RP-\Lambda}$ is larger than that of $\triangle f(N)_{w-\Lambda}$ for all cell radii and redshifts. Generally, as redshift increases, the variance of the difference tends to first decrease and then increase, although this trend is not seen for all cell radii. The skewness of the residuals $\triangle f(N)_{RP-\Lambda}$ is generally smaller than that of $\triangle f(N)_{w-\Lambda}$. Both values approach zero at higher redshifts and in larger cells except for cell radius $R=2h^{-1}$Mpc. The kurtosis of the differences $\triangle f(N)_{RP-\Lambda}$ and $\triangle f(N)_{w-\Lambda}$ are not strongly distinguished by dark energy cosmological model; the kurtosis of these differences generally decreases for larger cells and higher redshifts. It increases very slightly from $z=2.3$ to $z=4$. Figure 4 shows that the variance, skewness and kurtosis of the residuals $\triangle f(N)_{RP-\Lambda}$ and $\triangle f(N)_{w-\Lambda}$ on $2h^{-1}$Mpc scale at $z>2.3$ are the largest and the most different among RPCDM, $w$CDM and $\Lambda$CDM. The skewness and kurtosis of the residuals on $4-6h^{-1}$Mpc scales are larger at $z<0.65$. At scales larger than $8h^{-1}$Mpc, the moments of the residuals are not significantly far from zero at $0<z<4$.

The percentage differences of the lowest moments of the measured CiC $f(N)$ in other dark energy models compared to the moment values obtained in $\Lambda$CDM are shown in Figure 5. Overall, RPCDM clearly shows larger amplitude and opposite signs in percentage differences than $w$CDM when compared with $\Lambda$CDM. The percentage difference of the CiC $f(N)$ mean between RPCDM and $\Lambda$CDM decreases almost linearly as a function of redshift regardless of cell radius from about -8.0% at $z=0$ to about -67.7% at $z=4$. The percentage difference of the CiC mean between $w$CDM and $\Lambda$CDM increases almost linearly as a function of redshift regardless of cell radius from about 3.9% at $z=0$ to about 25.5% at $z=4$. The colored lines for different cell radii overlap almost exactly in the upper left panel of Figure 5. This shows that the percentage difference of the CiC means between RPCDM/$w$CDM and $\Lambda$CDM depends only on redshift, not scale. The trend for variance is similar, but shows slightly larger percentage differences in smaller cells up to a maximum $\sim-80\%$ for RPCDM and about 30% for $w$CDM at $z=4$. The percentage difference of skewness as a function of redshift shows only smooth trends for cell radii up to 10$h^{-1}$Mpc, but has erratic and non-monotonic variation in larger cells. The percentage difference in kurtosis between RPCDM and $\Lambda$CDM increases up to more than 60% and that between $w$CDM and $\Lambda$CDM more than -10%. The percentage difference changes sign and crosses 0 at $z\sim 0.25$. Smaller cells show larger percentage differences in kurtosis. Overall, Figure 5 shows that in the redshift range of $0<z<4$, the mean, variance, skewness and kurtosis of the CiC PDFs are more promising in showing differences among RPCDM, $w$CDM and $\Lambda$CDM at high redshifts than at low redshifts. The percentage differences in mean and variance of the CiC PDFs do not depend strongly on scale between $2-25h^{-1}$Mpc. The percentage difference in skewness is a smooth function of redshift on $2-10h^{-1}$Mpc scales and increases as a function of scale. The percentage differences in kurtosis are more pronounced on small scales between $2-6h^{-1}$Mpc.

The parameters $\bar{N}$ and $b$ obtained from the best fit of the GQED model (Equation 1) to the measured CiC PDF $f(N)$ are plotted in Figure 6 as a function of cell radius and redshift; this Figure shows the fitted parameter values for the RPCDM and wCDM cosmologies and their percentage difference relative to the same fitted parameters for the $\Lambda$CDM cosmology. Comparable plots for best-fit NBD parameters $n$ and $g$ (Equation 2), PLN parameters $\bar{N}$ and $\omega$ (Equation 3), and the PLNB parameters ${\bar{N}}$, $b$, and $C_{b}$ (Equation 5) are shown in Figures 7, 8, and 9 respectively.

Examination of these figures shows that in all four fitted models, the parameters related to the mean halo count (either $\bar{N}$ or $n$) all show very similar trends in percentage differences with redshift relative to $\Lambda$CDM as the first-moment (mean) of the CiC $f(N)$ computed directly and shown in Figure 5. The best-fit $b$ in GQED increases and then decreases as a function of redshift, with larger $b$ in large cells; we remind the reader that the $b$ parameter in the GQED model (Equation 1) is a clustering rather than a bias parameter. The percentage difference of best-fit $b$ between RPCDM and $\Lambda$CDM grows positively but decreases strongly to negative amplitudes as redshift increases. In Figure 6 the smaller cells with $R=2h^{-1}$Mpc show more than -18% difference and the largest cells with $R=25h^{-1}$Mpc show about 5% difference at $z=4$. The percentage differences in $b$ between RPCDM and $\Lambda$CDM at $z<2.3$ and $w$CDM and $\Lambda$CDM at $z<4$ are all below 2%. Only RPCDM at very high redshift shows significant percentage difference compared to $\Lambda$CDM. The best-fit $g$ in the NBD is smaller for larger cells and decreases as a function of redshift. The percentage difference of $g$ between RPCDM and $\Lambda$CDM increases from about 8% to about 50-80% and from a few percent negative to more than -10% between $w$CDM and $\Lambda$CDM. Larger cells generally show larger percentage differences in $g$. The trend for $\omega$ in PLN is very similar with the percentage differences roughly halved. In PLNB, the best-fit bias parameter $b$ generally increases then decreases as a function of redshift. The bias and the percentage differences in bias between alternative dark energy models and $\Lambda$CDM is generally larger in larger cells. The percentage difference in bias in RPCDM and $w$CDM have opposite signs relative to $\Lambda$CDM, ranging in absolute magnitude from zero to 40%. The best-fit variance parameter $C_{b}$ in PLNB decreases as a function of increasing redshift and cell radius. RPCDM and $\Lambda$CDM show significant percentage difference (more than 100%) in the variance parameter $C_{b}$, whereas $w$CDM and $\Lambda$CDM show less than 50% absolute magnitude difference in $C_{b}$. The percentage differences for the bias and the variance parameters in PLNB are greater for large scales up to $25h^{-1}$Mpc at $1<z<4$.

The bias and the variance parameters in the PLNB fit were found to be effectively mathematically degenerate during model fitting. The MPFIT algorithm could not then converge to a global minimum for large cells at high redshifts. To confirm the degeneracy, we fixed $\bar{N}$ and evaluated the sums of squares for different bias and variance parameters on a grid. A two-dimensional grid search for $b$ and $C_{b}$ to locate the global minimum showed that the smallest sums of squares fell on a monotonic curve relating these two parameters. An arbitrarily small variance and a large bias (and vice versa) can therefore yield the same sum of squares. As a result, when the fitted bias was larger than 10 or reached an upper bound of 10, we took this as evidence of strong degeneracy and we used radial basis functions to perform linear and cubic interpolation over existing converged solutions to obtain the bias and variance parameters instead. The residuals of the linearly interpolated parameters were found to be smaller and to vary more uniformly. The linearly interpolated variance could be negative in this approach however, which is not physical, so we adopted the linearly interpolated bias $b$ and fixed its value when fitting the mean count $\overline{N}$ and variance $C_{b}$. The interpolated valued of $b$ and the fitted $C_{b}$ are included with other successful three-parameter fits in Figure 9. For the parameter range where degeneracy is not dominant we found the Levenberg-Marquardt algorithm to be robust to choices of initial parameter values. The cases of PLNB with linearly interpolated bias $b$ and fitted $C_{b}$ and $\bar{N}$ give the smallest least squares among all four analytical models considered.

Figures 10 and 11 show the scaled least squares (Equation 9) for all four models as a function of redshift, cell radius and dark energy model. The best-fit model for the halo CiC distribution is the NBD model for $R<6h^{-1}$Mpc and the PLNB model for $R>6h^{-1}$Mpc. The choice of dark energy model does not appear to influence significantly the choice of the best-fit analytical model. For $R=6h^{-1}$Mpc, the best-fit model transitions from NBD to PLNB for $z>0.25$. The quality of fits are generally the best at high redshifts for small cells $R<10h^{-1}$Mpc and at medium redshifts for large cells $R>15h^{-1}$Mpc. This improves monotonically as a function of redshift for $R>6h^{-1}$Mpc. For $R<6h^{-1}$Mpc, the scaled least squares have a peak as a function of redshift. In the cell radius range $2\leqslant R\leqslant 25h^{-1}$Mpc, the PLN model always has the lowest residual least squares compared to the other three models. For $R>8h^{-1}$Mpc, the ranking of the best-fit model is always PLNB, NBD, GQED, and PLN. The GQED ranks second for small cells $R<4h^{-1}$Mpc at low redshifts after NBD as the best-fit. The CiC distribution $f(N)$ is most accurately modeled by the GQED, NBD and PLN for cells $R=25h^{-1}$Mpc at $z=1$, $R=20h^{-1}$Mpc at $z=1.5$ and $R=15h^{-1}$Mpc at $z=2.3$. It is most accurately modeled by PLNB for cells $R=25h^{-1}$Mpc at $z=0.65$, $R=20h^{-1}$Mpc at $z=1$ and $R=15h^{-1}$Mpc at $z=2.3$. In Figure 11, for cell $R=2h^{-1}$Mpc, RPCDM is best fit by the four models, but for larger cells RPCDM produces worse fits than $\Lambda$CDM and $w$CDM. Only cell radius $R=2h^{-1}$Mpc shows significant differences in the quality of fits between various dark energy models. The differences in the quality of fits between various dark energy models are smaller than the differences in the quality of fits between the analytical models. The best-fits of models for small cells $R=2h^{-1}$Mpc are at low and high redshifts. The best-fits of models for cells $R=4-10h^{-1}$Mpc are at high redshifts. The best-fits of models for large cells $R=15-25h^{-1}$Mpc are at medium redshift ($1<z<2.3$). These results are discussed further in Section 4.

Consider the Friedmann equation for $\Lambda$CDM (Peacock, 1999)

where $H$ is the Hubble parameter, $\rho$ is the energy density and $\Lambda$ is the cosmological constant. The timescale for the dark matter halos to merge in an expanding universe can be approximated by the gravitational timescale

where $\Omega_{\Lambda}$ is the density parameter for the cosmological constant. Because $\Omega_{\Lambda}>0$, the gravitational timescale $\tau_{grav}$ should always be longer than the universe’s expansion timescale $\tau_{Hubble}=1/H$. In the thermodynamic description of galaxy clustering that gives rise to the GQED (Saslaw, 2000), the assumption for quasi-equilibrium evolution holds if the macroscopic timescale, $\tau_{macro}$, exceeds the microscopic timescale $\tau_{micro}\approx\frac{L}{v}$, where $L$ and $v$ are characteristic microscopic lengths and velocities respectively. Here

where $H$ is the Hubble parameter and $b$ is the clustering parameter in Equation 1. The microscopic timescale $\tau_{micro}\ll\tau_{Hubble}$ in the nonlinear regime and $\tau_{micro}\approx\tau_{Hubble}$ in the linear regime (Saslaw, 2000). The CiC distribution of dark matter halos is a result of halo formation, clustering, merging and the expansion of the universe. Merging halos reduces the number of halos and decreases the two-point correlation function of halos. The effect of merging is more prominent when the halo number density is higher at larger redshifts. The expansion of the universe also reduces the number density of halos and reduces halo correlation functions, but its effect is only dominant at low redshifts after the expansion accelerates. The value of $b$ in the GQED model increases as halos cluster gravitationally, but also decreases because of the competing effects of halo merging and the expansion of the universe. Assuming $\Omega_{\Lambda}=1-\Omega_{m}=1-0.2573=0.7427$, then $\tau_{macro}>\tau_{grav}$ if $b>0.746$ and $\tau_{macro}>\tau_{Hubble}$ if $b>0.5$. In Figure 6, we see $b>0.746$ for cell radii $R>6h^{-1}$Mpc. This means at scales larger than $6h^{-1}$Mpc, the degree of gravitational clustering for a many-body system in $\Lambda$CDM evolves slower than the expansion of the universe and slower than halo formation through merging. In contrast, $b<0.746$ for cell radius $R=2h^{-1}$Mpc at $0<z<4$. At this non-linear scale, local gravitational collapse and mergers happen faster than the global evolution of gravitational clustering. When a halo forms through merging, the number of halos changes for the whole thermodynamic system, so the system changes from one system in quasi-equilibrium to a new system with a different number of halos. Gravitational clustering takes longer to respond to abrupt changes due to mergers, so at small scales the deviation from equilibrium is larger. As a result, the GQED tends to fit less accurately at smaller scales than at large scales in Figure 10.

The CiC distribution measured by our oversampling CiC algorithm has small uncertainties from uncertain counts near cell boundaries. The small jackknife errors indicate that the CiC distribution function measurements are quite precise. The overlapping cells sample the same volume in the simulation box many times and produce very smooth CiC distributions. However, the overlapping cells are inherently correlated, so the deleted jackknife subsamples are not independent and may underestimate the errors. Also, excluding incomplete cells near the edges of the simulation box may introduce additional errors, because the simulation box is then no longer uniformly sampled. This edge effect is larger at high redshifts when a significant fraction ($\sim 40\%$) of cells are excluded. The fact that the resultant CiCs at high redshifts are still quite smooth means that the simulation box is large enough. By using the multi-precision library for model evaluation, the uncertainty due to evaluating complex functions should be negligible. The fitting parameters obtained are robust against varied initial guesses as long as the parameters are not degenerate. Only the PLNB model has degenerate parameters in certain parameter ranges and requires appropriate numerical treatment as described above.

Among the four analytical models we fitted to the measured CiC PDF $f(N)$, there is no model that fits the CiC best universally in all ranges of cell radii and redshifts. The CiC is best described by the NBD model for cells with $R<6h^{-1}$Mpc and by the PLNB model for cells $R>8h^{-1}$Mpc. Based on comparisons of the residual least squares (Figure 11), these individual dark energy models do not favor any particular analytical model of CiC. We find that the CiC $f(N)$ is best modeled for all four analytical models with large cells at medium redshifts. Figure 11 indicates that various dark energy models show significant differences only at a scale of $R=2h^{-1}$Mpc.

All these various models have ranges of the counts-in-cells distribution function where they more or less agree with the simulations and with each other, as judged by their least squares differences. Yet as Tufte (1983) discusses, statistical agreement, even (or especially) with many parameters, does not convey the full import of data. Therefore we would be cautious in claiming that small differences in distribution functions are meaningful.

To help overcome this problem it is useful to relate these statistics to the fundamental physical properties of halo clustering. If the statistics need to be modified, so do the physical properties. This is easier to understand in the case of the GQED since its physical origin is much less obscure than for the NBD, PLN and PLNB distributions. This may be because in the GQED only one force, gravity, dominates the evolution.

An important factor concerning the differences shown by different dark energy models lies in the cosmological parameters chosen for this set of simulations (Table 1). Bouillot et al. chose degenerate parameters the mean cosmic mater density, $\Omega_{m}$, the dark energy equation of state, $w$, and the root-mean-square fluctuation amplitude of the density contrast at 8 $h^{-1}$ Mpc, $\sigma_{8}$ within the marginalized $\sim$ 1 $\sigma$ contour that fits the cosmic microwave background anisotropy power spectra and the luminosity distances measurements to supernova Type Ia (Bouillot et al., 2015). The Bouillot et al. simulation aims to break the degeneracy between these cosmological parameters by testing non-linear clustering. The accelerated expansion in $\Lambda$CDM is stronger than in RPCDM, so to fit the observations as well as in the $\Lambda$CDM, the smaller acceleration in RPCDM is compensated by a larger dark energy density than is required. This is reflected in the smaller matter density, $\Omega_{m}$, in RPCDM (Alimi et al., 2010). The opposite is true for $w$CDM. As a result, the CiC mean, as well as the mean number count in the fitted analytical models, all show fewer halos in RPCDM and more halos in $w$CDM throughout their cosmic history.

Figure 3 shows that the CiC distributions measured with the same physical cell radius R have an increasing average halo count and an increasing most probable halo count from $z=0$ to $z=4$ for all three cosmologies. This is mainly because the matter density at a higher redshift is higher and a spherical cell with the same physical volume encloses a larger comoving volume at a higher redshift. In the same comoving volume, more bound halos should have formed at a lower redshift, but this increase in the number of halos is clearly offset by the expansion of the universe at $z\leqslant 4$. At sufficiently high $z$, such as $z>4$, large halos above our minimum mass criterion are very few, hence reversing the trend of increasing halo count with increasing redshift. The differences in halo counts between RPCDM or $w$CDM and $\Lambda$CDM are larger at higher redshifts.

Apart from the effect of various dark energy models on the number density of halos through changing the expansion history, a more subtle effect on the degree of clustering as a function of redshift was also shown in the halo CiC distributions. The best-fit GQED $b$ parameter as an indicator of the degree of clustering is found to be closer to unity for larger cells. Because $b$ is an integral of the two-point correlation function within the scale of the cell, the values of $b$ closer to 1 in larger cells are expected. The absolute difference of the best-fit $b$ between $w$CDM and $\Lambda$CDM is very small across all redshifts and cell radii. The largest percentage difference of $b$ between $w$CDM and $\Lambda$CDM is found to be about 2% with a cell radius of $2h^{-1}$Mpc at $z=4$. Therefore, $b$ is not very useful for distinguishing $w$CDM and $\Lambda$CDM. The difference in $b$ between RPCDM and $\Lambda$CDM is more obvious. The absolute difference between $b_{RPCDM}$ and $b_{\Lambda CDM}$ is larger with smaller cells and at higher redshifts. It is up to about 18%. It is interesting to note that the percentage differences between $b_{RPCDM}$ or $b_{wCDM}$ compared to $b_{\Lambda CDM}$ change sign at $z\sim 1.7$, (except for RPCDM with $R=2h^{-1}$Mpc cells). In other words, at the turning point $z\sim 1.7$, all three dark energy models have nearly identical best-fit values of $b$ for all CiC with $R=4$ to $25h^{-1}$Mpc. At lower $z<1.7$, $b_{RPCDM}$ is most distinguishable at $z=0.65$ in cells $R=4h^{-1}$Mpc with a percentage difference of about 2%. In Figure 6 it is worth noting that the turning point for $b$ occurs at a higher redshift for $w$CDM and at a lower redshift for RPCDM compared to $\Lambda$CDM. Our method makes a clear prediction for the redshift at which $b$ reaches the maximum for different dark energy models. CiC in galaxy surveys should be able to test our prediction.

The NBD parameter $g$ can be interpreted as another clustering parameter that is approximately equivalent to the two-point correlation (Yang & Saslaw, 2011). The absolute difference between $g_{RPCDM}$ and $g_{\Lambda CDM}$ is larger than the absolute difference between $g_{wCDM}$ and $g_{\Lambda CDM}$ for a given cell radius and redshift. Both tend to be more significant for smaller cell radii across all redshifts. In the last panel of Figure 7 the percentage differences of $g_{RPCDM}$ and $g_{wCDM}$ compared to $g_{\Lambda CDM}$ are both larger for large cell radii and increase to their maxima at $z\sim 2.3$ to about 48% difference for RPCDM and 18% for $w$CDM. The substantial differences in the evolution history of clustering indicated by $b$ and $g$ as a result of various dark energy models suggest that the clustering parameters as a function of redshift are promising probes for distinguishing RPCDM and $w$CDM from $\Lambda$CDM.

The comparison between PLN and PLNB clearly shows that if the matter density is modeled by a log-normal distribution, bias must be considered in order to model the CiC accurately. The bias parameter in PLNB does not take the scale dependence of bias and nonlinear bias into account, so a better model of bias would be necessary to improve the accuracy of the PLN or the PLNB. The PLNB, being the best-fit model to the CiC in large cells, shows about 40% differences in $b$ and larger than 100% differences in $C_{b}$ between RPCDM and $\Lambda$CDM. The variance $C_{b}$ in the PLNB increases from high redshift to low redshift as matter density fluctuations grow. The bias parameter clearly shows a scale dependence, but does not vary much as a function of redshift at small scales. The bias at large scales tends to decrease towards lower redshifts. The PLN and the PLNB both show larger percentage differences in larger cells at higher redshifts. The trends of the parameters in PLNB are not as smooth as in the other models, which may make clearer predictions.

Our CiC methods and model fitting can be applied to current and future galaxy surveys. The trends of CiC and their best-fit models are expected to be similar for galaxies, but whether the large percentage differences would be measured in the CiC of galaxies might be less certain for a number of reasons. Firstly, the galaxy surveys may not have large enough samples of galaxies at $z>2.3$ or $z>4$ to measure the very large percentage differences between dark energy models precisely. At very high redshifts, the galaxy samples are not complete samples because very faint galaxies are below detection limits.

Secondly, the percentage differences between dark energy models may be smaller when cylindrical cells are used for galaxy surveys, which are subject to photometric redshift uncertainties and redshift-space distortions. These cylindrical cells are circular cells elongated in the redshift direction. The CiC measured with the cylindrical cells are redshift-averaged, whereas our measurements of simulations are snapshots at particular redshifts. CiC measurements done in 3D cells and 2D projected cells are not directly comparable. We will explore 2D CiC in redshift-space in further detail in future work. The redshift-averaged CiC will only be close to our simulation measurements if the redshift thickness is small. However, the galaxies in a thin redshift bin may not be a large enough sample to give reliable results. Furthermore, there are larger uncertainties associated the photometric redshifts that make galaxy samples in very thin redshift bins even more uncertain. The count of objects per steradian per unit increment in redshift can be approximated by $\frac{dN}{dz}=\frac{n_{0}}{H_{0}}\frac{r^{2}(z)}{\sqrt{\Omega_{m}(1+z)^{3}+\Omega_{\Lambda}}}$, where $r(z)=\frac{1}{H_{0}}\int_{0}^{z}\frac{dz^{\prime}}{\sqrt{\Omega_{m}(1+z^{\prime})^{3}+\Omega_{\Lambda}}}$ (Peebles, 1993). Assuming that photometric redshift uncertainties and redshift-space distortions produce $\sigma_{z}=0.002(1+z)$, the fractional galaxy count uncertainty $dN/N$ for a redshift bin $0.4<z<0.6$ is about 2.8%. Assuming a larger $\sigma_{z}=0.02(1+z)$ for a redshift bin $0.8<z<1.2$, $dN/N\approx 20\%$. Assuming a small $\sigma_{z}=0.002(1+z)$ at high redshift $3.8<z<4.2$, $dN/N\approx 5\%$. The relative uncertainty in the galaxy count for the DES Science Verification catalog at $0.1<z<0.5$ is about 5% (Clerkin et al., 2017). The requirement for Year 10 LSST large-scale structure analysis is that the systematic uncertainty in the mean redshift of each tomographic bin shall not exceed $0.003(1+z)$ and the systematic uncertainty in the photometric redshift scatter $\sigma_{z}$ shall not exceed $0.03(1+z)$ (The LSST Dark Energy Science Collaboration et al., 2018). The average number of galaxies in our fitting models should have similar fractional uncertainty as the above values of $dN/N$. The uncertainties and percentage differences of the clustering parameters of our fitting models for different cosmologies must be found by fitting the CiC of galaxies. Our predicted large differences between dark energy models at $z>1$ should still be detectable if the galaxy counts have 5% uncertainties due to photometric redshift uncertainty. Further redshift-space CiC studies with simulations that sample the cosmological parameter space more densely would better forecast the detectability of the differences from different dark energy models in future surveys and are a target of our future research.

Lastly, the process of removing satellite galaxies may introduce more uncertainties. Once a halo merges into a larger halo, its count in our measurement will become zero. This means that for galaxy CiC, satellite galaxies in galaxy clusters that reside in smaller halos that previously merged into larger halos would not be counted. From a theoretical standpoint, not considering subhalos or satellite galaxies is easier, because the stochasticity of the halo occupation number is no longer a concern. However, for observations, it is harder to decide if a galaxy belongs to a galaxy cluster, especially without the spectroscopic redshifts of the galaxies in the cluster.

CiC for galaxy surveys do not need accurate mass determination, especially when compared to the CiC distribution of galaxy-sized halos in simulations, because the observable galaxies are above the minimum observable mass for galaxy-mass halos and dwarf galaxies. The large sample of galaxies and the easily obtained discrete number of galaxies make CiC measurements easier to implement than other probes of dark energy mentioned in Section 1.