COSMOS2020: The galaxy stellar mass function. The assembly and star formation cessation of galaxies at $0.2< z \le 7.5$

J. R. Weaver, I. Davidzon, S. Toft, O. Ilbert, H. J. Mccracken, K. M. L.

Gould, C. K. Jespersen, C. Steinhardt, C. D. P. Lagos, P. L. Capak, et al.

► To cite this version:

J. R. Weaver, I. Davidzon, S. Toft, O. Ilbert, H. J. Mccracken, et al., COSMOS2020: The galaxy stellar mass function. The assembly and star formation cessation of galaxies at $0.2< z \le 7.5$. Astronomy & Astrophysics, 2023, 677, 10.1051/0004-6361/202245581. insu-04505563

HAL Id: insu-04505563

https://insu.hal.science/insu-04505563

Submitted on 15 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# COSMOS2020: The galaxy stellar mass function

## The assembly and star formation cessation of galaxies at $0.2 < z \le 7.5^*$

J. R. Weaver¹,²,³,*, I. Davidzon¹,², S. Toft¹,², O. Ilbert⁴, H. J. McCracken⁵, K. M. L. Gould¹,², C. K. Jespersen⁶, C. Steinhardt¹,², C. D. P. Lagos⁷,⁸,¹, P. L. Capak¹, C. M. Casey⁹, N. Chartab¹⁰, A. L. Faisst¹¹, C. C. Hayward¹², J. S. Kartaltepe¹³, O. B. Kauffmann⁴, A. M. Koekemoer¹⁴, V. Kokorev¹,², C. Laigle⁵, D. Liu¹⁵, A. Long⁹,**, G. E. Magdis¹,¹⁶,², C. J. R. McPartland¹,², B. Milvang-Jensen¹,², B. Mobasher¹⁷, A. Moneti⁵, Y. Peng¹⁸,¹⁹, D. B. Sanders²⁰, M. Shuntov⁵, A. Sneppen¹,², F. Valentino¹,², L. Zalesky²⁰, and G. Zamorani²¹

*(Affiliations can be found after the references)*

Received 29 November 2022 / Accepted 10 April 2023

## ABSTRACT

**Context.** How galaxies form, assemble, and cease their star formation is a central question within the modern landscape of galaxy evolution studies. These processes are indelibly imprinted on the galaxy stellar mass function (SMF), and its measurement and understanding is key to uncovering a unified theory of galaxy evolution.

**Aims.** We present constraints on the shape and evolution of the galaxy SMF, the quiescent galaxy fraction, and the cosmic stellar mass density across 90% of the history of the Universe from $z = 7.5 \rightarrow 0.2$ as a means to study the physical processes that underpin galaxy evolution.

**Methods.** The COSMOS survey is an ideal laboratory for studying representative galaxy samples. Now equipped with deeper and more homogeneous near-infrared coverage exploited by the COSMOS2020 catalog, we leverage the large 1.27 deg² effective area to improve sample statistics and understand spatial variations (cosmic variance) – particularly for rare, massive galaxies – and push to higher redshifts with greater confidence and mass completeness than previous studies. We divide the total stellar mass function into star-forming and quiescent subsamples through *NUVrJ* color-color selection. The measurements are then fit with single- and double-component Schechter functions to infer the intrinsic galaxy stellar mass function, the evolution of its key parameters, and the cosmic stellar mass density out to $z = 7.5$. Finally, we compare our measurements to predictions from state-of-the-art cosmological simulations and theoretical dark matter halo mass functions.

**Results.** We find a smooth, monotonic evolution in the galaxy stellar mass function since $z = 7.5$, in general agreement with previous studies. The number density of star-forming systems have undergone remarkably consistent growth spanning four decades in stellar mass from $z = 7.5 \rightarrow 2$ whereupon high-mass systems become predominantly quiescent ("downsizing"). Meanwhile, the assembly and growth of low-mass quiescent systems only occurred recently, and rapidly. An excess of massive systems at $z \approx 2.5-5.5$ with strikingly red colors, with some being newly identified, increase the observed number densities to the point where the SMF cannot be reconciled with a Schechter function.

**Conclusions.** Systematics including cosmic variance and/or active galactic nuclei contamination are unlikely to fully explain this excess, and so we speculate that they may be dust-obscured populations similar to those found in far infrared surveys. Furthermore, we find a sustained agreement from $z \approx 3-6$ between the stellar and dark matter halo mass functions for the most massive systems, suggesting that star formation in massive halos may be more efficient at early times.

**Key words.** galaxies: evolution – galaxies: statistics – galaxies: luminosity function, mass function – galaxies: high-redshift

## 1. Introduction

The galaxy stellar mass function (SMF) is defined as the number density of galaxies $\Phi(M, z)$ in bins of stellar mass $\Delta M$ at each redshift $z$, and is a fundamental cosmological observable in the study of the statistical properties of galaxies. Understanding its shape and evolution ultimately informs us about the growth of the baryonic content of the Universe, and can help infer the star formation activity and the availability of cold molecular gas across cosmic time. Its integral over $M$ is the galaxy stellar mass density (SMD) $\rho_*(z)$, which describes the cumulative mass assembled by a given epoch.

Remarkable progress has been made since the first mass-selected measurements of the local $z \approx 0$ SMF by Cole et al. (2001). The intervening years have been marked by order-of-magnitude increases in sample size, photometric preci-

sion, and redshift accuracy which have enabled more precise determinations of the local SMF as well as ground-breaking extensions to increasingly higher redshifts (e.g., Fontana et al. 2006; Marchesini et al. 2009; Drory et al. 2009; Peng et al. 2010; Ilbert et al. 2010, 2013; Pozzetti et al. 2010; Mortlock et al. 2011; Vulcani et al. 2013; Muzzin et al. 2013; Grazian et al. 2015; Mortlock et al. 2015; Duncan et al. 2014; Song et al. 2016; Bernardi et al. 2017; Davidzon et al. 2017; Wright et al. 2018; Santini et al. 2021, 2022; Stefanon et al. 2021; Adams et al. 2021; Thorne et al. 2021; McLeod et al. 2021).

A coherent picture has emerged from these studies. For instance, the shape of the SMF for star-forming galaxies at $z < 3$ is well described by the empirically constructed Schechter function (Schechter 1976), which features an exponential downturn at a characteristic stellar mass $M^*$, with a low-mass end that declines with a slope $\alpha$; both seem to remain constant out to $z \approx 2$ (Marchesini et al. 2009; Ilbert et al. 2013; Tomczak et al. 2014; Davidzon et al. 2017; Adams et al. 2021). This constancy

<footnote>* Data files containing sample IDs and key measurements are available for download: https://doi.org/10.5281/zenodo.7808832</footnote>

<footnote>** NASA Hubble fellow.</footnote>
# COSMOS2020: The galaxy stellar mass function

## The assembly and star formation cessation of galaxies at $0.2 < z \le 7.5^*$

J. R. Weaver¹,²,³,*, I. Davidzon¹,², S. Toft¹,², O. Ilbert⁴, H. J. McCracken⁵, K. M. L. Gould¹,², C. K. Jespersen⁶, C. Steinhardt¹,², C. D. P. Lagos⁷,⁸,¹, P. L. Capak¹, C. M. Casey⁹, N. Chartab¹⁰, A. L. Faisst¹¹, C. C. Hayward¹², J. S. Kartaltepe¹³, O. B. Kauffmann⁴, A. M. Koekemoer¹⁴, V. Kokorev¹,², C. Laigle⁵, D. Liu¹⁵, A. Long⁹,**, G. E. Magdis¹,¹⁶,², C. J. R. McPartland¹,², B. Milvang-Jensen¹,², B. Mobasher¹⁷, A. Moneti⁵, Y. Peng¹⁸,¹⁹, D. B. Sanders²⁰, M. Shuntov⁵, A. Sneppen¹,², F. Valentino¹,², L. Zalesky²⁰, and G. Zamorani²¹

*(Affiliations can be found after the references)*

Received 29 November 2022 / Accepted 10 April 2023

## ABSTRACT

**Context.** How galaxies form, assemble, and cease their star formation is a central question within the modern landscape of galaxy evolution studies. These processes are indelibly imprinted on the galaxy stellar mass function (SMF), and its measurement and understanding is key to uncovering a unified theory of galaxy evolution.

**Aims.** We present constraints on the shape and evolution of the galaxy SMF, the quiescent galaxy fraction, and the cosmic stellar mass density across 90% of the history of the Universe from $z = 7.5 \rightarrow 0.2$ as a means to study the physical processes that underpin galaxy evolution.

**Methods.** The COSMOS survey is an ideal laboratory for studying representative galaxy samples. Now equipped with deeper and more homogeneous near-infrared coverage exploited by the COSMOS2020 catalog, we leverage the large 1.27 deg² effective area to improve sample statistics and understand spatial variations (cosmic variance) – particularly for rare, massive galaxies – and push to higher redshifts with greater confidence and mass completeness than previous studies. We divide the total stellar mass function into star-forming and quiescent subsamples through *NUVrJ* color-color selection. The measurements are then fit with single- and double-component Schechter functions to infer the intrinsic galaxy stellar mass function, the evolution of its key parameters, and the cosmic stellar mass density out to $z = 7.5$. Finally, we compare our measurements to predictions from state-of-the-art cosmological simulations and theoretical dark matter halo mass functions.

**Results.** We find a smooth, monotonic evolution in the galaxy stellar mass function since $z = 7.5$, in general agreement with previous studies. The number density of star-forming systems have undergone remarkably consistent growth spanning four decades in stellar mass from $z = 7.5 \rightarrow 2$ whereupon high-mass systems become predominantly quiescent ("downsizing"). Meanwhile, the assembly and growth of low-mass quiescent systems only occurred recently, and rapidly. An excess of massive systems at $z \approx 2.5-5.5$ with strikingly red colors, with some being newly identified, increase the observed number densities to the point where the SMF cannot be reconciled with a Schechter function.

**Conclusions.** Systematics including cosmic variance and/or active galactic nuclei contamination are unlikely to fully explain this excess, and so we speculate that they may be dust-obscured populations similar to those found in far infrared surveys. Furthermore, we find a sustained agreement from $z \approx 3-6$ between the stellar and dark matter halo mass functions for the most massive systems, suggesting that star formation in massive halos may be more efficient at early times.

**Key words.** galaxies: evolution – galaxies: statistics – galaxies: luminosity function, mass function – galaxies: high-redshift

## 1. Introduction

The galaxy stellar mass function (SMF) is defined as the number density of galaxies $\Phi(M, z)$ in bins of stellar mass $\Delta M$ at each redshift $z$, and is a fundamental cosmological observable in the study of the statistical properties of galaxies. Understanding its shape and evolution ultimately informs us about the growth of the baryonic content of the Universe, and can help infer the star formation activity and the availability of cold molecular gas across cosmic time. Its integral over $M$ is the galaxy stellar mass density (SMD) $\rho_*(z)$, which describes the cumulative mass assembled by a given epoch.

Remarkable progress has been made since the first mass-selected measurements of the local $z \approx 0$ SMF by Cole et al. (2001). The intervening years have been marked by order-of-magnitude increases in sample size, photometric preci-

sion, and redshift accuracy which have enabled more precise determinations of the local SMF as well as ground-breaking extensions to increasingly higher redshifts (e.g., Fontana et al. 2006; Marchesini et al. 2009; Drory et al. 2009; Peng et al. 2010; Ilbert et al. 2010, 2013; Pozzetti et al. 2010; Mortlock et al. 2011; Vulcani et al. 2013; Muzzin et al. 2013; Grazian et al. 2015; Mortlock et al. 2015; Duncan et al. 2014; Song et al. 2016; Bernardi et al. 2017; Davidzon et al. 2017; Wright et al. 2018; Santini et al. 2021, 2022; Stefanon et al. 2021; Adams et al. 2021; Thorne et al. 2021; McLeod et al. 2021).

A coherent picture has emerged from these studies. For instance, the shape of the SMF for star-forming galaxies at $z < 3$ is well described by the empirically constructed Schechter function (Schechter 1976), which features an exponential downturn at a characteristic stellar mass $M^*$, with a low-mass end that declines with a slope $\alpha$; both seem to remain constant out to $z \approx 2$ (Marchesini et al. 2009; Ilbert et al. 2013; Tomczak et al. 2014; Davidzon et al. 2017; Adams et al. 2021). This constancy

<footnote>* Data files containing sample IDs and key measurements are available for download: https://doi.org/10.5281/zenodo.7808832</footnote>

<footnote>** NASA Hubble fellow.</footnote>
in the shape of the SMF implies that star formation activity has acted consistently to transform baryons from cold gas into stars, and appears to have formed >75% of the total stellar mass of the Universe in ≈10 Gyr. Only the normalization of the SMF, $\Phi^*(z)$, of star-forming galaxies has been confidently shown to evolve with cosmic time, and its rapid evolution at earlier times points directly to enhanced rates of mass growth at earlier epochs (Popesso et al. 2023, and references therein).

Which physical processes are responsible for the behavior of the SMF, and the extent of their influence at different cosmic epochs, are not well understood even in the low redshift Universe ($z < 1$). Examples include feedback from supermassive black holes, supernovae, galaxy-galaxy mergers, stellar and gas dynamics, and gaseous inflows and outflows, all of which act on varying scales both physically and in time. As a result, significant effort has been undertaken to build detailed cosmological simulations to identify and understand the role of physical processes that underpin the observed shape and evolution of the SMF (e.g., Somerville et al. 2015; Furlong et al. 2015; Lagos et al. 2018; Laigle et al. 2019; Pillepich et al. 2018; Davé et al. 2019; Lovell et al. 2021). Hence, precise measurements of the SMF and SMD are key observables utilized by most (but not all, see e.g., Dubois et al. 2014) large-scale cosmological simulations to tune input parameters such that the SMF of the local Universe ($z \approx 0$) is recovered. Comparisons of predicted and observed SMFs at earlier cosmic ages can point to the physical processes at play and/or suggest recipes in need of refinement (Torrey et al. 2014; Furlong et al. 2015).

Owing to the challenges associated with large-scale simulations, purely analytical, data-driven models have enjoyed great popularity with the introduction of the first large galaxy samples with photometric redshifts (photo-$z$). For example, Peng et al. (2010) constructed a phenomenological model to explain the bimodality of galaxy types (star-forming and quiescent) as seen from their mass functions, megaparsec-scale environment, and star formation activity. That is, as a consequence of two mechanisms driving star formation cessation (often referred to as “quenching”): massive galaxies cease forming stars irrespective of environment (“mass quenching”), and galaxies in dense environments cease forming stars irrespective of mass (“environmental quenching”). Several discrete mechanisms have been proposed to explain the latter, such as ram pressure stripping (e.g., Gunn & Gott 1972), gas strangulation (e.g., Larson et al. 1980; Balogh et al. 2000), and dynamical heating of gas within haloes (e.g., Gabor & Davé 2015). Proposed mechanisms for the former must, by definition, involve secular processes. This includes star formation cessation due to structural changes, or heating and/or ejection of gas by central super-massive black holes for the most massive galaxies – radiative feedback from active galactic nuclei (AGN) – or by supernovae for less massive systems as they are more weakly self-gravitating (e.g., Shankar et al. 2006). Peng et al. (2010) hypothesized that while environmental effects reproduce the single Schechter shape observed for star-forming, blue galaxies in the local Universe, it is through a combination of both environmental and mass effects that the two-component Schechter shape observed for quiescent, red galaxies is reproduced. A wide variety of extensions have been applied to this model to directly incorporate other measurements such as the gas fraction (Bouché et al. 2010), and wholly new models continue to be developed (e.g., Peng et al. 2015; Belli et al. 2019; Suess et al. 2021; Varma et al. 2022).

The observed distributions of stellar mass and redshift are not in themselves the intrinsic, underlying SMF. Instead, they are bridged by statistical inferences between carefully selected

samples and a vast, practically unobservable parent population from which a sample is taken. Insufficient control of biases and systematics obscure such inferences by creating unrepresentative samples, and hence weaken comparisons with simulation and analytical models (Marchesini et al. 2009; Fontanot et al. 2009). Great effort has been expended over the past 20 years to strip away a number of these biases, and their obscuring effects (Cole et al. 2001; Pozzetti et al. 2010; Marchesini et al. 2009; Bernardi et al. 2017; Davidzon et al. 2017; Leja et al. 2019a; Adams et al. 2021). Examples of these biases include sample incompleteness (including Malmquist biases and mass completeness), the so-called “cosmic” variance relating to sampling over- and under-dense regions of large scale structure, and Eddington bias which acts to overestimate the number density of the most massive galaxies (Malmquist 1920, 1922; Eddington 1913). While many studies incorporate only poisson noise (Fontanot et al. 2009, although this is steadily improving with time), other notable uncertainties exist including sample variance and effective volume size, as well as uncertainties on photo-$z$ and stellar mass; the latter two items can introduce significant bias. Since the ultimate goal of any survey is to generalize the observed properties of a specific sample to that of the entire population, ignoring any of these important biases or uncertainties severely complicate this effort.

Studies of the SMF at increasingly higher redshift have been made possible due to advances in facilities and continued investment in deep, primarily photometric surveys of the distant Universe. Building on work by Songaila et al. (1994) and Glazebrook et al. (1995), Cowie et al. (1996) secured the first mass-complete samples at $z \approx 2-3$, selected in the near-infrared ($K$) to directly constrain the Balmer continuum at $\sim 0.6$ µm. Sampled nearer the stellar bulk at $\sim 1-2$ µm, their mass-selected samples enabled a higher degree of mass completeness, and in turn provided a more representative view of the high-$z$ Universe compared to optically selected forerunner surveys. More recently, precise mass estimates from similarly selected samples have been obtained by fitting observed spectral energy distributions (Tomczak et al. 2014; Martis et al. 2016; Straatman et al. 2015), from the ground (with VISTA, UKIRT) and space (HST/WFC3, *Spitzer*/IRAC). Although limited in area, samples recovered by HST have continued to pose new and significant challenges to existing paradigms by highlighting a series of stark changes in the shape of the SMF that indicates earlier galaxy populations were fundamentally different from those in the present day Universe. Additional studies at these early times ($z \ge 2$) by degree-scale surveys capable of finding the rarest sources have revealed the existence of surprisingly mature, massive quiescent galaxies whose mass has already been assembled by $z \approx 4$ (e.g., Ilbert et al. 2013), challenging the assumed timeline typical for galaxy assembly (Steinhardt et al. 2016; Behroozi & Silk 2018). Limited numbers of them have been confirmed by detailed spectroscopic follow-up (e.g., Glazebrook et al. 2017; Schreiber et al. 2018a; Tanaka et al. 2019; Valentino et al. 2020; Forrest et al. 2020a,b). Likewise, several studies have placed constraints on the massive end of the SMF at the highest-redshifts ($z > 6$) from samples selected by their rest-frame UV luminosity (e.g., Song et al. 2016; Harikane et al. 2016), who utilized empirically calibrated $L_{UV}-M$ relations to convert directly measured UV luminosity into indirect estimates of stellar mass.

Although tantalizing, deriving mass estimates from UV luminosity involves significant uncertainties; only deep infrared photometry can directly measure the rest-frame stellar bulk ($\lambda > 4000$ Å, although ideally $\sim 1-2$ µm) required to more directly
in the shape of the SMF implies that star formation activity has acted consistently to transform baryons from cold gas into stars, and appears to have formed >75% of the total stellar mass of the Universe in ≈10 Gyr. Only the normalization of the SMF, $\Phi^*(z)$, of star-forming galaxies has been confidently shown to evolve with cosmic time, and its rapid evolution at earlier times points directly to enhanced rates of mass growth at earlier epochs (Popesso et al. 2023, and references therein).

Which physical processes are responsible for the behavior of the SMF, and the extent of their influence at different cosmic epochs, are not well understood even in the low redshift Universe ($z < 1$). Examples include feedback from supermassive black holes, supernovae, galaxy-galaxy mergers, stellar and gas dynamics, and gaseous inflows and outflows, all of which act on varying scales both physically and in time. As a result, significant effort has been undertaken to build detailed cosmological simulations to identify and understand the role of physical processes that underpin the observed shape and evolution of the SMF (e.g., Somerville et al. 2015; Furlong et al. 2015; Lagos et al. 2018; Laigle et al. 2019; Pillepich et al. 2018; Davé et al. 2019; Lovell et al. 2021). Hence, precise measurements of the SMF and SMD are key observables utilized by most (but not all, see e.g., Dubois et al. 2014) large-scale cosmological simulations to tune input parameters such that the SMF of the local Universe ($z \approx 0$) is recovered. Comparisons of predicted and observed SMFs at earlier cosmic ages can point to the physical processes at play and/or suggest recipes in need of refinement (Torrey et al. 2014; Furlong et al. 2015).

Owing to the challenges associated with large-scale simulations, purely analytical, data-driven models have enjoyed great popularity with the introduction of the first large galaxy samples with photometric redshifts (photo-$z$). For example, Peng et al. (2010) constructed a phenomenological model to explain the bimodality of galaxy types (star-forming and quiescent) as seen from their mass functions, megaparsec-scale environment, and star formation activity. That is, as a consequence of two mechanisms driving star formation cessation (often referred to as “quenching”): massive galaxies cease forming stars irrespective of environment (“mass quenching”), and galaxies in dense environments cease forming stars irrespective of mass (“environmental quenching”). Several discrete mechanisms have been proposed to explain the latter, such as ram pressure stripping (e.g., Gunn & Gott 1972), gas strangulation (e.g., Larson et al. 1980; Balogh et al. 2000), and dynamical heating of gas within haloes (e.g., Gabor & Davé 2015). Proposed mechanisms for the former must, by definition, involve secular processes. This includes star formation cessation due to structural changes, or heating and/or ejection of gas by central super-massive black holes for the most massive galaxies – radiative feedback from active galactic nuclei (AGN) – or by supernovae for less massive systems as they are more weakly self-gravitating (e.g., Shankar et al. 2006). Peng et al. (2010) hypothesized that while environmental effects reproduce the single Schechter shape observed for star-forming, blue galaxies in the local Universe, it is through a combination of both environmental and mass effects that the two-component Schechter shape observed for quiescent, red galaxies is reproduced. A wide variety of extensions have been applied to this model to directly incorporate other measurements such as the gas fraction (Bouché et al. 2010), and wholly new models continue to be developed (e.g., Peng et al. 2015; Belli et al. 2019; Suess et al. 2021; Varma et al. 2022).

