Quantifying the tension between cosmological models and JWST red candidate massive galaxies

In the standard $\Lambda$ cold dark matter ($\Lambda$CDM) picture, baryon decouples from radiation after recombination and tracks the gravitational potential of dark matter. Before the formation of early galaxies, baryon and dark matter are well mixed. The fraction of baryonic mass in a massive dark matter halo therefore approximately equals to the cosmological mean $f_{b}\equiv\Omega_{b}/\Omega_{m}$, where $\Omega_{b}$ and $\Omega_{m}$ are the baryon and matter abundance parameters, respectively. If the fraction of stellar mass in baryonic matter, namely the star formation efficiency (SFE) $\epsilon$ is known, we can connect the stellar mass $M_{\star}$ to the halo mass $M_{\rm halo}$ via $M_{\star}=\epsilon f_{b}M_{\rm halo}$. This provides a way to test cosmology by estimating stellar masses of massive galaxies in a given volume. In particular, at high redshift where massive halos are predicted to be very rare, detection of very massive galaxies may pose a stringent constraint on cosmology. The red candidate massive galaxies released by the James Webb Space Telescope (JWST), for instance, have inspired a debate whether a beyond-$\Lambda$CDM cosmology is on demand (Labbé et al. 2023; Boylan-Kolchin 2023; Lovell et al. 2023; Haslbauer et al. 2022; Chen et al. 2023; Ferrara et al. 2023; Qin et al. 2023; Wang et al. 2023c). A plenty of beyond-$\Lambda$CDM models including alternative dark energy (Wang et al. 2023b; Adil et al. 2023; Menci et al. 2022), non-standard dark matter (Bird et al. 2023; Lin et al. 2023; Maio & Viel 2023; Dayal & Giri 2023), cosmic string (Jiao et al. 2023), primordial black holes (Yuan et al. 2023; Hütsi et al. 2023; Gouttenoire et al. 2023; Guo et al. 2023; Huang et al. 2023), and many others (Melia 2023; John & Babu Joseph 2023; Lovyagin et al. 2022; Lei et al. 2023) have been confronted with the JWST high-redshift candidate galaxies. Modification of primordial conditions may also lead to an overabundance of dark matter halos and hence more massive galaxies at high redshift (Huang 2019; Parashari & Laha 2023; Padmanabhan & Loeb 2023; Keller et al. 2023; Forconi et al. 2023; Su et al. 2023; Wang & Liu 2022; Tkachev et al. 2023).

The major obstacle to an accurate comparison between a theory and observed galaxies is the uncertainty in the star formation efficiency $\epsilon$ at high redshift. Galaxy formation models and low-redshift observation typically suggest $\epsilon\lesssim 0.2$ (Guo et al. 2020; Boylan-Kolchin 2023), though $\epsilon$ at high redshift may be significant different (Zhang et al. 2022; Qin et al. 2023). Because the UV radiation from star formation can ionize the neutral hydrogen at the reionization epoch, a larger star formation efficiency tends to accelerate the reionization process and finish it at an earlier time. Observations of the Lyman-$\alpha$ emitter luminosity function suggest that reionization is still on-going at $z\gtrsim 7$ (Wold et al. 2022; Goto et al. 2021), consistent with the $\epsilon\lesssim 0.2$ picture (Minoda et al. 2023; Mobina Hosseini et al. 2023). It is however difficult to derive an upper-bound of $\epsilon$ solely from the reionization history, due to the degeneracy between $\epsilon$ and other astrophysical parameters.

The common practice is then to study the cosmological implication for an assumed star formation efficiency. Although most of the studies qualitatively agree, quantitative discrepancies still exist. For example, while Boylan-Kolchin (2023) and Chen et al. (2023) find that $\epsilon\sim 0.5$ suffices to explain the Labbé et al. (2023) data in the standard $\Lambda$CDM cosmology, Menci et al. (2022) that uses the same data claims that $\Lambda$CDM is ruled out at $2\sigma$ level even with $\epsilon=1$. This is not only due to different choices of mass cut and volume cut, but also because it is hard to give an accurate description of the systematic errors of redshift and stellar mass obtained with spectral energy distribution (SED) fitting. By using seven individual redshift and stellar mass measurements with different SED-fitting codes/templates, Labbé et al. (2023) estimated extremes of stellar mass and redshift of each candidate galaxy. These extremes however are insufficient for quantitative analysis, which requires precise knowledge of the probability density functions (PDFs) of the systematic errors.

