The mass, colour, and structural evolution of today's massive galaxies since z~5

TL;DR

This study traces the mass, colour, and structural evolution of massive galaxies from redshift 5 to the present, revealing early core formation, ongoing outer growth, and a transition in growth mechanisms around redshift 1.5.

Contribution

It provides a detailed analysis of the evolution of massive galaxies using stacking and a novel progenitor selection method, highlighting the transition in growth processes over cosmic time.

Findings

Most stellar mass in centers was in place by z~2.

Outer regions experience largest fractional mass increase.

A transition at z~1.5 from in-situ star formation to minor mergers.

Abstract

In this paper, we use stacking analysis to trace the mass-growth, colour evolution, and structural evolution of present-day massive galaxies ($\log(M_{*}/M_{\odot})=11.5$) out to $z=5$. We utilize the exceptional depth and area of the latest UltraVISTA data release, combined with the depth and unparalleled seeing of CANDELS to gather a large, mass-selected sample of galaxies in the NIR (rest-frame optical to UV). Progenitors of present-day massive galaxies are identified via an evolving cumulative number density selection, which accounts for the effects of merging to correct for the systematic biases introduced using a fixed cumulative number density selection, and find progenitors grow in stellar mass by $\approx1.5~\mathrm{dex}$ since $z=5$. Using stacking, we analyze the structural parameters of the progenitors and find that most of the stellar mass content in the central regions was in place by $z\sim2$, and while galaxies continue to assemble mass at all radii, the outskirts experience the largest fractional increase in stellar mass. However, we find evidence of significant stellar mass build up at $r<3~\mathrm{kpc}$ beyond $z>4$ probing an era of significant mass assembly in the interiors of present day massive galaxies. We also compare mass assembly from progenitors in this study to the EAGLE simulation and find qualitatively similar assembly with $z$ at $r<3~\mathrm{kpc}$. We identify $z\sim1.5$ as a distinct epoch in the evolution of massive galaxies where progenitors transitioned from growing in mass and size primarily through in-situ star formation in disks to a period of efficient growth in $r_{e}$ consistent with the minor merger scenario.

Tables

Table 1. Table 1 Number of galaxies in each redshift range by catalog

$z$-range	$UVISTA$	$3DHST$
	16	0
	56	5
	96	22
	166	31
	276	79
	466	104
	160	69
	…	110
	…	154**We are incomplete in mass for this point

Equations

$$r_e=a{M}_{\ast }^b$$

$$r_e=\gamma M_{\ast }^{\alpha }{\left(1+\frac{M_{\ast }}{M_0}\right)}^{\alpha -\beta }$$

Full text

1. Introduction

The mass growth and structural evolution of today’s most massive galaxies is an important tracer of galaxy assembly at early times. These systems are host to the oldest stars, suggesting they were the first galaxies to assemble. Because they are the oldest systems, their progenitors can theoretically be traced to higher redshifts than their low mass counterparts and can be studied from the onset of re-ionization to give a complete history of galactic evolution. Additionally, the most massive systems tend to be the most luminous, and they are the easiest to observe at high redshift with high fidelity. Massive galaxies also provide important constraints on the physics involved in cosmological simulations, as they impose upper limits on growth as well as the efficiency of various feedback mechanisms such as active galactic nuclei, mergers and supernovae.

Today’s massive ($\log{\mathrm{M_{*}/M_{\odot}}}\sim 11.5$) galaxies, to first order, are a uniform population. They are homogeneous in morphology and star formation, appearing spheroidal and have low specific star formation rates, and high quiescent fractions (e.g., Thomas et al., 2005; Gallazzi et al., 2005; Kuntschner et al., 2010; Thomas et al., 2010; Cappellari et al., 2011; Mortlock et al., 2013; Moustakas et al., 2013; Ilbert et al., 2013; Muzzin et al., 2013b; Davis et al., 2014; McDermid et al., 2015). In contrast to today’s massive galaxies, massive galaxies at high redshift show increasing diversity (e.g., Franx et al., 2008; van Dokkum et al., 2011). With increasing redshift, massive galaxies become increasingly star forming (e.g., Papovich et al., 2006; Kriek et al., 2008; van Dokkum et al., 2010; Brammer et al., 2011; Bruce et al., 2012; Ilbert et al., 2013; Muzzin et al., 2013b; Patel et al., 2013; Stefanon et al., 2013; Barro et al., 2014; Duncan et al., 2014; Marchesini et al., 2014; Toft et al., 2014; van Dokkum et al., 2015; Barro et al., 2016; Man et al., 2016; Tomczak et al., 2016), and the massive galaxies which are identified as quiescent at high redshift are structurally distinct from their low redshift counterparts as seen in their small effective radii ($r_{e}$) and more centrally concentrated stellar-mass density profiles (Daddi et al., 2005; Trujillo et al., 2006; Toft et al., 2007; Cimatti et al., 2008; van Dokkum et al., 2008; Damjanov et al., 2009; Newman et al., 2010; Szomoru et al., 2010; Williams et al., 2010; van de Sande et al., 2011; Bruce et al., 2012; Muzzin et al., 2012; Oser et al., 2012; Szomoru et al., 2012, 2013; McLure et al., 2013; van de Sande et al., 2013; Newman et al., 2015; Straatman et al., 2015; Hill et al., 2016).

Although the central regions of massive galaxies contain a higher fraction of the total mass at high redshift, their central stellar densities show remarkably little evolution between $z\approx 2-3$ and $z=0$ (e.g., Bezanson et al., 2009; van Dokkum et al., 2010; Toft et al., 2012; van de Sande et al., 2013; Patel et al., 2013; Belli et al., 2014a; van Dokkum et al., 2014; Williams et al., 2014; Whitaker et al., 2016) with the majority of stellar-mass build-up occurring in the outer regions (with galaxies growing in an ‘inside-out’ fashion). This mass assembly is thought to occur via minor, dissipation-less mergers; a scenario which is able to account for the size growth, while leaving the interior regions relatively undisturbed (e.g., Bezanson et al., 2009; Naab et al., 2009; Hopkins et al., 2010; Trujillo et al., 2011; Newman et al., 2012; Hilz et al., 2013; McLure et al., 2013). The aims of the present study are to determine whether these trends continue to high redshifts and to identify the epoch when galaxies’ central regions assemble their mass.

Obtaining a census of massive galaxies across a broad redshift range is technically challenging, as they have low number densities on the sky (Cole et al., 2001; Bell et al., 2003; Conselice et al., 2005; Marchesini et al., 2009; Bezanson et al., 2011; Caputi et al., 2011; Baldry et al., 2012; Ilbert et al., 2013; Muzzin et al., 2013b; Duncan et al., 2014; Tomczak et al., 2014; Caputi et al., 2015; Stefanon et al., 2015; Huertas-Company et al., 2016) and their rest-frame optical emission shifts into the near-infrared (NIR) at intermediate redshifts. To study the evolution of massive galaxies across cosmic time, as a population, necessitates deep and wide NIR surveys to both probe large volumes and obtain rest-frame optical emission to significant signal-to-noise (S/N).

In this study we use stacking analysis to obtain high-fidelity profiles of the progenitors of massive galaxies out to significant radii (at low $z$, $r>60~{}\mathrm{kpc}$). We take advantage of the unparalleled combination of depth and area in the third data release of the UltraVISTA survey (McCracken et al., 2012) to study the structural evolution of massive galaxies out to $z=3.5$. Due to incompleteness in UltraVISTA at the highest redshifts considered in this study, we also use the deeper CANDELS F160W data from the 3DHST photometric catalogs (Brammer et al., 2012; Skelton et al., 2014; Momcheva et al., 2016) to extend the redshift coverage to $z=5$. This is a significant gain in redshift over previous studies, and provides the most extensive redshift range over which the profiles of massive galaxies have been traced.

