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

CrossMark

# Constraints on Quenching of $Z \lesssim 2$ Massive Galaxies from the Evolution of the Average Sizes of Star-forming and Quenched Populations in COSMOS

A. L. Faisst¹,², Ç. M. Carollo¹, P. L. Capak²,³, S. Tacchella¹, A. Renzini⁴, O. Ilbert⁵, H. J. McCracken⁶, and N. Z. Scoville³

¹Institute for Astronomy, Swiss Federal Institute of Technology, 8093 Zürich, Switzerland; afaisst@ipac.caltech.edu

²Infrared Processing and Analysis Center, California Institute of Technology, Pasadena, CA 91125, USA

³Cahill Center for Astronomy and Astrophysics, California Institute of Technology, Pasadena, CA 91125, USA

⁴NAF—Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, I-35122 Padova, Italy

⁵Aix Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille) UMR 7326, F-13388, Marseille, France

⁶Institut d'Astrophysique de Paris, CNRS & UPMC, UMR 7095, 98 bis Boulevard Arago, F-75014, Paris, France

Received 2016 July 15; revised 2017 March 22; accepted 2017 March 23; published 2017 April 17

## Abstract

We use $>9400 \log(m/M_\odot) > 10$ quiescent and star-forming galaxies at $z \lesssim 2$ in COSMOS/UltraVISTA to study the average size evolution of these systems, with focus on the rare ultra-massive population at $\log(m/M_\odot) > 11.4$. The large 2 square degree survey area delivers a sample of $\sim 400$ such ultra-massive systems. Accurate sizes are derived using a calibration based on high-resolution images from the *Hubble Space Telescope*. We find that at these very high masses, the size evolution of star-forming and quiescent galaxies is almost indistinguishable in terms of normalization and power-law slope. We use this result to investigate possible pathways of quenching massive $m > M^*$ galaxies at $z < 2$. We consistently model the size evolution of quiescent galaxies from the star-forming population by assuming different simple models for the suppression of star formation. These models include an instantaneous and delayed quenching without altering the structure of galaxies and a central starburst followed by compaction. We find that instantaneous quenching reproduces the observed mass-size relation of massive galaxies at $z > 1$ well. Our starburst+compaction model followed by individual growth of the galaxies by minor mergers is preferred over other models without structural change for $\log(m/M_\odot) > 11.0$ galaxies at $z > 0.5$. None of our models is able to meet the observations at $m > M^*$ and $z < 1$ without significant contribution of post-quenching growth of individual galaxies via mergers. We conclude that quenching is a fast process in galaxies with $m \gtrsim 10^{11} M_\odot$, and that major mergers likely play a major role in the final steps of their evolution.

**Key words:** galaxies: evolution – galaxies: fundamental parameters – galaxies: structure

## 1. Introduction

Quiescent (or quenched) galaxies—here defined to be galaxies that have heavily suppressed specific star formation rates (specific SFRs) relative to the star-forming “main-sequence” galaxies (e.g., Daddi et al. 2007; Noeske et al. 2007)—host about half of the mass in stars in the local Universe (Baldry et al. 2004) and have been observed in substantial numbers as early as $z \sim 2$ (e.g., Ilbert et al. 2013; Muzzin et al. 2013; Davidzon et al. 2017). Understanding the dominant processes responsible for the shutdown of their star formation (often referred to as “quenching”) as well as the connection between these processes and galaxy structure are key for understanding the evolution of the whole galaxy population over cosmic time.

Suppressed specific SFRs are not the only intriguing property of quiescent galaxies. At least in terms of “light”, they have on average substantially larger spheroidal components and smaller half-light radii ($R_c$) than their star-forming counterparts at a given stellar mass and redshift (e.g., Shen et al. 2003; Szomoru et al. 2012; Cibinel et al. 2013b). At least in part this size difference is likely contributed by “nurture”, in particular post-quenching “fading” of stellar populations at large radii (Tacchella et al. 2015; Carollo et al. 2016), but it is also possible that part of the difference may be imprinted by “nature”, i.e., different formation processes for spheroids and disks. It is also intriguing that the population-averaged sizes of quiescent galaxies of a given mass have increased by a factor of $\sim 3$ since $z = 2$. This average size growth is similar to that of star-forming

disk galaxies, which are expected and observed to increase their individual (disk) sizes more or less proportional to $(1+z)^{-1}$ through continuous accretion of gas from their halos (e.g., Mo et al. 1998; Oesch et al. 2010a; Mosleh et al. 2012; Newman et al. 2012; Huang et al. 2013; Shibuya et al. 2015; Mosleh et al. 2017). Individual quiescent galaxies by definition form no new stars, however, and thus their only channel for individual mass and size growth is provided by gas-poor mergers. Averaging over a large mass range, several studies suggest indeed that mergers are important contributors to the size growth of quiescent galaxies (e.g., Toft et al. 2007; Buitrago et al. 2008; Franx et al. 2008; Stockton et al. 2008; Kriek et al. 2009; Williams et al. 2010; Newman et al. 2012; Oser et al. 2012; Whitaker et al. 2012; Belli et al. 2014, 2015).

Analyses in thinner bins of stellar mass, however, suggest a threshold mass—roughly around $m < M^* \sim 10^{11} M_\odot$, the characteristic mass of the Schechter (Schechter 1976) fit to galaxy mass functions<fnref n="7" />—below and above which different mechanisms may be responsible for the average size growth of quiescent galaxies. In particular, at $m \lesssim M^*$, a number of studies indicate that the growth in average size of the quiescent population is dominated by the addition of larger galaxies at later times, as a result of the continuous addition of newly quenched galaxies to the large-size-end of the size function (Carollo et al. 2013a; Cassata et al. 2013; Saracco et al. 2014).

<footnote n="7">⁷ The value of $M^*$ is remarkably constant for star-forming and quiescent galaxies and at all epochs since $z \sim 4$; (see e.g., Ilbert et al. 2013; Muzzin et al. 2014; Davidzon et al. 2017).</footnote>
This picture is substantiated by the stellar ages of compact (older) and large (younger) quiescent galaxies at a given stellar mass and epoch (Saracco et al. 2011; Onodera et al. 2012; Carollo et al. 2013a; Belli et al. 2015; Fagioli et al. 2016; Williams et al. 2017). It is only above $M^*$ that dissipationless mergers are expected to be important (e.g., Peng et al. 2010b), and all studies of galaxy sizes indeed agree on them playing the dominant role in leading to the growth of individual quiescent galaxies in mass and size (Carollo et al. 2013a; Poggianti et al. 2013; Belli et al. 2014, 2015).

The above results may indicate that different quenching mechanisms could be at work below and above $M^*$. Theoretically, there are many candidate mechanisms for quenching (e.g., Birnboim & Dekel 2003; Croton et al. 2006; Birnboim et al. 2007; Bournaud et al. 2007; Kawata & Mulchaey 2008; Martig et al. 2009; Peng et al. 2010b; Feldmann et al. 2011; De Lucia et al. 2012; Hearin et al. 2013; Cen 2014; Dekel & Burkert 2014; Mandelker et al. 2014; Schaye et al. 2015; Tacchella et al. 2016a, 2016b), and observationally identifying the correct mechanisms is a non-negligible challenge—not least since, as discussed in Carollo et al. (2013b) and demonstrated in Lilly & Carollo (2016), correlations of observed quantities do not necessarily indicate a causal relation between them.

Different quenching mechanisms are expected to act on different timescales and result in different morphological transformations of galaxies. Constraining them therefore is an important step toward understanding the dominant processes that lead to galaxy quiescence in these populations. For example, the cutoff of gas inflow onto star-forming galaxies is expected to lead to the exhaustion of star formation over long timescales, which is set by the time needed for star formation to consume the gas reservoir of a galaxy. It is likely that the star formation ceases smoothly over the galaxy disk, thereby not significantly changing its observed morphology. In contrast, a gas-rich major merger might lead to a starburst and thus to a fast consumption of gas on dynamical timescales of order 100–200 Myr. Furthermore, a substantial change in the morphology of the galaxies is expected, with an apparent compaction in light induced by the centrally confined starburst (e.g., Barro et al. 2013; Zolotov et al. 2015). A number of studies have focused their attention on quenching timescales and their dependence on galaxy properties. For relatively massive galaxies, and in particular for satellites in groups and clusters at low redshifts, there is growing evidence that the transition from active star formation to quiescence takes of order 2–4 Gyr (von der Linden et al. 2010; De Lucia et al. 2012; Cibinel et al. 2013a; Mok et al. 2013; Trinh et al. 2013; Wetzel et al. 2013; Hirschmann et al. 2014; Muzzin et al. 2014; Schawinski et al. 2014; Taranu et al. 2014; Peng et al. 2015). At redshifts of $z \sim 2$ and for massive galaxies, Tacchella et al. (2015, 2016a) show that suppression of star formation starts at the center of galaxies and slowly progresses outward on timescales of 1–3 Gyr. It is unclear at this point, however, whether the observed centrally suppressed specific SFRs are the outcome of a “quenching mechanism” (e.g., through central gas and stellar “compaction”; Dekel & Burkert 2014; Tacchella et al. 2016b) or if they are the natural outcome of inside-out galaxy formation (e.g., Lilly & Carollo 2016).

In this paper we make a new attempt to constrain the processes that quench massive $m > M^*$ star-forming galaxies

at $z < 2$ via studying the timescales and morphological changes, using as diagnostic tool the size evolution of both the star-forming and quiescent galaxy populations at $z < 2$. We show that this further enables us to set constraints on the amount and properties of mergers in this massive population.

