However, MCMC cannot provide an uncertainty on these maximum likelihood parameters, limiting its use. Worse, the model uncertainty cannot be derived from the posterior widths, as they also are covariant and so the resulting error envelopes will be meaningless. The most obvious choice then is to compute the medians and widths of posterior distributions of the chains after burn-in. However, the curve traced by the median is not guaranteed to reflect the actual model function, and the 68% confidence envelopes may not contain the maximum likelihood nor the median posterior. We therefore strongly advise that models are reconstructed using the maximum likelihood parameters, and that caution should be exercised when interpreting best-fit parameters from median posteriors. Thankfully, the situation is less severe for symmetric posteriors, which for the total SMF are generally symmetric at $z \le 6$. At $6.5 < z \le 7.5$, the posteriors are highly skewed as the relatively weak constraints produce widespread parameter covariance. The maximum likelihood results in all cases appear reliable. The evolution of the fitted Schechter models for both the median posterior and maximum likelihood is shown in Fig. 10.

In addition to the fits at fixed redshift, we also use the method of Leja et al. (2019a) to fit both the shape and evolution of the SMF simultaneously. This aptly named “continuity model” is fitted on the unbinned distribution of mass-complete sources in $z$ and $M^8$. Double Schechter parameters $M^*$, $\Phi_1$, and $\Phi_2$ are treated as continuous functions over time described by second order polynomial expansions in $z$. The slopes $\alpha_1$ and $\alpha_2$ are assumed to be constant. These 11 parameters are joined by a minor parameter $\sigma_{\text{eff}}$ which attempts to capture the cosmic variance. The effects of Eddington bias are incorporated by resampling the entire catalog by 10 random draws from their $\mathcal{L}(M|z)$. Only the mass-complete points are fitted, which allows sources near the mass limit to scatter in and out of the SMF realizations. Importantly, this method accounts for the significant covariances in the Schechter function between adjacent redshift bins which is neglected when binning in redshift, and exploits it to provide generally tighter parameter constraints (see Leja et al. 2019a for details).

We employ EMCEE to determine the posterior distributions and maximum likelihood fit of the continuity model, noting the aforementioned caveats. Initially we tried to fit the entire SMF out of $z = 7.5$, but the fits were inaccurate due to the limited range of evolution which can be described by a second order expansion. Already known from Davidzon et al. (2017), $\Phi$ rises quickly before slowing down toward $z \sim 0$. The second order expansion in $z$ can either fit the quick rise or the slow down, but not both. We therefore only consider galaxies at $z \le 3$ in our continuity model fit, as Leja et al. (2019a) do, and leave modifications to the fitting functions to future work. The expansion shapes are determined not by their coefficients, but rather by the amplitudes of three anchor points at fixed redshifts. The location of these anchor points is largely arbitrary, and so we follow (Leja et al. 2019a) in choosing $z = 0.2, 1.6$, and $3.0$. Given the larger parameter space, we randomly initialize 100 walkers which explore the space until converged as before. We cannot apply the $1/V_{\text{max}}$ correction directly to the continuity model, and so to remain conservative we exclude sources in the lowest 0.25 dex from the $M_{\text{lim}}$ given in Eq. (3) for the total sample.

The evolution of the Schechter parameters inferred from the star-forming and quiescent SMFs are shown in the rightmost panels in Fig. 11. Tables C.2 and C.3 record the best-fit parameters,

with fits shown in Figs. C.2 and C.3, respectively. Adapting the continuity model to the star-forming and quiescent SMFs is left to future work.

The treatment for the star-forming SMF fit is similar to that of the total form. We begin with a double Schechter form at $0.2 < z \le 0.5$, and transition to a single Schechter form at $2.5 < z \le 5.5$ with $\alpha_1$ fixed. The best-fit model describes the observed SMF reasonably well until $z \approx 3.5-5.5$ where an excess of high-mass star-forming systems is observed that cannot be described by a single Schechter form. Possible reasons for this are discussed in Sect. 6.3.

The quiescent SMF behaves noticeably differently from that of the star-forming sample and therefore requires a different strategy. We begin with a double Schechter form at $0.2 < z \le 1.5$ but fix $\alpha_2$ at $1.1 < z \le 1.5$ as the constraints on the low-mass regime deteriorate. With no significant low-mass constraints, we transition to a single Schechter at $1.5 < z \le 5.5$. We note that this disregards a possible low-mass contribution and so the extrapolation to low-masses is likely an underestimate at $z > 1.5$, depending on the growth of low-mass quiescent galaxies beyond our mass limit (see Santini et al. 2022). The normalization of the quiescent SMF continues to steadily decline with $z$. Given the uncertainty about sample completeness at $4.5 < z \le 5.5$, we necessarily fix $\alpha_1$ to ensure a reliable fit of the secure measurements at high-$M$. Consequently, we find a mild rise in the expected $M^*$ from $z \approx 3 \to 5$.

The evolution of the inferred total, star-forming, and quiescent SMFs resulting from the fixed redshift and continuity model fits are shown in Fig. 10. Furthermore, the evolution of the Schechter function parameters inferred from each sample are shown in Fig. 11. Although the evolutionary physics inferred from the fitting will be discussed in Sect. 6.1, we briefly describe the immediate result here. In general, parameters derived at fixed redshift from the $\chi^2$, maximum likelihood, and median posterior agree well for the total, star-forming, and quiescent SMFs, with few exceptions.

The characteristic mass is found to be $M^* \approx 10^{10.7-11.0} M_{\odot}$ with very little significant evolution until $z \approx 3$ when it begins to decrease with increasing $z$, with the continuity model suggesting a potentially increasing value with $z$. This contrasts with results from Davidzon et al. (2017), who find an initial decrease in $M^*$ which increases at $z > 3$. The results of the continuity model of Leja et al. (2019a) agree well with the lack of evolution suggested by our measurements, although in some tension with our continuity model results. $M^*$ of the star-forming sample similar in shape but generally larger than that of the quiescent sample. The likeness to the total $M^*$ at fixed redshift is proportional to the sub sample that dominates around the knee in that $z$ interval.

The low-mass slope of the low-mass Schechter component $\alpha_1$ is roughly constant with time, although may experience a maximum at $z \approx 1$, in agreement with Davidzon et al. (2017). However, whereas Davidzon et al. (2017) finds a steeply decreasing $\alpha_1$ at $2 < z \le 5$, we find a constant evolution which we then fix at $z > 2.5$ to avoid mounting degeneracies as the low-$M$ constraints are lost; this is also where a single Schechter becomes the preferred form. The value of $\alpha_1$ in good agreement with that found from our continuity model fit, as well as that of Leja et al. (2019a) who formed a composite SMF with significantly lower-$M$ and hence $\alpha_1$ constraints. While $\alpha_1$ slope of the star-forming sample is similar to that of the total sample, the $\alpha_1$ slope of the quiescent sample is poorly constrained and appears to decrease with redshift (i.e. steepen with time), consistent with the rapid appearance of low-$M$ quiescent systems at late times.

<footnote n="8">⁸ This is not formally an STY method, as the normalization is constrained as a fitted parameter.</footnote>
However, MCMC cannot provide an uncertainty on these maximum likelihood parameters, limiting its use. Worse, the model uncertainty cannot be derived from the posterior widths, as they also are covariant and so the resulting error envelopes will be meaningless. The most obvious choice then is to compute the medians and widths of posterior distributions of the chains after burn-in. However, the curve traced by the median is not guaranteed to reflect the actual model function, and the 68% confidence envelopes may not contain the maximum likelihood nor the median posterior. We therefore strongly advise that models are reconstructed using the maximum likelihood parameters, and that caution should be exercised when interpreting best-fit parameters from median posteriors. Thankfully, the situation is less severe for symmetric posteriors, which for the total SMF are generally symmetric at $z \le 6$. At $6.5 < z \le 7.5$, the posteriors are highly skewed as the relatively weak constraints produce widespread parameter covariance. The maximum likelihood results in all cases appear reliable. The evolution of the fitted Schechter models for both the median posterior and maximum likelihood is shown in Fig. 10.

In addition to the fits at fixed redshift, we also use the method of Leja et al. (2019a) to fit both the shape and evolution of the SMF simultaneously. This aptly named “continuity model” is fitted on the unbinned distribution of mass-complete sources in $z$ and $M^8$. Double Schechter parameters $M^*$, $\Phi_1$, and $\Phi_2$ are treated as continuous functions over time described by second order polynomial expansions in $z$. The slopes $\alpha_1$ and $\alpha_2$ are assumed to be constant. These 11 parameters are joined by a minor parameter $\sigma_{\text{eff}}$ which attempts to capture the cosmic variance. The effects of Eddington bias are incorporated by resampling the entire catalog by 10 random draws from their $\mathcal{L}(M|z)$. Only the mass-complete points are fitted, which allows sources near the mass limit to scatter in and out of the SMF realizations. Importantly, this method accounts for the significant covariances in the Schechter function between adjacent redshift bins which is neglected when binning in redshift, and exploits it to provide generally tighter parameter constraints (see Leja et al. 2019a for details).

We employ EMCEE to determine the posterior distributions and maximum likelihood fit of the continuity model, noting the aforementioned caveats. Initially we tried to fit the entire SMF out of $z = 7.5$, but the fits were inaccurate due to the limited range of evolution which can be described by a second order expansion. Already known from Davidzon et al. (2017), $\Phi$ rises quickly before slowing down toward $z \sim 0$. The second order expansion in $z$ can either fit the quick rise or the slow down, but not both. We therefore only consider galaxies at $z \le 3$ in our continuity model fit, as Leja et al. (2019a) do, and leave modifications to the fitting functions to future work. The expansion shapes are determined not by their coefficients, but rather by the amplitudes of three anchor points at fixed redshifts. The location of these anchor points is largely arbitrary, and so we follow (Leja et al. 2019a) in choosing $z = 0.2, 1.6$, and $3.0$. Given the larger parameter space, we randomly initialize 100 walkers which explore the space until converged as before. We cannot apply the $1/V_{\text{max}}$ correction directly to the continuity model, and so to remain conservative we exclude sources in the lowest 0.25 dex from the $M_{\text{lim}}$ given in Eq. (3) for the total sample.

The evolution of the Schechter parameters inferred from the star-forming and quiescent SMFs are shown in the rightmost panels in Fig. 11. Tables C.2 and C.3 record the best-fit parameters,

with fits shown in Figs. C.2 and C.3, respectively. Adapting the continuity model to the star-forming and quiescent SMFs is left to future work.

The treatment for the star-forming SMF fit is similar to that of the total form. We begin with a double Schechter form at $0.2 < z \le 0.5$, and transition to a single Schechter form at $2.5 < z \le 5.5$ with $\alpha_1$ fixed. The best-fit model describes the observed SMF reasonably well until $z \approx 3.5-5.5$ where an excess of high-mass star-forming systems is observed that cannot be described by a single Schechter form. Possible reasons for this are discussed in Sect. 6.3.

The quiescent SMF behaves noticeably differently from that of the star-forming sample and therefore requires a different strategy. We begin with a double Schechter form at $0.2 < z \le 1.5$ but fix $\alpha_2$ at $1.1 < z \le 1.5$ as the constraints on the low-mass regime deteriorate. With no significant low-mass constraints, we transition to a single Schechter at $1.5 < z \le 5.5$. We note that this disregards a possible low-mass contribution and so the extrapolation to low-masses is likely an underestimate at $z > 1.5$, depending on the growth of low-mass quiescent galaxies beyond our mass limit (see Santini et al. 2022). The normalization of the quiescent SMF continues to steadily decline with $z$. Given the uncertainty about sample completeness at $4.5 < z \le 5.5$, we necessarily fix $\alpha_1$ to ensure a reliable fit of the secure measurements at high-$M$. Consequently, we find a mild rise in the expected $M^*$ from $z \approx 3 \to 5$.

The evolution of the inferred total, star-forming, and quiescent SMFs resulting from the fixed redshift and continuity model fits are shown in Fig. 10. Furthermore, the evolution of the Schechter function parameters inferred from each sample are shown in Fig. 11. Although the evolutionary physics inferred from the fitting will be discussed in Sect. 6.1, we briefly describe the immediate result here. In general, parameters derived at fixed redshift from the $\chi^2$, maximum likelihood, and median posterior agree well for the total, star-forming, and quiescent SMFs, with few exceptions.

The characteristic mass is found to be $M^* \approx 10^{10.7-11.0} M_{\odot}$ with very little significant evolution until $z \approx 3$ when it begins to decrease with increasing $z$, with the continuity model suggesting a potentially increasing value with $z$. This contrasts with results from Davidzon et al. (2017), who find an initial decrease in $M^*$ which increases at $z > 3$. The results of the continuity model of Leja et al. (2019a) agree well with the lack of evolution suggested by our measurements, although in some tension with our continuity model results. $M^*$ of the star-forming sample similar in shape but generally larger than that of the quiescent sample. The likeness to the total $M^*$ at fixed redshift is proportional to the sub sample that dominates around the knee in that $z$ interval.

The low-mass slope of the low-mass Schechter component $\alpha_1$ is roughly constant with time, although may experience a maximum at $z \approx 1$, in agreement with Davidzon et al. (2017). However, whereas Davidzon et al. (2017) finds a steeply decreasing $\alpha_1$ at $2 < z \le 5$, we find a constant evolution which we then fix at $z > 2.5$ to avoid mounting degeneracies as the low-$M$ constraints are lost; this is also where a single Schechter becomes the preferred form. The value of $\alpha_1$ in good agreement with that found from our continuity model fit, as well as that of Leja et al. (2019a) who formed a composite SMF with significantly lower-$M$ and hence $\alpha_1$ constraints. While $\alpha_1$ slope of the star-forming sample is similar to that of the total sample, the $\alpha_1$ slope of the quiescent sample is poorly constrained and appears to decrease with redshift (i.e. steepen with time), consistent with the rapid appearance of low-$M$ quiescent systems at late times.

<footnote n="8">⁸ This is not formally an STY method, as the normalization is constrained as a fitted parameter.</footnote>
The normalization of the low-mass Schechter component $\Phi_1$ has little evolution at $0.2 < z \le 1.0$, afterwards rapidly decreasing until it appears to approach zero asymptotically. Although Davidzon et al. (2017) finds slightly larger values of $\Phi_1$ at low-$z$, both measurements generally agree on the rapid decline of the low-mass normalization. $\Phi_1$ derived from the continuity model finds still lower values and hence a more gradual evolution<fnref n="9" />. Yet, $\Phi_1$ derived by Leja et al. (2019a) agrees well with our fixed redshift measurements. As with $\alpha_1$, the $\Phi_1$ normalization of the star-forming sample is similar to that of the total sample and $\Phi_1$ of the quiescent sample remains significantly smaller and decreases with redshift.

The low-mass slope of the high-mass Schechter component $\alpha_2$ is statistically consistent with being constant with $z$, but may rise slightly toward $z \approx 2$. This is broadly comparable with Davidzon et al. (2017) within the stated uncertainties. $\alpha_2$ derived from our continuity model agrees well with the fixed redshift measurements, and with that of Leja et al. (2019a). Interestingly, at $z \le 3$ we find that $\alpha_2-\alpha_1 \approx 1$, which is predicted by Peng et al. (2010) as an indicator of mass quenching. The $\alpha_2$ slope of the star-forming sample is similar to that of the total sample, with $\alpha_2$ of the quiescent sample changing from negative at $z \le 1$ to positive at $1 \le z \le 2$ where it appears to stabilize.

Lastly, the normalization of the high-mass Schechter component $\Phi_2$ generally declines over $z = 0.2 \rightarrow 3.0$. Interestingly, the continuity model finds an initial increase in $\Phi_2$ toward $z \approx 1$, declining afterwards in agreement with the fixed redshift estimates. We find lower values compared to those of Davidzon et al. (2017) and Leja et al. (2019a), which are in general agreement in part because they both use COSMOS2015. This suggests that the different evolution in $\Phi_2$ between this work and the literature may be driven by subtle differences between the catalogs, sample selection, assumed Eddington kernels, or a combination of the three. The $\Phi_2$ normalization of the star-forming sample is subdominant to that of the quiescent sample with both decreasing steadily with time, consistent with the large fraction of quiescence found in massive systems out to $z \approx 2.5$.

Overall we find good agreement with the most directly comparable literature measurements from Davidzon et al. (2017) and Leja et al. (2019a). The implications for these evolutionary trends, and their relation to the growth of massive quiescent galaxies, are explored in detail in Sects. 6.3 and 6.4.

# 6. Discussion

In this section we focus on some of the open issues in galaxy evolution that can be addressed from a statistical perspective, using the results of Sect. 5 as a starting point. Besides the assembly of $z \approx 3-5$ massive galaxies and their quiescent fractions, we also discuss derived measurements such as the cosmic stellar mass density $\rho_*$, as well as the relation between stellar and dark matter halo mass functions. We include comparisons with simulations to provide physical insight and also to highlight areas for improvement from both sides.

## 6.1. Cosmic stellar mass density evolution

Galaxy mass assembly and growth is inextricably related to star formation. As such, reconciling direct measurements of

galaxy mass growth with the behavior of the star formation main sequence (e.g., Brinchmann et al. 2004; Noeske et al. 2007; Daddi et al. 2007; Salim et al. 2007; Whitaker et al. 2012, 2014) is of great interest. Meaningful assessment of empirical models of galaxy growth that link the two (e.g., Peng et al. 2010; Behroozi et al. 2019) depends on unbiased, accurate measurements of the integrated mass density $\rho_*$ and its evolution with $z$.

We integrate the SMF measurements presented in Sect. 5 to derive an estimate of the stellar mass density ($\rho_*$) for each bin of $z$. Although definitions vary, $\rho_*$ is commonly integrated down to $10^8 M_{\odot}$. Since our $M_{\lim}$ at all redshifts is larger than $10^8 M_{\odot}$, we integrate into the extrapolated low-$M$ regime of our (unconvolved) Schechter models. The resulting $\rho_*$ for the total sample in Fig. 10 is compared with other literature measurements in Fig. 12, converting to a Chabrier IMF where relevant. All SMF-based literature measurements of $\rho_*$ have been reintegrated consistently to the same mass limit. We also show $\rho_*$ derived from integrating the star formation rate density (SFRD) function of Madau & Dickinson (2014), assuming a 41% return fraction corresponding to Chabrier’s IMF (see Sect. 6.1 of Ilbert et al. 2010).

We find remarkably good agreement with previous observational studies. At $z \le 3$, the agreement seems to be limited by systematics as these are generally well measured, secure samples. However, $\rho_*$ at $z \ge 3$ is dominated by significantly less certain measurements, due both to the size of the samples and their typically noisy photometry. Our measurements place $\rho_*$ near the midpoint of the scatter, in closest agreement with Davidzon et al. (2017, $z \le 5$) and Grazian et al. (2015, $z \le 7$). The relatively large uncertainties on $\rho_*$ beyond $z \approx 4$ stem from the increasing degeneracy of $M^*$ and $\phi$ with decreasing low-$M$ constraints ($\alpha$ being fixed at $z > 1.5$; see Fig. C.1).

At $z \ge 5$, our estimated median $\rho_*$ as well as that of Grazian et al. (2015) is lower than the SFRD-derived predictions of Madau & Dickinson (2014). This discrepancy could suggest that while star formation is high, the mass growth is lagging behind. As intriguing as this is, we stress that the measurements are consistent within $1\sigma$, that the SFRD constructed in Madau & Dickinson (2014) are constrained by only two surveys at $z > 5$ available at the time (Bouwens et al. 2012a,b; Bowler et al. 2012), and that our measurements assume a fixed low-mass slope $\alpha$. Further work utilizing larger samples complete to lower masses will be required to confidently evaluate this divergence at $z > 7$.

In addition to $\rho_*$ of the total sample, the lower panel of Fig. 12 includes $\rho_*$ for the star-forming and quiescent samples. While stellar mass is overwhelmingly concentrated in star-forming systems from $z \approx 7 \rightarrow 3$, the mass density of quiescent systems grows rapidly until flattening at $z \approx 1$, consistent with Fig. 8. From $z = 1 \rightarrow 0$, there is no growth in $\rho_*$ for either star-forming and quiescent systems where the former remains larger than the latter.

While $\rho_*$ of the star-forming sample is well constrained out to lower masses, the same cannot be said about $\rho_*$ of the quiescent sample as its SMF is less mass complete in comparison. Since we cannot directly determine its low-$M$ shape, we assume that the low-$M$ Schechter component effectively vanishes at $z \ge 1.5$. If this is not the case, then we underestimate $\rho_*$ for the quiescent sample. This may help explain the differences observed with respect to UNIVERSE MACHINE from Behroozi et al. (2019) who report generally larger values of $\rho_*$ for their quiescent sample, which features a significant quiescent population at low $M$ (see Fig. B.1). Although in agreement at $0.2 < z \le 0.5$, (Behroozi et al. 2019) report higher quiescent fractions at $<10^{10.7} M_{\odot}$ with redshift, finding more

<footnote n="9">⁹ We note that a second order expansion cannot describe an exponential tail as seen here, which justifies $z \sim 3$ as the rightmost domain limit reachable by our 12 parameter continuity model.</footnote>
The normalization of the low-mass Schechter component $\Phi_1$ has little evolution at $0.2 < z \le 1.0$, afterwards rapidly decreasing until it appears to approach zero asymptotically. Although Davidzon et al. (2017) finds slightly larger values of $\Phi_1$ at low-$z$, both measurements generally agree on the rapid decline of the low-mass normalization. $\Phi_1$ derived from the continuity model finds still lower values and hence a more gradual evolution<fnref n="9" />. Yet, $\Phi_1$ derived by Leja et al. (2019a) agrees well with our fixed redshift measurements. As with $\alpha_1$, the $\Phi_1$ normalization of the star-forming sample is similar to that of the total sample and $\Phi_1$ of the quiescent sample remains significantly smaller and decreases with redshift.

The low-mass slope of the high-mass Schechter component $\alpha_2$ is statistically consistent with being constant with $z$, but may rise slightly toward $z \approx 2$. This is broadly comparable with Davidzon et al. (2017) within the stated uncertainties. $\alpha_2$ derived from our continuity model agrees well with the fixed redshift measurements, and with that of Leja et al. (2019a). Interestingly, at $z \le 3$ we find that $\alpha_2-\alpha_1 \approx 1$, which is predicted by Peng et al. (2010) as an indicator of mass quenching. The $\alpha_2$ slope of the star-forming sample is similar to that of the total sample, with $\alpha_2$ of the quiescent sample changing from negative at $z \le 1$ to positive at $1 \le z \le 2$ where it appears to stabilize.

Lastly, the normalization of the high-mass Schechter component $\Phi_2$ generally declines over $z = 0.2 \rightarrow 3.0$. Interestingly, the continuity model finds an initial increase in $\Phi_2$ toward $z \approx 1$, declining afterwards in agreement with the fixed redshift estimates. We find lower values compared to those of Davidzon et al. (2017) and Leja et al. (2019a), which are in general agreement in part because they both use COSMOS2015. This suggests that the different evolution in $\Phi_2$ between this work and the literature may be driven by subtle differences between the catalogs, sample selection, assumed Eddington kernels, or a combination of the three. The $\Phi_2$ normalization of the star-forming sample is subdominant to that of the quiescent sample with both decreasing steadily with time, consistent with the large fraction of quiescence found in massive systems out to $z \approx 2.5$.

Overall we find good agreement with the most directly comparable literature measurements from Davidzon et al. (2017) and Leja et al. (2019a). The implications for these evolutionary trends, and their relation to the growth of massive quiescent galaxies, are explored in detail in Sects. 6.3 and 6.4.

# 6. Discussion

In this section we focus on some of the open issues in galaxy evolution that can be addressed from a statistical perspective, using the results of Sect. 5 as a starting point. Besides the assembly of $z \approx 3-5$ massive galaxies and their quiescent fractions, we also discuss derived measurements such as the cosmic stellar mass density $\rho_*$, as well as the relation between stellar and dark matter halo mass functions. We include comparisons with simulations to provide physical insight and also to highlight areas for improvement from both sides.

## 6.1. Cosmic stellar mass density evolution

Galaxy mass assembly and growth is inextricably related to star formation. As such, reconciling direct measurements of

galaxy mass growth with the behavior of the star formation main sequence (e.g., Brinchmann et al. 2004; Noeske et al. 2007; Daddi et al. 2007; Salim et al. 2007; Whitaker et al. 2012, 2014) is of great interest. Meaningful assessment of empirical models of galaxy growth that link the two (e.g., Peng et al. 2010; Behroozi et al. 2019) depends on unbiased, accurate measurements of the integrated mass density $\rho_*$ and its evolution with $z$.

