# Stellar Mass Growth of Brightest Cluster Galaxy Progenitors in COSMOS Since $z \sim 3$

Kevin C. Cooke, Jeyhan S. Kartaltepe, K. D. Tyler, Behnam Darvish, Caitlin

M. Casey, Olivier Le Fèvre, Mara Salvato, Nicholas Scoville

► To cite this version:

Kevin C. Cooke, Jeyhan S. Kartaltepe, K. D. Tyler, Behnam Darvish, Caitlin M. Casey, et al. Stellar Mass Growth of Brightest Cluster Galaxy Progenitors in COSMOS Since $z \sim 3$. The Astrophysical Journal, 2019, 881, 10.3847/1538-4357/ab30c9. insu-03666736

HAL Id: insu-03666736

https://insu.hal.science/insu-03666736v1

Submitted on 13 Aug 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Stellar Mass Growth of Brightest Cluster Galaxy Progenitors in COSMOS Since $z \sim 3$

Kevin C. Cooke¹, Jeyhan S. Kartaltepe¹, K. D. Tyler¹, Behnam Darvish², Caitlin M. Casey³, Olivier Le Fèvre⁴, Mara Salvato⁵,⁶, and Nicholas Scoville²

¹School of Physics and Astronomy, Rochester Institute of Technology, Rochester, NY 14623, USA; kcc7952@rit.edu, jeyhan@astro.rit.edu

²Cahill Center for Astrophysics, California Institute of Technology, 1216 East California Boulevard, Pasadena, CA 91125, USA

³Department of Astronomy, The University of Texas at Austin, 2515 Speedway Boulevard, Stop C1400, Austin, TX 78712, USA

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

⁵MPE, Giessenbachstrasse 1, Garching D-85748, Germany

⁶Cluster of Excellence, Boltzmann Strasse 2, D-85748, Germany

Received 2019 May 8; revised 2019 June 29; accepted 2019 July 8; published 2019 August 22

## Abstract

We examine the role of environment on the in situ star formation (SF) hosted by the progenitors of the most massive galaxies in the present-day universe, the brightest cluster galaxies (BCGs), from $z \sim 3$ to present in the COSMOS field. Progenitors are selected from the COSMOS field using a stellar mass cut motivated by the evolving cumulative comoving number density of progenitors within the Illustris simulation, as well as the Millennium-II simulation and a constant comoving number density method for comparison. We characterize each progenitor using far-ultraviolet–far-infrared observations taken from the COSMOS field and fitting stellar, dust, and active galactic nucleus components to their spectral energy distributions. Additionally, we compare the SF rates of our progenitor sample to the local density maps of the COSMOS field to identify the effects of environment. We find that BCG progenitors evolve in three stages, starting with an in situ SF-dominated phase ($z > 2.25$). This is followed by a phase until $z \sim 1.25$ where mass growth is driven by in situ SF and stellar mass deposited by mergers (both gas rich and poor) on the same order of magnitude independent of local environment. Finally, at low redshift dry mergers are the dominant stellar mass generation process. We also identify this final transition period as the time when progenitors quench, exhibiting quiescent NUVrJ colors.

**Key words:** galaxies: abundances – galaxies: clusters: general – galaxies: elliptical and lenticular, cD – galaxies: star formation – galaxies: structure

## 1. Introduction

The evolutionary history of today's most massive galaxies is an important component in understanding the evolution of large-scale structure in the universe and galaxy populations in general, placing vital limits on the speed that baryonic matter can assemble into gravitationally bound structures. These are massive elliptical galaxies, called brightest cluster galaxies (BCGs), and are normally located at the center of their cluster in position/velocity space (Lauer et al. 2014). They are composed of an old stellar population and a hot gaseous component in an extended envelope larger than predicted by the distribution of neighboring non-BCG ellipticals (Oemler 1976). While most are quiescent, those that inhabit relaxed, massive halos may exhibit atypically large ($>1 M_{\odot} \text{ yr}^{-1}$) star formation rates (SFRs). This is due to the surrounding intra-cluster medium (ICM) cooling and precipitating into the BCGs, forming cool-core galaxy clusters (e.g., Fabian 1994; O'Dea et al. 2008, 2010; Rawle et al. 2012; Fogarty et al. 2019), or due to the accretion of gas-rich satellites (e.g., Tremblay et al. 2014; Webb et al. 2015a; Cooke et al. 2016). BCG stellar mass and SFRs may also be affected by the mass of their cluster halos at low redshift (e.g., Cooke et al. 2018; Gozaliasl et al. 2018). Their atypically high mass, quiescence, and extended light profiles collectively point to a formation mechanism unique with regard to their satellites, in which most of their stars formed at $z > 2$, and conglomerated into a single BCG through many major and minor mergers (De Lucia & Blaizot 2007). With this complex model in mind, a full identification of BCG progenitors requires the characterization

of all galaxies at high redshift that are likely to merge and contribute their stellar mass onto proto-BCG “seeds.”

BCGs inhabit the high-mass end of the galaxy stellar mass function (SMF); at $z \sim 0.5$, BCGs in massive, relaxed clusters are typically $\sim 10^{11.5} M_{\odot}$ or greater (e.g., Lidman et al. 2012; Lin et al. 2013; Chiu et al. 2016; Cooke et al. 2016). At this stellar mass, the population is dominated by passive (elliptical) galaxies, both BCGs and non-BCG ellipticals, that are three orders of magnitude more common than star-forming (spiral) galaxies (Davidzon et al. 2017), particularly in dense environments (e.g., Scoville et al. 2013; Darvish et al. 2016). By starting with this low-redshift population, and then identifying the progenitor populations responsible, we examine the role of internal versus external processes at a time when proto-clusters have not finished merging into the virialized clusters observed at present day. In situ star formation (SF) refers to any SF processes within a galaxy, such as those triggered through secular (disk instability; Livermore et al. 2012, 2015) or external stimuli (interactions or gas-rich mergers; Sanders & Mirabel 1996; Kartaltepe et al. 2012). Ex situ stellar mass growth is the accretion of stellar mass formed outside of the galaxy, such as the old stellar populations composing the secondary galaxy in a galaxy merger (e.g., Pillepich et al. 2015).

The current hierarchical model of galaxy formation relies on a combination of stellar mass growth through internal processes and mergers that build galaxies from several individual galaxies. Mergers directly inject stellar mass into the primary galaxy concurrent with in situ SF triggered if either galaxy is gas rich. The relative roles of internal processes and mergers have been simulated at length in large-scale cosmological models (e.g., the Millennium-II Simulation (MS-II) and the
Illustris project; Springel et al. 2005; Vogelsberger et al. 2014) and zoom-in models (Narayanan et al. 2010, 2015; Ragone-Figueroa et al. 2018). While there is still tension with observations of properties such as morphology (e.g., Genel et al. 2014), many predictions from simulations, such as merger rates, remain within observational errors (Lotz et al. 2011; Rodriguez-Gomez et al. 2017; Duncan et al. 2019). The total stellar mass growth of a BCG, or any massive galaxy, can be summarized by two drivers with multiple contributors.

1. In situ star formation

(a) Self-sourced—a galaxy's cold gas supply condenses via local or gravitational interaction/flyby.

(b) ICM precipitation—cooling gas from the cluster ICM precipitates down the potential well and forms stars.

(c) Wet mergers—stars form due to the shocks and greater gas supply provided during a gas-rich merger.

2. External delivery of stellar mass

(a) Wet and dry mergers—immediate contribution of an existing, older stellar population.

(b) Local interaction—capture of stellar mass from an interaction with an unbound secondary galaxy.