It is now well established that studies of $z \sim 2$ galaxies crucially need imaging in the near-infrared in order to measure their rest-frame optical properties (especially sizes); the near-infrared images furthermore need to be quite deep in order to detect the faint and fading stellar populations of quiescent systems. Much progress has been made using data from the Hubble Space Telescope (HST) CANDELS survey (Grogin et al. 2011; Koekemoer et al. 2011). Very massive galaxies are rare, however, and increasingly so at increasingly higher redshifts (with number densities lower than $10^{-4}$ per Mpc$^3$ at $z > 1$; Ilbert et al. 2013). Assembling a sufficiently large number of such galaxies to enable a statistical study requires imaging over a large area of the sky. With its 2 square-degree area coverage, the Cosmological Evolution Survey (COSMOS,<fnref n="8" /> Scoville et al. 2007) enables us to assemble a sample of more than 400 ultra-massive galaxies (UMGs, $\log(m/M_\odot) > 11.4$) in the redshift range $0.2 < z < 2.5$. Another advantage of COSMOS is its $>30$ passband coverage from UV to IR wavelengths, which enables the derivation of very accurate stellar masses and photometric redshifts. Last but not least, the deep near-IR data of the UltraVISTA survey on COSMOS (McCracken et al. 2012; Laigle et al. 2016) allow an accurate separation of star-forming and quiescent galaxies across this entire redshift range (e.g., Ilbert et al. 2013). A drawback of the UltraVISTA data is their seeing-limited point-spread function (PSF), whose full width at half maximum (FWHM) is typically about 0.8 and therefore hampers the measurement of reliable galaxy sizes. To overcome this limitation, we correct the UltraVISTA size measurements using as a calibration reference the ~3% of the COSMOS area that is covered by the CANDELS/COSMOS legacy survey.

The paper is organized as follows. In Section 2 we describe the data sets that we have used in this work. In Section 3 we describe the selection criteria for separating star-forming and quiescent UMGs, and in Section 4 the procedure that we have followed to measure the galaxy sizes. In the same Section we also present the calibration of the UltraVISTA sizes that we have performed using the HST CANDELS size measurements for the ~9000 galaxies for which both data sets are available. The final (calibrated) size measurements are presented and discussed in Section 5. In Section 6 we present our model that we use to predict the average size evolution of quiescent galaxies through the redshift range of our analysis. The model predictions are compared with the observed size evolutions in Section 7, where we furthermore describe the additional modifications to the predicted trends that are introduced by galaxy mergers. We summarize our main results in Section 8.

Note that all magnitudes are given in the AB system (Gunn et al. 1986); stellar masses ($m$) are scaled to a Chabrier (2003) initial mass function (IMF); we assume a flat cosmology with $\Omega_\Lambda = 0.7$, $\Omega_m = 0.3$, and $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$.

<footnote n="8">⁸ http://cosmos.astro.caltech.edu/</footnote>
## 2. Data

### 2.1. UltraVISTA Near-IR Imaging Data

As mentioned in the previous section, near-IR data on a large area are crucial for the study of massive galaxies at high redshifts. Therefore the backbone of this work is the UltraVISTA survey carried out on the 4.1-meter Visible and Infrared Survey Telescope for Astronomy (VISTA) located at the Paranal observatory in Chile. This survey covers $1.5 \text{ deg}^2$ of the COSMOS field in the near-IR bands Y, J, H, and Ks. Specifically, we use the UltraVISTA data release (DR) 2 imaging data. Compared to DR1, this release has an improvement in H band by up to 1 mag in the ultra-deep stripes (covering roughly 50% of the field) and ~0.2 mag on the deep stripes. The typical exposure times per pixel are between 53 and 82 hr, leading to $5\sigma$ sensitivities of 25.4AB, 25.1AB, 24.7AB, and 24.8AB in Y, J, H, and Ks band within a $2''$ aperture. The reduction of the imaging data is similar to DR1 (see McCracken et al. 2012) and is briefly outlined in the following: the data were taken in three complete observing seasons between 2009 December and 2012 May. The individual science frames are visually inspected to remove bad frames (e.g., due to loss of auto-guiding). Each frame is sky subtracted before stacking, which leads to a very flat combined image with a very small variation in background flux. The combined frames have an *average H-band* seeing of $0.75 \pm 0.10$. The final photometric calibration is done by using nonsaturated stars from the Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006) sample, leading to an absolute photometric error smaller than 0.2 mag.

### 2.2. Photometric Redshift and Stellar Mass Catalog

Our galaxy selection (see below) is based on the public COSMOS/UltraVISTA catalog in which galaxies are selected from a combined $YJHK_s$ image (Ilbert et al. 2013). This has advantages compared to purely optical selected catalogs as it is more sensitive to galaxies with red colors, e.g., dusty star-forming galaxies or quiescent galaxies with old stellar populations. The catalog comprises photometric redshifts, stellar masses, and other physical quantities derived from SED fitting on $>30$ passbands from UV to IR (PSF homogenized) for more than 250,000 galaxies on COSMOS (see e.g., Capak et al. 2007; Ilbert et al. 2013). The photometric redshifts in that catalog are derived using Le Phare (Arnouts et al. 2002; Ilbert et al. 2006), employing different templates including a range of galaxy types from elliptical to young and star-forming. These redshifts have been verified to have a precision of $\sigma_{\Delta z/(1+z)} = 0.01$ up to $z=3$ by comparison to a sample of more than ~10,000 spectroscopically confirmed star-forming and quiescent galaxies. Physical quantities (mass, SFR, etc.) are fitted by Le Phare at fixed photometric redshift using a library of synthetic composite stellar population models based on Bruzual & Charlot (2003). These models include different dust extinctions (following a Calzetti et al. (2000) dust extinction law), metallicities, and star formation histories (following exponentially declining $\tau$ models). Emission line templates are also included. The emission line flux is derived from the observed UV light using empirical relations. All these parameters have been verified by a number of other fitting routines, including ZEBRA (Feldmann et al. 2008) and its upgraded version ZEBRA+ (Oesch et al. 2010b; Carollo et al. 2013a). The typical uncertainties in masses are on the

order of 0.3 dex. All quantities are computed for a Chabrier (2003) IMF. The stellar masses are defined as the integral of the star formation histories of the galaxies, thus representing the total galaxy mass of a galaxy rather than its mass in active stars. In the following, the stellar masses quoted by other studies are converted into total masses if necessary. These corrections, calculated using Bruzual & Charlot (2003) models with solar metallicity and exponentially declining as well as constant star formation histories, can be up to 0.2 dex for quiescent galaxies with ages of one billion years and older, while they are less substantial for star-forming galaxies.

### 2.3. CANDELS/COSMOS Near-IR Imaging Data

To calibrate the sizes measured on the ground based on UltraVISTA imaging data, we make use of the overlap between UltraVISTA and the HST-based CANDELS/COSMOS survey (Grogin et al. 2011; Koekemoer et al. 2011). The latter covers $0.06 \text{ deg}^2$ on sky (roughly 1/25 of the total UltraVISTA field) in the WFC3/IR F160W passband, similar to the UltraVISTA H band, but at a much higher resolution (the PSF is more than eight times smaller). We use the latest publicly available data release of the COSMOS/F160W mosaic (by 2013 February) with a total exposure time of 3200 s and a sensitivity of 26.9 AB ($5\sigma$ for a point source).

## 3. The Sample

In the following, we describe the selection of massive galaxies at $\log(m/M_\odot) > 11.4$ from which we built our main galaxy sample, as well as less massive galaxies ($10.0 < \log(m/M_\odot) < 11.4$) that we use for the calibration of the ground-based size measurements. Furthermore, we split this sample into quiescent and star-forming galaxies.

### 3.1. High- and Low-mass Galaxies

The selection of the high- and low-mass galaxy sample is based on the near-IR COSMOS/UltraVISTA photometric catalog (as described above), which allows for the selection of dusty star-forming and quiescent galaxies.

We select a total sample of 403 massive galaxies satisfying $\log(m/M_\odot) > 11.4$ and $0.2 < z_{phot} < 2.5$ (green hatched region in Figure 1). We have visually verified these galaxies to be real (i.e., not artifacts or stars). The exact value of this mass limit has been chosen to correspond to the 90% completeness limit at a H-band magnitude of 21.5 AB at $z < 2.5$, which allows us to provide reliable size measurements for these galaxies (see Section 4). For the estimation of the mass completeness we have used the identical method as described in Pozzetti et al. (2010). With this mass cut, we select the most massive observable galaxies with a number density lower than $10^{-4} \text{ Mpc}^{-3}$ and $10^{-5} \text{ Mpc}^{-3}$ at $z \sim 0.5$ and $z \sim 2$. These galaxies may be the progenitors of today's most massive galaxies, assuming these most massive galaxies keep their ranking through cosmic time. This is verified by more complicated methods of progenitor selections, including the selection of galaxies at a constant galaxy number density (Marchesini et al. 2014), or using semi-empirical models that take galaxy mergers into account (Behroozi et al. 2013).

The (mass-complete) low-mass galaxy control/calibration sample is selected in a similar way to have $10.0 < \log(m/M_\odot) < 11.4$ and $H < 21.5$ AB. The mass completeness limit at $H = 21.5$AB as a function of redshift is shown in
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Sample selection. Our sample of UMGs ($\log(m/M_{\odot}) > 11.4$) at $0.2 < z < 2.5$ is shown with large symbols. The sample is split into quiescent (red circles) and star-forming (blue squares) according to their location on the rest-frame (NUV $-r$) vs. ($r-J$) diagram (see also Figure 2). Other galaxies with $\log(m/M_{\odot}) < 11.4$ and $H < 21.5AB$ are shown in gray. The dark red line shows the 90% mass completeness for star-forming (solid) and quiescent (dashed) galaxies as described in the text. The low-mass sample therefore consists of the galaxies shown in gray scale that are located in the dark red hatched region.</figcaption></figure>

Figure 1 by the red line (solid for star-forming and dashed for quiescent galaxies). The low-mass control sample (9000 galaxies in total) is consequently selected to be above the combined completeness limit of the star-forming and quiescent galaxies and satisfies three stellar mass bins of $10.0 < \log(m/M_{\odot}) < 10.5$, $10.5 < \log(m/M_{\odot}) < 11.0$, and $11.0 < \log(m/M_{\odot}) < 11.4$ with the corresponding redshift ranges $0.2 < z < 0.45$, $0.2 < z < 0.75$, and $0.2 < z < 1.25$.

## 3.2. Selection of Quiescent and Star-forming Galaxies

