Massive Dark Matter Halos at High Redshift: Implications for Observations in the JWST Era

Throughout this paper, we use ELUCID (Wang et al., 2016), a constrained N-body simulation obtained using the code L-Gadget, a memory-optimized version of Gadget-2 (Springel, 2005). The simulation uses a WMAP5 (Dunkley et al., 2009) cosmology with the following parameters: $\Omega_{\rm K,0}=0$, $\Omega_{\rm M,0}=0.258$, $\Omega_{\rm B,0}=0.044$, $\Omega_{\rm\Lambda,0}=0.742$, $H_{0}=100\ h\,{\rm km\,s^{-1}\,Mpc^{-1}}$ with $h=0.72$, and a spectral index of $n=0.96$ with an amplitude specified by $\sigma_{8}=0.80$ for the Gaussian initial density field. A total of 100 snapshots, from redshift $z=18.4$ to $0$, are saved. Halos are identified using the friends-of-friends (FoF) algorithm (Davis et al., 1985) with a scaled linking length of $0.2$. Subhalos are identified using the Subfind algorithm (Springel et al., 2001; Dolag et al., 2009), and subhalo merger trees are constructed using the SubLink algorithm (Springel, 2005; Boylan-Kolchin et al., 2009). ELUCID has a simulation box with a side length of $500\,h^{-1}{\rm{Mpc}}$ and uses a total of $3072^{3}$ particles to trace the cosmic density field. The mass of each dark matter particle is $3.08\times 10^{8}\,h^{-1}{\rm M_{\odot}}$, and the mass resolution limit of FoF halos is about $10^{10}\,h^{-1}{\rm M_{\odot}}$. As we only use halos well above this limit to model massive galaxies, the relatively low resolution of ELUCID does not affect our conclusions. Additionally, because all the uncertainties considered here are significant, the slight difference in cosmological parameters between WMAP5 and the more recent Planck data does not have a significant impact on our findings.

The initial condition of ELUCID is reconstructed from real galaxy groups identified from the SDSS (Yang et al., 2007; Yang et al., 2012) by a sequence of numerical methods detailed in Wang et al. (2016). Analysis based on cross-matching the simulated and observed halos at $z\sim 0$ showed that more than $95\%$ of massive halos with halo mass $M_{\rm halo}\geqslant 10^{14}\,h^{-1}{\rm M_{\odot}}$ can be recovered in the constrained simulation with a $0.5\,{\rm dex}$ tolerance on the error of halo mass and a $4\,h^{-1}{\rm{Mpc}}$ tolerance on the error of spatial location. These features of ELUCID enable not only the statistical study of halo and galaxy populations within it, but also a halo-by-halo comparison to massive clusters in the real universe. In §3, we use ELUCID to statistically infer uncertainties in measuring the abundance of extremely massive galaxies at $z=7\sim 10$, and in §4, we use ELUCID to study the assembly history of several real massive clusters in the local universe. The bottom-right panel of Fig. 1 shows some examples of the matched clusters in a slice of the simulation volume, including two massive clusters, Coma and Leo, and a number of Abell clusters (Abell et al., 1989; Struble & Rood, 1999).

Here we take ELUCID as the input, populate dark matter halos with galaxies and estimate the uncertainties in counting galaxies for surveys at high $z$. Our goal is to estimate the uncertainty of a given summary statistic, $s_{*}$, of galaxies in a sub-volume whose geometry is consistent with the survey in question, starting from the sample $S_{\rm halo}$ of all halos at a given redshift $z$ in ELUCID. To achieve this, we decompose the transformation from $S_{\rm halo}$ to $s_{*}$ into a chain of four steps based on our understanding of galaxy formation in the $\Lambda$CDM paradigm:

where $\mathcal{T}_{\rm CV}$ denotes the volume-based sub-sampling of halos, $\mathcal{T}_{M_{\rm halo}}$ denotes the determination of the mass $M_{\rm halo}$ for each halo in the sub-sample, $\mathcal{T}_{M_{*}}$ denotes the halo-to-galaxy mapping, and $\mathcal{T}_{\rm stat}$ denotes the statistical function that extracts the summary statistic from the galaxy sample. We will describe numerical implementations of these steps in detail later.

As discussed in §1, there are multiple uncertainties that are critical in estimating $s_{*}$. Our decomposition strategy described above ensures that each type of uncertainty can be physically modeled in the relevant step, and that the forward incorporation of all steps naturally propagates all uncertainties into the total uncertainty of $s_{*}$. In addition, this also allows us to quantify contributions of individual uncertainties by incrementally adding them into the chain.

In principle, the counts of halos are noisy because massive objects are rare in a small volume by definition. This sampling uncertainty is further enhanced by other sources of uncertainty in the halo-to-galaxy mapping due to the increased steepness of the Schechter function toward the high-mass end. Additionally, as the sample size goes below $\sim 100$, the discreteness in galaxy counts and its skewness must be carefully taken into account (e.g., Trenti & Stiavelli, 2008). To address these issues, we deliberately choose the cumulative cosmic stellar mass density, $s_{*}=\rho_{*}(>M_{*})$, as the summary statistic when comparing the theoretical prediction with observational data (see also Labbé et al., 2023; Lovell et al., 2022; Boylan-Kolchin, 2023, for the same usage). We describe the uncertainty of $s_{*}$ using the probability distribution, $P(s_{*})$, instead of summary statistics with compressed information, such as the average, variance, and quantile. If desired, these compressed quantities can be estimated from $P(s_{*})$.

In the following, we specify each of the transformations involved in the mapping from halos to galaxy statistic and describe uncertainties injected into them.