We integrate the SMF measurements presented in Sect. 5 to derive an estimate of the stellar mass density ($\rho_*$) for each bin of $z$. Although definitions vary, $\rho_*$ is commonly integrated down to $10^8 M_{\odot}$. Since our $M_{\lim}$ at all redshifts is larger than $10^8 M_{\odot}$, we integrate into the extrapolated low-$M$ regime of our (unconvolved) Schechter models. The resulting $\rho_*$ for the total sample in Fig. 10 is compared with other literature measurements in Fig. 12, converting to a Chabrier IMF where relevant. All SMF-based literature measurements of $\rho_*$ have been reintegrated consistently to the same mass limit. We also show $\rho_*$ derived from integrating the star formation rate density (SFRD) function of Madau & Dickinson (2014), assuming a 41% return fraction corresponding to Chabrier’s IMF (see Sect. 6.1 of Ilbert et al. 2010).

We find remarkably good agreement with previous observational studies. At $z \le 3$, the agreement seems to be limited by systematics as these are generally well measured, secure samples. However, $\rho_*$ at $z \ge 3$ is dominated by significantly less certain measurements, due both to the size of the samples and their typically noisy photometry. Our measurements place $\rho_*$ near the midpoint of the scatter, in closest agreement with Davidzon et al. (2017, $z \le 5$) and Grazian et al. (2015, $z \le 7$). The relatively large uncertainties on $\rho_*$ beyond $z \approx 4$ stem from the increasing degeneracy of $M^*$ and $\phi$ with decreasing low-$M$ constraints ($\alpha$ being fixed at $z > 1.5$; see Fig. C.1).

At $z \ge 5$, our estimated median $\rho_*$ as well as that of Grazian et al. (2015) is lower than the SFRD-derived predictions of Madau & Dickinson (2014). This discrepancy could suggest that while star formation is high, the mass growth is lagging behind. As intriguing as this is, we stress that the measurements are consistent within $1\sigma$, that the SFRD constructed in Madau & Dickinson (2014) are constrained by only two surveys at $z > 5$ available at the time (Bouwens et al. 2012a,b; Bowler et al. 2012), and that our measurements assume a fixed low-mass slope $\alpha$. Further work utilizing larger samples complete to lower masses will be required to confidently evaluate this divergence at $z > 7$.

In addition to $\rho_*$ of the total sample, the lower panel of Fig. 12 includes $\rho_*$ for the star-forming and quiescent samples. While stellar mass is overwhelmingly concentrated in star-forming systems from $z \approx 7 \rightarrow 3$, the mass density of quiescent systems grows rapidly until flattening at $z \approx 1$, consistent with Fig. 8. From $z = 1 \rightarrow 0$, there is no growth in $\rho_*$ for either star-forming and quiescent systems where the former remains larger than the latter.

While $\rho_*$ of the star-forming sample is well constrained out to lower masses, the same cannot be said about $\rho_*$ of the quiescent sample as its SMF is less mass complete in comparison. Since we cannot directly determine its low-$M$ shape, we assume that the low-$M$ Schechter component effectively vanishes at $z \ge 1.5$. If this is not the case, then we underestimate $\rho_*$ for the quiescent sample. This may help explain the differences observed with respect to UNIVERSE MACHINE from Behroozi et al. (2019) who report generally larger values of $\rho_*$ for their quiescent sample, which features a significant quiescent population at low $M$ (see Fig. B.1). Although in agreement at $0.2 < z \le 0.5$, (Behroozi et al. 2019) report higher quiescent fractions at $<10^{10.7} M_{\odot}$ with redshift, finding more

<footnote n="9">⁹ We note that a second order expansion cannot describe an exponential tail as seen here, which justifies $z \sim 3$ as the rightmost domain limit reachable by our 12 parameter continuity model.</footnote>
<figure><img src="image_12.png" /><figcaption>**Fig. 12.** $0.2 < z \le 7.5$ cosmic stellar mass density. Upper: evolution of the cosmic stellar mass density of the total sample computed from the best-fit likelihood models (blue) integrated above $10^8 M_{\odot}$. Literature results from observational studies of mass-selected samples (Grazian et al. 2015; Tomczak et al. 2014; Caputi et al. 2015, 2011; Ilbert et al. 2013; Muzzin et al. 2013; Santini et al. 2012; Adams et al. 2021; McLeod et al. 2021; Wright et al. 2018; Thorne et al. 2021) and mass inferred from rest UV measurements (González et al. 2011). By integrating their SFRD functions, we can plot $\rho_*$ from Behroozi et al. (2013) and Madau & Dickinson (2014). In both cases we assume a return fraction of 41% (based on Chabrier's IMF, see Sect. 6.1 of Ilbert et al. 2010). For Madau & Dickinson (2014), we include a shaded area based on return fractions between 25 and 50% (the latter value is similar to the one given by Salpeter's IMF). Lower: evolution of the cosmic stellar mass density of the total (gray, repeated from above), star-forming (light blue), and quiescent (orange) samples compared to literature measurements (Behroozi et al. 2019; Santini et al. 2021; McLeod et al. 2021).</figcaption></figure>

than $10 \times$ as many $10^9 M_{\odot}$ by $z \approx 1$. UNIVERSE MACHINE primarily defines quiescence as sSFR $< 10^{-11} \text{ yr}^{-1}$, which according to Davidzon et al. (2018) is comparable our *NUVrJ*-selection, and so this overproduction of quiescent galaxies may indeed be at odds with our observations. Although not readily available, a consistently *NUVrJ*-selected sample from UNIVERSE MACHINE would clarify this sensitive comparison.

Comparisons with the results of Santini et al. (2021) and McLeod et al. (2021), although generally in agreement, are also complicated by differences in selection. Whereas Santini et al. (2021) adopts a novel SFR-driven selection, McLeod et al. (2021) adopts a *UVJ* selection following Carnall et al. (2018). Therefore, we stress that these estimates of $\rho_*$ for our quiescent sample are particular to *NUVrJ*-selected quiescent systems, and our assumption of a single Schechter description at $z > 1.5$ may drive our relatively low estimates. Finally, despite differences in selection and constraining power, these observational studies collectively indicate a general agreement as to the evolution of $\rho_*$ for quiescent systems.

## 6.2. Comparison to simulations

Observational constraints on the shape and evolution of the SMF may be useful in their own right, but we can also gain meaningful insight by comparing them with SMFs produced by galaxy formation simulations. Figure 13 shows the SMF constraints from this work and the inferred SMF derived by the binned MCMC fits evaluated at maximum likelihood, with the flagship or reference versions of simulated hydrodynamical SMFs overlaid: TNG100 of the ILLUSTRISTNG project<fnref n="10" />, EAGLE, SIMBA and FLARES (Pillepich et al. 2018; Furlong et al. 2015; Davé et al. 2019; Lovell et al. 2021, respectively). It is worth noting that FLARES uses the same code and model as EAGLE (specifically they adopted the strong AGN feedback EAGLE model introduced in Crain et al. 2015), but sample regions that span a wider dynamical range, and therefore has an much larger effective volume than EAGLE. We also include the semi-analytical results from SHARK (Lagos et al. 2018) as well as those of SANTA CRUZ (Yung et al. 2019a,b, see also Somerville et al. 2015). To note, we do not apply any artificial normalization (aside from *h* considerations) to any model or observation such that direct comparisons are possible from the figures alone.

Since our observed SMF measurements are affected by Eddington bias, it is preferable to compare the SMF directly predicted by the simulations assuming no errors to our inferred SMF with considerations as to the assumptions made in our Eddington bias correction. We acknowledge that each simulation has multiple flavors with variations to their physical recipes. Although these variations provide additional insight, care must be taken when making these more complicated comparisons. As the goal of this work is to provide constraints on the observed and inferred SMF, we reserve comparisons to the variations of simulated SMFs for future work.

Overall, FLARES, SHARK, and EAGLE perform best below $M^*$, with SIMBA, FLARES, and the SANTA CRUZ SAM performing best above $M^*$. At high-$z$ ($z \ge 4$), we find the best agreement with SIMBA, ILLUSTRIS TNG100, and the SANTA CRUZ SAM. At low-$z$ ($z \le 1.5$), however, there is a slight preference toward SIMBA as ILLUSTRIS TNG100 underestimates the number densities at $M^*$ for $z \le 1.5$. The situation is different at $1.5 < z \le 3.0$ where SIMBA significantly overestimates the high-$M$ end. While this could be explained by a systematic bias in observed masses or a missing high-$M$ population (see Sect. 6.3), it could also be due to SIMBA potentially over-grouping several separate dark matter haloes into one massive halo and/or insufficient AGN feedback at early times. Over this same range, it is apparent that EAGLE suffers from volume limitations and thereby does not contain the most massive galaxies. This is most clear

<footnote n="10">¹⁰ We choose TNG100 as a compromise between resolution and volume, see Pillepich et al. (2018), Donnari et al. (2021a,b).</footnote>
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** $0.2 < z \le 7.5$ cosmic stellar mass density. Upper: evolution of the cosmic stellar mass density of the total sample computed from the best-fit likelihood models (blue) integrated above $10^8 M_{\odot}$. Literature results from observational studies of mass-selected samples (Grazian et al. 2015; Tomczak et al. 2014; Caputi et al. 2015, 2011; Ilbert et al. 2013; Muzzin et al. 2013; Santini et al. 2012; Adams et al. 2021; McLeod et al. 2021; Wright et al. 2018; Thorne et al. 2021) and mass inferred from rest UV measurements (González et al. 2011). By integrating their SFRD functions, we can plot $\rho_*$ from Behroozi et al. (2013) and Madau & Dickinson (2014). In both cases we assume a return fraction of 41% (based on Chabrier's IMF, see Sect. 6.1 of Ilbert et al. 2010). For Madau & Dickinson (2014), we include a shaded area based on return fractions between 25 and 50% (the latter value is similar to the one given by Salpeter's IMF). Lower: evolution of the cosmic stellar mass density of the total (gray, repeated from above), star-forming (light blue), and quiescent (orange) samples compared to literature measurements (Behroozi et al. 2019; Santini et al. 2021; McLeod et al. 2021).</figcaption></figure>

than $10 \times$ as many $10^9 M_{\odot}$ by $z \approx 1$. UNIVERSE MACHINE primarily defines quiescence as sSFR $< 10^{-11} \text{ yr}^{-1}$, which according to Davidzon et al. (2018) is comparable our *NUVrJ*-selection, and so this overproduction of quiescent galaxies may indeed be at odds with our observations. Although not readily available, a consistently *NUVrJ*-selected sample from UNIVERSE MACHINE would clarify this sensitive comparison.

Comparisons with the results of Santini et al. (2021) and McLeod et al. (2021), although generally in agreement, are also complicated by differences in selection. Whereas Santini et al. (2021) adopts a novel SFR-driven selection, McLeod et al. (2021) adopts a *UVJ* selection following Carnall et al. (2018). Therefore, we stress that these estimates of $\rho_*$ for our quiescent sample are particular to *NUVrJ*-selected quiescent systems, and our assumption of a single Schechter description at $z > 1.5$ may drive our relatively low estimates. Finally, despite differences in selection and constraining power, these observational studies collectively indicate a general agreement as to the evolution of $\rho_*$ for quiescent systems.

## 6.2. Comparison to simulations

Observational constraints on the shape and evolution of the SMF may be useful in their own right, but we can also gain meaningful insight by comparing them with SMFs produced by galaxy formation simulations. Figure 13 shows the SMF constraints from this work and the inferred SMF derived by the binned MCMC fits evaluated at maximum likelihood, with the flagship or reference versions of simulated hydrodynamical SMFs overlaid: TNG100 of the ILLUSTRISTNG project$^{10}$, EAGLE, SIMBA and FLARES (Pillepich et al. 2018; Furlong et al. 2015; Davé et al. 2019; Lovell et al. 2021, respectively). It is worth noting that FLARES uses the same code and model as EAGLE (specifically they adopted the strong AGN feedback EAGLE model introduced in Crain et al. 2015), but sample regions that span a wider dynamical range, and therefore has an much larger effective volume than EAGLE. We also include the semi-analytical results from SHARK (Lagos et al. 2018) as well as those of SANTA CRUZ (Yung et al. 2019a,b, see also Somerville et al. 2015). To note, we do not apply any artificial normalization (aside from $h$ considerations) to any model or observation such that direct comparisons are possible from the figures alone.

Since our observed SMF measurements are affected by Eddington bias, it is preferable to compare the SMF directly predicted by the simulations assuming no errors to our inferred SMF with considerations as to the assumptions made in our Eddington bias correction. We acknowledge that each simulation has multiple flavors with variations to their physical recipes. Although these variations provide additional insight, care must be taken when making these more complicated comparisons. As the goal of this work is to provide constraints on the observed and inferred SMF, we reserve comparisons to the variations of simulated SMFs for future work.

Overall, FLARES, SHARK, and EAGLE perform best below $M^*$, with SIMBA, FLARES, and the SANTA CRUZ SAM performing best above $M^*$. At high-$z$ ($z \ge 4$), we find the best agreement with SIMBA, ILLUSTRIS TNG100, and the SANTA CRUZ SAM. At low-$z$ ($z \le 1.5$), however, there is a slight preference toward SIMBA as ILLUSTRIS TNG100 underestimates the number densities at $M^*$ for $z \le 1.5$. The situation is different at $1.5 < z \le 3.0$ where SIMBA significantly overestimates the high-$M$ end. While this could be explained by a systematic bias in observed masses or a missing high-$M$ population (see Sect. 6.3), it could also be due to SIMBA potentially over-grouping several separate dark matter haloes into one massive halo and/or insufficient AGN feedback at early times. Over this same range, it is apparent that EAGLE suffers from volume limitations and thereby does not contain the most massive galaxies. This is most clear

<footnote>$^{10}$ We choose TNG100 as a compromise between resolution and volume, see Pillepich et al. (2018), Donnari et al. (2021a,b).</footnote>
<figure><img src="image_13.png" /><figcaption>**Fig. 13.** Comparison of observed and inferred galaxy stellar mass function (gray points, and colored curve with $1 \text{ and } 2\sigma$ envelopes, respectively) to the reference flavors of four simulations: TNG100 of the ILLUSTRISTNG project (Pillepich et al. 2018), EAGLE (Furlong et al. 2015), SHARK (Default; Lagos et al. 2018), SANTA CRUZ (Yung et al. 2019a,b), and SIMBA (Flagship; Davé et al. 2019). To note, we do not apply any artificial normalization to any model or observation. Upper limits for empty bins are shown by the horizontal gray line with an arrow. Mass incomplete measurements are not shown.</figcaption></figure>

from the comparison between EAGLE and FLARES, which are based on the same physical recipes, but the latter sample much rarer overdensities, and hence captures the high-$M$ population. Both the SHARK and Santar Cruz SAMs fare better as their semi-analytical prescriptions are able to produce high-$M$ systems. In some $z$-bins, both SHARK and EAGLE assemble low-$M$ systems too early relative to the observed SMF where there are direct constraints (most apparent at $1 \le z \le 2$). We find a lesser degree of agreement at $z > 3.5$: while both SHARK and EAGLE under-produce the number of galaxies at all $M$, there is significant scatter between ILLUSTRIS TNG, SIMBA, FLARES, and SANTA CRUZ. Meanwhile the volume limitations of ILLUSTRIS TNG100 become significant at $z > 6$ (see Pillepich et al. 2018 for comparisons with TNG300).

It is remarkable that all of the simulations<fnref n="11" /> reproduce the $0.2 < z < 0.5$ SMF, and yet at higher redshifts display severe disagreement with our observations and each other with number densities differing by more than a factor of ten. While this should not come as a surprise given that simulations are typically tuned to reach an end state in agreement with the local Universe (despite known disagreements with the observed sSFR history and merger rates, see Popesso et al. 2023; Conselice et al. 2022), it suggests that their initial conditions and early evolu-

tionary behavior are strikingly different such that variously tuned physical recipes and initial conditions can produce the same end state. This highlights the need for continued observations to improve constraints on the SMF at high-$z$ where simulations can be critically tested, and their physical thereby recipes refined. In addition to the treatment of AGN and associated outflows (e.g., Debuhr et al. 2012; Richardson et al. 2016; Mitchell et al. 2020), modeling of dust attenuation at high-$z$ is a particularly relevant concern at $z \ge 4$ where the dust content of galaxies has been largely unknown, and in turn makes stellar mass estimates based on rest-frame UV/optical light increasingly uncertain. This, coupled with uncertainties in dust production mechanisms at high-$z$, can lead to simulations diverging significantly (from each other and observations), and has often been cited as a major contributor to such discrepancies (see discussions in e.g., Kokorev et al. 2021). Although future observations will benefit from improved more comprehensive mass uncertainties, our current constraints indicate that while it appears that simulations are able to reproduce the observed abundance of high-$M$ galaxies and so the production of early low-$M$ systems needs improvement.

### 6.3. Abundant massive galaxies at $z \sim 3-5$

One of the most striking results highlighted in Fig. 6 is the high number density of massive $M > 10^{11} M_\odot$ galaxies not only at $z > 3.5$. Although few in number, their identification

<footnote n="11">¹¹ With the exception FLARES and the SANTA CRUZ SAM which are limited to $z > 4.5$ and $z > 4.0$, respectively.</footnote>
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Comparison of observed and inferred galaxy stellar mass function (gray points, and colored curve with 1 and $2\sigma$ envelopes, respectively) to the reference flavors of four simulations: TNG100 of the ILLUSTRISTNG project (Pillepich et al. 2018), EAGLE (Furlong et al. 2015), SHARK (Default; Lagos et al. 2018), SANTA CRUZ (Yung et al. 2019a,b), and SIMBA (Flagship; Davé et al. 2019). To note, we do not apply any artificial normalization to any model or observation. Upper limits for empty bins are shown by the horizontal gray line with an arrow. Mass incomplete measurements are not shown.</figcaption></figure>

from the comparison between EAGLE and FLARES, which are based on the same physical recipes, but the latter sample much rarer overdensities, and hence captures the high-$M$ population. Both the SHARK and Santar Cruz SAMs fare better as their semi-analytical prescriptions are able to produce high-$M$ systems. In some $z$-bins, both SHARK and EAGLE assemble low-$M$ systems too early relative to the observed SMF where there are direct constraints (most apparent at $1 \le z \le 2$). We find a lesser degree of agreement at $z > 3.5$: while both SHARK and EAGLE under-produce the number of galaxies at all $M$, there is significant scatter between ILLUSTRIS TNG, SIMBA, FLARES, and SANTA CRUZ. Meanwhile the volume limitations of ILLUSTRIS TNG100 become significant at $z > 6$ (see Pillepich et al. 2018 for comparisons with TNG300).

It is remarkable that all of the simulations<fnref n="11" /> reproduce the $0.2 < z < 0.5$ SMF, and yet at higher redshifts display severe disagreement with our observations and each other with number densities differing by more than a factor of ten. While this should not come as a surprise given that simulations are typically tuned to reach an end state in agreement with the local Universe (despite known disagreements with the observed sSFR history and merger rates, see Popesso et al. 2023; Conselice et al. 2022), it suggests that their initial conditions and early evolu-

tionary behavior are strikingly different such that variously tuned physical recipes and initial conditions can produce the same end state. This highlights the need for continued observations to improve constraints on the SMF at high-$z$ where simulations can be critically tested, and their physical thereby recipes refined. In addition to the treatment of AGN and associated outflows (e.g., Debuhr et al. 2012; Richardson et al. 2016; Mitchell et al. 2020), modeling of dust attenuation at high-$z$ is a particularly relevant concern at $z \ge 4$ where the dust content of galaxies has been largely unknown, and in turn makes stellar mass estimates based on rest-frame UV/optical light increasingly uncertain. This, coupled with uncertainties in dust production mechanisms at high-$z$, can lead to simulations diverging significantly (from each other and observations), and has often been cited as a major contributor to such discrepancies (see discussions in e.g., Kokorev et al. 2021). Although future observations will benefit from improved more comprehensive mass uncertainties, our current constraints indicate that while it appears that simulations are able to reproduce the observed abundance of high-$M$ galaxies and so the production of early low-$M$ systems needs improvement.

### 6.3. Abundant massive galaxies at $z \sim 3-5$

One of the most striking results highlighted in Fig. 6 is the high number density of massive $M > 10^{11} M_\odot$ galaxies not only at $z > 3.5$. Although few in number, their identification

<footnote n="11">¹¹ With the exception FLARES and the SANTA CRUZ SAM which are limited to $z > 4.5$ and $z > 4.0$, respectively.</footnote>
in COSMOS2020 is a direct consequence of utilizing the larger 1.27 deg² now accessible with deep, homogeneous NIR coverage. While no $M > 10^{11}$ M☉ galaxies are observed at $z > 5$, their growth since then has been similar to galaxies at other masses, as shown in Fig. 8. The majority of $M > 10^{11}$ M☉ systems are star-forming at $z > 1$, as shown by Fig. 7 and quantified in Fig. 9, with only $z < 1$ systems at $M \approx 10^{10.5-11}$ M☉ being equally divided between star-forming and quiescent states. We find evidence for a sustained population of massive quiescent systems at $z > 2$, but their number densities are dwarfed by that of star-forming systems by a factor of ~10. The existence of these massive quiescent systems seems to defy the timescales expected for mass assembly (Steinhardt et al. 2016; Faisst et al. 2017; Schreiber et al. 2018a; Cecchi et al. 2019), and so their tremendously early formation and growth are a topic of great interest (Caputi et al. 2011; Toft et al. 2017; Hill et al. 2017; Carnall et al. 2018; Tanaka et al. 2019; Valentino et al. 2020; Whitaker et al. 2021; Akhshik et al. 2023; Marsan et al. 2022), including a more focused investigation into the origins of similarly selected massive galaxies from COSMOS2020 by Gould et al. (2023).