2. Sample Selection

2.1. Number-density selection

Linking the progenitors of present day galaxies to their high redshift counterparts is challenging, as the merger and star formation history (SFH) of any individual galaxy is not well constrained. One way to circumvent these issues is to assume that galaxies maintain rank-order across cosmic time (i.e., the most massive galaxies today will have been the most massive galaxies yesterday, cosmologically speaking). This assumption predicts a constant co-moving number-density with redshift, an outcome used by van Dokkum et al. (2010) to trace the mass and size growth of galaxies from $z=2$ (corresponding to $n=2\times 10^{-4}~{}\mathrm{Mpc^{-3}dex^{-1}}$). Subsequent studies have used the same assumptions to select progenitors based on a constant cumulative number density (e.g., Bezanson et al., 2011; Brammer et al., 2011; Papovich et al., 2011; Fumagalli et al., 2012; van Dokkum et al., 2013; Patel et al., 2013; Ownsworth et al., 2014; Morishita et al., 2015), which has the advantage over its non-cumulative counterpart of being single valued in mass.

The selection of progenitors and their descendants at a constant cumulative number density implicitly assumes that mergers and in-situ star-formation do not broadly effect rank-order, an assumption which has been shown to result in systematically biased progenitor selection (Behroozi et al., 2013; Leja et al., 2013; Torrey et al., 2015). To account for the affects of mergers on the progenitor mass, we utilize an evolving cumulative number density selection following the prescription of Behroozi et al. (2013), who use halo-abundance matching within a $\Lambda\mathrm{CDM}$ cosmology to connect progenitors and their descendants. It is important to note, that we have used the prescription to trace progenitors of low redshift massive galaxies, not the descendants of high redshift massive galaxies, of which the former yields a steeper evolution in cumulative number density due to the shape of the halo mass function, and scatter in mass accretion histories (see Behroozi et al., 2013; Leja et al., 2013).

2.2. The implied stellar mass growth of the progenitors of massive galaxies since $z\sim 5$

In Fig. 1 we show the integrated Schecter fits of the mass functions of Muzzin et al. (2013b) between $0.2<z<3.0$, and Grazian et al. (2015) between $3.5<z<5.5$. These mass functions are based on photometric redshifts determined via ground and space based NIR imaging from the UltraVISTA and CANDELS surveys respectively. In the left-panel of Fig. 1, we show our evolving cumulative number density selection based on the abundance matching of Behroozi et al. (2013). The masses implied from a fixed-cumulative number density selection are also shown to illustrate the effect of the bias when the effects of mergers are ignored in the selection. In the right-panel of Fig. 1 the implied progenitor masses from the left-panel are plotted for both the fixed and evolving cumulative number density selection, as a function of redshift. The error bars are the uncertainties from the mass functions, which take into account the uncertainties in the photometric redshifts, SFHs, and cosmic variance. The solid grey region represents the scatter in the number densities from the abundance matching of Behroozi et al. (2013), and the hatched regions illustrate an estimate of the mass completeness which is discussed in detail in Sec. 2.3.

Below $z=2$, Fig. 1 shows that both constant and evolving cumulative number density selections yield progenitor masses which are consistent within the uncertainties in the mass-functions. However beyond $z=2$, the bias in the fixed cumulative number density becomes significant, and over-predicts the median progenitor mass. Using the abundance matching technique, we see an overall increase in stellar mass of $1.5~{}\mathrm{dex}$ since $z\sim 5$. Our fractional mass growth out to $z=3$ is consistent within the uncertainties with Marchesini et al. (2014) who use the same abundance matching selection for ultra-massive $\log(M_{*}/M_{\odot})\sim 11.8$) descendants, and with Ownsworth et al. (2014), who use a constant cumulative number density selection which is corrected for major mergers to trace progenitors. Using their correction, they find $75\pm 9\%$ of the descendant mass is assembled after $z=3$, which is consistent with $\sim 80\%$ which we find in the current study.

We note that in Fig. 1 we have selected a progenitor mass for a redshift bin between $3.0<z<3.5$ (orange point), even though we have indicated no mass-function for this redshift. The mass-function from Muzzin et al. (2013b) for this redshift range proved to be unreliable for the mass considered due to incompleteness from UltraVISTA DR1 (the source catalog used in generating the mass functions). However, with the deeper exposures from the third data release (DR3) of UltraVISTA, we are complete to the progenitor masses considered out to $z=3.5$. To calculate the expected progenitor mass between $3.0<z<3.5$, we linearly interpolated the mass between adjacent redshift bins. We also observe a trend of the uncertainties in the mass function monotonically increasing from low to high redshift. Thus, we similarly linearly interpolated the uncertainties to estimate the uncertainty in mass for $3.0<z<3.5$ due to uncertainties in photo-$z$, SFH and cosmic variance. We also use the uncertainties in the progenitor mass selection as the upper and lower mass bounds for the galaxies that contribute to the resulting stack, thus we select a larger range of masses at higher redshift, than at lower redshift.

It has been shown that the Behroozi et al. (2013) prescription for selecting progenitors performs well in terms of recovering the average stellar mass of the progenitors of present-day high-mass galaxies, however this method fails in capturing the diversity in mass of all progenitors as implied by simulations (e.g., Torrey et al., 2015; Clauwens et al., 2016; Wellons & Torrey, 2016), which also predict that the scatter in progenitor masses tends to increase with redshift. Given this large scatter, there is no guarantee that the evolution of other galaxy properties, such as size, will follow from the Behroozi et al. (2013) selection. However, in an upcoming paper (Clauwens et al., in prep) we will show that for the property of interest in our study (i.e. the average radial build-up of stellar mass for the progenitors of massive galaxies), the Behroozi et al. (2013) selection yields average agreement with progenitors within the EAGLE simulation.

2.3. Data

2.3.1 UltraVISTA

In order to study the evolution of the average properties of massive galaxies, it was necessary to utilize both wide field ground-based, and deep space-based imaging for our stacking analysis. Massive galaxies ($\log(M_{*}/M_{\odot})\sim 11$) are exceedingly rare objects, with low number densities ($\sim 10^{-5}~{}\mathrm{Mpc^{-3}}$) on the sky (e.g., Cole et al., 2001; Bell et al., 2003; Baldry et al., 2012; Muzzin et al., 2013b; Ilbert et al., 2013; Tomczak et al., 2014; Caputi et al., 2015; Stefanon et al., 2015), and require wide-field surveys to characterize a significant population. To that end, we utilize the NIR imaging from the DR3 of the UltraVISTA survey (McCracken et al., 2012) for our stacking analysis.

The DR3 UltraVISTA catalog (Muzzin et al., in prep) is a K-selected, multi-band catalog constructed from the UltraVISTA survey. Briefly, the survey covers the COSMOS field with a total area of $1.7~{}\mathrm{deg^{2}}$, with deep imaging in the $Y,J,H$ and $Ks$ bands. The survey also contains ultra-deep stripes with longer exposures which cover a $0.75~{}\mathrm{deg^{2}}$ area, and also includes imaging in the VISTA NB118 NIR filter (Milvang-Jensen et al., 2013). The newest data release is constructed with the same techniques as the DR1 30-band catalog (Muzzin et al., 2013a), with the inclusion of new and higher-quality data to determine photo-z’s, and stellar population parameters. The DR3 survey depths in the ultra-deep stripes are $\sim 1.4$ magnitudes deeper than DR1 (with $5\sigma$ limiting magnitudes in the ultra-deep regions of 25.7, 25.4, 25.1, and 24.9 in $Y,J,H$ and $Ks$).