We split our sample into quiescent and star-forming galaxies by making use of the rest-frame (NUV $-r$) versus ($r-J$) color diagnostics (see Williams et al. 2009; Ilbert et al. 2010; Carollo et al. 2013a; Ilbert et al. 2013). In Figure 2 we show the rest-frame (NUV $-r$) versus ($r-J$) diagram for six different redshift bins with our main sample of massive galaxies. The black line in each panel divides the quiescent (upper left) from the star-forming (lower right) galaxy population. Our $\log(m/M_{\odot}) > 11.4$ galaxies are shown with large symbols color-coded by their specific star formation rate (sSFR $\equiv$ SFR/m, the inverse of the mass-doubling timescale) derived from SED fitting. All the other galaxies at lower stellar masses in the same redshift bin and $H < 21.5AB$ are shown in gray scale. We find that the color-color diagram efficiently isolates quiescent galaxies with $\log(sSFR/Gyr) \sim -1$ to $-2$ (depending on redshift, as expected). We note that this color selection is very similar to the widely used ($U-V$) versus ($V-J$) selection, but it is a slightly better indicator of the current versus past star formation activity (e.g., Arnouts et al. 2007; Martin et al. 2007). We have verified that other selections of quiescent and star-forming galaxies (e.g., by sSFR or ($U-V$) versus ($V-J$)) do not change the results of this paper.

# 4. Size Measurements and Calibration

As we have already discussed in the introduction to this paper, we are investigating the quenching process in massive galaxies via the average size evolution of star-forming and quiescent galaxies. Reliable size measurements are therefore crucial. We denote with “size” the observed semimajor axis half-light radius, $R_e$. While we benefit from the large area of the COSMOS/UltraVISTA survey to select very massive galaxies, its poor resolution and PSF hampers the accurate measurement of galaxy structure parameters.

In this section, we lead in detail through (i) the determination of a spatially varying PSF, (ii) the basic measurement of galaxy sizes, and (iii) our two-step size-calibration procedure using simulated galaxies and the HST-based CANDELS imaging. Finally, we outline how we correct for the band-shifting across redshift in our sample.

## 4.1. Determination of the Spatially Varying PSF

Galaxy sizes are measured by the use of GALFIT, which takes the effect of PSF into account (Peng et al. 2010a). Understanding the PSF size (FWHM), shape, and spatial variation is therefore crucial. We represent the two-dimensional PSF at a given position $(x, y)$ by a Moffat profile (Moffat 1969):

$$F(x, y) = \frac{\beta - 1}{\pi \alpha^2} \left[ 1 + \left( \frac{(x - \mu_x)^2 + (y - \mu_y)^2}{\alpha^2} \right)^{-\beta} \right], \quad (1)$$

where $\mu_x, \mu_y, \alpha$, and $\beta$ are free fitting parameters. The FWHM of a PSF in this parametrization is given by

$$\text{FWHM}(\alpha, \beta) = 2\alpha\sqrt{2^{1/\beta} - 1}. \quad (2)$$

This has been shown to be a good approximation for ground-based PSFs and has the advantage over a pure Gaussian as it represents the wings of the PSF better (e.g., Trujillo et al. 2001). In order to create a spatially comprehensive PSF map, we select unsaturated stars between 16 AB and 21 AB from the HST-based COSMOS/ACS $I_{F814W}$-band catalog (Leauthaud et al. 2007). We select them according to their SExtractor stellarity parameter (larger than 0.9) and using diagnostic diagrams such as color versus color and magnitude versus size. Furthermore, we inspect the stars visually and verify that there are no close companion stars (or galaxies) visible on the ACS images.

For each of these more than 3000 stars, we extract a $10'' \times 10''$ image stamp from the UltraVISTA H-band mosaic on which we fit the PSF. We note small shifts of the center of the stars between ACS and UltraVISTA data of a few tenths of arcseconds (likely caused by small differences in the coordinate systems, the large differences in the PSF size, and differences in the resolution of the images), which we correct for. We then fit the selected stars according to the above parametrization $F(x, y|\mu_x, \mu_y, \alpha, \beta)$. The accuracy and robustness of the fitting method was verified by generating stars with random FWHM between $0''2 < \text{FWHM} < 1''2$, add noise taken from real background images, and fit them in the same way as the real data. This test shows that we are able to recover the FWHM with an accuracy of better than $0''05$. As a last cut, we require less than 5% difference between the model and data in the enclosed flux up to 1.5 times the PSF FWHM. We end up with $\sim 800$ PSF models across COSMOS/UltraVISTA. The PSFs show variations in their FWHM between $0''65$ and $0''80$.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Selection of quiescent (upper panels, left of the solid line) and star-forming (lower panels, right of the solid line) UMGs on the rest-frame (NUV − *r*) vs. (*r* − *J*) diagram for six different redshift bins. The UMGs are color-coded by their sSFR. This shows that the color-color cut efficiently separates quiescent galaxies with log(sSFR/Gyr) ~ −1 to −2. The gray background shows less massive galaxies with *H* < 21.5AB in the same redshift bins.</figcaption></figure>

assign to each galaxy an average PSF model created from the
stars within 6', which we use for GALFIT.

4.2. *Guess-parameters for Surface Brightness Fitting*

In this section, we describe the determination of the initial
values that are fed to GALFIT. In order to have consistency
between the initial values and the actual images on which we
run GALFIT, we do not use the values given in the public
COSMOS/UltraVISTA catalog, but we rerun Source Extractor
(SExtractor, version 2.5.0, Bertin & Arnouts 1996) on the
DR2 UltraVISTA H-band images. We run SExtractor with
two different values of the DEBLEND_MINCONT for a better
deblending of galaxies next to brighter galaxies or stars. The
SExtractor input parameters are tuned manually in order to
optimize the source extraction. We mask (using Weight-
watcher; Marmo & Bertin 2008) each star identified on the
HST-based COSMOS/ACS *I*<sub>F814W</sub>-band images by a circle
with a maximal radius *r*<sub>σ</sub> at which its flux decays to the
background flux level. This maximal radius (which depends on
the magnitude of the star) is determined by fitting *r*<sub>σ</sub> as a
function of magnitude for several different stars in a broad
magnitude range. Furthermore, we match our catalog to the
public UltraVISTA catalog and compare the measured
magnitudes, which we find to be in excellent agreement.
Finally, we extract each of our galaxies from our

SExtractor catalog to use the measured galaxy position
(X_IMAGE and Y_IMAGE), magnitude (MAG_AUTO), half-
light radius (FLUX-radius), axis ratio (ratio of A_IMAGE
and B_IMAGE), and position angle (THETA_IMAGE) as initial
parameters for GALFIT.

4.3. Uncalibrated Size Measurements