Volume sampling operator: The transformation $\mathcal{T}_{\rm CV}$, implemented by drawing a sub-volume from the periodic box of ELUCID and filtering out all halos from $S_{\rm halo}$ outside the selected volume, constitutes the volume sampling operator. The uncertainty in the sampling is naturally incorporated due to the fluctuation of the density field seeded by the initial condition. Throughout this paper, we refer to the uncertainty introduced by this sampling step collectively as the cosmic variance (CV). We do not subtract the Poisson shot noise from this uncertainty because its exact effect on $\rho_{*}(>M_{*})$ is not analytically traceable, and any inaccuracy in the approximation can change the extreme statistics significantly. When estimating CV for a real survey, the sub-volume must be chosen to have a consistent shape with the survey, especially for a survey with a pencil-beam-like space coverage. This is because a beam-shaped volume passes many different environments and is thus less likely to be affected by a single over- or under-density region in comparison to a cube-shaped volume, as discussed by Trenti & Stiavelli (2008) and Moster et al. (2011).

Halo mass estimator: The operator $\mathcal{T}_{\rm M_{\rm halo}}$ obtains halo masses from the sample of halos obtained from the previous step. In this paper, we define the halo mass, $M_{\rm halo}$, as the total mass of all the dark matter particles bound to the central subhalo when mapping the halo to its central galaxy. We have verified that the definition of halo mass has a negligible effect on our results for halos with $M_{\rm halo}\geqslant 10^{10.5}\,h^{-1}{\rm M_{\odot}}$ in comparison to other more significant uncertainties. In the $\Lambda$CDM paradigm, the combination of halo mass and cosmic baryon fraction, $f_{\rm b}=\Omega_{B,0}/\Omega_{M,0}$, sets a natural upper limit on the baryon mass that can be converted into stars. However, backsplash halos (e.g., Balogh et al., 2000; Wang et al., 2009; Wetzel et al., 2014; Diemer, 2021; Wang et al., 2023a), whose dark matter has left the host halo while its gas may remain in the host halo, can increase the baryon-to-dark matter ratio in the host halo, thereby lifting the upper limit. To model this uncertainty, we add the mass of backsplash halos to the $M_{\rm halo}$ of the host halo, and we use this new mass as an “effective halo mass” of the host to estimate the upper limit on the available baryons. Our tests show that the added mass is, on average, about $50\%$ of the original mass for the most massive halos at $z=7\sim 10$, and less significant for less massive halos.

Stellar mass estimator: The transformation $\mathcal{T}_{\rm M_{*}}$ assigns each halo a stellar mass, $M_{*}$, for its central galaxy, based on the halo mass obtained from the previous step. Assuming a constant star formation efficiency $\epsilon_{*}$, we model this mapping as:

The simplification using a constant star formation efficiency is motivated by observational measurements from a HST+Spitzer sample by Stefanon et al. (2021), where the star formation efficiency is found to be flat at the high-halo-mass end and shows little evolution in $z\sim 6-10$. Recent empirical models (e.g., Behroozi et al., 2019), semi-analytical models (e.g., Yung et al., 2023), and hydrodynamic simulations (e.g., Kannan et al., 2023) do not seem capable of producing such high-mass galaxies at high $z$, suggesting a high value of $\epsilon_{*}$ must be assumed to match the observational data.

However, in real observations, the measurement of $M_{*}$ suffers from systematic and random errors. Each assumption on the IMF, star formation history, dust attenuation, and photometric redshift can introduce significant amounts of uncertainty into the SED fitting (see Conroy, 2013; Behroozi et al., 2010; Stanway, 2020, for a comprehensive review of techniques and uncertainties). Indeed, Harikane et al. (2023) and Finkelstein et al. (2023) found that these assumptions can introduce significant uncertainty in the results of JWST galaxies. van Mierlo et al. (2023) performed a test on a $z\sim 7$ massive galaxy with different SED fitting codes, and found a $0.76\,{\rm dex}$ systematic offset in stellar mass among the five codes adopted. Similarly, Labbé et al. (2023) performed SED fittings with seven different methods, and found $\sim 0.1-0.3\,{\rm dex}$ random error in stellar mass by using a given method, and $\sim 0.2-1\,{\rm dex}$ systematic difference between different methods. For some galaxies, the systematic difference was found to be as large as $2\,{\rm dex}$. To incorporate such uncertainties, we introduce two free parameters in the mapping from $M_{\rm halo}$ to $M_{*}$: a constant bias factor $b_{M_{*}}$ and a normally distributed random error with zero mean and a standard deviation of $\sigma_{M_{*}}$. A typical value for these parameters is about $0.1$ to $0.5$ dex in low-redshift measurements, depending on galaxy sample in use and the model choices made (e.g., Li & White, 2009; Conroy, 2013), and they are expected to be at least as large for the high-$z$ galaxies concerned here.

Statistical operator: The final transformation, $\mathcal{T}_{\rm stat}$, is a statistical function that compresses the sample of $M_{*}$ into $s_{*}=\rho_{*}(>M_{*})$. This is simply a cumulative histogram weighted by $M_{*}$. No error should be added in this step.

By applying all four steps, we obtain an estimate for the stellar mass density that mimics real observations. To obtain a sample of $\rho_{*}$, we repeatedly select sub-samples from the entire volume of ELUCID and obtain $\rho_{*}$ from each of them. From this sample of $\rho_{*}$, we numerically compute the probability distribution $P(\rho_{*})$ and other summary statistics based on $\rho_{*}$. Finally, we make a model-observation comparison using $P(\rho_{*})$ to obtain the probability of observing a particular value of $\rho_{*}$.

In some of our analyses, we use halo and galaxy properties other than those defined above. When linking halos across redshifts, we are more concerned about properties of a halo as a whole, rather than stellar contents within it. In this case, we adopt the “tophat” halo mass, calculated using dark matter particles enclosed in a virial radius within which the mean density is equal to that of the spherical collapse model (Bryan & Norman, 1998). When demonstrating the absolute number of halos/galaxies in bins of given masses, we use the un-weighted, un-normalized histogram and its cumulative version, instead of the distribution density. We will clarify their usage in the description of our results.

