The JWST Extragalactic Mock Catalog: Modeling Galaxy Populations from the UV through the Near-IR over 13 Billion Years of Cosmic History

Here we describe the procedure we follow to generate galaxy number counts (i.e., the expected number of galaxies of a given mass, at fixed redshift and on-sky area). We first model the evolution of stellar mass functions across cosmic time using continuously evolving Schechter functions for both star-forming and quiescent galaxies. We then generate the expected number of star-forming or quiescent galaxies for a given redshift bin over a given survey area by integrating their respective model mass function, multiplying by the co-moving volume, and drawing from a Poisson distribution with this mean. By computing the cumulative distribution function (CDF) of the stellar mass function, we can then effectively draw this number of galaxies from the mass function using inverse transform sampling.

For star-forming galaxies at z ≳ 4, stellar mass function measurements become increasingly difficult and uncertain. With current facilities, this epoch represents a transition to rest-frame UV selections tracing young stars (with HST) instead of rest-frame optical selections that trace stellar mass (which would require Spitzer/IRAC whose sensitivity is lower). At z ≳ 4 the UV luminosity function becomes much more easily measurable than the stellar mass function with current facilities. Therefore, at redshifts z ≳ 4 we rely on observed UV luminosity functions to constrain the number densities of star-forming galaxies. (Quiescent galaxies at z > 4 are instead based on an informed extrapolation, which is discussed in Section 4).

While generating galaxy counts from a stellar mass function is straightforward, using a UV luminosity function requires modeling a theoretical or an empirical relation linking a galaxy's UV luminosity, or ${M}_{{\rm{UV}}}$, to its stellar mass—hereafter ${M}_{\star }=\mathrm{log}(M/{M}_{\odot })$. The observed connection between UV luminosity and stellar mass is not a simple monotonic relationship; galaxies exhibit a range of UV luminosities at fixed stellar mass that likely depends on other galaxy properties, including stellar population age and metallicity, dust, and gas content. We will describe this distribution in terms of a Gaussian scatter about an average ${M}_{{\rm{UV}}}$–${M}_{\star }$ relation, where the standard deviation of ${M}_{{\rm{UV}}}$ at fixed ${M}_{\star }$ is given by ${\sigma }_{{uv}}$ (assumed independent of ${M}_{\star }$). We can then express the probability of a galaxy of stellar mass ${M}_{\star }$ to have a given ${M}_{{\rm{UV}}}$ as

$$\begin{eqnarray}&&\displaystyle \frac{{dP}}{{{dM}}_{{\rm{UV}}}}({M}_{\star },z)={ \mathcal N }[{M}_{{\rm{UV}}},{\bar{M}}_{{\rm{UV}}}({M}_{\star },z),{\sigma }_{{uv}}],\end{eqnarray} \tag{ 1 }$$

where the mean relationship between UV luminosity at a given stellar mass and redshift is

$$\begin{eqnarray}&&{\bar{M}}_{\mathrm{UV}}({M}_{\star },z)\equiv \int {M}_{{\rm{UV}}}\displaystyle \frac{{dP}}{{{dM}}_{{\rm{UV}}}}({M}_{\star },z){{dM}}_{{\rm{UV}}}.\end{eqnarray} \tag{ 2 }$$

Once such a relation for ${M}_{{\rm{UV}}}$–${M}_{\star }$ and its scatter has been adopted (see Section 3.2), the observed UV luminosity function ${\rm{\Phi }}({M}_{{\rm{UV}}},z)$ can be modeled as the convolution of the stellar mass function ${\rm{\Phi }}({M}_{\star },z)$ with the distribution of ${M}_{{\rm{UV}}}$

$$\begin{eqnarray}&&{\rm{\Phi }}({M}_{{\rm{UV}}},z)\\ &&\qquad =\,{\displaystyle \int }_{0}^{\infty }{\rm{\Phi }}({M}_{\star },z)\,{ \mathcal N }[{M}_{{\rm{UV}}},{\bar{M}}_{{\rm{UV}}}({M}_{\star },z),{\sigma }_{{uv}}]\,{{dM}}_{\star },\end{eqnarray} \tag{ 3 }$$

where ${\rm{\Phi }}({M}_{{\rm{UV}}},z)$ represents the number of galaxies per co-moving volume with UV absolute magnitude ${M}_{{\rm{UV}}}$ as a function of redshift. Therefore to calculate star-forming galaxy counts at z ≳ 4, where we have the best constraints from the UV luminosity function, we forward model the continuously evolving stellar mass function, convolved with an empirical characterization of ${ \mathcal N }[{M}_{{\rm{UV}}},{\bar{M}}_{{\rm{UV}}}({M}_{\star },z),{\sigma }_{{uv}}]$ in order to fit with observed UV luminosity functions over the range $4\lesssim z\lesssim 10$.¹⁸

This procedure enables us to produce one continuously evolving stellar mass function that, when sampled randomly as outlined above, produces star-forming galaxy counts that follow observed stellar mass functions at z ≲ 4, the ${M}_{{\rm{UV}}}$–${M}_{\star }$ distribution given by ${ \mathcal N }[{M}_{{\rm{UV}}},{\bar{M}}_{{\rm{UV}}}({M}_{\star },z),{\sigma }_{{uv}}]$, and the observed UV luminosity functions at z ≳ 4. We will describe the characterization of the empirical ${M}_{{\rm{UV}}}$–${M}_{\star }$ distribution, ${ \mathcal N }[{M}_{{\rm{UV}}},{\bar{M}}_{{\rm{UV}}}({M}_{\star },z),{\sigma }_{{uv}}]$, that we use to forward model the stellar mass function in Section 3.2. In Section 3.1 we will describe our procedure to fit the observed evolving stellar mass function over $0.2\lt z\lesssim 4$ and forward model by convolving the mass function with ${ \mathcal N }[{M}_{{\rm{UV}}},{\bar{M}}_{{\rm{UV}}}({M}_{\star },z),{\sigma }_{{uv}}]$ at z ≳ 4, to produce galaxy counts.

For each object generated in the mock galaxy population, we use redshift and stellar mass to assign other integrated galaxy properties including UV absolute magnitude ${M}_{{\rm{UV}}}$ and continuum slope β for star-forming galaxies, as well as type-dependent structural parameters. To assign the integrated properties we use empirical relations, plus appropriate scatter, to generate smoothly redshift-evolving distributions of ${M}_{{\rm{UV}}}$–${M}_{\star }$ (Section 3.2), β–${M}_{{\rm{UV}}}$ (Section 3.3), size–mass (z < 4), and size-UV luminosity (at z > 4; see Section 5). Figures 1 and 2 provide more details on this procedure. The integrated properties inform the assignment of a fully consistent SED to the mock galaxies, from which we derive JWST, HST, and Spitzer filter photometry. These SEDs are created using beagle and span a range of physical properties, as described in the following section.

beagle (Chevallard & Charlot 2016; C16) is a new-generation tool for the modeling and interpretation of spectro-photometric galaxy SEDs based on a self-consistent approach to describe stellar emission and its transfer through the interstellar (ISM) and intergalactic (IGM) media. In this section, we first describe the general characteristics of beagle and the models integrated therein, and then summarize our two methods for assigning SEDs to mock galaxies according to whether or not the realized properties overlap with those of observed galaxies from current surveys.

In beagle, the emission from simple stellar populations of different ages, t', and metallicities, $Z$ (the mass fraction of all elements heavier than Helium), is described by the latest version of the Bruzual & Charlot (2003) population synthesis code. Stellar emission is computed using the MILES stellar library (Sánchez-Blázquez et al. 2006) and includes new prescriptions for the evolution of massive stars (Bressan et al. 2012; Chen et al. 2015) and their spectra (Hamann & Gräfener 2004; Leitherer et al. 2010). We account for the (continuum+line) emission of gas photoionized by young stars by considering the large grid of photoionization models of Gutkin et al. (2016). These are based on the standard photoionization code cloudy (version 13.3; Ferland et al. 2013) and assume "ionization bounded" nebulae (i.e., a zero escape fraction of H-ionizing photons). The models are described in terms of "effective" (i.e., galaxy-wide) parameters following the prescription of Charlot & Longhetti (2001). Adjustable model parameters include the ionization parameter $\mathrm{log}{U}_{{\rm{S}}}$, which sets the ratio of H-ionizing photons to H atoms at the edge of the Strömgren sphere, the interstellar metallicity ${Z}_{{\rm{ISM}}}$, and the dust-to-metal (mass) ratio ${\xi }_{{\rm{d}}}$, which traces metal depletion onto dust grains. Since the gas density ${n}_{{\rm{H}}}$ and depletion factor ${\xi }_{{\rm{d}}}$ do not significantly affect emission line ratios at sub-solar metallicities (see Figures 3 and 5 of Gutkin et al. 2016), and most of our galaxies exhibit $\mathrm{log}(Z/{Z}_{\odot })\,\lesssim -0.5$ (see Figure 12), we fix ${n}_{{\rm{H}}}={10}^{2}\,{\mathrm{cm}}^{-3}$, the typical value measured in ${\text{}}z\sim 2\mbox{--}3$ galaxies (e.g., Sanders et al. 2016; Strom et al. 2017), and ${\xi }_{{\rm{d}}}=0.3$, a value similar to that measured in the Solar neighborhood (although, see Section 6.5). We account for attenuation by dust of the emission from stars and photoionized gas using the two-component model of Charlot & Fall (2000), parametrized in terms of the total attenuation optical depth ${\hat{\tau }}_{{\rm{V}}}$, and the fraction of this arising in the diffuse ISM μ. The mean effects of intergalactic medium absorption are included following the model of Inoue et al. (2014).

For mock galaxies with properties that are observable using current facilities, we use beagle to generate a distribution of model SEDs consistent with the observations and assign these SEDs to the mock objects. To achieve this, we fit SED models from beagle to the multi-band photometry of galaxies in two CANDELS fields using the 3D-HST catalog (Skelton et al. 2014). When performing parameter estimation, beagle employs the nested sampling algorithm (Skilling 2006) as implemented in MultiNest (Feroz et al. 2009). This procedure creates a range of statistically acceptable SED fits for each observed galaxy in a subset of the 3D-HST sources (see Sections 3.4.2 and 4.2, while for more detail of the beagle output, see C16, Section 3.3), which are then used to produce a parent catalog. This parent catalog is used to assign SEDs to mock objects with high stellar mass (i.e., those with mass above $\mathrm{log}(M/{M}_{\odot })\gt 8$, or above the mass completeness of the 3D-HST catalog if larger in that redshift bin) and low redshift (z < 4), where the $\lambda \lesssim 4.5\,\mu {\rm{m}}$ photometry provides firm constraints on stellar mass. The SEDs are assigned by finding the closest match in stellar mass and redshift for each mock galaxy within the parent catalog, allowing us to encapsulate the observed diversity of galaxy SEDs at z < 4 with relatively few assumptions.

For mock galaxies with realized properties that extend beyond current measurements of real sources, we can leverage the capabilities of beagle to produce theoretical SEDs and generate model spectra for the mock objects. In this second method, we generate a parent catalog built of theoretical SEDs covering a range of model parameters that can be matched to mock galaxy stellar mass, redshift, and, for star-forming galaxies, ${M}_{{\rm{UV}}}$, and β (see Sections 3.4.3 and 4.2). We use this method at low stellar masses $[\mathrm{log}(M/{M}_{\odot })\lt 8]$, where current galaxy survey sampling of the population is less complete, and at $z\geqslant 4$, where SED coverage in the rest-frame optical is only available from imaging taken with IRAC, the 3.6–8 μm camera on Spitzer (Fazio et al. 2004).

In generating a mock galaxy catalog, we aim to reproduce measurements of the star-forming galaxy stellar mass functions at low redshift (z ≲ 4) and the UV luminosity function at high redshift (z ≳ 4). Our primary mass function constraints come from Tomczak et al. (2014, hereafter T14), while our UV luminosity function constraints are adopted from Bouwens et al. (2015) at $4\lesssim z\lesssim 8$ and the newest z ∼ 10 estimate presented in Oesch et al. (2017).

T14 provide measurements of the stellar mass function of star-forming and quiescent galaxies in eight redshift bins in the range $0.2\lt z\lt 3$. They employed imaging data from the FourStar galaxy evolution (ZFOURGE) survey (Straatman et al. 2016) covering the CDFS, COSMOS, and UDS fields with five near-IR medium-bandwidth filters spanning the J and H bands, as well as broad-band K_S imaging. Specifically they used the regions that also overlap with CANDELS J₁₂₅ and H₁₆₀ imaging (to ∼26.5 depth to 5σ), covering a total area of ∼316 arcmin². Additionally, imaging from NEWFIRM Medium-band Survey (Whitaker et al. 2011) was used in the AEGIS and COSMOS fields, employing the same filter sets as the ZFOURGE survey to shallower depths but wider area to leverage better constraints of the high-mass end of the mass function. Each of the fields also benefit from further imaging that allows comprehensive sampling of galaxy SEDs over the wavelength range 0.3–8 μm, with the field-specific filter sets and imaging programs summarized in Section 2.4 of Straatman et al. (2016).

T14 inferred photometric redshifts and rest-frame colors (used to separate galaxies into star-forming or quiescent based on the UVJ diagram of Whitaker et al. 2011) using the template-based eazy code (Brammer et al. 2008), while stellar masses were estimated using fast (Kriek et al. 2009). Within fast, they used the original Bruzual & Charlot (2003) population synthesis code at fixed solar metallicity, employing a Chabrier (2003) IMF, and a declining exponential star formation history. The 80% mass completeness limits of their sample increase from $\mathrm{log}(M/{M}_{\odot })\sim 7.75$ at z ∼ 0.5 to $\mathrm{log}(M/{M}_{\odot })\sim 9.25$ at z ∼ 3. T14 fit their resulting stellar mass functions with a sum of two Schechter (1976) functions:

$$\begin{eqnarray}{\rm{\Phi }}({M}_{\star }){\text{}}{{dM}}_{\star } & = & \,{{\rm{\Phi }}}_{1}({M}_{\star }){\text{}}{{dM}}_{\star }+{{\rm{\Phi }}}_{2}({M}_{\star })\,{\text{}}{{dM}}_{\star }\\ & = & \,\mathrm{ln}10{\phi }_{1,{\rm{M}}}^{* }\,{10}^{({M}_{\star }-{M}_{1,{\rm{M}}}^{* })(1+{\alpha }_{1,{\rm{M}}})}\\ & & \times {\rm{\exp }}(-{10}^{{M}_{\star }-{M}_{1,{\rm{M}}}^{* }}){\text{}}{{dM}}_{\star }\\ & & +\mathrm{ln}10{\phi }_{2,{\rm{M}}}^{* }\,{10}^{({M}_{\star }-{M}_{2,{\rm{M}}}^{* })(1+{\alpha }_{2,{\rm{M}}})}\\ & & \times {\rm{\exp }}(-{10}^{{M}_{\star }-{M}_{2,{\rm{M}}}^{* }}){\text{}}{{dM}}_{\star },\end{eqnarray} \tag{ 4 }$$

where ${M}_{\star }=\mathrm{log}(M/{M}_{\odot })$, as defined in Section 2.1, ${\rm{\Phi }}({M}_{\star })$ indicates the number of galaxies per Mpc³ with stellar masses between ${M}_{\star }$ and ${M}_{\star }$ + d${M}_{\star }$, and ${M}_{1,{\rm{M}}}^{* }$, ${M}_{2,{\rm{M}}}^{* }$, ${\phi }_{1,{\rm{M}}}^{* }$, ${\phi }_{2,{\rm{M}}}^{* }$, ${\alpha }_{1,{\rm{M}}}$, and ${\alpha }_{2,{\rm{M}}}$ are the six free parameters of the function.¹⁹ In a single Schechter function, ${M}_{{\rm{M}}}^{* }$ is the mass at the turnover, or "knee" of the mass function; ${\phi }_{{\rm{M}}}^{* }$ is the characteristic number density of galaxies at the turnover; and ${\alpha }_{{\rm{M}}}$ is the low-mass slope. In the double-Schechter function used in T14, they explicitly set ${M}_{1,{\rm{M}}}^{* }={M}_{2,{\rm{M}}}^{* }={M}_{{\rm{M}}}^{* }$, meaning that they fit with a single "knee" but the different normalizations and faint-end slopes of each function enable them to fit the observed steepening of the mass function to low masses (see Figure 4).

At z > 4 stellar masses become progressively less well constrained from measurements, in part because the rest-frame optical SED (a key region containing the Balmer break at ∼3600 Å, and the 4000 Å break) shifts into the infrared where current facilities have low sensitivity. Additionally, high equivalent width (EW) emission lines can add to the flux in the reddest photometric bands, leading to an overprediction of galaxy stellar masses (Schaerer & de Barros 2010; Curtis-Lake et al. 2013; Stark et al. 2013; de Barros et al. 2014). As a result, relative uncertainties on stellar mass measurements are high (e.g., 0.4 dex at ${10}^{10}\,{M}_{\odot }$ at z = 4, increasing with redshift and decreasing mass; Grazian et al. 2015; see also Mobasher et al. 2015) and may contribute to the large scatter of mass function measurements in the literature (nearly ∼1 dex in counts; see Figure 9 in Song et al. 2016, and Figure 11 in Davidzon et al. 2017). Therefore, to generate galaxy counts at z > 4, we leverage the constraints provided by the observed UV luminosity function in the range of $4\lesssim z\lesssim 8$ from Bouwens et al. (2015) with luminosity function measurements with mean redshifts at $\langle z\rangle =[3.8,4.9,5.9,6.8,7.9]$ using data from the HST Legacy Fields, as well as the z ∼ 10 luminosity function of Oesch et al. (2017). The binned UV luminosity function measurements we use for this work are overall consistent with many other results in the literature at ${M}_{{\rm{UV}}}$ < −17 (e.g., McLure et al. 2013; Atek et al. 2015; Finkelstein et al. 2015; Laporte et al. 2015; Castellano et al. 2016; Bouwens et al. 2017b; Livermore et al. 2017; Ono et al. 2018; Yue et al. 2017).

