On the (Lack of) Evolution of the Stellar Mass Function of Massive Galaxies from z = 1.5 to 0.4

For this study, we use a new K-band-selected catalog for SHELA. This catalog is based on imaging from the NEWFIRM-HETDEX Survey (Stevans et al. 2020), which covers 17.5 deg² to a median 5σ depth of K_s = 22.7 AB mag (2'' diameter apertures). The advantage of using the K-band catalog is that for our galaxies of interest, this filter samples the rest-frame near-IR (∼1 μm) and is therefore more sensitive to the older stellar populations which dominate the light of many systems. We have also verified that including the K-band photometry improves the quality of our photometric redshifts and stellar population parameter fits (including stellar masses). A full description of the NEWFIRM K-band imaging, catalog construction, and derived properties (photometric redshifts and stellar masses) is provided by Stevans et al. (2020). For reference, we use the NEWFIRM K-band photometry measured in the FLUX_AUTO apertures from SExtractor (this is relevant for the next subsection).

In addition to the NEWFIRM K-band data, the SHELA survey includes ugriz imaging from the Dark Energy Camera over 17.5 deg² (DECam; Wold et al. 2019) and deep 3.6 and 4.5 μm imaging from Spitzer/IRAC (Papovich et al. 2016). The riz-band-selected DECam catalogs reach a 5σ depth of ∼24.5 AB mag. However, we do not use the DECam catalogs and the DECam-selected Spitzer/IRAC forced photometry (Wold et al. 2019). In our analysis, we use our new K-band image for detection for the reason we have described in the previous section. Following our work in Kawinwanichakij (2018) and Wold et al. (2019), we perform "forced photometry" to derive optimal fluxes in the ugriz + [3.6] and [4.5] data for sources detected in the K_s catalog. We use the code "The Tractor" (Lang et al. 2016a, 2016b) for this process. This allows us to measure flux densities (or accurate limits) even for sources below the 5σ depth threshold of the original DECam or IRAC imaging or for sources blended at the resolution of Spitzer/IRAC. We follow identical procedures to those described in Kawinwanichakij (2018) and Wold et al. (2019), except we use the NEWFIRM K-band image for detection. With this technique, our forced photometric IRAC catalog is 80% complete to a limiting magnitude of 22.6 AB mag in both 3.6 and 4.5 μm bands, in contrast to 22.3 AB mag for the original published Spitzer/IRAC-selected catalog (Papovich et al. 2016).

For all of the analyses here, we estimate the stellar mass, ${M}_{* ,\mathrm{iSED}}$, using the Bayesian iSEDfit¹¹ package (Moustakas 2017) presented in Moustakas et al. (2013). The iSEDfit code performs a refined grid search of the posterior distributions of stellar mass and enables priors with nonflat probability distributions. The advantage of using iSEDfit is that it allows the use of different assumptions in the SFHs (including "bursts" of star formation) and in the SPS models, including those of the Flexible Stellar Population Synthesis¹² (FSPS; Conroy & Gunn 2010a, 2010b), Bruzual & Charlot (2003), and Maraston (2005) models. In this way, we are able to test for systematic uncertainties in the stellar masses resulting from differences in the underlying assumptions, without systematic uncertainties resulting from different stellar-population fitting codes. We summarize the key aspects of the iSEDfit code below.

For iSEDfit, we adopt fiducial prior parameters for ${M}_{* ,\mathrm{iSED}}$ from Moustakas et al. (2013). The basic set of iSEDfit priors is based only on a set (randomly generated for each run of iSEDfit) of SFHs using 10,000 declining exponential models, where star formation rate (SFR) $\propto \,\exp (-t/\tau )$, for age t and e-folding timescale τ. The parameters for each iSEDfit model vary independently, and it is therefore important to sample the entire the range of each prior. We also allow the iSEDfit age t (time since the onset of star formation) of each model to have uniform probability from 0.1 to 13 Gyr; however, we disallow ages older than the age of the universe at the redshift of each galaxy. We draw the e-folding time (τ) from the linear range 0.1–5 Gyr. We assume a uniform prior on stellar metallicity, Z, in the range of Z = 0.004–0.03 (roughly 20%–150% times the solar metallicity; Asplund et al. 2009). For the stellar masses based on the Bruzual & Charlot (2003) and FSPS (Conroy & Gunn 2010a, 2010b) models, we assume the Chabrier (2003) IMF. For the stellar masses derived from the Maraston (2005) models, we assume the Kroupa (2001) IMF (these choices of IMF produce systematic shifts in the derived stellar masses of ≃0.04 dex). Finally, we adopt the time-dependent dust-attenuation curve of Charlot & Fall (2000), in which stellar populations older than 10 Myr are attenuated by a factor of μ times less than younger stellar populations. Following Moustakas et al. (2013), we draw μ from an order four Gamma distributions that range from zero to unity centered on a typical value ($\left\langle \mu \right\rangle =0.3$).

We consider both smooth SFHs and superpositions of smooth SFHs with "bursts" of star formation with varying strengths and duration. Following Bundy et al. (2017) and Moustakas et al. (2013), we add stochastic bursts randomly to the SFHs. For every 2 Gyr interval over the lifetime of a given model, the cumulative probability that a burst occurs is P_burst = 0.2. The SFH of each burst is modeled as a Gaussian as a function of time, with an amplitude of F_b, defined as the total amount of stellar mass formed in the burst divided by the underlying mass of the smooth SFH at the burst's peak time. The values of F_b are drawn from the range of 0.03–4.0. The allowed burst duration (or the width of a Gaussian distribution characterizing the SF burst) uniformly ranges from 0.03 to 0.3 Gyr.

Additionally, we allow the time for the onset of star formation of each model to range with equal probability from 0.1 to 13 Gyr (Salim et al. 2007; Wild et al. 2009), but as above, we restrict these times to be less than the age of the universe at the redshift of each galaxy.

We then apply iSEDfit to the galaxies in our catalog with difference assumptions for the SFHs, stellar population libraries, and models. Table 1 describes the details of each set of runs. We adopt fiducial prior parameters for ${M}_{* ,\mathrm{iSED}}$ from Moustakas et al. (2013), marginalizing over all stellar population parameters, to produce posterior probability distribution functions (PDFs). We then adopt the median PDF as the reported value for each quantity and derive 68% confidence intervals by taking the values that correspond to the PDF integrated between 0.16 and 0.84. We refer to the stellar mass derived for each of these runs using the "name" listed in the first column of Table 1. We discuss how the different model assumptions impact both the stellar mass estimates and the derived SMF in Section 6.

Table 1.  Stellar Mass Estimates

			Prior
Name	Models	Star Formation History	Metallicity (Z)	Decay Timescale (τ/Gyr)	Age (t/Gyr)
(1)	(2)	(3)	(4)	(5)	(6)
${M}_{* ,\mathrm{iSED}}^{\mathrm{FSPS},\mathrm{burst}}$	FSPS	Exponentially	[0.004, 0.03]	[0.5–1]	[0.1–13]
	(Conroy & Gunn 2010b)	declining with Pburst = 0.2
${M}_{* ,\mathrm{iSED}}^{\mathrm{FSPS},\mathrm{no}\,\mathrm{burst}}$	FSPS	Exponentially	[0.004, 0.03]	[0.5–1]	[0.1–13]
	(Conroy & Gunn 2010b)	declining
${M}_{* ,\mathrm{iSED}}^{\mathrm{BC}03}$	Bruzual & Charlot (2003)	Exponentially	[0.004, 0.03]	[0.5–1]	[0.1–13]
		declining
${M}_{* ,\mathrm{iSED}}^{\mathrm{Ma}05}$	Maraston (2005)	Exponentially	[0.004, 0.03]	[0.5–1]	[0.1–13]
		declining

Note. (1) Name of the stellar mass estimate, (2) stellar population synthesis (SPS) model, (4) the prior on metallicity, (5) the prior on the exponentially declining star formation timescale, τ, in units of Gyr, (6) the prior on the time for the onset of star formation, t, in units of Gyr. Priors of the form [A, B] are flat, with the minimum and maximum given by A and B. P_burst denotes the cumulative probability that a burst occurs. We assume a Chabrier (2003) initial mass function.

Download table as:  ASCIITypeset image

In Figure 2, we show representative examples of SHELA galaxies with $\mathrm{log}({M}_{* }/{M}_{\odot })\gt 11$ in three redshift bins: 0.5 < z < 0.75, 0.75 < z < 1.0, and 1.0 < z < 1.5. The best-fit SEDs and photometry are based on FSPS models without stochastic bursts. Figure 2 also shows false-color images of the massive galaxies in the DECam z- (red color) combined with i- (green color) and g-band (blue color) images. By inspection, our massive galaxies are typically spheroidal, or reddened, bulge-dominated disks.

Zoom In Zoom Out Reset image size