Despite nearly identical sample selections, we still find greater number densities of massive galaxies compared to Davidzon et al. (2017) at $3 < z \le 5.5$. By comparing to COSMOS2015 (matched within 0.6"), we find that 78% of $M > 10^{11}$ M☉ galaxies in our sample are also found in COSMOS2015 where they are exclusively high-$z$ (91% agree within $\Delta z \pm 0.5$) and high-mass (78% agree within $\Delta M \pm 0.5$ M☉). The remaining 23% are new sources found only in COSMOS2020. As shown by Fig. 14, they are predominantly faint, having a median $K_s$ magnitude $\approx 24.2$ AB compared to the median sample brightness $K_s \approx 23.4$ AB with a photometric uncertainty 0.05–0.1 AB (see Fig. 10 of Weaver et al. 2022). They are generally as massive as the already known sources from COSMOS2015 and follow the same $M$ distribution. They are also remarkably red, having a median $H - K_s$ color of $\approx 1.2$ AB compared to that of the sample overall ($\approx 0.8$ AB). It seems unlikely that they are now found because of the better de-blending by The Farmer as visual inspection shows that they are generally isolated sources and are detected by Source Extractor in the CLASSIC version of the catalog. Therefore the most probable explanation for the new faint, red sources is that the deeper UltraVISTA $YJHK_s$ (as well as HSC $iz$) have enabled the detection of these faint red sources which in COSMOS2015 could not be detected. Interestingly, we find that their distribution in stellar mass is consistent with that of the total sample such that number densities at all masses $M > 10^{11}$ M☉ are proportionately represented.

Visual inspection of photometry and SED fits indicate that nearly all of the $3 < z \le 5.5$ $M > 10^{11}$ M☉ galaxies have red colors. About 75% of them are selected as star-forming by NUVrJ, $\ge 80\%$ of which are likely attenuated by a thick screen of dust ($A_V > 1$; compared to only ~10% across entire sample). The remaining ~25% are classified as quiescent. The red colors appear to be genuine and not driven by blends, as confirmed by visually inspecting the imaging for each galaxy. Their broad-band photometry lacks the strong spectral features that contribute to a secure photo-$z$: there is no detectable flux contamination by nebular emission lines and both the Lyman and Balmer breaks are weak. However, they are surprisingly well fit by LePhare (typical $\chi^2_N \approx 1.5$). Their likelihood redshift distributions $L(z)$ are narrow as $> 68\%$ of the probability is typically contained within $\Delta z \approx 1$, with similarly well constrained $L(M|z)$. Recently, Lower et al. (2020) has shown how the relatively simple parametric star formation histories assumed by

most template-based SED fitting codes are susceptible to biases on the order of 0.5 dex in mass, which suggests that these uncertainties are likely underestimated (see also Michałowski et al. 2012, 2014).

One possible explanation for smooth SEDs with a power-law slope is contamination by AGN. The COSMOS2020 results computed by LePhare include classifications for sources with strong X-ray detections<fnref n="12" /> determined from crossmatching to Civano et al. (2016) sources within 0.6" radius. In our sample we only use those sources identified as galaxies, which excludes these X-ray detections. Their inclusion would only serve to increase these already surprisingly large number densities. We do not attempt to quantify this, as the expected accretion disc light from Type 1 Seyferts and quasars make estimates of photo-$z$ and $M$ unreliable and susceptible to catastrophic failures. See Weaver et al. (2022) and Salvato et al. (2011) for details.

Identification of AGN (including X-ray faint AGN) from broadband colors has been explored in the literature (Stern et al. 2005; Donley et al. 2012; Hviding et al. 2022), and their discussion is a standard component of SMF studies at these redshifts (see Grazian et al. 2015; Davidzon et al. 2017), although a consensus has yet to be reached. In general, AGN selection criteria rely on colors derived from (near-)infrared wavelengths from most notably Spitzer/IRAC. While the Donley et al. (2012) criteria have been used successfully at $z < 3$, they require constraints in all four IRAC bands, which is not the case for COSMOS2020 where channel 3 and channel 4 lack sufficient depth to detect

<figure><img src="image_14.png" /><figcaption>**Fig. 14.** $K_s$ magnitude and $H - K_s$ colors of 125 new massive $M > 10^{11}$ M☉ galaxy candidates $3.0 < z \le 5.5$ found in COSMOS2020 (orange), relative to the total sample used in this work (gray) and the 444 galaxies found in COSMOS2015 (blue). Measurements are taken from COSMOS2020, although they are similar to those from COSMOS2015, where matched. Number densities are summarised using gaussian kernel density estimators to produce smoothed distributions and contours for each sample. Median $K_s$ and $H - K_s$ values for each sample are shown by colored dotted lines.</figcaption></figure>

<footnote n="12">¹² Sources with X-ray counterparts (i.e., AGN) are classified in Weaver et al. (2022) by lp_type=2.</footnote>
in COSMOS2020 is a direct consequence of utilizing the larger 1.27 deg² now accessible with deep, homogeneous NIR coverage. While no $M > 10^{11}$ M☉ galaxies are observed at $z > 5$, their growth since then has been similar to galaxies at other masses, as shown in Fig. 8. The majority of $M > 10^{11}$ M☉ systems are star-forming at $z > 1$, as shown by Fig. 7 and quantified in Fig. 9, with only $z < 1$ systems at $M \approx 10^{10.5-11}$ M☉ being equally divided between star-forming and quiescent states. We find evidence for a sustained population of massive quiescent systems at $z > 2$, but their number densities are dwarfed by that of star-forming systems by a factor of ~10. The existence of these massive quiescent systems seems to defy the timescales expected for mass assembly (Steinhardt et al. 2016; Faisst et al. 2017; Schreiber et al. 2018a; Cecchi et al. 2019), and so their tremendously early formation and growth are a topic of great interest (Caputi et al. 2011; Toft et al. 2017; Hill et al. 2017; Carnall et al. 2018; Tanaka et al. 2019; Valentino et al. 2020; Whitaker et al. 2021; Akhshik et al. 2023; Marsan et al. 2022), including a more focused investigation into the origins of similarly selected massive galaxies from COSMOS2020 by Gould et al. (2023).

Despite nearly identical sample selections, we still find greater number densities of massive galaxies compared to Davidzon et al. (2017) at $3 < z \le 5.5$. By comparing to COSMOS2015 (matched within 0.6"), we find that 78% of $M > 10^{11}$ M☉ galaxies in our sample are also found in COSMOS2015 where they are exclusively high-$z$ (91% agree within $\Delta z \pm 0.5$) and high-mass (78% agree within $\Delta M \pm 0.5$ M☉). The remaining 23% are new sources found only in COSMOS2020. As shown by Fig. 14, they are predominantly faint, having a median $K_s$ magnitude $\approx 24.2$ AB compared to the median sample brightness $K_s \approx 23.4$ AB with a photometric uncertainty 0.05–0.1 AB (see Fig. 10 of Weaver et al. 2022). They are generally as massive as the already known sources from COSMOS2015 and follow the same $M$ distribution. They are also remarkably red, having a median $H - K_s$ color of $\approx 1.2$ AB compared to that of the sample overall ($\approx 0.8$ AB). It seems unlikely that they are now found because of the better de-blending by The Farmer as visual inspection shows that they are generally isolated sources and are detected by Source Extractor in the CLASSIC version of the catalog. Therefore the most probable explanation for the new faint, red sources is that the deeper UltraVISTA $YJHK_s$ (as well as HSC $iz$) have enabled the detection of these faint red sources which in COSMOS2015 could not be detected. Interestingly, we find that their distribution in stellar mass is consistent with that of the total sample such that number densities at all masses $M > 10^{11}$ M☉ are proportionately represented.

Visual inspection of photometry and SED fits indicate that nearly all of the $3 < z \le 5.5$ $M > 10^{11}$ M☉ galaxies have red colors. About 75% of them are selected as star-forming by NUVrJ, $\ge 80\%$ of which are likely attenuated by a thick screen of dust ($A_V > 1$; compared to only ~10% across entire sample). The remaining ~25% are classified as quiescent. The red colors appear to be genuine and not driven by blends, as confirmed by visually inspecting the imaging for each galaxy. Their broad-band photometry lacks the strong spectral features that contribute to a secure photo-$z$: there is no detectable flux contamination by nebular emission lines and both the Lyman and Balmer breaks are weak. However, they are surprisingly well fit by LePhare (typical $\chi^2_N \approx 1.5$). Their likelihood redshift distributions $\mathcal{L}(z)$ are narrow as $> 68\%$ of the probability is typically contained within $\Delta z \approx 1$, with similarly well constrained $\mathcal{L}(M|z)$. Recently, Lower et al. (2020) has shown how the relatively simple parametric star formation histories assumed by

most template-based SED fitting codes are susceptible to biases on the order of 0.5 dex in mass, which suggests that these uncertainties are likely underestimated (see also Michałowski et al. 2012, 2014).

One possible explanation for smooth SEDs with a power-law slope is contamination by AGN. The COSMOS2020 results computed by LePhare include classifications for sources with strong X-ray detections<fnref n="12" /> determined from crossmatching to Civano et al. (2016) sources within 0.6" radius. In our sample we only use those sources identified as galaxies, which excludes these X-ray detections. Their inclusion would only serve to increase these already surprisingly large number densities. We do not attempt to quantify this, as the expected accretion disc light from Type 1 Seyferts and quasars make estimates of photo-$z$ and $M$ unreliable and susceptible to catastrophic failures. See Weaver et al. (2022) and Salvato et al. (2011) for details.

Identification of AGN (including X-ray faint AGN) from broadband colors has been explored in the literature (Stern et al. 2005; Donley et al. 2012; Hviding et al. 2022), and their discussion is a standard component of SMF studies at these redshifts (see Grazian et al. 2015; Davidzon et al. 2017), although a consensus has yet to be reached. In general, AGN selection criteria rely on colors derived from (near-)infrared wavelengths from most notably Spitzer/IRAC. While the Donley et al. (2012) criteria have been used successfully at $z < 3$, they require constraints in all four IRAC bands, which is not the case for COSMOS2020 where channel 3 and channel 4 lack sufficient depth to detect

<figure><img src="image_15.png" /><figcaption>**Fig. 14.** $K_s$ magnitude and $H - K_s$ colors of 125 new massive $M > 10^{11}$ M☉ galaxy candidates $3.0 < z \le 5.5$ found in COSMOS2020 (orange), relative to the total sample used in this work (gray) and the 444 galaxies found in COSMOS2015 (blue). Measurements are taken from COSMOS2020, although they are similar to those from COSMOS2015, where matched. Number densities are summarised using gaussian kernel density estimators to produce smoothed distributions and contours for each sample. Median $K_s$ and $H - K_s$ values for each sample are shown by colored dotted lines.</figcaption></figure>

<footnote n="12">¹² Sources with X-ray counterparts (i.e., AGN) are classified in Weaver et al. (2022) by `lp_type=2`.</footnote>
these sources. Even if deeper IRAC images could be taken, the selection criteria rely on the four bands sampling the continuum shape at restframe 2–10 µm but at $z > 3$ instead sample the rest-frame stellar bulk at $\le 2$ µm. While the MIPS 24 µm data from S-COSMOS (Sanders et al. 2007) is an attractive solution, it also is too shallow (20.2 AB at $3\sigma$, Jin et al. 2018) to fully constrain the restframe 2–10 µm continuum at $3 < z \le 5$: only galaxies $H \le 20$ with large AGN fractions can be positively identified, and both low-fraction AGN at $H \approx 20$ and normal galaxies at $H < 20$ cannot be classified with MIPS. Full spectral fitting including far-infrared measurements is challenging without the constraints from channel 3, channel 4, and 24 µm, and attempts to gain further insight using STARDUST (Kokorev et al. 2021) were broadly unsuccessful with tentative contamination found to be on the order of 10–30%. We note, importantly, that removing potentially contaminated sources provides no significant change to the SMF at these epochs.

Despite these challenges, we are able to leverage the elementary AGN template fitting from LePhare to statistically assess the spectral similarity between the best-fit AGN and galaxy templates for each source, limited to rest-frame UV/optical light. While only 10% of galaxies in the total sample are best-fit with AGN templates, this fraction grows to 30% for all mass complete galaxies $3.5 < z < 5.5$ and then to 50% for those with $M > 10^{10.8} M_\odot$ with broad wings, having >80% of these sources statistically indistinguishable as either galaxy or AGN ($|\Delta\chi^2| < 0.5$). Turning to further infrared data of individual sources, we find 15% of the sample are detected by VLA-COSMOS (Smolčić et al. 2017). While not conclusive, we find 5% of these massive $M > 10^{11} M_\odot$ galaxies at $3 < z \le 5$ are detected in the ALMA maps of A$^3$COSMOS (Liu et al. 2019), which currently covers ~5% of the COSMOS field, highlighting the need for further observations. While similar sample statistics for individual sources have been extrapolated in the literature (e.g., Santini et al. 2021), anticipated surveys such as the TolTEC Ultra-Deep Galaxy Survey (Pope et al. 2019) will expand and deepen the far-infrared/(sub)mm coverage in COSMOS so that these samples can be more conclusively investigated.

A recent X-ray and radio stacking analysis of similarly selected $z > 1.5$, $M > 10^{10} M_\odot$ COSMOS2020 galaxies by Ito et al. (2022; selected from independent SED fitting with MIZUKI, Tanaka 2015) revealed that low-luminosity AGN are likely ubiquitous in massive quiescent galaxies from $1.5 < z \le 5$, even if individually they are not X-ray or radio detected. They also estimate the AGN contribution to optical/NIR continuum and find that at these redshifts, e.g., the rest-frame B-band luminosities of their quiescent galaxies are a factor of 30× larger than the expected rest-frame AGN luminosity. Ito et al. (2022) also find that the AGN luminosities for quiescent galaxies are significantly larger than those of star-forming galaxies. Together, their findings provide further confidence that our redshifts and stellar mass estimates for our X-ray undetected sample are not strongly contaminated by AGN light.

The lack of obvious systematics or likely only weak AGN contamination increases our confidence that these sources, or at least a part of them, are truly massive at $z \ge 3$. In agreement with Ito et al. (2022), we also find that >60% of $M > 10^{11} M_\odot$ galaxies appear to be star-forming, and it is likely that at least some of them are also dust obscured (DSFGs; Casey et al. 2014; Zavala et al. 2021, and references therein). Since there are a number of sources found in the deeper NIR images of COSMOS2020, it makes sense that we are now more sensitive to fainter red sources than ever before (see Fig. 14). It

is precisely this class of galaxy which are efficiently captured by infrared facilities such as Spitzer, ALMA, and JWST until now being optically “dark” (e.g., Schreiber et al. 2018b; Wang et al. 2019; Gruppioni et al. 2020; Sun et al. 2021; Fudamoto et al. 2021; Shu et al. 2022). If genuine, their existence not only points at an incompleteness particularly in the massive end of SMFs reported in previous studies lacking the necessarily deep infrared data over degree scales, but also highlights the sensitive interplay between the shape of their SEDs and the selection function consisting of the bands, their depths, and the detection methodology (see Fig. 3 of Weaver et al. 2022). Although detailed simulations will be explored in future work, qualitatively this could explain why the excess of sources relative to a Schechter at $z \approx 3-5$ diminishes at $z \approx 6$ as similarly red sources become too faint to still be detected. Optically dark galaxies selected from 2 mm ALMA observations in COSMOS from Casey et al. (2021) and Manning et al. (2022) are constrained to similar redshifts, stellar masses, and number densities. Importantly, Manning et al. (2022) studied two systems with star formation rates above $200 M_\odot \text{ yr}^{-1}$ but a gas depletion timescale <1 Gyr, suggesting rapid star formation cessation and a potential transition to massive quiescent galaxies by $z \sim 3-4$. Casey et al. (2021) find that DSFGs are responsible for ~30% of the integrated star formation rate density at $3 < z < 6$ and that 2 mm selection is an efficient way to identify larger samples in future surveys (see also Cooper et al. 2022). In an empirically based numerical model, Long et al. (2023) uses known DSFG population properties to predict the DSFG stellar mass function out to $z \sim 5$. They find a shallow power law extension beyond the knee of traditional star-forming Schechter SMFs, and posit that 60–80% of massive ($M > 10^{11} M_\odot$) star-forming galaxies at $z \approx 3-6$ are significantly dust-obscured and not captured in previous SMF studies (see also Martis et al. 2016). This is in line with our aforementioned results on significantly attenuated NUVrJ-selected massive star-forming galaxies at similar redshifts.

As introduced in Sect. 6.2 and shown in Fig. 13, there is generally good agreement between number densities of massive galaxies found in several simulations and in this work. This may be surprising given the range of physical recipes utilized in these simulations, and may suggest that several different tuning of physical recipes (mainly those of radiative and mechanical AGN feedback) can reproduce observations. However, from $2.5 < z \le 5.5$, we observe galaxies with mass $1.8 \times$ larger than the most massive galaxies in any of these simulations. While bias and underestimated mass uncertainties may contribute, the simulations are usually considered to be limited by their volume. Figure 15 compares the observed number densities of massive galaxies (including quiescent) to the upper limits set by volumes of different surveys as well as simulations. The fully hydrodynamical codes (EAGLE, TNG100, and SIMBA) have the smallest volumes comparable to CANDELS and are similar to the observed number densities of $M > 10^{11} M_\odot$ galaxies. However, the fact that they contain a large enough volume to catch the most massive galaxy seen in our observations may suggest that their volumes are adequate, but that their DM halo growth physics are not providing the large halos from which they can grow. Alternatively, recent observations of a highly star-forming galaxy at $z = 6.9$ (SFR ≈ 2900 $M_\odot \text{ yr}^{-1}$, Marrone et al. 2018) place close to a maximally massive DM halo not seen in current simulations – populations expected to evolve into the most massive systems at $z \approx 2-4$ (Lower et al. 2023). This is especially true of SHARK, whose volume is similar to that of COSMOS itself. More massive galaxies, if they exist, are not likely to be found in COSMOS, and so larger volumes (such as that probed by the two
these sources. Even if deeper IRAC images could be taken, the selection criteria rely on the four bands sampling the continuum shape at restframe 2–10 µm but at $z > 3$ instead sample the rest-frame stellar bulk at $\le 2$ µm. While the MIPS 24 µm data from S-COSMOS (Sanders et al. 2007) is an attractive solution, it also is too shallow (20.2 AB at $3\sigma$, Jin et al. 2018) to fully constrain the restframe 2–10 µm continuum at $3 < z \le 5$: only galaxies $H \le 20$ with large AGN fractions can be positively identified, and both low-fraction AGN at $H \approx 20$ and normal galaxies at $H < 20$ cannot be classified with MIPS. Full spectral fitting including far-infrared measurements is challenging without the constraints from channel 3, channel 4, and 24 µm, and attempts to gain further insight using STARDUST (Kokorev et al. 2021) were broadly unsuccessful with tentative contamination found to be on the order of 10–30%. We note, importantly, that removing potentially contaminated sources provides no significant change to the SMF at these epochs.

Despite these challenges, we are able to leverage the elementary AGN template fitting from LePhare to statistically assess the spectral similarity between the best-fit AGN and galaxy templates for each source, limited to rest-frame UV/optical light. While only 10% of galaxies in the total sample are best-fit with AGN templates, this fraction grows to 30% for all mass complete galaxies $3.5 < z < 5.5$ and then to 50% for those with $M > 10^{10.8} M_\odot$ with broad wings, having >80% of these sources statistically indistinguishable as either galaxy or AGN ($|\Delta\chi^2| < 0.5$). Turning to further infrared data of individual sources, we find 15% of the sample are detected by VLA-COSMOS (Smolčić et al. 2017). While not conclusive, we find 5% of these massive $M > 10^{11} M_\odot$ galaxies at $3 < z \le 5$ are detected in the ALMA maps of A$^3$COSMOS (Liu et al. 2019), which currently covers ~5% of the COSMOS field, highlighting the need for further observations. While similar sample statistics for individual sources have been extrapolated in the literature (e.g., Santini et al. 2021), anticipated surveys such as the TolTEC Ultra-Deep Galaxy Survey (Pope et al. 2019) will expand and deepen the far-infrared/(sub)mm coverage in COSMOS so that these samples can be more conclusively investigated.

A recent X-ray and radio stacking analysis of similarly selected $z > 1.5$, $M > 10^{10} M_\odot$ COSMOS2020 galaxies by Ito et al. (2022; selected from independent SED fitting with MIZUKI, Tanaka 2015) revealed that low-luminosity AGN are likely ubiquitous in massive quiescent galaxies from $1.5 < z \le 5$, even if individually they are not X-ray or radio detected. They also estimate the AGN contribution to optical/NIR continuum and find that at these redshifts, e.g., the rest-frame B-band luminosities of their quiescent galaxies are a factor of 30× larger than the expected rest-frame AGN luminosity. Ito et al. (2022) also find that the AGN luminosities for quiescent galaxies are significantly larger than those of star-forming galaxies. Together, their findings provide further confidence that our redshifts and stellar mass estimates for our X-ray undetected sample are not strongly contaminated by AGN light.

The lack of obvious systematics or likely only weak AGN contamination increases our confidence that these sources, or at least a part of them, are truly massive at $z \ge 3$. In agreement with Ito et al. (2022), we also find that >60% of $M > 10^{11} M_\odot$ galaxies appear to be star-forming, and it is likely that at least some of them are also dust obscured (DSFGs; Casey et al. 2014; Zavala et al. 2021, and references therein). Since there are a number of sources found in the deeper NIR images of COSMOS2020, it makes sense that we are now more sensitive to fainter red sources than ever before (see Fig. 14). It

is precisely this class of galaxy which are efficiently captured by infrared facilities such as Spitzer, ALMA, and JWST until now being optically “dark” (e.g., Schreiber et al. 2018b; Wang et al. 2019; Gruppioni et al. 2020; Sun et al. 2021; Fudamoto et al. 2021; Shu et al. 2022). If genuine, their existence not only points at an incompleteness particularly in the massive end of SMFs reported in previous studies lacking the necessarily deep infrared data over degree scales, but also highlights the sensitive interplay between the shape of their SEDs and the selection function consisting of the bands, their depths, and the detection methodology (see Fig. 3 of Weaver et al. 2022). Although detailed simulations will be explored in future work, qualitatively this could explain why the excess of sources relative to a Schechter at $z \approx 3-5$ diminishes at $z \approx 6$ as similarly red sources become too faint to still be detected. Optically dark galaxies selected from 2 mm ALMA observations in COSMOS from Casey et al. (2021) and Manning et al. (2022) are constrained to similar redshifts, stellar masses, and number densities. Importantly, Manning et al. (2022) studied two systems with star formation rates above $200 M_\odot \text{ yr}^{-1}$ but a gas depletion timescale <1 Gyr, suggesting rapid star formation cessation and a potential transition to massive quiescent galaxies by $z \sim 3-4$. Casey et al. (2021) find that DSFGs are responsible for ~30% of the integrated star formation rate density at $3 < z < 6$ and that 2 mm selection is an efficient way to identify larger samples in future surveys (see also Cooper et al. 2022). In an empirically based numerical model, Long et al. (2023) uses known DSFG population properties to predict the DSFG stellar mass function out to $z \sim 5$. They find a shallow power law extension beyond the knee of traditional star-forming Schechter SMFs, and posit that 60–80% of massive ($M > 10^{11} M_\odot$) star-forming galaxies at $z \approx 3-6$ are significantly dust-obscured and not captured in previous SMF studies (see also Martis et al. 2016). This is in line with our aforementioned results on significantly attenuated NUVrJ-selected massive star-forming galaxies at similar redshifts.

As introduced in Sect. 6.2 and shown in Fig. 13, there is generally good agreement between number densities of massive galaxies found in several simulations and in this work. This may be surprising given the range of physical recipes utilized in these simulations, and may suggest that several different tuning of physical recipes (mainly those of radiative and mechanical AGN feedback) can reproduce observations. However, from $2.5 < z \le 5.5$, we observe galaxies with mass $1.8 \times$ larger than the most massive galaxies in any of these simulations. While bias and underestimated mass uncertainties may contribute, the simulations are usually considered to be limited by their volume. Figure 15 compares the observed number densities of massive galaxies (including quiescent) to the upper limits set by volumes of different surveys as well as simulations. The fully hydrodynamical codes (EAGLE, TNG100, and SIMBA) have the smallest volumes comparable to CANDELS and are similar to the observed number densities of $M > 10^{11} M_\odot$ galaxies. However, the fact that they contain a large enough volume to catch the most massive galaxy seen in our observations may suggest that their volumes are adequate, but that their DM halo growth physics are not providing the large halos from which they can grow. Alternatively, recent observations of a highly star-forming galaxy at $z = 6.9$ (SFR ≈ 2900 $M_\odot \text{ yr}^{-1}$, Marrone et al. 2018) place close to a maximally massive DM halo not seen in current simulations – populations expected to evolve into the most massive systems at $z \approx 2-4$ (Lower et al. 2023). This is especially true of SHARK, whose volume is similar to that of COSMOS itself. More massive galaxies, if they exist, are not likely to be found in COSMOS, and so larger volumes (such as that probed by the two
<figure><img src="image_15.png" /><figcaption>**Fig. 15.** Comparison of observed number densities and upper limits on the rarity probed by various observed (gray steps) and simulated (colored lines) volumes. Number densities correspond to $10^{11.0} < M < 10^{11.5} M_{\odot}$ and $10^{11.5} M < 10^{12.0} M_{\odot}$ from the total sample (light and dark gray, respectively), and $10^{11.5} M < 10^{12.0} M_{\odot}$ from the quiescent sample (orange) from Fig. 8. $1\sigma$ upper limits are computed following Gehrels (1986), which for the observed volumes are dependent on widths of redshift bins.</figcaption></figure>

Euclid Deep Fields North and Fornax, as well as *Roman* will be required, see McPartland et al., in prep.).

Naturally, massive galaxies in the $z > 3$ Universe are of great interest as targets for JWST. The widest deep-field ERS program is the 100 arcmin² Cosmic Evolution Early Release Science Survey (CEERS; Finkelstein et al. 2017), and based on our estimates it stands to find approximately 22, 16, 5 for $M > 10^{10.5} M_{\odot}$, and 6, 4, 2 for $> 10^{11.0} M_{\odot}$) galaxies at $3 < z \le 3.5$, $3.5 < z \le 4.5$, and $4.5 < z \le 5.5$, respectively, although the small area will equate to a stronger density bias from cosmic variance. The widest galaxy survey field of any GO program is the 0.6 deg² COSMOS-Web (Kartaltepe et al. 2021), which is expected to find 493 (137), 340, (94), 103 (40) galaxies for the same redshift ranges (and masses). Although JWST will be instrumental in studying these sources in exquisite spatial resolution and with efficient spectroscopy (Barrufet et al. 2021; Carnall et al. 2021; Glazebrook et al. 2021), ground-based NIR observations that can efficiently identify them in wide-area surveys will retain their importance.

## 6.4. Rise of quiescent galaxies

As shown by Fig. 7, low-$z$ quiescent galaxies are well described by a two component Schechter function whose low-$M$ component diminishes rapidly from $z \approx 0.2 \rightarrow 1$ (see Fig. 11). Simultaneously our sample becomes less mass complete with redshift, doing so more rapidly for quiescent systems due to their red color characteristic of their high mass-to-light ratios. Despite the considerable uncertainties on the completeness of our low-$M$ quiescent sample outlined in Sect. 3.3, it seems likely that the shape of the quiescent stellar mass function in this work

and in previous literature in fact does turn down at low-masses, with selection effects playing a comparably minor role (see also Ilbert et al. 2010). However, there are hints of $M \lesssim 10^9 M_{\odot}$ quiescent systems at $z \ge 1.5$ (beyond our mass limit) reported by Santini et al. (2022), albeit tentative from a ~130× smaller effective area from the 33.4 arcmin² Hubble Frontier Field parallels.

Given their low apparent brightness and rarity, the present work is the first to consistently quantify the number densities of $M \approx 10^{9.5} M_{\odot}$ low-mass quiescent systems at $1.5 < z \le 4.5$ at high signal-to-noise, based on the deepest NIR data taken over such a homogeneously measured degree-scale area required to find them in sufficient numbers. As seen in Fig. 8, the rate of growth in the number density of low-mass ($M \lesssim 10^{10} M_{\odot}$) quiescent galaxies has been seemingly rapid over the past ~10 billion years. Still, they constitute only a minor fraction of low-mass sources by $z \sim 0.2$, and by extrapolating the number densities from $z \approx 1-2$ one may expect to find none within COSMOS by $z \sim 3$ (see Fig. 8). This is typically interpreted by the phenomenological model of Peng et al. (2010) to mean that the processes which act to halt star formation cessation in low-$M$ systems are inefficient at early times. For example the lack of virialized at $z > 2-3$ structure makes influence from environmental effects unlikely.

As shown by Fig. 8, the apparent lack of growth in the abundance of massive systems at $z < 1$ is the result of a decline in star-forming galaxies simultaneous with an increase in quiescent number densities. As noted by Ilbert et al. (2013), this decrease in the star-forming population is consistent with star formation cessation becoming extremely efficient, to the extent that massive star-forming galaxies are becoming quiescent faster than they can be replaced. Therefore, the mass assembly of massive star-forming systems at $z < 1$ is slower than the cessation of their star formation. However, there is also a slowing down in the rate of growth in the number density of massive quiescent systems themselves. While this may suggest high incidences of dry mergers or even rejuvination, it must also relate to evolving demographics: from $z \approx 1.5 \rightarrow 0.3$ there are fewer and fewer high-$M$ star-forming galaxies available to become quiescent. While the number density growth of massive systems seems to have stalled out, that of lower-mass systems continue to grow; this is the so-called phenomenon of “downsizing with time” (Cowie et al. 1996; Neistein et al. 2006; Fontanot et al. 2009).

The quiescent mass function at $M > 10^{11} M_{\odot}$ does not change much from $z \approx 4.5 \rightarrow 2.0$, with surprisingly elevated quiescent fractions (Fig. 9) being above 20% at $z < 5$. Figure 8 shows why: while the number density of star-forming galaxies at $z \approx 5$ is lower than at $z \approx 2$, quiescent galaxies are roughly constant in density over this same range. However, we caution that their number densities are only marginally above what should be possible to find within the nominal volume of COSMOS at these redshifts, and so it is plausible that they are dominated by noise and/or photo-$z$ bias. The question of whether or not this behavior is genuine can only be explored in future large volume surveys, as demonstrated by Fig. 15. Even though the most massive galaxies at $z \ge 3-4$ are typically too faint to obtain continuum features, spectroscopic follow-up and supporting SED fitting will continue to provide valuable insights (Gobat et al. 2012; Schreiber et al. 2018a; Valentino et al. 2020; Glazebrook et al. 2017).

Figure 16 shows the fraction of quiescent galaxies in bins of mass for three epochs: $z \approx 0.3, 1.3$, and $2.7$. A key advantage of examining fractions is that they are less sensitive to the overall normalization of the simulation and biases from observations (e.g., in simulations: overproduction of all galaxy masses; in
<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Comparison of observed number densities and upper limits on the rarity probed by various observed (gray steps) and simulated (colored lines) volumes. Number densities correspond to $10^{11.0} < M < 10^{11.5} M_{\odot}$ and $10^{11.5} M < 10^{12.0} M_{\odot}$ from the total sample (light and dark gray, respectively), and $10^{11.5} M < 10^{12.0} M_{\odot}$ from the quiescent sample (orange) from Fig. 8. $1\sigma$ upper limits are computed following Gehrels (1986), which for the observed volumes are dependent on widths of redshift bins.</figcaption></figure>

*Euclid Deep Fields North and Fornax, as well as Roman will be required, see McPartland et al., in prep.).*

Naturally, massive galaxies in the $z > 3$ Universe are of great interest as targets for JWST. The widest deep-field ERS program is the 100 arcmin² Cosmic Evolution Early Release Science Survey (CEERS; Finkelstein et al. 2017), and based on our estimates it stands to find approximately 22, 16, 5 for $M > 10^{10.5} M_{\odot}$, and 6, 4, 2 for $> 10^{11.0} M_{\odot}$) galaxies at $3 < z \le 3.5$, $3.5 < z \le 4.5$, and $4.5 < z \le 5.5$, respectively, although the small area will equate to a stronger density bias from cosmic variance. The widest galaxy survey field of any GO program is the 0.6 deg² COSMOS-Web (Kartaltepe et al. 2021), which is expected to find 493 (137), 340, (94), 103 (40) galaxies for the same redshift ranges (and masses). Although JWST will be instrumental in studying these sources in exquisite spatial resolution and with efficient spectroscopy (Barrufet et al. 2021; Carnall et al. 2021; Glazebrook et al. 2021), ground-based NIR observations that can efficiently identify them in wide-area surveys will retain their importance.

## 6.4. Rise of quiescent galaxies

As shown by Fig. 7, low-$z$ quiescent galaxies are well described by a two component Schechter function whose low-$M$ component diminishes rapidly from $z \approx 0.2 \rightarrow 1$ (see Fig. 11). Simultaneously our sample becomes less mass complete with redshift, doing so more rapidly for quiescent systems due to their red color characteristic of their high mass-to-light ratios. Despite the considerable uncertainties on the completeness of our low-$M$ quiescent sample outlined in Sect. 3.3, it seems likely that the shape of the quiescent stellar mass function in this work

and in previous literature in fact does turn down at low-masses, with selection effects playing a comparably minor role (see also Ilbert et al. 2010). However, there are hints of $M \lesssim 10^9 M_{\odot}$ quiescent systems at $z \ge 1.5$ (beyond our mass limit) reported by Santini et al. (2022), albeit tentative from a ~130× smaller effective area from the 33.4 arcmin² Hubble Frontier Field parallels.

Given their low apparent brightness and rarity, the present work is the first to consistently quantify the number densities of $M \approx 10^{9.5} M_{\odot}$ low-mass quiescent systems at $1.5 < z \le 4.5$ at high signal-to-noise, based on the deepest NIR data taken over such a homogeneously measured degree-scale area required to find them in sufficient numbers. As seen in Fig. 8, the rate of growth in the number density of low-mass ($M \lesssim 10^{10} M_{\odot}$) quiescent galaxies has been seemingly rapid over the past ~10 billion years. Still, they constitute only a minor fraction of low-mass sources by $z \sim 0.2$, and by extrapolating the number densities from $z \approx 1-2$ one may expect to find none within COSMOS by $z \sim 3$ (see Fig. 8). This is typically interpreted by the phenomenological model of Peng et al. (2010) to mean that the processes which act to halt star formation cessation in low-$M$ systems are inefficient at early times. For example the lack of virialized at $z > 2-3$ structure makes influence from environmental effects unlikely.

As shown by Fig. 8, the apparent lack of growth in the abundance of massive systems at $z < 1$ is the result of a decline in star-forming galaxies simultaneous with an increase in quiescent number densities. As noted by Ilbert et al. (2013), this decrease in the star-forming population is consistent with star formation cessation becoming extremely efficient, to the extent that massive star-forming galaxies are becoming quiescent faster than they can be replaced. Therefore, the mass assembly of massive star-forming systems at $z < 1$ is slower than the cessation of their star formation. However, there is also a slowing down in the rate of growth in the number density of massive quiescent systems themselves. While this may suggest high incidences of dry mergers or even rejuvination, it must also relate to evolving demographics: from $z \approx 1.5 \rightarrow 0.3$ there are fewer and fewer high-$M$ star-forming galaxies available to become quiescent. While the number density growth of massive systems seems to have stalled out, that of lower-mass systems continue to grow; this is the so-called phenomenon of “downsizing with time” (Cowie et al. 1996; Neistein et al. 2006; Fontanot et al. 2009).

The quiescent mass function at $M > 10^{11} M_{\odot}$ does not change much from $z \approx 4.5 \rightarrow 2.0$, with surprisingly elevated quiescent fractions (Fig. 9) being above 20% at $z < 5$. Figure 8 shows why: while the number density of star-forming galaxies at $z \approx 5$ is lower than at $z \approx 2$, quiescent galaxies are roughly constant in density over this same range. However, we caution that their number densities are only marginally above what should be possible to find within the nominal volume of COSMOS at these redshifts, and so it is plausible that they are dominated by noise and/or photo-$z$ bias. The question of whether or not this behavior is genuine can only be explored in future large volume surveys, as demonstrated by Fig. 15. Even though the most massive galaxies at $z \ge 3-4$ are typically too faint to obtain continuum features, spectroscopic follow-up and supporting SED fitting will continue to provide valuable insights (Gobat et al. 2012; Schreiber et al. 2018a; Valentino et al. 2020; Glazebrook et al. 2017).

Figure 16 shows the fraction of quiescent galaxies in bins of mass for three epochs: $z \approx 0.3, 1.3$, and $2.7$. A key advantage of examining fractions is that they are less sensitive to the overall normalization of the simulation and biases from observations (e.g., in simulations: overproduction of all galaxy masses; in
<figure><img src="image_16.png" /><figcaption>**Fig. 16.** Quiescent mass fractions. Upper: fraction of quiescent galaxies as a function of mass for three redshift ranges compared with predictions from EAGLE (Furlong et al. 2015) and SHARK (Lagos et al. 2018) as well as measurements from the empirically calibrated UNIVERSE MACHINE (Behroozi et al. 2019). Lower: fractional difference between this work and literature predictions/measurements. See Appendix B for more details.</figcaption></figure>

observations: systematics in effective survey volume), and provide additional insight which is obscured by simply comparing mass functions. Although comparisons to UNIVERSE MACHINE have been already discussed in Sect. 6.2, we introduce two samples of quiescent galaxies which we selected from EAGLE and SHARK with an *NUVrJ* selection consistent with our methodology. While EAGLE underproduces quiescent systems by 10–20% at all masses, SHARK overproduces them in all but the most massive bins at *z* ≈ 0.3–1.3 but underproduces them at *z* ≈ 2.7. Roughly, the fraction of quiescent galaxies in SHARK at *z* ≈ 1.3 matches the observed fractions at *z* ≈ 0.3. This may suggest that SHARK’s physical recipes that halt star formation in lower mass galaxies are too aggressive. These same systems are also seemingly overproduced (see Fig. 13), and so may be assembling and maturing too early. For additional figures and details, see Appendix B.

## 6.5. Dark matter halo connection

The mass assembly of a galaxy is inherently connected to the dark matter halo in which it formed and grew (see Wechsler & Tinker 2018 for a review). Yet, stellar masses *M* have been observed to be <2% of their halo mass *M*<sub>h</sub>, which point to galaxy formation as a strikingly inefficient process (Mandelbaum et al. 2006; Conroy & Wechsler 2009;

Behroozi et al. 2010). This has led to investigations into the role of dark matter halo mass in influencing star formation cessation, including promoting thermal heating and/or gas expulsion by AGN (Shankar et al. 2006; Main et al. 2017), as well as virial shock heating of in-falling cold gas whose Jeans mass inhibits the formation of star-forming molecular clouds (so-called hot-halo mode, Birnboim & Dekel 2003). Hence, the gravitational influence of the halo mass on the cold gas reservoir regulates the ability of a galaxy to form stars, and hence the stellar-to-halo mass relationship (SHMR). Therefore it is no surprise that there is a similar evolution between the specific mass increase rate of the halo by accretion (sMIR ≡ $M_h^{-1}∂M_h/∂t$, Neistein & Dekel 2008; Saintonge et al. 2013) and the instantaneous mass growth by star formation (sSFR; see discussions in Lilly et al. 2013). While in this work we restrict ourselves to comparisons with theoretical HMFs, another paper in this series computes a self consistent SHMR based on measuring the halo occupation distribution directly from angular correlation functions and SMFs of COSMOS2020 galaxies (Shuntov et al. 2022). For an investigation into the SHMR split by star-forming and quiescent samples, see Cowley et al. (2019), which is also based on galaxies from the COSMOS field.

As shown in Fig. 17, we compare our observed and inferred SMFs to the halo mass function (HMF) of Tinker et al. (2008)<fnref n="13" /> from *z* = 1.5 → 7.5. We choose not to show *z* < 1.5 as these comparisons have been thoroughly explored by previous investigations (e.g., Davidzon et al. 2017; Legrand et al. 2019) and we do not observe any significant differences. The effects of feedback can be seen in the first panel of Fig. 17 at 1.5 < *z* ≤ 2.0 that explains the relatively lower number densities of both low- and high-*M* systems, with those around *M*<sup>*</sup> being most similar to HMF. This apparent tension is a well known feature and lies at the foundation of the contemporary galaxy evolution paradigm, involving halting star formation by secular (inter nal) and/or environmental (external) action on the gas reservoir such as thermal heating, dynamical turbulence, and/or removal. The growth of low-mass *M* < *M*<sup>*</sup> galaxies can be impeded by secular (e.g., supernovae and stellar winds) and environmental processes (e.g., ram-pressure stripping, thermal evaporation). Similarly, secular and environmental processes can also impede the growth of massive *M* > *M*<sup>*</sup> galaxies, albeit from different driving forces such as AGN and major mergers, respectively (see Peng et al. 2010, 2012, 2015; Wechsler & Tinker 2018; Förster Schreiber & Wuyts 2020). In this context, the characteristic knee at *M*<sup>*</sup> is the result of a build-up of massive galaxies which can no longer sustain mass growth. At *z* ≈ 3, the low-mass end is still considerably lower than the scaled HMF but the number density of massive systems comes into agreement. Although the stellar mass function slightly lies above the SHMR-scaled HMF at some points, we caution that this should not be taken as a challenge to theory as it assumes that the SHMR at *z* = 0 is appropriate at higher-*z* (which is unlikely; see Legrand et al. 2019) and small modifications can reconcile this difference. We note that Fig. 18 of Davidzon et al. (2017) finds no such offset

<footnote n="13">¹³ The Tinker et al. (2008) HMF is computed according to our cosmology ($\sigma_8$ = 0.82) at the mid-point in each z-bin using COLUSSUS (Diemer 2018), which explicitly takes into account that these mass functions derived were originally derived from spherical overdensities, which is not universal with redshift (see Eqs. (3)–(8) of Tinker et al. 2008). We adopt the definition for a halo as the DM mass contained in a region 200× the mean matter density, and found no significant differences from using other definitions (e.g., friends-of-friends, spherical overdensity, virial radius).</footnote>
<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Quiescent mass fractions. Upper: fraction of quiescent galaxies as a function of mass for three redshift ranges compared with predictions from EAGLE (Furlong et al. 2015) and SHARK (Lagos et al. 2018) as well as measurements from the empirically calibrated UNIVERSE MACHINE (Behroozi et al. 2019). Lower: fractional difference between this work and literature predictions/measurements. See Appendix B for more details.</figcaption></figure>

observations: systematics in effective survey volume), and provide additional insight which is obscured by simply comparing mass functions. Although comparisons to UNIVERSE MACHINE have been already discussed in Sect. 6.2, we introduce two samples of quiescent galaxies which we selected from EAGLE and SHARK with an *NUVrJ* selection consistent with our methodology. While EAGLE underproduces quiescent systems by 10–20% at all masses, SHARK overproduces them in all but the most massive bins at *z* ≈ 0.3–1.3 but underproduces them at *z* ≈ 2.7. Roughly, the fraction of quiescent galaxies in SHARK at *z* ≈ 1.3 matches the observed fractions at *z* ≈ 0.3. This may suggest that SHARK’s physical recipes that halt star formation in lower mass galaxies are too aggressive. These same systems are also seemingly overproduced (see Fig. 13), and so may be assembling and maturing too early. For additional figures and details, see Appendix B.