Several other datasets have also been added since the first data release including 5 CFHTLS filters, $u^{*}g^{\prime}r^{\prime}i^{\prime}z^{\prime}$, as well as 2 new Subaru narrow bands (NB711, NB816). Most importantly for this analysis we also include the latest data from SPLASH (Capak et al., 2012) and SMUVS (PI Caputi; Ashby et al., in prep). These are post-cryo Spitzer-IRAC observations that improve the $\mathrm{[3.6]}$ and $\mathrm{[4.5]}$ depth from 23.9 to 25.3. Overall this is a $38$-band catalog (compared to 30 in Muzzin et al. 2013a), and the substantial increase in depth in the $Y,J,H,Ks$, $\mathrm{[3.6]}$ and $\mathrm{[4.5]}$ bands make it a powerful dataset for studying massive galaxies at intermediate and high redshifts.

In the right panel of Fig. 1 we have indicated our estimated mass completeness limits with the filled hatched regions. To estimate our mass completeness at $z<4$, we used the limits on the mass functions from Muzzin et al. (2013b) (which were derived using UltraVISTA DR1), and adjusted the mass limit according to the gain in K-band depth (the K-band limit is 1.5 magnitudes deeper between DR1 and DR3) assuming a constant mass-to-light ratio. Since galaxy mass-to-light ratios decrease with redshift (e.g. van de Sande et al., 2015), this likely represents a conservative estimate of the limiting mass at high redshifts.

2.3.2 CANDELS

As UltraVISTA DR3 is only mass complete for our selection out to $z=3.5$, we use the reddest band available from CANDELS in order to explore redshifts unobtainable through UltraVISTA. We select galaxies using the photometric data products from the 3DHST survey (Brammer et al., 2012; Skelton et al., 2014) from all 5 CANDELS fields. As an estimate of our mass completeness in CANDELS, we adopt the limiting mass derived from the $75\%$ magnitude completeness limit ($F160W=25.9$) in the shallower pointings in the GOODS-S and UDS fields as described in Grazian et al. (2015). They estimated their mass completeness using the technique of Fontana et al. (2004), which assumes the distribution of mass-to-light ratios immediately above the magnitude limit holds at slightly lower fluxes, and compute the fraction of objects lost due to large mass-to-light ratios. The estimated completeness for CANDELS is indicated in the right panel of Fig. 1 as the grey cross-hatched region.

Although the aforementioned estimates of mass completeness take into account galaxies with varied mass-to-light ratios, it is worth stressing inherent uncertainties when determining mass limits at high redshift. At $z>3.5$, we increasingly rely on photometric redshifts, as high-fidelity spectroscopic redshifts are fewer in number (Grazian et al., 2015). In addition, sub-mm galaxies (SMGs) likely account for at least a fraction of the progenitors of massive galaxies at high redshift (e.g., Toft et al., 2014), and they have been shown to have high optical extinction (e.g., Swinbank et al., 2010; Couto et al., 2016). As the progenitors selected at $z>3.5$ of this study tend to be less massive than a typical SMG, we do not expect that they will form a significant fraction of the sample. However, we cannot rule out a tail of less, but still obscured sources to lower masses in the distribution of SMGs. This would have the effect of biasing our high redshift progenitor selection to bluer, less-obscured sources.

Table 1 provides a summary of the number of galaxies in the given redshift range, at the implied mass as determined from our evolving cumulative number density selection (see Sec. 2.1) from both the UltraVISTA and 3DHST catalogs. In order to boost the number of galaxies in UltraVISTA, we have used galaxies from both the deep (DR1) and ultra-deep (DR3) catalog out to $2.0<z<2.5$ where we are complete in mass for the shallower catalog (DR1). For the $3.0<z<3.5$ bin, we have only utilized the DR3 catalog, as we are incomplete in DR1. As evident from Table 1, UltraVISTA has a larger population of massive galaxies at low redshift, while there are 0 galaxies in all 5 CANDELS fields which are massive ($\log(M_{*}/M_{\odot})\sim 11.5$) at $z=0.35$, and only 5 galaxies in the next highest redshift bin. However CANDELS is crucial to continue the progenitor selection beyond $z>3.5$ as we are mass incomplete in this region with UltraVISTA. Additionally, as galaxies had smaller $r_{e}$ at high redshift (see discussion in Sec. 1 and references therein), the space-based seeing of CANDELS is necessary to properly map the density profiles at these epochs. Thus we utilize both data sets in our analysis.

3. Rest-Frame Color Evolution

Cumulative number density selection is a method which selects solely on stellar mass, and is therefore blind to other galaxy properties such as levels of star-formation activity. A simple, but effective way to establish star-forming activity in a population of galaxies is to observe where they are located in rest-frame $U-V$ and $V-J$ color space, commonly referred to as a $UVJ$-diagram. First proposed by Labbé et al. (2005), it is observed that galaxies exhibit a bi-modality in rest-frame $UVJ$ colour space which is correlated with the level of obscured and unobscured star formation. Actively star-forming and quiescent galaxies separate into a ‘blue’ and ‘red’ sequence in the $UVJ$-diagram (e.g., Williams et al., 2009, 2010; Whitaker et al., 2011; Fumagalli et al., 2014; Yano et al., 2016).

In Fig. 2, we plot the rest-frame $U-V$ and $V-J$ colours for all redshift bins to provide a diagnostic of star-formation activity within each stack. Each of the nine panels represents a different redshift range, with galaxy masses selected according to their expected evolving cumulative number density (see Fig. 1). The first seven panels are galaxies from UltraVISTA DR3, and the last two panels contain galaxies from the 3DHST photometric catalog. It is important to note that we are mass incomplete for the $4.5<z<5.5$ bin (see Fig. 1). However we have chosen to include it as part of our analysis, with the caveat that we are likely biased towards bluer galaxies. Overlaid in each panel are the colour selections used by Muzzin et al. (2013b) to separate quiescent and star forming sequences.

As one progresses in redshift, it becomes apparent from Fig. 2 that the number of galaxies selected dramatically increases. This is a result of two competing effects. The first, is that the size of our mass range becomes progressively larger with redshift, as seen in the error bars on the right-panel of Fig. 1. By selecting in a wider mass range, we will inevitably select more galaxies. The second effect is that as the number densities of progenitors increases with redshift, we are progressing towards the lower mass end of the mass-functions (Ilbert et al., 2013; Muzzin et al., 2013b; Grazian et al., 2015). Thirdly, at low redshift, the probed co-moving volume is also smaller than at high redshift. The combined effect is to have our lowest redshift, and least populated stack contain only 16 galaxies, whereas our most populated stack at $2.5<z<3.0$ contains 276 objects (Table 1).

The most prominent trend in Fig. 2 comes in the colour evolution of the progenitors across redshift. They begin very blue in both $U-V$ and $V-J$ in the lower-left of the star-forming sequence and progress red-ward along the star-forming sequence to the upper-right until $2.5<z<3.0$ before reddening in $V-J$ and joining the quiescent sequence. Assuming our number density selection is valid, this represents a true evolution in $UVJ$ colour.

Fig. 3 show the average $UVJ$ colour evolution for each redshift bin, separated into star-forming and quiescent progenitors, and highlights explicitly the trends observed in Fig. 2. In this figure, we see most of the early ($z>3$) colour evolution is driven by the star-forming progenitors. At $z<3$, star-forming progenitors are beginning to quench in large numbers and the two tracks are broadly parallel until $z<1$ where the quiescent progenitor fractions are high, and $UVJ$ colour evolution is driven by the quiescent progenitors. This seems to indicate that massive galaxies begin their existence as star forming galaxies, which progress along the blue sequence (via aging of the stellar populations, and increase in stellar mass through star formation), before quenching and joining the red sequence.