Finally, we note that our stellar mass estimates (${M}_{* ,\mathrm{iSED}}$) refer to the stellar mass implied by the visible flux from the living stellar population within a galaxy and not the total living plus stellar remnants (i.e., white dwarf, neutron stars, black holes, etc.). Stellar remnants can make an important contribution to the total stellar mass of a galaxy and the SMFs for massive galaxies at low redshift (e.g., Shimizu & Inoue 2013; Bernardi et al. 2016). For example, Shimizu & Inoue (2013) found a weak correlation between the remnant mass fraction and the total stellar mass of galaxies, and the remnant fraction can be regarded as a redshift-dependent constant. Also, the shape of the SMF is almost unchanged, but simply shifts horizontally depending on the inclusion or omission of the remnant mass. This shift in the SMF is larger at lower redshift (∼0.05 dex at z = 3 and ∼0.15 dex at z = 0). However, we find that this difference is small and comparable to the SHELA stellar mass uncertainty, which we already take into account using a forward-modeling method. Therefore, the choice of including or excluding stellar remnants in the "stellar mass" does not significantly impact our inferred evolution of the SMFs.

The left panel of Figure 3 compares the stellar masses of FSPS, no-burst SFHs (${M}_{* ,\mathrm{iSED}}^{\mathrm{FSPS},\mathrm{no}\,\mathrm{burst}}$) to models with bursts (${M}_{* ,\mathrm{iSED}}^{\mathrm{FSPS},\mathrm{burst}}$). For these two cases, the median offset is ∼–0.02 dex, and the scatter is tight (0.01 dex). There is no measurable dependence on stellar mass for the range (${M}_{* }={10}^{9-12}\,{M}_{\odot }$).

Zoom In Zoom Out Reset image size

The middle of Figure 3 shows the comparison between the FSPS ${M}_{* ,\mathrm{iSED}}^{\mathrm{FSPS},\mathrm{no}\,\mathrm{burst}}$ stellar masses and those based on the Bruzual & Charlot (2003) SPS model (${M}_{* ,\mathrm{iSED}}^{\mathrm{BC}03}$). In both cases we compare models with smoothly varying SFHs. For these two cases, the median offset is −0.02 dex, with a weak dependence on stellar mass. The scatter is likewise small with σ = 0.02 dex. This result confirms our expectations that adding bursts primarily modifies the bluer bandpasses of a galaxy, while leaving the redder wavelengths, which count the accumulated stellar mass of a galaxy, relatively unaffected.

The right panel of Figure 3 compares the stellar masses of FSPS ${M}_{* ,\mathrm{iSED}}^{\mathrm{FSPS},\mathrm{no}\,\mathrm{burst}}$ models to those based on the Maraston (2005) (${M}_{* ,\mathrm{iSED}}^{\mathrm{Ma}05}$) stellar population models. These SPS models exhibit larger mass-dependent offsets, with a median of −0.15 dex and a scatter substantially larger than in the previous comparisons, σ = 0.09 dex. This is most likely due to the different treatment of the thermally pulsing asymptotic giant branch (TP-AGB) phase; the Maraston (2005) models have a much larger flux at wavelengths longward of ∼0.7–0.8 μm.

In summary, the stellar masses for the galaxies in our samples are fairly robust to changes in SFH or model stellar population (with the exception of the Maraston models). For this reason, we average the results from the SMFs derived from different mass estimates (defined as "assumption-averaged" SMF; see Section 6.1) to infer biases associated with the assumptions of the SFH and SPS models. We adopt the assumption-averaged SMF as our fiducial measurement. However, in Section 6.3, we further discuss how the differences in stellar masses from the different model assumptions impact the SMFs.

We estimated our stellar mass completeness limits by using Bruzual & Charlot (2003) models to generate a series of solar metallicity Simple Stellar Population (SSP) models with a formation redshift of z_f = 4. We then and used EZGal (Mancone & Gonzalez 2012) to infer the models' observables and compared these values to our Bruzual & Charlot (2003) models without bursts. The result is shown in Figure 4 as a function of redshift using the limits defined from our limiting NEWFIRM K-band magnitude of 22.7 AB mag (5σ depth in 2'' diameter apertures). We adopt this stellar mass completeness for the following analysis. We also computed the stellar mass completeness limits using other SFHs, SPS models, metallicities (Z), and exponential decay timescales (τ) spanning the parameters we used. We find a less than 0.2 dex change in the derived stellar mass completeness limits, consistent with the results from the previous section.

Zoom In Zoom Out Reset image size

Figure 4 also shows that the K band provides the deepest stellar mass limit at all redshifts we consider here, i.e., z < 1.5. For z < 1.0, the IRAC depth corresponds to slightly lower stellar mass limits. However, at 1.0 < z < 1.5, the K band provides a deeper stellar mass limit and is well matched to the DECam imaging. For this reason, we use the limit derived from the K-band data, which provide a galaxy sample "complete" to $\mathrm{log}{M}_{* }/{M}_{\odot }=10.0$ at z = 1 and $\mathrm{log}{M}_{* }/{M}_{\odot }$ = 10.3 at z = 1.5.

We measure SFRs for the galaxies in our sample using iSEDfit. Similar to the stellar mass estimates, we adopt the median of the posterior PDF marginalized over the SFR as the best estimate. We estimate SFRs for the four different sets of SFH and SPS models with the same parameters as we did for our stellar mass estimates (Table 1).

We classify the galaxy population as either star-forming or quiescent based on whether they lie on or below the so-called star-forming main sequence (Noeske et al. 2007). The star-forming (SF) main sequence is the correlation between the SFR and stellar mass of star-forming galaxies that has been observed out to z ∼ 2.5 (e.g., Rodighiero et al. 2011; Wuyts et al. 2011; Whitaker et al. 2012; Shivaei et al. 2015; Schreiber et al. 2016). In Figure 5, we plot SFR versus stellar mass in 11 redshift bins from z = 0.4–1.5 for our SHELA sample. The figure demonstrates the existence of a well-defined SF main sequence whose amplitude increases smoothly toward higher redshift, similar to other results (see, e.g., Speagle et al. 2014). Additionally, we find a distinct population of quiescent galaxies that lie below the SF main sequence at a given stellar mass.

Zoom In Zoom Out Reset image size

To classify the galaxy population, we use an evolving threshold of specific SFR (sSFR), computed as $\mathrm{sSFR}\,=\mathrm{SFR}/{M}_{* }$, to trace the lower envelope of the SF main sequence in each redshift bin. Specifically, we plot the distribution of SFR in each stellar mass and redshift bin, and model the bimodal distribution of each bin's SFR as the sum of two normal Gaussian functions (e.g., Strateva et al. 2001; Baldry et al. 2004). We then measure the mean (μ_SF) and dispersion (σ_SF) of the distribution of the star-forming population. We define all galaxies whose SFRs lie below 3σ_SF from μ_SF to be quiescent. This results in the threshold for sSFRs evolving from 10⁻¹¹ yr⁻¹ at z = 0.4 to 10^−10.2 yr⁻¹ at z = 1.5. We adopt the evolving threshold of sSFR to classify the galaxy population because the star-forming main sequence evolves with redshift—the main sequence as a whole moves to higher SFR as redshift increases (e.g., Noeske et al. 2007). However, in Section 7.4, we discuss the effect of using a nonevolving threshold of sSFR to classify galaxy population on the evolution of number densities of quiescent and star-forming populations.

We also must emphasize an important caveat of our derived SFRs. At z < 0.5, the u-band photometry samples the rest-frame wavelength of galaxies longer than near-ultraviolet (NUV), and we need NUV or far-UV observations to probe the recent star formation of galaxies. As a result, in our lowest redshift bins, we may have less accurate SFRs and also less reliable separation between star-forming and quiescent populations. In Section 7.2, we show that our interpretation of the number density and stellar mass density evolution may be impacted by this limitation.

Finally, in Figure 4, we show the distributions of redshift and stellar mass for our full SHELA sample and the subsamples of quiescent and star-forming galaxies. Over 0.4 < z < 1.5, the population of massive galaxies ($\mathrm{log}{M}_{* }/{M}_{\odot }\gt 11$) are dominated by quiescent systems. In Section 6.5, we will quantify how the evolution in the SMF of each subsample accounts for the evolution in the SMF for all massive galaxies.

