BULGE GROWTH AND QUENCHING SINCE z = 2.5 IN CANDELS/3D-HST

The HST CANDELS observations used for this study comprise high-resolution imaging in five distinct fields: GOODS-South, GOODS-North, COSMOS, UDS, and EGS, covering a total area of 625 arcmin². Typical limiting depths in H₁₆₀ are 27.0 mag for CANDELS/wide and 27.7 mag for CANDELS/deep (the central halves of the GOODS fields), respectively. Additional data used for this work include pre-existing Advanced Camera for Surveys imaging in the GOODS, EGS, and COSMOS fields (Giavalisco et al. 2004; Davis et al. 2007; Koekemoer et al. 2007), plus WFC3 imaging in the Early Science Release (ERS; Windhorst et al. 2011). The available passbands are B₄₃₅, V₆₀₆, i₇₇₅, z₈₅₀, J₁₂₅, H₁₆₀ for the GOODS fields, and V₆₀₆, I₈₁₄, J₁₂₅, H₁₆₀ for the remaining fields. Additionally, Y₀₉₈ imaging was available for ERS as part of GOODS-South and Y₁₀₅ for the regions with CANDELS/Deep coverage.

All imaging used in our analysis was drizzled to a 006 pixel scale. For details on the observations and data reduction, we refer the reader to Koekemoer et al. (2011) and Grogin et al. (2011).

The main focus of this paper is on the structural shape of stellar mass distributions in a mass-selected sample of galaxies in the redshift range 0.5 < z < 2.5. To this end, we apply a mass cut M > 10¹⁰ M_☉ to select our galaxies, which is well above the mass completeness limit of the H₁₆₀-selected catalogs for all five fields. Our sample consequently spans a wide range of SFRs—from normal SFGs on the MS to QGs below—that already formed the bulk of their stars. For selection, we use the SED modeled galaxy-integrated mass estimates, but note that they nicely match the masses obtained by summing the resolved stellar mass distributions (see also Wuyts et al. 2012).

In order to determine redshifts for the galaxies in our sample, we use ground-based spectroscopic redshift information whenever available. Otherwise, redshifts are fitted to the combination of 3D-HST grism data and broadband photometry.

The total sample comprises 6764 galaxies, of which 3839 and 2925 lie within the redshift range 0.5 < z < 1.5 and 1.5 < z < 2.5, respectively. In the following, these two redshift ranges are referred to as z ∼ 1 and z ∼ 2, respectively.

A detailed description of the resolved SED modeling can be found in Wuyts et al. (2012). Here, we review only the key steps involved and additional processing steps.

First, all images in the available wavelength bands are brought to the WFC3 H₁₆₀ resolution (018) by using the IRAF PSFMATCH algorithm. Next, a Voronoi-binning scheme from Cappellari & Copin (2003) is applied in order to ensure a minimum signal-to-noise ratio level of 10 for each bin in the corresponding H-band image. The binned multi-wavelength images are then fit with stellar population synthesis models from Bruzual & Charlot (2003), assuming a Chabrier (2003) IMF, a uniform solar metallicity, and a Calzetti et al. (2000) reddening law with visual extinctions in the range 0 < A_v < 4. The adopted SFHs are exponentially declining, allowing for e-folding timescales down to 300 Myr. This choice is simplistic and may not be representative of the real SFHs (e.g., for SFGs, see Maraston et al. 2010). However, we are interested in the stellar mass, which, among SED-derived properties, is most robust against variations (and uncertainties) in model assumptions (e.g., Papovich et al. 2001; Shapley et al. 2005; Maraston et al. 2010).

In order to conduct further structural analysis using the (binned) output of this SED modeling technique, we constructed pixelized (instead of Voronoi binned) images of the stellar mass surface density distribution. First, we construct a M/L map that, within the galaxy's H₁₆₀ segmentation map, has uniform values for pixels belonging to the same Voronoi bins. Pixels outside the segmentation map (i.e., containing sky noise and possibly faint wings of the galaxy extending below the signal-to-noise threshold) are assigned the average M/L of the nearest three Voronoi bins. The resulting expanded M/L map is then combined with the H₁₆₀ image to construct a final mass map at full H₁₆₀ resolution. This ensures a smooth transition of the galaxy's mass profile from the brighter central parts to the faintest regions.

In order to conduct a structural analysis of our galaxy sample, we employ the GALFIT (Peng et al. 2010) morphology fitting code to fit two-dimensional (2D) parametric models to the stellar mass distribution. For comparison purposes, we model the 2D H₁₆₀ surface brightness distributions as well, following identical procedures.

As part of our set of models, we use single Sérsic models to parameterize the galaxy's shape with the effective radius (R_e) and Sérsic index n. R_e is the galactocentric distance containing half the total light/mass, and n is a measure of the cuspiness of the overall light/mass profile. We allow n to vary within the range 0.2 < n < 8.