The summary statistics of the statistical errors of stellar mass are often represented as $M^{+\Delta M_{\rm sup}}_{-\Delta M_{\rm inf}}$, where $M$ is the median mass. The upper and lower errors ($\Delta M_{\rm sup}$ and $\Delta M_{\rm inf}$) are defined with $16$th and $84$th percentiles. The errors in stellar mass, and similarly in redshift if only photometry information is available, are typically asymmetric and difficult to model analytically. Accurate observation-theory comparison usually requires an end-to-end comparison between high-resolution simulations and observations. On the other hand, while the summary statistics (and the extremes from different SED-fitting methods, if available) contain limited information, they are important science products of many galaxy surveys. The aim of the present work is to develop a tool that performs a preliminary scan of the theory space based on these simple products, before planning subsequent spectroscopic observation and running costly simulations.

This paper is organized as follows. Section 2 describes the physical model of the abundance of high-redshift massive galaxies and the workflow of observation-theory comparison. Section 3 uses the tool to analyze the cosmological implication of the candidate galaxies in Labbé et al. (2023). Section 4 concludes.

Figure 1 summarizes our algorithm to do an observation-theory comparison. We consider galaxies beyond a stellar-mass threshold ($M_{\star}>M_{\rm\star,cut}$) and in a fixed comoving volume defined by the sky coverage $f_{\rm sky}$ and redshift range ($z_{\min}<z<z_{\max}$). Because both the stellar mass and redshift from SED-fitting have significant uncertainties, whether an observed galaxy is above the stellar-mass threshold and in the volume is probabilistic. The total observed galaxies above the stellar-mass threshold and in the volume, denoted as $N_{\rm obs}$, is therefore a random variable. Due to the cosmic variance and the Poisson shot noise, the theoretical prediction of the total number of galaxies above the stellar-mass threshold and in the volume, denoted as $N_{\rm th}$, is also probabilistic.

A model is ruled out when and only when $N_{\rm th}<N_{\rm obs}$. The tension between a cosmological model and the observed galaxies is therefore quantified by the probability $\mathrm{P}\left(N_{\rm th}<N_{\rm obs}\right)$, or more conveniently, the smallness of $\mathrm{P}\left(N_{\rm th}\geq N_{\rm obs}\right)=1-\mathrm{P}\left(N_{\rm th}<N_{\rm obs}\right)$. For instance, $\mathrm{P}\left(N_{\rm th}\geq N_{\rm obs}\right)=0.01$ indicates that the model under consideration is ruled at $99\%$ confidence level. Here and in what follows, $\mathrm{P}\left(\cdot\right)$ stands for probability (for discrete variable or an event) or probability density function (for continuous variable).

For a given cosmology, $\mathrm{P}\left(N_{\rm th}\geq N_{\rm obs}\right)$ is evaluated in three steps.

• 1.

  For $N_{\rm obs}=0,1,2,\ldots$ compute the distribution function $\mathrm{P}\left(N_{\rm obs}\right)$ from summary statistics of observed galaxies and user-specified assumptions about the functional form of probability density functions of stellar mass and redshift.

• 2.

  For $N_{\rm th}=0,1,2,\ldots$ compute the distribution function $\mathrm{P}\left(N_{\rm th}\right)$ from the cosmological model and a user-specified star formation efficiency.

• 3.

  Compute

  $$\mathrm{P}\left(N_{\rm th}\geq N_{\rm obs}\right)=\sum_{N_{\rm obs}=0}^{\infty}\mathrm{P}\left(N_{\rm obs}\right)\sum_{N_{\rm th}=N_{\rm obs}}^{\infty}\mathrm{P}\left(N_{\rm th}\right).$$

Explicit evaluation of $\mathrm{P}\left(N_{\rm th}\geq N_{\rm obs}\right)$ requires perfect knowledge about the PDFs of statistical/systematic errors of the stellar masses and redshifts of the candidate galaxies. While the PDFs of statistical errors can be recovered by going back to the SED-fitting process, it is not easy to give an accurate description of the PDFs of systematic errors, which often dominate the error budget. Without further knowledge about the systematic errors, we simply model the PDFs with a few generic distribution functions: skew normal, skew triangular, and skew rectangular, whose definitions are given in the Appendix. The code is organized in such a way that the users can easily add their own distribution functions. We also allow the users to turn off systematic errors by using Dirac delta distribution. In the case when systematic errors are given in the form of extremes from different SED-fitting methods, we match the extremes with exact upper/lower bounds of skew triangular/rectangular distributions or $2\sigma$ bounds of skew normal distribution. Finally, the stellar mass $M_{\star}$ and redshift $z$ from different SED-fitting methods are typically positively correlated. We test the impact of $M_{\star}$-$z$ correlation by forcing $(M_{\star}-M_{\rm\star,median})(z-z_{\rm median})\geq 0$ in the joint distribution of $M_{\star}$ and $z$, and find that the correlation does not make a qualitative difference.