The progression in the $UVJ$-diagram between $0.2<z<3.0$ is qualitatively similar to Marchesini et al. (2014) who tracked the progenitors of local ultra-massive ($\log(M_{*}/M_{\odot})\sim 11.8$) galaxies, with the main difference being that this study contain galaxies which are bluer than those of Marchesini et al. (2014). The origin of this difference is rooted in the fact that we select progenitors for a lower local mass galaxy ($\log(M_{*}/M_{\odot})\sim 11.5$). Our galaxies in the higher-$z$ bins are also bluer than the sample of Ownsworth et al. (2016), who select progenitors of massive galaxies based on fixed-cumulative number density. As previously discussed, a fixed cumulative number density selection will yield progenitors which are systematically more massive and thus, redder in $U-V$ and $V-J$ colours, and the inconsistencies in galaxy properties between the samples is likely attributed to differences in stellar mass.

The progression of galaxies between different redshift bins within Fig. 2 and Fig. 3 already provides clues as to the structure of the galaxies within them. Numerous studies find that galaxies in the quiescent region of the $UVJ$-diagram tend to have higher $n$ and smaller $r_{e}$ (e.g., Williams et al., 2010; Patel et al., 2012; Yano et al., 2016). However, those analyses were for galaxies at fixed masses and did not connect progenitor to descendent, and therefore do not make a direct evolutionary link. In the next section we examine the size and structural evolution of the galaxies selected using the cumulative number density method.

4. Evolution in Far-Infrared Star Formation Rates

In Fig. 2 and Fig. 3, we see evidence that the evolution of massive galaxies can be broadly separated into two epochs. At $z>1.5$, galaxies have colours which are consistent with growth mainly through in-situ star formation. At $z<1.5$, galaxy colours are consistent with quenched systems, with mergers becoming the dominant mechanism for growth. We can estimate this epoch more directly by comparing star formation rates to the mass assembly implied from the evolving cumulative number density selection.

Fig. 4 shows the SFR plotted against the derivative of the progenitor mass growth from the right panel of Fig. 1. The SFRs are calculated from far-infrared (FIR) luminosities, which are derived from stacks which include Spitzer $24~{}\micron$, and Herschel PACS and SPIRE bands. For each UltraVISTA stack, FIR stacks were generated in the same manner as described in Schreiber et al. (2015). From Fig. 4, of Schreiber et al. (2015), we see that we do not have sufficient numbers of galaxies at $z>3.5$ with the CANDELS data to expect a FIR detection. Thus we only calculate SFRs out to $z=3.5$. The FIR luminosities were converted to SFRs via the relation from Kennicutt (1998), with a factor of 1.6 correction to convert between the Salpeter IMF used in Kennicutt (1998), to the Chabrier IMF used for the DR3 catalog.

In order to more directly compare the net stellar mass growth as implied from the abundance matching technique to the stellar mass growth from star formation, a $50\%$ conversion factor has been applied to the SFR to account for stellar mass which is lost in outflows from stellar winds (see van Dokkum et al., 2008, 2010). From Fig. 4, we see that SF is able to account for all of the stellar mass growth at $z>1.5$, with little to no contribution from mergers. In contrast, the SFR at $z<1.5$ are insufficient to explain the mass growth, suggesting stellar mass is accreted via mergers.

Between $1.5<z<2.5$, the stellar mass growth predicted from star formation is greater than what is found from the abundance matching techniques by $0.1-0.2~{}dex$. This discrepancy is also seen in model and observation comparisons (see Somerville & Davé 2015; Madau & Dickinson 2014), with potential for the FIR SFRs to be over estimated during this epoch (see Madau & Dickinson 2014 and discussion therein). In spite of this, the FIR SFR support the notion that massive galaxies grow via star-formation until $z\sim 1.5$, where merger driven growth dominates, consistent with the rest-frame $UVJ$ colours, and what is found in the literature (see Sec. 1 and references therein).

5. Analysis

5.1. Stacked Images

For galaxies at $z<3.5$, images were stacked using $48^{\prime\prime}\times 48^{\prime\prime}$ cutouts taken from the UltraVISTA mosaics, which contain both deep and ultra-deep stripes. For each cut-out, SEDs were generated using the ancillary data available in the UltraVISTA and CANDELS source catalogs. These SEDs were used to flag potential active galactic nuclei (AGN), which were removed from the resultant stack. The individual cut-outs were also visually inspected to remove objects which were identified as doubles, or triples (i.e., were not separated by SExtractor; Bertin & Arnouts 1996) or in close proximity to saturated stars to maintain image fidelity. In total, $<4\%$ of the entire sample was discarded.

Cutouts were centered using coordinates taken from the UltraVISTA DR1 (only deep stripes) and DR3 (only ultra-deep stripes) catalogs, with cubic spline interpolation performed for sub-pixel shifting. For galaxies at $z>3.5$, images were stacked using $24^{\prime\prime}\times 24^{\prime\prime}$ cutouts, taken from the 5 CANDELS fields (AEGIS, COSMOS, GOODS-S, GOODS-N and the UDS), with images centred using the coordinates from the 3DHST photometric catalogs (Brammer et al., 2012; Skelton et al., 2014) with sub-pixel shifting also performed using cubic spline interpolation.

From these cutouts, bad-pixel masks were also constructed using SExtractor segmentation maps. These bad-pixel masks were also used to construct a weight-map for the final stack, by summing the bad-pixel masks (in a similar manner to van Dokkum et al. 2010).

For the UltraVISTA stacks, the ultra-deep and deep cutouts were weighted differently in the final stack as the ultra-deep stripes have an exposure time a factor $\sim 10\times$ greater than the deep stripes. The images are weighted by the expected S/N gain, based on the exposure time (i.e. an image with a factor of $\sim 10$ more exposure time, will result in a S/N gain of $\sim 3$). The exact exposures varied between the Y, J, H and Ks bands, with the relative weights between the deep and ultra-deep also changing slightly.

The cutouts were normalized to the sum of the flux contained in the central $1.5^{\prime\prime}\times 1.5^{\prime\prime}$ (corresponding to $10\times 10~{}\mathrm{pixel}$ for UltraVISTA images and $25\times 25~{}\mathrm{pixel}$ for CANDELS images). A weighted sum was performed on the masked cutouts, with the cutouts contained in the ultra-deep stripes given a heavier weight than those in the deep stripes. This summed image was divided by the weight-map to provide the final stack.

For the UltraVISTA stacks, PSFs were generated similarly to the stacked-galaxy images. Stars within a magnitude range were chosen such that the stars had sufficiently high S/N without being saturated ($\approx 16.5$ $Ks$-band magnitude). The stars were treated in the same manner as the stacks of the galaxies (i.e. normalized and averaged). To account for variations in the PSF across the mosaic, 12 different PSFs for each band were generated corresponding to 12 different regions of the mosaic ultra-deep stripes, and 9 for the deep stripes. A final PSF for the relevant band was generated from a weighted average of the 12/9 PSFs, with the weights corresponding to the number of galaxies from each field that went into the making of the stack. Thus, each stack has a uniquely generated PSF.

For the CANDELS F160W stacks, PSFs for each of the 5 fields were taken from the 3DHST-CANDELS data release (Grogin et al., 2011; Koekemoer et al., 2011; Skelton et al., 2014). In a similar manner to the UltraVISTA stacks, the PSF for the relevant band was generated from a weighted average of the PSF from each field, with the weights corresponding to the number of galaxies from each field that contributed to the final stack, thus each F160W stack similarly has a uniquely generated PSF.