While the Sérsic index n is often used and indeed serves as an approximate measure of the contribution of the bulge, it is important to note that n does not translate one to one to the bulge-to-total ratio, B/T (see Appendix A; also Andredakis et al. 1995; de Jong et al. 2004; Cibinel et al. 2013; Bruce et al. 2012; Bruce et al. 2014). For a given bulge + disk composite profile, the best-fit Sérsic index can be increased both by increasing the B/T, and by leaving the B/T constant while growing the extent of the disk relative to that of the bulge (i.e., lowering R_{e, B}/R_{e, D}). For this reason, and to facilitate a more direct comparison to (semi-analytic) models that are parameterized in terms of bulge and disk units, we furthermore decompose the galaxies by fitting two-component (bulge + disk) models. In addition, we also performed one-component pure disk and pure bulge fits by forcing n = 1 and n = 4, respectively.

For the bulge-to-disk decomposition, we adopt a procedure similar to that implemented in Bruce et al. (2012) for the bulge+disk decomposition of H₁₆₀ light profiles, where we fixed the Sérsic indices to be n = 1 for the disk and n = 4 for the bulge.¹⁵ The centers are left free, but we restrict the relative distance between the bulge and disk centers to be less than 2 pixels. All other parameters defining the two components (R_e, the axial ratio b/a, and the total magnitude/mass of both components) are allowed to vary independently.

In two-component modeling, the higher number of degrees of freedom increases the odds of the fit being trapped in a local χ² minimum. In order to mitigate this risk, we ran GALFIT using a grid of initial starting values. Our grid was constructed by using a range of size ratios between bulge and disk (R_{e, B}/R_{e, D}), ranging from 0.1 to 1, in steps of 0.1. For each initial guess of R_{e, B}/R_{e, D}, the corresponding initial guess on B/T was then set, such that the Sérsic index matching this initial configuration matches the one measured in the single-component fit (see Appendix A). Likewise, the initial magnitudes and absolute values of the initial size guesses for bulge and disk were set such that the total magnitude and half-light/mass radius of the composite profile matches the respective values determined from Sérsic fitting.

In cases where GALFIT yields solutions with implausibly small bulge sizes (<0.1 pixel, corresponding to ≲ 1/30 of the resolution, or ≲50 pc at z ∼ 2) and flags the outcome as potentially not converged and unphysical, we excluded that run from the grid. Also, solutions yielding a disk smaller in size than the bulge (i.e., R_{e, B}/R_{e, D} > 1) were not included.¹⁶

After performing the fits for each point of the grid with initial guesses, we assigned a final B/T ratio for each object as the solution of the fit with the lowest $\chi _{{\rm red}}^2$. Here, the pure disk and pure bulge fits were also included. Their solutions generally show higher $\chi _{{\rm red}}^2$ than the two-component fits, but are occasionally preferred over those in a $\chi _{{\rm red}}^2$ sense (oftentimes, these are bulgeless systems with n < 1). The two-component decompositions are statistically preferred over the single Sérsic models (as based on both their $\chi _{{\rm red}}^2$ values and the Akaike information criterion, AIC¹⁷) for ∼2/3 of the total sample. Those systems for which the single Sérsic model yields a lower value of $\chi _{{\rm red}}^2$ and AIC typically feature shallow profiles with n < 1.

In both single- and two-component fitting, we use an automated scheme that pre-determines neighboring sources that need to be masked or fitted simultaneously, and passes initial guesses of fitting parameters to GALFIT. Those include estimates on size, total magnitude, and center. The initial guess on size is based on the distance from the center to the radius at which the curve of growth reaches 50% of the galaxy's total flux. The mass-weighted center of the galaxy derived within its segmentation map is adopted as initial estimate for the center. GALFIT takes into account the convolution of the model with the point-spread function (PSF). Both for fitting light and mass, we use a PSF that is a combination of stacked stars and a Tiny-Tim (Krist 1995) model PSF. For a more detailed description of the used PSF, see van der Wel et al. (2012).

We emphasize that throughout this paper our working definition adopted for bulge and disk components is based on the above bulge+disk modeling of stellar mass or light maps, as there are currently no empirical constraints on whether or not stars assigned to a bulge/disk component are dynamically hot/cold. As the SINS survey of z ∼ 2 galaxy kinematics demonstrated that high-redshift SFGs are dynamically (Förster Schreiber et al. 2009) and morphologically (Förster Schreiber et al. 2011) distinct from local spiral galaxies, we caution that bulge fractions as derived by our decompositions may have a somewhat different meaning than they would have in the local universe.

In addition, we estimated a typical measurement uncertainty on B/T by setting up an array of model galaxies with a range of mass, size, B/T ratio, and R_{e, B}/R_{e, D}, which is similar to our data. The grid of model galaxies consisted of 5, 3, 11, and 10 grid points in mass, R_e, B/T, and R_{e, B}/R_{e, D}, respectively. These were then inserted in multiple empty sky regions of the CANDELS UDS field to mimic the typical background noise.¹⁸ We next ran GALFIT on all 8250 mock galaxies using our two-component fitting scheme. The measurement error on B/T of each galaxy in our sample given its mass, radius, and profile shape is then finally assigned as the scatter among the recovered B/T ratios of the corresponding model galaxy.