We use GALFIT to fit single Sérsic profiles (parametrized by
the half-light radius R<sub>e</sub>, Sérsic index *n*, total magnitude *M*<sub>tot</sub>,
axis ratio *b*/*a*, and position angle θ) to the observed surface
brightness of our galaxies. As described in the previous section,
we use the *SExtractor* values measured on the DR2
COSMOS/UltraVISTA images as initial parameters. For the
Sérsic index, which is not known a priori, we assume *n* = 2
(and let it vary between 0 < *n* < 8 during the fitting process).
The size of the image cutout on which GALFIT is run is
variable between 71 × 71 and 301 × 301 pixels. The size is
set to optimize the estimate of local sky background and to
minimize the running time of GALFIT, and it is defined such
that the cutout contains three times more sky pixels than pixels
attributed to galaxy detections. Companion galaxies on the
image cutout are fit simultaneously with the main galaxy if they
are brighter than 25AB in *H* band. All other detections of
fainter objects are masked out and are not taken into account in
the χ<sup>2</sup> minimization. To access the stability of the fits, we run
GALFIT in two different configurations: In the first config-
uration (referred to as “VARPOS”) we let GALFIT fit the center
of the galaxy within ±10 pixels of the SExtractor input. In
the second configuration (referred to as “FIXPOS”) we fix the
galaxy position to its initial SExtractor value.

We select good fits (either from the *FIXPOS* or *VARPOS*
run) by comparing the results from the two configurations. We
require that (i) *R*<sub>e</sub> > 0.1 px, (ii) the fitted position differ by less
than √2/2 times the PSF FWHM from the *SExtractor*
input, (iii) the *R*<sub>e</sub> of the two configurations agree better than
50%, and (iv) the total magnitude does not differ by more than
0.5 from the *SExtractor* total magnitude. Roughly 70% of
our total sample galaxies satisfy these criteria and are used in
the following to access the size evolution as a function of
cosmic time. Owing to their brightness and relatively large
size, the above criteria result in a negligible cut for our massive
log(*m*/*M*<sub>☉</sub>) > 11.4 but in principle could affect the following
results and conclusions. We have investigated this in depth and
find that mostly unresolved galaxies are affected by this,
without any clear relation with redshift. However, adding this
small amount of galaxies to our sample at log(*m*/*M*<sub>☉</sub>) > 11.4
(keeping their small sizes as lower limits) impacts the median
size versus redshift relations by less than 5% compared to the
general systematic uncertainties of the ground-based sizes of up
to 50%. Furthermore, star-forming and quiescent galaxies are
equally affected, and therefore we do not expect significant
impacts on our results.

## 4.4. Correcting for Measurement Biases using Simulated Galaxies

The measurement of galaxy structure is prone to many biases,
as discussed by several authors (Cameron & Driver 2007;
Carollo et al. 2013a; Cibinel et al. 2013b). Small and compact
 galaxies are affected by the PSF (leading to an overestimation of
 $R_e$); large and extended galaxies suffer surface brightness
dimming in the outskirts (leading in underestimation of $R_e$).
Although GALFIT does take into account the effects of PSF and
therefore partially cures these problems, it has its limits. It is
therefore important to investigate possible biases and correct for
them by using simulated galaxies. In the following, we outline
this first step in our two-step calibration process in more detail.

4.4.1. Simulating Galaxies

We use GALFIT to create ~1.5 million model galaxies on a
grid in (R<sub>e</sub>, M<sub>tot</sub>, n, b/a)in parameter space: 0.2 < n < 10,
15 mag < M<sub>tot</sub> < 26 mag, 0.2 < b/a < 1, and 0.5 < R<sub>e</sub> < 15
pixels (corresponding to 0.075 < R<sub>e</sub> < 2.250). The model
 galaxies are subsequently convolved with a PSF, equipped with
 Poisson noise, and added onto realistic sky backgrounds. For
 the latter, we account for the fact that the sky background noise
 (σ<sub>sky</sub>) varies across the COSMOS/UltraVISTA field by a factor
 2 or more (mainly between the deep and ultra-deep stripes). We
 compute σ<sub>sky</sub> automatically in rectangles of ~0.1 × 0.1
degrees across the field. For this end, we use the SExtractor
catalog (see Section 4.2) to mask out all detections and fit σ<sub>sky</sub>
to the remaining non-masked pixels by assuming a Gaussian
noise distribution. In order to ensure that we remove all the
light of galaxies and stars, we increase their semimajor and
semiminor axis as given by SExtractor by a factor of 10.
We verify this procedure by manually measuring σ<sub>sky</sub> at

random positions. To take the variations in PSF and σ<sub>sky</sub> into
account, we simulate galaxies in four different representations,
which will be interpolated in the end. We use two bracketing
PSFs (FWHM = 0.65 and 0.85) as well as two bracketing
σ<sub>sky</sub> (5.5 × 10<sup>-6</sup> and 2.0 × 10<sup>-5</sup> counts/s). On each of these
model galaxies we run SExtractor and GALFIT in the
same manner as for the real galaxies (as described in
Section 4.3) to obtain (R<sub>e</sub>, M<sub>tot</sub>, n, b/a)<sub>out</sub>. This allows us to
derive a correction function and discuss possible measurement
biases as outlined below.

4.4.2. Correction Function

We obtain a correction function, *S*(*R*<sub>e</sub>, *M*<sub>tot</sub>, *n*, *b/a*), in an
identical fashion as in Carollo et al. (2013a), and we refer the
reader to this paper for additional details. We construct *S* such
that it returns a four-dimensional median correction vector
(*ΔR*<sub>e</sub>, *ΔM*<sub>tot</sub>, *Δn*, *Δb/a*) for each point in measured
(*R*<sub>e</sub>, *M*<sub>tot</sub>, *n*, *b/a*)<sub>meas</sub> parameter space. The median correction
vector is constructed as the difference between the median
of the 50 closest (*R*<sub>e</sub>, *M*<sub>tot</sub>, *n*, *b/a*)<sub>out</sub> (with respect to
(*R*<sub>e</sub>, *M*<sub>tot</sub>, *n*, *b/a*)<sub>meas</sub>) and the median of their true values
(*R*<sub>e</sub>, *M*<sub>tot</sub>, *n*, *b/a*)<sub>in</sub>. We obtain this correction vector for each
combination of PSF and σ<sub>sky</sub>. The final correction vector is then
obtained by an interpolation of the grid at the PSF and σ<sub>sky</sub>
attributed to the galaxy for which the correction is computed.

Because of our imposed magnitude cut of bright
H = 21.5 AB, the correction in size (usually overestimated)
is on the order of smaller than 20%. The simulations also show
that the detection rate of galaxies is 100% in the worst case up
to half-light sizes of at least 3" at H = 21.5 AB, corresponding
to a surface brightness limit of ~25.2 mag arcsec⁻². This size
corresponds to ~25 kpc (~20 kpc) at z ~ 2 (z ~ 0.5).

The correction function allows an assessment of detection
limits and a first correction for measurement biases. However,
the simulated galaxies are ideal cases. The overlap between
UltraVISTA and CANDELS is ideal to perform a more
thorough calibration of our size measurement.

4.5. Final Calibration of Size Measurements using CANDELS

The second step of our calibration process consists of the
comparison of our measured (and corrected with S) sizes with
HST-based structural measurements on COSMOS/CANDELS,
which has an overlap of 3% with the central part of COSMOS.
Because of the 2.5 times higher resolution and 4 times smaller
PSF of the HST images, we consider the HST-based size
measurements to reflect the true galaxy sizes. We first measure
the sizes of galaxies on the publicly available CANDELS
F160W mosaic as these match the UltraVISTA H-band data
most closely. For this end, we use SExtractor in order to
extract the sources and to obtain the initial parameters for
GALFIT in the same manner as described above for the
UltraVISTA-based measurements. Subsequently, we run
GALFIT for the extracted sources in the two configurations
FIXPOS and VARPOS, thereby applying the same selection
criteria for good fits as described in Section 4.3. Furthermore,
we apply a correction function S as done before, but with
PSF and σsky matching those of the COSMOS/CANDELS
images. In turn, we find corrections smaller than 5% for
galaxies at H < 21.5 AB. As a further check, we compare the
size measurements to the publicly available COSMOS/
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Final calibration of sizes (in arcseconds) using the CANDELS imaging data. Shown are two different magnitude bins below $H = 21.5AB$. The left panels show the comparison of semimajor half-light radii measured on CANDELS and ground-based UltraVISTA images after correction of systematic biases and calibration. Orange and blue points show galaxies with $0.1 < n < 2.5$ and $2.5 < n < 8.0$, the gray regions show 20%, 50%, and 100% discrepancies, and the black lines show a running median with scatter (dashed black lines). Furthermore, the dashed horizontal (vertical) line shows the CANDELS (UltraVISTA) PSF size. The right panels show the normalized histograms of $log(R_{\text{candels}}/R_{\text{ultravista}})$ for uncalibrated (empty) and calibrated (filled/hatched) UltraVISTA size measurements. Without calibrations, the UltraVISTA sizes are substantially underestimated (tail toward large $log(R_{\text{candels}}/R_{\text{ultravista}})$). This is true in particular for compact ($n > 2.5$) galaxies. The histograms indicate an uncertainty for calibrated sizes of ~50%. The gray regions show 20%, 50%, and 100% discrepancies.</figcaption></figure>

CANDELS size catalog by van der Wel et al. (2012) and find excellent agreement.

The comparison between the *HST*-based ($R_{\text{candels}}$) and ground-based ($R_{\text{ultravista}}$) galaxy sizes and their calibration is shown in Figure 3. Shown are galaxies with Sérsic indices $n < 2.5$ (blue) and $n > 2.5$ (orange) measured on the ground-based images in two magnitude bins at $H < 21.5AB$ (top and bottom row). Looking at the *empty* histograms (showing the log-ratio of the sizes) in the right panels, we see an underestimation of $R_{\text{ultravista}}$ by a factor 3 and more, which we find occurs preferentially for galaxies smaller than the (UltraVISTA) PSF radius (~0.′′3) and with large Sérsic $n$ (i.e., compact light distribution). Furthermore, an overestimation of

galaxy sizes preferentially occurs for large galaxies ($R_e > 2′′$) with small Sérsic $n$.

We calibrate our ground-based size measurements by constructing a calibration function $C(R_e, M_{\text{tot}}, n, b/a)$ in a similar way as described in Section 4.4.2. Returning to Figure 3, the measurements with the calibration function applied are shown in the filled and hatched histograms in the right two panels (for different magnitudes and $n$). Furthermore, the left panels show the one-to-one comparison of the size measurements with a running median with $1\sigma$ scatter (dashed). The comparison of the fully calibrated sizes with the *HST*-based size measurements show that we are able to recover $R_e$ on UltraVISTA to an accuracy of better than 50% ($1\sigma$ scatter).
As shown in Figure 3, the uncertainty of the calibrated sizes of galaxies close to the resolution limit of UltraVISTA can be up to a factor of three. We note that fewer than 5% of our massive $log(m/M_\odot) > 11.4$ are unresolved and thus could have much larger uncertainties.

### 4.6. Correction for Internal Color Gradients

Mostly negative internal color gradients are ubiquitously measured in star-forming galaxies up to at least $z \sim 3$, whereas this effect is much weaker in quiescent galaxies (Cassata et al. 2010; Bond et al. 2011; Szomoru et al. 2011; Wuys et al. 2012; Cibinel et al. 2013a; Pastrav et al. 2013; Bond et al. 2014; Hemmati et al. 2014; Vulcani et al. 2014). The observed color gradients are caused by different stellar populations and dust attributed to inside-out growth of galaxies and therefore depend on galaxy age, stellar mass, redshift, and star formation activity. Such color gradients cause the observed size to change as a function of wavelengths. Vice versa, at a fixed observed wavelength the observed size of galaxies changes as a function of redshift since the rest-frame wavelengths shifts. The effect of color gradients may introduce artificial effects in the size evolution across redshift. Several studies have constrained this effect using observations at different wavelengths for different types of galaxies and stellar masses at various redshifts (e.g., Kelvin et al. 2012; van der Wel et al. 2014; Lange et al. 2015). Typical gradients for galaxies at $log(m/M_\odot) = 10$ are on the order of $|\Delta\log R/\Delta\log \lambda| = 0.1-0.3$ depending on data quality, resolution, and redshift. This leads to corrections in size of 10%–50% over a wavelength range of rest-frame 0.5–1.0 $\mu$m. In the following, we use the parameterization by Lange et al. (2015) to correct our size measurements for internal color gradients. However, other parametrizations (e.g., van der Wel et al. 2014) result in similar corrections and do not change the results of this paper.

### 4.7. Verification of the Accuracy of the Size Measurement

Because our measurements at $log(m/M_\odot) > 11.4$ are unique so far, we cannot directly check whether they are reasonable.

In the following, we use our (fully calibrated and mass-complete) low-mass control samples at $10.0 < \log(m/M_\odot) < 11.4$ (see Section 3.1) to investigate possible systematics in our size measurement.

In panels (B) through (D) of Figure 4 we compare our measured size evolution of quiescent (open, color) and star-forming (filled, color) galaxies to measurements taken from the literature (gray lines and symbols; Carollo et al. 2013a; van der Wel et al. 2014). The latter are based on high-resolution HST imaging and corrected for color gradients in the same way as we do here. We find a very good agreement with our measurements.

In panel (A) we compare our final size evolution at $log(m/M_\odot) > 11.4$ to spectroscopically confirmed quiescent galaxies at the same stellar mass in two redshift bins from the literature as black circles (Krogager et al. 2014; Onodera et al. 2015; Belli et al. 2015). These galaxies reside well within the $1-2\sigma$ scatter of our measurements (indicated by the thin error bar), although at the lower end. This can be explained by the higher success rate of spectroscopic surveys for compact galaxies with high surface brightness.

To conclude, we do not expect any severe systematic biases in our measurements.

## 5. Results: Size Evolution of Very Massive Galaxies

### 5.1. Size Evolution of Massive Galaxies

In Figure 4 (panel (A)) we show the final median size evolution with cosmic time of our massive $log(m/M_\odot) > 11.4$ quiescent (red, open) and star-forming (red, filled) galaxies. These are compared to literature measurements at lower masses (Carollo et al. 2013a; van der Wel et al. 2014, gray lines and symbols) and spectroscopically confirmed quiescent galaxies at $z > 1$ (Krogager et al. 2014; Onodera et al. 2015; Belli et al. 2015, median in two redshift bins, black open points). The dashed and solid lines show fits to the size evolution of quiescent and star-forming galaxies, respectively, parametrized as $R_e = B \times (1+z)^{-\beta}$. We find a slope $\beta = 1.22 \pm 0.20$ and $\beta = 1.18 \pm 0.15$ for quiescent and star-forming galaxies with $log(m/M_\odot) > 11.4$, respectively. Note that this slope is statistically identical, in contrast to lower masses, where quiescent galaxies show a faster size increase with cosmic time than star-forming galaxies on average. In addition, very massive star-forming galaxies are only ~20% larger on average at a fixed redshift and stellar mass, whereas at lower masses, the difference can be as large as a factor of two (see gray lines).

### 5.2. The Stellar Mass versus Size Relation

The relation between stellar mass and size (MR relation) has been measured so far on statistically large samples at $log(m/M_\odot) < 11.0$. Our measurement on a large sample of galaxies at $log(m/M_\odot) > 11.4$ enables us to provide an additional data point at high masses. In Figure 5 we show the MR relation in three redshift bins measured over two orders of magnitudes in stellar mass. Shown are our data at $log(m/M_\odot) > 11.4$ (large filled dots) for quiescent (red) and star-forming (blue) galaxies as well as measurements at lower masses. The latter include the 3D-HST survey (cloud of thin points, van der Wel et al. 2014), spectroscopically confirmed quiescent galaxies at $z > 1$ (crosses, asterisks, and pluses, Krogager et al. 2014; Onodera et al. 2015; Belli et al. 2015), and galaxies at $z < 0.1$ with measurements in g band from the Galaxy and Mass Assembly survey (small points with error bars, Driver et al. 2011; Lange et al. 2015). The large blue and red symbols show the median size of star-forming and quiescent galaxies in different stellar mass bins. The lines show the corresponding log-linear fits ($R_e(m) \propto m^\alpha$, see Table 1) to the medians with errors from bootstrapping.

The MR relation of quiescent galaxies is much steeper than for star-forming galaxies. The average sizes of quiescent and star-forming galaxies are comparable at $log(m/M_\odot) \sim 11.5$ independent of redshift. The logarithmic slope of the MR relation ($\langle\alpha_{qu}\rangle \sim 0.6$ for quiescent and $\langle\alpha_{sf}\rangle \sim 0.2$ for star-forming galaxies) does not evolve significantly with cosmic time. It is also very consistent with the measurements in the local universe ($z < 0.1$; Shen et al. 2003; Lange et al. 2015), finding values between $\alpha_{qu} \sim 0.35-0.60$ and $\alpha_{sf} \sim 0.15-0.25$ for quiescent and star-forming galaxies, respectively (see also Figure 5). The study by van der Wel et al. (2014) finds steeper slopes for quiescent galaxies ($\alpha \sim 0.75$), likely due to the authors missing very massive quiescent galaxies at $log(m/M_\odot) > 11.5$. The constant slope of the MR relation is indicative of a constant relation between the growth of galaxies in size (e.g., due to accretion) and stellar mass over cosmic time.
<figure><img src="image_5.png" /></figure>

**Figure 4.** Size evolution as a function of redshift of star-forming and quiescent galaxies at different masses. Panel (A): for quiescent (open, red) and star-forming (filled, red) UMGs at $log(m/M_\odot) > 11.4$ with fits (dashed, solid red lines). Spectroscopically confirmed quiescent UMGs from the literature (Krogager et al. 2014; Onodera et al. 2015; Belli et al. 2015) are shown as black circles (median in two redshift bins with scatter). Clearly, star-forming UMGs are systematically offset to larger observed sizes at all redshifts. However, this offset is smaller compared to lower masses (shown in gray from the literature). The thin error bars show the $1\sigma$ scatter and the thick error bars show the error on the median. Panels (B) through (D): comparison of our size measurements (color points, up to redshift where the mass is complete) with measurements in the literature in gray for $11.0 < log(m/M_\odot) < 11.4$, $10.5 < log(m/M_\odot) < 11.0$, and $10.0 < log(m/M_\odot) < 10.5$, respectively. The lines (dashed: quiescent, solid: star-forming) are fits to the size evolution by van der Wel et al. (2014). The filled (open) symbols show measurements of star-forming (quiescent) galaxies from Carollo et al. (2013a) (diamonds) and van der Wel et al. (2014) (triangles). The good agreement to our measurements shows that our measurement are not biased.

Finally, we note that a recent study by Peng et al. (2015) suggests that the bulk of star-forming $z \sim 0$ galaxies at $log(m/M_\odot) < 11$ are being quenched via strangulation<fnref n="9" /> within ~4 Gyr. We would therefore expect the $m - R_e$ relation of star-forming galaxies at $z \sim 0.5$ and $log(m/M_\odot) < 11.0$ to be similar to the relation of the quiescent galaxies at $z \sim 0$ if the observed sizes of the galaxies do not change during or after quenching. This is, however, not seen from the left panel of Figure 5, which shows that star-forming galaxies ~4 Gyr ago are significantly (up to a factor two) larger than local quiescent galaxies at $log(m/M_\odot) < 11.0$. This contradiction can be alleviated by post-quenching disk-fading, which would substantially decrease the observed sizes of quiescent galaxies and is shown to be at work at low redshifts (Carollo et al. 2016) and most likely also at $z \sim 2$ (Tacchella et al. 2015, 2016a). In

addition to this, at high redshifts, morphological transformation as a result of quenching cannot be ruled out.

**6. Model for the Size Evolution of Massive Quiescent Galaxies**

The similar sizes at a given redshift of star-forming and quiescent galaxies at $log(m/M_\odot) > 11.0$ at all redshifts $z < 2$ suggest a very close connection of these galaxies. This might have important implications on the process that quenches these galaxies. In this section, we investigate this further by modeling the size evolution of quiescent galaxies, thereby applying different assumptions on the quenching process. A quantitative visualization of our model is shown in Figure 6.

**6.1. Evolution of Star-forming Model Galaxies**

Our model assumes that galaxies—as long as they are forming stars—evolve along the star-forming main-sequence (MS) spanned by stellar mass and SFR. In addition, we assign

<footnote n="9">⁹ Strangulation means that the supply of cold gas is halted and thus star formation is shut down.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Mass vs. size relation for quiescent (red) and star-forming (blue) galaxies for the combined sample of 3D HST at log(*m*/M<sub>☉</sub>) < 11.5 (thin dots) and our sample of massive galaxies with log(*m*/M<sub>☉</sub>) > 11.4 (large points). Also included are spectroscopically confirmed quiescent galaxies from Onodera et al. (2015, pluses), Krogager et al. (2014, crosses), and Belli et al. (2015, asterisks) at various masses, as well as *z* < 0.1 galaxies from Lange et al. (2015, small dots with error bars). The medians in mass bins (large points with error bars) are fit linearly in log-space for three redshift bins using a least-squares method (red and blue points with error bars indicating the error on the median). The logarithmic slope α (*R*<sub>e</sub> ∝ *m*<sup>α</sup>) is indicated for star-forming and quiescent galaxies. We note the large dispersion of sizes at high masses, consistent with the findings of Mancini et al. (2010).</figcaption></figure>

