Selection of Massive Evolved Galaxies at 3 ≤ z ≤ 4.5 in the CANDELS Fields

This method uses the rest-frame U, V, and J bands to select the quiescent population, as the rest-frame U − V probes the prominent Balmer break (3646 Å) seen in poststarbursts. The V − J color is used to break the degeneracy between dusty star-forming and quiescent galaxies. We can classify galaxies using a color cut that separates the quiescent and star-forming regions in the UVJ plane. This method has been developed and extensively used for classifying galaxies in photometric surveys (e.g., Labbé et al. 2005; Wuyts et al. 2007; Williams et al. 2009; Arnouts et al. 2013; Whitaker et al. 2013; Barro et al. 2014; Straatman et al. 2014; Fumagalli et al. 2016; Siudek et al. 2017; Fang et al. 2018). However, measuring UVJ colors for galaxies at high redshifts is challenging, making the selection less reliable. A typical method to infer rest-frame UVJ colors is from the best-fit model SED. The model SEDs are generally built assuming a τ model star formation history (SFH), which also becomes less reliable at high redshifts, since the galaxy model needs more time to evolve into the quiescent regions (Merlin et al. 2018) and more generally suffers from other uncertainties associated with SED fits used to infer rest-frame colors. Therefore, the boundaries of the UVJ criteria are modified at different redshift bins (Whitaker et al. 2013). However, this can introduce contamination from the dusty star-forming galaxies, given the similarity of their rest-frame U − V color. Nevertheless, the overlap in the wavelength coverage, including medium-band data, can ameliorate some of these problems when selecting high-z galaxies (e.g., Spitler et al. 2014).

To measure the rest-frame UVJ colors, we employ the prescription described by Wolf et al. (2003) for the COMBO-17 survey, where the rest-frame colors are estimated from the best-fit template SED. Rest-frame UVJ colors can also be measured by interpolating the observed bands that track the redshifted UVJ (e.g., Rudnick et al. 2003; Taylor et al. 2009; Williams et al. 2009). In the latter case, data from many bands outside of the UVJ region of the spectra are left unused, while in constraining the best-fit model SED, these data are beneficial (Brammer et al. 2011). We will further investigate the effect of the SFHs on the UVJ colors from the best-fit SEDs when analyzing the effects of the photometric uncertainties in Section 4.

We combine the selection criteria on the rest-frame UVJ plane employed in Muzzin et al. (2013), Whitaker et al. (2012), and Straatman et al. (2014) in Figure 1. We define the grading scheme by giving a weight to individual galaxies based on their distance from the boundary defined in Table 7. The combination of these selections is done by finding the likelihood that each galaxy is a member of the union of selected populations from different criteria. To find the rest-frame UVJ colors of CANDELS galaxies, we use the Le Phare SED fitting code (Arnouts et al. 1999; Ilbert et al. 2006) with the standard libraries defined in Table 3 (Chartab et al. 2020).

Zoom In Zoom Out Reset image size

In addition to the criteria on the UVJ plane, we impose a mild detection constraint on the observed J band (S/N > 2), which probes blueward of the Balmer break for the highest redshift bins. This reassures that the break lies redward of the J band, reducing the contamination from low-z interlopers. The results are presented in Figure 1.

The most commonly used methods for identifying the population of high-redshift galaxies are variations of the dropout technique, based on the observed colors of galaxies. This uses well-known features in the galaxy SEDs and follows them as they move to redder passbands when the galaxy is redshifted. Examples of this are the Lyman-break features used for the selection of UV-bright Lyman-break galaxies (Steidel & Hamilton 1993; Steidel et al. 1995; Steidel et al. 2003; Reddy et al. 2005; Stark et al. 2010; Bouwens et al. 2011; Bouwens et al. 2014; Roberts-Borsani et al. 2016; Oesch et al. 2016) and evolved systems using Balmer-break features (Cimatti et al. 2002; Roche et al. 2002; Franx et al. 2003; van Dokkum et al. 2003; Daddi et al. 2004, 2007; Reddy et al. 2005; Mobasher et al. 2005; Lane et al. 2007; Rodighiero et al. 2007; Wiklind et al. 2007; Caputi et al. 2012; Nayyeri et al. 2014; Girelli et al. 2019).