Typical measurement errors in B/T for star-forming and QGs are, on average, ∼0.05 and ∼0.06 at z ∼ 1 and ∼0.1 and ∼0.13 at z ∼ 2, respectively. The distribution of errors peaks below the median (≲ 0.05), and shows a tail toward higher B/T errors. Two alternative methods to estimate the uncertainty in B/T, namely, the formal random uncertainty reported by GALFIT and re-fitting the observed galaxies after applying additional background and Poisson noise, generally lead to lower estimated uncertainties (by a factor of ∼2 in the case of GALFIT). In the remainder of the paper, we therefore adopt the most conservative error estimates inferred from our analysis of the inserted mock galaxies.

When exploiting the bulge–disk decompositions of the stellar mass maps derived for our sample of massive galaxies, we first evaluate the average mass budget of disks and bulges. The following is an example of picking a random star out of the sample of massive galaxies above 10¹⁰ M_☉. At 1.5 < z < 2.5, the probability that this star belongs to a bulge component is 46%. Increasing the mass limit to log (M) = 10.5 or 11 yields a higher probability for the star to be associated to the bulge (49% and 54%, respectively). Perhaps somewhat surprisingly, the fraction of stars residing in a bulge component rises only slightly, to 0.5 < z < 1.5 (47%, 50%, and 56%) for galaxies more massive than log (M) = 10, 10.5, and 11, respectively.

The formal uncertainties to the above stated probabilities, including sample variance and typical measurement errors on B/T, are limited to a few percent. The total error budget is likely dominated by systematics, for example, related to the assumptions made in stellar population modeling (see Section 3.1). We note, however, that only M/L uncertainties with a differential impact on bulges and disks will affect the above numbers. Even if the M/L ratio of bulges were systematically under- or overestimated by 0.2 (0.3) dex with respect to those of disks, the change would be limited to ∼7(10)%.

Since bulges, unlike stellar disks, can be considered sinks in the continuous assembly of a galaxy's stellar component, the rising mass density of stars in bulges (by a factor of ∼1.8 from the higher to the lower redshift bin) seems to be compensated largely by the continuing assembly of new disks. Splitting our sample in finer redshift bins, we do find the fraction of stellar mass in bulges to increase more significantly, by a factor ∼1.5 over the entire 6 Gyr timespan sampled in our study.

Overall, the bulge mass fractions are higher than what would be inferred from fits to the H-band surface brightness profiles, as the median mass-to-light ratio of disk components is 0.2 dex lower than the bulge components. The above numbers address the evolving mass budget of disks and bulges for a mass-limited sample, including both star-forming and QGs. In the remainder of the paper, we will break down our sample by star formation activity.

In recent years, several HST-based studies have investigated the structural differences between star-forming and QGs at high redshift. In all of these analyses star-forming systems are found to have significantly larger rest-optical sizes than their quiescent counterparts at the same mass and redshift (e.g., Toft et al. 2009; van der Wel et al. 2014). In addition, their surface brightness profile shapes tend to be shallow (n ∼ 1), whereas QGs feature cuspier light profiles (e.g., Wuyts et al. 2011; Bell et al. 2012; Cheung et al. 2012, although see also Bruce et al. 2012 for decomposed light profiles at high stellar masses and redshifts).

At the same time, the same multi-wavelength, high-resolution look-back surveys have also established that substantial mass-to-light ratio variations in the rest-optical can occur, not only between but also within galaxies (see, e.g., Wuyts et al. 2012; Guo et al. 2012; S. Boada et al. in preparation). Since new stars typically tend to form from gas settled in a disk configuration, M/L ratio variations are anticipated to give rise to a composite light profile in which the disk component has a relatively larger weight (per unit mass) than the bulge. It is therefore important to address to what degree the above structural distinction between the two classes of galaxies is intrinsic to their distribution of stellar mass, or, conversely, can be attributed to stellar population effects. The answer is of immediate relevance to our understanding of quenching, as the latter scenario could imply that compact quiescent systems can evolve from the star-forming MS by simple fading, without invoking an associated morphological transition. Kriek et al. (2009) investigate this scenario for a spectroscopically confirmed sample of massive z ∼ 2 galaxies, finding that three out of six massive star-forming systems have dense cores, and thus may passively evolve into compact galaxies due to the fading of the outer star-forming regions. Szomoru et al. (2010), on the other hand, exploit the exquisite depth of the Hubble Ultra Deep Field to probe the surface brightness profile of a massive compact QG at z = 1.91, ruling out the existence of a faint extended envelope or disk around the observed galaxy. Another argument against fading comes from Cheung et al. (2012), who derived the stellar masses of bulges for both star-forming and QGs at 0.5 ⩽ z < 0.8. They found that the bulges of SFGs are half as massive as those of similar-mass QGs, implying they cannot simply fade onto the red sequence without structural evolution.