<table><caption>**Table 1**
Power-law Slope (*R*<sub>e</sub>(*m*) ∝ *m*<sup>α</sup>) of the Stellar Mass vs. Size Relation at *z* > 0.5 (this Work including 3D HST and Spectroscopically Confirmed Quiescent Galaxies) as well as Integrated Over Cosmic Time at Lower Redshifts (see References)</caption><thead><tr><th>Redshift range</th><th>Star-forming<br>α<sub>sf</sub></th><th>Quiescent<br>α<sub>qu</sub></th><th>References</th></tr></thead><tbody><tr><td>z ~ 0</td><td>0.14 to 0.39</td><td>0.56</td><td>(1)</td></tr><tr><td>z < 0.1</td><td>0.19 ± 0.02</td><td>0.41 ± 0.06</td><td>(2)</td></tr><tr><td>0.5 < z < 1.0</td><td>0.30 ± 0.10</td><td>0.55 ± 0.05</td><td>This work</td></tr><tr><td>1.0 < z < 1.5</td><td>0.22 ± 0.08</td><td>0.62 ± 0.09</td><td>This work</td></tr><tr><td>1.5 < z < 2.0</td><td>0.14 ± 0.06</td><td>0.59 ± 0.15</td><td>This work</td></tr><tr><td>z < 3</td><td>0.22 ± 0.05</td><td>0.75 ± 0.05</td><td>(3)</td></tr></tbody></table>

**Notes.** (1) Shen et al. (2003). For star-forming galaxies they fit α = 0.14 at log(*m*/M<sub>☉</sub>) < 10.6 and α = 0.39 at log(*m*/M<sub>☉</sub>) > 10.6. (2) Lange et al. (2015). (3) van der Wel et al. (2014). Report no significant change in slope over 0 < z < 3.

our model galaxies a half-light radius *R*<sub>e</sub> and a gas fraction *f*<sub>gas</sub> using empirical relations and observations. The galaxies eventually become quenched in concordance with the observed quiescent fraction observed as a function of redshift and stellar mass. As fiducial model for quenching we assume an instantaneous (on-the-spot) quenching process without any change in the structure (i.e., half-light size *R*<sub>e</sub>) of the galaxies. We complement this model with two additional models featuring a structural change (compaction due to starburst) as well as a delayed quenching. These different assumptions of quenching processes are explained in more detail later on.

The main steps of our empirical model are the following.

1. Our model starts at *z* = 2.5 and uses the observed stellar mass function by Ilbert et al. (2013) as initial condition for the mass distribution of the 100,000 simulated star-forming galaxies with stellar masses between 7 < log (*m*/M<sub>☉</sub>) < 12. The initial mass distribution and fraction of quiescent galaxies is derived from the quiescent fraction *f*<sub>q</sub>(*m*, *z*) at a given redshift and stellar mass (Figure 7).