The top-left panel of Fig. 1 shows a sample of halos with $M_{\rm halo}\geqslant 2\times 10^{11}\,h^{-1}{\rm M_{\odot}}$ at $z=8$ in a slice of the simulation box. At this redshift, there are already a large number of massive halos. Assuming a large but reasonable star formation efficiency, $\epsilon\sim 0.5$, these halos are capable of hosting massive galaxies with $M_{*}\gtrsim 10^{10}\,h^{-1}{\rm M_{\odot}}$ that have been identified by JWST (e.g., Rodighiero et al., 2022; Labbé et al., 2023). The two panels in the first row of Fig. 2 show the cumulative number density, $N(>M_{*})$, and the cumulative cosmic stellar mass density, $\rho_{*}(>M_{*})$, respectively, as functions of $M_{*}$ for galaxies in the entire simulation volume of ELUCID at snapshots from $z=7$ to $10$, with $M_{*}$ obtained by assuming a constant star formation efficiency $\epsilon_{*}=0.5$. More than 1000 galaxies with $M_{*}\geqslant 10^{10}\,h^{-1}{\rm M_{\odot}}$ can be found in the entire volume of ELUCID. In the following section, we examine whether or not the number of massive objects observed by JWST can be accommodated by the simulation.

Using the halo-to-galaxy transformations defined in the previous section, we now quantify the effect of their uncertainties on the estimate of the stellar mass density, $\rho_{*}(>M_{*})$. We first describe the effect of each uncertainty separately by zeroing out the other uncertainties. We then combine them to see their synergistic effect.

The second and third rows of Fig. 2 show the effect of mass enhancements introduced by backsplash halos on the transformation $\mathcal{T}_{M_{\rm halo}}$, where the “effective halo mass” is used to predict the stellar mass in it. To compare the results at the massive end, where the number/mass density drops down to the limit of the simulation, we fit the histogram with a Schechter function and extrapolate it to higher mass. Compared with the results without backsplash, the $N(>M_{*})$ or $\rho_{*}(>M_{*})$ for galaxies with $M\lesssim 10^{9}\,h^{-1}{\rm M_{\odot}}$ do not change. However, the effect of backsplash is more significant for galaxies of higher mass. The density can reach $\sim 200\%$ of its original value at $M_{*}\sim 10^{10}\,h^{-1}{\rm M_{\odot}}$ for all redshifts shown, and go beyond $1000\%$ at the high-mass end. This mass-dependency of the backsplash enhancement is the consequence of two effects. The first is that more massive halos are, on average, embedded in denser environments where high-speed close encounters are more frequent and the fraction of mass brought in by backsplash halos is higher. Second, the Schechter function for halo mass distribution declines exponentially at the high-mass end, so that the change in the halo mass is magnified in the halo mass function. It is not yet clear whether or not the baryon component of a backsplash halo can be effectively acquired by the host halo after the dark matter component is ejected, and whether or not the acquired gas can form stars effectively. Future detailed hydrodynamic simulations are needed to clarify the ambiguity.

Fig. 3 shows the effect of cosmic variance (CV) on $\rho_{*}(>M_{*})$ introduced by the sub-sampling operator $\mathcal{T}_{\rm CV}$:

Here, $p\in[0,1]$ is the target fraction of the probability mass centralized at the median, $\rho_{*,0.5(1+p)}$ is the quantile at the percentile $0.5(1+p)$, and $\rho_{*,0.5}$ is the median. We deliberately avoid using the average and standard deviation in the measurement of error, as they are not stable nor informative for a highly skewed distribution.

Throughout this paper, we use 1, 2 and 3$\sigma$ to denote cases where $p=0.68$, $0.95$ and $0.99$ are used for the quantiles, respectively. The distribution of $\rho_{*}$ is obtained numerically by sub-sampling the periodic box of ELUCID with a given volume, and the quantiles are estimated from the distribution. Here, we take the $z=8$ snapshot as an example and show CV for different volumes and for different thresholds of $M_{*}$. Solid, dashed, and dotted colored lines in Fig. 3 show 1, 2, and 3$\sigma$ CVs for a cubic volume of a given size. For low-mass galaxies with $M_{*}\sim 10^{9}\,h^{-1}{\rm M_{\odot}}$, the $1\sigma$ CV is estimated to be about $10\%$ for a volume of $(100\,h^{-1}{\rm{Mpc}})^{3}$ and increases monotonically to $\sim 70\%$ for a volume of $(25\,h^{-1}{\rm{Mpc}})^{3}$. This is consistent with the results obtained from the cosmic variance estimator given by Chen et al. (2019) for low-$z$ samples. Black lines show the CV expected for JWST CEERS at $z=7-8.5$, assuming a $38\,{\rm arcmin}^{2}$ effective sky area as used by (Labbé et al., 2023). The sub-sampling in ELUCID is conducted using beam-shaped volumes with a size of $(36\,h^{-1}{\rm{Mpc}})^{3}$ and a tangential-to-normal aspect ratio of $12/342$. The 1, 2, and 3$\sigma$ CVs are estimated to be $25\%$, $60\%$, and $80\%$ for low-mass galaxies. The moderate effect of CV on this mass scale indicates that the density of $M_{*}\lesssim 10^{9}\,h^{-1}{\rm M_{\odot}}$ galaxies can be estimated reliably in JWST CEERS if other uncertainties are well controlled.

However, the effect of CV has a significant dependence on stellar mass, with uncertainty reaching $100\%$ at higher stellar masses and eventually becoming divergent at the highest-mass end, regardless of the survey volume. The point of divergence shifts rightward as the volume increases, as it emerges when the median of $\rho_{*}(>M_{*})$ approaches zero. In Fig. 4, the vertical lines represent the masses at which $\rho_{*}(>M_{*})$ is measured by Labbé et al. (2023) using the four massive galaxies from JWST CEERS. Three of these lines intersect with the black solid line at CV $\sim 100\%$, suggesting that the CV is a significant effect for this small sample. The rightmost line, obtained from a galaxy with the most extreme stellar mass of $10^{10.89}{\rm M_{\odot}}$ at $z\sim 7.48$, is located in the divergent region of CV, indicating an incredible effect of uncertainty in the statistics drawn from this single galaxy. Since the divergence of CV is tightly related to the vanishing of the median $\rho_{*}(>M_{*})$ by its definition, it is more informative to model $\rho_{*}(>M_{*})$ forwardly, by directly predicting the distribution $P[\rho_{*}(>M_{*})]$ and the probability of observing a value of $\rho_{*}(>M_{*})$. We will demonstrate this later in this section.