We choose to model the redshift evolution of the six mass function parameters across the entire redshift range of the mock (i.e., $0.2\lt z\lt 15$). This ensures a smooth evolution in number counts across the transition from mass to luminosity function-based constraints. At z < 3.8 (the mean redshift of the B-dropout sample used to produce the Bouwens et al. (2015) z ∼ 4 luminosity function), we use the measured mass functions of T14 to directly constrain their redshift evolution, while at $z\geqslant 3.8$ we use our model of the redshift-evolving ${M}_{{\rm{UV}}}$–${M}_{\star }$ relation (see Section 3.2) to fit to the observed luminosity functions with mass function parameters. However, it is important to note that this is not a direct prediction of the shape or evolution of the z ≳ 4 mass functions that we expect to measure with JWST. Our z ≳ 4 mass functions are dependent on our model of the ${M}_{{\rm{UV}}}$–${M}_{\star }$ relation, and additionally we do not yet know how incomplete the current ${M}_{{\rm{UV}}}$-selected samples at z ≳ 4 may be.

To determine a suitable form for the redshift evolution of the Schechter function parameters, we first need to know what mass function parameters can reproduce the observed UV luminosity functions at z ≳ 4. The details of this fitting are given in Appendix A.1, and we plot the Schechter (mass) function parameters that best fit the z ≳ 4 luminosity function observations in Figure 3, as well as the individual maximum-likelihood estimates of T14 at z < 4. The estimates of ${M}_{{\rm{M}}}^{* }$ derived from the measured luminosity functions at z ≳ 4 are significantly lower than the T14 measurements. However, if we fit the evolution of a double-Schechter function with different "knees," as in Equation (4), we can use ${M}_{1,{\rm{M}}}^{* }$ to fit to the high-mass end of the z < 3 mass functions, while ${M}_{2,{\rm{M}}}^{* }$ (plus the fast evolution of ${\phi }_{1,{\rm{M}}}^{* }$) can be used to account for the rapid evolution in the bright end required to fit to the z ≳ 4 luminosity functions. We therefore choose to set the $z\geqslant 3.8$ evolution of ${M}_{2,{\rm{M}}}^{* }$, ${\alpha }_{2,{\rm{M}}}$, and ${\phi }_{2,{\rm{M}}}^{* }$ using a weighted least-squares linear regression to the luminosity function fits. We extrapolate the linear fit of ${M}_{2,{\rm{M}}}^{* }$ to z < 3.8 but re-fit the T14 measured mass functions, allowing the other five Schechter function parameters to vary. In fact, this choice of ${M}_{2,{\rm{M}}}^{* }$ evolution somewhat under-estimates the high-mass end of the $2\lt z\lt 3$ mass functions (see Figure 4). It is entirely possible that the reason for the strong evolution in ${M}_{{\rm{M}}}^{* }$ seen between the z < 3.8 and $z\geqslant 3.8$ samples is due to the ${M}_{{\rm{UV}}}$-selected samples missing a population of dusty, high-mass star-forming galaxies. If they exist, these objects will be revealed by JWST, but currently we lack firm constraints on their number density evolution. We are basing this mock catalog on current observational constraints, and so choose to favor the fit to the z ∼ 4 luminosity function over the $2\lesssim z\lesssim 3$ mass function at the high-mass end, as it allows us to produce a model with number counts that vary relatively smoothly with redshift. As such, a caveat of our model is that we are not modeling the dusty star-forming galaxies currently missed in UV-selected samples, and mildly under-represent the high-mass end of the $2\lesssim z\lesssim 3$ mass functions. A model that simultaneously fits the z ∼ 2.75 T14 data and the z ∼ 4 Bouwens et al. (2015) luminosity function would require a strong gradient discontinuity in ${M}_{1,{\rm{M}}}^{* }$ that would lead to a step discontinuity in the number counts of galaxies at high stellar masses.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When fixing the evolution of the Schechter function parameters at z ≳ 4, we use a weighted least-squares linear regression excluding the point at z ∼ 10, which has noticeably lower number densities than can be accounted for by a simple linear relation in all three parameters. In fact, the exact form of the redshift evolution of the UV luminosity function, and associated CSFRD above z ∼ 8, has been an area of active debate in the literature—for example, with McLeod et al. (2016) presenting measurements of the z ∼ 9–10 luminosity function that are consistent with a smooth decline in the CSFRD. For our fiducial mock catalog we choose to base the model on the Oesch et al. (2017) results in order to provide a conservative limit on z ≳ 8 galaxy number counts likely to be detected with JWST. We defer further discussion of this issue to Section 7. The constraints at z ∼ 10 are not strong enough to constrain the likely evolution in ${M}_{2,{\rm{M}}}^{* }$, ${\phi }_{2,{\rm{M}}}^{* }$ and ${\alpha }_{2,{\rm{M}}}$. We thus choose to re-fit the z ∼ 10 luminosity function with ${\alpha }_{2,{\rm{M}}}$ and ${M}_{2,{\rm{M}}}^{* }$ fixed to the values defined by the extrapolated linear fits at z = 10, giving $\mathrm{log}({\phi }_{2,{\rm{M}}}^{* }/{{\rm{Mpc}}}^{3}\,{{\rm{dex}}}^{-1})=-4.67\pm 0.3$. We then require the gradient of the evolution in $\mathrm{log}({\phi }_{2,{\rm{M}}}^{* })$ to decrease further at z > 8, so that this value is reached by the relation at z = 10.

At z < 3.8 we re-fit to the T14 mass functions using a Bayesian multi-level modeling approach (see Appendix A.2), which allows us to derive the best-fit redshift evolution of the Schechter function parameters by fitting to the mass function measurements in each redshift bin simultaneously. This approach is more powerful than fitting a functional form to the published Schechter parameter estimates, as it accounts for parameter covariance self-consistently. At z < 3.8, we choose a functional form for the redshift evolution for ${\alpha }_{2,{\rm{M}}}$ and ${\phi }_{2,{\rm{M}}}^{* }$ that asymptotically approaches the value of the best-fit linear relation at z = 3.8, but decreases rapidly at the lowest redshifts. Without this dip to low redshifts, the mass function is too shallow, with too high a normalization at low masses. We accept a mildly discontinuous evolution at z = 3.8 because allowing the functional form in either ${\alpha }_{2,{\rm{M}}}$ or ${\phi }_{2,{\rm{M}}}^{* }$ to increase and turn over by z = 3.8 (to give smooth evolution at z = 3.8) produces mass functions that cross over at low masses, a situation that we are trying to avoid by requiring that our model is monotonically increasing at given mass with decreasing redshift.

The redshift evolution of each Schechter function parameter is summarized below:

$$\begin{eqnarray}&&{M}_{1,{\rm{M}}}^{* }(z)={a}_{1}\,\,\,\,\,\,\,\,\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\mathrm{log}[{\phi }_{1,{\rm{M}}}^{* }(z)]={b}_{1}+{b}_{2}\,z+{b}_{3}\,{z}^{2}\,\,\,\,\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\alpha }_{1,{\rm{M}}}(z)={c}_{1}+{c}_{2}\,z\,\,\,\,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{M}_{2,{\rm{M}}}^{* }(z)={D}_{1}+{D}_{2}\,z\,\,\,\,\,\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}\mathrm{log}[{\phi }_{2,{\rm{M}}}^{* }(z)] & = & {e}_{1}\,[1-\exp (-z)]+{e}_{2} & \,z\lt 3.8\\ & = & {E}_{1}+{E}_{2}z & 3.8\leqslant z\lt 8\\ & = & {E}_{1}^{{\prime} }+{E}_{2}^{{\prime} }z & z\geqslant 8\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}\,{\alpha }_{2,{\rm{M}}}(z) & = & {f}_{1}\,[1-\exp (-z)]+{f}_{2} & & & z\lt 3.8\\ & = & {F}_{1}+{F}_{2} & & & \,zz\geqslant 3.8,\end{eqnarray} \tag{ 10 }$$

where the parameters D₁, D₂, E₁, E₂, F₁, and F₂ are all determined from the linear regression to the forward-modeled luminosity function fitting, and ${E}_{1}^{{\prime} }$ and ${E}_{2}^{{\prime} }$ are chosen to fit to the z ∼ 10 luminosity function while maintaining continuous evolution in log[${\phi }_{2,{\rm{M}}}^{* }(z)]$ at z = 8. The parameters e₂ and f₂ are fixed to the values required to produce continuous evolution at z = 3.8 with ${e}_{2}={E}_{1}+3.8\,{E}_{2}-{e}_{1}[1-\exp (-3.8)]$ and ${f}_{2}\,={F}_{1}+3.8\,{F}_{2}-{f}_{1}[1-\exp (-3.8)]$. The remaining free parameters a₁, b₁, b₂, b₃, c₁, c₂, e₁, and f₁, are then constrained using multi-level modeling (see Appendix A.2) to the published T14 star-forming mass functions (their Table 1).

Table 1.  The Values of the Parameters Used in Our Model of the Mass Function Evolution, as Described in Equations (5)–(10)

Parameter	Median	1σ Uncertainty	Prior/Source of Fits
a1	10.69	0.04	${ \mathcal N }(0,50)$
b1	−2.68	0.16	${ \mathcal N }(0,50)$
b2	0.06	0.24	${ \mathcal N }(0,50)$
b3	−0.19	0.08	${ \mathcal N }(0,50)$, $\in [-\infty ,0]$
c1	−1.02	0.16	${ \mathcal N }(0,50)$
c2	0.29	0.13	${ \mathcal N }(0,50)$
D1	10.30	0.10	Linear fitting 4 ≲ z < 8
D2	−0.15	0.02	Linear fitting 4 ≲ z < 8
e1	0.73	0.26	${ \mathcal N }(0,50)$, $\in [0,\infty ]$
e2	−3.60	⋯	$={E}_{1}+3.8{E}_{2}-{e}_{1}[1-\exp (-3.8)]$
E1	−2.03	0.41	Linear fitting 4 ≲ z < 8
E2	−0.23	0.09	Linear fitting 4 ≲ z < 8
${E}_{1}^{{\prime} }$	−0.67	⋯	${\phi }_{2,{\rm{M}}}^{* }$ fit to z ∼ 10 LF
${E}_{2}^{{\prime} }$	−0.40	⋯	${\phi }_{2,{\rm{M}}}^{* }$ fit to z ∼ 10 LF
f1	0.41	0.17	${ \mathcal N }(0,50)$, $\in [0,\infty ]$
f2	−1.82	⋯	$={F}_{1}+3.8{F}_{2}-{f}_{1}[1-\exp (-3.8)]$
F1	−1.16	0.10	Linear fitting 4 ≲ z < 8
F2	−0.07	0.02	Linear fitting 4 ≲ z < 8

Note. For those parameters determined using the multi-level model fitting to z < 4 mass functions, we report the median of the posterior distribution function, its 1σ confidence interval, as well as the prior used in the fitting.

Download table as:  ASCIITypeset image

We report the median values and associated uncertainties, along with the values for the model parameters defined by linear fits to the z ≳ 4 individual mass function estimates in Table 4. The chosen redshift evolution of each parameter is plotted as the orange lines in Figure 3. The resulting mass function comparisons to the T14 measurements at z < 4 are plotted in Figure 4, and the luminosity function comparisons are shown in Figure 5.

Zoom In Zoom Out Reset image size

In this section, we describe our method to characterize the relation (slope, intercept, scatter) between ${M}_{{\rm{UV}}}$ and ${M}_{\star }$ of galaxies at redshifts $0.2\leqslant z\leqslant 15$. Hereafter, we use the definition of ${M}_{{\rm{UV}}}$ adopted, for example, in Robertson et al. (2013), as the average magnitude at rest-frame wavelength within a flat filter centered at 1500 Å and with a width of 100 Å, which is the definition adopted by beagle (Chevallard & Charlot 2016). This definition of ${M}_{{\rm{UV}}}$ differs slightly from that used to measure the UV luminosity functions in Bouwens et al. (2015), which define ${M}_{{\rm{UV}}}$ to be at rest-frame wavelength of 1600 Å. We calculate the typical color correction based on the mean β as a function of ${M}_{{\rm{UV}}}$ and redshift presented in Bouwens et al. (2014) and find that the typical difference in magnitudes between 1500 and 1600 Å is negligible ($| \delta {M}_{{\rm{UV}}}| \lesssim 0.05$). This correction is significantly smaller than the k-correction applied to estimate ${M}_{{\rm{UV}}}$ at 1600 Å from broad-band photometry in the first place ($| \delta {M}_{{\rm{UV}}}| \lesssim 0.1$), and so we apply no conversion between rest-frame 1500 and 1600 Å ${M}_{{\rm{UV}}}$ values.

The ${M}_{{\rm{UV}}}$–${M}_{\star }$ distribution and its evolution are critical components of our underlying phenomenological model, and are required to statistically assign UV luminosities to mock galaxies generated from our continuously evolving stellar mass function model. However, we note the following uncertainties to this procedure. At all redshifts, galaxies exhibit a diversity of mass-to-light ratios, which depend on the stellar population properties (age, metallicity), star formation history, and dust content of a galaxy. As a result, the exact form of the relation between ${M}_{{\rm{UV}}}$ and ${M}_{\star }$ and its dependency on galaxy properties are largely unknown. In general, brighter galaxies at UV wavelengths correspond to more massive objects (e.g., Stark et al. 2009; González et al. 2011; Lee et al. 2011), and this holds out to z ∼ 7 (Duncan et al. 2014; Grazian et al. 2015; Salmon et al. 2015; Song et al. 2016). Although the relation between ${M}_{{\rm{UV}}}$ and ${M}_{\star }$ follows a general trend of decreasing ${M}_{{\rm{UV}}}$ with ${M}_{\star }$ out to $\mathrm{log}(M/{M}_{\odot })\sim 10$, at higher masses the average ${M}_{{\rm{UV}}}$ becomes fainter due to the appearance of a population of fainter objects. This trend could be attributable to several effects, such as increased dust content and older average stellar ages among massive galaxies (e.g., Spitler et al. 2014). Characterizing the relationship is further complicated by the difficulty of measuring stellar mass owing to emission line contamination at high redshift (Labbé et al. 2013; Stark et al. 2013) and at low stellar masses (Whitaker et al. 2014), and the lack of a direct photometric probe of ${M}_{{\rm{UV}}}$ at intermediate redshifts ($0.6\lt z\lt 1.5$). The procedure we outline here has a direct impact on the resulting UV luminosity functions (see Sections 2.1, 3.1, and 6.1). We have therefore developed a straightforward description of the ${M}_{{\rm{UV}}}$–${M}_{\star }$ distribution and its evolution that is designed to encapsulate the diversity of real galaxies.

3.2.1. Characterizing the Evolution of ${M}_{{\rm{UV}}}\mbox{--}{M}_{\star }$ from Observations

We characterize the ${M}_{{\rm{UV}}}$–${M}_{\star }$ relationship at z ≲ 4 using measurements from the 3D-HST catalog (using SED fitting with beagle; see description in Section 3.4.1). As discussed extensively in Stefanon et al. (2017a), selection effects can heavily influence the observed shape of the ${M}_{{\rm{UV}}}$–${M}_{\star }$ distribution. Therefore we avoid including observed galaxies whose ${M}_{{\rm{UV}}}$ or ${M}_{\star }$ measurements are poorly constrained by the beagle fits. Specifically we only use galaxies with $\delta \mathrm{log}(M/{M}_{\odot })\lt 1$, $\delta z\lt 1$, and $\delta {M}_{{\rm{UV}}}\lt 1$ (where, e.g., $\delta z$ is the 68% credibility interval on redshift). The limits imposed were chosen to avoid biasing the characterization of the ${M}_{{\rm{UV}}}$–${M}_{\star }$ distribution with overly strict ${M}_{{\rm{UV}}}$ or ${M}_{\star }$ cuts, which we discuss further below.

In Figure 6, we plot the ${M}_{{\rm{UV}}}$–${M}_{\star }$ distributions for the 3D-HST galaxies with well-constrained ${M}_{{\rm{UV}}}$, ${M}_{\star }$, and redshift measurements. As discussed above, these distributions show a trend of increasing stellar mass with decreasing ${M}_{{\rm{UV}}}$ at low stellar mass. At high stellar mass the ${M}_{{\rm{UV}}}$ values tend to be fainter than the linear relation, as observed in Spitler et al. (2014). Rather than attempting to fully model this mass-dependent behavior, especially given the Malmquist biases that begin to affect the higher-redshift bins, we adopt the following two-step procedure to describe the ${M}_{{\rm{UV}}}$–${M}_{\star }$ distributions at z ≲ 4. We fit the observed ${M}_{{\rm{UV}}}$–${M}_{\star }$ distribution under the simplest assumption of a linear relationship to extrapolate to low masses, while at higher masses ($\mathrm{log}(M/{M}_{\odot })$ ≳ 8–8.5, depending on the mass limit at a given redshift), we assign ${M}_{{\rm{UV}}}$ values by sampling from real galaxies of the same mass. The matching procedure allows us to maintain the observed flattening of the distribution at high masses, and is fully described in Section 3.4.2 below.