The observed distributions of stellar mass and redshift are not in themselves the intrinsic, underlying SMF. Instead, they are bridged by statistical inferences between carefully selected

samples and a vast, practically unobservable parent population from which a sample is taken. Insufficient control of biases and systematics obscure such inferences by creating unrepresentative samples, and hence weaken comparisons with simulation and analytical models (Marchesini et al. 2009; Fontanot et al. 2009). Great effort has been expended over the past 20 years to strip away a number of these biases, and their obscuring effects (Cole et al. 2001; Pozzetti et al. 2010; Marchesini et al. 2009; Bernardi et al. 2017; Davidzon et al. 2017; Leja et al. 2019a; Adams et al. 2021). Examples of these biases include sample incompleteness (including Malmquist biases and mass completeness), the so-called “cosmic” variance relating to sampling over- and under-dense regions of large scale structure, and Eddington bias which acts to overestimate the number density of the most massive galaxies (Malmquist 1920, 1922; Eddington 1913). While many studies incorporate only poisson noise (Fontanot et al. 2009, although this is steadily improving with time), other notable uncertainties exist including sample variance and effective volume size, as well as uncertainties on photo-$z$ and stellar mass; the latter two items can introduce significant bias. Since the ultimate goal of any survey is to generalize the observed properties of a specific sample to that of the entire population, ignoring any of these important biases or uncertainties severely complicate this effort.

Studies of the SMF at increasingly higher redshift have been made possible due to advances in facilities and continued investment in deep, primarily photometric surveys of the distant Universe. Building on work by Songaila et al. (1994) and Glazebrook et al. (1995), Cowie et al. (1996) secured the first mass-complete samples at $z \approx 2-3$, selected in the near-infrared ($K$) to directly constrain the Balmer continuum at $\sim 0.6$ µm. Sampled nearer the stellar bulk at $\sim 1-2$ µm, their mass-selected samples enabled a higher degree of mass completeness, and in turn provided a more representative view of the high-$z$ Universe compared to optically selected forerunner surveys. More recently, precise mass estimates from similarly selected samples have been obtained by fitting observed spectral energy distributions (Tomczak et al. 2014; Martis et al. 2016; Straatman et al. 2015), from the ground (with VISTA, UKIRT) and space (HST/WFC3, *Spitzer*/IRAC). Although limited in area, samples recovered by HST have continued to pose new and significant challenges to existing paradigms by highlighting a series of stark changes in the shape of the SMF that indicates earlier galaxy populations were fundamentally different from those in the present day Universe. Additional studies at these early times ($z \gtrsim 2$) by degree-scale surveys capable of finding the rarest sources have revealed the existence of surprisingly mature, massive quiescent galaxies whose mass has already been assembled by $z \approx 4$ (e.g., Ilbert et al. 2013), challenging the assumed timeline typical for galaxy assembly (Steinhardt et al. 2016; Behroozi & Silk 2018). Limited numbers of them have been confirmed by detailed spectroscopic follow-up (e.g., Glazebrook et al. 2017; Schreiber et al. 2018a; Tanaka et al. 2019; Valentino et al. 2020; Forrest et al. 2020a,b). Likewise, several studies have placed constraints on the massive end of the SMF at the highest-redshifts ($z > 6$) from samples selected by their rest-frame UV luminosity (e.g., Song et al. 2016; Harikane et al. 2016), who utilized empirically calibrated $L_{UV}-M$ relations to convert directly measured UV luminosity into indirect estimates of stellar mass.

Although tantalizing, deriving mass estimates from UV luminosity involves significant uncertainties; only deep infrared photometry can directly measure the rest-frame stellar bulk ($\lambda > 4000$ Å, although ideally $\sim 1-2$ µm) required to more directly
assess stellar mass at these redshifts. Only recently have deep ($K_s \approx 25$) near-infrared photometric studies enabled the first mass-selected samples at $z > 2$, although with mass uncertainties increasing with redshift (Retzlaff et al. 2010; Fontana et al. 2014; Laigle et al. 2016). Importantly, photometric surveys of galaxies have – and are expected to remain – the primary means by which these measurements will be made; obtaining even elementary parameters for large ($N > 100 000$) samples with deep spectroscopy, while precise, is simply too expensive even when utilizing highly multiplexed instrumentation. For the time being, spectroscopy of the distant Universe will remain a follow-up exercise to strengthen larger photometrically selected samples.

Of these deep photometric surveys, those with large areas have a key advantage: they probe a wider range of environments (i.e., density) compared to narrower surveys. As such, they provide more representative samples that are significantly less likely to be affected by field-to-field cosmic variances (and, for the same reason, are more suited to find rare, massive galaxies). Although one may resort to combining disparate data sets into a single analysis to combat this (e.g., Moster et al. 2013; Henriques et al. 2015; Volonteri & Reines 2016; McLeod et al. 2021), the unknown systematics between survey selection functions and depths make the interpretation of these joined samples more uncertain.

Stitching together surveys of low- and high-$z$ samples carries similar concerns. A lack of uniform selection owing to different survey areas, detection bands, and photo-$z$ determination strategies (i.e., dropout selection) can likewise complicate interpretations arising from such composite samples (e.g., Leja et al. 2019a; McLeod et al. 2021; Adams et al. 2021; Santini et al. 2022). Worse, these photometric samples may have been processed with different spectral energy distribution (SED) fitting codes that assume different templates, emission line recipes, and dust attenuation laws. Differing choices of cosmology and initial mass function, while reversible, nonetheless add complexity.

Accurate estimates of the galaxy stellar mass function at increasingly higher redshifts requires complementary deep observations from near-infrared selected samples to ensure both mass completeness and data reliability. Currently the widest field with near-infrared coverage of the requisite depth to probe large samples of $z \gg 3$ galaxies is the Cosmic Evolution Survey (COSMOS; Scoville et al. 2007). This work takes advantage of the latest NIR-selected catalog of the COSMOS field, COSMOS2020 (Weaver et al. 2022). We adopt the “THE FARMER” catalog whose photometry was consistently measured across all bands by fitting galaxy profiles while taking into account PSF variations, depth, and source crowding. From the resulting photo-$z$ and stellar masses we construct a consistently measured SMF from $z = 0.2–7.5$, identify quiescent and star-forming systems, and study the build-up and assembly of stellar mass over 10 Gyr of cosmic history. This includes a detailed study of key moments in the development of galaxy populations from the Era of Reionization ($z > 6$) to the peak of star formation activity at Cosmic Noon ($z \sim 2$), up to the rich bimodality of star-forming and quiescent galaxies observed at more recent times ($z \le 1$).

This paper is organized as follows. Section 2 introduces the dataset chosen for this analysis from which samples are drawn and possible uncertainties discussed in Sect. 3. Section 4 provides a brief overview of the Schechter function. Results are presented in Sect. 5, including the presentation and fitting of the total, star-forming, and quiescent mass functions. These are compared to literature measurements in Sect. 6, whereupon fur-

ther discussion is had with regards to galaxy assembly, star formation cessation, and connections to dark matter halos. This work concludes in Sect. 7.

These results are computed adopting a standard $\Lambda$CDM cosmology with $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega_{m,0} = 0.3$, and $\Omega_{\Lambda,0} = 0.7$ throughout, such that the dimensionless Hubble parameter $h_{70} \equiv H_0/(70 \text{ km s}^{-1} \text{ Mpc}^{-1}) = 1$. Galaxy stellar masses, when derived from SED fitting, scale as the square of the luminosity distance (i.e., $D_L^2$); hence a factor of $h_{70}^{-1}$ is retained implicitly for all relevant measurements, unless explicitly noted otherwise (see Croton 2013 for an overview of $h$ and best practices). Estimates of stellar mass (henceforth $M$) assume an Chabrier (2003) initial mass function (IMF). All magnitudes are expressed in the AB system (Oke 1974), for which a flux $f_v$ in $\mu$Jy ($10^{-29} \text{ erg cm}^{-2} \text{ s}^{-1} \text{ Hz}^{-1}$) corresponds to $AB_v = 23.9–2.5 \log_{10}(f_v/\mu$Jy).

## 2. Data: COSMOS2020

We utilize the most recent release of the COSMOS catalog, COSMOS2020 (Weaver et al. 2022), comprised of ~100000 galaxies selected from a $izYJHK_s$ coadd with near-infrared depths approaching 26 AB ensuring a mass-selected sample complete down to $10^9 M_{\odot}$ at $z \approx 3$. This sample is complemented by extensive supporting photometry ranging from the UV to 8 $\mu$m over a total area of $2 \text{ deg}^2$, making it the widest near-infrared-selected multiwavelength catalog to these depths. This is made possible by the latest release of the UltraVISTA survey (DR4; McCracken et al. 2012; Moneti et al. 2023), which is the longest running near-infrared survey to date, and complemented by even deeper optical *grizy* data provided by Subaru's Hyper Suprime-Cam instrument (HSC PDR2; Aihara et al. 2019). Flanking this core imaging are $u$ band measurements made by the CLAUDS survey from the Canada-France-Hawaii Telescope (Sawicki et al. 2019) and deep *Spitzer*/IRAC imaging in all four channels, reprocessed to include nearly every exposure taken in COSMOS for use in this catalog (Euclid Collaboration 2022, and references therein). As in the previous COSMOS catalog by Laigle et al. (2016), several intermediate and narrow bands from both Subaru/Suprime-Cam and VISTA (Milvang-Jensen et al. 2013) are used to provide precise determinations of photometric redshifts.

In addition to the comprehensive set of data, photometry is performed by two methods to produce two separate catalogs. The catalog called CLASSIC contains aperture photometry extracted with Source Extractor (Bertin & Arnouts 1996) for the optical/near-infrared and with IRACLEAN (Hsieh et al. 2012) for the infrared (as in Laigle et al. 2016). The catalog called THE FARMER contains profile-fitting photometry extracted with The Farmer (Weaver et al. 2023), which uses The Tractor (Lang et al. 2016) to construct and apply parametric models to estimate source fluxes and does so consistently across more than 30 bands, including all four IRAC channels. Each photometric catalog is paired with photometric redshift and physical parameter estimates from both LePhare (Arnouts et al. 2002; Ilbert et al. 2006) and EAzY (Brammer et al. 2008) resulting in four measurements of photo-$z$ and $M$ per source. The combination of comprehensive deep images, independent photometry techniques, and independent SED fitting enables an unparalleled control of systematics and uncertainties, and hence further insight when inferring the underlying true nature of galaxy populations. See Weaver et al. (2022) for details.
assess stellar mass at these redshifts. Only recently have deep ($K_s \approx 25$) near-infrared photometric studies enabled the first mass-selected samples at $z > 2$, although with mass uncertainties increasing with redshift (Retzlaff et al. 2010; Fontana et al. 2014; Laigle et al. 2016). Importantly, photometric surveys of galaxies have – and are expected to remain – the primary means by which these measurements will be made; obtaining even elementary parameters for large ($N > 100 000$) samples with deep spectroscopy, while precise, is simply too expensive even when utilizing highly multiplexed instrumentation. For the time being, spectroscopy of the distant Universe will remain a follow-up exercise to strengthen larger photometrically selected samples.

Of these deep photometric surveys, those with large areas have a key advantage: they probe a wider range of environments (i.e., density) compared to narrower surveys. As such, they provide more representative samples that are significantly less likely to be affected by field-to-field cosmic variances (and, for the same reason, are more suited to find rare, massive galaxies). Although one may resort to combining disparate data sets into a single analysis to combat this (e.g., Moster et al. 2013; Henriques et al. 2015; Volonteri & Reines 2016; McLeod et al. 2021), the unknown systematics between survey selection functions and depths make the interpretation of these joined samples more uncertain.

Stitching together surveys of low- and high-$z$ samples carries similar concerns. A lack of uniform selection owing to different survey areas, detection bands, and photo-$z$ determination strategies (i.e., dropout selection) can likewise complicate interpretations arising from such composite samples (e.g., Leja et al. 2019a; McLeod et al. 2021; Adams et al. 2021; Santini et al. 2022). Worse, these photometric samples may have been processed with different spectral energy distribution (SED) fitting codes that assume different templates, emission line recipes, and dust attenuation laws. Differing choices of cosmology and initial mass function, while reversible, nonetheless add complexity.

Accurate estimates of the galaxy stellar mass function at increasingly higher redshifts requires complementary deep observations from near-infrared selected samples to ensure both mass completeness and data reliability. Currently the widest field with near-infrared coverage of the requisite depth to probe large samples of $z \gg 3$ galaxies is the Cosmic Evolution Survey (COSMOS; Scoville et al. 2007). This work takes advantage of the latest NIR-selected catalog of the COSMOS field, COSMOS2020 (Weaver et al. 2022). We adopt the “THE FARMER” catalog whose photometry was consistently measured across all bands by fitting galaxy profiles while taking into account PSF variations, depth, and source crowding. From the resulting photo-$z$ and stellar masses we construct a consistently measured SMF from $z = 0.2–7.5$, identify quiescent and star-forming systems, and study the build-up and assembly of stellar mass over 10 Gyr of cosmic history. This includes a detailed study of key moments in the development of galaxy populations from the Era of Reionization ($z > 6$) to the peak of star formation activity at Cosmic Noon ($z \sim 2$), up to the rich bimodality of star-forming and quiescent galaxies observed at more recent times ($z \le 1$).

This paper is organized as follows. Section 2 introduces the dataset chosen for this analysis from which samples are drawn and possible uncertainties discussed in Sect. 3. Section 4 provides a brief overview of the Schechter function. Results are presented in Sect. 5, including the presentation and fitting of the total, star-forming, and quiescent mass functions. These are compared to literature measurements in Sect. 6, whereupon fur-

ther discussion is had with regards to galaxy assembly, star formation cessation, and connections to dark matter halos. This work concludes in Sect. 7.

These results are computed adopting a standard $\Lambda$CDM cosmology with $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega_{m,0} = 0.3$, and $\Omega_{\Lambda,0} = 0.7$ throughout, such that the dimensionless Hubble parameter $h_{70} \equiv H_0/(70 \text{ km s}^{-1} \text{ Mpc}^{-1}) = 1$. Galaxy stellar masses, when derived from SED fitting, scale as the square of the luminosity distance (i.e., $D_L^2$); hence a factor of $h_{70}^{-1}$ is retained implicitly for all relevant measurements, unless explicitly noted otherwise (see Croton 2013 for an overview of $h$ and best practices). Estimates of stellar mass (henceforth $M$) assume an Chabrier (2003) initial mass function (IMF). All magnitudes are expressed in the AB system (Oke 1974), for which a flux $f_v$ in $\mu$Jy ($10^{-29} \text{ erg cm}^{-2} \text{ s}^{-1} \text{ Hz}^{-1}$) corresponds to $AB_v = 23.9–2.5 \log_{10}(f_v/\mu$Jy).

## 2. Data: COSMOS2020

We utilize the most recent release of the COSMOS catalog, COSMOS2020 (Weaver et al. 2022), comprised of ~100000 galaxies selected from a $izYJHK_s$ coadd with near-infrared depths approaching 26 AB ensuring a mass-selected sample complete down to $10^9 M_{\odot}$ at $z \approx 3$. This sample is complemented by extensive supporting photometry ranging from the UV to 8 $\mu$m over a total area of $2 \text{ deg}^2$, making it the widest near-infrared-selected multiwavelength catalog to these depths. This is made possible by the latest release of the UltraVISTA survey (DR4; McCracken et al. 2012; Moneti et al. 2023), which is the longest running near-infrared survey to date, and complemented by even deeper optical *grizy* data provided by Subaru's Hyper Suprime-Cam instrument (HSC PDR2; Aihara et al. 2019). Flanking this core imaging are $u$ band measurements made by the CLAUDS survey from the Canada-France-Hawaii Telescope (Sawicki et al. 2019) and deep *Spitzer*/IRAC imaging in all four channels, reprocessed to include nearly every exposure taken in COSMOS for use in this catalog (Euclid Collaboration 2022, and references therein). As in the previous COSMOS catalog by Laigle et al. (2016), several intermediate and narrow bands from both Subaru/Suprime-Cam and VISTA (Milvang-Jensen et al. 2013) are used to provide precise determinations of photometric redshifts.

In addition to the comprehensive set of data, photometry is performed by two methods to produce two separate catalogs. The catalog called CLASSIC contains aperture photometry extracted with Source Extractor (Bertin & Arnouts 1996) for the optical/near-infrared and with IRACLEAN (Hsieh et al. 2012) for the infrared (as in Laigle et al. 2016). The catalog called THE FARMER contains profile-fitting photometry extracted with The Farmer (Weaver et al. 2023), which uses The Tractor (Lang et al. 2016) to construct and apply parametric models to estimate source fluxes and does so consistently across more than 30 bands, including all four IRAC channels. Each photometric catalog is paired with photometric redshift and physical parameter estimates from both LePhare (Arnouts et al. 2002; Ilbert et al. 2006) and EAzY (Brammer et al. 2008) resulting in four measurements of photo-$z$ and $M$ per source. The combination of comprehensive deep images, independent photometry techniques, and independent SED fitting enables an unparalleled control of systematics and uncertainties, and hence further insight when inferring the underlying true nature of galaxy populations. See Weaver et al. (2022) for details.
# 3. Sample selection, completeness, and uncertainties

## 3.1. Selection function

Galaxies in COSMOS2020 are selected from a near-infrared izYJHKs CHI-MEAN coadded detection image (Szalay et al. 1999; Bertin 2010). The deepest band is *i*, with a $3\sigma$ sensitivity limit at ~27 mag, and the shallowest is $K_s$ at ~25 mag$^1$. The reddest selection band, $K_s$ delivers mass-selected samples to $z \le 4.5$. Furthermore, the comparably deeper $K_s$ imaging (by 0.5–0.8 mag relative to COSMOS2015) translates directly into a higher degree of completeness where low-$M$ samples are more mass complete down to lower masses. Taken together with IRAC, these images provide reliable detections of the reddest, oldest, and heavily dust obscured sources not seen in COSMOS2015.

Although the nominal area of the COSMOS survey is 2 deg$^2$, the most secure region comprises 1.279 deg$^2$ after removing contamination due to bright star halos and requiring an intersection of the deep near-infrared UltraVISTA imaging and the Subaru Suprime-Cam intermediate bands, where the photo-z performance is generally best. Compared to Davidzon et al. (2017) who use UltraVISTA DR2 (McCracken et al. 2012) imaging to probe the near-IR, here we rely on UltraVISTA DR4 (Moneti et al. 2023), reaching greater and simultaneously uniform NIR depth across the entire field (especially $K_s$). These non-uniform depths restricted the SMFs of Davidzon et al. (2017) to only the four deepest stripes in UltraVISTA spanning 0.62 deg$^2$. The gains leveraged in this work are hence five-fold. Firstly, the total number of recovered sources at all apparent magnitudes increases due to the additional area, improving statistical margins especially of rare sources. Secondly, the increased optical and near-infrared depths permit the recovery of fainter sources, improving flux completeness. Thirdly, the increased depth reduces photometric uncertainties leading to more precise, less biased photo-z and $M$. Fourth, the consistent depth and improved imaging performance from Hyper Suprime-Cam provides significantly better photometry relative to Suprime-Cam. Lastly, the wider field makes it less likely to become biased due to specific structures along the line-of-sight (e.g., clusters, voids) and instead probes a greater variety of environments affecting the evolution of galaxies within them.

As presented in Fig. 13 of Weaver et al. (2022), the best photo-z performance at $i > 22.5$ AB is achieved by the combination of The Farmer photometry and photo-z from LePhare with a $\delta z/(1+z)$ precision of <1% at $i \approx 20$ and <4% at $i \approx 26$ AB. For this reason, and to expedite comparisons with the most similar work in the literature (Ilbert et al. 2013; Davidzon et al. 2017), this work adopts the photometry from The Farmer paired with the photo-z, $M$, and rest-frame magnitudes estimated by LePhare. This combination will be henceforth referred to as COSMOS2020, unless explicitly stated otherwise. Our photo-z estimates correspond to the column LP_zPDF, defined as the median of the redshift likelihood distribution $L(z)$ (often referred to as $p(z)$). We choose to adopt $M$ estimates from MASS_MED as the median of the mass likelihood distributions are generally less susceptible to template fitting systematics than MASS_BEST which is taken at the minimum $\chi^2_{SED}$. We find that the two agree at all redshifts within 0.01 dex with a narrow scatter whose 68% range is below 0.05 dex.

Identification of stellar contaminants proceeds as described in Sect. 5.1 of Weaver et al. (2022). Briefly, LePhare computes the best-fit stellar template to each SED from a range of templates, including white dwarf and brown dwarf stars. Sources that achieve a $\chi^2$ from a stellar template fit lower than any galaxy template are removed. An additional criterion on morphology is used, identifying likely stars as point-like sources from the COSMOS HST/ACS mosaics (Koekemoer et al. 2007) and Subaru/HSC (PDR2, Aihara et al. 2019) for the $i < 23$ and $i < 21.5$ AB, respectively. Bright, resolved sources with SED shapes similar to stars are not removed.