Figure  4 illustrates the impact of the uncertainty in the stellar mass estimator, $\mathcal{T}_{M_{*}}$. To demonstrate the effect, we take all halos at the $z=8$ snapshot of ELUCID as an example. In the left panel, we show the change of $\rho_{*}(M_{*})$ when different levels of systematic error i.e. bias, are added to the stellar mass predicted by Eq. 2. A negative (positive) bias naturally decreases (increases) $\rho_{*}(M_{*})$, as it is equivalent to a lower-left (upper-right) shift of the distribution function. However, when considering the ratio of the biased distribution to the unbiased one, a stellar-mass dependency is observed. This is similar to the mass dependency of CV, where a small perturbation to the transformation is magnified significantly owing to the steepness of the Schechter function at the high-stellar-mass end. The overall outcome is an increased effect of the uncertainty in the stellar mass estimate at the high-mass end, and a divergence occurs when the unbiased galaxy number $N(>M_{*})$ approaches zero.

In the right panel of Fig. 4, we show the effect of random error, known as the Eddington bias (Eddington, 1913), in the stellar mass estimate. Here, we model the error as a Gaussian random variate with zero mean and varying standard deviation, $\sigma_{M_{*}}$. Lines with different colors represent results with different $\sigma_{M_{*}}$, and are compared to the result assuming no error. Once again, to enable a comparison in the full mass range, we fit the histogram to a Schechter function and extrapolate it to the uncovered mass range. As the distribution of $M_{*}$ has a negative derivative, a symmetric noise in $M_{*}$ produces an asymmetric effect on the number counts of galaxies. As the derivative becomes more negative toward the high-stellar-mass end, the effect of this type of noise becomes more significant. We do see this expected behavior in the panel, where the histogram is lifted everywhere by a nonzero error of $M_{*}$. With a moderate random error of 0.2 dex (0.3 dex), the increase of $\rho_{*}(M_{*})$ is $\sim 100\%$ at $M_{*}\sim 10^{9.2}\,h^{-1}{\rm M_{\odot}}$ ($10^{9.7}\,h^{-1}{\rm M_{\odot}}$), and becomes divergent when the number of galaxies $N_{>M_{*}}$ approaches zero.

As suggested by Lovell et al. (2022) and Boylan-Kolchin (2023), an assumption of $\epsilon_{*}\sim 1$ still results in a $\sim 3\sigma$ tension in $\rho_{*}(>M_{*})$ with recent JWST observations, especially when the massive galaxies at $z=7-10$ from the JWST CEERS sample (Labbé et al., 2023) are included. With all the uncertainties introduced above, and the fact that all of them have a strong effect on the statistics drawn from samples of small size, it is likely that the observational results can be reproduced by the theory with much less tension, as demonstrated in the following.

Fig. 5 shows the cumulative cosmic stellar mass density $\rho_{*}(M_{*})$ as a function of the limiting stellar mass $M_{*}$ for two different redshifts, $z\sim 8$ and $z\sim 9$, where six extremely massive galaxy candidates are identified and $\rho_{*}(>M_{*})$ estimated by Labbé et al. (2023). We demonstrate effects of different sources of uncertainties by incrementally adding them into the transformations, $\mathcal{T}_{\rm CV}$, $\mathcal{T}_{\rm M_{halo}}$, and $\mathcal{T}_{\rm M_{*}}$, respectively, in the mapping from halo sample, $S_{\rm halo}$, to the galaxy statistic, $\rho_{*}(>M_{*})$, as introduced in §2.2. The JWST detections are overplotted for comparison.

Black lines with gray shading in the left column of Fig. 5 show the median and different $\sigma$ ranges derived from the probability distribution $P[\rho_{*}(>M_{*})|M_{*}\,,\epsilon_{*}=0.5\,,{\rm CV}]$ when cosmic variance is incorporated by volume sub-sampling and a star formation efficiency $\epsilon_{*}=0.5$ is adopted to convert halo mass to stellar mass. The results shown here are consistent with those obtained by Lovell et al. (2022) and Boylan-Kolchin (2023) that some data points are outside the $\sim 3\sigma$ range. The exact tension level is slightly different, likely because of their simplification in the error calculation with analytical approximations, the difference in the cosmologies adopted, and the difference in the version of galaxy data. The point obtained by us from a $z\sim 8$ extreme is the only one that lies outside the $3\sigma$ range, while the $z\sim 9$ extreme marginally touches the $2\sigma$ boundary. Other points are all contained within the $2\sigma$ boundary.

As demonstrated in §3.1 and §4, Eddington bias is able to lift the galaxy number by more than an order of magnitude, which may help explain the outliers we see here. The middle column of Fig. 5 shows the probability distribution $P[\rho_{*}(>M_{*})|M_{*},\,\epsilon_{*}=0.5,\,{\rm CV},\,{\rm\sigma_{M_{*}}}=0.3\,{\rm dex}]$ when a Gaussian random error with a conservative 0.3 dex standard deviation is added to the estimate of $\log\,M_{*}$. A significant shift of the median $\rho_{*}(>M_{*})$ and a significant broadening of the quantile ranges can be seen after the incorporation of this uncertainty. The $z=8$ extreme is now safely contained within the 2$\sigma$ range of $\rho_{*}(>M_{*})$, and all the $z=9$ points are contained by the $1\sigma$ range. Thus, once the cosmic variance and the Eddington bias are included, there is no significant tension between the $\Lambda$CDM paradigm and the JWST observations, if the star formation efficiency can reach $0.5$ at these redshifts.