Zoom In Zoom Out Reset image size

Figure 6 illustrates the substantial scatter in the observed ${M}_{{\rm{UV}}}$–${M}_{\star }$ distributions at z ≲ 4. Owing to the large scatter, the best-fitting slope will depend strongly on the uncertainties on the data points, and the size of the uncertainties may depend on ${M}_{{\rm{UV}}}$, ${M}_{\star }$, and also plausibly on redshift. Indeed, we find that when fitting with both slope and normalization as free parameters, neither parameter is well constrained, and the best-fitting slope is highly variable between redshift bins. Therefore we adopt a fixed slope for the ${M}_{{\rm{UV}}}$–${M}_{\star }$ relation at all redshifts and fit only the intercept at each redshift. This procedure essentially fits the average redshift-dependent mass to light ratio, which has lower uncertainty and is less dependent on the error on individual galaxy measurements and stellar mass-dependent systematics. Several studies have reported the measurement of constant slope for UV-selected galaxies, with normalization evolving in redshift (Duncan et al. 2014; Grazian et al. 2015; Salmon et al. 2015; Song et al. 2016; Stefanon et al. 2017a), and find a reasonable description of the data. The blue dashed line in Figure 6 shows our best-fit relation to the ${M}_{{\rm{UV}}}$–${M}_{\star }$ distribution in each redshift bin, where the slope is fixed to a value of −1.66. We find excellent agreement with the observed distribution at all redshifts. For reference we also indicate the stellar mass limits in each redshift bin, above which we assign ${M}_{{\rm{UV}}}$ values by sampling from real galaxies (red dashed lines). Fitting the ${M}_{{\rm{UV}}}$–${M}_{\star }$ distribution only above these mass limits instead has a negligible effect on the result at z < 3. At z ∼ 3.75, where there are fewer well-constrained measurements, fitting above this mass limit would increase the ${M}_{{\rm{UV}}}$–${M}_{\star }$ intercept by ∼0.1 mag, an indication that fitting only at the high-mass end biases the characterization of the ${M}_{{\rm{UV}}}$–${M}_{\star }$ due to the high-mass end flattening. Therefore we choose to proceed using all galaxies with well-characterized stellar mass, redshift, and ${M}_{{\rm{UV}}}$.

To set the full redshift evolution of the ${M}_{{\rm{UV}}}$–${M}_{\star }$ relation, we combine the intercept values for the best-fit relations with fixed slope at each redshift z ≲ 4, with measurements of the ${M}_{{\rm{UV}}}$–${M}_{\star }$ intercept at z > 4. We use the average observed value of stellar mass for bright (${M}_{{\rm{UV}}}$ = −20) galaxies at $4\lt z\leqslant 7$ to set the overall normalization in each z > 4 redshift bin, while assuming the same constant ${M}_{{\rm{UV}}}$–${M}_{\star }$ slope. We utilize the high-redshift stellar mass measurements shown in Figure 7 of Stark et al. (2013), where the measured stellar masses were fit while including the contribution to the SED from nebular emission lines. The normalization value at ${M}_{{\rm{UV}}}=-20$ shows an overall decline in the range of $4\lesssim z\lt 7$, indicating a decrease in the average mass to light ratios of galaxies with increasing redshift. The measured values for the ${M}_{{\rm{UV}}}$–${M}_{\star }$ intercepts at all redshift bins, evaluated at ${M}_{{\rm{UV}}}=-20$, are shown as points in Figure 7.

Zoom In Zoom Out Reset image size

3.2.2. Continuous Redshift Evolution of ${M}_{{\rm{UV}}}\mbox{--}{M}_{\star }$

Our model of the redshift evolution of the ${M}_{{\rm{UV}}}$–${M}_{\star }$ relation is defined by fitting to the evolving intercept values shown in Figure 7. These intercepts show a rapid decline from z ∼ 0–3, with a shallower decline at z > 4. We find that at $z\leqslant 4$ the intercept measurements are adequately described by a quadratic function, and a linear function at z > 4. To ensure these two functions remain continuous at z ∼ 4, we include the boundary condition that the derivatives of the two functions are equal at the center of the highest redshift bin where we use the 3D-HST data (z = 3.75). This constraint results in the following function to describe the intercept of the ${M}_{{\rm{UV}}}$–${M}_{\star }$ relation, ${M}_{{M}_{{\rm{UV}}}\mbox{--}{M}_{\star }}^{0}(z)$, evaluated at ${M}_{{\rm{UV}}}=-20$ and $z\leqslant 3.75$:

$$\begin{eqnarray}{M}_{{M}_{{\rm{UV}}}\mbox{--}{M}_{\star }}^{0}(z) & = & a\,{(z-3.75)}^{2}\\ & & +[b-2\,a(z-3.75)]\,(z-3.75)+c,\end{eqnarray} \tag{ 11 }$$

where a = 0.12, b = 0.08, and c = 9.41. This evolution of the ${M}_{{\rm{UV}}}$–${M}_{\star }$ intercept is shown as the black solid curve in Figure 7, along with the observed data (blue points). The resulting linear ${M}_{{\rm{UV}}}$–${M}_{\star }$ relations in each redshift bin, according to the smoothly evolving intercept function defined in Equation (11), are also shown as black solid lines in Figure 6. At z > 3.75, the ${M}_{{\rm{UV}}}$–${M}_{\star }$ intercept evaluated at ${M}_{{\rm{UV}}}=-20$ evolves approximately linearly with the following form:

$$\begin{eqnarray}{M}_{{M}_{{\rm{UV}}}\mbox{--}{M}_{\star }}^{0}(3.75\lt z\lt 8) & = & -0.12\times z+9.88\\ {M}_{{M}_{{\rm{UV}}}\mbox{--}{M}_{\star }}^{0}(z\geqslant 8) & = & 8.92.\end{eqnarray} \tag{ 12 }$$

We use Equation (11) to assign ${M}_{{\rm{UV}}}$ values to galaxies of a given stellar mass and redshift at $z\leqslant 3.75$, and Equation (12) to assign ${M}_{{\rm{UV}}}$ values at z > 3.75, assigned randomly within the scatter observed in 3D-HST data in Figure 6. We have characterized this scatter in both stellar mass, ${M}_{{\rm{UV}}}$, and redshift, and find that the scatter in ${M}_{{\rm{UV}}}$ is remarkably constant in both stellar mass and redshift, with an average value of ${\sigma }_{{uv}}\sim 0.7$ magnitudes. Therefore, at all redshifts, we assign ${M}_{{\rm{UV}}}$ values randomly according to this relation and that assumed Gaussian scatter ${\sigma }_{{uv}}$. This direct assignment of ${M}_{{\rm{UV}}}$ applies only to low-redshift low-mass star-forming mock galaxies ($z\leqslant 4$ and have $\mathrm{log}(M/{M}_{\odot })\leqslant 8$), or high-redshift star-forming galaxies at z > 4. Massive, low-redshift mock galaxies (with $\mathrm{log}(M/{M}_{\odot })\gt 8$ and $z\leqslant 4$) are assigned ${M}_{{\rm{UV}}}$ values according to their matched 3D-HST realization. The ${M}_{{\rm{UV}}}$–${M}_{\star }$ relation and scatter, as described here, are used as priors when drawing realizations from fits to 3D-HST galaxies to ensure a smooth transition in the mock catalog ${M}_{{\rm{UV}}}$–${M}_{\star }$ relation at $\mathrm{log}(M/{M}_{\odot })=8$. This procedure is detailed in Section 3.4.2.

3.2.3. Theoretical Limits on the Mass to Light Ratios of Galaxies

As mass to light ratios continue to decrease with increasing redshift and decreasing stellar mass according to our model, the mass to light ratios approach a theoretical limit of the stellar population models we generate with beagle. This limit is set by the UV luminosities of individual massive stars, and represents the minimum mass to light ratio possible for an instantaneous burst of star formation for any given IMF (with no dust attenuation, the lowest metallicity, and corresponding nebular continuum emission). Any IMF choice will result in such a limit in the possible ${M}_{{\rm{UV}}}$ given a stellar mass, with more top-heavy or bottom-light IMFs allowing for brighter limiting ${M}_{{\rm{UV}}}$ and bottom-heavy IMFs producing fainter limiting ${M}_{{\rm{UV}}}$. The Chabrier IMF that we use in this work is relatively bottom-light and has a larger mass to light ratio parameter space than a more bottom heavy IMF (e.g., Salpeter 1955; Kroupa 2001). For a Chabrier IMF with our assumed high-mass cutoff of 100 ${M}_{\odot }$, we find that the stellar plus nebular continuum emission results in a theoretical mass to light ratio given by ${M}_{{\rm{UV}}}\approx -2.45\,\mathrm{log}(M/{M}_{\odot })-1.3$.

To accommodate this theoretical limit in the ${M}_{{\rm{UV}}}$–${M}_{\star }$ evolution of our phenomenological model, we truncate the Gaussian distribution that we use to assign ${M}_{{\rm{UV}}}$ values. For galaxies at the detection limit of future blank surveys (e.g., apparent magnitude ${m}_{\mathrm{app}}\sim 31$) this truncation has a negligible effect on the overall ${M}_{{\rm{UV}}}$–${M}_{\star }$ distribution at z < 4. At z > 4, we find that scatter using the truncated Gaussian at the limit changes the overall shape of the ${M}_{{\rm{UV}}}$–${M}_{\star }$ distribution by steepening the low-mass end. The effect becomes significant by z ∼ 8. We therefore additionally halt the redshift evolution of the mean ${M}_{{\rm{UV}}}$–${M}_{\star }$ relation parametrized in Equation (12) at z = 8. A demonstration of the beagle mass to light ratio limit is shown in Figure 8, compared with the projected evolution of the mean ${M}_{{\rm{UV}}}$–${M}_{\star }$ relation allowed to evolve past z ∼ 8.

Zoom In Zoom Out Reset image size

Although we make every effort to choose reasonable constraints where available, the overall shape of the ${M}_{{\rm{UV}}}$–${M}_{\star }$ distribution at z > 4 is still an extrapolation that impacts various evolutionary relations at high redshift in the model, including the UV luminosity function and sSFR. We will discuss these issues in depth in Section 6.

The spectral slope of the UV continuum (β where ${f}_{\lambda }\propto {\lambda }^{\beta }$) of galaxies is sensitive to the properties of stellar populations (e.g., metallicity and age), star formation history, and dust attenuation. Population studies of star-forming galaxies indicate well-characterized relationships between β and UV luminosity. Bright, massive galaxies tend to have red (i.e., shallower) UV continua, which likely owes to a combination of old stars, higher stellar metallicities, and a larger dust content. The bluer (i.e., steeper) UV continua of lower luminosity galaxies are often associated with younger, less metal-rich stellar populations, and less dust attenuation (e.g., Stanway et al. 2005; Labbé et al. 2007, 2010; Bouwens et al. 2009, 2012b; Rogers et al. 2013, 2014). The detailed relations between β and ${M}_{{\rm{UV}}}$, scatter, and evolution out to redshift z ∼ 8 are still areas of active research, but many studies are consistent with a linear relationship between β and ${M}_{{\rm{UV}}}$ (e.g., Bouwens et al. 2012b, 2014; Alavi et al. 2014; Kurczynski et al. 2014; Rogers et al. 2014) with an average evolutionary trend toward bluer β with increasing redshift (e.g., Labbé et al. 2007; Bouwens et al. 2009; Wilkins et al. 2011; Castellano et al. 2012; Finkelstein et al. 2012b; Bouwens et al. 2014).

3.3.1. Mean $\beta \mbox{--}{M}_{{\rm{UV}}}$ Relation across Cosmic Time

We use a compilation of measured β–${M}_{{\rm{UV}}}$ relations and their scatter across redshifts to assign the rest-frame UV SEDs of mock galaxies. Following several studies (see the previous section), we model the average β–${M}_{{\rm{UV}}}$ relation with a linear function, where the slope $d\beta (z)/{{dM}}_{{\rm{UV}}}$ and intercept $\beta ({M}_{{\rm{UV}}}=-19.5,z)$ of the function vary with redshift. We consider sets of β–${M}_{{\rm{UV}}}$ relations at $1\leqslant z\leqslant 8$ obtained from HST/ACS (Alavi et al. 2014; Kurczynski et al. 2014; Mehta et al. 2017) and HST/WFC3 imaging (Bouwens et al. 2009, 2014). The relationships describing β–${M}_{{\rm{UV}}}$ at the high redshifts we model are broadly consistent with fits measured in other studies (e.g., Bouwens et al. 2012b; Rogers et al. 2014; Finkelstein et al. 2015).

We find that the slope of the β–${M}_{{\rm{UV}}}$ relation shows little evolution at redshifts $1\leqslant z\leqslant 8$, as already found in previous works (e.g., Bouwens et al. 2012b, 2014; Kurczynski et al. 2014). The intercept of the relation increases significantly from redshift z ∼ 1–8, reflecting the evolutionary trend that galaxies have older ages and higher metallicities at later cosmic times (e.g., Labbé et al. 2007). We perform least-squares linear fits to measurements of both $d\beta (z)/{{dM}}_{{\rm{UV}}}$ and $\beta ({M}_{{\rm{UV}}}=-19.5,z)$ and their errors from the literature to produce a mean relation that smoothly evolves with redshift (see Figure 9), described by

$$\begin{eqnarray}d\beta (z)/{{dM}}_{{\rm{UV}}} & = & -0.007\,z-0.09\\ \beta ({M}_{{\rm{UV}}}=-19.5,z) & = & -0.09\,z-1.49.\end{eqnarray} \tag{ 13 }$$

For mock galaxies at $z\leqslant 1$, we extrapolate this relationship to lower redshifts to assign β values. At z > 8 we use the β–${M}_{{\rm{UV}}}$ relationship at z = 8 to assign β to mock galaxies. There exist several motivations to curb the evolution of β–${M}_{{\rm{UV}}}$, with the foremost being the existence of a theoretical limit on the steepness of the UV spectrum emitted by non-Pop. III stars (see discussion in Section 3.3.2 below). Further, the data do not yet constrain evolutionary trends at the highest redshifts currently accessible (${\text{}}z\sim 5\mbox{--}8$). While evolutionary trends with redshift are observed in most analyses (Finkelstein et al. 2012b; Bouwens et al. 2014), the combination of high-redshift color selections with flux boosting from noise in the filters used to measure β are still likely causing statistical studies to be biased against redder β measurements at z > 5 (Dunlop et al. 2012; Rogers et al. 2013). At the very highest redshifts currently accessible (${\text{}}z\sim 7\mbox{--}10$), the data do not provide strong evidence for or against any evolutionary trend in β with redshift or ${M}_{{\rm{UV}}}$ (Dunlop et al. 2013; Wilkins et al. 2016), although the dynamic range in ${M}_{{\rm{UV}}}$ is relatively small at such early times. While evolution cannot be excluded by current data, deep imaging surveys with JWST will enable more robust characterization of the evolution beyond z ∼ 8.

Zoom In Zoom Out Reset image size

In our model, all mock galaxies are assigned a β according to the mean relationship and intrinsic, photometric bias-corrected scatter described in Section 3.3.3. We now detail the theoretical predictions for the UV continuum slope of our combined stellar population and photoionization model.

3.3.2. Theoretical Limits on UV Continuum Slopes

Model predictions on the shape of the UV continuum emission of galaxies depend on the assumed properties of stellar populations and ISM (gas and dust). Young, massive stars show blue UV continua, which become bluer with decreasing metallicity. The effects of stellar age and metallicity on the UV continuum emission over time in our spectral evolution model are illustrated in Figure 10, which shows that the bluest UV spectra are obtained for very young ($\lesssim 1$ Myr) stellar populations with sub-solar metallicities (dashed lines). Dust reddens the emitted stellar spectrum, and leads to a relation between the attenuation suffered by a galaxy and its UV slope β (Meurer et al. 1999). Recombination-continuum from ionized hydrogen also reddens the UV continuum emission emerging from a galaxy. This effect is shown by the solid lines in Figure 10, which illustrates how our combined stellar population + photoionization model predicts redder β slopes than a model accounting only for stellar emission.

Zoom In Zoom Out Reset image size

In order to avoid unphysical values for the β slopes associated to our mock galaxies through the redshift-dependent β–${M}_{{\rm{UV}}}$ relation presented in Section 3.3.1 above, we impose a limit ${\beta }_{\min }=-2.6$ for the bluest possible value. This limit corresponds, approximately, to the bluest β value obtained with beagle for a model with constant SFH, sub-solar metallicity, low depletion factor (${\xi }_{{\rm{d}}}=0.1$), and low ionization parameter ($\mathrm{log}{U}_{{\rm{S}}}=-3.5$). We note that models with non-zero escape fraction of H-ionizing photons can reach bluer values.

3.3.3. Scatter in the $\beta -{M}_{{\rm{UV}}}$ Relation across Cosmic Time

We assign a set of spectral properties to each mock galaxy, allowing us to provide filter photometry as well as a full spectrum for each object. The general method is to produce a parent catalog of spectra that can be matched to galaxies in the mock. Where possible we produce this parent catalog from the results of SED fitting to galaxies in the 3D-HST catalog, allowing the observed photometry to provide the diversity of observed SEDs at given stellar mass and redshift (see Section 2.3). We limit the use of these empirical SEDs to $z\leqslant 4$ galaxies, as beyond this redshift the rest-frame optical is only sampled by the Spitzer IRAC bands where the poor resolution leads to a significant confusion of sources at faint magnitudes.