An initial sample of 636 567 galaxies$^2$ in the range $0.2 < z \le 6.5$ is selected from the contiguous 1.27 deg$^2$ COMBINED region defined as the union of the UltraVISTA and Subaru/SC footprints but removing regions around bright stars$^3$. However, the source density at $6.5 < z \le 7.5$ becomes noticeably greater in the ultra-deep stripes and so we restrict sources in this particular $z$-bin to only the $0.716$ deg$^2$ of the ultra-deep stripes also covered by Subaru/SC but away from bright stars, finding 1327 galaxies. The construction of a robust stellar mass function requires precise redshifts and accurate stellar masses. Since an infrared detection is required to accurately measure galaxy stellar masses at $z \ge 2-3$, where this study seeks to cover new ground, 227 305 sources with $m_{ch1} > 26$ AB (i.e., $S/N \le 5$) are removed; $\approx 93\%$ of which fall below the mass completeness limit introduced in Sect. 3.3. An additional 13 644 sources with ambiguous redshifts are removed, quantified for simplicity by having >32% of their redshift probability lying outside $z_{phot} \pm 0.5$; $\approx 70\%$ of which fall below the mass completeness limit. Finally, we remove a further 1850 sources with unreliable SED fits with reduced $\chi^2 > 10$. These three necessary cuts create a final sample of 395 095 galaxies with precise redshifts and accurate stellar masses. We note that many excluded sources fail more than one of these three criteria (e.g., 90% are undetected in the infrared); Fig. A.1 shows their cross-section. Although the rejected objects are potentially interesting, their deeply uncertain nature does not permit us to reliably incorporate them into our measurements without a substantially more complex approach (e.g., to model potential outliers, see Hogg et al. 2010). Overall the excluded sources are predominantly faint (median $K_s \sim 25.7$ AB) such that their best-fit masses are generally below our adopted mass completeness limit and their potential for bias minimal within the mass ranges examined in this work.

## 3.2. Quiescent galaxy selection

Accurate identification of distant and/or unresolved quiescent galaxies has been a longstanding problem (see e.g., Leja et al. 2019b; Shahidi et al. 2020; Steinhardt et al. 2020). While quiescent galaxies can be identified by their low star formation rates, SFR estimates from fitting broad-band photometry are subject to large uncertainties and are dependent on the particular SED template library (Pacifici et al. 2015; Leja et al. 2019a; Carnall et al. 2019). The other most successful approach is that of rest-frame colors, as the lack of blue O- and B-type stars in quiescent populations implies a dearth of UV emission and a red optical continuum. Similar photometric UV to NIR SEDs are measured for heavily dust-enshrouded but otherwise star-forming systems, making the challenge of color-color selection that of distinguishing truly quiescent galaxies from dusty star-forming ones, in the absence of FIR measurements. Popular rest-frame color-color

<footnote n="1">¹ Computed at $3\sigma$ from 3'' apertures randomly placed in regions without detected sources; see Table 1 of Weaver et al. (2022).</footnote>

<footnote n="2">² Galaxies are selected by $lp\_type=0$, see catalog documentation.</footnote>

<footnote n="3">³ The COMBINED region corresponds to `FLAG_COMBINED=0`.</footnote>
# 3. Sample selection, completeness, and uncertainties

## 3.1. Selection function

Galaxies in COSMOS2020 are selected from a near-infrared izYJHKs CHI-MEAN coadded detection image (Szalay et al. 1999; Bertin 2010). The deepest band is *i*, with a $3\sigma$ sensitivity limit at ~27 mag, and the shallowest is $K_s$ at ~25 mag$^1$. The reddest selection band, $K_s$ delivers mass-selected samples to $z \le 4.5$. Furthermore, the comparably deeper $K_s$ imaging (by 0.5–0.8 mag relative to COSMOS2015) translates directly into a higher degree of completeness where low-$M$ samples are more mass complete down to lower masses. Taken together with IRAC, these images provide reliable detections of the reddest, oldest, and heavily dust obscured sources not seen in COSMOS2015.

Although the nominal area of the COSMOS survey is 2 deg$^2$, the most secure region comprises 1.279 deg$^2$ after removing contamination due to bright star halos and requiring an intersection of the deep near-infrared UltraVISTA imaging and the Subaru Suprime-Cam intermediate bands, where the photo-z performance is generally best. Compared to Davidzon et al. (2017) who use UltraVISTA DR2 (McCracken et al. 2012) imaging to probe the near-IR, here we rely on UltraVISTA DR4 (Moneti et al. 2023), reaching greater and simultaneously uniform NIR depth across the entire field (especially $K_s$). These non-uniform depths restricted the SMFs of Davidzon et al. (2017) to only the four deepest stripes in UltraVISTA spanning 0.62 deg$^2$. The gains leveraged in this work are hence five-fold. Firstly, the total number of recovered sources at all apparent magnitudes increases due to the additional area, improving statistical margins especially of rare sources. Secondly, the increased optical and near-infrared depths permit the recovery of fainter sources, improving flux completeness. Thirdly, the increased depth reduces photometric uncertainties leading to more precise, less biased photo-z and $M$. Fourth, the consistent depth and improved imaging performance from Hyper Suprime-Cam provides significantly better photometry relative to Suprime-Cam. Lastly, the wider field makes it less likely to become biased due to specific structures along the line-of-sight (e.g., clusters, voids) and instead probes a greater variety of environments affecting the evolution of galaxies within them.

As presented in Fig. 13 of Weaver et al. (2022), the best photo-z performance at $i > 22.5$ AB is achieved by the combination of The Farmer photometry and photo-z from LePhare with a $\delta z/(1+z)$ precision of <1% at $i \approx 20$ and <4% at $i \approx 26$ AB. For this reason, and to expedite comparisons with the most similar work in the literature (Ilbert et al. 2013; Davidzon et al. 2017), this work adopts the photometry from The Farmer paired with the photo-z, $M$, and rest-frame magnitudes estimated by LePhare. This combination will be henceforth referred to as COSMOS2020, unless explicitly stated otherwise. Our photo-z estimates correspond to the column LP_zPDF, defined as the median of the redshift likelihood distribution $L(z)$ (often referred to as $p(z)$). We choose to adopt $M$ estimates from MASS_MED as the median of the mass likelihood distributions are generally less susceptible to template fitting systematics than MASS_BEST which is taken at the minimum $\chi^2_{SED}$. We find that the two agree at all redshifts within 0.01 dex with a narrow scatter whose 68% range is below 0.05 dex.

Identification of stellar contaminants proceeds as described in Sect. 5.1 of Weaver et al. (2022). Briefly, LePhare computes the best-fit stellar template to each SED from a range of templates, including white dwarf and brown dwarf stars. Sources that achieve a $\chi^2$ from a stellar template fit lower than any galaxy template are removed. An additional criterion on morphology is used, identifying likely stars as point-like sources from the COSMOS HST/ACS mosaics (Koekemoer et al. 2007) and Subaru/HSC (PDR2, Aihara et al. 2019) for the $i < 23$ and $i < 21.5$ AB, respectively. Bright, resolved sources with SED shapes similar to stars are not removed.

An initial sample of 636 567 galaxies$^2$ in the range $0.2 < z \le 6.5$ is selected from the contiguous 1.27 deg$^2$ COMBINED region defined as the union of the UltraVISTA and Subaru/SC footprints but removing regions around bright stars$^3$. However, the source density at $6.5 < z \le 7.5$ becomes noticeably greater in the ultra-deep stripes and so we restrict sources in this particular $z$-bin to only the $0.716$ deg$^2$ of the ultra-deep stripes also covered by Subaru/SC but away from bright stars, finding 1327 galaxies. The construction of a robust stellar mass function requires precise redshifts and accurate stellar masses. Since an infrared detection is required to accurately measure galaxy stellar masses at $z \ge 2-3$, where this study seeks to cover new ground, 227 305 sources with $m_{ch1} > 26$ AB (i.e., $S/N \le 5$) are removed; $\approx 93\%$ of which fall below the mass completeness limit introduced in Sect. 3.3. An additional 13 644 sources with ambiguous redshifts are removed, quantified for simplicity by having >32% of their redshift probability lying outside $z_{phot} \pm 0.5$; $\approx 70\%$ of which fall below the mass completeness limit. Finally, we remove a further 1850 sources with unreliable SED fits with reduced $\chi^2 > 10$. These three necessary cuts create a final sample of 395 095 galaxies with precise redshifts and accurate stellar masses. We note that many excluded sources fail more than one of these three criteria (e.g., 90% are undetected in the infrared); Fig. A.1 shows their cross-section. Although the rejected objects are potentially interesting, their deeply uncertain nature does not permit us to reliably incorporate them into our measurements without a substantially more complex approach (e.g., to model potential outliers, see Hogg et al. 2010). Overall the excluded sources are predominantly faint (median $K_s \sim 25.7$ AB) such that their best-fit masses are generally below our adopted mass completeness limit and their potential for bias minimal within the mass ranges examined in this work.

## 3.2. Quiescent galaxy selection

Accurate identification of distant and/or unresolved quiescent galaxies has been a longstanding problem (see e.g., Leja et al. 2019b; Shahidi et al. 2020; Steinhardt et al. 2020). While quiescent galaxies can be identified by their low star formation rates, SFR estimates from fitting broad-band photometry are subject to large uncertainties and are dependent on the particular SED template library (Pacifici et al. 2015; Leja et al. 2019a; Carnall et al. 2019). The other most successful approach is that of rest-frame colors, as the lack of blue O- and B-type stars in quiescent populations implies a dearth of UV emission and a red optical continuum. Similar photometric UV to NIR SEDs are measured for heavily dust-enshrouded but otherwise star-forming systems, making the challenge of color-color selection that of distinguishing truly quiescent galaxies from dusty star-forming ones, in the absence of FIR measurements. Popular rest-frame color-color

<footnote n="1">¹ Computed at $3\sigma$ from 3'' apertures randomly placed in regions without detected sources; see Table 1 of Weaver et al. (2022).</footnote>

<footnote n="2">² Galaxies are selected by $lp\_type=0$, see catalog documentation.</footnote>

<footnote n="3">³ The COMBINED region corresponds to `FLAG_COMBINED=0`.</footnote>
selections include (U - V) vs. (V - J), (NUV - r) vs. (r - J), (NUV - r) vs. (r - K), and (FUV - r, r - J) (Williams et al. 2009; Ilbert et al. 2010; Arnouts et al. 2013; Leja et al. 2019b, respectively), in addition to observed-frame methods such as (B - z) vs. (z - K) for sources at 1.4 ≤ z ≤ 2.5 (Daddi et al. 2004). However, each implies a different definition of quiescence (see Leja et al. 2019b; Shahidi et al. 2020). These pairs of rest-frame colors may be measured directly from the observed photometry (where available) or derived from the best-fit templates that should largely agree with observed-band photometry, and can also be meaningfully extrapolated in cases where appropriate observed-frame data are not available. While the former may be adversely affected by uncertainties and systematics propagated from the observed photometry, the latter assumes that the SED template library or basis contains a suitable model. It is worth noting that efforts are being made to employ machine learning techniques to lessen the dependence and bias impact from model assumptions (e.g., Steinhardt et al. 2020), although no such method has gained comparable use as yet.

Recent work by Leja et al. (2019b) has clarified the efficacy of these color-color selections, proposing to extend the typically adopted U - V vs. V - J baseline further into the UV on the basis of stronger correlation with intrinsic sSFR as measured in simulated photometry (see also Arnouts et al. 2013). As such, catalogs containing reliable and deep FUV and NUV photometry have the advantage of being able to identify low-z quiescent populations which though extrapolation can be consistently identified up to the highest redshifts even when the observed-frame U band no longer corresponds to rest-frame U. At redder wavelengths the rest-frame J band is crucial for measuring the rest-frame stellar bulk, however, it becomes easily redshifted out of most deep surveys with IRAC photometry by z ≈ 2–3 which renders quiescent galaxy selections highly dependent on template libraries even at these intermediate redshifts.

In order to provide broad comparisons to other mass functions, in particular those measured in COSMOS (e.g., Ilbert et al. 2013; Davidzon et al. 2017), this work selects quiescent galaxies following the criteria introduced by Ilbert et al. (2013):

$$ (NUV - r) > 3(r - J) + 1 \text{ and } (NUV - r) > 3.1 \quad (1) $$

whereby the slant line runs perpendicular to the direction of increasing sSFR and parallel to an increase in dust attenuation to separate truly quiescent systems from otherwise dusty star-forming contaminants. For clarity, we henceforth refer to this (NUV - r) vs. (r - J) selection as NUVrJ. Generally, this selection has been found to approximate a cut in sSFR ≲ 10⁻¹¹ yr⁻¹ (Davidzon et al. 2018).

It is important to note that LePhare computes rest-frame magnitudes by attempting to find an observed frame band which directly probes the rest-frame band, which for the reasons listed above is less model-dependent than simply adopting the rest-frame colors of the best-fit galaxy template (or sSFR) and hence conveys a view of the true variety of observed galaxies that is less biased by our assumptions. Although the rest-frame photometry are estimated using a K-correction and color-term, they are most strongly dependent on the best-fit template at redshifts where observed photometry are not well matched to the rest-frame band and are then effectively extrapolated from the best-fit template (see Appendix 1 of Ilbert et al. 2005 for details).

Unbiased sample selection becomes increasingly difficult with redshift. By z ≈ 3, rest-frame J-band fluxes are no longer directly measured even by IRAC channel 2⁴ and so become

increasingly model-dependent and uncertain. While rest-frame NUV remains constrained even at z > 6, rest-frame r must be extrapolated by z ≈ 5 at which point differentiation between quiescent and dusty star-forming galaxies becomes statistically impossible. Because of this, and the expectation that quiescence is unlikely at such early times, selection of quiescent galaxies in this work is limited to z ≤ 5.5, noting that selection between 3 < z ≤ 5.5 is subject to a significantly higher degree of uncertainty. Advancement in this selection will only be made possible with deeper infrared imaging from facilities such as JWST (Rigby et al. 2023).

Selection of quiescent galaxies is presented in Fig. 1, showing star-forming and quiescent galaxies whose masses are above their respective mass limits (see Sect. 3.3). Error bars on the rest-frame colors are estimated from the quadrature addition of the observed-frame filter nearest to that of the rest-frame. Given that even rest-frame photometric uncertainties are generally inversely proportional to mass at given redshift, these error bars are most representative for median mass systems but are overestimated for bright, massive galaxies. In addition, uncertainties become increasingly underestimated at z > 3 as the extrapolation based on the best-fit template is not propagated. However, it is already obvious that the dividing power of the slant line is significantly diminished at z > 3 where the rest-frame J band is no longer directly observed. Both the quiescent and dusty star-forming regions appear to be devoid of sources by z ≈ 6, such that no intrinsically quiescent (misclassified or not) are found at these early times. Whether this is an intrinsic feature of galaxy populations or merely a selection effect is discussed below in Sect. 3.3. For these reasons this work restricts distinction between quiescent and star-forming galaxies to z ≤ 5.5. Galaxies at z > 5.5 in this work should be considered star-forming.

### 3.3. Completeness

Accurate estimates of completeness are crucial when inferring general properties about a population from an otherwise incomplete sample. While advancements in near-infrared facilities have enabled breakthroughs in selecting representative sample of galaxies by measuring their stellar bulk, samples are still mass limited and these mass limits evolve with z, owing to the faintness of increasingly distant galaxies. As such, mass limits are highly dependent on accurate estimates of survey depths and their impact on the selection function as it relates to the detection of the lowest mass galaxies.

There are three known populations which are expected to be missed by a near-infrared selection whose deepest images are on the blue-end of the selection function. Firstly, the faintest star-forming blue galaxies in the local Universe (e.g., dwarf irregular starbursts) may have detectable fluxes blueward of i, but they will not be included in existing izYJHKs selections as their near-infrared emission is too faint. However, their contribution is expected to be limited to only the low-mass end of the galaxy stellar mass function, and can henceforth be well characterized in large numbers by deeper small field surveys (e.g., CANDELS: Grogin et al. 2011; Koekemoer et al. 2011).

Secondly, quiescent galaxies are more difficult to detect than star-forming galaxies owing to a lack of rest-frame UV and blue optical emission, meaning their detection relies upon deep near-infrared imaging at both low and high-redshifts. For this reason, an izYJHKs detection will be insufficient to detect the quiescent systems compared to star-forming ones of the same mass, meaning that the quiescent sample will be as complete as the star-forming sample at higher masses. We compute their respective

<footnote n="4">⁴ Channel 3 and 4 are too shallow to provide useful constraints at z ≥ 3.</footnote>
selections include (U - V) vs. (V - J), (NUV - r) vs. (r - J), (NUV - r) vs. (r - K), and (FUV - r, r - J) (Williams et al. 2009; Ilbert et al. 2010; Arnouts et al. 2013; Leja et al. 2019b, respectively), in addition to observed-frame methods such as (B - z) vs. (z - K) for sources at 1.4 ≤ z ≤ 2.5 (Daddi et al. 2004). However, each implies a different definition of quiescence (see Leja et al. 2019b; Shahidi et al. 2020). These pairs of rest-frame colors may be measured directly from the observed photometry (where available) or derived from the best-fit templates that should largely agree with observed-band photometry, and can also be meaningfully extrapolated in cases where appropriate observed-frame data are not available. While the former may be adversely affected by uncertainties and systematics propagated from the observed photometry, the latter assumes that the SED template library or basis contains a suitable model. It is worth noting that efforts are being made to employ machine learning techniques to lessen the dependence and bias impact from model assumptions (e.g., Steinhardt et al. 2020), although no such method has gained comparable use as yet.

Recent work by Leja et al. (2019b) has clarified the efficacy of these color-color selections, proposing to extend the typically adopted U - V vs. V - J baseline further into the UV on the basis of stronger correlation with intrinsic sSFR as measured in simulated photometry (see also Arnouts et al. 2013). As such, catalogs containing reliable and deep FUV and NUV photometry have the advantage of being able to identify low-z quiescent populations which though extrapolation can be consistently identified up to the highest redshifts even when the observed-frame U band no longer corresponds to rest-frame U. At redder wavelengths the rest-frame J band is crucial for measuring the rest-frame stellar bulk, however, it becomes easily redshifted out of most deep surveys with IRAC photometry by z ≈ 2–3 which renders quiescent galaxy selections highly dependent on template libraries even at these intermediate redshifts.

In order to provide broad comparisons to other mass functions, in particular those measured in COSMOS (e.g., Ilbert et al. 2013; Davidzon et al. 2017), this work selects quiescent galaxies following the criteria introduced by Ilbert et al. (2013):

$$ (NUV - r) > 3(r - J) + 1 \text{ and } (NUV - r) > 3.1 \quad (1) $$

whereby the slant line runs perpendicular to the direction of increasing sSFR and parallel to an increase in dust attenuation to separate truly quiescent systems from otherwise dusty star-forming contaminants. For clarity, we henceforth refer to this (NUV - r) vs. (r - J) selection as NUVrJ. Generally, this selection has been found to approximate a cut in sSFR ≲ 10⁻¹¹ yr⁻¹ (Davidzon et al. 2018).

It is important to note that LePhare computes rest-frame magnitudes by attempting to find an observed frame band which directly probes the rest-frame band, which for the reasons listed above is less model-dependent than simply adopting the rest-frame colors of the best-fit galaxy template (or sSFR) and hence conveys a view of the true variety of observed galaxies that is less biased by our assumptions. Although the rest-frame photometry are estimated using a K-correction and color-term, they are most strongly dependent on the best-fit template at redshifts where observed photometry are not well matched to the rest-frame band and are then effectively extrapolated from the best-fit template (see Appendix 1 of Ilbert et al. 2005 for details).

Unbiased sample selection becomes increasingly difficult with redshift. By z ≈ 3, rest-frame J-band fluxes are no longer directly measured even by IRAC channel 2⁴ and so become

increasingly model-dependent and uncertain. While rest-frame NUV remains constrained even at z > 6, rest-frame r must be extrapolated by z ≈ 5 at which point differentiation between quiescent and dusty star-forming galaxies becomes statistically impossible. Because of this, and the expectation that quiescence is unlikely at such early times, selection of quiescent galaxies in this work is limited to z ≤ 5.5, noting that selection between 3 < z ≤ 5.5 is subject to a significantly higher degree of uncertainty. Advancement in this selection will only be made possible with deeper infrared imaging from facilities such as JWST (Rigby et al. 2023).

Selection of quiescent galaxies is presented in Fig. 1, showing star-forming and quiescent galaxies whose masses are above their respective mass limits (see Sect. 3.3). Error bars on the rest-frame colors are estimated from the quadrature addition of the observed-frame filter nearest to that of the rest-frame. Given that even rest-frame photometric uncertainties are generally inversely proportional to mass at given redshift, these error bars are most representative for median mass systems but are overestimated for bright, massive galaxies. In addition, uncertainties become increasingly underestimated at z > 3 as the extrapolation based on the best-fit template is not propagated. However, it is already obvious that the dividing power of the slant line is significantly diminished at z > 3 where the rest-frame J band is no longer directly observed. Both the quiescent and dusty star-forming regions appear to be devoid of sources by z ≈ 6, such that no intrinsically quiescent (misclassified or not) are found at these early times. Whether this is an intrinsic feature of galaxy populations or merely a selection effect is discussed below in Sect. 3.3. For these reasons this work restricts distinction between quiescent and star-forming galaxies to z ≤ 5.5. Galaxies at z > 5.5 in this work should be considered star-forming.

### 3.3. Completeness

Accurate estimates of completeness are crucial when inferring general properties about a population from an otherwise incomplete sample. While advancements in near-infrared facilities have enabled breakthroughs in selecting representative sample of galaxies by measuring their stellar bulk, samples are still mass limited and these mass limits evolve with z, owing to the faintness of increasingly distant galaxies. As such, mass limits are highly dependent on accurate estimates of survey depths and their impact on the selection function as it relates to the detection of the lowest mass galaxies.