Errors in the stellar mass estimates (M_*) introduce a bias into the derived galaxy SMF, as random errors cause more objects to have "upscattered" stellar masses than "downscattered" values. This is a form of Eddington (1913) bias and is especially problematic on the exponential (high-mass) end of the SMF, due to the steep decline in the number of galaxies. Even a small fraction of these upscattered low-mass galaxies can dominate the number densities at high stellar masses (Ilbert et al. 2013; Moster et al. 2013; Caputi et al. 2015; D'Souza et al. 2015; Grazian et al. 2015; Davidzon et al. 2017).¹³ This Eddington-type bias in stellar mass may increase with increasing redshift because of the decrease in signal-to-noise ratio with increasing redshift. Some previous studies have shown that Eddington bias can affect the interpretation of the evolution of the galaxy SMF (see discussion in, e.g., Fontanot et al. 2009; Moster et al. 2013; Bundy et al. 2017).

Here, we account for a varying Eddington bias using a forward-model method. This method requires that we assume an intrinsic shape of a galaxy SMF, to which we then apply measurement (and systematic) uncertainties. For our model, we first assume that the galaxy SMF is well described by a double Schechter function (Baldry et al. 2008) of the form

$$\begin{eqnarray}\begin{array}{rcl}\phi ({M}_{* }) & = & (\mathrm{ln}10)\exp \left[-\displaystyle \frac{{M}_{* }}{{M}^{* }}\right]\times \left\{{\phi }_{1}{10}^{({\alpha }_{1}+1)(\mathrm{log}{M}_{* }-\mathrm{log}{M}^{* })}\right.\\ & & +\left.{\phi }_{2}{10}^{({\alpha }_{2}+1)(\mathrm{log}{M}_{* }-\mathrm{log}{M}^{* })}\right\},\end{array}\end{eqnarray} \tag{ 2 }$$

where α₂ < α₁ are the faint-end (power-law) slopes of the SMF (where the second term dominates at the low-mass end). We denote the "knee" of the SMF with the characteristic stellar mass as M*, which marks the stellar mass above which the SMF declines exponentially. We require that both terms in the double Schechter function have the same M*.

We then construct a series of mock galaxy samples with the stellar mass distribution following this double Schechter function. We generate 500,000 mock data sets that sample various parameter ranges of the double Schechter function, then we perturb the stellar mass estimates using uncertainties drawn from a Gaussian distribution, with the width of the distribution equal to the 1σ uncertainty in the stellar masses based on the SED fitting, ${\sigma }_{{M}_{* }}(\mathrm{log}{M}_{* },z)$. This uncertainty is both stellar mass and redshift dependent. Figure 6 shows the stellar mass uncertainties that we used, including their dependence on redshift.

Zoom In Zoom Out Reset image size

We do not directly constrain the normalization ϕ₁ and ϕ₂ of the double Schechter function. Instead, we define the parameter λ_mix that is varied between between 0 and 1 to indicate the relation between the first term and the second term of the double Schechter form, such that λ_mix ∝ ϕ₁ and (1 − λ_mix) ∝ ϕ₂. We then evaluate the overall normalization factor C of the SMF of the mock sample $\phi ({M}_{* ,\mathrm{mock}})$such that the comoving number density of mock galaxies, n_mock,

$$\begin{eqnarray}&&{n}_{\mathrm{mock}}=C{\int }_{\mathrm{log}{M}_{* \min }}^{\mathrm{log}{M}_{* \max }}\phi ({M}_{* ,\mathrm{mock}})d\mathrm{log}{M}_{* },\end{eqnarray} \tag{ 3 }$$

equals the observed comoving number density of SHELA galaxies, n_SHELA, and ${M}_{* \min }={M}_{\mathrm{lim}}$ and ${M}_{* \max }={10}^{12.5}{M}_{\odot }$, where M_lim denotes the stellar mass completeness limit at a given redshift bin (see Table 3). We then bin the mock samples identically to the data.

Finally, we constrain the double Schechter parameters by performing a grid search over ranges of parameters (M*, α₁, α₂, λ_mix) and comparing the mock SMFs with the observed one in each stellar mass bin. We demonstrate our method of forward modeling the galaxy SMF in Figure 7.

Zoom In Zoom Out Reset image size

In addition to the scatter in stellar mass estimates, we also consider the effects of Poissonian (counting) uncertainties and "cosmic variance." The latter are large-scale fluctuations in the spatial distribution of the number of galaxies in the universe. While all of these effects contribute to the uncertainties in the SMF, cosmic variance is most significant for small fields and highly biased objects, i.e., objects with strong spatial clustering, such as massive galaxies. One advantage to our SHELA survey is that, as it spans such a large volume (1.4 × 10⁸ Mpc³ = 0.14 Gpc³), the effects of cosmic variance on the number density of massive galaxies are mitigated, and we can estimate cosmic variance from the data itself. For example, using the formalism set by Moster et al. (2011), the relative cosmic variance, σ_v, of galaxies more massive than 10¹¹ M_⊙ at z = 0.35 is 36% (∼0.1 dex) for COSMOS (≈2 deg²), but only ∼17% (∼0.07 dex) for SHELA.

Bundy et al. (2017) estimated the cosmic variance in S82-MGC (140 deg²) using bootstrap resampling. In the 0.3 < z < 0.65 bin, their estimated 1σ error due to the cosmic variance is ∼0.01 dex (corresponding to σ_v of ∼2%) at $\mathrm{log}({M}_{* }/{M}_{\odot })\sim 11.0$, and 0.02–0.05 dex (${\sigma }_{v}\sim 5 \% \mbox{--}10 \%$) at $\mathrm{log}({M}_{* }/{M}_{\odot })\sim 11.6$.

We adopt the method of Bundy et al. (2017) to estimate the cosmic variance in our SHELA samples. We divide the SHELA survey into 150 roughly equal-area regions and recompute SMFs after resampling with replacement (see Appendix A). For low-mass galaxies ($\mathrm{log}({M}_{* }/{M}_{\odot })\lt 11.5$), the bootstrap resampling yields a cosmic variance, σ_v, of 5%–12% at 0.3 < z < 1.5. At higher masses ($\mathrm{log}({M}_{* }/{M}_{\odot })\gt 11.5$), the cosmic variance rises to 6%–12% (Figure 23 in Appendix A). In the redshift range where we overlap with Bundy et al. (2017), our cosmic variance is larger. This is expected as the area of S82-MGC is eight times larger than that of the SHELA.

To test for other sources of systematic uncertainty, we compare our SMF from SHELA to that in S82-MGC (Bundy et al. 2017) for all galaxies between 0.3 < z < 0.65, where our samples overlap (see Section 7 and Figure 19). The two estimates, which both use the forward-modeling method, are in good agreement. The agreement is particularly good for the massive galaxies ($\mathrm{log}({M}_{* }/{M}_{\odot })\gtrsim 11.0$) between 0.4 < z < 0.6. This suggests that we are not significantly affected by either cosmic variance or other systematics at this redshift range. However, between 0.3 < z < 0.4, the normalization of the SHELA SMF is lower than that of S82-MGC. In this redshift bin, the estimated relative cosmic variance (Moster et al. 2011) in our SHELA survey is 20% (∼0.1 dex), compared to the 7% cosmic variance in the S82-MGC survey. Therefore, our results in this (smallest-volume) redshift bin may be affected by a larger than typical cosmic variance. In the following analysis, we omit this lowest redshift bin and only study the evolution of the SMF for galaxies at 0.4 < z < 1.5. We also note that we do not include the effect of cosmic variance in our forward-modeling method because the effects are correlated in stellar mass, and all mass bins should be equally affected by the same large-scale fluctuation. Consequently, the measurement of the galaxy SMF should be mainly affected by the random errors in the stellar mass estimates rather than cosmic variance.

To estimate the contamination from QSOs and AGNs on our SHELA galaxy SMFs, we first cross-match our sample with the SDSS spectral catalog (DR13; Albareti et al. 2017) to find 760 QSOs in our galaxy sample. Second, we cross-match the SHELA sample with the 31 deg² Stripe 82X X-ray Catalog (LaMassa et al. 2016). The flux limits of this Stripe 82 X-ray survey are 8.7 × 10⁻¹⁶ erg s⁻¹ cm⁻², 4.7 × 10⁻¹⁵ erg s⁻¹ cm⁻², and 2.1 × 10⁻¹⁵ erg s⁻¹ cm⁻² in the soft (0.5–2 keV), hard (2–10 keV), and full X-ray bands (0.5–10 keV), respectively (LaMassa et al. 2016). We found 1253 matched X-ray sources in our sample. In total, these sources account for less than 1% of galaxies more massive than 10¹¹M_⊙.¹⁴ We then exclude these QSOs and X-ray sources from our galaxy sample and repeat our measurement of galaxy SMFs. Because these sources account for only a few percent of the total number of massive galaxies, the resulting SMFs and the inferred redshift evolution are not significantly affected by those contaminations. Therefore, we conclude that our measurements are not adversely affected by the presence of AGNs. For the rest of this paper, we exclude SDSS QSOs and X-ray sources from our galaxy sample.

We define the "assumption-averaged" SMF, i.e., the average of the SMFs derived using the different methods listed in Table 1, as our fiducial measurement. In practice, we take the average number density in each bin of stellar mass using the different assumptions in the SFH and SPS models.¹⁵ In this way, our results are a statistical mean, and we can study the variance in the SMF. These four M_* estimates encompass the range of M_* values obtained by adopting currently uncertain priors. Therefore, the assumption-averaged result represents a compromise among differing approaches.

In Figure 10, we show the observed galaxy SMF in each redshift bin with shaded regions corresponding to the Poisson errors computed by taking the square root of the number of galaxies in the stellar mass bin. We indicate the stellar mass completeness limit at each redshift bin with a vertical dotted line. We present our measurements of the assumption-averaged SMF and the number of galaxies in each redshift bin in Tables 2 and 3, respectively.

Zoom In Zoom Out Reset image size

Table 2.  Assumption-averaged Stellar Mass Functions for All Galaxies

	0.4 < z < 0.5	0.5 < z < 0.6	0.6 < z < 0.7	0.7 < z < 0.8	0.8 < z < 0.9
$\mathrm{log}({M}_{* }/{M}_{\odot })$	$\mathrm{log}(\phi /{\mathrm{Mpc}}^{-3}\,{\mathrm{dex}}^{-1}$)	$\mathrm{log}(\phi /{\mathrm{Mpc}}^{-3}\,{\mathrm{dex}}^{-1}$)	$\mathrm{log}(\phi /{\mathrm{Mpc}}^{-3}\,{\mathrm{dex}}^{-1}$)	$\mathrm{log}(\phi /{\mathrm{Mpc}}^{-3}\,{\mathrm{dex}}^{-1}$)	$\mathrm{log}(\phi /{\mathrm{Mpc}}^{-3}\,{\mathrm{dex}}^{-1}$)
9.53	−2.22 ± 0.01	⋯	⋯	⋯	⋯
9.68	−2.27 ± 0.01	−2.31 ± 0.01	⋯	⋯	⋯
9.83	−2.32 ± 0.01	−2.36 ± 0.01	−2.49 ± 0.01	⋯	⋯
9.98	−2.38 ± 0.01	−2.43 ± 0.01	−2.53 ± 0.01	−2.59 ± 0.01	−2.54 ± 0.01
10.13	−2.45 ± 0.01	−2.49 ± 0.01	−2.56 ± 0.01	−2.61 ± 0.01	−2.57 ± 0.01
10.28	−2.50 ± 0.01	−2.57 ± 0.01	−2.59 ± 0.01	−2.62 ± 0.01	−2.59 ± 0.01
10.43	−2.52 ± 0.01	−2.62 ± 0.01	−2.59 ± 0.01	−2.62 ± 0.01	−2.60 ± 0.01
10.58	−2.53 ± 0.01	−2.64 ± 0.01	−2.60 ± 0.01	−2.61 ± 0.01	−2.64 ± 0.01
10.73	−2.58 ± 0.01	−2.64 ± 0.01	−2.60 ± 0.01	−2.64 ± 0.01	−2.65 ± 0.01
10.88	−2.66 ± 0.01	−2.69 ± 0.01	−2.69 ± 0.01	−2.69 ± 0.01	−2.70 ± 0.01
11.03	−2.78 ± 0.01	−2.81 ± 0.01	−2.81 ± 0.01	−2.81 ± 0.01	−2.81 ± 0.01
11.18	−2.99 ± 0.01	−2.99 ± 0.01	−3.04 ± 0.01	−3.03 ± 0.01	−3.02 ± 0.01
11.33	−3.32 ± 0.02	−3.25 ± 0.02	−3.35 ± 0.02	−3.35 ± 0.02	−3.31 ± 0.01
11.48	−3.80 ± 0.04	−3.73 ± 0.03	−3.77 ± 0.03	−3.82 ± 0.03	−3.76 ± 0.02
11.63	−4.43 ± 0.07	−4.37 ± 0.06	−4.41 ± 0.05	−4.28 ± 0.04	−4.28 ± 0.04
11.78	−5.07 ± 0.14	−5.08 ± 0.12	−5.00 ± 0.10	−5.16 ± 0.11	−5.02 ± 0.09
11.93	−5.59 ± 0.22	⋯	−5.92 ± 0.24	−5.84 ± 0.21	−5.94 ± 0.22
12.08	−5.94 ± 0.30	−6.07 ± 0.30	−6.28 ± 0.33	−6.14 ± 0.28	−6.20 ± 0.28
0.9 < z < 1.0	1.0 < z < 1.1	1.1 < z < 1.2	1.2 < z < 1.3	1.3 < z < 1.4	1.4 < z < 1.5
$\mathrm{log}(\phi /{\mathrm{Mpc}}^{-3}\,{\mathrm{dex}}^{-1}$)	$\mathrm{log}(\phi /{\mathrm{Mpc}}^{-3}\,{\mathrm{dex}}^{-1}$)	$\mathrm{log}(\phi /{\mathrm{Mpc}}^{-3}\,{\mathrm{dex}}^{-1}$)	$\mathrm{log}(\phi /{\mathrm{Mpc}}^{-3}\,{\mathrm{dex}}^{-1}$)	$\mathrm{log}(\phi /{\mathrm{Mpc}}^{-3}\,{\mathrm{dex}}^{-1}$)	$\mathrm{log}(\phi /{\mathrm{Mpc}}^{-3}\,{\mathrm{dex}}^{-1}$)
⋯	⋯	⋯	⋯	⋯	⋯
⋯	⋯	⋯	⋯	⋯	⋯
⋯	⋯	⋯	⋯	⋯	⋯
⋯	⋯	⋯	⋯	⋯	⋯
−2.52 ± 0.01	−2.53 ± 0.01	−2.57 ± 0.01	⋯	⋯	⋯
−2.55 ± 0.01	−2.57 ± 0.01	−2.63 ± 0.01	−2.67 ± 0.01	−2.70 ± 0.01	⋯
−2.58 ± 0.01	−2.62 ± 0.01	−2.70 ± 0.01	−2.74 ± 0.01	−2.74 ± 0.01	−2.73 ± 0.01
−2.59 ± 0.01	−2.65 ± 0.01	−2.77 ± 0.01	−2.82 ± 0.01	−2.83 ± 0.01	−2.81 ± 0.01
−2.60 ± 0.01	−2.65 ± 0.01	−2.77 ± 0.01	−2.87 ± 0.01	−2.91 ± 0.01	−2.90 ± 0.01
−2.65 ± 0.01	−2.67 ± 0.01	−2.79 ± 0.01	−2.92 ± 0.01	−2.97 ± 0.01	−2.95 ± 0.01
−2.76 ± 0.01	−2.80 ± 0.01	−2.91 ± 0.01	−3.02 ± 0.01	−3.06 ± 0.01	−3.05 ± 0.01
−2.97 ± 0.01	−2.99 ± 0.01	−3.14 ± 0.01	−3.23 ± 0.01	−3.25 ± 0.01	−3.23 ± 0.01
−3.27 ± 0.01	−3.28 ± 0.01	−3.40 ± 0.01	−3.57 ± 0.02	−3.56 ± 0.02	−3.55 ± 0.01
−3.70 ± 0.02	−3.73 ± 0.02	−3.83 ± 0.02	−3.94 ± 0.02	−4.01 ± 0.03	−4.01 ± 0.02
−4.32 ± 0.04	−4.31 ± 0.04	−4.43 ± 0.04	−4.50 ± 0.04	−4.53 ± 0.04	−4.57 ± 0.05
−5.01 ± 0.08	−5.16 ± 0.10	−5.12 ± 0.09	−5.19 ± 0.09	−5.08 ± 0.08	−5.19 ± 0.09
−6.16 ± 0.26	−6.20 ± 0.26	−5.97 ± 0.20	−5.89 ± 0.19	−5.70 ± 0.15	−5.77 ± 0.16
−6.16 ± 0.26	⋯	−7.01 ± 0.48	⋯	−6.46 ± 0.30	−6.23 ± 0.24

Download table as:  ASCIITypeset image

Table 3.  Intrinsic Mass Function Shape Parameters from Forward Modeling for All Galaxies

Redshifts	$\mathrm{log}({M}_{\mathrm{lim}}/{M}_{\odot })$	Ngal	$\mathrm{log}({\phi }_{1}/{\mathrm{Mpc}}^{-3}\,{\mathrm{dex}}^{-1})$	$\mathrm{log}({\phi }_{2}/{\mathrm{Mpc}}^{-3}\,{\mathrm{dex}}^{-1})$	$\mathrm{log}({M}^{* }/{M}_{\odot })$	${\alpha }_{1}$	${\alpha }_{2}$
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
(0.3,0.4)	9.22	33848	−2.03	−2.98	${10.70}_{-0.03}^{+0.01}$	−0.45 ± 0.17	−1.70 ± 0.10
(0.4,0.5)	9.44	36283	−2.15	−2.75	${10.76}_{-0.03}^{+0.02}$	−0.12 ± 0.25	−1.50 ± 0.20
(0.5,0.6)	9.60	36389	−2.25	−2.72	${10.82}_{-0.00}^{+0.02}$	0.00 ± 0.16	−1.45 ± 0.05
(0.6,0.7)	9.72	36200	−2.32	−2.92	${10.76}_{-0.04}^{+0.03}$	−0.10 ± 0.21	−1.35 ± 0.11
(0.7,0.8)	9.83	36114	−2.39	−3.00	${10.73}_{-0.01}^{+0.01}$	0.00 ± 0.03	−1.25 ± 0.08
(0.8,0.9)	9.92	39318	−2.36	−3.31	${10.87}_{-0.02}^{+0.02}$	−0.57 ± 0.05	−1.27 ± 0.03
(0.9,1.0)	9.98	45584	−2.33	−3.29	${10.82}_{-0.03}^{+0.02}$	−0.45 ± 0.09	−1.30 ± 0.05
(1.0,1.1)	10.05	43303	−2.37	−3.32	${10.84}_{-0.02}^{+0.02}$	−0.47 ± 0.07	−1.27 ± 0.03
(1.1, 1.2)	10.11	34655	−2.45	−3.40	${10.86}_{-0.01}^{+0.01}$	−0.65 ± 0.03	−1.25 ± 0.08
(1.2, 1.3)	10.16	28914	−2.55	−3.50	${10.82}_{-0.03}^{+0.02}$	−0.55 ± 0.05	−1.25 ± 0.08
(1.3, 1.4)	10.23	25776	−2.63	−3.58	${10.78}_{-0.00}^{+0.02}$	−0.35 ± 0.10	−1.25 ± 0.08
(1.4, 1.5)	10.30	24697	−2.70	−3.65	${10.74}_{-0.02}^{+0.03}$	−0.10 ± 0.13	−1.25 ± 0.08

Note. (1) Redshift range used for the SMF, (2) $\mathrm{log}({M}_{\mathrm{lim}}/{M}_{\odot })$ denotes the stellar mass completeness limit of each redshift bin, (3) N_gal denotes the number of galaxies with stellar mass above the stellar mass completeness limit at each redshift bin, (4) the normalization, ${\phi }_{1}$, of a double Schechter function, (5) the normalization, ϕ₂, (6) the characteristic stellar mass, (7) the power-law slope of the high-mass end, and (8) the power-law slope of the low-mass end. The uncertainties of the double Schechter parameters are from the forward-modeling fits to the assumption-averaged SMF. The normalization of a double Schechter function (ϕ₁ and ϕ₂) does not have error bars because we do not directly constrain the normalization. Instead, we define the parameter λ_mix that is varied between 0 and 1 to indicate the relation between of the first term and the second term of the double Schechter form (Equation (2)). We then scale the stellar mass function of the mock sample so that its comoving number density is equal to the observed comoving number density of SHELA galaxies in each redshift bin (see Section 5.1).

Download table as:  ASCIITypeset image

We further perform forward modeling on the observed galaxy SMF in each redshift bin as described in Section 5. Figure 10 compares the measured SMF, which includes the measurement uncertainties, and the fitted intrinsic SMF. The forward modeling involves random draws from the estimated error distributions of stellar masses (${\sigma }_{{M}_{* }}$); as a result, the intrinsic models can vary from run to run with a scatter consistent with the error bars indicated on the observed SMFs.

To quantify the evolution in the SMF, we first present the evolution of the characteristic stellar mass (M*) resulting from the forward modeling of the assumption-averaged mass function (Figure 11). Within the systematic uncertainty due to the different stellar mass estimators, we detect no redshift evolution in M* (≲0.1 dex) from z = 1.5 to z = 0.4 even after accounting for the Eddington bias caused by random errors in stellar mass measurement.

Zoom In Zoom Out Reset image size

Second, we derive the cumulative number density of galaxies with stellar mass greater than 10¹¹ M_⊙ by integrating the intrinsic SMF inferred from the forward modeling,

$$\begin{eqnarray}&&n(\gt {M}_{* 11})={\int }_{{M}_{* 11}}^{\infty }\phi ({M}_{* }){{dM}}_{* }.\end{eqnarray} \tag{ 4 }$$

In practice, we use an upper limit of the integral of M_*max = 10^12.5 M_⊙, because our catalog contains no objects at higher stellar mass.

We note that the cumulative number density is less sensitive to the degeneracy between the characteristic stellar mass, M*, and other derived Schechter parameters. In Figure 11, we plot the cumulative number density of galaxies with stellar mass greater than 10¹¹ M_⊙ and the corresponding 68% range. We find no significant evolution (≲0.1 dex) in these densities out to z < 1. In contrast, at higher redshift, the cumulative number density of massive galaxies increases by ≲0.3 dex from z = 1.5 to z = 1.0.

We further use the abundance matching technique to identify galaxy cumulative number densities with dark matter halo cumulative number densities and estimate the evolution in the median cumulative number density of the progenitors of 10¹¹ M_⊙ galaxies at z = 0.4. We specifically implement the Number Density Redshift Evolution Code (NDE;¹⁶ Behroozi et al. 2013a, 2013b) to convert the cumulative number density of 10¹¹ M_⊙ galaxies at z = 0.4 (resulting from integrating the modeled intrinsic assumption-averaged SMF) to number densities at higher redshifts. Similarly, we estimate the evolution in median cumulative number density of the descendants of 10¹¹ M_⊙ galaxies at z = 1.5. The predicted evolution in the cumulative number density of the progenitors of 10¹¹ M_⊙ galaxies at z = 0.4 is consistent with that found from forward modeling the SHELA SMF (Figure 11) at z < 1. However, the predicted evolution in the cumulative number density of the progenitors of 10¹¹ M_⊙ galaxies stays constant out to z = 1.5, which we do not observe.

The discrepancy between the predicted evolution in the cumulative number density of the progenitors of ${10}^{11}\,{M}_{\odot }$ galaxies at z = 0.4 and the observed evolution, particularly at z > 1, could arise from assumptions of the NDE. This code ignores scatter in mass accretion and galaxy–galaxy mergers histories. This can lead to errors when comparing the evolution of galaxies over large redshift ranges (Δz > 1). This is evidenced by comparing the evolution over the redshift range of 0.4 < z < 1.5 (but the predicted evolution in the number density of the progenitors and the observed evolution is more consistent at 0.4 < z < 1). The scatter in stellar mass at fixed halo mass will influence the inferred 1σ range of cumulative number densities for galaxy progenitors and descendants. The NDE assumes that the growth in the differences in the ranked order of galaxy stellar mass are the same as the growth in the differences in the ranked order of halo mass as a function of time. In reality, this would be the case if the star formation efficiency in individual galaxies depends much more on halo mass than on cosmic time or environment (Behroozi et al. 2013c).

To better compare the relative evolutionary trend, we normalized the cumulative number density of galaxies with mass $\mathrm{log}{M}_{* }/{M}_{\odot }\gt 11.0$ at each redshift bin to that at z = 0.4 (Figure 12). This is the lowest bin where the comparison between SHELA SMFs and S82-MGC (Bundy et al. 2017, 140 deg²) shows that the SHELA sample is not significantly affected by cosmic variance (see Figure 19). The assumption-average mass function suggests no more than a ≲0.1 dex increase in the cumulative number density of these massive galaxies since z = 1.0 relative to those at z = 0.4. On the other hand, the number density of galaxies increases by ≲0.3 dex from z = 1.5 to z = 1 relative to those at z = 0.4. In the following section, we further explore the impact of the different SFH priors and SPS models on the SMFs and their redshift evolution.

Zoom In Zoom Out Reset image size

In the previous section, we derive galaxy SMFs by averaging all different sets of M_* measurements that include our various assumptions for the SFH and SPS models. Within the systematic uncertainty due to the different stellar mass estimators, we detect no redshift evolution in either the characteristic stellar mass (M*) or the cumulative number density of massive galaxies ($\mathrm{log}({M}_{* }/{M}_{\odot })\gt 11$) over 0.4 < z < 1.0 (even after accounting for the Eddington bias using the forward-modeling method). In this section, we further investigate the evolution of the SMFs derived using specific sets of M_* measurements.

We first consider how the SMF changes based on different assumptions for the SFH. In Figure 8, we show galaxy SMFs derived using stellar masses from different SFHs, including the effects of bursts.

It is clear that the effects of SFH, at least among the set of stellar mass estimates used here, on the SMF are minor. In Figure 11, we illustrate this by comparing the SMFs based on FSPS models with bursts (${M}_{* ,\mathrm{iSED}}^{\mathrm{FSPS},\mathrm{burst}}$) with those with no bursts (${M}_{* ,\mathrm{iSED}}^{\mathrm{FSPS},\mathrm{no}\,\mathrm{burst}}$). While there is a slight (≲0.1 dex) increase in characteristic mass derived from SFHs that allow bursts, this falls within the range of uncertainties, and we do not consider it significant. In both cases where we measure the SMF from stellar masses with SFHs that allow and disallow bursts, we find that the characteristic mass shows little evolution from z = 0.4 to 1.5.

In Figure 12, we show the cumulative comoving number density of massive galaxies (${M}_{* }\gt {10}^{11}{M}_{\odot }$) derived using stellar masses with different SFHs. In both cases, we normalize the results to the measurement at z = 0.4. Over 0.4 < z < 1.0, the normalized comoving number density of massive galaxies is approximately constant regardless of SFH. Similarly, the number densities drop by ≲0.3 dex from z = 1.0 to z = 1.5 relative to that at z = 0.4 if we use stellar masses from SFHs that allow bursts. We conclude that the systematic uncertainties arising from the choice of SFH contribute <0.1 dex to the error budget in the growth of the characteristic stellar mass of massive galaxies, which we determined from the combined assumption-average mass function.

In Figure 9 we evaluate how three choices for the stellar population models underlying iSEDfit M_* estimates impact the derived SMFs and constraints on the growth of massive galaxies. In all cases, we compare only models with smoothly varying SFHs (e.g., no bursts). We show again the FSPS ${M}_{* ,\mathrm{iSED}}^{\mathrm{FSPS},\mathrm{no}\,\mathrm{burst}}$ SMF in the left panel. We compare these to the SMFs based on Bruzual & Charlot (2003) masses (${M}_{* ,\mathrm{iSED}}^{\mathrm{BC}03}$, middle panel) and Maraston (2005) masses (${M}_{* ,\mathrm{iSED}}^{\mathrm{Ma}05}$, right panel). The different SPS models lead to different trends in terms of the redshift of the evolution of both the characteristic mass (M*) of the SMF and comoving number density of galaxies with stellar mass >10¹¹M_⊙.

There are competing claims as to the ability of the treatment of TP-AGB stars in the Maraston (2005) models to reproduce the colors of galaxies and star clusters. Kriek et al. (2010) found that the Maraston (2005) models could not simultaneously reproduce the rest-frame optical and near-IR portions of galaxy SEDs. Similarly, Conroy & Gunn (2010a) showed that the Maraston (2005) models produce redder colors at intermediate ages inconsistent with the colors of star clusters in the Magellanic Clouds. However, Capozzi et al. (2016) argued that the Maraston (2005) models fit better the SEDs of a sample of high-redshift galaxies in COSMOS with spectroscopic redshifts and that therefore the contribution from TP-AGB stars remains an important component in galaxy models. These points illustrate that uncertainties in stellar population models (in particular the treatment of TP-AGB stars) is an important component of the total error budget in the evolution of the SMF. We therefore include the results from Maraston (2005) fits with those from the FSPS and Bruzual & Charlot (2003) models in our assumption-averaged SMF (see below) to estimate the effect of uncertainties in the stellar population models to our results.

The SMFs based on all SPS models we are using in this study exhibit a ≲0.1 dex change in characteristic stellar mass from z = 1.0 to z = 0.4. On the other hand, at z > 1, the SMF based on Maraston (2005) masses (${M}_{* ,\mathrm{iSED}}^{\mathrm{Ma}05}$) exhibits a ∼0.2 dex increase in characteristic mass (Figure 13) from z = 1.5 to z = 1.0. This evolution is milder for the stellar masses based on the Bruzual & Charlot (2003) models (${M}_{* ,\mathrm{iSED}}^{\mathrm{BC}03}$) and those based on FSPS (${M}_{* ,\mathrm{iSED}}^{\mathrm{FSPS},\mathrm{no}\,\mathrm{burst}}$).

Zoom In Zoom Out Reset image size

We derive the cumulative number density of massive galaxies (>10¹¹M_⊙) normalized to that at z = 0.4 for each SPS model (Figure 12). Similar to the observed evolution in the characteristic stellar mass, at z < 1.0 we find almost no evolution in the cumulative number density of massive galaxies with stellar mass >10¹¹M_⊙ based on the stellar masses from any of the stellar population models. On the other hand, at z > 1.0, Maraston (2005) models exhibit a 0.4 dex increase in the number density of massive galaxies since z = 1.5.

The larger evolution of the SMF based on Maraston (2005) models likely arises from the different prescriptions for the TP-AGB stars in these models (compared to the assumptions used by the Bruzual & Charlot 2003 and Conroy & Gunn 2010b FSPS models). For galaxies at z > 0.8, the fits using the Maraston (2005) models have ages near 0.5–2 Gyr where the effects of the TP-AGB stars are most pronounced. This lowers the stellar M/L ratios of the models (Maraston et al. 2006). As a result, SED fitting with these models produce fits with lower stellar masses (${M}_{* ,\mathrm{iSED}}^{\mathrm{Ma}05}$) than those of the other stellar population models (see Figure 3). This reduces the number density of massive galaxies at the high-mass end, yielding increased evolution in both characteristics stellar mass (M*) and the number density. Additionally, ${M}_{* ,\mathrm{iSED}}^{\mathrm{Ma}05}$ has a larger spread relative to the other mass estimates, and this may contribute to the evolution in the SMF.

In our analysis of the evolution of the SMF, we average the results from the different stellar population codes (Conroy & Gunn 2010b, FSPS; Bruzual & Charlot 2003, and Maraston 2005). Given the different treatments in the prescription of the TP-AGB phases (which leads to a stronger redshift evolution in the number density of massive galaxies), this highlights how our uncertainties in the details of the later stages of stellar evolution propagate into uncertainties on measures of galaxy evolution such as the galaxy SMF.

We conclude from the combined assumption-average SMF in Figure 13 (left panel) that systematic uncertainties arising from the choice of SPS model contribute ≲0.1 dex to the error budget in the growth of the characteristic stellar mass of massive galaxies at z < 1 and 0.2 dex at 1.0 < z < 1.5. In addition, at least among the set of stellar mass estimates used here, the difference in SPS models is more important and lead to a larger variance in the implied number density evolution than the assumptions in the SFH.

We further quantify how the evolution in SMF depends on different SFH priors in various stellar mass bins. In Figures 14 and 15, we plot the number density of galaxies versus redshift in four bins with stellar mass between 10^10.4 and 10^12.5 M_⊙. Each redshift bin is normalized to that at z = 0.4. These figures show the evolutionary trend for a set of different priors in SFH and a set of different SPS models, respectively. The shaded region shows the 68th percentile range allowed by the different model assumptions in deriving the stellar masses. The number density of galaxies increases at different rates, with a dependence on stellar mass.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This leads to one of the main conclusions in this work: the number density of galaxies more massive than 10¹¹ M_⊙ does not change significantly over 0.4 < z < 1.0 where our sample is complete. Galaxies with masses between 10^10.4 and 10¹¹ M_⊙ show weak evidence for number density evolution. In contrast, there is a decline in the number density of these massive galaxies of ≃0.25–0.50 dex from z = 1 to 1.5. In each stellar mass bin, the individual evolutionary trends are consistent with one another at the ±1σ level for the different choices of SFH priors and SPS models used in this study (with the exception of the results using the Maraston 2005 models). We will discuss the implications of this in the next section.

In the previous subsections we measured the evolution of the SMF for the global population of galaxies from z = 0.4 to 1.5. We found no significant change in both the characteristic stellar mass and the cumulative number density of galaxies more massive than 10¹¹ M_⊙ at z < 1.0. At these redshifts, the evolution at lower stellar masses 10^10.4–10^10.7 M_⊙ is largest (≲0.1 dex) relative to the higher mass bins. One explanation for this difference is that a higher fraction of the lower mass galaxies are still star-forming (and therefore the number density of galaxies at fixed stellar mass grows with time). We therefore compare the evolution in the SMF for galaxies that are star-forming and those that are quiescent. To make this classification, we use median of the SFR posteriors reported by iSEDfit and compute the specific SFR as sSFR = SFR/M_*. We then divide the sample into quiescent and star-forming galaxies using the evolving sSFR threshold described in Section 4.4.

We recompute the SMFs for the quiescent galaxies and star-forming galaxies using our forward-modeling method. We plot the results in Figures 16 and 17 and present our measurements in Appendix B. Additionally, we compare the evolution of forwarded-modeled intrinsic SMFs for both populations in Figure 18. For massive quiescent galaxies with stellar mass ≳10¹¹ M_⊙, we do not detect growth (≲0.1 dex) in the characteristic stellar mass from z = 1.5 to z = 0.4. However, there is strong evolution in the number density of lower mass quiescent galaxies, similar to that seen in other studies (e.g., Tomczak et al. 2014; Moutard et al. 2016). This build-up in low-mass quiescent galaxies is expected to occur from the quenching of satellites (see discussion in Kawinwanichakij et al. 2017; Papovich et al. 2018).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the star-forming population, we find moderate growth of ∼0.1–0.2 dex in the characteristic stellar mass from z = 1.5 to z = 0.4, and this growth in M_* is larger at higher redshift. This is also consistent with previous studies (e.g., Tomczak et al. 2014; Moutard et al. 2016). However, in this paper, we focus on the evolution of the total SMF, and we save a detailed comparison of the SMF as a function of star formation activity for a future paper.

Figures 16–18 also show that, over 0.4 < z < 1.5, the population of massive galaxies ($\mathrm{log}{M}_{* }/{M}_{\odot }\gt 11$) is dominated by quiescent systems. It is therefore the evolution of this population that must account for the (lack of) evolution in the SMF of all massive galaxies. The evolution of quiescent galaxies will involve mass losses from stellar evolution processes (Girardi et al. 2000), and we expect additional mass growth either by mergers and/or the quenching of massive star-forming galaxies, as there appear to be too few of the latter at z < 1. In the next section, we discuss the implication of our finding on the rate of mass growth for these massive galaxies through merging.

The lack of evolution in the SMF of massive galaxies places constraints on models of galaxy growth and evolution. According to the two-phase formation scenario for the formation of massive galaxies (e.g., Oser et al. 2010, 2012), mass assembly at late times is dominated by minor mergers (e.g., Bédorf & Portegies Zwart 2013; Hilz et al. 2013; Laporte et al. 2013; Oogi & Habe 2013). We also expect some mass loss of quiescent galaxies from stellar evolutionary processes and dynamical processes in clusters. The lack of observed evolution in the SMF of massive galaxies could be a balance between the build-up of stellar mass through mergers and mass loss due to stellar evolution processes.

We can estimate the rate at which these massive galaxies grow by mergers at this late epoch using the results of Moster et al. (2013), who provided a parameterization for the SFH and mass accretion for galaxies of arbitrary present-day stellar mass. We integrate the Moster et al. (2013) fitting functions with respect to time and account for mass losses from passive stellar evolution (see Moster et al. 2013, their Equation (16)) to derive the expected stellar mass evolution of galaxies. For systems with a present-day stellar mass of $\mathrm{log}{M}_{* }/{M}_{\odot }=11$, the fraction of stellar mass loss relative to the stellar mass growth (both due to star formation and mass accretion) from z = 1.0 to z = 0.4 is 39%–45%. To be consistent with our measurements and to account for the lack of evolution of the SMF of massive galaxies, the upper limit on the amount of mass growth from mergers from z = 1.0 to z = 0.4 must be ∼45% (≃0.16 dex).

Our estimate of mass growth by mergers is in good agreement with a study by van Dokkum et al. (2010), who used a stacking analysis to study the growth of massive galaxies with a constant number density of 2 × 10⁻⁴ Mpc⁻³, corresponding to a galaxy with a stellar mass of 3 × 10¹⁰ M_⊙. At 0.6 < z < 0.1, van Dokkum et al. found ∼0.1 dex mass growth for these massive systems. In addition, Marchesini et al. (2014) used the UltraVISTA catalogs to investigate the evolution of the progenitors of local ultramassive galaxies ($\mathrm{log}({M}_{* }/{M}_{\odot })\approx 11.8;$ hereafter UMGs). They selected progenitors using the semiempirical approach of abundance matching and found a growth in stellar mass of ${0.27}_{-0.12}^{+0.08}$ dex from z = 1 to z = 0 after including the scatter in the progenitor's number density in the error budget. Marchesini et al. also found that half of the assembled stellar mass of local UMGs is formed primarily by merging over this redshift range. Our infer stellar mass growth and that of Marchesini et al. (2014) is consistent within the range of the uncertainties.

On the other hand, Ownsworth et al. (2014) presented a study on the stellar mass growth for the progenitors of galaxies with ${M}_{* }={10}^{11.24}{M}_{\odot }$ at z = 0.3 and showed that these massive galaxies have grown by a factor of ∼1.8 (∼0.25 dex) in total stellar mass since z = 1.0. They also found that, on average, major and minor mergers account for ∼17% and ∼34% of the mass assembled to galaxies at z = 0.3, respectively. In contrast, the process of star formation accounts for ∼24% of the total stellar mass. We observe a lower rate of mass growth from mergers compared to that from Ownsworth et al. (2014), and this discrepancy may result from the different SED-modeling assumptions and the manner with which we have estimated the effects of the Eddington bias.

Our finding can be directly compared to the recent study by Bundy et al. (2017), who followed a similar analysis as we have here. Bundy et al. (2017) detected no growth (with an uncertainty of 9%) in the characteristic stellar mass of massive galaxies ($\mathrm{log}({M}^{* }/{M}_{\odot })\gt 11.2$) from z = 0.65 to z = 0.3 in S82-MGC. We reproduce their observed SMFs and find a consistent intrinsic SMF after accounting for the Eddington bias (Figure 19). Our number density of massive galaxies ($\mathrm{log}({M}_{* }/{M}_{\odot })\gt 11.5$) is lower than that of Bundy et al. but because of the larger volume probed by the S82-MGC, particularly at 0.3 < z < 0.4, and the possible effects of cosmic variance, our results are still consistent.

Capozzi et al. (2017) studied the evolution of the galaxy SMF since z = 1, using ∼155 deg² of the Dark Energy Survey. In good agreement with our finding, Capozzi et al. find that the number densities of galaxies with $\mathrm{log}({M}^{* }/{M}_{\odot })\gt 11$ are constant from z ∼ 1 to z ∼ 0.2. In addition, these authors also find a mass dependence of the galaxy number density—less massive galaxies exhibit larger evolution in the number density compared to more massive galaxies. Again, this is qualitatively consistent with our finding (see Section 6.4), and we are able to verify the robustness of these results by fully accounting for the statistical and systematic uncertainties on stellar mass estimates. Also, the deep mid-infrared photometry from Spitzer/IRAC allows us to better constrain stellar masses and improve the uncertainties.

Moutard et al. (2016) presented an analysis of the evolution of the SMF from redshift z = 0.2 to z = 1.5 of a K_s < 22 mag-selected sample, over an effective area of ∼22.4 deg² of the VIPERS Multi-Lambda Survey. To account for scatter in the stellar mass measurements, Moutard et al. (2016) corrected the SMF during their fitting procedure by convolving the parametric form of the SMF with the stellar mass uncertainty (Ilbert et al. 2013). Moutard et al. showed that the number density of the most massive galaxies ($\mathrm{log}({M}_{* }/{M}_{\odot })\gt 11.5$) increases by a factor of ∼2 from z ∼ 1 to z ∼ 0.3. The higher number density of massive galaxies inferred by Moutard et al. compared to our finding could arise from the different methods used to account for the scatter in stellar mass measurement and our inclusion of Spitzer/IRAC mid-IR measurements in our stellar mass determinations. Moutard et al. also demonstrated that the quiescent population largely dominates the massive galaxies population since z ∼ 1; this agrees with our result that the massive galaxies assemble their stellar masses through mergers.

In this section, we compute the evolution of the cumulative number density of galaxies as a function of stellar mass threshold and compare it with other studies. Here, we are particularly interested in the redshift evolution measured within our SHELA field, as it mitigates against systematic uncertainties in the analysis between other surveys (and our survey is one of the few that attempt to forward model the SMF).

We integrate our best-fitting intrinsic SHELA SMF (i.e., the best-fit SMF derived from the forward modeling of the assumption-averaged SMF) for stellar masses greater than $\mathrm{log}{M}_{* }/{M}_{\odot }$ = 10, 10.5, 11, and 11.5. In Figure 20, we show the cumulative number density of galaxies in the four mass bins. For galaxies at all masses ($\mathrm{log}{M}_{* }/{M}_{\odot }\gt 10$), the cumulative number density in our sample is consistent with little evolution (as we have showed above). We see here that for galaxies more massive than $\mathrm{log}{M}_{* }/{M}_{\odot }\gt 10.5$, this (lack of) evolution in number density is primarily due to quiescent galaxies, which show a constant number density (out to z < 1). In contrast, star-forming galaxies show an increase in cumulative number density at all redshifts and stellar masses. This is consistent with the scenario that these galaxies continue to form stars (and stellar mass) and build up their number densities at later cosmic times.

Zoom In Zoom Out Reset image size

For all $\mathrm{log}{M}_{* }/{M}_{\odot }\gt 10.5$ galaxies, quiescent systems dominate by number density for all redshifts considered here (z ≲ 1.5). For galaxies more massive than $\mathrm{log}{M}_{* }/{M}_{\odot }\gt 11$, quiescent galaxies at fixed mass outnumber star-forming galaxies by a factor of ∼3:1. Quenching of star-forming galaxies can therefore at most contribute to roughly a 33% increase in the number density of quiescent galaxies at z < 1.5.

We can gain additional insight by comparing our results to other studies. These cumulative number densities are plotted in Figure 20 at the redshifts and stellar masses where they overlap with our study. We begin by comparing our data to that of S82-MGC (Bundy et al. 2017, ∼140 deg²). Because both their study and our study use the forward-modeling method to account for the Eddington bias, we can directly take their best-fit intrinsic SMF and integrate it to compute the cumulative number density, at least to >10^11.2 M_⊙, where Bundy et al. are complete.

Given that our result is in excellent agreement with that of the larger survey area of S82-MGC, this strongly suggests that our 17.5 ${\deg }^{2}$ SHELA survey comoving volume of ∼0.15 Gpc³ in the redshift range of 0.4 < z < 1.5 is sufficient to mitigate the effects of cosmic variance, even for very rare galaxy populations. We can thus put strong constraints on the evolution of the galaxy SMF down to a stellar mass of 10^10.3 M_⊙.

The lack of evolution (≲0.1 dex) seen in our cumulative number density of all galaxies more massive than M_* > 10¹⁰ M_⊙ is also in agreement with the result of Moustakas et al. (2013; Figure 20). By integrating the observed SMF of the 5.5 deg² PRIsm MUlti-object Survey (PRIMUS; Coil et al. 2011), these authors find a ≲10% change in the number density of M_* > 10¹¹ M_⊙ galaxies since z ≈ 1. It is interesting that when Moustakas et al. (2013) divided their sample into quiescent and star-forming galaxies, they found that the number density of quiescent galaxies with 10¹¹–10^11.5M_⊙ has changed relatively little since z = 1, and the decline in the number density of massive star-forming galaxies is ≲0.2 dex. This is in agreement with our finding. We do observe slight offsets in the normalization of the SMF compared to our result, which may be a result of systematics in the sample selection or SED analysis.

In contrast, the cumulative number density evolution observed in the UltraVISTA (Muzzin et al. 2013, 1.62 deg²) survey finds a ∼0.2–0.4 dex growth in the number density of galaxies with stellar mass above 10¹⁰, 10¹¹, and 10^11.5 M_⊙, respectively, from z = 1 to z = 0.2. While on the surface this runs counter to our results, an inspection of Figure 20 shows that we are in agreement (within the uncertainties) for regions where both surveys are complete: higher mass galaxies with $\mathrm{log}{M}_{* }/{M}_{\odot }\gt 11$ and redshifts z > 0.4. At higher masses, $\mathrm{log}{M}_{* }/{M}_{\odot }\gt 11.5$, it is likely that UltraVISTA is limited by the cosmic variance. While at z < 0.4, we have already argued that SHELA may be incomplete (in comparison to Bundy et al.) and this may be true as well for the much smaller-area UltraVISTA survey. Additionally, as we already noted in Section 4.4, our interpretation of the number density and stellar mass density evolution may be impacted by less accurate SFRs at z < 0.5 because the u band does not sample the rest-frame NUV (see Moutard et al. 2016, who perform a similar analysis including UV photometry with the Galaxy Evolution Explorer satellite).

Figure 20 compares the results from our study to those from Capozzi et al. (2017), using galaxies from at 0.1 < z < 1 observed in DES. Capozzi et al. (2017) find that the cumulative number density evolves grows by ∼0.4 (0.3) dex over this redshift range for galaxies with stellar mass above 10¹⁰ (10^10.5 M_⊙). The larger evolution in number density could arise from a difference in survey selection: the Capozzi et al. sample is optical i-band-selected (i < 23 mag), and this could impact their stellar mass completeness limit particularly at lower masses and higher redshifts. This is where we observed the greatest discrepancy in the evolution of the cumulative number density. Additionally, the discrepancy could arise from the different methods used to account for the statistical and systematic uncertainties in stellar mass measurements. Interestingly, the number density of galaxies with M_* > 10^11.5 of Capozzi et al. (2017) show no evolution since z ∼ 1, in good agreement with our finding and the result at z < 0.06 from GAMA (covering 143 deg², Baldry et al. 2012; see also Figure 15 of Capozzi et al. 2017).

Figure 20 also shows number densities from the integrated SMF from BOSS (Maraston et al. 2013) for galaxies with M_* > 10^11.5, out to z ∼ 0.7 where their survey is complete in stellar mass. The lack of evolution between 0.45 < z < 0.7 in the cumulative number density of BOSS galaxies in the highest mass bin (M_* > 10^11.5) is in good agreement with our finding here. There is a slight offset in the normalization of the cumulative number density of Maraston et al. (2013) compared to our results, but that may arise from the systematics in the sample selection or SED analysis.

Finally, our results are broadly in agreement with (smaller-area) surveys selected using deep near-IR data. Mortlock et al. (2015) find little evolution in the characteristic stellar mass (M*) of the SMF from z = 3 to z = 0.3 in their analysis of the combination of deep near-IR data from the Ultra Deep Survey (UDS) and the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS) UDS and GOODS-S fields (see Figure 11). Additionally, Conselice et al. (2007) presented the evolution of massive galaxies at z ∼ 0.4–2 by combining wide and deep NIR imaging from the Palomar telescope with DEEP2 spectroscopy. They found that galaxies with $\mathrm{log}({M}_{* }/{M}_{\odot })\gt 10.5$ exhibit no significant evolution in number density since z < 1, but show an increase in the number density of galaxies with $11\lt \mathrm{log}{M}_{* }/{M}_{\odot }\lt 11.5$ from z ∼ 1 to 1.5 (see their Figure 4). At higher masses, $\mathrm{log}{M}_{* }/{M}_{\odot }\gt 11.5$, it is likely that their sample is limited by the cosmic variance. These results are consistent with our findings.