For galaxies at redshifts z > 4, or $z\leqslant 4$ but below the mass completeness limits of 3D-HST, we rely on extrapolations of observed relationships between ${M}_{{\rm{UV}}}$–${M}_{\star }$ and β–${M}_{{\rm{UV}}}$ to provide constraints on the galaxy SEDs. These constraints are used to match to a parent catalog, produced using beagle in mock catalog mode. When generating this parent catalog we have the full parameter space of the stellar and nebular templates to choose from (described in Section 2.3), and so we use observed trends in galaxy physical parameters, albeit with large scatter, to restrict this parameter space. Specifically we use three observed relations: ${M}_{\star }\mbox{--}Z\mbox{--}\psi$, where $\psi$ is the SFR of the object (the "fundamental metallicity relation"); $\psi \mbox{--}Z\mbox{--}{\hat{\tau }}_{{\rm{V}}}$, to provide physically motivated constraints on dust attenuation (where ${\hat{\tau }}_{{\rm{V}}}$ is the effective V band optical depth); and Z−$\mathrm{log}{U}_{{\rm{S}}}$.

3.4.1. SED Fitting to 3D-HST Catalogs

The photometric catalogs produced by the 3D-HST team (Skelton et al. 2014) are selected from the noise-equalized combination of HST/WFC3 J₁₂₅, ${{JH}}_{140}$, and H₁₆₀ images taken from an extensive set of publicly available imaging data over five fields (AEGIS, COSMOS, GOODS-North, GOODS-South, and the UKIDSS UDS) covering $\simeq 900$ arcmin². From these catalogs we use the data in the deeper regions of the CANDELS (Grogin 2011; Koekemoer et al. 2011) GOODS-South and GOODS-North fields to provide a parent catalog of redshift and mass-dependent SEDs that can be assigned to our mock catalog galaxies.

We use version 4.1 of the 3D-HST photometric catalogs (Momcheva et al. 2016). These catalogs include selection from mosaics that include HUDF-09 (11563; PI:Illingworth) and HUDF-12 (12498; PI: Ellis) WFC3 imaging in the HUDF and parallels (11563; PI: Illingworth) that was performed as part of release 3.0.²⁰ The catalogs do not include deeper HUDF ACS imaging of Beckwith et al. (2006), and so ACS photometry across the GOODS-South is at the depths of the original GOODS imaging (Giavalisco et al. 2004).

From the GOODS-South and GOODS-North 3D-HST catalogs we fit to the broad-band HST fluxes (B₄₃₅, V₆₀₆, i₇₇₅, z₈₅₀, J₁₂₅, JH₁₄₀, and H₁₆₀), as well as the Spitzer/IRAC Channel 1 (3.6 μm) and Channel 2 (4.5 μm) imaging from SEDS (Ashby et al. 2013) to provide constraints in the rest-frame optical at high redshifts. We also use a subset of the ground-based filters that required small photometric zeropoint corrections in the SED fitting analysis of Skelton et al. (2014) compared to the H₁₆₀ band (see their Table 11). In the GOODS-South field we use photometry from VLT/ISAAC J, H, and K_s band imaging from the ESO/GOODS and FIREWORKS surveys (Wuyts et al. 2008; Retzlaff et al. 2010), and VLT/VIMOS U band imaging from the ESO/GOODS survey (Nonino et al. 2009). In the GOODS-North field we use Subaru/MOIRCS imaging in J, H, and K_s bands from the MODS survey (Kajisawa et al. 2011).²¹ We do not apply the zeropoint offsets reported in Skelton et al. (2014), Table 11, after verifying that applying these corrections does not improve the accuracy of photometric redshifts output by beagle.

We fit the broad-band photometry using beagle (see Section 2.3 for details on the model). We use a delayed star formation history $\psi (t)\propto t\exp (-t/{\tau }_{{\rm{SFR}}})$, where ${\tau }_{{\rm{SFR}}}$ is the star formation timescale and t the age of the galaxy, taken to lie between 10⁷ years and the maximum time allowed since the onset of star formation at the galaxy redshift. This parametrization gives a star formation history that rises at early times and declines exponentially at later times. This star formation history is shown to better reproduce the colors and mass-to-light ratios of galaxies in the smoothed particle hydrodynamics (SPH) simulations of Simha et al. (2014) than the widely used exponentially decreasing star formation histories. Additionally, simulations have been shown to predict that high-redshift galaxies have rising star formation histories (Finlator et al. 2011), a scenario that is naturally achieved using this parametrization. To ensure that galaxies are not fitted with models that are older than the age of the universe, we set an upper limit of ${z}_{\mathrm{form}}^{\max }=15$ to the redshift of onset of star formation. We further employ a weakly informative Gaussian prior on $\mathrm{log}t$, with mean at $\mathrm{log}(t/\mathrm{year})=9.3$ and $\sigma =0.7$.²² We approximate the distributions of stellar and interstellar metallicities in a galaxy with a single metallicity ${Z}_{{\rm{ISM}}}=Z$. We use an exponential prior for ${\hat{\tau }}_{V}$ and fix μ = 0.4. The free parameters in the model fitting are summarized in Table 2, and see Section 2.3 for a general overview of the individual parameters.

Table 2.  Parameters Allowed to Vary in the beagle Fitting to Galaxies in the 3D-HST Catalog with Their Priors

Parameter	Prior	Description
z	Uniform $\in \,[0,15]$	Redshift
$\mathrm{log}({M}_{\mathrm{tot}}/{M}_{\odot })$ a	Uniform $\in \,[7,13]$	Integrated SFH
$\mathrm{log}(t/\mathrm{year})$	Gaussian ${ \mathcal N }(9.3;0.7)$ truncated $\in \,[7,10.15]$	Age of oldest stars in the galaxy
$\mathrm{log}({\tau }_{{\rm{SFR}}}/\mathrm{year})$	Uniform $\in \,[7,12]$	Timescale of star formation
$\mathrm{log}(Z/{Z}_{\odot })$	Uniform $\in \,[-2.2,0.24]$	Stellar (and interstellar) metallicity ($Z$ = ${Z}_{{\rm{ISM}}}$)
${\hat{\tau }}_{{\rm{V}}}$	Exponential exp(−${\hat{\tau }}_{{\rm{V}}}$) truncated $\in \,[0,4]$	V band attenuation optical depth
μ	Fixed 0.4	Fraction of attenuation arising in the diffuse ISM
$\mathrm{log}{U}_{{\rm{S}}}$	Uniform $\in \,[-4,-1]$	Effective gas ionization parameter
${\xi }_{{\rm{d}}}$	Fixed 0.3	Dust-to-metal mass ratio

Note.

^abeagle samples over the integral of the past star formation history of the galaxy (${M}_{\mathrm{tot}}$). It returns the stellar mass (${M}_{\star }$), which accounts for the mass returned by evolved stars to the ISM.

Download table as:  ASCIITypeset image

3.4.2. Parent SED Catalog of Galaxies Based on 3D-HST Catalog for $z\leqslant 4$, $\mathrm{log}({M}_{\star }/{M}_{\odot })\gt 8$

To generate our parent SED catalog for high-mass galaxies at z < 4, we use fits to the broad-band photometry of star-forming galaxies in the 3D-HST catalog using beagle. Since T14 use rest-frame U − V versus V − J colors to separate galaxies into star-forming and quiescent galaxies before measuring the type-dependent mass functions, we select star-forming galaxies from the 3D-HST catalog using a similar star-forming/quiescent (SF/Q) classification scheme. Specifically, we select objects in the star-forming region of U − V versus V − J color space as defined by Whitaker et al. (2011; see their Figure 17 and Equations (14) and (25)) using the reported U, V, and J band absolute rest-frame magnitudes supplied in the 3D-HST catalog.

As described in Section 2.3, beagle uses MultiNest (Feroz et al. 2009) to sample the parameter space, and records the associated SEDs and MultiNest weights. This information can be used to produce samples drawn from the corresponding posterior probability distribution, and we use these samples to populate a parent catalog with a range of statistically acceptable SED fits (plus associated physical parameters) for each object.

At low masses, constraints on the rest-frame UV, metallicity, and dust attenuation can suffer from poor photometric constraints. For each of these parameters we therefore impose conditional priors on ${M}_{{\rm{UV}}}$–${M}_{\star }$, ${M}_{\star }\mbox{--}Z\mbox{--}\psi$, and $\psi \mbox{--}Z\mbox{--}{\hat{\tau }}_{{\rm{V}}}$. These are additional priors to those already set in the SED fitting (listed in Table 2) that we must apply after the fact, as beagle does not accept conditional priors. Luckily it is relatively simple to apply these priors as one effectively needs to adjust the MultiNest weights derived from the initial beagle fitting. This method was presented in Chevallard et al. (2017; see their Section 2.3, where they describe drawing SEDs that match the observed main sequence of star-forming galaxies). Essentially we perform weighted draws from the MultiNest output for each galaxy fitted to in the 3D-HST catalog such that the probability of a galaxy entering into the parent catalog follows:

$$\begin{eqnarray}&&P({M}_{{\rm{UV}}},Z,{\hat{\tau }}_{{\rm{V}}}| {M}_{\star },\psi ,z)\\ &&\qquad \propto \,P({M}_{{\rm{UV}}}| {M}_{\star },z)\,P(Z| {M}_{\star },\psi )\,P({\hat{\tau }}_{{\rm{V}}}| Z,\psi )\,P({\boldsymbol{\Theta }}| {\boldsymbol{D}}),\end{eqnarray} \tag{ 14 }$$

where $P({\boldsymbol{\Theta }}| {\boldsymbol{D}})$ is the posterior probability of a set of parameters sampled over in the fitting (${\boldsymbol{\Theta }}=[z$, ${M}_{\mathrm{tot}}$, $t$, ${\tau }_{{\rm{SFR}}}$, $Z$, ${\hat{\tau }}_{{\rm{V}}}$, μ, $\mathrm{log}{U}_{{\rm{S}}}$, ${\xi }_{{\rm{d}}}]$) given the data, D. The shape of the prior in ${M}_{{\rm{UV}}}$–${M}_{\star }$ is set by the observational constraints detailed in Section 3.2, and for $P(Z| {M}_{\star },\psi )$ and $P({\hat{\tau }}_{{\rm{V}}}| Z,\psi )$ we use the same priors imposed when constructing the parent catalog not based on 3D-HST photometry (as described in the following sub-section; specifically Equations (15)–(16) describe the prior in $P(Z| {M}_{\star },\psi )$ and Equations (19)–(22) describe the prior in $P({\hat{\tau }}_{{\rm{V}}}| Z,\psi )$).²³ The prior on ${M}_{{\rm{UV}}}$–${M}_{\star }$ ensures a smooth transition in observables between galaxies assigned from the 3D-HST catalog constraints and those at lower masses or higher redshift that lie outside the parameter space covered by the observed galaxies. For objects in the 3D-HST catalog with firm ${M}_{{\rm{UV}}}$ constraints, the prior changes the sampling very little, allowing the final mock ${M}_{{\rm{UV}}}$–${M}_{\star }$ relation to display the wide diversity of SEDs seen in the observed population.

When drawing SEDs from the fits to 3D-HST galaxies, we account for the varying sensitivity limits of the adopted data sets. The HUDF field and parallels provide the faintest H₁₆₀ objects in the catalogs, while the sensitivity in CANDELS varies between the Deep and Wide regions. To avoid favoring SEDs from the numerous H₁₆₀-bright objects found in shallow, wider survey areas, we limit the shallower area to that in the CANDELS deep regions within GOODS-N and GOODS-S. Figure 11 displays the residual difference in depths between the deeper HUDF and parallel fields and the CANDELS deep regions of GOODS-S and GOODS-N (plotted together). The number counts at given total H₁₆₀ magnitude start to turn over at H₁₆₀ ∼ 27 for the GOODS-N+GOODS-S region, while the number counts for the HUDF and parallels remain flat to ${H}_{160}\sim 27.5$. We apply a conservative cut of H₁₆₀ = 26.8, above which only objects in the deeper HUDF and parallels are included in the parent catalog and below which only objects from GOODS-N+GOODS-S are included. We also require that the transition in number counts across this magnitude limit is smooth, ensuring that the distribution of H₁₆₀ magnitudes at a given mass is not more heavily weighted by the more numerous objects in the shallower region. The ratio between cumulative number counts at ${H}_{160}\lt 26.7$ for the two regions is equal to 7.5, which essentially accounts for the difference in area between the two regions. (This factor has already been applied to the number counts in the HUDF+parallels in the figure.) We therefore draw 45 realizations per object from the HUDF+parallel regions, while only 5 random draws are included in the parent catalog for each object in the shallower region.

Zoom In Zoom Out Reset image size

3.4.3. Parent Catalog for Galaxies at $z\leqslant 4$ and with $\mathrm{log}({M}_{\star }/{M}_{\odot })\lt 8$, and for All z > 4 Galaxies

For mock galaxies outside of the stellar mass and redshift range of 3D-HST galaxies, we generate a parent catalog by using beagle to produce SEDs covering a wide range of parameter values (e.g., z, ${M}_{\star }$, τ, ${\hat{\tau }}_{{\rm{V}}}$, $\mathrm{log}{U}_{{\rm{S}}}$, $Z$, and t). We impose constraints to avoid populating the parent catalog with objects in unphysical regions of parameter space, such as galaxies with high metallicity and high $\mathrm{log}{U}_{{\rm{S}}}$, or high-mass galaxies with extremely small metallicities. To constrain these parameters, we utilize the distributions of ${M}_{\star }\mbox{--}Z\mbox{--}\psi$ and $Z$−$\mathrm{log}{U}_{{\rm{S}}}$ inferred from observations.

We use delayed SFHs to produce the SEDs, where the τ and t values are chosen to ensure that the galaxies would be classified as star-forming, with $\mathrm{log}({\psi }_{{\rm{S}}}/{\mathrm{yr}}^{-1})\gt -10$. Here we allow $\mathrm{log}(t/\mathrm{year})$ to vary between six and the age of the universe. The lower limit in age is lower than that introduced in the prior on $\mathrm{log}(t/\mathrm{year})$ used in the fitting to objects in the 3D-HST catalog (see Table 2), as this parent catalog is going to be used to match to lower mass/higher redshift objects. beagle provides the stellar masses accounting for mass returned to the ISM as stars evolve and die, as well as the current SFR, which can be used to assign metallicity, ionization parameter, and V band optical depth due to attenuation by dust.

To constrain the metallicities of galaxies, we use the fundamental metallicity relation between ${M}_{\star }\mbox{--}Z\mbox{--}\psi$ measured by Hunt et al. (2016) from a compilation of ∼1000 galaxies covering a wide range in $\psi$ and stellar mass, with oxygen abundance estimates derived from consistent calibrations and, crucially, covering a wide range in redshifts up to z ∼ 3.7. The redshift distribution of Hunt et al. (2016) can be viewed in their Figure 1, which can be compared with the fundamental metallicity relation measured in Mannucci et al. (2010) from galaxies in the Sloan Digital Sky Survey for objects with redshifts between 0.07 and 0.3. A fit to this fundamental metallicity relation is given by

$$\begin{eqnarray}12+\mathrm{log}({\rm{O}}/{\rm{H}}) & = & -0.14\,\mathrm{log}(\psi /{M}_{\odot }\,{\mathrm{yr}}^{-1})\\ & & +0.37\,\mathrm{log}(M/{M}_{\odot })+4.82.\end{eqnarray} \tag{ 15 }$$

This relation is a fit to the gas-phase oxygen abundance, while we require a relation based on the nebular metallicity ${Z}_{{\rm{ISM}}}=Z$ associated with our models. From the grid of ${\xi }_{{\rm{d}}}=0.3$ models used here, we infer the approximate relation $12+\mathrm{log}({\rm{O}}/{\rm{H}})\,\simeq \mathrm{log}({Z}_{{\rm{ISM}}}/{Z}_{\odot })+8.7$. While this approximation is not suitable to determine accurately $12+\mathrm{log}({\rm{O}}/{\rm{H}})$ from the output metallicity for individual nebular models (see Gutkin et al. 2016, Table 2, for values of $12+\mathrm{log}({\rm{O}}/{\rm{H}})$ for different model metallicities), the errors it generates are much smaller than the scatter we introduce below in the fundamental metallicity relation.

It is important to highlight here that the Hunt et al. (2016) relation only models a linear dependence between oxygen abundance and stellar mass, whereas we know the mass–metallicity relation is not linear at high stellar masses (e.g., Tremonti et al. 2004). In fact the upper limit in metallicity of $Z/{Z}_{\odot }=0.24$ (also the upper limit of the prior in metallicity employed when fitting to objects in the 3D-HST catalogs) introduces a turnover in the catalog mass–metallicity relation at low redshifts and high masses (see Section 6.4).

We wish to instill a broad diversity in our parent catalog SEDs and avoid an over-representation of unphysical parameters in our resulting mock catalogs. We therefore apply a broad scatter to the fundamental metallicity relation, and do not attempt to predict the form of the ${M}_{\star }$−$Z$ relation or ${M}_{\star }\mbox{--}Z\mbox{--}\psi$ plane to high redshifts. The broad range of spectral parameters will enable investigations of selection effects in future observations, especially for redshift and magnitude regimes where current measurements cannot yet reach. We characterize the scatter with a student's t-distribution:

$$\begin{eqnarray}&&f(x)=\displaystyle \frac{{\rm{\Gamma }}(\tfrac{\nu +1}{2})}{\sqrt{\nu \pi }{\rm{\Gamma }}(\tfrac{\nu }{2})}{\left(1+\displaystyle \frac{{x}^{2}}{\nu }\right)}^{-\tfrac{\nu +1}{2}},\end{eqnarray} \tag{ 16 }$$