2. We evolve the stellar mass and SFR of star-forming galaxies along their MS, for which we use the

parameterization by Schreiber et al. (2015) compiled from deep Herschel observations. We verified that the use of other parameterizations of the MS does not change the results and conclusions of this work. Furthermore, we assign to each of the star-forming model galaxies gas fractions *f*<sub>gas</sub>(*m*, *z*) from our compilation of the literature, as outlined in the Appendix, as well as sizes according to the measured size versus stellar mass relation *R*<sub>e</sub>(*m*, *z*) for star-forming galaxies (including our new measurements at log(*m*/M<sub>☉</sub>) > 11.0). When drawing values from the above empirical relations, we also include the observed scatter, which we characterized by a Gaussian centered on the median. The typical scatter for the MS and MR relation is assumed to be ~0.3 dex.

3. At each redshift, we quench galaxies in mass bins randomly, such that the model quiescent fraction reproduces the observed *f*<sub>q</sub>(*m*, *z*). After a galaxy is quenched, we set its gas fraction and SFR to zero. The remaining gas is added instantaneously to the stellar mass under the simple assumption that the gas is fully converted into stars and not stripped. Depending on the quenching model (see below), the new stars are either distributed evenly in the galaxy disk or in a central region of 1 kpc. We also do not implement the rejuvenation of galaxies once they are quiescent.

For visualization purposes, example tracks of our model galaxies with different initial stellar masses at *z* = 2.5 as well as the fraction of quiescent galaxies (*f*<sub>q</sub>) are shown in Figure 6.

## 6.2. Quenching of Model Galaxies

We implement three simple models that should bracket different pathways of quenching processes. In the following, we describe these in more detail.

**Instantaneous/no structural change.** This is our fiducial model in which galaxies quench instantaneously without any structural change. A physical scenario could be the cutoff of cold gas inflow by the heating of the gas in massive dark matter halos above 10<sup>12</sup> M<sub>☉</sub> (e.g., Croton et al. 2006). The galaxy then
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Tracks of our model galaxies on the SFR vs. stellar mass (left) and SFR vs. redshift (right) plane. We show the evolution of three sets of galaxies with different initial stellar masses at $z = 2.5$ (log(*m*/$M_\odot$) = 9, 10, and 11 shown in green circles, orange squares, and blue diamonds, respectively). The large symbols show the SFR in redshift bins with $\Delta z \sim 0.3$. The redshifts are indicated as numbers in the symbols in the left panel, and the median stellar mass in each redshift bin is shown as numbers in the lower part of the right panel. The error bars represent $1\sigma$ errors. The lines show five randomly chosen individual galaxies in each of the sets with different initial masses. The insets show the quiescent fraction $f_q$ as a function of stellar mass marginalized over all redshifts (left panel) and as a function of redshift marginalized over all stellar masses (right panel).</figcaption></figure>

consumes its remaining gas according to its SFR on the star-forming MS and evenly increases its mass in its disk.

**Instantaneous/compaction.** In this model, the galaxy decreases its overall size (compaction) due to an increase of its surface density instantaneously after the shutdown of star formation. We assume that this compaction is triggered by a starburst in the inner 1 kpc region of the galaxy, which may be induced by a major merger event (e.g., Barro et al. 2013). We compute the decrease in overall half-light radius after the starburst by adding the gas of the galaxy in a 1 kpc bulge component characterized by a $n = 4$ Sérsic profile to the disk dominated ($n = 1$) star-forming galaxy. For simplicity we assume that all of the gas mass is turned into stars in the bulge component. Furthermore, we assume that the bulge component has the same mass-to-light ratio as the disk, such that the ratio in luminosity of the bulge component and the disk is proportional to the ratio of stellar mass added to the bulge and stellar mass in the disk.

**Delayed/no structural change.** This model is similar to our fiducial model, with the difference that the quenching does not occur instantaneously, but with a delay. A possible scenario could be the slow consumption of gas off the star-forming MS after the gas supply onto the galaxy is cut off. We assume the delay (i.e., the time the galaxy spends in the green valley) to be 50% of the cosmic time between the start of quenching and $z = 2.5$ (~800 Myr at $z = 1.5$ and ~3 Gyr at $z = 0.5$).

## 7. Discussion

We now compare the predicted size evolution of quiescent galaxies from our simple empirical models with observations to investigate possible processes that quench massive galaxies at $z < 2$. Figure 8 shows the MR relation of our quiescent model galaxies (symbols) together with the observed relations for star-forming (blue, hatched and dashed line) and quiescent (red, hatched and solid line), respectively. The width of the hatched bands and the error bars on the points represent the scatter in

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Fraction of quiescent galaxies, $f_q(m, z)$, as a function of redshift for different bins in stellar mass (color). The lines show continuous fractions derived from the ratio of the stellar mass function of quiescent and star-forming galaxies from Ilbert et al. (2013). The symbols show the fraction derived from mass-complete samples in UltraVISTA.</figcaption></figure>

the observed as well as modeled relations. The different quenching models are shown in different colors and symbols, as indicated in the legend.

Focusing on galaxies above the characteristic knee of the stellar mass function (log($M_*$) ~ 10.8, e.g., Ilbert et al. 2013; Davidzon et al. 2017), we note the following.

(i) The instantaneous quenching model without altering the structure of the galaxies (green circles) predicts the quiescent MR relation well at all masses at $z > 0.5$.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Observed MR relation in three redshift bins of star-forming (blue, hatched) and quiescent (red, hatched) galaxies together with our three different quenching models (symbols). The filled symbols show the models without mergers, and the open symbols assume that 90% of the quiescent galaxies above $log(M_*) \sim 10.8$ experience ten 1:10 minor mergers after being quenched. The models suggest a fast quenching of massive ($log(m/M_\odot) > 11$) galaxies at $z > 1$. In order to reproduce the MR relation at lower redshifts, a series of minor mergers is necessary. The same is true for quenching and compaction via major mergers.</figcaption></figure>

However, this model underpredicts the sizes of the most massive galaxies ($log(m/M_\odot) \sim 11.5$) at $z < 0.5$.

(ii) The instantaneous quenching model followed by a compaction triggered by a starburst within a 1 kpc central region (orange squares) predicts the sizes of $log(m/M_\odot) < 11$ galaxies well down to $z \sim 0.5$, but underpredicts the sizes at later times as well as at higher masses.

(iii) The delayed quenching model without structural change (blue diamonds) is only able to explain the MR relation at $log(m/M_\odot) > 11$ at the highest redshifts, but underpredicts the sizes at lower redshifts by factors of 2–3. It reproduces the relations at $log(m/M_\odot) < 11$ and $z < 1$ well.

We explain and interpret these findings in more detail in the following subsections.

## 7.1. Slow versus Fast Quenching at $m > M^*$

The star formation in very massive galaxies can be shut down without significant structural change of the light profile by cutting off the gas supply onto the galaxies. In current theoretical models and simulations, this can be achieved in galaxies with massive dark matter halos of $m_{DM} > 10^{12} M_\odot$, which cause the infalling gas to be heated (e.g., Croton et al. 2006). This may result in a uniform decrease of the star formation in the galaxy disk without significantly altering its structure. Taking the above results at face value suggests that if there is no net structural change after the turnoff of star formation, massive galaxies ($log(m/M_\odot) > 11$) have to transition from star-forming to quiescent on relatively short timescales. This is suggested by the fact that our fiducial model (instantaneous quenching) is able to reproduce the sizes of galaxies at these stellar masses reasonably well, at least in the two upper redshift bins at $z \gtrsim 0.5$. An instantaneous quenching might be too much of a simplification, and a non-zero quenching time is suggested by recent observational studies (e.g., Schawinski et al. 2014; Peng et al. 2015). Our delayed quenching model works well for redshifts $z \gtrsim 1$ and $log(m/M_\odot) > 11$, where the delay times are shorter than 1–1.5 Gyr according to our definition (50% of the difference in cosmic time between the quenching event and $z = 2.5$). Note

that this is compatible with the time a galaxy on the star-forming MS needs to consume all of its gas given its main-sequence gas fraction and SFR: less than 1–2 Gyr for a galaxy at $log(m/M_\odot) > 11$ and $z > 0.5$ (e.g., Tacconi et al. 2017). Note that our delayed model overpredicts the sizes of galaxies at $m < M^*$ at $z > 1$. This would suggest that the delay as defined here is not long enough, and instead a longer delay (2 Gyr or more) is favored. This mass dependence of the quenching time (i.e., slow versus fast quenching) is also strongly suggested by recent simulations (e.g., Hahn et al. 2016) and could be explained by different quenching processes taking place at different stellar masses and as a function of environment the galaxies are living in.

## 7.2. Merger-induced Starbursts and Compaction

It is suggested that mergers play an important role in shaping galaxies at high stellar masses. Thus a smooth quenching without significantly altering the structure of a galaxy is likely too simplistic. Our model of a merger-triggered compact starburst inducing a fast consumption of gas and quiescence might therefore be a better approach to characterize the quenching mechanism at high redshifts and high stellar masses. As shown in Figure 8, such a scenario underpredicts the sizes of massive quiescent galaxies by factors of two or more at all redshifts. If such a scenario is the dominant way of quenching massive galaxies, then the galaxies have to grow individually to meet the observed MR relation. This is similar to the fast-track quenching mechanism proposed by Barro et al. (2013) (see also Zolotov et al. 2015), in which galaxies experience a compaction phase with subsequent growth due to minor and major mergers. Note that changing the parameters of this particular model does not significantly change this conclusion. For example, assuming a 2 kpc central starburst would only increase the sizes by ~50% and still leads to a significant underprediction.

## 7.3. Post-quenching Growth through Mergers in Massive Galaxies

Figure 8 shows that all of our bracketing models in some cases severely underpredict the sizes of quiescent galaxies
above $M^*$ and $z \lesssim 1$. One possible way to bring the models in agreement with observations is to introduce a series of minor and/or major mergers following the quenching event. We investigate this further by assuming a simple toy model in which 90% of the quiescent galaxies above $\log(M_*) \sim 10.8$ experience ten 1:10 minor mergers during their lives after being quenched. We choose this case because minor mergers are more common and are dominantly increasing the size of galaxies and less so their stellar mass. For the implementation of this model, we assume that the virial condition holds for gas-poor ellipticals and compute the resulting size increase ($\Delta R_e$) as a function of the merger mass fraction ($\Delta m$) and change in velocity dispersion ($\Delta \sigma$) during the merger event as