Using the stellar mass maps reconstructed for our mass-selected sample of 6764 galaxies at 0.5 < z < 2.5 with log (M) > 10, we are now able to draw statistically significant conclusions on the structural distinction between high-z galaxies prior to and after quenching. In Figure 1, we compare the shape of the stellar mass distributions (i.e., corrected for spatial M/L variations) of star-forming and QGs, and study their dependence on the total galaxy stellar mass. We consider profile parameters based on single-component (i.e., Sérsic) fits as well as two-component (bulge + disk) decompositions, and show the results for two separate redshift intervals: 0.5 < z < 1.5 and 1.5 < z < 2.5. In both cases, we identified a galaxy as quiescent if its specific SFR (sSFR) satisfied sSFR < 1/3t_Hubble, where t_Hubble is the Hubble time at the redshift of the galaxy, and as star-forming otherwise. We tested that a definition of quiescence based on the location of a galaxy in the UVJ diagnostic diagram (Wuyts et al. 2007; Williams et al. 2009) yields effectively identical results.

Zoom In Zoom Out Reset image size

Figure 1 immediately highlights that the distinct structural appearance of star-forming and QGs is intrinsic to its internal distribution of stellar mass, and not just driven by stellar population or obscuration effects. In fact, a comparison to the equivalent plots based on H-band surface brightness profiles rather than mass maps (see Appendix B) indicates that stellar population effects (when measuring at rest-optical wavelengths) only induce a modest, albeit non-negligible shift. At all masses, QGs feature cuspier stellar mass distributions (i.e., higher n) than star-forming systems. Their typical best-fit Sérsic index is also an increasing function of galaxy mass. Interestingly, among SFGs the profile shape is not independent of stellar mass. An increase in n is apparent above 10¹¹ M_☉, both at z ∼ 1 and at z ∼ 2. A similar trend of increasing cuspiness at the tip of the MS was noted by Wuyts et al. (2011, Figure 1; see also E. J. Nelson et al. in preparation).

Next, it is worthwhile to reflect on what is measured when fitting Sérsic profiles. Appendix A illustrates that, when considering galaxies as superpositions of bulge and disk components, a given best-fit Sérsic index does not necessarily correspond one-to-one to a unique B/T value, even though it is often interpreted as such. Given a bulge+disk system with an associated best-fit n, one can increase its n either by boosting B/T, or, alternatively, by growing the extent of the disk with respect to that of the bulge without any change to B/T.

The right-hand panel of Figure 1 explores the B/T ratio as a function of galaxy mass, for SFGs and QGs separately. Again, we find a clear anti-correlation between star-formation activity and bulge prominence. Focusing on the star-forming population, the median B/T is limited to below 30% for intermediate mass SFGs (10 < log (M) < 11), while typical bulge mass fractions rise to 40%–50% above 10¹¹ M_☉. We note that there is a significant scatter in the distribution of individual B/T values around the median for both quiescent and SFGs. We investigated the variation in median trends when varying the binning intervals, finding negligible changes at lower masses, while the median B/T of the most massive (log (M) ∼ 11.3) SFG bin changes by ±0.1, depending on the applied binning intervals.¹⁹

We note that measurements on the H band yield bulge mass fractions among SFGs that are lower by ∼30% (on average), as can be understood from a disk component composed of a younger, lower M/L stellar population than the bulge.

From the two-component fits, we infer a typical R_{e, B}/R_{e, D} size ratio of ∼0.2, albeit with significant scatter (see Figure 6). The median size ratio shows little dependence on star-formation activity or mass, over the range probed by our sample. Given the enhanced B/T values in QGs, and the fact that bulges have smaller half-mass radii than disks, one might wonder if the difference in total size between SFGs and QGs can be completely accounted for by a redistribution of stellar material from the disk to the bulge, without changing the extent of each of the components individually. Our analysis confirms that the change in the B/T of SFGs prior to or during quenching is, to a large extent, responsible for the size difference between the quiescent and star-forming population. However, some fraction of the shrinking size is still attributed to the individual components being smaller. In Figure 2, we show the total sizes as well as the sizes of the individual components for star-forming and QGs, as measured on the mass maps. While the total sizes of SFGs and QGs are noticeably different (by a factor ∼3 at log (M_*) ∼ 10.5), the difference in the size of the bulge and the disk components between SFGs and QGs, respectively, is smaller, typically by a factor ∼1.5.

Zoom In Zoom Out Reset image size

With the morphological parameters of the mass maps for our entire galaxy sample in hand, we now proceed to relate those with galaxy-integrated star formation properties.

The three panels of Figure 3 show, from left to right, the fraction f_quench of quenched galaxies as a function of total stellar mass, bulge mass, and disk mass, respectively. Here, we again define galaxies as quenched/quiescent when sSFR < 1/3*t_Hubble(z), and as star-forming otherwise.

Zoom In Zoom Out Reset image size

The uncertainties in f_quench are derived via a bootstrapping method and represent the 68% confidence levels. They include sample variance and the typical measurement errors on B/T, which are derived as described in Section 3.2. For the bootstrapping, we computed f_quench for 1000 samples, which are randomly drawn from the original sample, with replacement. For each bootstrap iteration, we displace the B/T values for each galaxy by the typical measurement error in B/T, given the galaxy's magnitude, size, and measured profile shape.