There are three known populations which are expected to be missed by a near-infrared selection whose deepest images are on the blue-end of the selection function. Firstly, the faintest star-forming blue galaxies in the local Universe (e.g., dwarf irregular starbursts) may have detectable fluxes blueward of i, but they will not be included in existing izYJHKs selections as their near-infrared emission is too faint. However, their contribution is expected to be limited to only the low-mass end of the galaxy stellar mass function, and can henceforth be well characterized in large numbers by deeper small field surveys (e.g., CANDELS: Grogin et al. 2011; Koekemoer et al. 2011).

Secondly, quiescent galaxies are more difficult to detect than star-forming galaxies owing to a lack of rest-frame UV and blue optical emission, meaning their detection relies upon deep near-infrared imaging at both low and high-redshifts. For this reason, an izYJHKs detection will be insufficient to detect the quiescent systems compared to star-forming ones of the same mass, meaning that the quiescent sample will be as complete as the star-forming sample at higher masses. We compute their respective

<footnote n="4">⁴ Channel 3 and 4 are too shallow to provide useful constraints at z ≥ 3.</footnote>
<figure><img src="image_1.png" /><figcaption>Fig. 1. Galaxies are classified as either star-forming or quiescent in bins of redshift by selection in rest frame NUV - r and r - J colors, following Ilbert et al. (2013). The rest-frame J band is not directly measured at z ≥ 3.5 and so the sloped part of the selection box is dashed to indicate that the r - J colors are strongly model dependent. Typical uncertainties for the rest-frame colors of the star-forming and quiescent subsamples are estimated as medians of the nearest observed-frame bands, consistent with the derivation of the rest-frame colors. See the text for details.</figcaption></figure>

mass completeness limits separately and apply them consistently throughout.

Lastly, and similar to quiescent galaxies, the most heav-
ily dust obscured star-forming galaxies (A<sub>V</sub> ≫ 5) at high
redshift (z ≥ 2) will not present appreciable optical or near-
infrared fluxes to be detected in COSMOS2020, but unlike
quiescent galaxies are ubiquitously and efficiently detected
in far-infrared, submillimeter, and radio surveys (e.g.,
Schreiber et al. 2018b; Wang et al. 2019; Sun et al. 2021;
Jin et al. 2019; Fudamoto et al. 2020, 2021; Casey et al. 2021;
Shu et al. 2022). The nature and extent of this recently dis-
covered population remains difficult to quantify, owing to a
combination of complex selection functions and serendipitous
detections. The majority are likely to be missed by an izYJHK<sub>s</sub>
selection function even at the depths of COSMOS2020 (see
discussions in Sect. 6.3).

Although mass completeness estimates are presented in
Weaver et al. (2022), they are derived for a comparably less
secure sample than used in this work. For this reason the mass
completeness limits are rederived in identical fashion following
the method of Pozzetti et al. (2010). A critical advantage of this
method is that it does not rely on the theoretical mass distribu-
tion of galaxies fainter than the magnitude limit, but assumes
that those just above the threshold share similar properties
with the undetected ones, modulo a rescaling factor as detailed
below. This contrasts with studies that estimate mass complete-
ness either through injection-recovery of simulated sources, or
extrapolation of the observed distribution itself below the magni-
tude limit. The latter approach, in particular, may underestimate

the stellar mass limit if the galaxy sample is sparse due purely
to astrophysical reasons and not genuine incompleteness. In our
case, the sample in each z-bin is first cleaned by discounting the
1% worst-fit sources via χ². The stellar masses M of the 30%
faintest galaxies in channel 1 are then rescaled such that their
observed channel 1 apparent magnitude m<sub>ch1</sub> matches the IRAC
sensitivity limit m<sub>lim</sub> = 26:

$$
\mathrm{Log}_{10}(\mathcal{M}_{\mathrm{resc}}) = \mathrm{Log}_{10}(\mathcal{M}) + 0.4 (m_{\mathrm{ch1}} - 26.0). \quad (2)
$$

The limiting mass M<sub>lim</sub> is then taken to be the 95th percentile
of the M<sub>resc</sub> distribution. Finally 2nd order expansions in (1 + z)
are fitted to each M<sub>lim</sub> per z-bin up to z = 7 for the total sample,
and z = 5 for star-forming and quiescent samples, to produce a
smoothly evolving limiting mass. Limits above which samples
are ~70% mass complete (see justification below) are derived
consistently for the total sample from 0.2 < z ≤ 7.5 and for the
star-forming, and quiescent samples from 0.2 < z ≤ 5.5:

$$
\text{Total: } -3.23 \times 10^7 (1+z) + 7.83 \times 10^7 (1+z)^2 \quad (3)
$$

$$
\text{Star-forming: } -5.77 \times 10^7 (1+z) + 8.66 \times 10^7 (1+z)^2 \quad (4)
$$

$$
\text{Quiescent: } -3.79 \times 10^8 (1+z) + 2.98 \times 10^8 (1+z)^2 \quad (5)
$$

and are shown in Fig. 2. Despite the more conservative selec-
tion adopted in this work, the derived mass completeness limits
are essentially identical to those derived in Weaver et al. (2022),
which indicates the robustness of these limits against sample
selections. For reference, at z ≈ 1 the quiescent sample is com-
plete at ~0.25 dex in mass above the star-forming sample, grow-
ing to ~0.45 dex above by z ≈ 3–5.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Galaxies are classified as either star-forming or quiescent in bins of redshift by selection in rest frame NUV - r and r - J colors, following Ilbert et al. (2013). The rest-frame J band is not directly measured at z ≥ 3.5 and so the sloped part of the selection box is dashed to indicate that the r - J colors are strongly model dependent. Typical uncertainties for the rest-frame colors of the star-forming and quiescent subsamples are estimated as medians of the nearest observed-frame bands, consistent with the derivation of the rest-frame colors. See the text for details.</figcaption></figure>

mass completeness limits separately and apply them consistently throughout.

Lastly, and similar to quiescent galaxies, the most heav-
ily dust obscured star-forming galaxies (A<sub>V</sub> ≫ 5) at high
redshift (z ≥ 2) will not present appreciable optical or near-
infrared fluxes to be detected in COSMOS2020, but unlike
quiescent galaxies are ubiquitously and efficiently detected
in far-infrared, submillimeter, and radio surveys (e.g.,
Schreiber et al. 2018b; Wang et al. 2019; Sun et al. 2021;
Jin et al. 2019; Fudamoto et al. 2020, 2021; Casey et al. 2021;
Shu et al. 2022). The nature and extent of this recently dis-
covered population remains difficult to quantify, owing to a
combination of complex selection functions and serendipitous
detections. The majority are likely to be missed by an izYJHK<sub>s</sub>
selection function even at the depths of COSMOS2020 (see
discussions in Sect. 6.3).

Although mass completeness estimates are presented in
Weaver et al. (2022), they are derived for a comparably less
secure sample than used in this work. For this reason the mass
completeness limits are rederived in identical fashion following
the method of Pozzetti et al. (2010). A critical advantage of this
method is that it does not rely on the theoretical mass distribu-
tion of galaxies fainter than the magnitude limit, but assumes
that those just above the threshold share similar properties
with the undetected ones, modulo a rescaling factor as detailed
below. This contrasts with studies that estimate mass complete-
ness either through injection-recovery of simulated sources, or
extrapolation of the observed distribution itself below the magni-
tude limit. The latter approach, in particular, may underestimate

the stellar mass limit if the galaxy sample is sparse due purely
to astrophysical reasons and not genuine incompleteness. In our
case, the sample in each z-bin is first cleaned by discounting the
1% worst-fit sources via χ². The stellar masses M of the 30%
faintest galaxies in channel 1 are then rescaled such that their
observed channel 1 apparent magnitude m<sub>ch1</sub> matches the IRAC
sensitivity limit m<sub>lim</sub> = 26:

$$
\mathrm{Log}_{10}(\mathcal{M}_{\mathrm{resc}}) = \mathrm{Log}_{10}(\mathcal{M}) + 0.4 (m_{\mathrm{ch1}} - 26.0). \quad (2)
$$

The limiting mass M<sub>lim</sub> is then taken to be the 95th percentile
of the M<sub>resc</sub> distribution. Finally 2nd order expansions in (1 + z)
are fitted to each M<sub>lim</sub> per z-bin up to z = 7 for the total sample,
and z = 5 for star-forming and quiescent samples, to produce a
smoothly evolving limiting mass. Limits above which samples
are ~70% mass complete (see justification below) are derived
consistently for the total sample from 0.2 < z ≤ 7.5 and for the
star-forming, and quiescent samples from 0.2 < z ≤ 5.5:

$$
\text{Total: } -3.23 \times 10^7 (1+z) + 7.83 \times 10^7 (1+z)^2 \quad (3)
$$

$$
\text{Star-forming: } -5.77 \times 10^7 (1+z) + 8.66 \times 10^7 (1+z)^2 \quad (4)
$$

$$
\text{Quiescent: } -3.79 \times 10^8 (1+z) + 2.98 \times 10^8 (1+z)^2 \quad (5)
$$

and are shown in Fig. 2. Despite the more conservative selec-
tion adopted in this work, the derived mass completeness limits
are essentially identical to those derived in Weaver et al. (2022),
which indicates the robustness of these limits against sample
selections. For reference, at z ≈ 1 the quiescent sample is com-
plete at ~0.25 dex in mass above the star-forming sample, grow-
ing to ~0.45 dex above by z ≈ 3–5.
<figure><img src="image_2.png" /><figcaption>**Fig. 2.** Stellar mass completeness limits. The three panels show, from left to right, the density of the total, star-forming, and quiescent samples with their respective mass completeness limits. The limits are derived in discrete z-bins following Pozzetti et al. (2010) with magnitude limits adopted from IRAC channel 1 (filled circles), then interpolated with Eqs. (3)–(5), respectively, resulting in the solid lines. For quiescent galaxies, consistent estimates using $K_s$ are also shown (empty circles). To verify these estimates we compute conservative mass limits from simulated spectra following a delta-like burst at $z = 15$ which are normalized separately to the $K_s$ and channel 1 magnitude limits; they are shown by the dotted and solid red lines, respectively. See the text for details.</figcaption></figure>

There remains an additional incompleteness arising from the fact that mass completeness is derived from IRAC channel 1 photometry, and yet our *izYJHKs* selection function does not include sources which are only identified in IRAC. As discussed in Davidzon et al. (2017), despite this drawback it is also disad-vantageous to estimate the mass completeness from any one of the six *izYJHKs* detection bands either. For instance, the red-desst, $K_s$, samples the rest-frame stellar bulk only out to $z \lesssim 2$, and will tend to underestimate stellar masses at $z \ge 2$ until $z \approx 4.5$ when it becomes a UV tracer. Thankfully, it is possi-ble to estimate the impact of this additional incompleteness by examining a sample of galaxies common to this work and those of the comparably deeper ~200 arcmin² CANDELS-COSMOS catalog of Nayyeri et al. (2017). This analysis is performed and discussed at length in Weaver et al. (2022). In summary, 75% of CANDELS sources at $M_{\text{lim}}$ are recovered by both THE FARMER and CLASSIC COSMOS2020 catalogs, which combined with the choice of $M_{\text{lim}}$ being the 95th percentile of $M_{\text{resc}}$ implies that $M_{\text{lim}}$ of the total sample corresponds to a mass completeness of ~70%.

As shown in Fig. 2, the mass limit of the total sample is almost identical with that of the star-forming sample at all red-shifts. This is unsurprising as star-forming galaxies generally dominate galaxy demographics. As in our *NUVrJ* selection from Fig. 1, we see the lack of quiescent systems at $z > 5.5$ and so we consider the star-forming subsample to be statistically equivalent to the total sample at $z > 5.5$.

Deriving a consistent mass completeness limit for quiescent galaxies is more challenging. It is well known that the predom-inantly older, redder stellar populations of quiescent systems imply higher mass-to-light ratios than found in star-forming galaxies which in turn imply a lower degree of mass complete-ness at the same flux limit. Figure 2 shows that our channel 1 derived mass completeness falls $\ge 0.2$ dex lower than the bulk of the quiescent sources. Taken at face value, this would appear to indicate a lack of quiescent systems at low-intermediate masses at $z > 2$. Again, we have cause for concern as channel 1 is not included in our *izyJHKs* detection strategy, and hence our selection function is liable to miss IRAC-bright sources that are not present in bluer bands. Furthermore, due to the predomi-nantly red optical spectral slopes of quiescent galaxies, contin-

uum emission in $K_s$ should be lower than in channel 1, implying that one would need a deeper $K_s$ magnitude limit to detect the same sources from shallower channel 1 imaging. In other words, it is expected that there are red sources visible only in channel 1 that are not included in our *izYJHKs* detected catalog. Further-more, by $z \approx 4.5$ $K_s$ no longer traces mass but rather UV con-tinuum (channel 1 does so by $z \approx 8$). Worse, the UV continuum in these systems is expected to be faint (even given ‘frostings’ of star formation in *NUVrJ*-selected post-starburst systems). There-fore, while deriving a mass completeness from channel 1 mag-nitude limits may be suboptimal due to possible selection effects with regards to a *izYJHKs*-selected sample, we cannot turn to $K_s$ as it is no longer a mass indicator at $z \ge 4.5$.

One solution is to modify our selection function by incor-porating *izYJHKs* undetected IRAC-only sources into our sam-ple. While a systematic search for IRAC only sources is ongo-ing, their identification is made difficult not only because of the significantly lower resolution and consequently higher source crowding in necessarily deep IRAC images, but these sources also lack optical/NIR data which is, by definition, insufficiently deep to identify low-z interlopers. Deep mid-infrared (MIR) data redward of IRAC over sufficient degree-scale areas do not cur-rently exist, making a determination of redshift, mass, and qui-escent nature of these IRAC-detected sources problematically uncertain.

For the time being, whether or not this absence of intermediate-low mass quiescent systems is astrophysical can-not be determined. Literature measurements do not yield more mass-complete quiescent samples as the UltraVISTA DR4 NIR depths are now similar to those from even the deepest small-field NIR imaging, $H_{160} \approx 25.9$ (e.g., CANDELS, Tomczak et al. 2014), and even so, comparisons are hampered by field-to-field and photometric systematics. Comparisons with other stellar mass functions measured more consistently in COSMOS (e.g., Davidzon et al. 2017; Ilbert et al. 2013) are still affected by sys-tematics from comparably less certain measurements from pre-vious, shallower data, despite the lessened impact of cosmic variance. Additionally, comparisons with Davidzon et al. (2017) are complicated by the fact that they refit the photometry of Laigle et al. (2016) to produce new $z$ and $M$ estimates. Modi-fications include expanding the redshift baseline from $z \le 6$ to
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Stellar mass completeness limits. The three panels show, from left to right, the density of the total, star-forming, and quiescent samples with their respective mass completeness limits. The limits are derived in discrete z-bins following Pozzetti et al. (2010) with magnitude limits adopted from IRAC channel 1 (filled circles), then interpolated with Eqs. (3)–(5), respectively, resulting in the solid lines. For quiescent galaxies, consistent estimates using $K_s$ are also shown (empty circles). To verify these estimates we compute conservative mass limits from simulated spectra following a delta-like burst at $z = 15$ which are normalized separately to the $K_s$ and channel 1 magnitude limits; they are shown by the dotted and solid red lines, respectively. See the text for details.</figcaption></figure>

There remains an additional incompleteness arising from the fact that mass completeness is derived from IRAC channel 1 photometry, and yet our *izYJHKs* selection function does not include sources which are only identified in IRAC. As discussed in Davidzon et al. (2017), despite this drawback it is also disad-vantageous to estimate the mass completeness from any one of the six *izYJHKs* detection bands either. For instance, the red-desst, $K_s$, samples the rest-frame stellar bulk only out to $z \lesssim 2$, and will tend to underestimate stellar masses at $z \ge 2$ until $z \approx 4.5$ when it becomes a UV tracer. Thankfully, it is possi-ble to estimate the impact of this additional incompleteness by examining a sample of galaxies common to this work and those of the comparably deeper ~200 arcmin² CANDELS-COSMOS catalog of Nayyeri et al. (2017). This analysis is performed and discussed at length in Weaver et al. (2022). In summary, 75% of CANDELS sources at $M_{\text{lim}}$ are recovered by both THE FARMER and CLASSIC COSMOS2020 catalogs, which combined with the choice of $M_{\text{lim}}$ being the 95th percentile of $M_{\text{resc}}$ implies that $M_{\text{lim}}$ of the total sample corresponds to a mass completeness of ~70%.

As shown in Fig. 2, the mass limit of the total sample is almost identical with that of the star-forming sample at all red-shifts. This is unsurprising as star-forming galaxies generally dominate galaxy demographics. As in our *NUVrJ* selection from Fig. 1, we see the lack of quiescent systems at $z > 5.5$ and so we consider the star-forming subsample to be statistically equivalent to the total sample at $z > 5.5$.

Deriving a consistent mass completeness limit for quiescent galaxies is more challenging. It is well known that the predom-inantly older, redder stellar populations of quiescent systems imply higher mass-to-light ratios than found in star-forming galaxies which in turn imply a lower degree of mass complete-ness at the same flux limit. Figure 2 shows that our channel 1 derived mass completeness falls $\ge 0.2$ dex lower than the bulk of the quiescent sources. Taken at face value, this would appear to indicate a lack of quiescent systems at low-intermediate masses at $z > 2$. Again, we have cause for concern as channel 1 is not included in our *izyJHKs* detection strategy, and hence our selection function is liable to miss IRAC-bright sources that are not present in bluer bands. Furthermore, due to the predomi-nantly red optical spectral slopes of quiescent galaxies, contin-

uum emission in $K_s$ should be lower than in channel 1, implying that one would need a deeper $K_s$ magnitude limit to detect the same sources from shallower channel 1 imaging. In other words, it is expected that there are red sources visible only in channel 1 that are not included in our *izYJHKs* detected catalog. Further-more, by $z \approx 4.5$ $K_s$ no longer traces mass but rather UV con-tinuum (channel 1 does so by $z \approx 8$). Worse, the UV continuum in these systems is expected to be faint (even given ‘frostings’ of star formation in *NUVrJ*-selected post-starburst systems). There-fore, while deriving a mass completeness from channel 1 mag-nitude limits may be suboptimal due to possible selection effects with regards to a *izYJHKs*-selected sample, we cannot turn to $K_s$ as it is no longer a mass indicator at $z \ge 4.5$.

One solution is to modify our selection function by incor-porating *izYJHKs* undetected IRAC-only sources into our sam-ple. While a systematic search for IRAC only sources is ongo-ing, their identification is made difficult not only because of the significantly lower resolution and consequently higher source crowding in necessarily deep IRAC images, but these sources also lack optical/NIR data which is, by definition, insufficiently deep to identify low-z interlopers. Deep mid-infrared (MIR) data redward of IRAC over sufficient degree-scale areas do not cur-rently exist, making a determination of redshift, mass, and qui-escent nature of these IRAC-detected sources problematically uncertain.

For the time being, whether or not this absence of intermediate-low mass quiescent systems is astrophysical can-not be determined. Literature measurements do not yield more mass-complete quiescent samples as the UltraVISTA DR4 NIR depths are now similar to those from even the deepest small-field NIR imaging, $H_{160} \approx 25.9$ (e.g., CANDELS, Tomczak et al. 2014), and even so, comparisons are hampered by field-to-field and photometric systematics. Comparisons with other stellar mass functions measured more consistently in COSMOS (e.g., Davidzon et al. 2017; Ilbert et al. 2013) are still affected by sys-tematics from comparably less certain measurements from pre-vious, shallower data, despite the lessened impact of cosmic variance. Additionally, comparisons with Davidzon et al. (2017) are complicated by the fact that they refit the photometry of Laigle et al. (2016) to produce new $z$ and $M$ estimates. Modi-fications include expanding the redshift baseline from $z \le 6$ to
<figure><img src="image_3.png" /><figcaption>Fig. 3. Adopted estimates for the total uncertainty $\sigma_{\Phi}$ (solid) as a function of stellar mass at several redshifts (by color), for mass-complete samples only. Contributions include uncertainties from Poisson noise $\sigma_N$ (dashed), Cosmic Variance $\sigma_{CV}$ (dotted), and SED fitting $\sigma_{SED}$ (dash-dot). 10% and 100% levels are indicated for reference.</figcaption></figure>

$z \le 8$ and adding in additional SED templates to probe star-bursting galaxies with rising star formation histories as well as dust-obscured systems. As demonstrated by a recent comparison by Lustig et al. (2023), these modifications produce lower number densities of massive galaxies compared to the original measurements.

We turn to theoretical frameworks to investigate this further. Namely, we use Flexible Stellar Population Synthesis (FSPS; Conroy et al. 2009; Conroy & Gunn 2010) to estimate the stellar masses produced by extraordinarily old stellar population produced by a delta-burst evaluated at $z = 15$. We assume a Chabrier (2003) IMF with $\log_{10}(Z/Z_{\odot}) = -0.3$. At each epoch (i.e., $z$) we normalize the evolved model spectrum to match our channel 1 observed frame magnitude limits (26.0 mag) and estimate the corresponding stellar mass shown by the solid red curve in Fig. 2. This exercise is repeated consistently with the $K_s$ magnitude limit (25.5 mag) shown by the dashed red curve.

We caution however, that the assumption that these models accurately describe real high-$z$ quiescent galaxies is becoming increasingly dubious. Such a system formed in a monolithic delta-burst just following the big bang cannot reach quiescence (as defined by NUVrJ) above $z \approx 5.3$, and yet remarkably mature systems at $z \approx 4-5$ have already been reported in the literature (e.g., Schreiber et al. 2018a; Tanaka et al. 2019; Valentino et al. 2020). Worse, the spread of mass-to-light ratios found in quiescent systems means that a single mass completeness limit for all quiescent galaxies at a given redshift is ill-defined even for consistently selected (e.g., UVJ, NUVrJ, sSFR) samples, resulting in a non negligible selection effects.

