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 a constant Λ and w. Scalar-field models introduce quintessence (see Ratra & Peebles 1988; Wetterich 1988; Frieman et al. 1995; Copeland et al. 2006, for reviews) and allow a time-varying equation of state parameterized by a Taylor expansion w(a) = w₀ + (1 − a)w_a, where a is the scale factor (Planck Collaboration et al. 2018). A cosmological constant Λ corresponds to w₀ = −1 and w_a = 0. A constant Λ and w = −1 is consistent with CMB+BAO+SNe+H₀ measurements within approximately 10% uncertainty with 95% confidence (Planck Collaboration et al. 2018). However, w₀ and w_a for a dynamical dark energy equation of state are less constrained by combinations of different probes (Figure 4 in Planck Collaboration et al. 2016b). The case where 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 parameterized by the slope parameter _s, whose negative values correspond to phantom models. The value of _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 (Munshi et al. 2011; Smith & Zaldarriaga 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 nonrandom 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., the mean, of the CiC distribution with a high mass cut at ∼10¹⁴ M_⊙. 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 CiC distributions of halos is a crucial first step to measuring 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 (Benson et al. 2000; Weinberg et al. 2004) indicates 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 Λ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 nonzero skewness compared to a Poisson distribution. These indicate nonlinear gravitational clustering (Wall & Jenkins 2012). The second-order moment σ₈ 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 (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 Λ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 PDFs 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⁻¹ Mpc and 2048³ dark matter particles (Alimi et al. 2012) are selected for comparison between three dark energy models, namely, the cosmological constant (ΛCDM), Ratra–Peebles potential for quintessence scalar field (RPCDM), and phantom dark energy model (wCDM) with a 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 wCDM were chosen by DEUS to be statistically consistent with the best-fit flat ΛCDM model to the WMAP-7 yr data while being indistinguishable from ΛCDM at the significance level of the WMAP-7 yr data. In our work, we use the simulations as provided by DEUS. It is worth noting that the values of Ω_m, w₀, and w₁ chosen at the time of DEUS are not fully consistent with recent constraints combining multiple observational probes (Abbott et al. 2019a, 2019b). 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 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 × 10¹¹ M_⊙ for all three cosmologies. The halos are from nine 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 × 10¹¹ M_⊙ inside spherical cells of a given radius. We make 216 CiC measurements using the three halo catalogs for eight different cell sizes of radii R = 2, 4, 6, 8, 10, 15, 20, 25 h⁻¹ Mpc, at nine 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³ cubes are the centers of the 512³ spherical cells used for the CiC measurements (see the Appendix). Because 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 the Appendix). 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 a 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.

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

$$\begin{eqnarray}\begin{array}{rcl}{f}_{\mathrm{GQED}}(N) & = & \displaystyle \frac{\overline{N}(1-b)}{N!}\\ & & \times \,{\left(\overline{N}(1-b)+{Nb}\right)}^{N-1}{e}^{-\overline{N}(1-b)-{Nb}},\end{array}\end{eqnarray} \tag{ 1 }$$

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),

$$\begin{eqnarray}&&{f}_{\mathrm{NBD}}(N)=\displaystyle \frac{{\rm{\Gamma }}\left(N+\tfrac{1}{g}\right)}{{\rm{\Gamma }}\left(\tfrac{1}{g}\right)N!}\displaystyle \frac{{n}^{N}{\left(\tfrac{1}{g}\right)}^{\tfrac{1}{g}}}{{\left(n+\tfrac{1}{g}\right)}^{N+\tfrac{1}{g}}},\end{eqnarray} \tag{ 2 }$$

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

$$\begin{eqnarray}&&{f}_{\mathrm{PLN}}(N)={\int }_{-1}^{\infty }\displaystyle \frac{{\bar{N}}^{N}{\left(1+{\delta }_{g}\right)}^{N}}{N!}{e}^{-\bar{N}(1+{\delta }_{g})}f({\delta }_{g})d{\delta }_{g},\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray}&&f({\delta }_{g})d{\delta }_{g}=\displaystyle \frac{1}{w\sqrt{2\pi }}\exp \left(\displaystyle \frac{-{x}^{2}}{2{w}^{2}}\right){dx},\end{eqnarray} \tag{ 4 }$$

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

$$\begin{eqnarray}&&{f}_{\mathrm{PLNB}}({\rm{\Delta }})=\displaystyle \frac{1}{\sqrt{2\pi {H}_{0}}}\displaystyle \frac{\exp (-\tfrac{1}{2}\tfrac{{y}^{2}}{{H}_{0}})}{{\rm{\Delta }}+b-1},\end{eqnarray} \tag{ 5 }$$

where