The top left panel of Figure 3 illustrates that the fraction of quenched galaxies increases with increasing mass, from ∼0.1 at around 10¹⁰ M_☉ to ∼0.5 at 2 × 10¹¹ M_☉. A second conclusion is that the fraction of quenched galaxies is overall higher for the lower redshift bin, by, on average, a factor of ∼2. Both of these results are well established in the literature. The rising mass function of QGs over cosmic time has most recently been quantified on a firm statistical footing by Muzzin et al. (2013) and Ilbert et al. (2013), both of which exploit the wide-area UltraVISTA survey. What CANDELS lacks in number statistics compared to UltraVISTA, it adds in depth and high resolution. Exploiting these key strengths, we now turn to the dependence of the quenched fraction on galaxy subcomponents: the mass of their bulge (middle panel) and disk (right-hand panel). Clearly, the dependence of f_quench on the bulge mass is much stronger than on the disk mass, which does not show any significant correlation with f_quench above log (M_Disk) ∼ 9.5. Toward lower disk masses, f_quench increases rapidly, but we point out that this trend is entirely driven by the (total) stellar mass limit of our sample (log (M_*) > 10; that is, the galaxies occupying the lowest M_Disk bins are necessarily heavily bulge-dominated systems, which tend to form relatively few stars). If lower mass galaxies were to be included, less massive, disk-dominated SFGs would likely outnumber these massive spheroids with small residual disks in the low M_Disk bins, producing a flat relation of f_quench with M_Disk over the full range probed. Above respective masses of 10¹⁰ M_☉, f_quench increases more rapidly with bulge mass than with total stellar mass in both redshift ranges (∼0.35 per dex of M_bulge compared to ∼0.3 per dex M_*, or ∼0.1 per dex of M_Disk).

With the bulge-to-disk decompositions in hand, we next split the galaxy sample in bins of B/T, and explore second parameter dependencies. Considering first the dependence of f_quench on the total stellar mass, it is apparent that, at a given total mass, f_quench is increasing significantly with increasing B/T ratio. The middle panels of Figure 3 illustrate that, when considering the dependence of f_quench on bulge mass, the different B/T bins align along a much tighter locus. In contrast, a large spread is seen as a function of disk mass (right-hand panels of Figure 3). In order to quantify these trends, we compute the Spearman's rank correlation coefficient (r_s) for the relations of f_quench with log (M_*), log (M_Disk), and log (M_bulge) for respective masses log (M) > 10. We find that r_s is indeed significantly higher for the relation f_quench versus log (M_bulge) (r_s ∼ 0.68), than for both f_quench versus log (M_*) (r_s ∼ 0.32) and f_quench versus log (M_Disk) (r_s ∼ −0.05), as measured for z ∼ 1. Consistent results are found for z ∼ 2.

We investigated the impact of defining quiescence based on a UVJ color–color criterion instead of a sSFR cut. When applying a UVJ-based selection of QGs, we find an overall good agreement with the trends presented in Figure 3. Quantitatively, small changes occur, with f_quench increasing by ∼7% for the entire 0.5 < z < 2.5 sample integrated over all masses. The good agreement is not surprising, since the precise threshold in sSFR used to select QGs (sSFR < 1/3*t_Hubble(z)) was chosen to yield maximum overlap with the UVJ selection criterion.

Taken together, this demonstrates that the build-up of a bulge seems to play a critical role in the quenching process of galaxies, whereas the disk does not. The amount of stars in the disk component of a galaxy has little to no predictive power regarding its star-forming or quenched state, unless B/T (and hence the bulge mass) is also known. We find a qualitatively similar behavior at z ∼ 2 as at z ∼ 1, but note that the cosmic evolution in the quiescent fraction cannot solely be attributed to continuing bulge growth over time, as galaxies in the same M_* and B/T bin at z ∼ 1 are more likely to be quenched than those at z ∼ 2. Appendix B illustrates how the equivalent diagrams composed from fits to the H-band surface brightness, rather than the stellar mass distribution, exhibit a larger spread from low to high B/T bins. This generic behavior can be understood from a physical picture, where the disk component has a relatively larger weight in light than in mass.

Our work is in agreement with, and takes the next step beyond, previous reports that the inner stellar mass density is better related to the SFH than the total stellar mass (Franx et al. 2008; Bell et al. 2012), as inferred from rest-optical imaging of smaller samples of high-redshift galaxies (see Kauffmann et al. 2003 and Fang et al. 2013 for a local universe reference, and Cheung et al. 2012 for intermediate redshifts z < 0.8)

Importantly, the same behavior explored here over the redshift range 0.5 < z < 2.5 extends in a strikingly similar fashion all the way to the present day, as demonstrated by Bluck et al. (2014) who exploit the large number statistics of SDSS.