Nonetheless, we stress that the difference in completeness between the effectively flux-complete $K_s$-derived mass limit and the mass limit derived from channel 1 magnitudes is only ~0.3 dex, which is typically less than a single bin in our analysis. In light of this considerable uncertainty, we adopt the optimistic quiescent galaxy mass completeness limits derived via Pozzetti et al. (2010) from channel 1 magnitudes with the caveat that the lowest mass bins in each measurement at $z \ge 2$ are potentially incomplete.

We indicate the $K_s$ mass limit for quiescent samples throughout. Also, we note that our quiescent mass limit, although consistently determined in Davidzon et al. (2017), is more conservative despite the deeper NIR data. As will be discussed in Sect. 5.2, we attribute this to an overestimate of the quiescent galaxy mass completeness by Davidzon et al. (2017).

The difference in mass completeness between the star-forming and quiescent samples presents an additional complication. Because the star-forming galaxies can be reliably detected to lower masses than quiescent galaxies from our izYJHKs selection function, the low-mass regime of the total SMF at a given epoch, as defined in this work, does not actually include contributions from the lowest-mass quiescent systems. Therefore the total SMF is simply the star-forming SMF at masses below the quiescent mass limit. The effect from the undetected low-mass quiescent systems is almost certainly negligible: a simple extrapolation from the evolution in the number density of $z \approx 0.2-2.5$ low-mass quiescent galaxies predict a number density at least $10 \times$ lower than that of the star-forming galaxies at $z > 3$. Furthermore, given that general uncertainties are still above the 10% level at low-masses even at $z \approx 2$ (Fig. 3), any bias arising in the total number density from undetected low-mass quiescent systems is lower than other sources of error.

### 3.4. Derivation of the $1/V_{max}$ correction

Intrinsically faint galaxies at any given redshift are more likely to be missed by survey selection functions compared to brighter sources. For a NIR-selected sample, this equates to a mass bias by which low-luminosity, low-mass galaxies can only be detected in smaller volumes (i.e. out to lower $z$) relative to brighter, more massive ones which could be detected if they were at higher $z$. This is the well known Malmquist bias (Malmquist 1920, 1922).

The most straight-forward approach to correct for such a bias is the $1/V_{max}$ method of Schmidt (1968), which has enjoyed significant popularity owing to its simplicity. Briefly, the $1/V_{max}$ method statistically corrects for selection incompleteness by weighting each detected object by the maximum comoving volume in which it can be observed, given the characteristics of the telescope survey. The $V_{max}$ estimate per individual object is computed after finding the maximum redshift $z_{max}$ by which the best-fit SED would no longer be observable<fnref n="5" /> because of the survey’s flux limit. On the other hand, the minimum redshift ($z_{min}$) should be the one at which the source would become too bright and saturate the camera, although in practice is the lower boundary of the $z$ bin in which the considered galaxy lies ($z_{low}$). Therefore, for the $i$th galaxy inside the bin $z_{low} < z < z_{high}$, the maximum observable volume is

$$V_{max,i} = \frac{4\pi}{3} \frac{\Omega}{\Omega^{\text{sky}}} (D_{\text{cov}}(\min(z_{\text{max},i}, z_{\text{high}}))^{3} - D_{\text{cov}}(\max(z_{\text{min},i}, z_{\text{low}}))^{3}), \quad (6)$$

where $\Omega$ is the solid angle subtended by the sample, $\Omega^{\text{sky}} \equiv 41\,253 \, \text{deg}^2$ is the solid angle of a sphere, and $D_{\text{cov}}(z)$ is the comoving distance at $z$ (Hogg 1999). If $z_{max}$ exceeds the upper boundary of the redshift bin, the latter is used instead, meaning that the brightest sources are often assigned a weight that corresponds to the full volume of that redshift slice. As such, this correction is expected to be significant for only the faintest, lowest-mass sources in a given $z$-bin. While it is nonparametric and does not assume a functional form of the SMF, the $1/V_{max}$ technique does assume that samples are drawn from a uniform

<footnote n="5">⁵ $z_{max}$ estimates are computed with ALF (Ilbert et al. 2005).</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Adopted estimates for the total uncertainty $\sigma_{\Phi}$ (solid) as a function of stellar mass at several redshifts (by color), for mass-complete samples only. Contributions include uncertainties from Poisson noise $\sigma_{N}$ (dashed), Cosmic Variance $\sigma_{CV}$ (dotted), and SED fitting $\sigma_{SED}$ (dash-dot). 10% and 100% levels are indicated for reference.</figcaption></figure>

$z \le 8$ and adding in additional SED templates to probe star-bursting galaxies with rising star formation histories as well as dust-obscured systems. As demonstrated by a recent comparison by Lustig et al. (2023), these modifications produce lower number densities of massive galaxies compared to the original measurements.

We turn to theoretical frameworks to investigate this further. Namely, we use Flexible Stellar Population Synthesis (FSPS; Conroy et al. 2009; Conroy & Gunn 2010) to estimate the stellar masses produced by extraordinarily old stellar population produced by a delta-burst evaluated at $z = 15$. We assume a Chabrier (2003) IMF with $\log_{10}(Z/Z_{\odot}) = -0.3$. At each epoch (i.e., $z$) we normalize the evolved model spectrum to match our channel 1 observed frame magnitude limits (26.0 mag) and estimate the corresponding stellar mass shown by the solid red curve in Fig. 2. This exercise is repeated consistently with the $K_s$ magnitude limit (25.5 mag) shown by the dashed red curve.

We caution however, that the assumption that these models accurately describe real high-$z$ quiescent galaxies is becoming increasingly dubious. Such a system formed in a monolithic delta-burst just following the big bang cannot reach quiescence (as defined by NUVrJ) above $z \approx 5.3$, and yet remarkably mature systems at $z \approx 4-5$ have already been reported in the literature (e.g., Schreiber et al. 2018a; Tanaka et al. 2019; Valentino et al. 2020). Worse, the spread of mass-to-light ratios found in quiescent systems means that a single mass completeness limit for all quiescent galaxies at a given redshift is ill-defined even for consistently selected (e.g., UVJ, NUVrJ, sSFR) samples, resulting in a non negligible selection effects.

Nonetheless, we stress that the difference in completeness between the effectively flux-complete $K_s$-derived mass limit and the mass limit derived from channel 1 magnitudes is only ~0.3 dex, which is typically less than a single bin in our analysis. In light of this considerable uncertainty, we adopt the optimistic quiescent galaxy mass completeness limits derived via Pozzetti et al. (2010) from channel 1 magnitudes with the caveat that the lowest mass bins in each measurement at $z \ge 2$ are potentially incomplete.

We indicate the $K_s$ mass limit for quiescent samples throughout. Also, we note that our quiescent mass limit, although consistently determined in Davidzon et al. (2017), is more conservative despite the deeper NIR data. As will be discussed in Sect. 5.2, we attribute this to an overestimate of the quiescent galaxy mass completeness by Davidzon et al. (2017).

The difference in mass completeness between the star-forming and quiescent samples presents an additional complication. Because the star-forming galaxies can be reliably detected to lower masses than quiescent galaxies from our izYJHKs selection function, the low-mass regime of the total SMF at a given epoch, as defined in this work, does not actually include contributions from the lowest-mass quiescent systems. Therefore the total SMF is simply the star-forming SMF at masses below the quiescent mass limit. The effect from the undetected low-mass quiescent systems is almost certainly negligible: a simple extrapolation from the evolution in the number density of $z \approx 0.2-2.5$ low-mass quiescent galaxies predict a number density at least $10 \times$ lower than that of the star-forming galaxies at $z > 3$. Furthermore, given that general uncertainties are still above the 10% level at low-masses even at $z \approx 2$ (Fig. 3), any bias arising in the total number density from undetected low-mass quiescent systems is lower than other sources of error.

### 3.4. Derivation of the $1/V_{max}$ correction

Intrinsically faint galaxies at any given redshift are more likely to be missed by survey selection functions compared to brighter sources. For a NIR-selected sample, this equates to a mass bias by which low-luminosity, low-mass galaxies can only be detected in smaller volumes (i.e. out to lower $z$) relative to brighter, more massive ones which could be detected if they were at higher $z$. This is the well known Malmquist bias (Malmquist 1920, 1922).

The most straight-forward approach to correct for such a bias is the $1/V_{max}$ method of Schmidt (1968), which has enjoyed significant popularity owing to its simplicity. Briefly, the $1/V_{max}$ method statistically corrects for selection incompleteness by weighting each detected object by the maximum comoving volume in which it can be observed, given the characteristics of the telescope survey. The $V_{max}$ estimate per individual object is computed after finding the maximum redshift $z_{max}$ by which the best-fit SED would no longer be observable<fnref n="5" /> because of the survey’s flux limit. On the other hand, the minimum redshift ($z_{min}$) should be the one at which the source would become too bright and saturate the camera, although in practice is the lower boundary of the $z$ bin in which the considered galaxy lies ($z_{low}$). Therefore, for the $i$th galaxy inside the bin $z_{low} < z < z_{high}$, the maximum observable volume is

$$V_{max,i} = \frac{4\pi}{3} \frac{\Omega}{\Omega^{\text{sky}}} (D_{\text{cov}}(\min(z_{\text{max},i}, z_{\text{high}}))^{3} - D_{\text{cov}}(\max(z_{\text{min},i}, z_{\text{low}}))^{3}), \quad (6)$$

where $\Omega$ is the solid angle subtended by the sample, $\Omega^{\text{sky}} \equiv 41\,253 \, \text{deg}^2$ is the solid angle of a sphere, and $D_{\text{cov}}(z)$ is the comoving distance at $z$ (Hogg 1999). If $z_{max}$ exceeds the upper boundary of the redshift bin, the latter is used instead, meaning that the brightest sources are often assigned a weight that corresponds to the full volume of that redshift slice. As such, this correction is expected to be significant for only the faintest, lowest-mass sources in a given $z$-bin. While it is nonparametric and does not assume a functional form of the SMF, the $1/V_{max}$ technique does assume that samples are drawn from a uniform

<footnote n="5">⁵ $z_{max}$ estimates are computed with ALF (Ilbert et al. 2005).</footnote>
spatial distribution, which is not accurate in the case of over- or under-dense environments (Efstathiou & Bond 1999). However, the assumption of a uniform spatial distribution is expected to be problematic only at $z < 1$, where large-scale structures have fully assembled, or in narrower surveys that can be biased by structures at smaller scales. Other methods exist which do not make this assumption such as STY (Sandage et al. 1979) and SWML (Efstathiou et al. 1988), a parametric and nonparametric maximum likelihood method, respectively. Already, Davidzon et al. (2017) found that the constraints provided by COSMOS2015 were sufficiently strong for the $1/V_{max}$ method as well as more complex methods (e.g., STY, SWML) to essentially converge. With even stronger constraints and larger effective area provided by COSMOS2020, we can expect even better agreement, with minimal advantages to the more complex methods. More extensive discussions on strengths and weaknesses of the various approaches can be found in the literature (e.g., Ilbert et al. 2004; Binggeli & Jerjen 1998; Johnston 2011; Takeuchi 2000; Weigel et al. 2016).

### 3.5. Further considerations of uncertainty

We adopt a statistical error budget on the SMF number density $\Phi$ consisting of the quadrature addition of Poisson noise ($\sigma_N$), cosmic variance fluctuations ($\sigma_{cv}$), and uncertainties on masses induced by SED fitting ($\sigma_{SED}$) such that $\sigma_{\Phi} = (\sigma_N^2 + \sigma_{cv}^2 + \sigma_{SED}^2)^{1/2}$. Figure 3 shows the composition of the total error budget from $z = 1.1$ to 6.5 as a function of stellar mass for mass-complete bins.

#### 3.5.1. Poisson noise

Poisson noise arises from processes wherein the abundance of a discrete quantity (or counts) is measured. Although in the limit of many events a Poisson process becomes indistinguishable with that of a Gaussian, in small number counts it can be the dominant source of uncertainty. Here we compute the Poisson error $\sigma_N$ for each mass bin as $\sqrt{N}$ where $N$ is the number of objects in that bin. These values are recomputed for the star-forming and quiescent subsamples separately. As shown in Fig. 3, the fractional contribution from $\sigma_N$ increases with mass and redshift with the largest contribution at $M > 10^{11.5} M_\odot$.

The discrete nature of Poisson processes allows us to also provide upper limits on bins containing zero detected galaxies. Following Table 1 of Gehrels (1986), the statistical upper limit on $\Phi(N=0)$ for a given observed volume $V$ is $\sigma_{N,limit} = 1.841/V$. See Ebeling (2003) and Weigel et al. (2016) for details and further discussions.

#### 3.5.2. Cosmic variance

It is well established that galaxy properties are correlated with environmental density (i.e., clustering). Galaxy clusters, while being an important laboratory for galaxy evolution, are not typical of galaxy environments. Because of their density, they impart a higher overall normalization to the stellar mass function. More noticeably, they tend to inflate the massive end of the mass function as they preferentially contain the most massive systems. This environmental field-to-field bias (so-called “cosmic variance”) is a topic of intense study, and is a key component to accurately assessing sample uncertainties when trying to infer universal or intrinsic properties of galaxies.

There are many published methods to estimate cosmic variance, based on numerical simulations (Bhowmick et al. 2020; Ucci et al. 2021), analytical models calibrated to observations solved either using linear theory (Moster et al. 2011; Trapp & Furlanetto 2020) or on forward simulation corrections to linear theory (Steinhardt et al. 2021), and observationally (e.g., Driver & Robotham 2010). Trenti & Stiavelli (2008) combines results from cosmological simulations with analytical predictions.

Cosmic variance $\sigma_{cv}$ is estimated following Steinhardt et al. (2021), who adapt the methods of Moster et al. (2011) which, importantly, scale with stellar mass (up to $10^{11.25} M_\odot$) and are commonly adopted for use in $0.1 \le z \le 3.5$ measurements as that is the redshift range in which the Moster et al. (2011) calculator was devised. However, above $z \approx 3.5$ these estimates become increasingly underestimated, so Steinhardt et al. (2021) use linear perturbation theory to extend this work more reasonably to the early Universe, while maintaining agreement at $z < 3$. Although environmental density has known covariance with star formation (e.g., Davidzon et al. 2016; Bolzonella et al. 2010), we assume cosmic variance is equivalent between star-forming and quiescent subsamples. As shown in Fig. 3, the fractional contribution from $\sigma_{cv}$ dominates the error budget for low-$M$ systems at all redshifts and becomes progressively more important at high-$M$ with increasing $z$.

#### 3.5.3. SED fitting uncertainties and bias

Another consideration is the uncertainty on the stellar mass estimate provided by the SED fit. Figure 4 shows the likelihood distributions on stellar mass at fixed redshifts and masses produced by LePhare sampling at $\Delta \log_{10}(M/M_\odot) = 0.025$ dex. Trends with width of the likelihood distributions indicate that mass is best constrained for low-redshift, massive (i.e., bright) sources. Although there is nonzero skewness and kurtosis in individual cases, the overall median distribution is symmetric. This is expected, as the uncertainty on stellar mass is essentially a measurement of the range of allowable templates and their normalization in the fitting procedure, and thus $\sigma_{SED}$ scales with photometric uncertainties. However, these likelihood distributions on stellar mass should be treated as lower limits as they do not take into account any covariance with redshift.

While these typical per-bin distributions can be valuable, especially for injecting noise into measurements from simulations, attempting to compute its contribution to the SMF, $\sigma_{SED}$, using the typical width in a given bin is suboptimal as the wings of neighboring mass-bins contribute asymmetrically. To address this, we use the individual mass likelihood distributions to draw 1000 independent realizations of the galaxy stellar mass function and thereby directly estimate the variance produced by the mass uncertainties, which we take as the 68% range about the median number density per bin of mass. As shown in Fig. 3, the contribution from $\sigma_{SED}$ become dominant only at $M > 10^{11.5} M_\odot$, in some cases becoming larger than unity. They are comparable to contributions from $\sigma_N$ across the entire mass range.

It is important to note that this does not account for systematic biases arising from SED fitting, such as assumptions of the stellar initial mass function$^6$, potential photo-z offsets, and assumptions as to the star formation histories that all propagate into the stellar mass determination. However, concerning the latter case, we show in Weaver et al. (2022) that the combination of

<footnote>$^6$ For example, $M$ computed assuming Salpeter (1955) is on average 0.24 dex larger than those computed assuming Chabrier (2003).</footnote>
spatial distribution, which is not accurate in the case of over- or under-dense environments (Efstathiou & Bond 1999). However, the assumption of a uniform spatial distribution is expected to be problematic only at $z < 1$, where large-scale structures have fully assembled, or in narrower surveys that can be biased by structures at smaller scales. Other methods exist which do not make this assumption such as STY (Sandage et al. 1979) and SWML (Efstathiou et al. 1988), a parametric and nonparametric maximum likelihood method, respectively. Already, Davidzon et al. (2017) found that the constraints provided by COSMOS2015 were sufficiently strong for the $1/V_{max}$ method as well as more complex methods (e.g., STY, SWML) to essentially converge. With even stronger constraints and larger effective area provided by COSMOS2020, we can expect even better agreement, with minimal advantages to the more complex methods. More extensive discussions on strengths and weaknesses of the various approaches can be found in the literature (e.g., Ilbert et al. 2004; Binggeli & Jerjen 1998; Johnston 2011; Takeuchi 2000; Weigel et al. 2016).

### 3.5. Further considerations of uncertainty

We adopt a statistical error budget on the SMF number density $\Phi$ consisting of the quadrature addition of Poisson noise ($\sigma_N$), cosmic variance fluctuations ($\sigma_{cv}$), and uncertainties on masses induced by SED fitting ($\sigma_{SED}$) such that $\sigma_{\Phi} = (\sigma_N^2 + \sigma_{cv}^2 + \sigma_{SED}^2)^{1/2}$. Figure 3 shows the composition of the total error budget from $z = 1.1$ to 6.5 as a function of stellar mass for mass-complete bins.

#### 3.5.1. Poisson noise

Poisson noise arises from processes wherein the abundance of a discrete quantity (or counts) is measured. Although in the limit of many events a Poisson process becomes indistinguishable with that of a Gaussian, in small number counts it can be the dominant source of uncertainty. Here we compute the Poisson error $\sigma_N$ for each mass bin as $\sqrt{N}$ where N is the number of objects in that bin. These values are recomputed for the star-forming and quiescent subsamples separately. As shown in Fig. 3, the fractional contribution from $\sigma_N$ increases with mass and redshift with the largest contribution at $M > 10^{11.5} M_\odot$.

The discrete nature of Poisson processes allows us to also provide upper limits on bins containing zero detected galaxies. Following Table 1 of Gehrels (1986), the statistical upper limit on $\Phi(N=0)$ for a given observed volume $V$ is $\sigma_{N,limit} = 1.841/V$. See Ebeling (2003) and Weigel et al. (2016) for details and further discussions.

#### 3.5.2. Cosmic variance

It is well established that galaxy properties are correlated with environmental density (i.e., clustering). Galaxy clusters, while being an important laboratory for galaxy evolution, are not typical of galaxy environments. Because of their density, they impart a higher overall normalization to the stellar mass function. More noticeably, they tend to inflate the massive end of the mass function as they preferentially contain the most massive systems. This environmental field-to-field bias (so-called “cosmic variance”) is a topic of intense study, and is a key component to accurately assessing sample uncertainties when trying to infer universal or intrinsic properties of galaxies.

There are many published methods to estimate cosmic variance, based on numerical simulations (Bhowmick et al. 2020; Ucci et al. 2021), analytical models calibrated to observations solved either using linear theory (Moster et al. 2011; Trapp & Furlanetto 2020) or on forward simulation corrections to linear theory (Steinhardt et al. 2021), and observationally (e.g., Driver & Robotham 2010). Trenti & Stiavelli (2008) combines results from cosmological simulations with analytical predictions.

Cosmic variance $\sigma_{cv}$ is estimated following Steinhardt et al. (2021), who adapt the methods of Moster et al. (2011) which, importantly, scale with stellar mass (up to $10^{11.25} M_\odot$) and are commonly adopted for use in $0.1 \le z \le 3.5$ measurements as that is the redshift range in which the Moster et al. (2011) calculator was devised. However, above $z \approx 3.5$ these estimates become increasingly underestimated, so Steinhardt et al. (2021) use linear perturbation theory to extend this work more reasonably to the early Universe, while maintaining agreement at $z < 3$. Although environmental density has known covariance with star formation (e.g., Davidzon et al. 2016; Bolzonella et al. 2010), we assume cosmic variance is equivalent between star-forming and quiescent subsamples. As shown in Fig. 3, the fractional contribution from $\sigma_{cv}$ dominates the error budget for low-$M$ systems at all redshifts and becomes progressively more important at high-$M$ with increasing $z$.

#### 3.5.3. SED fitting uncertainties and bias

Another consideration is the uncertainty on the stellar mass estimate provided by the SED fit. Figure 4 shows the likelihood distributions on stellar mass at fixed redshifts and masses produced by LePhare sampling at $\Delta \log_{10}(M/M_\odot) = 0.025$ dex. Trends with width of the likelihood distributions indicate that mass is best constrained for low-redshift, massive (i.e., bright) sources. Although there is nonzero skewness and kurtosis in individual cases, the overall median distribution is symmetric. This is expected, as the uncertainty on stellar mass is essentially a measurement of the range of allowable templates and their normalization in the fitting procedure, and thus $\sigma_{SED}$ scales with photometric uncertainties. However, these likelihood distributions on stellar mass should be treated as lower limits as they do not take into account any covariance with redshift.