$$ \Delta R_e = \Delta m \left( \frac{1}{\Delta \sigma} \right)^2, \qquad (3) $$

where $\Delta R_e = \frac{R_{e,post}}{R_{e,pre}}$, $\Delta m = \frac{m_{post}}{m_{pre}}$, and $\Delta \sigma = \frac{\sigma_{post}}{\sigma_{pre}}$ are the ratios of quantities before (“pre”) and after (“post”) the merging event. We assume that the change in the velocity dispersion is negligible during the merger event, i.e., $\Delta \sigma \sim 1$ (e.g., Hopkins et al. 2009; Oser et al. 2012).

The open symbols in Figure 8 show the impact of post-quenching mergers on our previous results. The addition of a series of minor mergers to our instantaneous quenching + compaction model (orange open squares) leads indeed to a good agreement with the observed MR relation at $z > 0.5$ at all stellar masses probed here. We note, however, that the sizes of massive $m > M^*$ galaxies are still underestimated at $z < 0.5$. It is therefore likely that if the compaction model holds, these galaxies must experience more minor mergers cumulatively than anticipated in our simple merger toy model. Alternatively, massive galaxies at later cosmic times might be quenched via other paths that do not include a compaction phase, such as heating of cold gas alone. Such a possibility is shown by our fiducial model with post-quenching minor mergers (green open circles), which is able to predict the sizes of massive galaxies at $z < 0.5$ well (however, it fails at high redshifts). We note that Equation (3) describes the effect of size growth by major mergers. Instead, strictly speaking, the size growth due to minor mergers is expected to be steeper ($\Delta R \propto \Delta m^\alpha$ with $\alpha > 1$, e.g., Bezanson et al. 2009; Naab et al. 2009), which would decrease the number of mergers needed in our model. For example, assuming $\alpha = 2$, we find that only $\sim 30\%$ of the galaxies at $\log(m/M_\odot) > 10.8$ are needed to experience a 1:10 merger in order to meet the observations.

Finally, we note that mergers for low-redshift ($z < 1$) less massive galaxies ($m < M^*$) are not needed to bring our models in agreement with observations. This is in line with the idea that the size evolution of quiescent galaxies below $\log(m/M_\odot) \sim 11$ with cosmic time is mainly driven by the addition of newly quenched galaxies, while at higher masses it is more dominated by individual growth due to mergers (e.g., Carollo et al. 2016; Belli et al. 2015).

## 8. Summary and Conclusions

We use the size evolution of massive star-forming and quiescent galaxies as an independent diagnostic tool to investigate the process of quenching at $\log(m/M_\odot) > 11$ and $z \lesssim 2$. To this end, we measure the half-light size evolution of a large sample of very massive star-forming and quiescent galaxies at $\log(m/M_\odot) \gtrsim 11.4$ on the 2 square degree survey field of

COSMOS/UltraVISTA. We find the size evolution of both populations of galaxies at $\log(m/M_\odot) > 11.4$ to be similar in slope and normalization and to be consistent with the extrapolation of the mass versus size relation from lower masses.

In order to investigate different quenching mechanisms and the impact of mergers, we predict the MR relation of massive $m > M^*$ quiescent galaxies within our simple empirical models as a function of redshift. Our main results are the following.

1. Massive galaxies quench fast. Models with instantaneous quenching or a short delay of up to $\sim 1$ Gyr are able to predict the sizes of quenching galaxies at $z > 1$ and $m > M^*$. Longer quenching times are more favored at lower masses and redshifts.

2. A more realistic model incorporating a compaction phase (e.g., due to a merger-triggered central starburst within 1 kpc) followed by quiescence and subsequent individual growth by mergers is able to reproduce the observed MR relation of massive $m > M^*$ quiescent galaxies at all redshifts.

3. None of our models is able to predict the size evolution of $m > M^*$ galaxies at low redshifts ($z \lesssim 1$). We show that with 1:10 minor mergers for 90% of the quiescent galaxies at $m > M_*$, the models can be brought into agreement with observations. In contrast, no mergers are needed at lower stellar masses, in agreement with the size evolution being driven by the addition of larger newly quenched galaxies.

It is important to note that we are not able to distinguish the dominant pathways of quenching of massive quiescent galaxies with our simple models as these yield very similar predictions for the size evolution. Nonetheless, our study suggests that quenching is likely a fast process at the stellar masses probed here, with a significant involvement of mergers in the post-quenching growth of massive galaxies. To further distinguish these models, more information on the (resolved) structural properties of the galaxies is necessary. This will be possible with high-resolution imaging and spectroscopy of massive quiescent galaxies by the HST or the James Webb Space Telescope.

We would like to thank Dan Masters, Charles Steinhardt, Behnam Darvish, and Bahram Mobasher for valuable discussions. Furthermore, we would like to thank the referee for valuable feedback that greatly improved this manuscript. A.F. acknowledges support from the Swiss National Science Foundation. Based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under ESO programme ID 179.A-2005 and on data products produced by TERAPIX and the Cambridge Astronomy Survey Unit on behalf of the UltraVISTA consortium. This work is based on observations taken by the CANDELS Multi-Cycle Treasury Program with the NASA/ESA HST, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555.

## Appendix

### The Gas Fraction $f_{gas}(m, z)$

We use studies from the literature to fit an empirical relation $f_{gas}(m, z)$, which is used in our models. The data used include PHIBSS at $z \sim 1-1.5$ (Tacconi et al. 2013) and COLDGASS at $z \sim 0$ (Saintonge et al. 2011), as well as data from lensed and
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Gas fraction as a function of redshift in the four different stellar mass bins (different panels). The black and gray symbols show literature values measured from individual galaxies (black squares: PHIBSS at $z \sim 1-1.5$. Tacconi et al. (2013); black diamonds: COLDGASS at $z \sim 0$, Saintonge et al. (2011); gray points: lensed and other star-forming galaxies from Dessauges-Zavadsky et al. (2014) and references therein). The hatched region is a fit ($f_{gas}(m, z)$) to the measured individual galaxies including the uncertainty as described in the text. For comparison, the color symbols show $f_{gas}$ derived from the Kennicutt-Schmidt relation for the UltraVISTA galaxies. Note that these are not used in the fitting.</figcaption></figure>

other star-forming galaxies from Dessauges-Zavadsky et al. (2014) and references therein.

The result is shown in Figure 9 for four different bins in stellar mass. The PHIBSS and COLDGASS data are shown in black, the other measurements are shown in gray. We also show in color $f_{gas}$ derived from our galaxies (UMGs and lower mass control sample) using the Kennicutt-Schmidt relation (KS relation, Schmidt 1959; Kennicutt 1998), relating $\Sigma_{gas} \propto \Sigma_{SFR}^N$, where we take $N = 1.31$ (Krumholz et al. 2012). Note that these derivations are not used to fit the parametrization for $f_{gas}(m, z)$.

We derive $f_{gas}(m, z)$ and its uncertainty (95% CLs) by fitting the observed data as follows. We first perform a linear fit forced through the COLDGASS data point at $z = 0$ in order to determine

the slope. The error on the slope is derived from the systematic error of the fit and the uncertainty of the data points by bootstrapping, which we both add in quadrature. In a second fit, we fix the slope to the one determined before and fit for the intercept including error. The resulting uncertainty region as shown in Figure 9 as hatched region is then derived by the unification of the errors of the two fits. To obtain a continuous function for $f_{gas}$, we interpolate between the four stellar mass bins.

We compared our fit to the recent work by Genzel et al. (2015). We find that their $f_{gas}(m, z)$ parametrization has a slightly steeper redshift dependence, resulting in 10%–30% larger gas fractions at the highest redshifts. We have verified that our results do not change when we use the Genzel et al. (2015) parametrization for $f_{gas}(m, z)$.
References

Arnouts, S., Moscardini, L., Vanzella, E., et al. 2002, MNRAS, 329, 355

Arnouts, S., Walcher, C. J., Le Fèvre, O., et al. 2007, A&A, 476, 137

Baldry, I. K., Glazebrook, K., Brinkmann, J., et al. 2004, ApJ, 600, 681

Barro, G., Faber, S. M., Pérez-González, P. G., et al. 2013, ApJ, 765, 104

Behroozi, P. S., Marchesini, D., Wechsler, R. H., et al. 2013, ApJL, 777, L10

Belli, S., Newman, A. B., & Ellis, R. S. 2015, ApJ, 799, 206

Belli, S., Newman, A. B., Ellis, R. S., & Konidaris, N. P. 2014, ApJL, 788, L29

Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393

Bezanson, R., van Dokkum, P. G., Tal, T., et al. 2009, ApJ, 697, 1290

Birnboim, Y., & Dekel, A. 2003, MNRAS, 345, 349

Birnboim, Y., Dekel, A., & Neistein, E. 2007, MNRAS, 380, 339

Bond, N. A., Gardner, J. P., de Mello, D. F., et al. 2014, ApJ, 791, 18

Bond, N. A., Gawiser, E., & Koekemoer, A. M. 2011, ApJ, 729, 48

Bournaud, F., Elmegreen, B. G., & Elmegreen, D. M. 2007, ApJ, 670, 237

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Buitrago, F., Trujillo, I., Conselice, C. J., et al. 2008, ApJ, 687, L61

Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682

Cameron, E., & Driver, S. P. 2007, MNRAS, 377, 523

Capak, P., Aussel, H., Ajiki, M., et al. 2007, ApJS, 172, 99

Carollo, C. M., Bschorr, T. J., Renzini, A., et al. 2013a, ApJ, 773, 112

Carollo, C. M., Cibinel, A., Lilly, S. J., et al. 2013b, ApJ, 776, 71

Carollo, C. M., Cibinel, A., Lilly, S. J., et al. 2016, ApJ, 818, 180

Cassata, P., Giavalisco, M., Guo, Y., et al. 2010, ApJL, 714, L79