SAMs have a rich history of trying to reproduce galaxy scaling relations and abundances, with the goal of guiding our interpretation of the observational results. Here, we focus specifically on the SAM developed by Somerville et al. (2008) and further updated by Somerville et al. (2012) and L. A. Porter et al. (2014, in preparation), which is rooted in the Bolshoi cosmological dark matter simulation (Klypin et al. 2011).²⁰ As is generic to all SAMs, the model relies on simplified analytic prescriptions for the dynamical and astrophysical processes down from entire galaxy scales, rather than on kiloparsec to parsec scales (the resolution below which state-of-the-art cosmological and zoom-in hydro-simulations resort to subgrid physics, respectively). This limitation, however, yields the enhanced flexibility of a relatively inexpensive runtime, allowing the straightforward generation of statistically significant model galaxy populations, and the tuning of parameters to observational constraints, such as mass functions and scaling relations (only empirical constraints from the nearby universe were used in tuning the parameters of the model considered here). The fact that SAMs conceptually are formulated in units of bulge and disk components furthermore makes them suitable for a direct and meaningful comparison to the diagnostics explored in this paper.

A detailed description of the prescriptions for cooling, star formation, feedback, and structural growth is provided by Somerville et al. (2008, 2012), with extensions and applications to the CANDELS data set presented by L. A. Porter et al. (2014, in preparation). The input to and output from the model is further contrasted to that of other SAMs by Lu et al. (2013).

For a detailed discussion of the physical recipes of this SAM and the resulting output in the context of a larger set of SAMs, we refer the reader to Lu et al. (2013). For the sake of the comparison presented here, we emphasize that none of the relations investigated in this paper were part of the set of observational constraints to which the model parameters were tuned. Model parameters of the SAM were tuned to (approximately) match the global stellar mass function, the stellar mass function of early- and late-type galaxies, the gas fraction as a function of stellar mass for disks, and the mass–metallicity relation for stars.

Also of particular relevance is the fact that, in the model, stars form either in the disk following a Kennicutt–Schmidt law (Kennicutt 1998), where the disk scalelength is set by a similar methodology as Mo et al. (1998), or, in the event of a merger during a starburst. Bulge formation, as well as feeding of the central supermassive black hole, can happen through two channels: mergers or disk instabilities (see L. A. Porter et al. 2014, in preparation). The starburst, black hole accretion, and morphological transformation induced by mergers depends on the mass ratio and gas fraction, as calibrated using a large suite of binary merger simulations (Somerville et al. 2008, 2012; Hopkins et al. 2009). Star formation is moderated by heating the gas by supernovae as well as through AGN feedback. In addition to the quasar mode, during which AGN can drive powerful outflows, black holes also grow more gradually over longer timespans through the so-called radio mode (i.e., suppression of cooling via radio jets). No explicit connection between the bulge mass and quenching (as may for example be expected from the Toomre Q stability criterion in a gravitational quenching scenario; see Section 6.2) was built into the model.

In Figure 4, we show the model equivalent of Figure 3, describing how the fraction of galaxies that are quenched depends on the total stellar mass, the mass of stars in the bulge, and the mass of stars in the disk component, respectively. As for the observations, we define the threshold for a galaxy to be quenched with a ruler moving with redshift: sSFR < 1/(3t_Hubble). This definition, rather than an application of the UVJ diagnostic, remains closer to the direct output of the SAM, avoiding a translation to mock spectral energy distributions, which would introduce additional assumptions and uncertainties.²¹

Zoom In Zoom Out Reset image size

At first glance, the SAM features several of the characteristic trends noted earlier for the CANDELS galaxies: f_quench rises toward later cosmic times, increases with the total stellar mass (more steeply so with bulge mass), and shows no appreciable correlation with disk mass above log (M_Disk) = 9.5. In the interval 9 ≲ log (M_Disk) ≲ 9.5, a sharp drop with increasing M_Disk is noted, as also shown in the observations. Given the M > 10¹⁰ M_☉ threshold of our sample selection, the latter objects are necessarily heavily bulge-dominated. Without imposing such a mass limit, less massive, disk-dominated SFGs would outnumber massive early-type galaxies with small residual disks in the lower M_Disk bins. Despite the qualitative success, quantitative differences in the quenched fractions of model galaxies with respect to those observed are clearly present. The discrepancy is most severe at z ∼ 2, where modeled f_quench values are, on average over the whole displayed mass range, of order a factor ∼2–3 short of the observed, hinting at an underestimated quenching rate and/or inefficiency to prevent quiescent systems from rejuvenating.²²