$$\begin{eqnarray}&&{\rm{\Delta }}=N/\bar{N},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{H}_{0}=\mathrm{log}(1+{C}_{b}),\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&y=\mathrm{log}\left({\rm{\Delta }}+b-1)\displaystyle \frac{\sqrt{1+{C}_{b}}}{b}\right),\end{eqnarray} \tag{ 8 }$$

The Appendix 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 halo and subhalo 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:

$$\begin{eqnarray}&&x=\sum _{N=0}^{{N}_{\max }}{\left(\displaystyle \frac{{f}_{\mathrm{model}}(N)-{f}_{\mathrm{CiC}}(N)}{E}\right)}^{2},\end{eqnarray} \tag{ 9 }$$

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 roundoff errors during fitting. We use MPFIT³ (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⁴ (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⁵ (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⁻¹ Mpc and R = 10h⁻¹ Mpc at z = 0 in ΛCDM. The sample of 512³ 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

$$\begin{eqnarray}&&{\sigma }_{\mathrm{JK}}=\sqrt{\displaystyle \frac{n-1}{n}\sum _{i=1}^{n}{\left({\delta }_{i}-{\delta }_{\mathrm{JK}}\right)}^{2}},\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray}&&{\delta }_{\mathrm{JK}}=\displaystyle \frac{1}{n}{\delta }_{i}\end{eqnarray} \tag{ 11 }$$

and δ_i denotes the rms error when all the cells except the ith are used. The sum of the squared jackknife errors for a CiC distribution may vary with the delete fraction. For the two test cases, we choose n between 2 and 2048.

For two test cases, cells of radii R = 2h⁻¹ Mpc and R = 10h⁻¹ Mpc at z = 0 in ΛCDM, the sums of jackknife errors for a given delete fraction are similar (Figure 1). Except when half or a quarter 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.

Zoom In Zoom Out Reset image size

The measured CiC distributions f(N) shown in Figures 2 and 3 are smooth unimodal PDFs. 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 toward 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 4, we examine the lowest-order statistical moments of the residuals between quintessence and ΛCDM and phantom dark energy and ΛCDM, denoted by $\bigtriangleup f{(N)}_{\mathrm{RP}-{\rm{\Lambda }}}\,=f{(N)}_{\mathrm{RPCDM}}-f{(N)}_{{\rm{\Lambda }}\mathrm{CDM}}$ and $\bigtriangleup f{(N)}_{w-{\rm{\Lambda }}}=f{(N)}_{\mathrm{wCDM}}\,-f{(N)}_{{\rm{\Lambda }}\mathrm{CDM}}$ respectively. The mean, variance, skewness, and kurtosis of the residuals of f(N) between both alternative dark energy models and ΛCDM are larger in smaller cells. The mean of the residuals $\bar{\bigtriangleup }f{(N)}_{\mathrm{RP}-{\rm{\Lambda }}}$ and $\bar{\bigtriangleup }f{(N)}_{w-{\rm{\Lambda }}}$ decrease at higher redshifts. The residual means are computed over all N (from 0 to the largest N with nonzero count), so their values are vanishingly small and do not distinguish well between dark energy models due to this averaging. The variance of the residuals $\bigtriangleup f{(N)}_{\mathrm{RP}-{\rm{\Lambda }}}$ is larger than that of $\bigtriangleup f{(N)}_{w-{\rm{\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 $\bigtriangleup f{(N)}_{\mathrm{RP}-{\rm{\Lambda }}}$ is generally smaller than that of $\bigtriangleup f{(N)}_{w-{\rm{\Lambda }}}$. Both values approach zero at higher redshifts and in larger cells except for cell radius R = 2 h⁻¹ Mpc. The kurtosis of the differences $\bigtriangleup f{(N)}_{\mathrm{RP}-{\rm{\Lambda }}}$ and $\bigtriangleup f{(N)}_{w-{\rm{\Lambda }}}$ is not strongly distinguished by the 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 $\bigtriangleup f{(N)}_{\mathrm{RP}-{\rm{\Lambda }}}$ and $\bigtriangleup f{(N)}_{w-{\rm{\Lambda }}}$ on 2 h⁻¹ Mpc scale at z > 2.3 are the largest and the most different among RPCDM, wCDM, and ΛCDM. The skewness and kurtosis of the residuals on 4–6 h⁻¹ Mpc scales are larger at z < 0.65. At scales larger than 8 h⁻¹ Mpc, the moments of the residuals are not significantly far from zero at 0 < z < 4.

Zoom In Zoom Out Reset image size

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 ΛCDM are shown in Figure 5. Overall, RPCDM clearly shows a larger amplitude and opposite signs in percentage differences than wCDM when compared with ΛCDM. The percentage difference of the CiC f(N) mean between RPCDM and Λ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 wCDM and Λ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/wCDM and Λ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 of ∼−80% for RPCDM and about 30% for wCDM at z = 4. The percentage difference of skewness as a function of redshift shows only smooth trends for cell radii up to 10 h⁻¹ Mpc, but has erratic and non-monotonic variation in larger cells. The percentage difference in kurtosis between RPCDM and ΛCDM increases up to more than 60% and that between wCDM and ΛCDM more than −10%. The percentage difference changes sign and crosses 0 at z ∼ 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, wCDM, and ΛCDM at high redshifts than at low redshifts. The percentage differences in the mean and variance of the CiC PDFs do not depend strongly on scales between 2 and 25 h⁻¹ Mpc. The percentage difference in skewness is a smooth function of redshift on 2–10 h⁻¹ Mpc scales and increases as a function of scale. The percentage differences in kurtosis are more pronounced on small scales between 2 and 6 h⁻¹ Mpc.

Zoom In Zoom Out Reset image size

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 ΛCDM cosmology. Comparable plots for best-fit NBD parameters n and g (Equation (2)), PLN parameters $\bar{N}$ and ω (Equation (3)), and the PLNB parameters $\bar{N}$, b, and C_b (Equation (5)) are shown in Figures 7–9, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 Λ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 a 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 the best-fit b between RPCDM and ΛCDM grows positively but decreases strongly to negative amplitudes as redshift increases. In Figure 6, the smaller cells with R = 2 h⁻¹ Mpc show more than a −18% difference and the largest cells with R = 25 h⁻¹ Mpc show about a 5% difference at z = 4. The percentage differences in b between RPCDM and ΛCDM at z < 2.3 and wCDM and ΛCDM at z < 4 are all below 2%. Only RPCDM at the very high redshift shows a significant percentage difference compared to Λ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 ΛCDM increases from about 8% to about 50%–80% and from a few percent negative to more than −10% between wCDM and ΛCDM. Larger cells generally show larger percentage differences in g. The trend for ω 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 ΛCDM is generally larger in larger cells. The percentage differences in bias in RPCDM and wCDM have opposite signs relative to Λ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 ΛCDM show a significant percentage difference (more than 100%) in the variance parameter C_b, whereas wCDM and Λ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 25 h⁻¹ 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 values 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 < 6 h⁻¹ Mpc and the PLNB model for R > 6 h⁻¹ Mpc. The choice of dark energy model does not appear to influence significantly the choice of the best-fit analytical model. For R = 6 h⁻¹ Mpc, the best-fit model transitions from NBD to PLNB for z > 0.25. The quality of fits is generally the best at high redshifts for small cells R < 10 h⁻¹ Mpc and at medium redshifts for large cells R > 15 h⁻¹ Mpc. This improves monotonically as a function of redshift for R > 6 h⁻¹ Mpc. For R < 6 h⁻¹ Mpc, the scaled least squares have a peak as a function of redshift. In the cell radius range 2 ≤ R ≤ 25 h⁻¹ Mpc, the PLN model always has the lowest residual least squares compared to the other three models. For R > 8 h⁻¹ Mpc, the ranking of the best-fit model is always PLNB, NBD, GQED, and PLN. The GQED ranks second for small cells R < 4 h⁻¹ 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 = 25 h⁻¹ Mpc at z = 1, R = 20 h⁻¹ Mpc at z = 1.5, and R = 15 h⁻¹ Mpc at z = 2.3. It is most accurately modeled by PLNB for cells R = 25 h⁻¹ Mpc at z = 0.65, R = 20 h⁻¹ Mpc at z = 1, and R = 15 h⁻¹ Mpc at z = 2.3. In Figure 11, for cell R = 2 h⁻¹ Mpc, RPCDM is best fit by the four models, but for larger cells, RPCDM produces worse fits than ΛCDM and wCDM. Only cell radius R = 2 h⁻¹ 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 = 2 h⁻¹ Mpc are at low and high redshifts. The best fits of models for cells R = 4–10 h⁻¹ Mpc are at high redshifts. The best fits of models for large cells R = 15–25 h⁻¹ Mpc are at medium redshift (1 < z < 2.3). These results are discussed further in Section 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Consider the Friedmann equation for ΛCDM (Peacock 1999),

$$\begin{eqnarray}&&{H}^{2}=\displaystyle \frac{8\pi G}{3}\rho +\displaystyle \frac{{\rm{\Lambda }}{c}^{2}}{3},\end{eqnarray} \tag{ 12 }$$

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

$$\begin{eqnarray}&&{\tau }_{\mathrm{grav}}^{-2}\equiv \displaystyle \frac{8\pi G}{3}\rho ={H}^{2}-\displaystyle \frac{{\rm{\Lambda }}{c}^{2}}{3}={H}^{2}(1-{{\rm{\Omega }}}_{{\rm{\Lambda }}}),\end{eqnarray} \tag{ 13 }$$

where Ω_Λ is the density parameter for the cosmological constant. Because Ω_Λ > 0, the gravitational timescale τ_grav should always be longer than the universe's expansion timescale τ_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, τ_macro, exceeds the microscopic timescale ${\tau }_{{\rm{micro}}}\approx L/v$, where L and v are characteristic microscopic lengths and velocities respectively. Here,

$$\begin{eqnarray}&&{\tau }_{\mathrm{macro}}^{-1}=2H(1-b),\end{eqnarray} \tag{ 14 }$$

where H is the Hubble parameter and b is the clustering parameter in Equation (1). The microscopic timescale τ_micro ≪ τ_Hubble in the nonlinear regime and τ_micro ≈ τ_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 Ω_Λ = 1 − Ω_m = 1 − 0.2573 = 0.7427, then ${\tau }_{\mathrm{macro}}\gt {\tau }_{\mathrm{grav}}$ if b > 0.746 and τ_macro > τ_Hubble if b > 0.5. In Figure 6, we see b > 0.746 for cell radii R > 6 h⁻¹ Mpc. This means at scales larger than 6 h⁻¹ Mpc, the degree of gravitational clustering for a many-body system in Λ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 = 2 h⁻¹ Mpc at 0 < z < 4. At this nonlinear 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 entire 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 (∼40%) of cells is excluded. The fact that the resultant CiC at high redshifts are still quite smooth means that the simulation box is large enough. By using the multiprecision 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 < 6 h⁻¹ Mpc and by the PLNB model for cells with R > 8 h⁻¹ 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 = 2 h⁻¹ Mpc.

All these various models have ranges of the CiC distribution function where they more or less agree with the simulations and with each other, as judged by their least-squares differences. As Tufte (1983) discusses, however, statistical agreement, even (or especially) with many parameters, does not convey the full import of data. Therefore, we should 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 as 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 the degenerate parameters the mean cosmic matter density, Ω_m; the dark energy equation of state, w; and the rms fluctuation amplitude of the density contrast at 8 h⁻¹ Mpc, σ₈ within the marginalized ∼1σ contour that fits the CMB anisotropy power spectra and the luminosity distance measurements to SNe Ia (Bouillot et al. 2015). The Bouillot et al. simulation aims to break the degeneracy between these cosmological parameters by testing nonlinear clustering. The accelerated expansion in ΛCDM is stronger than in RPCDM, so to fit the observations as well as in the Λ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, Ω_m, in RPCDM (Alimi et al. 2010). The opposite is true for wCDM. 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 wCDM 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 wCDM and Λ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 wCDM and ΛCDM is very small across all redshifts and cell radii. The largest percentage difference of b between wCDM and ΛCDM is found to be about 2% with a cell radius of 2 h⁻¹ Mpc at z = 4. Therefore, b is not very useful for distinguishing between wCDM and ΛCDM. The difference in b between RPCDM and ΛCDM is more obvious. The absolute difference between b_RPCDM and b_Λ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}_{{\rm{\Lambda }}\mathrm{CDM}}$ change sign at z ∼ 1.7, (except for RPCDM with R = 2 h⁻¹ Mpc cells). In other words, at the turning point z ∼ 1.7, all three dark energy models have nearly identical best-fit values of b for all CiC with R = 4 to 25 h⁻¹ Mpc. At lower z < 1.7, b_RPCDM is most distinguishable at z = 0.65 in cells R = 4 h⁻¹ 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 wCDM and at a lower redshift for RPCDM compared to ΛCDM. Our method makes a clear prediction for the redshift at which b reaches the maximum for different dark energy models. The 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_ΛCDM is larger than the absolute difference between g_wCDM and g_Λ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_ΛCDM are both larger for large cell radii and increase to their maxima at z ∼ 2.3 to about 48% difference for RPCDM and 18% for wCDM. 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 wCDM from ΛCDM.

The comparison between PLN and PLNB clearly shows that if the matter density is modeled by a lognormal 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 Λ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 toward 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. First, 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.

Second, 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 with 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 $\tfrac{{dN}}{{dz}}=\tfrac{{n}_{0}}{{H}_{0}}\tfrac{{r}^{2}(z)}{\sqrt{{{\rm{\Omega }}}_{m}{\left(1+z\right)}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}}}$, where $r(z)=\tfrac{1}{{H}_{0}}{\int }_{0}^{z}\tfrac{{{dz}}^{{\prime} }}{\sqrt{{{\rm{\Omega }}}_{m}{\left(1+{z}^{{\prime} }\right)}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}}}$ (Peebles 1993). Assuming that photometric redshift uncertainties and redshift-space distortions produce σ_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 σ_z = 0.02(1 + z) for a redshift bin 0.8 < z < 1.2, dN/N ≈ 20%. Assuming a small σ_z = 0.002(1 + z) at high redshifts 3.8 < z < 4.2, dN/N ≈ 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 σ_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 to 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.

The 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.