While these typical per-bin distributions can be valuable, especially for injecting noise into measurements from simulations, attempting to compute its contribution to the SMF, $\sigma_{SED}$, using the typical width in a given bin is suboptimal as the wings of neighboring mass-bins contribute asymmetrically. To address this, we use the individual mass likelihood distributions to draw 1000 independent realizations of the galaxy stellar mass function and thereby directly estimate the variance produced by the mass uncertainties, which we take as the 68% range about the median number density per bin of mass. As shown in Fig. 3, the contribution from $\sigma_{SED}$ become dominant only at $M > 10^{11.5} M_\odot$, in some cases becoming larger than unity. They are comparable to contributions from $\sigma_N$ across the entire mass range.

It is important to note that this does not account for systematic biases arising from SED fitting, such as assumptions of the stellar initial mass function$^6$, potential photo-z offsets, and assumptions as to the star formation histories that all propagate into the stellar mass determination. However, concerning the latter case, we show in Weaver et al. (2022) that the combination of

<footnote>$^6$ For example, $M$ computed assuming Salpeter (1955) is on average 0.24 dex larger than those computed assuming Chabrier (2003).</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 4. Likelihood distributions of galaxy stellar mass are derived from SED fitting with LePhare and assume a fixed redshift. Individual distributions (gray) are summarized by a median stack (blue) grouped by redshift and stellar mass (indicated in ranges of $log_{10}(M/M_{\odot})$). Estimates of standard deviation $\sigma$ are shown. The size of a typical mass bin used in this work is 0.25 dex, indicated by the pair of dotted gray lines in each panel.</figcaption></figure>

The Farmer and LePhare achieves a subpercent photo-z bias even for faint, high-z sources (-0.004 $\delta z/1 + z$ at $25 < i < 27$) improving over other works including COSMOS2015. Systematic errors cannot be combined with random errors, and so additionally complicate measurements of the SMF. Given this indication of relatively low bias arising from the photo-z and significantly better constrained SEDs relative to previous measurements, we omit these considerations in the present work. See Marchesini et al. (2009) and Davidzon et al. (2017) for detailed discussions of various sources of bias and their effect on the SMF.

### 3.5.4. Eddington bias

The number of low-mass galaxies is orders of magnitude larger than the highest-mass systems, and so a randomly chosen galaxy is overwhelmingly likely to be lower-mass. If even a small fraction of such truly low-mass systems scatter toward high-mass (owing to a $M$ overestimate) it can significantly change the poisson-dominated high-mass number density estimate. The converse situation, while depleting the high-mass end, would have virtually no effect on the low-mass estimates. This is the well known Eddington bias (Eddington 1913). While generally understood to mean that there is a net bias leading to the overestimation of the density of massive galaxies, a small, but highly asymmetric uncertainty on the mass of low-mass systems can similarly generate such a bias that affects the shape of the low-mass regime. See Grazian et al. (2015) for further discussions.

Effectively correcting for Eddington biases has been a leading point of discussion in recent literature, generally favoring the convolution of the fitting function with a kernel that describes the uncertainty in stellar mass (e.g., Davidzon et al. 2017; Ilbert et al. 2013). Recently, Adams et al. (2021) compared the effect of using three different forms for the convolution kernel finding a significant difference in the inferred intrinsic SMF. Alternative approaches have been proposed (e.g.,

Leja et al. 2019a), developed a nonparametric formalism for incorporating $\sigma_{\Phi}$ into an unbinned Likelihood fitting, whereby multiple realizations of the parent catalog are made, each time sampling stellar mass from the mass likelihood distributions of each galaxy. In this work we primarily adopt the traditional convolution kernel method to estimate Eddington bias. At the same time, we also fit the mass function using the method of Leja et al. (2019a), and so follow their approach where relevant.

## 4. The Schechter function

Galaxy luminosity and stellar mass functions can be described empirically by the parametric formulation first introduced by Schechter (1976) in the context of the local Universe. Since then, the Schechter function has been adopted ubiquitously in statistical studies of galaxy mass assembly. It is more convenient to express the number density of galaxies per logarithmic mass bin $d \log M$ as given by Weigel et al. (2016):

$$ 
\begin{aligned}
n_{\text{gal}} &= \Phi d \log M \\
&= \ln(10) \Phi^* e^{-10^{\log M - \log M^*}} \times (10^{\log M - \log M^*})^{\alpha+1} d \log M,
\end{aligned} 
\quad (7) 
$$

which describes a power law of slope $\alpha$ at masses smaller than a characteristic stellar mass $M^*$, whereupon the function is cut off by a high-mass exponential tail. The overall normalization is set by $\Phi^*$, which corresponds to the number density at $M^*$. Although $\Phi(M)$ is expected to evolve smoothly with redshift such that $\Phi(M, z)$ can be mapped from secondary functions $M^*(z)$, $\Phi^*(z)$, $\alpha(z)$, most literature applications fit for each parameter independently at each redshift. A notable exception is Leja et al. (2019a).

Many studies have since found evidence that the total galaxy population at low-$z$ is better described as the coaddition of two Schechter functions (e.g., Peng et al. 2010) whereby the low-mass and high-mass end acquire individual normalizations ($\Phi_1^*$ and $\Phi_2^*$) and low-mass slopes ($\alpha_1$ and $\alpha_2$) but retaining a single characteristic stellar mass $M^*$. This so-called ‘Double’ Schechter form is as follows:

$$ 
\begin{aligned}
\Phi d \log M = \ln(10) e^{-10^{\log M - \log M^*}} \\
\times [\Phi_1^* (10^{\log M - \log M^*})^{\alpha_1 + 1} \\
+ \Phi_2^* (10^{\log M - \log M^*})^{\alpha_2 + 1}]
d \log M.
\end{aligned} 
\quad (8) 
$$

The Double Schechter function has been found to be as suitable description in most of the studies in the local Universe and up to $z \approx 2$. At higher redshifts, it is less clear whether this is a better description of $\Phi$ than a single Schechter function. Moreover, deviations have been recently been reported at $z > 7$ (e.g., Stefanon et al. 2021) in which the observed stellar mass function is better described by a power law than by any Schechter profile.

## 5. Results

Now having established the selections and methods adopted in this work, we present the resulting measurement of the SMF for the total, star-forming, and quiescent samples. We investigate also the evolution of number densities and quiescent fractions at fixed mass, and then fit the SMF for each sample with several methods to derive the evolution of the best-fit Schechter parameters.
<figure><img src="image_5.png" /><figcaption>Fig. 4. Likelihood distributions of galaxy stellar mass are derived from SED fitting with LePhare and assume a fixed redshift. Individual distributions (gray) are summarized by a median stack (blue) grouped by redshift and stellar mass (indicated in ranges of $log_{10}(M/M_{\odot})$). Estimates of standard deviation $\sigma$ are shown. The size of a typical mass bin used in this work is 0.25 dex, indicated by the pair of dotted gray lines in each panel.</figcaption></figure>

The Farmer and LePhare achieves a subpercent photo-z bias even for faint, high-z sources ($-0.004\delta z/1 + z$ at $25 < i < 27$) improving over other works including COSMOS2015. Systematic errors cannot be combined with random errors, and so additionally complicate measurements of the SMF. Given this indication of relatively low bias arising from the photo-z and significantly better constrained SEDs relative to previous measurements, we omit these considerations in the present work. See Marchesini et al. (2009) and Davidzon et al. (2017) for detailed discussions of various sources of bias and their effect on the SMF.

### 3.5.4. Eddington bias

The number of low-mass galaxies is orders of magnitude larger than the highest-mass systems, and so a randomly chosen galaxy is overwhelmingly likely to be lower-mass. If even a small fraction of such truly low-mass systems scatter toward high-mass (owing to a $M$ overestimate) it can significantly change the poisson-dominated high-mass number density estimate. The converse situation, while depleting the high-mass end, would have virtually no effect on the low-mass estimates. This is the well known Eddington bias (Eddington 1913). While generally understood to mean that there is a net bias leading to the overestimation of the density of massive galaxies, a small, but highly asymmetric uncertainty on the mass of low-mass systems can similarly generate such a bias that affects the shape of the low-mass regime. See Grazian et al. (2015) for further discussions.

Effectively correcting for Eddington biases has been a leading point of discussion in recent literature, generally favoring the convolution of the fitting function with a kernel that describes the uncertainty in stellar mass (e.g., Davidzon et al. 2017; Ilbert et al. 2013). Recently, Adams et al. (2021) compared the effect of using three different forms for the convolution kernel finding a significant difference in the inferred intrinsic SMF. Alternative approaches have been proposed (e.g.,

Leja et al. 2019a), developed a nonparametric formalism for incorporating $\sigma_{\Phi}$ into an unbinned Likelihood fitting, whereby multiple realizations of the parent catalog are made, each time sampling stellar mass from the mass likelihood distributions of each galaxy. In this work we primarily adopt the traditional convolution kernel method to estimate Eddington bias. At the same time, we also fit the mass function using the method of Leja et al. (2019a), and so follow their approach where relevant.

## 4. The Schechter function

Galaxy luminosity and stellar mass functions can be described empirically by the parametric formulation first introduced by Schechter (1976) in the context of the local Universe. Since then, the Schechter function has been adopted ubiquitously in statistical studies of galaxy mass assembly. It is more convenient to express the number density of galaxies per logarithmic mass bin $d \log M$ as given by Weigel et al. (2016):

$$ 
\begin{aligned}
n_{\text{gal}} &= \Phi d \log M \\
&= \ln(10) \Phi^* e^{-10^{\log M - \log M^*}} \times (10^{\log M - \log M^*})^{\alpha+1} d \log M,
\end{aligned} 
\quad (7) 
$$

which describes a power law of slope $\alpha$ at masses smaller than a characteristic stellar mass $M^*$, whereupon the function is cut off by a high-mass exponential tail. The overall normalization is set by $\Phi^*$, which corresponds to the number density at $M^*$. Although $\Phi(M)$ is expected to evolve smoothly with redshift such that $\Phi(M, z)$ can be mapped from secondary functions $M^*(z)$, $\Phi^*(z)$, $\alpha(z)$, most literature applications fit for each parameter independently at each redshift. A notable exception is Leja et al. (2019a).

Many studies have since found evidence that the total galaxy population at low-$z$ is better described as the coaddition of two Schechter functions (e.g., Peng et al. 2010) whereby the low-mass and high-mass end acquire individual normalizations ($\Phi_1^*$ and $\Phi_2^*$) and low-mass slopes ($\alpha_1$ and $\alpha_2$) but retaining a single characteristic stellar mass $M^*$. This so-called ‘Double’ Schechter form is as follows:

$$ 
\begin{aligned}
\Phi d \log M = \ln(10) e^{-10^{\log M - \log M^*}} \\
& \times [\Phi_1^* (10^{\log M - \log M^*})^{\alpha_1 + 1} \\
& + \Phi_2^* (10^{\log M - \log M^*})^{\alpha_2 + 1}] d \log M.
\end{aligned} 
\quad (8) 
$$

The Double Schechter function has been found to be as suitable description in most of the studies in the local Universe and up to $z \approx 2$. At higher redshifts, it is less clear whether this is a better description of $\Phi$ than a single Schechter function. Moreover, deviations have been recently been reported at $z > 7$ (e.g., Stefanon et al. 2021) in which the observed stellar mass function is better described by a power law than by any Schechter profile.

## 5. Results

Now having established the selections and methods adopted in this work, we present the resulting measurement of the SMF for the total, star-forming, and quiescent samples. We investigate also the evolution of number densities and quiescent fractions at fixed mass, and then fit the SMF for each sample with several methods to derive the evolution of the best-fit Schechter parameters.
<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Evolution of the galaxy stellar mass function over 12 redshift bins ($0.2 < z \le 7.5$) for the total, star-forming, and quiescent samples. Mass incomplete bins based on the channel 1 limiting magnitude are not shown. For the quiescent sample, bins that are fully incomplete based on the $K_s$ limiting magnitude are indicated by cross hatching.</figcaption></figure>

## 5.1. The total galaxy stellar mass function

We measure the SMF for our total (star-forming and quiescent) sample divided in 12 redshift bins from $z = 0.2$ to 7.5 in fixed bins of stellar mass $\Delta M = 0.25$ above the mass limit for that $z$-range[^7]. Shown in the left panel of Fig. 5, the shape and normalization of the SMF changes considerably over the ~10 billion years corresponding to this redshift range. At $z \le 3$, the SMF features a smooth, monotonically decreasing number density at low-$M$ that flattens before falling off steeply at $M \approx 10^{11} M_{\odot}$, around the so-called ‘knee’ of the function. Its overall shape and normalization remains roughly constant until $z \sim 1.5$ indicating a lack of mass growth at recent times, consistent with the decline of the cosmic star formation rate (Madau & Dickinson 2014). However, at $z \sim 1.5$ the normalization decreases dramatically. The number density at the knee is only ~1% of its $z \approx 0.5$ value with the fastest growth occuring on the low mass end, consistent with mass “downsizing in time” (Cowie et al. 1996; Neistein et al. 2006; Fontanot et al. 2009). The slope of the low mass end steeps with time and is expected to be driven by physical processes including outflows, star formation efficiency, slope of the main sequence, and mergers (see discussions in Peng et al. 2014). At the same time, the knee itself becomes difficult to distinguish as the SMF takes the form of a smooth power law, and we become increasingly unable to constrain the low-mass end of the SMF due to worsening mass completeness (mass incomplete points are omitted in the figure). As expected from Fig. 3, the overall uncertainty grows significantly with increasing redshift and mass. We note that at a few epochs (e.g., $z \approx 5$), the SMF is not strictly monotonic, likely driven by systematics rather than a real physical phenomenon (see discussions in Sect. 6).

The evolution of the total SMF measured in this work is compared with literature results in Fig. 6. We begin at $z \approx 0.2$ by comparing with two SMFs previously measured in the same field: Ilbert et al. (2013) and Davidzon et al. (2017); both use photo-$z$ and $M$ computed with LePhare over COSMOS, out to $z \approx 4.0$ and $z \approx 5.5$, respectively, with the nearly same redshift binning scheme as we use. We note one exception that (Ilbert et al. 2013) bins sources at $3.0 < z < 4.0$ whereas (Davidzon et al. 2017) uses $3.0 < z \le 3.5$ and $3.5 < z \le 4.5$, and so we opt to follow the scheme of Davidzon et al. (2017) as it allows a comparison up to higher redshift and omit the compari-

son with the highest-$z$ measurement of Ilbert et al. (2013). Additionally, Davidzon et al. (2017) only considered sources in the ultra-deep regions of COSMOS2015 (corresponding to UltraVISTA DR2) as the spatial inhomogeniety of the NIR bands implies a significant variation in selection function and mass completeness between the deep and ultra-deep regions. Thankfully, COSMOS2020 (corresponding to UltraVISTA DR4) contains nearly uniform NIR coverage ($\Delta z \approx 0.4$ mag, Moneti et al. 2023) and so can leverage an area almost 2× larger out to at least $z \approx 6.5$, beyond which the source density becomes clearly different between the deep and ultra-deep stripes. Thus, for the $6.5 < z \le 7.5$ bin we use the $0.72 \text{ deg}^2$ subset of our primary area covered by the UltraVISTA deep stripes.

At $z < 2.5$, in the range that all three can be directly compared, we find excellent agreement with Ilbert et al. (2013) and Davidzon et al. (2017). This is unsurprising, as measurements from Davidzon et al. (2017) at $z < 2.5$ are adopted directly from COSMOS2015 and computed nearly identically as Ilbert et al. (2013) but with deeper imaging. Our measurements are similar, but with visibly less structure around the knee especially between $0.8 < z \le 1.1$ where the volume density of sources is slightly higher, with the greatest increase at $M \sim 10^{10} M_{\odot}$. However, at $z \approx 2.5$ Davidzon et al. (2017) predicts a lower volume density than either Ilbert et al. (2013) or our measurements, which lie in agreement. This is not surprising, as Davidzon et al. (2017) expanded the template library of LePhare to include both starbursting and dust-obscured systems (see Sect. 3 of Davidzon et al. 2017) and so our work is naturally more similar to that of Ilbert et al. (2013) – a trend made clearer when looking at quiescent systems in the Sect. 5.2. At $z > 3.0$, we observe a significantly higher volume densities of massive galaxies compared to Davidzon et al. (2017), although the general shape of the low-mass end of the SMF remains similar. However, if we similarly limit our SMF to only the Ultra-Deep region we find it decreases to about the existing lower $1\sigma$ limit at $M > 10^{11} M_{\odot}$. We speculate that this may be due to the presence of (proto-)clusters preferentially in the Deep region at $3 < z \le 4$ (see Brinch et al. 2023), one of which has been recently spectroscopically confirmed by McConachie et al. (2022). Constraints from Grazian et al. (2015) can be introduced at $3.5 < z \le 4.5$, although at a significantly higher degree of caution as Grazian et al. (2015) uses a combination of smaller CANDELS fields (GOODS-South, UDS, HUDF) with an overall area of ~12× smaller than the present study which implies a higher degree of uncertainty from cosmic variances and Poisson noise,

<footnote>[^7]: The range of the lowest-$M$ bin at each $z$ is adjusted to include its respective mass limit at its lowest extent. Bins of width <0.05 dex are discarded.</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Evolution of the galaxy stellar mass function over 12 redshift bins ($0.2 < z \le 7.5$) for the total, star-forming, and quiescent samples. Mass incomplete bins based on the channel 1 limiting magnitude are not shown. For the quiescent sample, bins that are fully incomplete based on the $K_s$ limiting magnitude are indicated by cross hatching.</figcaption></figure>

## 5.1. The total galaxy stellar mass function

We measure the SMF for our total (star-forming and quiescent) sample divided in 12 redshift bins from $z = 0.2$ to 7.5 in fixed bins of stellar mass $\Delta M = 0.25$ above the mass limit for that $z$-range[^7]. Shown in the left panel of Fig. 5, the shape and normalization of the SMF changes considerably over the ~10 billion years corresponding to this redshift range. At $z \le 3$, the SMF features a smooth, monotonically decreasing number density at low-$M$ that flattens before falling off steeply at $M \approx 10^{11} M_{\odot}$, around the so-called ‘knee’ of the function. Its overall shape and normalization remains roughly constant until $z \sim 1.5$ indicating a lack of mass growth at recent times, consistent with the decline of the cosmic star formation rate (Madau & Dickinson 2014). However, at $z \sim 1.5$ the normalization decreases dramatically. The number density at the knee is only ~1% of its $z \approx 0.5$ value with the fastest growth occuring on the low mass end, consistent with mass “downsizing in time” (Cowie et al. 1996; Neistein et al. 2006; Fontanot et al. 2009). The slope of the low mass end steeps with time and is expected to be driven by physical processes including outflows, star formation efficiency, slope of the main sequence, and mergers (see discussions in Peng et al. 2014). At the same time, the knee itself becomes difficult to distinguish as the SMF takes the form of a smooth power law, and we become increasingly unable to constrain the low-mass end of the SMF due to worsening mass completeness (mass incomplete points are omitted in the figure). As expected from Fig. 3, the overall uncertainty grows significantly with increasing redshift and mass. We note that at a few epochs (e.g., $z \approx 5$), the SMF is not strictly monotonic, likely driven by systematics rather than a real physical phenomenon (see discussions in Sect. 6).

The evolution of the total SMF measured in this work is compared with literature results in Fig. 6. We begin at $z \approx 0.2$ by comparing with two SMFs previously measured in the same field: Ilbert et al. (2013) and Davidzon et al. (2017); both use photo-$z$ and $M$ computed with LePhare over COSMOS, out to $z \approx 4.0$ and $z \approx 5.5$, respectively, with the nearly same redshift binning scheme as we use. We note one exception that (Ilbert et al. 2013) bins sources at $3.0 < z < 4.0$ whereas (Davidzon et al. 2017) uses $3.0 < z \le 3.5$ and $3.5 < z \le 4.5$, and so we opt to follow the scheme of Davidzon et al. (2017) as it allows a comparison up to higher redshift and omit the comparison

with the highest-$z$ measurement of Ilbert et al. (2013). Additionally, Davidzon et al. (2017) only considered sources in the ultra-deep regions of COSMOS2015 (corresponding to UltraVISTA DR2) as the spatial inhomogeniety of the NIR bands implies a significant variation in selection function and mass completeness between the deep and ultra-deep regions. Thankfully, COSMOS2020 (corresponding to UltraVISTA DR4) contains nearly uniform NIR coverage ($\Delta z \approx 0.4$ mag, Moneti et al. 2023) and so can leverage an area almost 2× larger out to at least $z \approx 6.5$, beyond which the source density becomes clearly different between the deep and ultra-deep stripes. Thus, for the $6.5 < z \le 7.5$ bin we use the $0.72 \text{ deg}^2$ subset of our primary area covered by the UltraVISTA deep stripes.

At $z < 2.5$, in the range that all three can be directly compared, we find excellent agreement with Ilbert et al. (2013) and Davidzon et al. (2017). This is unsurprising, as measurements from Davidzon et al. (2017) at $z < 2.5$ are adopted directly from COSMOS2015 and computed nearly identically as Ilbert et al. (2013) but with deeper imaging. Our measurements are similar, but with visibly less structure around the knee especially between $0.8 < z \le 1.1$ where the volume density of sources is slightly higher, with the greatest increase at $M \sim 10^{10} M_{\odot}$. However, at $z \approx 2.5$ Davidzon et al. (2017) predicts a lower volume density than either Ilbert et al. (2013) or our measurements, which lie in agreement. This is not surprising, as Davidzon et al. (2017) expanded the template library of LePhare to include both starbursting and dust-obscured systems (see Sect. 3 of Davidzon et al. 2017) and so our work is naturally more similar to that of Ilbert et al. (2013) – a trend made clearer when looking at quiescent systems in the Sect. 5.2. At $z > 3.0$, we observe a significantly higher volume densities of massive galaxies compared to Davidzon et al. (2017), although the general shape of the low-mass end of the SMF remains similar. However, if we similarly limit our SMF to only the Ultra-Deep region we find it decreases to about the existing lower $1\sigma$ limit at $M > 10^{11} M_{\odot}$. We speculate that this may be due to the presence of (proto-)clusters preferentially in the Deep region at $3 < z \le 4$ (see Brinch et al. 2023), one of which has been recently spectroscopically confirmed by McConachie et al. (2022). Constraints from Grazian et al. (2015) can be introduced at $3.5 < z \le 4.5$, although at a significantly higher degree of caution as Grazian et al. (2015) uses a combination of smaller CANDELS fields (GOODS-South, UDS, HUDF) with an overall area of ~12× smaller than the present study which implies a higher degree of uncertainty from cosmic variances and Poisson noise,