Fig. 5 displays the results from the stacking analysis. Each panel contains a $24\times 24^{\prime\prime}$ display of one of the UltraVISTA bands (either Y,J,H, or Ks), except the last two panels which are stacks of the CANDELS F160W data. The UltraVISTA stacks are all displayed at the same colour scale to highlight the differences in background, which increases with increasing $z$. The F160W stacks were plotted at a different colour scale for clarity, as the background is much higher.

In addition to the stacks in Fig. 5, 100 bootstrapped images were also generated for each stack to constrain uncertainties in the structural parameters determination (see Sec. 5.2). Each bootstrapped image also comes with its own unique PSF that reflects the proportion of galaxies from various fields in the same manner as the original stacked images.

5.2. Sersic Profile Fitting

Sersic fitting (Sersic, 1968) of the stacked and bootstrapped images was performed using GALFIT (Peng et al., 2010), with the only constraints imposed on the fits being to restrict the value of the Sersic indices to between $1<n<6$. Fig. 6 shows the best fit $r_{e}$ and $n$ for each band, and each redshift bin, with the uncertainty derived from the $1\sigma$ distribution of the bootstrapped fits. The UltraVISTA derived values are in black, with each symbol corresponding to a different band. The F160W values are indicated with red diamonds, with the last symbol plotted unfilled to mark where we are incomplete. The HWHM of the UltraVISTA and CANDELS PSFs are indicated on the top panel in light-grey and dark-grey regions respectively. As seen from the top panel of Fig. 6, we resolve the stacked images to within an effective radius for UltraVISTA below $z=2$, and the $r_{e}$ is fully resolved for CANDELS in all redshift bins. Additionally, for the redshift bins for which we have stacks for CANDELS and UltraVISTA, the derived sizes and Sersic indices are roughly consistent with one another, suggesting our ground-based structural parameters are reliable.

Absent from Fig. 6 are best-fit values for $r_{e}$ and $n$ below $z=1$ for the F160W band. In these redshift bins, at the mass ranges considered, there were no galaxies present in the catalogue to contribute to a stack (see Table 1). Similarly, best-fit values for the UltraVISTA bands are not present for all redshift bins with the Y, J, H and Ks dropping out at $z=2,2.5,3$ and $3.5$ respectively, due to insufficient signal-to-noise in the resultant stack (and that we are incomplete in UltraVISTA at $z>3.5$).

5.3. Evolution in $r_{e}$

Due to the progression of redshift between the stacks, the $r_{e}$ and $n$ are measured at varying rest-frame wavelengths. In order to measure as closely as possible the same rest-frame wavelength, we have measured how $r_{e}$ and $n$ change with wavelength. In the left panel of Fig. 7 we have plotted $r_{e}$ as a function of rest-frame wavelength. Different colours correspond to different redshift bins, and different symbols demarcate the observed band (with the same symbol convention as in Fig. 6). The desired rest-frame wavelength of $0.5~{}\micron$ was chosen to minimize extrapolation, as well as still be red-ward of the optical-break.

At $z<3$, we have measurements in multiple bands, and find the effective radii decrease with increasing rest-frame wavelength which is consistent with results from previous studies (e.g., Cassata et al., 2011; Kelvin et al., 2012; van der Wel et al., 2014; Lange et al., 2015). However between $2<z<3$, the uncertainties are consistent with little to no evolution in $r_{e}$ with rest-frame wavelength. When considering the evolving properties of the progenitors with redshift, this result is also consistent with the literature. van der Wel et al. (2014) who measured the sizes of galaxies from CANDELS at $0<z<3$, found the size-gradient with rest-frame wavelength was steepest for galaxies at high-mass, and low-redshift, and flatter for low-mass galaxies. As the progenitors decrease in mass with redshift, we expect a flattening of this gradient. The difference in size-gradients is also seen in local populations. Kelvin et al. (2012) found size-gradients to be flatter for late-type galaxies in the GAMA survey. Because we only have measurements in one band for $z>3.5$, and we are dominated by late-type galaxies at high redshift, we have not extrapolated $r_{e}$ between $3.5<z<5.5$ and assume the measurement is representative of the $r_{e}$ at $0.5~{}\micron$. This is assuming that the size-gradient will be flat for low-mass, late-type galaxies at high redshift.

In the right panel of Fig. 7, we have plotted $r_{e}$ at $0.5~{}\micron$ as a function of redshift. It is clear from both panels of Fig. 7 that the $r_{e}$ decreases out to $z=5$ which is consistent with previous results using a diverse set of methods to select progenitors (e.g., van Dokkum et al., 2010; Williams et al., 2010; Damjanov et al., 2011; Mosleh et al., 2011; Oser et al., 2012; Barro et al., 2013; Patel et al., 2013; Straatman et al., 2015; Ownsworth et al., 2016). In spite of the different choice of progenitor selection, below $z=2$ our measurements fall broadly on the same relation found by van Dokkum et al. (2010) (our values are systematically larger, but this is likely a reflection of our slightly higher mass selection). This result is not surprising, and the consistency is reflected in the right panel of Fig. 1, where at $z<2$, the mass of the progenitors chosen using a fixed vs. evolving cumulative number density are within the uncertainties in both the mass function and the semi-analytic models. Although we measure a slightly steeper relation than van Dokkum et al. (2010), it is surprising how well the relation is extrapolated at $z>2$ given that we are selecting galaxies which are distinct in mass from the fixed-cumulative number density selection.

In Fig. 8 we investigate the evolution of the mass-size plane. We have taken the values of $r_{e}$ from the right panel of Fig. 7, and plotted them against their respective progenitor masses, with the highest mass associated with the lowest redshift bin. For comparison purposes, we have over-plotted the mass-size relations from Shen et al. (2003) and Lange et al. (2015) for both early and late type galaxies. For Lange et al. (2015), who investigate the mass-size relations as a function of rest-frame wavelength, we use their $g$-band relations which corresponds most closely to a rest-frame wavelength of $0.5~{}\mathrm{\micron}$. The measured $r_{e}$ from our stacking analysis fall on the SDSS and GAMA mass-size relations for early-type galaxies at $z<0.1$. However for all other redshift bins, our galaxies fall below the local-mass size relation, consistent with van Dokkum et al. (2010). Also plotted in Fig. 8 are simple single (dash-dot line) and double (solid line) power-law fits of the form

[TABLE]

[TABLE]

From Fig. 8, we see that the double power-law is a more appropriate fit for our data, with the parameters $\gamma=2.9\times 10^{-4}$, $\alpha=0.35$, $\beta=2.1\times 10^{3}$ and $\log(M_{0}/M_{\odot})=14.77$. Continuing with the plotting convention in previous figures, our $4.5<z<5.5$ point has been plotted as an open-face symbol to highlight incompleteness issues within that bin. It is interesting to note that this point has not been included in any of the power-law fits, and the fact that it falls on on the extrapolation of the double power-law is not designed.

The evolution in the mass-size relation in Fig. 8 can be broadly separated into two phases. At $z>1.5$, the mass-size evolution is relatively linear (in log-log space), with most points falling along a single power-law. At $z<1.5$, the size growth becomes more efficient, and no longer follows the same single power law as before. This is broadly consistent with patterns we have seen in Fig. 3 and Fig. 7 i.e., that star-formation and mergers are dominating mass and size growth at different epochs, with this changeover occurring at $z\sim 1-2$. Before this time, mass was primarily added via star formation, which has been shown to be ineffective at altering the structure of massive galaxies (Ownsworth et al., 2012). At these redshifts, we see a marked increase in the quiescent fraction of the progenitors. As star formation is no longer an available pathway to mass growth, the growth is dominated by minor mergers, which efficiently increases the $r_{e}$ (see Sec. 1 and references therein).

5.4. Evolution in $n$