The evaluation of $\mathrm{P}\left(N_{\rm th}\right)$ involves the extended Press-Schechter ellipsoidal collapse model (Press & Schechter 1974; Sheth et al. 2001; Sheth & Tormen 2002), where the comoving halo number density per mass interval, namely the halo mass function (at a fixed redshift $z$) is given by

The $\rho_{m}$ on the right-hand side is the comoving average background matter density and $\sigma$ is the mass fluctuation at scale $R=\left(\frac{3M}{4\pi\rho_{m}}\right)^{1/3}$, given by

where $k$ is the comoving wavenumber, and $\mathrm{P}_{m}(k,z)$ is the linear matter power spectrum at redshift $z$. For the filter function $\mathrm{W}(kR)$, we take the form of top-hat

The simulation calibrated factor $f(\nu)$ is given by

with $A=0.322$, $a=0.707$, $q=0.3$ and $\nu=\frac{\delta_{c}}{\sigma(R,z)}$, where parameter $\delta_{c}=1.686$ corresponds to the critical linear overdensity. We ignore the tiny difference in $\delta_{c}$ for different cosmological models, as we are only interested in the models that are close to the concordance $\Lambda$CDM model.

For a massive galaxy with stellar mass $M_{\star}>M_{\rm\star,cut}$, the mass of its host halo must exceed $M_{\rm cut}=M_{\rm\star,cut}/(\epsilon f_{b})$. The halo mass function Eq. (1) allows to compute the ensemble average of the mean number of such halos

where $\frac{\mathrm{d}V}{\mathrm{d}z\mathrm{d}\Omega}$ is comoving volume per redshift interval per solid angle. Here we have equalized the number of massive galaxies and the number of their host halos, assuming each massive halo has a central galaxy.

For the purpose of using rare objects to study cosmological tensions, we are only interested in the $\langle N_{\rm th}\rangle\lesssim O(1)$ case, where $N_{\rm th}$ approximately follows a Poisson distribution

where $\lambda\approx\langle N_{\rm th}\rangle$ is given in Eq. (5).

To obtain a more accurate result, we integrate $\lambda$ on the right-hand side of Eq. (6) over a Gaussian distribution, giving

The beyond-Poisson cosmic variance from the large scale structure, $\sigma_{\lambda}$, can be read out from the online emulator at https://www.ph.unimelb.edu.au/~mtrenti/cvc/CosmicVariance.html (Trenti & Stiavelli 2008). For a pencil-like survey volume, we find the cosmic variance approximately equal to the product of the linear result, which equals to the product of linear halo bias and matter cosmic variance, and a fudge factor $\approx 1.4$ that accounts for the nonlinear correction. Since the cosmic-variance correction is a subdominant effect¹¹1Our result does not contradict that of Chen et al. (2023), which claims that cosmic variance taken from their simulations is a dominant factor. Our “cosmic variance” here, as well as that defined in Trenti & Stiavelli (2008), refers to the variance due to power spectrum uncertainties caused by finite sampling and does not include Poisson shot noise, while “cosmic variance” taken from simulations include both effects. , this approximation suffices and allows us to make our tool self-contained.

The algorithm is realized with a Python tool publicly available at http://zhiqihuang.top/codes/massivehalo.htm.

We now apply the tool to the 13 candidate galaxies in Labbé et al. (2023). The summary statistics of stellar mass and redshift are given in the Extended Data Table 2 of Labbé et al. (2023), as well as in the online repository given in Sec. LABEL:sec:data. Stellar mass and redshift values include two uncertainties: $\pm(\text{ran})\pm(\text{sys})$, where the random uncertainty (ran) corresponds to the $16$th and $84$th percentiles of the combined posterior distribution, and the systematic uncertainty (sys) corresponds to the extreme values from seven different methods (EAZY, Prospector, and Bagpipes with 5 variations, including dust, SFH, age prior, and SNR limit). In this study, we update the data of three candidate galaxies (L23-13050, L23-35300 and L23-39575) with the information from spectroscopic follow-up observations (Kocevski et al. 2023; Pérez-González et al. 2023; Fujimoto et al. 2023). The updated data is also available in the online repositary given in Sec. LABEL:sec:data.