<footnote>[^7]: The range of the lowest-$M$ bin at each $z$ is adjusted to include its respective mass limit at its lowest extent. Bins of width <0.05 dex are discarded.</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 6.** Evolution of the galaxy stellar mass function across 12 redshift bins ($0.2 < z \le 7.5$). The $0.2 < z \le 0.5$ SMF from the first redshift bin is repeated in each panel for reference shown by the purple dotted curve. Two other COSMOS stellar mass functions from Ilbert et al. (2013) and Davidzon et al. (2017) are shown for comparison, along with Grazian et al. (2015) from UDS/GOODS-S/HUDF and the recent work of Stefanon et al. (2021) from GREATS at $z > 6$. Mass incomplete measurements are not shown. Upper limits for empty bins are shown by the horizontal gray line with an arrow.</figcaption></figure>

as well as a reduced constraint on rare, high-mass systems found only in larger volumes. They are, however, completely independent as Grazian et al. (2015) do not include the CANDELS-COSMOS field. Where a comparison is possible ($3.5 < z \le 7.5$), our results are consistent with those of Grazian et al. (2015) within the stated uncertainties. Similarly at $z > 5.5$, we are able to compare with the recent measurements of the low-to-intermediate mass regime from Stefanon et al. (2021), with which the present work also is consistent within the stated uncertainties.

We do not include comparisons to stellar mass functions inferred from $L_{UV}$-selected samples with stellar masses estimated from their UV luminosities (e.g., Song et al. 2016; Harikane et al. 2016) as UV-bright sources are only a portion of the general galaxy population and conversations between $L_{UV}$ and $M$ are hazardous. Moreover, such studies often rely on color-color selections which are less certain than photo-z, in addition to being susceptible to a number of systematic effects when translating a UV luminosity function to a SMF by means of a $L_{UV}-M$ relation. We refer the reader to Davidzon et al. (2017) for basic comparisons and discussion.

While the smaller area surveys used by Grazian et al. (2015) and Stefanon et al. (2021) are more mass complete at lower masses, COSMOS2020 contains the widest contiguous NIR imaging at these depths and consequently provides a larger volume than any previous study (including Davidzon et al. 2017) and so introduces new constraints on the rarest systems at all redshifts $z \le 7.5$ such as massive galaxies. Yet it appears that we are quickly approaching the statistical limit beyond which a larger survey is needed to find more massive systems, if they exist (e.g., Euclid, see discussions in McPartland et al., in prep.). The nature of these sources, their authenticity, and their potential Eddington bias is further assessed in Sect. 6.3.

## 5.2. The star-forming and quiescent galaxy stellar mass functions

In Sect. 3.2, we used an *NUVrJ* color-color selection to distinguish quiescent systems from star-forming ones and then estimated their respective mass completeness limits, which will differ since their $M/L$ ratios are not the same. The corresponding SMFs are shown in Fig. 5 and compared to literature in Fig. 7.
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Evolution of the galaxy stellar mass function across 12 redshift bins ($0.2 < z \le 7.5$). The $0.2 < z \le 0.5$ SMF from the first redshift bin is repeated in each panel for reference shown by the purple dotted curve. Two other COSMOS stellar mass functions from Ilbert et al. (2013) and Davidzon et al. (2017) are shown for comparison, along with Grazian et al. (2015) from UDS/GOODS-S/HUDF and the recent work of Stefanon et al. (2021) from GREATS at $z > 6$. Mass incomplete measurements are not shown. Upper limits for empty bins are shown by the horizontal gray line with an arrow.</figcaption></figure>

as well as a reduced constraint on rare, high-mass systems found only in larger volumes. They are, however, completely independent as Grazian et al. (2015) do not include the CANDELS-COSMOS field. Where a comparison is possible ($3.5 < z \le 7.5$), our results are consistent with those of Grazian et al. (2015) within the stated uncertainties. Similarly at $z > 5.5$, we are able to compare with the recent measurements of the low-to-intermediate mass regime from Stefanon et al. (2021), with which the present work also is consistent within the stated uncertainties.

We do not include comparisons to stellar mass functions inferred from $L_{UV}$-selected samples with stellar masses estimated from their UV luminosities (e.g., Song et al. 2016; Harikane et al. 2016) as UV-bright sources are only a portion of the general galaxy population and conversations between $L_{UV}$ and $M$ are hazardous. Moreover, such studies often rely on color-color selections which are less certain than photo-z, in addition to being susceptible to a number of systematic effects when translating a UV luminosity function to a SMF by means of a $L_{UV}-M$ relation. We refer the reader to Davidzon et al. (2017) for basic comparisons and discussion.

While the smaller area surveys used by Grazian et al. (2015) and Stefanon et al. (2021) are more mass complete at lower masses, COSMOS2020 contains the widest contiguous NIR imaging at these depths and consequently provides a larger volume than any previous study (including Davidzon et al. 2017) and so introduces new constraints on the rarest systems at all redshifts $z \le 7.5$ such as massive galaxies. Yet it appears that we are quickly approaching the statistical limit beyond which a larger survey is needed to find more massive systems, if they exist (e.g., Euclid, see discussions in McPartland et al., in prep.). The nature of these sources, their authenticity, and their potential Eddington bias is further assessed in Sect. 6.3.

## 5.2. The star-forming and quiescent galaxy stellar mass functions

In Sect. 3.2, we used an *NUVrJ* color-color selection to distinguish quiescent systems from star-forming ones and then estimated their respective mass completeness limits, which will differ since their $M/L$ ratios are not the same. The corresponding SMFs are shown in Fig. 5 and compared to literature in Fig. 7.
<figure><img src="image_7.png" /><figcaption>**Fig. 7.** Evolution of the star-forming (blue) and quiescent (orange) galaxy components of the galaxy stellar mass function in nine bins of redshift (0.2 < z < 4.5). Uncertainty envelopes correspond to 1 and 2σ limits. For comparison, we show other literature studies in COSMOS from Ilbert et al. (2013) and Davidzon et al. (2017) that adopt similar selections and methodologies. For reference, the total SMF is shown in each bin (solid gray) and the 0.2 < z < 0.5 total SMF is repeated in each panel (purple dotted). Mass incomplete measurements as defined by the channel 1 limiting magnitude are not shown, with the mass limits corresponding to the $K_s$ limiting magnitude shown by orange arrows. Upper limits for empty bins are shown by the horizontal gray line and arrow.</figcaption></figure>

Fitting results are discussed in Sect. 5.3 and details can be found in Appendix C; their corresponding Schechter parameters are reported in Tables C.2 and C.3 for the star-forming and quiescent samples, respectively. Fractional uncertainties on mass and cosmic variance are adopted from the total sample, leaving only the impact of poisson noise differentiated between the star-forming and quiescent samples (see Sect. 3.5).

We can follow the development of quiescent systems out to $z \approx 5$, although with significant uncertainty at $z \gtrsim 4$. As evi-

denced by Fig. 1, only ~200 of quiescent systems are found at $z \gtrsim 4$, dropping precipitously to only one candidate by $z \approx 6$. This is partly driven by mass completeness. At $z < 3$ the difference between the IRAC channel 1 mass limit and the effective mass completeness dictated by our *izYJHKs* selection function is ≤0.2 dex in mass. This difference grows to ~0.4 dex at $z \sim 4.5$ when $K_s$ falls blueward of the Balmer break causing the mass completeness to shift upwards to higher masses, indicated by the hatched region of the SMF of $4.5 < z \le 5.5$ quiescent systems
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Evolution of the star-forming (blue) and quiescent (orange) galaxy components of the galaxy stellar mass function in nine bins of redshift (0.2 < z < 4.5). Uncertainty envelopes correspond to 1 and 2σ limits. For comparison, we show other literature studies in COSMOS from Ilbert et al. (2013) and Davidzon et al. (2017) that adopt similar selections and methodologies. For reference, the total SMF is shown in each bin (solid gray) and the 0.2 < z < 0.5 total SMF is repeated in each panel (purple dotted). Mass incomplete measurements as defined by the channel 1 limiting magnitude are not shown, with the mass limits corresponding to the $K_s$ limiting magnitude shown by orange arrows. Upper limits for empty bins are shown by the horizontal gray line and arrow.</figcaption></figure>

Fitting results are discussed in Sect. 5.3 and details can be found in Appendix C; their corresponding Schechter parameters are reported in Tables C.2 and C.3 for the star-forming and quiescent samples, respectively. Fractional uncertainties on mass and cosmic variance are adopted from the total sample, leaving only the impact of poisson noise differentiated between the star-forming and quiescent samples (see Sect. 3.5).

We can follow the development of quiescent systems out to $z \approx 5$, although with significant uncertainty at $z \gtrsim 4$. As evi-

denced by Fig. 1, only ~200 of quiescent systems are found at $z \gtrsim 4$, dropping precipitously to only one candidate by $z \approx 6$. This is partly driven by mass completeness. At $z < 3$ the difference between the IRAC channel 1 mass limit and the effective mass completeness dictated by our *izYJHKs* selection function is $\le 0.2$ dex in mass. This difference grows to ~0.4 dex at $z \sim 4.5$ when $K_s$ falls blueward of the Balmer break causing the mass completeness to shift upwards to higher masses, indicated by the hatched region of the SMF of $4.5 < z \le 5.5$ quiescent systems
<figure><img src="image_8.png" /><figcaption>**Fig. 8.** Evolution of the volume number density $\Phi$ of the total, star-forming, and quiescent samples at fixed mass, showing only mass complete bins. $1\sigma$ uncertainties are indicated by the colored envelopes. Gray shaded regions indicate upper limits corresponding to $\Phi(N = 0)$ at each redshift bin. The lower panel shows the evolution of the fractional change in the volume number density at fixed mass, relative to $z_0 \equiv 0.2 < z \le 0.5$ of the total sample.</figcaption></figure>

in Fig. 5. At this point the identification of quiescent systems is reliant on only a few bands, and is dependent on the particular SED templates (as suggested by the two distinct clusters in the upper right corner of the $4.5 < z \le 5.5$ NUVrJ diagram), making such a distinction hazardous and susceptible to interlopers.

Our measurements of the star-forming and quiescent SMFs are generally in good agreement with other literature measurements in COSMOS, namely Davidzon et al. (2017) and Ilbert et al. (2013) as they are the only NIR-selected, mass-complete samples from which star-forming and quiescent sub-samples are identified by NUVrJ. Other selections may induce additional systematics, and other separation methods (e.g., UVJ, BzK, sSFR) implicitly adopt a different criteria for quiescence (see Davidzon et al. 2018; Leja et al. 2019b). As shown in Fig. 7, our measurements indicate similar high-$M^*$ quiescent number densities compared to Ilbert et al. (2013) and Davidzon et al. (2017), but find significantly more low-mass lowest-$M^*$ quiescent galaxies compared to Davidzon et al. (2017; which are not measurable from the data used by Ilbert et al. 2013). Since our sample is derived from comparably deeper NIR data, it is expected to be complete down to lower masses relative to Davidzon et al. (2017). Given the increased quiescent galaxy number densities near the mass completeness limit in our work, we conclude that the 70–80% completeness threshold of Davidzon et al. (2017) is underestimated by a factor of ~2–10 across this $z$-range, and is more likely only 15–35% complete. We note, however, that this in agreement with worst-case sce-

nario discussed by the authors (see Sect. 4.2 of Davidzon et al. 2017).

The SMFs of star-forming and quiescent galaxies have remarkably different shapes and evolutionary histories. As shown by Fig. 5, star-forming galaxies at $0.2 < z \le 0.5$ follow a double Schechter form (Eq. (7)) with a characteristic stellar mass $M^* \approx 10^{10-11} M_\odot$ and yet by $z \approx 3$ it appears to flatten into a smooth power-law-like form with a lower overall normalization. Meanwhile the SMF of $0.2 < z \le 0.5$ quiescent galaxies follows the form of a double Schechter (Eq. (8)) with a low-mass upturn and a similarly positioned $M^*$ (see Moutard et al. 2016) but the form appears to lose its lower-mass Schechter component around $z \sim 2$. Although this may be physical, the fact that the quiescent sample is less mass complete at a given $z$ means that the low-mass end of the total sample reflects only the contribution from star-forming systems. However, the contribution from low-mass quiescent systems, if they exist, can be expected to be <1%.

This evolutionary picture is perhaps more easily understood by Fig. 8. Here galaxies have been binned by mass instead of redshift, allowing for a more direct comprehension of the growth, or lack thereof, in the number density of galaxies at fixed mass. We first notice a growth rate in the number densities that is strikingly consistent across all masses from $z = 7.5 \rightarrow 1$ (i.e., similar slopes in each mass bin), constituting ~5 Gyr or ~36% of the history of the Universe, consistent with recent findings by Wright et al. (2018). This consistent growth is made especially
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Evolution of the volume number density $\Phi$ of the total, star-forming, and quiescent samples at fixed mass, showing only mass complete bins. $1\sigma$ uncertainties are indicated by the colored envelopes. Gray shaded regions indicate upper limits corresponding to $\Phi(N = 0)$ at each redshift bin. The lower panel shows the evolution of the fractional change in the volume number density at fixed mass, relative to $z_0 \equiv 0.2 < z \le 0.5$ of the total sample.</figcaption></figure>

in Fig. 5. At this point the identification of quiescent systems is reliant on only a few bands, and is dependent on the particular SED templates (as suggested by the two distinct clusters in the upper right corner of the $4.5 < z \le 5.5$ NUVrJ diagram), making such a distinction hazardous and susceptible to interlopers.

Our measurements of the star-forming and quiescent SMFs are generally in good agreement with other literature measurements in COSMOS, namely Davidzon et al. (2017) and Ilbert et al. (2013) as they are the only NIR-selected, mass-complete samples from which star-forming and quiescent sub-samples are identified by NUVrJ. Other selections may induce additional systematics, and other separation methods (e.g., UVJ, BzK, sSFR) implicitly adopt a different criteria for quiescence (see Davidzon et al. 2018; Leja et al. 2019b). As shown in Fig. 7, our measurements indicate similar high-$M^*$ quiescent number densities compared to Ilbert et al. (2013) and Davidzon et al. (2017), but find significantly more low-mass lowest-$M^*$ quiescent galaxies compared to Davidzon et al. (2017; which are not measurable from the data used by Ilbert et al. 2013). Since our sample is derived from comparably deeper NIR data, it is expected to be complete down to lower masses relative to Davidzon et al. (2017). Given the increased quiescent galaxy number densities near the mass completeness limit in our work, we conclude that the 70–80% completeness threshold of Davidzon et al. (2017) is underestimated by a factor of ~2–10 across this $z$-range, and is more likely only 15–35% complete. We note, however, that this in agreement with worst-case sce-

nario discussed by the authors (see Sect. 4.2 of Davidzon et al. 2017).

The SMFs of star-forming and quiescent galaxies have remarkably different shapes and evolutionary histories. As shown by Fig. 5, star-forming galaxies at $0.2 < z \le 0.5$ follow a double Schechter form (Eq. (7)) with a characteristic stellar mass $M^* \approx 10^{10-11} M_\odot$ and yet by $z \approx 3$ it appears to flatten into a smooth power-law-like form with a lower overall normalization. Meanwhile the SMF of $0.2 < z \le 0.5$ quiescent galaxies follows the form of a double Schechter (Eq. (8)) with a low-mass upturn and a similarly positioned $M^*$ (see Moutard et al. 2016) but the form appears to lose its lower-mass Schechter component around $z \sim 2$. Although this may be physical, the fact that the quiescent sample is less mass complete at a given $z$ means that the low-mass end of the total sample reflects only the contribution from star-forming systems. However, the contribution from low-mass quiescent systems, if they exist, can be expected to be <1%.

This evolutionary picture is perhaps more easily understood by Fig. 8. Here galaxies have been binned by mass instead of redshift, allowing for a more direct comprehension of the growth, or lack thereof, in the number density of galaxies at fixed mass. We first notice a growth rate in the number densities that is strikingly consistent across all masses from $z = 7.5 \rightarrow 1$ (i.e., similar slopes in each mass bin), constituting ~5 Gyr or ~36% of the history of the Universe, consistent with recent findings by Wright et al. (2018). This consistent growth is made especially
<figure><img src="image_9.png" /><figcaption>Fig. 9. Evolution of the fraction of quiescent galaxies as a function of redshift in four bins spanning $10^{9.5} < M \le 10^{11.5} M_{\odot}$. Uncertainty envelopes correspond to 1 and $2\sigma$ limits. Points which are mass incomplete are not shown.</figcaption></figure>

clear in the lower left panel where the growth is relative to that of the total sample at $z_0 \equiv 0.2 < z \le 0.5$. Although the uncertainties are considerable, there may be a hint that the growth is fastest (i.e., the slope is higher) for systems at $M \sim 10^{10.5-11.0} M_{\odot}$ or at the ‘knee’ of the mass function where star formation is hypothesized to be the most efficient (Moster et al. 2013; Behroozi et al. 2013; Wechsler & Tinker 2018). However, we note that Eddington bias remains uncorrected and so a growth in the bias strength with redshift may produce this apparently slow evolution (as suggested by Fig. 4). While the least massive systems are always more common than more massive ones, this is simply a consequence of the monotonic shape of the SMF.

Another interesting feature is that the number density of $M = 10^{9.5-10.0}$ and $M = 10^{10.0-10.5} M_{\odot}$ systems appear to decrease mildly at $z \sim 2.5$. There are fewer such systems in this work relative to both Davidzon et al. (2017) and Ilbert et al. (2013), and so this may be indicative of an incompleteness in our selection, or alternatively a bias in $z$ and/or $M$. Determining the origin of this difference is nontrivial, and so we simply caution against over-interpretation of this feature.

The highest-$z$ constraints are difficult to interpret due to the incompleteness of low-mass systems (which are omitted) and the rarity of the most massive sources. The fact that the constraints at $z \ge 5$ overlap with the $\Phi(N=0)$ upper limit region in gray indicates that even the comparably large, deep volume of COSMOS is insufficient to provide robust measurements of the number density of massive galaxies at these early times. Future constraints will be derived from even larger volumes, which are expected to be delivered soon by near-infrared wide-area deep surveys such as *Roman* and *Euclid*.

This nearly uniform growth in the number density of sources slows down at $z \le 1$ in all but the least massive bin, with almost no growth for $M > 10^{10.5} M_{\odot}$ systems. To understand this further, we examine the growth of the number densities of the star-forming and quiescent subsamples. As seen in the middle column of Fig. 8, star-forming systems maintain their monotonic mass-ranked order. In addition, their number density grows rapidly until $z \approx 1$, when it starts to decrease. As shown in the

rightmost column, quiescent galaxies follow a different pattern. Instead of a monotonic, mass ranked growth in number density, quiescent systems around the knee at $10^{10.5} M_{\odot}$ are always the most numerous, with a slower rate of growth compared to the least massive and most massive bins. While the evidence for massive quiescent systems at $z > 2$ is hampered by the limited volume of COSMOS, we can confidently see that they appear to grow quickly in number at $z < 2$. By $z \approx 1$ the number density of the most massive ones ($M = 10^{11.5-12.0} M_{\odot}$) stabilizes at $\approx 3 \times 10^{-5} Mpc^{-3} \text{ dex}^{-1}$. Reasons for this are discussed in Sect. 6.3.

Figure 9 shows the evolution of the fractional contribution of quiescent systems in four 0.5 dex bins at fixed mass and one wider 2 dex bin. The quiescent fraction of low mass systems increases with time monotonically, and at a higher rate of growth than those of the highest mass, at least for $z < 2$ where comparisons can be made. The quiescent fraction of the most massive systems at any redshift is more affected by both Eddington bias and poisson noise compared to low-mass systems and so should be taken as an upper limit. Nevertheless, it seems that between 20 and 30% of $M > 10^{11} M_{\odot}$ galaxies are quiescent from $z \approx 5 \rightarrow 1$, growing above 50% and plateauing at $z < 1$.

### 5.3. Inferring the intrinsic mass function

So far we have accounted for three key sources of statistical uncertainty in Sect. 3.5. However, Eddington bias remains a source of systematic uncertainty that has not yet been removed.

To do so, we fit the observed total, star-forming, and quiescent SMFs with Schechter functions convolved with a kernel which attempts to describe the mass uncertainty for a given mass and redshift interval. This approach is common in the literature (Ilbert et al. 2013; Davidzon et al. 2017; Adams et al. 2021). We adopt a two component parametric kernel of the form introduced by Ilbert et al. (2013):

$$ \mathcal{D}(M) = \frac{1}{2\pi} \exp\left(\frac{-M^2}{2\sigma_{\text{Edd}}^2}\right) \times \frac{\tau_{\text{Edd}}}{2\pi} \frac{1}{(\tau_{\text{Edd}}/2)^2 + M^2} \quad (9) $$