This technique uses the fact that magnitudes and colors are sensitive to redshift and the shape of their SEDs, which follows the physical properties of their stellar population and interstellar medium. For example, for poststarburst galaxies, we can use Balmer-break features from the 3648 Å Balmer limit. Therefore, the observed colors can help us find the candidate galaxies directly from photometric measurements using a few photometric bands. For the redshift range of interest here, Balmer-break features redshift toward NIR wavelengths.

We use the observed colors to select candidates in the CANDELS fields with the color selection criteria from Nayyeri et al. (2014) in Figure 2. This selection was initially developed by finding the color criteria in the color–color space for a sample of old (low-z) quiescent, dusty starburst, and poststarburst populations using Bruzual & Charlot (2003a) stellar population synthesis evolutionary tracks and considering intergalactic medium (IGM) absorption following Madau (1995). In the redshift range of interest ($3\leqslant z\leqslant 4.5$), H − K color will constrain the presence of the Balmer break, accompanied by a Y − J or J − H constraint that can discriminate between dusty starburst and quiescent galaxies. Also, we impose a nondetection requirement on the U and B bands, since they fall blueward of the break. By doing this, we reduce the contamination from the low-redshift interlopers. Another likely source of contamination is due to the absence of nebular emission lines. These affect the selection based on broadband photometry because they can mimic a Balmer-break feature, which can mislead the classification scheme (see Nayyeri et al. 2014; Straatman et al. 2014; Merlin et al. 2018). To minimize the effect of these lines, we perform the SED fitting procedure described in Section 3.3 using libraries with and without the nebular line emissions.

Zoom In Zoom Out Reset image size

We assign likelihoods to the selected sources based on their proximity to the selection boundaries to reduce the dependence on the chosen selection boundaries (more details are presented in the Appendix). The final likelihood will be the average likelihood for all of the realizations of a galaxy within its error bars, which we call Balmer Break Galaxy (BBG) likelihood (Section 4.2).

Here we discuss another method used for finding quiescent candidates based on the inferred physical properties measured from their SEDs. The advantage of this method is that it makes use of all of the photometric data available and, therefore, imposes stronger constraints on the selection process. Furthermore, it predicts the physical parameters for each galaxy. The disadvantage is that the method is model-dependent and based on the SFH and extinction used to generate SED templates (Grazian et al. 2007; Fontana et al. 2009; Pacifici et al. 2016; Carnall et al. 2018; Merlin et al. 2018, 2019; Carnall et al. 2019b). We fit template SEDs for the candidates in the range 2.8 < z_phot < 5.4 with masses larger than 10¹⁰ M_⊙ using CANDELS photometric and spectroscopic redshifts (where available) and a catalog of the physical properties. We rely only on the mass measurements, since the stellar masses are less susceptible to the parameters chosen in the template libraries of galaxies, particularly the SFH used for SED fitting (compared to the inferred SFRs; Papovich et al. 2001; Shapley et al. 2001; Wuyts et al. 2009; Muzzin et al. 2009; Mobasher et al. 2015). However, the presence of the nebular lines can affect stellar mass measurements. Hence, we treat the stellar mass as a free parameter (although we made a prior mass selection on the subsample) when fitting the subsample photometric measurements to keep the stellar mass uncertainty of this type limited to the initial selection.

We find the physical properties using the Bayesian Analysis of Galaxies for Physical Inference and Parameter EStimation (Bagpipes), a Bayesian SED fitting code (see Carnall et al. 2018). Bagpipes uses the 2016 version of Bruzual & Charlot (2003a) and Multinest (Feroz & Hobson 2008; Feroz et al. 2009) for a multimodal nested sampling algorithm (Skilling 2006). The multimodal nested sampling algorithm employed is a huge improvement over a simple χ² fit, which is incapable of producing a reliable error estimate. Moreover, the Markov Chain Monte Carlo (MCMC) algorithm employed in some SED fitting procedures for finding the posterior of the physical properties can be problematic when sampling from a multimodal posterior. These include models with a large degeneracy, which is often the case when modeling a complex system such as a galaxy under large photometric uncertainties.