The right column of Fig. 5 shows the results for $\rho_{*}(M_{*})$ when the third source of uncertainty, the backsplashed mass, is taken into account and added into the halo mass estimator, $\mathcal{T}_{M_{\rm halo}}$. Unlike cosmic variance and random error in the stellar mass estimate, this effect always increases halo mass and thus provides a larger effective baryon mass for star formation. The outcome is obvious: all the data points at $z\sim 8$ and $9$ are now below the upper $2\sigma$ quantile once a star formation efficiency $\epsilon_{*}=0.5$ is assumed. At $z\sim 8$, three of the four data points actually lie close to the lower $2\sigma$ quantile, while the most massive is close to the upper $1\sigma$ quantile.

The colored curves in Fig. 5 show the results obtained by considering lower star formation efficiencies, $\epsilon_{*}=0.2$, $0.1$, and $0.02$, which are more realistic expected from low-z observations. Each curve is overplotted with a series of horizontal error bars, which indicate the $1\sigma$, $2\sigma$, and $3\sigma$ ranges at the stellar mass density corresponding to the leftmost observational data point in the same panel. When solely considering cosmic variance and assuming $\epsilon_{*}=0.2$, none of the observational data points fall within the $3\sigma$ ranges. However, when introducing a standard deviation of 0.3 dex to the estimate of $\log,M_{*}$, the leftmost data points for both the $z=8$ and $z=9$ observations fall within the $3\sigma$ range. Furthermore, when incorporating the backsplash effect, these observations shift into the $2\sigma$ ranges. These findings suggest that a star formation efficiency of $\epsilon_{*}=0.2$ can marginally account for the high stellar mass densities obtained from the CEERS sample. Nevertheless, when considering more common values for the low-z Universe, for example, $\epsilon_{*}=0.1$ and $\epsilon_{*}=0.02$ (e.g., Yang et al., 2003; Moster et al., 2010; Yang et al., 2012; Reddick et al., 2013; Behroozi et al., 2013; Moster et al., 2013; Birrer et al., 2014; Lu et al., 2014; Rodríguez-Puebla et al., 2017; Shankar et al., 2017; Moster et al., 2018; Behroozi et al., 2019), the modeled stellar mass densities do not reach the observational values. This indicates that a simple extrapolation of low-z observations is insufficient in explaining the CEERS observations.

All the results presented above highlight the importance of fully modeling all the uncertainties when interpreting the observational data. Since each of the three sources of uncertainty considered above has a positive effect on $\rho_{*}(>M_{*})$, a lower star formation efficiency $\epsilon_{*}$ is needed to explain the observed data when the uncertainty involved is larger. Conversely, a greater $\epsilon_{*}$ is necessary to explain the observational data if the uncertainty is smaller. Since the exact level of uncertainty in each transformation step is not known a priori, we provide a general prediction for the relationship between the level of uncertainty and the probability to accommodate the observational data. For a given observed density, $\rho_{*,{\rm obs}}(>M_{*})$, we define the probability to observe it as the remaining probability mass of having $\rho_{*}(>M_{*})>\rho_{*,{\rm obs}}(>M_{*})$:

where ${\rm CDF}[\rho_{*}(>M_{*})]$ is the cumulative distribution function of $\rho_{*}(>M_{*})$ derived from $P[\rho_{*}(>M_{*})]$. This probability can be used to obtain the level of uncertainties and the star formation efficiency required to explain the observation.

Fig. 6 displays the predicted $p_{\rm obs}$ for the six data points obtained by Labbé et al. (2023) as a function of star formation efficiency $\epsilon_{*}$, taking into account different sources of uncertainty. The upper left panel shows the case that does not include the random error in the stellar mass estimator. Without mass added by backsplash halos, four of the samples have $p_{\rm obs}>10\%$ for $\epsilon_{*}=0.5$. However, one extreme at $z\sim 8$ has $p_{\rm obs}<10\%$ even for $\epsilon_{*}=1$, suggesting a possible tension. With the backsplash effect included, the $p_{\rm obs}$ value is three times the original value at $\epsilon_{*}=1$ for this extreme, and it is non-zero at $\epsilon_{*}=0.5$. A random error with $\sigma_{\rm M_{*}}$ ranging from $0.15\,{\rm dex}$ to $0.3\,{\rm dex}$ on the stellar mass estimate increases the probability further. A $\sigma_{M_{*}}$ of 0.5 dex drastically changes the trend of $p_{\rm obs}$; even $\epsilon_{*}=20\%$ is sufficient to move all the observed samples to within the $1\sigma$ range ($p_{\rm obs}\geqslant 0.16$).

It is important to note that the probability prediction presented here is based solely on empirical modeling of uncertainties in the $\Lambda$CDM cosmology. Therefore, the results are quite general and can be extended to other surveys, regardless of the details of the survey strategy, data processing and statistical model. However, further work is needed to confirm conjectures presented here. Surveys with large areas are useful in suppressing cosmic variance and thus in reducing the uncertainty in volume sampling. Hydrodynamic simulations and semi-analytical modeling are needed to verify or rule out the proposed effect of backsplash, as well as to understand the conversion of baryons to stars. Deeper imaging, high S/N spectroscopy, and precise stellar population synthesis are critical in order to narrow down the posterior parameter space in the estimates of redshift and stellar mass estimate. All these, together, may eventually resolve the tension or provide definitive evidence for the need of new cosmology and new physics.

In order to investigate the properties of the descendants of high-redshift massive galaxies, we use subhalo merger trees in the ELUCID simulation. This enables us to establish a connection between the high-redshift galaxies, as predicted by our empirical transformations in §2.2, and low-redshift halos which host the descendants of these galaxies. For convenience these low-$z$ halos are referred to as “ancient descendant halos”, or “ADH” for short. Note that ADH are defined at a given redshift, $z$, for galaxies modeled at a higher redshift, $z_{\rm g}$. We study both the spatial and mass distribution of these ADH. It is important to note that the density field of ELUCID is constrained by real observations, as outlined in §2.1. Thus, the connection established can also be used to understand the formation history of real massive clusters in the local universe back to a time when the universe was only approximately 500 million years old.