Unless otherwise specified, we choose the stellar-mass threshold to be $M_{\rm\star,cut}=10^{10}M_{\odot}$ and the redshift range to be $z_{\min}=7$ and $z_{\rm max}=10$, and we assume the statistical errors follow a skew normal distribution. The data covers a sky area $38\,\mathrm{arcmin}^{2}$, which translates to $f_{\rm sky}=2.56\times 10^{-7}$. We take the cosmology to be the Planck best-fit $\Lambda$CDM model (Aghanim et al. 2020) with $\Omega_{m}=0.3158$, $\Omega_{b}=0.049389$, primordial spectrum index $n_{s}=0.96605$, matter fluctuation parameter $\sigma_{8}=0.812$, and Hubble constant $H_{0}=67.32\,\mathrm{km/s/Mpc}$. Figure 2 shows the dependence of $\mathrm{P}\left(N_{\rm th}\geq N_{\rm obs}\right)$ on the star formation efficiency and the model of the PDFs of systematic errors. For $\epsilon\gtrsim 0.3$ that is fairly possible at the early epoch of galaxy formation (Zhang et al. 2022; Qin et al. 2023), the data is well consistent with $\Lambda$CDM model even when the systematic errors are ignored (dotted blue line in Figure 2). For the most conservative $\epsilon\sim 0.1$ case, $\Lambda$CDM model is slightly disfavored by the data. For different forms of the PDFs of the systematic errors, the statistical significance of the tension varies between $\sim 2\sigma$ and $\sim 3\sigma$.

It has been shown that dynamic dark energy model with equation of state $w(a)=w_{0}+w_{a}(1-a)$ (Chevallier & Polarski 2001; Linder 2003), where $a$ is the scale factor, may alleviate the tension between cosmology and the Labbé et al. (2023) data. Since the tension is only significant when the star formation efficiency $\epsilon$ is $\lesssim 0.1$, we hereafter fix $\epsilon=0.1$ and investigate whether the dynamic dark energy model is well consistent with the data. For simplicity, we also assume that the systematic errors follow the skew triangular distribution. We take cosmological parameters from the Markov chains that are sampled with Planck+BAO+SNe likelihood. Here Planck refers to the TTTEEE+lowl+LowE likelihood of the cosmic microwave background maps measured by the Planck satellite (Aghanim et al. 2020); BAO refers to the baryon acoustic oscillations data that are summarized in  Alam et al. (2017); SNe stands for the Pantheon supernova catalog (Scolnic et al. 2018). We download the chains from Planck Legacy Archive (https://pla.esac.esa.int/pla/) and, for better visual effect, thin the chains by a factor of 5. For each set of cosmological parameters of the thinned chains, we compute $\mathrm{P}\left(N_{\rm th}\geq N_{\rm obs}\right)$ value and show it in Figure 3 as a dot, whose color represents the value of $\mathrm{P}\left(N_{\rm th}\geq N_{\rm obs}\right)$ and coordinates indicate the dark energy parameters $(w_{0},w_{a})$. The variation of the dynamic dark energy equation of state in the range that is allowed by Planck+BAO+SNe does lead to a variation of $\mathrm{P}\left(N_{\rm th}\geq N_{\rm obs}\right)$, however, not at a very significant level that would make the result qualitatively different.

In this work, we present and make publicly available an efficient numeric tool that can quickly estimate the tension between cosmological models and observed high-redshift massive galaxies. Our method is based on the statistics of the total numbers of massive halos, whereas the previous works are all based on the statistics of the cumulative stellar mass $\rho_{\star}(>M_{\star})$ (Labbé et al. 2023; Menci et al. 2022; Boylan-Kolchin 2023; Wang et al. 2023b; Chen et al. 2023). The advantage of using the number of halos as a fundamental variable is that the dominating Poisson shot noise contribution can be explicitly formulated.

We update the parameters of some of the candidate galaxies in Labbé et al. (2023) with recent spectroscopic follow-up observations. For a low star formation efficiency $\epsilon\sim 0.1$, the data remains in $\gtrsim 2\sigma$ tension with the concordance $\Lambda$CDM cosmology, and the tension is insensitive to the dark energy equation of state and other cosmological parameters, provided that the cosmological parameters are confined by the Planck+BAO+SNe likelihood. We also verify that the updated spectroscopic redshift information slightly reduces the tension between $\Lambda$CDM model and the JWST candidate galaxies. For instance, assuming a skew normal distribution of both systematic and statistical errors and star formation efficiency $\epsilon=0.1$, we find Planck best-fit $\Lambda$CDM model is ruled out at $2.59\sigma$ level with the original data from Labbé et al. (2023) and $2.49\sigma$ level with the updated data.

Although the algorithm in Figure 1 captures the main physics, there are subtleties beyond the scope of the present work. First, it is known that the extended Press-Schechter halo mass functions at high redshift are 20%-50% higher than those obtained from N-body simulations (Reed et al. 2003; Despali et al. 2016; Shirasaki et al. 2021; Wang et al. 2022). Secondly, backsplash halos, whose dark matter has left the host halo while its gas may remain in the host halo (Balogh et al. 2000; Wang et al. 2009; Wetzel et al. 2014; Diemer 2021; Wang et al. 2023a), may make the interpretation of star formation efficiency ambiguous in the real world and may alter the statistical significance of the cosmological tension (Chen et al. 2023). Finally, we have ignored catastrophic redshift errors which, despite being rare, might have a significant impact on the counting of rare objects. We leave further studies on these nontrivial subjects as our future work.