where ν is the number of degrees of freedom, Γ is the gamma function, and

$$\begin{eqnarray}&&x=\displaystyle \frac{\mathrm{log}({Z}_{{\rm{ISM}}}/{Z}_{\odot })-\mathrm{log}({\bar{Z}}_{{\rm{ISM}}}/{Z}_{\odot })}{{\sigma }_{x}}.\end{eqnarray} \tag{ 17 }$$

We set ${\sigma }_{x}=0.3$ and ν = 3, where ν = 3 has been chosen to provide a distribution with more weight in the tails compared to a Gaussian.

To constrain the ionization parameter of these galaxies, we use a linear fit between metallicity and $\mathrm{log}{U}_{{\rm{S}}}$ measurements at low redshift from Carton et al. (2017, see their Figure 2):

$$\begin{eqnarray}&&\mathrm{log}{U}_{{\rm{S}}}(Z)=-0.8\,\mathrm{log}(Z/{Z}_{\odot })-3.58.\end{eqnarray} \tag{ 18 }$$

We again use the student's t-distribution with three degrees of freedom to introduce scatter in this relation.

We account for dust attenuation by using an approach commonly featured in semi-analytic models of galaxy formation (e.g., Guiderdoni & Rocca-Volmerange 1987; De Lucia & Blaizot 2007; Fontanot et al. 2009). Following Devriendt et al. (1999, their Equation (6)), we estimate the V band, face-on attenuation optical depth ${\hat{\tau }}_{{\rm{V}}\perp }$ using the relation

$$\begin{eqnarray}&&{\hat{\tau }}_{{\rm{V}}\perp }={\left(\displaystyle \frac{{Z}_{{\rm{ISM}}}}{{Z}_{\odot }}\right)}^{1.6}\left(\displaystyle \frac{{N}_{{\rm{H}}}}{2.1\times {10}^{21}\,{\mathrm{cm}}^{-2}}\right),\end{eqnarray} \tag{ 19 }$$

where ${Z}_{{\rm{ISM}}}$ is the interstellar metallicity and ${N}_{{\rm{H}}}$ the mean hydrogen column density. As in Devriendt et al. (1999), we compute ${N}_{{\rm{H}}}$ from the (cold) gas fraction

$$\begin{eqnarray}&&{N}_{{\rm{H}}}=6.8\times {10}^{21}\,\displaystyle \frac{{M}_{\mathrm{gas}}}{{M}_{\star }+{M}_{\mathrm{gas}}},\end{eqnarray} \tag{ 20 }$$

where the cold gas mass ${M}_{\mathrm{gas}}$ is computed by inverting the Schmidt–Kennicutt relation (Schmidt 1959; Kennicutt 1998). In practice, we consider the size and SFR of each galaxy to compute the SFR density ${{\rm{\Sigma }}}_{\psi }=\psi /(\pi {r}^{2})$, in units of ${M}_{\odot }\,{\mathrm{yr}}^{-1}\,{\mathrm{pc}}^{-2}$, where r is the galaxy effective radius. We can then compute the cold gas surface density ${{\rm{\Sigma }}}_{\mathrm{gas}}$ from the Schmidt–Kennicutt relation, and estimate the cold gas mass as ${M}_{\mathrm{gas}}={{\rm{\Sigma }}}_{\mathrm{gas}}\,\pi {r}^{2}$.

Equation (19) provides us with the face-on attenuation optical depth, from which we can derive an angle-averaged attenuation optical depth by assuming a spatial distribution of dust and stars. Following Devriendt et al. (1999), we approximate our galaxies as oblate ellipsoids where dust and stars are homogeneously mixed. The V band attenuation optical depth averaged over all galaxy inclinations i can then be written as

$$\begin{eqnarray}&&\langle {\hat{\tau }}_{{\rm{V}}}{\rangle }_{i}=-2.5\,\mathrm{log}\left(\displaystyle \frac{{a}_{{\rm{V}}}}{1-{\omega }_{{\rm{V}}}+{\omega }_{{\rm{V}}}\,{a}_{{\rm{V}}}}\right)\,\displaystyle \frac{1}{1.086},\end{eqnarray} \tag{ 21 }$$

where ${\omega }_{{\rm{V}}}=0.87$ is the albedo at 5500 Å for dust grains with properties as in the Small Magellanic Cloud (Pei 1992), and ${a}_{{\rm{V}}}$ is computed as

$$\begin{eqnarray*}&&{a}_{{\rm{V}}}=\displaystyle \frac{3}{4\,{\hat{\tau }}_{{\rm{V}}\,\perp }^{{\prime} }}\left[1-\displaystyle \frac{1}{2\,{\hat{\tau }}_{{\rm{V}}\,\perp }^{{\prime} \,2}}+\left(\displaystyle \frac{1}{{\hat{\tau }}_{{\rm{V}}\,\perp }^{{\prime} }}+\displaystyle \frac{1}{2\,{\hat{\tau }}_{{\rm{V}}\,\perp }^{{\prime} \,2}}\right)\exp (-2\,{\hat{\tau }}_{{\rm{V}}\,\perp }^{{\prime} })\right],\end{eqnarray*}$$

where ${\hat{\tau }}_{{\rm{V}}\,\perp }^{{\prime} }=2.62\,{\hat{\tau }}_{{\rm{V}}\perp }$.

Equations (19)–(21) enable us to associate to each mock galaxy a physically motivated value for the angle-averaged V band attenuation optical depth, which depends on the galaxy SFR, size, and metallicity. We then account for the effect of galaxy inclination on the V band attenuation optical depth ${\hat{\tau }}_{{\rm{V}}}$ by randomly drawing ${\hat{\tau }}_{{\rm{V}}}$ from a Gaussian distribution centered at $\langle {\hat{\tau }}_{{\rm{V}}}{\rangle }_{i}$ and truncated at ${\hat{\tau }}_{{\rm{V}}}=0$. The width of the Gaussian used to draw the ${\hat{\tau }}_{{\rm{V}}}$ value is chosen to be dependent on $\langle {\hat{\tau }}_{{\rm{V}}}{\rangle }_{i}$ and $Z$, according to

$$\begin{eqnarray}&&\sigma =0.2+Z/{Z}_{\odot }-\langle {\hat{\tau }}_{{\rm{V}}}{\rangle }_{i}.\end{eqnarray} \tag{ 22 }$$

This function ensures that at low metallicities, there is a smaller scatter in ${\hat{\tau }}_{{\rm{V}}}$, limiting range of dust attenuation in the regime where we expect low dust-to-gas ratios, while the minimum value of 0.2 prevents the values of ${\hat{\tau }}_{{\rm{V}}}$ being too constrained at the lowest metallicities. We note that non-negligible attenuation by dust even at very low metallicities is not unreasonable, as even if they have low gas-to-dust ratios they may be gas-rich, allowing for non-negligible dust-to-stellar mass ratios (e.g., da Cunha et al. 2010). Our choice of the negative dependence of σ on $\langle {\hat{\tau }}_{{\rm{V}}}{\rangle }_{i}$ mimics results obtained from radiative transfer calculations of dust attenuation in galaxies (e.g., Pierini et al. 2004; Tuffs et al. 2004), which show that galaxies with low angle-averaged attenuation optical depths $\langle {\hat{\tau }}_{{\rm{V}}}{\rangle }_{i}\lesssim 0.1$ exhibit a larger fractional range of inclination-dependent attenuations than galaxies with larger $\langle {\hat{\tau }}_{{\rm{V}}}{\rangle }_{i}$.

3.4.4. Matching Mock Galaxies to the Parent Catalog

To assign SEDs from the parent catalog to mock galaxies, we find the closest match between the mock galaxy parameters and the physical parameters of the parent catalog galaxies. Which properties are matched differ when assigning SEDs drawn from fits to 3D-HST galaxies or from the wide grid of galaxies produced by beagle. For those objects with z < 4 and stellar mass higher than $\mathrm{log}(M/{M}_{\odot })=8$ or the mass completeness limit of the 3D-HST catalog ($\mathrm{log}(M/{M}_{\odot })\,\gt 0.7z+0.63$;²⁴ red vertical lines shown in Figure 6), we find the closest match in ${M}_{\star }$ and redshift, as we rely on observed broad-band photometry to constrain the expected SED shape of an object at a given stellar mass. For objects at z > 4 or with lower stellar masses at $z\leqslant 4$, the expected SED shape is based on extrapolations of the ${M}_{{\rm{UV}}}$–${M}_{\star }$ and β–${M}_{{\rm{UV}}}$ relations, as described in Sections 3.2 and 3.3, respectively. These relations are used to assign ${M}_{{\rm{UV}}}$ and β values to each galaxy in the star-forming galaxy mock catalog. Each mock object is then assigned an SED based on the closest match in redshift, stellar mass, ${M}_{{\rm{UV}}}$, and UV slope in the parent catalog.

For all matches between mock galaxies and parent catalog, it is then possible to shift the redshift of the SED to the exact redshift of the mock galaxy. Figure 12 displays the distributions of physical parameters assigned to all star-forming galaxies in the mock catalog.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

An overview of our method of generating mock quiescent galaxies is outlined in Figure 13, and as with star-forming galaxies, is based on a model for the redshift evolution of the stellar mass function of quiescent galaxies. This model is based on observed stellar mass functions that have been measured in the redshift range $0.2\leqslant z\leqslant 3.0$ by T14. These authors use the redshift-dependent UVJ color selection from Whitaker et al. (2011) to select quiescent galaxies from the ZFOURGE medium-band photometric survey (Straatman et al. 2016) and data from the CANDELS survey (see Section 3.1). T14 find that the quiescent galaxy stellar mass function is best fitted by a double-Schechter function (Equation (4) above) in the redshift range $0.2\leqslant z\leqslant 1.5$, and by a single-Schechter function at higher redshifts, $1.5\leqslant z\leqslant 3.0$. The differing functional forms were chosen to match the observed upturn in the stellar mass function below $\mathrm{log}(M/{M}_{\odot })\leqslant 9.5$ at $z\leqslant 1$, a result in line with observations by Santini et al. (2012), Muzzin et al. (2013), and Ilbert et al. (2013).

In this work, we use both the observed binned stellar mass functions and the fitted Schechter parameters in bins of redshift from T14 to construct a continuous model for the redshift evolution of the quiescent stellar mass function. To produce smooth evolution at all redshifts, we choose to adopt the double-Schechter function description for the mass function at all redshifts, even at z > 1.5 where observations find consistency with a single Schechter function. This double-Schechter function has five parameters: ${M}_{{\rm{M}}}^{* }$, ${\alpha }_{1,{\rm{M}}}$, ${\phi }_{1,{\rm{M}}}^{* }$, ${\alpha }_{2,{\rm{M}}}$, and ${\phi }_{2,{\rm{M}}}^{* }$. We fit the double-Schechter parameters from T14 at z < 1.5, but substituted the single Schechter parameters (M^*, α, and ${\phi }^{* }$) for the double-Schechter parameters (${M}_{{\rm{M}}}^{* }$, ${\alpha }_{1,{\rm{M}}}$, ${\phi }_{1,{\rm{M}}}^{* }$) at z > 1.5. In a double-Schechter function, ${\alpha }_{2,{\rm{M}}}$ and ${\phi }_{2,{\rm{M}}}^{* }$ control the slope and overall normalization at low masses.

We choose to fix ${M}_{{\rm{M}}}^{* }$ to a value of $\mathrm{log}(M/{M}_{\odot })\sim 10.6$ at all redshifts, owing to the lack of significant shift in the observed evolution in ${M}_{{\rm{M}}}^{* }$ with redshift from the T14 observations, and we fit the ${\alpha }_{1,{\rm{M}}}$ and ${\phi }_{1,{\rm{M}}}^{* }$ evolution with a cubic function with redshift. We fit a linear function to the ${\alpha }_{2,{\rm{M}}}$ evolution and a quadratic function to the ${\phi }_{2,{\rm{M}}}^{* }$ evolution.

Given the chosen forms for the evolution in the parameters, in order to prevent the cubic functions from diverging at low redshift and high redshift, we stop the evolution for some of the parameters. For ${\alpha }_{1,{\rm{M}}}$, ${\phi }_{1,{\rm{M}}}^{* }$, ${\alpha }_{2,{\rm{M}}}$, and ${\phi }_{2,{\rm{M}}}^{* }$, we stop the redshift evolution at z < 0.5, z < 0.75, z < 0.5, and z < 0.5, respectively. We require ${\phi }_{1,{\rm{M}}}^{* }$ to stay constant below z < 0.75 because the cubic fit to the T14 Schechter parameters would otherwise over-predict by an order of magnitude the number of objects, as compared to observations at low redshift. We fix ${\alpha }_{1,{\rm{M}}}$ at z > 1.75 to the value of the cubic fit at that redshift, and at z > 1.88 we impose a linear decline on the evolution of ${\phi }_{1,{\rm{M}}}^{* }$ to extrapolate the quiescent galaxy counts to z > 3, beyond the range probed by the observations. This projected decline in ${\phi }_{1,{\rm{M}}}^{* }$ predicts a number density of 2.16 ± 0.8 $\times \,{10}^{-5}$ Mpc⁻³ quiescent galaxies at $3.4\leqslant z\leqslant 4.2$, consistent with the number density of quiescent galaxies identified in Straatman et al. (2014, 2015, 1.8 ± 0.8 × 10⁻⁵ Mpc⁻³) to their stellar mass limit of $\mathrm{log}(M/{M}_{\odot })\sim 10.6$ and the ZFOURGE survey area of ∼363 arcmin². Thus our extrapolation is in broad agreement with the few constraints available on counts of quiescent galaxies at z > 3.5.

The resulting quiescent galaxy stellar mass function is compared with the T14 data in Figure 14, and we provide the detailed functional fit parameters in Table 3. Our model for the stellar mass function evolution broadly agrees within the uncertainties of the binned stellar mass function observations from T14. We note that the observations imply more low-mass $[\mathrm{log}(M/{M}_{\odot })\leqslant 9.5]$ quiescent galaxies at $0.5\leqslant z\leqslant 0.75$ than below $z\leqslant 0.5$. The continuously varying stellar mass function evolution model we describe here monotonically increases in counts with decreasing redshift, and therefore it does not replicate this rise and fall of the low-mass end that is seen in the observations (Figure 14).

Zoom In Zoom Out Reset image size

Table 3.  Quiescent Mass Function Double-Schechter Parameters (Figure 14) and Their Evolution

Parameter	Redshift Range	Functional Form
${M}_{{\rm{M}}}^{* }$	$z\geqslant 0.2$	${M}_{{\rm{M}}}^{* }=10.617$
${\alpha }_{1,{\rm{M}}}$	z < 0.5	${\alpha }_{1,{\rm{M}}}=-0.225$
	$0.5\leqslant z\lt 1.75$	${\alpha }_{1,{\rm{M}}}={b}_{0}\times {z}^{3}+{b}_{1}\times {z}^{2}+{b}_{2}\times z+{b}_{3}$
		${b}_{0}=0.43$, ${b}_{1}=-2.33$, b2=3.49, ${b}_{3}=-1.44$
	$z\geqslant 1.75$	${\alpha }_{1,{\rm{M}}}=-0.150$
${\phi }_{1,{\rm{M}}}^{* }$	z < 0.75	$\mathrm{log}({\phi }_{1,{\rm{M}}}^{* })=-2.67$
	$0.75\leqslant z\lt 1.877$	$\mathrm{log}({\phi }_{1,{\rm{M}}}^{* })={c}_{0}\times {z}^{3}+{c}_{1}\times {z}^{2}+{c}_{2}\times z+{c}_{3}$
		c0 = −0.35, c1 = 1.92, c2 = −3.77, c3 = −0.78
	$z\geqslant 1.877$	$\mathrm{log}({\phi }_{1,{\rm{M}}}^{* })={c}_{4}\times z+{c}_{5}$
		c4 = − = 0.43, c5 = −2.59
${\alpha }_{2,{\rm{M}}}$	z < 0.5	${\alpha }_{2,{\rm{M}}}=-1.83$
	z ≥ 0.5	${\alpha }_{2,{\rm{M}}}={d}_{0}\times z+{d}_{1}$
		d0 = 1.15, d1 = −2.41
${\phi }_{2,{\rm{M}}}^{* }$	z < 0.5	$\mathrm{log}({\phi }_{2,{\rm{M}}}^{* })=-4.71$
	z ≥ 0.5	$\mathrm{log}({\phi }_{2,{\rm{M}}}^{* })={e}_{0}\times {z}^{2}+{e}_{1}\times z+{e}_{2}$
		e0 = −0.59, e1 = 1.93, e2 = −5.52

Note. The parameter is shown on the left column, the functional form of the parameter (with associated constants) on the right column, and the redshift range of that functional form in the middle column.

Download table as:  ASCIITypeset image

We apply the same SED assignment technique for quiescent galaxies as that for star-forming galaxies, described in Section 3.4, whereby we produce a large parent catalog of SEDs and associated physical parameters that can then be used to associate an SED with each mock galaxy. Where possible, we use observed galaxies to guide the allowed range of mock galaxy SEDs at a given stellar mass. Outside of the parameter space covered by 3D-HST objects, we produce theoretical mock SEDs using beagle.

To produce a parent catalog from observed objects, we use beagle fits to 3D-HST galaxies classified as quiescent using the rest-frame U, V, and J band absolute magnitudes and the Whitaker et al. (2011) U − V versus V − J color criteria. Following the method described in Section 3.4.2, we draw SEDs (and associated physical parameters) from the beagle fits to produce the parent catalog. Specifically, the stellar mass and associated 68% central credible interval are estimated for each galaxy, and the SEDs are randomly drawn from the BEAGLE output files if their stellar mass lies within this interval, without further weighting based on other physical parameters.