When splitting the SAM galaxies at each M_*, M_bulge, and M_Disk in bins of B/T (middle and bottom panels of Figure 4), we reproduce a similar behavior as found for the real universe in Section 4.3. Namely, the total stellar mass acts as a poorer predictor of the quenched state of a galaxy. This situation can be remedied if in addition to M_*, B/T (and hence M_bulge) is also known. Quantitatively, the correlation between f_quench and M_bulge is measured to be the strongest (r_s ∼ 0.46), whereas M_* and M_Disk only show weak correlation with f_quench (r_s ∼ 0.21 and r_s ∼ −0.1, respectively). The quoted values of r_s are measured for respective masses of log (M) > 10 and at z ∼ 1. The values for r_s at z ∼ 2 are similar, with the correlation between f_quench and log (M_*), as well as between f_quench and log (M_Disk), being somewhat stronger. At 0.5 < z < 1.5, less than 20% of all massive galaxies with B/T < 0.2 are classified as quiescent. Conversely, the majority of galaxies in the upper B/T bin (with B/T > 0.8) have low sSFR. These inferences are in agreement with the SAM and the observations. Also in agreement is the fact that M_bulge serves as a better predictor of f_quench than the total stellar mass, with different B/T bins being more (albeit not perfectly) aligned along a single locus in the f_quench versus M_bulge diagram.

At 1.5 < z < 2.5, the model predictions are skewed toward too low f_quench values, as noted earlier. However, in relative terms the same generic behavior, as a function of bulge, prominence is notable.

We note that most of the trends in Figure 4 are driven by physical prescriptions in the SAM affecting central galaxies rather than satellites, as centrals account for 80% (90%) of the model galaxy population above log(M) = 10 (11). Those massive galaxies classified as satellites are further subjected to additional environmental quenching processes, resulting in a higher f_quench(M_*, M_bulge) for this particular subpopulation.

Given the qualitative agreement between the model and observations, we can now pose the question how, in the context of the Somerville SAM, the relation between structure and stellar populations could be interpreted. To this end, we consider the dependence of f_quench on two physical properties of galaxies in the SAM that are not observationally accessible for our CANDELS sample: the halo mass M_Halo and the mass of the central supermassive black hole M_BH. The top panels of Figure 5 illustrate that the probability of a galaxy being quenched increases toward high M_Halo and high M_BH. In detail, however, the dependencies on the two look different. The increase in f_quench with M_Halo is gradual over nearly two orders of magnitude. In contrast, a much sharper upturn of f_quench emerges above log (M_BH) = 7.5. This behavior is especially notable at z ∼ 1, but a rise above the same threshold is also present in the z ∼ 2 population. Breaking the model galaxy population down by its structural properties, we find a wide spread in f_quench for different B/T at a given M_Halo. A different, less scattered behavior would be expected if halo mass quenching were the sole and dominant mechanism (Birnboim et al. 2007; Dekel & Birnboim 2008; Dekel et al. 2009). However, as a function of M_BH, a similar upturn in f_quench is present for all B/T bins equally above log (M_BH) = 7.5.

Zoom In Zoom Out Reset image size

The increased scatter in the relation between f_quench and M_bulge, compared to the tight correlation of f_quench and M_BH can be explained by the SAM's M_BH–M_bulge relation. For a given B/T bin, the scatter in the M_BH–M_bulge relation is significant compared to the dynamic range in the bulge masses plotted. At a given bulge mass, the scatter stems in part from an anti-correlation between the black hole mass and the level of star formation activity (sSFR).

Our analysis illustrates that in the Somerville model, which includes feedback from both supernovae and AGN, the central supermassive black hole acts as the primary agent of quenching in massive galaxies, and its accumulated mass (i.e., the integral over past accretion activity) is tightly related to the probability of finding a galaxy in a quenched state. Since the physical processes giving rise to bulge and black hole growth are the same in the SAM (mergers and disk instabilities), the stronger relation of f_quench with M_bulge than with M_*, present in both observations and model predictions, is not surprising.

We stress that in the Somerville model, no direct causal link between the presence of a bulge and quenching is implemented. The bulge is simply the accessible observable that correlates most tightly with the actual agent of quenching in this particular model: the supermassive black hole. Observationally, there is increasing evidence for AGN-driven outflows of massive z ∼ 1–3 galaxies (Förster Schreiber et al. 2014; Nesvadba et al. 2011; Harrison et al. 2012; Cano-Díaz et al. 2012). Given the shortcomings of the SAM in a quantitative sense, notably its under-prediction of the quiescent population at z ∼ 2, additional or other quenching processes may be at play in the real universe. One such process, which is causally linked to the presence of a bulge, could be morphological quenching (Martig et al. 2009; Genzel et al. 2014). In such a scenario, the high central stellar density provided by a bulge stabilizes the gas disk and prevents it from forming stars. While plausibly only a temporary measure (as no gas is expelled, or stopped from accreting through this mechanism), it could potentially contribute to suppressing star formation in z ∼ 2 galaxies more efficiently and/or prevent them from returning to the star-forming branch in the SAM. In the local universe, star formation efficiencies of bulge-dominated systems are reduced by factors of ∼2–3 compared to disk-dominated galaxies (Saintonge et al. 2012; Martig et al. 2013).

In order to study the morphological differences between SFGs and QGs, and to draw conclusions on the possible structural changes of SFGs as they move along the MS, we first examined the mass dependence of the profile shape of SFGs and QGs as traced by the Sérsic index and B/T ratio. We have shown that QGs are structurally distinct from the star-forming population as seen by overall higher Sérsic indices and B/T ratios at a given stellar mass. SFGs show rising trends of their median Sérsic index and B/T ratio with increasing stellar mass, with the latter rising up to ∼40%–50% above 10¹¹ M_☉. These findings give insights about the link between the structural evolution of SFGs and the quenching process as they move along the MS.