## 6.5. Dark matter halo connection

The mass assembly of a galaxy is inherently connected to the dark matter halo in which it formed and grew (see Wechsler & Tinker 2018 for a review). Yet, stellar masses *M* have been observed to be <2% of their halo mass *M*<sub>h</sub>, which point to galaxy formation as a strikingly inefficient process (Mandelbaum et al. 2006; Conroy & Wechsler 2009;

Behroozi et al. 2010). This has led to investigations into the role of dark matter halo mass in influencing star formation cessation, including promoting thermal heating and/or gas expulsion by AGN (Shankar et al. 2006; Main et al. 2017), as well as virial shock heating of in-falling cold gas whose Jeans mass inhibits the formation of star-forming molecular clouds (so-called hot-halo mode, Birnboim & Dekel 2003). Hence, the gravitational influence of the halo mass on the cold gas reservoir regulates the ability of a galaxy to form stars, and hence the stellar-to-halo mass relationship (SHMR). Therefore it is no surprise that there is a similar evolution between the specific mass increase rate of the halo by accretion (sMIR ≡ $M_h^{-1}∂M_h/∂t$, Neistein & Dekel 2008; Saintonge et al. 2013) and the instantaneous mass growth by star formation (sSFR; see discussions in Lilly et al. 2013). While in this work we restrict ourselves to comparisons with theoretical HMFs, another paper in this series computes a self consistent SHMR based on measuring the halo occupation distribution directly from angular correlation functions and SMFs of COSMOS2020 galaxies (Shuntov et al. 2022). For an investigation into the SHMR split by star-forming and quiescent samples, see Cowley et al. (2019), which is also based on galaxies from the COSMOS field.

As shown in Fig. 17, we compare our observed and inferred SMFs to the halo mass function (HMF) of Tinker et al. (2008)<fnref n="13" /> from *z* = 1.5 → 7.5. We choose not to show *z* < 1.5 as these comparisons have been thoroughly explored by previous investigations (e.g., Davidzon et al. 2017; Legrand et al. 2019) and we do not observe any significant differences. The effects of feedback can be seen in the first panel of Fig. 17 at 1.5 < *z* ≤ 2.0 that explains the relatively lower number densities of both low- and high-*M* systems, with those around *M*<sup>*</sup> being most similar to HMF. This apparent tension is a well known feature and lies at the foundation of the contemporary galaxy evolution paradigm, involving halting star formation by secular (inter nal) and/or environmental (external) action on the gas reservoir such as thermal heating, dynamical turbulence, and/or removal. The growth of low-mass *M* < *M*<sup>*</sup> galaxies can be impeded by secular (e.g., supernovae and stellar winds) and environmental processes (e.g., ram-pressure stripping, thermal evaporation). Similarly, secular and environmental processes can also impede the growth of massive *M* > *M*<sup>*</sup> galaxies, albeit from different driving forces such as AGN and major mergers, respectively (see Peng et al. 2010, 2012, 2015; Wechsler & Tinker 2018; Förster Schreiber & Wuyts 2020). In this context, the characteristic knee at *M*<sup>*</sup> is the result of a build-up of massive galaxies which can no longer sustain mass growth. At *z* ≈ 3, the low-mass end is still considerably lower than the scaled HMF but the number density of massive systems comes into agreement. Although the stellar mass function slightly lies above the SHMR-scaled HMF at some points, we caution that this should not be taken as a challenge to theory as it assumes that the SHMR at *z* = 0 is appropriate at higher-*z* (which is unlikely; see Legrand et al. 2019) and small modifications can reconcile this difference. We note that Fig. 18 of Davidzon et al. (2017) finds no such offset

<footnote n="13">¹³ The Tinker et al. (2008) HMF is computed according to our cosmology ($\sigma_8$ = 0.82) at the mid-point in each z-bin using COLUSSUS (Diemer 2018), which explicitly takes into account that these mass functions derived were originally derived from spherical overdensities, which is not universal with redshift (see Eqs. (3)–(8) of Tinker et al. 2008). We adopt the definition for a halo as the DM mass contained in a region 200× the mean matter density, and found no significant differences from using other definitions (e.g., friends-of-friends, spherical overdensity, virial radius).</footnote>
<figure><img src="image_17.png" /><figcaption>**Fig. 17.** Comparison of the inferred $z > 1.5$ galaxy stellar mass function (colored curve with 1 and $2\sigma$ envelopes) fitted to observed measurements (gray points). We scale the Tinker et al. (2008) halo mass function (HMF; lower bound of the gray shaded region) by 0.018 corresponding to the stellar-to-halo mass relation at $z=0$ and $M_h = M_h^*$ (Behroozi et al. 2013) to produce an idealized SMF (gray curve). Upper limits on the SMF at each redshift are derived from the HMF assuming a fixed baryon fraction (0.166). Upper limits for empty bins are shown by the horizontal gray line with an arrow. Mass incomplete measurements are not shown.</figcaption></figure>

using the same scaling, but we are unable to reproduce their precise result.

We derive an upper limit for the baryonic matter distribution by rescaling the HMF by $\Omega_b/\Omega_m = 0.166$, which for our adopted cosmology is redshift independent, and assume a 100% efficient baryon-to-stellar mass conversion. This is the maximum SMF physically allowed under our simple assumptions. This upper limit becomes relevant especially at $z > 3.5$ where our observed number densities exceed those inferred by the Schechter model. While a large Eddington bias or selection systematic could explain this excess (see Sect. 6.3), we stress that our inferred SMF assumes the applicability of Schechter formalism, which cannot accurately describe the observed number densities at $3.5 < z \le 5.5$. Nevertheless, the inferred stellar mass function agrees well with the SHMR-scaled HMF up to $z \approx 7$. This suggests that the most efficient haloes during these early epochs are not around $M^*$, but rather the most massive ones, and with little observed evolution consistent with the findings from Stefanon et al. (2021). This phenomenon has also been observed in some local massive spiral galaxies (e.g., Di Teodoro et al. 2023), which have seemingly matured without significant star formation cessation. One interpretation is that this high star formation efficiency in massive haloes is the result of diminished feedback from AGN, with stellar mass growing similarly to the host halo at early times. Indeed, this is consistent with findings of inefficient radiative AGN feedback from simulations (Kaviraj et al. 2017; Laigle et al. 2019; Roos et al. 2015; Bieri et al. 2017; Habouzit et al. 2022), as well as FIR/radio observations of AGN activity at $z > 3$ (Maiolino et al. 2012; Cicone et al. 2014, 2015; Padovani et al. 2015; Vito et al. 2018).

At no point do our mass functions, observed or inferred, exceed this upper limit. Therefore we do not report evidence of “impossibly early galaxies” introduced by Steinhardt et al. (2016) who point out that there appears to be too many massive

galaxies at $z > 4$ compared to the dark matter haloes that should host them. However at $z > 6.5$, where Steinhardt et al. (2016) predicts that the effect will be most obvious, we report observed number densities approaching this upper limit and in clear excess of the SHMR-scaled HMF. We caution that these sources are the most vulnerable to misclassification and bias, being susceptible to blending in addition to being constrained by only a handful of NIR bands yielding proportionately uncertain Schechter fits. Extrapolation to $M \ge 10^{11.5} M_\odot$ would place their number density below that which can be probed in a volume contained by the $0.716 \text{ deg}^2$ area of the Ultra-Deep region, and so COSMOS2020 is unlikely to find them if they exist. While they are not “impossibly early galaxies”, their surprising abundance hint that explorations of $z > 7$ with future deep, large-volume surveys may provide the evidence necessary to firmly challenge theoretical frameworks.

## 7. Summary and conclusions

Following on from COSMOS2020 photo-$z$ catalog (Weaver et al. 2022), we study the shape and evolution of the galaxy stellar mass function. Our measurements span three decades in mass (at $z \approx 0.2$) across 10 billion years of cosmic history ($z = 7.5 \rightarrow 0.2$), including the most mass complete sample of quiescent galaxies at $z > 2$ enabled by our unprecedented, homogeneous NIR depths across an effective 1.27 deg$^2$. We probe a volume nearly 2x that of Davidzon et al. (2017) which not only improves sample statistics but also finds new galaxies of still greater mass at all redshifts. Complementary deep IRAC coverage has allowed us to directly measure the stellar bulk, and hence galaxy stellar mass, in a less biased way and to higher redshifts compared to $K_s$-based measurements. We developed a robust, mass-dependent error budget with contributions from poisson, stellar mass, and cosmic variance, and account for the effects of
<figure><img src="image_18.png" /><figcaption>**Fig. 17.** Comparison of the inferred $z > 1.5$ galaxy stellar mass function (colored curve with 1 and $2\sigma$ envelopes) fitted to observed measurements (gray points). We scale the Tinker et al. (2008) halo mass function (HMF; lower bound of the gray shaded region) by 0.018 corresponding to the stellar-to-halo mass relation at $z=0$ and $M_h = M_h^*$ (Behroozi et al. 2013) to produce an idealized SMF (gray curve). Upper limits on the SMF at each redshift are derived from the HMF assuming a fixed baryon fraction (0.166). Upper limits for empty bins are shown by the horizontal gray line with an arrow. Mass incomplete measurements are not shown.</figcaption></figure>