The relative role of each of these contributors is a rich field of research, with many works specializing in one or more facet of this picture (e.g., O'Dea et al. 2008; Naab et al. 2009; Loubser & Sánchez-Blázquez 2012; Webb et al. 2015a, 2017; Contini et al. 2018). Another factor is the role of local environment on the stellar mass assembly processes described above. A correlation between local environment and SF is observed at low redshift (e.g., Oemler 1976; Dressler 1980; Kauffmann et al. 2004; Darvish et al. 2018), but is still a subject of debate at high redshift (e.g., Elbaz et al. 2007; Scoville et al. 2013; Darvish et al. 2016). Active galactic nucleus (AGN) activity as a function of environment at high redshift is likewise under investigation (e.g., Ehlert et al. 2013; Martini et al. 2013; Umehata et al. 2015; Alberts et al. 2016; Yang et al. 2018; Macuga et al. 2019).

In this paper, we seek to understand the roles of in situ SF and external delivery toward the total growth of BCG progenitors as a function of environment and redshift. However, identifying progenitors to specific populations in the present-day universe remains one of the most difficult, but critical, elements of galaxy evolution studies. Due to the long timescales involved, we require a method that connects galaxies across redshift and is independent of SFR characteristics such as color. This is often implemented as a stellar mass selection as a function of redshift and descendent stellar mass.

A relatively straightforward method for identifying the progenitors of massive galaxies is to assume that the most massive galaxies in clusters and groups remain the most massive throughout their history (van Dokkum et al. 2010). One implementation of this concept is the constant comoving number density method (Leja et al. 2013; Patel et al. 2013; van Dokkum et al. 2013), which assumes that the number of galaxies within a given comoving volume of the universe is constant. As the galaxies age within the volume, the population within each stellar mass bin will grow in stellar mass with time, and inhabit the new stellar mass bin at an equal density as their old bin. By selecting galaxies at high redshift that inhabit the universe at an equal cumulative number density, one selects the direct progenitors of the low-redshift population. This method can be corrected to include

merger partners by increasing the cumulative comoving number density cutoff (i.e., lowering the stellar mass cutoff), and this evolving cutoff has been investigated using semi-empirical approaches or cosmological simulations (e.g., Marchesini et al. 2014; Torrey et al. 2015).

We use the evolving comoving cumulative number density selection method to identify the progenitors of BCGs out to redshift $z \sim 3$. We select a descendant stellar mass cutoff of $10^{11.5} M_{\odot}$ at $z < 0.35$ to ensure we track progenitors to BCGs in the low-redshift universe. To characterize each progenitor, we fit the spectral energy distributions (SEDs) of each target from the far-ultraviolet (FUV) to far-infrared (FIR) to ensure sensitivity to low and obscured SFRs. Additionally, we build upon previous studies (e.g., Hill et al. 2017; Torrey et al. 2017) by correlating our progenitor characteristics to the local galaxy environment using redshift-sensitive adaptive kernel density maps of the COSMOS field (Darvish et al. 2015, 2017).

In Section 2 we discuss the data we used to construct the SEDs of our sample. In Section 3 we describe how we implemented an evolving comoving number density selection function to acquire our BCG progenitor sample. The median SED fits of our sample are discussed in Section 4. We consider the physical implications of our results in Section 5 and we summarize our findings in Section 6. We use the $\Lambda$CDM standard cosmological parameters of $H_0 = 70 \text{ Mpc}^{-1} \text{ km s}^{-1}$, $\Omega_M = 0.3$, and $\Omega_{\text{vac}} = 0.7$ and a Chabrier (2003) initial mass function.

## 2. COSMOS Multiwavelength Data

We first require a large contiguous area and volume to examine progenitors across the full range of cosmological environments. Additionally, we need well-sampled SEDs from the FUV to the FIR to constrain the stellar mass and SFRs of the massive galaxies within these environments. These requirements are uniquely satisfied by the COSMOS field (Scoville et al. 2007). COSMOS is a 2 sq. deg. field centered at R.A.(J2000) = 10:00:28.600, decl.(J2000) = +02:12:21.00 observed from the X-ray to the radio. We use the public COSMOS photometric data catalog, COSMOS2015 (Laigle et al. 2016), to select progenitors and construct their SEDs from the FUV to FIR. Below, we concisely summarize each of the COSMOS data sets used in this study; for further details of COSMOS2015 see Laigle et al. (2016).

### 2.1. Photometric and Spectroscopic Redshifts

The quality of our SED fits and environment estimates depends on accurate distance measurements. To that end, we start with COSMOS photometric redshift measurements (Ilbert et al. 2008) in the COSMOS2015 catalog (Laigle et al. 2016). We then crossmatch the COSMOS2015 progenitors' coordinates with the spectroscopic detections in the COSMOS spectroscopic catalog to the closest match within 1" of the original COSMOS2015 location. Our BCG progenitor sample includes spectroscopic redshifts available from the 3D-Hubble Space Telescope (HST) Survey (Brammer et al. 2012; Momcheva et al. 2016), Keck DEIMOS 10 K Spectroscopic Survey (Hasinger et al. 2018), FMOS-COSMOS Survey (Kashino et al. 2013; Kartaltepe et al. 2015b; Silverman et al. 2015), the Gemini GMOS-S spectra of Balogh et al. (2011), COSMOS AGN Spectroscopic Survey (Trump et al. 2007, 2009), the Keck LRIS spectra of Casey et al. (2017),
MOIRCS Deep Survey (Yoshikawa et al. 2010), the Keck MOSFIRE spectra of Trakhtenbrot et al. (2016), the Keck DEIMOS spectra of Capak et al. (2011) and M. Salvato et al. (2019, in preparation), the SINFONI spectra of Perna et al. (2015), the zCOSMOS Survey (Lilly et al. 2007), the VIMOS Ultra-Deep Survey (Le Fèvre et al. 2015), and the HST/WFC3 grim spectra of Krogager et al. (2014).

## 2.2. GALEX

We use FUV and near-ultraviolet (NUV) band point-spread function (PSF)-fit photometric magnitudes from the *Galaxy Evolution Explorer (GALEX*; Martin et al. 2005) to constrain the degree of unobscured SF. *GALEX* observations in the COSMOS2015 catalog were originally reduced by Zamojski et al. (2007). We correct the *GALEX* FUV-Spitzer IRAC4 observations for Milky Way foreground extinction using a galactic reddening of $R_v = 3.1$ (Morrissey et al. 2007) and $E(B-V)$ values from the dust maps of Schlegel et al. (1998).

## 2.3. Canada-France-Hawaii Telescope

Canada-France-Hawaii Telescope (CFHT)/MegaPrime (Aune et al. 2003; Boulade et al. 2003) $u^*$ observations of the COSMOS field were taken with a consistent depth of $m_{u^*} \sim 26.4$ (Capak et al. 2007).

## 2.4. Subaru

To constrain our stellar mass estimates, we need sensitive optical continuum measurements across as many bands as possible. Therefore, we include Subaru/Suprime-Cam optical observations using five broadband filters (B, V, R, i+, z++) and 11 medium-band filters (IA427, IA464, IA484, IA505, IA527, IA574, IA624, IA679, IA738, IA767, and IA827), observed to a $3\sigma$ depth of $m_{AB} \sim 25.2$ or deeper. Out of these filters, IA464 suffers the worst resolution with a PSF FWHM of $1.89$ (Taniguchi et al. 2007, 2015).

## 2.5. Vista

We include Vista/VIRCAM (Sutherland et al. 2015) J-, H-, and K-band observations from the UltraVISTA-DR2 survey (McCracken et al. 2012). These have limiting magnitudes $m_{AB}$ of 24.7, 24.3, and 24.0 respectively.

## 2.6. Spitzer Observations

Further IR observations of the COSMOS field are available from the SPLASH (Steinhardt et al. 2014) and S-COSMOS surveys (Sanders et al. 2007) using the *Spitzer Space Telescope* (*Spitzer*; Werner et al. 2004). Data from *Spitzer's* Infrared Array Camera (IRAC)'s 3.6, 4.5, 5.7, and 8.0 µm channels (for more information, see Fazio et al. 2004) have PSF widths $1.6''$, $1.6''$, $1.8''$, and $1.9''$, respectively, and were observed down to a $3\sigma$ depth of $m_{AB}$ of 25.5, 25.5, 23.0, and 22.9, respectively. We also use MIPS (Rieke et al. 2004) 24 µm measurements, which act as an important constraint to the total IR fit. MIPS 24 µm was observed to a $5\sigma$ depth of $71$ µJy.

## 2.7. Herschel Observations

To accurately estimate the obscured AGN and SF activity in our sample, we use FIR observations to constrain the FIR peak.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Example of the constant number density selection method, using the assumption of a constant number of galaxies in the universe (total box). A population of low-redshift, quiescent galaxies (shaded red spirals) on the right correspond to their high-redshift progenitors on the left (shaded yellow spirals) by the density they inhabit the universe, irrespective of parameters such as color or star formation activity.</figcaption></figure>

COSMOS FIR observations include the *Herschel Space Observatory* (Herschel; Pilbratt et al. 2010) Photoconductor Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) 100 and 160 µm bands and Spectral and Photometric Imaging Receiver (SPIRE) 250, 350, and 500 µm bands. The PACS Evolutionary Probe (Lutz et al. 2011) data used in this study were observed down to a $3\sigma$ depth of 5 and 10.2 mJy for 100 and 160 µm bands, respectively. Herschel Multi-tiered Extra-galactic Survey (Oliver et al. 2012) SPIRE 250 µm, 350 µm, and 500 µm observations reach a $3\sigma$ depth of 8.1, 10.7 mJy, 15.4 mJy, respectively.

Blending is possible at the longest wavelengths due to the large beam size ($18''$1 at 250 µm, $24''$9 at 350 µm, $36''$6 at 500 µm), therefore de-blending was performed band by band using the next most-resolved observation starting with *Spitzer* IRAC observations as priors (Lee et al. 2010). Blending is the largest issue at 500 µm, so to ensure that the 500 µm detections correspond to our progenitor sample, we use the “clean index” from Elbaz et al. (2011). In brief, a progenitor’s 500 µm data are rejected if there is more than one bright neighbor within $1.1 \times$ the FWHM at 24 µm or there are bright neighbors within $1.1 \times$ the FWHM of the respective *Herschel* FIR band, with a bright neighbor defined as $S_{\text{Neighbor}}/S_{\text{Target}} > 0.5$.

# 3. Methods

## 3.1. Cumulative Number Density Selection Function

Due to the long timescales involved in an individual galaxy's evolution, large samples over many epochs are needed to probe the growth of stellar mass in the universe. Selecting a population of progenitors at high redshift is difficult because the entire point of the exercise is that an important characteristic of the progenitor population remains unknown. The constant comoving cumulative number density method introduced by van Dokkum et al. (2010) was a solution to this problem. This method assumes a constant number of galaxies in the past and present universe. To identify the progenitors of massive galaxies in the present day, one would measure the density they inhabit in a volume element of the present-day universe and identify a corresponding population at high redshift that occupied the universe at the same volume density (illustrated in Figure 1).

This method has been successfully used to identify some of the progenitors for today's population of massive galaxies (van Dokkum et al. 2010; Patel et al. 2013; Morishita et al. 2015), but it is biased toward higher-mass progenitors. This bias arises because the galaxies that merge onto the selected progenitor
galaxy are not included. Secondary galaxies are often lower mass (of higher number density) and therefore their exclusion results in a bias toward lower number density and higher-mass progenitors (Mundy et al. 2015). Progenitors that evolved via a merger-dominant evolutionary path would not be selected, as these low-mass progenitors would inhabit a lower stellar mass (and higher number density) parameter space than the rest of the progenitor population. The disparity between predicted and “actual” progenitor stellar masses could reach a factor of two at high redshift (Leja et al. 2013; Mundy et al. 2015). To better correct for the effect of mergers, an evolving number density method was developed by tracking the median progenitor mass in cosmological simulations (Mundy et al. 2015; Torrey et al. 2015; Jaacks et al. 2016).

An evolving number density selection is derived by identifying a population in a low-redshift slice of the universe and tracing back in redshift the number density that corresponds to the median stellar mass of their most massive progenitors. This method has been used with multiple simulations, and our analysis includes the Illustris simulation number density evolution tracks (Torrey et al. 2015, 2017), with comparison tracks from MS-II (Boylan-Kolchin et al. 2009) and a constant number density cut as discussed in Section 5.

The Illustris simulation is a cosmological hydrodynamic simulation of volume $1.2079 \times 10^6$ Mpc$^3$ that includes dark and baryonic physics (Vogelsberger et al. 2014). Baryonic processes include radiative cooling, SF and supernova feedback, and AGN feedback. Baryonic mass is resolved with particles of mass $1.3 \times 10^6 M_\odot$. The inclusion of baryonic physics enables Illustris to more accurately model the evolution of SF at late times in the universe as AGNs and outflows heat and push out the gas supply within each galaxy.