Fig. 9 is analogous to Fig. 7, except we investigate how the Sersic index $n$ changes with rest-frame wavelength in place of $r_{e}$. From the left panel of Fig. 9 we see little to no evolution in $n$ with wavelength at any redshift. We have therefore taken an average $n$ weighted by the bootstrapped uncertainty in each band to measure a representative $n$ for each redshift bin. At $z>3$ where we only have one measurement for each stack, the measurement was considered representative. The resulting values are plotted in the right panel of Fig. 9.

Fig. 9 shows a clear downward trend of $n$ with redshift, consistent with previous findings out $z=2$ (e.g., van Dokkum et al., 2010). This trend is also expected given the evolution in the quiescent fraction. Actively star forming galaxies tend to have lower $n$/are more centrally concentrated than their quiescent counterparts (e.g., Lee et al., 2013; Freeman, 1970; Lange et al., 2015; Mortlock et al., 2015), thus at $z>1$ the decrease in $n$ is likely driven by morphological changes between each redshift bin which we also see reflected in the evolution of the mass-size relations (Fig. 8). van Dokkum et al. (2010) also found $n$ to decrease with redshift, although their relation is steeper than the one measured in the current study. However the $n-z$ relation from van Dokkum et al. (2010) was derived from galaxies at $z<2$, where the slopes are comparable, but where we measure systematically higher $n$.

5.5. Mass Assembly

Equipped with measurements of $r_{e}$ and $n$, we can investigate surface-density profiles, and mass assembly as a function of radius. To generate these profiles, we have assumed that the mass-to-light ratio is constant across the profile, and that all the mass can be found within a radius of $75~{}\mathrm{kpc}$. Given these assumptions, and that the integrated mass within $75~{}\mathrm{kpc}$ must equal the total mass found in the right panel of Fig 1 (i.e., the same constraints used in van Dokkum et al. 2010), we have generated stellar-mass density profiles which can be found in Fig. 10.

Fig. 10 shows the Sersic fits using the values of $r_{e}$ and $n$ for each redshift bin in the right panels of Fig. 7 and Fig. 9 respectively. The transition between the solid and dashed lines for each profile marks the point when the error in the background becomes significant. Since many of the values of $r_{e}$ and $n$ are either interpolated, or averaged (see Sec. 5.3 and Sec. 5.4), the profiles from which this transition point was determined were the closest to the rest-frame wavelength of $0.5~{}\micron$ (these are the same bands which are displayed in Fig. 5).

Fig. 10 illustrates that the majority of mass build-up in galaxies since $z=4$ occurs in the outskirts, consistent with previous findings and the inside-out growth paradigm for massive galaxies (van Dokkum et al., 2010; Toft et al., 2012; Bezanson et al., 2013; van de Sande et al., 2013; Belli et al., 2014a; Margalef-Bentabol et al., 2016). It is only at $z=5$ that we begin to see significant growth in the inner regions. Important to note is that as we are incomplete in that redshift bin, we will be biased towards blue, and possibly diskier galaxies which would likely have lower values of $n$. However given the trend of Sersic index with redshift found in Fig. 9, this does not seem to be an unreasonable depiction of the progenitors.

In Fig. 11, we have divided the surface mass density profile for each redshift bin from Fig. 10 by the surface mass density profile at $0.2<z<0.5$. In this way, we are able to trace the fractional mass assembly as a function of radius. At the highest redshift bins, we see the central regions are the first to form, with very little of the stellar mass beyond $3~{}\mathrm{kpc}$ in place at $z\sim 5$. Between $3.0<z<4.5$, we see rapid growth, with the fraction of mass assembled in the inner regions more than doubling. It is also in this redshift interval that a not insignificant fraction of stellar mass is assembled between $3$ and $10~{}\mathrm{kpc}$. We can trace the redshift of formation as a function of radius by tracing the horizontal dashed-line in Fig. 11, which marks the point at which half of the stellar mass was assembled. As you trace from small to large radii, the dashed line crosses different coloured regions, indicating that the interior regions were the first to assemble, with the outer regions assembling at later and later times, indicative of ‘inside-out’ growth.

We can trace this growth quantitatively by considering the total mass in and outside the $3~{}\mathrm{kpc}$ boundary. We have de-projected the surface density profiles of Fig. 10, and separated the mass growth into stellar mass assembly that is within $r<3~{}\mathrm{kpc}$, and exterior to $r>3~{}\mathrm{kpc}$. The total mass assembly is indicated in black, and is the same mass assembly seen in the right panel of Fig. 1. From the red line, we see continuous, albeit decelerating, mass assembly from $z=5$ to $z=0$. This is inconsistent with previous works such as van Dokkum et al. (2010) and Patel et al. (2013) who found the interior regions are consistent with no assembly since $z=2$, oft cited to be evidence of ‘inside-out’ growth, although it depends on precisely what is meant by this term.

It is important to note from Fig. 11 and Fig. 12 that even though the regions outside $3~{}\mathrm{kpc}$ experience a greater growth rate than the inner regions, there is still significant mass build-up from $z=5$ to $z=0$ in the interior. Although the growth between the inner and outer regions is not self-similar, the growth is not necessarily ‘inside-out’ as described in previous works (e.g., van Dokkum et al., 2010; van de Sande et al., 2013), especially when considering the mass assembly at $z>3$. At these redshifts, significant stellar mass is assembled at all radii (although mass accretion is concentrated in the central regions).

5.6. Comparisons with simulations

There have been many comparisons between the mass growth of galaxies in extra-galactic surveys (i.e. mass functions) to hydrodynamical galaxy simulations (e.g., Vogelsberger et al., 2014; Schaye et al., 2015). In fact, the EAGLE simulation has been calibrated to reproduce the galaxy stellar mass function at $z=0$. However, there remain few examples (e.g, Snyder et al., 2015; Wellons et al., 2015) in the literature which explicitly compare the evolution of structure in simulations to observations . In this section, we endeavour to make such a comparison.

In Fig 13, we see how the mass assembly as implied by our observations compares to the EAGLE simulation (Schaye et al., 2015). In Fig. 13, we see the total mass assembly (black line), the mass assembly within a $3~{}\mathrm{kpc}$ aperture and the mass assembly outside a $3~{}\mathrm{kpc}$ aperture. Also plotted in Fig. 13 are the de-projected aperture masses from this study for comparison. The progenitors in EAGLE are defined as the ‘true’ progenitors, and are selected in a similar method to the dark-matter halo mergers trees from Behroozi et al. (2013) which inform the abundance matching technique, i.e. only the most massive progenitor from the precursors of a merger is considered. The progenitors were traced from all galaxies within the EAGLE simulation that have a stellar mass within $0.1~{}\mathrm{dex}$ of $\log(M_{*}/M_{\odot})\sim 11.5$ (i.e. chosen to match the starting point of this study), which amounted to 24 galaxies. The aperture masses from EAGLE quoted above are averages from the progenitors of these 24 galaxies.

A qualitative comparison between the simulations and the observations show remarkable agreement. For the mass within $3~{}\mathrm{kpc}$, the agreement is always within a factor of 2, which is within the uncertainty associated with the assumptions made when determining stellar masses from photometry (Conroy et al., 2009). Both methods predict the same overall trend i.e. that there is a steady build up of stellar mass within $3~{}\mathrm{kpc}$, and rapid assembly at later times at radii larger than $3~{}\mathrm{kpc}$. The main difference between the simulations and observations is that EAGLE predicts a more rapid assembly of the progenitors. The progenitors in EAGLE must assemble more mass in the same period of time in order to result in the final stellar mass of $\mathrm{\log{M_{*}/M_{\odot}}}=11.5$ at $z\sim 0.3$. This offset is not entirely unexpected, given differences between the evolution of the observed and simulated galaxy stellar mass functions at high-$z$ in the mass ranges considered for this study ($\sim 10^{10}-10^{11}~{}\mathrm{M_{\odot}}$, Furlong et al. 2015).