using the same scaling, but we are unable to reproduce their precise result.

We derive an upper limit for the baryonic matter distribution by rescaling the HMF by $\Omega_b/\Omega_m = 0.166$, which for our adopted cosmology is redshift independent, and assume a 100% efficient baryon-to-stellar mass conversion. This is the maximum SMF physically allowed under our simple assumptions. This upper limit becomes relevant especially at $z > 3.5$ where our observed number densities exceed those inferred by the Schechter model. While a large Eddington bias or selection systematic could explain this excess (see Sect. 6.3), we stress that our inferred SMF assumes the applicability of Schechter formalism, which cannot accurately describe the observed number densities at $3.5 < z \le 5.5$. Nevertheless, the inferred stellar mass function agrees well with the SHMR-scaled HMF up to $z \approx 7$. This suggests that the most efficient haloes during these early epochs are not around $M^*$, but rather the most massive ones, and with little observed evolution consistent with the findings from Stefanon et al. (2021). This phenomenon has also been observed in some local massive spiral galaxies (e.g., Di Teodoro et al. 2023), which have seemingly matured without significant star formation cessation. One interpretation is that this high star formation efficiency in massive haloes is the result of diminished feedback from AGN, with stellar mass growing similarly to the host halo at early times. Indeed, this is consistent with findings of inefficient radiative AGN feedback from simulations (Kaviraj et al. 2017; Laigle et al. 2019; Roos et al. 2015; Bieri et al. 2017; Habouzit et al. 2022), as well as FIR/radio observations of AGN activity at $z > 3$ (Maiolino et al. 2012; Cicone et al. 2014, 2015; Padovani et al. 2015; Vito et al. 2018).

At no point do our mass functions, observed or inferred, exceed this upper limit. Therefore we do not report evidence of “impossibly early galaxies” introduced by Steinhardt et al. (2016) who point out that there appears to be too many massive

galaxies at $z > 4$ compared to the dark matter haloes that should host them. However at $z > 6.5$, where Steinhardt et al. (2016) predicts that the effect will be most obvious, we report observed number densities approaching this upper limit and in clear excess of the SHMR-scaled HMF. We caution that these sources are the most vulnerable to misclassification and bias, being susceptible to blending in addition to being constrained by only a handful of NIR bands yielding proportionately uncertain Schechter fits. Extrapolation to $M \ge 10^{11.5} M_\odot$ would place their number density below that which can be probed in a volume contained by the $0.716 \text{ deg}^2$ area of the Ultra-Deep region, and so COSMOS2020 is unlikely to find them if they exist. While they are not “impossibly early galaxies”, their surprising abundance hint that explorations of $z > 7$ with future deep, large-volume surveys may provide the evidence necessary to firmly challenge theoretical frameworks.

## 7. Summary and conclusions

Following on from COSMOS2020 photo-z catalog (Weaver et al. 2022), we study the shape and evolution of the galaxy stellar mass function. Our measurements span three decades in mass (at $z \approx 0.2$) across 10 billion years of cosmic history ($z = 7.5 \rightarrow 0.2$), including the most mass complete sample of quiescent galaxies at $z > 2$ enabled by our unprecedented, homogeneous NIR depths across an effective 1.27 deg$^2$. We probe a volume nearly 2x that of Davidzon et al. (2017) which not only improves sample statistics but also finds new galaxies of still greater mass at all redshifts. Complementary deep IRAC coverage has allowed us to directly measure the stellar bulk, and hence galaxy stellar mass, in a less biased way and to higher redshifts compared to $K_s$-based measurements. We developed a robust, mass-dependent error budget with contributions from poisson, stellar mass, and cosmic variance, and account for the effects of
Eddington bias by fitting a kernel-convolved Schechter function to our observed SMF. We use three fitting techniques, including the continuity model of Leja et al. (2019a), finding good agreement with literature measurements with smaller bin-to-bin variance with $z$. We stress that our derived parameters are dependent on the assumed Eddington correction, and while the inferred SMF evaluated at maximum likelihood and associated parameters are robust, parameters (and their uncertainties) derived from the median of posterior distributions can be unreliable if constraints are weak. To make these results more transparent and to facilitate comparisons, we have made the object IDs and key measurements presented in this work available to download<fnref n="14" />.

Although literature comparisons on the shape of the SMF at fixed redshift show good agreement, the novel advantage of COSMOS2020 is the extended historical baseline over which the mass functions (as well as many other properties) can be consistently measured. Not only do we examine the evolution of the integrated mass density $\rho_*$ over this time, we also examine the remarkably consistent rate of growth in the number densities of systems of differen masses from $z \approx 7 \rightarrow 1$, whereupon the most massive star-forming galaxies become quiescent faster than they can be replaced. Similarly, we find evidence for the sharp rise in low-mass quiescent systems consistent with the phenomenological model of Peng et al. (2010) probed to $z \approx 2.5$ where our sample becomes incomplete. Although tentative, we also find evidence for a sustained >20% fraction of high-mass quiescent systems from $z \approx 5 \rightarrow 2$.

Furthermore, we highlight three main results:

- Comparisons with several hydrodynamical simulations and semi-analytical models indicate an exceptional degree of agreement for the most massive galaxies out to high-z. This comes despite the surprisingly high number densities of massive galaxies at $z \approx 3-5$ in excess of a Schechter function, suggesting that existing physical recipes are assembling massive $M \approx 10^{10-11.5} M_\odot$ systems in sufficient quantity at early times. In order to explore star formation cessation and feedback modes, we identify quiescent galaxies out to $z \approx 5.5$ by means of a *NUVrJ* selection and compare them to consistently selected quiescent samples produced by two simulation codes, finding evidence for delayed assembly of low-mass quiescent systems in *EAGLE*, and too rapid assembly in the *SHARK*.

- A closer examination of these massive systems reveals that a quarter are not found in COSMOS2015. Not only are they $K_s$-faint, but their extremely red colors challenge SED fitting templates. We find no strong evidence for AGN contamination, although we stress the need for future infrared facilities with deep surveys capable of measuring the rest-frame MIR light at $z \approx 3-5$. Recent findings of optically dark galaxies from IRAC and ALMA suggest that previous studies have missed contributions from dust-obscured star-forming galaxies. Their brightness, redshifts, mass, and number densities are consistent with our findings, suggesting that the $K_s \sim 26$ depth of UltraVISTA DR4 may indeed be sufficient to reach out into the tail end of this population missed by previous optical-NIR selections. Further work is required to conclusively establish the nature of these massive galaxy candidates.

- Lastly, we investigate the connection to dark matter halos by comparing both our observed and inferred SMFs to constraints provided by the HMF. While we confirm the divergence of the SMF from the HMF at both low- and high-mass

regimes which has been historically interpreted as evidence for feedback processes, the massive end of the inferred, Schechter-fit SMF comes into agreement with the HMF at $z \ge 2$. While we find no evidence of tension which would challenge theoretical models, our observed number densities at $z \approx 3-5$ approach the upper limit for fully efficient star formation in the most massive halos. Larger volume surveys containing even more massive systems, if they exist, may be able to challenge these models, especially at $z \ge 6-7$.

The launch of JWST has opened the door on a new era, and it will soon be flanked by efficient survey facilities from space (Euclid, Roman) and the ground (Rubin). While massive quiescent systems may exist at $z \sim 5$ and perhaps even at earlier times (Mawatari et al. 2016, 2020), their identification is beyond the reach of COSMOS2020. They may be identified soon by deep degree-scale JWST surveys i.e., COSMOS-Web (Kartaltepe et al. 2021; Casey et al. 2023) and possibly by narrower ones including the JWST Advanced Deep Extragalactic Survey (JADES, Eisenstein et al. 2017; Ferruit 2017), the Cosmic Evolution Early Release Science Survey (CEERS, Finkelstein et al. 2017), the Next Generation Deep Extragalactic Exploratory Public Survey (NGDEEP, Finkelstein et al. 2021), and Ultra-deep NIRCam and NIRSpec Observations Before the Epoch of Reionization (UNCOVER, Labbé et al. 2021; Bezanson et al. 2022) to name a few.

While low mass quiescent systems at $z \ge 2.5$ may be found by JWST (with hints emerging from Marchesini et al. 2023), a truly precise quantification of their number density and demographics is likely to remain a future objective. By extrapolating these observations, <1% of $M \approx 10^{9.5-10.0} M_\odot$ are expected to be quiescent by $z \sim 2.5$, become even rarer at earlier times. While identifying even one quiescent, low-mass system at $z > 2$ in the absence of virialized structure would present a significant challenge to the paradigm of Peng et al. (2010), performing a statistically meaningful survey of them will require incredibly deep, degree-scale NIR surveys, placing them out of reach by current facilities. While the deepest degree-scale surveys from JWST (COSMOS-Web), *Roman*, and *Euclid* stand to establish the rarity of these systems at $z \approx 2-3$, it seems that no currently planned survey will be able to definitively quantify their contribution at $z \ge 3$.

While explorations of mass-selected samples at $z > 7$ are being made possible by JWST, the most UV-luminous sources at these epochs will likely be missed by small area programs and yet are crucial to a comprehensive study of the ionizing UV budget (Kauffmann et al. 2022; Donnan et al. 2023). Following up known samples with JWST will not only allow us to measure the star formation rate and dust content of the first ultra-luminous galaxies from deep within the epoch of reionization, but also to directly identify the progenitors of $z \sim 3-4$ massive galaxies while still in their formation stage (e.g., Weaver et al. 2021). Yet despite the incredible promises of highly resolved IR imaging and spectroscopy from JWST, identifying the rarest and potentially most informative populations that can challenge and thereby improve galaxy formation models will remain the domain of wide-field surveys.

*Acknowledgements.* The authors thank Vadim Ruskov, Sidney Lower, Desika Narayanan, Stephen Wilkins, William Roper, Lukas Furtak, and Pratika Dayal for helpful discussions. We are also grateful for the many helpful and constructive comments from the anonymous referee. The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation under grant No. 140. S.T. and J.W. acknowledge support from the European Research Council (ERC) Consolidator Grant funding scheme (project ConTExt, grant

<footnote n="14">¹⁴ https://doi.org/10.5281/zenodo.7808832</footnote>
Eddington bias by fitting a kernel-convolved Schechter function to our observed SMF. We use three fitting techniques, including the continuity model of Leja et al. (2019a), finding good agreement with literature measurements with smaller bin-to-bin variance with $z$. We stress that our derived parameters are dependent on the assumed Eddington correction, and while the inferred SMF evaluated at maximum likelihood and associated parameters are robust, parameters (and their uncertainties) derived from the median of posterior distributions can be unreliable if constraints are weak. To make these results more transparent and to facilitate comparisons, we have made the object IDs and key measurements presented in this work available to download<fnref n="14" />.

Although literature comparisons on the shape of the SMF at fixed redshift show good agreement, the novel advantage of COSMOS2020 is the extended historical baseline over which the mass functions (as well as many other properties) can be consistently measured. Not only do we examine the evolution of the integrated mass density $\rho_*$ over this time, we also examine the remarkably consistent rate of growth in the number densities of systems of differen masses from $z \approx 7 \rightarrow 1$, whereupon the most massive star-forming galaxies become quiescent faster than they can be replaced. Similarly, we find evidence for the sharp rise in low-mass quiescent systems consistent with the phenomenological model of Peng et al. (2010) probed to $z \approx 2.5$ where our sample becomes incomplete. Although tentative, we also find evidence for a sustained >20% fraction of high-mass quiescent systems from $z \approx 5 \rightarrow 2$.

Furthermore, we highlight three main results:

- Comparisons with several hydrodynamical simulations and semi-analytical models indicate an exceptional degree of agreement for the most massive galaxies out to high-z. This comes despite the surprisingly high number densities of massive galaxies at $z \approx 3-5$ in excess of a Schechter function, suggesting that existing physical recipes are assembling massive $M \approx 10^{10-11.5} M_\odot$ systems in sufficient quantity at early times. In order to explore star formation cessation and feedback modes, we identify quiescent galaxies out to $z \approx 5.5$ by means of a *NUVrJ* selection and compare them to consistently selected quiescent samples produced by two simulation codes, finding evidence for delayed assembly of low-mass quiescent systems in *EAGLE*, and too rapid assembly in the *SHARK*.

- A closer examination of these massive systems reveals that a quarter are not found in COSMOS2015. Not only are they $K_s$-faint, but their extremely red colors challenge SED fitting templates. We find no strong evidence for AGN contamination, although we stress the need for future infrared facilities with deep surveys capable of measuring the rest-frame MIR light at $z \approx 3-5$. Recent findings of optically dark galaxies from IRAC and ALMA suggest that previous studies have missed contributions from dust-obscured star-forming galaxies. Their brightness, redshifts, mass, and number densities are consistent with our findings, suggesting that the $K_s \sim 26$ depth of UltraVISTA DR4 may indeed be sufficient to reach out into the tail end of this population missed by previous optical-NIR selections. Further work is required to conclusively establish the nature of these massive galaxy candidates.

- Lastly, we investigate the connection to dark matter halos by comparing both our observed and inferred SMFs to constraints provided by the HMF. While we confirm the divergence of the SMF from the HMF at both low- and high-mass

regimes which has been historically interpreted as evidence for feedback processes, the massive end of the inferred, Schechter-fit SMF comes into agreement with the HMF at $z \ge 2$. While we find no evidence of tension which would challenge theoretical models, our observed number densities at $z \approx 3-5$ approach the upper limit for fully efficient star formation in the most massive halos. Larger volume surveys containing even more massive systems, if they exist, may be able to challenge these models, especially at $z \ge 6-7$.

The launch of JWST has opened the door on a new era, and it will soon be flanked by efficient survey facilities from space (Euclid, Roman) and the ground (Rubin). While massive quiescent systems may exist at $z \sim 5$ and perhaps even at earlier times (Mawatari et al. 2016, 2020), their identification is beyond the reach of COSMOS2020. They may be identified soon by deep degree-scale JWST surveys i.e., COSMOS-Web (Kartaltepe et al. 2021; Casey et al. 2023) and possibly by narrower ones including the JWST Advanced Deep Extragalactic Survey (JADES, Eisenstein et al. 2017; Ferruit 2017), the Cosmic Evolution Early Release Science Survey (CEERS, Finkelstein et al. 2017), the Next Generation Deep Extragalactic Exploratory Public Survey (NGDEEP, Finkelstein et al. 2021), and Ultra-deep NIRCam and NIRSpec Observations Before the Epoch of Reionization (UNCOVER, Labbé et al. 2021; Bezanson et al. 2022) to name a few.

While low mass quiescent systems at $z \ge 2.5$ may be found by JWST (with hints emerging from Marchesini et al. 2023), a truly precise quantification of their number density and demographics is likely to remain a future objective. By extrapolating these observations, <1% of $M \approx 10^{9.5-10.0} M_\odot$ are expected to be quiescent by $z \sim 2.5$, become even rarer at earlier times. While identifying even one quiescent, low-mass system at $z > 2$ in the absence of virialized structure would present a significant challenge to the paradigm of Peng et al. (2010), performing a statistically meaningful survey of them will require incredibly deep, degree-scale NIR surveys, placing them out of reach by current facilities. While the deepest degree-scale surveys from JWST (COSMOS-Web), *Roman*, and *Euclid* stand to establish the rarity of these systems at $z \approx 2-3$, it seems that no currently planned survey will be able to definitively quantify their contribution at $z \ge 3$.

While explorations of mass-selected samples at $z > 7$ are being made possible by JWST, the most UV-luminous sources at these epochs will likely be missed by small area programs and yet are crucial to a comprehensive study of the ionizing UV budget (Kauffmann et al. 2022; Donnan et al. 2023). Following up known samples with JWST will not only allow us to measure the star formation rate and dust content of the first ultra-luminous galaxies from deep within the epoch of reionization, but also to directly identify the progenitors of $z \sim 3-4$ massive galaxies while still in their formation stage (e.g., Weaver et al. 2021). Yet despite the incredible promises of highly resolved IR imaging and spectroscopy from JWST, identifying the rarest and potentially most informative populations that can challenge and thereby improve galaxy formation models will remain the domain of wide-field surveys.

*Acknowledgements.* The authors thank Vadim Ruskov, Sidney Lower, Desika Narayanan, Stephen Wilkins, William Roper, Lukas Furtak, and Pratika Dayal for helpful discussions. We are also grateful for the many helpful and constructive comments from the anonymous referee. The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation under grant No. 140. S.T. and J.W. acknowledge support from the European Research Council (ERC) Consolidator Grant funding scheme (project ConTExt, grant

<footnote n="14">¹⁴ https://doi.org/10.5281/zenodo.7808832</footnote>
No. 648179). O.I. acknowledges the funding of the French Agence Nationale de la Recherche for the project iMAGE (grant ANR-22-CE31-0007). H.J.McC. acknowledges support from the PNCG. I.D. has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 896225. This work used the CANDIDE computer system at the IAP supported by grants from the PNCG, CNES, and the DIM-ACAV and maintained by S. Rouberol. B.M.J. is supported in part by Independent Research Fund Denmark grant DFF - 7014-00017. G.E.M. acknowledges the Villum Fonden research grant 13160 "Gas to stars, stars to dust: tracing star formation across cosmic time". D.R. acknowledges support from the National Science Foundation under grant numbers AST-1614213 and AST-1910107. Y.P. acknowledges National Science Foundation of China (NSFC) Grant No. 12125301, 12192220, 12192222, and the science research grants from the China Manned Space Project with No. CMS-CSST-2021-A07. The authors wish to recognize and acknowledge the very significant cultural role and relevance that the summit of Mauna Kea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. This work is based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under ESO program ID 179.A-2005 and on data products produced by CALET and the Cambridge Astronomy Survey Unit on behalf of the UltraVISTA consortium. This work is based in part on observations made with the NASA/ESA Hubble Space Telescope, obtained from the Data Archive at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. Some of the data presented herein were obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W.M. Keck Foundation. This research is also partly supported by the Centre National d'Etudes Spatiales (CNES). These data were obtained and processed as part of the CFHT Large Area U-band Deep Survey (CLAUDS), which is a collaboration between astronomers from Canada, France, and China described in Sawicki et al. (2019). CLAUDS is based on observations obtained with MegaPrime/MegaCam, a joint project of CFHT and CEA/DAPNIA, at the CFHT which is operated by the National Research Council (NRC) of Canada, the Institut National des Sciences de l'Univers of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. CLAUDS uses data obtained in part through the Telescope Access Program (TAP), which has been funded by the National Astronomical Observatories, Chinese Academy of Sciences, and the Special Fund for Astronomy from the Ministry of Finance of China. CLAUDS uses data products from TERAPIX and the Canadian Astronomy Data Centre (CADC) and was carried out using resources from Compute Canada and Canadian Advanced Network For Astrophysical Research (CANFAR).

## References

Adams, N. J., Bowler, R. A. A., Jarvis, M. J., Häußler, B., & Lagos, C. D. P. 2021, MNRAS, 506, 4933

Aihara, H., AlSayyad, Y., Ando, M., et al. 2019, PASJ, 71, 114

Akhshik, M., Whitaker, K. E., Leja, J., et al. 2023, ApJ, 943, 179

Arnouts, S., Moscardini, L., Vanzella, E., et al. 2002, MNRAS, 329, 355

Arnouts, S., Le Floc'h, E., Chevallard, J., et al. 2013, A&A, 558, A67

Balogh, M. L., Navarro, J. F., & Morris, S. L. 2000, ApJ, 540, 113

Barrufet, L., Oesch, P., Fudamoto, Y., et al. 2021, JWST Proposal, Cycle 1, 2198

Behroozi, P., & Silk, J. 2018, MNRAS, 477, 5382

Behroozi, P. S., Conroy, C., & Wechsler, R. H. 2010, ApJ, 717, 379

Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57

Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143

Belli, S., Newman, A. B., & Ellis, R. S. 2019, ApJ, 874, 17

Bernardi, M., Meert, A., Sheth, R. K., et al. 2017, MNRAS, 467, 2217

Bertin, E. 2010, Astrophysics Source Code Library [record ascl:1010.068]

Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393

Bezanson, R., Labbé, I., Whitaker, K. E., et al. 2022, ApJ, submitted [arXiv:2121.04026]

Bhowmick, A. K., Somerville, R. S., Di Matteo, T., et al. 2020, MNRAS, 496, 754

Bieri, R., Dubois, Y., Rosdahl, J., et al. 2017, MNRAS, 464, 1854

Binggeli, B., & Jerjen, H. 1998, A&A, 333, 17

Birnboim, Y., & Dekel, A. 2003, MNRAS, 345, 349

Bolzonella, M., Kovač, K., Pozzetti, L., et al. 2010, A&A, 524, A76

Bouché, N., Dekel, A., Genzel, R., et al. 2010, ApJ, 718, 1001

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2012a, ApJ, 754, 83

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2012b, ApJ, 752, L5

Bowler, R. A. A., Dunlop, J. S., McLure, R. J., et al. 2012, MNRAS, 426, 2772

Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503

Brinch, M., Greve, T. R., Weaver, J. R., et al. 2023, ApJ, 943, 153

Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151

Caputi, K. I., Cirasuolo, M., Dunlop, J. S., et al. 2011, MNRAS, 413, 162

Caputi, K. I., Ilbert, O., Laigle, C., et al. 2015, ApJ, 810, 73

Carnall, A. C., McLure, R. J., Dunlop, J. S., & Davé, R. 2018, MNRAS, 480, 4379

Carnall, A. C., Leja, J., Johnson, B. D., et al. 2019, ApJ, 873, 44

Carnall, A., Begley, R. A., Cimatti, A., et al. 2021, JWST Proposal, Cycle 1, 2285

Casey, C. M., Narayanan, D., & Cooray, A. 2014, Phys. Rep., 541, 45

Casey, C. M., Zavala, J. A., Manning, S. M., et al. 2021, ApJ, 923, 215

Casey, C. M., Kartaltepe, J. S., Drakos, N. E., et al. 2023, ApJ, 954, 31

Cecchi, R., Bolzonella, M., Cimatti, A., & Girelli, G. 2019, ApJ, 880, L14

Chabrier, G. 2003, PASP, 115, 763

Cicone, C., Maiolino, R., Sturm, E., et al. 2014, A&A, 562, A21

Cicone, C., Maiolino, R., Gallerani, S., et al. 2015, A&A, 574, A14

Civano, F., Marchesi, S., Comastri, A., et al. 2016, ApJ, 819, 62

Cole, S., Norberg, P., Baugh, C. M., et al. 2001, MNRAS, 326, 255

Conroy, C., & Gunn, J. E. 2010, ApJ, 712, 833

Conroy, C., & Wechsler, R. H. 2009, ApJ, 696, 620

Conroy, C., Gunn, J. E., & White, M. 2009, ApJ, 699, 486

Conselice, C. J., Mundy, C. J., Ferreira, L., & Duncan, K. 2022, ApJ, 940, 168

Cooper, O. R., Casey, C. M., Zavala, J. A., et al. 2022, ApJ, 930, 32

Cowie, L. L., Songaila, A., Hu, E. M., & Cohen, J. G. 1996, AJ, 112, 839

Cowley, W. I., Caputi, K. I., Deshmukh, S., et al. 2019, ApJ, 874, 114

Crain, R. A., Schaye, J., Bower, R. G., et al. 2015, MNRAS, 450, 1937

Croton, D. J. 2013, PASA, 30, e052

Daddi, E., Cimatti, A., Renzini, A., et al. 2004, ApJ, 617, 746

Daddi, E., Dickinson, M., Morrison, G., et al. 2007, ApJ, 670, 156

Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827

Davidzon, I., Cucciati, O., Bolzonella, M., et al. 2016, A&A, 586, A23

Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, A&A, 605, A70

Davidzon, I., Ilbert, O., Faisst, A. L., Sparre, M., & Capak, P. L. 2018, ApJ, 852, 107

Debuhr, J., Quataert, E., & Ma, C.-P. 2012, MNRAS, 420, 2221

Di Teodoro, E. M., Posti, L., Fall, S. M., et al. 2023, MNRAS, 518, 6340

Diemer, B. 2018, ApJS, 239, 35

Donley, J. L., Koekemoer, A. M., Brusa, M., et al. 2012, ApJ, 748, 142

Donnan, C. T., McLeod, D. J., Dunlop, J. S., et al. 2023, MNRAS, 518, 6011

Donnari, M., Pillepich, A., Joshi, G. D., et al. 2021a, MNRAS, 500, 4004

Donnari, M., Pillepich, A., Nelson, D., et al. 2021b, MNRAS, 506, 4760

Driver, S. P., & Robotham, A. S. G. 2010, MNRAS, 407, 2131

Drory, N., Bundy, K., Leauthaud, A., et al. 2009, ApJ, 707, 1595

Dubois, Y., Pichon, C., Welker, C., et al. 2014, MNRAS, 444, 1453

Duncan, K., Conselice, C. J., Mortlock, A., et al. 2014, MNRAS, 444, 2960

Ebeling, H. 2003, MNRAS, 340, 1269

Eddington, A. S. 1913, MNRAS, 73, 359

Efstathiou, G., & Bond, J. R. 1999, MNRAS, 304, 75

Efstathiou, G., Ellis, R. S., & Peterson, B. A. 1988, MNRAS, 232, 431

Eisenstein, D. J., Ferruit, P., & Rieke, M. J. 2017, JWST Proposal, Cycle 1, 1181

Euclid Collaboration (Moneti, A., et al.) 2022, A&A, 658, A126

Faisst, A. L., Carollo, C. M., Capak, P. L., et al. 2017, ApJ, 839, 71

Ferruit, P. 2017, JWST Proposal, Cycle 1, 1210

Finkelstein, S. L., Dickinson, M., Ferguson, H. C., et al. 2017, JWST Proposal, 1345, Cycle 0 Early Release Science

Finkelstein, S. L., Papovich, C., Pirzkal, N., et al. 2021, JWST Proposal, Cycle 1, 2079

Fontana, A., Salimbeni, S., Grazian, A., et al. 2006, A&A, 459, 745

Fontana, A., Dunlop, J. S., Paris, D., et al. 2014, A&A, 570, A11

Fontanot, F., De Lucia, G., Monaco, P., Somerville, R. S., & Santini, P. 2009, MNRAS, 397, 1776

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Forrest, B., Annunziatella, M., Wilson, G., et al. 2020a, ApJ, 890, L1

Forrest, B., Marsan, Z. C., Annunziatella, M., et al. 2020b, ApJ, 903, 47

Förster Schreiber, N. M., & Wuyts, S. 2020, ARA&A, 58, 661

Fudamoto, Y., Oesch, P. A., Faisst, A., et al. 2020, A&A, 643, A4

Fudamoto, Y., Oesch, P. A., Schouws, S., et al. 2021, Nature, 597, 489

Furlong, M., Bower, R. G., Theuns, T., et al. 2015, MNRAS, 450, 4486

Gabor, J. M., & Davé, R. 2015, MNRAS, 447, 374

Gehrels, N. 1986, ApJ, 303, 336

Glazebrook, K., Peacock, J. A., Miller, L., & Collins, C. A. 1995, MNRAS, 275, 169

Glazebrook, K., Schreiber, C., Labbé, I., et al. 2017, Nature, 544, 71

Glazebrook, K., Nanayakkara, T., Esdaile, J., et al. 2021, JWST Proposal, Cycle 1, 2565

Gobat, R., Strazzullo, V., Daddi, E., et al. 2012, ApJ, 759, L44

González, V., Labbé, I., Bouwens, R. J., et al. 2011, ApJ, 735, L34
No. 648179). O.I. acknowledges the funding of the French Agence Nationale de la Recherche for the project iMAGE (grant ANR-22-CE31-0007). H.J.McC. acknowledges support from the PNCG. I.D. has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 896225. This work used the CANDIDE computer system at the IAP supported by grants from the PNCG, CNES, and the DIM-ACAV and maintained by S. Rouberol. B.M.J. is supported in part by Independent Research Fund Denmark grant DFF - 7014-00017. G.E.M. acknowledges the Villum Fonden research grant 13160 "Gas to stars, stars to dust: tracing star formation across cosmic time". D.R. acknowledges support from the National Science Foundation under grant numbers AST-1614213 and AST-1910107. Y.P. acknowledges National Science Foundation of China (NSFC) Grant No. 12125301, 12192220, 12192222, and the science research grants from the China Manned Space Project with No. CMS-CSST-2021-A07. The authors wish to recognize and acknowledge the very significant cultural role and relevance that the summit of Mauna Kea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. This work is based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under ESO program ID 179.A-2005 and on data products produced by CALET and the Cambridge Astronomy Survey Unit on behalf of the UltraVISTA consortium. This work is based in part on observations made with the NASA/ESA Hubble Space Telescope, obtained from the Data Archive at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. Some of the data presented herein were obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W.M. Keck Foundation. This research is also partly supported by the Centre National d'Etudes Spatiales (CNES). These data were obtained and processed as part of the CFHT Large Area U-band Deep Survey (CLAUDS), which is a collaboration between astronomers from Canada, France, and China described in Sawicki et al. (2019). CLAUDS is based on observations obtained with MegaPrime/MegaCam, a joint project of CFHT and CEA/DAPNIA, at the CFHT which is operated by the National Research Council (NRC) of Canada, the Institut National des Sciences de l'Université of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. CLAUDS uses data obtained in part through the Telescope Access Program (TAP), which has been funded by the National Astronomical Observatories, Chinese Academy of Sciences, and the Special Fund for Astronomy from the Ministry of Finance of China. CLAUDS uses data products from TERAPIX and the Canadian Astronomy Data Centre (CADC) and was carried out using resources from Compute Canada and Canadian Advanced Network For Astrophysical Research (CANFAR).