The green histograms in Fig. 7 depict the distribution of mass, $M_{\rm halo,host}$, for the ADH of high-$z_{\rm g}$ massive galaxies. In this section, we use the top-hat mass for halos, since it is easier to infer from observation and since we are conducting statistics for a halo as a whole (e.g., Yang et al., 2007; Yang et al., 2012; Lim et al., 2017; Tinker, 2020; Yan et al., 2020; Wang et al., 2020; Hung et al., 2021). We select high-$z_{\rm g}$ galaxies as the central galaxies in the most massive halos with $M_{\rm halo,host}\geqslant 10^{11}\,h^{-1}{\rm M_{\odot}}$ at $z_{\rm g}=8$ in the entire ELUCID volume. From $z=8$ to $z=0$, the descendant halos of these massive high-redshift galaxies are found preferentially in the massive end of the mass distribution of the total population (the black histogram), although the distribution can extend to masses as low as $\sim 10^{12}\,h^{-1}{\rm M_{\odot}}$. At the massive end of the histogram, almost all halos are ADH of massive galaxies at $z_{\rm g}=8$. The minimum mass of the descendant halo is approximately $10^{11.5}\,h^{-1}{\rm M_{\odot}}$ ($10^{12}\,h^{-1}{\rm M_{\odot}}$) at $z=2$ ($z=0$), while the median is around $10^{12.8}\,h^{-1}{\rm M_{\odot}}$ ($13.5\,h^{-1}{\rm M_{\odot}}$). However, among the total population of halos at $z=0$ with $M_{\rm halo}\sim 10^{13.5}\,h^{-1}{\rm M_{\odot}}$, only about $10\%$ are ADH of massive galaxies at $z_{\rm g}=8$, and the fraction drops exponentially towards lower halo masses.

To quantify the precise fraction of high-$z_{\rm g}$ massive galaxies that end up in low-$z$ halos of different masses, we present the ratio between the ADH mass histogram with the unconditional halo-mass histogram at different redshift, $z$, in Fig. 8 using the same green color as in Fig. 7. The shaded area attached to each curve represents the standard deviation estimated from 100 bootstrapped samples. A halo with a mass of $10^{14}\,h^{-1}{\rm M_{\odot}}$ ($10^{15}\,h^{-1}{\rm M_{\odot}}$) at $z=2$ ($z=0$) almost certainly contains at least one descendant of the most massive galaxies at $z_{\rm g}=8$, while a halo with a mass of $10^{13.5}\,h^{-1}{\rm M_{\odot}}$ ($10^{14.5}\,h^{-1}{\rm M_{\odot}}$) has a probability of $\gtrsim 70\%$ to do so. For halos with a mass of $10^{13}\,h^{-1}{\rm M_{\odot}}$ ($10^{13.5}\,h^{-1}{\rm M_{\odot}}$), the probability drops to $\lesssim 30\%$. Note that we do not differentiate between centrals and satellites in descendant halos.

The brown histograms labeled as “ADH-C” in Fig. 7 represent the number of ADH that contain the descendant galaxies as their central galaxies. The lines in Fig. 8 show the fraction of these halos among the total population at the same redshift. The similar values of the green (ADH) and brown (ADH-C) histograms/lines suggest that most of the high-$z_{\rm g}$ massive galaxies eventually become the dominant galaxies or parts of the dominant galaxies of massive halos at low-$z$. Due to hierarchical formation of galaxies in the $\Lambda$CDM model, some of the high-$z_{\rm g}$ massive galaxies become parts of more dominating centrals at low-$z$ while others remain as the dominating centrals. To demonstrate this, the purple histograms and lines in Fig. 7 and Fig. 8, marked as “ADH-C(1st)”, represent the number and fraction, respectively, for low-$z$ halos whose central subhalos are the first descendants of the $z_{\rm g}=8$ massive galaxies. Here, a subhalo A is considered to be the first progenitor of another subhalo B in the subsequent snapshot, if A has the greatest bound mass among all of B’s progenitors. In this case, we refer to B as the first descendant of A, only to make clear the special status of A. In addition, the chain of first progenitors (i.e., the main branch) of B can be determined by recursively identifying the first progenitor along the links in the tree towards higher redshift. Thus, ADH-C(1st) at a given $z$ are halos at $z$ that each has a massive progenitor galaxy at $z_{\rm g}$ in the main branch of its central subhalo. The results indicate that approximately $30\%-40\%$ of the central galaxies in the most massive halos at low-$z$ halos are direct descendants of the most massive, high-$z_{\rm g}$ galaxies.

The findings presented here provide a compelling guidance regarding the sites to identify remnants of the old stellar population formed at high-$z$. For instance, since central galaxies in the majority of massive clusters with $M_{\rm halo}\sim 10^{14}\,h^{-1}{\rm M_{\odot}}$ ($10^{15}\,h^{-1}{\rm M_{\odot}}$) at $z=2$ ($z=0$) are highly likely to contain ancient stars from $z\sim 8$, targeting the central galaxies of these potential ADH with infrared spectroscopic observations might reveal the presence of this ancient stellar population. This, in turn, would facilitate comparisons with stellar populations that formed later to investigate star formation in high-$z$ galaxies.