which describes a core Gaussian component with a constant width $\sigma_{\text{Edd}}$ in product with a Lorenztian component which provides wide wings that are crucial to capturing catastrophic errors in SED measurements. The other free parameter $\tau_{\text{Edd}}$ is redshift dependent such that $\tau_{\text{Edd}} = \tau_c(1+z)$, widening the wings with increasing $z$. Instead of fitting the kernel directly to the $\mathcal{L}(M|z)$ distributions (as in Davidzon et al. 2017; Adams et al. 2021), we determine the kernel parameters from directly fitting the total SMF leaving the kernel free to vary independently at each redshift. However, both $\sigma_{\text{Edd}}$ and $\tau_c$ are nuisance parameters and leaving them free to vary is may produce over-fitting of any one SMF at a given $z$. As expected, the resulting parameter distributions with respect to $z$ contain outliers on either side of the median particularly where the kernel shape is most degenerate with that of the intrinsic SMF; we adopt the median values $\sigma_{\text{Edd}} = 0.2$ and $\tau_c = 0.1$. The resulting kernel is significantly wider than the $\mathcal{L}(M|z)$ distributions shown earlier in Fig. 4. This is because the latter assume perfect redshifts and therefore underestimate the full $M$ uncertainty as there is likely to be a covariance with $z$. In addition, there is a mild evolution of $\mathcal{L}(M|z)$ with $M$, which we have omitted from our kernel form for the sake of simplicity and to avoid overfitting. Interestingly, while $\sigma_{\text{Edd}} \approx 0.2$ is smaller than that found by Davidzon et al. (2017)
<figure><img src="image_10.png" /><figcaption>Fig. 9. Evolution of the fraction of quiescent galaxies as a function of redshift in four bins spanning $10^{9.5} < M \le 10^{11.5} M_{\odot}$. Uncertainty envelopes correspond to 1 and $2\sigma$ limits. Points which are mass incomplete are not shown.</figcaption></figure>

clear in the lower left panel where the growth is relative to that of the total sample at $z_0 \equiv 0.2 < z \le 0.5$. Although the uncertainties are considerable, there may be a hint that the growth is fastest (i.e., the slope is higher) for systems at $M \sim 10^{10.5-11.0} M_{\odot}$ or at the ‘knee’ of the mass function where star formation is hypothesized to be the most efficient (Moster et al. 2013; Behroozi et al. 2013; Wechsler & Tinker 2018). However, we note that Eddington bias remains uncorrected and so a growth in the bias strength with redshift may produce this apparently slow evolution (as suggested by Fig. 4). While the least massive systems are always more common than more massive ones, this is simply a consequence of the monotonic shape of the SMF.

Another interesting feature is that the number density of $M = 10^{9.5-10.0}$ and $M = 10^{10.0-10.5} M_{\odot}$ systems appear to decrease mildly at $z \sim 2.5$. There are fewer such systems in this work relative to both Davidzon et al. (2017) and Ilbert et al. (2013), and so this may be indicative of an incompleteness in our selection, or alternatively a bias in $z$ and/or $M$. Determining the origin of this difference is nontrivial, and so we simply caution against over-interpretation of this feature.

The highest-$z$ constraints are difficult to interpret due to the incompleteness of low-mass systems (which are omitted) and the rarity of the most massive sources. The fact that the constraints at $z \ge 5$ overlap with the $\Phi(N=0)$ upper limit region in gray indicates that even the comparably large, deep volume of COSMOS is insufficient to provide robust measurements of the number density of massive galaxies at these early times. Future constraints will be derived from even larger volumes, which are expected to be delivered soon by near-infrared wide-area deep surveys such as *Roman* and *Euclid*.

This nearly uniform growth in the number density of sources slows down at $z \le 1$ in all but the least massive bin, with almost no growth for $M > 10^{10.5} M_{\odot}$ systems. To understand this further, we examine the growth of the number densities of the star-forming and quiescent subsamples. As seen in the middle column of Fig. 8, star-forming systems maintain their monotonic mass-ranked order. In addition, their number density grows rapidly until $z \approx 1$, when it starts to decrease. As shown in the

rightmost column, quiescent galaxies follow a different pattern. Instead of a monotonic, mass ranked growth in number density, quiescent systems around the knee at $10^{10.5} M_{\odot}$ are always the most numerous, with a slower rate of growth compared to the least massive and most massive bins. While the evidence for massive quiescent systems at $z > 2$ is hampered by the limited volume of COSMOS, we can confidently see that they appear to grow quickly in number at $z < 2$. By $z \approx 1$ the number density of the most massive ones ($M = 10^{11.5-12.0} M_{\odot}$) stabilizes at $\approx 3 \times 10^{-5} \text{ Mpc}^{-3} \text{ dex}^{-1}$. Reasons for this are discussed in Sect. 6.3.

Figure 9 shows the evolution of the fractional contribution of quiescent systems in four 0.5 dex bins at fixed mass and one wider 2 dex bin. The quiescent fraction of low mass systems increases with time monotonically, and at a higher rate of growth than those of the highest mass, at least for $z < 2$ where comparisons can be made. The quiescent fraction of the most massive systems at any redshift is more affected by both Eddington bias and poisson noise compared to low-mass systems and so should be taken as an upper limit. Nevertheless, it seems that between 20 and 30% of $M > 10^{11} M_{\odot}$ galaxies are quiescent from $z \approx 5 \rightarrow 1$, growing above 50% and plateauing at $z < 1$.

### 5.3. Inferring the intrinsic mass function

So far we have accounted for three key sources of statistical uncertainty in Sect. 3.5. However, Eddington bias remains a source of systematic uncertainty that has not yet been removed.

To do so, we fit the observed total, star-forming, and quiescent SMFs with Schechter functions convolved with a kernel which attempts to describe the mass uncertainty for a given mass and redshift interval. This approach is common in the literature (Ilbert et al. 2013; Davidzon et al. 2017; Adams et al. 2021). We adopt a two component parametric kernel of the form introduced by Ilbert et al. (2013):

$$ \mathcal{D}(M) = \frac{1}{2\pi} \exp\left(\frac{-M^2}{2\sigma_{\text{Edd}}^2}\right) \times \frac{\tau_{\text{Edd}}}{2\pi} \frac{1}{(\tau_{\text{Edd}}/2)^2 + M^2} \quad (9) $$

which describes a core Gaussian component with a constant width $\sigma_{\text{Edd}}$ in product with a Lorenztian component which provides wide wings that are crucial to capturing catastrophic errors in SED measurements. The other free parameter $\tau_{\text{Edd}}$ is redshift dependent such that $\tau_{\text{Edd}} = \tau_c(1+z)$, widening the wings with increasing $z$. Instead of fitting the kernel directly to the $\mathcal{L}(M|z)$ distributions (as in Davidzon et al. 2017; Adams et al. 2021), we determine the kernel parameters from directly fitting the total SMF leaving the kernel free to vary independently at each redshift. However, both $\sigma_{\text{Edd}}$ and $\tau_c$ are nuisance parameters and leaving them free to vary is may produce over-fitting of any one SMF at a given $z$. As expected, the resulting parameter distributions with respect to $z$ contain outliers on either side of the median particularly where the kernel shape is most degenerate with that of the intrinsic SMF; we adopt the median values $\sigma_{\text{Edd}} = 0.2$ and $\tau_c = 0.1$. The resulting kernel is significantly wider than the $\mathcal{L}(M|z)$ distributions shown earlier in Fig. 4. This is because the latter assume perfect redshifts and therefore underestimate the full $M$ uncertainty as there is likely to be a covariance with $z$. In addition, there is a mild evolution of $\mathcal{L}(M|z)$ with $M$, which we have omitted from our kernel form for the sake of simplicity and to avoid overfitting. Interestingly, while $\sigma_{\text{Edd}} \approx 0.2$ is smaller than that found by Davidzon et al. (2017)
<figure><img src="image_10.png" /><figcaption>**Fig. 10.** Summary of inferred galaxy stellar mass function evolution. Best-fit parameters are estimated by a Markov chain Monte Carlo fitting with a Schechter function at fixed redshift convolved with a redshift dependent kernel from which the inferred intrinsic mass function is recovered with the unconvolved Schechter fit (colored solid curves) for the total (leftmost), star-forming (center left), and quiescent (center right) samples. Estimates are shown corresponding to both the median posterior parameter values including $1\sigma$ envelopes (upper row) and the single set of maximum likelihood parameters (bottom row). The rightmost panel shows the results of fitting the Double Schechter continuity model of Leja et al. (2019a) to the total SMF between $0.2 < z \le 3.0$ where the fits are reliable; see later in Sect. 5.3 for details.</figcaption></figure>

(0.35), $\tau_c \approx 0.1$ is conspicuously larger (0.04, same as in Ilbert et al. 2013). If we instead fit to the $L(M|z)$ distributions directly, we find $\sigma_{Edd} \approx 0.05$ and $\tau_c \approx 0.04$. These are much smaller than (Davidzon et al. 2017), which may be explained in part by the slightly different approach they used to estimate $L(M)$ that more directly incorporates uncertainties on photo-$z$ (see Sect. 4.1 of Davidzon et al. 2017). Nevertheless, we opt for the pessimistic case and fix the kernel shape and evolution with $z$ to $\sigma_{Edd} \approx 0.2$ and $\tau_c \approx 0.1$ for the subsequent analysis of the total, as well as the star-forming and quiescent SMFs; their shapes are shown in the inset panels in Figs. C.1–C.3, respectively.

At the present, precisely determining the correct, intrinsic shape of the Eddington bias, and its evolution with $z$ as well as $M$ is problematically difficult. In addition, the results stemming from these kernel-convolved fits suffer a degree of degeneracy with the kernel parameters (fixed or unfixed). We are not alone in issuing this caution; although the SMF measurements of Davidzon et al. (2017) is similar to those of Grazian et al. (2015), their different choice of convolution kernel caused their inferred Schechter parameters to differ. Recently, Adams et al. (2021) explored the impact of the assumed shape of the kernel, as well as various other systematics at $z < 2$, highlighting the full extent of the problem. Pushing measurements of the SMF to higher redshifts, where fewer constraints are available, naturally increases the influence of the kernel shape, and so the results derived here for $\alpha, M^*$, and $\Phi^*$ (Tables C.1–C.3) should be taken in conjunction with our assumed kernel and its evolution with $z$.

When fitting each sample (total, star-forming, and quiescent) we assume a double Schechter parametric form (Eq. (8)) and move to a single Schechter form (Eq. (7)) when it achieves a lower $\chi^2$ per degree of freedom. The point at which this change occurs, and when various parameters are fixed, are different

for each of the three samples as it depends on not only their shape but also $M_{lim}$. The fitted points follow from before: we account for minor incompleteness on the low-mass end using the $1/V_{max}$ correction and include only mass complete $M$-bins adopting the uncertainty budget from Sect. 3.5. We proceed in two stages, fitting first using a simple and efficient $\chi^2$ minimization routine (MINUIT, James & Roos 1975) whose resulting best-fit parameters are used to set the initial positions of walkers in a secondary and longer Markov chain Monte Carlo (MCMC) optimization (EMCEE, Foreman-Mackey et al. 2013). It is generally appropriate to initialize walkers at well defined initial positions, assuming each chain is linearly independent (in this case scattering them by 1% of the initial parameter values) and allowed to converge well past its respective autocorrelation length such that it is not sensitive to those initial conditions (Hogg et al. 2010). Flat priors are used for each parameter, with generous limits that are typically not encountered during a given fit. 500 MCMC walkers are initialized to seek the state of maximum likelihood derived similarly to the minimum $\chi^2$, and do so following the “stretch move” (Goodman & Weare 2010) until converged, defined as having every chain iterate for at least $10 \times$ their autocorrelation length, and every autocorrelation length having changed by <1% of their previous value. We do not see any significant differences by using a different move style (e.g., Red-Blue, Metropolis-Hastings), suggesting that the data provide sufficient constraining power. Although we include the original $\chi^2$ results throughout, we focus on the MCMC results which provide full posterior sampling, which are taken from the last 90% of each chain to avoid only the initial burn-in when the parameters are only beginning to converge. The MCMC and $\chi^2$ methods generally find consistent results. Fitting solutions corresponding to both the median posterior and maximum likelihood for the total, star-forming, and quiescent samples is shown in Fig. 10.
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Summary of inferred galaxy stellar mass function evolution. Best-fit parameters are estimated by a Markov chain Monte Carlo fitting with a Schechter function at fixed redshift convolved with a redshift dependent kernel from which the inferred intrinsic mass function is recovered with the unconvolved Schechter fit (colored solid curves) for the total (leftmost), star-forming (center left), and quiescent (center right) samples. Estimates are shown corresponding to both the median posterior parameter values including $1\sigma$ envelopes (upper row) and the single set of maximum likelihood parameters (bottom row). The rightmost panel shows the results of fitting the Double Schechter continuity model of Leja et al. (2019a) to the total SMF between $0.2 < z \le 3.0$ where the fits are reliable; see later in Sect. 5.3 for details.</figcaption></figure>

(0.35), $\tau_c \approx 0.1$ is conspicuously larger (0.04, same as in Ilbert et al. 2013). If we instead fit to the $\mathcal{L}(M|z)$ distributions directly, we find $\sigma_{Edd} \approx 0.05$ and $\tau_c \approx 0.04$. These are much smaller than (Davidzon et al. 2017), which may be explained in part by the slightly different approach they used to estimate $\mathcal{L}(M)$ that more directly incorporates uncertainties on photo-$z$ (see Sect. 4.1 of Davidzon et al. 2017). Nevertheless, we opt for the pessimistic case and fix the kernel shape and evolution with $z$ to $\sigma_{Edd} \approx 0.2$ and $\tau_c \approx 0.1$ for the subsequent analysis of the total, as well as the star-forming and quiescent SMFs; their shapes are shown in the inset panels in Figs. C.1–C.3, respectively.

At the present, precisely determining the correct, intrinsic shape of the Eddington bias, and its evolution with $z$ as well as $M$ is problematically difficult. In addition, the results stemming from these kernel-convolved fits suffer a degree of degeneracy with the kernel parameters (fixed or unfixed). We are not alone in issuing this caution; although the SMF measurements of Davidzon et al. (2017) is similar to those of Grazian et al. (2015), their different choice of convolution kernel caused their inferred Schechter parameters to differ. Recently, Adams et al. (2021) explored the impact of the assumed shape of the kernel, as well as various other systematics at $z < 2$, highlighting the full extent of the problem. Pushing measurements of the SMF to higher redshifts, where fewer constraints are available, naturally increases the influence of the kernel shape, and so the results derived here for $\alpha, M^*$, and $\Phi^*$ (Tables C.1–C.3) should be taken in conjunction with our assumed kernel and its evolution with $z$.

When fitting each sample (total, star-forming, and quiescent) we assume a double Schechter parametric form (Eq. (8)) and move to a single Schechter form (Eq. (7)) when it achieves a lower $\chi^2$ per degree of freedom. The point at which this change occurs, and when various parameters are fixed, are different

for each of the three samples as it depends on not only their shape but also $M_{lim}$. The fitted points follow from before: we account for minor incompleteness on the low-mass end using the $1/V_{max}$ correction and include only mass complete $M$-bins adopting the uncertainty budget from Sect. 3.5. We proceed in two stages, fitting first using a simple and efficient $\chi^2$ minimization routine (MINUIT, James & Roos 1975) whose resulting best-fit parameters are used to set the initial positions of walkers in a secondary and longer Markov chain Monte Carlo (MCMC) optimization (EMCEE, Foreman-Mackey et al. 2013). It is generally appropriate to initialize walkers at well defined initial positions, assuming each chain is linearly independent (in this case scattering them by 1% of the initial parameter values) and allowed to converge well past its respective autocorrelation length such that it is not sensitive to those initial conditions (Hogg et al. 2010). Flat priors are used for each parameter, with generous limits that are typically not encountered during a given fit. 500 MCMC walkers are initialized to seek the state of maximum likelihood derived similarly to the minimum $\chi^2$, and do so following the “stretch move” (Goodman & Weare 2010) until converged, defined as having every chain iterate for at least $10 \times$ their autocorrelation length, and every autocorrelation length having changed by <1% of their previous value. We do not see any significant differences by using a different move style (e.g., Red-Blue, Metropolis-Hastings), suggesting that the data provide sufficient constraining power. Although we include the original $\chi^2$ results throughout, we focus on the MCMC results which provide full posterior sampling, which are taken from the last 90% of each chain to avoid only the initial burn-in when the parameters are only beginning to converge. The MCMC and $\chi^2$ methods generally find consistent results. Fitting solutions corresponding to both the median posterior and maximum likelihood for the total, star-forming, and quiescent samples is shown in Fig. 10.
<figure><img src="image_11.png" /><figcaption>**Fig. 11.** Schechter parameter evolution with redshift. Left: evolution of best-fit Schechter parameters computed from fits to the total SMF from $\chi^2$ regression (red diamonds), as well as likelihood methods using fixed redshifts bins (blue; median posteriors as filled points with 68% error bars, maximum likelihood parameter set as unfilled squares), and our continuity model fit (maroon curve; median posterior). Error bars are not shown when parameters are fixed. Center: same points as on the left panels compared to literature values from Davidzon et al. (2017, gray points with hatching) and Leja et al. (2019a, gray dashed curves). Right: evolution of the best-fit Schechter parameters computed from fits to the star-forming (blue) and quiescent (orange) SMFs from the $\chi^2$ regression (unfilled colored diamonds) and likelihood methods using fixed redshift bins (median posterior as filled points with $\pm$34% envelope, maximum likelihood as unfilled squares).</figcaption></figure>

An alternative version binned in redshift (as opposed to total, star-forming, or quiescent sample) is included for reference in Fig. C.4 and discussed in Appendix C.

The evolution of the Schechter parameters inferred from the total SMF are shown in the leftmost panels Fig. 11, and compared with (Davidzon et al. 2017) in the center panels. Table C.1 records the best-fit parameters, with detailed fits shown in Fig. C.1.

As evidenced by the low $\chi^2$ values, the double Schechter well describes the observed SMF at $z \le 3$. No parameters are fixed other than those of the kernel. While a single Schechter finds a better fit at $z > 3$, the increasing mass completeness limit weakens the constraints on $\alpha_1$. To avoid overfitting, we fix $\alpha_1$ to the value from the best-fit value at $2.5 < z \le 3.0$ and do so for the $\chi^2$ and MCMC methods independently but to excellent agreement. All of the unfixed search parameters are generally bounded such that the chains are not likely to encounter

them. However, for the sake of physicality especially at high-z, we choose to bound the search space for the normalization $\Phi(z > 0.5)$ (both $\Phi_1$ and $\Phi_2$ where applicable) to below the upper 68% on the posterior distribution of $\max(\Phi_1, \Phi_2)$ of the previous mass bin.

It is important to note that from $\chi^2$ minimization we obtain a single set of parameters that have minimized the $\chi^2$ as well as their symmetric, Gaussian uncertainties. On the other hand, MCMC provides only posterior distributions that can be used to estimate parameter uncertainties but do not imply a unique definition of “best-fit parameters”. Commonly, the median is the best-fit statistic of choice, bounded by a 68% one-parameter confidence interval (which ignores covariance). However, for highly skewed posteriors the median may lie out in the wings and is therefore not a typical value for that parameter. In this case, the most obvious choice may be the parameters corresponding to the model which has found the maximum likelihood.
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Schechter parameter evolution with redshift. Left: evolution of best-fit Schechter parameters computed from fits to the total SMF from $\chi^2$ regression (red diamonds), as well as likelihood methods using fixed redshifts bins (blue; median posteriors as filled points with 68% error bars, maximum likelihood parameter set as unfilled squares), and our continuity model fit (maroon curve; median posterior). Error bars are not shown when parameters are fixed. Center: same points as on the left panels compared to literature values from Davidzon et al. (2017, gray points with hatching) and Leja et al. (2019a, gray dashed curves). Right: evolution of the best-fit Schechter parameters computed from fits to the star-forming (blue) and quiescent (orange) SMFs from the $\chi^2$ regression (unfilled colored diamonds) and likelihood methods using fixed redshift bins (median posterior as filled points with $\pm$34% envelope, maximum likelihood as unfilled squares).</figcaption></figure>

An alternative version binned in redshift (as opposed to total, star-forming, or quiescent sample) is included for reference in Fig. C.4 and discussed in Appendix C.

The evolution of the Schechter parameters inferred from the total SMF are shown in the leftmost panels Fig. 11, and compared with (Davidzon et al. 2017) in the center panels. Table C.1 records the best-fit parameters, with detailed fits shown in Fig. C.1.

As evidenced by the low $\chi^2$ values, the double Schechter well describes the observed SMF at $z \le 3$. No parameters are fixed other than those of the kernel. While a single Schechter finds a better fit at $z > 3$, the increasing mass completeness limit weakens the constraints on $\alpha_1$. To avoid overfitting, we fix $\alpha_1$ to the value from the best-fit value at $2.5 < z \le 3.0$ and do so for the $\chi^2$ and MCMC methods independently but to excellent agreement. All of the unfixed search parameters are generally bounded such that the chains are not likely to encounter

them. However, for the sake of physicality especially at high-z, we choose to bound the search space for the normalization $\Phi(z > 0.5)$ (both $\Phi_1$ and $\Phi_2$ where applicable) to below the upper 68% on the posterior distribution of $\max(\Phi_1, \Phi_2)$ of the previous mass bin.

It is important to note that from $\chi^2$ minimization we obtain a single set of parameters that have minimized the $\chi^2$ as well as their symmetric, Gaussian uncertainties. On the other hand, MCMC provides only posterior distributions that can be used to estimate parameter uncertainties but do not imply a unique definition of “best-fit parameters”. Commonly, the median is the best-fit statistic of choice, bounded by a 68% one-parameter confidence interval (which ignores covariance). However, for highly skewed posteriors the median may lie out in the wings and is therefore not a typical value for that parameter. In this case, the most obvious choice may be the parameters corresponding to the model which has found the maximum likelihood.