Taken together, our analysis, combined with results in the literature, paint a picture where the number density of massive galaxies, $\mathrm{log}{M}_{* }/{M}_{\odot }\gtrsim 11$, is roughly constant out to z ∼ 1. Additionally, this evolution is dominated by the number density of quiescent galaxies.

We compute the total stellar mass density by integrating our best-fitting intrinsic (assumption-averaged) SHELA SMF for stellar masses greater than 10⁹ M_⊙ and compare our results with other studies. Again, we are particularly interested in the internal redshift evolution, as there are significant discrepancies in the normalization of the stellar mass density between the different studies.

The stellar mass density we derive from SHELA shows an overall increase from z = 1.5 to z = 0.4. However, most of this evolution occurs before z ∼ 1. This is consistent with our previous statements that the number density evolution at z < 1 is dominated by quiescent galaxies. Figure 21 illustrates this: quiescent galaxies show no measurable growth in stellar mass density at z < 1. Most evolution occurs in star-forming galaxies, which show a continuous increase in the stellar mass density from z = 1.5 to z = 0.4.

Zoom In Zoom Out Reset image size

We compare our stellar mass density with those from UltraVISTA (Muzzin et al. 2013) and PRIMUS (Moustakas et al. 2013). The values from UltraVISTA are measured by integrating the maximum-likelihood Schechter function fits of the observed SMF, down to a limit of 10⁹ M_⊙; those from PRIMUS are derived by integrating the observed SMF and also the best-fit Schechter function down to a limit of 10^9.5 M_⊙. Overall, the evolution of the stellar mass density for all galaxy populations is consistent, with a ∼0.2 dex growth in the stellar mass density from z ∼ 1.5 to z ∼ 1 (but our data provide the most detailed accounting of uncertainties). At lower redshifts, our result and those from UltraVISTA and PRIMUS are generally consistent with ≲0.1 dex growth z ∼ 1 to z ∼ 0.4. However, we note that our measurements encompass the largest area, and include in the error budget the effects from forward modeling, and uncertainties associated with the derivation of stellar masses.