In Fig. 9, we partition the $z_{\rm g}=8$ galaxy sample, whose descendants are centrals at some lower redshift (indicated by orange symbols in Fig. 7 and Fig. 8), into subsets based on $N_{\rm merge}$, the number of mergers they underwent before entering the central galaxies of their ADH. Note that here we only include merger events in which the subhalo under consideration was the less-massive one in the encountering subhalos. The fractions of galaxies in these subsets as a function of host halo mass are depicted by lines and shadings of different colors for different redshifts below $z=8$. The fraction of galaxies with $N_{\rm merge}=0$ consistently declines with increasing host halo mass at all redshifts. This is because a host halo with greater mass harbors more satellites, thereby increasing the likelihood of close encounters among its satellites. The proportion of galaxies with $N_{\rm merge}=1$ or $2$ steadily increases with host halo mass at all redshifts. At the high-mass extreme, galaxies with $N_{\rm merge}=1$ emerge as the dominant population, indicating that most central galaxies in the most massive halos at lower redshifts have undergone at least one merger event in their history. As redshift decreases, the $N_{\rm merge}=2$ cases become increasingly prevalent, eventually outnumbering the $N_{\rm merge}=0$ case.

Given that the galaxies we have chosen at $z_{\rm g}=8$ are the most massive ones at that time, it is highly probable that the merger events in which they entered the central galaxies of their ADH were all major ones. These major mergers can cause angular momentum loss, black hole growth, morphology transformation, and quenching of galaxies. Our findings thus imply that the central galaxies in local massive halos underwent one or two bursts in star formation. These outcomes may serve as valuable priors for Bayesian spectral synthesis models (e.g., Zhou et al., 2019; Johnson et al., 2021), where one or more burst components can be incorporated into the star formation history.

The significant spread of the $M_{\rm halo,host}$ distribution for ADH at $z\lesssim 5$ shown in Fig. 7, coupled with the diverse merger histories of the descendants of high-$z_{\rm g}$ massive galaxies shown in Fig. 9, implies a disruption of the rank order for halo/stellar mass during the evolution. Consequently, abundance matching between progenitors and descendants using only halo and stellar mass may not be accurate, thus presenting a problem for rank-based empirical methods to link galaxies between different redshifts (e.g., Zheng et al., 2007; Behroozi et al., 2019; Wang et al., 2023b; Zhang et al., 2022a). Our findings suggest that these methods need to be improved, and that the inclusion of secondary halo properties might be needed. We also note that some efforts have been made to identify proto-clusters at intermediate redshifts ($1\lesssim z\lesssim 3$) and use them to establish a statistical link between galaxies at different redshifts (e.g., Chiang et al., 2013, 2014; Cai et al., 2016, 2017; Wang et al., 2021). With ongoing and future surveys, these methods can be extended and applied to high-$z$ data. This will open a new avenue to study the evolution of the galaxies across the history of the universe.

The diversity of the descendant distribution highlights the need to investigate the evolution of individual galaxy clusters in the low-$z$ universe. In this section, we present the assembly histories of selected halos that can be matched to real clusters of galaxies, capitalizing on ELUCID with its nature of being constrained by the galaxy distribution in the SDSS survey.

Coma is a massive cluster located in the nearby universe, positioned at $({\rm RA}$, ${\rm Dec}$, ${\rm z_{\rm CMB}})$ $=(194.95^{\circ}$, $27.98^{\circ}$, $0.0231)$. According to Yang et al. (2007); Yang et al. (2012), its estimated halo mass is about $9.3\times 10^{14}\,h^{-1}{\rm M_{\odot}}$. In ELUCID, the matched halo is located at $(194.81^{\circ}$, $27.91^{\circ}$, $0.0240)$ and has a halo mass of $9.4\times 10^{14}\,h^{-1}{\rm M_{\odot}}$. In Fig. 1, we follow a halo with $M_{\rm halo}$ slightly larger than $10^{11}\,h^{-1}{\rm M_{\odot}}$ at $z=8$ (top-left panel), which joins the Coma cluster at $z=0$ (lower-right panel). Note that this halo is actually not the most massive one at $z=5$, $3$ and $1.5$ among all the halos displayed in the same panels. Only at lower redshifts, $z=1.5$ and $0$, does Coma become the most massive cluster. This suggests an unusual characteristic of Coma, which formed relatively late, perhaps via violent mergers during later epochs.

To understand the evolution of Coma in more detail, Fig. 10 displays the mass distribution reconstructed by ELUCID for its progenitor halos, with a red arrow indicating the host halo mass, at each snapshot, of the first progenitor subhalo of the central subhalo of Coma at $z=0$. For comparison, we depict the median mass distribution, as well as the 1 and $2\sigma$ ranges, for progenitors of halos at $z=0$ with $10^{14.75}\leq M_{\rm halo}/(\,h^{-1}{\rm M_{\odot}})\leq 10^{15.25}$. At $z=8$, The most massive progenitor of Coma has $M_{\rm halo}\sim 2\times 10^{11}\,h^{-1}{\rm M_{\odot}}$, barely touches the high-mass end, as shown in the top-left panel of Fig. 7. At $z=0$, Coma has $M_{\rm halo}\approx 10^{15}\,h^{-1}{\rm M_{\odot}}$, which is at the high-mass end, as demonstrated in the lower-right panel of Fig. 7. At $z=5$, progenitors of Coma are less massive than the average for similarly massive halos at $z=0$. At $z=3$, the first progenitor of Coma forms, and at $z=2$, it merges with another halo and becomes a satellite. This suggests that the birthplace of Coma was in a crowded environment. At $z\leq 1.5$, the massive end of the distribution of Coma’s progenitor mass exceeds the average, and the overall amplitude of the distribution also goes above the average (as seen from the red and black solid histograms). At $z\leq 1$, several massive progenitors form and eventually merge with the first progenitor. All these suggest that violent mergers at $z\lesssim 3$, particularly at $z\leq 1$, played a critical role for Coma to assemble a large amount of mass by $z=0$.