We fixed redshifts to their photo-z and, where available, spec-z values. We built separate model libraries based on three prescriptions for SFHs: exponentially declining, top-hat, and double power-law forms (Behroozi et al. 2013b) assuming a Chabrier (2003) initial mass function (IMF); dust attenuation based on Calzetti et al. (2000); and IGM absorption from Inoue et al. (2014), which is a revised version of the Madau (1995) model. It is known that the presence of the nebular lines can mimic a Balmer break-like feature in broadband photometry. Therefore, for controlling their effects on the physical properties and consequently, the selection made based on them, we run the code with and without the nebular emission (since the target population is expected to have little to no nebular emission). The priors used for the physical parameters in the SED fitting are listed in Table 4. We then find the posterior distribution for each model parameter. Following Carnall et al. (2018), we define ψ_SFR as the ratio of the SFR at any given time (SFR(t)) to the average SFR over the age of a given galaxy ($\langle \mathrm{SFR}(t)\rangle$),

$$\begin{eqnarray}&&{\psi }_{{\rm{SFR}}}=\displaystyle \frac{{\rm{SFR}}(t)}{\langle {\rm{SFR}}\rangle (t)}=\displaystyle \frac{{\rm{SFR}}(t)}{\tfrac{1}{t}{\int }_{0}^{t}{\rm{SFR}}(t^{\prime} ){dt}^{\prime} },\end{eqnarray} \tag{ 1 }$$

and we define the quiescent population as those with ${\psi }_{\mathrm{SFR}}$ less than 0.1 at the observed age of the galaxy. In other words, the quiescent galaxies are those with an average SFR over the last 100 Myr of less than 10% of the average SFR over its lifetime.

Table 4.  Priors Used for Different Free Physical Parameters and the Fixed Parameters Used in the Fit

SFH	Free Parameter	Prior	Limits	Fixed Parameter	Value
	AV1	Uniform	(0, 2)	SPS models	BC03
Double power law	${\mathrm{log}}_{10}({M}_{\mathrm{formed}}\ /\ {{\rm{M}}}_{\odot })$ 2	Uniform	(1, 13)	IMF	Chabrier (2003)
	${ \mathcal Z }\ /\ {{ \mathcal Z }}_{\odot }$ 3	Uniform	(0.2, 2.5)	zobs	${z}_{\mathrm{phot}/\mathrm{spec}}$
${\rm{SFR}}{(t)=C{[(t/\tau )}^{\alpha }+{(t/\tau )}^{-\beta }]}^{-1}$	$\tau \ /\ \mathrm{Gyr}$ 4	Uniform	(0, $t({z}_{\mathrm{obs}}))$	${\mathrm{log}}_{10}(U)$ 5	−3
	α6	Logarithmic	(10−2, 103)
	β7	Logarithmic	(10−2, 103)
	AV	Uniform	(0, 2)	SPS models	BC0312
Exponentially declining	${\mathrm{log}}_{10}({M}_{\mathrm{formed}}\ /\ {{\rm{M}}}_{\odot })$	Uniform	(1, 13)	IMF	Chabrier (2003)
	${ \mathcal Z }\ /\ {{ \mathcal Z }}_{\odot }$	Uniform	(0.2, 2.5)	zobs	${z}_{\mathrm{phot}/\mathrm{spec}}$
${\rm{SFR}}(t)=C\exp (-t/\tau )$	$\tau \ /\ \mathrm{Gyr}$ 8	Uniform	(0.05, 10)	${\mathrm{log}}_{10}(U)$	−3
	Age	Uniform	(0, $t({z}_{\mathrm{obs}}))$
	AV	Uniform	(0, 2)	SPS models	BC03
Top hat	${\mathrm{log}}_{10}({M}_{\mathrm{formed}}\ /\ {{\rm{M}}}_{\odot })$	Uniform	(1, 13)	IMF	Chabrier (2003)
	${ \mathcal Z }\ /\ {{ \mathcal Z }}_{\odot }$	Uniform	(0.2, 2.5)	zobs	${z}_{\mathrm{phot}/\mathrm{spec}}$