Lastly, we consider how our results would be impacted using different thresholds to separate quiescent from star-forming galaxies. In this study, we have adopted the evolving threshold of $\mathrm{log}(\mathrm{sSFR}/{\mathrm{yr}}^{-1})$ that ranges from −11 at z = 0.4 to −10.2 at z = 1.5 (see Section 4.4 and Figure 5). We tested how our results would change if we used instead a nonevolving (fixed) threshold of log (sSFR/yr⁻¹) of −11 at all redshifts. With the unevolving sSFR selection, we find no significant evolution (≲0.1 dex) in the cumulative number density of massive quiescent galaxies (with $\mathrm{log}{M}_{* }/{M}_{\odot }\gt 11$) from z = 1.0 to z = 0.4, consistent with the behavior derived using the evolving threshold of sSFR (Figure 20). We do observe a difference when we consider the number density of quiescent galaxies down to more moderate masses of $\mathrm{log}({M}_{* }/{M}_{\odot })\gt 10$, which now exhibits a 0.2 dex increase from z = 1.0 to z = 0.4. Of course, this comes with a trade-off in the evolution of star-forming galaxies, which now show ≲0.1 dex evolution in the cumulative number density from z = 1.0 to z = 0.4 in all stellar mass bins.

The choice of sSFR threshold to differentiate between star-forming and quiescent galaxies clearly impacts our interpretation of evolution (and this is true for other studies in the literature). We favor a redshift-dependent sSFR cut because we wish to study galaxies that are star-forming (with a current SFR higher than or equal to their past average) and quiescent (with a current SFR that is much less than their past average). This is only achieved by using a redshift-dependent sSFR for the reason that the SFR–stellar mass relation itself evolves (see Figure 5) and because there is less time available at higher redshift. In contrast, using a fixed sSFR threshold would change the definition of quiescence and would even include galaxies that are still on the star-forming main sequence at z ∼ 0.3. Our results specifically describe the evolution of quiescent and star-forming galaxies as defined here. Using a different selection threshold to identify different populations (and describe the evolution of different populations of galaxies) would impact the interpretation. It is therefore important to take into account the definition of quiescent and star-forming galaxies when comparing results in the literature.