Uncertainties associated with the star-forming/quiescent separation have to be dealt with carefully when producing this parent catalog. More details are given in Appendix B, where we define a redshift-dependent H₁₆₀ magnitude cut of ${H}_{160}\lt 24.5+z$, brighter than which quiescent objects from the 3D-HST catalog can be used to provide SEDs for the parent catalog. This limit lies well above the sensitivity limit of the CANDELS deep region of the 3D-HST catalog (and so we do not encounter the situation we found for the star-forming galaxies, where the faintest objects were only drawn from the smaller area HUDF portion of the catalog). We therefore skip the areal correction described in Section 3.4.2, and populate the 3D-HST parent catalog of quiescent galaxies with five random draws from the beagle fit output for each quiescent galaxy in the 3D-HST catalog satisfying the ${H}_{160}\lt 24.5+z$ criterion. The H₁₆₀ limit translates to an approximate mass limit of $\mathrm{log}(M/{M}_{\odot })\gt 8.7+0.4z$. Thus only objects in the mock catalog with $\mathrm{log}(M/{M}_{\odot })\gt 8.7+0.4z$ are paired with the closest match in redshift and stellar mass among parent catalog galaxies.

For objects with $\mathrm{log}(M/{M}_{\odot })\lt 8.7+0.4z$, we produce a parent catalog of SEDs and associated physical parameters using beagle. For quiescent galaxies, we do not need to assign nebular H ii-region parameters (e.g., ${\xi }_{{\rm{d}}}$, $\mathrm{log}{U}_{{\rm{S}}}$), and we neglect dust attenuation. We therefore vary only the galaxy age t, star formation timescale ${\tau }_{{\rm{SFR}}}$, and metallicity $Z$ to generate SEDs for the parent galaxy catalog. Using a delayed SFH (Section 3.4.1) with $\mathrm{log}({\tau }_{{\rm{SFR}}})\lt 1.11\times \mathrm{log}(t)-2.02$ ensures that the sSFR of objects be less than $\mathrm{log}({\psi }_{{\rm{S}}}/{\mathrm{yr}}^{-1})\,=-10$. The parameters t and ${\tau }_{{\rm{SFR}}}$ are assigned to each mock catalog galaxy from uniform distributions. The parameter t is allowed to vary between 30 Myr and the age of the universe at the redshift of the object. We allow ${\tau }_{{\rm{SFR}}}$ to vary between 10 Myr and the maximum value required to produce $\mathrm{log}({\psi }_{{\rm{S}}}/{\mathrm{yr}}^{-1})\lt -10$.

Measuring the stellar metallicities of quiescent galaxies and their evolution is technically challenging, requiring deep rest-frame optical spectra to measure stellar absorption-line indices. (JWST will provide new opportunities to probe stellar metallicities of quiescent galaxies at high redshift.) However, with existing data, most stellar metallicity measurements for quiescent galaxies exist only for high-mass field galaxies ($\mathrm{log}(M/{M}_{\odot })\gtrsim 9.5$) to moderate redshifts (e.g., Gallazzi et al. 2006, z ∼ 0.1, and Gallazzi et al. 2014, z ∼ 0.7) or cluster galaxies (e.g., Sánchez-Blázquez et al. 2009; Jørgensen & Chiboucas 2013). The highest mass galaxies in our mock have their physical properties assigned from SED fits to 3D-HST galaxies, and so, with the lack of current constraints on the stellar metallicities of low-mass and high-redshift quiescent galaxies, we assign metallicities with a uniform weight between the limits of our templates ($-2.2\lt \mathrm{log}(Z/{Z}_{\odot })\lt 0.24$).

After each galaxy in the mock catalog is assigned t, ${\tau }_{{\rm{SFR}}}$, and $Z$ values, beagle computes the fraction of mass $1-{M}_{\star }/{M}_{\mathrm{tot}}$ returned by evolved stars to the ISM, and hence the scaling required to generate an SED with the corresponding stellar mass assigned to the mock catalog object. beagle then generates the SEDs in "mock mode" using the assigned parameters. As with the star-forming catalog, the SEDs assigned to each of the quiescent mock catalog galaxies are used to generate NIRCam filter fluxes.

We aim to generate a continuously evolving model in stellar mass, UV luminosity, and redshift, using the observed size–mass relationships that have been measured for galaxies at $0\leqslant z\leqslant 4$. For this purpose, we use the relationships measured in van der Wel et al. (2014) for both star-forming and quiescent galaxies using CANDELS data, and extrapolate their behavior down to $\mathrm{log}(M/{M}_{\odot })$ ∼ 6. These relationships have also been shown to agree with the measured size–mass relation of local galaxies in SDSS (Shen et al. 2003; Guo et al. 2009).

van der Wel et al. (2014) parametrize the redshift evolution of the half-light semimajor radius in kpc, ${R}_{\mathrm{eff},\mathrm{maj}}$, as a power-law function of the Hubble parameter at a given redshift, H(z), in bins of stellar mass. The parametrization has the form

$$\begin{eqnarray}&&{R}_{{\rm{eff,maj}}}(z)={B}_{H}{\left(\displaystyle \frac{H(z)}{{H}_{0}}\right)}^{{\beta }_{H}},\end{eqnarray} \tag{ 23 }$$

where both B_H and ${\beta }_{H}$ (respectively, the coefficient and power-law slope of the redshift evolution) vary with stellar mass. To generate a smoothly evolving model, we generalize this evolution of ${R}_{\mathrm{eff},\mathrm{maj}}(z)$ between stellar mass bins by fitting both B_H and ${\beta }_{H}$ as functions of stellar mass, to produce one smooth function in both redshift and stellar mass.

The behavior of B_H and ${\beta }_{H}$ with stellar mass differs between star-forming galaxies and quiescent, and therefore we model the stellar mass dependence differently between the two samples. For star-forming galaxies, both B_H and ${\beta }_{H}$ appear linear with ${M}_{\star }=\mathrm{log}(M/{M}_{\odot })$. The best-fitting linear relationships are

$$\begin{eqnarray}{B}_{H}({M}_{\star }) & = & 0.23{M}_{\star }-1.61\\ {\beta }_{H}({M}_{\star }) & = & -0.08{M}_{\star }+0.25,\end{eqnarray} \tag{ 24 }$$

and the shape of the resulting function ${R}_{\mathrm{eff},\mathrm{maj}}(z,{M}_{\star })$ for star-forming galaxies is shown in the left panel of Figure 15. We extrapolate the relationship out to z ∼ 4, and assign sizes to all star-forming galaxies in this redshift range using this method.

Zoom In Zoom Out Reset image size

For quiescent galaxies, the behavior of both B_H and ${\beta }_{H}$ with stellar mass are not linear. We find that B_H (${M}_{\star }$) is well described by either a quadratic or exponentially declining function of mass; however, we choose to parametrize B_H using the exponential to avoid the undesirable quadratic feature that galaxy size increases unphysically to low masses. We find that ${\beta }_{H}$ is well fit by an exponentially increasing function with decreasing stellar mass; however, increasingly large values of this power-law exponent at low masses produce unphysical size evolution at low mass. Therefore we fix the value of ${\beta }_{H}$ below $\mathrm{log}(M/{M}_{\odot })\lt 9.75$. The relationships we use for quiescent galaxies are

$$\begin{eqnarray}{B}_{H}({M}_{\star }) & = & 3.8{\rm{e}}-4{e}^{0.71{M}_{\star }}-0.11\quad (\mathrm{at}\ \mathrm{all}\ \mathrm{masses})\\ {\beta }_{H}({M}_{\star }) & = & 1.38{\rm{e}}12{e}^{-2.87{M}_{\star }}-1.21\quad ({M}_{\star }\geqslant 9.75)\\ {\beta }_{H}({M}_{\star }) & = & -0.19\qquad \qquad ({M}_{\star }\lt 9.75).\end{eqnarray} \tag{ 25 }$$

The resulting function ${R}_{\mathrm{eff},\mathrm{maj}}(z,{M}_{\star })$ for quiescent galaxies is shown in Figure 15. The significant size evolution among massive quiescent galaxies ($\mathrm{log}(M/{M}_{\odot })$ ≳ 10) is in agreement with other studies (e.g., Cassata et al. 2013). The flattening of the quiescent galaxy size–mass relation evolution at lower stellar mass is consistent with the expectations of environmental effects due to satellite quenching (e.g., van der Wel et al. 2010; Kawinwanichakij et al. 2017) and the observations that quiescent galaxies in high-density regions typically have larger sizes (e.g., Cooper et al. 2012; Delaye et al. 2014).

As we have outlined in Section 4, our model extrapolates the evolution of the stellar mass function for quiescent galaxies at z > 4, and therefore mock catalogs will contain samples of quiescent galaxies at redshifts beyond where current surveys can identify them or measure their morphologies. Although the expected number of quiescent galaxies at z > 4 in mock surveys would be small, we estimate that mock catalogs of comparable area to one GOODS field (∼150 arcmin²) will contain quiescent galaxies out to z ∼ 6. Therefore we note here that sizes for such objects come from an extrapolation of the relationship presented in Equation (25), which can be calculated for arbitrarily large redshift. We show the extrapolation out to z ∼ 6 in the right panel of Figure 15.

The data presented in van der Wel et al. (2014) indicate that the scatter in galaxy sizes within quiescent and star-forming galaxy samples is approximately uniform in both redshift and stellar mass. The size distributions within redshift and mass bins for quiescent galaxies are log-normal with σ = 0.16 dex, whereas the sizes of star-forming galaxies follow a skewed distribution better fit by a Gauss–Hermite polynomial expansion in ${\mathrm{log}}_{10}{R}_{\mathrm{eff},\mathrm{maj}}$ (van der Marel & Franx 1993) with fixed dispersion of 0.18 dex and skewness h3 = −0.15. The skewed distribution is motivated by the observation that the star-forming galaxy size distribution has a tail of small-sized galaxies (e.g., Williams et al. 2014, 2015; van Dokkum et al. 2015). The distribution has the form