To estimate stellar masses of progenitors observed in COSMOS, we use Davidzon et al. (2017)'s SMF to identify the cumulative number density of galaxies with stellar masses above $10^{1.5} M_\odot$ below $z \sim 0.35$, which corresponds to the regime of low-redshift BCGs (e.g., Lidman et al. 2012; Liu et al. 2012; Fraser-McKelvie et al. 2014; Cooke et al. 2016). This high stellar mass selection specifically examines high-mass BCGs commonly hosted in clusters with total mass $>10^{14} M_\odot$ (e.g., Cavagnolo et al. 2009; Postman et al. 2012; Lin et al. 2013) to better understand how the most massive galaxies acquired their stellar populations. We choose to use the SMF of the total population and not only the passive population, as the latter outnumbers the active population by a factor of ten at this mass and redshift and we do not wish to bias our selection toward the most passive progenitors. The results shown in Section 5 are consistent within errors if we exclude or include active descendants at low redshift. Once the density of the total descendant population is identified, we determine the cumulative comoving number density of BCG progenitors in redshift slices out to $z \sim 3$ using an Illustris-derived number density evolution relation (Equation (A4) of Torrey et al. 2017),

$$ \langle N(z) \rangle' = N_0 + \Delta z (A'_0 + A'_1 N_0 + A'_2 N_0^2) + \Delta z^2 (B'_0 + B'_1 N_0 + B'_2 N_0^2), $$

where $N_0$ represents the cumulative comoving number density of BCGs below $z \sim 0.35$, and the constants $A'_{0,1,2}$ and $B'_{0,1,2}$ are best-fit parameters to the backward median number density evolution.

<figure><img src="image_4.png" /><figcaption>Figure 2. Stellar mass selection function of our sample via constant comoving number density (orange), and the evolving comoving number density predicted via Illustris (red) and MS-II (teal). We also identify the stellar mass cutoff for the massive galaxy population that deposits 50% of the merger delivered stellar mass at that redshift (blue dotted) and plot the galaxies from COSMOS2015 that satisfy this lower mass cut (gray circles). This lower cut enables us to identify galaxies with initial COSMOS estimates below our Illustris cutoff that may be refit to higher masses. The mass incompleteness limit is hatched green (Laigle et al. 2016), and is defined as the stellar mass above which 90% of galaxies are detected given the $K_s$ limiting magnitude of the COSMOS field observations (Pozzetti et al. 2010; Davidzon et al. 2013; Ilbert et al. 2013; Moustakas et al. 2013).</figcaption></figure>

## 3.2. Application to Observations

To construct our final sample, we take the cumulative comoving number density predictions estimated above and calculate the stellar mass corresponding to the progenitor cumulative number density in that redshift slice (red line in Figure 2) given the observed COSMOS SMF (Davidzon et al. 2017). To maintain sample completeness, we fit the SEDs of an expanded sample down to a lower stellar mass limit than required by our Illustris-derived selection function. By fitting COSMOS galaxies within a factor of 3 of the Illustris selection function, we hope to recover any galaxies with underestimated masses in their initial fit in the COSMOS2015 catalog. After fitting our initial sample of prospective progenitors (dotted blue line in Figure 2), we select galaxies that have a fit stellar mass (described in Section 4) above the Illustris evolving comoving number density method mass cut for that given redshift bin (red line in Figure 2). Our final sample includes 1444 BCG progenitors from $0.5 < z < 3.0$.

## 3.3. Comparison Sample Selection Methods

We also compare our results based on the Illustris selection function to other commonly used selection methods. MS-II (Boylan-Kolchin et al. 2009) is a dark-matter-only simulation and does not include baryonic physics. Number density is thus derived from dark matter halos characteristic of our low-redshift massive elliptical sample. For MS-II halos, we identify the stellar mass of a given halo mass’s baryonic counterpart using Guo et al.’s (2013) Millennium WMAP 7 (MR7) semi-analytic galaxy models scaled to the MS-II simulation. We start similarly to our Illustris method by using Davidzon et al. (2017)’s SMF to identify the cumulative number density of galaxies with stellar masses above $10^{11.5} M_\odot$ at $z \sim 0.35$. Then we identify the dark matter halo mass corresponding to this number density cutoff. We identify the most massive halo
progenitor in each higher-redshift bin, and find the number
density corresponding to the median of the most massive halo.
This number density is then used with the corresponding
Davidzon et al. (2017) SMF at that redshift to identify the galaxy
progenitor population from Guo et al.'s (2013) Millennium-II
mock galaxy catalog. We make the assumption that a single
BCG progenitor lies in each halo. The MS-II derived evolution
track is shown as the teal solid line in Figure 2. We discuss how
the MS-II-selected sample results compare with our primary
Illustris-selected results in Section 5.6.

Finally, we wish to compare our results to those derived
from the original constant number density method. Our
constant number density comparison sample is chosen by
identifying the cumulative comoving number density of
ellipticals with stellar masses $10^{11.5} M_{\odot}$ at $z \sim 0.35$ in
Davidzon et al. (2017)'s SMF, and identifying the stellar mass
of the cumulative population above that density in all redshift
bins out to $z \sim 3$. This stellar mass and redshift combination
corresponds to the massive BCG population in the low- to
intermediate-redshift universe (e.g., Lidman et al. 2012; Liu
et al. 2012; Fraser-Mcelvie et al. 2014; Cooke et al. 2016).
Our constant number density evolution track is shown as the
orange solid line in Figure 2. A comparison between the
constant number density selection method and the evolving
number density method used for the main results of this work is
included in Section 5.7.

3.4. Quantifying Local Environment

With progenitor stellar mass selection criteria calculated, our
next step is to select galaxies likely to merge onto the BCG
"seeds" over time given their local environment. The sample
selection discussed in Section 3.1 identified progenitors,
irrespective of environment, that are predicted to relocate to
the overdense regions where present day BCGs ultimately
reside. Therefore we use density maps of the COSMOS field to
identify progenitors already in overdense regions (very likely to
merge into a BCG within a Hubble time) and underdense
regions (less likely to merge into a BCG).

We quantify local environment using density field maps of
the COSMOS field produced using the “weighted” adaptive
kernel smoothing approach of Darvish et al. (2015). The
density field is evaluated in a series of overlapping redshift
slices. Darvish et al. assign galaxy weights that are proportional
to the likelihood of a galaxy belonging to a redshift. They are
estimated by measuring what fraction of a galaxy’s redshift
distribution is within a particular redshift slice. Given these
weights, the density field is then adaptively smoothed with a
Gaussian kernel whose “global” width is 0.5 Mpc. The
adaptive widths are evaluated around this global width,
depending on how sparse or dense the local neighborhood of
each galaxy is (see Darvish et al. 2015 for details).

3.5. Identification of AGNs

We also quantify the contribution of the AGN emission to
the IR regime of progenitor SEDs. This step is important, as a
dominant AGN could cause the SED fit solution to over-
estimate stellar mass, and an IR-bright AGN could result in a fit
with an overestimated SFR. We crossmatch our sample with
the Chandra X-ray observations of the COSMOS field (Elvis
et al. 2009; Civano et al. 2016; Marchesi et al. 2016; Lanzuisi
et al. 2017) to identify X-ray-bright AGNs. We also use the

*Spitzer* IRAC photometry available for the COSMOS field to identify obscured AGNs through their mid-infrared (MIR) colors using the Donley et al. (2012) IRAC color criteria. Radio AGNs are identified by crossmatching our progenitor sample with the radio AGNs identified in the VLA-COSMOS 3 GHz Large Project (Delvecchio et al. 2017; Smolčić et al. 2017, 2008). Our full progenitor and AGN sample selection totals and statistics are listed in Table 1.

4. SED Fitting

4.1. Software and Methodology

To determine the SFR and total stellar mass of each
progenitor, we require a way to reliably fit thousands of galaxy
SEDs from rest-frame FUV to FIR. To do so, we choose the
SED fitting software MAGPHYS (da Cunha et al. 2008) with
50,000 stellar population models generated from Bruzual &
Charlot (2003) and 50,000 infrared dust models from da Cunha
et al. (2008).

MAGPHYS fits the SED of each galaxy in two stages. First, it estimates the best-fit stellar model between 0.0912 µm ~ λ ~ 10 µm from the Bruzual & Charlot (2003) stellar population models while taking into account dust obscuration. Once the optical component is fit, MAGPHYS fits the infrared from 2.5 to 1000 µm while taking the total luminosity of obscured starlight found during the optical fitting as a prior. FIR characteristics are calculated between 8 and 1000 µm. This ensures energy balance between the energy absorbed and emitted by the dust component.

At time of writing, the public version of MAGPHYS does
not include an AGN component. This is an important
consideration, as additional luminosity across the optical and
near-IR from an AGN accretion disk and the dusty torus
surrounding it may skew stellar mass estimates to higher
values. Additionally, any UV or IR emission may also cause
SFRs to be overestimated. Therefore to identify potential
overestimations due to a lack of an AGN fitting component, we
investigated several goodness-of-fit metrics to identify which
progenitors require an AGN component. We find that a
fractional difference cutoff of 0.4 between the MAGPHYS fit
8.0 µm (Spitzer IRAC4) flux and the observed value worked
most consistently to identify sources that require an AGN
component in the MIR (for further details, see K. D. Tyler et al.
2019, in preparation). This parameter measures how poorly the
stellar and dust models fit the slope of the MIR, which could
indicate the requirement of a third (AGN) component. We find
that ~1% of our progenitor sample fits have MAGPHYS
fractional 8.0 µm residuals >0.4, and we refit this subsample
using the package SED3FIT. By choosing to select and fit
AGN components to only the galaxies that have poor IR fits
without an AGN, we prevent the overfitting of our sample and
the possible underestimation of SFRs that an overestimation of
AGN activity and frequency would induce.

SED3FIT (Berta et al. 2013) works similarly to MAGPHYS
by fitting a series of optical models first, and then using the
optical results as a prior during the fitting of the IR component.
SED3FIT, however, includes AGN torus and accretion disk
emission (Feltre et al. 2012) across both components while
maintaining energy balance. For the IR-bright AGNs identified
through their 8 µm residuals, we use the SFR and stellar mass
estimates from SED3FIT instead of MAGPHYS in the
following results.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Left: original 10 Feltre et al. (2012) AGN models included in SED3FIT with viewing angles of 0° or 90°. Right: our selection of 240 AGN torus+disk models from Feltre et al. equally sampling viewing angle.</figcaption></figure>

<table><caption>**Table 1**

Illustris-selected BCG Progenitor Sample Distribution</caption><thead><tr><th>z</th><th>Total Progenitors</th><th>Total Photo-z</th><th>Total Spec-z</th><th>Total AGNs<sup>a</sup></th><th>X-ray AGNs</th><th>MIR AGNs</th><th>Radio AGNs</th><th>SED3FIT AGNs<sup>b</sup></th></tr></thead><tbody><tr><td>0.5–0.8</td><td>26</td><td>2</td><td>24</td><td>2 (7.7%)</td><td>1</td><td>1</td><td>2</td><td>0</td></tr><tr><td>0.8–1.1</td><td>180</td><td>74</td><td>106</td><td>14 (7.8%)</td><td>1 2</td><td>6</td><td>0</td><td>5</td></tr><tr><td>1.1–1.5</td><td>216</td><td>156</td><td>60</td><td>8 (3.7%)</td><td>4</td><td>4</td><td>1</td><td>3</td></tr><tr><td>1.5–2.0</td><td>401</td><td>371</td><td>30</td><td>29 (7.2%)</td><td>22</td><td>15</td><td>1</td><td>12</td></tr><tr><td>2.0–2.5</td><td>347</td><td>308</td><td>41</td><td>32 (9.2%)</td><td>15</td><td>24</td><td>1</td><td>9</td></tr><tr><td>2.5–3.0</td><td>274</td><td>254</td><td>20</td><td>31 (11.3%)</td><td>10</td><td>26</td><td>0</td><td>8</td></tr></tbody></table>

Notes. Galaxies selected from the COSMOS2015 catalog using the stellar mass cut predicted via the median BCG progenitor population in the Illustris simulation.

<sup>a</sup> Total AGN fractions represent the fraction of the AGN population identified through any method out of the total progenitor sample in that redshift bin. An AGN may satisfy multiple criteria, and are counted in any classification count that they satisfy.

<sup>b</sup> Defined as AGNs by a fractional residual >0.4 at 8.0 µm (Spitzer IRAC4) after fitting with a MAGPHYS two-component model (see Section 4).

## 4.2. AGN Template Library

The original AGN template library included in the public SED3FIT distribution included 10 AGN templates from the Feltre et al. (2012) AGN template library. These templates were selected at extreme viewing angles, effectively probing archetypical Type-1 and Type-2 AGNs. These models assume either negligible or strong extinction with respect to line of sight. While satisfactory for first-order fits, the model of the obscuring torus surrounding the supermassive black hole in an AGN has changed in recent years to a clumpy torus with regions of high and low obscuration at high and low viewing angles (e.g., Krolik & Begelman 1988; Shi et al. 2006; Markowitz et al. 2014). Therefore, we expand our AGN template library to include 240 models (K. D. Tyler et al. 2019, in preparation) equally probing between 0° and 90°. Shown in Figure 3, this new library contains models that fit AGNs at intermediate obscurations and viewing angles.

In Figure 4, we plot the median SED of the progenitors found in three environments in six redshift bins with $1\sigma$ confidence intervals (shaded regions). We classify each environment as dense ($\frac{\delta}{\delta_{median}} > 2$), intergroup ($1 < \frac{\delta}{\delta_{median}} < 2$), or field ($\frac{\delta}{\delta_{median}} < 1$).

In the top panel, we find that progenitors in dense regions have SEDs similar to those in other environments at high redshift ($z > 1.1$), but have ~$\frac{1}{5}$ the MIR emission (10–20 µm) at $z < 1.1$ in comparison to progenitors at lower density at low redshift. This indicates a divergence in evolution due to environment at low redshift, where massive progenitors in high-density environments undergo a more efficient removal of their warm gas and dust supplies.

## 5. Results and Discussion

To characterize our total sample of 1444 BCG progenitors, we use the stellar parameters fit by MAGPHYS for 1407 progenitors and the estimates from SED3FIT for 37 progenitors due to their significant 8.0 µm residuals in their original MAGPHYS fits (hereafter referred to as SED3FIT AGNs). The SED3FIT AGN fraction is <1% in any redshift bin, with an overall AGN fraction found through X-ray, MIR, and radio methods of ~8% (Table 1). For comparison, previous Sunyaev–Zeldovich–selected studies have found X-ray cavities indicating past AGN activity in ~7% of low-redshift BCGs (Hlavacek-Larrondo et al. 2015).
<figure><img src="image_6.png" /></figure>

**Figure 4.** Median spectral energy distributions (solid lines) with $1\sigma$ confidence intervals (shaded regions) of massive progenitors between $0.5 < z < 3$ for dense regions (top), intergroup regions (middle), and field regions (bottom). Dense environments are defined as $\frac{\delta}{\delta_{\text{median}}} > 2$, where $\delta$ is the density at that location and $\delta_{\text{median}}$ represents the median density at that redshift. Intergroup is defined as $1 < \frac{\delta}{\delta_{\text{median}}} < 2$. Field environments are densities $\frac{\delta}{\delta_{\text{median}}} < 1$. Color coding corresponds to redshift bins and is identical in all three subplots. We find that across all environments, FIR luminosity decreases with decreasing redshift. However, MIR luminosity decreases with decreasing redshift at a much higher rate in dense environments in comparison to intergroup or field starting in our $0.8 < z < 1.1$ redshift bin. This may be due to either a lack of warm dust radiating in the MIR or a lack of SFR illuminating the dust supply in progenitors residing in dense environments. The latter is more likely as the SFR is also decreasing with decreasing redshift.

## 5.1. SF and Stellar Mass Estimates

Our sample spans an order of magnitude in stellar mass from $10^{11}$ to $10^{12} M_{\odot}$ at $z \sim 3$ (Figure 5), motivated by the stellar mass selection function shown in Figure 2. Our sample below $z \sim 0.8$ is subject to small-number statistics due to the low volume available in COSMOS at this distance ($\sim 4.4 \times 10^6 \text{ Mpc}^3$).

Shown in Figure 5, both SFR (averaged over $10^8$ yr) and specific SFR (sSFR) span three orders of magnitude from $0.5 < z < 3$. Progenitors do not span this parameter space equally, with a concentration of highly star-forming progenitors with sSFR $\sim 10^{-9.25} \text{ yr}^{-1}$ at $z > 2$, that transition to a lower sSFR of $10^{-10.75} \text{ yr}^{-1}$ at $z < 2$. This is partly driven by our stellar mass selection function, which is limiting the stellar mass parameter space to higher stellar masses at low redshift. Our observed relation at the bottom of Figure 5 is systematically an order of magnitude below the sSFR–z correlation found for less-massive ($10^{10.5}$) galaxies found in the COSMOS field (Davidzon et al. 2018), consistent with galactic “downsizing” where less-massive galaxies continue to form stars while massive galaxies are quenching (Cowie et al. 1996; Pérez-González et al. 2008). This will yield descendants dominated by the old stellar populations, just as observed in BCGs in the nearby universe (Loubser et al. 2009).

## 5.2. Evolution of Stellar Mass Growth

Our goal is to investigate what role in situ SF plays in the total growth of stellar mass in BCG progenitors. To do so, we need to compare the growth solely due to SF to a measurement of total stellar mass growth. We define our total growth rate to be the change in median stellar mass between each redshift bin

in units of $M_{\odot} \text{yr}^{-1}$. We then perform a polynomial fit to these growth rates at each bin boundary and plot the result as the solid green line in Figure 6.

We determine the growth rate of each progenitor due to in situ SF by taking the SFRs from our SED fits and applying a reduction of 50%. This scale factor used by Brinchmann et al. (2004) and van Dokkum et al. (2010) corrects for stellar mass loss over a Gyr timescale. This correction to SFR produces a stellar mass growth measurement that considers the net growth influenced by SF and star destruction processes. Stellar mass loss occurs as the population of stars born during an SF episode lose their massive members over time as they age off the main sequence. The final in situ growth rate represents the total stellar mass per year generated in a galaxy that will remain on multi-Gyr timescales. After we apply this correction, we calculate the in situ growth rate for progenitors in dense ($\frac{\delta}{\delta_{\text{median}}} > 2$) and field ($\frac{\delta}{\delta_{\text{median}}} < 1$) environments and plot the mean stellar mass growth rates due to SF for each sample as red and blue lines, respectively, in Figure 6. $\delta$ is defined as the density at that location while $\delta_{\text{median}}$ is the median density of the redshift slice.

Finally, we wish to understand how the stellar mass directly delivered via gas-rich and gas-poor mergers compares to the above values. Therefore we integrate the stellar mass deposition rate ($M_{\odot} \text{yr}^{-1}$) from mergers of secondary galaxies with mass ratios of 1:1 to $1:10^{-5}$ of the progenitor’s stellar mass in that redshift bin using the merger rates provided by the Illustris simulation (Rodríguez-Gómez et al. 2017). This is plotted as the dashed gray line in Figure 6.

Following Hill et al. (2017), we plot the above values as $dM/dt (M_{\odot} \text{yr}^{-1})$ versus redshift from SF, mergers, and all sources in Figure 6. From this comparison, we identify three epochs of BCG progenitor growth.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Top: stellar mass vs. redshift of massive BCG progenitors between 0.5 < z < 3 fit using MAGPHYS and SED3FIT. This plot includes all sources with stellar mass estimates above the stellar mass selection cut described in Section 3. Middle: SFR vs. redshift of massive BCG progenitors between 0.5 < z < 3. Bottom: specific SFR vs. redshift of massive BCG progenitors between 0.5 < z < 3.</figcaption></figure>

### 5.2.1. The SF-Dominated Epoch ($z > 2.25$)

At $z > 2.25$, the in situ SF estimated in progenitors in all environments is consistent with the total growth rate of the median stellar mass of our sample. This consistency indicates that progenitor growth is dominated by active SF with contributions of direct stellar mass injection and gas via gas-rich and gas-poor mergers. The total stellar mass growth rate shown in Figure 6 is motivated by the stellar mass selection function that was used to originally select our sample; however, we may still conclude that the active SFRs via secular or wet merging processes are high enough to be the primary method of individual progenitor growth. Dry mergers are unnecessary, but may still contribute to the highest-mass cases.

### 5.2.2. The Transitionary Epoch ($1.25 < z < 2.25$)

Between $1.25 < z < 2.25$, progenitors are forming their in situ stellar mass at a rate ~60%–75% of the observed

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Observed stellar mass growth rate for massive galaxy progenitors in dense regions (red) and field regions (blue) compared with the median growth rate of our sample's median stellar masses (green). The galaxies selected in dense environments are those most likely to form the BCGs we see today, while those in the field sample are candidates to merge into a BCG during cluster mergers at later times. We overplot the stellar mass deposition rate of mergers in the Illustris simulation (gray dashed-dotted) with a stellar mass capture efficiency of 100%. The hashed blue region in our lowest-redshift bin ($0.5 < z < 0.8$) indicates that only one progenitor was in a field environment. We identify three epochs of stellar mass growth; an SF-dominated phase at $z > 2.25$, a phase between $1.25 < z < 2.25$ where SF is responsible for ~50% of the total stellar mass growth and gas-rich and poor mergers are required to match the total (green), and finally a dry merger-dominated phase at $z < 1.25$ where SF is insignificant to the total stellar mass growth.</figcaption></figure>

median total mass growth rate. In this era, progenitors are still forming their mass predominantly through in situ SF independent of environment; however, additional methods of stellar mass generation are required. It is during this time that the total stellar mass delivered by mergers becomes important to BCG progenitor evolution, as the mass formed via SF (red/ blue lines in Figure 6) is insufficient to account for the observed total stellar mass growth rate (green line in Figure 6).

### 5.2.3. The Dry Merger Epoch ($z < 1.25$)

Finally, at $z < 1.25$ the in situ SF of BCG progenitors in dense environments sharply declines with redshift, while massive field galaxies maintain SF activity down to lower redshifts. In our lowest-redshift bin ($z = 0.55$), the limited area of the COSMOS field introduces small-number statistics. In this bin, we observe only one galaxy above our mass cut in a field environment, compared with 16 in dense environments. We see a general trend where SF is more than an order of magnitude insufficient to account for the total stellar mass growth rate. This era requires ex situ mass delivery systems, specifically dry mergers, to be the dominant growth mechanism.

## 5.3. The Onset of Quenching

We contextualize the diminishment of in situ SF toward progenitor growth by examining how rest-frame colors and morphological parameters such as the Sérsic index evolve with redshift. Shown in Figure 7, we plot the rest-frame NUV-r and r-J colors (Ilbert et al. 2013) of our progenitor sample in each of our six redshift bins with Sérsic index from the COSMOS Zurich Morphological Catalog (Sargent et al. 2007; Scarlata et al. 2007) as the colorbar. We use the rest-frame absolute magnitudes provided by the COSMOS2015 catalog.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Evolution of the Illustris-selected BCG progenitor sample in rest-frame NUVrJ space. Quiescent galaxies are selected by the top left quadrant indicated by the black lines. Targets outside this region are considered star-forming. Gray points are used for progenitors identified with a junk flag, used for spurious measurements or for targets too faint to fit a Sérsic index ($I_{AB} > 22.5$; Sargent et al. 2007). The median Sérsic index of the overall population (bottom left of each panel) is constant within errors across all redshift bins with more than one measured index. We find that the majority of progenitors selected by Illustris are star-forming until $z \sim 1.5$, after which our progenitor sample predominantly resides in a quiescent colorspace.</figcaption></figure>

Rest-frame magnitudes were calculated using the observed nearest filter magnitude and a k-correction estimated from the SED fitting performed by Laigle et al. (2016). Our SFR estimates are consistent with their target’s NUVrJ classification, with the star-forming galaxy subpopulation hosting median SFRs 10 times the SFR of the quiescent population in our lowest-redshift bin, and five times the SFR of the quiescent population in our highest-redshift bin. A caveat remains that this diagnostic may misclassify star-forming galaxies as quiescent, and contaminate >10% of quiescent selected galaxies (Ilbert et al. 2013). Therefore we compare the total number of galaxies classified as quiescent to the number of quiescent galaxies above the star-forming main sequence, the correlation between SF and stellar mass indicating the median star-forming activity for a galaxy of a given stellar mass and redshift (e.g., Schreiber et al. 2015). We find that ~3%–5% of progenitors classified as quiescent lie above the star-forming main sequence, potentially the misclassified star-forming galaxies found in Ilbert et al. (2013). This contamination rate is low enough that the results shown in Figure 7 are not affected.

Sérsic indices were estimated using the HST/ACS F814W observations of the COSMOS field. NUVrJ is more sensitive to recent SF evolution than UVJ-derived color measurements (Martin et al. 2007; Davidzon et al. 2017), an important consideration in determining when a galaxy has begun quenching. Position in this color space has also been shown to correlate with sSFR as the galaxy’s star-forming efficiency changes (Martin et al. 2007; Arnouts et al. 2013; Ilbert et al. 2015). We find that the stellar populations of our BCG progenitors transition to a quiescent state between our $0.8 < z < 1.1$ and $1.1 < z < 1.5$ redshift bins, including those without morphological classifications. This is consistent with our results from Figure 6, where in situ SF becomes a negligible contributor to a BCG’s evolution between $1.0 < z < 1.3$. Cotemporal with the transition to a quiescent stellar population, the morphology of the progenitor population also changes. The median Sérsic index in each redshift bin (Figure 7) is constant within the scatter, centered at $\langle n \rangle \sim 2.5$. These median Sérsic indices are consistent with those found for

S0 galaxies in cluster environments (e.g., D’Onofrio et al. 2015), suggesting that galaxies in our progenitor sample have hosted a composite bulge+disk structure since $z \sim 2$.

5.4. Visual Classification of Morphology

To provide context to the light profile evolution observed in Section 5.3, we classify the morphological state of all progenitors imaged in the COSMOS field ACS mosaic (Koekemoer et al. 2007), totaling 1372 of our original sample of 1444. Progenitors not included are outside the bounds of the ACS mosaic. The COSMOS ACS mosaic includes HST/ACS F814W broadband filter observations with a final pixel scale of 0.003 per pixel and a $3\sigma$ surface brightness density of $m_{AB} \sim 25.1$ mag arcsec$^{-2}$ (Wen & Zheng 2016). This spatial resolution corresponds to 183 pc per pixel at $z \sim 0.5$ and 231 pc per pixel at $z \sim 3$. We also supplement this observation set with any publicly available images from HST/WFC3 F105W, F110W, F125W, F140W, and F160W. This subset of 381 targets is most useful in our high-redshift bins as the peak of the stellar continuum is redshifted into the WFC3 IR broad bands. The following results shown in this work are consistent independent of the inclusion or exclusion of targets with WFC3 observations. To ensure consistency of classification, we utilize the morphological classifications scheme and user interface implemented by the CANDELS team (Kartaltepe et al. 2015a; Kocevski 2015) to classify each galaxy’s morphology (spheroid, disk, irregular) as well as interaction class (merger, interaction, non-interacting companion).

The user examines the HST image and labels the progenitor with as many flags as necessary to classify it. Our sample was classified by three authors (K.C.C., J.S.K., K.D.T.), and two of three authors must agree on a designation to assign a galaxy to any of the classification bins. For galaxies that lack the organized, symmetric morphology of a disk and/or spheroid, the irregular flag is used. To provide context to the primary classification, interaction classes are also used for the apparent stage of an interaction. The merger flag is used in highly disturbed cases where the primary and secondary galaxies of the mergers are nearly indistinguishable. The interaction flag is
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Distribution of morphological types (spheroid, disk, irregular, unclassifiable) for our sample of BCG progenitors. To be classified as spheroid, disk, irregular, or unclassifiable two out of three classifiers must agree on the designation and each subsample may include members of others if two labels are used by two classifiers.</figcaption></figure>

used for disturbances such as tidal tails where the interacting galaxies remain independently resolved. For more details on this classification scheme see Kartaltepe et al. (2015a). The distribution of morphologies for our BCG progenitor sample is shown in Figure 8. Our sample is predominantly spheroidal at low redshift, with a growing fraction of unclassifiable progenitors toward high redshift, where the rest-frame optical stellar continuum is redshifted out of the observed band. We see the majority of progenitors are spheroidal, indicating that a significant fraction of our total progenitor sample has already developed the spheroidal structure exhibited by present day BCGs. This appears to be consistent with recent work that identifies a population of massive spheroidals at high redshift (e.g., Cassata et al. 2013).

To compare the environmental effects discussed in Sections 5.2.1–5.2.3, we plot the net classification, or “diskiness” of progenitors in our three environment bins (top half of Figure 9). Our diskiness parameter is defined by the averaging of disk flags (valued as +1) and spheroid flags (valued as −1). For example, a target with two spheroid flags and one disk flag will have a diskiness value of −0.33. Flags indicating disk- or bulge-dominated features are also included as values of ±0.5, respectively. These flags are used to refine the classification of a galaxy with multiple primary flags and identify objects such as bulge-dominated S0 galaxies. We find that low-redshift progenitors have a decreasing likelihood of being labeled a disk with increasing local density. This trend is muddled in our highest-redshift bins, where we do not observe evidence for an environmental dependence. Next, we plot the diskiness parameters of progenitors identified as star-forming and quiescent through their NUVrJ colors (see Section 5.3) in the bottom half of Figure 9. We find a clear offset of average net classification for these two subsamples, where quiescent galaxies are more spheroidal on average than their star-forming counterparts below $z \sim 2.25$.

We next consider how often progenitors are classified as an irregular with redshift, and how this flag is used in combination with other flags such as disk or spheroid. We include example HST/ACS F814W images of the possible irregular combinations in Figure 10. In Figure 11, we plot the fraction of the total

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Top: median diskiness parameter, described in Section 5.4, for our progenitor sample binned by field, intergroup, and dense environments. We identify a morphological dependence on environment at low redshift. However, it is difficult to classify fainter, more dust-obscured galaxies at high redshift. Bottom: median diskiness parameter with redshift for progenitors identified as quiescent (magenta) and star-forming (teal) using their NUVrJ colors. We find the quiescent population evolves to a spheroid-dominated population with time, while star-forming galaxies remain a combination of spheroid and disk. We caution that the ACS F814W images predominantly used here will be dominated by the rest-frame UV past $z > 1$.</figcaption></figure>

sample that has at least two irregular flags per target, identical to the confident irregular subsample in Figure 9. This total irregular sample is divided into six redshift bins with subsamples that include at least one spheroid, disk, or composite (spheroid+disk) flag. As these are confident irregulars (two or more irregular flags), we lower the threshold of additional flags to one to examine how often other characteristics such as disks or spheroids are assigned to this population. The frequency of irregulars in our sample declines with decreasing redshift. Also with redshift, we see the irregular population transition from a diverse range of morphologies to an irregular-spheroid and irregular-composite population toward present day. One result of interest is the lowest-redshift bin, in which disk components are identified in all the irregular candidates, a rare feature (≤5%) in the total sample at $z < 1$.

Finally, we examine the distribution of interaction classifications (Figure 12). To classify the state of interactions, we use flags to identify any ongoing “mergers” in which the structure of both galaxies have been disturbed and combined. An “interaction” flag is used when an interacting pair, identified using tidal features, are visibly distinct from one another. We note that these classifications do not distinguish between major and minor merger pairs. We observe a slowly decreasing interaction fraction with redshift (~25%–15%), with a roughly constant merger fraction of 5% of the total sample. We find a $1\sigma$
COSMOS 728559

z = 0.87

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

COSMOS 429644

z = 1.00

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

COSMOS 258757

z = 0.95

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

COSMOS 832597

z = 0.85

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

**Figure 10.** Example irregular galaxies from the HST/ACS F814W broadband filter mosaic of the COSMOS field, identified by their ID in the COSMOS2015 catalog and photometric redshift (Laigle et al. 2016). Top left: an irregular spheroid with a tail to the northwest of the nucleus. Top right: an irregular disk with a predominantly disk-like structure and a tidal tail to the west of the disk, potentially indicating an ongoing merger. Bottom left: a composite disk+spheroid irregular, with two spiral arms to the south of a spheroid component. Bottom right: a purely irregular galaxy with an asymmetric appearance devoid of spheroidal core or spiral arms.

Poisson confidence interval of ±5% between 0.5 < z < 2.0 and
±3% between 2.0 < z < 3.0. Our results are consistent within
errors of those calculated by Man et al. (2016) for the general
galaxy population; they find a total interaction frequency of
~20% at z ∼ 1 using a lower-mass sample of pairs down to
$10^{10.8} M_{\odot}$. The lack of observed mergers in our lowest-redshift
bin is consistent with a low-redshift, minor merger-dominated
model for BCG growth (e.g., Edwards & Patton 2012). It is also

consistent with the results from Duncan et al. (2019), who found
a merger fraction of <10% for their massive subsample
($10^{10.3} M_{\odot}$) in the CANDELS fields. The lack of interactions
in our highest-redshift bin may be due to the difficulty in
visually identifying low surface brightness features present in
interactions at z > 2. Overall this slowly evolving fraction of
interactions is consistent with a model with fewer mergers at
later times.
<figure><img src="image_16.png" /><figcaption>**Figure 11.** Fraction of our visually classified progenitor sample which were labeled as irregulars, using “only irregular” (pink), “irregular+disk” (blue), “irregular+spheroid” (red), or all three component (green) flags.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 12.** Fraction of our visually classified progenitor sample which included any interaction flags, using “interaction” (green), and “merger” (pink) flags. Mergers are differentiated from interactions by the degree to which the observer can distinguish the two interacting galaxies from each other.</figcaption></figure>

### 5.5. Comparison with the Literature

Several works on massive galaxy evolution and the comparison of progenitor selection methods have been published in recent years (e.g., Zhao et al. 2016; Hill et al. 2017; Torrey et al. 2017), with results consistent with ours. Hill et al. (2017) selected massive galaxy progenitors using the number density evolution functions derived from the Bolshoi simulation (Klypin et al. 2011). Due to the dark-matter-only nature of the Bolshoi simulation, a comparison between our work and that done by Hill et al. (2017) helps illustrate how stellar mass selection functions derived using different combinations of dark and baryonic physics affect measured evolution.

Hill et al. (2017) used a lower stellar mass progenitor selection than ours. Since lower-mass galaxies build their stellar mass over longer timescale on average (mass downsizing; Cowie et al. 1996; Pérez-González et al. 2008), Hill et al.’s selection led to an in situ SF-dominated phase which persisted to lower redshift than ours, until $z \sim 1.75$. However, they found a fast transition to a quenched state from $z < 1.75$ (Figure 1 of Hill et al. 2017). Dark-matter-only simulations producing a lower-mass selection function is hypothesized to

be due to the lack of baryonic physics feedback effects that limit the speed of stellar mass growth, enabling the inclusion of faster evolving, lower-mass progenitors.

The median stellar mass of our sample increases by a factor of ~80% over $0.35 < z < 2.00$, and a factor of two when extrapolated to the present day from $z=2$. This rate of BCG stellar mass growth is slower than the factors of two from $z=1$ shown in BCGs residing in clusters selected from *Spitzer* observations (Lidman et al. 2012; Lin et al. 2013; Shankar et al. 2015). However, our results are in rough agreement with Zhao et al. (2016), who found that BCG progenitors grow by a factor of 2 from $z \sim 2$ to present day. The Zhao et al. sample was selected using a hybrid selection function which selected the most massive galaxies in the 38 most overdense regions identified in the CANDELS UDS field. The Zhao et al. sample also exhibited dry merger-dominated growth at $z < 1$. Moreover, our work is consistent with additional works that indicate a slow rate of growth at low redshift (e.g., Tonini et al. 2012; Hill et al. 2017), shown in BCGs selected via the Maraston (2005) stellar population models or comoving number density cuts. Our observation of a transition to a dry merger-dominated epoch below $z \sim 1$ is also consistent with the *Spitzer* Wide-Area Infrared Extragalactic (SWIRE) Survey (Lonsdale et al. 2003) samples of Webb et al. (2015b). However, above $z \sim 1$, our sample includes progenitors with SFRs an order of magnitude lower than those of Webb et al. This is believed to be primarily due to the difference in depth of the 24 µm imaging in COSMOS versus SWIRE (71 µJy versus ~150 µJy, respectively).

### 5.6. Comparison with the Millennium-II Simulation

As there is not yet a perfect cosmological simulation of the universe, it is worthwhile to examine how choosing a different simulation may change our results. Dark-matter-only simulations such as MS-II and Bolshoi, for example, have a different merger tree history than Illustris, which changes the rate ex situ stellar mass is accumulated. To test this, we retrieve the median cumulative number density of progenitor halos from MS-II (Boylan-Kolchin et al. 2009). To identify BCG progenitors simulated in MS-II, we identify the halo population inhabiting a comoving volume at the same density as observed low-redshift BCGs (Davidzon et al. 2017). We then identify the most massive progenitor halo corresponding to each descendant halo in ascending redshift slices to $z \sim 3$ and measure the cumulative comoving number density corresponding to the median most massive halo in that redshift slice. This number density is used to select mock galaxies from Guo et al. (2013), which represent the baryonic content of the MS-II BCG progenitor halos. The median Guo et al. (2013) sample stellar mass is used as to select observed galaxies in COSMOS based on their stellar mass (teal line in Figure 2).

We fit the SEDs of MS-II-selected progenitors identified in the COSMOS2015 catalog in the same manner as the Illustris-selected progenitors, as described in Section 3. Shown in Figure 13, the in situ growth of massive galaxy progenitors due to SF is at least half of the total growth until $z \sim 1.5$, lower than our findings using Illustris-selected progenitors. However, the total growth rate track has a much higher normalization than Illustris, indicating a later, faster growth history in MS-II. To be physically consistent between SFR and total growth in MS-II-selected progenitors, a much richer merger history is required to deliver mass onto the simulated BCG seeds in addition to their intrinsic SF. These differences between
<figure><img src="image_18.png" /><figcaption>**Figure 13.** Observed stellar mass growth rate for MS-II-selected massive galaxy progenitors in dense regions (BCG progenitors; red) and field regions (blue), and the median growth rate of our sample's stellar mass (green).</figcaption></figure>

Illustris and MS-II may have several causes. One is the difference in feedback physics between the MS-II (Guo et al. 2013) semi-analytical models assigned to their halos and the baryonic physics in Illustris. A less efficient feedback process in MS-II would prevent the slow- or shut-down of highly star-forming progenitors. This faster evolution produces a mass selection function with a larger slope, including galaxies with correspondingly lower SFRs.

## 5.7. Comparison with the Constant Density Method

We also compare our results to the constant cumulative number density selection method (van Dokkum et al. 2010). Following the constant number density selection function plotted in Figure 2, only the most massive progenitors are selected. As discussed in Section 3, a constant number density selection only selects progenitors that inhabit the universe at an equal density, an assumption that the number of galaxies in the universe is constant. As this ignores the effect of mergers, the median progenitor mass is not pulled down by the lower-mass secondaries. We use the same fits described in Section 3 to estimate the average stellar mass accretion via in situ SF in comparison with the evolution of the median stellar mass of the progenitor sample in Figure 14. We find that the combination of slow predicted stellar mass evolution and selection of very massive galaxies that are likely to be star-forming result in a model where SF is dominant until $z < 1$, after which gas poor mergers are necessary.

# 6. Summary and Conclusions

We examine the in situ SF and ex situ stellar deposition components of the total stellar mass growth of progenitors for the present day BCG population. By using FUV-FIR photometric data in the COSMOS field, we fit the SEDs for several samples of BCG progenitors and determine their SFR, sSFR, and stellar mass. Once corrected for stellar mass loss, we identify multiple epochs of BCG evolution based on the prominence, or lack thereof, of in situ SF toward stellar growth. Our estimates of in situ SF include episodes fueled by a progenitor's own gas supply as well as those triggered by gas-rich mergers. Total stellar mass growth also includes the direct contribution of stellar mass by dry mergers. These behaviors

<figure><img src="image_19.png" /><figcaption>**Figure 14.** Observed stellar mass growth rate for massive galaxy progenitors selected by constant number density in dense regions, i.e., BCG progenitors (red) and field regions (blue), and the median growth rate of our sample's stellar mass (green).</figcaption></figure>

are then compared with local environment to test for the presence of any dependence. Our results are as follows.

1. BCG progenitor growth is dominated by in situ SF, either secular or merger driven, until approximately $z \sim 2.25$, after which galaxy mergers (gas rich and gas poor) contribute increasing fractions of ex situ stellar mass growth with time until $z \sim 1.25$. After this time, dry mergers become the dominant growth mechanism.

2. Any dependence on local environment at $z > 1$ is indistinguishable within our errors. The SFR-local density correlation observed today is observed in our $z \sim 0.8$ bin; however, volume effects limit our ability to measure this effect at lower redshifts. We also find that spheroidal progenitors become more common in dense environments ($\frac{\delta}{\delta_{\text{median}}} > 2$) than field environments ($\frac{\delta}{\delta_{\text{median}}} < 1$), which most often host spheroid+disk composites.

3. The era of quenching, either through merger or AGN activity, is encountered by most progenitors by $z \sim 1.00-1.25$. This era is characterized by a shut-down of in situ SF and the transition from a disk+spheroid morphology to an overall spheroidal type.

4. Across all redshifts, progenitors requiring an AGN-fitting component represent a small fraction (~1%) of the total progenitor population. SED3FIT AGNs are most often observed at $z > 1.5$, with 29/37 of our SED3FIT AGN sample observed at high redshift.

5. Our progenitor sample consists of spheroid+disk composites until $z \sim 1.3$, after which we identify an environmental dependence on morphology. Progenitors in dense environments are labeled with a more spheroidal net classification as redshift decreases. This is primarily driven by the change in diskiness of the quiescent population in dense regions. Progenitors in intergroup regions ($1 < \frac{\delta}{\delta_{\text{median}}} < 2$) retain their composite structure. The few galaxies in field environments that meet our stellar mass cut increase in diskiness. Mergers and interactions are no more common than the general
population (<25%), with irregular galaxies exhibiting a greater range of morphological classifications with higher redshift.

We thank the referee for the constructive feedback that improved the clarity of this work. Support for this work was provided by NASA through grants HST-GO-13657.010-A and HST-AR-14298.004-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. Support and observations made with the NASA/ESA Hubble Space Telescope, obtained from the Data Archive at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., were also provided by NASA through grant NNX16AB36G as part of the Astrophysics Data Analysis Program. Spectral energy distribution fitting was performed using the scientific computing resources of the Rochester Institute of Technology ION cluster. B.D. acknowledges financial support from NASA through the Astrophysics Data Analysis Program (ADAP), grant No. NNX12AE20G, and the National Science Foundation, grant No. 1716907. We thank Paul Torrey for his council on evolving number density methodology. We thank the generous advice of Elisabete da Cunha and Stefano Berta on the use of MAGPHYS and SED3FIT, respectively. This research made use of Astropy,<fnref n="7" /> a community-developed core Python package for Astronomy (Robitaille et al. 2013; Price-Whelan et al. 2018). This research made use of the iPython environment (Perez & Granger 2007) and the Python packages SciPy (Jones et al. 2001), NumPy (van der Walt et al. 2011), and Matplotlib (Hunter 2007).

Based on observations made with the NASA *Galaxy Evolution Explorer*. *GALEX* is operated for NASA by the California Institute of Technology under NASA contract NAS5-98034. MegaPrime/MegaCam, a joint project of CFHT and CEA/DAPNIA, at the Canada-France-Hawaii Telescope (CFHT) which is operated by the National Research Council (NRC) of Canada, the Institut National des Sciences de l'Univers of the Centre National de la Recherche Scientifique of France, and the University of Hawaii. Based (in part) on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan. Based on observations obtained with WIRCam, a joint project of CFHT, Taiwan, Korea, Canada, and France. This work is based (in part) on observations made with the *Spitzer Space Telescope*, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA.

This work was supported by a NASA Keck PI Data Award, administered by the NASA Exoplanet Science Institute. Data presented herein were obtained at the W. M. Keck Observatory from telescope time allocated to the National Aeronautics and Space Administration through the agency's scientific partnership with the California Institute of Technology and the University of California. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Mauna Kea has always

had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.

## ORCID iDs

Kevin C. Cooke https://orcid.org/0000-0002-2200-9845

Jeyhan S. Kartaltepe https://orcid.org/0000-0001-9187-3605

Behnam Darvish https://orcid.org/0000-0003-4919-9017

Caitlin M. Casey https://orcid.org/0000-0002-0930-6466

Olivier Le Fèvre https://orcid.org/0000-0001-5891-2596

Mara Salvato https://orcid.org/0000-0001-7116-9303

Nicholas Scoville https://orcid.org/0000-0002-0438-3323

## References

Alberts, S., Pope, A., Brodwin, M., et al. 2016, ApJ, 825, 72

Arnouts, S., Le Floc'h, E., Chevallard, J., et al. 2013, A&A, 558, A67

Aune, S., Boulade, O., Charlot, X., et al. 2003, Proc. SPIE, 4841, 513

Balogh, M. L., McGee, S. L., Wilman, D. J., et al. 2011, MNRAS, 412, 2303

Berta, S., Lutz, D., Santini, P., et al. 2013, A&A, 551, A100

Boulade, O., Charlot, X., Abbon, P., et al. 2003, Proc. SPIE, 4841, 72

Boylan-Kolchin, M., Springel, V., White, S. D. M., Jenkins, A., & Lemson, G. 2009, MNRAS, 398, 1150

Brammer, G. B., van Dokkum, P. G., Franx, M., et al. 2012, ApJS, 200, 13

Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151

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

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

Capak, P., Mobasher, B., Scoville, N. Z., et al. 2011, ApJ, 730, 68

Casey, C. M., Cooray, A., Killi, M., et al. 2017, ApJ, 840, 101

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

Cavagnolo, K. W., Donahue, M., Voit, G. M., & Sun, M. 2009, ApJS, 182, 12

Chabrier, G. 2003, PASP, 115, 763

Chiu, I., Mohr, J., McDonald, M., et al. 2016, MNRAS, 455, 258

Civano, F., Marchesi, S., Comastri, A., et al. 2016, ApJ, 819, 62

Contini, E., Yi, S. K., & Kang, X. 2018, MNRAS, 479, 932

Cooke, K. C., Fogarty, K., Kartaltepe, J. S., et al. 2018, ApJ, 857, 122

Cooke, K. C., O'Dea, C. P., Baum, S. A., et al. 2016, ApJ, 833, 224

Cowie, L. L., Songaila, A., Hu, E. M., & Cohen, J. G. 1996, AJ, 112, 839

da Cunha, E., Charlot, S., & Elbaz, D. 2008, MNRAS, 388, 1595

D'Onofrio, M., Marziani, P., & Buson, L. 2015, FrASS, 2, 4

Darvish, B., Martin, C., Gonçalves, T. S., et al. 2018, ApJ, 853, 155

Darvish, B., Mobasher, B., Martin, D. C., et al. 2017, ApJ, 837, 16

Darvish, B., Mobasher, B., Sobral, D., Scoville, N., & Aragon-Calvo, M. 2015, ApJ, 805, 121

Darvish, B., Mobasher, B., Sobral, D., et al. 2016, ApJ, 825, 113

Davidzon, I., Bolzonella, M., Coupon, J., et al. 2013, A&A, 558, A23

Davidzon, I., Ilbert, O., Faisst, A. L., Sparre, M., & Capak, P. L. 2018, ApJ, 852, 107

Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, A&A, 605, A70

De Lucia, G., & Blaizot, J. 2007, MNRAS, 375, 2

Delvecchio, I., Smolčić, V., Zamorani, G., et al. 2017, A&A, 602, A3

Donley, J. L., Koekemoer, A. M., Brusa, M., et al. 2012, ApJ, 748, 142

Dressler, A. 1980, ApJ, 236, 351

Duncan, K., Conselice, C. J., Mundy, C., et al. 2019, ApJ, 876, 110

Edwards, L. O. V., & Patton, D. R. 2012, MNRAS, 425, 287

Ehlert, S., der Linden, A. v., Allen, S. W., et al. 2013, MNRAS, 437, 1942

Elbaz, D., Daddi, E., Le Borgne, D., et al. 2007, A&A, 468, 33

Elbaz, D., Dickinson, M., Hwang, H. S., et al. 2011, A&A, 533, A119

Elvis, M., Civano, F., Vignali, C., et al. 2009, ApJS, 184, 158

Fabian, A. C. 1994, ARA&A, 32, 277

Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004, ApJS, 154, 10

Feltre, A., Hatziminaoglou, E., Fritz, J., & Franceschini, A. 2012, MNRAS, 426, 120

Fogarty, K., Postman, M., Li, Y., et al. 2019, arXiv:1905.01377

Fraser-McKelvie, A., Brown, M. J. I., & Pimbblet, K. A. 2014, MNRAS, 444, L63

Genel, S., Vogelsberger, M., Springel, V., et al. 2014, MNRAS, 445, 175

Gozaliasl, G., Finoguenov, A., Khosroshahi, H. G., et al. 2018, MNRAS, 475, 2787

Guo, Q., White, S., Angulo, R. E., et al. 2013, MNRAS, 428, 1351

Hasinger, G., Capak, P., Salvato, M., et al. 2018, ApJ, 858, 77

<footnote n="7">7 http://www.astropy.org</footnote>
Hill, A. R., Muzzin, A., Franx, M., et al. 2017, ApJ, 837, 147

Hlavacek-Larrondo, J., McDonald, M., Benson, B. A., et al. 2015, ApJ, 805, 35

Hunter, J. D. 2007, CSE, 9, 90

Ilbert, O., Arnouts, S., Le Floc'h, E., et al. 2015, A&A, 579, A2

Ilbert, O., Capak, P., Salvato, M., et al. 2008, ApJ, 690, 1236

Ilbert, O., McCracken, H. J., Le Fevre, O., et al. 2013, A&A, 556, A55

Jaacks, J., Finkelstein, S. L., & Nagamine, K. 2016, ApJ, 817, 174

Jones, E., Oliphant, T., & Peterson, S. P. 2001, SciPy: Open Source Scientific Tools for Python, http://www.scipy.org

Kartaltepe, J. S., Dickinson, M., Alexander, D. M., et al. 2012, ApJ, 757, 23

Kartaltepe, J. S., Mozena, M., Kocevski, D., et al. 2015a, ApJS, 221, 11

Kartaltepe, J. S., Sanders, D. B., Silverman, J. D., et al. 2015b, ApJL, 806, L35

Kashino, D., Silverman, J. D., Rodighiero, G., et al. 2013, ApJL, 777, L8

Kauffmann, G., White, S. D. M., Heckman, T. M., et al. 2004, MNRAS, 353, 713

Klypin, A. A., Trujillo-Gomez, S., & Primack, J. 2011, ApJ, 740, 102

Kocevski, D. D. 2015, MorphGUI, Zenodo, doi:10.5281

Koekemoer, A. M., Aussel, H., Calzetti, D., et al. 2007, ApJS, 172, 196

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

Krolik, J. H., & Begelman, M. C. 1988, ApJ, 329, 702

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

Lanzuisi, G., Delvecchio, I., Berta, S., et al. 2017, A&A, 602, A123

Lauer, T. R., Postman, M., Strauss, M. A., Graves, G. J., & Chisari, N. E. 2014, ApJ, 797, 82

Le Fèvre, O., Tasca, L. A. M., Cassata, P., et al. 2015, A&A, 576, A79

Lee, N., Le Floc'h, E., Sanders, D. B., et al. 2010, ApJ, 717, 175

Leja, J., van Dokkum, P., & Franx, M. 2013, ApJ, 766, 33

Lidman, C., Suherli, J., Muzzin, A., et al. 2012, MNRAS, 427, 550

Lilly, S. J., Fevre, O. L., Renzini, A., et al. 2007, ApJS, 172, 70

Lin, Y.-T., Brodwin, M., Gonzalez, A. H., et al. 2013, ApJ, 771, 61

Liu, F. S., Mao, S., & Meng, X. M. 2012, MNRAS, 423, 422

Livermore, R. C., Jones, T., Richard, J., et al. 2012, MNRAS, 427, 688

Livermore, R. C., Jones, T. A., Richard, J., et al. 2015, MNRAS, 450, 1812

Lonsdale, C. J., Smith, H. E., Rowan-Robinson, M., et al. 2003, PASP, 115, 897

Lotz, J. M., Jonsson, P., Cox, T. J., et al. 2011, ApJ, 742, 103

Loubser, S. I., & Sánchez-Blázquez, P. 2012, MNRAS, 425, 841

Loubser, S. I., Sánchez-Blázquez, P., Sansom, A. E., & Soechting, I. K. 2009, MNRAS, 398, 133

Lutz, D., Poglitsch, A., Altieri, B., et al. 2011, A&A, 532, A90

Macuga, M., Martini, P., Miller, E. D., et al. 2019, ApJ, 874, 54

Man, A. W. S., Zirm, A. W., & Toft, S. 2016, ApJ, 830, 89

Maraston, C. 2005, MNRAS, 362, 799

Marchesi, S., Civano, F., Elvis, M., et al. 2016, ApJ, 817, 34

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

Markowitz, A. G., Krumpe, M., & Nikutta, R. 2014, MNRAS, 439, 1403

Martin, D. C., Fanson, J., Schiminovich, D., et al. 2005, ApJL, 619, L1

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

Martini, P., Miller, E. D., Brodwin, M., et al. 2013, ApJ, 768, 1

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

Momcheva, I. G., Brammer, G. B., van Dokkum, P. G., et al. 2016, ApJS, 225, 27

Morishita, T., Ichikawa, T., Noguchi, M., et al. 2015, ApJ, 805, 34

Morrissey, P., Conrow, T., Barlow, T. A., et al. 2007, ApJS, 173, 682

Moustakas, J., Coil, A. L., Aird, J., et al. 2013, ApJ, 767, 50

Mundy, C. J., Conselice, C. J., & Ownsworth, J. R. 2015, MNRAS, 450, 3696

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

Narayanan, D., Dey, A., Hayward, C. C., et al. 2010, MNRAS, 407, 1701

Narayanan, D., Turk, M., Feldmann, R., et al. 2015, Natur, 525, 496

O'Dea, C. P., Baum, S. A., Privon, G., et al. 2008, ApJ, 681, 1035

O'Dea, K. P., Quillen, A. C., O'Dea, C. P., et al. 2010, ApJ, 719, 1619

Oemler, A., Jr. 1976, ApJ, 209, 693

Oliver, S. J., Bock, J., Altieri, B., et al. 2012, MNRAS, 424, 1614

Patel, S. G., van Dokkum, P. G., Franx, M., et al. 2013, ApJ, 766, 15

Perez, F., & Granger, B. E. 2007, CSE, 9, 21

Pérez-González, P. G., Rieke, G. H., Villar, V., et al. 2008, ApJ, 675, 234

Perna, M., Brusa, M., Salvato, M., et al. 2015, A&A, 583, A72

Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, 518, L1

Pillepich, A., Madau, P., & Mayer, L. 2015, ApJ, 799, 184

Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2

Postman, M., Coe, D., Benítez, N., et al. 2012, ApJS, 199, 25

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

Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M., et al. 2018, AJ, 156, 123

Ragone-Figueroa, C., Granato, G. L., Ferraro, M. E., et al. 2018, MNRAS, 479, 1125

Rawle, T. D., Edge, A. C., Egami, E., et al. 2012, ApJ, 747, 29

Rieke, G. H., Young, E. T., Engelbracht, C. W., et al. 2004, ApJS, 154, 25

Robitaille, T. P., Tollerud, E. J., Greenfield, P., et al. 2013, A&A, 558, A33

Rodriguez-Gomez, V., Sales, L. V., Genel, S., et al. 2017, MNRAS, 467, 3083

Sanders, D. B., & Mirabel, I. F. 1996, ARA&A, 34, 749

Sanders, D. B., Salvato, M., Aussel, H., et al. 2007, ApJS, 172, 86

Sargent, M. T., Carollo, C. M., Lilly, S. J., et al. 2007, ApJS, 172, 434

Scarlata, C., Carollo, C. M., Lilly, S., et al. 2007, ApJS, 172, 406

Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525

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

Scoville, N., Arnouts, S., Aussel, H., et al. 2013, ApJS, 206, 3

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

Shankar, F., Buchan, S., Rettura, A., et al. 2015, ApJ, 802, 73

Shi, Y., Rieke, G. H., Hines, D. C., et al. 2006, ApJ, 653, 127

Silverman, J. D., Kashino, D., Sanders, D., et al. 2015, ApJS, 220, 12

Smolčić, V., Novak, M., Delvecchio, I., et al. 2017, A&A, 602, A6

Smolčić, V., Schinnerer, E., Scodeggio, M., et al. 2008, ApJS, 177, 14

Springel, V., White, S. D. M., Jenkins, A., et al. 2005, Natur, 435, 629

Steinhardt, C. L., Speagle, J. S., Capak, P., et al. 2014, ApJL, 791, L25

Sutherland, W., Emerson, J., Dalton, G., et al. 2015, A&A, 575, A25

Taniguchi, Y., Kajisawa, M., Kobayashi, M. A. R., et al. 2015, PASJ, 67, 104

Taniguchi, Y., Scoville, N., Murayama, T., et al. 2007, ApJS, 172, 9

Tonini, C., Bernyk, M., Croton, D., Maraston, C., & Thomas, D. 2012, ApJ, 759, 43

Torrey, P., Wellons, S., Ma, C.-P., Hopkins, P. F., & Vogelsberger, M. 2017, MNRAS, 467, 4872

Torrey, P., Wellons, S., Machado, F., et al. 2015, MNRAS, 454, 2770

Trakhtenbrot, B., Civano, F., Urry, C. M., et al. 2016, ApJ, 825, 4

Tremblay, G. R., Gladders, M. D., Baum, S. A., et al. 2014, ApJL, 790, L26

Trump, J. R., Impey, C. D., Elvis, M., et al. 2009, ApJ, 696, 1195

Trump, J. R., Impey, C. D., McCarthy, P. J., et al. 2007, ApJS, 172, 383

Umehata, H., Tamura, Y., Kohno, K., et al. 2015, ApJL, 815, L8

van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, CSE, 13, 22

van Dokkum, P. G., Leja, J., Nelson, E. J., et al. 2013, ApJL, 771, L35

van Dokkum, P. G., Whitaker, K. E., Brammer, G., et al. 2010, ApJ, 709, 1018

Vogelsberger, M., Genel, S., Springel, V., et al. 2014, MNRAS, 444, 1518

Webb, T., Noble, A., DeGroot, A., et al. 2015a, ApJ, 809, 173

Webb, T. M. A., Lowenthal, J., Yun, M., et al. 2017, ApJL, 844, L17

Webb, T. M. A., Muzzin, A., Noble, A., et al. 2015b, ApJ, 814, 96

Wen, Z. Z., & Zheng, X. Z. 2016, ApJ, 832, 90

Werner, M. W., Roellig, T. L., Low, F. J., et al. 2004, ApJS, 154, 1

Yang, G., Brandt, W. N., Darvish, B., et al. 2018, MNRAS, 480, 1022

Yoshikawa, T., Akiyama, M., Kajisawa, M., et al. 2010, ApJ, 718, 112

Zamojski, M. A., Schiminovich, D., Rich, R. M., et al. 2007, ApJS, 172, 468

Zhao, D., Conselice, C. J., Aragón-Salamanca, A., et al. 2016, MNRAS, 464, 1393