In Fig. 11, we present the projected 2-D spatial positions of Coma’s progenitors. A star is used to mark the spatial location of the first progenitor subhalo, at each snapshot, of the central subhalo of Coma at $z=0$, shown in red if the progenitor subhalo is a central and blue if it is a satellite. The abundance of progenitors at $z\lesssim 5$ indicates the dense environment of Coma, where frequent mergers are taking place. This is in agreement with observations, such as those from SDSS (Abazajian et al., 2004), which report a substantial fraction of red galaxies in Coma (e.g., de Propris et al., 1998; Eisenhardt et al., 2007; Jenkins et al., 2007; Adami et al., 2009; Miller et al., 2009; Mahajan et al., 2010; Hammer et al., 2010; Mahajan et al., 2011; De Propris, 2017). At $z=0.5$, two massive structures are visible among the progenitors, consistent with the observed presence of two massive elliptical galaxies near the center of Coma. It is noteworthy that overall mass distribution of Coma is not spherically symmetric and that it is surrounded by filamentary structures that contribute to its continuous mass assembly. Future X-ray observations of these filaments may provide additional constraints on the environment and assembly history of the Coma cluster.

According to the extended Press-Schechter formalism (Lacey & Cole, 1993), the progenitor mass distribution is expected to depend on the mass of descendant halos. Fig. 12 depicts the progenitor mass distributions at various redshift for present-day halos in the mass range $10^{14.2}\leq M_{\rm halo}/(\,h^{-1}{\rm M_{\odot}})\leq 10^{14.4}$. The median, $1\sigma$, and $2\sigma$ ranges are displayed by solid, dashed and dotted lines, respectively. In addition, we also show the progenitor distributions for two $z=0$ halos matched with Abell 1630 and Abell 1564.

Abell 1630 is a galaxy cluster located at $({\rm RA}$, ${\rm Dec}$, ${\rm z_{\rm CMB}})$ $=(192.93^{\circ}$, $4.56^{\circ}$, $0.064)$, and its estimated halo mass is $1.78\times 10^{14}\,h^{-1}{\rm M_{\odot}}$ by Yang et al. (2007); Yang et al. (2012). In the ELUCID simulation, its position is simulated at $(193.18^{\circ},4.42^{\circ},0.064)$, and its halo mass is $2.01\times 10^{14}\,h^{-1}{\rm M_{\odot}}$. Abell 1630 is quite unique in that its first progenitor significantly dominates the progenitor population. The mass gap between the first and other progenitors is evident from $z=8$ and continues throughout its history until $z=0$. This suggests that Abell 1630 has assembled most of its mass through continuous accretion or relatively minor mergers, which is very different from Coma, where late-time major mergers dominate the mass assembly. This cluster has a massive progenitor at $z=8$, which lies well outside the $2\sigma$ range of the mass distribution and is even more massive than the most massive progenitor of Coma at the same redshift. These suggest that Abell 1630 should have a dominating member at $z=0$ near the center. Indeed, a luminous red galaxy with an absolute magnitude of $M_{\rm r}^{0.1}=-22.17$ and a color index of $(g-r)^{0.1}=0.93$ has been identified in this cluster, consistent with the expectation of ELUCID.

The spatial locations of the progenitors of Abell 1630 at different snapshots are displayed in Fig. 13. A comparison with the corresponding plot for Coma in Fig. 11 reveals that Abell 1630 formed in an environment with only a moderate matter density, which may be the cause of its assembly history distinctive from that of Coma.

Abell 1564, the final example, is observed at $(188.74^{\circ}$, $1.84^{\circ}$, $0.0780)$. Its estimated halo mass is $2.87\times 10^{14}\,h^{-1}{\rm M_{\odot}}$, and it is reconstructed by ELUCID at $(188.96^{\circ}$, $2.02^{\circ}$, $0.0784)$ with a halo mass of $1.99\times 10^{14}\,h^{-1}{\rm M_{\odot}}$. The progenitor mass and spatial distributions are shown in Fig. 12 (the red histograms) and in Fig. 14, respectively.

The assembly history of Abell 1564 falls between the two extremes shown above in that it has gone through only one major merger in its entire history. The two merging progenitors at $z=0.5$ have nearly equal mass and, as a result, the first progenitor of Abell 1564 became a satellite subhalo during a short period of time at $z\sim 0.2$ after the merger. This late-time major merger implies that the cluster at $z=0$ may not yet be fully virialized. The non-spherical light distribution of Abell 1564 in real observations provides supports to our conclusion.

The three examples presented here illustrate three distinct assembly modes of massive clusters: one with multiple major mergers like Coma, one with only one major merger like Abell 1564, and one with continuous accretion (or minor mergers) like Abell 1630. To illustrate the differences between the modes, we present the mass assembly histories of the three clusters as a function of redshift in Fig. 15 for $M_{\rm halo,main\ branch}$, which is the mass of the main branch progenitor, $M_{\rm halo,total}$, which is the total mass contained in all progenitor halos with $M_{\rm halo}\geqslant 10^{10}\,h^{-1}{\rm M_{\odot}}$, and the ratio between them. Among the three clusters, Coma is the most massive at $z=0$, but its main progenitor is the least massive at $z\gtrsim 3$. Three faster growth stages at $z\sim 5$, $2$ and $0.5$ respectively, are clearly seen in the growth of Coma’s main branch, indicating violent merger events in its history. From Fig. 11, it is evident that these three fast accretion stages result from the richness of Coma’s progenitor halos and their asymmetric spatial distribution. The other extreme, Abell 1630, shows continuous growth of its main branch. It was massive enough to host a bright galaxy at $z\gtrsim 8$, but its main-branch mass assembly history is smooth and slow, culminating in a final halo less massive than Coma. The in-between case, Abell 1564, shows only one main-branch mass jump at $z\sim 0.5$, which results from the major merger event seen in Fig. 14. During most of its history, Abell 1564 has a main-branch-to-total ratio lying between the two extremes, Coma and Abell 1630. Note that in some cases during mergers, the curves of total mass even decrease. This is an artificial effect, due to the fact that some particles linked by the FoF algorithm are not enclosed in the "tophat" filter that defines the halo mass. The diversity of the formation mode, as represented by the main branch growth rate and the number of merger events, may be responsible for the significant variation observed in the descendant mass distribution in §4.1. Such variation must be carefully taken into account when constructing models to link galaxies over a wide range of redshifts.