## References

Adams, N. J., Bowler, R. A. A., Jarvis, M. J., Häußler, B., & Lagos, C. D. P. 2021, MNRAS, 506, 4933

Aihara, H., AlSayyad, Y., Ando, M., et al. 2019, PASJ, 71, 114

Akhshik, M., Whitaker, K. E., Leja, J., et al. 2023, ApJ, 943, 179

Arnouts, S., Moscardini, L., Vanzella, E., et al. 2002, MNRAS, 329, 355

Arnouts, S., Le Floc'h, E., Chevallard, J., et al. 2013, A&A, 558, A67

Balogh, M. L., Navarro, J. F., & Morris, S. L. 2000, ApJ, 540, 113

Barrufet, L., Oesch, P., Fudamoto, Y., et al. 2021, JWST Proposal, Cycle 1, 2198

Behroozi, P., & Silk, J. 2018, MNRAS, 477, 5382

Behroozi, P. S., Conroy, C., & Wechsler, R. H. 2010, ApJ, 717, 379

Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57

Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143

Belli, S., Newman, A. B., & Ellis, R. S. 2019, ApJ, 874, 17

Bernardi, M., Meert, A., Sheth, R. K., et al. 2017, MNRAS, 467, 2217

Bertin, E. 2010, Astrophysics Source Code Library [record ascl:1010.068]

Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393

Bezanson, R., Labbé, I., Whitaker, K. E., et al. 2022, ApJ, submitted [arXiv:2121.04026]

Bhowmick, A. K., Somerville, R. S., Di Matteo, T., et al. 2020, MNRAS, 496, 754

Bieri, R., Dubois, Y., Rosdahl, J., et al. 2017, MNRAS, 464, 1854

Binggeli, B., & Jerjen, H. 1998, A&A, 333, 17

Birnboim, Y., & Dekel, A. 2003, MNRAS, 345, 349

Bolzonella, M., Kovač, K., Pozzetti, L., et al. 2010, A&A, 524, A76

Bouché, N., Dekel, A., Genzel, R., et al. 2010, ApJ, 718, 1001

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2012a, ApJ, 754, 83

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2012b, ApJ, 752, L5

Bowler, R. A. A., Dunlop, J. S., McLure, R. J., et al. 2012, MNRAS, 426, 2772

Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503

Brinch, M., Greve, T. R., Weaver, J. R., et al. 2023, ApJ, 943, 153

Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151

Caputi, K. I., Cirasuolo, M., Dunlop, J. S., et al. 2011, MNRAS, 413, 162

Caputi, K. I., Ilbert, O., Laigle, C., et al. 2015, ApJ, 810, 73

Carnall, A. C., McLure, R. J., Dunlop, J. S., & Davé, R. 2018, MNRAS, 480, 4379

Carnall, A. C., Leja, J., Johnson, B. D., et al. 2019, ApJ, 873, 44

Carnall, A., Begley, R. A., Cimatti, A., et al. 2021, JWST Proposal, Cycle 1, 2285

Casey, C. M., Narayanan, D., & Cooray, A. 2014, Phys. Rep., 541, 45

Casey, C. M., Zavala, J. A., Manning, S. M., et al. 2021, ApJ, 923, 215

Casey, C. M., Kartaltepe, J. S., Drakos, N. E., et al. 2023, ApJ, 954, 31

Cecchi, R., Bolzonella, M., Cimatti, A., & Girelli, G. 2019, ApJ, 880, L14

Chabrier, G. 2003, PASP, 115, 763

Cicone, C., Maiolino, R., Sturm, E., et al. 2014, A&A, 562, A21

Cicone, C., Maiolino, R., Gallerani, S., et al. 2015, A&A, 574, A14

Civano, F., Marchesi, S., Comastri, A., et al. 2016, ApJ, 819, 62

Cole, S., Norberg, P., Baugh, C. M., et al. 2001, MNRAS, 326, 255

Conroy, C., & Gunn, J. E. 2010, ApJ, 712, 833

Conroy, C., & Wechsler, R. H. 2009, ApJ, 696, 620

Conroy, C., Gunn, J. E., & White, M. 2009, ApJ, 699, 486

Conselice, C. J., Mundy, C. J., Ferreira, L., & Duncan, K. 2022, ApJ, 940, 168

Cooper, O. R., Casey, C. M., Zavala, J. A., et al. 2022, ApJ, 930, 32

Cowie, L. L., Songaila, A., Hu, E. M., & Cohen, J. G. 1996, AJ, 112, 839

Cowie, W. I., Caputi, K. I., Deshmukh, S., et al. 2019, ApJ, 874, 114

Crain, R. A., Schaye, J., Bower, R. G., et al. 2015, MNRAS, 450, 1937

Croton, D. J. 2013, PASA, 30, e052

Daddi, E., Cimatti, A., Renzini, A., et al. 2004, ApJ, 617, 746

Daddi, E., Dickinson, M., Morrison, G., et al. 2007, ApJ, 670, 156

Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827

Davidzon, I., Cucciati, O., Bolzonella, M., et al. 2016, A&A, 586, A23

Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, A&A, 605, A70

Davidzon, I., Ilbert, O., Faisst, A. L., Sparre, M., & Capak, P. L. 2018, ApJ, 852, 107

Debuhr, J., Quataert, E., & Ma, C.-P. 2012, MNRAS, 420, 2221

Di Teodoro, E. M., Posti, L., Fall, S. M., et al. 2023, MNRAS, 518, 6340

Diemer, B. 2018, ApJS, 239, 35

Donley, J. L., Koekemoer, A. M., Brusa, M., et al. 2012, ApJ, 748, 142

Donnan, C. T., McLeod, D. J., Dunlop, J. S., et al. 2023, MNRAS, 518, 6011

Donnari, M., Pillepich, A., Joshi, G. D., et al. 2021a, MNRAS, 500, 4004

Donnari, M., Pillepich, A., Nelson, D., et al. 2021b, MNRAS, 506, 4760

Driver, S. P., & Robotham, A. S. G. 2010, MNRAS, 407, 2131

Drory, N., Bundy, K., Leauthaud, A., et al. 2009, ApJ, 707, 1595

Dubois, Y., Pichon, C., Welker, C., et al. 2014, MNRAS, 444, 1453

Duncan, K., Conselice, C. J., Mortlock, A., et al. 2014, MNRAS, 444, 2960

Ebeling, H. 2003, MNRAS, 340, 1269

Eddington, A. S. 1913, MNRAS, 73, 359

Efstathiou, G., & Bond, J. R. 1999, MNRAS, 304, 75

Efstathiou, G., Ellis, R. S., & Peterson, B. A. 1988, MNRAS, 232, 431

Eisenstein, D. J., Ferruit, P., & Rieke, M. J. 2017, JWST Proposal, Cycle 1, 1181

Euclid Collaboration (Moneti, A., et al.) 2022, A&A, 658, A126

Faisst, A. L., Carollo, C. M., Capak, P. L., et al. 2017, ApJ, 839, 71

Ferruit, P. 2017, JWST Proposal, Cycle 1, 1210

Finkelstein, S. L., Dickinson, M., Ferguson, H. C., et al. 2017, JWST Proposal, 1345, Cycle 0 Early Release Science

Finkelstein, S. L., Papovich, C., Pirzkal, N., et al. 2021, JWST Proposal, Cycle 1, 2079

Fontana, A., Salimbeni, S., Grazian, A., et al. 2006, A&A, 459, 745

Fontana, A., Dunlop, J. S., Paris, D., et al. 2014, A&A, 570, A11

Fontanot, F., De Lucia, G., Monaco, P., Somerville, R. S., & Santini, P. 2009, MNRAS, 397, 1776

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Forrest, B., Annunziatella, M., Wilson, G., et al. 2020a, ApJ, 890, L1

Forrest, B., Marsan, Z. C., Annunziatella, M., et al. 2020b, ApJ, 903, 47

Förster Schreiber, N. M., & Wuyts, S. 2020, ARA&A, 58, 661

Fudamoto, Y., Oesch, P. A., Faisst, A., et al. 2020, A&A, 643, A4

Fudamoto, Y., Oesch, P. A., Schouws, S., et al. 2021, Nature, 597, 489

Furlong, M., Bower, R. G., Theuns, T., et al. 2015, MNRAS, 450, 4486

Gabor, J. M., & Davé, R. 2015, MNRAS, 447, 374

Gehrels, N. 1986, ApJ, 303, 336

Glazebrook, K., Peacock, J. A., Miller, L., & Collins, C. A. 1995, MNRAS, 275, 169

Glazebrook, K., Schreiber, C., Labbé, I., et al. 2017, Nature, 544, 71

Glazebrook, K., Nanayakkara, T., Esdaile, J., et al. 2021, JWST Proposal, Cycle 1, 2565

Gobat, R., Strazzullo, V., Daddi, E., et al. 2012, ApJ, 759, L44

González, V., Labbé, I., Bouwens, R. J., et al. 2011, ApJ, 735, L34
Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Comput. Sci., 5, 65

Gould, K. M. L., Brammer, G., Valentino, F., et al. 2023, AJ, 165, 248

Grazian, A., Fontana, A., Santini, P., et al. 2015, A&A, 575, A96

Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35

Gruppioni, C., Béthermin, M., Loiacono, F., et al. 2020, A&A, 643, A8

Gunn, J. E., & Gott, J. R., III 1972, ApJ, 176, 1

Habouzit, M., Somerville, R. S., Li, Y., et al. 2022, MNRAS, 509, 3015

Harikane, Y., Ouchi, M., Ono, Y., et al. 2016, ApJ, 821, 123

Henriques, B. M. B., White, S. D. M., Thomas, P. A., et al. 2015, MNRAS, 451, 2663

Hill, A. R., Muzzin, A., Franx, M., et al. 2017, ApJ, 837, 147

Hogg, D. W. 1999, ArXiv e-prints [arXiv:astro-ph/9905116]

Hogg, D. W., Bovy, J., & Lang, D. 2010, ArXiv e-prints [arXiv:1008.4686]

Hsieh, B.-C., Wang, W.-H., Hsieh, C.-C., et al. 2012, ApJS, 203, 23

Hviding, R. E., Hainline, K. N., Rieke, M., et al. 2022, AJ, 163, 224

Ilbert, O., Tresse, L., Arnouts, S., et al. 2004, MNRAS, 351, 541

Ilbert, O., Tresse, L., Zucca, E., et al. 2005, A&A, 439, 863

Ilbert, O., Arnouts, S., McCracken, H. J., et al. 2006, A&A, 457, 841

Ilbert, O., Salvato, M., Le Floc'h, E., et al. 2010, ApJ, 709, 644

Ilbert, O., McCracken, H. J., Le Fèvre, O., et al. 2013, A&A, 556, A55

Ito, K., Tanaka, M., Miyaji, T., et al. 2022, ApJ, 929, 53

James, F., & Roos, M. 1975, Comput. Phys. Commun., 10, 343

Jin, S., Daddi, E., Liu, D., et al. 2018, ApJ, 864, 56

Jin, S., Daddi, E., Magdis, G. E., et al. 2019, ApJ, 887, 144

Johnston, R. 2011, A&ARv, 19, 41

Kartaltepe, J., Casey, C. M., Bagley, M., et al. 2021, JWST Proposal, Cycle 1, 1727

Kauffmann, O. B., Ilbert, O., Weaver, J. R., et al. 2022, A&A, 667, A65

Kaviraj, S., Laigle, C., Kimm, T., et al. 2017, MNRAS, 467, 4739

Koekemoer, A. M., Aussel, H., Calzetti, D., et al. 2007, ApJS, 172, 196

Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36

Kokorev, V. I., Magdis, G. E., Davidzon, I., et al. 2021, ApJ, 921, 40

Labbé, I., Bezanson, R., Atek, H., et al. 2021, JWST Proposal, Cycle 1, 2561

Lagos, C. D. P., Tobar, R. J., Robotham, A. S. G., et al. 2018, MNRAS, 481, 3573

Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24

Laigle, C., Davidzon, I., Ilbert, O., et al. 2019, MNRAS, 486, 5104

Lang, D., Hogg, D. W., & Mykytyn, D. 2016, Astrophysics Source Code Library [record ascl:1604.0008]

Larson, R. B., Tinsley, B. M., & Caldwell, C. N. 1980, ApJ, 237, 692

Legrand, L., McCracken, H. J., Davidzon, I., et al. 2019, MNRAS, 486, 5468

Leja, J., Carnall, A. C., Johnson, B. D., Conroy, C., & Speagle, J. S. 2019a, ApJ, 876, 3

Leja, J., Tacchella, S., & Conroy, C. 2019b, ApJ, 880, L9

Lilly, S. J., Carollo, C. M., Pipino, A., Renzini, A., & Peng, Y. 2013, ApJ, 772, 119

Liu, D., Lang, P., Magnelli, B., et al. 2019, ApJS, 244, 40

Long, A. S., Casey, C. M., Lagos, C. D. P., ApJ, 953, 11

Lovell, C. C., Vijayan, A. P., Thomas, P. A., et al. 2021, MNRAS, 500, 2127

Lower, S., Narayanan, D., Leja, J., et al. 2020, ApJ, 904, 33

Lower, S., Narayanan, D., Li, Q., & Davé, R. 2023, ApJ, 950, 94

Lustig, P., Strazzullo, V., Remus, R. S., et al. 2023, MNRAS, 518, 5953

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Mair, R. A., McNamara, B. R., Nulsen, P. E. J., Russell, H. R., & Vantyghem, A. N. 2017, MNRAS, 464, 4360

Maiolino, R., Gallerani, S., Neri, R., et al. 2012, MNRAS, 425, L66

Malmquist, K. G. 1920, Meddelanden fran Lunds Astronomiska Observatorium Serie I, 96, 1

Malmquist, K. G. 1922, Meddelanden fran Lunds Astronomiska Observatorium Serie I, 100, 1

Mandelbaum, R., Seljak, U., Kauffmann, G., Hirata, C. M., & Brinkmann, J. 2006, MNRAS, 368, 715

Manning, S. M., Casey, C. M., Zavala, J. A., et al. 2022, ApJ, 925, 23

Marchesini, D., van Dokkum, P. G., Förster Schreiber, N. M., et al. 2009, ApJ, 701, 1765

Marchesini, D., Brammer, G., Morishita, T., et al. 2023, ApJ, 942, L25

Marrone, D. P., Spilker, J. S., Hayward, C. C., et al. 2018, Nature, 553, 51

Marsan, Z. C., Muzzin, A., Marchesini, D., et al. 2022, ApJ, 924, 25

Martis, N. S., Marchesini, D., Brammer, G. B., et al. 2016, ApJ, 827, L25

Mawatari, K., Yamada, T., Fazio, G. G., Huang, J.-S., & Ashby, M. L. N. 2016, PASJ, 68, 46

Mawatari, K., Inoue, A. K., Hashimoto, T., et al. 2020, ApJ, 889, 137

McConachie, I., Wilson, G., Forrest, B., et al. 2022, ApJ, 926, 37

McCracken, H. J., Milvang-Jensen, B., Dunlop, J., et al. 2012, A&A, 544, A156

McLeod, D. J., McLure, R. J., Dunlop, J. S., et al. 2021, MNRAS, 503, 4413

Michałowski, M. J., Dunlop, J. S., Cirasuolo, M., et al. 2012, A&A, 541, A85

Michałowski, M. J., Hayward, C. C., Dunlop, J. S., et al. 2014, A&A, 571, A75

Milvang-Jensen, B., Freudling, W., Zabl, J., et al. 2013, A&A, 560, A94

Mitchell, P. D., Schaye, J., Bower, R. G., & Crain, R. A. 2020, MNRAS, 494, 3971

Moneti, A., McCracken, H. J., Rouberol, S., et al. 2023, VizieR On-line Data Catalog: II/373

Mortlock, A., Conselice, C. J., Bluck, A. F. L., et al. 2011, MNRAS, 413, 2845

Mortlock, A., Conselice, C. J., Hartley, W. G., et al. 2015, MNRAS, 447, 2

Moster, B. P., Somerville, R. S., Newman, J. A., & Rix, H.-W. 2011, ApJ, 731, 113

Moster, B. P., Naab, T., & White, S. D. M. 2013, MNRAS, 428, 3121

Moutard, T., Arnouts, S., Ilbert, O., et al. 2016, A&A, 590, A103

Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJ, 777, 18

Nayyeri, H., Hemmati, S., Mobasher, B., et al. 2017, ApJS, 228, 7

Neistein, E., & Dekel, A. 2008, MNRAS, 388, 1792

Neistein, E., van den Bosch, F. C., & Dekel, A. 2006, MNRAS, 372, 933

Noeske, K. G., Weiner, B. J., Faber, S. M., et al. 2007, ApJ, 660, L43

Oke, J. B. 1974, ApJS, 27, 21

Pacifi, C., da Cunha, E., Charlot, S., et al. 2015, MNRAS, 447, 786

Padovani, P., Bonzini, M., Kellermann, K. I., et al. 2015, MNRAS, 452, 1263

Peng, Y.-J., Lilly, S. J., Kovač, K., et al. 2010, ApJ, 721, 193

Peng, Y.-J., Lilly, S. J., Renzini, A., & Carollo, M. 2012, ApJ, 757, 4

Peng, Y.-J., Lilly, S. J., Renzini, A., & Carollo, M. 2014, ApJ, 790, 95

Peng, Y., Maiolino, R., & Cochrane, R. 2015, Nature, 521, 192

Pillepich, A., Nelson, D., Hernquist, L., et al. 2018, MNRAS, 475, 648

Pope, A., Arestxaga, I., Hughes, D., Wilson, G., & Yun, M. 2019, Am. Astron. Soc. Meet. Abstr., 233, 363.20

Popesso, P., Concas, A., Cresci, G., et al. 2023, MNRAS, 519, 1526

Pozzetti, L., Bolzella, M., Zucca, E., et al. 2010, A&A, 523, A13

Retzlaff, J., Rosati, P., Dickinson, M., et al. 2010, A&A, 511, A50

Richardson, M. L. A., Scannapieco, E., Devriendt, J., et al. 2016, ApJ, 825, 83

Rigby, J., Perrin, M., McElwain, M., et al. 2023, PASP, 135, 048001

Roos, O., Juneau, S., Bournaud, F., & Gabor, J. M. 2015, ApJ, 800, 19

Saintonge, A., Lutz, D., Genzel, R., et al. 2013, ApJ, 778, 2

Salim, S., Rich, R. M., Charlot, S., et al. 2007, ApJS, 173, 267

Salpeter, E. E. 1955, ApJ, 121, 161

Salvato, M., Ilbert, O., Hasinger, G., et al. 2011, ApJ, 742, 61

Sandage, A., Tammam, G. A., & Yahil, A. 1979, ApJ, 232, 352

Sanders, D. B., Salvato, M., Aussel, H., et al. 2007, ApJS, 172, 86

Santini, P., Fontana, A., Grazian, A., et al. 2012, A&A, 538, A33

Santini, P., Castellano, M., Merlin, E., et al. 2012, A&A, 652, A30