Analyzing the Schechter functional forms of the SFG and QG mass function as a function of redshift, Peng et al. (2010) conclude that the quenching rate of galaxies climbing the MS rises proportionally to the SFR (and given the near-linear MS slope therefore also proportionally to the stellar mass, hence their terminology "mass quenching"). This corresponds to a survival probability on the MS that drops exponentially with mass, implying that, while nearly all low-mass SFGs are destined to continue growing along the MS, toward the high-mass end the MS becomes progressively more dominated by nearly dead SFGs. In fact, the sub-unity slope of the MS, and possible flattening at the high-mass end (Whitaker et al. 2012), may well be interpreted in this context—the typical SFG above 10¹¹ M_☉ is already undergoing some level of quenching, thereby deviating from the projected path along a SFR–mass relation of slope unity, which could be expected from cosmological accretion rates in the absence of quenching. Tying in our observational results on galaxy structure, the deviation toward high median n and B/T at the massive end reflects the typical structure of nearly dead star-formers that account for the bulk of SFGs above 10¹¹ M_☉. The fact that they are structurally distinct implies that the morphological transition happens first, to be followed later by the departure from the MS. Bulge growth precedes quiescence. Such a morphological change prior to quenching is in line with qualitative predictions based on a toy model by Dekel & Burkert (2013). In the latter study, about half of the star-forming disk galaxies at z ∼ 2 are predicted to evolve into compact star-forming "blue nuggets" due to violent disk instabilities before they are quenched into compact QGs ("red nuggets"). An observed population of "blue nuggets" has been proposed by Barro et al. (2013a, Barro et al. 2013b) to represent an evolutionary link, originating from extended disk galaxies and evolving into compact quiescent systems.

This does not refute that galaxies also undergo further structural evolution after they are quenched. At least part of the size growth (Cassata et al. 2013; van der Wel et al. 2014) and evolution toward rounder axial ratios (Chang et al. 2013) has been attributed to (minor and/or major) dry mergers, and it is conceivable that similar processes contribute to the observed trend of increasing B/T toward the massive end for the quiescent population.

We have demonstrated that the bulge mass of a system is well correlated with its quenched state and has a stronger predictive power of quiescence than the total stellar mass. The observed trends of f_quench with total stellar mass, bulge mass, and disk mass as viewed among galaxies in different B/T bins are in good qualitative agreement with predictions from the Somerville et al. SAM. In the context of this model, the growth of the central supermassive black hole, which is the primary quenching agent for massive galaxies in this SAM, is tightly coupled with the growth of bulges through both merging and disk instabilities.

If a black hole—bulge scaling relation is in place during the peak of cosmic star formation as it is in the present-day universe (Häring & Rix 2004), our observational results together with the model comparison could therefore hint at the bulge not being the causal link to quenching, but rather the most accessible observational proxy for the AGN acting as the quenching agent.²³ In detail, however, there are quantitative differences between the SAM and our observations, most severely in the highest redshift bin (1.5 < z < 2.5), where the observed quenched fraction exceeds the value predicted by the SAM by a factor of ∼3.5. The latter difference could hint at a need for more frequent, efficient, or lasting quenching, a possible mechanism that we speculate about below. We also note that the same behavior is not necessarily a generic feature to all SAMs (see Appendix C).

It is tempting to draw connections between the emerging bulges in massive MS galaxies out to z ∼ 2.5 revealed by our analysis, and recent observational results based on deep adaptive-optics-assisted integral field data sets and grism spectroscopy over the same redshift range. Förster Schreiber et al. (2014) found a high prevalence of powerful nuclear outflows in log (M) > 11 galaxies driven by AGN, which appear to be absent in galaxies at lower masses. Along with star formation driven winds in the outer parts of the galaxies, such outflows could efficiently remove gas from the galaxies and, in this way, contribute to the quenching process.

Meanwhile, the 3D-HST and CANDELS legacy programs have yielded evidence for nuclear depressions in the Hα equivalent width in z ∼ 1 SFGs (Nelson et al. 2012, 2013; Wuyts et al. 2013). At the highest stellar surface mass densities, star formation no longer appears to proceed in lockstep with the assembled stellar mass. Likewise, Genzel et al. (2014) report on ring-shaped Hα distributions in z ∼ 2 SFGs, surrounding a more quiescent center where the dynamically inferred Toomre Q parameter significantly exceeds unity, owing to the emergence of a stellar bulge. Thus the Toomre stability criterion is satisfied in the central galaxy regions, which consequently could prevent the gas reservoir, if present there, from fragmenting and forming stars. While this result suggests that some causal connection between bulge growth and quenching may be at play, it should be noted (as is done also by Genzel et al. 2014) that gravitational quenching by itself does not expel the gas present, or stop the accumulation of a larger gas reservoir by continuing cosmological accretion. Additional maintenance modes might be required for a long-term shut down of further gas supply.