${\rm{SFR}}(t)=\left\{\begin{array}{lc}C & t\in [{\mathrm{Age}}_{\min },{\mathrm{Age}}_{\max }]\\ 0 & \mathrm{otherwise}\end{array}\right.$	Age${}_{\min }$ 9	Uniform	(0, $t({z}_{\mathrm{obs}}))$	${\mathrm{log}}_{10}(U)$	−3
	Age${}_{\max }$ 10	Uniform	(0, $t({z}_{\mathrm{obs}})$)

Notes.

¹A_v is the attenuation at 5500 Å. ²M_formed is the mass formed. ³ ${ \mathcal Z }$ is the metallicity. ⁴τ is the peak time for double power-law SFH. ⁵ ${\mathrm{log}}_{10}U$ is the ionization parameter. ⁶α is the rising power. ⁷β is the falling power. ⁸τ is the exponential decay timescale in τ model SFH. ⁹Age_min is the initial time for the top-hat SFH. ¹⁰Age_max is the final time for the top-hat SFH, and 11C is the normalization constant. ¹²BC03 is the stellar population synthesis at the resolution of 2003 (Bruzual & Charlot 2003b).

Download table as:  ASCIITypeset image

Carnall et al. (2018) showed that the selection criterion mentioned above is consistent with the definition proposed in Pacifici et al. (2016; see Figure 4), which is a criterion on the sSFR evolving with redshift to define the quiescent population at a given age of the universe.

The Bayesian nature of the SED fitting code allows us to have the posterior distribution for each galaxy and all parameters associated with it. Then we can apply the selection on the sample from the posterior and count the number of selected samples from the posterior versus the total count. We assign this ratio as the probability of being selected given the posterior sample. We apply the same prescription for the rest of the models employed.

We use the Bayesian evidence (marginal likelihood) on different models used for SED fitting. By doing so, we find the relative evidence for different models according to the data. The definition of the Bayesian evidence is

$$\begin{eqnarray}&&\zeta =p(D| H)=\int p(D| \theta ,H)p(\theta | H)d\theta ,\end{eqnarray} \tag{ 2 }$$

where D, H, and θ represent the data, hypothesis (model), and model parameters, respectively. The $p(D| H)$ is the probability of getting the data D under the assumption of the validity of a specific hypothesis H, which is the Bayesian evidence for H.

We adopt the Jeffreys criterion (Jeffreys 1961; also see Kass & Raftery 1995) for interpretation of the relative evidence (Bayes factor) calculated for each galaxy and different models. We use the same prior for fixed parameters such as redshift, and for the rest of the parameters, we use noninformative priors (Kass & Wasserman 1996) across different models. By doing so, the Bayes factor represents the posterior odds of the models, ($\displaystyle \frac{p({H}_{2}| D)}{p({H}_{1}| D)}=\tfrac{p(D| {H}_{1})}{p(D| {H}_{2})}$). In other words, we impose no prior judgment about the validity of a certain model. Results presented in Figure 3 express that statistically; none of the models show substantially more evidence. However, the Bayes factors support models that include nebular lines and generally favor double power-law over τ models. This is consistent with Reddy et al. (2012) regarding the inability of the τ models in reproducing SFRs from UV+IR measurements in the star-forming population at higher redshifts. Pacifici et al. (2016) showed that the double power law is the best model to describe the SFHs of their samples. Carnall et al. (2019a) also confirmed that a double power-law model could produce relatively stronger evidence. As Figure 3 depicts, similar to Belli et al.'s (2019) investigation of the effect of SFH on age measurements, top-hat models generally fail to produce a comparable level of evidence from data. Therefore, we do not use the results from the top-hat models for making a selection.

Zoom In Zoom Out Reset image size

We cross-matched the potential candidates for high-redshift massive evolved systems with the publicly available Chandra X-ray catalogs and used the Spitzer MIPS 24 μm to detect possible dusty active galactic nuclei (AGNs). We excluded any candidates with a counterpart (shows detection) in either the X-ray or MIPS. For the X-ray catalogs, we looked for any counterparts closer than 5'', which is 10 times the resolution of the Chandra X-ray Observatory (cross-matched with Laird et al. 2009; Evans et al. 2010; Nandra et al. 2015; Xue et al. 2016; Cappelluti et al. 2016; Civano et al. 2016; Masini et al. 2018), and for MIR, we used the catalog published in Barro et al. (2019) for finding any counterpart and flux measurements in the mentioned MIR bands. About 30% of the candidates have X-ray or MIR counterparts, which were removed from the final sample.

To each galaxy, we assigned a likelihood measure that is the median of the likelihoods estimated for that galaxy from each of the three techniques: the UVJ selection (Section 3.1), observed color selection (Section 3.2), and SED fitting method with and without the nebular emission contribution (Section 3.3). We combine the results from the SED fitting under different assumptions using the weighted average of the likelihoods for each galaxy using the marginal likelihood calculated for a particular model as their corresponding weights. By doing this, we make sure that we have put more emphasis on likelihoods calculated from the models with higher marginal likelihoods. Then we take the median of all of these methods as our final indicator. We select the final sample as those with a median likelihood higher than 0.5. However, we use the 0.3 and 0.7 criteria as our least and most conservative samples, respectively. By using the median indicator, we limit our sample to those galaxies that were assigned by at least two of the methods mentioned to have a high likelihood of being a massive evolved galaxy. Figure 8 shows how this threshold changes the number density measurements of each selection method, as well as the median value that is taken as the final indicator. The final selected galaxies with their assigned likelihoods and estimated physical parameters are listed in Table 5. The galaxies that were selected with a median likelihood higher than 0.3 constitute the most inclusive sample and the one used for finding the upper bounds on the number/mass densities. Also, we control the possible contamination in our final sample, since we rely on the composite indicator compared to a single measure. Figure 4 shows a comparison between different selection methods and how they relate to each other (all highlighted points are galaxies in our final sample). For example, the selection in the sSFR versus M_s plane is generally consistent across different methods, since the objects with higher BBG/UVJ likelihood are close to or inside the selection boundary. For selection in the UVJ colors, we have several candidates that are far from the boundary but show a much higher BBG/SED likelihood. This shows the sensitivity of the final results to the libraries used, and by changing the SED fitting code and/or the libraries used, we are not necessarily searching through the same part of the models' parameter space. In terms of the selection based on the observed colors, we tend to have a higher UVJ/BBG likelihood, but there are a couple of candidates that show higher UVJ/SED likelihoods but fall outside of the criteria. This shows the importance of using information from other bands as well.

Zoom In Zoom Out Reset image size

Here we consider the effect of the photometric uncertainties in finding UVJ colors and how that could affect the final selected sample. We first choose a subset of the selected galaxies, construct 50 realizations of each, and run SED fitting on the resampled galaxies. Then we quantify how these uncertainties propagate into our inferred rest-frame colors and how they change the result from the UVJ selection in Section 3.1.

Figure 5 shows the effect of photometric uncertainty on the UVJ selection. It is clear how the choice of the SFH and whether to include the nebular emissions is pivotal in selecting the quiescent population. We find the UVJ colors from the Bagpipes SED fitting code, since we are interested in understanding the effect of the photometric uncertainty on the inference from the SED, which makes a Bayesian posterior estimate more reliable given the clear uncertainty definition based on the posterior. For all but a small subset of galaxies fitted with τ models and nebular lines, the posterior UVJ of the resampled galaxies falls around the posterior from the true photometry with considerable scatter. For a small subset of galaxies (five out of 28) fitted with τ models, the locus of the posterior varies dramatically and with higher scatter than the posterior UVJ from the true photometry. The double power-law model is found to be more immune to the photometric uncertainties, as is also the case for the physical parameters used in SED-based selection.

Zoom In Zoom Out Reset image size

To find the effect of photometric uncertainties on the selection based on the observed colors, we generate 10⁵ realizations for every galaxy and assign a likelihood based on the number of realizations that fall into the hard selection criteria to the total number of realizations. The results for the GOODS-South galaxies are shown in Figure 6, where they are colored and resized based on the assigned likelihood. Few candidates fall outside of the selection criteria, but when we include the photometric uncertainties, we find realizations that fall in the selected region. On the other hand, few galaxies fall inside the criteria, but because of their photometric uncertainties, they do not show a clear indication of being quiescent.

Zoom In Zoom Out Reset image size

With having a posterior distribution of every physical parameter, we can estimate the uncertainties of these parameters. We use the resampling method mentioned above to get even more reliable error estimates (more details on the effect of the sampling error can be found in Higson et al. 2018).

Here we used the same sample from the SED fitting results for finding the sensitivity of the UVJ colors on photometric uncertainties. We follow the effect on the posterior ψ_SFR, which is the critical parameter in the SED-based selection discussed in Section 3.3. As Figure 7 shows, for most of the sample used, the τ-model SFH is more susceptible to photometric uncertainties than the double power law, which has more stable results under photometric perturbations similar to the UVJ colors. Also, including the models' nebular emission leads to models with UVJ colors that are further away from the quiescent region and hence reduce the number densities of the quiescent population.

Zoom In Zoom Out Reset image size

In this section, we found the effect of the photometric uncertainties on each method's assigned likelihood. As mentioned in Section 3.2, for finding the likelihood based on the observed colors for all of the galaxies in the CANDELS fields, we assign a likelihood to every photometric realization of a galaxy based on their proximity to the selection boundary, then we take the average of these as the observed color likelihood for that galaxy. By doing this, we can take into account both the photometric uncertainties and the closeness to the selection boundary. Since we only calculated the effect of the photometric uncertainties on likelihoods based on the UVJ- and SED-based selection for a subsample of galaxies, we did not incorporate the photometric uncertainties in the final assigned likelihood of these methods.

To estimate the completeness of the sample used, we follow the prescription adopted in Pozzetti et al. (2010). We divide the CANDELS galaxies into the two redshift bins mentioned above. We find the 20% of the galaxies with the faintest apparent magnitude in the H band, with sSFRs lower than 75% of the galaxies in that redshift bin. We calculate the M_lim, defined as the mass that the galaxy would have if it were at the survey limiting magnitude (H_lim = 26; assuming a constant mass-to-light ratio); in other words, we find $\mathrm{log}({M}_{\mathrm{lim}}/{M}_{* })=0.4(H-{H}_{\mathrm{lim}})$. Then we define the fraction of the sample with $\mathrm{log}({M}_{\mathrm{lim}}/{M}_{\odot })\gt 10$ as the completeness of the quiescent galaxies with that mass range. The effective volume of the survey at each redshift bin is defined to be V_eff = V_co η(z), where V_co is the comoving volume at each redshift bin and η(z) is the completeness of the sample for $\mathrm{log}({M}_{* }/{M}_{\odot })\gt 10$ at each redshift bin of the quiescent population.

In this section, we find the target population predicted from the several existing models listed below and compare their number and stellar mass densities with the sample found in this study.

• 1.  
  Simulated infrared dusty extragalactic sky. We use the simulated catalog presented in Béthermin et al. (2017), which uses the abundance matching technique for occupying dark matter halos from the Bolshoi–Planck simulation (Klypin et al. 2016; Rodríguez-Puebla et al. 2016) and an updated version of the two star formation modes (2SFM) galaxy evolution model (Sargent et al. 2012; Béthermin et al. 2012) from which a light cone covering 2 deg² was produced.
• 2.  
  Universe machine. We employ the available CANDELS light cones produced by the universe machine semiempirical model, which uses the same Bolshoi–Planck simulation for halo properties and mass assembly histories. There are eight realizations for each field in CANDELS, and we use all of these realizations (Behroozi et al. 2019).
• 3.  
  Eagle simulation. We use the Eagle cosmological hydrodynamical simulations of a box with L = 50 and 100 comoving megaparsecs (cMpc) on each side. We use RefL0050N0752, RefL0100N1504, and AGNdT9L0050N0752, in which the first two are the reference physical model used and the last one has a higher AGN heating temperature and lower subgrid black hole accretion disk viscosity. All of these simulations have the same mass resolutions (Crain et al. 2015; Schaye et al. 2015; McAlpine et al. 2016; The EAGLE team 2017).
• 4.  
  IllustrisTNG simulation. We use the latest IllustrisTNG cosmological hydrodynamical simulations of a box with L = 100 and 300 cMpc on each side. There are three resolutions for each simulation box, and here we used the highest resolution (Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2018, 2019; Pillepich et al. 2018; Springel et al. 2018).

To identify the quiescent galaxies in the simulations and be consistent with the observation of the massive evolved galaxies at z > 2, which have been known to be quite compact, ${R}_{e}\leqslant 2.5\,\mathrm{kpc}$ (Daddi et al. 2005; Trujillo et al. 2006; Trujillo et al. 2007; van Dokkum et al. 2008; Newman et al. 2010; Damjanov et al. 2011; van der Wel et al. 2014), we use the stellar mass and SFRs within twice the half-mass radius. Donnari et al. (2019) showed that the star formation main sequence in the IllustrisTNG100 for z > 2 is quite identical when assuming physical properties within an aperture of 5 kpc compare to the twice the half-mass radius definition for SFR measurements averaged over 10 and 1000 Myr. Also, following Merlin et al. (2019) and staying consistent in our comparison within the simulations, and since the physical properties of the Eagle simulations are available within certain apertures in the range 1–100 kpc, we use the physical properties within an aperture that is closest to four times the half-mass radius. For all of the models, we select galaxies with stellar masses $\mathrm{log}({M}_{* }/{M}_{\odot })\gt 10$ and sSFRs within the half-mass radius lower than the evolving sSFR constraint employed in Pacifici et al. (2016), as sSFR_lim = 0.2/t_U(z), where t_U(z) in Gyr is the age of the universe at redshift z. Carnall et al. (2019a) showed this to be consistent with the ψ_SFR measure used in this work across different redshifts.

Figure 9 shows the comparison between the observed number density of the quiescent galaxies within the corresponding comoving volumes, corrected for their completeness. The observed and predicted counts from the numerical simulations and empirical models of galaxy evolution are in generally good agreement. However, at the higher end of the redshift distribution, there are more observed candidates than predicted by models by a factor of ∼5. Nevertheless, it is somewhat consistent with the uncertainties shown in Figure 9 (especially in the case of the universe machine). Additionally, the presence of possible contaminants in the sample can ameliorate the tension.

Zoom In Zoom Out Reset image size

Figure 9 further shows the stellar mass density for our sample compared to the simulations and empirical models. We used the mass catalog from CANDELS to measure the stellar mass density. The tension between the stellar mass density for our sample and the simulations seems to be worse; however, the argument above for the number density applies here as well. The stellar mass measurements from the SED fitting results when including nebular line emission can also mitigate the tension, as the presence of the nebular lines can reduce the contribution from the continuum part of the SED, which can result in lower inferred stellar masses. Also, looking at the mass distribution of the selected samples from models and our CANDELS candidates (Figure 10) shows that the CANDELS candidates are relatively more massive than their counterparts in models, except for those from the TNG simulations. The TNG100-1 sample is shown to be more massive, and, using the fact that the stellar mass density at the highest redshift bins is lower than what we found here, we argue that the sample galaxies from the TNG tend to build their stellar mass later and at a higher rate than what is suggested from the CANDELS sample, similar to what Fontana et al. (2009) found. Also, by looking at the merger tree information of the selected sample (Rodriguez-Gomez et al. 2015) and following their history back to their formation epoch, we find the mass assembly and SFHs of the selected sample to check whether it is consistent with their selection as quiescent galaxies. Figure 11 shows the median and 20th and 80th percentiles of the mass assembly and SFR history of the sample selected. The SFHs were calculated from two SFR measures, one within the half-mass radius and one within twice the half-mass radius. Figure 11 reveals that the TNG sample at the lowest redshift bin is consistent with quenching even within twice the half-mass radius. Although the higher-redshift sample shows modest evidence of quenching within the half-mass radius, the star formation seems to continue outside of the half-mass radius.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size