Santini, P., Castellano, M., Fontana, A., et al. 2022, ApJ, 940, 135

Sawicki, M., Arnouts, S., Huang, J., et al. 2019, MNRAS, 489, 5202

Schechter, P. 1976, ApJ, 203, 297

Schmidt, M. 1968, ApJ, 151, 393

Schreiber, C., Glazebrook, K., Nanayakkara, T., et al. 2018a, A&A, 618, A85

Schreiber, C., Labbé, I., Glazebrook, K., et al. 2018b, A&A, 611, A22

Scoville, N., Aussel, H., Brusa, M., et al. 2007, ApJS, 172, 1

Shahidi, A., Mobasher, B., Nayyeri, H., et al. 2020, ApJ, 897, 44

Shankar, F., Lapi, A., Salucci, P., De Zotti, G., & Danese, L. 2006, ApJ, 643, 14

Shu, X., Yang, L., Liu, D., et al. 2022, ApJ, 926, 155

Shuntov, M., McCracken, H. J., Gavazzi, R., et al. 2022, A&A, 664, A61

Smolčić, V., Novak, M., Bondi, M., et al. 2017, A&A, 602, A1

Somerville, R. S., Popping, G., & Trager, S. C. 2015, MNRAS, 453, 4337

Song, M., Finkelstein, S. L., Ashby, M. L. N., et al. 2016, ApJ, 825, 5

Songaila, A., Cowie, L. L., Hu, E. M., & Gardner, J. P. 1994, ApJS, 94, 461

Stefanon, M., Labbé, I., Oesch, P. A., et al. 2021, ApJS, 257, 68

Steinhardt, C. L., Capak, P., Masters, D., & Speagle, J. S. 2016, ApJ, 824, 21

Steinhardt, C. L., Weaver, J. R., Maxfield, J., et al. 2020, ApJ, 891, 136

Steinhardt, C. L., Jespersen, C. K., & Linzer, N. B. 2021, ApJ, 923, 8

Stern, D., Eisenhardt, P., Gorjian, V., et al. 2005, ApJ, 631, 163

Straatman, C. M. S., Labbé, I., Spitler, L. R., et al. 2015, ApJ, 808, L29

Suess, K. A., Kriek, M., Price, S. H., & Barro, G. 2021, ApJ, 915, 87

Sun, F., Egami, E., Pérez-González, P. G., et al. 2021, ApJ, 922, 114

Szalay, A. S., Connolly, A. J., & Szokoly, G. P. 1999, AJ, 117, 68

Takeuchi, T. T. 2000, Ap&SS, 271, 213

Tanaka, M. 2015, ApJ, 801, 20

Tanaka, M., Valentino, F., Toft, S., et al. 2019, ApJ, 885, L34

Thorne, J. E., Robotham, A. S. G., Davies, L. J. M., et al. 2021, MNRAS, 505, 540

Tinker, J., Kravsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709

Toft, S., Zabl, J., Richard, J., et al. 2017, Nature, 546, 510

Tomczak, A. R., Quadri, R. F., Tran, K.-V. H., et al. 2014, ApJ, 783, 85

Torrey, P., Vogelsberger, M., Genel, S., et al. 2014, MNRAS, 438, 1985

Trapp, A. C., & Furlanetto, S. R. 2020, MNRAS, 499, 2401

Trenti, M., & Stiavelli, M. 2008, ApJ, 676, 767

Ucci, G., Dayal, P., Hutter, A., et al. 2021, MNRAS, 506, 202

Valentino, F., Tanaka, M., Davidzon, I., et al. 2020, ApJ, 889, 93

Varma, S., Huertas-Company, M., Pillepich, A., et al. 2022, MNRAS, 509, 2654
Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Comput. Sci., 5, 65

Gould, K. M. L., Brammer, G., Valentino, F., et al. 2023, AJ, 165, 248

Grazian, A., Fontana, A., Santini, P., et al. 2015, A&A, 575, A96

Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35

Gruppioni, C., Béthermin, M., Loiacono, F., et al. 2020, A&A, 643, A8

Gunn, J. E., & Gott, J. R., III 1972, ApJ, 176, 1

Habouzit, M., Somerville, R. S., Li, Y., et al. 2022, MNRAS, 509, 3015

Harikane, Y., Ouchi, M., Ono, Y., et al. 2016, ApJ, 821, 123

Henriques, B. M. B., White, S. D. M., Thomas, P. A., et al. 2015, MNRAS, 451, 2663

Hill, A. R., Muzzin, A., Franx, M., et al. 2017, ApJ, 837, 147

Hogg, D. W. 1999, ArXiv e-prints [arXiv:astro-ph/9905116]

Hogg, D. W., Bovy, J., & Lang, D. 2010, ArXiv e-prints [arXiv:1008.4686]

Hsieh, B.-C., Wang, W.-H., Hsieh, C.-C., et al. 2012, ApJS, 203, 23

Hviding, R. E., Hainline, K. N., Rieke, M., et al. 2022, AJ, 163, 224

Ilbert, O., Tresse, L., Arnouts, S., et al. 2004, MNRAS, 351, 541

Ilbert, O., Tresse, L., Zucca, E., et al. 2005, A&A, 439, 863

Ilbert, O., Arnouts, S., McCracken, H. J., et al. 2006, A&A, 457, 841

Ilbert, O., Salvato, M., Le Floc'h, E., et al. 2010, ApJ, 709, 644

Ilbert, O., McCracken, H. J., Le Fèvre, O., et al. 2013, A&A, 556, A55

Ito, K., Tanaka, M., Miyaji, T., et al. 2022, ApJ, 929, 53

James, F., & Roos, M. 1975, Comput. Phys. Commun., 10, 343

Jin, S., Daddi, E., Liu, D., et al. 2018, ApJ, 864, 56

Jin, S., Daddi, E., Magdis, G. E., et al. 2019, ApJ, 887, 144

Johnston, R. 2011, A&ARv, 19, 41

Kartaltepe, J., Casey, C. M., Bagley, M., et al. 2021, JWST Proposal, Cycle 1, 1727

Kauffmann, O. B., Ilbert, O., Weaver, J. R., et al. 2022, A&A, 667, A65

Kaviraj, S., Laigle, C., Kimm, T., et al. 2017, MNRAS, 467, 4739

Koekemoer, A. M., Aussel, H., Calzetti, D., et al. 2007, ApJS, 172, 196

Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36

Kokorev, V. I., Magdis, G. E., Davidzon, I., et al. 2021, ApJ, 921, 40

Labbé, I., Bezanson, R., Atek, H., et al. 2021, JWST Proposal, Cycle 1, 2561

Lagos, C. D. P., Tobar, R. J., Robotham, A. S. G., et al. 2018, MNRAS, 481, 3573

Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24

Laigle, C., Davidzon, I., Ilbert, O., et al. 2019, MNRAS, 486, 5104

Lang, D., Hogg, D. W., & Mykytyn, D. 2016, Astrophysics Source Code Library [record ascl:1604.008]

Larson, R. B., Tinsley, B. M., & Caldwell, C. N. 1980, ApJ, 237, 692

Legrand, L., McCracken, H. J., Davidzon, I., et al. 2019, MNRAS, 486, 5468

Leja, J., Carnall, A. C., Johnson, B. D., Conroy, C., & Speagle, J. S. 2019a, ApJ, 876, 3

Leja, J., Tacchella, S., & Conroy, C. 2019b, ApJ, 880, L9

Lilly, S. J., Carollo, C. M., Pipino, A., Renzini, A., & Peng, Y. 2013, ApJ, 772, 119

Liu, D., Lang, P., Magnelli, B., et al. 2019, ApJS, 244, 40

Long, A. S., Casey, C. M., Lagos, C. D. P., et al. 2023, ApJ, 953, 11

Lovell, C. C., Vijayan, A. P., Thomas, P. A., et al. 2021, MNRAS, 500, 2127

Lower, S., Narayanan, D., Leja, J., et al. 2020, ApJ, 904, 33

Lower, S., Narayanan, D., Li, Q., & Davé, R. 2023, ApJ, 950, 94

Lustig, P., Strazzullo, V., Remus, R. S., et al. 2023, MNRAS, 518, 5953

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Mair, R. A., McNamara, B. R., Nulsen, P. E. J., Russell, H. R., & Vantyghem, A. N. 2017, MNRAS, 464, 4360

Maiolino, R., Gallerani, S., Neri, R., et al. 2012, MNRAS, 425, L66

Malmquist, K. G. 1920, Meddelanden fran Lunds Astronomiska Observatorium Serie I, 96, 1

Malmquist, K. G. 1922, Meddelanden fran Lunds Astronomiska Observatorium Serie I, 100, 1

Mandelbaum, R., Seljak, U., Kauffmann, G., Hirata, C. M., & Brinkmann, J. 2006, MNRAS, 368, 715

Manning, S. M., Casey, C. M., Zavala, J. A., et al. 2022, ApJ, 925, 23

Marchesini, D., van Dokkum, P. G., Förster Schreiber, N. M., et al. 2009, ApJ, 701, 1765

Marchesini, D., Brammer, G., Morishita, T., et al. 2023, ApJ, 942, L25

Marrone, D. P., Spilker, J. S., Hayward, C. C., et al. 2018, Nature, 553, 51

Marsan, Z. C., Muzzin, A., Marchesini, D., et al. 2022, ApJ, 924, 25

Martis, N. S., Marchesini, D., Brammer, G. B., et al. 2016, ApJ, 827, L25

Mawatari, K., Yamada, T., Fazio, G. G., Huang, J.-S., & Ashby, M. L. N. 2016, PASJ, 68, 46

Mawatari, K., Inoue, A. K., Hashimoto, T., et al. 2020, ApJ, 889, 137

McConachie, I., Wilson, G., Forrest, B., et al. 2022, ApJ, 926, 37

McCracken, H. J., Milvang-Jensen, B., Dunlop, J., et al. 2012, A&A, 544, A156

McLeod, D. J., McLure, R. J., Dunlop, J. S., et al. 2021, MNRAS, 503, 4413

Michałowski, M. J., Dunlop, J. S., Cirasuolo, M., et al. 2012, A&A, 541, A85

Michałowski, M. J., Hayward, C. C., Dunlop, J. S., et al. 2014, A&A, 571, A75

Milvang-Jensen, B., Freudling, W., Zabl, J., et al. 2013, A&A, 560, A94

Mitchell, P. D., Schaye, J., Bower, R. G., & Crain, R. A. 2020, MNRAS, 494, 3971

Moneti, A., McCracken, H. J., Rouberol, S., et al. 2023, VizieR On-line Data Catalog: II/373

Mortlock, A., Conselice, C. J., Bluck, A. F. L., et al. 2011, MNRAS, 413, 2845

Mortlock, A., Conselice, C. J., Hartley, W. G., et al. 2015, MNRAS, 447, 2

Moster, B. P., Somerville, R. S., Newman, J. A., & Rix, H.-W. 2011, ApJ, 731, 113

Moster, B. P., Naab, T., & White, S. D. M. 2013, MNRAS, 428, 3121

Moutard, T., Arnouts, S., Ilbert, O., et al. 2016, A&A, 590, A103

Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJ, 777, 18

Neistein, E., & Dekel, A. 2008, MNRAS, 388, 1792

Neistein, E., van den Bosch, F. C., & Dekel, A. 2006, MNRAS, 372, 933

Noeske, K. G., Weiner, B. J., Faber, S. M., et al. 2007, ApJ, 660, L43

Oke, J. B. 1974, ApJS, 27, 21

Pacifici, C., da Cunha, E., Charlot, S., et al. 2015, MNRAS, 447, 786

Padovani, P., Bonzini, M., Kellermann, K. I., et al. 2015, MNRAS, 452, 1263

Peng, Y.-J., Lilly, S. J., Kovač, K., et al. 2010, ApJ, 721, 193

Peng, Y.-J., Lilly, S. J., Renzini, A., & Carollo, M. 2012, ApJ, 757, 4

Peng, Y.-J., Lilly, S. J., Renzini, A., & Carollo, M. 2014, ApJ, 790, 95

Peng, Y., Maiolino, R., & Cochrane, R. 2015, Nature, 521, 192

Pillepich, A., Nelson, D., Hernquist, L., et al. 2018, MNRAS, 475, 648

Pope, A., Arestxaga, I., Hughes, D., Wilson, G., & Yun, M. 2019, Am. Astron. Soc. Meet. Abstr., 233, 363.20

Popesso, P., Concas, A., Cresci, G., et al. 2023, MNRAS, 519, 1526

Pozzetti, L., Bolzolina, M., Zucca, E., et al. 2010, A&A, 523, A13

Retzlaff, J., Rosati, P., Dickinson, M., et al. 2010, A&A, 511, A50

Richardson, M. L. A., Scannapieco, E., Devriendt, J., et al. 2016, ApJ, 825, 83

Rigby, J., Perrin, M., McElwain, M., et al. 2023, PASP, 135, 048001

Roos, O., Juneau, S., Bournaud, F., & Gabor, J. M. 2015, ApJ, 800, 19

Saintonge, A., Lutz, D., Genzel, R., et al. 2013, ApJ, 778, 2

Salim, S., Rich, R. M., Charlot, S., et al. 2007, ApJS, 173, 267

Salpeter, E. E. 1955, ApJ, 121, 161

Salvato, M., Ilbert, O., Hasinger, G., et al. 2011, ApJ, 742, 61

Sandage, A., Tammam, G. A., & Yahil, A. 1979, ApJ, 232, 352

Sanders, D. B., Salvato, M., Aussel, H., et al. 2007, ApJS, 172, 86

Santini, P., Fontana, A., Grazian, A., et al. 2012, A&A, 538, A33

Santini, P., Castellano, M., Merlin, E., et al. 2012, A&A, 652, A30

Santini, P., Castellano, M., Fontana, A., et al. 2022, ApJ, 940, 135

Sawicki, M., Arnouts, S., Huang, J., et al. 2019, MNRAS, 489, 5202

Schechter, P. 1976, ApJ, 203, 297

Schmidt, M. 1968, ApJ, 151, 393

Schreiber, C., Glazebrook, K., Nanayakkara, T., et al. 2018a, A&A, 618, A85

Schreiber, C., Labbé, I., Glazebrook, K., et al. 2018b, A&A, 611, A22

Scoville, N., Aussel, H., Brusa, M., et al. 2007, ApJS, 172, 1

Shahidi, A., Mobasher, B., Nayyeri, H., et al. 2020, ApJ, 897, 44

Shankar, F., Lapi, A., Salucci, P., De Zotti, G., & Danese, L. 2006, ApJ, 643, 14

Shu, X., Yang, L., Liu, D., et al. 2022, ApJ, 926, 155

Shuntov, M., McCracken, H. J., Gavazzi, R., et al. 2022, A&A, 664, A61

Smolčić, V., Novak, M., Bondi, M., et al. 2017, A&A, 602, A1

Somerville, R. S., Popping, G., & Trager, S. C. 2015, MNRAS, 453, 4337

Song, M., Finkelstein, S. L., Ashby, M. L. N., et al. 2016, ApJ, 825, 5

Songaila, A., Cowie, L. L., Hu, E. M., & Gardner, J. P. 1994, ApJS, 94, 461

Stefanon, M., Labbé, I., Oesch, P. A., et al. 2021, ApJS, 257, 68

Steinhardt, C. L., Capak, P., Masters, D., & Speagle, J. S. 2016, ApJ, 824, 21

Steinhardt, C. L., Weaver, J. R., Maxfield, J., et al. 2020, ApJ, 891, 136

Steinhardt, C. L., Jespersen, C. K., & Linzer, N. B. 2021, ApJ, 923, 8

Stern, D., Eisenhardt, P., Gorjian, V., et al. 2005, ApJ, 631, 163

Straatman, C. M. S., Labbé, I., Spitler, L. R., et al. 2015, ApJ, 808, L29

Suess, K. A., Kriek, M., Price, S. H., & Barro, G. 2021, ApJ, 915, 87

Sun, F., Egami, E., Pérez-González, P. G., et al. 2021, ApJ, 922, 114

Szalay, A. S., Connolly, A. J., & Szokoly, G. P. 1999, AJ, 117, 68

Takeuchi, T. T. 2000, Ap&SS, 271, 213

Tanaka, M. 2015, ApJ, 801, 20

Tanaka, M., Valentino, F., Toft, S., et al. 2019, ApJ, 885, L34

Thorne, J. E., Robotham, A. S. G., Davies, L. J. M., et al. 2021, MNRAS, 505, 540

Tinker, J., Kravsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709

Toft, S., Zabl, J., Richard, J., et al. 2017, Nature, 546, 510

Tomczak, A. R., Quadri, R. F., Tran, K.-V. H., et al. 2014, ApJ, 783, 85

Torrey, P., Vogelsberger, M., Genel, S., et al. 2014, MNRAS, 438, 1985

Trapp, A. C., & Furlanetto, S. R. 2020, MNRAS, 499, 2401

Trenti, M., & Stiavelli, M. 2008, ApJ, 676, 767

Ucci, G., Dayal, P., Hutter, A., et al. 2021, MNRAS, 506, 202

Valentino, F., Tanaka, M., Davidzon, I., et al. 2020, ApJ, 889, 93

Varma, S., Huertas-Company, M., Pillepich, A., et al. 2022, MNRAS, 509, 2654
Vito, F., Brandt, W. N., Yang, G., et al. 2018, MNRAS, 473, 2378

Volonteri, M., & Reines, A. E. 2016, ApJ, 820, L6

Vulcani, B., Poggianti, B. M., Oemler, A., et al. 2013, A&A, 550, A58

Wang, T., Schreiber, C., Elbaz, D., et al. 2019, Nature, 572, 211

Weaver, J. R., Brammer, G., Casey, C. M., et al. 2021, JWST Proposal, Cycle 1, 2659

Weaver, J. R., Kauffmann, O. B., Ilbert, O., et al. 2022, ApJS, 258, 11

Weaver, J. R., Zalesky, L., Kokorev, V., et al. 2023, ApJS, submitted

Wechsler, R. H., & Tinker, J. L. 2018, ARA&A, 56, 435

Weigel, A. K., Schawinski, K., & Bruderer, C. 2016, MNRAS, 459, 2150

Whitaker, K. E., van Dokkum, P. G., Brammer, G., & Franx, M. 2012, ApJ, 754, L29

Whitaker, K. E., Franx, M., Leja, J., et al. 2014, ApJ, 795, 104

Whitaker, K. E., Williams, C. C., Mowla, L., et al. 2021, Nature, 597, 485

Williams, R. J., Quadri, R. F., Franx, M., van Dokkum, P., & Labbé, I. 2009, ApJ, 691, 1879

Wright, A. H., Driver, S. P., & Robotham, A. S. G. 2018, MNRAS, 480, 3491

Yung, L. Y. A., Somerville, R. S., Finkelstein, S. L., Popping, G., & Davé, R. 2019a, MNRAS, 483, 2983

Yung, L. Y. A., Somerville, R. S., Popping, G., et al. 2019b, MNRAS, 490, 2855

Zavala, J. A., Casey, C. M., Manning, S. M., et al. 2021, ApJ, 909, 165

<footnote n="1">¹ Cosmic Dawn Center (DAWN), Copenhagen N, Denmark
e-mail: john.weaver.astro@gmail.com</footnote>

<footnote n="2">² Niels Bohr Institute, University of Copenhagen, Jagtvej 128, 2200 Copenhagen N, Denmark</footnote>

<footnote n="3">³ Department of Astronomy, University of Massachusetts, Amherst, MA 01003, USA</footnote>

<footnote n="4">⁴ Aix Marseille Univ., CNRS, CNES, LAM, 13000 Marseille, France</footnote>

<footnote n="5">⁵ Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, 75014 Paris, France</footnote>

<footnote n="6">⁶ Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA</footnote>

<footnote n="7">⁷ International Centre for Radio Astronomy Research (ICRAR), M468, University of Western Australia, 35 Stirling Hwy, Crawley, WA 6009, Australia</footnote>

<footnote n="8">⁸ ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Stromlo, Australia</footnote>

<footnote n="9">⁹ Department of Astronomy, The University of Texas at Austin, 2515 Speedway Blvd Stop C1400, Austin, TX 78712, USA</footnote>

<footnote n="10">¹⁰ The Observatories of the Carnegie Institution for Science, 813 Santa Barbara St., Pasadena, CA 91101, USA</footnote>

<footnote n="11">¹¹ Caltech/IPAC, MS314-6, 1200 E. California Blvd, Pasadena, CA 91125, USA</footnote>

<footnote n="12">¹² Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA</footnote>

<footnote n="13">¹³ Laboratory for Multiwavelength Astrophysics, School of Physics and Astronomy, Rochester Institute of Technology, 84 Lomb Memorial Drive, Rochester, NY 14623, USA</footnote>

<footnote n="14">¹⁴ Space Telescope Science Institute, 3700 San Martin Dr., Baltimore, MD 21218, USA</footnote>

<footnote n="15">¹⁵ Max-Planck-Institut für Extraterrestrische Physik (MPE), Giessen-bachstr. 1, 85748 Garching, Germany</footnote>

<footnote n="16">¹⁶ DTU-Space, Technical University of Denmark, Elektrovej 327, 2800 Kgs. Lyngby, Denmark</footnote>

<footnote n="17">¹⁷ Physics and Astronomy Department, University of California, 900 University Avenue, Riverside, CA 92521, USA</footnote>

<footnote n="18">¹⁸ Kavli Institute for Astronomy and Astrophysics, Peking University, 5 Yiheyuan Road, Beijing 100871, PR China</footnote>

<footnote n="19">¹⁹ Department of Astronomy, School of Physics, Peking University, 5 Yiheyuan Road, Beijing 100871, PR China</footnote>

<footnote n="20">²⁰ Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu, HI 96822, USA</footnote>

<footnote n="21">²¹ Istituto Nazionale di Astrofisica – Osservatorio di Astrofisica e Scienza dello Spazio, Via Gobetti 93/3, 40129 Bologna, Italy</footnote>
Vito, F., Brandt, W. N., Yang, G., et al. 2018, MNRAS, 473, 2378

Volonteri, M., & Reines, A. E. 2016, ApJ, 820, L6

Vulcani, B., Poggianti, B. M., Oemler, A., et al. 2013, A&A, 550, A58

Wang, T., Schreiber, C., Elbaz, D., et al. 2019, Nature, 572, 211

Weaver, J. R., Brammer, G., Casey, C. M., et al. 2021, JWST Proposal, Cycle 1, 2659

Weaver, J. R., Kauffmann, O. B., Ilbert, O., et al. 2022, ApJS, 258, 11

Weaver, J. R., Zalesky, L., Kokorev, V., et al. 2023, ApJS, submitted

Wechsler, R. H., & Tinker, J. L. 2018, ARA&A, 56, 435

Weigel, A. K., Schawinski, K., & Bruderer, C. 2016, MNRAS, 459, 2150

Whitaker, K. E., van Dokkum, P. G., Brammer, G., & Franx, M. 2012, ApJ, 754, L29

Whitaker, K. E., Franx, M., Leja, J., et al. 2014, ApJ, 795, 104

Whitaker, K. E., Williams, C. C., Mowla, L., et al. 2021, Nature, 597, 485

Williams, R. J., Quadri, R. F., Franx, M., van Dokkum, P., & Labbé, I. 2009, ApJ, 691, 1879

Wright, A. H., Driver, S. P., & Robotham, A. S. G. 2018, MNRAS, 480, 3491

Yung, L. Y. A., Somerville, R. S., Finkelstein, S. L., Popping, G., & Davé, R. 2019a, MNRAS, 483, 2983

Yung, L. Y. A., Somerville, R. S., Popping, G., et al. 2019b, MNRAS, 490, 2855

Zavala, J. A., Casey, C. M., Manning, S. M., et al. 2021, ApJ, 909, 165

<footnote n="1">¹ Cosmic Dawn Center (DAWN), Copenhagen N, Denmark
e-mail: john.weaver.astro@gmail.com</footnote>

<footnote n="2">² Niels Bohr Institute, University of Copenhagen, Jagtvej 128, 2200 Copenhagen N, Denmark</footnote>

<footnote n="3">³ Department of Astronomy, University of Massachusetts, Amherst, MA 01003, USA</footnote>

<footnote n="4">⁴ Aix Marseille Univ., CNRS, CNES, LAM, 13000 Marseille, France</footnote>

<footnote n="5">⁵ Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, 75014 Paris, France</footnote>

<footnote n="6">⁶ Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA</footnote>

<footnote n="7">⁷ International Centre for Radio Astronomy Research (ICRAR), M468, University of Western Australia, 35 Stirling Hwy, Crawley, WA 6009, Australia</footnote>