Cassata, P., Giavalisco, M., Williams, C. C., et al. 2013, ApJ, 775, 106

Cen, R. 2014, ApJ, 781, 38

Chabrier, S. 2003, Archives de Pédiatrie, 10, 168

Cibinel, A., Carollo, C. M., Lilly, S. J., et al. 2013a, ApJ, 777, 116

Cibinel, A., Carollo, C. M., Lilly, S. J., et al. 2013b, ApJ, 776, 72

Croton, D. J., Springel, V., White, S. D. M., et al. 2006, MNRAS, 365, 11

Daddi, E., Dickinson, M., Morrison, G., et al. 2007, ApJ, 670, 156

Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, arXiv:1701.02734

De Lucia, G., Weinmann, S., Poggianti, B. M., Aragón-Salamanca, A., & Zaritsky, D. 2012, MNRAS, 423, 1277

Dekel, A., & Burkert, A. 2014, MNRAS, 438, 1870

Dessauges-Zavadsky, M., Zamojski, M., Schaerer, D., et al. 2014, arXiv:1408.0816

Driver, S. P., Hill, D. T., Kelvin, L. S., et al. 2011, MNRAS, 413, 971

Fagioli, M., Carollo, C. M., Renzini, A., et al. 2016, ApJ, 831, 173

Feldmann, R., Carollo, C. M., & Mayer, L. 2011, ApJ, 736, 88

Feldmann, R., Carollo, C. M., Porciani, C., Lilly, S. J., & Oesch, P. 2008, arXiv:0801.3275

Franx, M., van Dokkum, P. G., Schreiber, N. M. F., et al. 2008, ApJ, 688, 770

Genzel, R., Tacconi, L. J., Lutz, D., et al. 2015, ApJ, 800, 20

Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35

Gunn, J. E., Hoessel, J. G., & Oke, J. B. 1986, ApJL, 306, 30

Hahn, C., Tinker, J. L., & Wetzel, A. R. 2016, arXiv:1609.04398

Hearin, A. P., Zentner, A. R., Berlind, A. A., & Newman, J. A. 2013, MNRAS, 433, 659

Hemmati, S., Miller, S. H., Mobasher, B., et al. 2014, ApJ, 797, 108

Hirschmann, M., De Lucia, G., Wilman, D., et al. 2014, MNRAS, 444, 2938

Hopkins, P. F., Hernquist, L., Cox, T. J., Keres, D., & Wuyts, S. 2009, ApJ, 691, 1424

Huang, K.-H., Ferguson, H. C., Ravindranath, S., & Su, J. 2013, ApJ, 765, 68

Ilbert, O., Arnouts, S., McCracken, H. J., et al. 2006, A&A, 457, 841

Ilbert, O., McCracken, H. J., Le Fèvre, O., et al. 2013, A&A, 556, A55

Ilbert, O., Salvato, M., Le Floc'h, E., et al. 2010, ApJ, 709, 644

Kawata, D., & Mulchaey, J. S. 2008, ApJL, 672, L103

Kelvin, L. S., Driver, S. P., Robotham, A. S. G., et al. 2012, MNRAS, 421, 1007

Kennicutt, R. C., Jr. 1998, ApJ, 498, 541

Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36

Kriek, M., van Dokkum, P. G., Franx, M., Illingworth, G. D., & Magee, D. K. 2009, ApJL, 705, L71

Krogager, J.-K., Zirm, A. W., Toft, S., Man, A., & Brammer, G. 2014, ApJ, 797, 17

Krumholz, M. R., Dekel, A., & McKee, C. F. 2012, ApJ, 745, 69

Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24

Lange, R., Driver, S. P., Robotham, A. S. G., et al. 2015, MNRAS, 447, 2603

Leauthaud, A., Massey, R., Kneib, J.-P., et al. 2007, ApJS, 172, 219

Lilly, S. J., & Carollo, C. M. 2016, ApJ, 833, 1

Mancini, C., Daddi, E., Renzini, A., et al. 2010, MNRAS, 401, 933

Mandelker, N., Dekel, A., Ceverino, D., et al. 2014, MNRAS, 443, 3675

Marchesini, D., Muzzin, A., Stefanon, M., et al. 2009, ApJ, 794, 65

Marmo, C., & Bertin, E. 2008, in ASP Conf. Ser. 394, Astronomical Data Analysis Software and Systems XVII, ed. R. W. Argyle, P. S. Bunclark, & J. R. Lewis (San Francisco, CA: ASP), 619

Martig, M., Bournaud, F., Teyssier, R., & Dekel, A. 2009, ApJ, 707, 250

Martin, D. C., Wyder, T. K., Schiminovich, D., et al. 2007, ApJS, 173, 342

McCracken, H. J., Milvang-Jensen, B., Dunlop, J., et al. 2012, A&A, 544, A156

Mo, H. J., Mao, S., & White, S. D. M. 1998, MNRAS, 295, 319

Moffat, A. F. J. 1969, A&A, 3, 455

Mok, A., Balogh, M. L., McGee, S. L., et al. 2013, MNRAS, 431, 1090

Mosleh, M., Tacchella, S., Renzini, A., et al. 2017, ApJ, 837, 2

Mosleh, M., Williams, R. J., Franx, M., et al. 2012, ApJL, 756, L12

Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJ, 777, 18

Muzzin, A., van der Burg, R. F. J., McGee, S. L., et al. 2014, ApJ, 796, 65

Naab, T., Johansson, P. H., & Ostriker, J. P. 2009, ApJ, 699, L178

Newman, A. B., Ellis, R. S., Bundy, K., & Treu, T. 2012, ApJ, 746, 162

Noeske, K. G., Weiner, B. J., Faber, S. M., et al. 2007, ApJL, 660, L43

Oesch, P. A., Bouwens, R. J., Carollo, C. M., et al. 2010a, ApJL, 709, L21

Oesch, P. A., Carollo, C. M., Feldmann, R., et al. 2010b, ApJL, 714, L47

Onodera, M., Carollo, C. M., Renzini, A., et al. 2015, ApJ, 808, 161

Onodera, M., Carollo, C. M., et al. 2012, ApJ, 755, 26

Oser, L., Naab, T., Ostriker, J. P., & Johansson, P. H. 2012, ApJ, 744, 63

Pastrav, B. A., Popescu, C. C., Tuffs, R. J., & Sansom, A. E. 2013, A&A, 553, A80

Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2010a, AJ, 139, 2097

Peng, Y., Maiolino, R., & Cochrane, R. 2015, Natur, 521, 192

Peng, Y.-j., Lilly, S. J., Kovač, K., et al. 2010b, ApJ, 721, 193

Poggianti, B. M., Moretti, A., Calvi, R., et al. 2013, ApJ, 777, 125

Pozzetti, L., Bolzonella, M., Zucca, E., et al. 2010, A&A, 523, A13

Saintonge, A., Kauffmann, G., Kramer, C., et al. 2011, MNRAS, 415, 32

Saracco, P., Casati, A., Gargiulo, A., et al. 2014, A&A, 567, A94

Saracco, P., Longhetti, M., & Gargiulo, A. 2011, MNRAS, 412, 2707

Schawinski, K., Urry, C. M., Simmons, B. D., et al. 2014, MNRAS, 440, 889

Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521

Schechter, P. 1976, ApJ, 203, 297

Schmidt, M. 1959, ApJ, 129, 243

Schreiber, C., Pannella, M., Elbaz, D., et al. 2015, A&A, 575, A74

Scoville, N., Aussel, H., Brusa, M., et al. 2007, ApJS, 172, 1

Shen, S., Mo, H. J., White, S. D. M., et al. 2003, MNRAS, 343, 978

Shibuya, T., Ouchi, M., & Harikane, Y. 2015, ApJS, 219, 15

Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163

Stockton, A., McGrath, E., Canalizo, G., Iye, M., & Maihara, T. 2008, ApJ, 672, 146

Szomoru, D., Franx, M., Bouwens, R. J., et al. 2011, ApJL, 735, L22

Szomoru, D., Franx, M., & van Dokkum, P. G. 2012, ApJ, 749, 121

Tacchella, S., Carollo, C. M., Renzini, A., et al. 2015, Sci, 348, 314

Tacchella, S., Dekel, A., Carollo, C. M., et al. 2016a, MNRAS, 458, 242

Tacchella, S., Dekel, A., Carollo, C. M., et al. 2016b, MNRAS, 457, 2790

Taconi, L. J., Genzel, R., Saintonge, A., et al. 2017, arXiv:1702.01140

Taconi, L. J., Neri, R., Genzel, R., et al. 2013, ApJ, 768, 74

Taranu, D. S., Hudson, M. J., Balogh, M. L., et al. 2014, MNRAS, 440, 1934

Toft, S., van Dokkum, P., Franx, M., et al. 2007, ApJ, 671, 285

Trinh, C. Q., Barton, E. J., Bullock, J. S., et al. 2013, MNRAS, 436, 635

Trujillo, J., Aguerri, J. A. L., Cepa, J., & Gutiérrez, C. M. 2001, MNRAS, 328, 977

van der Wel, A., Bell, E. F., Häussler, B., et al. 2012, ApJS, 203, 24

van der Wel, A., Franx, M., van Dokkum, P. G., et al. 2014, ApJ, 788, 28

von der Linden, A., Wild, V., Kauffmann, G., White, S. D. M., & Weinmann, S. 2010, MNRAS, 404, 1231

Vulcani, B., Bamford, S. P., Häußler, B., et al. 2014, MNRAS, 441, 1340

Wetzel, A. R., Tinker, J. L., Conroy, C., & van den Bosch, F. C. 2013, MNRAS, 432, 336

Whitaker, K. E., Kriek, M., van Dokkum, P. G., et al. 2012, ApJ, 745, 179

Williams, C. C., Giavalisco, M., Bezanson, R., et al. 2017, ApJ, 838, 94

Williams, R. J., Quadri, R. F., Franx, M., et al. 2010, ApJ, 713, 738

Williams, R. J., Quadri, R. F., Franx, M., van Dokkum, P., & Labbé, I. 2009, ApJ, 691, 1879

Wuyts, S., Förster Schreiber, N. M., Genzel, R., et al. 2012, ApJ, 753, 114

Zolotov, A., Dekel, A., Mandelker, N., et al. 2015, MNRAS, 450, 2327