The progenitors in EAGLE must assemble more mass in the same period of time in order to come to the same descendant mass by $z\sim 0.3$; and given the agreement with observations at $r<3~{}\mathrm{kpc}$, nearly all of this mass growth must occur in the outer regions. This suggests the progenitors in EAGLE are more centrally concentrated than observed, except at $z>4$. Between $4<z<5$, the fraction of stellar mass outside a $3~{}\mathrm{kpc}$ aperture is broad agreement with the observations, which does not follow the trend at $z<4$. One possible reason for this is the effective radius at these redshifts is close to $1~{}\mathrm{kpc}$, which suggests nearly all of the total bound mass in the galaxy would be within $3~{}\mathrm{kpc}$, which is not true at lower redshifts.

Some caveats that could affect the above comparison are some assumptions that were made in the observations, in particular the assumption of a constant mass-to-light ratio for our surface mass density profiles. If there is a strong gradient of stellar age with radius in the progenitors, and the interiors are older (which would be consistent with what we see in Fig. 11), then we would over-predict the fraction of the total stellar mass which is located at large radii, bringing us closer to agreement with the simulations. A similar effect would be expected if there are also strong gradients in dust. An analysis of forthcoming virtual observations from EAGLE with the effects of dust and inter-cluster light taken into account would be a better comparison, the investigation of which is beyond the scope of this paper.

6. Discussions and Conclusions

6.1. Mass and size growth at $z<2$

In this paper, we have selected the progenitors of today’s massive galaxies through an evolving cumulative number density technique, and have made image stacks to infer their evolution with redshift. Based on rest-frame $U-V$ and $V-J$ colours, we find the progenitors of massive galaxies become increasingly star forming out to higher redshift, and by assuming Sersic profiles for the mass distribution, we find the progenitors decrease in both $r_{e}$ and $n$. These trends are qualitatively similar to previous studies which select based on fixed (e.g., van Dokkum et al., 2010; Patel et al., 2013; Ownsworth et al., 2014), and evolving (Marchesini et al., 2014) cumulative number densities at $z\lesssim 2$

Although the qualitative trends are consistent with the literature, there are quantitative differences, especially in regards to the evolution of the central mass densities with redshift. Previous works (e.g, van Dokkum et al., 2010; Patel et al., 2013; van Dokkum et al., 2014) have found little to no mass assembly in the inner regions ($r<\mathrm{2~{}kpc}$), and find that mass assembly occurs in an inside-out fashion with the majority of mass growth since $z\sim 2$ occurring at $r>\mathrm{2~{}kpc}$ (although in van Dokkum et al. 2010 at $\sim 1~{}\mathrm{kpc}$, there is a spread in mass-density of at least $0.1~{}\mathrm{dex}$ since $z=2$ suggesting modest mass growth). In this study, we find the central regions have accumulated $\approx 50\%$ of their mass between $2.0<z<2.5$, but continue to experience mass growth out to $z=0.2$, albeit at a lower rate (i.e. we find $\sim 90\%$ of the mass within $2~{}\mathrm{kpc}$ was in place by $z\sim 1$).

The suspected cause of this discrepancy is the differences which arise between a fixed vs. evolving cumulative number density selection. By using a fixed cumulative number density selection, one is biased towards the most massive progenitors (e.g., Clauwens et al., 2016; Wellons & Torrey, 2016). This is a result of the fact that an abundance matching technique (i.e Behroozi et al. 2013) predicts higher number densities with increasing redshift, whereas a fixed cumulative number density will select galaxies at a steeper point in the mass function which is inhabited by higher mass galaxies. We have tested this hypothesis by re-measuring the surface mass density profiles for a fixed cumulative number density selection (see Fig. 1 for the mass assembly history), and do find the redshift evolution in central regions of the stellar surface mass density profiles is considerably weaker than for an evolving number density selection (Fig. 14). Details of this analysis can be found in an attached appendix.

In contrast to a fixed cumulative number density selection, van Dokkum et al. (2014) selected galaxies based on their stellar surface mass density within $1~{}\mathrm{kpc}$ (i.e. ’dense cores’), and found evidence that the interiors are formed first, with the outer radii forming around them. This inconsistency can also be attributed to selection, and the progenitors van Dokkum et al. (2014) select are likely a subpopulation of the progenitors of massive galaxies. Since they are selected on central stellar density, and central stellar density is correlated with quiescence, they will not select star forming progenitors. This is evidenced by the differences in quiescent fraction at $2.0<z<2.5$; van Dokkum et al. (2014) find a quiescent fraction of $57\%$, whereas the selection of the current study has a quiescent fraction of $23\%$ in the same redshift range.

The most massive progenitors are likely to host older stellar populations, have less star formation, and more compact configurations due to rapid early assembly . As these progenitors would have assembled first, they experience more passive evolution in their central regions between $z=2$ and today (e.g., van de Sande et al., 2013). The star forming progenitors however, still must quench, and might involve more violent events, such as disk instabilities which result in compaction, i.e. the driving of mass towards smaller radii (Dekel & Burkert, 2014; Barro et al., 2014). By averaging these populations, one would expect modest gains in stellar mass density in the central regions, which is what is seen in our analysis.

An important caveat to consider when selecting progenitors at systematically higher number densities is the effect of a lower normalization to the mass profiles. Our profiles are designed such that $100\%$ of the stellar mass, as determined from the mass functions as outlined in Fig. 1, is contained within $75~{}\mathrm{kpc}$. If at each $z$ step we have a slightly lower mass selection than studies based on a fixed cumulative number density selection, the normalization of the profile will trend to lower values which imposes sustained mass growth in the central regions (see Fig. 14, and discussion in the appendix).

Although we find the progenitors continue to assemble mass at all radii, the growth rate at small and large radii is not self-similar. The fractional growth rate is higher at larger radii, consistent with the idea that minor mergers play a dominant role in the mass assembly at $z<1.5$, and especially at $z<1$ as found by Newman et al. (2012); Whitaker et al. (2012); Belli et al. (2014b, 2015); Vulcani et al. (2016).

This is also in agreement with our quiescent fractions, which are $>90\%$ at $z<1$, suggesting that the majority of the mass growth cannot be from star formation. However between $1<z<2$ our star-forming fraction exceeds $50\%$, suggesting the increasing importance of star formation in mass assembly, which is in broad agreement with Vulcani et al. (2016) who find star-formation and minor mergers play equal roles in mass growth during this epoch. Additionally, $H{\alpha}$ maps of massive star-forming galaxies between $0.7<z<1.5$ reveal that the disk scale lengths are larger in $H{\alpha}$ than in the stellar continuum, suggesting that star formation also contributes to the mass build-up at large radii (Nelson et al., 2016), and not just in the inner regions.

6.2. Mass and size growth at $z>2$

In addition to comparisons with other works, which are largely limited to $z<2$, we have selected progenitors, and generated stacks for galaxies out to $z=5.5$. In this regime we see a continuation of the trends at $z<2$, i.e. progenitors are smaller, and have Sersic indices which imply more disk-like configurations than spheroidal. This is consistent with the evolution of our quiescent fraction which continues to decrease with increasing $z$, suggesting the progenitors are dominated by star forming galaxies which also tend to have disk-like morphology, which is observed in massive galaxies at high redshift (e.g., van der Wel et al., 2011; Wuyts et al., 2011; Bruce et al., 2012; Newman et al., 2015). This is in agreement with the prediction of Patel et al. (2013) who posited that the progenitors of massive galaxies at $z>3$ will continue the trend towards smaller sizes.