<footnote n="8">⁸ ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Stromlo, Australia</footnote>

<footnote n="9">⁹ Department of Astronomy, The University of Texas at Austin, 2515 Speedway Blvd Stop C1400, Austin, TX 78712, USA</footnote>

<footnote n="10">¹⁰ The Observatories of the Carnegie Institution for Science, 813 Santa Barbara St., Pasadena, CA 91101, USA</footnote>

<footnote n="11">¹¹ Caltech/IPAC, MS314-6, 1200 E. California Blvd, Pasadena, CA 91125, USA</footnote>

<footnote n="12">¹² Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA</footnote>

<footnote n="13">¹³ Laboratory for Multiwavelength Astrophysics, School of Physics and Astronomy, Rochester Institute of Technology, 84 Lomb Memorial Drive, Rochester, NY 14623, USA</footnote>

<footnote n="14">¹⁴ Space Telescope Science Institute, 3700 San Martin Dr., Baltimore, MD 21218, USA</footnote>

<footnote n="15">¹⁵ Max-Planck-Institut für Extraterrestrische Physik (MPE), Giessen-bachstr. 1, 85748 Garching, Germany</footnote>

<footnote n="16">¹⁶ DTU-Space, Technical University of Denmark, Elektrovej 327, 2800 Kgs. Lyngby, Denmark</footnote>

<footnote n="17">¹⁷ Physics and Astronomy Department, University of California, 900 University Avenue, Riverside, CA 92521, USA</footnote>

<footnote n="18">¹⁸ Kavli Institute for Astronomy and Astrophysics, Peking University, 5 Yiheyuan Road, Beijing 100871, PR China</footnote>

<footnote n="19">¹⁹ Department of Astronomy, School of Physics, Peking University, 5 Yiheyuan Road, Beijing 100871, PR China</footnote>

<footnote n="20">²⁰ Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu, HI 96822, USA</footnote>

<footnote n="21">²¹ Istituto Nazionale di Astrofisica – Osservatorio di Astrofisica e Scienza dello Spazio, Via Gobetti 93/3, 40129 Bologna, Italy</footnote>
**Appendix A: Sample selection**

As discussed in Section 3, we select sources from the 1.27 deg² COMBINED region for sources 0.2 < z ≤ 6.5, and impose an additional restriction of the ultra-deep stripe for sources 6.5 < z ≤ 7.5. The three remaining criteria are aimed to ensure secure photo-z and *M* based on the IRAC ch1 magnitude, SED χ² fit quality, and *L*(z). As shown in Fig. A.1, there is considerable overlap between these three criteria, with m<sub>ch1</sub} > 26 accounting for 94% of the sample removed by combining these three.

<figure><img src="image_18.png" /><figcaption>**Fig. A.1.** Fraction of sources removed by a single criterion that are also removed by another, after applying the primary criteria on area and photo-z. E.g., 71% of χ²<sub>SED</sub> > 10 sources have m<sub>ch1</sub> > 26.</figcaption></figure>
**Appendix A: Sample selection**

As discussed in Section 3, we select sources from the 1.27 deg² COMBINED region for sources 0.2 < z ≤ 6.5, and impose an additional restriction of the ultra-deep stripe for sources 6.5 < z ≤ 7.5. The three remaining criteria are aimed to ensure secure photo-z and *M* based on the IRAC ch1 magnitude, SED χ² fit quality, and *L*(z). As shown in Fig. A.1, there is considerable overlap between these three criteria, with m<sub>ch1</sub} > 26 accounting for 94% of the sample removed by combining these three.

<figure><img src="image_19.png" /><figcaption>**Fig. A.1.** Fraction of sources removed by a single criterion that are also removed by another, after applying the primary criteria on area and photo-z. E.g., 71% of χ²<sub>SED</sub> > 10 sources have m<sub>ch1</sub> > 26.</figcaption></figure>
**Appendix B: Further comparisons with simulations**

Here we include additional figures comparing the total, star-
forming, and quiescent mass functions (Fig. B.1) and quies-

cent fractions (Fig. B.2) from this work to those of UNIVERSE
MACHINE Behroozi et al. (2019, sSFR < 10⁻¹¹ M☉ yr⁻¹) as well
as those of consistently NUVrJ-selected quiescent galaxies from
EAGLE (Furlong et al. 2015) and SHARK (Lagos et al. 2018).

<figure><img src="image_19.png" /><figcaption>**Fig. B.1.** Galaxy stellar mass functions for total (gray), star-forming (blue), and quiescent (orange) samples from this work compared to simulations. We include those from UNIVERSE MACHINE Behroozi et al. (2019, sSFR < 10⁻¹¹ M☉ yr⁻¹) shown by dotted lines with hatched uncertainty envelopes, as well as consistently NUVrJ-selected samples from EAGLE (Furlong et al. 2015) and SHARK (Lagos et al. 2018) shown by the dashed and dot-dashed lines, respectively. Upper limits for empty bins are shown by the horizontal gray line with an arrow. Mass incomplete measurements are not shown.</figcaption></figure>
**Appendix B: Further comparisons with simulations**

Here we include additional figures comparing the total, star-
forming, and quiescent mass functions (Fig. B.1) and quies-

cent fractions (Fig. B.2) from this work to those of UNIVERSE
MACHINE Behroozi et al. (2019, sSFR < 10⁻¹¹ M☉ yr⁻¹) as well
as those of consistently NUVrJ-selected quiescent galaxies from
EAGLE (Furlong et al. 2015) and SHARK (Lagos et al. 2018).

<figure><img src="image_20.png" /><figcaption>**Fig. B.1.** Galaxy stellar mass functions for total (gray), star-forming (blue), and quiescent (orange) samples from this work compared to simulations. We include those from UNIVERSE MACHINE Behroozi et al. (2019, sSFR < 10⁻¹¹ M☉ yr⁻¹) shown by dotted lines with hatched uncertainty envelopes, as well as consistently NUVrJ-selected samples from EAGLE (Furlong et al. 2015) and SHARK (Lagos et al. 2018) shown by the dashed and dot-dashed lines, respectively. Upper limits for empty bins are shown by the horizontal gray line with an arrow. Mass incomplete measurements are not shown.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Fig. B.2.** Fraction of star-forming (blue) and quiescent (orange) galaxy samples as a function of mass at for three redshift ranges compared to simulations. We include measurements from UNIVERSE MACHINE Behroozi et al. (2019, sSFR < $10^{-11}$ M$_\odot$ yr$^{-1}$) shown by dotted lines, as well as those from consistently *NUVrJ*-selected samples from EAGLE (Furlong et al. 2015) and SHARK (Lagos et al. 2018).</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Fig. B.2.** Fraction of star-forming (blue) and quiescent (orange) galaxy samples as a function of mass at for three redshift ranges compared to simulations. We include measurements from UNIVERSE MACHINE Behroozi et al. (2019, sSFR < $10^{-11}$ M$_{\odot}$ yr$^{-1}$) shown by dotted lines, as well as those from consistently *NUVrJ*-selected samples from EAGLE (Furlong et al. 2015) and SHARK (Lagos et al. 2018).</figcaption></figure>
Appendix C: Fitting Results

Here we include the inferred total, star-forming, and quiescent mass functions (Figs. C.1, C.2, and C.3) optimized by both $\chi^2$ and Likelihood maximization; their derived parameters are contained in Tables C.1, C.2, and C.3, respectively. These results are summarized together in Fig. C.4, along with quiescent frac-

tions as a function of mass (defined three ways). Shown in the same figure are the SMFs realized by evaluation of the continuity model (based on Leja et al. 2019a), with related parameters contained in Table C.4. For those interested in parameter covariances, we have made the full Markov chains available to download<fnref n="15" />.

<figure><img src="image_21.png" /><figcaption>**Fig. C.1.** Results of fitting a double ($z < 3$) and single Schechter ($z \ge 3$) functional forms to the observed galaxy stellar mass function. Measurements are binned in redshift and inferred from both $\chi^2$ minimization (red) and likelihood methods (evolving colors) taken at the median of each parameter posterior distribution. In both cases, Schechter functions are convolved with parameterized, redshift-dependent kernels in $\delta M$ (upper right subpanels) to account for Eddington bias (dashed lines) which is then removed in the corrected fit (solid lines). For completeness, the corrected maximum likelihood solutions are also indicated (dash-dot colored curves). In the lower panels, fitting residuals between the uncorrected kernel convolved fits and the data are shown by both the relative fractional difference (dashed curve) and the difference weighted by the uncertainty (square points). Corresponding parameters are shown in Table C.1.</figcaption></figure>

<footnote n="15">¹⁵ https://doi.org/10.5281/zenodo.7808832</footnote>
Appendix C: Fitting Results

Here we include the inferred total, star-forming, and quiescent mass functions (Figs. C.1, C.2, and C.3) optimized by both $\chi^2$ and Likelihood maximization; their derived parameters are contained in Tables C.1, C.2, and C.3, respectively. These results are summarized together in Fig. C.4, along with quiescent frac-

tions as a function of mass (defined three ways). Shown in the same figure are the SMFs realized by evaluation of the continuity model (based on Leja et al. 2019a), with related parameters contained in Table C.4. For those interested in parameter covariances, we have made the full Markov chains available to download<fnref n="15" />.

<figure><img src="image_22.png" /><figcaption>Fig. C.1. Results of fitting a double ($z < 3$) and single Schechter ($z \ge 3$) functional forms to the observed galaxy stellar mass function. Measurements are binned in redshift and inferred from both $\chi^2$ minimization (red) and likelihood methods (evolving colors) taken at the median of each parameter posterior distribution. In both cases, Schechter functions are convolved with parameterized, redshift-dependent kernels in $\delta M$ (upper right subpanels) to account for Eddington bias (dashed lines) which is then removed in the corrected fit (solid lines). For completeness, the corrected maximum likelihood solutions are also indicated (dash-dot colored curves). In the lower panels, fitting residuals between the uncorrected kernel convolved fits and the data are shown by both the relative fractional difference (dashed curve) and the difference weighted by the uncertainty (square points). Corresponding parameters are shown in Table C.1.</figcaption></figure>

<footnote n="15">¹⁵ https://doi.org/10.5281/zenodo.7808832</footnote>
<figure><img src="image_22.png" /><figcaption>**Fig. C.2.** Results of fitting a double ($z \le 3.0$) and single Schechter ($z > 3.0$) functional forms to the star-forming sample. The figure follows the same layout as Fig. C.1. Corresponding parameters are shown in Table C.2.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. C.3.** Results of fitting a double ($z \le 1.5$) and single Schechter ($z > 1.5$) functional forms to the quiescent sample. The figure follows the same layout as Fig. C.1. Corresponding parameters are shown in Table C.3.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Fig. C.2.** Results of fitting a double ($z \le 3.0$) and single Schechter ($z > 3.0$) functional forms to the star-forming sample. The figure follows the same layout as Fig. C.1. Corresponding parameters are shown in Table C.2.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Fig. C.3.** Results of fitting a double ($z \le 1.5$) and single Schechter ($z > 1.5$) functional forms to the quiescent sample. The figure follows the same layout as Fig. C.1. Corresponding parameters are shown in Table C.3.</figcaption></figure>
<figure><img src="image_24.png" /></figure>

**Fig. C.4.** Best fit Schechter models corresponding to maximum likelihood for the total, star-forming, and quiescent samples (solid colored curves) based on the observed data (colored points). Fits for the median posterior will be similar for symmetric parameter posterior distributions, which are found in all but the last bin shown here. Kernel convolved fits are shown by the dashed curves. Fits measured using the continuity model are shown in gray solid curves. Lower panels indicate the fraction of quiescent galaxies $\mathcal{F}_{\text{QG}} \equiv \Phi_{\text{QG}}/(\Phi_{\text{SF}} + \Phi_{\text{QG}})$ as a function of $\log_{10}(M/M_{\odot})$ for the data (colored points) and maximum likelihood fits (solid curve). Since the star-forming and quiescent fits are not required to sum to the total fit, changing the definition to $\mathcal{F}_{\text{QG}} \equiv \Phi_{\text{QG}}/\Phi_{\text{Total}}$ as reported by Ilbert et al. (2013) and Davidzon et al. (2017) produces the dotted curves with noticeably lower $\mathcal{F}_{\text{QG}}$ at high $M$.
<figure><img src="image_25.png" /></figure>

**Fig. C.4.** Best fit Schechter models corresponding to maximum likelihood for the total, star-forming, and quiescent samples (solid colored curves) based on the observed data (colored points). Fits for the median posterior will be similar for symmetric parameter posterior distributions, which are found in all but the last bin shown here. Kernel convolved fits are shown by the dashed curves. Fits measured using the continuity model are shown in gray solid curves. Lower panels indicate the fraction of quiescent galaxies $\mathcal{F}_{\text{QG}} \equiv \Phi_{\text{QG}}/(\Phi_{\text{SF}} + \Phi_{\text{QG}})$ as a function of $\log_{10}(M/M_{\odot})$ for the data (colored points) and maximum likelihood fits (solid curve). Since the star-forming and quiescent fits are not required to sum to the total fit, changing the definition to $\mathcal{F}_{\text{QG}} \equiv \Phi_{\text{QG}}/\Phi_{\text{Total}}$ as reported by Ilbert et al. (2013) and Davidzon et al. (2017) produces the dotted curves with noticeably lower $\mathcal{F}_{\text{QG}}$ at high $M$.
<table><caption>**Table C.1.** Double ($z \le 3$) and single ($z > 3$) Schechter parameters derived for the total mass complete sample.</caption><thead><tr><th>z-bin</th><th>Log<sub>10</sub> M*<br>(M<sub>☉</sub>)</th><th>α<sub>1</sub></th><th>Φ<sub>1</sub> × 10<sup>-3</sup><br>(Mpc<sup>-3</sup> dex<sup>-1</sup>)</th><th>α<sub>2</sub></th><th>Φ<sub>2</sub> × 10<sup>-3</sup><br>(Mpc<sup>-3</sup> dex<sup>-1</sup>)</th><th>ρ<sub>*</sub>(M ≤ 10<sup>8</sup>) × 10<sup>7</sup><br>(M<sub>☉</sub> Mpc<sup>-3</sup>)</th></tr></thead><tbody><tr><td colspan="7"><strong>Likelihood Fit<sup>a</sup></strong></td></tr><tr><td>0.2 &lt; z ≤ 0.5</td><td>10.89<sup>+0.14</sup><sub>-0.14</sub>[10.92]</td><td>-1.42<sup>+0.05</sup><sub>-0.06</sub>[-1.46]</td><td>0.73<sup>+0.25</sup><sub>-0.27</sub>[0.59]</td><td>-0.46<sup>+0.50</sup><sub>-0.46</sub>[-0.66]</td><td>1.09<sup>+0.50</sup><sub>-0.54</sub>[1.12]</td><td>15.91<sup>+1.50</sup><sub>-1.43</sub>[16.03]</td></tr><tr><td>0.5 &lt; z ≤ 0.8</td><td>10.96<sup>+0.10</sup><sub>-0.10</sub>[10.98]</td><td>-1.39<sup>+0.05</sup><sub>-0.07</sub>[-1.45]</td><td>0.66<sup>+0.22</sup><sub>-0.27</sub>[0.48]</td><td>-0.61<sup>+0.46</sup><sub>-0.39</sub>[-0.79]</td><td>0.83<sup>+0.37</sup><sub>-0.43</sub>[0.95]</td><td>15.49<sup>+1.25</sup><sub>-1.17</sub>[15.44]</td></tr><tr><td>0.8 &lt; z ≤ 1.1</td><td>11.02<sup>+0.08</sup><sub>-0.09</sub>[11.01]</td><td>-1.32<sup>+0.04</sup><sub>-0.06</sub>[-1.35]</td><td>0.84<sup>+0.20</sup><sub>-0.31</sub>[0.71]</td><td>-0.63<sup>+0.48</sup><sub>-0.44</sub>[-0.72]</td><td>0.66<sup>+0.34</sup><sub>-0.42</sub>[0.74]</td><td>17.51<sup>+1.33</sup><sub>-1.27</sub>[16.80]</td></tr><tr><td>1.1 &lt; z ≤ 1.5</td><td>11.00<sup>+0.07</sup><sub>-0.11</sub>[11.00]</td><td>-1.33<sup>+0.04</sup><sub>-0.05</sub>[-1.34]</td><td>0.72<sup>+0.15</sup><sub>-0.23</sub>[0.68]</td><td>-0.51<sup>+0.62</sup><sub>-0.62</sub>[-0.57]</td><td>0.34<sup>+0.30</sup><sub>-0.25</sub>[0.34]</td><td>12.64<sup>+0.86</sup><sub>-0.82</sub>[12.27]</td></tr><tr><td>1.5 &lt; z ≤ 2.0</td><td>10.86<sup>+0.07</sup><sub>-0.08</sub>[10.89]</td><td>-1.48<sup>+0.07</sup><sub>-0.09</sub>[-1.52]</td><td>0.29<sup>+0.11</sup><sub>-0.11</sub>[0.23]</td><td>-0.43<sup>+0.37</sup><sub>-0.31</sub>[-0.56]</td><td>0.64<sup>+0.13</sup><sub>-0.16</sub>[0.65]</td><td>7.62<sup>+0.62</sup><sub>-0.63</sub>[7.63]</td></tr><tr><td>2.0 &lt; z ≤ 2.5</td><td>10.78<sup>+0.16</sup><sub>-0.14</sub>[10.79]</td><td>-1.46<sup>+0.05</sup><sub>-0.06</sub>[-1.47]</td><td>0.27<sup>+0.09</sup><sub>-0.08</sub>[0.26]</td><td>0.07<sup>+0.58</sup><sub>-0.60</sub>[0.07]</td><td>0.27<sup>+0.12</sup><sub>-0.12</sub>[0.31]</td><td>4.25<sup>+0.58</sup><sub>-0.54</sub>[4.52]</td></tr><tr><td>2.5 &lt; z ≤ 3.0</td><td>10.97<sup>+0.06</sup><sub>-0.07</sub>[10.98]</td><td>-1.46</td><td>0.24<sup>+0.03</sup><sub>-0.02</sub>[0.23]</td><td>—</td><td>—</td><td>3.51<sup>+0.33</sup><sub>-0.33</sub>[3.54]</td></tr><tr><td>3.0 &lt; z ≤ 3.5</td><td>10.83<sup>+0.11</sup><sub>-0.09</sub>[10.83]</td><td>-1.46</td><td>0.21<sup>+0.03</sup><sub>-0.03</sub>[0.21]</td><td>—</td><td>—</td><td>2.22<sup>+0.24</sup><sub>-0.28</sub>[2.22]</td></tr><tr><td>3.5 &lt; z ≤ 4.5</td><td>10.46<sup>+0.09</sup><sub>-0.06</sub>[10.44]</td><td>-1.46</td><td>0.20<sup>+0.03</sup><sub>-0.03</sub>[0.20]</td><td>—</td><td>—</td><td>0.87<sup>+0.08</sup><sub>-0.08</sub>[0.87]</td></tr><tr><td>4.5 &lt; z ≤ 5.5</td><td>10.30<sup>+0.10</sup><sub>-0.10</sub>[10.32]</td><td>-1.46</td><td>0.14<sup>+0.04</sup><sub>-0.03</sub>[0.13]</td><td>—</td><td>—</td><td>0.42<sup>+0.07</sup><sub>-0.06</sub>[0.42]</td></tr><tr><td>5.5 &lt; z ≤ 6.5</td><td>10.14<sup>+0.10</sup><sub>-0.12</sub>[10.18]</td><td>-1.46</td><td>0.06<sup>+0.03</sup><sub>-0.02</sub>[0.05]</td><td>—</td><td>—</td><td>0.13<sup>+0.03</sup><sub>-0.03</sub>[0.13]</td></tr><tr><td>6.5 &lt; z ≤ 7.5</td><td>10.18<sup>+0.37</sup><sub>-0.27</sub>[10.49]</td><td>-1.46</td><td>0.03<sup>+0.03</sup><sub>-0.02</sub>[0.01]</td><td>—</td><td>—</td><td>0.06<sup>+0.03</sup><sub>-0.02</sub>[0.06]</td></tr><tr><td colspan="7"><strong>χ<sup>2</sup> Fit<sup>b</sup></strong></td></tr><tr><td>0.2 &lt; z ≤ 0.5</td><td>10.93 ± 0.09(0.14)</td><td>-1.45 ± 0.05(0.08)</td><td>0.60 ± 0.21(0.35)</td><td>-0.66 ± 0.33(0.54)</td><td>1.07 ± 0.32(0.51)</td><td>16.15<sup>+7.52</sup><sub>-5.44</sub></td></tr><tr><td>0.5 &lt; z ≤ 0.8</td><td>10.98 ± 0.08(0.11)</td><td>-1.44 ± 0.08(0.12)</td><td>0.48 ± 0.29(0.41)</td><td>-0.79 ± 0.36(0.50)</td><td>0.92 ± 0.27(0.37)</td><td>16.73<sup>+8.09</sup><sub>-5.98</sub></td></tr><tr><td>0.8 &lt; z ≤ 1.1</td><td>11.01 ± 0.09(0.11)</td><td>-1.35 ± 0.07(0.09)</td><td>0.70 ± 0.35(0.47)</td><td>-0.72 ± 0.45(0.61)</td><td>0.77 ± 0.34(0.46)</td><td>18.89<sup>+9.13</sup><sub>-7.72</sub></td></tr><tr><td>1.1 &lt; z ≤ 1.5</td><td>11.01 ± 0.10(0.13)</td><td>-1.34 ± 0.04(0.06)</td><td>0.69 ± 0.18(0.25)</td><td>-0.56 ± 0.62(0.84)</td><td>0.32 ± 0.24(0.33)</td><td>11.97<sup>+6.51</sup><sub>-4.12</sub></td></tr><tr><td>1.5 &lt; z ≤ 2.0</td><td>10.88 ± 0.08(0.08)</td><td>-1.52 ± 0.10(0.11)</td><td>0.23 ± 0.12(0.13)</td><td>-0.55 ± 0.34(0.37)</td><td>0.67 ± 0.15(0.16)</td><td>7.72<sup>+3.32</sup><sub>-2.82</sub></td></tr><tr><td>2.0 &lt; z ≤ 2.5</td><td>10.80 ± 0.14(0.15)</td><td>-1.46 ± 0.06(0.06)</td><td>0.26 ± 0.09(0.10)</td><td>0.04 ± 0.63(0.66)</td><td>0.29 ± 0.11(0.11)</td><td>4.48<sup>+4.07</sup><sub>-2.04</sub></td></tr><tr><td>2.5 &lt; z ≤ 3.0</td><td>11.00 ± 0.06(0.06)</td><td>-1.46</td><td>0.22 ± 0.02(0.02)</td><td>—</td><td>—</td><td>3.69<sup>+0.50</sup><sub>-0.47</sub></td></tr><tr><td>3.0 &lt; z ≤ 3.5</td><td>10.85 ± 0.14(0.12)</td><td>-1.46</td><td>0.20 ± 0.04(0.04)</td><td>—</td><td>—</td><td>2.36<sup>+1.20</sup><sub>-0.95</sub></td></tr><tr><td>3.5 &lt; z ≤ 4.5</td><td>10.46 ± 0.18(0.10)</td><td>-1.46</td><td>0.20 ± 0.07(0.04)</td><td>—</td><td>—</td><td>0.81<sup>+0.50</sup><sub>-0.36</sub></td></tr><tr><td>4.5 &lt; z ≤ 5.5</td><td>10.33 ± 0.12(0.10)</td><td>-1.46</td><td>0.13 ± 0.04(0.03)</td><td>—</td><td>—</td><td>0.46<sup>+0.18</sup><sub>-0.16</sub></td></tr><tr><td>5.5 &lt; z ≤ 6.5</td><td>10.19 ± 0.05(0.09)</td><td>-1.46</td><td>0.05 ± 0.01(0.02)</td><td>—</td><td>—</td><td>0.13<sup>+0.03</sup><sub>-0.03</sub></td></tr><tr><td>6.5 &lt; z ≤ 7.5</td><td>10.50 ± 0.14(0.36)</td><td>-1.46</td><td>0.01 ± 0.00(0.01)</td><td>—</td><td>—</td><td>0.05<sup>+0.05</sup><sub>-0.02</sub></td></tr></tbody></table>

<footnote><sup>a</sup> For the Likelihood fit, values are shown for the median posterior distributions with even-tailed 68% range and the values corresponding to the maximum likelihood solution in brackets. <sup>b</sup> For the $\chi^2$ regression fit, uncertainties on parameter values are shown multiplied by $\sqrt{\chi_N^2}$ with formal uncertainties following in brackets.</footnote>