$$\begin{eqnarray}P({R}_{\mathrm{mock}}| M* ,\,z) & \propto & \,{e}^{-\tfrac{1}{2}{y}^{2}(1+h3(\tfrac{1}{\sqrt{6}}(2\sqrt{2}{y}^{3}-3\sqrt{(}2)y)}\\ y & = & \displaystyle \frac{1}{\sigma }({\mathrm{log}}_{10}{R}_{\mathrm{mock}}-{\mathrm{log}}_{10}{R}_{\mathrm{eff},\mathrm{maj}}).\end{eqnarray} \tag{ 26 }$$

To assign sizes to mock galaxies of each type, we draw a size ${R}_{\mathrm{mock}}$ from random distributions of these forms, where ${R}_{\mathrm{eff},\mathrm{maj}}$ is the median size in kpc for each mock galaxy's redshift and stellar mass from Equations (24) and (25) (Figure 15).

At z > 3.5 the WFC3 H band moves out of the rest-frame optical. At higher redshifts, the current best imaging data set for morphological measures, the CANDELS H band imaging, will then become a probe of the rest-frame UV morphology of high-redshift galaxies. For mock star-forming galaxies at z > 4, we therefore assign morphologies according to the redshift evolution of the observed M_UV-size relationships presented in Shibuya et al. (2015). (As mentioned in the last section, quiescent galaxies at z > 4 follow the extrapolated relation presented in Equation (25).) The size evolution in that work is parametrized by circularized half-light radius, defined as

$$\begin{eqnarray}&&{R}_{\mathrm{eff},\mathrm{circ}}={R}_{\mathrm{eff},\mathrm{maj}}\sqrt{b/a},\end{eqnarray} \tag{ 27 }$$

where b/a is the axis ratio (ratio of the semiminor to semimajor axis size). Shibuya et al. (2015) find that median ${R}_{\mathrm{eff},\mathrm{circ}}$ of galaxies at z > 4 is correlated with UV luminosity, and at fixed UV luminosity, galaxy size significantly decreases with increasing redshift. Sizes at z > 4 in Shibuya et al. (2015) are measured from imaging at $1500\lt {\lambda }_{{\rm{rest}}}\lt 3000$ Å. We generate sizes according to their size parametrizations with UV luminosity:

$$\begin{eqnarray}&&{R}_{\mathrm{eff},\mathrm{circ}}={r}_{o}{\left(\displaystyle \frac{{L}_{\mathrm{UV}}}{{L}_{o}}\right)}^{0.27},\end{eqnarray} \tag{ 28 }$$

where ${r}_{o}=6.9{(1+z)}^{-1.2}$ represents the effective radius in kpc at a characteristic UV luminosity L_o (corresponding to ${M}_{{\rm{UV}}}=-21$). This relation is generally consistent with other measurements at z ∼ 6–7 (Oesch et al. 2010; Grazian et al. 2012; Ono et al. 2013; Holwerda et al. 2015; Kawamata et al. 2015; Curtis-Lake et al. 2016; Bowler et al. 2017). The observed distribution of sizes is log-normal at all redshifts, with a mean size according to the above parametrization and a constant scatter with redshift and luminosity, ${\sigma }_{\mathrm{ln}{R}_{\mathrm{eff},\mathrm{circ}}}=0.5$.

In addition to sizes, galaxy shapes are characterized by the following additional properties: projected axis ratio (semiminor half-light size/semimajor half-light size; b/a), Sérsic index n (setting the overall concentration of the light profile), and position angle on the sky. The distribution of observed axis ratios and Sérsic indices are known to correlate with star formation (e.g., Franx et al. 2008; Bell et al. 2012; Mortlock et al. 2013). Star-forming galaxies often exhibit extended, disk-like light profiles characterized by $n\lt 2$ and lower values of b/a, whereas quiescent galaxies exhibit centrally concentrated light profiles with $n\gt 2$ and higher b/a. These properties are known to evolve within each star formation sub-class, in that the n and b/a of quiescent galaxies tend to decrease with increasing redshift (e.g., van der Wel et al. 2011), while at high-redshift star-forming galaxies exhibit more clumpy concentrated morphologies in contrast to the extended disks seen at low redshift (e.g., Elmegreen et al. 2009; Förster Schreiber et al. 2011; Lee et al. 2013; Conselice 2014; Guo et al. 2015). To capture the differing evolution within each class of galaxy, we use the observed axis ratios and Sérsic indices that have been measured using GALFIT (Peng et al. 2002) from the deep HST/WFC3 F160W CANDELS imaging in both GOODS-North and GOODS-South (van der Wel et al. 2012). To produce the observed distribution of these parameters for each galaxy class and binned in redshift (Δz = 1), we match the real CANDELS galaxies in those morphology catalogs with robust morphology measurements ("flag = 0") to the photometric redshifts and rest-frame UVJ colors published by 3D-HST (Skelton et al. 2014). The resulting evolution with redshift of the distribution of parameters for star-forming and quiescent galaxies are shown in Figure 16.

Zoom In Zoom Out Reset image size

We assign axis ratios and Sérsic indices to mock star-forming and quiescent galaxies at a given redshift by generating random variates from those redshift-dependent distributions. We choose to draw directly from the binned distribution, rather than fit with an assumed functional shape for the evolution of the distributions. We first bin axis ratio measurements in bins of Δ(b/a) = 0.1 and Sérsic index in bins of Δn = 0.5, and treat these distributions as probability distribution functions. Then we draw random variates from these distributions to assign to mock galaxies, assuming uniform probability within each bin. Quiescent galaxies at z > 4 are assigned morphologies according to the $3\leqslant z\leqslant 4$ distributions, and star-forming galaxies at z > 6 follow the $5\leqslant z\leqslant 6$ distributions.

The resulting axis ratio for each mock galaxy is used to convert the semimajor axis at $z\leqslant 4$ into both semiminor axis and circularized half-light radius, and to convert circularized half-light radius into semimajor and minor axes at z > 4. Position angles for all galaxies are assigned randomly from a uniform distribution.

We compare the UV luminosity function at $z\leqslant 4$ computed from the mock catalog with measurements from the literature, as this enables us to test the adopted evolutionary model of both the star-forming galaxy stellar mass function (Section 3.1) and the ${M}_{{\rm{UV}}}$–${M}_{\star }$ relation (Section 3.2).

To compare with observations of the UV luminosity function, we use the compilation of literature measurements analyzed in Parsa et al. (2016) from $0.4\lt z\lt 4$. Parsa et al. (2016) provide a compilation of Schechter function parameters with error bars but do not quantify the covariance(s) between these parameters, which are known to be strong. Because of the degeneracies among the parameters, we avoid comparing directly to individually measured Schechter functions, and rather convert these into stepwise binned luminosity functions with equal magnitudes and bin widths. We then average the binned luminosity functions to produce a mean stepwise luminosity function at each redshift. We quantify the scatter in the literature as the standard deviation of the galaxy counts in each luminosity bin, divided by the square root of the number of measured luminosity functions contributing to the average. The binned averages and scatter from the literature are presented in Table 4.

Table 4.  Stepwise Average Binned Luminosity Functions from the Literature in the Redshift Bins Presented in Figure 17

${M}_{{\rm{UV}}}$	z ∼ 0.5	z ∼ 0.8	z ∼ 1.25	z ∼ 1.75	z ∼ 2.25	z ∼ 2.75	z ∼ 3.75
	Log Φa	Log Φ	Log Φ	Log Φ	Log Φ	Log Φ	Log Φ
−22.75	−14.59	−12.02 ± 1.34	−10.38 ± 0.88	−9.06 ± 0.21	−6.55 ± 0.70	−6.27 ± 0.30	−6.44 ± 0.17
−22.25	−10.42	−8.96 ± 1.10	−7.91 ± 0.69	−6.92 ± 0.17	−5.34 ± 0.53	−5.15 ± 0.24	−5.16 ± 0.11
−21.75	−7.74	−6.98 ± 0.85	−6.27 ± 0.50	−5.35 ± 0.15	−4.49 ± 0.38	−4.35 ± 0.19	−4.28 ± 0.07
−21.25	−6.02	−5.68 ± 0.62	−5.15 ± 0.34	−4.26 ± 0.15	−3.88 ± 0.27	−3.77 ± 0.17	−3.66 ± 0.05
−20.75	−4.88	−4.78 ± 0.43	−4.37 ± 0.21	−3.54 ± 0.16	−3.43 ± 0.18	−3.33 ± 0.17	−3.23 ± 0.05
−20.25	−4.13	−4.11 ± 0.27	−3.78 ± 0.11	−3.07 ± 0.17	−3.08 ± 0.12	−3.01 ± 0.17	−2.92 ± 0.05
−19.75	−3.61	−3.58 ± 0.16	−3.34 ± 0.04	−2.76 ± 0.17	−2.8 ± 0.09	−2.76 ± 0.17	−2.69 ± 0.05
−19.25	−3.25	−3.17 ± 0.08	−3.0 ± 0.03	−2.54 ± 0.16	−2.59 ± 0.10	−2.56 ± 0.17	−2.51 ± 0.06
−18.75	−2.98	−2.86 ± 0.03	−2.75 ± 0.05	−2.39 ± 0.14	−2.41 ± 0.12	−2.39 ± 0.17	−2.36 ± 0.07
−18.25	−2.76	−2.63 ± 0.01	−2.54 ± 0.06	−2.27 ± 0.12	−2.27 ± 0.15	−2.23 ± 0.17	−2.22 ± 0.08
−17.75	−2.59	−2.45 ± 0.01	−2.37 ± 0.07	−2.16 ± 0.10	−2.14 ± 0.17	−2.09 ± 0.18	−2.10 ± 0.10
−17.25	−2.44	−2.3 ± 0.00	−2.23 ± 0.09	−2.06 ± 0.08	−2.03 ± 0.21	−1.95 ± 0.20	−1.99 ± 0.12

Note.

^aThe z ∼ 0.5 binned luminosity function only includes one measurement from the literature; therefore the observed scatter in that bin is zero.

Download table as:  ASCIITypeset image

Figure 17 shows the comparison of the $0.4\lt z\lt 4$ mock galaxy UV luminosity functions (black circles) with the average literature measurements (blue squares).²⁶ At z > 1.5, the mock catalog exhibits excellent agreement with the observations, while at z < 1.5 the agreement is less robust. We note, however, that the mock catalog never overpredicts the measured number density of galaxies by more than ∼0.36 dex (at ${M}_{{\rm{UV}}}\sim -19$ at $0.6\lt z\lt 1$ and $1\lt z\lt 1.5$). This overprediction is likely caused by the poor observational constraints on the rest-frame UV emission of galaxies at z < 1.5, where photometry in the 3D-HST catalog no longer provides coverage at 1500 Å (rest wavelength). This affects both the characterization of the ${M}_{{\rm{UV}}}$–${M}_{\star }$ relation, which we use to directly assign ${M}_{{\rm{UV}}}$ values to low-mass mock galaxies at these redshifts, as well as the ${M}_{{\rm{UV}}}$ values of high-mass galaxies that are derived directly from fits to the 3D-HST photometry. The same uncertainties affect the observed UV luminosity functions compiled from the literature between z ∼ 0, where GALEX data are available, and z ∼ 1, where HST near-UV bands probe the rest-frame 1500 Å flux (e.g., Oesch et al. 2010; Windhorst et al. 2011; Teplitz et al. 2013; Rafelski et al. 2015).

Zoom In Zoom Out Reset image size

In Figure 18, we compare the CSFRD of mock galaxies with the average CSFRD evolution presented in Madau and Dickinson (2014, based on a uniform analysis of luminosity function measurements in the literature), converted to a Chabrier IMF. To approximate the same UV luminosity limits imposed by Madau & Dickinson (2014), we estimate the limiting ${M}_{{\rm{UV}}}$ corresponding to $0.03\,{L}^{* }$ in each redshift bin. These values are taken directly from the Schechter UV luminosity function parameter ${M}_{{\rm{}}{UV}}^{* }$, which is used to estimate the equivalent limit in ${M}_{{\rm{UV}}}$ from $0\lt z\lt 10$. The equivalent of $0.03\,{L}^{* }$, as published in the following studies, corresponds to ${M}_{{\rm{UV}}}\sim -14.5$ at $1\lt z$ (Cucciati et al. 2012); ${M}_{{\rm{UV}}}\sim -15.5$ at $1\lt z\lt 2$ (Cucciati et al. 2012); ${M}_{{\rm{UV}}}\sim -16.89$ at $1.9\lt z\lt 2.7$ (Reddy & Steidel 2009); and ${M}_{{\rm{UV}}}\sim -17$ at z > 2.7 (Bouwens et al. 2015, 2016b; Finkelstein et al. 2015; Oesch et al. 2017).

Zoom In Zoom Out Reset image size

To calculate the CSFRD of mock galaxies in each redshift bin, we identify all mock galaxies above the limiting ${M}_{{\rm{UV}}}$ corresponding to $0.03\,{L}^{* }$ and take the sum of their star formation rates (averaged over the past 100 Myr) per co-moving volume. We find that the evolution of the CSFRD estimated from the mock catalog qualitatively reproduces the overall shape of the observed CSFRD; however, at $4\lt z\lt 6$ mock galaxies show a slight excess with respect to the results of Madau & Dickinson (2014), who derive their z > 4 points exclusively from the UV luminosity functions of Bouwens et al. (2012a, 2012b). Bouwens et al. (2012a) provided measurements of the z ∼ 6 luminosity function, along with a compilation of earlier measurements from at z ∼ 4–5 and z ∼ 7–8 (Bouwens et al. 2007, 2011). However, our model is based on the UV luminosity function of Bouwens et al. (2015), who find higher number counts than in previous works, in particular at the bright end at ${M}_{{\rm{UV}}}\lt -18$ and z ∼  4–5, and at ${M}_{{\rm{UV}}}\lt -20$ and z > 6, which likely explains our excess relative to Madau & Dickinson (2014). At $1\lt z\lt 3$, the CSFRD measured from mock galaxies is slightly lower than the one of Madau & Dickinson (2014). A reason for this may be that our phenomenological model does not explicitly include extremely dust-obscured galaxies, which may be missing from the HST-selected samples that we use to characterize the stellar mass function evolution, whereas Madau & Dickinson (2014) do incorporate far-infrared measurements of dust-obscured objects. To validate this explanation, we include individual measurements of the CSFRD from UV-selected samples at z < 3 and their uncertainties (Schiminovich et al. 2005; Reddy & Steidel 2009; Cucciati et al. 2012). These indicate a better agreement with our mock galaxies, within the scatter of binned observations, than implied by the Madau & Dickinson (2014) curve.

We additionally plot measurements of the CSFRD at z > 4 from Bouwens et al. (2015, 2016a), Finkelstein et al. (2015), McLeod et al. (2016), and Oesch et al. ( 2017), converted to Chabrier IMF when necessary. The CSFRD evolution of mock galaxies at z > 4 is in better agreement with these later measurements. The ${M}_{{\rm{UV}}}$ integration limits in these studies match those used to estimate the CSFRD from the mock catalog, with the exception of McLeod et al. (2016), whose integration limit is ${M}_{{\rm{UV}}}\sim -17.7$. The CSFRD of McLeod et al. (2016) would be higher if measured to the same limiting ${M}_{{\rm{UV}}}$ used for the mock catalog. The higher UV counts (and thus CSFRD) measured in McLeod et al. (2016) have been interpreted as evidence for a slower rate of decline in galaxy number counts with redshift, in contrast to the "accelerated" evolution that is seen in, for example, Bouwens et al. (2015) and Oesch et al. (2017). Since our model follows the more rapidly evolving luminosity function of Oesch et al. (2017), the CSFRD in the mock follows more closely the CSFRD measured in that work. We find that our mock realization generally reproduces the z > 4 evolution of the CSFRD based on the latest results from current extragalactic surveys, and in particular those we have used to develop our model for the evolution of the stellar mass function.

We compare the redshift evolution of the sSFR from the mock catalog with literature measurements at $0\lt z\lt 7$, the highest redshift at which observational constraints are available. We consider the median sSFR of galaxies with masses $8.8\lt \mathrm{log}(M/{M}_{\odot })\lt 10$, as computed by Noeske et al. (2007), Damen et al. (2009), Daddi et al. (2007), Reddy & Steidel (2009), Stark et al. (2013), González et al. (2014), Tasca et al. (2015), and Salmon et al. (2015). To calculate the median sSFR of mock galaxies, we use the beagle-assigned SFR averaged over the last 100 Myr divided by stellar mass of mock galaxies in the same stellar mass range as the observations. We plot this comparison in Figure 19, and show that the sSFR values for our mock are in excellent agreement with the observations within their uncertainties.

Zoom In Zoom Out Reset image size

As discussed in Section 3.2, the sSFR of galaxies at z ≳ 4 and the redshift dependence of the average sSFR are directly related to the extrapolation we adopted for the redshift evolution and normalization of the ${M}_{{\rm{UV}}}$–${M}_{\star }$ relation. This is in contrast to mock galaxy sSFRs at z < 4, where the ${M}_{{\rm{UV}}}$–${M}_{\star }$ relation is based directly on more reliable stellar mass measurements across a range of UV luminosities. Therefore the agreement between the mock galaxy sSFRs with measurements in the literature at z > 4 shown in Figure 19 is a validation of the adopted extrapolations of the ${M}_{{\rm{UV}}}$–${M}_{\star }$ relation at high redshift. We note, however, that the evolution of the sSFR of galaxies at z > 4 is still a matter of active research; most studies support an increasing average sSFR rather than a plateau at z > 4 (see the discussion in Stark 2016). An accurate characterization of the sSFR evolution at z > 4 will await improved measurements that will be made possible by JWST.

The mass–metallicity relation of galaxies is known to evolve, such that galaxies at higher redshifts are more metal-poor at a given stellar mass (e.g., Maiolino et al. 2008). When assigning SEDs (and associated physical properties) to each mock galaxy, we use a parent catalog that has a built-in redshift-independent prior linking ${M}_{\star }\mbox{--}Z\mbox{--}\psi$, based on the fundamental metallicity relation as measured by Hunt et al. (2016). The redshift-evolving ${M}_{{\rm{UV}}}$–${M}_{\star }$ relation, however, produces a redshift evolution of the ${M}_{\star }$−$\psi$ relation that, in turn, should generate an evolving mass–metallicity relation.

Figure 20 shows the mass–metallicity relation in two redshift bins: $0.2\leqslant z\leqslant 0.5$ and $2\leqslant z\leqslant 2.5$. Given that different metallicity indicators provide systematically different metallicity estimates (Kewley & Ellison 2008), we choose to estimate the mock galaxy metallicities using a metallicity calibration often used in the literature, based on the ratio between the lines $[{\rm{N}}\,{\rm{II}}]\lambda 6584$ and ${\rm{H}}\alpha$ (the Pettini & Pagel 2004 N2 metallicity calibration):

$$\begin{eqnarray}&&12+\mathrm{log}{({\rm{O}}/{\rm{H}})}_{{\rm{N}}2}=8.90+0.57\times \mathrm{log}([{\rm{N}}\,{\rm{II}}]\lambda 6584/{\rm{H}}\alpha ).\end{eqnarray} \tag{ 29 }$$

This allows us to compare directly to observationally derived mass–metallicity relations.

Zoom In Zoom Out Reset image size

Figure 20 indicates that the mass–metallicity relation derived from the mock catalog does show a turnover at high stellar masses, as seen in the observations, despite the Hunt et al. (2016) fundamental metallicity relation only being linearly dependent on $\mathrm{log}(Z/{Z}_{\odot })$. The turnover is caused by the presence of a maximum metallicity $\mathrm{log}(Z/{Z}_{\odot })=0.24$ that can be assigned to the mock galaxies (see Sections 3.4.2 and 3.4.3). Figure 20 demonstrates a good agreement among our mock galaxies and the mass–metallicity relation of Kewley & Ellison (2008; at z ∼ 0.07) and of Sanders et al. (2015; at z ∼ 2.3). As discussed above, the evolution of the mass–metallicity relation in our mock catalog is a result of the priors we impose on both ${M}_{\star }\mbox{--}Z\mbox{--}\psi$ and ${M}_{{\rm{UV}}}$–${M}_{\star }$.

Emission line diagnostic diagrams are commonly used in the literature to identify the ionization sources in galaxies and characterize their physical properties. Baldwin et al. (1981) (BPT hereafter) pioneered the use of the ratios of $[{\rm{O}}\,{\rm{III}}]\lambda 5007/{\rm{H}}\beta$ and $[{\rm{N}}\,{\rm{II}}]\lambda 6584/{\rm{H}}\alpha$ to separate star-forming galaxies and AGN-dominated galaxies into two distinct regions of the plot, with composite galaxies (those with both significant star formation and AGN contribution) spanning the region in between. Recent observations of rest-frame optical spectra at redshifts out to z ∼ 2.5 have revealed that the locus of star-forming galaxies is slightly shifted to higher $[{\rm{O}}\,{\rm{III}}]\lambda 5007/{\rm{H}}\beta$ ratios, at fixed $[{\rm{N}}\,{\rm{II}}]\lambda 6584/{\rm{H}}\alpha$, than in the local universe (e.g., Shapley et al. 2005; Brinchmann et al. 2008; Hainline et al. 2009; Kewley et al. 2013; Steidel et al. 2014; Kashino et al. 2017). This shift can been explained if the physical conditions of the ionized gas are evolving with redshift. For example, evolution in the ionization parameter (e.g., Kashino et al. 2017), the incident ionizing spectrum of stars (e.g., Steidel et al. 2014, 2016), gas-phase nitrogen abundance (e.g., Masters et al. 2016), and hydrogen density (e.g., Sanders et al. 2016) can all introduce a shift in observed $[{\rm{O}}\,{\rm{III}}]\lambda 5007/{\rm{H}}\beta$ ratio. However, it is also possible that selection effects play a role (e.g., Juneau et al. 2014). More likely, the observed evolution of the location of star-forming galaxies in the BPT diagram is caused by a combination of the effects above, but current data do not allow us to quantify their respective roles.

In Figure 21 we plot the BPT diagram of mock galaxies at low ($0.2\lt z\lt 0.5$) and high ($2\lt z\lt 2.5$) redshift. We compare the mock galaxies with SDSS galaxies (z ∼ 0.07) and with a sample of star-forming galaxies from the KBSS survey (Steidel et al. 2014) (z ∼ 2.5). For mock galaxies at $0.2\lt z\lt 0.5$, we adopt the lowest luminosity threshold of Juneau et al. (2014), and only plot objects with log luminosities $\gt 39.9$ erg s⁻¹ in the lines $[{\rm{O}}\,{\rm{III}}]\lambda 5007$, ${\rm{H}}\beta$, ${\rm{H}}\alpha$, $[{\rm{N}}\,{\rm{II}}]\lambda 6584$. As Figure 21 shows, the mock galaxies sit close to the main star-forming galaxy locus, despite being at higher redshift than the SDSS sample. At $2\lt z\lt 2.5$, we instead apply a flux limit close to that of Steidel et al. (2014; i.e., line flux $\gt 3\times {10}^{-18}$ erg s⁻¹ cm⁻²). At these higher redshifts, the mock galaxies cover similar regions to the Steidel et al. (2014) points at low $[{\rm{N}}\,{\rm{II}}]\lambda 6584/{\rm{H}}\alpha$, but not at high values of $[{\rm{N}}\,{\rm{II}}]\lambda 6584/{\rm{H}}\alpha$. This difference can be understood by appealing to the recent study of Hirschmann et al. (2017), where they self-consistently couple the nebular emission models of Gutkin et al. (2016) for star-forming galaxies and narrow-line AGN-driven models of Feltre et al. (2016) to cosmological zoom-in hydrodynamic simulations. In their work, they consistently tie $\mathrm{log}{U}_{{\rm{S}}}$ to the simulated galaxy properties, and fix ${\xi }_{{\rm{d}}}$ and ${n}_{{\rm{H}}}$ to the same values as fixed in this mock catalog. They reproduce the observed evolution of the line ratios (i.e., higher $[{\rm{O}}\,{\rm{III}}]\lambda 5007/{\rm{H}}\beta$ on average at higher redshifts), but their Figure 12 indicates that covering the elevated $[{\rm{O}}\,{\rm{III}}]\lambda 5007/{\rm{H}}\beta$ values at high $[{\rm{N}}\,{\rm{II}}]\lambda 6584/{\rm{H}}\alpha$ requires either a higher value of ${\xi }_{{\rm{d}}}$ and/or of ${n}_{{\rm{H}}}$, or some level of AGN activity. We add to the the right panel of Figure 21 two vectors indicating how variations of these two parameters modify the expected line ratios. Given that the possible reasons for the elevated $[{\rm{O}}\,{\rm{III}}]\lambda 5007/{\rm{H}}\beta$ ratio is far from determined, we choose to supply line ratios for mock galaxies with different values of ${\xi }_{{\rm{d}}}$ and ${n}_{{\rm{H}}}$. Oxygen is a refactory element, unlike nitrogen, and is thus heavily depleted onto dust grains. An elevated ${\xi }_{{\rm{d}}}$ ratio therefore acts to increase the relative abundance of gas-phase nitrogen to oxygen. Whether we expect an elevated nitrogen abundance to be due to dust depletion or the relative importance of primary versus secondary nitrogen production in stars at different epochs, the effect on the O iii Hβ and N ii Hα line ratios will be the same. We therefore supply a single mock realization with ξ_d = 0.5 (where the fiducial mock has ξ_d = 0.3). We also provide line fluxes and equivalent for galaxies in that realization with a higher hydrogen density. The Gutkin et al. (2016) models are evaluated in unit steps in $\mathrm{log}({n}_{{\rm{H}}})$, so we are able to supply this line information for ${n}_{{\rm{H}}}=1000$ cm⁻³ (where the fiducial mock has ${n}_{{\rm{H}}}\,=100$ cm⁻³).

Zoom In Zoom Out Reset image size

Strong nebular emission lines are now known to contaminate the broad-band photometry of galaxies. In particular, colors from Spitzer/IRAC have been extensively used to infer the EW of strong optical emission lines (Shim et al. 2011; Smit et al. 2014, 2015, 2016; Faisst et al. 2016; Mármol-Queraltó et al. 2016; Rasappu et al. 2016; Castellano et al. 2017). In Figure 22, we compare the IRAC 3.6–4.5 μm color of mock galaxies with measurements at z > 4 from Smit et al. (2014, 2015, 2016) and Rasappu et al. (2016). The redshift evolution of the 3.6–4.5 μm color predicted by the mock catalog, caused by the emission lines ${\rm{H}}\alpha$, $[{\rm{O}}\,{\rm{III}}]$, and ${\rm{H}}\beta$ entering and leaving the IRAC bands, is consistent with the observations. This suggests that the emission line strengths in the mock galaxies are compatible with the values inferred from IRAC observations. In Section 7.2 below, we will use the mock catalog to show how the sensitivity and wavelength coverage of JWST/NIRCam will constrain emission line EWs with greater accuracy than existing Spitzer data.

Zoom In Zoom Out Reset image size

In Figures 25 and 26, we show predictions for the high-redshift detections with NIRCam imaging in the JADES GTO survey. In these figures, we plot only those objects that are detected with at least $5\sigma$ (assuming point source detection limits) in two rest-frame UV photometric bands: that closest to 1500 Å, and the nearest band at a longer wavelength. The detection bands correspond to F115W and F150W at $6\lt z\lt 7;$ F150W and F200W at $7\lt z\lt 9.6;$ and F200W and F277W at $9.6\lt z\lt 13$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 25 shows the expected number of objects detected as a function of ${M}_{{\rm{UV}}}$ and redshift in the JADES survey based on our phenomenological model. We find that JADES should detect several thousands of galaxies at z > 6, and several 10s out to z > 10. These predictions are based on the galaxy number counts evolution of our phenomenological model, which follows observations supporting a more rapid ("accelerated") evolution of the UV luminosity function at z ≳ 8 (e.g., Bouwens et al. 2016b; Oesch et al. 2017), implying a factor of ∼10 decrease of galaxy counts between z ∼ 8 and z ∼ 10. This is substantially faster than the factor of ∼2 decrease of galaxy counts per unit redshift seen at $3\lt z\lt 8$ (e.g., Bouwens et al. 2015; Finkelstein 2016), which are similar to the decrease measured by McLeod et al. (2015, 2016) at z > 8. However, due to large uncertainties from small sample sizes and survey volumes at z ≳ 8 in current observations, there remain discrepancies in the literature of how fast the luminosity function evolves (Oesch et al. 2012, 2014, 2017; Zheng et al. 2012; Coe et al. 2013; Ellis et al. 2013; Bouwens et al. 2015, 2016b; McLeod et al. 2015, 2016; Calvi et al. 2016; Ishigaki et al. 2017; Stefanon et al. 2017b). Accurately measuring the evolution at z > 8 is an important goal for JWST because it will directly constrain the relationship between star formation efficiency and the evolution of the halo mass function at early times (Trenti et al. 2010; Tacchella et al. 2013; Mason et al. 2015; Liu et al. 2016; Mashian et al. 2016; Sun & Furlanetto 2016; Ceverino et al. 2017; Cowley et al. 2018), and is critical for our understanding of reionization.

Therefore, we have used our mock catalog tool to predict if JADES will distinguish between the constant versus "accelerated" evolutionary model at z ≳ 8. To do this, we have used the tool to produce a comparison mock realization based on the constant rate of evolution that was parametrized in Bouwens et al. (2015) at $4\lt z\lt 8$, which we extrapolate to z ≳ 8. This parametrization is in good agreement with the observations at ${\text{}}z\sim 9\mbox{--}10$ by McLeod et al. (2015, 2016). For this secondary set of realizations, we follow the same procedures for assigning physical properties and SEDs to mock galaxies as outlined previously in this work; the only difference is the surface density of counts of a given ${M}_{{\rm{UV}}}$ predicted by the UV luminosity functions. In Figure 26, we additionally plot the total predicted number of galaxies that will be detected per redshift bin, expected for the constant rate of evolution from Bouwens et al. (2015) from the JADES survey (red points) in addition to the accelerated model (black points). The error bars in this figure represent both Poisson uncertainties as well as estimates for the uncertainty due to cosmic variance (Trenti & Stiavelli 2008), assuming the JADES survey area and a Press-Schechter formalism. From this figure it is clear that the future JADES survey will discriminate between the two evolutionary scenarios.

The inclusion of nebular emission lines in the SEDs of mock galaxies enables analyses of the photometric contamination that is likely to be observed in future NIRCam imaging. Identifying the exact level of contamination in current surveys remains difficult, because there are only two deep Spitzer bands that provide rest-frame optical coverage of SEDs above z ∼ 4.5. These bands are broad and there are only limited, narrow redshift ranges where one or other band is free from emission line contamination, and therefore able to provide constraints on the stellar continuum. However, the extensive filter set of NIRCam, with three broad-band filters red-wards of current ground-based K band imaging as well as a large array of medium- and narrow-band filters, will offer the opportunity to better characterize both the typical and extreme emission line EWs, which in turn will provide much better constraints on stellar masses at z > 4. As summarized in Table 5, the JADES survey will provide imaging in two medium bands at λ > 3 μm (F335M and F410M). These filters provide continuum anchor points when the emission line is outside the medium-band filter, or direct emission line contamination estimates of the broad-band fluxes when the line sits within the filter.

In Figure 27, we plot the predicted redshift evolution of the NIRCam colors F335M–F356W and F410M–F444W of mock galaxies that would be detected at $\gt 5\sigma$ in the deep region of the JADES survey. The galaxies are color-coded by the rest-frame EW of the emission lines that are contributing to the observed photometric flux at that redshift: Hα at z < 5.2 and $[{\rm{O}}\,{\rm{III}}]\lambda 5007$ at z > 5.2 ($[{\rm{O}}\,{\rm{III}}]\lambda 5007$ dominates the EW of the $[{\rm{O}}\,{\rm{III}}]\lambda \lambda 4959,5007$ doublet). This figure demonstrates how strong emission lines can impact the observed infrared colors of galaxies in specific redshift ranges. The multiple overlapping NIRCam bands means that not all bands are contaminated at a given time, and stellar mass estimates are likely to be better constrained than currently possible at z ∼ 6, where both IRAC bands are affected by emission line contamination. Additionally, Figure 27 demonstrates how NIRCam colors can be used to select high-EW emission line objects at specific redshifts. Of those galaxies in the mock catalog with 5σ Deep JADES detections, a color cut of F335M–F356W >0.8 and <−0.5 will be successful at selecting a sample that is 99% composed of objects with rest-frame EW${}_{{\rm{H}}\alpha }\gtrsim 1000$ Å at z < 5.2 (the number density of this sample is 3.85 objects per square arcminute). Similarly, a color cut of F335M–F356W >0.8 and <−0.6 will be successful at selecting a sample that is 98% composed of objects with rest-frame EW${}_{[{\rm{O}}{\rm{III}}]\lambda 5007}\gtrsim 1000$ Å at z > 5.2 (the number density of this sample is 1.62 objects per square arcminute). Objects selected with such extreme EW emission lines would be obvious targets for follow-up spectroscopy with NIRSpec, and provide an independent probe of the Hα-derived star formation rates free from the effects of slit losses produced by the MSA on NIRSpec. They would additionally provide a consistent probe of the evolution of the number densities of extreme emission line galaxies from z ∼ 4.8–8.3.

Zoom In Zoom Out Reset image size

As described in Section 3.1, we re-fit the individual Bouwens et al. (2015) UV luminosity functions with mass function Schechter parameters. To do this we convolve a given mass function with our model of the evolving ${M}_{{\rm{UV}}}$–${M}_{\star }$ relation (described in Section 3.2) to produce the corresponding UV luminosity function (general procedure described in Section 2.1). We then use Markov-Chain Monte-Carlo (MCMC) sampling, employing the adaptive metropolis algorithm of Haario et al. (2001), to sample from the posterior probability distributions of the Schechter function parameters. At each iteration, Schechter function parameters are proposed and corresponding CDF function generated. The CDF function is used to randomly assign masses to a large population of objects, which can then be used to assign each object an ${M}_{{\rm{UV}}}$ value following our adopted redshift-evolving ${M}_{{\rm{UV}}}$–${M}_{\star }$ model. The resulting UV luminosity function is measured within the ${M}_{{\rm{UV}}}$ bins of the published measurements. The likelihood of those given Schechter function parameters is evaluated using the published ${\rm{\Phi }}({M}_{{\rm{UV}}},z)$ values and associated errors given in Bouwens et al. (2015), Table 5, and the modeled luminosity function values evaluated in the same ${M}_{{\rm{UV}}}$ bins.

The z ≳ 4 luminosity functions are well described by a single Schechter function; thus fitting these data with a double-Schechter function renders the parameters of ${{\rm{\Phi }}}_{1}({M}_{\star })$ unconstrained at z ≳ 4. We therefore constrain the values of each ${{\rm{\Phi }}}_{1}({M}_{\star })$ parameter while performing the luminosity function fits. The values of ${\alpha }_{1,{\rm{M}}}$ and ${\phi }_{1,{\rm{M}}}^{* }$ are fixed in each redshift bin using the extrapolated best-fit linear and quadratic relations to the published T14 maximum-likelihood measurements, and we also tie the value of ${M}_{1,{\rm{M}}}^{* }$ to that of ${M}_{2,{\rm{M}}}^{* }$ with ${M}_{1,{\rm{M}}}^{* }={M}_{2,{\rm{M}}}^{* }={M}_{{\rm{M}}}^{* }$. The exact form of the redshift evolution of the parameters of ${{\rm{\Phi }}}_{1}({M}_{\star })$ do not affect the results of the luminosity function fits, as long as the number density has fallen significantly by z ∼ 4.

The measurements are displayed in Figure 3.

To constrain the parameters a₁, b₁, b₂, b₃, c₁, c₂, e₁, and f₁ of our model of the redshift evolution of the mass function (Equations (5)–(10)), we fit to the T14 star-forming mass functions using a Bayesian multi-level modeling approach. Whereas the Bayesian fitting to each individual luminosity function entails sampling from the posterior distribution of each Schechter function parameter, the multi-level modeling involves sampling over the posterior distribution of the hyper-parameters (as we shall now refer to the parameters a₁, b₁, b₂, b₃, c₁, c₂, e₁, and f₁) and the conditional probability distributions for each Schechter function parameter in each redshift bin. The posterior distribution of the model free parameters can be expressed as

$$\begin{eqnarray}&&P({\boldsymbol{A}},{\boldsymbol{\Phi }}| {\boldsymbol{X}})\propto \displaystyle \prod _{i=1}^{8}P({A}_{i})\,\displaystyle \prod _{z=1}^{8}P({{\boldsymbol{x}}}_{z}| {{\boldsymbol{\phi }}}_{z})\,\displaystyle \prod _{z=1}^{8}P({{\boldsymbol{\phi }}}_{z}| {\boldsymbol{A}}),\end{eqnarray} \tag{ 30 }$$

where ${\boldsymbol{A}}=[{a}_{1},{b}_{1},{b}_{2},{b}_{3},{c}_{1},{c}_{2},{e}_{1},{f}_{1}]$ represents the free parameters describing the redshift evolution of the individual Schechter function parameters, $P({A}_{i})$ is the prior on the ith free parameter, ${\rm{\Phi }}=[{\phi }_{1},{\phi }_{2},\,\ldots ,\,{\phi }_{n}]$ represents the set of n mass function parameters, ${\boldsymbol{X}}=[{{\boldsymbol{x}}}_{0},{{\boldsymbol{x}}}_{1},\,\ldots ,\,{{\boldsymbol{x}}}_{8}]$ represents the measurements of the mass function in each redshift bin, and ${{\boldsymbol{\phi }}}_{z}=[{M}_{1,{\rm{M}},z}^{* },{\phi }_{1,{\rm{M}},z}^{* },{\alpha }_{1,{\rm{M}},z},{M}_{2,{\rm{M}},z}^{* },{\phi }_{2,{\rm{M}},z}^{* },{\alpha }_{2,{\rm{M}},z}]$ represents the Schechter function parameters in each redshift bin. We fit with weakly informative priors (broad Gaussian distributions) on each hyper-parameter, with the added constraint that b₃ be negative to ensure that the normalization of ${{\rm{\Phi }}}_{1}(M)$ decreases with redshift, and e₁ and f₁ be positive. All priors are reported in Table 4. We assume that the measurements represent true measurements of the mass function at the mid-point of each redshift bin (z = [0.35, 0.625, 0.875, 1.125, 1.375, 1.75, 2.25, 2.75]), and so do not include any modeling of the redshift distribution of the underlying sources that make up the mass bin volume density measurements.

The mass function estimates reported in T14 are supplied as $\mathrm{log}[{\rm{\Phi }}({M}_{\star })]$ (see their Table 1), with no information of the number of galaxies entering each bin. In fact, the associated errors adopted from T14 account for Poisson noise, cosmic variance, and the uncertainties arising from classifying galaxies as star-forming or quiescent and from the determination of stellar masses; thus they are not simply Poissonian errors. Therefore, in the absence of the information required to construct the correct function for our count distribution, we resort to a Gaussian assumption, giving:

$$\begin{eqnarray}&&P({{\boldsymbol{x}}}_{z}| {\phi }_{z})=\displaystyle \prod _{i}\,\displaystyle \frac{1}{\sqrt{2\pi }{\sigma }_{i}}\,\exp \left(-\displaystyle \frac{{({{\rm{x}}}_{i}-\hat{{x}_{i}}({{\boldsymbol{\phi }}}_{z}))}^{2}}{2{\sigma }_{i}^{2}}\right),\end{eqnarray} \tag{ 31 }$$

where x_i and ${\sigma }_{i}$ are the estimate and error of the mass function in bin i, respectively, and $\hat{{x}_{i}}({{\boldsymbol{\phi }}}_{z})$ is the model prediction.

We sample the posterior distribution of model parameters with a Metropolis-within-Gibbs sampler (see, e.g., Sharma 2017). Gibbs sampling involves drawing directly from the conditional distribution for each parameter in turn. However, when the conditional probability is tricky to derive, it is possible to use the Metropolis update step, where the next step in the chain is sampled from a proposal distribution, often a Gaussian, and accepted or rejected based on comparison of the posterior probability between the current and last steps of the chain. The advantage of using a Gibbs sampler for this problem is that the conditional distributions of the model parameters in each redshift bin are independent of each other, meaning that $P({{\boldsymbol{\phi }}}_{i}| \,...\,{{\boldsymbol{\phi }}}_{i-1},{{\boldsymbol{\phi }}}_{i+1},\,\ldots ,\,{\boldsymbol{A}},{\boldsymbol{X}})=P({{\boldsymbol{x}}}_{i}| {{\boldsymbol{\phi }}}_{i})\,P({{\boldsymbol{\phi }}}_{i}| {\boldsymbol{A}})$, and so the mass function parameters for each redshift bin can be sampled in turn, before updating the values of the hyper-parameters.

During each iteration of the chain, our multi-level modeling algorithm performs the following steps:

• 1.  
  Values for the double-Schechter function parameters are proposed for each individual redshift bin in turn, sampling from $P({{\boldsymbol{\phi }}}_{i}| \,...\,{{\boldsymbol{\phi }}}_{i-1},{{\boldsymbol{\phi }}}_{i+1},\,\ldots ,\,{\boldsymbol{A}},{\boldsymbol{X}})$, using a Metropolis update step.
• 2.  
  Once the mass function parameters of each redshift bin have been updated, new values for the hyper-parameters are proposed. They are updated simultaneously, sampling from $P({\boldsymbol{A}}| {\boldsymbol{\Phi }},{\boldsymbol{X}})$, again using a Metropolis update step.

After an initial burn-in stage of 10,000 iterations, the width of the Gaussian proposal distribution is set using the covariance matrix of the past history of the chains, following the adaptive Metropolis algorithm of Haario et al. (2001). We run two different chains with independent starting points for 40,000 iterations each, after which we test for convergence using the scale reduction factor, $\hat{R}$, which compares the within-chain and between-chain variance. We require $\hat{R}\lt 1.05$ in each free parameter.

After rejecting the initial 10,000 iterations, we estimate each parameter value as the median of the values in the chain and uncertainties as their standard deviation.