The trends in the evolution of the mass-size relation, the $r_{e}$, the Sersic index, the $UVJ$ colour evolution, and the FIR derived SFRs all corroborate the idea that $z\sim 1.5$ represents a transitional period in how the progenitors of massive galaxies assemble their mass. At $z>1.5$ the $UVJ$ colours and the FIR SFRs suggests the progenitors are actively forming stars, and the Sersic index suggests those stars are consistent with being distributed in an exponential disk. The change in power-law slope at $z\sim 1.5$ in the evolution of the mass-size plane suggests a change in assembly method; one in which the size evolves more efficiently with mass than at higher redshift, consistent with the minor merger scenario (see Sec. 1 and references therein). This is further corroborated by the fact that the FIR SFR is insufficient to account for the rate of stellar mass assembly at $z<1.5$ (Fig. 4).

This study supports the scenario that the progenitors of massive galaxies begin with a disk-like morphology with the disk forming concurrently with the central regions (i.e. the ‘bulge’). At some point, the disk morphology is destroyed, either by major mergers, or disk instabilities which may also be responsible for the increase in quiescent fraction. Evidence of disks (e.g., van Dokkum et al., 2008; van der Wel et al., 2011; Wuyts et al., 2011; Bruce et al., 2012; Bell et al., 2012) and rotation (Newman et al., 2015) in massive compact quenched galaxies are seen at intermediate ($1.5<z<3$) redshifts which confirms that at least some of the massive progenitors host/hosted disk-like morphology. By $z=1.5$, assembly is less violent, with mass growth dominated by minor mergers, and more passive quenching (i.e. gas exhaustion) until $z=0$.

The scenario that the progenitors of massive galaxies begin as disks have support in cosmological simulations. Fiacconi et al. (2016b) simulated the assembly of the main progenitor of a $z=0$ ultra-massive elliptical, and found the progenitor to be disk dominated, with an exponential brightness profile at $z>6$ which had experienced several major mergers at $z>9$. The ‘survival’, or more accurately, the reassembly of the disk after a major merger is feasible, provided the major mergers are sufficiently gas rich (e.g., Hopkins et al., 2009). Fiacconi et al. (2016b) also calculated the Toomre parameter for their simulated disk and found it to be stable against fragmentation for all resolved spatial scales, with the disk supported by a turbulent inter-stellar medium thought to be due to feedback from star-formation. They also predict, that gas-rich star forming disks at $z>5$ should not host a significant bulge, but is rather built up by mergers occurring at $2<z<4$ (Fiacconi et al., 2016a). This is consistent with our analysis which show the majority of the stellar mass in the central regions (i.e. $r<1~{}\mathrm{kpc}$, which we take as a proxy for the bulge) is assembled between $2.0<z<5.5$.

Stacking analysis is a useful tool to probe the average properties of low-surface brightness features of a population of galaxies. However, specific aspects of the morphology are lost in a stack. To verify our hypothesis about the nature of the progenitors of today’s massive galaxies will require resolution and sensitivity of spaced based observatories such as HST. At high-$z$, the rest-frame optical emission is shifted further into the infrared, of which future space observatories such as JWST which will observe at wavelengths beyond the $K$-band, will prove to be invaluable in determining the nature of ’regular’ galaxies at $z>2$.

7. Summary

To briefly summarize the paper, we have traced the stellar mass evolution of the median progenitors of $\log{M_{*}/M_{\odot}}=11.5$ galaxies at $z=0.35$ using abundance matching techniques. Using photometric data from the UltraVISTA and 3DHST surveys and their associated catalogs, we have used stacking analysis to trace the mass assembly of the progenitors out to $z=5.5$. By fitting the images stacks with 2D convolved Sersic profiles, we have found the following.

Selecting progenitors based on an evolving cumulative number density selection results in progenitors that are less massive than if selected based on a fixed cumulative number density selection. This discrepancy becomes significant at $z>2$. 2. 2.

The progenitors of massive galaxies become progressively more star forming, with star forming fractions exceeding 50% at $z>1.5$ as determined by their rest-frame $U-V$ and $V-J$ colours. 3. 3.

The progenitors decrease in both effective radius and Sersic index with increasing redshift, which is consistent with the picture that the progenitors of today’s massive galaxies began with disk-like morphology. 4. 4.

The progenitors continue to assemble mass at all radii until $z=0.35$, which suggests a more complex mass assembly then ‘inside-out’ growth. 5. 5.

Even though galaxies continue to assemble mass in their interiors to low redshift, the redshift at which half of the resultant stellar mass is assembled is higher for the interiors than the exterior regions, with $z_{f,r=3~{}kpc}\sim 2-3$, and $z_{f,r=10~{}kpc}\sim 1-2$. 6. 6.

A brief comparison between the implied mass assembly of this study to results from the EAGLE simulation show a very similar qualitative trend. However the results from simulations imply a more rapid assembly of the outer regions.

8. Acknowledgements

The authors would like to thank the anonymous referee, as well as Matthew Ashby, and Joop Schaye for the helpful comments and discussion which greatly improved this paper. The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Program (FP7/2007-2013)/ERC Grant agreement no. EGGS-278202. This work is 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 also made possible by the observations taken by the 3D-HST Treasury Program (GO 12177 and 12328) with the NASA/ESA HST, which is operated by the Association of University for Research in Astronomy, Inc., under NASA contract NAS5-26555. DM acknowledges the support of the Research Corporation for Science Advancement’s Cottrell Scholarship, and the National Science Foundation under Grant No. 1513473. KIC acknowledges funding from the European Research Council through the award of the Consolidator Grant ID 681627-BUILDUP. This research has made use of NASA’s Astrophysics Data System.

Appendix A The effects of a fixed cumulative number density selection on the stellar surface mass density profiles

In this appendix, we briefly explore the effects of mass selection on the stellar surface mass density profiles. A key finding of this study is that the central ($r<1-2~{}\mathrm{kpc}$) stellar surface mass densities evolve more strongly than observed in earlier works (e.g., van Dokkum et al., 2010; Patel et al., 2013; van Dokkum et al., 2014). It was suspected that this discrepancy was a result of the different number density selections (i.e. a fixed vs. evolving cumulative number density selection as discussed in Sec. 2.1), with a fixed number density selection yielding more massive progenitors (Fig. 1).

To properly investigate this, we repeated our analysis (as detailed in Sec. 5) for a fixed cumulative number density selected sample. In Fig. 14, we plot the resultant surface mass density profiles. A comparison of the right panels of Fig. 10 and Fig. 14 shows the new mass selection significantly alters the observed surface mass density profiles in the central regions. In Fig. 10, we see a difference of $\approx 1~{}\mathrm{dex}$ between the lowest and highest redshift bin at $r<2~{}\mathrm{kpc}$. In contrast, the inner profiles in Fig. 14 lie approximately on top of each other with most mass evolution occurring in the outskirts.

By choosing progenitors using the same methods as previous studies (e.g., van Dokkum et al., 2010; Patel et al., 2013; van Dokkum et al., 2014), we recover their trends, i.e., there is very little redshift evolution in the central stellar surface mass densities and that most mass evolution is occurring in the outskirts ($r>2~{}\mathrm{kpc}$). The effect of selection on the evolution of surface mass-density profiles is two-fold. First, a fixed cumulative number density selection yields higher mass progenitors, which will tend to be more spheroidal, and more centrally concentrated. Secondly, for an evolving cumulative number density selection, the mass evolves more steeply, with less-massive progenitors at high redshift. This will mean the normalization of surface-mass density profiles will also evolve more steeply which is reflected in the evolution of the central stellar surface mass density (as seen in Fig. 10).